\section{Introduction}
Despite strong observational evidence for non-baryonic dark matter (DM) which interacts gravitationally~\cite{Kolb:1990}, the
detection of dark matter interactions with the Standard Model (SM) remains elusive. Many extensions of the SM predict dark matter candidates
which should leave signatures in the direct and indirect dark matter detection experiments, and at hadron colliders.
If the particles mediating dark matter-Standard Model interactions are much heavier than the energy scales involved, then the
constraints on dark matter interactions arising from these disparate detection strategies can be related
to each other in a model independent fashion via a generalized effective field theory (EFT) framework,
in which the details of the ultraviolet (UV) physics have been integrated out of the
Lagrangian~\cite{Goodman:2011,Hooper:2009,Cao:2011,Kumar:2013iva,Rajaraman:2013,Dreiner:2013,Buckley:2013,Rajaraman:2013-2,
Busoni:2014,Buchmueller:2014,Crivellin:2014,Busoni:2014-2,Alves:2014,Fedderke:2014},
 and dark matter-Standard Model interactions occur through contact operators.

Weakly interacting massive particles (WIMPs) are stable dark matter candidates predicted by many models of
physics beyond the SM~\cite{Kolb:1990}. Although
recent hints of possible WIMP signals may be encouraging, the lack of clear and convincing evidence for the discovery of WIMP dark matter
motivates consideration of dark matter candidates which deviate from the expectations of the WIMP paradigm.  A well-motivated
example is light dark matter (LDM), a class of dark matter candidates with masses typically $\sim 10 \, {\rm MeV} - 10 \, {\rm GeV}$.  LDM would
elastically scatter at direct detection experiments, with nuclear recoil energies which are relatively small and may be below the experimental
threshold, rendering them undetectable.  In this case, other experimental means aside from direct detection would be required to probe
DM-SM interactions.  For example, complimentary bounds on LDM scattering can be inferred from collider monojet
searches~\cite{Feng:2005gj,Goodman:2010,Birkedal:2004,Hooper:2010,Bai:2010,Rajaraman:2011,Fox:2012,Bai:2013,Agrawal:2013,
Goodman:2011-2,Bai:2011,Papucci:2014}.
Nonresonant LDM production at low energy $e^+/e^-$colliders can also be used to set model independent limits on electron scattering and, if
the LDM couples universally, nucleon scattering~\cite{Essig:2013}.

In this work, we consider the prospects for probing LDM-quark interactions through bounds on invisible decays of heavy quarkonium states at
colliders~\cite{Fayet:2007,Fayet:2010,Cotta:2013,Schmidt-Hoberg:2013hba,McElrath:2007}.  Such bounds have already been considered in a variety of contexts, including
$B$ and $D$ meson decays~\cite{Badin:2010}, and $\Upsilon$ decays into scalar LDM~\cite{Yeghiyan:2010,McKeen:2009rm}.  However, the constraints which one
can obtain on dark matter-quark
interactions depend in detail on the quantum numbers of the heavy meson, as well on the choice of final state
(i.e., $\rightarrow invisible$ or $\rightarrow \gamma +invisible$).
The angular momentum and $C/P$ transformation properties of the initial state (as well as the presence or absence of a photon in the
final state) together determine
which of the possible dark matter-quark interaction structures can participate in the decay process (and can thus be bounded by constraints
on invisible decays).
We consider invisible decays of the heavy quarkonium states $\Upsilon (1S)$ and $J / \Psi $, mesons with $J^{PC} = 1^{--}$.
As the quark constituents annihilate in an $s$-wave, the dependence of the meson decay matrix element on
the associated nonrelativistic bound state wavefunction is very simple and can be determined experimentally, with relatively little uncertainty.
Moreover, since the mesons which we consider each have a quark and anti-quark of the same flavor, the DM-SM interactions which
we introduce are not constrained by bounds on flavor-violation and, given universal quark coupling, can contribute to nucleon scattering.

At quark level, the matrix element relevant for meson decay ($\bar q q \rightarrow \bar X X$) is also relevant for
monojet/photon/$W$,$Z$ searches at the LHC~\cite{CMS:2012, CMS:2012tea,ATLAS:2012ky,ATLAS:2012,ATLAS:2013}.  As a result, these searches
will share many features, and event rates will have the same
dependence on the energy of the process.  Moreover, bounds arising from both invisible meson decay rates and LHC monojet
searches do not weaken as the dark matter mass decreases, in notable contrast to direct detection searches.
However, if the particle mediating LDM interactions is light, then the interaction is poorly approximated by a contact
operator for the purposes of LHC mono-anything searches,
and model independent bounds on the interaction strength can no longer be obtained. Even given a mediating particle massive
enough to warrant the use of the contact interaction approximation, the masses of LDM particles would be difficult to resolve at LHC searches due to
the large center-of-mass energy of the beam.
The bounds obtained from meson decays thus provide a unique handle on some dark matter interaction models, which is complementary
to the information provided by other search strategies.

In this paper, we use the limits on bound state decay widths to constrain the coupling of scalar, fermion or vector LDM to Standard Model quarks through
all contact operators of dimension six or lower.
In section II, we review the relevant effective contact interactions and calculate the resulting meson
decay rates and dark matter annihilation and scattering cross sections.
In section III we present constraints on all of the relevant dark matter-quark interaction structures
arising from $\Upsilon (1S)$ decay and Fermi gamma-ray searches
of dwarf spheroidal galaxies, and relate these constraints to those arising from direct detection experiments and LHC searches.
We conclude in section IV with a discussion of our results.

\section{Framework and Constraints}

We will consider a framework in which dark matter-quark interactions can be parametrized by a four-point contact effective
operator.  Such a structure can be written as an appropriate Lorentz contraction of a Standard Model quark bilinear and a
dark matter bilinear.  We are only interested in effective operators which can yield a non-zero matrix element when acting
on a $1^{--}$ meson state, such as the $\Upsilon (1S)$ or $J/\psi$; for such operators, the quark bilinear must be either
$ \bar q \gamma^i q$ or $ \bar q \sigma^{0 i} q$, where $i$ is a spatial index~\cite{Kumar:2013iva}.
The angular momentum quantum numbers of the quark/anti-quark bound states of interest are $S=1$, $L=0$, $J=1$.
In general, the final state need not have the same $C$ and $P$ transformation properties as the initial state (and thus need not have
the same $S$ and $L$ quantum numbers), but must have the
same total angular momentum $J$.  The effective operators\footnote{Henceforth, we refer to spin-0, spin-1/2 and spin-1 dark matter
fields with $\phi$, $X$ and $B^\mu$, respectively.  When describing a dark matter particle of arbitrary spin, we will use $X$.}
of dimension 6 or less which can have a non-zero matrix element with either an $\Upsilon (1S)$ or
$J/\psi$ initial state are listed in Table~\ref{tab:structs}~\cite{Kumar:2013iva} (the operators are labeled using the conventions
of~\cite{Kumar:2013iva} and~\cite{Goodman:2010}).
Note that for Majorana fermion dark matter, F6 is the only non-vanishing contact operator.

\begin{table}[hear]
\centering
\begin{tabular}{|c|c|c|c|c|}
\hline
 Name & Interaction Structure &  Annihilation  & Scattering \\
\hline
 F5  & $ (1/\Lambda^2) \bar X \gamma^\mu X \bar q \gamma_\mu q$ & Yes & SI \\
\hline
 F6  & $ (1/\Lambda^2) \bar X \gamma^\mu \gamma^5 X \bar q \gamma_\mu q$ & No & No \\
\hline
 F9  & $ (1/\Lambda^2) \bar X \sigma^{\mu \nu} X \bar q \sigma_{\mu \nu} q$ & Yes & SD \\
\hline
 F10  & $ (1/\Lambda^2) \bar X \sigma^{\mu \nu} \gamma^5 X \bar q \sigma_{\mu \nu} q$ & Yes & No \\
\hline
 S3  & $ (1/\Lambda^2) \imath Im (\phi^\dagger\partial _\mu \phi ) \bar q \gamma^\mu q$ & No & SI \\
\hline
 V3  & $ (1/\Lambda^2 ) \imath Im (B_\nu^\dagger\partial _\mu B^\nu ) \bar q \gamma^\mu q$ & No & SI \\
\hline
 V5  & $ (1/\Lambda ) (B_\mu^\dagger B_\nu - B_\nu^\dagger B_\mu ) \bar q \sigma^{\mu \nu} q$ & Yes & SD \\
\hline
 V6  & $ (1/\Lambda ) (B_\mu^\dagger B_\nu - B_\nu^\dagger B_\mu ) \bar q \sigma^{\mu \nu} \gamma^5 q$ & Yes & No \\
\hline
 V7  & $ (1/\Lambda^2) B_\nu^{(\dagger)} \partial^\nu B_\mu \bar q \gamma^\mu q$ & No & No \\
\hline
 V9  & $ (1/\Lambda^2) \epsilon^{ \mu \nu \rho \sigma} B_\nu^{(\dagger)} \partial_\rho B_\sigma \bar q \gamma_\mu q$ & No & No \\
\hline
\end{tabular}
\caption{Effective contact operators which can mediate the decay of a $J^{PC} = 1^{--}$
quarkonium bound state.
We also indicate if the operator can permit an $s$-wave dark matter initial state to annihilate to a quark/anti-quark
pair; if so, then a bound can also be set by indirect observations of photons originating from
dwarf spheroidal galaxies. Lastly, we indicate if the effective operator can mediate velocity-independent nucleon scattering
which is either spin-independent (SI) or spin-dependent (SD).}
\label{tab:structs}
\end{table}

The measured limits on the invisible
decay branching fractions could possibly be contaminated by decays to an invisible final state as well as a soft photon.  However,
contributions to the processes $\Upsilon (1S), J/\psi \rightarrow \bar X X \gamma$ would not arise from the operators we consider.  A contribution could
arise from the other dark matter-quark four-point operators which we do not consider, but the rate would be suppressed by an additional
factor of $\alpha$; to constrain such operators, it would be more fruitful to search directly for $\gamma +invisible $ decays.
The contribution of soft photon processes to fake decays was also considered in~\cite{Aditya:2012im}, where
the contribution was also small.

The meson decay rate matrix element can be determined by convolving the quark/antiquark annihilation ($\bar q q \rightarrow \bar X X$)
matrix element with the meson bound state wavefunction.
Using crossing symmetry, one can relate the quark annihilation matrix element to the matrix elements for
either dark matter annihilation ($\bar X X \rightarrow \bar q q$) or dark matter-nucleon scattering ($XN \rightarrow XN$).

The signals which can be observed at various experiments
depend on which quarks appear in the effective operator.  The $b$- and $c$-quark couplings will be relevant for the decay of
$\Upsilon$ or $J/\psi$ states, respectively,  whereas $u$-, $d$- and $s$-quark couplings are the most relevant for either dark matter
annihilation, dark matter-nucleon scattering, or dark matter production at the LHC.

\subsection{Bound State Decays}

The matrix element for the decay of a bound state is given by the convolution of the nonrelativistic bound state wavefunction with the
annihilation matrix element for a free quark/anti-quark pair.  Since we consider $s$-wave meson bound states in the nonrelativistic
approximation, this convolution depends only
on the value of the spatial wavefunction at the origin, $\psi (0)$.  The wavefunction at the origin can then be determined
from the well-measured decay branching fraction to $e^+ e^-$, yielding
\begin{eqnarray}
{\cal B} ( \Upsilon (1S) \rightarrow e^{+} e^{-} ) = 16 \pi\alpha^2 Q_b^2 { | \psi_\Upsilon (0) |^2 \over \Gamma_\Upsilon M_\Upsilon^2 }
= 0.0238 \pm 0.0011 ,
\nonumber\\
{\cal B} ( J / \Psi  \rightarrow e^{+} e^{-} ) = 16 \pi \alpha^2 Q_c^2 { | \psi_{J/\psi} (0) |^2 \over \Gamma_{J/\psi} M_{J/\psi}^2 }
= 0.0594 \pm 0.0006 ,
\end{eqnarray}
with $M_\Upsilon = 9460.30 \pm 0.26 \, {\rm MeV}$, $\Gamma_\Upsilon = 54.02 \pm 1.25 \, {\rm keV}$, $M_{J/\psi} = 3096.916 \pm 0.011 \, {\rm MeV}$
and $\Gamma_{J/\psi} = 92.9 \pm 2.8 \, {\rm keV}$~\cite{PDG}. Note, we have ignored the contribution from $Z$,$h$-exchange; this
contribution is smaller than the uncertainties
in the measured branching fraction.
Searches for $\Upsilon$(1S) invisible decays have been performed by Belle~\cite{Belle:2007} and BaBar~\cite{BaBar:2009}
operating at the $\Upsilon (3S)$ resonance. They use the transition $\Upsilon(3S) \rightarrow\pi^{+}\pi^{-}\Upsilon (1S)$ to detect invisible
$\Upsilon (1S)$ decays and reconstruct the presence of the $\Upsilon$(1S) from the $\Upsilon(1S)$ peak in the recoil mass distribution, $M_{rec}$,
by tagging $\pi^{+}\pi^{-}$ pairs with kinematics
\begin{eqnarray}
M_{rec}^{2} \equiv s + M_{\pi \pi }^{2}-2\sqrt{s}E^{*}_{\pi \pi },
\end{eqnarray}
where $M_{\pi \pi}$ is the invariant mass of the dipion system, $E^{*}_{\pi \pi}$ is the energy of the dipion system in the
center-of-mass (CM) frame of the $\Upsilon(3S)$, and $\sqrt{s} = 10.3552 \, {\rm GeV} $ is $\Upsilon (3S)$ resonance energy.
Similar searches for invisible decays of $J/\Psi$ are based on the transition $\Psi(2S)\rightarrow\pi^{+}\pi^{-}J/\Psi$.
The 90\% CL constraints on branching fractions for invisible decays of $\Upsilon (1S)$ and $J / \Psi$, as
measured by BaBar and BES~\cite{BES:2008}, respectively, are
\begin{eqnarray}
{\cal B} ( \Upsilon (1S) \rightarrow invisible ) < 3.0 \times 10^{-4},
\nonumber\\
{\cal B} ( J / \Psi  \rightarrow invisible ) < 7.2 \times 10^{-4} .
\end{eqnarray}
There is a Standard Model contribution to invisible bound state decays, namely, the decay of
a meson to $\bar \nu \nu$ via a $Z$ boson.
But these partial widths have been calculated and are negligible~\cite{Chang:1998}:
\begin{eqnarray}
{ \cal B} ( \Upsilon (1S) \rightarrow \nu \bar \nu ) &=& 9.85 \times 10^{-6} ,
\nonumber\\
{ \cal B} ( J / \Psi \rightarrow \nu \bar \nu ) &=& 2.70 \times 10^{-8} .
\end{eqnarray}

