\section{Introduction}
\IEEEPARstart{T}{he} penetration of
new renewable sources of electrical energy is currently increasing in most electric power grids worldwide, as more and more traditional power plants are phased out.
A major concern is obviously that this substitution reduces the available inertia while it simultaneously induces
 larger fluctuations in power generation~\cite{Ulb14}. Both changes may jeopardize grid stability, either individually or taken together.
A key issue is accordingly to evaluate how much  power grids need to be adapted
to their resulting new modes of operation -- for instance through line extensions or deployment of resources providing ancillary services~\cite{AnnualEnOutlook}.
To ensure the stability of the grid and the safety of power supply, it is important to clarify the role of the generators dynamical parameters that will be affected by this transition, namely the
rotational
inertia and the frequency damping / droop control. Both are going to be globally reduced, moreover their geographical distribution will be modified. To try to identify
which of these dynamical parameters are most crucial, where they should be primarily deployed,  as well as when and in what operating conditions they are most needed, analytical and
numerical works have investigated the robustness of electric power grids under external disturbances.
The response of power grids to external disturbances has been investigated through
quadratic performance metrics~\cite{Bam13,Dor14a,Teg15,Sia16,Poo17,Pag17,Gru18,Col19,Pag19b, Yan19},
eigenvalue damping ratios and frequency overshoots~\cite{Mev16}, rate of change of frequency\cite{Pag17,Guo18} or
disturbance wave propagation~\cite{Tam18,Hae18}. Except for numerical results~\cite{Mev16,Tam18,Hae18},
these works considered disturbances with infinitely short time scales, such
as white noise power fluctuations or instantaneous power injection changes~\cite{Bam13,Teg15,Sia16,Poo17,Gru18,Col19,Pag19b, Yan19},
or with infinitely long time scales, such as step changes in power injection~\cite{Pag17,Guo18}. All of the analytical works
relied on one of the two homogeneity assumptions that
inertia and frequency damping are the same everywhere, or that their ratio is.

These homogeneity assumptions
are not representative of real power grids, where in particular, consumer nodes are inertialess but with small, albeit finite frequency damping~\cite{Ber81,Mac08}.
Often, this inconsistency  is circumvented by invoking a prior Kron reduction absorbing the inertialess nodes into an effective network.
One then measures the robustness of that reduced network, which may or may not be related to the robustness of the original one, because Kron reduction does not
capture the dynamics of the reduced, inertialess nodes.
To the best of our knowledge, Ref.~\cite{Pag19b} is the first work that tolerates deviations from homogeneity in an analytical calculation of a quadratic
performance metric. Its results suggest that grid robustness is crucially sensitive to
the geographic distribution of frequency control, while inertia has to be distributed rather evenly. This conclusion has to be revisited because,
first, Ref.~\cite{Pag19b} is based on an approximate method tolerating only small deviations from homogeneity and second, the only disturbances it considers
are long power losses.

In this manuscript, we investigate the response of power systems to colored noisy power fluctuations.
We quantify the response to these disturbances by a quadratic performance metric measuring the primary control effort
necessary to absorb the fault~\cite{Poo17}. Our analytical approach still relies on a homogeneity assumption.
Our results however emphasize the role played by the different time scales in the problem:
the performance metric depends on inertia only when the characteristic time scale of the disturbance is long compared to all other time scales in the system.
In that case we conjecture, and confirm numerically, that our analytical results also apply to heterogeneous systems.
In continental-size
high-voltage power grids, the network time scales are shorter than few seconds~\cite{Tyl18b}. Therefore, noise fluctuating on time scales of tens of second or
more induces transients whose amplitude, duration and oscillations are largely independent of inertia, and our analytical results directly apply to most
disturbances on realistic, inhomogeneous high voltage power grids.

The manuscript is organized as follows. Section \ref{notation} defines our mathematical notations.
Section \ref{section2} defines the power network model and gives an analytical expression for its linear response.  Section \ref{section3} introduces performance metrics and gives analytical expressions for them. Of particular interest are the short- and long-correlation time asymptotics. In Sec.~\ref{section5}, we numerically confirm our theory on both the IEEE 118-Bus test case and the PanTaGruEl model of the synchronous grid of continental Europe. We discuss time scales in such
high-voltage power grids and show that the corresponding performance metric is given by the long noise correlation time asymptotic limit.
Our conclusions are given in section \ref{section6}.

\section{Mathematical Notation}\label{notation}
Given a vector $\bm v\in \mathbb{R}^n$, we denote its transpose by $\bm v^\top$. We write   $\bm M = {\rm diag}\{m_i\}$
for the diagonal $n \times n$ matrix with $m_1,m_2,...,m_n \in \mathbb{R}$ on its diagonal. The $j$-th unit vector with a single nonzero component is
$(\hat{e}^j)_i=\delta_{ij}$\,. The scalar product of two vectors $\bm u, \bm v\in \mathbb{R}^n$ is written $\bm u^\top\bm v$ and the scalar product of a vector with itself
is $\bm v^2=\bm v^\top\bm v$. The statistical average of a random variable $x\in\mathbb{R}$ is $\overline{x}$\,. Finally, considering a diagonal matrix $\bm M$\,, we denote
its $p^{\rm th}$ power as $\bm M^p={\rm diag}\{m_i^p\}$\,.

\section{Power Grids and their Response to Fluctuating Power Injections}\label{section2}
\subsection{Swing dynamics near synchrony}
Transient dynamics in high-voltage power networks is commonly modelled by the swing equations which describe the dynamics of voltage angles assuming constant voltage amplitudes. In the lossless line approximation, appropriate to very high-voltages~\cite{Mac08},
they read
\begin{eqnarray}\label{eq:kuramoto}
m_i\,\dot{\omega}_i+d_i\,{\omega}_i=P_i-\sum_j b_{ij}\sin(\theta_i-\theta_j),
\end{eqnarray}
where each network node is labeled $i=1,...,n$ with a voltage angle $\theta_i$. Equation~\eqref{eq:kuramoto} is written
in a rotating frame, so that the frequency $\omega_i=\dot{\theta}_i$ refers to the deviation from the rated frequency of $50$ or $60$ Hz. Each node has inertia and damping control parameters $m_i$ and $d_i$ respectively, and an active power $P_i$ that is generated ($P_i>0$) or consumed ($P_i<0$). The coupling between node $i$ and $j$
is given by the susceptance $b_{ij} \ge 0$ of the corresponding power line. The operational state ${\bm \theta}^{(0)}$ is a synchronous
stationary solution to Eq.~(\ref{eq:kuramoto}).

Equation~(\ref{eq:kuramoto}) is governed by two sets of time scales. The first set is given by the ratio between inertia and damping coefficients $\gamma_i^{-1}\equiv m_i/d_i$. It corresponds to the local relaxation of synchronous machines. The second one is determined by the network characteristic time scales $d_i/\lambda_\alpha$
given by damping coefficients and the eigenvalues of the weighted Laplacian, see Eq.~(\ref{eq:laplacian}) below.
Depending on these two sets of time scales, perturbations are locally damped or spread across the network. In a synthetic synchronous grid of continental Europe
with constant damping and inertia corresponding to the average of their true values,
it has been found that
all modes are underdamped and propagate through the whole system with $d/\lambda_\alpha < \gamma^{-1}$, $\forall \alpha$~\cite{Tyl18b}.

