\section{Introduction}
\label{introduction}

Many astrophysical observations provide evidences for the existence
of non-baryonic dark matter from the galactic scale to the cosmological
scale~\cite{susydm,evidence,dm}.
Its existence in our universe has been revealed only
by its gravitational effects, while its mass and interactions to the standard
model (SM) particles are still unclear.
Among many scenarios that had been proposed to account for non-baryonic dark matter,
the weakly interacting massive particle (WIMP) scenario,
where the particle with weak scale mass and interactions
could be thermally frozen out and its abundance would be the observed
dark matter relic density~\cite{susydm,evidence},
is one of the most attractive scenarios because it may be produced and detected
at the current or future
colliders. Actually many works in the literature have dealt with this
issue~\cite{collider},
but they have mainly focused on the Large Hadron Colliders (LHC).

The dark matter search at colliders is based on the production of SM particle(s)
with missing energies. The typical signals are
mono-jet~\cite{effective1,effective2,monophoton},
mono-photon~\cite{monophoton}, mono-$W$~\cite{monow}
or mono-$Z$~\cite{monoz} with missing energies.
In principle, the particles which are produced together with the dark matter
particle(s) could be any SM fundamental particles or their composite states.
For example, the dark matter production associated
with a heavy-quark-and-heavy-anti-quark ($Q\bar{Q}$) pair
may also play an important role in searching for dark matter.
It is known that the dark matter pair production associated with the $Q\bar{Q}$
pair could be more effective in the case that
the dark matter interaction with the SM particles depends on
the mass of the relevant SM quark~\cite{monottbar}.

In this paper, we propose to use the process of
dark matter production associated with a heavy quarkonium
to search for dark matter at colliders.
A heavy quarkonium is a meson which is a bound state of
a $Q\bar{Q}$ pair.
Naively speaking, the heavy-quark-pair-associated production of dark matter would be
more effective for the dark matter detection than the light-quark-pair-associated
production if the effective coupling of dark matter to the quark pair
is proportional to the quark mass.
In a similar way,  the dark matter production associated with a heavy quarkonium
would play a more essential role in the probe of such operators
than dark matter production associated with a light meson.

The dark matter production associated with the SM particle(s) can be
investigated in any collider experiments.
In this work, we focus on the dark matter production
associated with a heavy quarkonium at $B$ factories,
whose center-of-mass (CM) energy is $\sqrt{s}=10.58$ GeV.
However, we note that this search will easily be
applied to the International Linear Collider (ILC) and Large Hadron Collider
(LHC)~\cite{future}.

There is no evidence for dark matter from the collider experiments so far.
However, there are some reports for dark matter candidate signals
from direct detection of dark matter,
which is carried out in underground experiments,
and also from indirect detection of dark matter from astrophysical observations,
in particular, in the dark matter mass region of $\sim O(10)$ GeV.
On the other hand, the dark matter mass region below a few GeV
has not been well investigated, especially, in the direct detection
experiments of dark matter because the very light dark matter cannot
hit a nucleon with significant recoil.
One of the strongest bounds on the dark matter signals
in the indirect detection of dark matter comes from the $\gamma$-ray
observation of  dwarf spheroidal galaxies
by the Fermi-LAT satellite~\cite{dwarf}.
However, the constraints become weakened below the dark matter mass region
less than about 4 GeV because of the uncertainty of hadronization
of final state particles in the dark matter annihilation
and the limit of the photon energy threshold of the detector in the Fermi Gamma-ray Space
Telescope~\cite{Hdecay}.

Since such light dark matter mass region could be investigated
at colliders~\cite{light}, the collider detection of dark matter would be complementary
to the direct detection of dark matter and the indirect detection
of dark matter.

At $B$ factories, dark matter cannot only be produced directly,
but also produced from the decays of a heavy quarkonium.
For example, the invisible decay of $\Upsilon(1S)$ can give constraints
on the properties of light dark matter through the decay of
$\Upsilon(3S) \to \pi^+ \pi^- \Upsilon(1S), \gamma\Upsilon(1S)$
followed by $\Upsilon(1S) \to $
nothing at $B$ factories~\cite{Hdecay}.
Also, the charm factories can contribute
to search for dark matter by investigating the invisible decay of a charmonium,
$J/\psi$~\cite{Hdecay}.
It is found that the $\Upsilon(1S)$ and $J/\psi$ decays at $B$ and charm
factories are more suitable to probe
light dark matter interactions than the mono-jet search at high energy hadron
colliders when the mass of mediator is not large~\cite{Hdecay}.
The dark matter production associated with a heavy quarkonium at $B$ factories
would also be complementary to the invisible decay of a heavy quarkonium
as well as the mono-jet searches at high energy colliders.

This paper is organized as follows.
In Sec.~\ref{model}, we describe the relevant effective operators
for the dark matter
interactions with the SM quarks.
As a candidate for dark matter, we take into account two cases: Dirac
fermionic dark matter and real scalar dark matter.
We adopt the effective field theory (EFT) approach for the interaction of dark matter
with the SM quarks.
In Sec.~\ref{calc}, we calculate the cross sections for
the dark matter production at $B$ factories
as well as for the SM backgrounds. We also interpret the constraints
on the properties of the dark matter interaction
to the upper limit on the dark matter-nucleon scattering cross section.
Finally, we summarize our results in Sec.~\ref{summary}.

\section{Effective operators}
\label{model}

For a model-independent search for dark matter, we assume that
$\chi$, which stands for either a Dirac fermion or real scalar, is
the only component of dark matter and singlet under the SM gauge group.
The extension to the Majorana fermion
or complex scalar dark matter would be straightforward and not so much
different from the Dirac fermion or real scalar studied here.
We adopt the effective field
theory approach, where the mediator connecting the dark sector to the SM
particles is heavy enough to be integrated
out~\cite{effective1,effective2,effective3}.
After symmetry breaking, the dark matter interactions to the SM fields
can be described by higher-dimensional operators.
Since we are interested in the dark matter production associated
with a heavy quarkonium, the relevant effective operators must contain
the heavy quark fields.

The dimension-6 operators in the Dirac fermion case which we will consider
in this paper are~\cite{effective2}
\begin{eqnarray}
O_{1(2)}^D = (i)\frac{m_Q}{\Lambda^3} \bar{\chi} (\gamma^5)\chi \bar{Q} Q,
&~~&
O_{3(4)}^D = (-i)\frac{i m_Q}{\Lambda^3} \bar{\chi}(\gamma^5) \chi \bar{Q}\gamma^5  Q,
\nonumber \\
O_{5(6)}^D = \frac{1}{\Lambda^2} \bar{\chi}\gamma_\mu (\gamma^5) \chi
\bar{Q}\gamma^\mu  Q,
&~~&
O_{7(8)}^D = \frac{1}{\Lambda^2} \bar{\chi}\gamma_\mu (\gamma^5) \chi
\bar{Q}\gamma^\mu \gamma^5 Q,
\nonumber \\
O_{9}^D = \frac{1}{\Lambda^2} \bar{\chi}\sigma_{\mu\nu} \chi
\bar{Q}\sigma^{\mu\nu}  Q,
&~~&
O_{10}^D = \frac{i}{\Lambda^2} \bar{\chi}\sigma_{\mu\nu}\gamma^5 \chi
\bar{Q}\sigma^{\mu\nu}  Q,
\label{effopfermion}%
\end{eqnarray}
where $m_Q$ is the heavy-quark mass and
$\Lambda$ represents the energy scale of new physics.
Effectively $\Lambda \sim M /\sqrt{g_1 g_2}$, where $M$ is the mass
of the mediator and $g_i$ are the coupling of the mediator to
the dark sector and SM sector, respectively.
For the (pseudo)scalar interaction operators $O_{1,2,3,4}^D$,
an additional helicity suppression factor is taken into account.

In the real scalar case, the relevant dimension-6 operators
are given by~\cite{effective2}
\begin{eqnarray}
O_{1}^R = \frac{m_Q}{2\Lambda^2} \chi^2 \bar{Q}Q,
&~~&
O_{2}^R = i\frac{m_Q}{2\Lambda^2} \chi^2 \bar{Q}\gamma^5 Q.
\label{effopscalar}%
\end{eqnarray}

