text:
\section{Introduction}

The composition of the electric gird is in state of flux~\citep{Milligan:2015ju}.
Motivated by the need of reducing carbon emissions, conventional synchronous combustion generators, with relatively large inertia, are being replaced with renewable energy sources with little (wind) or no inertia (solar) at all~\citep{Winter:2015dy}.
Alongside, the steady increase of power electronics in the demand side is gradually diminishing the load sensitivity to frequency variations~\citep{WoodWollenberg1996}.
As a result, rapid frequency fluctuations are becoming a major source of concern for several grid operators~\citep{Boemer:2010wa,Kirby:2005uy}. Besides increasing the risk of frequency instabilities, this dynamic degradation also places limits on the total amount of renewable generation that can sustained by the grid. Ireland, for instance, is already resourcing to wind curtailment --whenever wind becomes larger than $50\%$ of existing demand-- in order to preserve the grid stability.

One solution that has been proposed to mitigate this degradation is to use inverter-based generation to mimic synchronous generator behavior, i.e. implement virtual inertia~\citep{Driesen:ft,Driesen:fta}. However, while virtual inertia can indeed mitigate this degradation, it is unclear whether that particular choice of control is the most suitable for this task.
On the one hand, unlike generator dynamics that set the grid frequency, virtual inertia controllers estimate the grid frequency and its derivative using noisy and delayed measurements.
On the other hand, inverter-based control can be significantly faster than conventional generators. Thus using inverters to mimic generators behavior does not take advantage of their full potential.

Recently, a novel dynamic droop control (iDroop)~\citep{m2016cdc} has been proposed as an alternative to virtual inertia that seeks to exploit the added flexibility present in inverters. Unlike virtual inertia that is sensitive to noisy measurements (it has unbounded $\mathcal H_2$ norm~\citep{m2016cdc}), iDroop can improve the dynamic performance without unbounded noise amplification.
However, as more sophisticated controllers such as iDroop are deployed,  the dynamics of the power grid become more complex and uncertain, which makes the application of direct stability methods harder.
The challenge is therefore to design an inverter control architecture that takes advantage of the added flexibility, while providing stability guarantees. Such an architecture must take into account the effect of delays and measurement noise in the design. It must be robust to unexpected changes in the network topology. And must provide a ``plug-and-play'' functionality by yielding decentralized --yet not conservative-- stability certificates.

In this paper we leverage classical stability tools for the Lur'e problem~\citep{BW65} to develop a novel analysis framework for power systems that allows a decentralized design of inverter controllers that are robust to network changes, delay uncertainties, and heterogeneous network components.
More precisely, by modeling the power system dynamics as the feedback interconnection of input-out bus dynamics and network dynamics (Section \ref{sec:network-model}), we derive a decentralized stability condition that depends only on the individual bus dynamics and the aggregate susceptance of the transmission lines connected to it (Section \ref{sec:stability-criterion}). When no network information is available, our conditions reduces to the standard decentralized passivity sufficient condition for stability.
We illustrate the novelty and strength of our analysis framework by studying the design of inverter-based control laws in the presence of delays (Section \ref{sec:analysis-idroop}).

\section{Input-output Power Network Model}\label{sec:network-model}

In this section we describe the input-output representation of the power grid used to derive our decentralized stability results. We use $i\in V:=\{1\dots,n\}$ to denote the $i$th bus in the network and the unordered pair $\{i,j\}\in E$ to denote each transmission line.

As mentioned in the previous section, we model the power network as the feedback interconnection of two systems, $P:=\diag\funof{p_i,i\in{}V}$ and $N$ shown in Figure \ref{fig:1}. Each subsystem $p_i$ denotes the $i$th \emph{bus dynamics}, and
maps $u_{P,i}\mapsto y_{P,i}$, where $u_{P,i}$ denotes the $i$th bus exogenous real power injection, i.e., the real power incoming to bus $i$ from other parts of the network or due to unmodeled bus elements. The output $y_{P,i}:=\omega_i$ denotes the bus frequency deviation from steady state.
Similarly, the system $N$ denotes the \emph{network dynamics}, which maps the vector of system frequencies $u_N:=\omega = (\omega_i,\,i\in V)$ to the vector of electric power network demand $y_{N}= (y_{N,i},\,i\in V)$, i.e. $y_{N,i}$ is the total electric power at bus $i$ that is being drained by the network. Thus, if we let $d_P = (d_{P,i},\,,i\in V)$ denote the unmodeled bus power injection, then by definition we get $u_{P}=-y_N+d_P$.
\begin{figure}
    \centering
    \includegraphics[height=.45\columnwidth,width=.65\columnwidth]{Figs/network-model}\vspace{-2ex}
    \caption{Input-output Power Network Model}
    \label{fig:1}
\end{figure}

This input-output decomposition provides a general modeling framework for power system dynamics that encompasses several existing models as special cases. For example, it can include the standard swing equations~\citep{Shen:1954eo}, as well as several different levels of details in generator dynamics, including turbine dynamics, and governor dynamics, see e.g., \citep{zmbl2016pscc}.
The main implicit assumption in Figure \ref{fig:1}, which is standard in the literature and well justified in transmission networks~\citep{kundur_power_1994}, is that voltage magnitudes and reactive power flows are constant. However, a similar decomposition could be envisioned in the presence of voltage dynamics and reactive power flows. Such extension will be subject of future research.

In this paper we focus on power system models that satisfy the following assumptions:
\begin{assumption}\label{ass:1}
$P\!=\!\diag\funof{p_1,\ldots{},p_n}$, with $p_i\in\H^{1\times{}1}$.
\end{assumption}
\begin{assumption}\label{ass:2}
$N=\frac{1}{s}L_B$, where $L_B\in\ensuremath{\mathbb{R}}^{n\times n}$ is a weighted Laplacian matrix.
\end{assumption}
Here $\H^{m\times n}$ denotes the Hardy space of $m$ by $n$ complex matrix transfer functions analytic on the open right-half  plane ($\mathrm{Re}\cfunof{s} >0$) and bounded on the imaginary axis ($j\ensuremath{\mathbb{R}}$). We only consider such a general class of transfer functions to model delays. In almost all remaining cases we only use real rational transfer functions.

In the rest of this section we illustrate how different network components can be modeled using this framework as well as the implications of assumptions \ref{ass:1} and \ref{ass:2}. For concreteness, we make specific choices on  the models for generators and loads. We highlight however that our analysis framework can be extended to more complex models provided they satisfy assumptions \ref{ass:1} and \ref{ass:2}.

\subsection{Bus Dynamics}
We model the bus dynamics using the standard swing equations~\citep{Shen:1954eo}. Thus the frequency of each bus $i$ evolves according to
\begin{equation}\label{eq:bus-swing}
p_i~:\begin{array}{rl}
\quad M_i\dot\omega_i &= -D_i\omega_i + x_i + u_{P,i}\\
y_{P,i}&=\omega_i
\end{array},
\end{equation}
where the state $\omega_i$ represents the frequency deviation from nominal and $x_i$ denotes the power injected by inverter-based generation at bus $i$.
The parameter $M_i\geq0$ denotes the aggregate bus inertia.
For a generator bus $M_i>0$, and $D_i>0$ represents the damping coefficient. For a load bus $M_i=0$ and $D_i>0$ represents the load sensitivity to frequency variations~\citep{BergenHill-1981}.

We consider linear control laws for the inverter dynamics that  depend solely on the local frequency. Thus, we model the inverter dynamics as a negative feedback law of the form
\begin{equation}\label{eq:control}
\hat{x}_i\ensuremath{\left(s\right)}=-c_i\ensuremath{\left(s\right)}\hat{\omega}_i\ensuremath{\left(s\right)},
\end{equation}
where $\hat x_i\ensuremath{\left(s\right)}$ and $\hat \omega_i\ensuremath{\left(s\right)}$ denote the Laplace transform of $x_i(t)$ and $\omega_i(t)$, respectively.

Combining \eqref{eq:bus-swing} and \eqref{eq:control} we get the following input-output representation of the bus dynamics
\begin{equation}\label{eq:bus-dynamics}
p_i\ensuremath{\left(s\right)}=\frac{1}{M_is+D_i}\funof{1+\frac{c_i\ensuremath{\left(s\right)}}{M_is+D_i}}^{-1}.
\end{equation}
Whenever $c_i(s)=0$,  \eqref{eq:bus-dynamics} represents a bus without inverter control, and by choosing either $M_i>0$ or $M_i=0$, \eqref{eq:bus-dynamics} can model generator or load buses respectively. Thus, \eqref{eq:bus-dynamics} provides a compact and flexible representation of the different bus elements.
It is important to notice that Assumption \ref{ass:1} is rather mild and only requires that the transfer function \eqref{eq:bus-dynamics} of each bus is stable.

Finally,  the control law $c_i\ensuremath{\left(s\right)}$ defines a general modeling class for inverter-based decentralized controllers. A number of conventional architectures can be written in this form. For example, \emph{virtual inertia} inverters correspond to
\begin{equation}\label{eq:droop-control}
c_i\ensuremath{\left(s\right)}=K_i+K_i^\nu s
\end{equation}
where $K_i\geq0$ and  $K_i^\nu\geq0$ are the  droop and virtual inertia constants, and by setting $K_i^\nu=0$, we recover the standard \emph{droop control}.
Similarly, the \emph{iDroop} dynamic controller is given by
\begin{equation}\label{eq:idroop}
c_i\ensuremath{\left(s\right)}= \frac{K^\nu_is+ K_i^\delta K_i}{s+K_i^\delta}
\end{equation}
where $K_i^\nu\geq0$, $K_i^\delta\geq0$ are tunable parameters\citep{m2016cdc}.

\subsection{Network Dynamics}\label{ssec:network-dynamics}

The network dynamics, under Assumption \ref{ass:2},  are given by
\begin{equation}\label{eq:net-dynamics}
N:~\begin{array}{rl}
    \dot \theta &= u_N\\
    y_N &= L_B\theta
\end{array},
\end{equation}
where, as mentioned before, the matrix $L_B\in\ensuremath{\mathbb{R}}^{n\times n}$ is the $B_{ij}$-weighted Laplacian matrix that describes how the transmission network couples the dynamics of different buses.

Thus Assumption \ref{ass:2} is equivalent to the DC power flow approximation, where the parameter $B_{ij}$ usually denotes the susceptance of line $\{i,j\}$, although more generally represents the sensitivity of the power flowing through line $\{i,j\}$ due to changes on the phase difference, i.e., $B_{ij}=\frac{v_iv_j}{x_{ij}}\cos(\theta_i^*-\theta_j^*)$ where $v_i$ and $v_j$ are the (constant) voltage magnitudes, $x_{ij}$ is the line inductance, and $\theta_i^*-\theta_j^*$ is the steady state phase difference between buses $i$ and $j$, see e.g., \citep{zhao2013power}.

To simplify the exposition, we refer here to $B_{ij}$ as the transmission line susceptance. Therefore, $$[L_B]_{ii}=\sum_{j:\{i,j\}\in E}B_{ij}$$ denotes the aggregate susceptance of the lines connected to bus $i$. As we will soon see in the next section, an upper bound of this value, i.e. $[L_B]_{ii}\leq\frac{1}{\gamma_i}$, can be used to guarantee stability in a decentralized manner.