For each contact operator, one
can calculate the branching fraction for the bound state to decay to dark matter in terms of
the bound state mass, the mediation scale, the dark matter mass, and
the branching fraction to $e^{+} / e^{-}$  (assuming $q=b$ for $\Upsilon (1S)$ decay, or $q=c$ for $J/\psi$ decay):
\begin{eqnarray}
{\cal B}_{F 5}  (\bar X X )&=& { {\cal B} ( e^{+} e^{-} ) M^4 \over 16 \pi^2 \alpha^2 Q^2 \Lambda^4 }
\left( 1 - {4 m_X^2 \over M^2} \right)^{1/2} \left(1 + {2 m_X^2 \over M^2 } \right) ,
\nonumber\\
{\cal B}_{F 6}  (\bar X X )&=& { {\cal B} ( e^{+} e^{-} ) M^4 \over 16 \pi^2 \alpha^2 Q^2 \Lambda^4 } \left( 1 - {4 m_X^2 \over M^2} \right)^{3 / 2} ,
\nonumber\\
{\cal B}_{F 9}  (\bar X X )&=& { {\cal B} ( e^{+} e^{-} ) M^4 \over 8 \pi^2 \alpha^2 Q^2 \Lambda^4 }
\left( 1 - {4 m_X^2 \over M^2} \right)^{1/2} \left(1 + {8 m_X^2 \over M^2 } \right) ,
\nonumber\\
{\cal B}_{F 10}  (\bar X X )&=& { {\cal B} ( e^{+} e^{-}) M^4 \over 8 \pi^2 \alpha^2 Q^2 \Lambda^4 } \left( 1 - {4 m_X^2 \over M^2} \right)^{3 / 2} ,
\nonumber\\
{\cal B}_{S 3}  (\bar X X )&=& { {\cal B} ( e^{+} e^{-} ) M^4 \over 256 \pi^2 \alpha^2 Q^2 \Lambda^4 } \left( 1 - {4 m_X^2 \over M^2} \right)^{3 / 2} ,
\nonumber\\
{\cal B}_{V 3}  (\bar X X )&=& { {\cal B} ( e^{+} e^{-} ) M^4  \over 128 \pi^2 \alpha^2 Q^2 \Lambda^4 } \left( 1 - {4 m_X^2 \over M^2} \right)^{3/2}
\left(1 + { M^4 \over 8 m_X^4 } \left( 1 - {2 m_X^2 \over M^2} \right)^2 \right) ,
\nonumber\\
{\cal B}_{V 5}  (\bar X X )&=& { {\cal B} ( e^{+} e^{-} ) M^2 \over 16 \pi^2 \alpha^2 Q^2 \Lambda^2 } \left( 1 - {4 m_X^2 \over M^2} \right)^{3/2}
{M^2 \over m_X^2} \left( 1 + { M^2 \over 4 m_X^2} \right) ,
\nonumber\\
{\cal B}_{V 6}  (\bar X X )&=& { {\cal B} ( e^{+} e^{-} ) M^2 \over 16 \pi^2 \alpha^2 Q^2 \Lambda^2 }
\left( 1 - {4 m_X^2 \over M^2} \right)^{1/2} {M^2 \over m_X^2}  \left(1 + {2 m_X^2 \over M^2 } \right) ,
\nonumber\\
{\cal B}_{V 7}  (\bar X X )&=& {  {\cal B} ( e^{+} e^{-} ) M^4 \over 64 \pi^2 \alpha^2 Q^2 \Lambda^4 } \left( 1 - {4 m_X^2 \over M^2} \right)^{3/2}
{M^2 \over m_X^2} ,
\nonumber\\
{\cal B}_{V 9}  (\bar X X )&=& {{\cal B} ( e^{+} e^{-} ) M^4 \over 256 \pi^2 \alpha^2 Q^2 \Lambda^4 } \left( 1 - {4 m_X^2 \over M^2} \right)^{5/2}
{M^2 \over m_X^2} .
\end{eqnarray}
Note that for operator F6 we have written the branching fraction assuming that dark matter is a Dirac fermion.  If the dark matter were
instead a Majorana fermion, the branching fraction would be larger by a factor of $2$ (with a factor of $4$ arising from
the squared matrix element, and a factor of $1/2$ from the reduced phase space of the final state particles).
The result for operator S3, and the corresponding constraints, match that found in~\cite{Yeghiyan:2010}.

Note that if dark matter is spin-1,
the decay rates have terms which scale as $m_X^{-2}$ or $m_X^{-4}$.  These terms arise from final states if
either one or both of the dark matter particles is longitudinally polarized~\cite{Kumar:2013iva}.

In general, unitarity constrains the magnitude of the squared matrix element for both elastic and inelastic
scattering.  One may worry that, for small $m_X$, unitarity would require the presence of corrections which invalidate
the use of the contact operator approximation at tree-level; for example, if the dark matter is a
gauge boson which has become massive due to spontaneous breaking of a dark sector symmetry, then one may also need
to include additional diagrams involving the fields responsible for the spontaneous symmetry breaking.
However, the constraints from unitarity are
trivial when initial particles are at rest, because the elastic scattering cross-section is at threshold (for example,
see~\cite{Endo:2014mja}).  As a result, the contact operator approximation is consistent with unitarity for the
meson decay process in the non-relativistic bound state limit.  This approximation becomes more problematic as one
departs from the non-relativistic limit, but that analysis is beyond the scope of this work.

\subsection{Dark Matter Annihilation}

Applying crossing symmetry to the quark annihilation matrix elements yields the matrix element for dark matter annihilation
to quarks.  If the dark matter can annihilate from an $s$-wave initial state, then gamma ray observations can be used to set
limits on $\Lambda$.  But one should note that, although we have restricted ourselves to effective operators which have a
non-zero matrix element with an $s$-wave meson state, these operators need not necessarily have a non-zero matrix element
with an $s$-wave dark matter initial state.  Only five of these operators which we consider can permit unsuppressed dark matter
annihilation.  The corresponding annihilation cross-sections, at tree-level, are given by

\begin{eqnarray}
\langle \sigma_A^{F 5} v \rangle &=& { 3 \over 2 \pi \Lambda^4  } \left( 1 - { m_q^2 \over m_X^2} \right)^{1/2}
\left( 2 m_X^2 + m_q^2 \right) ,
\nonumber\\
\langle \sigma_A^{F 9} v \rangle &=& { 6 \over \pi \Lambda^4  } \left( 1 - { m_q^2 \over m_X^2} \right)^{1/2}
\left( m_X^2 + 2 m_q^2 \right) ,
\nonumber\\
\langle \sigma_A^{F 10} v \rangle &=& { 6\over \pi  \Lambda^4  } \left( 1 - { m_q^2 \over m_X^2} \right)^{ 3 / 2}
m_X^2 ,
\nonumber\\
\langle \sigma_A^{V 5} v \rangle &=& { 2 \over 3 \pi  \Lambda^2  } \left( 1 - { m_q^2 \over m_X^2} \right)^{ 3 / 2}
 ,
\nonumber\\
\langle \sigma_A^{V 6} v \rangle &=& { 2 \over 3 \pi \Lambda^2  } \left( 1 - { m_q^2 \over m_X^2} \right)^{1/2}
\left( 1 +  {2 m_q^2 \over m_X^2} \right) ,
\end{eqnarray}

These cross sections can be bounded by a stacked analysis of the number of photons
arriving from dwarf spheroidal galaxies~\cite{Fermi:2011,Fermi:2013,Geringer:2011,Geringer:2012,Essig:2009}.
The number of photons
expected to result from dark matter annihilation is the product of an astrophysics-dependent
factor and a particle-physics dependent factor.  The astrophysics-dependent factor depends on the density profile
of the dwarf spheroidal galaxies, and can be estimated from the rotation curves of visible matter.
The particle-physics dependent factor can be expressed as
\begin{eqnarray}
\Phi_{PP} &=& { \langle \sigma_A v \rangle \over 8 \pi m_X^2 } \int_{E_{thr}}^{m_X} \sum_f B_f {d N_f \over d E} d E ,
\label{eqn:gammaflux}
\end{eqnarray}
where $B_f$ is the branching ratio for dark matter to annihilate to a channel $f$ and $E_{thr} = 1 \, {\rm GeV}$ is
the photon energy analysis threshold. $ d N_f / d E$ is the photon spectrum for a given dark matter
annihilation channel. Note that our bounds from dark matter annihilation assume universal quark coupling, with the strongest
contribution to the photon spectrum coming from the $u$- and $d$- quark channels.

The $ 95 \% $ CL limit on $\Phi_{PP}$ arising from observations by Fermi-LAT of dwarf spheroidal galaxies is given by~\cite{Geringer:2011}
\begin{eqnarray}
\Phi_{PP} < 5.0_{-4.5}^{+4.3} \times 10^{-30} \, {\rm cm}^3 \, {\rm s}^{-1} \, {\rm GeV}^{-2},
\end{eqnarray}
where the asymmetric uncertainties are $95 \%$ CL systematic errors~\cite{Fermi:2011}, resulting from the uncertainty in
the mass density profiles of the various satellites.  This bound assumes that dark matter is its own anti-particle
(if the dark matter field is complex, as in the cases we consider,
this number increases by a factor of 2).
The associated bounds on the annihilation cross sections are plotted in Figure~\ref{fig:sigann}. For each
annihilation channel, the photon spectrum was produced by Pythia $6.403$~\cite{Pythia} and accounts for collinear photons
generated during the quark final state parton shower evolution.

The general shape of the bound on the annihilation cross-section is easily understood; for larger dark matter mass, the
bound strengthens with decreasing mass because of the resulting increase in the dark matter number density. At lower
masses the bound begins to sharply weaken because it is no longer possible to create
photons above the analysis threshold. We do not plot bounds on dark matter annihilation for $m_X \lower.7ex\hbox{$\;\stackrel{\textstyle<}{\sim}\;$} 4 \, {\rm GeV}$ since the energy of the final state quarks
becomes close to the hadronization scale.
Annihilation constraints for lighter dark matter would also be cut off by either, for the $\bar c c$ and $\bar s s$ channels,
the threshold to create quark/antiquark final states or, for the $\bar u u$ and $\bar d d$ channels, the analysis threshold.

Given any particle physics model for the dark matter annihilation cross section and final state branching
fractions, one can then determine if the model is consistent with data from Fermi-LAT~\cite{Kumar:2011}.
Note, however, that since all of the relevant operators are only non-vanishing when dark matter is
complex, it is consistent with this analysis for dark matter to be asymmetric.  In this case, only the
particle would be abundant at the current epoch, not the anti-particle; the bounds on dark
matter annihilation in dwarf spheroidal galaxies would thus necessarily be unconstraining.

\begin{figure}[hear]
\center\includegraphics[width=0.5 \textwidth]{Annihilation_from_Fermi_cut.pdf}
\caption{Bounds on the annihilation cross section, $\langle \sigma_A v \rangle$, for dark matter of mass $m_X$
annihilating to quarks in dwarf spheroidal galaxies. Note that the results for annihilation to the $u \bar u $ and $d \bar d$
channels are visually identical.}
\label{fig:sigann}
\end{figure}

