#  preprint April 2021; to appear open access in IEEE Transactions on Power
Systems
Extracting resilience metrics from distribution utility data using outage and
restore process statistics

Nichelle’Le K. Carrington, Ian Dobson, Zhaoyu Wang The authors are with ECpE
Dept., Iowa State University, Ames IA 50011 USA; emails
<EMAIL_ADDRESS>We gratefully acknowledge support in part
from NSF grants 1609080, 1735354, and 1549883. ©2021 Nichelle’Le Carrington
and Ian Dobson and Zhaoyu Wang
Creative Commons Attribution License CC BY 4.0.

###### Abstract

Resilience curves track the accumulation and restoration of outages during an
event on an electric distribution grid. We show that a resilience curve
generated from utility data can always be decomposed into an outage process
and a restore process and that these processes generally overlap in time. We
use many events in real utility data to characterize the statistics of these
processes, and derive formulas based on these statistics for resilience
metrics such as restore duration, customer hours not served, and outage and
restore rates. The formulas express the mean value of these metrics as a
function of the number of outages in the event. We also give a formula for the
variability of restore duration, which allows us to predict a maximum restore
duration with 95% confidence. Overall, we give a simple and general way to
decompose resilience curves into outage and restore processes and then show
how to use these processes to extract resilience metrics from standard
distribution system data.

###### Index Terms:

power distribution reliability, data analysis, statistics, resilience, power
system restoration

##  I Introduction

Under normal conditions, component outages in electric power distribution
systems occur at a low rate and are restored as they occur. However, when
there is severe weather or other extreme stresses, the component outages occur
at a high rate, and the outaged components accumulate until they are gradually
restored. This paper studies these resilience events in which outaged
components accumulate. We process 5 years of outage data recorded by a
distribution utility to obtain many such real events for statistical analysis.
Each event has a resilience curve that tracks the cumulative number of
components or customers outaged as they accumulate and then are restored.

It is customary to divide resilience curves into successive, non-overlapping
phases. For example, Nan [1] describes a disruptive outaging phase followed by
a recovery phase, while Panteli [2] and Ouyang [3] describe a resilience
trapezoid with the three phases of progressive disturbances, then a degraded
or assessment phase, then recovery in a resilience triangle. Instead of a
resilience triangle with a straight line hypotenuse, Reed [4] describes
restoration of power line outages with an exponential curve. Similarly, Yodo
[5] describes resilience curves in terms of successive unreliability,
disrupted, and recovery phases with triangles, trapezoids, and other curves.
Carrington [6] uses the nadir of resilience curves from utility data to
separate an outage phase from a recovery phase. Many other papers have similar
accounts of resilience phases.

These distinct phases of resilience are conceptually compelling. Moreover, the
dimensions, slopes, and areas of the resilience triangles and trapezoids
define standard resilience metrics of event duration, average rates of outage
or recovery, and overall impact.

However, working with our utility data suggests a different point of view in
which outage and restore processes routinely overlap in time. Indeed, the
average fraction of event duration for which outage and restore processes
overlap111that is, the average of $(o_{n}-r_{1})/(r_{n}-o_{1})$ in the
notation of section II is 0.61 for events with 10 to 20 outages, 0.89 for
events with 100 to 200 outages, and 0.95 for events with 1000 to 2000 outages.
That is, for practical processing of real distribution system outage data we
decompose the resilience curve not into successive phases but into outage and
restore processes that occur together for much of the time. Moreover, we will
show that standard resilience metrics can still be obtained with this
approach. Our point of view is general: outage and restore processes can
always be combined into a resilience curve, and any resilience curve can
always be decomposed into an outage process and a restore process.

We start by extracting from the utility data many resilience events and their
resilience curves and the outage and restore processes. These track the number
of components out during each event. Then we obtain the statistics of the
outage and restore processes and derive the resilience metrics. The outage
process is statistically characterized by the times between successive outages
or the outage rate. The restore process starts after a delay and is
statistically characterized by the times between successive restores or the
restore rate. The duration of the restore process and of the entire event are
then easily obtained standard metrics.

There is also a conventional customer resilience curve tracking the number of
customers out during the event that we also obtain from the utility data. This
customer resilience curve can also be decomposed into a customer outage
process and a customer restore process. We can measure the impact of an event
by the customer hours lost, which is the area under the customer resilience
curve and a well-known resilience metric [1, 2, 5]. We compute the mean
customer hours lost from the statistics of the customer outage and restore
processes.

Previous pioneering work on queueing models of reliability and resilience has
used outage and restore processes. Zapata [7] models distribution system
reliability with outages as a point process arriving at a queue that is
serviced by a repair process with multiple crews to produce an output that is
a restore process. Wei and Ji [8] analyze distribution system resilience to
particular severe hurricanes with an outage process arriving at a queue with a
repair process to produce a restore process. In [8], these processes vary in
both time and space as the hurricane progresses. Both [7] and [8]
statistically model the outage process and the repair process of components,
and then calculate the restore process. With our focus on studying the overall
system resilience with real data, we can model the restore process directly
from the data, and avoid the complexities of explicitly modeling the repair of
components and assuming an order in which they are repaired.

Our methods require a sufficient number of events for good statistics, so that
this paper addresses the more common, less extreme events. Therefore there is
little overlap of this paper with [8], which addresses individual instances of
the most extreme events (direct hits by a hurricane), for which there are few
events for a given utility.

There are several methods of estimating the number of outages in an
anticipated storm [9, 10, 11, 12, 13, 14, 15], including practical utility
application in [15]. Since some of our statistics show a dependence on the
number of outages, this capability to predict the number of outages will be
useful in applying the results in this paper to anticipated storms.

The utility data also yield estimates of the variability of the outage and
restore processes, enabling estimates of the variability of the restore
duration. Then, given the estimate of the number of outages, we can use our
restore time statistics to predict upper bounds of the restore duration of an
anticipated storm, such as its 95th percentile. The upper bound is intended to
help the utility predict when the restore process will be completed with more
confidence.

With a similar overall aim, previous work estimates individual component
restoration times from utility data in different ways. For example, Jaech [16]
predicts a gamma distribution of individual component outage restoration times
and customer hours lost with a neural network that processes utility records
and wind speed, outage time and date, and hence obtains upper bounds of
individual component restoration times. Chow [17] analyzes the contributions
of timing, faults, protection, outage types and weather to individual
component restoration times. Liu [18] fits generalized additive accelerated
failure time models to hurricane and ice storm utility data. The individual
outages were then combined to give system restoration curves at the county
level. Liu [18] and Maliszewski [19] model the spatial variation of average
outage times.

We summarize the innovations in the paper:

(1) We extract from 5 years of utility data many events in which multiple
outages accumulate before being restored and perform a novel statistical
analysis of these resilience events.

(2) We separate the events into overlapping outage and restore processes for
both components and customers, and obtain statistics of these processes. The
restore process is directly characterized from the data without modeling
repair of individual components.

(3) We derive formulas for the mean and standard deviations of restore and
event durations from the statistics of the outage and restore processes, and
illustrate their application to estimating an upper bound for the restore
duration. We derive a formula for the mean customer hours lost.

These innovations allow us to extract and compute standard metrics for
resilience events from practical utility data.

##  II Outage and restore processes

The resilience events of interest occur when outaged components accumulate
before being restored. Each event has a conventional resilience curve $C(t)$
for component outages. $C(t)$ is the negative of the cumulative number of
component outages as a function of time $t$. For example, the orange curve at
the bottom portion of Figure 2 shows $C(t)$ for an event with 10 component
outages. We now explain the outage and restore processes and how they relate
to the resilience curve.

Figure 1: A component resilience curve and its associated outage and restore
processes.

Figure 2: A customer resilience curve and its associated customer outage and
customer restore processes for the same event as Figure 2.

Figure 3: Component resilience curves (upper rows with one shaded curve) and
their corresponding decompositions into outage and restore processes (lower
rows with blue and red curves). The first example is an idealized case with
trapezoidal resilience curve and all the rest are examples from utility data.

###  II-A Examples of component outage and restore processes

We start with no components outaged. Then 10 components outage at times
$o_{1}\leq o_{2}\leq...\leq o_{10}$ as shown by the tick marks below the top
time line of Figure 2. The restore times are $r_{1}\leq r_{2}\leq...\leq
r_{10}$ as shown by the tick marks above the top time line of Figure 2. The
restore times are numbered in the time order that they occur. The cumulative
number of outages $O(t)$ at time $t$ and the cumulative number of restores
$R(t)$ at time $t$ are defined by counting $1$ for each outage or restore
before time $t$:

$\displaystyle O(t)$ $\displaystyle=\sum_{\displaystyle k\mbox{ with
}o_{k}\leq t}1$ (1) $\displaystyle R(t)$ $\displaystyle=\sum_{\displaystyle
k\mbox{ with }r_{k}\leq t}1$ (2)

Figure 2 shows the cumulative number of outages $O(t)$ and the cumulative
number of restores $R(t)$. In this case, $O(t)$ and $R(t)$ increase from zero
to the total number of outages 10.

The cumulative number of component outages at time $t$ is $O(t)-R(t)$. The
component resilience curve $C(t)$ is defined as the negative of the cumulative
number of component outages at time $t$ so that

$\displaystyle C(t)=R(t)-O(t)$ (3)

Figure 2 shows how the resilience curve $C(t)$ can be decomposed into the
restore process minus the outage process.

