# Detecting temporal scaling with modified diffusion entropy analysis

Garland Culbreth111Current affiliation: Institute for Health Metrics and
Evaluation, University of Washington, 3980 15th Ave. NE, Seattle, 98195, WA,
USA. Jacob Baxley  David Lambert

###### Abstract

We present a modification to the diffusion entropy analysis method for
detecting temporal scaling. Diffusion entropy analysis detects temporal
scaling in a data set by converting a time-series into a diffusion trajectory
and using the entropy of that trajectory to measure the temporal scaling in
the data. We modify this by performing an event detection step to construct
the diffusion trajectory. The new modified diffusion entropy analysis offers
substantial improvements over the original method, especially for noisy data.
We describe the method’s purpose, how it works step-by-step, its application,
and future development.

###### keywords:

diffusion entropy , anomalous diffusion , complex systems theory , time-series
analysis , stochastic processes

###### PACS:

02.50.Ey , 05.40.Ca , 05.40.Fb , 05.45.Tp , 05.70.Fh

###### MSC:

60G35 , 60G52 , 60H40 , 60K05 , 62M10 , 62P10 , 62P20 , 62P25

††journal: arXiv

[inst1] organization=Center for Nonlinear Science, University of North Texas,
addressline=1155 Union Cir, city=Denton, postcode=76203, state=TX, country=USA

## Introduction

Diffusion Entropy Analysis (DEA) is a time-series analysis method for
detecting temporal complexity in a dataset; such as heartbeat rhythm (Bohara
et al., 2017; Tuladhar et al., 2018; Jelinek et al., 2020), seismographs (Mega
et al., 2003), or financial markets (Cai et al., 2006). DEA converts the time-
series into a diffusion trajectory and uses the deviation of this diffusion
from that of ordinary Brownian motion as a measure of the temporal complexity
in the data. Several publications describe the method, mostly in broad strokes
and only with words, but none yet have provided a step by step breakdown of
what is happening with figures to illustrate.

This tutorial describes the modified DEA (MDEA) formally set forth by Culbreth
et al. (2019). This is sometimes referred to as DEA with stripes. Stripes
refers to lines or thresholds overlaid atop the time-series data at regularly
spaced intervals. An event is said to occur at any time index where the data
series crosses one of these stripes. MDEA uses the asymmetric diffusion method
described by Grigolini et al. (2001): events are taken to be positive, in our
case $+1$, and are cumulatively summed to generate a diffusion trajectory
based on the time-series. Once this diffusion trajectory is made a moving
window method is applied, as in the original DEA (Scafetta and Grigolini,
2002). The window is moved along the diffusion trajectory and many small
trajectories generated. Their displacements are binned in a histogram to
create a probability density, which is then used to compute the Shannon
entropy. Window length, $l$, is then increased and the process repeated to get
many values of the entropy, $S(l)$, one for each window length. The end goal
is to plot $S(l)$ vs. $ln(l)$ and to use the relation between them to extract
the scaling, $\delta$, of the time-series process.

## Method

To demonstrate the method, a random walk of 10,000 steps is put into the MDEA
program for analysis. Each step in the process is described in words and
illustrated with a figure. Since the full 10,000 step series is difficult to
visually digest, Figures 1–3 display only the first 100 steps and Figures 4–5
the first 1,000 steps. Figure 6 shows the result of MDEA when run over the
full 10,000 step series.

The first step in MDEA is to apply stripes to the data series under analysis.
The researcher must define a number of stripes to use. Rigorous rules for this
choice have yet to be formulated; a good practice is to choose a number
between 2 and 100. The stripes are applied to the data by finding the maximum
and minimum of the data, $x$, and computing the width, $w_{\mathrm{x}}$:

$\mathit{w_{\mathrm{x}}}=\left|\mathrm{max}(\mathit{x})-\mathrm{min}(x)\right|.$
(1)

The size of the stripes, $w_{\mathrm{s}}$, on the data width,
$w_{\mathrm{x}}$, are then computed by taking:

$\mathit{w_{\mathrm{s}}}=\frac{w_{\mathrm{x}}}{\mathit{N}},$ (2)

where $N$ is the number of stripes to be applied. The data is then uniformly
scaled so that the stripes have width 1 according to:

$\mathit{X}=\frac{x}{\mathit{w_{\mathrm{s}}}}.$ (3)

The result of this process is shown in Figure 1.

Figure 1: Original data and rounded data. (a) The original data series. (b)
The corresponding rounded data series, shifted to fit within evenly spaced
intervals according to the number of stripes chosen, with grid lines drawn at
stripe boundaries.

Once the data has been rounded to the stripe widths, events are recorded at
all times when the data crosses, or intersects, a stripe. To determine whether
a crossing occurs at time $t$, MDEA checks two conditions:

$X_{t}\stackrel{{\scriptstyle?}}{{<}}\mathrm{floor}(X_{t-1}+1),$ (4)
$X_{t}\stackrel{{\scriptstyle?}}{{>}}\mathrm{ceil}(X_{t-1}-1).$ (5)

If either of these conditions is false then a crossing must have occurred at
time $t$. If a crossing occurs at time $t$ a 1 is appended to the event array,
$E$, at index $t$; if not a 0 is appended. This process is shown in Figure 2.

Figure 2: Rounded data and events. (a) The rounded data. (b) Events being
recorded at every time index in the rounded data where either Equation 4 or
Equation 5 are false. All events have magnitude 1 and are positive (Grigolini
et al., 2001).

The next step is to construct a diffusion trajectory, $Y$, from the recorded
events. This is done by taking the cumulative sum of the event array, $E$:

$Y_{n}=\sum_{t=0}^{n}E_{t},$ (6)

where $n$ starts off at the first index of the array—0 in most programming
languages—then increases by 1 until $n$ equals the length of the array. This
is depicted in Figure 3.

Figure 3: Events and diffusion trajectory. (a) Event array. (b) Diffusion
trajectory constructed by taking cumulative sum of the events array.

DEA uses a moving window to examine this diffusion trajectory. The process is
similar to some other time-series analysis methods which use moving windows.
In particular, the moving window is used similarly to Detrended Fluctuation
Analysis (DFA) (Peng et al., 1994). The difference is that DFA uses windows
that do not overlap, while DEA uses window positions advancing one time index
by one with much overlap. This similarity and the rationale behind this
decision is discussed in Grigolini et al. (2001). For this demonstration only
one window length, $l=100$, is defined, which is placed at fixed positions
(end to end) on the diffusion trajectory. This makes the windows and snapshots
easier to visualize. In practice a set of window lengths is defined and each
stepped down the trajectory, shifting one index at a time. For this
demonstration only one possible window length is used and only a subset of the
window positions.

Figure 4: Moving window and trajectory slices. (a) The diffusion trajectory
with grid lines drawn at a set of window positions, using $l=100$. (b) The
slices of the diffusion trajectory across each of the window positions in (a),
laid atop each other and shifted to start from the same origin along the
$y$-axis. Here the first 1000 steps in the series are used.

Once these trajectory slices have been computed the next step is to make a
histogram of their displacements. For each window length $l$, stepping the
window along the trajectory generates many slices, each of which goes some
distance from the origin – along the y-axis in Figure 4. These displacements
are binned in a histogram, which is then normalized to turn it into a
probability distribution, $P(l)$. One such histogram is plotted in Figure 5.

Figure 5: Distribution of diffusion slices. Histogram of the displacements of
diffusion trajectory slices; normalized to form a probability distribution,
$P(l)$.

Information entropy, or Shannon entropy (Shannon, 1948), is defined:

$S(Y)=-\sum_{i}P(y_{i})\log(P(y_{i})),$ (7)

where $S$ is the entropy, $Y$ the random variable, $y_{i}$ the possible values
of $Y$, and $P(y_{i})$ the probabilities associated with each value. For DEA,
with $Y$ being the value of the diffusion trajectory at $t=l$ and with
probability distributions $P_{l}$ corresponding to window length $l$, this
goes:

$S_{l}(Y)=-\sum_{i}P_{l}(y_{i})\ln(P_{l}(y_{i})).$ (8)

DEA uses many window lengths, constructs probability distributions for each,
and calculates the entropy corresponding to each according to Equation 8.
Figure 6 illustrates an example of the result.

Figure 6: Principal result of MDEA. (a) Shannon entropy $S(l)$ vs. natural log
of window lengths $\ln(l)$. The choice of natural log of $l$ for axis scale
reveals any extant power-law behavior. (b) The slope of the linear fit,
$\delta$, is the scaling of the time-series process. For a random walk such as
this data, $\delta=0.5$ is the expected value for MDEA to measure.

The goal of DEA is to extract the scaling of the data series. This is done by
using the relation between $S$ and $l$:

$S(l)=\delta\ln(l)+C,$ (9)

where $\delta$ is the scaling of the process and $C$ is some constant
(Scafetta and Grigolini, 2002). This relation is linear, in log-scale.
Therefore, if a linear fit is performed in log-scale, as is shown in Figure 6,
the slope is the scaling, $\delta$, of the time-series.

Clauset et al. (2009) show that using least-squares to fit the slope for a
line of the form in equation 9 introduces systematic errors. MDEA therefore
includes support for using more robust linear fitting: the Siegel (1982) and
Theil (1950)-Sen (1968) methods. These methods fit an ensemble of lines
between pairs of data points and use the median slope of the ensemble as the
slope of the fit line. This makes them much more robust against outliers, and
avoids the systematic errors introduced by least-squares.

## Discussion

MDEA measures a kind of temporal complexity related to criticality. The
scaling, $\delta$, that MDEA outputs provides a measure of the complexity
present in the data set. The behavior of $\delta$ and its relation to
criticality and complexity parameters is discussed in depth in Grigolini et
al. (2001) and Scafetta and Grigolini (2002). For this tutorial it suffices to
understand the basic rules for interpreting the results. For a completely non-
complex process, such as a random walk, MDEA yields $\delta=0.5$. For a
process at criticality, MDEA yields $\delta=1$. Therefore, $\delta$ represents
a measure of the complexity present in the process: the closer $\delta$ is to
1 the closer the process is to criticality.

MDEA offers significant advantages over the original DEA. In DEA the original
time-series is used as the diffusion trajectory with no pre-processing. As a
result several kinds of noise are able to affect the measured scaling. For
example, if the time-series is characterized by fractional Brownian motion
(FBM) (Mandelbrot and Ness, 1968), DEA will measure a scaling equal to the
Hurst index of the FBM. Even in cases where FBM is only a component of the
data and not a complete characterization, the scaling measured by DEA is
significantly affected by the FBM contribution. MDEA filters out any existing
FBM contribution to the scaling, allowing examination and quantification of
other complexity sources in the data (Culbreth et al., 2019).

In particular, MDEA improves detection of events that occur at time intervals
distributed according to a power law. These events appear in several complex
processes, and they are conjectured to be a signature of a complex system
undergoing self-organization in time (Bohara et al., 2018).