We next investigate the response of the system to a time-dependent disturbance ${\bm{P}}(t) = {\bm{P}} + \delta {\bm{P}}(t)$ acting on
the operational state ${\bm \theta}^{(0)}$, following which angles become time-dependent,
${\bm{\theta}}(t) = {\bm{\theta}}^{(0)} + \delta {\bm{\theta}}(t)$. The small-signal response is governed by dynamical equations obtained after linearizing Eq.~(\ref{eq:kuramoto}) about
$\bm{\theta}^{(0)} $,
\begin{align}\label{eq:kuramoto_lin}
{\bm M}\,\dot{\bm \omega} + {\bm D}\, {\bm \omega} &= \delta {\bm P}(t) - {\mathbb L}({\bm\theta}^{(0)}) \, \delta {\bm \theta} \, ,
\end{align}
where we introduced inertia and damping matrices, ${\bm M}={\mathrm{diag}}\{m_i\}$ and ${\bm D}={\mathrm{diag}}\{d_i\}$ and the weighted Laplacian matrix ${\mathbb L}(\{ \theta_i^{(0)} \})$ with matrix elements
\begin{equation}\label{eq:laplacian}
{\mathbb L}_{ij} =
\left\{
\begin{array}{cc}
-b_{ij} \cos(\theta_i^{(0)} - \theta_j^{(0)}) \, , & i \ne j \, , \\
\sum_k b_{ik} \cos(\theta_i^{(0)} - \theta_k^{(0)}) \, , & i=j \, .
\end{array}
\right.
\end{equation}
This Laplacian is minus the stability matrix of the linearized dynamics, and
since we consider a stable synchronous state, it
is positive semidefinite, with a single vanishing eigenvalue $\lambda_1=0$ with eigenvector ${\bm u}_1=(1,1,...1)/\sqrt{n}$. All other eigenvalues are positive,
$\lambda_\alpha>0$, $\alpha=2,3,...n$.
Equation~(\ref{eq:kuramoto_lin}) can be integrated via a spectral decomposition provided either: (i) both $\bm M$ and $\bm D$ commute with $\mathbb{L}$, then Eq.~(\ref{eq:kuramoto_lin}) can be integrated in the eigenspace of $\mathbb{L}$, or (ii)  ${{\bm M}^{-1}{\bm D}}=\gamma\,\mathbb{I}$\,, in which case Eq.~(\ref{eq:kuramoto_lin}) can be integrated in the eigenspace of ${\bm D^{-1/2}}\,\mathbb{L}\,{\bm D^{-1/2}}$.

\subsection{Analytical solution for constant damping-to-inertia ratio}
We consider the case (ii) above of constant  inertia-to-damping ratio, $m_i/d_i=\gamma^{-1}$ $\forall i$\,. To calculate the response of the system,
we first perform a change of variable ${\delta\bm \varphi}={\bm D}^{1/2}\delta {\bm \theta}$
on Eq.~(\ref{eq:kuramoto_lin}). We obtain
\begin{align}\label{eq:kuramoto_lin3}
\gamma^{-1}\delta \ddot{{\bm \varphi}} + \delta \dot{{\bm \varphi}} &= {\bm D}^{-1/2}\delta {\bm P} -  {\bm D}^{-1/2}\,\mathbb{L}\,{\bm D}^{-1/2}\, \delta {\bm \varphi} \, .
\end{align}
We choose this normalization with $\bm D$, rather than with $\bm M$ as proposed in Ref.~\cite{Pag17}, because it allows us to treat
the inertialess case with $\gamma^{-1}=0$, from which we will extrapolate
the realistic case where consumer nodes are inertialess and generator nodes have nonhomogeneous inertia.
Equation~(\ref{eq:kuramoto_lin3}) can be solved by expanding angle deviations as $\delta{\bm \varphi}(t)=\sum_{\alpha}c_\alpha(t)\,{{\bm u}_\alpha^{D}}$,
over the eigenvectors ${\bm u}_{\alpha}^{D}$ of $\mathbb{L}^D={\bm D}^{-1/2}\,\mathbb{L}\,{\bm D}^{-1/2}$.
The matrix $\mathbb{L}^D$ is no longer Laplacian but it still has a zero-mode ${\bm u}_1^D=(\sqrt{d_1},...,\sqrt{d_n})/\sqrt{\sum_i d_i}$. Angle shifts of $\delta\bm \varphi$ along ${\bm u}_1^D$ do not modify the synchronous state because ${\bm u}_1\propto{\bm D}^{-1/2}{\bm u}_1^D$. Note also that, by orthogonality with ${\bm u}_1^D$, eigenvector components
must satisfy $\sum_i \sqrt{d_i}u_{\alpha,i}^D=0$ for $\alpha\ge 2$.

\begin{prop}\label{prop1}
\textit{The general solution to Eq.~(\ref{eq:kuramoto_lin3}) reads}
\begin{equation}\label{eq:calpha2}
\begin{aligned}
\delta\varphi_i(t)&=  \sum_\alpha \gamma e^{\frac{-\gamma-\Gamma_{\alpha}}{2}t}\int_0^{t}e^{{\Gamma_{\alpha}}t_1}  \\
&\times \int_{0}^{t_1}[{\bm D}^{-1/2}\delta {\bm{P}}(t_2)]^\top {\bm{u}}_{\alpha}^De^{\frac{\gamma-\Gamma_{\alpha}}{2}t_2}{\rm d}t_2{\rm d}t_1 \;u_{\alpha,i}^D \;,
\end{aligned}
\end{equation}
\textit{with $\Gamma_\alpha=\sqrt{\gamma^2-4\lambda_\alpha^D\gamma}$ where $\lambda_\alpha^D$ is the eigenvalue associated with the eigenvector ${\bm u}_\alpha^D$ of $\mathbb{L}^D$.}
\end{prop}
\begin{IEEEproof}
We first expand angle deviations over
the eigenbasis of $\mathbb{L}^D$ as $\delta\varphi_i(t)=\sum_{\alpha}c_\alpha(t)\, u_{\alpha,i}^D$. From the orthogonality of the
eigenbasis, $({\bm{u}}_\alpha^D)^\top {\bm{u}}_\beta^D=\delta_{\alpha \beta}$,  one straightforwardly rewrites Eq.~(\ref{eq:kuramoto_lin3}) as
\begin{eqnarray}\label{eq:ca}
\gamma^{-1}\ddot{c}_\alpha + \dot{c}_\alpha = ({\bm D}^{-1/2}\delta {\bm P})^\top {\bm{u}}_\alpha^D - \lambda_\alpha^D c_\alpha\;
\end{eqnarray}
for $\alpha=1,...,n$\,. The expansion coefficients $c_\alpha(t)$ can be read from
Eq.~(\ref{eq:calpha2}), and direct differentiation shows that they
solve this equation. This completes the proof.
\end{IEEEproof}

\section{Dynamical Parameters and Transient Excursions}\label{section3}
\subsection{Quantifying frequency excursions}\label{sec:perf}
To evaluate the global response of the system to an external disturbance, we use the following performance metric
\begin{align}\label{eq:c1}
\begin{split}
{\mathcal P}(T) &= T^{-1} \, \int_0^T \,({\bm \omega}^\top-\overline{{\bm \omega}}^\top){\bm D}({\bm \omega}-\overline{{\bm \omega}})\, {\rm d}t   \, ,
\end{split}
\end{align}
because it measures the primary control effort and therefore has a physical meaning~\cite{Poo17}. The quantity
$\overline{{\bm \omega}}^\top=(\dot\Delta,\dot\Delta,...,\dot\Delta)$ with $\dot\Delta (t)  = \sum_i d_i\omega_i(t)/\sum_i d_i$ gives the average frequency deviation over
all nodes in the system.
Because synchronous states are defined modulo any
homogeneous angle shift, the transformation $\theta_i^{(0)} \rightarrow \theta_i^{(0)} + C$ does not change the
synchronous state. Accordingly only phase and frequency shifts with $\sum_i \delta \theta_i(t)=0$ and $\sum_i d_i \omega_i(t)=0$ matter. This is included
in ${\mathcal P}$ by subtracting the average $\overline{{\bm \omega}}$. That ${\mathcal P}(T)$ is a performance metric is easily understood:
low values indicate that the system absorbs the perturbation with little fluctuations, while large values
indicate a large transient excursion around the initial synchronous state.
Using the above change of variables, $\delta\dot{{\bm{\varphi}}}={\bm D}^{1/2}{\bm \omega}=\sum_{\alpha}\dot{c}_\alpha(t)\, {\bm u}_{\alpha}^D$\,,
Eq.~(\ref{eq:c1}) becomes,
\begin{align}\label{eq:perf2}
\begin{split}
\mathcal{P}(T)&=T^{-1} \sum_{\alpha \ge 2}\int_0^T\dot{c}_\alpha^2(t){\rm d}t \; .
\end{split}
\end{align}
This can be calculated using  the explicit expression for $c_\alpha(t)$ from Eq.~(\ref{eq:calpha2}), once a perturbation $\delta {\bm P}(t)$ is given.\\
\begin{prop}\label{prop2}
\textit{Consider a noisy disturbance acting on $N_n$ of the network nodes. The noise ensemble is Gaussian and defined by its
vanishing first moments, $\overline{\delta P_i(t)}=0$, and its second moments $\overline{\delta P_i(t)}\overline{\delta P_j(t')}=\delta_{ij}\delta P_{0i}^2\exp[-|t-t'|/\tau_0]$ with the noise correlation time $\tau_0$. The performance metric $\overline{{\mathcal P}^\infty}$ for primary control effort averaged over this noise ensemble is given by}
\begin{align}\label{eq:P2}
\begin{split}
\overline{\mathcal{P}^\infty}&=\sum_{\alpha\ge 2}\frac{\sum_{i\in N_n} \delta P_{0i}^2{u_{\alpha,i}^D}^2d_i^{-1}}{\lambda_\alpha^D\tau_0+1+\gamma^{-1}\tau_0^{-1}} \; .
\end{split}
\end{align}
\end{prop}
\begin{cor}
\textit{In the limit of short noise correlation time, $\tau_0\ll \gamma^{-1},{\lambda_\alpha^D}^{-1}$ one has}
\begin{align}
\begin{split}\overline{\mathcal{P}^\infty}&=\tau_0 {\sum_{i\in N_n} \delta P_{0i}^2  \, ( 1/m_i - 1 / \sum_j m_j )}\; .\label{eq:P2s}
\end{split}
\end{align}
\end{cor}
\begin{cor}
\textit{In the opposite asymptotic limit, $\tau_0\gg \gamma^{-1},{\lambda_\alpha^D}^{-1}$ one has}
\begin{align}
\begin{split}
\overline{\mathcal{P}^\infty}&=\tau_0^{-1}\sum_{\alpha\ge 2}\frac{\sum_{i\in N_n} \delta P_{0i}^2{u_{\alpha,i}^D}^2d_i^{-1}}{{\lambda_\alpha^D}}  \; .\label{eq:P2l}
\end{split}
\end{align}
\end{cor}
\begin{IEEEproof} Inserting the time derivative of Eq.~\eqref{eq:calpha2} into Eq.~(\ref{eq:perf2}) and taking the average over the noise ensemble
with $\overline{\delta P_i(t)}\overline{\delta P_j(t')}=\delta_{ij}\delta P_{0i}^2\exp[-|t-t'|/\tau_0]$
gives Eq.~\eqref{eq:P2}, with few straightforwardly calculated exponential integrals. The two asymptotic results \eqref{eq:P2s} and \eqref{eq:P2l}
are easily obtained by a Taylor expansion, keeping only the first non-vanishing term. To obtain Eq.~\eqref{eq:P2s},
we also used $\sum_{\alpha\ge 2}{u_{\alpha,i}^D}^2 = \sum_{\alpha\ge 1}{u_{\alpha,i}^D}^2 -  {u_{\alpha,1}^D}^2 =1-d_i/\sum_i d_i$.\\
\end{IEEEproof}

{\textit{Remark 1: The short correlation time asymptotic of Eq.~\eqref{eq:P2s} agrees with the result of~\cite{Poo17} obtained for either single-pulsed or
averaged white-noise perturbations.}}

{\textit{Remark 2: The noise correlators are defined either as time averages, $\overline{\delta P_i(t)}\overline{\delta P_j(t')} = {\rm lim}_{\tau \rightarrow \infty} \tau^{-1}
\int_0^\tau \delta P_i(t+\tau') \delta P_j(t'+\tau') {\rm d}\tau$ or as averages over different noise sequences.
 }}

{\textit{Remark 3: Finite-time corrections to Eqs.~\eqref{eq:P2}--\eqref{eq:P2l} disappear
as ${\cal O}(1/T)$ as $T \rightarrow \infty$.
 }}\\

\begin{prop}\label{prop3}
Under the same assumptions as Proposition~\ref{prop2}, the variance ${\rm var}\left[{\mathcal P}(T)\right]$ of the performance metric  for primary control effort
over the noise ensemble vanishes as $\sim T^{-1} + {\cal O}(1/T^2)$ as $T \rightarrow \infty$.
\end{prop}
The proof proceeds through direct calculation of ${\rm var}\left[{\mathcal P}(T)\right]$. It
is too long to fit in this article and here we only sketch it.
From Eq.~\eqref{eq:perf2} one has
\begin{align}\label{eq:varP}
\begin{split}
{\rm var}\left[{\mathcal P}(T) \right]&=
T^{-2} \sum_{\alpha,\beta \ge 2}\iint_0^T \overline{\dot{c}_\alpha^2(t)\dot{c}_\beta^2(t')}{\rm d}t {\rm d}t' - \overline{\mathcal{P}(T)}^2 \, .
\end{split}
\end{align}
From Eq.~\eqref{eq:calpha2},
each $\dot{c}_{\alpha,\beta}$ contains a noise term $\delta {\bm P}$. The noise average in the first term on the right-hand side of
Eq.~\eqref{eq:varP} therefore consists in pairings of four noise terms.
There are three such contributions. The first one pairs the two $\delta {\bm P}$'s in $\dot{c}_{\alpha}^2(t)$
and the two $\delta {\bm P}$'s in $\dot{c}_{\beta}^2(t)$. This contribution is cancelled by a similar pairing in $\overline{\mathcal{P}(T)}^2$. The other two contributions
pair $\delta {\bm P}$'s across indices $\alpha$ and $\beta$ and accordingly, they constrain the values that $t$ and $t'$ can take with respect to one another,
$|t-t'| \lesssim \tau_0$. Accordingly, the double time integral in Eq.~\eqref{eq:varP} gives a contribution $\sim T \tau_0$, instead of $\sim T^2$, resulting in
${\rm var}\left[{\mathcal P}(T)\right] \sim T^{-1}$.

{\textit{Remark 4: Proposition~\ref{prop3} means that for specific noisy disturbances satisfying the
assumption of Proposition~\ref{prop2} and for long enough observation times $T \gg \tau_0$,
$\mathcal{P}^\infty = \overline{\mathcal{P}^\infty} + {\cal O}(T^{-1/2})$.  The statistical average is therefore representative of specific noise disturbances for
sufficiently long observation time. The validity of Proposition~\ref{prop3} is illustrated numerically in Fig.~\ref{fig2} (c).}}

The two asymptotic limits of large and small $\tau_0$ are particularly interesting as they shed light on the influence of dynamical parameters,
in particular on the interplay between local disturbance absorption by inertia and long-range propagation through low-lying network modes.
First, in the short correlation time limit given by Eq.~(\ref{eq:P2s}),  $\overline{\mathcal{P}^\infty}$ explicitly depends on inertia but not on the
coupling network. This reflects the fact that, in the white-noise limit, the perturbation remains local and is easily absorbed, if
there is enough inertia. Second,
Eq.~(\ref{eq:P2l}) shows that, in the long correlation time limit, $\overline{\mathcal{P}^\infty}$ does not depend on
inertia. This suggests that changing inertia in any direction
will not change $\overline{\mathcal{P}^\infty}$ in the limit of long noise correlation time. This is a conjecture since Eq.~(\ref{eq:P2l}) has been derived
under the assumption of constant damping-to-inertia ratio, $\gamma=d_i/m_i$. Below we numerically confirm this conjecture.
Simultaneously, Eq.~(\ref{eq:P2l}) also shows that, in the long correlation time limit, $\overline{\mathcal{P}^\infty}$
is determined by the structure of the coupling network, with the modes with smallest eigenvalues having the largest influence. Those modes
are extended over the whole network in large power grids, as is illustrated in Fig.~\ref{fig3} (b). Accordingly, in the limit of long noise correlation time,
the disturbance is able to propagate over large distances in the network, and inertia has little influence on this large-scale propagation.
\begin{figure*}
\centering
\includegraphics[width=0.9\textwidth]{fig1.pdf}
\caption{Comparison between numerical calculations obtained by time-evolving Eq.~(\ref{eq:kuramoto}) and the theoretical result of Eq.~(\ref{eq:P2l}) for the primary control effort,
for short correlation time $\gamma\tau_0 = 4\times 10^{-3}$ (a) and long correlation time $\gamma\tau_0 = 40$ (b)
in the IEEE 118-Bus test case. Averages are made over $10$ noisy sequences and standard deviations are shown by barely visible vertical line. Blue crosses correspond to constant damping-to-inertia ratio, $d_i/m_i=\gamma=0.4s^{-1}$, while orange symbols correspond to inhomogeneous,
non-vanishing inertia on generator nodes (squares) and inertialess consumer nodes (crosses). On generation nodes one has $\gamma_i\tau_0\in [30 , 60]$\,. (c) Ratio between standard deviation and average of the primary control effort for the IEEE 118-Bus test case. For long enough $T$, the ratio scales as $T^{-1/2}$, confirming Prop.~\ref{prop3}. Averages are made over $40$ different noise sequences.}\label{fig2}
\end{figure*}

{\textit{Remark 5: It is important to realize that
fluctuations of renewable energy sources occur on time scales that are large compared to the intrinsic time scales of the power system~\cite{Mei11}. As a matter of fact, time scales in the synchronous grid of continental Europe have been found to
satisfy $\gamma^{-1} \simeq 2.5 s$ and  ${\lambda_\alpha^D}^{-1} \lesssim 0.5 s$~\cite{Tyl18b}.
Provided the above conjecture is corroborated, Eqs.~\eqref{eq:P2}--\eqref{eq:P2l} suggest that fluctuations from new renewables excite network modes and
are efficiently absorbed by optimizing the distribution of damping with little regard for inertia. This conjecture is numerically confirmed below.}

{\textit{Remark 6: Similar conclusions as in Eqs.~\eqref{eq:P2}--\eqref{eq:P2l} regarding local inertia absorption vs. large-scale mode propagation
are obtained in the case of step disturbances corresponding to sudden power losses, as a function of their duration $\tau_0$.}

{\textit{Remark 7: The lossless line approximation used in this paper does not account for ohmic dissipation. We expect that the latter enhances mode damping
and accordingly undermines disturbance propagation in the case of noise with long correlation time, but that it affects only marginally our result for short correlation time.
Investigations beyond the lossless line approximation would be very welcome but lie beyond the scope of the present paper.}

\section{Numerical Simulations}\label{section5}
\subsection{Dynamical parameters for simulations}\label{section5a}
We consider two different cases: (i) cases with homogeneous damping-to-inertia ratio $d_i/m_i=\gamma$; (ii) realistic heterogeneous cases, using the nonlinear swing dynamics of Eq.~(\ref{eq:kuramoto}).
In the homogeneous case, we used $d_i = \alpha |P_i^{(0)}|/\omega_0$, where $P_i^{(0)}$ is the produced/consumed power at nominal frequency $\omega_0$,
and set $m_i=\gamma^{-1}d_i$.
In the heterogeneous case, the inertia parameter vanishes on consumer nodes and is $m_i = 2H_i |P_i^{(0)}|/\omega_0$ on generator nodes,
where $H_i$ depends on the type of generator~\cite{Mac08}.
Damping is given by Eq.~(5.24) and Table 4.3 in Ref.~\cite{Mac08} for generators and by
$d_i = \alpha |P_i^{(0)}|/\omega_0$ for consumer nodes. In all
cases we use $\omega_0=2\pi\times 50Hz$. For the IEEE-118 Bus test case discussed in Section \ref{sec:ieee118}, $\alpha=1.5$ and $H_i=5s$. For the PanTaGruEl European model discussed in Section \ref{sec:panta}), $\alpha=1.5$ and $H_i$ varies according to the generator type as described in Ref.~\cite{Pag19c}.

\subsection{IEEE 118-Bus test case}\label{sec:ieee118}
The main prediction from Eqs.~\eqref{eq:P2}-\eqref{eq:P2l} is that, for noise correlation $\tau_0$ that is longer than the
other characteristic time scales in the system, the performance metric $\overline{\mathcal{P}^\infty}$ does not depend on  inertia parameters. This was conjectured from
Eq.~\eqref{eq:P2l} where inertia does not appear.
Fig.~\ref{fig2} shows performance metric obtained from individual noisy disturbance of a single node, repeating the operation for all nodes in the IEEE 118-Bus test case.
As disturbance, we take Gaussian noise with the same first two moments as in Proposition~\ref{prop2}.
First, the primary control effort $\overline{\mathcal{P}^\infty}$ is calculated for a constant damping-to-inertia ratio (blue crosses), then for a distribution of dynamical parameters where generator nodes have inertia while consumer nodes do not (orange squares and crosses). As predicted by Eq.~(\ref{eq:P2l}), for long correlation time of the noise [panel (b), $\gamma \tau_0=40$],
the two distributions of dynamical parameters give the same $\overline{\mathcal{P}^\infty}$, corroborating our conjecture that it does not depend on inertia. We have found (but do not show) that
the performance metric only depends on the damping distribution in that case.

For short noise correlation time, on the other hand,
Eq.~(\ref{eq:P2s}) explicitly depends on the damping-to-inertia ratio $d_i/m_i=\gamma$\,, and we expect that the numerical data will differ from the theoretical
prediction once this ratio is no longer constant. This is confirmed in
Fig.~\ref{fig2} [panel (a), $\gamma \tau_0=4\times 10^{-3}$] where with $d_i/m_i=\gamma$, numerical data fall on the theory (blue crosses). However, once
$d_i/m_i$ is no longer constant, numerical data and theoretical prediction differ significantly (orange symbols).
Quite interestingly, we found that for noisy perturbations on generator nodes
with inertia, the theory still gives a remarkably accurate estimate for the primary control effort $\mathcal{P}^\infty$. An understanding of this remarkable agreement would be highly welcome, particularly since it could justify dynamic performance analysis on Kron reduced networks.

\subsection{Time scales in high-voltage electric power grids}

We have shown that the primary control effort against fluctuating disturbances in the form of colored noise
behaves very differently depending on the position of the noise correlation time relative to the characteristic time scales in the system.
Furthermore the primary control effort is captured by our theory even for inhomogeneous dynamical parameters, when the noise correlation time is long enough.
It is therefore desirable to identify
what regime applies to a realistic high-voltage power grid subjected to fluctuating sources of power, in particular those generated by
new renewable sources of energy.
To that end, we consider in the next paragraph a realistic model of the synchronous grid of continental Europe~\cite{Tyl18b,Pag19c},
with the following time scales
\begin{subequations}\label{eq:time}
\begin{align}
{\lambda_\alpha^D}^{-1}&< 0.5s\;, \; {\text{for $\alpha=2,...,n$}\,,}\\
\gamma^{-1}&=\frac{\langle m_i \rangle}{\langle d_i \rangle}= 2.5s\, ,
\end{align}
\end{subequations}
where $\langle \;\ldots\; \rangle$ means that we take the average over all nodes in the grid. It is commonly accepted that
power fluctuations from renewable energy sources such as wind turbines or photovoltaic panels
fluctuate on time scales that are larger than both time scales in Eq.~(\ref{eq:time})~\cite{Mei11}. Therefore, the asymptotic limit
of large noise correlation time, corresponding to Eq.~(\ref{eq:P2l}) applies, and we expect that the primary control effort as measured by
Eq.~\eqref{eq:c1} is influenced only by damping, and not by inertia. Numerical results to be presented in the next paragraph corroborate this expectation.

As a side-remark we note that when the perturbation corresponds to the sudden disconnection of a power generator,
controls usually try to reconnect the bus several times quickly after the fault (typically within few AC cycles).
The typical time scale for such a perturbation is then less than several if not all system's time scales, and inertia obviously matters to absorb such sudden faults, as it is predicted by Eq.~(\ref{eq:P2s}).

\subsection{The PanTaGruEl European model}\label{sec:panta}
To further illustrate the influence of dynamical parameters on the primary control effort,
we numerically compute Eq.~(\ref{eq:c1}) on the large-scale PanTaGruEl model of the European high-voltage transmission grid~\cite{Pag19c}.
The model is shown on Fig.~\ref{fig3} (b). It has 3809 nodes and 4944 lines. More details can be found in Ref.~\cite{Pag19c}.
We considered three different cases, (i) a homogeneous situation that corresponds to today's grid in terms of its global amount of inertia, with $\gamma\tau_0=4$ and $\gamma=d_i/m_i$ constant,
(ii) a homogeneous situation where the inertia is reduced by a factor $10$\,, i.e. $\gamma\tau_0=40$, and $\gamma=d_i/m_i$ constant,
and (iii) a realistic situation with inhomogeneous damping parameters and where inertia vanishes on consumer nodes and is inhomogeneous on production nodes
as described in Section~\ref{section5a} and Ref.~\cite{Pag19c}, with
$\tau_0$ larger than all other time scales in the system.
Remarkably, Fig.~\ref{fig3} (a) shows that the primary control effort for all cases is well predicted by Eq.~(\ref{eq:P2l}). This confirms our main finding
that, for fluctuations with a correlation time longer than any other characteristic  time scale in the system, inertia does not affect the primary control effort, Eq.~(\ref{eq:P2l}).
Quite surprisingly, an overall reduction of the total available inertia by a factor of 10 does not affect the primary control
effort, Eq.~\eqref{eq:c1}.
\begin{figure*}
\centering
\includegraphics[width=0.95\textwidth]{fig2.pdf}
\caption{(a) Comparison between numerical calculations obtained by time-evolving Eq.~(\ref{eq:kuramoto}) and the theoretical result of Eq.~(\ref{eq:P2l}) for the PanTaGruEl model of the synchronous grid
of continental Europe~\cite{Pag19c} shown on panel (b).
Three cases are considered. The first two are grids with $\gamma=d_i/m_i$ constant,
with today's average inertia (orange crosses) and an inertia reduced by a factor of 10 (blue crosses). The third case corresponds to a realistic situations
as discussed in the text, with  $\gamma_i\tau_0\in [ 20,1600 ]$ on generation nodes.
In all cases, $\tau_0$ is the longest time scale,
consequently, the inertia-independent theoretical prediction of Eq.~\eqref{eq:P2l} accurately captures all numerical data. Averages are made over $10$ noisy sequences. (b) Network eigenmodes of $\mathbb{L}^D$
with the first two non-vanishing eigenvalues. These slow modes are extended over the whole network, with higher amplitudes on peripheral nodes.
Disturbances on the four buses highlighted in light blue correspond to the smallest primary control effort. These buses
lay in the center of the network where the slow modes have small amplitudes.
Disturbances on the four buses in dark red, on the other hand, have largest primary control effort. They are located at the periphery of the network
where the slow modes have large amplitudes~\cite{Tyl18b}. This shows that large primary control effort for noise with large correlation time correspond
to excitations of slow network modes, which in their turn propagate the disturbance over large distances in the network.}\label{fig3}
\end{figure*}

\section{Conclusion}\label{section6}
With the ongoing energy transition resulting in strongly increased penetrations of new renewable sources of electrical energy, a question of crucial importance is how
grid stability will evolve, given the resulting reduction of globally available rotational inertia and enhanced power fluctuations. We have shown that reduced inertia may pose
problems only for perturbation occurring/fluctuating on very short time scales, shorter than all other characteristic time scales in the system. In continental-size transmission grids,
these time scales are shorter than few seconds, consequently, power fluctuations from new renewables will not affect grid stability per se.

Inertia is of course important to absorb
sudden faults occurring on very short time scales such as line faults or disconnection/reconnection of large power plants and so forth. Simultaneously, our result
of Eq.~\eqref{eq:P2s} indicates that the resulting primary control effort is independent of the grid topology. Accordingly,
optimal inertia distribution needs to follow the distribution of potential faults, for instance being larger in regions with higher density of generators. A similar conclusion was
drawn in Refs.~\cite{Poo17} and \cite{Pag19b}.


\section{Introduction}

In this paper we will give new proofs for the following theorems of Schrijver. The first theorem is about perfect matchings of regular bipartite graphs.

\begin{Th}[Schrijver \cite{Sch2}] \label{matching}
Let $G=(A,B,E)$ be $d$--regular bipartite graph on $2n$ vertices. Let $\prm(G)$ denote the number of perfect matchings of $G$. Then
$$\prm(G)\geq \bigg( \frac{(d-1)^{d-1}}{d^{d-2}} \bigg)^{n}.$$
\end{Th}

The next theorem is about Eulerian orientations of regular graphs. Recall that an orientation of a graph $G$ is Eulerian if the in-degree and out-degree are equal at each vertex. (In particular, the degree of a vertex must be even.)

\begin{Th}[Schrijver \cite{Sch1}] \label{orientation}
Let $G$ be a graph, where the degree of the vertex $v$ is $d_v$. Suppose that $d_v$ is even for each $v$. Let $\e(G)$ denote the number of Eulerian orientations of the graph $G$. Then
$$\e(G)\geq \prod_{v\in V(G)}\frac{\binom{d_v}{d_v/2}}{2^{d_v/2}}.$$
\end{Th}

The main goal of this paper is to show that these two theorems have a common generalization. To spell out this generalization we will count the number of orientations in a graph with prescribed in-degree sequence.

\begin{Def}
Let $\underline{r}=(r_v)_{v\in V(G)}\in \mathbb{Z}^{V(G)}$. Let $\e_{\underline{r}}(G)$ denote the number of those orientations of the graph $G$, where the in-degree of the vertex $v$ is $r_v$.
\end{Def}

Observe that if we have a $d$--regular bipartite graph $G=(A,B,E)$, then the perfect matchings are in bijection with those orientations of the graph, where the in-degree of the vertices in $A$ is $1$, and is $d-1$ in case of the vertices of $B$. Indeed, simply orient each edge of a perfect matching towards $A$, and every other edge towards $B$. Clearly, if we have such an orientation, then the edges oriented towards $A$ form a perfect matching.
\bigskip

The following theorem might look technical, but it easily implies both  Theorem~\ref{matching} and \ref{orientation}.

\begin{Th} \label{general}
Let $G=(V,E)$ be a graph with degree $d_v$ at vertex $v$, and let
$\e_{\underline{r}}(G)$ denote the number of orientations of the graph $G$, where the in-degree of the vertex $v$ is $r_v$.  Then
$$\e_{\underline{r}}(G)\geq \prod_{v\in V(G)}\binom{d_v}{r_v}\bigg(\frac{r_v}{d_v}\bigg)^{r_v}\bigg(\frac{d_v-r_v}{d_v}\bigg)^{d_v-r_v}\cdot \inf_{x_u>0}\frac{\prod_{(u,v)\in E(G)}(x_u+x_v)}{\prod_{u\in V(G)}x_u^{r_u}}.$$
\end{Th}

We remark that the role of the multivariate polynomial $P_G(\underline{x}):=\prod_{(u,v)\in E(G)}(x_u+x_v)$ in the theorem comes from the fact that
$$\prod_{(u,v)\in E(G)}(x_u+x_v)=\sum_{\underline{r}}\e_{\underline{r}}(G)\prod_{u\in V(G)}x_u^{r_u}.$$
In other words, this theorem is about how to give a lower bound on a coefficient of a multivariate polynomial in terms of the polynomial. It turns out that this lower bound is possible, because the polynomial $P_G(\underline{x})$ is a real stable polynomial. The definition of real stability is the following.

\begin{Def}
A multivariate polynomial $P(x_1,\dots ,x_n)$ with complex coefficients is stable if $P(z_1,\dots ,z_n)\neq 0$ whenever $\mathrm{Im}(z_i)>0$ for $i=1,\dots ,n$. A polynomial is called real stable if it is stable and its coefficients are real.
\end{Def}

Note that a univariate polynomial with real coefficients is  stable if and only if it is real-rooted. So real stability is a generalization of real-rootedness for multivariate polynomials.

We remark that Gurvits \cite{Gur1} already gave a  proof of Theorem~\ref{matching} using real stable polynomials. He used the polynomial
$$Q(\underline{x})=\prod_{v\in B}\bigg(\sum_{u\in A}x_u\bigg)$$
in his proof. In this case the coefficient of $\prod_{u\in A}x_u$ is exactly $\prm(G)$. In fact, we will follow exactly the strategy of Gurvits.  This strategy is based on two concepts, the real stability and the capacity of a polynomial. The latter was invented by Gurvits himself. We will review these concepts in the next section.
\bigskip

\noindent \textbf{What is new in this paper?} The proof
of Theorem~\ref{coefficient theorem} is new, the theorem itself appeared in \cite{Gur3} in a slightly different form.
The use of the polynomial $P_G(\underline{x})$ in these proofs also seems to be new, although variants of this polynomial appeared in \cite{Gur3}, but never exactly this one. Proving Theorem~\ref{orientation} via stable polynomials is also new. On the other hand, the general strategy is not new at all. In fact, one of our main goals is to advertise this theory, so this paper can be considered as a mini survey.
\bigskip

\noindent \textbf{This paper is organized as follows.} In the next section we collect the basic facts about real stable polynomials and capacity. In Section 3 we prove Theorem~\ref{general} and derive Theorem~\ref{matching} and \ref{orientation} from it. In Section 4 we collected pointers to the literature.

\section{Real stability and capacity}

In this section we review the basic properties of real stability and capacity. To keep this paper  self-contained we will prove every result that we use apart from the Hermite-Sylvester criterion.

\subsection{Stability} Recall that a multivariate polynomial $P(x_1,\dots ,x_n)$ is real stable if it has real coefficients, and $P(z_1,\dots ,z_n)\neq 0$ whenever $\mathrm{Im}(z_i)>0$ for $i=1,\dots ,n$. We have seen that real stability is a generalization of univariate real-rooted polynomials. As the next lemma shows there is a more direct connection between the two concepts.

\begin{Lemma} \label{stable-real-rooted} A multivariate polynomial $P(z_1,\dots ,z_n)\in \R[z_1,\dots ,z_n]$ is stable if and only if for all $\underline{v}=(v_1,\dots, v_n)\in \R^n$ and $\underline{u}=(u_1,\dots ,u_n)\in \R^n_{>0}$, the univariate polynomial $g(t)=P(v_1+tu_1,\dots ,v_n+tu_n)$ is real-rooted.
\end{Lemma}

\begin{proof} Let $\mathbb{H}=\{z \ |\ \mathrm{Im}(z)>0\}$. First suppose that $g(t)$ is not real rooted, then it has a root $a+bi\in \C$, where $b\neq 0$. Since the coefficients of $g$ are real, $a\pm bi$ are both zeros of $g$. So we can assume that $b>0$. But then $\mathrm{Im}(v_j+u_j(a+bi))=bu_j>0$ shows that the numbers $z_j=v_j+u_j(a+bi) \in \HH$ for $j=1,\dots ,n$ and $P(z_1,\dots ,z_n)=0$, thus $P$ is not stable.

Next suppose that $P$ is not stable, it has some zero $(z_1,\dots ,z_n)\in \HH^n$. Let $z_j=a_j+b_ji$. Then $b_j>0$. Then the polynomial
$g(t)=P(a_1+tb_1,\dots ,a_n+tb_n)$ satisfies that $(a_1,\dots ,a_n)\in \R^n$ and $(b_1,\dots ,b_n)\in \R^n_{>0}$ and has a non-real zero, namely $i$.
\end{proof}

We will also need the following lemma.

\begin{Lemma} \label{real-rooted-convergence}
Given a polynomial $P(x)$ and suppose that $(P_n(x))_n$ is a sequence of real-rooted  polynomials such that
$\lim_{n\to \infty}P_n(x)=P(x)$ coefficientwise. Then $P(x)$ is real-rooted.
\end{Lemma}

Lemma~\ref{real-rooted-convergence} is a simple consequence of the well-known Hermite-Sylvester criterion. For a short, simple proof see that paper \cite{Nat}.

\begin{Lemma}[Hermite-Sylvester criterion] \label{HS-criterion}
Let $P(x)$ be a non-constant polynomial of degree $d$ with real coefficients and $\lambda_1,\dots ,\lambda_d$ not necessarily distinct zeros. Let $m_k=\sum_{j=1}^d\lambda_j^k$. Then $P(x)$ is real-rooted if and only if the $d\times d$ matrix $(m_{i+j-2})_{i,j=1,\dots ,d}$ is positive semi-definite.
\end{Lemma}

Note that $m_k$ can be computed from the coefficients by the Newton-Girard identities, so they are continuous functions of the coefficients. Hence the  Hermite-Sylvester criterion implies Lemma~\ref{real-rooted-convergence}.

Next we collect some operations that preserve stability. A general theory of stability preserver operations is developed by Borcea and Br\"and\'en \cite{BoBr}. We also recommend the paper of Choe, Oxley, Sokal and Wagner \cite{COSW} for a comprehensive list of operations that preserves stability.

\begin{Th}
Let $P(x_1,\dots ,x_n)$ be a real stable polynomial. Suppose that the degree of $x_1$ in $P$ is $d$. Then the following hold true.\\
(a) The polynomial $x_1^dP(-1/x_1,x_2,\dots x_n)$ is real stable.\\
(b) If $a\in \mathbb{R}$, then $P(a,x_2,\dots ,x_n)$ is real stable or the constant $0$ polynomial.\\
(c) The polynomial $\frac{\partial}{\partial x_1}P$ is real stable or the constant $0$ polynomial.\\
\end{Th}

\begin{proof}
As before let $\mathbb{H}=\{z \ |\ \mathrm{Im}(z)>0\}$. Then the first claim is trivial since $z\mapsto \frac{-1}{z}$ maps $\mathbb{H}$ to $\mathbb{H}$.
\medskip

Next we prove part (b). By Lemma~\ref{stable-real-rooted} the polynomial $P(a+\varepsilon t,v_2+tu_2,\dots ,v_n+tu_n)$ is real-rooted for every $\varepsilon,u_2,\dots u_n\in \R_{>0}$ and $a,v_2,\dots ,v_n\in \R$. Let $\varepsilon \to 0$, then by Lemma~\ref{real-rooted-convergence} we get that $P(a,v_2+tu_2,\dots ,v_n+tu_n)$ is real-rooted or the constant $0$ function for every $u_2,\dots ,u_n\in \R_{>0}$ and $v_2,\dots ,v_n$. Now using the  other direction of Lemma~\ref{stable-real-rooted} we get that $P(a,x_2,\dots ,x_n)$ is a real stable polynomial.

Next we prove part (c). Let
$$P(x_1,\dots ,x_n)=\sum_{k=0}^dP_k(x_2,\dots ,x_n)x_1^k.$$
First we show that $P_d(x_2,\dots ,x_n)$ is a real stable polynomial. By part (a)
$$R(x_1,\dots ,x_n):=x_1^dP(-1/x_1,x_2,\dots x_n)=\sum_{k=0}^dP_k(x_2,\dots ,x_n)(-1)^kx_1^{d-k}$$
is real stable. Then $R(0,x_2,\dots ,x_n)=(-1)^dP_d(x_2,\dots ,x_n)$ is real stable, and so $P_d(x_2,\dots ,x_n)$ is real stable.

Now let $Q=\frac{\partial}{\partial x_1}P$, and let $\underline{a}=(a_1,a_2,\dots, a_n)\in \mathbb{H}^n$. We show that if $Q\not\equiv 0$, then $\mathrm{Im}\left(\frac{Q(\underline{a})}{P(\underline{a})}\right)<0$. Note that here we use that $P$ is stable, so $P(\underline{a})\neq 0$, and we can divide with it. If $d=0$ then $Q\equiv 0$. We can assume that $d\geq 1$. Let
$$g(x):=P(x,a_2,\dots ,a_n)=\sum_{k=0}^dP_k(a_2,\dots ,a_n)x^k.$$
Note that $P_d(a_2,\dots ,a_n)\neq 0$ since $P_d$ is real-stable. So $g(x)$ has degree $d\geq 1$. Then
$g(x)=c\prod_{i=1}^d(x-\rho_i)$, and we have
$$\frac{g'(x)}{g(x)}=\sum_{i=1}^d\frac{1}{x-\rho_i}.$$
Note that $\mathrm{Im}(\rho_i)\leq 0$, otherwise $P(\rho_i,a_2,\dots ,a_n)=0$ would yield a zero in $\mathbb{H}^n$. Hence
$$\mathrm{Im}\left(\frac{Q(\underline{a})}{P(\underline{a})}\right)=\mathrm{Im} \left(\frac{g'(a_1)}{g(a_1)}\right)=\mathrm{Im} \left(
\sum_{i=1}^d\frac{1}{a_1-\rho_i}\right)<0.$$
In particular, this shows that $Q(\underline{a})\neq 0$. Hence $Q$ is stable. (Remark: we essentially repeated the proof of Gauss--Lucas theorem that asserts that the zeros of the derivative of a polynomial lie in the convex hull of the zeros of the polynomial.)
\medskip

\end{proof}

\subsection{Capacity of a polynomial} In this section we introduce the concept capacity.

\begin{Def}[Gurvits \cite{Gur1}]
Let $P(x_1,\dots ,x_n)$ be a multivariate polynomial with non-negative coefficients. Let $\underline{\alpha}=(\alpha_1,\dots ,\alpha_n)$ be a non-negative vector. Then the $\underline{\alpha}$-capacity of the polynomial $P(x_1,\dots ,x_n)$ is
$$\mathrm{cap}_{\underline{\alpha}}(P)=\inf_{x_1,\dots ,x_n>0}\frac{P(x_1,\dots ,x_n)}{\prod_{i=1}^nx_i^{\alpha_i}}.$$
\end{Def}

Note that if the numbers $\alpha_i$ are integers, then the capacity is an upper bound for the coefficient of the term $\prod_{i=1}^nx_i^{\alpha_i}$.

So in Theorem~\ref{general} we had $\mathrm{cap}_{\underline{r}}(P_G)$ in the statement.
In the definition of capacity we never used that $P$ is real stable, but it turns out that this concept is especially useful when we study stable polynomials. The main reason for this phenomenon is that one can often govern the capacity for stability preserver operators. A general theory of capacity preserver linear operators was developed in the  paper of Leake and Gurvits \cite{LeGu}.   The following theorem is a special case of their theory.

\begin{Th} \label{capacity preserving}
Let $P(x_1,\dots ,x_n)$ be a real stable polynomial with non-negative coefficients. Suppose that the degree of $x_1$ in $P$ is $d$. Let
$$Q=\frac{1}{r!}\bigg( \frac{\partial^r}{\partial x_1^r}P\bigg)\bigg|_{x_1=0}.$$
In other words, if we expand $P$ as a polynomial of $x_1$, then $Q$ is the coefficient of $x_1^r$. \\
Let $\underline{\alpha}=(\alpha_1,\dots ,\alpha_n)$, where $\alpha_1=r$, and $\underline{\alpha}'=(\alpha_2,\dots ,\alpha_n)$. Then
$$\mathrm{cap}_{\underline{\alpha}'}(Q)\geq \binom{d}{r}\bigg(\frac{r}{d}\bigg)^{r}\bigg(\frac{d-r}{d}\bigg)^{d-r}\mathrm{cap}_{\underline{\alpha}}(P).$$
\end{Th}

An immediate corollary of Theorem~\ref{capacity preserving} is the following theorem.

\begin{Th}[Coefficient lemma for  stable polynomials with non-negative coefficients] \label{coefficient theorem}
Let $P(x_1,\dots ,x_n)$ be a real stable polynomial with non-negative coefficients. Suppose that the degree of $x_i$ in $P$ is at most $d_i$ for $i=1,\dots ,n$. Let $\underline{r}=(r_1,\dots ,r_n)$ and $a_{\underline{r}}$ be the coefficient of $\prod_{i=1}^nx_i^{r_i}$ in $P$. Then
$$a_{\underline{r}}\geq \prod_{i=1}^n\binom{d_i}{r_i}\bigg(\frac{r_i}{d_i}\bigg)^{r_i}\bigg(\frac{d_i-r_i}{d_i}\bigg)^{d_i-r_i}\mathrm{cap}_{\underline{r}}(P).$$
\end{Th}

If we apply Theorem~\ref{coefficient theorem} to the real stable polynomial  $P_G(\underline{x})=\prod_{(u,v)\in E(G)}(x_u+x_v)$ we get Theorem~\ref{orientation}.

The key lemma to prove the above theorems is the following.

\begin{Lemma} \label{main lemma}
Let $p(z)=\sum_{k=0}^da_kz^k$ be a real-rooted polynomial with non-negative coefficients. Then
$$a_r\geq \binom{d}{r}\bigg(\frac{r}{d}\bigg)^r\bigg(\frac{d-r}{d}\bigg)^{d-r}\inf_{t>0}\frac{p(t)}{t^r}.$$
\end{Lemma}

We remark that this lemma for $r=1$ was the main ingredient of the proof of Gurvits for Theorem~\ref{matching}.

We will derive Lemma~\ref{main lemma} from the following statement.

\begin{Lemma} \label{main lemma2}
Let $p(z)=\sum_{k=0}^da_kz^k$ be a real-rooted polynomial with non-negative coefficients. Suppose that $p(1)=1$ and $p'(1)=r$ is an integer. Then
$$a_r\geq \binom{d}{r}\bigg(\frac{r}{d}\bigg)^r\bigg(\frac{d-r}{d}\bigg)^{d-r}.$$
\end{Lemma}

We remark that Lemma~\ref{main lemma2} is a special case of the following theorem of Hoeffding \cite{Hoe}. Nevertheless to keep our paper self-contained we will give a proof of Lemma~\ref{main lemma2}. The intuitive meaning of Hoeffding's theorem is that among probability distributions coming from real-rooted polynomials and fixed expected value the binomial distribution is the least concentrated around its expected value.

\begin{Th}[Hoeffding \cite{Hoe}] \label{Hoeffding} Let $p(z)=\sum_{k=0}^dp_kz^k$ be a real-rooted polynomial with $p_k\geq 0$, and $p(1)=1$, that is, $\sum_{k=0}^dp_k=1$. Let $s$ be defined by the equation $\sum_{k=0}^dkp_k=ds$.
Suppose that for non-negative integers $b$ and $c$ we have $b\leq ds\leq c$. Then
$$\sum_{k=b}^cp_k\geq \sum_{k=b}^c\binom{d}{k}s^k(1-s)^{d-k}.$$
\end{Th}

 It is  easy to see that Lemma~\ref{main lemma2} is a special case of Lemma~\ref{main lemma}, but as the following proof shows they are actually equivalent statements.

\begin{proof}[Proof of Lemma~\ref{main lemma} from Lemma~\ref{main lemma2}.]
Suppose that $p(z)=\sum_{k=m}^Ma_kz^k$, where $a_m,a_M>0$.
If $r<m$ or $r>M$, then $\inf_{t>0} \frac{p(t)}{t^r}=0$ so the claim is true in this case. If $r=m$, then $\inf_{t>0} \frac{p(t)}{t^r}=a_m$ so the claim is again true. If $r=M$, then  $\inf_{t>0} \frac{p(t)}{t^r}=a_M$ so we are again done.
Thus we can assume that $m<r<M$.
Observe that $\frac{tp'(t)}{p(t)}$ is monotone increasing, $\frac{tp'(t)}{p(t)}\big|_{t=0}=m$ and $\lim_{t\to \infty} \frac{tp'(t)}{p(t)}=M$. So we can choose $t_r>0$ in such a way that $\frac{t_rp'(t_r)}{p(t_r)}=r$. Let us consider the probability distribution $q_j=\frac{a_jt_r^j}{p(t_r)}$. Then $\sum_k kq_k=r$, and $\sum_{j=0}^dq_jz^j$ is still a real-rooted polynomial.
Next let us apply Lemma~\ref{main lemma2}. Then
$$\frac{a_rt_r^r}{p(t_r)}=q_r\geq \binom{d}{r}\left(\frac{r}{d}\right)^r\left(\frac{d-r}{d}\right)^{d-r}.$$
In other words,
$$a_r\geq \binom{d}{r}\left(\frac{r}{d}\right)^r\left(\frac{d-r}{d}\right)^{d-r}\frac{p(t_r)}{t_r^r}\geq \binom{d}{r}\left(\frac{r}{d}\right)^r\left(\frac{d-r}{d}\right)^{d-r}\inf_{t>0}\frac{p(t)}{t^r}.$$
\end{proof}

Before we prove Lemma~\ref{main lemma2} we show that Lemma~\ref{main lemma} indeed implies Thorem~\ref{capacity preserving}.

\begin{proof}[Proof of Theorem~\ref{capacity preserving}]
For fixed $a_2,\dots ,a_n$ consider the polynomial
$$g(x)=P(x,a_2,\dots ,a_{n}).$$
Clearly, $\frac{1}{r!}\frac{d^r}{dx^r}g(x)\Big|_{x=0}=Q(a_2,\dots ,a_{n})$. Since $P$ is stable, and we substituted $a_i\in \mathbb{R}$ into it, $g(x)$ is stable. In other words, it is real-rooted. Thus we can use Theorem~\ref{main lemma}:
$$a_r \geq  \binom{d}{r}\bigg(\frac{r}{d}\bigg)^r\bigg(\frac{d-r}{d}\Bigg)^{d-r}\inf_{x>0}\frac{g(x)}{x^r}.$$
Hence we have
\begin{align*}
\frac{Q(a_2,\dots ,a_{n})}{\prod_{i=2}^{n}a_i^{\alpha_i}}=\frac{\frac{1}{r!}\frac{d^r}{dx^r}g(x)\Big|_{x=0}}{\prod_{i=2}^{n}a_i^{\alpha_i}}&\geq \binom{d}{r}\bigg(\frac{r}{d}\bigg)^r\bigg(\frac{d-r}{d}\bigg)^{d-r}\frac{1}{\prod_{i=2}^{n}a_i^{\alpha_i}}\inf_{x>0}\frac{g(x)}{x^r}\\
                            &= \binom{d}{r}\bigg(\frac{r}{d}\bigg)^r\bigg(\frac{d-r}{d}\bigg)^{d-r}\inf_{x>0} \frac{P(x,a_2,\dots ,a_{n})}{x^r\cdot \prod_{i=2}^{n} a_i^{\alpha_i}}\\
							 &\geq \binom{d}{r}\bigg(\frac{r}{d}\bigg)^r\bigg(\frac{d-r}{d}\bigg)^{d-r}\mathrm{cap}_{\underline{\alpha}}(P).
\end{align*}
Taking infimum on the left side we get that
$\mathrm{cap}_{\underline{\alpha}'}(Q)\geq \binom{d}{r}\big(\frac{r}{d}\big)^r\big(\frac{d-r}{d}\big)^{d-r}\mathrm{cap}_{\underline{\alpha}}(P)$.
\end{proof}

\subsection{Proof of Lemma~\ref{main lemma2}}
In this section we prove Lemma~\ref{main lemma2}. The condition on the non-negativity of the coefficients and $p(1)=1$ implies that $p(z)$ can be written as follows:
$$p(z)=\prod_{i=1}^d(1-\alpha_i+\alpha_iz),$$
where $0\leq \alpha_i\leq 1$. Indeed, as all coefficients are non-negative, there can be no positive roots, thus using $p(1)=1$ the polynomial can be rewritten in the following way:
$$p(z)=a_n\prod_{i=1}^d(z+\rho_i)=\prod_{i=1}^d\left(\frac{z+\rho_i}{1+\rho_i}\right)=\prod_{i=1}^d(1-\alpha_i+\alpha_iz),$$
where $\alpha_i=\frac{1}{1+\rho_i}$.
Note that
$$r=p'(1)=\sum_{i=1}^d \alpha_i\prod_{j\not=i}(\alpha_j+1-\alpha_j)=\sum_{i=1}^d \alpha_i.$$
Consider the domain
$$D_{d,r}=\left\{(\alpha_1,\dots ,\alpha_d)\in \mathbb{R}^d\ \Big|\ 0\leq \alpha_i\leq 1\ (i=1,\dots ,d),\ \sum_{i=1}^d\alpha_i=r\right\}.$$
Clearly, the coefficient $a_r$ of $z^r$ in $p(z)$ can be expressed as
$$a_{r}=\sum_{\substack{K\subset[n]\\ |K|=r}}\prod_{j\in K}\alpha_j\prod_{j\not\in K}(1-\alpha_j).$$
So let us introduce the function
$$f_{d,r}(x_1,\ldots,x_d)=\displaystyle\sum_{\substack{K\subset[n]\\ |K|=r}}\prod_{j\in K}x_j\prod_{j\not\in K}(1-x_j).$$
Clearly, the statement of Lemma~\ref{main lemma2} is equivalent with
$$\min_{\underline{x}\in D_{d,r}}f_{d,r}(\underline{x})=f_{d,r}\left(\frac{r}{d},\dots ,\frac{r}{d}\right)=\binom{d}{r}\bigg(\frac{r}{d}\bigg)^r\bigg(\frac{d-r}{d}\bigg)^{d-r}.$$
We will prove this statement by induction on $d$. The case $d=1$ and $r=0$ or $r=1$ is trivial. In general, the case $r=0$ or $r=d$ is trivial since $D_{d,r}$ consists of only one point in this case. So we can always assume that $0<r<d$. First we prove that the statement is true for the boundary $\partial D_{d,r}$. Then we will prove that  if $\underline{x}\neq \left(\frac{r}{d},\dots ,\frac{r}{d}\right)$, then we either have a point $\underline{y}\in \partial D_{d,r}$ for which $f_{d,r}(\underline{x})\geq f_{d,r}(\underline{y})$ or there exists an $\underline{x}'\in D_{d,r}$ such that $f_{d,r}(\underline{x})> f_{d,r}(\underline{x}')$.
 The compactness of $D_{d,r}$ then implies that $\min_{\underline{x}\in D_{d,r}}f_{d,r}(\underline{x})=f_{d,r}\left(\frac{r}{d},\dots ,\frac{r}{d}\right)$.

So let us first prove that for $\underline{x}\in \partial D_{d,r}$ we have $f_{d,r}(\underline{x})\geq f_{d,r}\left(\frac{r}{d},\dots ,\frac{r}{d}\right)$.
Clearly, if $\underline{x}\in \partial D_{d,r}$ then one of its coordinates is $0$ or $1$. If we delete this coordinate, then the obtained vector $\underline{x}'$ is in $D_{d-1,r}$ in the first case, and in $D_{d-1,r-1}$ in the second case. Furthermore, $f_{d,r}(\underline{x})=f_{d-1,r}(\underline{x}')$ in the first case, and $f_{d,r}(\underline{x})=f_{d-1,r-1}(\underline{x}')$ in the second case. By induction we know that we have
$$f_{d-1,r}(\underline{x}')\geq \binom{d-1}{r}\bigg(\frac{r}{d-1}\bigg)^r\bigg(\frac{d-1-r}{d-1}\bigg)^{d-1-r}$$
for $\underline{x}'\in D_{d-1,r}$, and
$$f_{d-1,r-1}(\underline{x}')\geq \binom{d-1}{r-1}\bigg(\frac{r-1}{d-1}\bigg)^{r-1}\bigg(\frac{d-r}{d-1}\bigg)^{d-r}$$
for $\underline{x}'\in D_{d-1,r-1}$. Let us introduce the function $\ell_{n,k}=x^k(1-x)^{n-k}$. It is easy to see that it takes its maximum at the value $x=k/n$ in the interval $[0,1]$ as its derivative is $(k-nx)x^{k-1}(1-x)^{n-k-1}$. Hence
\begin{align*}
\binom{d}{r}\bigg(\frac{r}{d}\bigg)^r\bigg(\frac{d-r}{d}\bigg)^{d-r}&=\binom{d-1}{r}\bigg(\frac{r}{d}\bigg)^r\bigg(\frac{d-r}{d}\bigg)^{d-r-1}\\
&=\binom{d-1}{r}\ell_{d-1,r}\left(\frac{r}{d}\right)\\
&<\binom{d-1}{r}\ell_{d-1,r}\left(\frac{r}{d-1}\right)\\
&=\binom{d-1}{r}\bigg(\frac{r}{d-1}\bigg)^r\bigg(\frac{d-1-r}{d-1}\bigg)^{d-1-r}
\end{align*}
and
\begin{align*}
\binom{d}{r}\bigg(\frac{r}{d}\bigg)^r\bigg(\frac{d-r}{d}\bigg)^{d-r}&=\binom{d-1}{r-1}\bigg(\frac{r}{d}\bigg)^{r-1}\bigg(\frac{d-r}{d}\bigg)^{d-r}\\
&=\binom{d-1}{r-1}\ell_{d-1,r-1}\left(\frac{r}{d}\right)\\
&<\binom{d-1}{r-1}\ell_{d-1,r-1}\left(\frac{r-1}{d-1}\right)\\
&=\binom{d-1}{r-1}\bigg(\frac{r-1}{d-1}\bigg)^{r-1}\bigg(\frac{d-r}{d-1}\bigg)^{d-r}.
\end{align*}
Hence for $\underline{x}\in \partial D_{d,r}$ we have $f_{d,r}(\underline{x})\geq f_{d,r}\left(\frac{r}{d},\dots ,\frac{r}{d}\right)$.

Next we show that if $\underline{x}\neq \left(\frac{r}{d},\dots ,\frac{r}{d}\right)$, then either
we have a point $\underline{y}\in \partial D_{d,r}$ for which $f_{d,r}(\underline{x})\geq f_{d,r}(\underline{y})$ or
there exists an $\underline{x}'\in D_{d,r}$ such that $f_{d,r}(\underline{x})> f_{d,r}(\underline{x}')$. Since $\underline{x}\neq \left(\frac{r}{d},\dots ,\frac{r}{d}\right)$ there exists $i$ and $j$ such that $x_i\neq x_j$. Let $x_i+x_j=u$ and
$$\prod_{k\neq i,j}(1-x_k+x_kz)=\sum_{m=0}^{d-2}b_mz^m.$$
Then $f_{d,r}(\underline{x})$ is the coefficient of $z^r$ in the polynomial
$$\left(\sum_{m=0}^{d-2}b_mz^m\right)(1-x_i+zx_i)(1-(u-x_i)+z(u-x_i)).$$
A little computation shows that
$$f_{d,r}(\underline{x})=x_i(u-x_i)(b_{r-2}-2b_{r-1}+b_r)+(rb_{r-1}+(1-r)b_r).$$
If $b_{r-2}-2b_{r-1}+b_r<0$, then we get a strictly smaller value for $f_{d,r}(\underline{x})$ if we replace $(x_i,x_j)$ with $(u/2,u/2)$. If $b_{r-2}-2b_{r-1}+b_r\geq 0$, then we can replace $(x_i,x_j)$ with $(0,u)$ or $(1,u-1)$ depending on $u\in [0,1]$ or $[1,2]$ yielding a boundary point $\underline{y}$ for which $f_{d,r}(\underline{x})\geq f_{d,r}(\underline{y})$. This completes the proof.

\section{Capacity of the polynomial $P_G$}

In this section we prove Theorems~\ref{matching} and \ref{orientation} by computing the capacity of $P_G(\underline{x})=\prod_{(u,v)\in E}(x_u+x_v)$ with respect to various vectors $\underline{\alpha}$.

\begin{Lemma} \label{lemma-matching} Let $G=(A,B,E)$ be a $d$--regular bipartite graph on $2n$ vertices. Let $\underline{\alpha}$ be the vector that takes value $1$ at a vertex $u\in A$, and value $d-1$ at a vertex $v\in B$. Then
$$\capp_{\underline{\alpha}}(P_G) = \frac{d^{nd}}{(d-1)^{n(d-1)}}.$$
\end{Lemma}

\begin{proof} For sake of convenience let us denote the variables by $x_u$ if $u\in A$, and $y_v$ if $v\in B$. Then
$$\capp_{\underline{\alpha}} (P_G)=\inf_{x_u,y_v>0 \atop u\in A, v\in B} \frac{ \prod_{(u,v)\in E(G)}(x_u+y_v)}{\prod_{u\in A}x_u \cdot  \prod_{v\in B}y_v^{d-1}}=\inf_{x_u,y_v>0 \atop u\in A, v\in B}\prod_{(u,v)\in E(G)}\frac{x_u+y_v}{x_u^{1/d}y_v^{(d-1)/d}}.$$
Note that
$$x_u+y_v=\frac{1}{d}(dx_u)+\frac{d-1}{d}\left(\frac{dy_v}{d-1}\right)\geq (dx_u)^{1/d}\left(\frac{dy_v}{d-1}\right)^{(d-1)/d}$$
by weighted arithmetic-geometric mean inequality. In other words,
$$\frac{x_u+y_v}{x_u^{1/d}y_v^{(d-1)/d}}\geq \frac{d}{(d-1)^{(d-1)/d}}.$$
Hence
$$\capp_{\underline{\alpha}} (P_G)\geq \left(\frac{d}{(d-1)^{(d-1)/d}}\right)^{nd}=\frac{d^{nd}}{(d-1)^{n(d-1)}}.$$
Observe that if $x_u=1$ and $y_v=d-1$ for all $u\in A$ and $v\in B$ then this bound is sharp.
\end{proof}

\begin{proof}[Proof of Theorem~\ref{matching}]
As before let  $\underline{\alpha}$ be the vector that takes value $1$ at a vertex $u\in A$, and values $d-1$ at a vertex $v\in B$.
By Theorem~\ref{general} we have
$$\prm(G)\geq  \prod_{v\in A}\binom{d}{1}\bigg(\frac{1}{d}\bigg)^{1}\bigg(\frac{d-1}{d}\bigg)^{d-1}\cdot \prod_{v\in B}\binom{d}{d-1}\bigg(\frac{d-1}{d}\bigg)^{d-1}\bigg(\frac{1}{d}\bigg)^{1}\cdot \mathrm{cap}_{\underline{\alpha}}(P_G).$$
Then by Lemma~\ref{lemma-matching} we have
$$\prm(G)\geq \left(d\bigg(\frac{d-1}{d}\bigg)^{d-1}\right)^{2n}\frac{d^{nd}}{(d-1)^{n(d-1)}}=\bigg( \frac{(d-1)^{d-1}}{d^{d-2}} \bigg)^{n}.$$
\end{proof}

\begin{Lemma} \label{lemma-orientation} Let $\underline{\alpha}=(d_1/2,\dots ,d_n/2)$, where $d_j$ is the degree of the vertex $j$. Then
$\capp_{\underline{\alpha}} (P_G) = 2^{e(G)}$,
where  $e(G)$ is the number of edges.
\end{Lemma}

\begin{proof}
First we will prove that $\capp_{\underline{\alpha}}(P_G) \leq 2^{e(G)}$.
Indeed, if we substitute $\underline{1}$ as $\underline{x}$, then we get
$$\frac{P_G(1,1, \dots 1)}{\prod_{v\in V(G)}1^{d_v/2}}=\prod_{(i,j)\in E(G)}(1+1)=2^{e(G)}.$$
Next we will prove the other direction: $\capp _{\underline{\alpha}}(P_G) \geq 2^{e(G)}$.
Using that $x_i+x_j\geq 2\sqrt{x_ix_j}$ we get that
$$\frac{\prod_{(u,v)\in E}(x_u+x_v)}{\prod_{v\in V(G)}x_v^{d_v/2}} \geq \frac{\prod_{(i,j)\in E(G)} 2\sqrt{x_ix_j}}{\prod_{v\in V(G)}x_v^{d_v/2}}=2^{e(G)}.$$
Thus $\capp_{\underline{\alpha}} (P_G) = 2^{e(G)}$.
\end{proof}

\begin{proof}[Proof of Theorem~\ref{orientation}]
As before let $\underline{\alpha}=(d_1/2,\dots ,d_n/2)$, where $d_j$ is the degree of the vertex $j$.
By Theorem~\ref{general} we have
$$\e(G)\geq  \prod_{v\in V}\binom{d_v}{d_v/2}\bigg(\frac{d_v/2}{d_v}\bigg)^{d_v/2}\bigg(\frac{d_v/2}{d_v}\bigg)^{d_v/2}\cdot \mathrm{cap}_{\underline{\alpha}}(P_G).$$
Then by Lemma~\ref{lemma-orientation} we have
$$\e(G)\geq \prod_{v\in V}\frac{\binom{d_v}{d_v/2}}{2^{d_v}}\cdot 2^{e(G)}=\prod_{v\in V}\frac{\binom{d_v}{d_v/2}}{2^{d_v/2}}.$$
\end{proof}

\begin{Rem}
For $d$--regular graphs Las Vergnas \cite{Ver1} improved Theorem~\ref{orientation} as follows:
$$\e(G)\geq \frac{2^d}{\binom{d}{d/2}}\left(\frac{\binom{d}{d/2}}{2^{d/2}}\right)^n.$$
(At Proposition 5.3 of \cite{Ver1} there is a typo as it is pointed out in another paper of Las Vergnas  in the footnote of the first page of \cite{Ver2}.) This strengthening can be obtained by our method too: all we have to note that we only need to apply Theorem~\ref{capacity preserving} to $n-1$ variables corresponding to vertices of $P_G(\underline{x})$: at the very end we should get a polynomial of the form $cx_n^{d_n/2}$ since throughout the process we get homogeneous polynomials. For this univariate polynomial there is no need to apply Theorem~\ref{capacity preserving} once more, so we get a $2^d/\binom{d}{d/2}$ improvement. This argument applies to non-regular graphs too thereby saving a factor $2^{\Delta}/\binom{\Delta}{\Delta/2}$, where $\Delta$ is the largest degree.

\end{Rem}

\section{Beyond this paper}

In this section we collected some pointers to the literature.

Stable polynomials have a huge literature.
If someone is interested in a  comprehensive introduction to the this theory, then the paper of Choe, Oxley, Sokal and Wagner \cite{COSW} or Wagner's survey \cite{Wag} might be a good choice. Another excellent survey of the area is the paper of Vishnoi \cite{Vis}.  For capacity preserver operations the paper of  Gurvits and Leake \cite{LeGu} gives a treatment that is both very general and very readable. Statements on capacity often boils down to some statement about coefficients of univariate real-rooted polynomials often with a probabilistic flavour like Hoeffding's theorem \cite{Hoe}. A good source of such inequalities and results is Pitman's survey \cite{Pit}. The coefficient lemma, Theorem~\ref{coefficient theorem}, already appeared in the paper \cite{Gur3}. Interestingly this paper also considers various versions of the polynomial $P_G(\underline{x})=\prod_{(u,v)\in E(G)}(x_u+x_v)$, but never exactly this form. The coefficient lemma could have been easily derived from the work of  Gurvits and Leake \cite{LeGu} too.

There are many different proofs and generalizations of Theorem~\ref{matching}. The original proof of Schrijver \cite{Sch2} is elementary, but involved. The first proof based on stable polynomials and capacity is due to Gurvits \cite{Gur1}, his proof is simplified in the paper of Laurent and Schrijver \cite{LaSc}. Another proof based on the theory of graph covers is given by Csikv\'ari \cite{Csik}. (The relationship between the theory of graph covers and the theory of stable polynomials is not yet well-understood.) Theorem~\ref{matching} has a very natural generalization for permanents of non-negative matrices. This generalization was derived by Gurvits \cite{Gur2} from the original paper of Schrijver \cite{Sch2}. Subsequently, Anari and Oveis-Gharan \cite{A-OG} and Straszak and Vishnoi \cite{StVi} gave a proof that only relies on the theory of stable polynomials. Another possible generalization considers counting matchings of fixed size in bipartite graphs instead of perfect matchings. This question was treated in the papers Csikv\'ari \cite{Csik}, Lelarge \cite{Lela} and Gurvits and Leake \cite{LeGu} (this last one uses only stable polynomials).

Concerning Eulerian orientations, Theorem~\ref{orientation} has many different proofs either. The original proof of Schrijver is very elegant and simple. Las Vergnas \cite{Ver1}  gave another proof building on the theory of Martin's polynomial that gives a slightly stronger result. Borb\'enyi and Csikv\'ari \cite{BoCs} gave a proof using gauge transformation. The proof presented here is the first one using stable polynomials, but we remark that this result could have been easliy deduced from the paper of Straszak and Vishnoi \cite{StVi} too that uses stable polynomials.

Both Theorems~\ref{matching} and \ref{orientation} can be interpreted as a correlation inequality. If we divide by $2^{e(G)}=\prod_{v\in V(G)} 2^{d_v/2}$ in Theorem~\ref{orientation}, then the the right hand side of the inequality
$$\frac{\e(G)}{2^{e(G)}}\geq \prod_{v\in V(G)}\frac{\binom{d_v}{d_v/2}}{2^{d_v}}$$
is the probability that a random orientation is Eulerian, while on the left hand side the term $\frac{\binom{d_v}{d_v/2}}{2^{d_v}}$ is the probability that a random orientation is balanced at vertex $v$. So this inequality can be interpreted as a positive correlation inequality. Similarly, if $G=(A,B,E)$ is a $d$--regular bipartite graph, then we can consider the probability space where for each vertex $u$ in $A$ we pick exactly one of the edges incident to $u$ uniformly at random. This way we picked $n$ edges.  For a vertex $v\in B$ let $E_v$ be the event that we picked exactly one of the edges incident to $v$. Then $\mathbb{P}(\cap_{v\in B}E_v)=\frac{\prm(G)}{d^n}$ while $\mathbb{P}(E_v)=d\cdot \frac{1}{d}\left(1-\frac{1}{d}\right)^{d-1}=\left(\frac{d-1}{d}\right)^{d-1}$. Hence Theorem~\ref{matching} is equivalent with $\mathbb{P}(\cap_{v\in B}E_v)\geq \prod_{v\in B}\mathbb{P}(E_v)$. It seems that many results on capacity are indeed fuelled by positive correlation inequalities, often in a very disguised way. This connection is the most explicit in the paper \cite{StVi} where the authors introduce the iterated positive correlation property and connect it with stable polynomials. Surprisingly it is also possible to build out a theory of negative correlation based on stable polynomials, for details see the paper  \cite{Pem}.

One might wonder whether there is a deeper connection between Theorems~\ref{matching} and \ref{orientation}. It turns out that both theorems fall into a pattern that is about the so-called Bethe approximation. Hans Bethe  was a Nobel-prize laurate physicist. Among many other achievements he introduced the concept that is now known as Bethe--approximation. Originally this was a highly heuristic concept that approximates well quantities coming from counting objects with local constraints like perfect matchings and Eulerian orientations. For a long time it was overlooked by the mathematics community, then Bethe approximation showed up in two different lines of research. In the work of Dembo, Montanari and their coauthors \cite{DeMa,DeMa2,DMS,DMSS} about graph limit theory of sparse graphs, an appropriate version of Bethe approximation played the role of the limit value of certain graph parameters. Another line of research emerged from the work of Gurvits establishing inequalities between a graph parameter and its Bethe approximation. In many cases Bethe approxiation turns out to be a lower bound for the corresponding graph parameter. Theorem~\ref{matching} and \ref{orientation} belong to this line of research. This area almost exclusively relies on stable polynomials \cite{A-OG,StVi} and graph covers \cite{Csik,Ruo,Lela,Von}.
\bigskip

\noindent \textbf{Acknowledgment.} The first author thanks Jonathan Leake for the discussions on the topic of this paper.
The authors are very grateful to the anonymous referee for his/her suggestions that greatly improved the paper.