\subsection{Connection with the Swing Equations}
For sake of clarity we derive here the full state space model of the described models and show how the standard swing equation model fits into our framework.
Since $u_N=y_P=\omega$ and $u_P = d_P -y_N$, then using \eqref{eq:bus-dynamics} and \eqref{eq:net-dynamics}, the state space representation of the feedback interconnection in Figure \ref{fig:1} amounts to:
\begin{subequations}\label{eq:swing}
\begin{align}
\dot{ \theta_i} &= \omega_i\label{eq:swing-a}\\
M_i\dot{ \omega_i} & = - D_i\omega_i -\!\!\!\sum_{j:\{i,j\}\in E}\!\!\!B_{ij}(\theta_i-\theta_j)  +x_i+d_{P,i}\label{eq:swing-b}.
\end{align}
\end{subequations}
where the inverter-based power injection $x_i$ evolves according to either
\begin{align}
\quad x_i &= -K_i\omega_i - K_i^\nu\dot\omega_i\\
\intertext{for virtual inertia inverter-based control, or}
\quad \dot x_i &= -K_i^\delta(K_i\omega_i + x_i) -K_i^\nu\dot\omega_i
\end{align}
of the case of iDroop.

Equation \eqref{eq:swing} amounts to the standard swing equations in which $d_{P,i}$ represents the net constant power injection at bus $i$. One interesting observation, and perhaps also a peculiarity of our framework, is that while usually the phase of bus $i$ ($\theta_i$) is considered to be part of the bus dynamics, in our framework this state is part of the network. This sidesteps the need to define, for example, a reference bus to overcome the lack of uniqueness in these angles.

\section{A Scalable Stability Criterion for Power Systems}\label{sec:stability-criterion}

This section consists of two main parts. First we give a generalisation of a classical stability result for single-input-single-output systems that can be applied to the structured feedback interconnection in Figure~\ref{fig:1}. Second, we discuss how it can be given a plug and play interpretation, and used to guide controller design in the context of the electrical power system models from Section~\ref{sec:network-model}.

\subsection{A Generalisation of a Classical Stability Criterion}

\cite{BW65} gives a method for testing stability of a feedback interconnection in the face of an uncertain gain. More specifically it is shown that for a real rational, stable, strictly proper transfer function $p\ensuremath{\left(s\right)}$, the feedback interconnection of $p\ensuremath{\left(s\right)}$ and $k\in\ensuremath{\mathbb{R}}$ is stable for all
\[
0\leq{}k\leq{}\frac{1}{k^*}
\]
if and only if there exists a \gls{spr} $h\ensuremath{\left(s\right)}$ such that
\begin{equation}\label{eq:absstb}
h\ensuremath{\left(s\right)}\funof{k^*+p\ensuremath{\left(s\right)}}
\end{equation}
is \gls{spr} (see \cite{DA96} for this precise statement). The concept of an \gls{spr} transfer function is given by the following standard definition (e.g. \cite{BL+06}):
\begin{definition}
A  (not necessarily proper) transfer functions $g\ensuremath{\left(s\right)}$ is said to be \gls{pr} if:
\begin{enumerate}[(i)]
\item $g\ensuremath{\left(s\right)}$ is analytic in $\mathrm{Re}\cfunof{s}>0$;
\item $g\ensuremath{\left(s\right)}$ is real for all positive real $s$;
\item $\mathrm{Re}\cfunof{g\ensuremath{\left(s\right)}}\geq{}0$ for all $\mathrm{Re}\cfunof{s}>0$.
\end{enumerate}
If in addition there exists an $\epsilon>0$ such that $g\funof{s-\epsilon}$ is \gls{pr}, then $g\ensuremath{\left(s\right)}$ is said to be \gls{spr}.
\end{definition}

The following theorem extends this result to the feedback interconnection in Figure~\ref{fig:1}, where the Laplacian matrix $L_B$ plays the role of the uncertain gain. Our motivation for trying to extend this result is driven by a desire for decentralised stability conditions. Since `the gain' of a Laplacian matrix is dependent on the network topology, in order to obtain a result that does not require exact knowledge of the network structure \textit{it is necessary} to be able to handle a degree of uncertainty in this gain. The result in \cite{BW65} is then the natural candidate for extension, because it is the strongest possible result of this type in the scalar case.

\begin{theorem}\label{thm:abs}
Let $\gamma_1,\ldots{},\gamma_n,$ be positive constants, and $P\in\ensuremath{\mathbf{H}_{\infty}}^{n\times{}n}$ satisfy Assumption~\ref{ass:1}. If there exists an $h\in\ensuremath{\mathbf{H}_{\infty}}$ such that $sh$ is \gls{pr}, and for each $i\in\cfunof{1,\ldots{},n}$:
\begin{equation}\label{eq:maintest}
h\ensuremath{\left(s\right)}\funof{\frac{\gamma_i}{2}s+p_i\ensuremath{\left(s\right)}}
\end{equation}
is \gls{spr}, then for any
\[
N\in\cfunof{\frac{1}{s}L_B:L_B \text{meets Assumption~\ref{ass:2} and $\sqfunof{L_B}_{ii}\leq{}\frac{1}{\gamma_i}$}},
\]
the feedback interconnection of $P$ and $N$ has no poles in the closed right half plane.
\end{theorem}

\begin{proof}
Since $P\in\ensuremath{\mathbf{H}_{\infty}}^{n\times{}n}$, the interconnection of $P$ and $N$ is stable if and only if
\[
N\funof{I+PN}^{-1}\in\ensuremath{\mathbf{H}_{\infty}}^{n\times{}n}.
\]
We will now show that the conditions of the theorem are sufficient for the above. Let $D^{-1}=\text{diag}\funof{\frac{\gamma_1}{2},\ldots{},\frac{\gamma_n}{2}}$, and define
\[
A=\funof{D^{\frac{1}{2}}L_BD^{\frac{1}{2}}}.
\]
Hence $N=D^{-\frac{1}{2}}\funof{\frac{1}{s}A}D^{-\frac{1}{2}}$, and $A$ is a normalised weighted Laplacian matrix. Factorise $A$ as
\[
A=QXQ^*,
\]
where $X$ is positive definite with eigenvalues $\lambda(X)\leq 1$ (i.e., $I\geq{}X>0$), $Q\in\ensuremath{\mathbb{C}}^{n\times{}\funof{n-m}}$, $m>0$, $Q^*Q=I$. Hence $N\funof{I+PN}^{-1}$ equals
\[
\begin{aligned}
&D^{-\frac{1}{2}}Q^*XQD^{-\frac{1}{2}}\funof{sI+PD^{-\frac{1}{2}}Q^*XQD^{-\frac{1}{2}}}^{-1},\\
=&D^{-\frac{1}{2}}Q^*X\funof{sI+QD^{-1}PQ^*X}^{-1}QD^{-\frac{1}{2}}.
\end{aligned}
\]
Clearly then it is sufficient to show that
\begin{equation}\label{eq:intintint1}
\funof{sI+QD^{-1}PQ^*X}^{-1}\in\ensuremath{\mathbf{H}_{\infty}}^{\funof{n-m}\times{}\funof{n-m}}.
\end{equation}
The above can be immediately recognised as an eigenvalue condition:
\begin{equation}\label{eq:necsuf2}
-s\notin\lambda\funof{Q^*D^{-1}P\ensuremath{\left(s\right)}{}QX},\forall{}s\in\bar{\mathbb{C}}_+.
\end{equation}
By Theorem 1.7.6 of \cite{HJ91}, for any $s\in\ensuremath{\mathbb{C}}$:
\[
\lambda\funof{Q^*D^{-1}P\ensuremath{\left(s\right)}{}QX}\subset\text{Co}\funof{\sqfunof{D^{-1}}_{ii}p_i\ensuremath{\left(s\right)},i\in N}\times{}\sqfunof{\epsilon,1}.
\]
Therefore it is sufficient to show that
\begin{equation}\label{eq:convsets}
-s\notin\text{Co}\funof{\sqfunof{D^{-1}}_{ii}p_i\ensuremath{\left(s\right)},i\in n}\times{}(0,1],
\end{equation}
for all $s=\bar{\mathbb{C}}_+$ (observe that since each $p_i$ is bounded it is enough to check that this holds on a sufficiently large, but bounded, subset of $\bar{\mathbb{C}}_+$). This can be checked using the separating hyperplane theorem, applied pointwise in $s$. In particular, \eqref{eq:convsets} holds if and only if there exists a nonzero $h\in\ensuremath{\mathbb{C}}$, such that $\forall{}i\in\cfunof{i,\ldots{},n}$:
\begin{equation}\label{eq:prtest}
\text{Re}\cfunof
{h\funof{s+k\sqfunof{D^{-1}}_{ii}p_i\ensuremath{\left(s\right)}}}>0,\forall{}\;0<k\leq{}1.
\end{equation}
By a very minor adaptation of the argument in Theorem 2 of \cite{BW65}, we can use a \gls{pr} function to define this $h$ pointwise in $s$. From the conditions of the theorem, and the maximum modulus principle,
\begin{equation}\label{eq:finaltest}
\text{Re}\cfunof{h\ensuremath{\left(s\right)}\funof{s+\sqfunof{D^{-1}}_{ii}p_i\ensuremath{\left(s\right)}}}>0,\forall{s}\in\bar{\ensuremath{\mathbb{C}}}_+.
\end{equation}
Since $sh\ensuremath{\left(s\right)}$ is \gls{pr}, for all $k^*\geq{}0$, $$\text{Re}\cfunof{h\ensuremath{\left(s\right)}\funof{s+\sqfunof{D^{-1}}_{ii}p_i\ensuremath{\left(s\right)}}+k^*sh\ensuremath{\left(s\right)}}>0.$$ Dividing through by $\funof{1+k^*}$ shows that under these conditions
\[
\text{Re}\cfunof{h\ensuremath{\left(s\right)}\funof{s+\frac{1}{1+k^*}\sqfunof{D^{-1}}_{ii}p_i\ensuremath{\left(s\right)}}}>0,\forall{}s\in\bar{\mathbb{C}}_+.
\]
Therefore \eqref{eq:prtest} is satisfied on the entire right half plane for the required range of $k$ values. Consequently \eqref{eq:intintint1} is satisfied, and the interconnection is stable.\end{proof}

\subsection{A Plug and Play Interpretation}\label{subsec:plug-and-play}

In order to give Theorem~\ref{thm:abs} a plug and play interpretation, we have to introduce a little conservatism. In direct analogy with the approach pioneered in \cite{LV06}, we propose to make an a-priori choice of the function $h\ensuremath{\left(s\right)}$. Observe that once we have done this, Theorem~\ref{thm:abs} suggests the following entirely decentralised procedure for constructing the network:
\begin{enumerate}
    \item For each component, compute the smallest $\gamma_i$ such that \eqref{eq:maintest} is satisfied.
    \item Check that $\frac{1}{\gamma_i}\geq{}\sqfunof{L_B}_{ii}$. If so, the component can be connected to the network.
\end{enumerate}
The interpretation of the above is that it defines a \textit{network protocol}. Observe that the above process works completely independently of the size of the network, and is therefore highly scalable. In addition we have the following appealing properties:

\textit{Robustness:} Suppose we have an uncertain description of our model $p_i$, for example:
\[
p_i\ensuremath{\left(s\right)}\in\cfunof{p\ensuremath{\left(s\right)}:p\ensuremath{\left(s\right)}=p_{i0}\ensuremath{\left(s\right)}+\Delta{}\ensuremath{\left(s\right)},\norm{\Delta\ensuremath{\left(s\right)}}_\infty<1}.
\]
Then provided we compute the smallest $\gamma_i$ such that \eqref{eq:maintest} is satisfied for all $p_i$ in this set, we can guarantee stability in a way that is robust to this uncertainty. This has not compromised scalability, since this process remains entirely decentralised.

\textit{Controller Design:} Suppose $p_i$ has some controller parameters that can be specified, for example (as in \eqref{eq:control})
\[
p_i\ensuremath{\left(s\right)}=\frac{1}{M_is+D_i}\funof{1+\frac{c_i\ensuremath{\left(s\right)}}{M_is+D_i}}^{-1},
\]
where $c_i$ can be chosen. Then provided we can design $c_i$ so that $p_i$ satisfies the protocol, the component can be connected to the network. Observe that this is an \textit{entirely decentralised synthesis condition}, and there is no need to redesign the controller in response to buses being added and removed elsewhere in the network. This is an approach to synthesis that specifically addresses the need for scalability.

The above are just observations about the structure of stability tests. The fact that they are based on testing \gls{spr}ness also brings a number of advantages. More specifically, provided the transfer functions in \eqref{eq:maintest} are rational and proper, this condition can be efficiently tested using the \gls{kyp} lemma. Furthermore, the synthesis problem: design $c_i$ such that
\[
h\ensuremath{\left(s\right)}\funof{\frac{\gamma_i{}}{2}s+\frac{1}{M_is+D_i}\funof{1+\frac{c_i\ensuremath{\left(s\right)}}{M_is+D_i}}^{-1}}\in\text{\gls{spr}};
\]
is equivalent to an $H^\infty$ optimisation problem, with state space solution given in \cite{SKS94}. Finally if $h\ensuremath{\left(s\right)}$ has relative degree 1, then \eqref{eq:maintest} can be checked on the imaginary axis (for a rather complete list of the different characterisations of \gls{spr}ness, see \cite{Wen88}). That is \eqref{eq:maintest} is equivalent to
\begin{equation}\label{eq:freqdom}
\mathrm{Re}\cfunof{h\ensuremath{\left(j\omega\right)}\funof{\frac{\gamma_i}{2}j\omega{}+p_i\ensuremath{\left(j\omega\right)}}}>0,\forall\omega\in\ensuremath{\mathbb{R}}.
\end{equation}
This allows much classical intuition from frequency responses and Nyquist diagrams to brought to bear on the above synthesis and robustness problems.

While obviously having many appealing features, all of the above rests on one crucial decision: the a-priori choice of $h\ensuremath{\left(s\right)}$. We make a few observations:
\begin{enumerate}
    \item If an arbitrary choice is made, this approach can be very conservative.
    \item One need not be overly concerned with optimising the choice of $h$ too much, because the degrees of freedom opened up by the synthesis problem can typically be exploited to meet a wide range of different protocols.
    \item There exists an $h$ such that if
    \[
    p_i\in\gls{spr},
    \]
    then \eqref{eq:maintest} will be satisfied. That is standard passivity results for the interconnection in Figure~\ref{fig:1} are a special case of Theorem~\ref{thm:abs}.
\end{enumerate}

We recommend that $h\ensuremath{\left(s\right)}$ be chosen by defining a set of `expected models'
\[
\mathcal{P}=\cfunof{\bar{p}_1\ensuremath{\left(s\right)},\ldots{},\bar{p}_m\ensuremath{\left(s\right)}},
\]
and then designing $h\ensuremath{\left(s\right)}$ so that \eqref{eq:maintest} is satisfied for all the models in this set. Here $\bar{p}_k$ are transfer functions describing the dynamics of `typical' bus models, perhaps obtained by putting in average values of the model parameters. Provided $m$ is not too large, this problem should be resolvable. The idea is that an $h$ that is designed to work for these transfer functions is a sensible `a-priori' choice, since it guarantees that if the actual bus dynamics are of the `expected' type (or can be designed to be close to the expected type), they will satisfy the network protocol.