## Future Development

MDEA development is ongoing and several key areas remain where MDEA requires
improvement.222As of , the MDEA code repository is located at:
https://github.com/garland-culbreth/Diffusion-Entropy-Analysis Most urgent of
these is the determination of fitting interval. So far this interval has been
chosen on a case-by-case basis by the researcher looking at the plot of $S(l)$
vs. $\ln(l)$ and adjusting the fit to lie along the linear portion. This is
not rigorous and quite unsatisfactory. Using the Siegel (1982) and Theil
(1950)-Sen (1968) robust linear fitting methods removes the need to set a
lower bound for this fitting interval, but still require an upper bound. A
rigorous and standardized method of determining the correct fitting interval
is urgently needed, perhaps by performing iterated fits and adjusting the fit
interval based on the residuals. Alternatively, since equation 9 has a power
law form, it may be possible to use maximum likelihood estimation to obtain
the scaling (Clauset et al., 2009). Using maximum likelihood estimation would
remove the need to set a fit interval at all.

Secondly, the method for determining how many stripes ought to be used
requires improvement. So far this is another case-by-case choice on the part
of the researcher, usually done by testing several possible values and trying
to guess which one is best. The best guidance thus far developed is that the
correct number of stripes is usually between 2 and 70. However, there is
currently no rigorous method for determining the correct number, or even for
judging whether the researcher’s guess is good or bad. Further work is
urgently needed to develop a rigorous method for determining how many stripes
ought to be used, which should follow the direction established by Allegrini
et al. (2002).

Finally, the MDEA code requires further optimization, especially for run-time
speed. The MDEA code hosted on GitHub is written in Python and uses NumPy
(Harris et al., 2020), which supports using Fortran code modules to speed up
computations. Future development should investigate whether this functionality
can be implemented into MDEA, or should develop a bespoke integration using
modules written in other high-performance languages, such as Rust. Development
in this direction would provide substantial runtime improvements.

## References

  * Allegrini et al. (2002) Allegrini, P., Balocchi, R., Chillemi, S., Grigolini, P., Hamilton, P., Maestri, R., Palatella, L., Raffaelli, G., 2002\.  Real event detection and the treatment of congestive heart failure: an efficient technique to help cardiologists to make crucial decisions.  arXiv preprint cond-mat/0209038 .
  * Bohara et al. (2017) Bohara, G., Lambert, D., West, B.J., Grigolini, P., 2017\.  Crucial events, randomness, and multifractality in heartbeats.  Physical Review E 96, 062216\.
  * Bohara et al. (2018) Bohara, G., West, B.J., Grigolini, P., 2018.  Bridging waves and crucial events in the dynamics of the brain.  Frontiers in physiology 9, 1174\.
  * Cai et al. (2006) Cai, S.M., Zhou, P.L., Yang, H.J., Yang, C.X., Wang, B.H., Zhou, T., 2006\.  Diffusion entropy analysis on the scaling behavior of financial markets.  Physica A: Statistical Mechanics and its Applications 367, 337–344.
  * Clauset et al. (2009) Clauset, A., Shalizi, C.R., Newman, M.E., 2009.  Power-law distributions in empirical data.  SIAM review 51, 661–703.
  * Culbreth et al. (2019) Culbreth, G., West, B., Grigolini, P., 2019.  Entropic approach to the detection of crucial events.  Entropy 21, 178\.  URL: https://doi.org/10.3390/e21020178, doi:10.3390/e21020178.
  * Grigolini et al. (2001) Grigolini, P., Palatella, L., Raffelli, G., 2001.  Asymmetric anomalous diffusion: an efficient way to detect memory in time series.  Fractals 09, 439–449.  URL: https://doi.org/10.1142/s0218348x01000865, doi:10.1142/s0218348x01000865.
  * Harris et al. (2020) Harris, C.R., Millman, K.J., van der Walt, S.J., Gommers, R., Virtanen, P., Cournapeau, D., Wieser, E., Taylor, J., Berg, S., Smith, N.J., et al., 2020\.  Array programming with numpy.  Nature 585, 357–362.
  * Jelinek et al. (2020) Jelinek, H.F., Tuladhar, R., Culbreth, G., Bohara, G., Cornforth, D., West, B.J., Grigolini, P., 2020.  Diffusion entropy vs. multiscale and rényi entropy to detect progression of autonomic neuropathy.  Frontiers in Physiology 11.
  * Mandelbrot and Ness (1968) Mandelbrot, B.B., Ness, J.W.V., 1968\.  Fractional brownian motions, fractional noises and applications.  SIAM Review 10, 422–437.  URL: https://doi.org/10.1137/1010093, doi:10.1137/1010093.
  * Mega et al. (2003) Mega, M.S., Allegrini, P., Grigolini, P., Latora, V., Palatella, L., Rapisarda, A., Vinciguerra, S., 2003.  Power-law time distribution of large earthquakes.  Physical Review Letters 90, 188501\.
  * Peng et al. (1994) Peng, C.K., Buldyrev, S.V., Havlin, S., Simons, M., Stanley, H.E., Goldberger, A.L., 1994\.  Mosaic organization of dna nucleotides.  Physical review e 49, 1685\.
  * Scafetta and Grigolini (2002) Scafetta, N., Grigolini, P., 2002\.  Scaling detection in time series: Diffusion entropy analysis.  Phys. Rev. E 66, 036130\.  URL: https://link.aps.org/doi/10.1103/PhysRevE.66.036130, doi:10.1103/PhysRevE.66.036130.
  * Sen (1968) Sen, P.K., 1968.  Estimates of the regression coefficient based on kendall’s tau.  Journal of the American statistical association 63, 1379–1389.
  * Shannon (1948) Shannon, C.E., 1948.  A mathematical theory of communication.  Bell System Technical Journal 27, 379–423.  URL: https://doi.org/10.1002/j.1538-7305.1948.tb01338.x, doi:10.1002/j.1538-7305.1948.tb01338.x.
  * Siegel (1982) Siegel, A.F., 1982.  Robust regression using repeated medians.  Biometrika 69, 242–244.
  * Theil (1950) Theil, H., 1950.  A rank-invariant method of linear and polynomial regression analysis.  Indagationes mathematicae 12, 173\.
  * Tuladhar et al. (2018) Tuladhar, R., Bohara, G., Grigolini, P., West, B.J., 2018\.  Meditation-induced coherence and crucial events.  Frontiers in physiology 9, 626\.


# Multiwavelength Opportunities and Challenges in the Era of Public Fermi Data

D. J. Thompson  NASA Goddard Space Flight Center, Greenbelt, MD 20771 USA
on behalf of the Fermi Large Area Telescope Collaboration

###### Abstract

The gamma-ray survey of the sky by the Fermi Gamma-ray Space Telescope offers
both opportunities and challenges for multiwavelength and multi-messenger
studies. Gamma-ray bursts, pulsars, binary sources, flaring Active Galactic
Nuclei, and Galactic transient sources are all phenomena that can best be
studied with a wide variety of instruments simultaneously or
contemporaneously. Identification of newly-discovered gamma-ray sources is
largely a multiwavelength effort. From the gamma-ray side, a principal
challenge is the latency from the time of an astrophysical event to the
recognition of this event in the data. Obtaining quick and complete
multiwavelength coverage of gamma-ray sources can be difficult both in terms
of logistics and in terms of generating scientific interest. The Fermi LAT
team continues to welcome cooperative efforts aimed at maximizing the
scientific return from the mission through multiwavelength studies.

##  I Opportunities

During its first year, the Fermi Large Area Telescope has excelled in
producing scientific results using multiwavelength approaches. Some examples
include:

  * •

PSR J1741-2054 is a radio pulsar found based on gamma-ray timing PSR J1741 . A
bright Fermi LAT point source was the first step. An analysis of the LAT
timing discovered gamma-ray pulsations. A follow-up observation with the Swift
X-Ray Telescope (XRT) found an X-ray source that gave better position
information than could be determined from the LAT image. Using the LAT timing
information and the Swift location allowed archival analysis using Parkes
radio data and a deep search using the Green Bank Telescope that found the
radio pulsar.

  * •

PMN J0948+0022 is known as a narrow-line quasar or a radio-loud Narrow-Line
Seyfert 1 galaxy, a somewhat different class than the blazars that are
regularly seen in gamma rays. Contemporaneous observations combining the LAT
data with Swift (X-ray, UV, and optical) and Effelsberg (radio) revealed a
Spectral Energy Distribution that showed this source to be similar to a
blazar, indicating the presence of a relativistic jet NLS1 .

  * •