It should be noted that, for $m_X \lower.7ex\hbox{$\;\stackrel{\textstyle<}{\sim}\;$} {\cal O}(1)~\, {\rm GeV} $,  $s$-wave dark matter annihilation  can be strongly constrained by anisotropies
in the cosmic microwave background (CMB) arising from the energy deposition caused by annihilation during recombination~\cite{Essig:2013goa}.
But since we do not assume any coupling between dark matter and leptons, the only relevant annihilation channels are $\bar X X \rightarrow \bar q q$.
For these final states, if $m_X \lower.7ex\hbox{$\;\stackrel{\textstyle<}{\sim}\;$} 1~\, {\rm GeV}$, it is more difficult to determine CMB constraints on LDM because of the effects of hadronization.
A full treatment of these constraints is thus beyond the scope of this work.  But one can note that, in general,
the constraints on the scale $\Lambda$ arising from CMB studies (as with gamma-ray searches of dwarf spheroidal
galaxies) either tend to become weaker as $m_X$ decreases, for fermionc dark matter, or strengthen less, relative to collider-based constraints, for vector dark matter.

\subsection{Nuclear Scattering}

Several of the effective operators which we consider will also yield velocity independent terms in
the associated dark matter-nucleon scattering matrix element.

Operators with vector quark bilinears will yield spin-independent scattering, with associated cross sections given by
\begin{eqnarray}
\sigma_{SI}^{p,n} &=& { \mu_p^2 \over 32 \pi (2J_X+1) } \sum_{spins}
\left| \sum_q  { B_q^{p,n} \over m_X m_q} {\cal M}_{Xq \rightarrow Xq}  \right|^2   ,
\end{eqnarray}
where $\mu_p$ is the reduced mass of the dark matter-nucleon system and $J_X$ is
the dark matter spin.
The nucleon form factors associated with the vector quark bilinear are $B_u^p = B_d^n = 2$, $B_u^n = B_d^p = 1$, and
$B_{s,c,b,t}^{p,n}=0$~\cite{Ellis:2001}.

Similarly, operators with tensor quark bilinears will yield spin-dependent scattering, with associated cross sections given by
\begin{eqnarray}
\sigma_{SD}^{p,n} &=& { \mu_p^2 \over 32 \pi (2J_X +1)} \sum_{spins}
\left| \sum_q  {\delta_q^{p,n} \over m_X m_q}
{\cal M}_{Xq \rightarrow Xq}  \right|^2 ,
\end{eqnarray}
The nucleon spin form factors $\delta_q^{p,n}$ can be extracted from data and are roughly given by
$\delta_u^p = 0.54_{-0.22}^{+0.09}$ and $\delta_d^p = - 0.23_{-0.16}^{+0.09}$~\cite{Fan:2010,Anselmino:2009}.  These form factors are
in slight disagreement with lattice calculations~\cite{Belanger:2008,Aoki:1997,LHPC:2002}, although using alternative form
factors would only change our scattering
cross sections by a factor of order unity.

Note that for contact operators which we consider, only the coupling to first-generation quarks is relevant for scattering.  In general,
dark matter can exhibit velocity-independent spin-independent scattering from nucleons through interactions with heavy quarks, but only if
the dark matter bilinear is a scalar.  We do not consider such operators, however, because they cannot contribute to the decay of a
$1^{--}$ meson.  Thus, the bounds on dark matter-quark interactions arising from $\Upsilon (1S)$ or $J/\psi$ decay can only be related
to bounds on dark matter-nucleon scattering if one makes a particular choice for the relative strength of dark matter coupling to light and
heavy quarks.

\section{Results}

We assume that dark matter couples to quarks only through a single effective operator, but with equal coupling
to all quark flavors.  An example of a model which would yield this effective operator realization would be the
case where the mediating particle was a massive vector boson for a new $U(1)$ symmetry under which all quarks have
equal charge and under which the dark matter is also charged (for example, this $U(1)$ could be a linear combination
of $U(1)_{baryon}$ and another $U(1)$ symmetry under which the dark matter is charged).  But, of course, other choices
are possible and can be well-motivated by other UV completions.  We focus on the case of equal couplings simply as a benchmark.
We note also that the coefficient of the effective contact operator can exhibit RG-running between the energy scale relevant
for meson decay and the scale relevant for nuclear scattering; however, this is expected to be a relatively small effect~\cite{Haisch:2013uaa,Crivellin:2014qxa}.

\subsection{Mediator Scale}

In Figure~\ref{fig:lambda} we plot bounds on $\Lambda$, as a function of $m_X$,
arising from limits on invisible $\Upsilon (1S)$ decays and from
dark matter annihilation in dwarf spheroidal galaxies.
Note that the bounds arising from $\Upsilon (1S)$ decays are sensitive only to dark matter coupling to $b$-quarks.  But for values
of $m_X$ for which $\Upsilon (1S)$ decays are kinematically allowed, dark matter annihilation to $\bar b b$ is kinematically
forbidden\footnote{For the case where $m_X > m_b$, bounds arising from radiative $\Upsilon$ decays at the LHC
have been discussed in~\cite{Cotta:2013,Schmidt-Hoberg:2013hba}.}; instead, the bounds on dark matter annihilation arise from dark matter couplings to the lighter quarks.
Note that
the mediator scale, $\Lambda$, is always larger than $\sim 10~\, {\rm GeV}$, implying that the contact operator
approximation is valid at all of the colliders which are relevant for meson decay bounds.

\begin{figure}[hear]
\includegraphics*[width=0.49 \textwidth]{Lambda_vs_DM_mass.pdf}
\includegraphics*[width=0.49 \textwidth]{Lambda_vs_DM_mass_Vectors.pdf}
\caption{Bounds on the mediator scale, $\Lambda$, for fermionic dark matter (left panel) and bosonic dark matter (right panel)
of mass $m_X$ arising from constraints on $\Upsilon (1S) \rightarrow nothing$ decays, from constraints on dark matter annihilation
to light quarks in dwarf spheroidal galaxies, and from monojet/photon searches at ATLAS~\cite{ATLAS:2012ky,ATLAS:2012} and
CMS~\cite{CMS:2012, CMS:2012tea}, as indicated.}
\label{fig:lambda}
\end{figure}

Bounds on dark matter annihilation only arise for the operators which permit dark matter annihilation from an $s$-wave initial
state; dark matter annihilation through other operators is suppressed by factors of $v^2$, implying that current gamma ray observations
can provide no useful bounds.
In particular, if dark matter is spin-0, then dark matter annihilation is always $p$-wave if mediated by a contact operator
which can also mediate the decay of a $1^{--}$ meson, and it is thus unconstrained by gamma ray observations.
Bounds on dark matter annihilation in dwarf spheroidal galaxies are not constraining for $m_X \leq 1~\, {\rm GeV}$,
because very light dark matter cannot produce any photons above the analysis threshold.  Bounds on the dark matter annihilation
cross section have systematic uncertainties related to the dark matter density profile; these uncertainties can weaken these bounds
by up to a factor of 2 or strengthen them by up to a factor of 10.

We also plot bounds on the scale $\Lambda$ for each of the operators (for spin-1/2 dark matter) arising from searches
for monojet and monophoton production at the LHC~\cite{CMS:2012, CMS:2012tea,ATLAS:2012ky,ATLAS:2012} (similar bounds can also be
found for mono-$W$,$Z$ production~\cite{ATLAS:2013}).
These searches place bounds on the cross section for the process $pp \rightarrow XX+jet,\gamma$, where
dark matter interacts with Standard Model quarks through a contact operator.  As such, these bounds are very similar in spirit to
bounds on dark matter-quark couplings arising from meson decay, and in particular these bounds do not dramatically worsen for
very light dark matter.  However, it is important to note one significant difference.  These LHC monojet bounds are only valid if
the contact operator approximation is valid even at the energies of LHC processes, which requires the mediator mass to be ${\cal O}(\, {\rm TeV})$.
For lighter mediators, one cannot perform a model-independent operator analysis; although stringent bounds may be
possible~\cite{Busoni:2013lha,Papucci:2014iwa}, they depend on
the details of dark matter-quark interaction.  Of course, the bounds arising from meson decay are also only valid if the contact operator
approximation is valid, but in this case the relevant energy scale is the meson mass, which is ${\cal O}(1-10~\, {\rm GeV})$.  For mediators much
heavier than $\sim 10~\, {\rm GeV}$, the bounds which arise from meson decay will be valid.

Dark matter with $m_X \lower.7ex\hbox{$\;\stackrel{\textstyle>}{\sim}\;$} 10 \, {\rm MeV} $ is heavy enough to be consistent with the cold dark matter paradigm.
Note, however, that our analysis does not require that dark matter is a thermal relic or that dark matter
annihilation be predominantly to Standard Model final states.
Instead, our focus is on constraints on dark matter interactions arising from astrophysical observations of the annihilation products
in the current epoch, not on constraints on the mechanism of dark matter generation.
Even if the the cross section for dark matter annihilation to visible matter is
small, there could be a large branching fraction for annihilation to the dark sector, evading indirect detection constraints but
allowing for a thermal relic density which can satisfy observational constraints.

\subsection{Complementary Dark Matter Scattering}

Of the contact operators which we consider, the ones which permit velocity-independent scattering are
F5 (SI), F9 (SD), S3 (SI), V3 (SI), and V5 (SD).  The corresponding total dark matter-proton scattering
cross sections are given by
\begin{eqnarray}
\sigma_{SI}^{F5}
= { \mu_p^2 \over \pi \Lambda^4 } \left( B_u^p + B_d^p \right)^2 ,
\nonumber\\
\sigma_{SI}^{S3} = \sigma_{SI}^{V3}
= { \mu_p^2 \over 4 \pi \Lambda^4 } \left( B_u^p + B_d^p \right)^2 ,
\end{eqnarray}
\begin{eqnarray}
\sigma_{SD}^{F9} &=& { 12 \mu_p^2 \over \pi \Lambda^4 }
\left( \delta_u^p + \delta_d^p \right)^2 ,
\nonumber\\
\sigma_{SD}^{V5} &=& { 2 \mu_p^2 \over  \pi \Lambda^2 m_X^2 }
\left( \delta_u^p + \delta_d^p \right)^2 .
\end{eqnarray}

In Figure~\ref{fig:scattering}, we plot the bounds on spin-independent (left panel) and spin-dependent (right panel) scattering mediated by
each of the relevant operators.
We also plot 95\% CL bounds arising from Fermi-LAT constraints on dark matter annihilation in dwarf spheroidal galaxies
and 90\% CL bounds arising from monojet and monophoton searches (CMS~\cite{CMS:2012, CMS:2012tea} and ATLAS~\cite{ATLAS:2012ky,ATLAS:2012}).
The DAMA/LIBRA~\cite{Savage:2008er}, CRESST II (95\% CL)~\cite{Angloher:2011uu}, CoGeNT~\cite{Aalseth:2012if}, and CDMS II(Silicon)~\cite{Agnese:2013rvf}
90\% CL signal regions are also shown, as are the 90\% CL exclusion contours from SuperCDMS~\cite{CDMS:2014}, LUX~\cite{Akerib:2013tjd},
SIMPLE~\cite{Felizardo:2011uw}, PICASSO~\cite{Archambault:2012pm},
and COUPP~\cite{Behnke:2012ys}.

\begin{figure}[hear]
\includegraphics*[width=0.49 \textwidth]{SI_Cross_Section.pdf}
\includegraphics*[width=0.49 \textwidth]{SD_Cross_Section.pdf}
\caption{
Bounds on the dark matter-proton spin-independent (left panel) and spin-dependent (right panel) scattering cross section
for dark matter of mass $m_X$ coupling universally to quarks through the indicated effective contact operator.
The labeled exclusion contours indicate 90\% CL bounds arising from limits on invisible decays of $\Upsilon (1S)$,
95\% CL bounds arising from Fermi-LAT constraints on dark matter annihilation in dwarf spheroidal galaxies,
and 90\% CL bounds arising from monojet searches (CMS~\cite{CMS:2012, CMS:2012tea} and ATLAS~\cite{ATLAS:2012ky}).
The DAMA/LIBRA~\cite{Savage:2008er}, CRESST II (95\% CL)~\cite{Angloher:2011uu}, CoGeNT~\cite{Aalseth:2012if} and
CDMS II(Silicon)~\cite{Agnese:2013rvf} 90\% CL signal regions
are also shown, as are the 90\% CL exclusion contours from
SuperCDMS~\cite{CDMS:2014}, LUX~\cite{Akerib:2013tjd}, SIMPLE~\cite{Felizardo:2011uw}, PICASSO~\cite{Archambault:2012pm},
and COUPP~\cite{Behnke:2012ys}.}
\label{fig:scattering}
\end{figure}

Various experiments are able to set bounds on DM-electron scattering~\cite{Essig:2013,Essig:2012}, but some
assumption of universal dark matter coupling to quarks and leptons is required to, in turn, bound dark matter-nucleon scattering.

Similar bounds on dark matter interactions can be obtained from bounds on $J / \Psi$ invisible decays, and are presented in the appendix.
While these bounds are weaker than those obtained from invisible $\Upsilon (1S)$ decays, they are valid for a larger range of mediator masses
($\lower.7ex\hbox{$\;\stackrel{\textstyle>}{\sim}\;$} M_{J/\psi}$) and directly probe the couple of dark matter to $c$-quarks, thus providing non-trivial complementary constraints.
Note that the contact operator
approximation will begin to break down for mediator scales smaller than $\sim 10~\, {\rm GeV}$.