\begin{comment}

Before illustrating the above on the iDroop power system model, we will first connect the above theorem and approach to the most well known large scale network design method: the passivity based approach. In the context of the feedback interconnection in Figure~\ref{fig:standardfeedback}, loosely speaking, the passivity theorem states that if both $P$ and $C$ are positive real, then their interconnection is stable (though not in the sense used in this paper). Since under Assumption~\ref{ass:2} $N$ is positive real, this suggests that a decentralised criterion along the following lines would also be sufficient: for each $i$
\[
p_i\text{ is positive real}.
\]
This could also be used a network protocol, and is in fact entirely independent of the transmission network! The above has striking similarities with Theorem~\ref{thm:abs}. However one would also suspect that it is unduly conservative, because Assumptions~\ref{ass:1} and \ref{ass:2} endow the model with considerable structure beyond positive realness. It turns out that this is not the case. The following shows that not only does a passivity type theorem arise as a special case of Theorem~\ref{thm:abs}, but if one desires a network independent result it is not possible to do better. In order to make this more concrete, we need the following definition:
\begin{definition}
Define $\mathcal{P}\subset\ensuremath{\mathbf{H}_{\infty}}^{1\times{}1}$ to be the set of transfer functions that satisfy:
\begin{enumerate}
    \item $z\funof{0}>0$.
    \item $\text{Re}\cfunof{z\ensuremath{\left(j\omega\right)}}\geq{}0$, for all $\omega\in\funof{0,\infty}$.
    \item If for $\omega>0$ $\text{Re}\cfunof{z\ensuremath{\left(j\omega\right)}}=0$, then $\text{Im}\cfunof{z\ensuremath{\left(j\omega\right)}}\geq{}0$.
\end{enumerate}
\end{definition}
The above is closely related to, but not equivalent to, the notion of a strictly positive real (SPR) transfer function (it includes all the SPR transfer functions). The following theorem shows that this class can be used to define a network independent protocol. The most significant observation here is that this class \textbf{cannot be made larger} if one wishes to guarantee stability against any network meeting Assumption~\ref{ass:2}. The implication is then that without making assumptions beyond \ref{ass:1} and \ref{ass:2}, there is no strictly better network independent stability test.

\begin{theorem}[Passivity Theorem]\label{thm:pass}
If $\forall{i}\in\cfunof{1,\ldots{},n}$:
\[
p_i\in\mathcal{P},
\]
then for any
\[
N\in\cfunof{\frac{1}{s}L_B:L_B \text{satisfies Assumption~\ref{ass:2}}},
\]
the interconnection of $N$ and $P$ has no poles in the closed right half plane. Furthermore there exists no subset $\mathcal{P}'\subset\ensuremath{\mathbf{H}_{\infty}}$ that is strictly larger than $\mathcal{P}$ ($\mathcal{P}'\supset{}\mathcal{P}$) that may be used in place of $\mathcal{P}$ in the above.
\end{theorem}

\begin{proof}
The claims on the locations of the closed loop poles can be established by applying Theorem~\ref{thm:abs} with
\[
h\ensuremath{\left(s\right)}\equiv\frac{1}{s\epsilon_1+1},
\]
for sufficiently small $\epsilon_1>0$. Since this $h$ has relative degree one, we need only consider the frequency domain characterisation in \eqref{eq:freqdom}. The effect of this $h$ is to rotate the frequency responses of $\funof{p_i\ensuremath{\left(j\omega\right)}+j\omega}$ clockwise by up to $\frac{\pi}{2}$ (and stretch them by some arbitrary amount).  Now observe that since the $p_i's$ are bounded, there exists an $\omega_0>\norm{p_i}_\infty$. Therefore for all $\omega\geq\omega_0$ the frequency domain inequality will be satisfied, because the $h\ensuremath{\left(j\omega\right)}{}j\omega$ term will be sufficiently far enough into the right half plane to overcome the $h\ensuremath{\left(j\omega\right)}{}p_i\ensuremath{\left(j\omega\right)}$ term. We now turn our attention to the lower frequency range. Observe that since the $p_i$'s are bounded there exists an $\epsilon_2>0$ such that $\forall{i}\;\angle{p_i}\in\sqfunof{\frac{-\pi}{2}+\epsilon_2,\frac{\pi}{2}}$. Also note that by making $\epsilon_1$ sufficiently small, the clockwise rotation from $h\ensuremath{\left(j\omega\right)}$ can be made arbitrarily small over any frequency range $\funof{0,\omega^*}$. Therefore choosing $\epsilon_1$ small enough allows us to rotate any points in the frequency range $\funof{0,\omega_0}$ where the $p_i$'s touch the imaginary axis into the right half plane, without rotating any other point sufficiently to touch the negative imaginary axis. This means that $\text{Re}\cfunof{h\ensuremath{\left(j\omega\right)}{}p_i\ensuremath{\left(j\omega\right)}}>0$, except for points where $p\ensuremath{\left(j\omega\right)}=0$ where the inequality will not be strict. However here $\text{Re}\cfunof{h\ensuremath{\left(j\omega\right)}{}jw}>0$, hence \eqref{eq:freqdom} is guaranteed on this frequency range too. Finally, \eqref{eq:freqdom} holds at $\omega=0$ directly from the definition of $\mathcal{P}$.

We will now establish that there exists no $\mathcal{P}'\supset{}\mathcal{P}$ that may be used in place of $\mathcal{P}$. Assume the converse. We will now construct $p_i$'s and an $N$ from the allowable sets that contradict the claims of the theorem. Let $p_1\in\mathcal{P}'\setminus\mathcal{P}$, and let $\omega_0>0$ be a frequency such that $\angle{p_1\ensuremath{\left(j\omega\right)}}\notin\left(-\frac{\pi}{2},\frac{\pi}{2}\right]$.
Therefore $p_1\ensuremath{\left(j\omega_0\right)}$ must lie in the open left half plane or on the negative imaginary axis. Next we will construct a $p_2\in\mathcal{P}'$ such that
\[
p_1\ensuremath{\left(j\omega_0\right)}+p_2\ensuremath{\left(j\omega_0\right)}=-\alpha{}j\omega_0,
\]
for some $\alpha>0$. If $p_1\ensuremath{\left(j\omega_0\right)}$ is on the imaginary axis, pick $p_2=p_1$. Otherwise $-j\omega_0-p_2\ensuremath{\left(j\omega_0\right)}$ must lie in the open right half plane, and therefore all we must do is pick a function that interpolates this number at this frequency. This can always be done with a function in $\mathcal{P}$, for example, a function of the form
\[
k_1+\frac{k_2s}{1+k_3}\;\;\;\text{or}\;\;\;k_1+\frac{k_2s}{1+k_3s},
\]
where $k_1,k_2,k_3$ are strictly positive, depending on whether $-j\omega_0-p_2\ensuremath{\left(j\omega_0\right)}$ is in the first or fourth quadrant. Since $\mathcal{P}'\supset{}\mathcal{P}$, this function is also in $\mathcal{P}'$ as required. Now choose $L_B$ to be zero everywhere except with upper left corner given by
\[
\frac{1}{\alpha{}}\begin{bmatrix}
1&-1\\-1&1
\end{bmatrix},
\]
and $p_3,\ldots{},p_n\in\mathcal{P}'$. In order to satisfy the claim of Theorem~\ref{thm:pass} with $\mathcal{P}'$, it is necessary that the interconnection does not have an imaginary axis pole at $j\omega_0$ with the above choices of $P$ and $N$. However
\[
\begin{aligned}
\det\funof{I+P\ensuremath{\left(j\omega_0\right)}{}N\ensuremath{\left(j\omega_0\right)}}&=\det\funof{I+P\ensuremath{\left(j\omega_0\right)}{}L_B\frac{1}{\alpha{}j\omega_0}}\\
&=1+\frac{\funof{p_1\ensuremath{\left(j\omega_0\right)}{}+p_2\ensuremath{\left(j\omega_0\right)}{}}}{\alpha{}j\omega_0}=0.
\end{aligned}
\]
This implies that there is an imaginary axis pole at $j\omega_0$ with this $P$ and $N$. This is a contradiction, and hence there can exist no $\mathcal{P}'\supset\mathcal{P}$ that satisfies Theorem~\ref{thm:pass}.
\end{proof}