There are more dimension-6 operators which would be necessary
for complete study of the dark matter production.
For example, one can consider the operators which include
the gluonic field strengths ($G_{\mu\nu}G^{\mu\nu}$).
They may contribute to the production of dark matter particles associated
with a heavy quarkonium at the one-loop level so that the dark matter
production process would be suppressed.
For the probe of the operators which include the electromagnetic field strengths
($F_{\mu\nu}F^{\mu\nu}$) or a lepton pair,
the dark matter production associated with a photon would be more efficient
at the lepton colliders like $B$ factories or ILC.

It is well known that the EFT approach is valid
when the particle mediating the interaction is heavier than
the typical energy scale of the process, which is the maximum momentum transfer
in the process.
At high-energy hadron collider where the energy scale of the fundamental process
is not fixed, the EFT approach may break down unless the mediator is
not too heavy compared to the CM energy of the collider.
Then, the EFT approach may over- or under-estimate the result in a UV-complete theory~\cite{break}.
It is worthwhile to note that $B$ factories are relatively free from
this issue because the CM energy is fixed as well as not so high.

Some operators in Eqs.~(\ref{effopfermion}) and (\ref{effopscalar})
do not respect the SM gauge group symmetry~\cite{violate}.
In particular, the scalar operators $O^D_{1,2,3,4}$ require
a $s$- or $t$ channel scalar exchange for UV completion.
In Ref.~\cite{monojetsf}, it is shown that in the $s$-channel scalar-exchange
model, the EFT approach prediction is completely different from that
of the UV complete model.
Therefore, one should be more cautious about the interpretation of the bounds on such operators
obtained in the EFT approach.

\section{Calculation}
\label{calc}

\subsection{Standard Model Backgrounds}

The production of a heavy quarkonium is described by non-relativistic QCD
(NRQCD), which is an effective field theory of QCD~\cite{nrqcd}.
As usual, the factorization between short-distance and long-distance physics is assumed.
While the factorization between high-energy and low-energy physics
for the inclusive heavy-quarkonium production has not been proved, the factorization for its exclusive production at $B$ factories was given in~\cite{factorization}.

The schematic expression for the production cross section is
given by $d\sigma \sim \sum c_O \langle O \rangle_H$, where
$\langle O \rangle_H$ is the long-distance matrix element (LDME),
which is the probability of the evolution of a heavy-quark pair into a heavy
quarkonium $H$,
and $c_O$ is the corresponding short-distance coefficient, which is
responsible for the production of a heavy-quark and anti-quark pair.
The coefficient $c_O$ is obtained by integrating the squared amplitude over
the phase space after averaging the spin states of initial particles
and summing the spin states of final particles.
The amplitude for the production of $H$ can be obtained by
projecting the production amplitude of a $Q\bar{Q}$ pair
onto $H$ via the spin and color projection operators of $H$ (generically denoted by $\Pi_H$ here):
\begin{equation}
{\cal M} \sim \textrm{Tr}[{\cal A} \Pi_H],
\end{equation}
with ${\cal A}$ being a matrix that acts on spinors of
the $Q\bar{Q}$ pair in the amplitude of the $Q\bar{Q}$ pair production~\cite{BL1,resum}.

In NRQCD, a heavy-quark and anti-quark pair can be produced in a color-singlet
state or color-octet state. The color-octet state can evolve into
a color-singlet state by emitting or absorbing soft gluons.
Typically, the color-octet LDME is suppressed by $v^3$ or higher
in the heavy-quark velocity $v$ compared to the color-singlet LDME.
Therefore, unless there exists a kinematical enhancement factor in the short-distance
coefficient for the color-octet amplitude,
the color-singlet amplitude dominates.
In the $B$ factories, the initial colliding beams are the electron
and positron, which are neutral under $SU(3)_C$ gauge symmetry.
Therefore, the production of the heavy quarkonium in the color-singlet state
can occur only when the other particles are color-neutral or
in the color-singlet state.
However, the production of a heavy quarkonium in a color-octet state
must accompany with at least a color-octet particle,
for instance, an additional gluon or a pair of quark
and anti-quark in the color-octet state.
In the dark matter production or off-shell $Z$ boson production at $B$
factories, the color-singlet heavy quarkonium production is dominant.
The color-octet heavy quarkonium production has a wavefunction suppression
factor of $O(v^3)$ or higher.
Furthermore, the color-octet process has an additional suppression factor
due to the additional color-octet gluon or a pair of quark and anti-quark.
Therefore, for the dark matter production associated with a heavy quarkonium
at $B$ factories, it is sufficient to take into account
only the color-singlet contribution.

The color-singlet LDME is determined by the electromagnetic decay of the heavy
quarkonium. Explicitly, for charm quark mass $m_c=1.5$ GeV  and to leading-order
in $v$, we take the color-singlet LDMEs
to be~\cite{LDME,LDME2}
\begin{eqnarray}
\langle O \rangle_{\eta_c}  = 0.474~\textrm{GeV}^3,
&~~&
\langle O \rangle_{J/\psi}  = 0.436~\textrm{GeV}^3.
\end{eqnarray}

\begin{figure}[!t]
\begin{center}
{\epsfig{figure=diagramsm.eps,width=0.7\textwidth}}
\end{center}
\vspace{-0.5cm}
\caption{
The Feynman diagrams for $e^+ e^- \to H (J/\psi,\eta_c) \nu \bar{\nu}$ in the SM.
Other diagrams are obtained by reversing the flow of fermions.
}
\label{diagramsm}
\end{figure}

First we calculate the SM backgrounds.
In the SM, the Feynman diagrams contributing to the heavy quarkonium
production with missing energy have at least one weak-boson exchange,
which is connected to a pair of neutrinos.
There are also  diagrams with two weak-boson exchanges, but they
are suppressed by $O(s/M_Z^2)$ so that we ignore those diagrams.
Then, we left with seven diagrams in the SM, whose
representative diagrams are shown in Fig.~\ref{diagramsm}.
The irreducible background for the dark matter production
associated with a heavy quarkonium is the $e^+ e^- \to H Z^{\ast}$ process
followed by $Z^{\ast}\to \nu\bar{\nu}$, as depicted in the first two diagrams in Fig.~\ref{diagramsm}.
There is also a $W$ fusion diagram, which is required by gauge invariance,
as shown in the last diagram in Fig.~\ref{diagramsm}.
Other diagrams not shown are obtained by reversing the flow of fermions.

In the monojet production at hadron colliders, except the irreducible
SM background $pp\to Z+\textrm{jets}$ followed by $Z\to \nu\bar{\nu}$,
there is another type of the SM background, for example,
the process $pp\to W+\textrm{jets}, W\to l \nu$, where both $\nu$ and
$l$ are missed. We note that this kind of background can be ignored
at $B$ factories because the relevant process is
$e^+ e^- \to H W^{+\ast} W^{-\ast}$, which is highly suppressed.

We note that, at the $e^+ e^-$ machines, there may be continuum background coming from $e^+ e^- \to q \bar{q} \nu\bar{\nu}$.
We calculate this continuum background by choosing
the invariant mass  $m_{q\bar{q}}$ of the $q\bar{q}$ pair in the region
$m_{H} - 5 \Gamma_H \le m_{q\bar{q}} \le m_H + 5 \Gamma_H$, where
$m_H$ and $\Gamma_H$ are the mass and decay width
of $H$ respectively for $H=J/\psi$ or $\eta_c$. We find that the cross section for the continuum background is
less than $3.5\times 10^{-3}$ ab, which can be neglected at $B$ factories.
 There may also be a continuum background coming from process like $e^+ e^- \to \ell^+ \ell^- \nu\bar{\nu}$ ($\ell= e, \mu$)
 since $J/\psi$ may be detected experimentally via the decay $J/\psi \to \ell^+ \ell^-$.
Nevertheless the cross section from this continuum contribution is also negligible.