It is clear from (3) that any outage and restore processes $O(t)$ and $R(t)$
define a resilience curve $C(t)$. Moreover, any resilience curve $C(t)$ can be
uniquely decomposed as (3) into outage and restoration processes $O(t)$ and
$R(t)$ that increase from zero to the number of outaged components $n$. In
mathematics this decomposition is known as the Jordan decomposition [20] of
functions of bounded variation222The total variation of $C(t)$ is $2n$, which
is bounded. In our case the Jordan decomposition (3) is minimal and unique
since we require that $O(t)=R(t)=0$ for $t<o_{1}$ and $O(t)+R(t)=2n$ for
$t>r_{n}$ [21, defn. 2.4(2), thm. 2.5(2)], [22, sec. 9-4]. Since $C(t)$ is
minus the cumulative number of outages, if there are simultaneous restores and
outages, for example, $m_{r}$ restores and $m_{o}$ outages all occurring
precisely at time $t$, then only their difference $m_{r}-m_{o}$ contributes to
$C(t)$, and we assume that only $|m_{r}-m_{o}|$ contributes to the total
number $n$ of outages or restores. Note that $O(t),R(t),C(t)$ are right
continuous..

There is noticeable variety in the forms of the component resilience curves in
our utility data. The examples (except for the first example) in Figure 3 show
these curves and their decompositions into outage and restore processes. For
comparison, the first example in Figure 3 shows a conventional, idealized case
of a trapezoidal resilience curve.

Considering the outage and restore processes separately is useful because they
correspond to different aspects of system resilience: the outage process
results from individual component strengths under bad weather stress and the
restore process results from the control room, restoration plans, and the
number and performance of restoration crews.

###  II-B Extracting events from utility data

The historical outage data were recorded by one distribution utility from 2011
to 2016 over a territory including rural areas and small cities. 32 291
outages were recorded during this time. The start and end time of the outages
were recorded by fuse cards based on the loss of power at that location. The
number of customers out for each outage was also recorded.

We use our method from [6] to define events. The start of an event is defined
by an initial outage that occurs when all components are functional, and the
end of the same event is defined by the first subsequent time when all the
components are restored. That is, the event starts when the cumulative number
of failures $C(t)$ first changes from zero and ends when $C(t)$ returns to
zero.

In particular, we sort the combined component outage and restore times by
their order of occurrence and then calculate the cumulative number of outages
$C(t)$ at all the outage and restore times. Each restore at time $r_{n}$ for
which $C(r_{n})=0$ is the end of an event and the immediately following outage
is the start of the next event. Note that if the event has $n$ outages, then
it must have $n$ restores to allow the cumulative number of outage $C(t)$ to
return to zero at time $t=r_{n}$.

Applying this event processing to the historical data yields 2618 events. The
component and customer resilience curves for each event are decomposed into
outage and restore processes as explained in more detail in the next
subsection.

###  II-C Component outage and restore processes

This subsection explains the outage and restore processes in more detail and
shows how their statistics are extracted from the events in the utility data.

Figure 4: Mean and standard deviation of outage time difference empirical data
(dots) and fitted curve as a function of number of outages $n$.

Suppose that $o_{1}\leq o_{2}\leq...\leq o_{n}$ are the component outage times
in an event in order of occurrence and that $\Delta o_{k}=o_{k+1}-o_{k}$,
$k=1,...,n-1$ are the times between successive component outages. The outage
time differences $\Delta o_{k}$, $k=1,...,n-1$ can be regarded as independent
samples from a probability distribution of outage time differences $\Delta o$.
We want to find the mean and standard deviation of $\Delta o$ as a function of
$n$.

To do this, we combine the outage time differences $\Delta o_{k}$ for all the
events with $n$ outages and calculate their mean and standard deviation. Then
we fit the empirical mean and the standard deviation as functions of $n$ with
a linear combination of a constant and 2 exponential functions to smooth and
interpolate the data as shown in Figure 4. The constant plus exponentials
fitting function is chosen as a simple form to match the way the average data
decrease towards a positive value. The functional fits are

$\displaystyle\overline{\Delta o}$
$\displaystyle=7.45+23.3e^{-0.0388n}+32.2e^{-0.00391n}\quad{\rm min}$ (4)
$\displaystyle\sigma(\Delta o)$
$\displaystyle=25.6+19.5e^{-0.0375n}+30.9e^{-0.00153n}\quad{\rm min}$ (5)

Figure 5: Mean and standard deviation of restore time difference empirical
data (dots) and fitted curve as a function of number of outages $n$.

Suppose that $r_{1}\leq r_{2}\leq...\leq r_{n}$ are the component restore
times in order of occurrence. Note that the component outaged in the $k$th
outage can be different from the component restored in the $k$th restore. In
effect, we disregard which component is restored and only track that some
component is restored. (We call $r_{1},r_{2},...,r_{n}$ component restore
times to minimize any confusion with the restoration or repair times of
particular components.)

Let $\Delta r_{k}=r_{k+1}-r_{k}$, $k=1,...,n-1$ be the times between
successive component restores. The restore time differences $\Delta r_{k}$,
$k=1,...,n-1$ can be regarded as independent samples from a probability
distribution of restore time differences $\Delta r$. We extract the statistics
of $\Delta r$ from the utility data as a function of the number of outages $n$
similarly as $\Delta o$. Figure 5 plots the mean restore time difference
$\overline{\Delta r}$ and the standard deviation of the restore time
difference $\sigma(\Delta r)$ and the functions fitted. The functional fits
are

$\displaystyle\overline{\Delta r}$
$\displaystyle=7.64+30.8e^{-0.0514n}+33.8e^{-0.00391n}\quad{\rm min}$ (6)
$\displaystyle\sigma(\Delta r)$ $\displaystyle=35.3+43.7e^{-0.0224n}\quad{\rm
min}$ (7)

Let $\Delta r_{0}=r_{1}-o_{1}$ be the delay in the start of the restoration
process relative to the start of the event at time $o_{1}$. One factor
contributing to $\Delta{r_{0}}$ is utility inspection crews and clean-up crews
working to ensure the safety of the area and assess the damage needing repair.
There is no clear trend in the variation of $\Delta r_{0}$ with $n$, so we
combine the data for $\Delta r_{0}$ for all events with 2 or more outages and
compute

$\displaystyle\overline{\Delta r_{0}}=132\ {\rm
min}\qquad\mbox{and}\qquad\sigma(\Delta r_{0})=92.4\ {\rm min}$ (8)

The observed events become sparse as the number of outages $n$ increases as
shown in Table I. However, the number of samples of times between restores
$\Delta r$ and the times between outages $\Delta o$ available in the data for
number of outages $n$ is proportional to both $n-1$ and the observed number of
events with $n$ outages. Table I shows that there are more samples of $\Delta
r$ and $\Delta o$ available in the data both for the small outages with many
events and the largest outages with only one event. This improves the
estimation of $\Delta r$ and $\Delta o$ for the smallest and largest events.