\section{Conclusions}

We have presented bounds on dark matter-quark contact interactions which can be obtained from high luminosity $B$/charm-factories by
constraining decays
of the form $\Upsilon (1S), J/\psi \rightarrow \bar X X$.
These bounds on low mass dark matter probe a mass range significantly below the threshold of direct dark matter detection experiments
and complement bounds on dark matter interactions obtained from gamma ray searches of dwarf spheroidal galaxies and
from monojet/monophoton/mono-$W$,$Z$ searches at hadron colliders.

In particular, the effective interactions which permit decay of a $1^{--}$ meson state can also permit velocity-independent
dark matter-nucleon scattering (either spin-independent or spin-dependent).  For $m_X \sim 1-5~\, {\rm GeV}$, the bounds obtained from meson decay can
thus potentially complement those obtained from direct detection experiments.  For the case of spin-independent scattering, direct detection
experiments already place bounds which well exceed those obtained from meson decay.  However, for spin-dependent scattering, bounds arising
from $\Upsilon (1S)$ decay via the F9 and V5 operators are comparable to those obtained from direct detection experiments.  This is not surprising,
as direct detection experiments typically have much weaker sensitivity to spin-dependent scattering, due to the lack of constructive interference
in coherent scattering.

Moreover, bounds on spin-1 dark matter interactions improve dramatically as $m_X$ decreases, because of the enhancement in the matrix element
which arises when the dark matter particles are longitudinally polarized.  Relating the suppression scale $\Lambda$ to the
mediator mass scale $m_{med.}$ and
coupling $g$ by $\Lambda \sim m_{med.} /g$, this implies that interactions between spin-1 dark matter and
quarks can be constrained even if the coupling is very weak.

We have seen that invisible quarkonium decays probe the same parton-level process ($\bar q q \rightarrow \bar X X$) as
monojet/photon/$W$,$Z$ searches at hadron colliders.  However, quarkonium decays provide complementary information,
allowing robust probes of models with relatively light mediators ($\lower.7ex\hbox{$\;\stackrel{\textstyle>}{\sim}\;$} 10~\, {\rm GeV}$) for which the contact operator
approximation would fail at the LHC.  It is also worth noting that, since the heavy quarkonium bound states are non-relativistic,
searches based on invisible heavy quarkonium decays readily distinguish between DM-SM interactions which vanish in the limit
of non-relativistic quarks and those which do not.  This probe thus nicely complements LHC searches, in which the partons
are highly relativistic.

A similar analysis can be performed of heavy quarkonium decays to a photon and missing energy;
although the set of relevant contact operators would be
different for such an analysis, those branching fractions are much more tightly constrained.
Improved bounds on the $\Upsilon (1S) \rightarrow nothing$ decay
rate from Belle II~\cite{Browder} and a factor
of $\sim 14$ enhancement in sensitivity to $J / \Psi \rightarrow nothing $ decay rate from BESIII~\cite{Harris} will
allow for even tighter constraints
on the interactions of low mass dark matter with Standard Model particles.

{\bf Acknowledgements}

We are grateful to Y.~G.~Aditya, T.~Browder, F.~Harris, D.~Marfatia, A.~Petrov, A.~Rajaraman, X.~Tata, S.~Vahsen, D.~Walker and A.~Wijangco for useful discussions.
This work is supported in part by Department of Energy grant  DE-SC0010504.


\section{Quadratizations of pseudo-Boolean functions}
A {\em pseudo-Boolean function}  is a real-valued function $f(x)=f(x_1,x_2,\ldots,x_n)$ of $n$ binary
variables, a mapping from $\{0,1\}^n$ to ${\ensuremath{\mathbb{R}}}$. It is well-known that every pseudo-Boolean function can
be uniquely represented as a multilinear polynomial in its variables. \emph{Nonlinear binary optimization
problems}, or \emph{pseudo-Boolean optimization} (PBO) \emph{problems}, of the form
\begin{equation*}\label{eq:PBO}
  \min \bigl\{ f(x) : x \in \{0,1\}^n \bigr\},
\end{equation*}
where $f(x)$ is a pseudo-Boolean function, have attracted the attention of numerous researchers, and
they are notoriously difficult, as they naturally encompass a broad variety of models such as maximum
satisfiability, maximum cut, graph coloring, simple plant location, and so on; see, e.g.,~\cite{CH2011}.
In recent years, several authors have revisited an approach initially proposed by Rosenberg~\cite{Rosenberg75}.
This involves reducing PBO to its quadratic case (QPBO) by relying on the following concept.

\begin{definition}\label{def:quad}
  For a pseudo-Boolean function $f(x)$ on $\{0,1\}^n$, we say that $g(x,y)$ is a \emph{quadratization} of $f$
  if $g(x,y)$ is a quadratic polynomial depending on $x$ and on $m$ \emph{auxiliary} binary variables
  $y_1,y_2,\ldots,y_m,$ such that
  \begin{equation*}\label{eq:quad}
    f(x)= \min \bigl\{ g(x,y) : y \in \{0,1\}^m \bigr\} \quad \mbox{for all } x \in \{0,1\}^n.
  \end{equation*}
\end{definition}

Clearly, if $g(x,y)$ is a quadratization  of $f$, then
\[
  \min \bigl\{ f(x) : x \in \{0,1\}^n \bigr\} = \min \bigl\{ g(x,y) : x \in \{0,1\}^n, y \in \{0,1\}^m \bigr\},
\]
so that the minimization of $f$ is reduced through this transformation to the QPBO problem of minimizing
$g(x,y)$. We are also interested (see~\cite{ABCG}) in special types of quadratizations, which we call
\emph{$y$-linear quadratizations}, which contain no products of auxiliary variables. If $g(x,y)$ is a
$y$-linear quadratization, then $g$ can be written as
\[
  g(x,y)=q(x) + \sum_{i=1}^m a_i(x) y_i,
\]
where $q(x)$ is quadratic in $x$ and each $a_i(x)$ is a linear function of $x$. When minimizing $g$ over
$y$, each product $a_i(x)y_i$ takes the value $\min\{0,a_i(x)\}$. Thus, $y$-linear quadratizations can be
viewed as piecewise linear functions of the $x$-variables.

\begin{example}
  As an easy explicit example, consider the {\em negative monomial}
  \[
    -\prod_{i=1}^n x_i=-x_1x_2\cdots x_n.
  \]
  This elementary pseudo-Boolean function has the following \emph{standard quadratization} (Freedman and
  Drineas~\cite{FreedDrineas2005}):
  \[
    s_n(x_1,x_2,\ldots, x_n,y) = y\left(n-1-\sum_{i=1}^n x_i\right).
  \]
  The reason is as follows: unless all the $x_i$ are $1$, then the quantity in parentheses in the expression
  for $s_n$ is non-negative and the minimum value of $s_n$ is therefore $0$, obtained when $y=0$; and, if all
  $x_i$ are $1$, the expression equals $-y$, minimized when $y=1$, giving value $-1$. In both cases, the minimum
  value of $s_n$ is the same as the value of the negative monomial.
\end{example}

\begin{example}
  The {\em positive monomial} is the function
  \[
    \prod_{i=1}^n x_i=x_1x_2\cdots x_n.
  \]
  Ishikawa~\cite{Ishikawa2011} showed that it can be quadratized using $\left\lfloor \frac{n-1}{2}\right\rfloor$
  auxiliary variables, and this is currently the best available bound for positive monomials; see
  also~\cite{FGBZ2011,Ishikawa2011}.
\end{example}