We take the electromagnetic coupling to be $\alpha=1/132$ at the scale
$\mu=\sqrt{s}$ and the charm quark mass to be $m_c=1.5$ GeV.
Then, the cross sections for a heavy quarkonium production associated
with a pair of neutrinos in the SM are calculated as (1 ab = 1 attobarn = $10^{-42}$ cm$^2$)
\begin{eqnarray}
\sigma(e^+ e^- \to J/\psi \nu\bar{\nu}) &=& 0.81~\textrm{ab},
\\
\sigma(e^+ e^- \to \eta_c \nu\bar{\nu}) &=& 2\times 10^{-5}~\textrm{ab}.
\end{eqnarray}
The cross section for $\eta_c\nu\bar{\nu}$ production is much smaller than
that for the $J/\psi\nu\bar{\nu}$ production.
This is mainly because of the dominated contributions from photon-fragmentation diagrams, where
a virtual photon evolves into a $J/\psi$~\cite{resum,LDME2,photonfrag}.
All the diagrams in Fig.~\ref{diagramsm} contribute
to the $J/\psi \nu\bar{\nu}$ production
while only the first diagram contributes to the $\eta_c \nu\bar{\nu}$
production. However, the first diagram does not correspond
to the photon-fragmentation contribution and is suppressed by $O(m_c^2/s)$
as compared to other photon-fragmentation diagrams.

The expected event number for the $J/\psi\nu\bar{\nu}$ production is about
$1$ $(40)$ with an integrated luminosity of $1$ $(50)$ ab$^{-1}$
without any cut. Since the expected event number for the $\eta_c\nu\bar{\nu}$
production is entirely negligible, even one event of the $\eta_c$ plus missing energy
measured at $B$ factories would imply the existence of new physics.

\begin{figure}[!t]
\begin{center}
{\epsfig{figure=diagram.eps,width=0.3\textwidth}}
\end{center}
\vspace{-0.5cm}
\caption{
The Feynman diagram for $e^+ e^- \to J/\psi + \chi +\bar{\chi}$.
Another diagram is obtained by reversing the flow of the charm quark.
The bulb represents an effective operator.
}
\label{diagram}
\end{figure}

\subsection{The Dirac Fermion Dark Matter Case}

In this section, we calculate the cross section for the Dirac fermion dark matter production
associated with a heavy quarkonium.

There are two Feynman diagrams, one of
which is shown in Fig.~\ref{diagram}. The other diagram is obtained
by reversing the flow of the charm quark.
The bulb in the figure represents vertices of various effective operators.

The relevant operators to the heavy quarkonium production depend
on the quantum number of the heavy quarkonium,
where $J^{PC}=1^{--}$ for $J/\psi$ and $J^{PC}=0^{+-}$ for $\eta_c$.
For example, the operators $O_{1,2,3,4,7,8}^D$ contribute only to
the $e^+ e^- \to J/\psi \chi \bar{\chi}$ process,
while the operators $O_{5,6,9,10}^D$
only to the $e^+ e^- \to \eta_c \chi\bar{\chi}$.
We note that only the $O_{5,6,9,10}^D$ operators
can contribute to the invisible decays of the heavy quarkonium states
$\Upsilon(1S)$ and $J/\psi$~\cite{Hdecay}.
Thus, the dark matter production associated with a heavy
quarkonium could provide constraints on a larger set of operators, which
cannot be probed in the invisible decays of $\Upsilon(1S)$ and $J/\psi$.

\begin{figure}[!t]
\begin{center}
{\epsfig{figure=psixsec.eps,width=0.4\textwidth}}
{\epsfig{figure=etacxsec.eps,width=0.4\textwidth}}
\end{center}
\vspace{-0.5cm}
\caption{
The cross sections (a) for $e^+ e^- \to J/\psi  \chi \bar{\chi}$
and (b) for $e^+ e^- \to \eta_c \chi \bar{\chi}$
in unit of fb as a function of $m_\chi$ in unit of GeV.
The red and blue lines correspond to the $O_1^D (O_5^D)$
and $O_8^D (O_9^D)$
operators for the $J/\psi \chi \bar{\chi}$ ($\eta_c\chi\bar{\chi}$)
production, respectively.
The dashed, solid, and dotted lines
correspond to $\Lambda=20$, $100$, and $200$ GeV, respectively.
The green lines correspond to the $O_1^R$ operator for the real scalar dark
matter.
}
\label{figurexsec}
\end{figure}

Figure~\ref{figurexsec} shows the cross sections
(a) for $e^+ e^- \to J/\psi  \chi \bar{\chi}$ and
(b) for $e^+ e^- \to \eta_c  \chi \bar{\chi}$
as a function of the dark matter mass $m_\chi$ for $\Lambda=20$ GeV (dashed line),
$100$ GeV (solid line), and $200$ GeV (dotted line), respectively.
The effective operators for the red and blue lines
are $O_1^D$ ($O_5^D$) and $O_8^D$ ($O_9^D$) for the $J/\psi$ ($\eta_c$)
production, respectively.
The horizontal line in Fig.~\ref{figurexsec}(a) is the SM background
in the framework of NRQCD

In Fig.~\ref{figurexsec}(a), the cross sections
for the $J/\psi \chi\bar{\chi}$ production can reach
about $0.08$ fb for $\Lambda= 20$ GeV and about $0.01$ ab for $\Lambda=100$ GeV
respectively in the case of the operator $O_8^D$.
However, for the $O_1^D$ operator, the cross section
reaches about $0.4$ ab for $\Lambda=20$ GeV, but the production is negligible
for $\Lambda=100$ GeV due to a suppression factor $O(m_c/\Lambda)$.
The curves in the figure drop rapidly as the dark matter mass grows, indicating
the kinematical limit of the dark matter mass of $(\sqrt{s}-m_{J/\psi})/2$.
For other operators $O_{2,3,4,7}^D$, similar features can be obtained.

In Fig.~\ref{figurexsec}(b), the cross section for the $\eta_c\chi\bar{\chi}$
production can reach about $0.02$ fb for $\Lambda=20$ GeV
and about $10$ ab for $\Lambda=100$ GeV in the case of the $O_5^D$ operator.
For the $O_9^D$ operator, the cross section reaches about
$0.1$ fb for $\Lambda=20$ GeV and about $0.17$ ab for $\Lambda=100$ GeV,
respectively. Effects from the kinematical limit of the dark matter mass of
$(\sqrt{s}-m_{\eta_c})/2$ are also evident.
For operators $O_{6,10}^D$, similar features are obtained.

\begin{figure}[!t]
\begin{center}
{\epsfig{figure=psilam.eps,width=0.4\textwidth}}
{\epsfig{figure=etaclam.eps,width=0.4\textwidth}}
\end{center}
\vspace{-0.5cm}
\caption{
Bounds on $\Lambda$ (a) for $e^+ e^- \to J/\psi  \chi \bar{\chi}$
and (b) for $e^+ e^- \to \eta_c \chi \bar{\chi}$
as a function of $m_\chi$, determined by the signal-to-background ratio $R=5$ and 1 respectively.
The red and blue lines correspond to the $O_1^D (O_5^D)$ and $O_8^D (O_9^D)$
operators in the $J/\psi \chi \bar{\chi}$ ($\eta_c\chi\bar{\chi}$)
production, respectively.
The green lines correspond to the $O_1^R$ operator in the $J/\psi \chi \chi$
production where $\chi$ is a real scalar.
The dashed and solid lines
correspond to the integrated luminosity of $1$ ${\rm ab}^{-1}$ and $50$ ${\rm ab}^{-1}$, respectively.
}
\label{figurelam}
\end{figure}