We have labelled this theorem the passivity theorem, because, to all intensive purposes, the set $\mathcal{P}$ is the set of strictly positive real transfer functions. Certainly, any strictly positive real transfer function also lies in $\mathcal{P}$. The only difference is that we allow some frequency responses that touch the positive imaginary axis, for example
\begin{equation}\label{eq:wierd}
\frac{\funof{s+1}^2}{\funof{s+5}^2}-\frac{1}{50},
\end{equation}
or that do not meet the usual requirements at $\omega=\infty$, such as
\[
\frac{s+1}{s^2+s+1}.
\]
These differences are more or less meaningless from a practical perspective. Transfer functions such as \eqref{eq:wierd} are allowed because we not only know that $C$ is passive, but also that it only contains inductors. It is doubtful whether we would want to use this type of transfer function in practice. We are, however, obliged to include them in our definition to show the tightness of our result, since if we excluded them there would exist a $\mathcal{P}'\supset\mathcal{P}$ on which Theorem~\ref{thm:pass} would hold.

\begin{remark}
There do exist subsets $\mathcal{P}'\subset\ensuremath{\mathbf{H}_{\infty}}$ that are not strictly smaller than $\mathcal{P}$ that may be used in place of $\mathcal{P}$ in Theorem~\ref{thm:pass}. In other words passivity is not the only decentralised condition on the $p_i$'s that is independent of $C$. One interesting example is the set of `positive imaginary' systems
\[
\mathcal{P}'=\cfunof{z\in\ensuremath{\mathbf{H}_{\infty}}^{1\times{}1}:z\ensuremath{\left(j\omega\right)}-z^*\ensuremath{\left(j\omega\right)}>0,\forall\omega\in\funof{-\infty,\infty}}.
\]
It is possible to find such classes because $C$ has structure beyond passivity. We will not develop this further here, but note that for our application the set above is not useful since the $p_i$'s arising will additionally be strictly proper, and hence cannot lie in $\mathcal{P}'$.
\end{remark}

\end{comment}

\section{Scalable Stability Analysis of iDroop}\label{sec:analysis-idroop}

To motivate both the need for Theorem~\ref{thm:abs}, and to illustrate it's application, consider a two bus network. Assume that the buses have the same dynamics, given by
\begin{equation}\label{eq:fb11}
p_i=\frac{1}{s+0.1}\funof{1+\frac{1}{s+0.1}c_i\ensuremath{\left(s\right)}}^{-1},
\end{equation}
and the Laplacian matrix describing the transmission network is given by
\[
L_B=\begin{bmatrix}
1&-1\\-1&1
\end{bmatrix}.
\]
We now consider the problem of how to design an iDroop controller subject to delay. That is, how to design
\[
c_i\ensuremath{\left(s\right)}=e^{-s\tau_i}\frac{K_i^\nu{}s+ K^\delta_i K}{s+K_i^\delta},
\]
where $\tau_i>0$ is the delay, such that the interconnection is stable. If the delay equals zero, then passivity theory can be easily applied to answer the stability question. This is because \eqref{eq:fb11} can be rearranged as
\[
p_i=\funof{\funof{\frac{1}{s+0.1}}^{-1}+\funof{c_i\ensuremath{\left(s\right)}^{-1}}^{-1}}^{-1},
\]
from which we see that $p_i$ is the parallel interconnection of $c_i^{-1}$ and an \gls{spr} transfer function. This means that if $c_i$ is \gls{spr}, so is $p_i$, which in turn implies that the interconnection with the transmission network is stable (the transmission network is a \gls{pr} transfer function, and it is well known that the feedback interconnection of a \gls{pr} and \gls{spr} transfer function is stable). Hence we arrive at the rather surprising conclusion that any possible choice of iDroop controller will lead to a stable interconnection. This becomes even stranger when realise we can use the same argument in networks of any size. Is it really true that any possible interconnection of heterogeneous buses, equipped with arbitrarily chosen iDroop controllers, interconnected through any possible network, is stable?

The answer is, at least in any practical sense, no. To make this absolutely concrete, let us return to the two bus example, and suppose we choose the following parameters for both our iDroop controllers:
\[
K_i^\nu=1,\;K_i^\delta{}=5,\;K_i=30.
\]
A simple Nyquist argument shows that even a small value of $\tau_i$, for example 0.05, will destabilise this two bus network. It is of course not a limitation of iDroop that this can happen; the introduction of delays (and badly designed controllers) can similarly destabilise networks where the buses employ virtual inertia or droop controllers. Nor is it even that surprising, and the authors would claim no originality for the above argument. It does however serve to illustrate the importance of robust design in the network setting. We will now demonstrate that, as claimed in Section~\ref{subsec:plug-and-play}, Theorem~\ref{thm:abs} can be used to do this.

\begin{figure}
    \centering
    \includegraphics[width=.85\columnwidth]{Figs/Bode}\vspace{-2ex}
    \caption{Bode diagram of $p_i$ (solid curve), and first order approximation of the form in \eqref{eq:approxp}, with $a=1.37,b=1$ (dashed curve).}
    \label{fig:2}
\end{figure}

As discussed in Section~\ref{subsec:plug-and-play}, to obtain a plug and play criterion, we need to make an a-priori choice of $h$. It was also claimed that we need not be overly careful when making this choice, so let us just pick
\[
h\ensuremath{\left(s\right)}=\frac{1}{\frac{s}{\omega_0}+1},
\]
and proceed with the robust synthesis of $c_i\ensuremath{\left(s\right)}$. Conducting a classical lead-lag design (this is another advantage of the iDroop controller structure) allows us to achieve a $p_i$ which, frequency by frequency, is similar to a transfer function of the form
\begin{equation}\label{eq:approxp}
\frac{a_i}{s+b_i}.
\end{equation}
For a delay of $\tau_i=0.5$ and the model paramters in \eqref{eq:fb11}, this is illustrated in Figure~\ref{fig:2}, where the iDroop controller specified by
\[
K^\nu_i=1.3,\;K_i^\delta{}=8,\;K_i=0.65,
\]
has been chosen. To conclude stability using Theorem~\ref{thm:abs}, we are required to verify that
\begin{equation}\label{eq:newapprox}
\frac{1}{\frac{s}{\omega_0}+1}\funof{\frac{\gamma_i}{2}s+\frac{a_i}{s+b_i}+\Delta_i\ensuremath{\left(s\right)}}\in\textrm{\gls{spr}},
\end{equation}
for some $\gamma_i\leq{}1$. In the above $\Delta_i\ensuremath{\left(s\right)}$ is the difference between \eqref{eq:approxp} and the actual $p_i$. A simple way to do this is to instead verify the condition
\begin{equation}\label{eq:newnewapprox}
\frac{1}{\frac{s}{\omega_0}+1}\funof{\frac{\gamma_i}{2}s+\frac{a_i}{s+b_i}}-\epsilon_i{}\in\textrm{\gls{pr}}.
\end{equation}
Then provided
\[
\norm{\frac{1}{\frac{s}{\omega_0}+1}\Delta_i\ensuremath{\left(s\right)}}_\infty{}<\epsilon_i{},\;\Longleftrightarrow{}\;\abs{\Delta_i\ensuremath{\left(j\omega\right)}}<\epsilon{}_i\sqrt{1+\frac{\omega^2}{\omega_0^2}},
\]
\eqref{eq:newapprox} is also guaranteed. This relaxation is shown in Figure~\ref{fig:3}.

\begin{figure}
    \centering
    \includegraphics[width=.85\columnwidth]{Figs/Nyquist}\vspace{-3ex}

    \caption{Nyquist diagram of $p_i$ (black curve), and the uncertainty set on $\Delta_i$ with $\epsilon_i=0.08$ (shaded region). Since the Nyquist diagram lies completely inside the shaded region, satisfying \eqref{eq:newnewapprox} with this $\epsilon_i$ is sufficient for \eqref{eq:newapprox}.}
    \label{fig:3}
\end{figure}

In addition to guaranteeing further levels of robustness to unmodelled dynamics, \eqref{eq:newnewapprox} gives a simple way to characterises any allowable iDroop design in terms of only three parameters: $\funof{a_i,b_i,\epsilon_i{}}$. This is an advantage because the \gls{pr} condition on \eqref{eq:newnewapprox} can be rewritten in terms of inequalities written in terms of these parameters. This means that our network protocol can be understood in terms of (relatively) simple inequalities that only depend on coarse features of our iDroop design. After some simplifications, using the result of e.g.  \cite{FL63}, it can be shown that \eqref{eq:newnewapprox} is equivalent to the following two inequalities:
\[
\begin{aligned}
a_i-\epsilon{}_ib_i&\geq{}0\\
b_i\funof{\gamma{}_i\omega_0-2\epsilon{}_i}-2\epsilon_i\omega_0&\geq{}\ldots{}\\
\ldots{}\frac{\omega_0}{\funof{b_i+\omega_0}}&\funof{\sqrt{a_i-\epsilon{}_ib_i}-\sqrt{b_i\funof{\frac{\gamma_i\omega_0}{2}-\epsilon{}_i}}}^2
\end{aligned}
\]
These give expressions that the design of each individual iDroop controller must satisfy in order to (robustly) satisfy the scalable stability tests. We can easily check them for the two bus example. Substituting in $\funof{a_i,b_i,\epsilon_i}=\funof{1.37,1,0.08}$, and $\omega_0=30$ shows that they are satisfied for any
\[
\gamma_i\geq{}0.18,
\]
which is clearly sufficient to conclude stability. Note that if they were failed, we would then just have to go back and retune the design of our iDroop controller to obtain more favourable parameters (or test \eqref{eq:newapprox} instead). To stress the scalability of this procedure, observe that we did not use the fact that this was a two bus example at any stage. This same local method can be used for designing iDroop controllers for heterogeneous buses in networks of any size, and the above shows that this specific bus can be connected into any possible transmission network provided the local network susceptance is greater than 0.18.

redpajama_set_name: RedPajamaArXiv


text:
\section{Introduction:}
\label{intro.sec}
	Ultra-high energy cosmic ray particles with energies greater than
about $5 \times 10^{19}$ eV cannot traverse great distances through the
cosmos without
losing energy to pion photoproduction in collisions with the cosmic
microwave background. Such a particle will lose energy until it is
below threshold for this reaction. This process, referred to as the
GZK cutoff (Greisen, 1966, and Zatsepin \& Kuzmin, 1966), happens in
about 50-100 Mpc (Cronin, 1992 and Aharonian \& Cronnin, 1994).
 Hence for a $10^{20}$ eV particle to arrive at Earth, it must
originate in the local neighborhood. The observed arrival direction of
  events with energies greater than $5 \times 10^{19}$ eV  provide no
compelling
evidence for a nearby source. However there
is a possible indication of correlations in event arrival direction (Cronin
1999).
Further, there are no known nearby astrophysical objects which
are likely to accelerate particles to these energies (Hillas,
1984 and Biermann, 1997). This paradox
is one of the major current mysteries in astrophysics.

\section{NNODE and WALTA:}
\label{NNODE.sec}
We are forming a collaboration for the development of a
distributed detector network that will support measurements of air
showers from ultra-high energy cosmic rays and can also support a broad class
of other physical measurements. Other researchers have expressed interest in
using the framework, for example,
 as a seismograph network for
geophysics studies, a passive radar network for ionosphere studies, and air
sampling
for atmospheric studies.  We call the overall network NNODE (Northwest Network
for Operation of Distributed Experiments), and we call the cosmic ray
measurement
component WALTA (WAshington Large-area Time-coincidence Array).  The focus
of the efforts of our subgroup within NNODE is WALTA.
The project is to be a direct physical science
outreach program between faculty and students of the University of Washington
and the science teachers and students of Washington-area schools (grades 7-12).
 The detection techniques and philosophy
 are similar to those of the ALTA project in Alberta, Canada
(Pinfold). The groups in Washington and Alberta
have a loose collaboration.

Technological advances of the past decade have made such distributed
experiments practical. Powerful yet inexpensive desktop computers are
available and each detection site can have its own unit. The Global
Positioning System (GPS) can provide very good absolute timing (Berns
99) and positioning between distributed sites. Finally the Internet
provides convenient communication between each of the sites.  Each
WALTA/NNODE measurement module is envisioned to consist of a computer
with an Internet connection, a GPS timing system, and measuring
equipment.  The measuring devices used will depend on the science
investigation.  For the WALTA investigation, modules will consist of
scintillation paddles to detect distributed particle showers produced
by ultra-high energy cosmic rays. We plan to ultimately locate the
majority of the modules at public and private schools and to enlarge
the network as additional schools participate.

An aspect of this program is that the teachers, and through them their
students, will have the opportunity to directly contribute in cutting
edge scientific research.  We expect that this direct involvement will
motivate and energize teachers and result in their learning science at
a deeper level.  It may also spark the imaginations of their students,
encouraging their consideration of careers in the sciences.  The
teachers and their students will be able to learn about scientific
techniques, mathematical tools, and the latest measurement technology
and see them applied as new results emerge.  Many elements of the
project will provide a basis for classroom teaching units in science,
mathematics, and technology.  The classroom learning will be made more
immediate by the ongoing measurements of the systems and the emergence
of new results.

\section{WALTA and Cosmic Ray Studies:}
\label{WALTA.sec}
An ultra-high energy cosmic ray colliding with the upper atmosphere will
produce
a cascade of particles which at ground level contains
several billion particles and covers tens of km$^2$.  Sampling
the particles from such a shower requires a distributed detector with sites
separated on the order of a km.

To study the potential performance of WALTA, we chose 32 candidate sites in the
Seattle area. (See Fig. 1) These locations are associated with the
University of Washington,
high schools, middle schools, and other local universities and community
colleges.
The list is neither exhaustive nor confirmed. At each site we anticipate
placing
three or four 1-inch thick scintillators of approximate area 1 m$^{2}$
each. These detectors
would be separated by approximately 10-20 m. At this spacing, the
individual sites will measure
cosmic rays of energy near the knee ($\approx 10^{15}$ eV). This will
provide students with high rate
data similar in character to the low rate data from ultra high energy events.

Each site would trigger independently.
We would require that some number of detectors
at each site observe a signal corresponding to a minimum number of vertical
equivalent muons. The data from
each site would be transferred via the Internet to a central computing
location
to be analyzed off-line
for inter-site coincidences.

The 32-site array covers a large area, approximately 200 km$^{2}$. However,
the spacing
is not uniform. The mean nearest-neighbor spacing between these sites is
approximately
1.4 km but varies greatly. As a result, ``holes" in the array reduce the
trigger
efficiency significantly. However, the autonomous nature of the individual
sites makes
the array very easy
to expand. That is, holes can be filled and the array can easily be
extended into
the urban areas north, east, and south of Seattle. Thus the array can be
made more
efficient and larger.

Many of the site and array parameters must be optimized. This optimization
is being studied
through simulation using the CORSIKA code, version 5.624 (Heck, 1998).
Figure 1 shows the site
layout with a simulated vertical 10$^{19}$ eV event to provide some
indication of the array response.
In the figure, the squares
are centered on sites which have no response to the event. The 5 sites
responding to the
event are indicated by circles. The area of
each circle is proportional to the number of particles detected at the
site. The number
of detectors
and their geometry, the detector spacing, and the trigger requirements are
all parameters
which are currently
under study. Furthermore, site locations required to eliminate holes in the
effective
detector sensitive area are being identified. Subsequently, we will search
for new candidate sites to fill such holes.

One very dangerous possibility is that
of high-energy tails on the resolution function. The observation of the
highest energy
cosmic rays is the critical feature of the experiment. If low energy events
below the
GZK cutoff can be falsely identified as having a high energy, it will
seriously compromise
the conclusions one can draw from the measurements. Simulations to study
this possibility
are also underway.

Finally efforts are underway to finalize the hardware design. To ensure the
ease
of expansion, the cost of each site must be minimized. Prototypes of a site
system are
being assembled and characterized.

\section{Conclusion:}
\label{conclusion.sec}

The scientific activity of the WALTA project has a number of aspects which can
feed into teacher and student research experience.  The technological examples
provided by the network of computer systems, the display of data as it is
acquired, the use of GPS circuits for timing, and the electronics of the
detection system will serve as challenging examples for the technologically
inclined students and the teachers of technology-based classes.
The analysis and interpretation of the results from many of the measurements
will provide case studies for the application of statistics and mathematics
that can be linked to mathematics classes.  The physics that underlies the
research activity also has many potential connections that can enrich high
school physics classes.

We wish to acknowledge the support and encouragement of Dean David Hodge,
College of Arts and Sciences, University of Washington. We also thank Steve
Ellis and David Boulware for their encouragement and advice. These efforts are
supported in part by the Department of Energy through the University of
Washington
Nuclear Physics laboratory.

\begin{figure}[bp]
\centering
\epsfig{figure=seattle.eps,width=5in}
\vspace{0.7in}
\caption[WALTA Layout]
{The Seattle area indicating a number of suggested sites.}
\label{walta-layout}
\end{figure}
\vspace{1ex}
\begin{center}
{\Large\bf References}
\end{center}
Aharonian, F.A. \& Cronin, J.W., 1994, Phys. Rev. {\bf D50}, 1892.\\
Berns, H. G. \& Wilkes, R. J., 1999, Proc. 11th IEEE NPSS Conference.
Biermann, P., 1997, J. Phys. G {\bf 23}, 1.\\
Cronin, J.W., 1992, Nucl. Phys. B (Proc. Supp.) {\bf 28B}, 213.\\
Cronin, J.W., 1999, Rev. Mod. Phys. {\bf 71}, 165.\\
Greisen, K., 1966, Phys. Rev. Lett. 16, 748.\\
Heck, D., {\it et al.}, 1998, Karlsruhe Report FZKA 6019;
   http://www-ik3.fzk.de/~heck/corsika/Welcome.html \\
Hillas, A.M., 1984, Astron. Astrophys. {\bf 22}, 425.\\
Pinfold, J.L., {\it et al.}, see http://csr.phys.ualberta.ca/~alta for
proposal text.\\
Zatsepin, G.T., \& Kuzmin, V.A., 1966, JETP Lett. {\bf 4}, 78.\\

\end{document}


text:
\section{Transient alerts in multimessenger Astronomy}
The field of astrophysics is experiencing several fundamental changes, such as the increasing relevance of observations of transients, i.e., temporal variabilities of astronomical sources such as supernova explosions, fast radio bursts (FRBs), and gamma-ray bursts (GRBs). An increasing number of complementary cosmic messengers provide crucial information about these events. The detections of a diffuse flux of high-energy neutrinos and of gravitational waves from mergers of compact binaries have started to supplement electromagntic observations. Combining these complementary messengers can provide crucial insights on events in the transient sky. The study of these extreme (flaring) events in the multimessenger approach requires real-time and, to a large extent, fully automated communication between telescopes (see e.g.\,\citet{2019ICRC743M,2021JCAP...03..045A}) and researchers to respond instantaneously to events with follow-up observations to better understand the nature of these events through complementary observational data. To tackle these challenges, among others, an ecosystem of numerous complimentary services has been established over the last decades. For a review of high-energy alerts in the multimessenger era, see \citet{Dorner2021}.
\newpage

\section{Realtime alerts on mobile devices}
Whereas the time-critical decisions for follow-up observations must be automatic, especially the evaluation and the decision on subsequent observations require extensive processing and a graphical representation, easily comprehensible for humans, of the information available on the event in the context of already known sources and temporal and spatial correlations of other transient signals.
\articlefigure{X3-001_f1.pdf}{fig1}{The multimessenger pipeline from observations of transient events to the user notification and its visualization in the Astro-COLIBRI app. Most processes are fully automatized and with negligible latency and downtime. Exceptions are e.g. the manual creation of GCN circulars that deliver information in human-written reports about the observed event. Details of the Astro-COLIBRI pipeline and the message transport of the push alerts to the mobile devices are shown in Fig.~\ref{fig2}.}

Astro-COLIBRI \citep{Reichherzer2021} aims to be the top layer that connects existing subsystems into a large ecosystem, with a focus on optimized display for users through an interactive graphical user interface (GUI) and mobile apps that receive push notifications in real time, all usable in an intuitive way. The multimessenger pipeline from observations to the user is discussed in the following and sceteched in Fig.~\ref{fig1}.
\begin{enumerate}
    \item \textit{Monitoring satellites and observatories:} The foundation of multimessenger astronomy are satellites and observatories that constantly monitor the various messengers of the transient sky. Gravitational waves (Ligo and Virgo), photons (Fermi, HAWC, Swift, INTEGRAL, SVOM, etc.), as well as neutrinos (IceCube, ANTARES, and GVD), can be localized very rapidly, which makes them essential for follow-up observations. Alerts announcing new detections are emitted largely automatically when specific criteria are met.
    \item \textit{Alert creators / streams / broker:} Alert and correlation systems such as VOEvent alerts, Astrophysical Multimessenger Observatory Network alerts \citep{AMON_2020}, Transient Name Server (TNS) notifications\footnote{\url{https://www.wis-tns.org/}}, Gamma-ray Coordinates Network (GCN) circulars\footnote{\url{https://gcn.gsfc.nasa.gov/}}, brokers (e.g. via 4 Pi Sky; \citet{4pisky_2016}, FINK broker; \citet{FINK_2021}) provide alerts which differ both in their degree of machine or human-generated content and in the underlying communication technology. The technologies involve the use of web pages, emails, and brokers and require the end-user to make an effort to keep track of the (ever-increasing amount of) relevant information.
    \item \textit{Astro-COLIBRI pipeline:} The Astro-COLIBRI pipeline processes machine or human-generated alerts in real-time and provides the event in the context of persistent sources as well as transient events in the relevant phase space. To accommodate the growing range of services from the community, Astro-COLIBRI is modular in design and allows easy integration of new or modifying third-party services as shown in Fig.~\ref{fig2}.
    \item \textit{Push notifications:} Astro-COLIBRI informs users based on their preferences set in their Astro-COLIBRI app about these events in real-time (latency of a few 100\,ms) via push notifications on smartphones. Once opened the app, the interface provides a graphical representation with a summary of the relevant data to allow for the fast identification of interesting phenomena along with an assessment of observing conditions at a large selection of observatories around the world. See \citet{Schuessler2021} for a description of some use cases.
\end{enumerate}
\articlefigure{X3-001_f2.pdf}{fig2}{Architecture of Astro-COLIBRI comprises a REST API, both a static and a real-time database, a cloud-based alert system, as well as a website and apps for iOS and Android as clients. Public user interfaces are highlighted in green.}

\section{Outlook}
To enable a better exchange of relevant information on specific transient events, a comment function on individual events is planned. The collected information from the contributed comments can later be automatically processed and used in the context of new alerts. The automated evaluation of human-generated alerts will be accomplished in the future by means of deep learning-based natural language processing techniques within Astro-COLIBRI, thus enabling rapid extraction and synthesis of information from reports that will be associated with the respective transient events.\\\\
The Astro-COLIBRI development team welcomes comments and feedback from the community to further improve the platform.
\newpage
\acknowledgements This work was supported by the European Union's Horizon 2020 Programme under the AHEAD2020 project (grant agreement No. 871158). This work is supported by the ADI 2019 project funded by the IDEX Paris-Saclay, ANR-11-IDEX-0003-02 (PR). PR also acknowledges support from the German Academic Exchange Service and by the RUB Research School via the \textit{Project International} funding. JBT and PR acknowledge funding by the BMBF, grant number 05A20PC1.

\textit{Used software / packages / frameworks:} 4 Pi Sky~\citep{4pisky_2016}, Astropy \citep{2018AJ....156..123A}, Flutter\footnote{\url{https://flutter.dev/}}


text:
\section{\textbf{INTRODUCTION}}
During the past two decades since the successful isolation of graphene from graphite \cite{nov04,nov05},
the exploration of two-dimensional (2D) layered materials
has received much attention and is likely to remain one of the leading topics in material science into foreseeable future.
Although there have been a number of experimental successes and theoretical predictions on discovering 2D materials, including graphene, \textit{hexagonal} boron nitride (\emph{h}-BN), and transition metal dichalcogenides (TMDs) \cite{nov04,nov05,nagashima95,mak10,geim13}, their potentials in spintronic applications have been limited by the lacks of spin orders and possibilities of magnetization flip.
Hence, the current research in this field primarily focuses on seeking a novel 2D material with intriguing magnetic properties \cite{leht04,han13,gon16,od16},
which would offer great advantages for exploring the internal quantum degrees of freedom of electrons and their potentials in nanoscale-down spin-based electronic applications.

Zinc oxide (ZnO), as a 2D material, can possibly be a promising candidate in the same regard,
since it has been recently reported that a (0001) oriented ultrathin ZnO film can be transformed into a new form of stable graphite-like structure \cite{tusche07,freeman06,tu06,tu10,wu11}.
One to three or four atomic layers of ZnO(0001) have been successfully synthesized on a (111) surface of Ag \cite{tusche07} and Au substrates \cite{deng13,tumino16,lee16},
where Zn and O atoms are arranged in the lateral lattice as a honeycomb shape, similar with the \emph{h}-BN structure.
Theoretical studies have also predicted that thin films of ZnO prefer a graphite-like layered structure when the thickness is less than three atomic layers \cite{tu06,tu10},
and the stability of even thicker graphitic layers can be improved by an epitaxial strain \cite{wu11}.
More remarkably, an epitaxial growth of ZnO monolayer on graphene has been very recently demonstrated in atomic resolution transmission electron microscopy (ARTEM) measurements \cite{hong17}.

Numerous reports, on the other hand, have already shown that intrinsic defects and impurities can readily give rise to magnetism in ZnO allotropies, including monolayer, graphitic layers, thin films, nanowires and nanoribbons \cite{wang08,xu08,botello08,khalid09,topsakal09,kim09,schmidt10,xing11,guo12,tan15,fang16}.
For instance, nonmetal (C, B, N, and K)-\cite{guo12,fang16} and rare-earth metal (Ce, Eu, Gd, and Dy)-doped ZnO monolayer \cite{tan15} exhibits magnetic instability via first-principles calculations.
Moreover, it has been shown by both experimentally \cite{xu08,khalid09,xing11} and theoretically \cite{wang08,topsakal09,kim09} that the Zn vacancy induces ferromagnetism in either thin films \cite{xu08,khalid09,kim09,xing11} and graphitic layers \cite{wang08,topsakal09} of ZnO without the need of magnetic impurities.

In this article, we perform first-principles calculations to reveal the synergistic effect of biaxial strain on magnetism and its spin direction of a graphene-like single layer of ZnO induced by a Zn-vacancy.
Through the analyses of the spin-orbit Hamiltonian matrix elements and atom- and orbital-decomposed magnetic anisotropy energy (MAE),
we elucidate the underlying mechanism of magnetization reversal,
which occurs at the compressive biaxial strain of only $\sim$1--2\%,
in terms of the strain-induced changes in the spin-orbit coupled \emph{d} states of the vacancy-neighbored Zn atom.
Moreover, this anisotropic magnetism is identified to even be reversible upon the change in charge carrier density,
suggesting the possibilities of chemical doping or magnetoelectric control of magnetization reorientation.

\section{\textbf{COMPUTATIONAL DETAILS}}
Density-functional theory (DFT) calculations were performed using the Vienna \emph{ab initio} simulation package (VASP) with plane-wave basis set \cite{Kresse96a,Kresse96b}. The projector augmented wave (PAW) potentials \cite{Blochl994} were used to describe the core electrons, and the generalized gradient approximation (GGA) of Perdew, Burke, and Ernzernhof (PBE) was adopted for exchange-correlation energy \cite{PBE1996}.
As illustrated in Fig. 1(a), a 4$\times$4 lateral unit cell of the \emph{hexagonal} lattice was chosen as a model geometry of a single-layer ZnO.
We also consider an intrinsic defect of an oxygen vacancy [Fig. 1(b)] as well as a zinc vacancy [Fig. 1(c)] in ZnO monolayer,
which occurs notably often during sample preparation and treatments \cite{tuomisto03,tuomisto05,janotti07,chen07,zubiaga08,dong10,khan13}.
In order to avoid spurious interactions between periodic images of ZnO layer,
a vacuum spacing of no less than 15 {\AA} perpendicular to the plane was employed.
An energy cutoff of 500 eV and a 5$\times$5$\times$1 \emph{k} mesh were imposed for the ionic relaxation, where forces acting on atoms were less than 10$^{-2}$ eV/{\AA}.
The MAE is obtained based on the total energy difference when the
magnetization directions are in the \emph{xy} plane (\emph{E}$^\|$) and along the \emph{z} axis (\emph{E}$^\bot$),
MAE $=$ \emph{E}$^\|$ -- \emph{E}$^\bot$, where convergence with respect to the \emph{k} point sampling is ensured.
To obtain reliable values of MAE, the Gaussian smearing method with a smaller smearing of 0.05 and dense \emph{k} points of 15$\times$15$\times$1 were used in noncollinear calculations, where the spin-orbit coupling (SOC) term was included using the second-variation method employing the scalar-relativistic eigenfunctions of the valence states \cite{Koelling1977}.

\begin{figure} [t]
  \centering

   \includegraphics[width=1\columnwidth]{FIG1.png}
  \caption {Top views of the atomic structures for (a) a single-layer ZnO with (b) a O-vacancy and (c) a Zn-vacancy.
  The larger gray and smaller red spheres are Zn and O atoms, respectively. In (b) and (c), the red and black dotted-circles correspond to the O and Zn vacancies, respectively.}
\end{figure}

\begin{table} [b]
\caption{Optimized in-plane lattice \emph{a} (\AA),
formation energy $H_f$ (eV),
magnetic moment $\mu$ ($\mu_\textrm{B}$) and
charge $\rho$ (\emph{e}) of Zn and O atoms
in the nondefective ($V_{\textrm{free}}$)
and defective ZnO with an O- ($V_{\textrm{O}}$) and a Zn-vacancy ($V_{\textrm{Zn}}$).}
\begin{ruledtabular}
\begin{tabular}{ccccccccccc}
 & \emph{a} & $H_f$ & $\mu^{\textrm{Zn}}$ & $\mu^{\textrm{O}}$ &$\mu^{\textrm{Total}}$&
$\rho^{\textrm{Zn}}$ & $\rho^{\textrm{O}}$ \\
\hline
$V_{\textrm{free}}$ & 3.290 & -- &0.0 & 0.0 & 0.0 & 10.81 & 7.18 \\
$V_{\textrm{O}}$ & 3.248 & --3.75 & 0.0 & 0.0 & 0.0 & 10.86 & 7.15 \\
$V_{\textrm{Zn}}$ & 3.296 & --5.43 & 0.04 & 0.32 & 2.0 & 10.79 & 6.97 \\
\end{tabular}
\end{ruledtabular}
\end{table}

\section{\textbf{RESULTS AND DISCUSSION}}
We first optimize the lattice parameters of the nondefective and defective ZnO monolayer and the results are shown in Table I.
The lattice constant of the nondefective ZnO monolayer
is 3.29 \AA, which agrees pretty well with experimental values of $\sim$3.29--3.30 {\AA} \cite{lee16,tusche07,hong17}, superior to the previous theories ($\sim$3.20--3.28 {\AA}) \cite{tu06,topsakal09,tu10}.
While the lattice constant of the nondefective ZnO remains almost unchanged in the presence of a Zn-vacancy,
it is reduced by $\sim$1\% when an O-vacancy is introduced (in a 4$\times$4 unit cell of ZnO monolayer).
The former is due to the direct repulsive interaction between O--O anions without Zn mediation.
As presented in Table I, the further calculations of formation energy $H_f$, defined as $H_f = E_{\textrm{Nondef}}-E_{\textrm{Def}}-E_{\textrm{Zn/O}}$, where $E_{\textrm{Nondef}}$, $E_{\textrm{Def}}$, and $E_{\textrm{Zn/O}}$ represent the total energies of the nondefective and defective ZnO with a Zn/O-vacancy, and a Zn/O atom in its bulk/O$_2$-molecule form,
indicate the favorable formation of a vacancy in ZnO, particularly for the Zn-site vacancy, as reported in experiments \cite{tuomisto03,tuomisto05,janotti07,chen07,zubiaga08,dong10,khan13}.

\begin{figure} [t]
  \centering
 \includegraphics[width=1\columnwidth]{FIG2.png}
  \caption {Spin- and orbital-decomposed PDOS of (a) Zn and (b) O atoms of the nondefective ZnO monolayer. The same PDOS for the defective ZnO with (a) and (b) a O-vacancy and (e) and (f) a Zn-vacancy.
  The singlet $a_1$ ($d_{z^2}$) and doublet $e'$ ($d_{xy/x^2−y^2}$) and $e''$ ($d_{xz/yz}$) orbital states of Zn atom are shown in blue, black, and red lines, respectively. The $p_{x,y}$ and $p_{z}$ orbital states of O atom are shown in blue and red lines, respectively.
For reference, the total DOS per atom (shaded area) is also shown in each panel. The Fermi level is set to zero in energy.}
\end{figure}

The significance of induced magnetism, as fairly addressed in experimental \cite{xu08,khalid09,xing11} and previous theoretical studies \cite{wang08,topsakal09,kim09}, is then justified.
Note that the magnetic moments of the first-neighboring Zn and O atoms to the Zn-vacancy site
are the most prominent beyond those (negligibly small) of the other sites away from the vacancy.
We thus refer the results and discussion here and hereafter mainly to those corresponding to the first-neighboring Zn and O atoms.
As seen in Table I, the nondefective and defective ZnO with an O-vacancy behave nonmagnetic nature as in bulk \cite{hong07,chen13},
but the Zn and O atoms in the Zn-vacancy defected ZnO exhibit significant spin magnetic moments of 0.04 and 0.32 $\mu_\textrm{B}$, respectively.
This vacancy-induced magnetic instability
simply reproduces the results reported in previous \emph{ab initio} calculations, where the total magnetization in an unit cell is 2 $\mu_\textrm{B}$ \cite{wang08,topsakal09,kim09}.
Moreover, our results in ZnO monolayer are more and less comparable with those ($\mu^{\textrm{Zn}}$ $=$ $\sim$0.04 and $\mu^{\textrm{O}}$ $=$ $\sim$0.14--0.48 $\mu_\textrm{B}$) in ZnO thin films \cite{wang08,kim09}.
The vacancy defect driven magnetism can be understood in terms of the spin-polarized charge excess or deficit, where the total magnetization in an unit cell is determined by the unpaired electron counts.
Our additional calculations with an injected charge show that
such a charge deficit is compensated by an extra charge of +2$e^{-}$ per Zn-vacancy and thereby the induced magnetism disappears.

\begin{figure} [t]
  \centering

 \includegraphics[width=1\columnwidth]{FIG3.png}
  \centering
  \caption{(a) Magnetic moments $\mu$ of Zn (blue) and O (black) atoms and (b) MAE of a single-layer ZnO with a Zn-vacancy as a function of strain ($\varepsilon$). The negative and positive $\varepsilon$ correspond to the compressive and tensile strain, respectively.
  The insets in (b) show the schematics for magnetization switching from in-plane (bottom in negative MAE) to out-of-plane magnetization (top in positive MAE) at a compressive biaxial strain of $\sim$1--2\% .}
\end{figure}

To better understand the Zn-vacancy induced magnetism,
we show the projected density of states (PDOS) of the Zn \emph{d} and O \emph{p} orbital states in Figs. 2(e) and 2(f), respectively, in comparison to those
of the nondefective [Figs. 2(a) and 2(b)] and O-vacancy defected ZnO [Figs. 2(c) and 2(d)]. In hexagonal lattice, the five \emph{d}-orbital states split into one singlet
$a_1$ or $m=0$ ($d_{z^2}$) and two doublets $e'$ or $m=\pm2$ ($d_{xy/x^2−y^2}$) and $e''$ or $m=\pm1$ ($d_{xz/yz}$).
Both the nondefective and O-vacancy defected ZnO exhibit wide-bandgap insulating electronic features,
where the Zn \emph{d} orbital states are fully occupied while the conduction band is mainly the characteristics of the O \emph{p} states.
Furthermore, the spin subbands of the majority- and minority-spin states are essentially degenerate.
On the other hand, for the Zn-vacancy defected ZnO the majority- and minority-spin states of the Zn and O atoms reflect unequally (i.e., net exchange splitting); The minority-spin filled bands of the Zn $e'$ and O $p_{x,y}$ states move simultaneously upward across the Fermi level in energy relative to the corresponding majority-spin states,
thus exhibits a half-metallic behavior.
This is mainly due to the charge delocalization of the minority-spin electrons of the Zn and O atoms due to the broken charge quantization associated with the absence of a Zn atom,
which are most likely distributed toward the Zn-vacancy site.
From the charge analysis taken account within the Wigner-Seitz atomic radius (Table I), the charges of 0.02 and 0.21$e^{-}$ are depleted from the Zn and O atoms, respectively.

As the stability of the layered ZnO can be enhanced by an epitaxial strain \cite{wu11},
we next explore the impact of biaxial compressive and tensile strains ($\varepsilon$) on magnetism of the Zn-vacancy defected ZnO monolayer,
where $\varepsilon$ $=$ (\emph{a$_{\|}$} -- \emph{a$_0$})/\emph{a}$_0$$\times$100\%;
$a_0$ is the equilibrium lattice (3.296 \AA) of the Zn-vacancy defected ZnO and $a_{\|}$ is a variable.
We here would like to note that either in-plane lattice expansion and compression up to $\sim$3\% of graphitic ZnO layers
can be practically achieved by controlling the number of atomic layers \cite{lee16}.
The magnetic moments of Zn (blue symbol) and O (black symbol) atoms are shown in Fig. 3(a) for different strains, ranging from a compressive strain of --3\% to a tensile strain of $+$3\%.
The Zn moment decreases as strain increases (compressive $\rightarrow$ tensile), whereas it increases for the O case, which is a reflection of Vergard's law \cite{vegard21}.

\begin{table} [t]
\caption{The atomic contributions of Zn and O species to the total MAE (meV) of a single-layer ZnO with a Zn-vacancy for different strains $\varepsilon$ (\%),
where the corresponding in-plane lattice constants \emph{a} (\AA) are also listed.}
\begin{ruledtabular}
\begin{tabular}{cccccc}
$\varepsilon$ & \emph{a} & MAE(Zn) & MAE(O) & MAE(Total) \\
\hline
--2 & 3.230 & $+$0.05 & --0.01 & $+$1.06  \\
0 & 3.296 & --0.05 & --0.01 & --0.50  \\
2 & 3.361 & --0.03 & --0.01 & --0.29  \\
\end{tabular}
\end{ruledtabular}
\end{table}

Figure 3(b) shows the calculated MAE of the Zn-vacancy defected ZnO versus strain.
The result of the negative MAE (--0.5 meV) at zero strain indicates that the induced magnetism prefers an in-plane spin orientation parallel to the lateral lattice.
Both the compressive and tensile strains enhance MAE.
More interestingly, the sign of MAE is reversed into the positive and reaches up to $+$1 meV at --2\% strain.
This indicates that the magnetization direction can reorient
from in-plane to perpendicular orientation (PMA) upon an application of compressive strain.
Such a PMA remains at strain up to the largest strain of --3\% considered in our study.
On the other hand, the MAE tends to disappear (or reorient)
at the larger tensile strain beyond $+$3\% as the Zn moment reaches zero.

\begin{figure} [t]
  \centering

 \includegraphics[width=1\columnwidth]{FIG4.png}
  \centering
  \caption{Difference of the orbital-projected SOC energies, $\Delta E_{\textrm{soc}}$, between in- and out-of-plane magnetization orientation of Zn and O atoms in a single-layer ZnO with a Zn-vacancy for the (a) and (b) --2\% and (c) and (d) zero and (e) and (f) $+$2\% strain, respectively. The yellow and blue bars represent the positive and negative $\Delta E_{\textrm{soc}}$, respectively.}
\end{figure}

The contributions of different atomic species (Zn and O) to the total MAE are shown in Table II for the selected strains of $\varepsilon$ $=$ --2, 0, and $+$2\%,
from the difference of SOC energies between in- and out-of-plane magnetization, $\Delta E_{\textrm{soc}}$ $=$ $E_\textrm{soc}^\|$ -- $E_\textrm{soc}^\bot$.
Here, $E_{\textrm{soc}} = <\frac{\hbar^2}{2m^2c^2}\frac{1}{r}\frac{dV}{dr} \textbf{\emph{L}}\cdot\textbf{\emph{S}}>$,
where \emph{V}(\emph{r}) is the spherical part of the effective potential within the PAW sphere, and \textbf{\emph{L}} and \textbf{\emph{S}} are
orbital and spin operators, respectively.
The expectation value of $E_{\textrm{soc}}$ is twice the actual value of the total energy correction to the second-order in SOC, i.e., MAE $\approx$ 1/2$\Delta E_{\textrm{soc}}$ \cite{antropov14,zhang17}.
Our test calculations indicate that the second-order perturbation theory is a reasonable approximation as the total MAE overall agree within a few percent accuracy
with those obtained from the atom and orbital projected calculations.
The other 50\% of the SOC energy translates into the crystal-field energy and the formation of the unquenched orbital moment \cite{skomski11}.

\begin{figure} [t]
  \centering

 \includegraphics[width=1\columnwidth]{FIG5.png}
  \caption {Spin- and orbital-decomposed PDOS of Zn and O atoms of the defective ZnO monolayer with a Zn-vacancy for the (a) and (b) $\varepsilon=-2\%$ and (c) and (d) $\varepsilon=0$ and (e) and (f) $\varepsilon=+2\%$, respectively. The singlet $a_1$ ($d_{z^2}$) and doublet $e'$ ($d_{xy/x^2−y^2}$) and $e''$ ($d_{xz/yz}$) orbital states of Zn atom are shown in blue, black, and red lines, respectively. The $p_{x,y}$ and $p_{z}$ orbital states of O atom are shown in blue and red lines, respectively. The Fermi level is set to zero in energy.}
\end{figure}

As seen in Table II, the $\Delta E_{\textrm{soc}}$ of Zn atom exhibits a trend similar to the changes in sign of the total MAE under strain, whereas the contribution of O atom is rather small and strain-independent.
Thus, the dominant contribution to the MAE must come from the Zn atom, where the larger SOC effect of \emph{d} orbitals ($l = 2$) than \emph{p} orbitals ($l=1$) as a physics origin of anisotropic phenomena is well manifested.
The validity of this argument is further demonstrated in connection with the Bruno theory
$ \textrm{ MAE} = \frac{\xi}{4\mu_B}\Delta \mu_{\textrm{o}}$ \cite{Bruno1989},
where $\Delta \mu_{\textrm{o}}$ $=$ $\mu_\textrm{o}^\bot - \mu_\textrm{o}^\|$ is the orbital magnetic anisotropy (OMA) and $\xi$ is the SOC constant.
The calculated $\Delta \mu_{\textrm{o}}$ of the Zn and O atoms are --0.6$\times$ and 0.8$\times$10$^{-2}$ $\mu_{B}$ at zero strain, respectively.

To get more insights on the strain-induced magnetization reversal, we further decompose the $\Delta E_{\textrm{soc}}$
into the \emph{d}-orbital matrix of the Zn atom for --2, 0, and $+$2\% strains in Fig. 4.
For comparison, we also plot the same for the \emph{p}-orbital matrix of O atom (right panels in Fig. 4).
The O \emph{p}-orbital contributions of MAE are relatively small and insensitive,
compared with the Zn \emph{d}-orbital contributions,
to the sensitivity of MAE under strain effect, as discussed in atom-decomposed MAE shown in Table II.
We thus focus on the \emph{d}-orbital states of the Zn atom for the discussion of the electronic origin for magnetization reversal.

\begin{figure} [t]
  \centering

 \includegraphics[width=1\columnwidth]{FIG6.png}
  \caption {Energy levels of Zn $e'$ orbital states under the compressive ($\varepsilon$ $=$ $-$2\%) and expansive ($\varepsilon$ $=$ $+$2\%) strain.
  The upward red and downward blue arrows denote the majority-spin and
minority-spin states, respectively. The thin horizontal line indicates the Fermi level ($E_F$). The minority-spin $e'$ states are not much perturbed under strain.
The spin-orbit coupling pairs near the Fermi level between the occupied and unoccupied states are emphasized, where the spin up-up and up-down couplings
provide the positive (by $\hat{L}_{z}$) and negative contribution (by $\hat{L}_{x}$) to the MAE, respectively.
The atomic symbols at the top (crystallographic image) are the same as used in Fig. 1.}
\end{figure}

\begin{figure} [b]
  \centering

 \includegraphics[width=1\columnwidth]{FIG7.png}
  \caption {(a) Magnetic moments $\mu$ of Zn (blue) and O (black) atoms and (b) MAE as a function of the number of extra valence electrons $n_e$ ($e^{-}$) in a 4$\times$4 unit cell of ZnO monolayer with a Zn-vacancy.
  The negative and positive values in $n_e$ correspond to the charge deficit and addition, respectively.
    The insets in (b) show the schematics for magnetization switching from in-plane (bottom in negative MAE) to out-of-plane magnetization (top in positive MAE) at $n_e=-0.5e^{-}$ and $n_e=+1e^{-}$.}
\end{figure}

In the second-order perturbation theory, the MAE is determined by the SOC between occupied and unoccupied \emph{d}-orbital states \cite{wang93};
$
  \textrm{MAE}^{\sigma\sigma'} = \xi^2 \sum_{o,u}
  \frac{|\langle \Psi_{o, \sigma} | \hat{L}_{z} | \Psi_{u, \sigma'} \rangle |^2 - {|\langle \Psi_{o, \sigma} | \hat{L}_{x} | \Psi_{u, \sigma'} \rangle |^2}}
  { E_{u, \sigma'} - E_{o, \sigma} },
$
where $\Psi_{o}$ ($\Psi_{u}$) and $E_{o}$ ($E_{u}$) are the eigenstates and eigenvalues of occupied (unoccupied) states (band index \emph{n} and wave vector \textbf{k} are omitted for simplicity) for each spin state, $\sigma, \sigma'$ $=$ $\uparrow, \downarrow $, respectively, and $\hat{L}_{x(z)}$ is the \emph{x} (\emph{z}) component of the orbital angular momentum operator.
For $\sigma\sigma'$ $=$ $\uparrow \uparrow$ or $\downarrow \downarrow$, the positive (negative) contribution to MAE is determined by the
SOC interaction between the occupied and unoccupied states with the same (different by one) magnetic quantum number (\emph{m}) through the $\hat{L}_{z}$ ($\hat{L}_{x}$) operator.
Relative contributions of the nonzero $\hat{L}_{z}$ and $\hat{L}_{x}$ matrix elements are $ \langle xz|\hat{L}_{z}|yz\rangle $ $=$ 1, $ \langle xy|\hat{L}_{z}|x^2-y^2\rangle $ $=$ 2,
$ \langle xz,yz|\hat{L}_{x}|z^2\rangle $ $=$ $\sqrt{3}$, $ \langle xz,yz|\hat{L}_{x}|xy\rangle $ $=$ 1, and $ \langle xz,yz|\hat{L}_{x}|x^2-y^2\rangle $ $=$ 1 \cite{wang93}.
For $\sigma\sigma'$ $=$ $\uparrow \downarrow $, MAE has the opposite sign, so the positive (negative) contribution comes from the $\hat{L}_{x}$ ($\hat{L}_{z}$) coupling.

\begin{figure} [t]
  \centering

 \includegraphics[width=1\columnwidth]{FIG8.png}
  \caption {Spin- and orbital-decomposed PDOS of Zn and O of the defective ZnO monolayer with a Zn-vacancy for the (a) and (b) $n_e = -1e^{-}$ and (c) and (d) $n_e = 0$ and (e) and (f) $n_e = +1e^{-}$, respectively. The singlet $a_1$ ($d_{z^2}$) and doublet $e'$ ($d_{xy/x^2−y^2}$) and $e''$ ($d_{xz/yz}$) orbital states of Zn atom are shown in blue, black, and red lines, respectively. The $p_{x,y}$ and $p_{z}$ orbital states of O atom are shown in blue and red lines, respectively. The Fermi level is set to zero in energy.}
\end{figure}

Figures 5(a)--5(f) show the \emph{d} (\emph{p})-orbital PDOS of the Zn (O) atom at $\varepsilon = -2$, 0, and $+$2\%, respectively.
In Fig. 6, we also illustrate simple energy level diagram for the changes of Zn $e'$ peak states
near the Fermi level under strain.
The electronic structure analyses along with the \emph{d}-orbital decomposed MAE indicate that the underlying origin of the negative MAE at zero strain
is the spin-orbit coupling between the majority-spin occupied and the minority-spin unoccupied $e'$ states through the in-plane orbital angular momentum operator,
$ \langle xy _ \uparrow| \hat {L}_{x}|x^2-y^2 _ \downarrow \rangle $ (or $ \langle x^2-y^2  _ \uparrow| \hat {L}_{x}|xy_ \downarrow \rangle $).
This negative contribution weakens at $+$2\% strain because of the downward shift of the majority-spin filled $e'$ states away from the Fermi level.
For the compressive strain,
the positive MAE arises primarily from the same spin-channel ($\uparrow \uparrow$) matrix elements by
$ \langle xy _ \uparrow| \hat {L}_{z}|x^2-y^2 _ \uparrow \rangle $
and $ \langle x^2-y^2 _ \uparrow| \hat {L}_{z}| xy_ \uparrow \rangle $, as the majority-spin degenerate $e'$ states become partially unoccupied.
Such a modification of the energy landscape of the Zn $e'$ electronic state
around the Fermi level is the result of the changes in the strength of orbital hybridization with the in-plane O $p_{x/y}$ states under strain.

The engineering of MAE by Fermi level shifts further suggests to explore crucial effects of the external gating and chemical doping on magnetization reversal.
To anticipate this phenomenon, we have performed electrostatic doping calculations
by changing the number of valence electrons ($n_e$) in the whole system (4$\times$4 unit cell of ZnO monolayer)
from $n_e=-1.5$ to $+1.5e^{-}$.
This approach reflects to the charge carrier density that is confined over the metallic surface rather than exceeding to a vacuum region, which is analogues to the electric field effect by the positive and negative gating as well as chemical doping with the neighboring elements (e.g., Cu or Ga for Zn, and N or F for O site) \cite{zhang17}.

The calculated magnetic moments and MAE as a function of $n_e$ are shown in Figs. 7(a) and 7(b), respectively.
It is expected that the charge deficit enhances the magnetic moments of both Zn and O atoms
whereas the charge excess reduces because of the electron/hole doping.
Nevertheless, as shown in Fig. 7(b), both the deficit and excess valence charges promptly enhance the MAE
and change its sign from in-plane to perpendicular magnetization at $n_e=-0.5$ and $=+1e^-$, respectively.
These results are of considerable interest in the area of the voltage-controlled magnetic anisotropy (VCMA) and magnetoelectric phenomena.
In the present study, the VCMA can be determined by the curve slope of the electric-field or $n_e$ dependence of MAE.
The nonlinear parabola-like shape behavior of the VCMA within a small range of $n_e$ ($-1<n_e<+1$)
is quite typical in most magnetoelectric materials, as reported in previous calculations \cite{zhang17,odkhuu15,odkhuu16} and experiments \cite{raj13,noz13,noz16}.

In Figs. 8(a)--8(f), we plot the \emph{d} (\emph{p})-orbital PDOS of the Zn (O) atom for $n_e=-1$, 0, and $+1e^{-}$, respectively.
The entire band structures of both Zn and O atoms shift toward the low-energy region upon the presence of external electrons from $-1$ to $+1e^{-}$ due to the band filling effect.
Since the majority-spin Zn $e'$ peak exists at the Fermi level, the similar argument discussed for the --2\% strained ZnO applies
to the origin of PMA at $n_e=-1e^{-}$.
On the other hand, for $n_e=+1e^{-}$
the peak position of the minority-spin empty $e'$ bands appears just above the Fermi level while the majority-spin states are fully occupied, as shown in Fig. 8(e).
Thus the SOC pair between the minority-spin filled and empty $e'$ states around the Fermi level
leads to the positive MAE through
$ \langle xy _ \downarrow| \hat {L}_{z}|x^2-y^2 _ \downarrow \rangle $
and $ \langle x^2-y^2 _ \downarrow| \hat {L}_{z}| xy_ \downarrow \rangle $.

\section{\textbf{CONCLUSION}}
In summary, we have performed the first-principles electronic structure calculations to reveal the critical effects of strain and external gating on magnetism and magnetic anisotropy of a single-layer ZnO with a Zn-vacancy.
In contrast to the pristine and defective ZnO with an O-vacancy, the presence of a Zn-vacancy
induces significant magnetic moments to its first neighboring O and Zn atoms due to the charge deficit.
More remarkably, we show that a compressive biaxial strain of only $\sim$1--2\% or either negative and positive gate voltages
give rise to magnetization reorientation from in- to out-of-plane magnetization.
The further analyses provide the underlying mechanism for magnetization reversal,
which is the result of the strain- and electric field-induced changes in spin-orbit coupled Zn \emph{d}-orbital states hybridized with O \emph{p} states.
These predictions may inspire further experimental and theoretical explorations
of exploiting epitaxial biaxial strain and electric field to manipulate magnetism and magnetization direction for two-dimensional spintronic applications.

This work was supported by Postdoctoral Research Program (2017) through Incheon National University.


text:
\section{Supplementary material}
In this supplementary material, we provide more details regarding the qualitative understanding and the properties of the FFLO$_x$ state under the NIST scheme.

As we have discussed in the main text, the various FFLO$_x$ states can be qualitatively understood as shifts of the local minima to finite $Q_x$ induced by Fermi surface asymmetries along the $x$-axis. This can be more easily understood by comparing the $Q=0$ phase diagram (see Fig. \ref{figsupp1}(a)) with the phase diagram in the main text (Fig. 3(a) therein). In the presence of a transverse Zeeman field, up to three local minima can be found in the thermodynamic potential for $Q=0$, corresponding to the normal, the nSF2 and the SF states, respectively (see Fig. \ref{figsupp1}(b)). This is reflected in the various first-order boundaries in Fig. \ref{figsupp1}(a). When we consider finite center-of-mass momentum pairing, these local minima will be shifted into the plane of finite $Q_x$. As a result, the first-order boundaries in Fig. 3(a) of the main text are reminiscent of the first-order boundaries for $Q=0$ in Fig. \ref{figsupp1}(a), but are shifted slightly. Of course, the phase diagram would be further modified by the competition from the FFLO$_y$ states, and the various local minima along $\mathbf{Q}=Q_x \hat{x}$ may merge as the parameters are tuned, which generates the end points we see in Fig. 3(a) of the main text. As an example, we plot in Fig. \ref{figsupp1}(c) a typical thermodynamic potential with $Q=0$ and in Fig. \ref{figsupp1}(d) the corresponding thermodynamic potential landscape when finite $Q_x$ pairing is considered. It appears that the original two local minima are both shifted into the finite $Q_x$ plane.

\begin{figure}
\includegraphics[width=6cm]{figsupp1a.pdf}
\includegraphics[width=6cm]{figsupp1b.pdf}
\includegraphics[width=6cm]{figsupp1c.pdf}
\includegraphics[width=6cm]{figsupp1d.pdf}
\caption{(Color online) (a) Phase diagram on the $\mu$--$\alpha$ plane for $E_b/h=0.2$, $h_x/h=0.2$. Only zero center-of-mass momentum pairing states are considered here. The solid curves are first-order boundaries, while the dashed curves represent phase boundaries of continuous phase transitions.
(b) Typical thermodynamic potentials with up to three local minima, with the parameters:
$\alpha k_h/h=0.33$, $\mu/h=4.53$ (solid curve); $\alpha k_h/h=0.32$, $\mu/h=4.53$ (dotted curve); $\alpha k_h/h=0.42$, $\mu/h=4.2$ (dash-dotted curve); $\alpha k_h/h=0.33$, $\mu/h=4.6$ (dashed curve); while $E_b/h=0.2$, $h_x/h=0.2$ for all cases.
(c) Thermodynamic potential for $Q=0$ and $E_b/h=0.2$, $h_x/h=0.2$, $\alpha k_h/h\sim2.12$, $\mu/h\sim-0.98$. The local minima are located at $\Delta/h\sim 0.3$ and $\Delta/h=0$. (d) Thermodynamic potential as a function of $\Delta_Q$ and $Q_x$ for $E_b/h=0.2$, $h_x/h=0.2$, $\alpha k_h/h\sim2.12$, $\mu/h\sim-0.98$. The local minima are located at ($\Delta_{Q}/h\sim 0.29$, $Q_x/k_h\sim -0.11$) and ($\Delta_{Q}/h\sim 0.16$, $Q_x/k_h\sim -0.35$).}
\label{figsupp1}
\end{figure}

A remarkable feature of the phase diagram in the presence of a transverse Zeeman field and SOC is the enhanced stability region of the FFLO$_x$ state which seems to replace the SF state in the $Q=0$ phase diagram. As we have discussed in the main text, this is due to the Fermi surface deformation. As the Fermi surface deformation becomes smaller with increasing chemical potential or with increasing SOC strength, it is natural to expect that the energy gain by adopting a finite center-of-mass momentum pairing in the FFLO$_x$ state should become less. This can be seen by a direct scaling of the pairing parameters $Q_x$ and $\Delta_Q$ of the FFLO$_x$ state with increasing SOC strength (see Fig. 4(a) in the main text) or with increasing chemical potential (see Fig. \ref{figsupp2}(a)). Apparently, while the pairing order parameter $\Delta_Q$ keeps increasing with chemical potential and SOC strength, the magnitude of the center-of-mass momentum of the FFLO$_x$ state becomes exponentially small in the large $\mu$ and/or large $\alpha$ limit. Furthermore, the difference in thermodynamic potential between the SF state and the FFLO$_x$ state also becomes exponentially small. Therefore in real experiments, the parameter region where the FFLO$_x$ states can be observed should feature moderate SOC and moderate chemical potential.

\begin{figure}
\includegraphics[width=6cm]{figsupp2a.pdf}
\includegraphics[width=6cm]{figsupp2b.pdf}
\caption{(Color online) (a) Evolution of the center-of-mass momentum and the pairing order parameter (inset) for the ground state as the chemical potential increases. (b) Evolution of the difference in thermodynamic potential between the SF state and the FFLO$_x$ state with increasing chemical potential, with $\delta\Omega=\Omega_{\rm SF}-\Omega_{\rm FFLO_x}$. The parameters are: $E_b/h=0.2$, $h_x/h=0.2$, and $\alpha k_h/h=1.9$. }
\label{figsupp2}
\end{figure}

\end{appendix}
\end{widetext}

\end{document}