TABLE I: Numbers of samples of time differences and events range of $n$ | mean # samples of | mean # events | total # events
---|---|---|---
($n$ = # outages) | $\Delta r\,\&\,\Delta o$ for each $n$ | for each $n$ | in range of $n$
2–10 | 738 | 238.8 | 2149
11–20 | 367 | 27.7 | 277
21–30 | 173 | 7.2 | 72
31–40 | 112 | 3.3 | 30
41–50 | 106 | 2.4 | 17
51–60 | 127 | 2.3 | 14
61–70 | 74 | 1.2 | 7
71–80 | 88 | 1.2 | 6
81–90 | 118 | 1.4 | 7
91–100 | 92 | 1.0 | 2
101–120 | 110 | 1.0 | 6
121–200 | 158 | 1.1 | 20
201–1000 | 364 | 1.0 | 7
1001–2000 | 1703 | 1.0 | 3

###  II-D Customer outage and restore processes

The utility data record the number of customers outaged for each outage,
allowing us to similarly analyze resilience curves tracking the number of
customers out to obtain outage and restore processes for the number of
customers. Generalizing (1) and (2), the cumulative numbers of customers out
$O^{\rm cust}(t)$ and customers restored $R^{\rm cust}(t)$ at time $t$ are

$\displaystyle O^{\rm cust}(t)$ $\displaystyle=\sum_{\displaystyle k\mbox{
with }o_{k}\leq t}c_{k}^{\rm out}$ (9) $\displaystyle R^{\rm cust}(t)$
$\displaystyle=\sum_{\displaystyle k\mbox{ with }r_{k}\leq t}c_{k}^{\rm res}$
(10)

The customer resilience curve $C^{\rm cust}(t)$ is now obtained similarly to
the component resilience curve (3) as

$\displaystyle C^{\rm cust}(t)=R^{\rm cust}(t)-O^{\rm cust}(t).$ (11)

Figure 2 shows an example of customer processes $O^{\rm cust}(t)$ and $R^{\rm
cust}(t)$ and the resilience curve $C^{\rm cust}(t)$. The processes $O^{\rm
cust}(t)$ and $R^{\rm cust}(t)$ increase from zero to the total number of
customers out, which is 181 in this example.

When an outage that disconnected customers is restored, the same number of
customers are restored. However, outages are not necessarily restored in the
order that the outages occurred. Therefore the numbers of customers restored
$c_{1}^{\rm res},c_{2}^{\rm res},...,c_{n}^{\rm res}$ are a permutation of the
numbers of customers out $c_{1}^{\rm out},c_{2}^{\rm out},...,c_{n}^{\rm
out}$.

The numbers of customers out $c_{1}^{\rm out},c_{2}^{\rm out},...,c_{n}^{\rm
out}$ can be regarded as $n$ independent samples from a distribution $c$ of
number of customers out. The customers restored $c_{1}^{\rm res},c_{2}^{\rm
res},...,c_{n}^{\rm res}$ can also be regarded as $n$ independent samples from
$c$. We combine the data for customers out for all events3332% of the customer
data are blank entries that we replaced with 0. and compute the mean and
standard deviation of the number of customers:

$\displaystyle\overline{c}=54.0\qquad\mbox{and}\qquad\sigma(c)=180.$ (12)

The customers out for each outage are determined by the location of the outage
in the network, and the distribution of the number of customers out is
determined by the overall network design and its vulnerabilities.

##  III Resilience metrics

We express event durations in terms of the time differences of the restore
process and then derive formulas for the mean and standard deviations of the
restore duration and the event duration. We also combine the time differences
with the customers outage statistics to derive a formula for the mean customer
hours lost. The average outage and restore rates are obtained.

The restore process starts at time $r_{1}$ and ends at time $r_{n}$, so the
restore duration is

$\displaystyle D_{R}$
$\displaystyle=r_{n}-r_{1}=(r_{2}-r_{1})+(r_{3}-r_{2})+...+(r_{n}-r_{n-1})$
$\displaystyle=\Delta r_{1}+\Delta r_{2}+...+\Delta r_{n-1}$ (13)

The mean restore duration is then

$\displaystyle\overline{D_{R}}$ $\displaystyle=(n-1)\overline{\Delta r}$ (14)

Assuming that $\Delta r_{1}$, $\Delta r_{2}$, …, $\Delta r_{n-1}$ are
independent, we obtain $\sigma^{2}(D_{R})=(n-1)\,\sigma^{2}(\Delta r)$ and

$\displaystyle\sigma(D_{R})$ $\displaystyle=\sqrt{n-1}\,\sigma(\Delta r)$ (15)

In (13), the restore duration $D_{R}$ is measured until the last restore of
the event. If it is preferred to measure the restore duration until, say, 95%
of the outages are restored, then this can easily be done by replacing $n-1$
in (14) and (15) by $\lceil 0.95\,n-1\rceil$, where the ceiling function
$\lceil\cdot\rceil$ rounds up to the nearest integer.

Each event starts at the first outage time $o_{1}$ and ends at the last
restore time $r_{n}$. Then the event duration is

$\displaystyle D_{E}$
$\displaystyle=r_{n}-o_{1}=(r_{1}-o_{1})+(r_{n}-r_{1})=\Delta r_{0}+D_{R}$
(16)

Using (14), the mean event duration is

$\displaystyle\overline{D_{E}}$ $\displaystyle=\overline{\Delta
r_{0}}+(n-1)\overline{\Delta r}$ (17)

and, since $\Delta r_{0}$ and $D_{R}$ are independent, we use (15) to obtain
$\sigma^{2}(D_{E})=\sigma^{2}(\Delta r_{0})+(n-1)\sigma^{2}(\Delta r)$ and

$\displaystyle\sigma(D_{E})$ $\displaystyle=\sqrt{\sigma^{2}(\Delta
r_{0})+(n-1)\sigma^{2}(\Delta r)}$ (18)

Figure 6: Resilience metrics (durations and rates) for the component outage
and restore processes.

The restore and outage rates during events444The outage rate measured over a
year is much lower than the outage rate during events because it accounts for
the time between events. are

$\displaystyle\lambda_{R}$ $\displaystyle=(\Delta r)^{-1}$ (19)
$\displaystyle\lambda_{O}$ $\displaystyle=(\overline{\Delta o})^{-1}$ (20)

In (19) we obtain the restore rate $\lambda_{R}$ from $\overline{\Delta r}$,
which is a quantity averaged over events. This restore rate $\lambda_{R}$
should be distinguished from the instantaneous restore rate $\lambda_{R}^{\rm
inst}(t)$, which has been observed in [8] to vary with time for the largest
events. The restore rate $\lambda_{R}$ averaged over events can usefully apply
even as the instantaneous restore rate varies555Consider the idealized case of
a nonhomogeneous Poisson recovery process. Suppose there are 3 restoring
processes on the time intervals $D^{(1)}$, $D^{(2)}$, $D^{(3)}$, each of
duration $T$ and with $n$ restores. Then the expected value of $\int_{t\in
D^{(i)}}\lambda_{R}^{\rm inst}(t)dt$ is $n-1$ and the duration
$T=\sum_{k=1}^{n-1}\Delta r_{k}^{(i)}=(n-1)\overline{\Delta r}^{(i)}$ for
$i=1,2,3$. Then we estimate $\lambda_{R}=\frac{1}{3T}\sum_{i=1}^{3}\int_{t\in
D^{(i)}}\lambda_{R}^{\rm inst}(t)dt$ as $\frac{3(n-1)}{(n-1)(\overline{\Delta
r}^{(1)}+\overline{\Delta r}^{(2)}+\overline{\Delta
r}^{(3)})}=(\overline{\Delta r})^{-1}$. .

The restore process depends largely on the restoration capability available to
the utility. The restore process metrics are $\overline{\Delta r_{0}}$ and
$\overline{D_{R}}$ or $\lambda_{R}$. Increasing the number of utility crews or
their effectiveness would decrease $\overline{\Delta r_{0}}$ and
$\overline{D_{R}}$, and increase $\lambda_{R}$.

The outage process depends on a combination of the weather impact and the
condition and strength of the grid components. The number of outages $n$ will
vary with the weather and the condition of the grid components, increasing if
the weather is more extreme or more prolonged, or if the grid components are
weaker.

Figure 7:  Area $A$ under resilience curve is the customer hours metric and is
equal to the area $A$ between the outage and restore processes.

We use the customer hours lost $A$ to quantify the customer impact of an
event. $A$ is the area under the customer resilience curve:

$\displaystyle A=-\int_{\displaystyle o_{1}}^{\displaystyle r_{n}}C^{\rm
cust}(t)dt$ (21)

The minus sign in (21) makes $A$ a positive area. Using (11), $A$ is also the
area between the customer outage and restore curves:

$\displaystyle A=\int_{\displaystyle o_{1}}^{\displaystyle r_{n}}\big{[}O^{\rm
cust}(t)-R^{\rm cust}(t)\big{]}dt$ (22)

The two interpretations of area $A$ are illustrated in Figure 7.

Consider the rectangle indicated by the dashed lines and the axes in Figure 7.
Let $A_{R}$ be the area in the rectangle above the restore curve and let
$A_{O}$ be the area in the rectangle above the outage curve. Then, since the
customer hours lost $A$ is the area between the outage curve and the restore
curve,

$\displaystyle A$ $\displaystyle=A_{R}-A_{O},$ (23) where $\displaystyle
A_{R}$ $\displaystyle=\bigg{(}\sum_{j=1}^{n}c_{j}^{\rm res}\bigg{)}\Delta
r_{0}+\sum_{i=2}^{n}\sum_{j=i}^{n}c_{j}^{\rm res}\Delta r_{i-1}$ (24)
$\displaystyle A_{O}$
$\displaystyle=\sum_{k=2}^{n}\sum_{\ell=k}^{n}c_{\ell}^{\rm out}\Delta
o_{k-1}$ (25)

Since the restore and outage curves are piecewise constant, the expressions
for the areas $A_{R}$ and $A_{O}$ in (24) and (25) were obtained by summing
the areas of the rectangles under the curves (recall the Riemann sum
definition of an integral).

Since the customers out are independent of the times between restores or
outages (so that expectations of products are products of expectations) and
using $\sum_{i=2}^{n}\sum_{j=i}^{n}1=\textstyle{\frac{1}{2}}n(n-1)$, we take
expectations of (24) and (25) to get

$\displaystyle\overline{A_{R}}$
$\displaystyle=\bigg{(}\sum_{j=1}^{n}\overline{c}\bigg{)}\overline{\Delta
r_{0}}+\sum_{i=2}^{n}\sum_{j=i}^{n}\overline{c}\,\overline{\Delta r}$
$\displaystyle=n\overline{c}\,\overline{\Delta
r_{0}}+{\textstyle\frac{1}{2}}n(n-1)\overline{c}\,\overline{\Delta r}$ (26)
$\displaystyle\overline{A_{O}}$
$\displaystyle=\sum_{k=2}^{n}\sum_{\ell=k}^{n}\overline{c}\,\overline{\Delta
o}={\textstyle\frac{1}{2}}n(n-1)\overline{c}\,\overline{\Delta o}$ (27)

Hence,

$\displaystyle\overline{A}$
$\displaystyle=\overline{A_{R}}-\overline{A_{O}}=n\overline{c}\,\overline{\Delta
r_{0}}+{\textstyle\frac{1}{2}}n(n-1)\overline{c}(\overline{\Delta
r}-\overline{\Delta o})$ (28)

An alternative expression for (28) can be obtained using (17):

$\displaystyle\overline{A}$
$\displaystyle=n\overline{c}\overline{D_{E}}-{\textstyle\frac{1}{2}}n(n-1)\overline{c}(\overline{\Delta
r}+\overline{\Delta o})$ (29)

The terms in (26) and (27) can be understood by examining the corresponding
areas in Figure 8. For example, the area $\overline{A_{R}}$ above the average
restore curve is the area of the rectangle with sides $\overline{\Delta
r_{0}}$ and $n\overline{c}$ plus the area of the triangle with sides
$(n-1)\overline{\Delta r}$ and $n\overline{c}$. And the area
$\overline{A_{O}}$ above the average outage curve is the area of the triangle
with sides $(n-1)\overline{\Delta o}$ and $n\overline{c}$.

Figure 8: Averaged dimensions and customer outage and restoring processes
shown to calculate the customer hours $\overline{A}$.

It is useful to describe the outage and restore processes with separate
parameters and separate metrics because they respond to different resilience
investments. For example, a program of renewing or strengthening components
would affect the outage process whereas an increased number of repair crews
would affect the restoration process. In more detail, (17) shows the effects
on the average event duration of reducing the number of outages $n$ (by
hardening the infrastructure) reducing $\overline{\Delta r_{0}}$ (by deploying
more inspection crews) and reducing the average time between restores
$\overline{\Delta r}$ (by deploying more repair crews). For a larger event,
$n$ is larger and reducing $\overline{\Delta r}$ will have a larger effect
because $\overline{\Delta r}$ is multiplied by $n-1$. Formula (28) shows the
corresponding effects on the customer hours $\overline{A}$. Hardening the
upstream system or installing more reclosers can reduce the average customers
disconnected per outage $\overline{c}$ and proportionally reduce
$\overline{A}$. Reducing the number of outages $n$ or $\Delta r$ has an even
greater effect for larger events because the second term of (28) grows like
$n^{2}$.

##  IV Results

This section gives numerical results illustrating the application of the
formulas for the statistics of the metrics.

We can evaluate restore duration mean $\overline{D_{R}}$ and standard
deviation $\sigma(D_{R})$ for a given number of outages from (6) and (7). For
example, if there are $n=10$ outages then the restore duration has mean 527
min and standard deviation 211 min. If there are $n=100$ outages, then the
restore duration has mean 3038 min and standard deviation 397 min. If these
formulas are to be used for predicting restoration duration for an incoming
storm, then the number of outages $n$ can be predicted by a number of methods
as reviewed in the introduction. As well as estimating the mean, it is useful
in applying the restore duration to compute its variability with its standard
deviation.

The event duration $D_{E}$ is the restore duration $D_{R}$ plus the delay
until the first restore $\Delta r_{0}$. From (8), $\Delta r_{0}$ has mean 132
min and standard deviation 92.4 min. Then we can evaluate event duration mean
and standard deviation from (17) and (18). For example, if there are $n=10$
outages, then the event duration has mean 660 min and a standard deviation of
230 min. If there are $n=100$ outages, then the event duration has mean 3171
min and a standard deviation of 408 min.

When an outage event occurs, the primary question the customers want an answer
for is “How long will the power be out?”. To help determine what should be
announced to the public to answer this question, it is useful to estimate an
upper bound on the restore duration that will be satisfied with a specified
confidence level.

For each given value of number of outages $n$, the restore duration $D_{R}$
approximately follows a gamma distribution. We can estimate from (6) and (7)
the mean and standard deviation of $D_{R}$ and then calculate the gamma
distribution with that mean and standard deviation. This method of estimation
is called the method of moments. We now give more details. The gamma
distribution has probability density function

$\displaystyle f(x)=\frac{\beta^{\alpha}}{\Gamma(\alpha)}x^{\alpha-1}e^{\beta
x},\quad x\geq 0,$ (30)

where $\Gamma(\cdot)$ is the gamma function, $\alpha$ is the shape parameter
and $\beta$ is the rate parameter. The gamma distribution mean is
$\alpha/\beta$ and the standard deviation is $\sqrt{\alpha}/\beta$. For each
$n$, we can evaluate (6) and (7) and solve

$\displaystyle\overline{D_{R}}$ $\displaystyle=\frac{\alpha}{\beta}$ (31)
$\displaystyle\sigma(D_{R})$ $\displaystyle=\frac{\sqrt{\alpha}}{\beta}$ (32)

for $\alpha$ and $\beta$. These parameters specify a gamma distribution
$f(D_{R})$ for each value of $n$. Then we can easily evaluate the 95th
percentile $D_{R95}$ of the gamma distribution that satisfies
$\int_{0}^{D_{R95}}f(x)dx=0.95$. $D_{R95}$ estimates an upper bound on the
restore duration that is satisfied by the actual restore duration with
probability 0.95.

The curves in Figure 9 show an increasing and initially decelerating increase
of mean restore duration as the number of outages increase, and a similar
increase in the 95th percentile of restore duration. The dots in Figure 9 are
the restore durations for the events in the data; they show how the mean and
95th percentile of restore duration calculated from the estimated gamma
distribution summarize the empirical data. The event data become sparser as
the number of outages increase.

Figure 9: Curves show mean and 95th percentile of restore duration $D_{R}$
versus number of outages. Dots show the restore durations of events in the
data.

Figure 10 shows the outage rate $\lambda_{O}$ and the restore rate
$\lambda_{R}$ obtained from (20) and (19) as the number of outages varies.
Both rates increase significantly as the number of outages increase. The
outage rate results from the interaction of the weather with the grid, and
depends on the design margin, age, and maintenance of the grid components. The
restore rate results from the utility crews, restoration plans, and control
room procedures. For up to 250 outages, the restore rate is quite close to the
outage rate. The restore rate slightly lags the outage rate, showing the
extent to which the utility restore process succeeds in keeping up with the
outage process.

Figure 10: Outage rate $\lambda_{O}$ and restore rate $\lambda_{R}$ versus
number of outages.

The curve in Figure 11 shows the mean customer hours $\overline{A}$ calculated
from (28) increasing as a function of the number of outages. The dots in
Figure 11 show the customer hours $A$ for each event in the data. There is
considerable variability in the customer hours in the data for more than 100
outages. Future work could aim to analyze and quantify this variability.

Although sections II-C and II-D fit the data for the full range of our data up
to 2000 outages, the data for the events with more than 250 outages become
sparse and more variable. In order to be cautious in our conclusions, we limit
all the presented results in this section to events with up to 250 outages.
Future work with more data or with more elaborate statistical methods might
well extend the range of prediction.

Figure 11: Curve shows mean customer hours $\overline{A}$ calculated from (28)
versus number of outages. Dots show customer hours $A$ of the events in the
data.

##  V Discussion

This section discusses the contributions of the paper in the larger context of
power system resilience.

Our analysis is for distribution systems since the data are from a
distribution system. However, the decomposition of resilience curves into
outage and restore processes is general and could be applied to resilience
curves derived from data in transmission systems or in microgrids. (Resilience
curves for transmission systems are discussed in many papers, including [1, 2,
23, 24, 3], and resilience curves for microgrids are discussed in [25, 26,
27].) Note that the decomposition applies to resilience curves that track any
quantity on the vertical axis; this paper is restricted to tracking numbers of
components and customers simply because these are the quantities available in
our data.

Detailed outage data for transmission systems usually track number of
components, and the corresponding outage and restore processes and associated
metrics could be extracted. However, the details of the event definition and
the specific forms of the processes could differ. For example, we would expect
the outage and restore processes for transmission systems to have some
reduction in overlap, and their restore processes to often have an exponential
recovery form. There are also other approaches to quantifying resilience in
transmission systems based on detailed outage data. Resilience metrics are
obtained by sampling transmission system outage data statistics in [28], and
resilience is assessed by sampling from a Markovian influence graph driven by
transmission system outage data in [29].

There is extensive work with simulation, mitigation, optimization, and
calculation of metrics based on detailed models of phases of resilience; for
example, [32, 33, 31, 30, 34, 3]. We see this detailed modeling as
complementary to methods driven by historical data. For example, to reduce
customer hours lost, one could use formula (25) to identify how the mean time
between restores influences mean customer hours lost and then use a detailed
simulation model of the restoration to show how mean time between restores
could be reduced.

Our approach analyses resilience events with metrics at a systems level. The
mechanisms of outage and restoration of individual components are complicated,
and how individual outages and restorations and component design and
maintenance and weather and restoration procedures interact is also
complicated. The outage and restore processes describe the aggregated outcomes
of these complexities by tracking the number of components or customers
outaged or restored. Note that which individual component is outaged or
restored or the order of outage or restore is not tracked by the outage and
restore processes. We find that it is simpler and more direct to avoid
description of individual component repair processes, and this approach is
made possible by the utility data that directly describe the outage and
restore processes.

Instead of analyzing individual resilience events, we analyze the statistics
of many resilience events, as do [6, 17, 7, 18, 16]. It is certainly useful to
analyze individual resilience events (e.g., [4, 35]) and respond by proposing
upgrades that would mitigate similar events. However, an exclusive focus on
individual events may skew analysis and mitigation towards the last major
event; quantifying resilience more broadly across many events is a useful
complement. Improving resilience is inherently a hard problem since the
largest events with highest impact and significant risk are rare, and we
conclude that analyses of both individual events and the statistics of all
events is warranted.

Smaller distribution utilities may not have enough outages over a period of
years to extract good statistics. In that case, data could be combined from
nearby distribution utilities with similar design practices and weather (for
outage process metrics) or similar or shared restoration crews and procedures
(for restore process metrics). Also, government or reliability authorities
could aggregate utility data to quantify resilience over a region to inform
investments in resilience.

While the practical advantages such as realism and lack of modeling
assumptions of analyzing historical data are clear, the limitations should
also be noted. The historical data are what did happen and a pragmatic sample
of what could happen, but cannot have many samples of all the possible rare
events. Moreover, some of the extremes of weather driving resilience events
will become gradually more frequent and more severe as the planet warms [36],
so that the severity of some types of resilience events will increase. This
underlines the importance of developing multiple approaches to quantifying and
improving resilience, including the approach developed in this paper.
Extracting resilience metrics from utility data and thus quantifying aspects
of resilience can help to justify resilience improvements to mitigate
disasters.

##  VI Conclusions

We process 5 years of distribution system outage data to extract and study
many resilience events in which outages accumulate and are restored. As
appropriate for quantifying resilience, we focus only on the resilience
events, and do not analyze the frequency of these events or the times between
events that are of interest in other kinds of reliability analysis.

It is usual to separate resilience curves for events into successive non-
overlapping phases in time such as outage and recovery. However, our
distribution utility data show that outage and restore processes typically
occur together for most of the event. Therefore, instead of using successive
phases, we show how resilience curves tracking the number of outaged
components or customers can be easily decomposed into outage and restore
processes that can occur at the same time. These outage and restore processes
describe the same information as the resilience curve, but usefully correspond
to different aspects of resilience: the outages are caused by weather
interacting with the strength of the components, whereas the restores are done
by utility crews and control rooms. The decomposition of the resilience curve
into outage and restore processes is known as the Jordan decomposition in
mathematics.

We compute some basic statistics of the outage and restore processes. In
particular, we estimate fits for the mean and standard deviations of the times
between outages and the times between restores as functions of the number of
the outages. The function fitting has the effect of smoothing and
interpolating the noisy data. We also estimate the means and standard
deviations of the number of customers outaged and of the delay until the
restore process starts.

Then, given the predicted number of outages, which is estimated in several
ways by previous work [9, 10, 11, 12, 13, 14, 15], we obtain formulas for the
means of standard resilience metrics, such as restore and event durations,
restore and outage rates, and the customer hours lost. These formulas for
standard resilience metrics usefully quantify the resilience processes and
show how the metrics depend on the number of outages, the delay before
restoration starts, the average time between restores, and the average number
of customers disconnected per outage. The outage rate quantifies the overall
grid fragility under weather stress and the restore rate quantifies the
overall performance of utility crews. This quantification can help inform
investments that improve these metrics. We also estimate the standard
deviations of the restore duration and the event duration. This leads to
estimates of probable upper bounds of restore durations. These credible upper
bounds based on past performance should be useful to utilities for informing
customers about their outage duration with confidence. The utility data
available to us limited our resilience processes to counting outages of
components or customers. If available, quantities such as power outaged could
be similarly analyzed to obtain useful metrics such as energy not served.

Our approach models the outage and restore processes directly from utility
data and avoids the complexities of modeling individual component restoration
or repair times and their order of restoration. That is, since the utility
data itself incorporate the detailed complexities of resilience, we can give a
high-level description and quantification of resilience. This high-level data-
driven approach is very much a useful complement to the detailed modeling of
the restoration processes by other authors.

In summary, we extract and separate the outage and restore processes from
distribution utility outage data in a new way, estimate the statistics of
times between successive restores or outages, and then show how standard
resilience metrics can be derived from these statistics. The overall effect is
to compute some useful resilience metrics from practical utility data.

## References

  * [1] C. Nan and G. Sansavini, “A quantitative method for assessing resilience of interdependent infrastructures,” Reliab. Eng. Syst. Safety, vol. 157, 2017, pp. 35–53.
  * [2] M. Panteli, D. N. Trakas, P. Mancarella, and N. D. Hatziargyriou, “Power systems resilience assessment: hardening and smart operational enhancement,” Proc. IEEE, vol. 105, no. 7, May 2017, pp. 1202–1213.
  * [3] M. Ouyang, L. Dueñas-Osorio, and X. Min, “A three-stage resilience analysis framework for urban infrastructure systems,” Structural Safety, vol. 36, 2012, pp. 23–31.
  * [4] D. A. Reed, K. C. Kapur, and R. D. Christie, “Methodology for assessing the resilience of networked infrastructure,” IEEE Syst. J., vol. 3, no. 2, Jun. 2009, pp. 174–180.
  * [5] N. Yodo and P. Wang, “Engineering resilience quantification and system design implications: a literature survey,” ASME J. Mechanical Design, vol. 138, Nov. 2016, 111408, pp. 1–13.
  * [6] N. K. Carrington, S. Ma, I. Dobson, and Z. Wang, “Extracting resilience statistics from utility data in distribution grids,” in Proc. IEEE PES General Meeting, Montreal, Canada, Aug. 2020.
  * [7] C. J. Zapata, S. C. Silva, H. I. Gonzalez, O. L. Burbano, and J. A. Hernandez, “Modeling the repair process of a power distribution system,” in Proc. IEEE/PES T&D Conf. & Exp.: Latin America, 2008.
  * [8] Y. Wei, C. Ji, F. Galvan, S. Couvillon, G. Orellana, and J. Momoh, “Non-stationary random process for large-scale failure and recovery of power distributions,” Applied Mathematics, vol. 7, Feb. 2016, pp. 233–249.
  * [9] H. Liu, R. A. Davidson, and T. V. Apanasovich, “Statistical forecasting of electric power restoration times in hurricanes and ice storms,” IEEE Trans. Power Systems, vol. 22, 2007, pp. 2270–2279.
  * [10] H. Liu, R. A. Davidson, D. V. Rosowsky, and J. R. Stedinger, “Negative binomial regression of electric power outages in hurricanes,” J. Infrastructure Systems, vol. 11, no. 4, 2005, pp. 258–267.
  * [11] P. Kankanala, S. Das, and A. Pahwa, “AdaBoost+: An ensemble learning approach for estimating weather-related outages in distribution systems,” IEEE Trans. Power Systems, vol. 29, no. 1, Jan. 2014, pp. 359–367.
  * [12] D. Zhu, D. Cheng, R. P. Broadwater, and C. Scirbona, “Storm modeling for prediction of power distribution system outages,” Elect. Power Syst. Res., vol. 77, no. 8, 2007, pp. 973–979.
  * [13] Y. Zhou, A. Pahwa, and S. S. Yang, “Modeling weather-related failures of overhead distribution lines,” IEEE Trans. Power Systems, vol. 21, no. 4, Nov. 2006, pp. 1683–1690.
  * [14] K. Alvehag and L. Söder, “A reliability model for distribution systems incorporating seasonal variations in severe weather,” IEEE Trans. Power Delivery, vol. 26, no. 2, Apr. 2011, pp. 910–919.
  * [15] H. Li, L. A. Treinish, and J. R. M. Hosking, “A statistical model for risk management of electric outage forecasts,” IBM J. Research and Development, vol. 54, no. 3, paper 8, May/Jun. 2010.
  * [16] A. Jaech, B. Zhang, M. Ostendorf, and D. S. Kirschen, “Real-time prediction of the duration of distribution system outages,” IEEE Trans. Power Systems vol. 34, Jan. 2019, pp. 773–781.
  * [17] M-Y. Chow, L. S. Taylor, and M-S. Chow, “Time of outage restoration analysis in distribution systems,” IEEE Trans. Power Delivery, vol. 11, no. 3, Jul. 1996, pp. 1652–1658.
  * [18] H. Liu, R. A. Davidson, and T. V. Apanasovich, “Spatial generalized linear mixed models of electric power outages due to hurricanes and ice storms,” Reliab. Eng. Syst. Safety, vol. 93, no. 6, 2008, pp. 897–912.
  * [19] P. J. Maliszewski and C. Perrings, “Factors in the resilience of electrical power distribution infrastructures,” Appl. Geogr., vol. 32, no. 2, 2012, pp. 668–679.
  * [20] Jordan decomposition (of a function). Encyclopedia of Mathematics. [Online]. Available: https://
encyclopediaofmath.org/wiki/Jordan_decomposition_(of_a_function)

  * [21] T. Jafarikhah and K. Weihrauch, “Computable Jordan decomposition of linear continuous functionals on C[0;1],” Logical Methods in Computer Science, vol. 10, Sept. 2014, pp. 1–13.
  * [22] A.E. Taylor, General theory of functions and integration, Dover NY 1985.
  * [23] M. Panteli, P. Mancarella, D. N. Trakas, E. Kyriakides, and N. D. Hatziargyriou, “Metrics and quantification of operational and infrastructure resilience in power systems,” IEEE Trans. Power Systems, vol. 32, no. 6, Nov. 2017, pp. 4732–4742.
  * [24] J. Lu, J. Guo, Z. Jian, Y. Yang, and W. Tang, “Resilience assessment and its enhancement in tackling adverse impact of ice disasters for power transmission systems,” Energies, vol. 11, no. 9, Aug. 2018, 2272. [Online]. Available: https://www.mdpi.com/1996-1073/11/9/2272
  * [25] Q. Zhou, M. Shahidehpour, A. Paaso, S. Bahramirad, A. Alabdulwahab, and A. Abusorrah, “Distributed control and communication strategies in networked microgrids,” IEEE Communications Surveys & Tutorials, vol. 22, no. 4, fourth quarter 2020, pp. 2586–2633.
  * [26] W. W. J. Anderson, “Resilience assessment of islanded renewable energy microgrids,” PhD dissertation, Naval Postgraduate School, Monterey CA, USA, 2020. [Online]. Available: https://calhoun.nps.edu/handle/10945/66574
  * [27] S. Ma, N. K. Carrington, A. Arif, and Z. Wang, “Resilience assessment of a self-healing distribution system under extreme weather events,” in Proc. IEEE PES General Meeting, Atlanta GA USA, Aug. 2019.
  * [28] M. R. Kelly-Gorham, P. D. H. Hines, K. Zhou, and I. Dobson, “Using utility outage statistics to quantify improvements in bulk power system resilience,” Elect. Power Syst. Res., vol. 189, Dec. 2020, 106676.
  * [29] K. Zhou, I. Dobson, and Z. Wang, “Can the Markovian influence graph simulate cascading resilience from historical outage data?,” in Proc. Prob. Methods Applied to Power Sys., Liège, Belgium, Aug. 2020.
  * [30] N. Bhusal, M. Abdelmalak, M. Kamruzzaman, and M. Benidris, “Power system resilience: current practices, challenges, and future directions,” IEEE Access, vol. 8, Jan. 2020, pp. 18064–18086.
  * [31] Z. Bie, Y. Lin, G. Li, and F. Li, “Battling the extreme: A study on the power system resilience,” Proc. IEEE, vol. 105, no. 7, Jul. 2017, pp. 1253–1266.
  * [32] S. Ma, L. Su, Z. Wang, F. Qiu, and G. Guo, “Resilience enhancement of distribution grids against extreme weather events,” IEEE Trans. Power Systems, vol. 33, no. 5, Sept. 2018, pp 4842–4853.
  * [33] J. Xu, R. Yao, and F. Qiu, “Mitigating cascading outages in severe weather using simulation-based optimization,” IEEE Trans. Power Systems, vol. 36, no. 1, Jan. 2021, pp. 204–213.
  * [34] Y. Wang, C. Chen, J. Wang, and R. Baldick, “Research on resilience of power systems under natural disasters—a review,” IEEE Trans. Power Systems, vol. 31, Mar. 2016, pp. 1604–1613.
  * [35] Y. Wei, C. Ji, F. Galvan, S. Couvillon, G. Orellana, and J. Momoh, “Learning geotemporal nonstationary failure and recovery of power distribution,” IEEE Trans. Neural Networks and Learning Systems, Jan. 2014, pp. 229–240.
  * [36] Intergovernmental Panel on Climate Change [Core Writing Team, R. K. Pachauri and L. A. Meyer (eds.)], IPCC, 2014: Climate Change 2014: Synthesis Report. Contribution of Working Groups I, II and III to the Fifth Assessment Report of the Intergovernmental Panel on Climate Change. IPCC, Geneva, Switzerland, 2014.

|  Nichelle’Le K. Carrington (S‘18) received the B.S. degree in electrical
engineering in 2014 from North Carolina Agricultural and Technical State
University, Greensboro NC, USA. She is currently pursuing the Ph.D. in
electrical engineering with a focus in power and energy systems at Iowa State
University. Her current research interests include distribution system
resilience, analysis and assessment of real utility data to extract resilience
metrics, data analytics in smart grids, and industry applications leveraging
power system smart meter data.
---|---
Ian Dobson (F’06) received the B.A. degree in mathematics from Cambridge
University and the Ph.D. degree in electrical engineering from Cornell
University. He previously worked for British industry and the University of
Wisconsin-Madison, and is currently Sandbulte Professor of Electrical
Engineering at Iowa State University. His interests are blackout risk, power
system stability, synchrophasors, complex systems, and nonlinear dynamics.
---
Zhaoyu Wang (S’13-M’15-SM’20) is the Harpole-Pentair Assistant Professor with
Iowa State University. He received the B.S. degree from Shanghai Jiaotong
University, and Ph.D. degree from Georgia Institute of Technology. Dr. Wang
has been a recipient of the National Science Foundation (NSF) CAREER Award and
the IEEE PES Outstanding Young Engineer Award. His research interests include
optimization and data analytics in power distribution systems and microgrids.
---


# The Adiabatic Phase Mixing and Heating of Electrons in Buneman Turbulence

H. Che1, J. F. Drake2, M. Swisdak2, M. L. Goldstein1 1:Goddard Space Flight
Center, NASA, Greenbelt, MD, 20771, USA  2: IREAP, University of Maryland,
College Park, MD, 20742, USA

###### Abstract

The nonlinear development of the strong Buneman instability and the associated
fast electron heating in thin current layers with $\Omega_{e}/\omega_{pe}<1$
are explored. Phase mixing of the electrons in wave potential troughs and a
rapid increase in temperature are observed during the saturation of the
instability. We show that the motion of trapped electrons can be described
using a Hamiltonian formalism in the adiabatic approximation. The process of
separatrix crossing as electrons are trapped and de-trapped is irreversible
and guarantees that the resulting electron energy gain is a true heating
process.

††preprint: APS/123-QED

The exploration of how waves and particles interact in strong turbulence has
been an important challenge in plasma physics. Kadomtsev (1965); O’Neil
(1965); Dupree (1966); Sagdeev and Galeev (1969); Galeev  _et al._ (1975);
Goldman (1984); Krommes (2002); Yampolsky and Fisch (2009); Bénisti, Morice,
and Gremillet (2012) Using particle-in-cell simulations we explore the
nonlinear development and nonlinear wave-particle interactions of the Buneman
instability to reveal how particle acceleration and heating take place. The
Buneman instabilityBuneman (1958) is driven by the relative drift between ions
and electrons. Its quasi-linear theory is well understood, but strong Buneman
turbulence is still a subject with open questions though it has been widely
discussedDavidson  _et al._ (1970); Ishihara, Hirose, and Langdon (1981);
Hirose, Ishihara, and Langdon (1982); Cargill and Papadopoulos (1988). The
previous work either did not consider the trapping regime (where the wave
electric field is large enough to trap thermal particles) or treated it under
the assumption that the particle heating growth rate was slow compared to the
instability. We investigate the regime in which rapid electron heating takes
place near the saturation of the Buneman instability when the particle’s
bounce rate in the wave potential is far larger than the growth rate of the
instability. As a consequence, the trapped particle’s motion is approximately
adiabatic. Heating is thus a consequence of coherent trapping, phase mixing
and de-trapping of the particles. Our simulations also demonstrate the
difference between the nonlinear development of the Buneman instability and an
idealized adiabatically-growing single sine wave, which supports that the
heating can be achieved by adiabatic motion and de-trapping .

Electron heating as a result of the Buneman instability is associated with the
intense electron current layers formed during magnetic reconnectionDrake  _et
al._ (2003); Che  _et al._ (2010); Khotyaintsev  _et al._ (2010),
shocksCargill and Papadopoulos (1988); Riquelme and Spitkovsky (2009);
Matsumoto, Amano, and Hoshino (2012) and turbulent energy cascades to sub-
proton scalesAlexandrova  _et al._ (2009); Sahraoui  _et al._ (2010). In
particular, understanding how kinetic turbulence transfers momentum and energy
is important for revealing the role of anomalous resistivity in magnetic
reconnection, which has been a long-standing puzzleKulsrud  _et al._ (2005);
Che, Drake, and Swisdak (2011).

We propose a new mechanism that is responsible for extremely fast electron
heating, in a few tens of electron plasma periods, during the nonlinear
evolution of the Buneman instability. The dynamics is dominated by the
coherent trapping and de-trapping of streaming electrons (with drift $v_{de}$)
in the nearly non-propagating electric field from the instability. The wave
amplitude grows until nearly all of the streaming electrons have been trapped.
Thus the electrostatic potential at saturation is approximately given by
$e\phi\sim m_{e}v_{de}^{2}$. The bounce frequency
$\omega_{b}=k_{0}\sqrt{e\phi/m_{e}}\sim\omega_{pe}$ of electrons trapped in
the potential $\phi$ greatly exceeds the characteristic growth rate of the
wave $\gamma\sim(m_{e}/m_{i})^{1/3}\omega_{pe}$. As a result, the electrons
trapped in the growing potential behave adiabatically, preserving their phase
space area as the wave amplitude slowly changes in time. Phase mixing of the
electrons in the wave potential troughs guarantees that, as the wave amplitude
decreases following saturation, the de-trapping of electrons leaves a
distribution of particles that forms a velocity-space plateau over the
interval ($-v_{de},v_{de}$). The process of separatrix crossing as electrons
are trapped and de-trapped is irreversible and guarantees that the resulting
electron energy gain $\sim m_{e}v_{de}^{2}/2$ is a true heating process.

We carry out 3D PIC simulations with strong electron drifts in an
inhomogeneous current-carrying plasma with a guide field. We apply no external
perturbations to initiate reconnection, and consequently reconnection does not
develop during relatively short duration of the simulation. We specify the
initial magnetic field to be $B_{x}/B_{0}=\tanh[(y-L_{y}/2)/w_{0}]$, where
$B_{0}$ is the asymptotic amplitude of $B_{x}$, and $w_{0}$ and $L_{y}$ are
the half-width of the initial current sheet and the box size in the $y$
direction, respectively. The guide field $B_{z}^{2}=B^{2}-B_{x}^{2}$ is chosen
so that the total field $B$ is constant. In our simulation, we take
$w_{0}=0.1d_{i}$ and $B=\sqrt{26}B_{0}$, where $d_{i}=c/\omega_{pi}$ and
$\omega_{pi}$ is the ion plasma frequency. The initial temperature is
$T_{e0}=T_{i0}=0.04m_{i}c_{A0}^{2}$, the mass ratio is $m_{i}/m_{e}=100$, and
$\Omega_{e}\sim 0.625\omega_{pe}$, where $c_{A0}=B_{0}/(4\pi
n_{0}m_{i})^{1/2}$ is the asymptotic ion Alfvén wave speed. The simulation
domain has dimensions $L_{x}\times L_{y}\times L_{z}=0.5\times 1\times 4$
$d_{i}$ with periodic boundaries in $x$ and $z$ and a conducting boundary in
$y$. The initial electron drift along $z$, $v_{de}\sim 10c_{A0}$, is above the
threshold for triggering the Buneman instability.

Figure 1: a: The time evolution of $d\langle T_{ezz}\rangle/dt$ (black solid
line), $|dW/dt|=-\dfrac{m_{e}}{e^{2}}d\langle j_{ez}^{2}/n_{e}^{2}\rangle/dt$
(black dot-dot-dashed line), and the parallel electric field variance
$Var(E_{z})=\langle\sqrt{(E_{z}-\langle E_{z}\rangle)^{2}}\rangle$ (blue dot-
dashed line), where $n_{e}$ is the electron density and $j_{ez}$ is the
$z$-component of the electron current density. b: Electron velocity
distribution functions in the current sheet $f(v_{ez})$ at $\Omega_{e}t=25.5$,
38.2, 51, 63.7 and 102 plotted as black dot-dashed, red short-dashed, orange
long-dashed, blue dot-dot-dashed and black solid lines respectively. The
narrow ion velocity distribution function $f(v_{iz})$ at $\Omega_{e}t=102$ is
reduced by 5 times. c: Electric field $E_{z}$ parallel to $B$ at
$\Omega_{e}t=$25.5 (blue dot-dashed), 38.2 (black solid), 63.7 (red dashed).

During the simulation the Buneman instability onsets around $\Omega_{e}t\simeq
25.5$ with a wave vector that is aligned along the magnetic field. The
instability reaches its peak around $\Omega_{e}t\simeq 40$ and ceases around
$\Omega_{e}t\simeq 102$, as indicated by the turbulence strength $Var(E_{z})$
in Fig.1a (blue dot-dashed line). The electric field $E_{z}$ parallel to
$B_{z}$ abruptly increases from a few $E_{0}=c_{A0}B_{0}/c$ to $E_{z}\sim
40-60E_{0}$ at $\Omega_{e}t\sim 40$ and then falls to a value $20E_{0}$ at
$\Omega_{e}t\sim 64$ (Fig.1c). At the same time, the average parallel
component of the electron temperature, $\langle T_{ezz}\rangle$ sharply
increases, from 0.04 to 0.5, by more than a factor of 10 while the ion
temperature increases only slightly. $\langle\rangle$ denotes an average over
the mid-plane of the current sheet at $y=L_{y}/2$. The electron drift velocity
decreases from $9c_{A0}$ to $\sim 6c_{A0}$. It is noteworthy in Fig. 1a that
the increase of $\langle T_{ezz}\rangle$ nearly matches the damping rate of
the electron parallel kinetic energy $W$, which implies that the streaming
kinetic energy is nearly fully converted into thermal energy, i.e.
$\dfrac{m_{e}}{2}\triangle\langle
v_{de}^{2}\rangle\approx\dfrac{1}{2}\triangle\langle T_{ezz}\rangle$, where
the Boltzmann constant has been absorbed into T. Panel b in Fig. 1 shows the
electron velocity distribution function in the current sheet at
$\Omega_{e}t=25.5$, 38.3, 51, 63.7 and 102. The narrow ion distribution
function at $\Omega_{e}t=102$ is shown with a solid line. We can see that the
electron velocity distribution functions $f(v_{ez})$ become flatter and
broader at late times, but the significant change takes place during
$\Omega_{e}t\sim 38-51$. The electron distribution functions reveal that a few
electrons are accelerated to very high velocity, which is a consequence of the
inductive electric field $E_{z}$ that maintains the integrated current. The
ion velocity distribution function is slightly affected by the turbulence.

Figure 2: Electron distribution functions in phase space $(z,v_{ez})$ at
$\Omega_{e}t=$25.5, 38.2, 51 and 63.7.

To fully understand the fast electron thermalization, we plot the electron
distribution function in phase space $(z,v_{ez})$ at times $\Omega_{e}t=$
25.5, 38.2, 51, and 63.7 in Fig. 2. We see that at $\Omega_{e}t=38.2$,
electrons begin their trapped orbits in the potential wells and by
$\Omega_{e}t=51$ they have phase mixed along their trapped orbits. The
localized electric field structure remains intact. Phase mixing does not
change significantly after $\Omega_{e}t=63.7$. The period of fast phase mixing
is coincident with that of the flattening of electron distribution functions
shown in Fig. 1. The phase mixing occurs near the time of saturation when the
change in the electric field is small. We now explore the physical mechanism
that produces the fast phase mixing and electron heating.

To reveal the physics behind this phase mixing and electron heating, we need
to identify the properties of the instability driving the turbulence. We use a
double drifting-Maxwellian kinetic model to trace the evolution of the
instabilities in the simulation in the following way Che (2009); Che  _et al._
(2009, 2010). We fit the ion and electron distribution functions averaged over
$z$ at $x=0.025$ in the mid-plane $x-z$ of the current sheet at $\Omega_{e}t=$
25.5, 38.2, 51, and 63.7, and then substitute the fits into the local
dispersion relation derived from a double drifting-Maxwellian kinetic model
for waves with $\Omega_{i}\ll\omega\ll\Omega_{e}$:

$1+\frac{2\omega_{pi}^{2}}{k^{2}v^{2}_{ti}}[1+\zeta_{i}Z(\zeta_{i})]+\frac{2\delta\omega^{2}_{pe}}{k^{2}v^{2}_{zte1}}[1+I_{0}(\lambda)e^{-\lambda}\zeta_{e1}Z(\zeta_{e1})]\\\
+\frac{2(1-\delta)\omega^{2}_{pe}}{k^{2}v^{2}_{zte2}}[1+I_{0}(\lambda)e^{-\lambda}\zeta_{e2}Z(\zeta_{e2})]=0,$
(1)

where $\zeta_{i}=(\omega-k_{z}v_{di})/kv_{ti}$, $\zeta_{e1}=(\omega-
k_{z}v_{de1})/k_{z}v_{zte1}$, $\zeta_{e2}=(\omega-
k_{z}v_{de2})/k_{z}v_{zte2}$, $\lambda=k^{2}_{x}v^{2}_{xte}/2\Omega^{2}_{e}$,
$Z$ is the plasma dispersion function and $I_{0}$ is the modified Bessel
function of the first kind with order zero. The thermal velocity of species
$j$ is defined by $v^{2}_{tj}=2T_{tj}/m_{j}$ and the drift speed by $v_{dj}$
which is parallel to the $z$ direction. The electron temperature takes on
different values parallel and across the magnetic field while the ions are
assumed to be isotropic. $\delta$ is the weight of the high velocity drifting
Maxwellian.

We numerically solve the dispersion relation in Eq. 1 and obtain the unstable
modes at $\Omega_{e}t=$ 25.5, 38.2, 51, 63.7 and 102. We find that the Buneman
instability dominates. The growth rate of the fastest growing Buneman mode
decreases with time from $\gamma\sim 0.12\omega_{pe}$ (close to the linear
value given by Ishihara et al Ishihara, Hirose, and Langdon (1981)
$\gamma\sim\sqrt{3}/2(m_{e}/2m_{i})^{1/3}(1-(m_{e}/2m_{i})^{1/3}/2)\omega_{pe}\sim
0.13\omega_{pe}$) at $\Omega_{e}t=25.5$ to $\gamma\sim 0.06\omega_{pe}$ at
$\Omega_{e}t=38.2$ and $\gamma\sim 0.006\omega_{pe}$ at $\Omega_{e}t=102$. The
frequency of the fastest growing mode is about $\omega_{0}\sim
0.013\omega_{pe}$ and its wavenumber $k_{0}d_{i}$ decreases from 90 to 75. The
phase speed $v_{p}$ increases slightly with time and has a value of $v_{p}\sim
0.05c_{A0}$. A transient two stream instability with growth rate $\gamma\sim
0.006\omega_{pe}$ appears at $\Omega_{e}t=51$ and is stable by
$\Omega_{e}t=63.7$. An oblique lower hybrid instability develops with growth
rate $\gamma\sim 0.02\omega_{pe}$ after $\Omega_{e}t=51$.

Figure 3: Panel (a) and (b) are electron distributions in phase space
$z,v_{ez}$ at the middle and the end time of the test particle simulation with
$E_{max}\sim 80$. The black dashed lines in Panel (c) and (d) are the initial
test electron velocity distribution functions. Panel (c) is at the middle time
of the test particle simulations and panel (d) is at the end time of the
simulations. The black solid lines are for the test particle simulation with
$E_{max}\sim 80$ and the red dot-dashed line are for the test particle
simulation with $E_{max}\sim 40$.

It’s interesting to notice that during $\Omega_{e}t\sim 38-64$, the typical
parallel electric field is about $40E_{0}$ and the wavenumber of the fastest
mode is $k_{0}d_{i}\sim 90$. The corresponding bounce frequency is
$\omega_{b}=k_{0}v_{b}/\sqrt{2}\sim\omega_{pe}$, where
$v_{b}\sim\sqrt{2e\phi/m_{e}}\sim 10c_{A0}$. The bounce rate is more than ten
times larger than the growth rate. Thus the amplitude of the electric field
$E_{z}$ evolves slowly compared with $\omega_{b}^{-1}$ during this interval.
By assuming the slowly evolving and slowly propagating (i.e. $z-v_{p}t\approx
z$) wave potential is $\phi(z,t)$, the electron Hamiltonian can be
approximated as:

$H\approx\frac{m_{e}}{2}v_{ez}^{2}-e\phi(z,t).$ (2)

Eq. (2) shows that during $\Omega_{e}t=38-64$, it is possible to choose $z$
and $v_{ez}$ as two approximate Hamiltonian canonical coordinates so that the
area $S=\frac{1}{2\pi}\oint v_{ez}dz$ enclosed by the electron’s trajectory in
phase space $(z,v_{ez})$ is an adiabatic invariant for trapped electrons,
where $v_{ez}=\sqrt{2(W(t,z)+e\phi)/m_{e}}$ and $W(t,z)$ is the electron’s
total energy. With the slow variation of the electric field, the electron’s
trajectory in phase space $(z,v_{ez})$ becomes narrower in $z$ and longer in
$v_{ez}$ as $E_{z}$ increases and becomes wider in $z$ and shorter in $v_{ez}$
as $E_{z}$ decreases. The electrons are trapped when the electric field grows
and are de-trapped when the electric field decays. The trapping and de-
trapping are non-adiabatic due to the change of the phase area inside and
outside of the wave potential Cary, Escande, and Tennyson (1986). The final
electron velocity depends on whether it crosses the upper or lower separatrix
as it is de-trapped. The upper (lower) separatrix crossing results in a
positive (negative) velocity in the wave frame.

To investigate how the adiabatic process converts kinetic energy into thermal
energy through non-adiabatic separatrix crossings of the wave potential, we
perform two test particle simulations with 5000 electrons in one single
standing wave $E=E_{z0}e^{\gamma t}\sin kz$. We take $kd_{i}\sim 80\sim
k_{0}$; $E_{z0}$ and $|\gamma|$ are constant and small. They satisfy
$\omega_{b}/|\gamma|\sim 80$ at the peak value of $E_{z}$. $\gamma>0$ for the
first half of the total duration and $\gamma<0$ for the second half so that
$E$ grows and decays sufficiently slowly to guarantee that the motion of the
trapped particles is adiabatic during the entire duration. The duration is
determined by the peak value of $E_{max}$. We investigate the cases with
$E_{max}\sim 40$ and $E_{max}\sim 80$. The initial electron velocity
distribution is a Maxwellian with a drift $v_{de}\sim 9c_{A0}$ and
$T_{e}=0.04m_{i}c_{A0}^{2}$ and the electron density is uniform in space. The
value $E_{max}=40$ is similar to the peak value of $E_{z}$ observed when the
PIC simulations can trap electrons with velocity below $v_{de}$. $E_{max}=80$
is higher than the peak value of $E_{z}$ observed in the PIC simulations. The
test single wave with $E_{max}=80$ can trap almost all of the electrons. The
results are shown in Fig.3.

More and more electrons are trapped as the electric field slowly increases and
the most are trapped at the peak value of $E$. The slight energy difference
between two trapped electrons leads to a large separation in their phase angle
around their trapped orbit since their angular velocity depends on energy.
Thus, at the time of the maximum trapping the trapped electrons are nearly
uniformly distributed along their trajectories as shown in panel a and the
velocity distribution of trapped electrons become flat as shown in panel c. As
$E_{z}$ decreases, the electron energy gain during trapping reverses and the
electrons are eventually de-trapped with the same value of $W$. The total
energy $W$ is symmetric with respect to positive and negative velocity.
Therefore, at the end of the simulation, due to the same probability for de-
trapping at the positive and negative velocity (Fig.3 panel b), a dip appears
near zero velocity in the velocity distribution function shown in panel d. The
red lines in panel c and d are for $E_{max}\sim 40$ where $E_{max}$ is not
strong enough to trap all of the electrons. As a result, the distribution
functions are not completely symmetric around $v_{ez}=0$.

Figure 4: The contours of total electron energy in phase space $(z,v_{ez})$ at
$\Omega_{e}t=$ 40, 102 from the PIC simulations.

The red line in panel c is similar to the electron velocity distribution
function at the saturation stage of the Buneman instability displayed in
Fig.1b.

However, even at $\Omega_{e}t=102$ in Fig. 1b, the electron velocity
distribution function keeps a similar shape to that at the saturation stage
rather than show a dip near zero as seen in the test particle simulation.
There are two reasons for the missing dip in the PIC simulations. First, the
wave amplitude is not spatially uniform. When the Buneman instability enters
the nonlinear stage, strong wave-wave interactions cause the collapse of the
uniformly distributed waves in space, and form localized solitary waves. As a
result, the trapping/de-trapping is more complex than in the sample model. In
Fig. 4, we plot the constant energy contours of electrons in phase space at
$\Omega_{e}t=$ 40 and 102 which correspond, respectively, to the peak of the
Buneman instability and late time of turbulence. We can see that the islands
between $z\in[0,0.2]$ and $z\in[0.3,0.4]$ are longer than those between
$z\in[0.2,0.3]$ and $z\in[0.4,0.5]$ at $\Omega_{i}t=40$ due to the
corresponding variation of the electric field $E_{z}$ in $z$ as shown in Fig.
1c. The long islands in phase space correspond to weak electric field and weak
trapping. Second, after $\Omega_{e}t=51$, the electron two-stream and Buneman
instabilities remain unstable, albeit weaker, and trapping continues. At the
late stage of the turbulence $\Omega_{e}t=102$, the islands in phase space
(Fig. 4) still exist even though the islands become longer.

In the PIC simulation, the electron trapping velocity covers the range
$[-10,10]c_{A0}\sim[-v_{de},v_{de}]$. The heating stops when the source of
kinetic energy is completely drained, i.e. the electron distribution with
velocity below $v_{de}$ becomes flat, as shown in Fig. 1 b. Nearly half of the
kinetic energy is dissipated, and the final electron temperature can be
estimated by $\triangle T_{ezz}\sim\triangle W\sim m_{e}v_{de}^{2}/2\sim 0.4$
which is consistent with what is observed in the simulation (Fig. 1a).

In reconnection, current sheets shrink as reconnection evolves and the Buneman
instability might occur in a wider current sheet with a reduced drift but a
similar growth rate $\sim(m_{e}/2m_{i})^{1/3}\omega_{pe}$. The implications of
these results in reconnection are being explored.

This research was supported by the NASA Postdoctoral Program at NASA/GSFC
administered by Oak Ridge Associated Universities through a contract with
NASA. The simulations and analysis were partially carried out at NASA/Ames
High-End Computing Capacity, at the National Energy Research Scientific
Computing Center, and at the National Institute for Computation Sciences.

## References

  * Kadomtsev (1965) B. B. Kadomtsev,  _Plasma turbulence_ (New York: Academic Press, 1965, 1965)
  * O’Neil (1965) T. O’Neil, “Collisionless Damping of Nonlinear Plasma Oscillations,” Phys. Fluid 8, 2255–2262 (1965)
  * Dupree (1966) T. H. Dupree, “A Perturbation Theory for Strong Plasma Turbulence,” Phys. Fluid 9, 1773–1782 (1966)
  * Sagdeev and Galeev (1969) R. Z. Sagdeev and A. A. Galeev,  _Nonlinear Plasma Theory_ (Nonlinear Plasma Theory, New York: Benjamin, 1969)
  * Galeev  _et al._ (1975) A. A. Galeev, R. Z. Sagdeev, V. D. Shapiro, and V. I. Shevchenko, “Nonlinear effects in an inhomogeneous plasma,” Akademiia Nauk SSSR Otdelenie Obshchei Fiziki i Astronomii Nauchnaia Sessiia Moscow USSR Uspekhi Fizicheskikh Nauk 116, 546–548 (1975)
  * Goldman (1984) M. V. Goldman, “Strong turbulence of plasma waves,” Rev. Mod. Phys. 56, 709–735 (1984)
  * Krommes (2002) J. A. Krommes, “Fundamental statistical descriptions of plasma turbulence in magnetic fields,” Physics Reports 360, 1–4 (2002)
  * Yampolsky and Fisch (2009) N. A. Yampolsky and N. J. Fisch, “Simplified model of nonlinear Landau damping,” Phys. Plasma 16, 072104 (2009)
  * Bénisti, Morice, and Gremillet (2012) D. Bénisti, O. Morice, and L. Gremillet, “The various manifestations of collisionless dissipation in wave propagation,” Phys. Plasma 19, 063110 (2012), arXiv:1111.1391 [physics.plasm-ph]
  * Buneman (1958) O. Buneman, “Instability, Turbulence, and Conductivity in Current-Carrying Plasma,” Phys. Rev. Lett. 1, 8–9 (1958)
  * Davidson  _et al._ (1970) R. C. Davidson, N. A. Krall, K. Papadopoulos, and R. Shanny, “Electron Heating by Electron-Ion Beam Instabilities,” Phys. Rev. Lett. 24, 579–582 (1970)
  * Ishihara, Hirose, and Langdon (1981) O. Ishihara, A. Hirose, and A. B. Langdon, “Nonlinear evolution of Buneman instability,” Phys. Fluid 24, 452–464 (1981)
  * Hirose, Ishihara, and Langdon (1982) A. Hirose, O. Ishihara, and A. B. Langdon, “Nonlinear evolution of Buneman instability. II - Ion dynamics,” Phys. Fluid 25, 610–616 (1982)
  * Cargill and Papadopoulos (1988) P. J. Cargill and K. Papadopoulos, “A mechanism for strong shock electron heating in supernova remnants,” The Astrophsical Journal Letters 329, L29–L32 (1988)
  * Drake  _et al._ (2003) J. F. Drake, M. Swisdak, C. Cattell, M. A. Shay, B. N. Rogers, and A. Zeiler, “Formation of Electron Holes and Particle Energization During Magnetic Reconnection,” Science 299, 873–877 (2003)
  * Che  _et al._ (2010) H. Che, J. F. Drake, M. Swisdak, and P. H. Yoon, “Electron holes and heating in the reconnection dissipation region,” Geophys. Res. Lett. 37, 11105–+ (2010), arXiv:1001.3203
  * Khotyaintsev  _et al._ (2010) Y. V. Khotyaintsev, A. Vaivads, M. André, M. Fujimoto, A. Retinò, and C. J. Owen, “Observations of Slow Electron Holes at a Magnetic Reconnection Site,” Physical Review Letters 105, 165002 (2010)
  * Riquelme and Spitkovsky (2009) M. A. Riquelme and A. Spitkovsky, “Nonlinear Study of Bell’s Cosmic Ray Current-Driven Instability,” The Astrophsical Journal 694, 626–642 (2009), arXiv:0810.4565
  * Matsumoto, Amano, and Hoshino (2012) Y. Matsumoto, T. Amano, and M. Hoshino, “Electron Accelerations at High Mach Number Shocks: Two-dimensional Particle-in-cell Simulations in Various Parameter Regimes,” The Astrophsical Journal 755, 109 (2012), arXiv:1204.6312 [astro-ph.HE]
  * Alexandrova  _et al._ (2009) O. Alexandrova, J. Saur, C. Lacombe, A. Mangeney, J. Mitchell, S. J. Schwartz, and P. Robert, “Universality of Solar-Wind Turbulent Spectrum from MHD to Electron Scales,” Phys. Rev. Lett. 103, 165003 (2009), arXiv:0906.3236 [physics.plasm-ph]
  * Sahraoui  _et al._ (2010) F. Sahraoui, M. L. Goldstein, G. Belmont, P. Canu, and L. Rezeau, “Three Dimensional Anisotropic k Spectra of Turbulence at Subproton Scales in the Solar Wind,” Phys. Rev. Lett. 105, 131101 (2010)
  * Kulsrud  _et al._ (2005) R. Kulsrud, H. Ji, W. Fox, and M. Yamada, “An electromagnetic drift instability in the magnetic reconnection experiment and its importance to magnetic ,” Plasma Physics 12, 082301 (2005)
  * Che, Drake, and Swisdak (2011) H. Che, J. F. Drake, and M. Swisdak, “A current filamentation mechanism for breaking magnetic field lines during reconnection,” Nature 474, 184 (2011)
  * Che (2009) H. Che,  _Non-linear development of streaming instabilities in magnetic reconnection with a strong guide field_ , Ph.D. thesis, University of Maryland, College Park (2009)
  * Che  _et al._ (2009) H. Che, J. F. Drake, M. Swisdak, and P. H. Yoon, “Nonlinear Development of Streaming Instabilities in Strongly Magnetized Plasma,” Phys. Rev. Lett. 102, 145004–+ (2009), arXiv:0903.1311 [physics.space-ph]
  * Cary, Escande, and Tennyson (1986) J. R. Cary, D. F. Escande, and J. L. Tennyson, “Adiabatic-invariant change due to separatrix crossing,” Physical Review A 34, 4256–4275 (1986)