In Fig.~\ref{figurelam}, we plot bounds on the scale $\Lambda$
(a) for $e^+ e^- \to J/\psi \chi \bar{\chi}$
and (b) for $e^+ e^- \to \eta_c \chi \bar{\chi}$
in the cases of the integrated luminosity (${\cal L}$) of $1$ ${\rm ab}^{-1}$ (dashed line)
and $50$ ${\rm ab}^{-1}$ (solid line), respectively.
In the $J/\psi \chi\bar{\chi}$ production, we set the bound on $\Lambda$
to be the value where the signal-to-background ratio, which is defined by
$R={\cal L} \sigma_\textrm{sig}/\sqrt{{\cal L}\sigma_\textrm{bg}}$, is equal
to 5.
In the $\eta_c\chi\bar{\chi}$ production, we set the bound
to be the value where the number of events is one because the SM background
is negligible.
In the $J/\psi\chi\bar{\chi}$ production, the bound can reach
about $42$ ($67$) GeV with ${\cal L}=1$ ($50$) ${\rm ab}^{-1}$ for the $O_8^D$ operator.
For the $O_1^D$ operator, the bound is much less due to the chirality suppression
factor in the operator.
In the $\eta_c\chi\bar{\chi}$ production, the bound can reach
about $42$ ($111$) GeV for $O_5^D$ operator
and about $64$ ($169$) GeV for $O_9^D$ operator
with ${\cal L}=1$ ($50$) ${\rm ab}^{-1}$,  respectively.

A couple of comments are in order here.
First, as the dark matter mass approaches the kinematical limit,
the bound on $\Lambda$ which can be achieved at $B$ factories becomes quite tiny
as clearly shown in Fig.~\ref{figurelam}.
We note that for the region below some value of $\Lambda$, our limits become untrustworthy.
Recall that $\Lambda \sim M / \sqrt{g_1g_2}$. A very tiny $\Lambda$ can be achieved
only when $g_1 g_2$ becomes large since
the mediator mass $M$ should be larger than $\sqrt{s}$
for validity of EFT approach. On the other hand,
$g_{1,2}$ cannot be too large in order to perform perturbative calculation.
For $g_{1,2} \sim 1$, the EFT description breaks down
in the region $\Lambda \lesssim \sqrt{s}$.
Second, we note that the number of events is not large. So the best thing
one can hope for is to detect inclusive signals rather than exclusive ones.
That is, the heavy quarkonium would be identified by the invariant mass distribution of
all the final state particles in the decay products (except for missing energy)
instead of using the exclusive mode like $J/\psi \to \l^+l^-$.

\subsection{The Scalar Dark Matter Case}

In this section, we consider the real scalar dark matter production associated
with a heavy quarkonium at $B$ factories. The relevant Feynman diagrams
are the same as in the Dirac dark matter case with the replacement of
the dark matter fermion lines by the dark matter scalar lines.
The relevant effective operators are shown in Eq.~(\ref{effopscalar}).

In Fig.~\ref{figurexsec}(a), the cross section for the $J/\psi\chi\chi$
production in the $O_1^R$ operator case
is depicted for $\Lambda=20$ GeV (green dashed line) and $\Lambda=100$ GeV
(green solid line). The cross section can reach about $2$ ab for $\Lambda=20$ GeV,
but for $\Lambda=100$ GeV the cross section is below $O(10^{-2})$ ${\rm ab}^{-1}$.
Similar feature is observed for the $O_2^R$ operator.
We note that the cross section for the $\eta_c \chi\chi$ production
vanishes for real dark matter.

The bound on the scale $\Lambda$ for $e^+ e^- \to J/\psi \chi \chi$
is shown in Fig.~\ref{figurelam}(a) for the integrated luminosity
of 1 ${\rm ab}^{-1}$ (green dashed line) and 50 ${\rm ab}^{-1}$ (green solid line), respectively.
Again the bound is set to be the signal-to-background ratio $R=5$.
The bound on $\Lambda$ can reach about $17$ $(28)$ GeV
for ${\cal L}=1$ $(50)$ ${\rm ab}^{-1}$ respectively
for the operator $O^R_1$.

\begin{figure}[!t]
\begin{center}
{\epsfig{figure=si.eps,width=0.4\textwidth}}
{\epsfig{figure=sd.eps,width=0.4\textwidth}}
\end{center}
\vspace{-0.5cm}
\caption{
Experimental reach of (a) the spin-independent cross section
and (b) the spin-dependent cross section
for the dark matter-nucleon scattering in units of cm$^2$
as a function of $m_\chi$ in unit of GeV.
The red (blue) lines correspond to the $O_1^D (O_5^D)$ and $O_8^D (O_9^D)$
operators in the $J/\psi \chi \bar{\chi}$ ($\eta_c\chi\bar{\chi}$)
production, respectively.
The green lines in (a) correspond to the $O_1^R$ operator
in the $J/\psi \chi \bar{\chi}$ production.
The dashed and solid lines
correspond to the integrated luminosity of $1$ ${\rm ab}^{-1}$ and $50$ ${\rm ab}^{-1}$, respectively.
The black and gray lines on the lower right corners are 90\% C.L. exclusion limits
of second CDMSlite run~\cite{cdmslite2} and superCDMS~\cite{supercdms}, respectively.
The bounds from the monojet search at the LHC are shown in
the orange dotted (solid) lines for $O^D_5$ ($O^D_8$) from CMS,
the cyan dotted (solid) lines for $O^D_5$ ($O^D_8$) from ATLAS,
and the purple dotted (solid) lines for $O^D_1$ ($O^D_9$) from ATLAS,
respectively~\cite{cmsmonojet,atlasmonojet} from 1 to 10 GeV.
}
\label{figuredirect}
\end{figure}

\subsection{Conversion to the Dark Matter-Nucleon Cross Section}

The bound on $\Lambda$ deduced above may be converted into the bound for
the dark matter-nucleon scattering cross section for the direct detection of dark matter.
We note that, however, this conversion holds with the assumption
that the effective operators in Eqs.~(\ref{effopfermion}) and (\ref{effopscalar})
are either flavor blind or proportional to the quark masses.
In case that dark matter has charm-philic property, the bound from
the dark matter production associated with a heavy quarkonium cannot be
converted into the cross section for the dark matter-nucleon scattering.

For the effective operators, the dark matter-nucleon scattering
cross sections are shown in Refs.~\cite{belanger,Haisch:2012kf}:
\begin{eqnarray}
\label{dmnucleona}
\sigma_\textrm{SI}{(O^D_1)} &=&
\frac{\mu_\chi^2 m_n^2}{\pi\Lambda^6}f_n^2,
\\
\label{dmnucleonb}
\sigma_\textrm{SI}{(O^D_5)} &=&
\frac{\mu_\chi^2}{\pi \Lambda^4} \left(\sum_{q}f_{V_q}^n\right)^2,
\\
\label{dmnucleonc}
\sigma_\textrm{SI}{(O^R_1)} &=&
\frac{4 \mu_\chi^2 m_n^2}{\pi \Lambda^4 m_\chi^2} f_n^2,
\\
\label{dmnucleond}
\sigma_\textrm{SD}{(O^D_{8,9})} &=&
\frac{3 \mu_\chi^2}{\pi \Lambda^4} \left(\sum_{q} \Delta_q^n\right)^2,
\end{eqnarray}
where $\mu_\chi$ is the reduced mass of the dark matter and nucleon
and $m_n$ is the mass of the nucleon.
$\sigma_\textrm{SI,SD}$ stands for the spin-independent (SI)
and spin-dependent (SD) cross section for the dark matter-nucleon scattering.
The scalar form factor of the nucleon
\begin{equation}
f_n = \sum_{q=u,d,s} f_q^n + \frac{2}{27} \sum_{Q=c,b,t} f_Q^n,
\end{equation}
where $f_Q^n=1-f_u^n-f_d^n-f_s^n$. Here, we use
$f_d^p = 0.017$, $f_u^p = 0.023$, and $f_s^p=0.053$,
for the scalar form factors,
$f_{V_u}^p=2$, $f_{V_d}^p=1$, for the vector form factors, and
$\Delta_{u}^p=0.85$, $\Delta_{d}^p=-0.42$, and $\Delta_{s}^p=-0.08$
for the axial-vector form factors, respectively~\cite{formfactor,formfactor2}.
We note that if the dark matter particles interact only with the charm quarks,
the dark matter-nucleon scattering cross sections are negligible or significantly reduced.