Rosenberg~\cite{Rosenberg75} has proved that every pseudo-Boolean function $f(x)$ has a quadratization,
and that a quadratization can be efficiently computed from the polynomial expression of $f$. This also easily
follows from our foregoing observations that every monomial has a quadratization. (It is also the case ---
see~\cite{ABCG} --- that any pseudo-Boolean function has a $y$-linear quadratization.)  Of course, quadratic
PBO problems remain difficult in general, but this special class of problems has been thoroughly studied for
the last few decades, and much progress has been made in solving large instances of QPBO, either exactly or
heuristically. Quadratization has emerged in recent years as one of the most successful approach to the solution
of very large-scale PBO problems arising in computer vision applications. (See, for instance, Boykov, Veksler
and Zabih \cite{BVZ01}, Kolmogorov and Rother~\cite{KR2007}, Kolmogorov and Zabik~\cite{KZ2004}, Rother,
Kolmogorov, Lempitsky and Szummer~\cite{RKLS2007}, Boros and Gruber~\cite{BG2011}, Fix, Gruber, Boros and
Zabih~\cite{FGBZ2011}, Freedman and Drineas~\cite{FreedDrineas2005}, Ishikawa~\cite{Ishikawa2011}, Ramalingam,
Russell, Ladick\'{y} and Torr~\cite{RRLT2011}, Rother, Kohli, Feng and Jia~\cite{RKFJ}.)

In a related paper, the present authors~\cite{ABCG} initiated a systematic study of quadratizations of
pseudo-Boolean functions. We investigated the minimum number of auxiliary $y$-variables required in a
quadratization of an arbitrary pseudo-Boolean function. In this paper, our focus is on {\em symmetric}
pseudo-Boolean functions. A symmetric pseudo-Boolean function is one in which the value of the function
depends only on the weight of the input. More precisely, a pseudo-Boolean function $f: \{0,1\}^n \to \ensuremath{\mathbb{R}}$
is {\em symmetric} if there is a function $k: \{0,1,\ldots, n\} \to \ensuremath{\mathbb{R}}$ such that $f(x)=k(l)=k_l$ where
$l=|x|=\sum_{j=1}^n x_j$ is the Hamming weight (number of ones) of $x$. In another way, $f$ is symmetric
if it is invariant under any permutation of the coordinates $\{1,2,\ldots,n\}$ of its variables. Note,
for instance, that the positive and negative monomials are symmetric. Here, we investigate the number of
auxiliary variables required in quadratizations of such functions.

\subsection{Outline}
In Section~\ref{sec:rep}, we present a representation theorem and corollaries, which provide useful ways
of expressing symmetric pseudo-Boolean functions. In Section~\ref{sec:reptoquad} we explain how we can use
such representations to construct quadratizations, and we present the implications for upper bounds on the
number of auxiliary variables in Section~\ref{sec:numberaux}. Section~\ref{sec:lower} presents two types of
lower bounds on the number of auxiliary variables for quadratizations of symmetric functions: an existence
result, establishing that there are symmetric functions needing a rather large number of auxiliary variables;
and a concrete lower bound on the number of auxiliary variables in any $y$-linear quadratization of the parity
function.

\section{A representation theorem}\label{sec:rep}
We introduce a useful piece of notation: for any real number $a$, $\n{a}$ denotes $\min(a,0)$, the smaller
of $a$ and $0$. In this section, we present a result that will be key in our approach to obtaining quadratizations.
This is a `representation theorem' that expresses a symmetric pseudo-Boolean function on variables $x_1, x_2, \ldots, x_n$
as a linear combination of terms of the form $\n{a-\sum_{r=1}^n x_r}$, for a suitable range of values $a$.

Our main result, in its most general form, is as follows.

\begin{theorem} \label{thm:repepsilongeneral}
  Let  $0<\epsilon_i \leq 1$, for $i=0,1,\ldots,n$. Then every symmetric pseudo-Boolean function $f: \{0,1\}^n\to \ensuremath{\mathbb{R}}$
  can be represented uniquely in the form
  \begin{equation*}\label{eq:symrepresent}
    f(x)=\sum_{i=0}^{n}\alpha_i \n{i-\epsilon_i-\sum_{r=1}^n x_r}.
  \end{equation*}
\end{theorem}

When all the $\epsilon_i$ are equal, we can be more explicit about the coefficients in this representation.
Recall that $k_l=k(l)$ is the value of $f(x)$ in any point $x$ with Hamming weight equal to $l$. We set
$k_{-1}=0$ by convention.

\begin{theorem}\label{thm:repepsilon}
  Let $0<\epsilon \le 1$. Then every symmetric pseudo-Boolean function
  $f: \{0,1\}^n\to \ensuremath{\mathbb{R}}$ can be represented uniquely in the form
  \begin{equation*}\label{eq:symrepresent2}
    f(x)=\sum_{i=0}^{n}\alpha_i \n{i-\epsilon-\sum_{r=1}^n x_r}
  \end{equation*}
  where, for $j=0,1,\ldots, n$, the value of $\alpha_j$ is
  \begin{equation}\label{eq:alpha}
    \alpha_j=- \sum_{i=0}^{j-2} \frac{(\epsilon-1)^{j-i-2}}{\epsilon^{j-i+1}}  k_i
             + \left(\frac{1}{\epsilon} +\frac{1}{\epsilon^2}\right) k_{j-1} -\frac{1}{\epsilon} k_j.
  \end{equation}
\end{theorem}
(The first sum in~\eqref{eq:alpha} is, by usual convention, taken to be $0$ if $j<2$.)

\begin{proof}[Proof of Theorem~\ref{thm:repepsilongeneral}]
  When $\sum_{r=1}^n x_r =j$, we should have $f(x)=k_j$. So, to find a coefficient vector
  $\ve{\alpha}=(\alpha_0, \alpha_1, \ldots, \alpha_{n})^T \in \ensuremath{\mathbb{R}}^{n+1}$ which establishes the required
  representation, we need to find a solution to the following system of $n+1$ linear equations:
  \begin{equation}\label{eq:alpha2}
    k_j = \sum_{i=0}^{n}\alpha_i \n{i-\epsilon_i-j}
        = \sum_{i=0}^{j}\alpha_i \, (i-\epsilon_i-j) \quad \mbox{for } j=0,1,\ldots,n.
 \end{equation}
 The matrix underlying this system is the lower-triangular matrix
 \begin{align*}
   A &= \left(\n{q -\epsilon_q -p}\right)_{p,q=1,2,\ldots,n+1}\\
     &=
   \begin{pmatrix}
     -\epsilon_0 & 0 & 0 & 0 & \cdots & 0 \cr
     -1-\epsilon_0 & -\epsilon_1 & 0 &  0 & \cdots & 0 \cr
     -2-\epsilon_0 & -1-\epsilon_1 & -\epsilon_2 & 0 & \cdots & 0 \cr
     \vdots & \vdots & \vdots & \vdots & \ddots &  \vdots \cr
     -n-\epsilon_0 & -n+1-\epsilon_1 & -n+2-\epsilon_2 & -n+3-\epsilon_3 & \cdots & -\epsilon_n
   \end{pmatrix}\\
  \end{align*}
  Because $A$ is lower-triangular with nonzero diagonal entries $-\epsilon_{q}$ ($q=0,1,\ldots,n$), this
  system does indeed have a unique solution and therefore the representation exists and is unique.
\end{proof}

\begin{proof}[Proof of Theorem~\ref{thm:repepsilon}]
  We check that a solution (and hence {\em the} solution) of the system~\eqref{eq:alpha2} with
  $\epsilon_i=\epsilon$ for all $i=0,1,\ldots,n$, is given by~\eqref{eq:alpha} in the statement of the theorem.

  We proceed by induction on $j$. The case $j=0$ is easily verified, since the first equation
  in~\eqref{eq:alpha2} immediately yields $\alpha_0=-\frac{1}{\epsilon} k_0$, in agreement
  with~\eqref{eq:alpha}. Assume now that (\ref{eq:alpha}) is satisfied by the values of $\alpha_i$ up
  to $i=j-1$. Then, from (\ref{eq:alpha2}) and from the induction hypothesis,
  \begin{align}
    -\epsilon \alpha_j &= k_j + \sum_{i=0}^{j-1}\alpha_i \, (\epsilon+j-i) \nonumber\\
    &= k_j + \sum_{i=0}^{j-1}\alpha_i \, (\epsilon+j-1-i) + \sum_{l=0}^{j-1}\alpha_l \nonumber\\
    &= k_j - k_{j-1} + \sum_{l=0}^{j-1}\alpha_l. \label{eq:tech1}
  \end{align}
  Substituting (\ref{eq:alpha}) in the last term of~(\ref{eq:tech1}) yields
  \begin{align}
    \sum_{l=0}^{j-1}\alpha_l &= \sum_{l=0}^{j-1} \left[
      - \sum_{i=0}^{l-2} \frac{(\epsilon-1)^{l-i-2}}{\epsilon^{l-i+1}}  k_i
      + \left(\frac{1}{\epsilon} +\frac{1}{\epsilon^2}\right) k_{l-1} -\frac{1}{\epsilon} k_l \right] \nonumber\\
    &= - \sum_{i=0}^{j-3} k_i  \sum_{l=i+2}^{j-1} \frac{(\epsilon-1)^{l-i-2}}{\epsilon^{l-i+1}}
      + \left(\frac{1}{\epsilon} +\frac{1}{\epsilon^2}\right) \sum_{l=0}^{j-1} k_{l-1} -\frac{1}{\epsilon}
      \sum_{l=0}^{j-1} k_l \nonumber\\
    &= - \sum_{i=0}^{j-3} k_i \sum_{t=0}^{j-i-3} \frac{(\epsilon-1)^{t}}{\epsilon^{t+3}}
      - \frac{1}{\epsilon} k_{j-1} +  \frac{1}{\epsilon^2} \sum_{l=0}^{j-1} k_{l-1}\label{eq:penultimate}\\
    &= \sum_{i=0}^{j-3}    \frac{(\epsilon-1)^{j-i-2}}{\epsilon^{j-i}} \, k_i
      - \frac{1}{\epsilon} k_{j-1} +  \frac{1}{\epsilon^2}  k_{j-2} \label{eq:tech2}
  \end{align}
  where the last equality is obtained by summing the geometric series which appears in the first sum
  of equation~(\ref{eq:penultimate}).

  Combining (\ref{eq:tech1}) and (\ref{eq:tech2}), we find
  \[
    \alpha_j = - \sum_{i=0}^{j-3}  \frac{(\epsilon-1)^{j-i-2}}{\epsilon^{j-i+1}} \, k_i
    -\frac{1}{\epsilon^3} k_{j-2} + \left(\frac{1}{\epsilon} +\frac{1}{\epsilon^2}\right) k_{j-1} -\frac{1}{\epsilon} k_j,
  \]
  which is equivalent to (\ref{eq:alpha}).
\end{proof}

There are two special cases of Theorem~\ref{thm:repepsilon} that we will use in particular.
When $\epsilon=1/2$, Theorem~\ref{thm:repepsilon} yields:

\begin{corollary}\label{thm:rep2}
  Every symmetric pseudo-Boolean function $f: \{0,1\}^n\to \ensuremath{\mathbb{R}}$ can be represented uniquely in the form
  \[
    f(x)=\sum_{i=0}^{n}\alpha_i \n{i-\frac12-\sum_{r=1}^n x_r}
  \]
  where
  \[
    \alpha_i=-8 \sum_{j=0}^i (-1)^{i-j}k_j - 2k_{i-1}+6 k_i
  \]
  for $i=0,1,\ldots,n$, and  $k_{-1}=0$.
  \qed
\end{corollary}

Taking $\epsilon=1$ in Theorem~\ref{thm:repepsilon}, we obtain the following.

\begin{corollary}\label{thm:fix}
  Every symmetric pseudo-Boolean function $f: \{0,1\}^n\to \ensuremath{\mathbb{R}}$ can be represented in the form
  \[
    f(x)=k_0+(k_1-k_0)\sum_{r=1}^n x_r +\sum_{i=1}^{n-1}\left(-k_{i-1}+2k_i-k_{i+1}\right) \n{i-\sum_{r=1}^n x_r}.
  \]
  \qed
\end{corollary}

In fact, Corollary~\ref{thm:fix} follows from work of Fix~\cite{Fix}, and a simpler direct proof can be given.
As Fix observed, if  $\sum_{r=1}^n x_r=l$, then
\[
  f(x)=k(l)=\sum_{i=0}^n k(i)\delta_i(l)
\]
where $\delta_i(l)=1$ if $i=l$ and $\delta_i(l)=0$ otherwise. Then, it can be seen that
\[
  \delta_i(l)=-\n{i-1-l}+2\n{i-l}-\n{i+1-l}.
\]
From this, it follows that
\[
  f(x)=\sum_{i=0}^n k(i)\left(-\n{i-1-l}+2\n{i-l}-\n{i+1-l}\right).
\]
On simplification, this gives
\[
  f(x)=k_0+l(k_1-k_0) + \sum_{i=1}^{n-1} (-k_{i-1}+2k_i-k_{i+1})\n{i-l},
\]
as required.

\section{From representations to quadratizations}\label{sec:reptoquad}
In this section we explain how a representation of the type presented in the previous section can be
used to construct  quadratizations of pseudo-Boolean functions. One useful observation is that when a
coefficient $\alpha_i$ is non-negative, the corresponding term $\alpha_i \n{i-\epsilon_i-\sum_{r=1}^n x_r}$
in the representation of Theorem~\ref{thm:repepsilongeneral} of $f$ can be quadratized as
$\min_{y_i} \alpha_i y_i (i-\epsilon_i-\sum_{r=1}^n x_r)$. But this translation simply does not work if
$\alpha_i$ is negative. The strategy described in this section is to take an expression as given in
Theorem~\ref{thm:repepsilongeneral} (or one of its special cases) and add a quantity that is identically-$0$
and which will result in a final expression that has no terms with negative coefficients. The following
Lemma describes three possible such quantities. The first is going to be useful when working with
representations of the form given in Corollary~\ref{thm:fix}, and the second and third will be useful when
working with the representations from Corollary~\ref{thm:rep2}.

\begin{lemma}\label{lem:zero1}
  Let
  \begin{align*}
    E(l) &=l(l-1)+2\sum_{i=1}^{n-1} \n{i-l}, \\
    E'(l)&=\frac{l(l-1)}{2}+ 2 \sum_{\stackrel{i=2:}{i \; \text{even}}}^{n}\n{i-\frac12-l}, \label{eq:zero1}
  \end{align*}
  and
  \begin{equation*}\label{eq:zero2}
    E''(l)=\frac{l(l+1)}{2}+ 2 \sum_{\stackrel{i=1:}{i \; \text{odd}}}^{n}\n{i-\frac12-l}.
  \end{equation*}
  Then, for all $l=0,1,\ldots,n$, $E(l)= E'(l)=E''(l)=0$.
\end{lemma}

\begin{proof}
  First we show that $E(l)$ is identically-$0$. We have
  \begin{align*}
    E(l) &= l(l-1)+2\sum_{i=1}^{n-1} \n{i-l} \\
         &= l(l-1)+2\sum_{i=1}^{l-1} (i-l) \\
         &= l(l-1)-2\sum_{j=1}^{l-1} j \\
         &= l(l-1)-l(l-1)=0.
  \end{align*}
  We next show that $E'(l) = 0$ for all values of $l$. Fix $l$ and
  note that $i-\frac12-l \leq 0$ if and only if $i\leq l$. Hence,
  \begin{equation}\label{eqn:E}
      \sum_{\stackrel{i=2:}{i \; \text{even}}}^{n} \n{i-\frac12-l}
    = \sum_{\stackrel{i=2:}{i \; \text{even}}}^{l} \left(i-\frac12-l\right)
    = \sum_{\stackrel{i=2:}{i \; \text{even}}}^{l}  i - \left(\frac12+l\right) \left\lfloor \frac{l}{2}\right\rfloor.
  \end{equation}
  By considering separately the cases where $l$ is respectively even or odd, one can conclude
  that $E'(l)=0$ for all $l =0,\ldots,n$. For, if $l=2r$, then the expression on the right in
  equation~(\ref{eqn:E}) is $r/2-r^2=-l(l-1)/4$ and, if $l=2r+1$, it is $-r/2-r^2=-l(l-1)/4$.
  The identity $E''(l) =0$ (for all $l$) can be proved similarly, or can be deduced from the
  previous one by observing  that, for all $l=0,1,\ldots,n$,
  \[
    \sum_{i=1}^{n} \n{i-\frac12-l} = \sum_{i=1}^{l} \left(i-\frac12-l\right) = -\frac12 l^2.
  \]
  We then can note that $$E'(l)+E''(l)=l^2+2\sum_{i=1}^n \n{i-\frac{1}{2}-l}=l^2-l^2=0,$$ so that $E''=-E'=0$.
\end{proof}

We gave a direct, self-contained, proof of Lemma~\ref{lem:zero1}, but in fact these three identities follow
from Corollary~\ref{thm:rep2} and Corollary~\ref{thm:fix}. For, if we apply Corollary~\ref{thm:fix} to the
function $f(x)=\sum_{r=1}^n x_r \left(\sum_{r=1}^nx_r -1\right)$, we see that
\[
  f(x)=-2\sum_{i=1}^{n-1}\n{i-\sum_{r=1}^n x_i},
\]
which implies the first identity of Lemma~\ref{lem:zero1}. Applying Corollary~\ref{thm:rep2} to $f(x)$ shows
(after some calculation) that
\[
  f(x)=-4\sum_{\stackrel{i=2:}{i \; \text{even}}}^{n}\n{i-\frac12-l},
\]
giving the second identity (that $E'$ is identically-$0$). Applying Corollary~\ref{thm:rep2} to the function
$g(x)=\sum_{r=1}^n x_r \left(\sum_{r=1}^nx_r+1\right)$ yields the third identity.

\section{Upper bounds on number of auxiliary variables}\label{sec:numberaux}
\subsection{Any symmetric function}
We first have the following very general result, which provides an explicit construction of a quadratization
of any pseudo-Boolean function, using no more than $n-2$ auxiliary variables.

\begin{theorem}\label{thm:symm_quad}
  Every symmetric function of $n$ variables can be quadratized using $n-2$ auxiliary variables.
\end{theorem}

\begin{proof}
  Using Corollary~\ref{thm:rep2}, we can write any symmetric function $f$ as
  \[
    f(x)= -\alpha_0 \left(\frac12+\sum_{j=1}^n x_j\right) + \sum_{i=1}^{n}\alpha_i \n{i-\frac12-\sum_{j=1}^n x_j}.
  \]
  Let $\alpha_r = \min\{\alpha_i: i \mbox{ even}, i \geq2\}$ and $\alpha_s= \min\{\alpha_i: i \mbox{ odd}\}$.
  Now add to $f$ the expression
  \[
    -\frac{{\alpha_r}}{2}E'\left(\sum_{j=1}^n x_j\right) -\frac{{\alpha_s}}{2}E''\left(\sum_{j=1}^n x_j\right),
  \]
  which is identically-$0$. This results in an expression for $f$ of the form
  \[
    f(x)= a_0 + a_1\sum_{j=1}^n x_j +a_2 \sum_{1\leq i<j \leq n} x_ix_j
    + \sum_{i=1}^{n}\beta_i \n{i-\frac12-\sum_{j=1}^n x_j},
  \]
  where, for each $i$, if $i$ is even, $\beta_i=\alpha_i-{\alpha_r} \ge 0$, and if $i$ is odd,
  $\beta_i=\alpha_i-{\alpha_s} \ge 0$. So all the coefficients $\beta_i$ are non-negative. Furthermore,
  $\beta_r=\beta_s=0$, so we have an expression for $f$ involving no more than $n-2$   positive coefficients
  $\beta_i$. Then,
  \[
    g(x,y)= a_0 + a_1\sum_{j=1}^n x_j +a_2 \sum_{1\leq i<j \leq n} x_ix_j
    + \sum_{\stackrel{i=1:}{i \neq r,s}}^{n}\beta_i y_i\left(i-\frac12-\sum_{j=1}^n x_j\right)
  \]
  is a quadratization of $f$ involving at most $n-2$ auxiliary variables.
\end{proof}
(A construction in~\cite{Fix} shows an upper bound of $n-1$. This is obtained by adding a multiple of
$E(\sum_{r=1}^n x_r)$ to each term in the expression from Corollary~\ref{thm:fix}, rather than to the
expression as a whole, resulting in more complex quadratizations.)

Notice that the quadratization in the proof of Theorem~\ref{thm:symm_quad} is $y$-linear, so we have in fact shown:

\begin{theorem}\label{thm:upperlinear}
  Every symmetric function of $n$ variables has a $y$-linear quadratization involving at most
  $n-2$ auxiliary variables.
\end{theorem}

Furthermore, these quadratizations are also symmetric in the $x$-variables. Not every quadratization of a
symmetric function must itself be symmetric in the original variables. For example, consider the negative
monomial
\[
  -\prod_{i=1}^n x_i=-x_1x_2\cdots x_n.
\]
As we have seen, this has the quadratization $y\left(n-1-\sum_{j=1}^n x_j\right)$, which is symmetric.
However, it also has the quadratization
\[
  (n-2)x_{n}y - \sum_{i=1}^{n-1} x_i (y-{\overline x}_n),
\]
where ${\overline x}_n=1-x_n$, which is not symmetric in the $x$-variables.

\subsection{Monomials}
The quadratization of monomials (positive and negative) has been fairly well-studied. The \emph{standard quadratization}
of the negative monomial
\[
  f(x)=-\prod_{i=1}^n x_i=-x_1x_2\cdots x_n,
\]
is
\[
  s_n(x_1,x_2,\ldots, x_n,y)=y\left(n-1-\sum_{j=1}^n x_j\right).
\]
(A related paper by the present authors~\cite{ABCG} gives a complete characterization of all the quadratizations
of negative monomials involving one auxiliary variable and this is, in a sense, one of  the simplest.) If we apply
Corollary~\ref{thm:rep2} to the negative monomial, noting that $k_i=0$ for $i<n$ and $k_n=-1$, we obtain the
representation
\[
  f(x)=2\n{n-\frac{1}{2}-\sum_{r=1}^nx_r},
\]
which immediately leads to the quadratization
\[
  h=2y\left(n-\frac{1}{2}-\sum_{r=1}^nx_r\right),
\]
only slightly different from the standard one. We could, instead, apply Corollary~\ref{thm:fix}, which would show
that $f(x)=\n{n-1-\sum_{r=1}^nx_r},$ from which we immediately obtain the standard quadratization.

As we noted earlier, the best known result (smallest number of auxiliary variables) for positive monomials is that
they can be quadratized using $\left\lfloor \frac{n-1}{2} \right\rfloor$ auxiliary variables. This was shown by
Ishikawa~\cite{Ishikawa2011}. We can see that this many auxiliary variables suffice by using our representation
theorem, Corollary~\ref{thm:rep2}, together with the argument given in the proof of Theorem~\ref{thm:symm_quad}.

\begin{theorem} \label{thm:positiveterm}
  The positive monomial $P=\prod_{i=1}^n x_i$ can be quadratized using $\left\lfloor \frac{n-1}{2} \right\rfloor$
  auxiliary variables.
\end{theorem}

\begin{proof}
  Consider first the case where $n$ is even. By Corollary~\ref{thm:rep2}, noting that $k_i=0$ for $i<n$ and $k_n=1$,
  we have $P = -2 \n{n-\frac{1}{2}-l}$ where $l=\sum_{r=1}^n x_r$. By Lemma~\ref{lem:zero1},
  \begin{align*}
    P &= -2 \n{n-\frac{1}{2}-l} + E'(l) \\
      &= \frac{l(l-1)}{2}+ \sum_{\stackrel{i=2:}{i \; \text{even}}}^{n-2}2\n{i-\frac{1}{2}-l} \\
      &= \sum_{1\leq i<j \leq n} x_ix_j + \min_y \sum_{\stackrel{i=2:}{i \; \text{even}}}^{n-2}2 y_i\left(i-\frac{1}{2}-l\right).
  \end{align*}
  This provides the required quadratization using $\frac{n}{2}-1=\left\lfloor \frac{n-1}{2} \right\rfloor$ new
  variables.

  When $n$ is odd, one similarly derives the following from Lemma~\ref{lem:zero1}:
  \begin{align*}
    P &= -2 \n{n-\frac{1}{2}-l} + E''(l) \\
      &= \sum_{i=1}^n x_i + \sum_{1\leq i<j \leq n} x_ix_j
       + \min_y \sum_{\stackrel{i=1:}{i \; \text{odd}}}^{n-2}2 y_i\left(i-\frac{1}{2}-l\right).
  \end{align*}
\end{proof}

This quadratization of $P$ requires the same number of auxiliary variables as Ishikawa's construction. Both
quadratizations are, in fact, identical when $n$ is even, but appear to be different when $n$ is odd.

Note that an alternative approach to the case of odd $n$ would be as follows. Write
\[
  P = \prod_{i=1}^{n-1} x_i - \prod_{i=1}^{n-1} x_i {\overline x}_n,
\]
where ${\overline x}_n=1-x_n$. The first term can now be quadratized using $\frac{n-1}{2}-1$ new variables (since it
contains an even number of variables), and the second term, viewed as a negative monomial in $x_1,\ldots,x_{n-1},{\overline x}_n$,
has a standard quadratization requiring one further auxiliary variable. Thus, this leads again
to a quadratization of $P$ with $\frac{n-1}{2}=\left\lfloor \frac{n-1}{2}\right\rfloor$ new variables.
This quadratization is also different from Ishikawa's.

\subsection{$t$-out of $n$ and exact-$t$ functions}
Consider now the \emph{$t$-out-of-$n$} function defined by:
\[
  f_{t,n}(x)=1\qquad \text{if and only if}\qquad \sum_{i=1}^n x_i \geq t.
\]

The basic Boolean functions $\textrm{And}_n(x):=\prod_{i=1}^n x_i$ (a positive monomial) and
$\textrm{Or}_n(x):=1-\prod_{i=1}^n(1-x_i)$ are examples of $t$-out of $n$ functions with $t=n$
and $t=1$, respectively. Another popular example is the \emph{majority} function given by:
\[
  \textrm{Maj}_n(x) := \begin{cases}
    1 & \ \text{if}\ \sum_{i=1}^n x_i\geq \left\lceil n/2 \right\rceil,\\
    0 & \ \text{otherwise},
  \end{cases}
\]
which breaks ties in favor of ones when $n$ is even. In this case, $t=\left\lceil n/2 \right\rceil$.

\begin{corollary}\label{thm:k-out-of-n}
  The $t$-out-of-$n$ function $f_{t,n}$ can be quadratized using $\left\lceil n/2 \right\rceil$
  auxiliary variables.
\end{corollary}

\begin{proof}
  From Corollary~\ref{thm:rep2}, $f_{t,n}$ can be represented in the form
  \begin{equation}\label{eq:ftn}
    f_{t,n}(x)=\sum_{i=0}^{n}\alpha_i \n{i-\frac12-\sum_{j=1}^n x_j}
  \end{equation}
  where $\alpha_i=0$ when $i<t$, $\alpha_t=-2$, and $\alpha_i=4(-1)^{i-t-1}$ when $i>t$.

  Since the terms of $f_{t,n}$ alternate in sign when $i\geq t$, we can again use Lemma~\ref{lem:zero1}
  to make all coefficients non-negative by adding either $2E'(l)$ or $2E''(l)$ to~(\ref{eq:ftn}), depending
  on the parity of $t$. The resulting expression has $\left\lceil n/2 \right\rceil$ positive
  coefficients, and its remaining coefficients are zero. Thus, it can be quadratized with
  $\left\lceil n/2 \right\rceil$ auxiliary variables.
\end{proof}

A related function is the \emph{exact-$t$} (out of $n$) function, defined as $f^=_{t,n}(x)=1$ if and only if
the Hamming weight of $x$ equals $t$. Using Corollary~\ref{thm:rep2} again, we have that $f^=_{t,n}$
can be represented in the form given in \eqref{eq:ftn} with $\alpha_i=0$ when $i<t$, $\alpha_t=-2$,
$\alpha_{t+1}=6$, and $\alpha_i=0$ when $i>t+1$. Depending on the parity of $t$, we add $E'(l)$ or
$E''(l)$ to \eqref{eq:ftn} to obtain an expression with $\left\lfloor n/2\right\rfloor$ positive
coefficients, which can then be quadratized with $\left\lfloor n/2\right\rfloor$ auxiliary variables.
We have just proved the following:

\begin{corollary}\label{thm:exact-k}
  The exact-$t$ function $f^=_{t,n}$ can be quadratized using $\left\lfloor n/2 \right\rfloor$
  auxiliary variables.
  \qed
\end{corollary}

The positive monomial and the $\textrm{And}_n$ Boolean function are also special cases of exact-$t$
functions, both with $t=n$. It is apparent from the argument leading to Corollary~\ref{thm:exact-k}
that the reason the positive monomial (and hence, the $\textrm{And}_n$ function) requires
$\left\lfloor \frac{n-1}{2} \right\rfloor$ auxiliary variables instead of $\left\lfloor n/2 \right\rfloor$
is precisely because $t=n$.

\subsection{Parity and its complement}
The \emph{parity} function is the (pseudo-)Boolean function $\Pi(x)$ such that $\Pi(x)=1$ if the Hamming weight
of $x$ is odd, and $\Pi(x)=0$ otherwise. To derive a quadratization of this function, we will use
Corollary~\ref{thm:fix} rather than Corollary~\ref{thm:rep2}, and will make use of a variant of the argument
given to establish Theorem~\ref{thm:symm_quad}. By Corollary~\ref{thm:fix}, we can see that $\Pi$ has the
representation
\begin{equation}\label{parity}
  \Pi(x) = \sum_{j=1}^n x_j+2\sum_{i=1}^{n-1}(-1)^{i-1}\n{i-\sum_{j=1}^n x_j}.
\end{equation}

Let $E(l)$ be as in Lemma~\ref{lem:zero1}. By adding $E(\sum_{j=1}^n x_j)$ to this representation of $\Pi$,
we obtain a representation with non-negative coefficients, which leads to a quadratization with
$m=\left\lfloor n/2\right\rfloor$ auxiliary variables: $\Pi(x)=\min_{y \in \{0,1\}^m} g(x,y)$ where
\[
  g(x,y)=2 \sum_{i<j}x_ix_j+ \sum_{j=1}^n x_j + 4\,\sum_{\stackrel{i=1:}{i \;\text{odd}}}^{n-1} y_i \left(i-\sum_{j=1}^n x_j\right).
\]
(The terms with coefficient $-2$ in the expansion (\ref{parity}) disappear on the addition of $E$.)

The complement, $\overline{\Pi}$ of $\Pi$ can be represented as
\[
  \overline{\Pi}(x)= 1-\sum_{j=1}^n x_j +2\sum_{i=1}^{n-1} (-1)^{i} \n{i -\sum_{j=1}^n x_j},
\]
so, by adding $E(\sum_{j=1}^n x_j)$, to eliminate negative coefficients, we arrive at the following
quadratization involving $m=\left\lfloor \frac{n-1}{2}\right\rfloor$ auxiliary variables:
\[
  g'(x,y)=1+2 \sum_{i<j}x_ix_j -\sum_{j=1}^n x_j +4\,\sum_{\stackrel{i=2:}{i \;\text{even}}}^{n-1} y_i \left(i-\sum_{j=1}^n x_j\right).
\]
So we conclude the following:

\begin{theorem}
  The parity function of $n$ variables has a $y$-linear quadratization involving $\left\lfloor n/2\right\rfloor$
  auxiliary variables, and its complement has a $y$-linear quadratization involving $\left\lfloor \frac{n-1}{2}\right\rfloor$
  auxiliary variables.
\end{theorem}

\section{Lower bounds on the number of auxiliary variables}\label{sec:lower}
\subsection{Generic lower bounds}
The following result is inspired by (but is different and does not follow from) a transformation given in Siu,
Roychowdhury and Kailath~\cite{Siu95}, in the framework of the representation of Boolean functions by threshold
circuits. This result relates quadratizations of arbitrary (possibly non-symmetric) pseudo-Boolean functions to
the quadratization of symmetric functions on a larger, related, number of variables. We will then use a lower
bound result from~\cite{ABCG} in order to obtain a lower bound result for symmetric functions.

\begin{lemma}~\label{lem:arbtosym}
  Suppose that $n,m $ are positive integers and suppose that every symmetric pseudo-Boolean function
  $F(z)$ of $N=2^n-1$ variables (that is, every symmetric function $F:\{0,1\}^{2^n-1}\to \ensuremath{\mathbb{R}}$) has an
  $m$-quadratization. Then every  (arbitrary) pseudo-Boolean function $f(x)$ on $\{0,1\}^n$ also has an
  $m$-quadratization.
\end{lemma}

\begin{proof}
  Let $f(x)$ be an arbitrary pseudo-Boolean function of $n$ variables. We are going to construct a sequence of
  four functions $k$, $F$, $G$, $g$, such that $g$ is a quadratization of $f$. For this purpose, let $N=2^{n}-1$.
  \begin{enumerate}
  \item
    Let $k:\{0,1,\ldots, N\} \to \ensuremath{\mathbb{R}}$ be defined as follows: $k(w) := f(x)$ where $x$ is the binary representation
    of $w$, that is, $w=\sum_{i=1}^n 2^{i-1} x_i$.

  \item
    Let $F$ be the symmetric pseudo-Boolean function of $N$ variables defined by: for all $z \in \{0,1\}^N$,
    $F(z) := k(|z|)$, where $|z|$ is the Hamming weight of $z$. (This defines $F$ completely, given that it is
    symmetric.)

  \item
    Let $G(z,y)$ be an arbitrary quadratization of $F(z)$ using $m$ auxiliary variables. (The hypothesis of
    the theorem is that such quadratizations exist.)

  \item
    Finally, let $g(x,y)$ be the pseudo-Boolean function on $\{0,1\}^{n+m}$ that is obtained by identifying each of the
    variables $z_{2^{j-1}}, z_{2^{j-1}+1},\ldots,z_{2^{j}-1}$ with $x_{j}$ in $G(z,y)$, for $j=1,2,\ldots,n$; that is,
    \[
      g(x_1,x_2,x_3,\ldots,x_n,y) := G(x_1,x_2,x_2,x_3,x_3,x_3,x_3,\ldots,x_n,\ldots,x_n,y).
    \]
    (The unification makes sense since $2^{j-1}x_j = z_{2^{j-1}} + z_{2^{j-1}+1} +\cdots + z_{2^{j}-1}$, for all
    $j=1,2,\ldots,n$.)
  \end{enumerate}
  We claim that $g(x,y)$ is a quadratization of $f$. Indeed, $g$ is clearly quadratic, because $G$ is. Moreover,
  for every point $x \in \{0,1\}^n$,
  \begin{align}
    &\hspace*{-1truecm}\min \bigl\{ g(x,y) : y \in \{0,1\}^m \bigr\} \nonumber\\
    &=  \min \bigl\{ G(x_1,x_2,x_2,x_3,x_3,x_3,x_3,\ldots,x_n,y) : y \in \{0,1\}^m \bigr\} \label{eq:g1}\\
    &=  F(x_1,x_2,x_2,x_3,x_3,x_3,x_3,\ldots,x_n) \label{eq:g2}\\
    &=  k\left(\sum_{i=1}^n 2^{i-1} x_i\right) \label{eq:g3}\\
    &=  f(x) \label{eq:g4}
  \end{align}
  (where equality \eqref{eq:g1} is by definition of $g$, \eqref{eq:g2} is by definition of $G$, \eqref{eq:g3} is
  by definition of $F$, and \eqref{eq:g4} is by definition of $k$).
\end{proof}

We will now make use of the following result from~\cite{ABCG}.

\begin{theorem}\label{thm:lowergeneral}
  There are pseudo-Boolean functions of $n$ variables for which any quadratization must involve at least
  $\Omega(2^{n/2})$ auxiliary variables.
\end{theorem}

To be more concrete, the analysis in~\cite{ABCG} implies that  for any $n \ge 8$, there is a pseudo-Boolean
function on $n$ variables for which any quadratization will require at least $2^{n/2}/8$ auxiliary variables.

This leads to the following lower bound result for symmetric functions.

\begin{theorem}\label{thm:lb_sym1}
  There exist symmetric functions of $n$ variables for which any quadratization must involve
  at least $\Omega(\sqrt{n})$ auxiliary variables.
\end{theorem}

\begin{proof}
  Lemma~\ref{lem:arbtosym} shows that, if every symmetric function $F(z)$ on $\{0,1\}^N$, with $N=2^n-1$, has
  an $m$-quadratization, then every (arbitrary) function $f(x)$ on $\{0,1\}^n$ also has an $m$-quadratization.
  On the other hand, from Theorem~\ref{thm:lowergeneral}, we know that some pseudo-Boolean functions on $n$
  variables require $\Omega(2^{n/2})$ auxiliary variables. It follows that some symmetric functions on $N$
  variables must need $\Omega(\sqrt{N})$ auxiliary variables in every quadratization.
\end{proof}

We also have a similar lower bound result for $y$-linear quadratizations. It rests on the following result
from~\cite{ABCG}:

\begin{theorem}\label{thm:lowerbound2}
  There are pseudo-Boolean functions of $n$ variables for which any $y$-linear quadratization
  must involve at least $\Omega(2^n/n)$ auxiliary variables.
\end{theorem}

We then have the following.

\begin{theorem}\label{thm:lb_sym2}
  There exist symmetric functions of $n$ variables for which any $y$-linear quadratization must involve
  at least $\Omega(n/\log n)$ auxiliary variables.
\end{theorem}

\begin{proof}
  The proof is similar to the previous one: it suffices to observe that when $G(z,y)$ is $y$-linear, then so
  is $g(x,y)$, and to rely on the generic lower bound $\Omega(2^n/n)$ of Theorem~\ref{thm:lowerbound2} for the
  number of auxiliary variables required in every $y$-linear quadratization of certain pseudo-Boolean functions.
\end{proof}

Note that the lower bound in Theorem~\ref{thm:lb_sym2} for the number of auxiliary variables in $y$-linear
quadratizations comes within a factor $O(\log n)$ of the upper bound of $n-2$ from Theorem~\ref{thm:upperlinear}.

\subsection{A lower bound for the parity function}
The results just obtained prove the existence of symmetric pseudo-Boolean functions which require a significant
number of auxiliary variables to quadratize. Specifically, there exist functions needing $\Omega(\sqrt{n})$ auxiliary
variables in any quadratization, and functions needing $\Omega(n/\log n)$ auxiliary variables in any $y$-linear
quadratization. Those results do not, however, explicitly exhibit particular such functions. We next give a concrete
example of a function which needs a significant number of auxiliary variables in any $y$-linear quadratization.

\begin{theorem}\label{thm:lb_parity}
  Every $y$-linear quadratization of the parity function on $n$ variables must involve
  at least $\Omega(\sqrt{n})$ auxiliary variables.
\end{theorem}

\begin{proof}
  Let $g(x,y)$ be an arbitrary $y$-linear quadratization of the parity function. Then it can be written as
  \begin{equation}\label{g}
    g(x,y) = q(x) + \sum_{i=1}^m y_i (\ell_i(x) - b_i)
  \end{equation}
  where $q(x)$ is quadratic, and $\ell_1(x),\ldots,\ell_m(x)$ are linear functions of $x$ only.

  For each $i\in [m]=\{1,2,\ldots,m\}$, consider the regions
  \[
    R^+_i= \{x \in \ensuremath{\mathbb{R}}^n : \ell_i(x) \geq b_i \}, \;\;R^-_i= \{x \in \ensuremath{\mathbb{R}}^n : \ell_i(x) \leq b_i \},
  \]
  which are closed half-spaces defined by the linear functions $\ell_i$. For each $S \subseteq [m]$, let $R_S$
  denote the region $R_S= \bigl(\bigcap_{i \in S} R^-_i \bigr) \cap \bigl( \bigcap_{i \not\in S} R^+_i \bigr)$.
  This is one of the `cells' into which the $m$ hyperplanes defining the linear functions $\ell_i$ partition $\ensuremath{\mathbb{R}}^n$.

  On every cell $R_S$, the function $f(x)=\min \{ g(x,y) : y \in \{0,1\}^m \}$ is quadratic. Indeed, on $R(S)$,
  we have
  \[
    \min \bigl\{ g(x,y) : y \in \{0,1\}^m \bigr\} = q(x) + \sum_{i \in S} (\ell_i(x) - b_i).
  \]
  We now use a result from Saks~\cite{Saks93} and Impagliazzo, Paturi and Saks~\cite{IPS} (which was used to
  obtain lower bounds on the size of threshold circuits representing the parity function). Let us say that a
  set of hyperplanes {\em slices} all $r$-dimensional subcubes of the Boolean hypercube $\{0,1\}^n$ if for each
  subcube (or face) of $\{0,1\}^n$ of dimension $r$, there are two vertices of the subcube that lie on opposite
  sides of one of these hyperplanes. Then (Proposition~3.82 of~\cite{Saks93}), if a set of $m$ hyperplanes slices
  all $r$-dimensional subcubes, we have $m>\sqrt{n/(r+1)-1}$. In particular, therefore, any set of hyperplanes
  that slices every $3$-dimensional subcube of $\{0,1\}^n$ must contain more than $\sqrt{n/4 -1}$ planes.
  Suppose the hyperplanes defined by the linear functions $\ell_i$ do {\em not} slice all $3$-dimensional subcubes.
  Then there would be some cell $R_S$ containing a subcube of dimension $3$. The parity function restricted to that
  subcube would then be equal to the quadratic expression $q(x)+\sum_{i \in S}(\ell_i(x)-b_i)$. However, it is
  well-known (see, for instance~\cite{Saks93,MinskyPapert,WangWilliams}) that the parity function on a subcube
  of dimension $r$ cannot be represented as a pseudo-Boolean function of degree less than $r$ (and it cannot even
  be represented as the sign of a pseudo-Boolean function of degree less than $r$). So, we would then have a
  quadratic, degree-$2$, representation of parity on a cube of dimension $3$, which is not possible. It follows,
  therefore, that the set of hyperplanes in question must slice all $3$-dimensional subcubes and therefore has
  size $m >\sqrt{n/4 -1}$.
\end{proof}

\section{Conclusions}
In this paper, we have studied the number of auxiliary variables required in quadratizations (and $y$-linear
quadratizations) of symmetric pseudo-Boolean function. We have presented explicit general constructions of
quadratizations, via special types of representations of the functions. This shows that every such function
can be quadratized (with a $y$-linear quadratization, symmetric in the original variables) using at most $n-2$
auxiliary variables. We investigated in more detail the quadratizations of special functions (monomials,
$t$-out-of-$n$, exact-$t$, and parity functions), where it was possible to obtain quadratizations using significantly
fewer than $n-2$ auxiliary variables. By drawing on a general result from our related paper~\cite{ABCG} and
establishing a connection between quadratizations of general functions and of symmetric functions on a related
number of variables, we showed that there exist symmetric functions requiring $\Omega(\sqrt{n})$ auxiliary
variables in any quadratization, and that $y$-linear quadratization can require $\Omega(n/\log n)$ variables.
It would clearly be of interest to close the gaps between these lower bounds and the linear upper bound. We
established, further, that any $y$-linear representation of the parity function needs $\Omega(\sqrt{n})$
auxiliary variables. An open question is to determine whether a similar (or better) lower bound can be
obtained for any (not necessarily $y$-linear) quadratization of this, or another specific, symmetric function.
For instance, any example of a symmetric function where a non $y$-linear quadratization needs fewer variables
than the $y$-linear ones would be of interest, as would be any non constant lower bound on the number of
auxiliary variables for positive monomials. Furthermore, the number of positive quadratic terms in any known
quadratization of the positive monomial is at least $n-1$, but no lower bound on such quantity has been found
so far. Settling this question is also of great interest, as it is related to the quality of relaxations based
on quadratizations for PBO problems.

\vskip 0.5cm
{\bf Acknowledgements.}
We thank Gy\"{o}rgy Tur\'an for several discussions and references on Boolean circuits for
symmetric functions.
The second author thanks the National Science Foundation (Grant IIS-1161476) for partial support.
The third author was partially funded by the Interuniversity Attraction Poles Programme initiated by
the Belgian Science Policy Office (grant P7/36) and by a sabbatical grant from FNRS.
The fourth author thanks the joint CAPES (Brazil)/Fulbright (USA) fellowship process BEX-2387050/15061676
for partial suport.


\section{Supplementary material}
\subsection{Differential equations}
\label{app:diff}

In this section, we derive the equations of motion for the expectation values of the operators. As an example, we consider the expectation value of $\hat{\tau}_k$ in the Schr\"{o}dinger's picture:
\begin{equation}
 \tau_k(t)= \langle \psi(t)| \hat \tau_k |\psi(t)\rangle,
\end{equation}
where $|\psi(t)\rangle$ denotes the wave function at time $t$. The Schr\"odinger's equation gives:
\begin{equation}
i \frac{d |\psi(t)\rangle}{dt} = \hat H_{eff}(t) |\psi(t)\rangle.
\end{equation}
The equation of motion for $\tau_k(t)$ is then:
\begin{equation}
-i \frac{d \tau_k(t)}{dt}= \langle \psi(t)| \left[\hat H_{eff}(t), \hat \tau_k \right] |\psi(t)\rangle.
\end{equation}
The commutator on the right-hand side consists of six possible operators: $\{\hat \nu_{k\uparrow},\hat \nu_{k\downarrow},\hat n_{k\uparrow},\hat n_{k\downarrow},\hat \tau_k,\hat \sigma_k \}$. By calculating the derivative of their expectation values to $|\psi(t)\rangle$, we get a set of closed-equations for the expectation values of the operators:
\begin{equation}\label{diffequationgroup}
 \begin{split}
 -i \frac{d \nu_{k \uparrow} }{dt} = & 2\epsilon_{k\uparrow} \nu_{k\uparrow} + \Delta^* \sigma_{-k} - \Delta^* \sigma_k + \alpha_k \tau_{-k} + \alpha_k\tau_k \\
-i  \frac{d \nu_{k \downarrow} }{dt} = & 2\epsilon_{k\downarrow} \nu_{k\downarrow} + \Delta^* \sigma^*_{k} - \Delta^* \sigma^*_{-k} +  \alpha_k \tau_{-k} + \alpha_k\tau_k \\
-i \frac{d n_{k\uparrow}}{dt} = & - \Delta^* \tau_k^* + \Delta \tau_k -\alpha_k \sigma_k^* -\alpha_k \sigma_k \\
-i \frac{d n_{k\downarrow}}{dt} = & -\Delta^* \tau_{-k}^* + \Delta \tau_{-k} + \alpha_k \sigma^*_k + \alpha_k \sigma_k \\
-i \frac{d\tau_k}{dt} = & (\epsilon_{k\uparrow}+\epsilon_{k\downarrow}) \tau_k + \Delta^* n_{-k\downarrow} - \alpha_k \nu_{k\downarrow} + \Delta^* n_{k\uparrow} \\ & - \alpha_k \nu_{k\uparrow} - \Delta^* \\
-i \frac{d\sigma_k}{dt} = & (\epsilon_{k\uparrow}-\epsilon_{k\downarrow}) \sigma_k - \Delta^* \nu^*_{k\downarrow} - \Delta \nu_{k\uparrow} - \alpha_k n_{k\downarrow} \\ & + \alpha_k n_{k\uparrow}.
 \end{split}
\end{equation}
The pairing gap $\Delta(t)$ in Eq.~(\ref{diffequationgroup}) can be expressed as
\begin{equation}\label{defdeltaauxiliary}
\Delta(t) = U_f\sum_{k} \tau_k^{\ast}(t),
\end{equation}
which can be solved self-consistently from Eq.~(\ref{diffequationgroup}) and Eq.~(\ref{defdeltaauxiliary}).

\subsection{Calculation of tunneling conductances}
\label{app:method}

To calculate the tunneling conductance, we employ the numerical operator method first introduced in Ref.~\cite{pei12} for the real time Kondo problem, which was then extended to time-dependent systems~\cite{pei13a} and topological superconductors~\cite{pei13b}.

The full Hamiltonian describing the reservoir and the superfluid is expressed as
\begin{equation}
 \hat H_{tot}(t) = \hat H_{eff}(t) + \hat H_{res} + \hat H_V.
\end{equation}
We define the synonyms $\hat d_{j\sigma 1}=\hat c^\dag_{j\sigma}$ and $\hat d_{j\sigma 0} = \hat c_{j\sigma}$ and re-express the conductance as
\begin{equation}\label{defdiffcondapp}
 G(t)= \left( -2 g_V\textbf{Im}\sum_\sigma  \frac{d}{d\mu_R} \langle \hat d_{-1,\sigma, 1}(t) \hat d_{0,\sigma, 0}(t) \rangle \right)_{\mu_R=0}.
\end{equation}
We first calculate the operators $\hat d_{-1,\sigma, 1}(t)$ and $\hat d_{0,\sigma, 0}(t)$ in the Heisenberg picture defined as
\begin{equation}\label{defheisenbergoperator}
\begin{split}
 \hat d_{\alpha} (t)= & \lim_{\tau\to 0} e^{i\hat H_{tot}(t_0) \tau} \cdots e^{i\hat H_{tot} \left(t_{N-1} \right) \tau} \hat d_{\alpha} \\ & \times e^{-i\hat H_{tot} \left(t_{N-1} \right) \tau} \cdots e^{-i\hat H_{tot}(t_0) \tau},
\end{split}
 \end{equation}
where $\alpha$ is the abbreviation of $(j,\sigma, s)$ ($s=0,1$), $\hat H_{tot}$ is the time-dependent Hamiltonian, $N$ is the total number of steps, $t_0=0$ is the initial time, and the intermediate time $t_n=n\tau$, where the time step $\tau=t/N$.

The calculation of Eq.~(\ref{defheisenbergoperator}) is divided into $N$ steps with the time evolving from $t_{N-1}$ to $t_0$. In the {\it n}-th step, we undress the pair of operators $e^{i\hat H_{tot} \left(t_{N-n} \right) \tau} $ and $e^{-i\hat H_{tot} \left(t_{N-n} \right) \tau} $ from $\hat d_\alpha$. Since $\hat H_{tot}$ is quadratic, in principle, the intermediate results at each step can be written as
\begin{equation}\label{heisenbergsolution}
 \hat {\tilde d}_{\alpha} (t_n) = \sum_{\alpha'} W_{\alpha,\alpha'}(t_n) \hat d_{\alpha'},
\end{equation}
with $\hat {\tilde d}_{\alpha} (t_n) $ defined as
\begin{equation}\label{defintermediate}
\begin{split}
 \hat {\tilde d}_{\alpha} (t_n) = & e^{i\hat H_{tot}(t_{N-n}) \tau} \cdots e^{i\hat H_{tot} \left(t_{N-1} \right) \tau} \hat d_{\alpha} \\ & \times e^{-i\hat H_{tot} \left(t_{N-1} \right) \tau} \cdots e^{-i\hat H_{tot}(t_{N-n}) \tau},
\end{split}
 \end{equation}
 and $W_{\alpha,\alpha'}(t_n)$ denoting the propagator. At $n=N$, we get
\begin{equation}
\hat d_{\alpha} (t)=\sum_{\alpha'} W_{\alpha,\alpha'}(t_N) \hat d_{\alpha'}.
\end{equation}

To calculate the propagator $W_{\alpha,\alpha'}(t_N)$, we derive an iterative relation. According to Eq.~(\ref{defintermediate}), we have
\begin{equation}
\hat {\tilde d}_{\alpha} (t_n) =e^{i\hat H_{tot}(t_{N-n}) \tau} \hat {\tilde d}_{\alpha} (t_{n-1}) e^{-i\hat H_{tot}(t_{N-n}) \tau} .
\end{equation}
Substituting Eq.~(\ref{heisenbergsolution}) into the equation above, we find
\begin{equation}\label{iterativerelation}
\begin{split}
& \sum_{\alpha'} W_{\alpha,\alpha'}(t_n) \hat d_{\alpha'} \\ &
= \sum_{\alpha'} W_{\alpha,\alpha'}(t_{n-1})  \bigg( \hat d_{\alpha'}+ i\tau [\hat H_{tot}(t_{N-n}), \hat d_{\alpha'}] \\ & + \frac{(i\tau)^2}{2} [\hat H_{tot}(t_{N-n}),[\hat H_{tot}(t_{N-n}),\hat d_{\alpha'}]] + O(\tau^3) \bigg).
\end{split}
\end{equation}
The calculation of these commutators is straightforward with the results
\begin{equation}\label{commutatorexp}
[\hat H_{tot}(t_j), \hat d_{\alpha}]= \sum_{\beta} G_{\alpha,\beta}(t_j) \hat d_{\beta}.
\end{equation}
Substituting Eq.~(\ref{commutatorexp}) into Eq.~(\ref{iterativerelation}) and comparing the left and the right sides, we obtain the iterative relation
\begin{equation}\label{witerativerelation}
\begin{split}
 W_{\alpha,\alpha'}(t_n)= & W_{\alpha,\alpha'}(t_{n-1}) +i\tau \sum_{\beta} W_{\alpha,\beta}(t_{n-1}) G_{\beta,\alpha'}(t_{N-n}) \\ & -\frac{\tau^2}{2} \sum_{\beta,\beta'} W_{\alpha,\beta}(t_{n-1}) G_{\beta,\beta'}(t_{N-n}) G_{\beta',\alpha'}(t_{N-n})  .
\end{split}
\end{equation}
Here we drop terms of the order $O(\tau^3)$. In principle, the error caused by the discretization of time can be made arbitrarily small by letting $\tau\to 0$. In practice, we adaptively choose $\tau$ to make the discretization error negligible.

According to Eq.~(\ref{witerativerelation}), we start from $n=0$ when $W_{\alpha,\alpha'}(t_0)= \delta_{\alpha,\alpha'}$, and iteratively work out $W_{\alpha,\alpha'}(t_n)$ until $t_N=t$. This approach typically works well, due to the critical fact that the matrix $G_{\alpha,\beta}$ has only very few non-zero elements, such that the number of non-zero $W_{\alpha,\beta}$ which needs to be kept at each step increases slowly according to Eq.~(\ref{witerativerelation}). However, when $N$ is very large, which is inevitable since we must choose a small $\tau$ and at the same time a large $t$ to obtain the asymptotic limit of $G(t)$, there will be too many non-zero $W_{\alpha,\beta}$ which cannot be simultaneously kept. To work around this, we apply a truncation scheme and keep only a fixed number of non-zero propagators (the number is denoted by $M$) with the largest magnitudes at each step. This truncation scheme is necessary for obtaining the steady limit of $G(t)$. We decide the value of $M$ adaptively, i.e.,
set an original $M$ and
increase it until the desired precision is obtained. The logic in doing so is that the truncation error goes to zero in the limit $M\to \infty$.

Substituting the expression of $\hat d_{\alpha}(t)$ into Eq.~(\ref{defdiffcondapp}), we have
\begin{equation}
\begin{split}
 G(t)= & -2 g_V\textbf{Im}\sum_{\sigma\beta\beta'}  W_{-1\sigma 1,\beta}(t) W_{0\sigma 0,\beta'}(t) \\ & \times \left( \frac{d}{d\mu_R} \langle \hat d_{\beta} \hat d_{\beta'} \rangle \right)_{\mu_R=0}.
\end{split}
\end{equation}
Next we calculate the derivative of correlation functions. The correlation functions of the initial state are non-zero only if the two field operators are both in the superfluid or both in the reservoir, since the superfluid and the reservoir are decoupled at the initial time. Considering that the initial state of the superfluid located on sites from $0$ to $L-1$ is independent of $\mu_R$, we immediately get
\begin{equation}
 \frac{d}{d\mu_R} \langle \hat d_{j\sigma s} \hat d_{j'\sigma' s'} \rangle =  0,
\end{equation}
as $j,j' \geq 0$. As $j,j'< 0$, we notice that the reservoir is in the ground state of a free Fermi gas at the chemical potential $\mu_R$, and get
\begin{eqnarray}
\left\{ \begin{array}{c} \langle \hat d_{j\sigma 1} \hat d_{j'\sigma' 0} \rangle = \delta_{\sigma,\sigma'} \displaystyle \frac{\sin((j-j')\theta )}{\pi(j-j')} \\ \langle \hat d_{j\sigma 0} \hat d_{j'\sigma' 1} \rangle = \delta_{\sigma,\sigma'}\left( \delta_{j,j'}-\displaystyle \frac{\sin((j-j')\theta )}{\pi(j-j')} \right) \end{array} \right.
\end{eqnarray}
with $\theta=\arccos(-\mu_R/2g_l)$ and $g_l$ the hopping between two neighbor sites of the reservoir. Finally we find
\begin{eqnarray}
\begin{split}
\left(\displaystyle\frac{d}{d\mu_R} \langle \hat d_{j\sigma s} \hat d_{j'\sigma' s'} \rangle\right)_{\mu_R=0} =& \left[ (1-\delta_{s,s'})(2\delta_{s,1}-1) \right] \delta_{\sigma,\sigma'} \\ & \times \frac{\cos \left( (j-j')\displaystyle\frac{\pi}{2} \right)}{2g_l \pi}.
\end{split}
\end{eqnarray}
\end{document}