Using the public light curves made available by the LAT team (at the Fermi
Science Support Center Web site http://fermi.gsfc.nasa.gov/ssc/), Bonning et
al. 3C454.3 MW  studied simultaneous multiwavelength variability of blazar
3C454.3 using Small and Moderate Aperture Research Telescope System (SMARTS)
telescopes for optical and ultraviolet and X-ray data from the Swift
satellite. They found excellent correlation, with a time lag less than a day,
an important parameter for modeling this blazar in terms of an external
Compton model.

These few examples illustrate some of the ways Fermi results enhance
scientific understanding when combined with observations and analysis from
other wavelengths. Because the position uncertainties for gamma-ray sources
are still large compared to those at many other wavelengths, unidentified
gamma-ray sources are inherently subjects for multiwavelength studies,
depending on timing, spectral, and modeling to determine what objects produce
the gamma rays.

##  II Enabling Technologies

Multiwavelength opportunities like those described in the previous section
were not possible even a few years ago. Several developments have facilitated
such multiwavelength efforts:

  * •

Communication - The ubiquity of network connectivity has allowed rapid
exchange of data and ideas. Wireless Internet access and portable devices of
all sorts have accelerated the exchange of information. Campaigns that once
had to be organized by telephone and letter can now be arranged in a matter of
minutes or hours.

  * •

Facilities - Most parts of the electromagnetic spectrum (and several multi-
messenger fields) are now covered by ground-based and space-based
observatories. Fermi is just one of many facilities that produce prompt and
public results that can be used for multiwavelength study.

  * •

Consolidated Information Centers - Resources like ADS, NED, Simbad, ASDC,
HEASARC, and others facilitate rapid discoveries of existing coverage of
sources. Scientists can now almost instantaneously review archival results for
nearly any cataloged object.

##  III Challenges

Despite the tools and resources now available for multiwavelength studies, the
Fermi LAT presents some challenges in terms of making the best scientific use
of the gamma-ray data. Three of these issues are described in the sections
below.

###  III.1 Challenge 1:Time-Criticality of Response

Gamma-ray bursts (GRB), thanks to the the Gamma-ray bursts Coordinates Network
(GCN) http://gcn.gsfc.nasa.gov/, offer a paradigm for rapid response to
transient astrophysical events. The success of the GCN originates in part from
the intensity of GRB, which can be recognized automatically with high
confidence in satellite detectors. The time for disseminating initial
information about a burst is not generally governed by the actual detection
but rather by the speed of communication. Multiwavelength studies can often
begin within seconds of the initial detection.

The situation for other gamma-ray sources is more complicated, because none of
them approaches the instantaneous brightness of a GRB. Nevertheless, dramatic
changes in flux have been see on time scales of a day or less (e.g. PKS
1502+106, 1502 ). On-board analysis of such sources by the LAT is not
practical, and so the response depends on both communication and analysis. The
Fermi LAT data are stored onboard and transmitted through the Tracking and
Data Relay System to the ground in batches, not continuously. This process
introduces some delay, as does the data reduction to extract the gamma-ray
event candidates and compare the gamma-ray sky with its previous appearance.
Half a day can pass before a flare is discovered.

The LAT team has taken several steps to minimize the latency in reporting
events of multiwavelength interest:

  * •

Automation - Much of the data handling process is now automated, and efforts
continue to streamline the procedures. The analysis pipeline now produces
preliminary flux values for over 40 sources of interest routinely, and these
results are posted daily at
http://fermi.gsfc.nasa.gov/ssc/data/access/lat/msl_lc/.

  * •

Dedication - The LAT team has a group of scientists called Flare Advocates or
Skywatchers, who examine the automatically produced analysis results as soon
as they appear. By applying scientific expertise at this early stage, the
process optimizes the response to findings of astrophysical interest while
minimizing any reaction to statistical fluctuations of steady sources.

  * •

Communication - Flare advocates use three avenues to share results quickly
about activity in the gamma-ray sky. The first is the use of Astronomer’s
Telegrams http://www.astronomerstelegram.org/, over 40 of which have been
issued by the LAT team for quickest reporting of results. The second is a
multiwavelength mailing list, gammamw
https://lists.nasa.gov/mailman/listinfo/gammamw, which is used to contact
scientists directly about gamma-ray multiwavelength news. Anyone interested is
welcome to join this list. The third approach is the Fermi Sky Blog,
http://fermisky.blogspot.com/, which posts weekly summaries of the brightest
sources in the high-energy gamma-ray sky.

###  III.2 Challenge 2: Finding Enough Multiwavelength Coverage

Although the LAT is an all-sky, every-day monitor for high energy gamma rays,
most telescopes at other wavelengths have much smaller fields of view and sky
coverage. In addition, many telescopes have sun-angle constraints.
Multiwavelength coverage of an active gamma-ray source is not assured. Two
approaches are being used by the multiwavelength community to enhance the
coverage of the sky:

  * •

More all-sky or wide-field monitors are becoming available. The RXTE All-Sky
Monitor in X-rays http://xte.mit.edu/ has recently been complemented with the
Japanese MAXI all-sky X-ray monitor on the International Space Station
http://maxi.riken.jp/top/. In optical, the Palomar Quest program regularly
surveys a large area http://www.astro.caltech.edu/~george/pq/, and the Pan-
STARRS http://pan-starrs.ifa.hawaii.edu/public/ and Skymapper
http://msowww.anu.edu.au/skymapper/ programs will be surveying much of the
northern and southern hemispheres repeatedly. At longer wavelengths, Planck
http://www.rssd.esa.int/index.php?project=Planck, Herschel
http://herschel.esac.esa.int/, and WISE
http://www.nasa.gov/mission_pages/WISE/main/index.html are viewing the sky
with fairly long cadence.

  * •

Source monitoring programs have also emerged. In radio, many observers are
cooperating to provide multiple-insrument monitoring of many candidate gamma-
ray targets, particularly blazars. A summary of ongoing activity can be found
at http://pulsar.sternwarte.uni-erlangen.de/radiogamma/. Similarly, optical
programs like the one at the University of Arizona
(http://james.as.arizona.edu/~psmith/Fermi/), SMARTS
(http://www.astro.yale.edu/smarts/glast/), and the GLAST-AGILE Support Program
(GASP, http://www.to.astro.it/blazars/webt/gasp/homepage.html) observe many
gamma-ray sources regularly in the optical.

A useful collection of links to multiwavelength information can be found at
https://confluence.slac.stanford.edu/display/GLAMCOG/. The LAT team greatly
appreciates the ongoing cooperative activities of all these groups and
welcomes other telescope teams who participate in particular campaigns.

###  III.3 Challenge 3: Deciding When to Work with the LAT Team

With all the Fermi data public, along with software for analysis, anyone can
undertake multiwavelength studies incorporating gamma-ray results. There may
be times when contacting the LAT team could benefit such analysis, however.

Analysis of the LAT data does involve some important caveats (see
http://fermi.gsfc.nasa.gov/ssc/data/analysis/LAT_caveats.html for more
details):

  * •

The diffuse Galactic emission is bright and highly structured. The diffuse
model supplied by the LAT team has recently been updated and is likely to
continue to evolve. Separating weaker sources from the diffuse Galactic
emission is non-trivial. There are regions of the sky where the diffuse model
has deficiencies.

  * •

The LAT Instrument Response Functions (IRFs) have significant uncertainties at
energies near 100 MeV and a non-negligible charged particle background at
energies above 10 GeV. Improvements in the IRFs are expected but are not
imminent. Analysis of data below 100 MeV with the current IRFs is not
recommended

Some suggestions about when consulting the LAT team might be beneficial:

  * •

If you are searching for a source that is not in the LAT catalog, then it is
probably weak enough that a simple analysis will not be adequate.

  * •

If you need a detailed energy spectrum or are looking for particular spectral
features, especially at very low or very high energies, the LAT team has
experience with non-standard analysis.

  * •

If you are trying to analyze the Galactic Center region, you are strongly
advised not to go it alone!

  * •

If you are interested in the most complete multiwavelength coverage, consider
contacting the LAT team. The LAT team knows many cooperating groups across the
spectrum who may be interested in working with you (even if you do not include
the LAT team).

###### Acknowledgements.

The Fermi LAT Collaboration acknowledges generous ongoing support from a
number of agencies and institutes that have supported both the development and
the operation of the LAT as well as scientific data analysis. These include
the National Aeronautics and Space Administration and the Department of Energy
in the United States, the Commissariat à l’Energie Atomique and the Centre
National de la Recherche Scientifique / Institut National de Physique
Nucléaire et de Physique des Particules in France, the Agenzia Spaziale
Italiana and the Istituto Nazionale di Fisica Nucleare in Italy, the Ministry
of Education, Culture, Sports, Science and Technology (MEXT), High Energy
Accelerator Research Organization (KEK) and Japan Aerospace Exploration Agency
(JAXA) in Japan, and the K. A. Wallenberg Foundation, the Swedish Research
Council and the Swedish National Space Board in Sweden. Additional support for
science analysis during the operations phase is gratefully acknowledged from
the Istituto Nazionale di Astrofisica in Italy and the Centre National
d’Études Spatiales in France.

## References

  * (1) F. Camilo et al., “Radio Detection of LAT PSRs J1741-2054 and J2032+4127: No Longer Just Gamma-ray Pulsars”, ApJ 705, 1, 2009.
  * (2) E. W. Bonning et al., “Correlated Variability in the Blazar 3C 454.3”, ApJ 697, L81, 2009.
  * (3) A. A. Abdo et al., “Fermi/Large Area Telescope Discovery of Gamma-Ray Emission from a Relativistic Jet in the Narrow-Line Quasar PMN J0948+0022”, ApJ 699, 976, 2009.
  * (4) A. A. Abdo et al., “PKS 1502+106: A New and Distant Gamma-ray Blazar in Outburst Discovered by the Fermi Large Area Telescope”, ApJ in press, 2010.


# Microscopic structure of electromagnetic whistler wave damping by kinetic
mechanisms in hot magnetized Vlasov plasmas

Anjan Paul  Devendra Sharma  Institute for Plasma Research, Bhat, Gandhinagar,
India, 382428  Homi Bhabha National Institute, Training School Complex,
Anushaktinagar, Mumbai 400094, India

###### Abstract

The kinetic damping mechanism of low frequency transverse perturbations
propagating parallel to the magnetic field in a magnetized warm electron
plasma is simulated by means of electromagnetic (EM) Vlasov simulations. The
short-time-scale damping of the electron magnetohydrodynamic whistler
perturbations and underlying physics of finite electron temperature effect on
its real frequency are recovered rather deterministically, and analyzed. The
damping arises from an interplay between a global (prevailing over entire
phase-space) and the more familiar resonant-electron-specific kinetic damping
mechanisms, both of which preserve entropy but operate distinctly by leaving
their characteristic signatures on an initially coherent finite amplitude
modification of the warm electron equilibrium distribution. The net damping
results from a deterministic thermalization, or phase-mixing process, largely
supplementing the resonant acceleration of electrons at shorter time scales,
relevant to short-lived turbulent EM fluctuations. A kinetic model for the
evolving initial transverse EM perturbation is presented and applied to
signatures of the whistler wave phase-mixing process in simulations.

##  I Introduction

The electromagnetic turbulence prevails in collective excitations of charged
matter interacting with, and by means of, the electromagnetic field over a
vast range of spatiotemporal scales, usually terminated by dissipation at the
finer scales. The solar-wind spectrum, for example, shows that beyond a
frequency breakpoint a deviation exists from the inertial rage characterized
by exponent -5/3 of power law Borovsky (2012); Podesta  _et al._ (2007). In
one of the plousible scenarios, the whistler fluctuations can be the
fundamental mode and central means of dissipation in this weak turbulence
regime Gary  _et al._ (2010). A steepening present in the spectrum leading to
considerably high spectral exponents ($\sim$2-3) suggests presence of
considerable damping alongside to the intra-spectral energy transfer Gary  _et
al._ (2010). Besides by conversion to electrostatic modes Xu  _et al._ (2017),
damping by kinetic transverse wave-particle interaction must operate on the
short lived excitations Helliwell (1965); Xiao  _et al._ (2017); Thorne
(2010); Chen  _et al._ (2013) initiated by sponteneous field fluctuations.
Fresh perturbations, so triggered, excite warm plasma eigenmodes by leaving
long lasting remnants of their initially enforced phase-space structure in the
memory of nonthermal kinetic distributions Krall and Trivelpiece (1986). The
asymptotic long-time solutions of the collisionless kinetic formulation Landau
(1944, 1946) applied to them therefore have large scope of sophistication by
admitting a strong deterministic thermalization, or phase-mixing Krall and
Trivelpiece (1986), alongside the damping evaluated in usual time asymptotic,
$t\rightarrow 0$ limit.

The general kinetic evolution produced collisionless damping of
electromagnetic fluctuations Lutomirski (1970) involves a rather complicated
phase-space dynamics and is most accessible by deterministic Vlasov
simulations Fijalkow (1999); Mangeney  _et al._ (2002). Only a limited number
of studies have rather deterministically simulated the dynamics of the
transverse electromagnetic excitations and their damping/stability in a hot
collisionless magnetized plasma Palodhi  _et al._ (2010); Valentini  _et al._
(2007), even as the process remains critical for determining the operational
state of turbulence and the transport associated with it both in space plasmas
Gary  _et al._ (2010); Shuster  _et al._ (2021) as well as in modern magnetic
confinement fusion experiments Fülöp  _et al._ (2009).

In the collisionless limit, the modifications made to temperature, or width,
of an initially equilibrium warm electron velocity distribution produce a
higher order correction to the resonant wave particle interaction process. The
analytical model predicts a related downward shift in the whistler wave
frequency in collisionless plasmas with hotter electrons Chen  _et al._
(2013); Schreiner  _et al._ (2017); Schreiner (2016). The recovered strength
of damping due to wave particle interaction however remains underestimated in
comparison with that produced by the computer simulations implemented with
reasonably low collisionality.

This paper addresses above aspects of kinetic whistler damping mechanism,
subsequent to the recovery of general electromagnetic modes of a magnetized
plasma and their dispersive characterization in our flux-balance based Mandal
and Sharma (2016); Fijalkow (1999) Vlasov simulations. This is accompanied by
illustration of its detailed phase-spatiotemporal evolution. The interaction
of electromagnetic modes, propagating parallel to the magnetic field with the
resonant particles is studied, particularly recovering the damping of the
whistler waves via full kinetic mechanism and comparison of the simulation
results with those analytically prescribed in the linear Landau theory limit
Landau (1944); Krall and Trivelpiece (1986). Presented simulations and
analysis enter the finer regime of kinetic phase-mixing of the electromagnetic
mode uniquely achievable by Vlasov simulations. We qualitatively recover the
phase-mixing effects showing the dominance of frequency ${\bf v\cdot k}$ of
the ballistic term ($\propto\exp(ikvt)$) Krall and Trivelpiece (1986)
accounting for short time response, in addition to time asymptotic Landau
damping results that are routinely applied to turbulent electromagnetic
fluctuations of sufficiently short life time. First quantitative analysis of
the phase-spatiotemporal evolution of a 4D electron phase-space distribution
perturbation associated with the transverse electromagnetic whistler mode
simulated in a hot magnetized Vlasov plasma is presented.

The present paper is organized as follows. In Sec. II the electromagnetic
Vlasov-Maxwell model considered for the simulation is presented. Vlasov
simulation implementation and its characterization by dispersion of transverse
EM perturbation is presented in Sec. III. In Sec. IV the simulations of Landau
damping of transverse whistler wave perturbations in a warm electron plasma is
presented followed by their comparison with analytical results. Based on the
evolution of simulated perturbations, the signatures of phase-mixing
supplemented landau damping of the initial perturbation are presented in Sec.
V. The kinetic model for initial transverse EM perturbations is presented and
solved in Sec. V.2. Conclusions are presented in Sec. VI.

##  II The electromagnetic Vlasov plasma model

The electromagnetic magnetized plasma modes simulated in this paper follow
pure kinetic formulation and are well represented by the solutions of
collisionless, fully nonlinear Vlasov equation for species $\alpha$,

$\displaystyle\frac{\partial f_{\alpha}}{\partial t}+{\bf
v}\cdot\frac{\partial f_{\alpha}}{\partial{\bf
x}}+\frac{q_{\alpha}}{m_{\alpha}}\left({\bf E}+\frac{{\bf v}\times{\bf
B}}{c}\right)\frac{\partial f_{\alpha}}{\partial{\bf v}}=0.$ (1)

The evolution of electric field ${\bf E}$ and magnetic field ${\bf B}$ for the
electromagnetic processes follows the Maxwell’s equations,

$\displaystyle\nabla\times{\bf B}$ $\displaystyle=$
$\displaystyle\frac{1}{c}\frac{\partial{\bf E}}{\partial t}+\frac{4\pi}{c}{\bf
J},$ (2) $\displaystyle\nabla\times{\bf E}$ $\displaystyle=$
$\displaystyle-\frac{1}{c}\frac{\partial{\bf B}}{\partial t},$ (3)
$\displaystyle\nabla\cdot{\bf B}$ $\displaystyle=$ $\displaystyle 0$ (4)

where the current density ${\bf J}$ is described by

$\displaystyle{\bf J}=\sum_{\alpha}\int_{-\infty}^{\infty}d{\bf
v}~{}q_{\alpha}{\bf v}f_{\alpha}$ (5)

Consistent with the electron magnetohydrodynamic regime excitations, we use an
externally applied constant magnetic field ${\bf B}_{0}$ and infinitely
massive ions. The subscript for species $\alpha$ therefore only assumes value,
$e$, representing electrons (and henceforth omitted), which are the only
mobile species and contributing to the perturbation.

The full nonlinear kinetic model (1)-(5) is implemented for the case of waves
propagating parallel to an applied magnetic field (${\bf B_{0}\|k}$) where
principle modes are high and low frequency right and left handed circularly
polarized modes. This set up includes the whistler waves that are right handed
polarized and propagate in the frequency range
$\Omega_{ci}<\omega<\Omega_{ce}$, where $\Omega_{c\alpha}$ is the
gyrofrequency of the species $\alpha$.

##  III Kinetic simulations of electromagnetic waves in magnetized warm Vlasov
plasma

The simulations presented in this analysis progress by numerically evolving
the magnetized plasma (electron) velocity distribution function $f$ according
to the dynamics of the phase-fluid flow Krall and Trivelpiece (1986) which is
governed by the collisionless Vlasov equation (1), and associated Maxwell’s
equations (2)-(4). The 4-dimensional (4D $\equiv$ 1x-3v) phase-space
simulations are performed using an advanced flux balance technique Mandal  _et
al._ (2020); Mandal and Sharma (2016); Fijalkow (1999) generalized to simulate
the electromagnetic plasma modes in a large range of magnetization of plasma
species. The results of simulations are characterized first against the
analytical dispersion relation of left and right handed circularly polarized
low and high frequency modes and then against both analytic and numerical
evaluation of the linear Landau damping descriptions Krall and Trivelpiece
(1986); Gary (1993); Chen  _et al._ (2013).

###  III.1 The simulation set-up

In order to simulate the waves propagating parallel to an ambient magnetic
field ${\bf B}_{0}$ in a wide range of frequency $\omega$ and wave vector
${\bf k}$, we have assumed a setup where both the ${\bf B}_{0}$ and ${\bf k}$
are aligned to $z$-axis and the periodic boundary condition is used at both
the boundaries of the one-dimensional simulation zone located between $z=0$
and $L$. The setup therefore assumes symmetry along both $\hat{x}$ and
$\hat{y}$ directions with finite spread of electron velocities along these
dimensions, besides the direction $\hat{z}$. The three dimensional equilibrium
velocity distribution for the electron species is considered to be a
Maxwellian,

$\displaystyle f(z,{\bf v})=\left(\frac{1}{2\pi
v_{th}^{2}}\right)^{3/2}\prod_{j=x,y,z}\exp{\left\\{-\frac{(v_{j}-\langle
v_{j}\rangle)^{2}}{2v_{th}^{2}}\right\\}},$ (6)

where $v_{th}=(T_{e}/m_{e})^{1/2}$ is electron thermal velocity, $T_{e}$ is
electron temperature in energy units and $m_{e}$ is electron mass.

At time $t=0$, the equilibrium distribution is perturbed with sinusoidal
perturbations (having variation along $\hat{z}$) with initial amplitudes of
$\langle{\bf v}\rangle_{1}$, ${\bf E}_{1}$ and ${\bf B}_{1}$ ($\propto
exp{(ikz)}$) consistent with transverse electromagnetic (right or left handed
circularly polarized) linear magnetized plasma modes propagating with a
desired $k$ Chen (2012); Helliwell (1965).

###  III.2 Dispersion characterization of electromagnetic plasma modes

We present the dispersion relations recovered for the case of very small
electron thermal velocity ($v_{th}=0.001$ c) simulated and its comparison with
the corresponding analytical cold plasma dispersion relation Stix (1992); Chen
(2012).

Figure 1:  Dispersion of simulated frequency (dots) compared to the Right and
Left handed polarized branches of the analytical dispersion relation (solid
line).

Specific to simulation cases presented in this section, we have used a
3-dimensional velocity space grid of rather moderate size having $32\times
32\times 32$ grid points, in combination with spatial grid size also of $32$
grid points. In Fig. 1, the comparison is presented of the simulated values of
frequency $\omega$ of the perturbation plotted as function of $k$, with the
right and left handed circularly polarized (RHCP and LHCP) branches of the
analytical dispersion relation in the limit of infinitely massive ions
($\omega_{\rm pe}\gg\omega_{\rm pi}$),

$\displaystyle
k^{2}_{R,L}=\frac{\omega^{2}}{c^{2}}\left[1-\frac{\omega_{pe}^{2}}{\omega(\omega\mp\Omega_{ce})}\right].$
(7)

Considering the parameters, $v_{th}=0.001$ c, the ratio of electron cyclotron
frequency and plasma frequency $\Omega_{ce}/\omega_{pe}=0.1$ and sufficiently
small initial perturbation exclusively in electron average velocity amplitude
$\Delta V=\langle v\rangle_{\rm 1max}-\langle v\rangle_{\rm 1min}=2\times
10^{-3}c$, the high frequency ($\omega>\Omega_{ce}$) RHCP and LHCP mode phase
velocities are sufficiently high for these waves to stay out of resonance with
the cold electrons chosen for this case. Moreover, for the low frequency
whistler waves excited (dots on the red curve in the region
$0<\omega<\Omega_{ce}$ in Fig. 1), the electron thermal velocity
$v_{th}=0.001c$ chosen is still sufficiently low for no significant resonant
electron population to be available at the resonant velocity, $v_{z}=v_{\rm
res}=(\omega-\Omega_{ce})/k$. In rest of this analysis we exclusively
characterize this low frequency whistler branch of the perturbation for
relatively warmer electrons in order to analyze its resonant damping.

##  IV Landau damping of the transverse Whistler mode

An advanced simulation set up, with grid size of $64\times 64\times 64\times
64$ is implemented in the following sets of simulations, by adopting the
electron velocity range [-0.26, 0.26] c. The velocity distribution of
electrons as function of the parallel velocity $v_{\|}\equiv v_{z}$ is plotted
in Fig. 2 at $v_{x}=v_{y}=0$ for a range of electron thermal velocity
$v_{th}=0.001c$ (inner most profile) to $0.026c$ (outer most profile) explored
in the cases presented in this section. The resonant velocity $v_{\rm
res}=(\omega-\Omega_{ce})/k$ and phase velocity $v_{p}$ of the low frequency
whistler mode for parameters $k=1.0\omega_{p}/c$ and
$\Omega_{ce}/\omega_{p}=0.1$, are indicated by vertical dotted line. The
profiles for smaller values of $v_{th}$ show the population of resonant
electrons drops to negligible values such that no resonant damping is present
for these cases.

Figure 2:  The electron distribution function (normalized to its maximum
value) plotted as function of parallel velocity $v_{\|}\equiv v_{z}$ at
$v_{x}=v_{y}=0$. The thermal velocity of electrons for the curves with
increasing width ranges from $v_{th}=0.001c$ (inner most profile) to $0.026c$
(outer most profile) , respectively. The vertical dotted lines indicate
resonant velocity $v_{\rm res}$ and phase velocity $v_{p}$ of the low
frequency whistler mode at $k=1.0\omega_{p}/c$.

In Fig. 3(a), the time evolution of the amplitude of the velocity perturbation
$\Delta V$ is plotted in the simulations done for the range of electron
thermal velocity $v_{th}=0.001c$ to $0.026c$ for the value of
$k=0.8\omega_{pe}/c$. Two additional set of simulations done for the values of
$k=0.897$ and $k=1.0\omega_{pe}/c$ are presented in Fig. 3(b) and Fig. 3(c),
respectively.

Figure 3:  The decay in the velocity perturbation $\Delta v$ for
$v_{th}=0.001c$ (lowest damping) to $0.026c$ (highest damping) for (a)
$k=0.8\omega_{pe}/c$, (b) $k=0.897\omega_{pe}/c$ and (c) $k=1.0\omega_{pe}/c$.
The modulation seen in some cases are residual RHCP and LHCP excitations
($\omega_{R,L}>\Omega_{ce}$ satisfying (18)).

The time evolution of $\Delta V$ in all the cases above cover the initial
evolution of the wave amplitude only for the time duration $\Delta t\sim
3\omega_{p}^{-1}$ (i.e., about a fraction of one complete cycle of the
whistler wave with frequency $\omega<\Omega_{ce}$ while
$\Omega_{ce}/\omega_{p}=0.1$) which is the time duration in which the linear
Landau damping rate remains reasonable estimate for the wave damping. The
short time evolution in Fig. 3 has sufficiently low numerical widening of $f$
for resolving the effect of $T_{e}$ variation, which is varied with relatively
much larger increments in this study.

###  IV.1 Comparison of simulations with analytic whistler damping rates

The damping rates of whistler velocity perturbation amplitude can be compared
with the available analytic approximations of the Landau damping rate of
whistlers. For the purpose of comparison we have used the analytic results
from kinetic formulation in certain limiting cases. The numerical evaluation
of analytical expression is done in the limits where such approximations
become unavailable. We begin by using the kinetic dispersion relation for the
electromagnetic waves (Bittencourt (2013) and Sec. V.2) given by,

$\displaystyle
c^{2}k^{2}-\omega^{2}=\frac{8\pi^{2}q_{e}^{2}}{m}\int\limits_{-\infty}^{\infty}\int\limits_{0}^{\infty}$
$\displaystyle\frac{\left[(\omega-kv_{z})\frac{\partial f_{e0}}{\partial
v_{\perp}}+kv_{\perp}\frac{\partial f_{e0}}{\partial v_{z}}\right]}{(\omega-
kv_{z}\pm\Omega_{ce})}$ (8) $\displaystyle\times$ $\displaystyle
v_{\perp}^{2}dv_{\perp}dv_{z}.$

Assuming isotropy in the velocity space,

$\displaystyle\frac{\partial f_{0}}{\partial v_{\|}}=\frac{\partial
f_{0}}{\partial v_{\perp}},$ (9)

additionally considering the whistler wave frequency limit
$\Omega_{ci}\ll\omega\ll\Omega_{ce}$ and for wavelengths larger than electron
gyroradius, $\lambda>\rho_{ce}$, we recover Krall and Trivelpiece (1986),

$\displaystyle\gamma\approx-\frac{\omega_{pe}^{2}}{|k_{\|}|v_{th}}\frac{1}{1+k^{2}c^{2}/\omega_{r}^{2}}\exp\left(-\frac{\Omega_{ce}^{2}}{k_{\|}^{2}v_{th}^{2}}\right).$
(10)

However, for the cases explored in the present simulations have
$\omega\lesssim\Omega_{ce}$ and the expression (10) underestimates the damping
rate $\gamma$, yielding negligible values in comparison to what recovered in
the simulations. For example, for $k=1.0\omega_{pe}/c$, $v_{th}=0.026c$,
$\Omega_{ce}=0.1\omega_{pe}$, $\gamma\sim 10^{-8}\omega_{pe}$, although a
considerably higher damping ($\gamma\sim 10^{-3}\omega_{pe}$) is recovered in
the simulations. An estimate from (8) in the limit $\omega\lesssim\Omega_{ce}$
more relevant to simulations is prescribed by Gary Gary (1993) in the form of
a general expression for $\gamma$,

$\displaystyle\gamma\approx-\frac{\pi}{2k}\sum_{\alpha}\frac{\omega_{p\alpha}^{2}}{\sqrt{2\pi}v_{th\alpha}}\exp\left(-\frac{(\omega\pm\Omega_{c\alpha})^{2}}{2k_{\|}^{2}v_{th\alpha}^{2}}\right)$
(11)

which can be used as an alternate analytic estimate to compare with theory the
whistler damping rate recovered in the simulations. For an even improved
analytic estimate of damping rates we have also used the exact expression (8)
and evaluated $\gamma$ numerically, by the following procedure.

As in the simulation where ions are infinitely massive, considering electron
equilibrium distribution to be Maxwellian,

$\displaystyle
f_{0}=\left(\frac{1}{2\pi}\right)^{\frac{3}{2}}\frac{1}{v_{th\parallel
e}v_{th\perp e}^{2}}\exp\left(-{\frac{v_{\perp}^{2}}{2v_{th_{\perp
e}}^{2}}-\frac{v_{\parallel}^{2}}{2v_{th_{\parallel e}}^{2}}}\right),$

and substituting it in the dispersion (8), one obtains Chen  _et al._ (2013),

$\omega^{2}+\omega_{pe}^{2}I-c^{2}k^{2}=0,$ (13)

where,

$\displaystyle I=\frac{v_{th\perp e}^{2}}{v_{th\parallel
e}^{2}}-1+\left[\frac{v_{th\perp e}^{2}}{v_{th\parallel
e}^{2}}\frac{\omega\pm\Omega_{ce}}{\mp\Omega_{ce}}+1\right]\frac{\mp\Omega_{ce}}{\sqrt{2}kv_{th\parallel
e}}Z(\zeta),$ (14)

and $Z(\zeta)$ is the plasma dispersion function Abramowitz and Stegun (1964),

$Z(\zeta)=\frac{1}{\sqrt{\pi}}\int\limits_{-\infty}^{\infty}\frac{e^{-t^{2}}}{t-\zeta}dt,$
(15)

with the argument

$\zeta=\frac{\omega\mp\Omega_{ce}}{\sqrt{2}kv_{th\parallel e}}.$ (16)

The upper and lower sign in Eq. (14) are for RHCP and LHCP waves,
respectively. Applying dispersion (13) to whistlers which are right handed
polarized, and considering isotropy of the distribution
($v_{th\parallel}=v_{th\perp}=v_{th}$) we get (13) in the form,

$\omega^{2}+\omega\frac{\omega_{pe}^{2}}{\sqrt{2}kv_{th}}Z(\zeta)-c^{2}k^{2}=0.$
(17)

Note that under cold plasma approximation where imaginary part of the function
$Z(\zeta)=i\sqrt{\pi}e^{-\zeta^{2}}-\frac{1}{\zeta}$ is negligible and
$Z(\zeta)\sim-1/\zeta$, such that with
$\zeta=(\omega-\Omega_{ce})/\sqrt{2}kv_{th}$ one readily recovers from (17),
the cold plasma dispersion relation for the RHCP waves,

$\omega^{2}-\omega\frac{\omega_{pe}^{2}}{(\omega-\Omega_{ce})}-c^{2}k^{2}=0,$
(18)

which is the whistler dispersion relation in the limit
$\Omega_{ci}<\omega<\Omega_{ce}$. Considering $\omega=\omega_{r}+i\omega_{i}$,
in the warm plasma RHCP dispersion (17) we obtain for
$\gamma=\omega_{i}\ll\omega_{r}$ Chen  _et al._ (2013); Schreiner (2016),

$\gamma=-\frac{\omega_{r}\sqrt{\pi}e^{-\zeta^{2}}\sqrt{2}kv_{th}}{[\omega_{r}Z^{\prime}(\zeta)+\sqrt{2}kv_{th}Z(\zeta)+2\omega_{r}k^{2}v_{th}^{2}]},$
(19)

where $Z^{\prime}({\zeta})$ is the derivative of $Z({\zeta})$ with respect to
$\zeta$.

Comparison of the damping rate recovered in simulations with the analytical
approximations (10), (11) and (19) is presented in Fig. 4 for the above range
of $v_{th}$ values for which simulations are performed using
$k=0.897\omega_{p}/c$ or simulation box length $L=7c/\omega_{p}$.

Figure 4: Comparison of simulated damping rate $\gamma$ with that obtained
from analytical approximations, approx-1 (10), approx-2 (19) and approx-3 (11)
done for $k=1.0\omega_{p}/c$ or simulation box length $L=2\pi c/\omega_{p}$.

Note that for the nearly cold electron case, $v_{th}\rightarrow 0$, the
simulations (blue dots) duly recover a whistler propagation free of any
damping (or growth), confirming that the deviation from this undamped
propagation recovered at large electron temperature values (larger $v_{th}>0$)
represents pure kinetic effects in the simulations. The simulations however
show increasing damping from finite $v_{th}$ values while all the other
analytical approximations prescribe comparable damping only beyond
$v_{th}>1.0\times 10^{-2}c$, as discussed further below.

Figure 5: Simulated frequency and damping rate $\gamma$ as function of $k$ and
a range in electron thermal velocity values $v_{th}=0.001$ to $0.026$. The
vertical dashed lines indicate $k$ values for which sets of simulations are
performed by varying $v_{th}$.

Suitability of approximation (10) for the present range is ruled out by the
negligible resonant damping prescribed by (10) as plotted in Fig. 4 using
green line in comparison to other estimates of $\gamma$, namely, (11) and (19)
plotted with red and black lines, respectively.

Figure 6: Comparison of simulated frequencies (a) $\omega_{r}$ and (b)
$\omega_{r}$ with those obtained from analytical approximation (19) and
corresponding $\omega_{r}$ as functions of $k$ and $v_{th}$. Markers at
various $k$ values represent simulation data such that the markers and curves
from top to bottom correspond to $v_{th}=0.001$ to $0.026$.

The analytical approximation (11) very closely estimates the numerically
evaluated values of the exact expression (19). We note that despite no clear
approximations used in them (at least, in the numerically computed $\gamma$
from (19)) besides their linear origin, they still show finite discrepancy
from the simulated damping rates. In general, the simulations recover stronger
damping rates for the whistlers than what prescribed by the analytical damping
rate approximations duly obtained from the kinetic dispersion relation. We
explore the origin of this discrepancy for whistlers excitations in Sec. V
based on the deterministic evolution of the electron distribution function in
essential 4D phase-space set up available from the simulations. The missing
physics is located in the small time evolution of the initial perturbations as
asymptotic contributions are duly accounting for in estimating the kinetic
Landau damping rate. For turbulent situations, where the excited fluctuations
are short lived, the short time contributions need to be accounted for.
Moreover, they might dominate the asymptotic mechanisms in certain relevant
limits.

In Fig. 5 the analytically obtained whistler frequency $\omega$ and more
accurate numerically computed damping rate from (19) are plotted as function
of $k$ for several ($>10$) different values of electron thermal velocity in
range $v_{th}=0.001$ to $0.026$. In Fig 6, the $\omega_{r}$ and $\gamma$
values are separately compared with analytical values where the markers
represent simulated data and variation is done in both $k$ and $v_{th}$
values. As in Fig 6, simulations confirm the analytical prediction that the
$\omega_{r}$ slightly reduces at larger electron temperature, though the
analytical $\omega_{r}$ overestimates this reduction for larger $v_{th}$. The
comparison of $\gamma$ presented in Fig. 6 shows that while at
$k=0.8\omega_{p}/c$ the damping rate in simulation is recovered to be stronger
than the corresponding analytic values, at larger $k$ values,
$0.897~{}\omega_{p}/c$ and $1~{}\omega_{p}/c$, the overlap in the analytic and
simulated $\gamma$ values is better, especially at larger $v_{th}$ values. It
can also be noted that at smaller $k$, analytical expressions underestimate of
the damping and the spread in kinetic $\gamma$ values is larger than the
analytical value spread. At larger $k$, on the other hand, the analytical
approximation overestimates $\gamma$ and the kinetically obtained $\gamma$
values have a narrower spread.

##  V Phase-mixing supplemented damping of the EMHD perturbations

Relatively larger damping of the whistler perturbations at smaller $k$ in
comparison to the analytical estimates (10), (11) and (19) is interpreted
based on the deterministic evolution of the electron distribution function
governed by the Vlasov equation as simulated and analyzed in this section. In
warm electron plasmas the initial damping is found to have considerable
contribution from the corresponding deterministic phase-mixing process of the
fresh (initial) perturbations made to the plasma, confirming to a propagating
whistler wave. Note that the term phase-mixing used here refers to the varied
response of the electrons of different velocities (from a wider spread about
mean) to the propagating wave perturbation. The deterministic deformations of
the velocity distribution function that evolve antisymmetrically about a
reference velocity (distinct for electrostatic and electromagnetic waves) in
the phase-space however introduce incoherence between the substructures of
observable macroscopic spatio-temporal variations, or collective wave
perturbations. The resultant weakening of the perturbations translates in to
an additional damping (than what caused by the resonant electrons) detectable
in the deterministic Vlasov simulation data. At shorter time scales this
contribution supplements, and modifies, the damping of the (transverse) wave,
producing enhanced damping rates in comparison to the pure resonant damping
essentially estimated considering the long time limit Krall and Trivelpiece
(1986). Clear signatures of both resonant and non resonant processes are
visible in the simulated 3-dimensional electron velocity distribution function
evolving in the phase-space.

###  V.1 Phase-space evolution of whistler wave perturbation in the Vlasov
simulations

The plots of electron distribution function perturbation $f_{1}=f-f_{0}$, as
available from the simulations done with denser grid size, $64\times 160\times
160\times 160$, are analyzed in order to understand the kinetics of the
evolution of whistler wave perturbation in the plasma. Since the system is
periodic and transverse, analyzing the $f$ data from any chosen spatial
location is equivalent. The location $z=L/2$, or the center of the simulation
domain, is therefore chosen here as the representative location in the
simulation domain. For added accuracy, both $f$ and $f_{0}$ are evolved under
equivalent Vlasov simulation operation in order to obtain $f_{1}(t)$ presented
below. From the full 3-dimensional velocity distribution function available,
the evolution of $f_{1}(z=L/2,v_{z},v_{x},v_{y})$ at various stages of time is
presented in Fig. 7 from the cold electron case simulation, performed with
sufficiently small $T_{e}$ such that the $|v_{\rm res}|\gg v_{\rm
th}=0.001~{}c$. Because of the small electron temperature, both the resonant
and phase velocity of the whistler are far removed from the distribution width
(along $v_{z}\|{\bf k}$) and the $f_{1}$ plotted in three columns of frames
from left to right correspond to velocity values $v_{z}=-v_{\rm th}$, $0$ and
$v_{\rm th}$, respectively. The same data is plotted in Fig. 8 from simulation
of warmer electrons case, $|v_{\rm res}|\leq v_{\rm th}=0.001~{}c$, such that
sufficient population of resonant electrons is present with $v_{z}=v_{\rm
res}$, producing dominant and analyzable kinetic effects. The rows of frames
from top to bottom in both plots correspond to time from $t=0$ to 20
$\omega_{pe}^{-1}~{}(\equiv 2\Omega_{ce}^{-1})$.

Figure 7: Evolution of the electrons velocity distribution function
perturbation in the transverse velocity-space plane $v_{x}$-$v_{y}$ for the
cold electrons case, with $v_{th}=0.001c$. Frames from top to bottom
correspond to equal intervals over one whistler cycle ($\sim
20\omega_{pe}^{-1}$) from $t=0$, while those from left to right are for
$v_{z}=-v_{\rm th}$, 0, and $v_{\rm th}$, respectively. Red and black arrows
indicate, ${\bf B}_{1}$ and $\langle{\bf v}\rangle_{1}$, respectively.  Figure
8: Evolution of the electrons velocity distribution function perturbation in
the transverse velocity-space plane $v_{x}$-$v_{y}$ for the warm electrons
case, with $v_{th}=0.02c$. Frames from top to bottom correspond to equal
intervals over one whistler cycle ($\sim 20\omega_{pe}^{-1}$) from $t=0$,
while those from left to right are for $v_{z}=v_{\rm res}$, 0, and $v_{\rm
phase}$, respectively. Red and black arrows indicate, ${\bf B}_{1}$ and
$\langle{\bf v}\rangle_{1}$, respectively.

In both cold and warm electron cases plotted in top rows of Fig. 7 and Fig. 8,
respectively, the $f_{1}$ plotted at $t=0$, consistent with a whistler phase-
space perturbation in the velocity space, has a nonzero velocity perturbation
purely aligned to $\hat{v}_{y}$. The following rows cover evolution over a
period of $20\omega_{pe}^{-1}$, or one complete whistler cycle. We first
discuss evolution of perturbation in the cold electron case presented in Fig.
7. Consider the evolution of the perturbation component located at $v_{z}\sim
0$ plotted in the second column of the Fig. 7. Because of being drawn from
vanishing $v_{z}$, this component of $f_{1}$ does not have electrons that
leave to, or arrive from, the adjacent $z$ locations as the time progresses
(i.e., from top frame to bottom frame). Therefore, the participating electrons
are simply displaced to newer angular locations with (angular) velocity
$\Omega_{ce}$ on the same transverse velocity plane $v_{x}$-$v_{y}$ in the
subsequent frames, with respect to their original locations in the $t=0$
frame. However the polarization of the perturbation in this case rotates with
angular velocity of the wave $\omega$, dissimilar to $\Omega_{ce}$ ($>\omega$,
the reason for which is well encoded in the linear theory of waves in
magnetized plasmas). Drawn from the simulation data for ${\bf\langle
v\rangle}_{1}$ and ${\bf B}$, the black and red arrows represent phases of
right handed polarized velocity and magnetic field vectors, respectively, for
the reference. The rotation in perturbation polarization consistent with wave
frequency confirms that the initial perturbation efficiently excited a
whistler eigenmode in the simulated magnetized plasma right from $t=0$.

The kinetic effects because of resonant electrons are addressed further below
in warm electron case since for cold electrons the velocity distribution drops
to negligible values at the resonant velocity. The same is evident from the
identical evolution of the perturbation component at $v_{z}=|\pm
v_{th}|\ll|v_{\rm res}|$, plotted in the columns 1 and 3, respectively, which
show rotation of perturbation polarization with same frequency and a coherent
response of electrons of all $v_{z}$ values due to the distribution function
being nearly a delta function in velocities both along and across ${\bf k}$.
Therefore, the fluid theory results ($\gamma\rightarrow 0$, $\omega$ given by
(18)) remain suitably recoverable for this cold electron case as clear from
the data point corresponding to $v_{th}\sim 0.001c$ in Fig. 4.

We now show that an identical initial perturbation, that excited a nearly
ideal whistler eigenmode in the above cold electron case, evolves in a warm
electron magnetized plasma exhibiting a far richer kinetic response. This
initial evolution, covering $20\omega_{pe}^{-1}$ remains applicable, for
example, to the short lived transverse whistler mode fluctuations
characteristically running much of weak electromagnetic plasma turbulence Gary
_et al._ (2010). There evolution therefore necessitates invoking initial
phase-mixing description of analytical kinetic theory results in addition to
the routine estimates of damping by pure resonant electrons.

Considering the evolution for the perturbation component at $v_{z}\sim
0,z=L/2$ plotted in the middle column of the frames in Fig. 8 we note that its
polarization rotates nearly with the wave frequency $\omega$, however with a
visible shortfall in comparison to the cold electron case. The electrons being
sufficiently warm in this case, it is now possible to examine the $f_{1}$
evolution at the resonant velocity for this set of data, which is noted to be
much distinct as plotted in the first column of the Fig. 8.

Note that individual electrons streaming with resonant velocity $v_{z}$ would
see the wave frequency Doppler-shifted to $\Omega_{ce}$ and therefore will be
perpetually accelerated by the wave ${\bf E}_{1}$ in transverse the direction.
The resulting proliferation of electrons to higher $v_{\perp}$ is visible with
time in frames from (a)-(f). It is important to notice, though, that, unlike
the coherent response of the full distribution in cold electron case, the
polarization of the perturbation in the resonant electrons has rotated with a
considerably lower frequency, thereby covering a lesser angle as compared to
the bulk, or $v_{z}\sim 0$, electrons.

Quantitative estimation of this shift in frequency is possible based on the
fact that because of $\Omega_{ce}>\omega$, in time $\delta t$ the gyrating
electrons cover a far larger angle in the transverse plane than the rotating
polarization of the perturbation. Consequently, as the time progresses the
$v_{z}=v_{\rm res}$ electrons bring imprint of an increasingly older phase of
the wave, distinct from the instant phase of the bulk electron component
polarization (seen in the corresponding middle column frame plotted for
$v_{z}=0$). This is despite of resonant electrons arriving (because of their
negative streaming velocity $v_{\rm res}<0$) from a forward location $z=\delta
z>L/2$ having a slightly advanced phase of the wave.

The convected phase $\phi(z=L/2,t=\delta t)$ of the polarization for resonant
perturbation component (perturbation component with $v_{z}=v_{\rm res}$) at
time $t=\delta t$ can be calculated by using a transformation between present
and older phase of the polarization of the perturbation,

$\displaystyle\phi(L/2,\delta t)=\phi(L/2+\delta
z,0)=\phi_{0}+\delta\phi_{z}+\delta\phi_{\theta},$ (20)

where $\phi_{0}$ is instant phase of the polarization of the bulk perturbation
at $z=L/2$ (as visible in $v_{z}=0$, $t=\delta t$ frame), $\phi(z,t)$ is the
general phase of the polarization (ideally known from the cold electron
limit), $\delta\phi_{z}$ and $\delta\phi_{\theta}$ are the phase difference
arising from change in $z$ and $\theta$, respectively. For the resonant
component of the perturbation, we have,

$\displaystyle\delta\phi_{z}=\delta zk=v_{\rm res}(-\delta
t)k;~{}\delta\phi_{\theta}=(-\delta t)~{}\Omega_{ce}.$ (21)

Since $v_{\rm res}<0$ and $\Omega_{ce}>\omega=\Omega_{ce}-|v_{\rm res}|k$, we
have $\delta\phi_{z}>0$, $\delta\phi_{\theta}<0$ and
$|\delta\phi_{\theta}|>\delta\phi_{z}$, meaning that the phase of resonant
component of the perturbation relative to instant bulk electron phase will be
negative. These phase delays are in excellent agreement with the phase
difference in the frames presenting the resonant component of the perturbation
in the first column of Fig. 8.

An even distinct evolution is witnessed in the perturbation component at the
positive values of parallel velocity, $v_{z}>0$, as potted in third column of
Fig. 8 for $v_{z}=v_{\rm phase}$. Before explanation of this behavior, which
strongly relates to phase-mixing process of the whistler perturbation, it can
be clearly noted that the strength of macroscopic wave variables (especially
${\bf J}_{1}$) is going to be weakened with time if velocity integration is
performed at advanced times using $f_{1}(v_{x},v_{y})$ that has its peaks in
$v_{x}$-$v_{y}$ plane initially overlapping for all $v_{z}$, but dispersed, at
future times, over a wider range of angles for different $v_{z}$.

In the $f_{1}$ evolution plotted at $v_{z}=v_{\rm phase}$ in third column of
Fig. 8, the streaming (and simultaneously gyrating) electrons experience the
wave fields purely non-rotating. At each $\delta t$ they import, from their
original locations to $z=L/2$, the phase as determined by the transformation
(20). In this case, however, with each gyration completed, the electrons
additionally drift orthogonal to the phase of the wave electric field
experienced by them. Consequently, the imported $f_{1}$ is additionally
displaced away from $v_{\perp}=0$. With time progressing, the $f_{1}$ steadily
sampled at $z=L/2$ witnesses newer electrons arriving from backward locations.
An increasing drift is therefore witnessed in $f_{1}$, rotated at an
incremented angle in $v_{x}$-$v_{y}$ plane. This results in periodically
separating red and blue patches, absent initially. Note that the initial
transverse velocity perturbation at $t=0$ is introduced with no such
nonuniformity with respect to $v_{z}$. This nonuniformity clearly develops
over the time-scale of an electron gyroperiod.

###  V.2 Kinetic model for initial-value transverse electromagnetic
perturbations

Out of the several kinetic modifications of $f_{1}$ witnessed above, those
attributable to the wave fields and those to ballistic effects can be
distinguished and analyzed based on the linearized Laplace-Fourier transformed
Vlasov-Maxwell system of equations,

$\displaystyle(p$ $\displaystyle+$ $\displaystyle i{\bf k\cdot v})f_{kp1}({\bf
v})=f_{k1}(t=0)$ (22) $\displaystyle+$
$\displaystyle\frac{e}{m_{e}}\left(\frac{\bf v\times
B_{0}}{c}\right)\cdot\nabla_{v}f_{kp1}({\bf v})$ $\displaystyle+$
$\displaystyle\frac{e}{m_{e}}\left({\bf E}_{kp1}+\frac{{\bf v\times
B}_{kp1}}{c}\right)\cdot\nabla_{v}f_{0}({\bf v})$ $\displaystyle{\bf
B}_{kp1}={\bf B}_{k1}(t=0)+\frac{c}{ip}({\bf k\times E}_{kp1}).$ (23)

and

$\displaystyle{\bf E}_{kp1}={\bf E}_{k1}(t=0)-\frac{c}{ip}\left({\bf k\times
B}_{kp1}-\frac{4\pi}{c}{\bf j}_{kp1}\right),$ (24)

where $p=-i\omega_{r}+\omega_{i}$ and ${\bf E}_{kp1}$ and ${\bf B}_{kp1}$ are
the transformed wave electric and magnetic fields, respectively.

Using the well known identity for conjugate components of the velocity
variable, $v_{x}=v_{\perp}\cos{\phi}$ and $v_{y}=v_{\perp}\sin{\phi}$,

$\displaystyle\frac{\partial f_{pk1}({\bf
v})}{\partial\phi}=\left(-v_{y}\frac{\partial}{\partial
v_{x}}+v_{x}\frac{\partial}{\partial v_{y}}\right)=-({\bf
v\times\hat{z}})\cdot\nabla_{v}f_{pk1}({\bf v}),$

Eq. (22) can be rewritten as,

$\displaystyle\frac{(p+i{\bf k\cdot v})}{\Omega_{ce}}f_{kp1}({\bf
v})+\frac{\partial f_{pk1}({\bf
v})}{\partial\phi}=\frac{f_{k1}(t=0)}{\Omega_{ce}}$
$\displaystyle+\frac{e}{m_{e}\Omega_{ce}}\left({\bf E}_{kp1}+\frac{{\bf
v\times B}_{kp1}}{c}\right)\cdot\nabla_{v}f_{0}({\bf v}).$ (25)

Figure 9: Time evolution of the maximum amplitude, $f_{\rm 1max}$, of
electrons velocity distribution function perturbation $f_{1}(v_{x}$,$v_{y})$
plotted in Figs. 7 and 8 in full range of $v_{z}$, for the hot electrons case
($v_{th}=0.02c$). Plots (a) and (b) present relatively rotated views of the
same surface for clearer exposure of the both sides about the surface maximum.
Figure 10: Surface in the figure normalized to individual initial values at
each $v_{z}$.

Note that the wave field generated $v_{z}$ specific modifications, including
(i) the acceleration of resonant electrons and (ii) the slowly building
electron drift for $v_{z}\sim v_{\rm phase}$, can be resolved, approximately,
from the fundamentally kinetic ballistic modifications in the $f_{1}$ at
shorter times ($\Omega_{ce}^{-1}<t<\omega^{-1}$) when the latter dominate the
former. In order to do this, we first write the right handed polarized
component of the Vlasov equation Bittencourt (2013) applicable to the whistler
mode,

$\displaystyle\frac{d}{d\phi}\left\\{f_{1+}({\bf
v})\exp{\left[\frac{(p+ikv_{z})}{\Omega_{ce}}\phi\right]}\right\\}=$
$\displaystyle\frac{f_{k1+}(t=0)}{\Omega_{ce}}\exp{\left[\frac{(p+ikv_{z})}{\Omega_{ce}}\phi+i\phi\right]}$
$\displaystyle+F_{+}(v_{z},v_{\perp})\exp{\left[\frac{(p+ikv_{z})}{\Omega_{ce}}\phi+i\phi\right]}$
(26)

where $v_{\perp}=(v_{x}^{2}+v_{y}^{2})^{\frac{1}{2}}$, $f_{k1+}$ is right
handed circularly polarized part of the initial perturbation. The quantity
$F_{+}(v_{z},v_{\perp})$ is obtained by ignoring ${\bf B}_{k}(t=0)$ in Eq.
(23), as,

$\displaystyle
F_{+}(v_{z},v_{\perp})=\frac{e}{m_{e}\Omega_{ce}}\left[\left(1-\frac{kv_{z}}{\omega}\right)\frac{\partial
f_{0}}{\partial v_{\perp}}+\frac{kv_{z}}{\omega}\frac{\partial f_{0}}{\partial
v_{\perp}}\right]\frac{E_{+}}{\sqrt{2}}$

where $E_{+}=E_{x}-iE_{y}$. Both sides of Eq. (26) can be integrated with
respect to $\phi$ from $\phi=-\infty$ to instant $\phi$ by assuming very small
imaginary part in $\omega$, yielding the time asymptotic solution,

$\displaystyle f_{pk1+}({\bf v})$ $\displaystyle=$
$\displaystyle\frac{f_{k1+}(t=0)e^{i\phi}}{(p+ikv_{z}+i\Omega_{ce})}$ (28)
$\displaystyle+$
$\displaystyle\frac{\Omega_{ce}}{(p+ikv_{z}+i\Omega_{ce})}F_{+}(v_{z},v_{\perp})e^{i\phi}.$

This means that all the right handed circularly polarized transverse modes are
obtainable by accounting for the poles of the second term, which represents
the contribution of the electric and magnetic fields of the wave. For example,
the Whistler mode corresponds to $E_{+}$ obtained in terms of $B_{+}$ and
$J_{+}$ using (22)-(24) in the limit $\Omega_{pi}<\omega<\Omega_{ce}$ such
that the relevant pole appears where the dispersion function for an isotropic
electron distribution,

$\displaystyle
c^{2}k^{2}-\omega^{2}-\frac{8\pi^{2}q_{e}^{2}}{m}\int\limits_{-\infty}^{\infty}\int\limits_{0}^{\infty}\frac{\left[(\omega-
kv_{z}+kv_{\perp})\frac{\partial f_{e0}}{\partial v_{z}}\right]}{(\omega-
kv_{z}\pm\Omega_{ce})}$ $\displaystyle\times v_{\perp}^{2}dv_{\perp}dv_{z}$
(29)

vanishes, yielding both $\omega_{r}$ and $\gamma$ Chen  _et al._ (2013).
However, the solution (28) indicates that $f_{pk1+}({\bf v})$, being the
integrand of the Landau-like Laplace inversion integral producing
$f_{k1+}({\bf v},t)$ (which can, in turn, be compared to the simulation
outputs plotted in Fig. 7 and 8), has an additional pole at,

$\displaystyle p=-i(kv_{z}+\Omega_{ce}).$ (30)

Since $k$, $v_{z}$ and $\Omega_{ce}$ are real, there is no reduction expected
with time in the amplitude of perturbation components corresponding to various
$v_{z}$. Instead, $f_{k1+}(v_{z},t)$ must develop temporal oscillations with a
pure real frequency $\omega_{b}(v_{z})$ directly proportional to $v_{z}$. As
seen below, this temporal modulation is caused by the information of the wave
phases at remote locations being brought by the electrons at the location of
measurement, at different temporal rates for different $v_{z}$.

The net (kinetic) reduction in the wave strength is therefore caused by both,
(i) a phase-mixing process present throughout the velocity space and (ii) the
resonant damping which is dominant at the resonant velocity, accelerating only
the resonant electrons. While both these processes in a collisionless warm
plasma conserve the entropy despite diminishing the wave strength Krall and
Trivelpiece (1986), only the latter contributes to the routinely estimated
$\gamma$. The temporal damping caused by former involves pure frequencies and
cause a rather deterministic thermalization (reversible or entropy preserving
ergodization), or a phase-mixing, of the transverse distribution (and,
consequently of ${\bf J}_{+}$). As presented below, this phase-mixing is
clearly recoverable from the output of a deterministic simulation performed in
this study, duly described by the above analytical treatment.

Figure 11: The simulated $v_{z}$ dependence frequency $\omega_{b}$ (dotted
line, circular markers) estimated from the Fig. 10 compared with analytical
expression (30) (solid line) for frequency of amplitude oscillation resulting
from phase-mixing.

The maxima of $f_{1}(v_{x},v_{y})$ presented in Fig. 8, corresponding to hot
electron case, are plotted separately as a function of $v_{z}$ and $t$ in the
form of the surface $f_{\rm 1max}(v_{z},t)$ in Fig. 9. For the convenience,
two views of the surface are presented in (a) and (b) with relative rotation
to clearly highlight the visible temporal variation of $f_{\rm max}$ at
various $v_{z}$, including $v_{z}=v_{\rm res}$ as indicated by an arrow. The
perturbation amplitude in Fig. 9 shows that only at $v_{z}=v_{\rm res}$ a
nearly monotonic variation is present in the amplitude, whereas at velocities
away from $v_{\rm res}$ a slow temporal oscillation enters in the amplitude of
$f_{1}$ with its frequency increasing with $|v_{z}|$. The high frequency
oscillations in small $v_{z}$ region are remnants of high frequency
($\omega_{R,L}>\Omega_{ce}$) RHCP and LHCP waves (see Fig. 3) and should not
be confused with the discussed low frequency high amplitude oscillations in
the small $v_{z}$ region.

This additional frequency in the surface plot contours for larger $v_{z}$
$(>0)$ results from the streaming electrons importing information of
polarization from adjacent locations based on the transformation Eq. (20).
Note that because of the periodic set up, for a larger velocity, say
$v_{z}\sim n(k\delta t)^{-1}$, the information being imported would contain
$n$ cycles of modulation, corresponding to $n$ consecutive spatial periods
sampled by fast streaming electrons in the time interval $\delta t$ past
$t=0$. The frequency of temporal modulation in $f_{\rm 1max}(v_{z},t)$
therefore must increase with $v_{z}$ value being examined at a fixed location
$z=L/2$. This dependence on $v_{z}$ is indeed confirmed by plotting, in Fig.
10, the amplitude $f_{\rm 1max}(v_{z},t)$ normalized to its $t=0$ reference
value for each $v_{z}$, such that the modulations become clearly identifiable
even at the larger $v_{z}$ where the distribution function $f_{0}$ as well as
the amplitude of $f_{1}$ drop to significantly small values.

The perturbation information at $z=L/2$ arriving loaded on electrons streaming
with phase velocity is played faster than the rotation of bulk polarization
(as $\omega<\Omega_{ce}$) but imports an increasingly delayed phase (causing
the rotation appearing to stop or reversed over certain intervals).
Additionally, after certain time interval, it begins to display the signatures
of freshly developed ${\bf E}_{1}\times{\bf B}_{0}$ drift.

Finally, since electrons with $v_{z}>v_{\rm res}$ are effectively magnetized
to a higher degree on the wave time scale experienced by them, $(\omega-
kv_{z})^{-1}\gg\Omega_{ce}^{-1}$, they participate rather efficiently in the
collective wave propagation. However, a Doppler down-shifted frequency
experienced, and responded to, by this majority of electrons causes the wave
collective $\omega_{r}$ to reduce in accordance with dispersion (17). This
reduction in $\omega_{r}$ is duly captured by the present Vlasov simulation
results, as clear by copmparing Fig. 7 with Fig. 8.

##  VI Summary and Conclusions

The kinetic damping mechanism of low frequency transverse perturbation
propagating parallel to the magnetic field in a magnetized warm plasma is
simulated by means of Vlasov simulations. The impact of finite temperature on
the wave frequency and the underlying physics of the analytical approximations
prescribing reduction in the real frequency of the wave in a warm plasma is
addressed and resolved by full 4D phase-spatiotemporal evolution simulated of
the electron distribution function.

The analysis of the Vlasov simulation output illustrates and estimates the
evolution of initially coherent perturbations that undergo damping by mutually
separable initial kinetic evolution operating over the global and local
regions of the of electron velocity space. The net kinetic effects contributed
reduction in the strength of perturbations is caused by both, (i) a phase-
mixing process present throughout the velocity space and (ii) the resonant
damping which is dominant at the resonant velocity, accelerating only the
resonant electrons. The quantitative analysis of the amplitude oscillation
produced by the former is done and their interplay, at short times, with the
modifications of the initial perturbations by the resonant electrons is shown
to result in enhanced decay of the perturbation amplitude as compared to
typical long time estimates of the damping. It is emphasized that while both
these mechanisms in a collisionless warm plasma conserve the entropy, it is
only the resonant effects that contributes to the routinely estimated damping
rate $\gamma$. The analytically predicted reduction in the real wave frequency
$\omega_{r}$ is duly recovered by the simulations. This is shown to follow
from a greater effective magnetization of dominant $v_{z}>v_{\rm res}$
electron population, which, however, experiences frequency of the wave
Doppler-shifted to smaller values.

We conclude by discussing that access to rather self-consistent analytical
initial value solutions of the problem (22)-(24) remain of high relevance.
They should enable the output from the present fully kinetic electromagnetic
transverse wave Vlasov simulations to more carefully determine, for example,
the spectral boundaries of the electromagnetic collisionless plasma
turbulence. The associated particle and energy transport in the space and
fusion plasma, in turn, can also be approached more deterministically.

Acknowlegement: The simulations are performed on the IPR HPC supercomputing
cluster ANTYA.

## References

  * Borovsky (2012) J. E. Borovsky, Journal of Geophysical Research 117, 5104 (2012).
  * Podesta  _et al._ (2007) J. J. Podesta, D. A. Roberts, and M. L. Goldstein, The Astrophysical Journal 664, 543 (2007).
  * Gary  _et al._ (2010) S. P. Gary, S. Saito, and Y. Narita, The Astrophysical Journal 716, 1332 (2010).
  * Xu  _et al._ (2017) X. Xu, L. Chen, C. Zhou, X. Liu, Z. Xia, J. J. Simpson, and Y. Zhang, Journal of Geophysical Research: Space Physics 125, e2019JA027750 (2017).
  * Helliwell (1965) R. A. Helliwell,  _Whistlers and related ionospheric phenomena_ (Press, Stanford, Calif: Stanford Univ., 1965).
  * Xiao  _et al._ (2017) F. Xiao, S. Liu, X. Tao, Z. Su, Q. Zhou, C. Yang, Z. He, Y. He, Z. Gao, D. N. Baker, H. E. Spence, G. D. Reeves, H. O. Funsten, and J. B. Blake, Journal of Geophysical Research: Space Physics 122, 3201–3211 (2017).
  * Thorne (2010) R. M. Thorne, Geophys. Res. Lett. 37, L22107 (2010).
  * Chen  _et al._ (2013) L. Chen, R. M. Thorne, Y. Shprits, and B. Ni, Journal of Geophysical Research: Space Physics 118, 2185 (2013).
  * Krall and Trivelpiece (1986) N. A. Krall and A. W. Trivelpiece,  _Principles of Plasma Physics_ (San Francisco Press Inc., San Francisco, 1986).
  * Landau (1944) L. D. Landau, C. R. Acad. Sci. U. R. S. S. 44, 311 (1944).
  * Landau (1946) L. D. Landau, J. Phys. U.S.S.R. 10, 25 (1946).
  * Lutomirski (1970) R. Lutomirski, The Physics of Fluids 13, 149 (1970).
  * Fijalkow (1999) E. Fijalkow, Computer physics communications 116, 329 (1999).
  * Mangeney  _et al._ (2002) A. Mangeney, F. Califano, C. Cavazzonic, and P.Travnicek, Journal of Computational Physics 179, 495 (2002).
  * Palodhi  _et al._ (2010) L. Palodhi, F. Califano, and F. Pegoraro, Plasma Physics and Controlled Fusion 52, 095007 (2010).
  * Valentini  _et al._ (2007) F. Valentini, P. Travnicek, F. Califano, P. Hellinger, and A. Mangeney, Journal of Computational Physics 225, 753 (2007).
  * Shuster  _et al._ (2021) J. R. Shuster, N. Bessho, S. Wang, and J. Ng, Physics of Plasmas 28, 122902 (2021).
  * Fülöp  _et al._ (2009) T. Fülöp, H. Smith, and G. Pokol, Physics of Plasmas 16, 022502 (2009).
  * Schreiner  _et al._ (2017) C. Schreiner, P. Kilian, and F. Spanier, Communications in Computational Physics 21, 947 (2017).
  * Schreiner (2016) C. Schreiner,  _Numerical modelling of the microphysical foundation of astrophysical particle acceleration_ , Ph.D. thesis, North-West University (South Africa), Potchefstroom Campus (2016).
  * Mandal and Sharma (2016) D. Mandal and D. Sharma, Journal of Physics: Conference Series 759, 012068 (2016).
  * Mandal  _et al._ (2020) D. Mandal, D. Sharma, and H. Schamel, Phys. Plasmas 27, 022102 (2020).
  * Gary (1993) S. P. Gary,  _Theory of space plasma microinstabilities_ , 7 (Cambridge university press, 1993).
  * Chen (2012) F. F. Chen,  _Introduction to Plasma Physics_ (Springer Science & Business Media, 2012).
  * Stix (1992) T. H. Stix,  _Waves in plasmas_ (Springer Science & Business Media, 1992).
  * Bittencourt (2013) J. A. Bittencourt,  _Fundamentals of plasma physics_ (Springer Science & Business Media, 2013).
  * Abramowitz and Stegun (1964) M. Abramowitz and I. A. Stegun,  _Handbook of Mathematical Functions with Formulas, Graphs, and Mathematical Tables_ , ninth dover printing, tenth gpo printing ed. (Dover, New York, 1964).