In Fig.~\ref{figuredirect}, we show (a) the spin-independent cross section
and (b) the spin-dependent cross section versus the dark matter mass
with the bounds of $\Lambda$ extracted for each dark matter mass
from Fig.~\ref{figurelam} for the corresponding operators and luminosities
and applied to the appropriate formulas given in
Eqs.~(\ref{dmnucleona}) to  (\ref{dmnucleond}).
The black and gray lines at the lower right corners are 90\% C.L. exclusion limits
of CMDSlite~\cite{cdmslite2} and superCDMS~\cite{supercdms}, respectively.
Also at the lower right corners of these two plots in Fig.~\ref{figuredirect},
the orange dotted (solid) lines are 90\% C.L. exclusion limits
of CMS for $O^D_5$ ($O^D_8$)  from 1 to 10 GeV~\cite{cmsmonojet},
while the cyan dotted (solid) lines are for $O^D_5$ ($O^D_8$) from ATLAS
and the purple dotted (solid) lines are for $O^D_1$ ($O^D_9$) from ATLAS,
respectively~\cite{atlasmonojet}. Since the bounds from CMS and ATLAS
are obtained for the effective operators (\ref{effopfermion}),
they strongly depend on the UV completion of the operators
if other particles in the UV completion model are not heavy enough to be integrated
out at the LHC energy scale.

The red and blue lines correspond to the $O_1^D (O_5^D)$
and $O_8^D (O_9^D)$ operators in the $J/\psi \chi\bar{\chi}$
($\eta_c \chi\bar{\chi}$) production, respectively.
The green line corresponds to the $O_1^R$ operator in the $J/\psi \chi\chi$
production.
The dashed and solid lines represent the cross section corresponding
to ${\cal L}=1$ ${\rm ab}^{-1}$ and $50$ ${\rm ab}^{-1}$ at $B$ factories, respectively.
The bound on the dark matter-nucleon scattering cross section could reach about
$10^{-41}$ cm$^2$ in the spin-independent case and
$10^{-42}$ cm$^2$ in the spin-dependent case for the very light $m_\chi$.
For $m_\chi \sim 1 \, \textrm{GeV}$ and ${\cal L}=50$ ${\rm ab}^{-1}$ at $B$ factories,
the bound could be about $\sim 10^{-38}$ cm$^2$ in the spin-independent case and
about $\sim 10^{-39}$ cm$^2$ in the spin-dependent case.
The bounds obtained from the mono-jet search at the LHC might be stronger
than those from $B$ factories. However, the bounds from the LHC are restricted
to the large mass region of the mediator, but those from $B$ factories
would be applied to broader range of the mediator mass
($\gtrsim \sqrt{s}$).

As shown in Fig.~\ref{figuredirect}, the LHC bounds are more efficient
than the bounds obtained from the dark matter production associated
with a heavy quarkonium at $B$ factories by a few orders of magnitude
for $m_\chi > 1$ GeV. However, this search is meaningful to provide
complementary search in the lepton colliders for the universal couplings
of dark matter.
If the dark matter particle dominantly couples to the charm quark,
the LHC bound would become weaker for $O_{5,8,9}^D$
since the current bounds arise from the light quark interactions
in the $q\bar{q}$ or $qg$ ($q=u,d$) collisions.
Furthermore, the bounds from the direct detection could be negligible
since they do not give any constraints on such operators at leading order.
For the scalar operators $O_1^{D,R}$, the LHC bounds are still powerful,
but the direct detection bounds become weaker by $O(0.01)$.

We note that one should be careful for the comparison of the results of dark matter search
at colliders with those from direct detection.
As shown in Ref.~\cite{monojetsf}, the EFT approach may be
inadequate for the hadron collider dark matter search because the typical energy
scale in the process is not fixed.
In addition, some of the effective operators in Eqs.~(\ref{effopfermion})
and (\ref{effopscalar}) do not respect the SM gauge symmetry,
in particular, $SU(2)_L$ symmetry.
If one imposes the SM gauge symmetry for the effective operators,
the UV complete model, which is responsible for generating these operators,
may not approach the EFT even in the limit of infinite mass of the mediator.
This is true for the scalar operators in Eq.~(\ref{effopfermion})~\cite{monojetsf}.
The contact interaction for the scalar$\times$scalar operator cannot be realized
at the energy scale of the LHC because the SM Higgs boson which is needed to impose
the SM $SU(2)_L$ gauge symmetry is still a propagating mode.
In this case, the comparison in Fig.~\ref{figuredirect} should be
re-interpreted with the analysis in the full theory.
However, at $B$ factories, one can obtain the contact interaction
by integrating out both the heavy degree of freedom of the mediator and SM Higgs boson.
The effective scale $\Lambda$ will depend on the masses of the mediator as well as
the SM Higgs boson. The effective operators thus obtained would respect just
the gauge symmetries of QED and QCD, which is necessary.

\section{Summary}
\label{summary}

In this paper, we investigate light dark matter production associated with
a heavy quarkonium at $B$ factories. For the interaction of dark matter
to the SM sector, the effective field theory approach is adopted, but
it would be straightforward to extend this study to a realistic model
with UV completion.
So far the dark matter search at colliders has been focused on hadron colliders.
In this paper, we showed that the $B$ factories may play a role in searching
for dark matter.
Especially, for the light mediator mass, the $B$ factories would be more
effective.
We took into account the $J/\psi \chi\bar{\chi}$ and $\eta_c\chi\bar{\chi}$
production, but this would easily be extended to the charmonium in
the higher resonances like $\psi(2S)$, $\eta_c(2S)$, and $\chi_{cJ}$,
as well as the bottomonium production. Combining all the results, one may
obtain stronger bounds on the effective energy scale $\Lambda$,
which together with the dark matter mass determines the size of
the dark matter-nucleon scattering cross section for direct detection.

The EFT approach holds only when the mediator mass is larger than
the typical energy scale of the relevant processes.
Therefore the dark matter search via the EFT approach
would be reliable only for $M \gtrsim \sqrt{s}$.
For $M \lesssim \sqrt{s}$, one must consider a UV complete model
or another EFT model including the particles which cannot be integrated out
at the energy scale.
Another important point for the EFT approach is to impose
the full SM gauge symmetry~\cite{monojetsf}.
For example, the $SU(2)_L$ symmetry may require the SM Higgs boson
exchange between the dark sector and SM sector for the scalar$\times$scalar
interaction, which cannot be integrated out at the energy of the LHC.
Searches for dark matter at the $B$ factories are
relatively free from such complications.

The search for light dark matter via its production associated with a heavy quarkonium
at $B$ factories would provide both the alternative and complementary ways
to the hadron colliders and the heavy quarkonium invisible decays.

\acknowledgments
This work was supported in part by the Ministry of Science and Technology (MoST) of Taiwan
under grant number 101-2112-M-001-005-MY3.

\vspace{-1ex}


\section{Introduction}
\let\thefootnote\relax\footnote{This paper will be presented in IEEE PES General Meeting 2016.}
Frequency regulation service involves the injection or withdrawal of active power from the power grid to maintain the system frequency \cite{dvorkin_2014}. In the United States, frequency regulation is equivalent to secondary frequency control, while primary frequency control is more commonly known as frequency response.
The most common mean for a unit to provide frequency regulation is following a system operator's automatic generation control (AGC) signal, which computes the area control error (ACE) from frequency deviations and interchange power imbalances.

An appropriate response to the AGC signal requires dispatching flexibility within the assigned regulation capacity. Traditionally, the majority of frequency regulation capability has been provided by specially equipped generators. As technologies evolve, new types of flexibility resources emerge, such as battery and flywheel energy storage. These energy storage systems have significantly faster ramping capabilities compared to conventional generators, and their participation in frequency regulation can reduce the need to procure regulation capacity ~\cite{PNNL_Regulation}. To create an incentive for fast responding units to participate in frequency regulation, the Federal Energy Regulatory Commission (FERC) enacted FERC Order 755~\cite{FERC_755}, which requires system operators to add a performance payment with an accuracy adjustment to the capacity payment typically used in markets for ancillary services. This ``pay-for-performance" scheme has been implemented by most U.S. independent system operators (ISO) and regional transmission organizations (RTO). Under this new market rule, fast ramping regulation resources receive higher regulation payments because their performance is usually more accurate. A case study conducted by RES Americas has already shown that its 4~MW~/~2.6~MWh battery storage achieved a very high performance in following the regulation signal provided by PJM~\cite{RES_ES}.

To further utilize the fast responsive capability of energy storage systems (ESS) beyond the traditional AGC framework, some system operators, such as PJM, have introduced fast regulation. Units participating in fast regulation follow a regulation signal that changes much faster than the traditional AGC signal, and receive extra payments for doing so. The system operator also benefits from the improved regulation accuracy. A review by the California Energy Storage Alliance (CESA) shows that the combination of pay-for-performance and fast regulation has reduced the procurement requirement in the PJM regulation market by 30\%~\cite{CESA_Reg}.

However, the high investment cost of ESS has so far limited the installed energy storage capacity \cite{pandzic_2015}. The energy lost during charging and discharging and the energy imbalance in the regulation signal make it difficult for ESS operators to provide regulation services seamlessly over a relatively long time period (e.g. several hours). Some system operators have elected to lower this barrier by co-optimizing energy storage when dispatching regulation, i.e. taking into account the state of charge (SoC) and the cycling efficiency of ESS.

This paper compares market policies that affect energy storage in the frequency regulation market. It covers the ISOs and RTOs that have organized wholesale electricity markets, including the PJM Interconnection (PJM); the New York Independent System Operator (NYISO); the Midcontinent Independent System Operator (MISO); ISO New England (ISO-NE); California Independent System Operator (CAISO). Study results on pay-for-performance, fast regulation, and storage energy compensation are discussed in Section~II, Section~III, and Section~IV, respectively. Section~V analyzes the regulation payment in different markets. Section~VI summarizes the findings.

\section{Pay-for-Performance}

All ISO~/~RTOs in our study have complied with FERC Order 755 and implemented  pay-for-performance in their regulation markets, except for NYISO and MISO which pay for regulation capacity only in their  day-ahead regulation markets (DAM).
A typical pay-for-performance regulation market involves the following terms:

\begin{itemize}
	\item \emph{Regulation capacity ($C$)} is the maximum and minimum regulation power that a unit offers during a given period. The regulation capacity is usually symmetrical so that a unit provides the same positive and negative capacity (i.e. $\pm 5$ MW). With the exception that CAISO separates regulation into up and down components. The capacity price in $\mathrm{\$/MW\cdot h}$, represents the option price to reserve 1~MW of capacity to provide regulation in the future. Not to be confused with the exercise price of the energy required to follow the regulation signal
	\item \emph{Regulation mileage ($M$)} is the sum of the absolute values of the regulation control signal movements. The mileage of a regulation resource output $\{P_1,\dotsc,P_n\}$ over a time period with $n$ signal steps and a regulation capacity $P_{max}$ is
	\begin{align}\label{Eq:mileage}
	M=\sum_{i=1}^{n}|P_i-P_{i-1}|/P_{max
	\end{align}

	 The unit for mileage is $\mathrm{\Delta MW/MW}$, representing the mileage of 1~MW of regulation capacity. The unit for mileage price is $\mathrm{\$/\Delta MW}$.
	\item \emph{Performance factor ($\rho$)} is a score from 0 to 1 that indicates a unit's performance in following the regulation signal. Most ISO/RTOs model the performance factor as a unit's accuracy in signal following (integration of absolute errors). PJM uses a more complex measure of performance that involves delay, correlation and accuracy~\cite{PJM_M12}.
\end{itemize}
Table~\ref{tab:pfp} shows the names used for these quantities in each market.

\begin{table}[t]
	\begin{center}

		\centering

		\caption{Pay-for-performance market terms}
		\label{tab:pfp}

		\begin{tabular}{|c|c|c|c|c|}
			\cline{2-5}
			\multicolumn{1}{c|}{} & \multicolumn{2}{|c|}{Price}  & Regulation  & Performance  \Tstrut\\
			\multicolumn{1}{c|}{} &  Capacity &  Mileage & mileage &  factor \Bstrut\\
			\hline

			PJM 	& capability 	& performance 	& mileage 		& performance \Tstrut\\
			\cite{PJM_M11, PJM_M28}		& price 		& price			& 				& score\Bstrut\\
			\hline
			NYISO 	& capacity 		& movement 		& movement 	& performance\Tstrut\\
			\cite{NYISO_Bill}		& price		 	& price 		& 			& index\Bstrut\\
			\hline
			MISO 	& capacity 		& mileage 		& mileage 	& performance\Tstrut\\
			\cite{MISO_M5, chen2015development}		& price			& price		 	& 			& accuracy\Bstrut\\
			\hline
			ISO-NE  & capacity 		& service 		& movement 	& performance\Tstrut\\
			\cite{ISONE_Rule1}		& price			& price 		& 			& score\Bstrut\\
			\hline
			CAISO 	& capacity 		& mileage 		& mileage 	& accuracy\Tstrut\\
			\cite{CAISO_755}		& price			& price 		& 			& percentage\Bstrut\\
			\hline

		\end{tabular}

	\end{center}
	\vspace{-.5cm}
\end{table}

In a pay-for-performance market, participants offer a capacity price bid and a mileage price bid. The system operator credits providers of regulation using the market capacity clearing price ($\lambda_C$), the mileage clearing price ($\lambda_M$), and the performance factor of its unit during the procurement period. A unit must reach a minimum performance score to be qualified for receiving credits, and maintain this minimum performance over a certain period to be eligible to bid into the regulation market in the future. We generalize the method for awarding regulation credit ($\lambda_R$) as follows:
\begin{align}\label{Eq:reg_credit}
\lambda_R = C(\lambda_C+\rho M\lambda_M
\end{align}
Eq.~\ref{Eq:reg_credit} applies to CAISO and the real-time regulation market in NYISO. Eq.~\ref{Eq:reg_credit} also applies to MISO when a unit's performance is below 70\%, a unit receives full credit when above 70\% during each 5-minute performance test interval. At ISO-NE, the performance factor is used to adjust the capacity payment, thus Eq.~\ref{Eq:reg_credit} becomes
\begin{align}\label{Eq:reg_creditA}
\lambda_R' = \rho C(\lambda_C+M\lambda_M
\end{align}
PJM uses mileage ratio in its credit calculations. The mileage ratio is the ratio between the mileage of the assigned regulation signal and the mileage of the traditional regulation signal ($M_{RegA}$, see Section~\ref{Sec:FR:PJM}). The formula for calculating the credit becomes
\begin{align}\label{Eq:reg_creditB}
\lambda_R'' = \rho C(\lambda_C+\frac{M}{M_{RegA}}\lambda_M
\end{align}

\section{Fast Regulation}\label{Sec:FR}

A fast regulation signal is generated by applying a  high-pass filter to the AGC signal. Such signals have a fast ramping rate but are designed to have a zero-mean. They are therefore ideal for ESS interested in providing frequency regulation. Fast regulations have been implemented by PJM and ISO-NE, while MISO provides an alternative approach for utilizing fast responsive units.

A regulation signal's energy capacity requirement refers to the minimum energy capacity needed to follow perfectly the regulation signal, neglecting the efficiency and SoC constraints of an ESS.
The energy requirement of the signal in that hour is the difference between the maximum and minimum SoC during the operating hour.
For a regulation signal $\mathbf{S}_h$ at hour $h$ with a step size of $T_s$ and $N$ steps per hour, the hourly energy requirement $E_{h}$ is therefore
\begin{align}\label{Eq:energy_offset}
E_{h} = \max_{j=1\dotsc N}\left(\sum_{i=1}^{j}T_sS_i\right)-\min_{j=1\dotsc N}\left(\sum_{i=1}^{j}T_sS_i\right
\end{align}

The energy balance of a regulation signal is the initial SoC level that allows a storage with the minimum energy capacity ($E_h$) to follow the signal perfectly, ignoring its efficiency and SoC constraints. A perfectly balanced signal (at 50\%) allows the SoC level of the energy storage to return to its initial level at the end of the operating period. Since Eq.~\ref{Eq:energy_offset} assumes that the initial SoC is zero (i.e. it simulates an ESS with negative capacity), the energy balance $\sigma_{h}$ for the signal $\mathbf{S}_h$ is
\begin{align}
\sigma_h = -\frac{1}{E_h}\min_{j=1\dotsc N}\left(\sum_{i=1}^{j}T_sS_i\right
\end{align}

\begin{figure}[t]%
\centering
\vspace{-.3cm}
\subfloat[Regulation energy hourly requirement (1MW capacity).]{
	\includegraphics[trim = 05mm 5mm 10mm 0mm, clip,  width = .9\columnwidth]{energy_offset.eps}
	\label{Fig:energy_offset}%
}
\vspace{-.3cm}
\subfloat[Regulation energy hourly balance.]{

	\includegraphics[trim = 05mm 5mm 10mm 0mm, clip,  width = .9\columnwidth]{energy_balancing.eps}
	\label{Fig:energy_balancing}%
}
\caption{Comparison of PJM and ISO-NE regulation signals.}%
\label{Fig:fast_reg}%
\vspace{-.4cm}
\end{figure}

\subsection{PJM RegD}\label{Sec:FR:PJM}

PJM provides two regulation signals in its day-ahead regulation market: RegA, is a low-pass filtered ACE designed for traditional regulation units, while RegD, is a high-pass filtered ACE designed for fast regulation units~\cite{PJM_M12}. The RegD signal is controlled to have zero-mean within 15 minutes~\cite{PJM_RegHist}. Units can choose freely which signal they follow, but can only follow one signal at a time. Our analysis of these regulation signals from June 2013 to May 2014 shows that the mileage of RegA is about 5~$\mathrm{\Delta MW/MW}$ per hour, while the mileage of RegD is around $15~\mathrm{\Delta MW/MW}$ per hour, suggesting that RegD ramps much faster. In market settlements, the mileage ratio for RegA is therefore 1, and the mileage ratio for RegD is around 3. Section~\ref{Sec:Reg_Pro} provides a detailed comparison of regulation payments. The box plot in Fig.~\ref{Fig:energy_offset} summarizes the distribution of the hourly energy requirement and shows that in most cases PJM RegD has an energy requirement of less than 0.25~MWh or approximately 15 minutes of full charge/discharge. Fig.~\ref{Fig:energy_balancing} shows that RegD has an average energy balance of 50\%, but that the extreme cases can reach 0 and 100\%.

\subsection{ISO-NE energy-neutral}\label{Sec:FF:ISONE}

ISO-NE includes ESS in the Alternative Technology Regulating Resources (ATRR), and provides two fast regulation signals (energy-neutral continuous~(ENC), energy-neutral trinary~(ENT)), and a conventional regulation signal~\cite{ISONE_EN_AGC}. A resource can change once a month the signal that it follows. ENC is identical to PJM's RegD signal, while ENT is specifically designed for ATRRs and the dispatch is either zero, full charge power, or full discharge power.
Since each trinary dispatched resource is assumed to have a step response and hence equal
participation factors, the ENT algorithm dispatches the entire set of Trinary Dispatched ATRRs together.
Fig.~\ref{Fig:fast_reg} shows the energy requirement and the balance analysis using a one-day duration simulated signal from ISO-NE~\cite{ISONE_AGC} and suggest that these signals behave in a similar fashion as PJM's RegD and RegA . The mileage difference between the three simulated signals is not significant: the average mileages being 18, 16, and 14~$\mathrm{\Delta MW/MW}$ for ENC, ENT, and the conventional signal, respectively.

\subsection{MISO AGC-Enhancement}

In the MISO reserve market, the deployment interval for ESS (which are called stored energy resources or SER) is optimized for each ESS based on its initial condition~\cite{MISO_M2}. Taking things one step further, MISO is proposing a new regulation reserve deployment scheme called AGC-Enhancement~\cite{MISO_AGC}. The idea of this scheme is to deploy and un-deploy the fast ramping units first in the AGC dispatch, by creating a new fast ramping deployment group. MISO expects AGC-Enhancement to improve the utilization of fast ramping units.

\section{Storage Energy Offset Compensation}

Due to their less than perfect efficiency and the fact that the AGC signal has a non-zero mean, the SoC of an ESS is not maintained in the traditional regulation framework. While the fast frequency regulation alleviates  the effect of the signal energy imbalance, some ISOs have decided to directly compensate the energy offset of storage units in their regulation market through  dispatching. Such mechanisms include the Regulation Energy Management~(REM) of CAISO, and the Real-Time Dispatching~(RTD) of NYSIO.

\subsection{CAISO REM~\cite{CAISO_REM}}

CAISO treats ESS as Non-Generating Resources~(NGR), a set of resources that also includes responsive demands. NGR units have two options for participating in the regulation market: non-REM~(traditional), and REM. Non-REM units are subjected to the same requirements as traditional regulation units, and must meet the 60-minute continuous energy requirement (i.e., to provide 1~MW regulation capacity, a non-REM NGR unit must be able to deliver at least 1~MWh). The continuous energy requirement for NGR-REM units is 15 minutes. NGR-REM units can only participate in the CAISO regulation market, and the subscription to REM can only be changed monthly.

CAISO compensates the energy offset of NGR-REM units. By participating in REM, a NGR unit submits a preferred SoC set-point. CAISO will maintain this set-point by dispatching energy from the Real Time Market (RTM) for the next real time dispatch interval~\cite{CAISO_REMNGR}. This balancing energy is dispatched to each NGR-REM units together with the AGC signal by the Energy Management System (EMS). The EMS also takes into account a unit's efficiency in the dispatch. The REM energy compensation is guaranteed during normal operations, while in emergency cases the SoC is restored later.

In the CAISO regulation market settlement, the regulation energy charged and discharged from the storage and the dispatched REM energy are settled at the RTM locational marginal price (LMP), in addition to the capacity settlement and mileage settlement. Since REM guarantees to maintain the SoC of NGRs, the consequence of the energy settlement is that NGR-REM units must pay for their energy losses at their LMP. This settlement method is currently unique to CAISO.

\begin{figure}[t]
	\centering
	\vspace{-.5cm}
	\includegraphics[trim = 18mm 0mm 10mm 0mm, clip, width=.95\columnwidth]{AS_price_DA.eps}
	\caption{DAM regulation price comparison (outliers omitted).}
	\label{Fig:AS_price}
	\vspace{-5mm}
\end{figure}

\subsection{NYISO RTD~\cite{NYISO_AS}, \cite{NYISO_LESR}}

NYISO calls ESS Limited Energy Storage Resources (LESR). LESR can bid into the regulation DAM and RTM, but only receive mileage payments in the RTM settlement. In the day-ahead regulation market, LESR units are evaluated and scheduled on an hourly basis without considering potential energy limitations. In the real-time regulation market, the RTD manages the energy level of LESR to maintain regulation capability to the extent possible by charging and discharging the LESR as necessary. The RTD assigns regulation base points to an LESR every  5 minutes based solely on its SoC level. For example, at each 5-minute interval, a 20~MW regulation capacity LESR  will provide $\pm$20~MW regulation reserve at 0~MW base point when its SoC is within a dead-band centered at 50\% SoC. When its SoC is above this dead-band, it will provide regulation with reduced capacity and a positive (charging) base point; when its SoC is below this dead-band, it will provide regulation with reduced capacity and a negative (charging) base point. Thus the SoC of the LESR converges to 50\%.

As in CAISO, the RTD SoC management is suspended during emergencies, when it is called RTD-CAM (Corrective Action Mode). During RTD-CAM, all LESR providing regulation discharge to zero in order to make up for the loss of generation. The SoC level is restored after RTD-CAM terminates.

\section{Regulation Market Income}\label{Sec:Reg_Pro}

\subsection{Day-ahead Market Results}

This section compares regulation payments in different regulation markets using public regulation pricing data from each ISO/RTO. In pay-for-performance markets, we assume that the performance factor is an ideal 100\%, and use Eqs.~\ref{Eq:reg_credit}-\ref{Eq:reg_creditB} to calculate these payments. The box plot in Fig.~\ref{Fig:AS_price} shows the distribution of calculated regulation payments in the DAM for each ISO/RTO.

\subsubsection{PJM}
The payments for the two regulation products available in PJM are based on signals and market clearing prices from June~2013 to May~2014. The hourly regulation mileage for each signal is calculated using Eq.~\ref{Eq:mileage}.

\subsubsection{CAISO}

CAISO publishes all its historical market data on the CAISO Open Access Same-Time Information System (OASIS). Besides pricing data, CAISO also publishes the hourly reg-up and reg-down mileage data (system mileage multipliers) of its AGC signal. To be consistent with the calculations for the other markets, we assume that each unit offers the same regulation up and regulation down capacity, we calculate the payments separately using day-ahead market prices from June~2014 to May~2015, and add the two results together. The average calculated reg-up price is \$5.2/MWh, while the reg-down price is \$4.1/MWh.

\subsubsection{MISO}

MISO publishes ancillary service pricing data on its market report website~\cite{MISO_data}. We use the day-ahead ancillary service market prices from June~2014 to May~2015. This market does not have a mileage price.

\subsubsection{ISO-NE}

ISO-NE publishes ancillary service pricing data on its market data and information website~\cite{ISONE_data}. We use the final hourly clearing price data from June~2014 to May~2015. ISO-NE started mileage bidding on March 31st, 2015, so the clearing price for April and May contains a capacity price and a service (mileage) price. For these two months, we use mileage for the simulated ENT signal described in Section~\ref{Sec:FF:ISONE}, the payment for the other two signals are not calculated since the mileages are very close to those for ENT.

\subsubsection{NYISO}

NYISO publishes its ancillary service pricing data on its pricing data website~\cite{NYISO_data}. We use the DAM regulation clearing prices from June~2014 to May~2015. This market does not have a mileage price.

\subsection{Real-time Market}

CAISO, MISO, and NYISO have both RTM and DAM for operating reserves. The RTM is to procure additional reserves according to real-time system co-optimizations with updated system conditions~\cite{SANDIA_AS}. All three markets have mileage payments, however most mileage clearing prices are very close to zero and thus have little effect on the final price.
\subsubsection{CAISO}
CAISO determines its RTM regulation capacity in the hour-ahead scheduling process.
The weighted average RTM price for reg-up and reg-down is \$7.84/MWh and \$4.28/MWh, while for the RTM the prices are \$6.08/MWh and \$3.98/MWh, respectively~\cite{CAISO_Report2014}.
\subsubsection{MISO}
In MISO the RTM is cleared every five minutes to produce price and dispatch levels for regulation. Its annual average hourly RTM regulation prices are in general 10\% to 15\% higher than the DAM prices~\cite{MISO_Report2013}.

\subsubsection{NYISO}
NYISO calculates the RTM requirement over an RTD interval (1 hour) every 5~minutes. The average DAM regulation price is \$10.1/MWh in 2013 and \$12.9/MWh in 2014, while the average RTM regulation price is \$9.8/MWh in 2013 and \$13.8/MWh in 2014~\cite{NYISO_Report2014}.

\begin{table*}[t]
    \vspace{-12mm}
	\begin{center}

		\centering
		\footnotesize
		\caption{A Comparison of Regulation Market Policies on Energy Storage and Fast Responsive Units}
		\label{tab:rating}

			\begin{tabular}{|c|c|c|c|c|c|c|c|c|c|}
\cline{2-9}
\multicolumn{1}{c|}{} & \multicolumn{2}{c|}{AS Requirements} &	ES Offset & Fast Reg. & ES & Reg. Cap. & Reg. Cost & Gen. Avg. \Tstrut\\
\multicolumn{1}{c|}{} & \multicolumn{2}{c|}{1~MW Cap.(MWh)} & Control & Schemes &  Abbr. & (MW) & (M\$) & (GW) \Bstrut\\
\hline

CAISO \cite{CAISO_REM, CAISO_MO, CAISO_Report2014} & 0.25 (REM) & 1 (else) & REM & no &  NGR & $\sim325$ & 32  & $\sim28$ \Tstrut\Bstrut\\
\hline
PJM \cite{PJM_M11, PJM_M12, PJM_Report2014}& 0.25 (RegD) & 1 (RegA) & no & RegD &  - & $525-700$ & 254 & $\sim171$ \Bstrut\Tstrut\\
\hline
ISO-NE \cite{ISONE_Rule1, ISONE_OP14, ISONE_Report2014} & 0.25 (EN) & 1 (else)  & no & energy-neutral &  ATRR & $30-130$ &  29 & $\sim 12$ \Bstrut\Tstrut\\
\hline
MISO \cite{MISO_AGC, MISO_M2, MISO_Report2014} &  \multicolumn{2}{c}{-}  & \multicolumn{2}{|c|}{AGC-Enhancement}  & SER & $\sim 400$ & 20 & $\sim55$ \Bstrut\Tstrut\\
\hline
NYISO \cite{NYISO_AS, NYISO_Report2014} & \multicolumn{2}{|c|}{-} & RTD & no &  LESR & $150-250$ & 29 & $\sim 18$ \Bstrut\Tstrut\\
\hline

			\end{tabular}

	\end{center}
	\vspace{-8mm}
\end{table*}

	\vspace{-1mm}
\section{Conclusions}\label{Sec:Con}
	\vspace{-1mm}
Table~\ref{tab:rating} summarizes the findings of our study. In the column ancillary service (AS) requirements, we list the minimum energy capacity required to provide 1 MW of regulation reserve. This requirement is specified for CAISO, while for PJM and ISO-NE, it is derived from the simulations described in Section.~\ref{Sec:FR}. We were not able to find similar information for MISO and NYISO. The last three columns show the average required regulation capacity, the regulation cost, and the average generation capacity in each operator zone.

The comparison shows that PJM operates a very efficient regulation market in terms of the ratio between regulation capacity and generation, while the cost of regulation is similar to the other ISOs relative to the size of the interconnection. In particular, if we compare the ancillary market cost of PJM with CAISO, the historical cost of regulation in PJM per MWh load is $\sim$\$0.2-0.3/MWh, while this cost is $\sim$\$0.1-0.2/MWh for CAISO. On the other hand, the total per MWh ancillary service cost is much higher in PJM, at $\sim$\$2/MWh, than in CAISO, where it is $\sim$\$0.3/MWh~\cite{PJM_Report2014, CAISO_Report2014}. A possible explanation for this difference is that CAISO performs a co-optimization of energy and ancillary services, while PJM has a coupled co-optimization process with more complicated categories of ancillary services ~\cite{SANDIA_AS}. This latter approach is less efficient than CAISO's and results in a higher ancillary service cost. For example, sometimes the PJM RegD capacity is over-saturated, meaning that the capacity ratio of RegD is too high. Since the RegD logic is to maintain a 5-minute zero-mean energy, a high RegD capacity ratio can make the controller go against ACE correlation for extreme cases~\cite{PJM_RegHist}.

Our comparison of different markets shows that the regulation markets are becoming more suitable and profitable for ESS because the new pay-for-performance policy makes energy storage more competitive. New zero-mean  regulation policies  have lowered the barrier for ESS to provide regulation for longer periods. From the perspective of a system operator, the participation of energy storage under appropriate market designs makes the regulation more efficient. However, providing a zero-mean regulation signal can be challenging, and must be well coordinated with other ancillary service products so that the energy balancing signal does not perturb the system.

\vspace{-2mm}
\section*{Acknowledgment}
\vspace{-1mm}
The authors thank Dr. I. Gyuk and his colleagues at the US DOE Energy Storage Program for funding this research. Sandia National Laboratories is a multi-program laboratory managed and operated by Sandia Corporation, a wholly owned subsidiary of Lockheed Martin Corporation, for the U.S. Department of Energy's National Nuclear Security Administration under Contract DE-AC04-94AL85000.
\vspace{-2mm}
\bibliographystyle{IEEEtran}