\section{Introduction}

What is the origin of short gamma-ray bursts?
This is a long-standing problem for more than 40 years
\citep[e.g.,][]{Nakar07,Berger14}.
A merger of binary neutron stars has been thought to be
the most promising candidate
\citep[e.g.,][]{Paczynski86,Goodman86,Eichler+89},
although an unequivocal evidence was missing.
The detection of the gravitational wave (GW) event GW170817 \citep{GW170817}
and the associated electromagnetic counterparts \citep{GWEM170817}
revolutionized the situation.
In particular, the short gamma-ray burst sGRB 170817A
detected two seconds ($\sim 1.7$ s) after GW170817
\citep{GRB170817A,GRB170817A_GBM,GRB170817A_INT}
and the following afterglows in radio to X-ray
\citep[][]{Troja+17,Margutti+17,Haggard+17,Hallinan+17,Alexander+17,Lyman+18}
give the first direct clues to this problem.

The problem has not been solved yet
because sGRB 170817A was very weak
with an isotropic-equivalent energy
$E_{\gamma,\rm iso} \sim 5.35 \times 10^{46}$ erg,
which is many orders of magnitude smaller than ordinary values.
On the other hand,
the afterglow observations, in particular
of superluminal motion in radio \citep{Mooley+18b,Ghirlanda+18}
and the consistency between the spectral index
and the light curve slope after the luminosity peak
\citep{Troja+18b,Mooley+18c,Lamb+19},
strongly suggest that a relativistic jet is launched and
successfully breaks out the merger ejecta in this event
\citep{Nagakura+14,Murguia-Berthier+14}.

\citet{IN18} proposed that sGRB 170817A is faint because
the jet is off-axis to our line-of-sight
\citep[see also][]{GRB170817A,Granot+17,Granot+18,LK18}.
Most emission is beamed into the on-axis direction
via a relativistic effect
and an off-axis observer receives photons emitted outside the beaming cone.
Consequently the apparent energy of the off-axis jet becomes faint
\citep{IN01,Yamazaki+02,Yamazaki+18}.
The off-axis model is initially studied by using a top-hat jet
with uniform brightness and a sharp edge.
This zeroth-order approximation is useful to capture the essential features
of the off-axis model.
Taking the effect of a finite opening angle of the jet,
a typical jet with a typical viewing angle is found to be broadly
consistent with the observations.

However the simplest off-axis model seems to be difficult to explain
the spectral peak energy $\nu_{\rm peak} = 185 \pm 62$ keV
of the main pulse of GRB 170817A,
if the central value is adopted.
The de-beamed emission from an off-axis top-hat jet tends to have
a low peak energy $\nu_{\rm peak} \sim 10$ keV,
or the peak energy viewed on-axis is required to be very high $>10$ MeV
\citep{Kasliwal+17,IN18,Matsumoto+19}.
More precisely, if the jet is viewed on-axis,
the isotropic energy and the peak energy do not satisfy
the well-known relation, so-called Amati relation for sGRBs
\citep{Amati+02,Yonetoku+04,Tsutsui+13},
which is $\sim 100$ times dimmer than that for long GRBs
at the same peak energy $\nu_{\rm peak}$.
Furthermore the afterglow observations including
VLBI observations of superluminal motion
reveal a jet with energy $E_{\rm iso} > 10^{52}$ erg,
a narrow core $\theta_c \lesssim 5^{\circ}$ and
a viewing angle $\sim 14^{\circ}$--$28^{\circ}$
\citep{Mooley+18b,Ghirlanda+18}.
If we adopt these parameters for an off-axis top-hat jet,
the compactness problem is serious \citep{Matsumoto+19}.\footnote{
  The compactness problem is in principle avoidable
  if the spectral cutoff above the peak energy is extremely sharp,
  which is not known from the observations of sGRB 170817A \citep{IN18}.
  }
In the previous paper, we put the above spectral puzzle aside
partly because the error bar is large and
the puzzle is solved if we allow 3$\sigma$ error \citep{IN18}.
However the large error bar could mainly stem from
the spectral evolution during the burst
and the peak energy of the main pulse could be really high \citep{Veres+18}.
The another reason for neglecting the spectral puzzle is that the weak tail
with 34\% of the fluence of the main pulse
has a temperature $k_B T = 10.3\pm 1.5$ keV,
consistent with the simplest off-axis model.
The main pulse could be produced by a different mechanism,
such as scattering of an sGRB by a cocoon \citep{Kisaka+15,Kisaka+17,Kisaka+18}
and a cocoon breakout from ejecta
\citep{Kasliwal+17,Gottlieb+18,Nakar+18,KI+19}.

The approximation of a top-hat jet would be too simple
to discuss the spectral puzzle \citep{Meszaros+98,ZM02}.
The afterglow observation,
in particular the slowly-rising light curve,
is not consistent with a top-hat jet \citep{Mooley+18a},
but strongly suggests a structured jet
\citep{Troja+18,Ruan+18,Margutti+18,D'Avanzo+18,Lazzati+18,Lyman+18,Troja+18b,Ghirlanda+18,Lamb+19}.
The jet structure most likely affects
the prompt emission and its spectrum as well.

In this paper we reexamine
the spectrum of an off-axis jet in GW170817
taking the jet structure constrained by the afterglow observations into account.
We suggest a possible solution, in which
the off-axis emission from a typical sGRB jet
is compatible with sGRB 170817A.
In Sec.~\ref{sec:structure},
we summarize the jet structure inferred from the afterglow observations.
In Sec.~\ref{sec:formulation},
we generalize the formulation in \citet{IN01} to calculate
the off-axis emission from a structured jet.
In Sec.~\ref{sec:off-center},
we show that the off-axis emission mostly comes from the off-center jet,
neither the jet core nor the line-of-sight jet,
if the jet structure is exponentially faint outward,
and the off-center emission can solve the spectral puzzle.
Section~\ref{sec:summary} is devoted to summary and discussions.

\section{Jet structure from afterglow observations}\label{sec:structure}

Before discussing the off-axis emission,
we summarize the jet structure
inferred from the observations of the afterglows
in Fig.~\ref{fig:profile}. The solid and dashed   lines show
the gamma-ray energy of each model specified by different color
as a function of an angle $\theta$ from the jet center.
The afterglow of sGRB 170817A
was observed after $\sim$ 10 days from the prompt sGRB,
and the solid line part of the jet emission was not observed at $\sim$ 10 days
since only the part with $\theta >  \theta_v-\Gamma^{-1}$ has
a significant contribution to the afterglow radiation,
where $\Gamma$ and $\theta_v$ are the Lorentz factor of the emitting region
and the viewing angle, respectively.
Here we obtain  $\Gamma$
assuming that each fluid element at each angle follows
the spherical hydrodynamical evolution decelerated in the interstellar medium
because numerical simulations suggest that this is a good approximation
\citep{KG03,vanEerten+10}.
Note that the Lorentz factor $\Gamma$ in the afterglow phase is
basically different from that in the prompt phase.

As time goes, the angular size $\Gamma^{-1}$ expands.
Then we can progressively probe the jet structure outside
the angular size observable at $\sim 10$ days,
in particular the inner part close to the jet axis
because the inner jet dominates the energy and the flux.
However we can not infer the structure of the dashed region
because the emission from this region is integrated in the initial observation
at $\sim$ 10 days.
Thus the dashed lines are
the part unconstrained by the afterglow observations.
The observations of the afterglow start at $\sim 10$ days
and there is no information before that.
At this time, $\sim 10$ days, we are already observing
an angular size of $\sim 1/\Gamma$ around the line-of-sight in the jet.

In Fig.~\ref{fig:profile},
we assume a radiative efficiency $\epsilon_{\gamma}=10\%$,
which is typical for the observed sGRBs \citep[e.g.,][]{Fong+15}.
These structures are obtained
by assuming a functional form of the jet structure,
calculating the afterglow emission with the standard afterglow theory,
and fitting the afterglow observations
with model parameters such as the ambient density $n$,
the viewing angle of the jet $\theta_v$,
the total energy, and so on.
The different structures obtained by different authors
most likely reflect the difference of the assumed structure
and model parameters,
which are not completely determined
only by the light curve observations.

\begin{figure}
  \begin{center}
    \includegraphics[width=\linewidth]{fig1.pdf}
  \end{center}
  \caption{
    The jet structures
    inferred from the observations of the afterglow light curves
    ({\it solid lines}) are plotted,
    where a radiative efficiency $\epsilon_{\gamma}=10\%$ is assumed.
    The {\it dashed lines} show the part unconstrained
    by the afterglow observations.
    The isotropic energy of sGRB 170817A is also plotted
    at the viewing angles of the jet ({\it cross marks}).
    See the text for details.
  }
  \label{fig:profile}
\end{figure}

In Fig.~\ref{fig:profile}, we also plot
the isotropic energy of sGRB 170817A
at the viewing angle of the jet ({\it cross marks})
to compare it with the jet structure.
Note that the VLBI observations of superluminal motion
prefer a viewing angle of $14^{\circ} \lesssim \theta_v \lesssim 28^{\circ}$
\citep{Mooley+18b}.
Note also that in blue, cyan and purple line cases,
the gamma-ray energy exceeds
that of sGRB 170817A,
requiring smaller radiative efficiency than $\epsilon_{\gamma}=10\%$
at the viewing angle.

From Fig.~\ref{fig:profile},
we can find that
the central part should be much more energetic than
the observed sGRB 170817A,
regardless of the different structures obtained by the different authors.
In order not to exceed sGRB 170817A,
the isotropic gamma-ray energy of the jet
should decrease exponentially outward
(where it is not always Gaussian but
a sharply decreasing function).
This is a general property required from the afterglow and sGRB 170817A.
Therefore we adopt a fiducial case as
\begin{eqnarray}
  E_{\gamma}(\theta)=\epsilon_{\gamma} E_0 \exp(-\theta^2/2\theta_c^2),
  \label{eq:stjet}
\end{eqnarray}
with $E_0=10^{52.80}$ erg,
$\theta_c=0.059$,
$n=10^{-2.51}$ cm$^{-3}$,
and $\theta_v=0.38 \approx 22^{\circ}$ \citep{Troja+18b}.
Note that although we do not know whether the structure
reflects that of the jet energy or of the radiative efficiency,
it does not matter to the following discussions.
Note also that although the jet structure could be modified
after the prompt emission,
namely during the propagation in the interstellar medium,
it does not change the above conclusion that
the jet structure is exponentially faint outward.

\section{Formulation of off-axis emission}\label{sec:formulation}

To calculate the off-axis emission from a structured jet,
we generalize the formulation in \citet{IN01} in this section.
We consider an axisymmetric jet for simplicity.
Assuming that the distance to the source $d$
is much larger than the source size,
the observed flux $F_{\nu}$ at a frequency $\nu$
is obtained from volume integration
of the emission coefficient $j_{\nu}$ as
\begin{eqnarray}
  F_{\nu}
  &\simeq&\frac{1}{d^2} \int r^2 dr \sin\theta\, d\theta\, d\phi\, j_{\nu},
  \label{eq:Fnu}
\end{eqnarray}
where the jet has an origin at $r=0$ and an axis at $\theta=0$ in
the spherical coordinate $(r, \theta, \phi)$.
The jet axis has a viewing angle $\theta_v$
from the line-of-sight between the observer and the origin.

The Lorentz transformation of the emission coefficient and frequency
from the lab frame (i.e., source center frame) $j_{\nu}$
to the comoving frame $j'_{\nu'}$ is
\begin{eqnarray}
  j_{\nu}&=&\frac{j'_{\nu'}}{\Gamma^2(1-\beta \cos\theta_{\Delta})^2},
  \label{eq:jnu}

  \\
  \nu&=&\frac{\nu'}{\Gamma (1-\beta \cos\theta_{\Delta})},
\end{eqnarray}
respectively
where we assume that the jet moves in the radial direction and thereby
the angle $\theta_\Delta$ between the jet motion
and the line-of-sight direction is given by
that between the $(\theta,\phi)$ direction and
the line-of-sight direction as
\begin{eqnarray}
  \cos\theta_{\Delta}=\sin\theta \cos\phi \sin\theta_v + \cos\theta \cos\theta_v.
\end{eqnarray}

A single pulse of sGRBs is well approximated by
instantaneous thin-shell emission at time $t_0(\theta)$ and radius $r_0(\theta)$,
\begin{eqnarray}
  j'_{\nu'}=\frac{1}{(4\pi)^2 r^2} E'_{\gamma}(\theta)
  f(\nu',\theta) \delta[r-r_0(\theta)] \delta[t-t_0(\theta)],
  \label{eq:jnu_model}
\end{eqnarray}
where the angular structure of the jet is characterized by
the comoving radiation energy $E'_{\gamma}(\theta)$ [erg].
This is related with the radiation energy $E_{\gamma}(\theta)$
and total energy $E(\theta)$ in the lab frame as
\begin{eqnarray}
  \epsilon_{\gamma} E(\theta)
  =E_{\gamma}(\theta)
  =\Gamma E'_{\gamma}(\theta),
  \label{eq:E'}
\end{eqnarray}
where the Lorentz factor $\Gamma$
and the radiative efficiency $\epsilon_{\gamma}$
also have angular structures in general.
We adopt the spectral shape similar to the so-called Band function
\begin{eqnarray}
  f(\nu',\theta) = \frac{C}{\nu'_0(\theta)}
  \left(\frac{\nu'}{\nu'_0(\theta)}\right)^{1+\alpha_B}
  \left[1+\left(\frac{\nu'}{\nu'_0(\theta)}\right)^2\right]^{\frac{\beta_B-\alpha_B}{2}},
  \label{eq:f(nu',theta)}
\end{eqnarray}
with $\alpha_B\sim -1$ and $\beta_B\sim -2.5$
\citep{Kaneko+06}.
We take the constant $C$ so that $\int d\nu' f(\nu',\theta)=1$.
Note that the following discussions do not depend on
the exact shape of the spectrum as long as it has a peak.

The time in the lab frame $t$ is related with
the observed time $T$ as
\begin{eqnarray}
  t=T+\frac{r}{c}\cos\theta_{\Delta},
  \label{eq:times}
\end{eqnarray}
where the time is measured from the merger time and we neglect the cosmological effect.

The isotropic energy is obtained from Eqs.~(\ref{eq:Fnu}), (\ref{eq:jnu}), (\ref{eq:jnu_model}), (\ref{eq:E'}) and (\ref{eq:times})
by performing the integrals of the delta functions as
\begin{eqnarray}
  E_{\gamma,\rm iso}&=&\int dT \int d\nu\ 4\pi d^2 F_{\nu}
  \nonumber\\
  &=& \frac{1}{4\pi} \int \sin\theta\, d\theta\, d\phi\,
  \frac{E_{\gamma}(\theta)}{\Gamma^4(1-\beta \cos\theta_{\Delta})^3},
  \label{eq:Eiso1}
\end{eqnarray}
where the arbitrary functions $r_0(\theta)$ and $t_0(\theta)$
are integrated out.
We can further perform the $\phi$ integral,
\begin{eqnarray}
  E_{\gamma,\rm iso}=\int \frac{\sin\theta\, d\theta}{2}\,
  E_{\gamma}(\theta)\cdot \mathcal{B}(\theta),
  \label{eq:Eiso2}
\end{eqnarray}
where we call the last part as the beaming term,
\begin{eqnarray}
  \mathcal{B}(\theta)&\equiv &\int_{-\pi}^{\pi} \frac{d\phi}{2\pi}\,
  \frac{1}{\Gamma^4(1-\beta \cos\theta_{\Delta})^3}
  \nonumber\\
  &=&\frac{1}{2\Gamma^4}
  \frac{2\left(1-\beta \cos\theta \cos\theta_v\right)^2
    +\left(\beta \sin\theta \sin\theta_v\right)^2}
       {[1-\beta\cos(\theta_v+\theta)]^{5/2}
         [1-\beta\cos(\theta_v-\theta)]^{5/2}}.
       \label{eq:beam}
\end{eqnarray}
Note that we can explicitly show $E_{\gamma,\rm iso}=E_{\gamma}(\theta)$
if $E_{\gamma}(\theta)$ and $\Gamma(\theta)$ are isotropic
(where we can always put $\theta_v=0$ by changing a coordinate
in the integration).

The surface brightness (i.e., the isotropic energy per solid angle) is given by
\begin{eqnarray}
  \frac{dE_{\gamma,\rm iso}}{d\Omega}=
  \frac{1}{4\pi}\frac{E_{\gamma}(\theta)}{\Gamma^4(1-\beta \cos\theta_{\Delta})^3}.
  \label{eq:surface}
\end{eqnarray}

The spectral peak energy $\nu_{\rm peak}$ corresponds to the energy
at which $\nu dE_{\gamma,\rm iso}/d\nu$ takes a maximum value.
We can show
\begin{eqnarray}
  \frac{dE_{\gamma,\rm iso}}{d\nu}
  = \frac{1}{4\pi} \int \sin\theta\, d\theta\, d\phi\,
  \frac{E_{\gamma}(\theta) f(\nu,\theta,\phi)}{\Gamma^4(1-\beta \cos\theta_{\Delta})^3},
  \label{eq:spec}
\end{eqnarray}
where
\begin{eqnarray}
  f(\nu,\theta,\phi) = \frac{C}{\nu_0(\theta,\phi)}
  \left(\frac{\nu}{\nu_0(\theta,\phi)}\right)^{1+\alpha_B}
  \left[1+\left(\frac{\nu}{\nu_0(\theta,\phi)}\right)^2\right]^{\frac{\beta_B-\alpha_B}{2}},
\end{eqnarray}
and $\nu_0(\theta,\phi)=\nu'_0(\theta)/\Gamma(1-\beta\cos\theta_{\Delta})$.

\section{Off-axis emission comes from off-center}\label{sec:off-center}

\subsection{Dominance of off-axis emission}

In Fig.~\ref{fig:Eiso},
we calculate the off-axis emission ({\it red line}) from
a structured jet in Eq.~(\ref{eq:stjet}) ({\it blue line})
with Eqs.~(\ref{eq:Eiso2}) and (\ref{eq:beam}).
For the calculation we need the Lorentz factor,
which is not well constrained from observations.
As a fiducial Lorentz factor,
we adopt a profile decreasing outward
\begin{eqnarray}
  \Gamma=\frac{\Gamma_{\max}}{1+(\theta/\theta_c)^{\lambda}}.
  \label{eq:gamma}
\end{eqnarray}
We take $\Gamma_{\max}=2000$
since lower limits $\Gamma \gtrsim 1000$
are obtained for some sGRBs like sGRB 090510 detected by {\it Fermi}/LAT.
The index $\lambda$ is used
to match the isotropic energy with that of sGRB 170817A,
and is found to be $\lambda \approx 3.8$
for our fiducial case in Eq.~(\ref{eq:stjet}).
It is always possible to match the observed value.
Note that $\Gamma>1$ for our range of interest.
Even if the Lorentz factor profile is different,
the following discussions are similar as long as
it is smooth enough.

\begin{figure}
  \begin{center}
    \includegraphics[width=\linewidth]{fig2.pdf}
  \end{center}
  \caption{
    The isotropic gamma-ray energy of
    the off-axis emission ({\it red line}) from
    a structured jet in Eq.~(\ref{eq:stjet}) ({\it blue line})
    is plotted as a function of the viewing angle $\theta_v$ of the jet.
    The isotropic gamma-ray energy at the viewing angle of sGRB 170817A
    is also plotted ({\it green cross}).
    We adopt the Lorentz factor profile in Eq.~(\ref{eq:gamma}).
    The off-axis emission always dominates
    the line-of-sight emission in the outer region.
  }
  \label{fig:Eiso}
\end{figure}
As shown in Fig.~\ref{fig:Eiso},
the off-axis emission ({\it red line})
always dominates the line-of-sight emission ({\it blue line})
in the outer region.
This is general
irrespective of the uncertainty of the jet structure for GW170817
because the jet energy should decrease exponentially outward
in order to satisfy both the afterglow observation
(i.e., the large energy at $\theta=0$)
and the prompt sGRB observation (i.e., the small energy at $\theta=\theta_v$)
while the off-axis emission has a power-law profile
$\propto (\theta_v-\theta_c)^{-4}$
\citep{IN18}.
Therefore if sGRB 170817A arises from a jet,
it is most likely off-axis emission, not the line-of-sight emission.

\subsection{Off-center emission}

Which part of the jet makes a major contribution
to the observed off-axis emission?
In Fig.~\ref{fig:surface}, we plot
the surface brightness distribution in Eq.~(\ref{eq:surface})
with the same model parameters as in Fig.~\ref{fig:Eiso}.
It is remarkable that most emission comes from the off-center jet,
neither the jet core nor the line-of-sight jet but the middle.

\begin{figure}
  \begin{center}
    \includegraphics[width=\linewidth]{fig3.pdf}
  \end{center}
  \caption{
    The surface brightness distribution of the jet emission is plotted
    on the $(\theta \cos\phi, \theta \sin\phi)$ plane.
    The model parameters are the same as in Fig.~\ref{fig:Eiso}.
    Most emission comes from the off-center jet,
    neither the jet core nor the line-of-sight jet but the middle.
  }
  \label{fig:surface}
\end{figure}
The off-center emission is also
a general property of the off-axis emission
irrespective of the uncertainty of the adopted jet structure for GW170817.
As we can see from Eq.~(\ref{eq:Eiso2}),
the observed isotropic energy is determined by the product
of the jet structure $E_{\gamma}(\theta)$ and
the beaming term $\mathcal{B}(\theta)$.
These two functions are plotted in Fig.~\ref{fig:LF} ({\it lower panel}).
The jet energy ({\it green line}) should decrease exponentially outward
to satisfy both the observations of the afterglow and sGRB 170817A.
On the other hand the beaming term ({\it blue line}) increases outward
within $\theta \lesssim \theta_v - \Gamma^{-1}$
because the beaming cone approaches the line-of-sight.
The product of the decreasing function ($E_{\gamma}(\theta)$)
and the increasing function ($\mathcal{B}(\theta)$)
makes a peak in the middle,
neither the jet core nor the line-of-sight.\footnote{
  This is analogous to a Gamov peak.}

\begin{figure}
  \begin{center}
    \includegraphics[width=\linewidth]{fig4.pdf}
  \end{center}
  \caption{
    ({\it Lower panel}): The jet energy ($E_{\gamma}(\theta)$; {\it green line}),
    the beaming term in Eq.~(\ref{eq:beam})
    ($\mathcal{B}(\theta)$; {\it blue line})
    and their product ({\it red line}), which determines
    the isotropic energy in Eq.~(\ref{eq:Eiso1}),
    are plotted as a function of the polar angle $\theta$ of the jet.
    The model parameters are the same as in Fig.~\ref{fig:Eiso}.
    The product of the decreasing $E_{\gamma}(\theta)$ and
    the increasing $\mathcal{B}(\theta)$
    makes a peak in the off-center region,
    neither the jet core ($\theta_c=0.059=3.4^{\circ}$)
    nor the line-of-sight jet ($\theta_v=0.38$; {\it orange arrow}).
    ({\it Upper panel}): The Lorentz factor distribution
    in Eq.~(\ref{eq:gamma}) is plotted as a function of $\theta$.
  }
  \label{fig:LF}
\end{figure}
\subsection{Compactness problem}

The off-center emission is important
to avoid the compactness problem.
The isotropic energy of the emitting region is
$\sim 10^{48}$--$10^{49}$ erg
much smaller than that of the jet core $\sim 10^{52}$ erg.
The angular separation
between the emitting region and the line-of-sight
is also $\sim 5^{\circ}$--$10^{\circ}$ much smaller than
the top-hat case $\theta_v-\theta_c \sim 20^{\circ}$.
Taking into account that
the Lorentz factor at the emitting region $\Gamma \sim 10$
(see {\it upper panel} in Fig.~\ref{fig:LF})
is smaller than that of the jet core $\Gamma \sim 10^3$,
we can show that the compactness problem is not crucial
\citep{Matsumoto+19}.\footnote{
  The innermost part of the emitting region
  could suffer from the compactness problem.
}

\subsection{Spectral Amati relation}

The off-center emission could also solve
the seemingly inconsistency with the Amati relation
between the spectral peak energy and the isotropic energy.
Actually the off-center emission does not satisfy the Amati relation
as argued previously.
However the jet core can satisfy the Amati relation.
The point is that the off-center jet is different from the jet core.
The classical sGRBs are emitted by the on-axis jet core
and hence satisfy the Amati relation,
where the off-center region is too faint to observe.
The off-center region is observable only for an off-axis event
like sGRB 170817A,
in which the core emission is severely de-beamed and suppressed.
Usually an off-axis event is difficult to find
because they are faint.
The jet structure is exponentially faint outward
and is more or less top-hat-like.
Thus it is easy to understand that
the off-axis event is difficult to find
without any other signature like GWs.

In Fig.~\ref{fig:spec}, we construct a concrete example
of the off-axis emission where
the on-axis emission from the jet core satisfies the Amati relation
and the off-axis emission from the off-center jet
reproduces sGRB 170817A.
We adjust the comoving spectral energy $\nu'_0(\theta)$
in Eq.~(\ref{eq:f(nu',theta)})
as $\nu'_0(\theta)=0.15\, {\rm keV} [1+(\theta/\theta_c)^{3.4}]$
by changing the power-law index (3.4 here)
with the same model parameters as in Figs.~\ref{fig:Eiso}--\ref{fig:LF}.
Note that the product with Eq.~(\ref{eq:gamma})
  $\Gamma(\theta) \nu'_0(\theta)$ (i.e., lab-frame spectral energy) is
  a slowly decreasing function of $\theta$.
The Amati relation is satisfied
as long as the viewing angle is about the core size $\theta_v \sim \theta_c$,
where $\theta_c=0.059=3.4^{\circ}$ for our fiducial case.
As the viewing angle becomes large,
the spectral relation deviates from the Amati relation
and approaches the point of sGRB 170817A.
If this interpretation is correct,
future observations will find outliers of Amati relation,
although the event rate of the outliers is very small
because the jet structure is exponentially faint outward.

\begin{figure}
  \begin{center}
    \includegraphics[width=\linewidth]{fig5.pdf}
  \end{center}
  \caption{
    An example where
    the on-axis emission from the jet core satisfies the Amati relation
    and the off-axis emission from the off-center jet reproduces sGRB 170817A
    is constructed on the plane of the spectral peak energy $\nu_{\rm peak}$
    and isotropic energy $E_{\gamma,{\rm iso}}$.
    The color of the line shows the viewing angle of the jet $\theta_v$.
    The model is compared with the observations of sGRBs, which satisfies
    the Amati relation ({\it grey solid line}).
    Note that Amati relation for sGRBs is
    $\sim 100$ times dimmer than that for long GRBs at the same peak energy
    $\nu_{\rm peak}$,
    which is determined by the sGRBs used in \citet{Tsutsui+13}
    ({\it red filled squares}).
  }
  \label{fig:spec}
\end{figure}

\section{Summary and Discussions}\label{sec:summary}

In this paper we identify the spectral puzzle
and propose a natural solution
in the off-axis jet model for sGRB 170817A associated with GW170817.
Taking both the afterglow and prompt sGRB 170817A into account,
we first clarify that the jet structure should be exponentially faint outward,
so that the off-axis emission dominates the line-of-sight emission
at a large viewing angle.
Given that,
we generally show that the off-axis emission comes from an off-center jet,
neither the jet core nor the line-of-sight but the middle.
Because the off-center jet is much less energetic
and much closer to the line-of-sight than the jet core,
the compactness problem is avoidable.
In addition, the Amati relation may be violated as in sGRB 170817A
because the off-center jet is different from the jet core
that satisfies the Amati relation.

The off-center jet is too faint to find for ordinary observations
without any other messengers like GWs,
high-energy gamma-rays \citep{Murase+18} and neutrinos \citep{Kimura+18}.
Even if it is detected in gamma-rays, the redshift is not determined
without follow-up observations.
Since the jet structure is exponentially faint outward and hence top-hat-like,
the event rate of detecting the off-center jet is very small.
This explains why there are not so much outliers to the Amati relation.
Conversely the off-axis model predicts outliers to the Amati relation
once many events are observed.
It is reported that there are similar bursts with similar spectrum to
sGRB 170817A among previous sGRBs \citep{Burns+18,Troja+18a,vonKienlin+19}.
They could be off-axis events like sGRB 170817A.

The off-axis model predicts
brighter events at smaller viewing angles
than sGRB 170817A.
Taking the inclination dependence of the GW amplitude into account,
the probability of observing the viewing angle
$\theta_v \lesssim 20^{\circ}$, $15^{\circ}$ and $10^{\circ}$
is $\sim 20\%$, $10\%$ and $5\%$, respectively
\citep{LK17}.
Therefore future multi-messenger observations,
in particular with GWs and gamma-rays,
can test the off-axis jet model in GW170817.
However, the exact event rate of the off-axis event depends on
the outer structure of the jet,
which is not constrained by
the afterglow observations of GW170817
(see {\it dashed lines} in Figs.~\ref{fig:profile} and \ref{fig:Eiso}).
It is necessary to observe earlier afterglows than $\sim 10$ days
to probe the off-center jet that is relevant to the off-axis emission.
In addition the viewing angle of GW170817
still has an uncertainty $\theta_v \sim 14^{\circ}$--$28^{\circ}$
even after the VLBI observations of superluminal motion
because it depends on the jet core size and
the sideway expansion during the propagation in the interstellar medium,
which also depends on the unobserved off-center structure.
It is a future problem to encompass the uncertainty of the event rate.

It is not known how the jet structure is formed
in particular in the off-center region for sGRB 170817A.
There are several possibilities.
First the jet itself could have a structure.
Even if the jet has a weak outer structure,
it is cut off during the jet propagation through the merger ejecta.
However the baryon load could be angular dependent,
leading to a different angular expansion and hence
a jet structure after the jet breakout from the merger ejecta.
Note that baryon may be loaded at the base of the jet launch
as well as during the jet propagation.
Second the off-center region could be a cocoon,
which is produced by the shocked jet and shocked ejecta
during the jet propagation.
Although a typical velocity of the cocoon is sub-relativistic,
an imperfect mixing between the shocked jet and shocked ejecta
make a high-entropy cocoon surrounding a jet,
which could lead to a jet structure.
We did not also consider several effects such as
the temporal evolution,
the photosphere,
the high-energy spectral cutoff
and so on.
These are interesting future problems.

\section*{Acknowledgements}

The authors would like to thank
J.~Granot, H.~Hamidani, K.~Hotokezaka, K.~Kashiyama,
T. Kinugawa, S.~Kisaka, K. Kyutoku,
G.~P.~Lamb, A.~Levinson, T.~Matsumoto,
T.~Piran, G.~Ryan, M. Shibata, K.~Takahashi and R.~Yamazaki
for useful discussions.
This work is partly supported by
JSPS KAKENHI nos. 18H01215, 17H06357, 17H06362, 17H06131, 26287051
(KI) and 15H02087 (TN).

\bibliographystyle{mnras}


\section{Introduction}

The top quark, the most massive elementary particle discovered so far, is playing an important role in testing the standard model (SM) and searching for new physics beyond the SM. At a hadron collider, such as the large hadron collider (LHC), the top quarks can be produced in pairs via strong interaction or in association with a jet or a $W$ boson via weak interactions. The associated production with a $W$ boson offers a particular window to the weak interactions of the top quark and potentially can lead to a direct measurement of the CKM matrix element $V_{tb}$. Besides, $tW$ production is the second largest single top production channel and thus serves as an essential background in search for new physics. So far, the LHC has accumulated a large number of data, based on which the total and differential cross sections of this channel have been measured directly \cite{Aad:2012xca,Aad:2015eto,Aaboud:2016lpj,Aaboud:2017qyi,Chatrchyan:2012zca,Chatrchyan:2014tua,Sirunyan:2018lcp}.

On the other hand, the precise theoretical predictions provide a valid framework in which important information can be extracted from the experimental data. Since the leading-order (LO) cross section of the $W$ boson associated production is proportional to the strong coupling $\alpha_s$ and $|V_{tb}|^2$, it is crucial to understand the value of  $\alpha_s$ in order to extract $V_{tb}$. The inclusion of higher order QCD corrections can help to estimate the factorization and renormalization scale dependence of the cross sections. The next-to-leading order (NLO) QCD corrections have been calculated in refs.~\cite{Giele:1995kr,Zhu:2001hw,Cao:2008af, Kant:2014oha}, and the results including decays of the top quark and the $W$ boson are also available \cite{Campbell:2005bb}. The parton shower effects in this channel have been studied in refs.~\cite{Frixione:2008yi,Re:2010bp,Jezo:2016ujg}.

When considering the higher order corrections for $tW^-$ production, there is one subtlety to deal with. Due to the same $tW^-b$ final state in both the real correction to $tW^-$ and the $t{\bar t}$ production with the decay ${\bar t}\to W^-b$, one needs to find a way to differentiate the two processes or to define the $tW^-$ process properly beyond tree level. There are several proposals on the market \cite{Tait:1999cf, Belyaev:2000me, White:2009yt,Demartin:2016axk}. Here we point out that so far we discuss the $tW^-$ production in the five flavor scheme, i.e., the LO is $gb\to tW^-$. It is also possible to work in the four flavour scheme in which the LO is $gg\to tW^-\bar{b}$ \cite{Frederix:2013gra,Cascioli:2013wga}.

The next-to-next-to-leading order (NNLO) QCD corrections to $tW$ production are not available for now. Though the NNLO $N$-jettiness soft function, one of the ingredients for an NNLO calculation using $N$-jettiness subtraction method, has been computed in refs.~\cite{Li:2016tvb,Li:2018tsq}, the two-loop virtual correction is still the bottleneck because of its dependence on multiple scales. The soft gluon corrections near the threshold have been calculated up to NNNLO based on next-to-next-to-leading logarithms (NNLL) resummation \cite{Kidonakis:2006bu,Kidonakis:2010ux,Kidonakis:2016sjf}, which are considered as an approximation to the full higher order corrections.  The three-loop soft anomalous dimension for $tW$ production was calculated by ref.~\cite{Kidonakis:2019nqa} which can be used to study the full NNNLO threshold effects.

In this paper, we will present the soft gluon corrections to  $tW$ production using the soft-collinear effective theory \cite{Bauer:2000ew,Bauer:2000yr,Bauer:2001ct,Bauer:2001yt,Beneke:2002ph} (see \cite{Becher:2014oda} for a review) which separates the hard contributions with the large momentum transfer and the soft gluon corrections characterized as low energy contributions. Two different definitions of the soft limit are investigated. One is measured by the threshold variable $1-z=1-M_{tW}^2/\hat{s}\to 0$, while the other is given by $s_4=(p_1+p_2-p_t)^2-M_W^2\to 0$. In principle, these two definitions encode the same soft gluon physics in the threshold limits and they only differ by power suppressed corrections. The threshold contributions up to NNLO are obtained  and the resummation is achieved through solving the RG equations of the hard and soft functions. Our results could be taken as an important theoretical input in future experimental analyses.

The paper is organized as follows. In section \ref{sec:factorization}  we show the basic information about the kinematics in this process and the factorization formula of the cross section in the soft limit. The numerical results and relevant discussions are then presented in section \ref{sec:numerical}. We conclude in section \ref{conclu}. The evolution equation of parton distribution functions (PDFs) in the threshold limits and the analytic result of the soft function are given in the appendices.

\section{Factorization and resummation formalism} \label{sec:factorization}

We consider inclusive stable top quark and $W$ boson associated production at the LHC
\begin{align}
    p(P_1) + p(P_2) \to t(p_3) + W^-(p_4) + X(P_X),
\end{align}
where $X$ denotes all the other possible extra radiations in the final states. In the threshold limit at the leading power we only need to consider the partonic channel
\begin{align}
    b(p_1) + g(p_2) \to t(p_3) + W^-(p_4) + X(p_X).
\end{align}
The corresponding LO Feynman diagrams are shown in figure~\ref{fig:lo_dia}.
The partonic kinematic variables are defined to be
\begin{align}\label{kindef}
&\hat{s}=(p_1+p_2)^2,~~\hat{t}_1=(p_1-p_3)^2-m_t^2, ~~ \hat{u}_1=(p_2-p_3)^2-m_t^2,
\nonumber \\
&\hat{t}_1^{\,W}= (p_2-p_4)^2-M_W^2,  ~~~~ \hat{u}_1^W=(p_1-p_4)^2-M_W^2.
\end{align}
The corresponding variables at hadronic level are
\begin{align}
    s=\hat{s}/x_1/x_2 ,~~ t_1 =  \hat{t}_1 /x_1, ~~ u_1 =  \hat{u}_1 /x_2, ~~ t_1^{\, W} =  \hat{t}_1^{\,W}/x_2, ~~ u_1^W =  \hat{u}_1^W /x_1,
\end{align}
where $x_{1,2}$ are the Bjorken scaling variables.

\begin{figure}[t]
 \centering
 \includegraphics[scale=0.4,keepaspectratio=true]{figures/bornxsec.pdf}
 \caption{LO Feynman diagrams for $tW^-$ production.}
 \label{fig:lo_dia}
\end{figure}

In the soft limit or threshold limit, the real emissions are highly constrained, only soft gluons allowed in the final state $X$. This limit would be reached if the invariant mass of the final state $M=\sqrt{(p_3+p_4)^2}$ approaches the initial partonic center-of-mass energy $\sqrt{\hat{s}}$. As a result, the variable  $1-z\equiv 1-M^2/\hat{s}\to 0$ in the threshold limit\footnote{In the central-of-mass frame the energy of the soft radiation is $E_g \approx M(1-z)/2\sqrt{z}$. }, and the perturbative expansion of cross section contains a series of large logarithms $\alpha_s^n\left[\ln^{2n-i} (1-z)/(1-z)\right]_+ ~~(i=1,2,...,2n)$,  which might spoil the convergence of perturbative series. It is our purpose in this work to study the threshold behavior and resum such large logarithms to all orders. Since the soft limit is characterized by the final-state two particles' invariant mass, it is called the pair invariant mass (PIM) scheme. Besides, there is another scheme, called one particle inclusive (1PI) scheme, in which the soft limit is defined by partonic level $s_4\to 0$ with $s_4 \equiv \hat{s}+\hat{t}_1+\hat{u}_1+m_t^2-M_W^2 =(p_4+p_X)^2-M_W^2$. The two schemes measure the soft limit in different ways and the combination of the studies in two schemes provides more complete information on the structure.

In the rest part of this section we briefly show the factorization formula for $tW$ production, which can be derived in a similar way used in the other processes, such as the single top or top quark pair productions \cite{Ahrens:2010zv,Ahrens:2011mw,Zhu:2010mr,Wang:2010ue,Wang:2012dc}.
In the PIM scheme the cross section in the threshold limit can be factorized to a product of the hard and the soft function
\footnote{This factorization is carried out at the leading power of the threshold variable. At next-to-leading power, the threshold factorization becomes more complicated; see the recent paper \cite{Beneke:2018gvs}.}~\cite{Idilbi:2005ky,Becher:2007ty}, which describes the physics at two different scales, i.e., the large hard scale and the small soft scale, respectively.
They contain no large logarithms at their intrinsic scales ($\mu_h$ and $\mu_s$ respectively) as expected, since they depend only on a single scale there.
The resummation of all the large logarithms caused by soft gluon effects is achieved by evolving the two functions from the intrinsic scales to a common factorization scale using their renormalization group (RG) equations.
The RG improved differential cross section can be written as
\begin{align}\label{eq:pim}
    \frac{d^2\sigma^{\rm PIM}}{dM^2  d\cos\theta } = & \frac{\lambda^{1/2}}{32 \pi s M^2} \sum_{ij} \int_\tau^1 \frac{dz}{z} \int_z^1 \frac{dx}{x} f_{i/p}(x, \mu_f) f_{j/p}(z/x, \mu_f) H(\mu_h) U_{\rm PIM}(\mu_h,\mu_s,\mu_f)  \nonumber \\
    & \hspace{-1.2cm}\times \frac{z^{-\eta}}{(1-z)^{1-2\eta}} \tilde s_{\rm PIM} \left( \ln \frac{M^2(1-z)^2}{z\mu_s^2} + \partial_\eta,\mu_s \right)
     \frac{e^{-2\gamma_{\rm E} \eta}}{\Gamma(2\eta)} \Bigg |_{\eta = (C_A+C_F) a_{\gamma^{\rm cusp}} (\mu_s, \mu_f)},
\end{align}
where polar angle of top quark $\theta$ is defined in the center-of-mass frame, $\lambda = (1-m_t^2/\hat{s}-M_W^2/\hat{s} )^2-4 m_t^2 M_W^2/\hat{s}^2 $, $\tau = M^2/s$ and the other kinematic variables have already been defined in eq.~\eqref{kindef}.
For convenience we suppress the dependence of hard function $H$, soft function $\tilde s_{\rm PIM}$ and evolution factor $U_{\rm PIM}$ on the kinematic variables. $\tilde{s}_{\rm PIM}$ denotes the soft function defined in the Laplace space.
The evolution factor embodies RG running from hard and soft scale to factorization scale, which is expressed as
\begin{align}
U_{\rm PIM}(\mu_h, \mu_s, \mu_f) = &\left(\frac{M^2}{\mu_h^2}\right)^{(C_A+ C_F)\, a_{\gamma^{\rm cusp}}(\mu_s,\mu_h)} \exp\Big[ 2(C_A + C_F)S(\mu_h,\mu_s) \nonumber \\
 & \hspace{2.4cm} + \, a_{\gamma^h}(\mu_s,\mu_h)  + 2\,a_{\gamma^{\phi_q}}(\mu_s,\mu_f) + 2\,a_{\gamma^{\phi_g}}(\mu_s,\mu_f) \Big].
\end{align}
where the definitions of function $S$ and $a_\gamma$ and all relevant anomalous dimensions, e.g., $\gamma^{\phi_q}$, can be found in appendix A of ref.~\cite{Ahrens:2010zv}.
The hard anomalous dimension specific to the $tW$ process is given by~\cite{Ferroglia:2009ep,Ferroglia:2009ii}
\begin{align}
\gamma^h = 2\left( \gamma^Q + \gamma^q + \gamma^g \right) - \frac{C_A}{2} \gamma^{\rm cusp} \ln \frac{m_t^2 (-\hat s)}{\hat u_1^2} - \left( C_F - \frac{C_A}{2}\right) \gamma^{\rm cusp} \ln \frac{m_t^2 (-\hat s)}{\hat t_1^{\, 2}} .
\end{align}

In the 1PI scheme,  the soft radiations are characterized via the threshold variable $s_4 \approx  2 p_4 \cdot k $ with $k$ the sum of all the momenta of the soft final state~\cite{Laenen:1998qw,Becher:2009th,Ahrens:2011mw, Wang:2012dc}.  The analysis is simplified by going to the rest frame of the inclusive final state $W+X$ where $|\Vec{p}_W| = \mathcal{O}(s_4/M_W)$. In this frame the energy of the soft radiation is  $E_g \approx s_4/2 M_W$. In the 1PI scheme, the RG improved cross section is
\begin{align}\label{eq:1pi}
   \frac{d^2\sigma^{\rm 1PI}}{dp_T^2  dy } =& \frac{  1 }{ 16 \pi s} \sum_{ij}\int_{x_1^{\min}}^{1} \frac{dx_1}{x_1} \int_0^{s_4^{\max}}
   \frac{ds_4}{s_4+M_W^2-m_t^2-x_1 t_1}  f_{i/p}(x_1, \mu_f)f_{j/p}(x_2, \mu_f)  \nonumber \\ &
   \times H(\mu_h)  U_{\rm 1PI}(\mu_h,\mu_s,\mu_f) \tilde{s}_{\rm 1PI}(\partial_\eta,\mu_s) \frac{1}{s_4}\left(\frac{s_4}{M_W \mu_s}\right)^{2\eta} \frac{e^{-2\gamma_{\rm E} \eta}}{\Gamma(2\eta)},
\end{align}
where the integration range is defined as $ x_1^{\min} = (M_W^2-m_t^2-u_1)/(s+t_1)$, $s_4^{\rm max} = m_t^2 - M_W^2 + u_1 + x_1(s+t_1)$, and momentum fraction $x_2$ is defined via $x_1$ and $s_4$ as $x_2 = (s_4+M_W^2-m_t^2-x_1 t_1)/(u_1 + x_1 s)$. The Mandelstam variables are related to the top quark's transverse momentum $p_T$ and rapidity $y$ via $t_1=-\sqrt{s}\,m_{\bot}e^{-y}$, and $u_1=-\sqrt{s}\,m_{\bot}e^{y}$, with $m_{\bot}=\sqrt{\smash[b] p_T^2+m_t^2}$.
Here the 1PI evolution factor has the form as
\begin{align}
U_{\rm 1PI}(\mu_h, \mu_s, \mu_f) &=     \left(\frac{M^2}{\mu_h^2}\right)^{(C_A+C_F)a_{\gamma^{\rm cusp}}(\mu_s,\mu_h)}
  \exp\Bigg[ 2(C_A+C_F)S(\mu_h,\mu_s) + \, a_{\gamma^h}(\mu_s,\mu_h) \nonumber \\
 & \hspace{-2cm}   + 2\,a_{\gamma^{\phi_q}}(\mu_s,\mu_f) + 2\,a_{\gamma^{\phi_g}}(\mu_s,\mu_f)  +a_{\gamma^{\rm cusp}}(\mu_s,\mu_f)\Bigg(C_A \ln\frac{M_W^2\mu_s^2}{\left( \hat{t}_1^{\,W}\right)^2}+C_F \ln\frac{M_W^2\mu_s^2}{\left( \hat{u}_1^W\right)^2}\Bigg) \Bigg].
\end{align}
The difference between the RG factors in eq.~(\ref{eq:pim}) and eq.~(\ref{eq:1pi}) arises from the RG equation of the PDFs and soft function. We present the RG evolution of the PDFs in the PIM and 1PI schemes in appendix~\ref{app:pdfs}. The two-loop anomalous dimensions that govern the evolution of the hard and soft functions and thus determine the scale dependent part  are derived from the general structure of the anomalous dimension~\cite{Ferroglia:2009ep, Ferroglia:2009ii}.
The scale independent part of the hard function has been obtained at NLO using modified MadLoop~\cite{Hirschi:2011pa} which makes uses of Ninja~\cite{Peraro:2014cba}, CutTools~\cite{Ossola:2007ax} and OneLOop~\cite{vanHameren:2009dr} packages.
We have computed one-loop soft function analytically, which is shown in appendix~\ref{app:soft}. Combining all the ingredients together,  we have checked the RG invariance
\begin{align}
  \frac{d}{d\ln \mu} \left( f_{i/p}\otimes f_{j/p} \otimes H \otimes S_{\rm PIM, 1PI} \right) = 0
    \label{eq:sclindep}
\end{align}
in both of the kinematic schemes.

The NLO and NNLO leading power contributions are obtained by setting the scales  $\mu_h, \mu_s, \mu_f  $ in Eqs.~(\ref{eq:pim}) and (\ref{eq:1pi})  equal. In this way, for PIM scheme we capture all the threshold logarithms $\alpha_s \left[  \ln^n (1-z)/(1-z)\right]_+$ with $n=1,0$  at NLO and $\alpha_s ^2 \left[\ln^n (1-z)/(1-z)\right]_+$ with $n=3,2,1,0$ at NNLO, as well as the scale dependent logarithms predicted by eq.~(\ref{eq:pim}).
The similar procedure can be applied to obtain the threshold enhanced logarithms for 1PI scheme.  In the following calculations  the approximate NNLO (aNNLO) cross section is defined as
\begin{align}
    d\sigma({\rm aNNLO } ) = d\sigma({\rm NNLO\  leading } ) + d\sigma({\rm NLO }) - d\sigma({\rm NLO\  leading } ),
    \label{eq:annlo}
\end{align}
where the NLO power suppressed terms in $1-z$ or $s_4$ have been included to give more precise results. We can also match the resummed prediction to the fixed order result by
\begin{align}
    d\sigma({\rm NLO+NNLL } ) = d\sigma({\rm NNLL } ) + d\sigma({\rm NLO }) - d\sigma({\rm NLO\ leading } )
    \label{eq:nnllnlo}
\end{align}
with the NNLL result given by  eq.~(\ref{eq:pim}) and eq.~(\ref{eq:1pi}).

\section{Numerical results}
\label{sec:numerical}

To perform the numerical calculation,
the input parameters are set as $m_t=173.3$ GeV, $\Gamma_t=1.5$ GeV, $M_W=80.419$ GeV, $\alpha=1/132.5$ and the Fermi-constant $G_F=1.166390 \times 10^{-5} $ GeV$^{-2}$.
For the LO and NLO calculations we use the CT14 LO and NLO PDF sets~\cite{Dulat:2015mca} as provided by the LHAPDF library~\cite{Buckley:2014ana}, respectively. The aNNLO  and resummed  predictions are obtained using CT14 NNLO PDF sets.
For fixed-order calculations  the renormalization scale is set to be the same as the factorization scale.
It is natural to set the default hard scale to be the invariant mass of the top quark and $W$ boson, i.e., $\mu_h=M$, where the hard function contains no large logarithms.
The soft scale is chosen numerically according to the criterion that the perturbative series of the soft function are well behaved~\cite{Becher:2007ty}.
Explicitly, we find that the ratio of the soft and hard scale is $0.3 \sim 0.4$ in the PIM scheme and
$0.3 \sim 0.5$ in the 1PI scheme at the 8 TeV LHC.
The default factorization scale has been chosen to be $\mu_f=M$ and $\mu_f=m_t+M_W$ in the PIM and the 1PI schemes, respectively.
The final scale uncertainties are evaluated by varying these scales by a factor of two independently.

As discussed in the introduction there are several methods to deal with the problem of the interference between the real corrections to $tW$ production and $t\bar{t}$ calculation ~\cite{Tait:1999cf, Belyaev:2000me, Campbell:2005bb,Frixione:2008yi,White:2009yt,Demartin:2016axk} at NLO, such as Diagram Removal (DR), Diagram Subtraction (DS) and b-jet transverse momentum veto.
The differences between these schemes have been discussed a lot before; see e.g. refs.~\cite{Frixione:2008yi,Re:2010bp}.
Since they are only relevant in the power suppressed channels\footnote{The problem of the interference exists only in the $gg\to tW^-\bar{b}$ or $q\bar{q}\to tW^-\bar{b}$ channel which is at subleading power near the threshold.},
we will not repeat the discussion about their difference in this paper.
We notice that in refs.~\cite{Kidonakis:2010ux,Kidonakis:2016sjf} where the higher order threshold corrections were studied, the power corrections as well as certain leading power logarithmic independent terms are not taken into account.
In the rest of this section we only show the predictions  of the process $ p p \to tW^-$. The total cross section for $tW^-$ and $\bar{t}W^+$ can be obtained by doubling the results, as demonstrated in ref.~\cite{Frixione:2008yi}. The NLO cross sections in the b-jet veto scheme are evaluated using MCFM~\cite{Campbell:2005bb} with $p_T^{b-{\rm jet}}<50$ GeV. The cross sections in the DS and DR schemes are calculated by POWHEG-BOX~\cite{Alioli:2010xd,Re:2010bp}.

\begin{table}[t]
    \centering
    \resizebox{\textwidth}{!}{\begin{tabular}{l|c|c|c|c|c|c}
    \hline \hline
         [pb]   & \multicolumn{3}{c|}{PIM  }   & \multicolumn{3}{c}{1PI}  \\
    \hline
       &     b-veto     &  DS   &  DR  &   b-veto     &  DS   &  DR   \\
         \hline
          LO   & \multicolumn{3}{c|}{$6.96_{-6\%}^{+5\%}$   }   & \multicolumn{3}{c}{$7.21^{+5\%}_{-4\%}$}  \\
         \hline
          NLO $bg$   & \multicolumn{3}{c|}{$12.0_{+2\%}^{-6\%}$   }   & \multicolumn{3}{c}{$11.7_{+4\%}^{-8\%}$}  \\
    \hline
            NLO leading   & \multicolumn{3}{c|}{$11.2_{+2\%}^{-5\%}$   }   & \multicolumn{3}{c}{$11.5_{+0\%}^{-3\%}$}  \\
    \hline
     NLO   &  $9.31^{+0\%}_{-1\%}$    & $9.92^{+2\%}_{-2\%}$ &  $10.0^{+2\%}_{-2\%}$ & $9.32^{-1\%}_{+0\%}$    & $10.0^{+1\%}_{-2\%}$   & $10.2^{+2\%}_{-2\%}$
     \\
         \hline
     power corr.  & $-1.87_{-0.25}^{+0.5}$ & $-1.26_{-0.40}^{+0.75}$ & $-1.17_{-0.43}^{+0.76}$  &  $-2.14_{-0.05}^{+0.30}$ &   $-1.46_{-0.21}^{+0.53}$ &  $-1.29_{-0.27}^{+0.63}$
     \\
    \hline \hline
    \end{tabular}}
    \caption{The fixed-order  total cross section and the power corrections  for $tW^-$ production with $\sqrt{s}=8$ TeV. The power corrections are defined as $ d\sigma({\rm NLO }) - d\sigma({\rm NLO\ leading } ) $. The LO results are different in two schemes due to the different choice of the factorization scale. The scale uncertainties are shown.
    }
    \label{tab:nlosigma}
\end{table}

\begin{figure}
    \centering
    \includegraphics[width=0.6\textwidth]{figures/scale.pdf}
    \caption{The factorization scale dependence of the cross sections in the PIM scheme for $tW^{-}$ production with $\sqrt{s}=8$ TeV. The NLO result is obtained in the DS scheme. The plots are shown in the region $1/8<\mu_f/M<4$.}
    \label{fig:PIMscale}
\end{figure}

Before presenting the resummed result, we firstly investigate the contribution of the leading power terms.
From table \ref{tab:nlosigma} we can see that the NLO corrections are sizeable,  enhancing the LO result by $29\% \sim 44\%$ depending on
the different methods to isolate the $tW$ process.
These NLO corrections get contributions from all the $bg$, $gg$ and $qq^\prime$ channels, though at LO only $bg$ channel exists.
Among them,  the $bg$ channel dominates or even surpasses the NLO corrections, as indicated in table \ref{tab:nlosigma} too.
Moreover, the leading power terms of the $bg$ channel can approximate the total  result of the $bg$ channel very well,
the difference being only $7\%$ and $2\%$  in the PIM and 1PI schemes, respectively.
Since the leading power terms can be obtained from the resummed results as discussed above in eq.(\ref{eq:annlo}),
they can be calculated up to higher orders in $\alpha_s$, namely beyond NLO.
These make up a major part of the full NNLO corrections and can be taken as an approximation of the latter.
The quality of the approximation could be estimated by looking at the power corrections.
The NNLO results are still unavailable, so we study the NLO ones which are shown in  table \ref{tab:nlosigma}  as well.
It is ready to see that they are negative and around $-20\%\sim -12\%$ in PIM and $-23\%\sim -13\%$ in 1PI kinematic scheme
depending on the methods to deal with the interference problem.
The contributions of the higher order (in $\alpha_s$) power corrections
can be obtained by calculating the full NNLO QCD corrections or by making use of the next-to-leading power factorization and resummation,
both of which are difficult at the moment and beyond the scope of this paper.

Although the usual way to evaluate the scale uncertainty is to vary the scales by a factor of two, it is also
interesting to investigate the factorization scale dependence in a larger region.
From the figure~\ref{fig:PIMscale}, we can see that  the ratio of the NLO over LO result is insensitive to the factorization scale, always in the region $(1.37,1.53)$, when it is varied from $M/8$ to $4M$. This means that there is no clear choice of the factorization scale to ensure fastest convergence.
Moreover, we find that the $bg$ channel is very sensitive to the factorization scale when it is smaller than $M/2$. In order to avoid such a dependence, we have chosen the default factorization scale at $M$.
It can also be seen that the NLO leading power terms dominate the $bg$ channel over a large region.

\begin{figure}[t]
    \centering
    \includegraphics[width=0.49\textwidth]{figures/PIM8fig.pdf}~
    \includegraphics[width=0.49\textwidth]{figures/PIM13fig.pdf}
    \vspace{-1.5cm}
    \caption{Invariant mass distributions in the PIM scheme for $tW^-$ production.
    In the upper plots, the black lines represent the NLO cross section from $bg$ channel while the blue and red lines are the NLO leading and NNLO leading predictions, respectively. In the bottom plots, we show the ratio of NLO (NNLO) leading over NLO $bg$ by blue (red) lines. }
    \label{fig:PIM_mass}
\end{figure}

\begin{figure}[t]
    \centering
    \includegraphics[width=0.49\textwidth]{figures/OPI8fig.pdf}~
    \includegraphics[width=0.49\textwidth]{figures/OPI13fig.pdf}
    \vspace{-1.5cm}
    \caption{Top quark $p_T$ distributions in the 1PI scheme  for $tW^-$ production. The color scheme is the same as figure~\ref{fig:PIM_mass}.
    }
    \label{fig:1PI_pt}
\end{figure}

Then we turn to the differential cross sections.
We show the $tW$ invariant mass distributions in the PIM scheme in figure \ref{fig:PIM_mass}
and the top quark $p_T$ distributions in the 1PI scheme in figure~\ref{fig:1PI_pt}.
We show results at both the 8 TeV and 13 TeV LHC.
It can be seen that the leading power terms are dominant in all the invariant mass or the top $p_T$ regions, as in the case of total cross sections.
The NNLO leading terms increase the NLO leading cross section by about $10\%$ in  most of the region.

\begin{table}[b]
    \centering
    \begin{tabular}{l|c|c|c|c}
    \hline\hline
     [pb]   & \multicolumn{2}{c|}{PIM  }   & \multicolumn{2}{c}{1PI}  \\
     \hline
    $\sqrt{s}$ & 8 TeV & 13 TeV     &   8 TeV & 13 TeV  \\
     \hline
     LO  &  $7.0^{+5\%}_{-6\%}$ & $22.4^{+5\%}_{-2\%}$ & $7.2^{+5\%}_{-4\%}$  & $22.9^{+3\%}_{-1\%}$
      \\
            \hline
     NLO  &  $9.92^{+2\%}_{-2\%}$ & $32.8^{+1\%}_{-1\%}$ & $10.0^{+1\%}_{-2\%}$  & $33.0^{+1\%}_{-1\%}$
      \\
            \hline
     aNNLO  & $11.6^{+4\%}_{-5\%}$  & $37.1^{+5\%}_{-5\%}$  &  $11.2^{+6\%}_{-6\%}$  & $35.9^{+7\%}_{-6\%}$
      \\
            \hline
     NLO+NNLL & $11.4^{+7\%}_{-7\%}$& $36.7^{+7\%}_{-7\%}$ &$11.7^{+12\%}_{-17\%}$ & $37.3^{+16\%}_{-21\%}$
      \\
                  \hline
     aNNLO/NLO & 1.16 & 1.13 & 1.12 & 1.09
     \\
                  \hline
     (NLO+NNLL)/NLO & 1.15 & 1.12 & 1.17 & 1.13
      \\
     \hline \hline
    \end{tabular}
    \caption{Total cross sections for $tW^-$ production in PIM and 1PI schemes.  The NLO cross sections  are calculated using DS scheme. }
    \label{tab:sigma}
\end{table}

Now we present in table \ref{tab:sigma} the aNNLO and NLO+NNLL result defined in eq.~(\ref{eq:annlo}) and eq.~(\ref{eq:nnllnlo}), respectively.
The NLO+NNLL (aNNLO) predictions increase the NLO total cross section by  $12\%\sim 17\%$ ($9\%\sim 16\%$) depending on the collider energy and the threshold variable, but with larger scale uncertainties.
These large uncertainties are mainly from the variation of the factorization scale $\mu_f$. At first sight, this is unexpected since we have checked the scale independence near the threshold analytically in eq.~(\ref{eq:sclindep}).
However, this is based on the assumption $x_{1,2}\to1$ as discussed in appendix~\ref{app:pdfs}.
When the kinematics is far away from the threshold limit, this assumption is not valid.
The very small scale uncertainties of the NLO results seem like a coincidence because the NLO contributions from $gg$ and $qq^\prime$ channels are negative while the contributions from $bg$ channel are positive. Meanwhile they display an opposite behavior under the  scale variation; see table \ref{tab:nlosigma}.
Our resummed result or its expansion in $\alpha_s$ improves only the result in $bg$ channel.
It would be interesting to investigate whether the scale cancellation among different channels happens  at higher orders.
From table \ref{tab:sigma} we also find that the total cross sections in the PIM and 1PI scheme are compatible.
And the resummed cross sections in PIM kinematics have smaller scale uncertainties.

Lastly, we compare the theoretical results with the measurements of the total cross section for $tW^-$ and $\bar tW^+$ production at the LHC in figure~\ref{fig:measu}. After considering the large experimental uncertainties, the NLO+NNLL predictions are in good  agreement with the data at the 8 TeV and 13 TeV LHC.
We also give the predictions at the 14 TeV LHC.

\begin{figure}[t]
    \centering
    \hspace{-0.5cm}\includegraphics[width=0.7\textwidth]{figures/xsec.pdf}
    \caption{Comparison between measured cross section for $tW^-$ and $\bar{t}W^+$ production  at the LHC~\cite{Chatrchyan:2014tua,Aad:2015eto,Aaboud:2016lpj,Sirunyan:2018lcp} and RG-improved predictions.  }
    \label{fig:measu}
\end{figure}

\section{Conclusion}
\label{conclu}

We have investigated the soft-gluon resummation for $tW$ production in the framework of soft-collinear effective theory.
We considered the two different definitions of the threshold limit, $1-M^2/\hat{s} \to 0 $ and $ s_4 \to 0$, corresponding to the PIM and 1PI kinematic schemes, respectively.
We briefly discussed the factorization and resummation formalism in both kinematic schemes.
In addition, we have calculated the hard function and  soft function at NLO.
Expanding the resummed formula in $\alpha_s$ gives the leading power terms of the fixed-order results.
We found that the NLO leading power contribution is a good approximation to the  $bg$ channel at NLO not only for the total cross sections but also for the differential distributions.
After resumming the soft gluon effects to all orders using renormalization group equation,
we find that the NLO+NNLL results
increase the NLO cross sections by about 15(12)\% in PIM and 17(13)\% in 1PI scheme at the 8(13) TeV LHC, but with large uncertainties which is mostly generated by varying the factorization scale.
We compared with the data at the 8 and 13 TeV LHC and found good agreement within uncertainties.
We provide the prediction for the 14 TeV LHC.

In future, we can obtain more precise predictions for the $tW$ process by including higher order hard and soft functions in the resummation formalism or by calculating the full NNLO corrections.
The latter may be achieved making use of the $N$-jettiness subtraction method~\cite{Gao:2012ja,Gaunt:2015pea,Boughezal:2015dva}.
The NNLO beam function~\cite{Gaunt:2014xga,Gaunt:2014cfa} and $N$-jettiness soft function~\cite{Li:2018tsq} for this process have been computed.
The only missing part is the two-loop hard function, which requires  a huge amount of work. We defer this study to future work.

\section*{Acknowledgements}
We would like to thank Shi Ang Li for the contribution in the early stage of this work.  C.S.L. was supported by the National Nature  Science foundation of China, under Grants No. 11875072.  H.T.L. was supported by the Los Alamos National Laboratory LDRD program. J.W. was  supported  by the BMBF project No.  05H15WOCAA and 05H18WOCA1.


\section{Introduction}

Let $\bm{X}_{1},\dots,\bm{X}_{n}\in\mathbb{X}\subseteq\mathbb{R}^{d}$
be an independent and identically distributed (IID) random sample
from some data generating process (DGP), characterized by some probability
distribution $F$. Let $h:\mathbb{X}^{m}\rightarrow\mathbb{R}$ be
a symmetric function, in the sense that
\[
h\left(\bm{x}_{\pi_{1}\left(1\right)},\dots,\bm{x}_{\pi_{1}\left(m\right)}\right)=h\left(\bm{x}_{\pi_{2}\left(1\right)},\dots,\bm{x}_{\pi_{2}\left(m\right)}\right)\text{,}
\]
for all $\bm{\pi}_{i}^{\top}=\left(\pi_{i}\left(1\right),\dots,\pi_{i}\left(m\right)\right)\in\Pi_{m}$
($i\in\left\{ 1,2\right\} $), where $\Pi_{m}$ is the set of all
permutations of the first $m$ consecutive natural numbers. We say
that $h$ is an order $m$ symmetric kernel. Assuming that the parameter
\[
\theta\left(F\right)=\mathbb{E}_{F}h\left(\bm{X}_{1},\dots,\bm{X}_{m}\right)=\int_{\mathbb{X}}\cdots\int_{\mathbb{X}}h\left(\bm{x}_{1},\dots,\bm{x}_{m}\right)\text{d}F\left(\bm{x}_{1}\right)\dots\text{d}F\left(\bm{x}_{m}\right)\text{,}
\]
exists, we can unbiasedly estimate $\theta\left(F\right)=\theta$
via the so-called U-statistic
\[
U_{n}=U\left(\bm{X}_{1},\dots,\bm{X}_{n}\right)={n \choose m}^{-1}\sum_{\bm{\kappa}\in\mathbb{K}_{m}}h\left(\bm{X}_{\kappa\left(1\right)},\dots,\bm{X}_{\kappa\left(m\right)}\right)\text{,}
\]
where $\bm{\kappa}^{\top}=\left(\kappa\left(1\right),\dots,\kappa\left(m\right)\right)\in\mathbb{K}_{m}$
and $\mathbb{K}_{m}$ is the set of all $n!/\left[\left(n-m\right)!m!\right]$
distinct combinations of $m$ elements from the first $n$ consecutive
natural numbers.

The U-statistics were first studied in the landmark articles of \citet{Halmos:1946aa}
and \citet{Hoeffding:1948aa}. Since their introduction, a significant
body of work has been produced on the topic. Comprehensive treatments
of the topic can be found in \citet[Ch. 5]{Serfling:1980aa}, \citet{Lee:1990aa},
\citet{Koroljuk:1994aa}, and \citet{Bose:2018aa}.

In recent years, concentration inequalities have become an important
research theme in the machine learning, probability, and statistics
research, due to their range of practical and theoretical applications.
The current state of the literature is well-reported in the volumes
of \citet{Ledoux:2001aa}, \citet{Massart:2007aa}, \citet{Dubhashi:2009aa},
\citet{Boucheron:2013aa}, and \citet{Bercu:2015aa}.

The first results regarding the concentration of $U_{n}$ about its
mean value $\theta$ were those established in \citet{Hoeffding:1963aa}.
Assume that $h\in\left[a,b\right]$ (for $a,b\in\mathbb{R}$, such
that $a<b$), and denote the variance of $h$ by $\sigma^{2}=\mathbb{V}_{F}h\left(\bm{X}_{1},\dots,\bm{X}_{m}\right)$.
Then, for any $\epsilon>0$ and $m\le n$, \citet{Hoeffding:1963aa}
proved the one-sided inequalities
\begin{equation}
\Pr\left(U_{n}-\theta\ge\epsilon\right)\le\exp\left(-\frac{2\left\lfloor n/m\right\rfloor \epsilon^{2}}{\left(b-a\right)^{2}}\right)\text{ and}\label{eq: one side hoeffding}
\end{equation}
\begin{align}
\Pr\left(U_{n}-\theta\ge\epsilon\right) & \le\exp\left(-\frac{\left\lfloor n/m\right\rfloor \epsilon^{2}}{2\sigma^{2}+\left(2c/3\right)\epsilon}\right)\text{,}\label{eq: one side bernstein}
\end{align}
where $\left\lfloor z\right\rfloor =\max\left\{ \zeta\in\mathbb{Z}:\zeta\le z\right\} $
is the floor function (cf. \citealp[Prop. 2.3]{Arcones:1993aa}) and
$c=2\max\left\{ \left|a\right|,\left|b\right|\right\} $. It is procedural
to demonstrate that the right-hand sides (RHSs) of (\ref{eq: one side hoeffding})
and (\ref{eq: one side bernstein}) also upper bound $\Pr\left(\theta-U_{n}\right)$,
thus we have the absolute inequalities
\begin{equation}
\Pr\left(\left|U_{n}-\theta\right|\ge\epsilon\right)\le2\exp\left(-\frac{2\left\lfloor n/m\right\rfloor \epsilon^{2}}{\left(b-a\right)^{2}}\right)\text{ and}\label{eq: two side hoeffding}
\end{equation}
\begin{equation}
\Pr\left(\left|U_{n}-\theta\right|\ge\epsilon\right)\le2\exp\left(-\frac{\left\lfloor n/m\right\rfloor \epsilon^{2}}{2\sigma^{2}+\left(2c/3\right)\epsilon}\right)\text{,}\label{eq: two side bernstein}
\end{equation}
since $\Pr\left(\left|U_{n}-\theta\right|\ge\epsilon\right)=\Pr\left(U_{n}-\theta\ge\epsilon\right)+\Pr\left(\theta-U_{n}\ge\epsilon\right)$.

Since the establishment of (\ref{eq: two side hoeffding}) and (\ref{eq: two side bernstein}),
there had been little progress in the derivation of fundamentally
novel bounds for U-statistics. A major contribution in this direction
was due to \citet{Arcones:1995aa}. Assume that $h\in\left[0,1\right]$,
and that $\varsigma^{2}=\mathbb{V}_{F}\mathbb{E}_{F}\left[h\left(X_{1},\dots,X_{m}\right)|X_{1}\right]<\infty$
exists. Then, for any $\epsilon>0$ and $m\le n$, \citet{Arcones:1995aa}
proved that
\begin{equation}
\Pr\left(\left|U_{n}-\theta\right|\ge\epsilon\right)\le4\exp\left(-\frac{\left\lfloor n/m\right\rfloor \epsilon^{2}}{2m\varsigma^{2}+\left(2^{m+3}m^{m-1}+\left[2/3\right]m^{-2}\right)\epsilon}\right)\text{.}\label{eq: arcones}
\end{equation}

In practice, only bounds of $h$ tend to be known, regarding data
from any arbitrary DGP. Thus, without knowledge of the variances $\sigma^{2}$
or $\varsigma^{2}$, the Bernstein-type bounds (\ref{eq: one side bernstein}),
(\ref{eq: two side bernstein}), and (\ref{eq: arcones}) cannot be
used numerically. In such situations, only the Hoeffding-type bounds
(\ref{eq: one side hoeffding}) and (\ref{eq: two side hoeffding})
tend to see practical application.

Let $X_{1},\dots,X_{n}\in\left[a,b\right]$ be independent random
variables. If $\Sigma_{n}^{2}=n^{-1}\sum_{i=1}^{n}\mathbb{V}X_{i}$
and $\bar{X}_{n}=n^{-1}\sum_{i=1}^{n}X_{i}$, then for any $\epsilon>0$,
\citet{Bennett:1962aa} proved the classic Bernstein-type inequalities
\begin{equation}
\Pr\left(\bar{X}_{n}-\mathbb{E}\bar{X}_{n}\ge\epsilon\right)\le\exp\left(-\frac{n\epsilon^{2}}{\Sigma_{n}^{2}/2+\left(2c/3\right)\epsilon}\right)\label{eq: classic bernstein one-sided}
\end{equation}
and
\begin{equation}
\Pr\left(\left|\bar{X}_{n}-\mathbb{E}\bar{X}_{n}\right|\ge\epsilon\right)\le2\exp\left(-\frac{\epsilon^{2}}{\Sigma_{n}^{2}/2+\left(2c/3\right)\epsilon}\right)\text{,}\label{eq: classic bernstein two-sided}
\end{equation}
where $c=2\max\left\{ \left|a\right|,\left|b\right|\right\} $. See
also \citet[Thm. 2.28]{Bercu:2015aa}. Similar to the U-statistics
counterpart, the general inequalities (\ref{eq: classic bernstein one-sided})
and (\ref{eq: classic bernstein two-sided}) require knowledge of
the variance $\Sigma_{n}^{2}$ to be of practical use.

A primary application of concentration inequalities is the construction
of $\left(1-\delta\right)\times100\%$ confidence intervals (CIs)
for some estimator of interest, where $\delta\in\left(0,1\right)$.
Recently, empirical CIs based on the (\ref{eq: classic bernstein one-sided})
and (\ref{eq: classic bernstein two-sided}) forms, where $\Sigma_{n}^{2}$
is replaced by an estimator, have been constructed by \citet{Audibert:2009aa}
and \citet{Maurer:2009aa}. These results allow for the construction
of CIs around the mean of data with known bounds, but unknown variances.
Such bounds were applied to perform variance penalized estimation
in \citet{Maurer:2009aa} and for the analysis of multi-armed bandit
problems in \citet{Audibert:2009aa}. Other examples of applications
include racing-based online model selection \citep{Mnih:2008aa} and
classifier boosting \citep{Shivaswamy:2010aa}.

In \citet{Peel:2010aa}, analogous results to \citet{Audibert:2009aa}
and \citet{Maurer:2009aa} were obtained, for the specific context
of U-statistics. That is, empirical two-sided CIs based on the Bernstein-type
bounds (\ref{eq: two side bernstein}) and (\ref{eq: arcones}) were
constructed, where the variances $\sigma^{2}$ or $\varsigma^{2}$
were replaced by respective empirical estimators. These CIs have useful
applications in racing-based online model selection \citep{Peel:2010aa},
change detection \citep{Sakthithasan:2013aa}, and optimal treatment
allocation \citep{Liang:2018aa}. Since the work of \citet{Peel:2010aa},
\citet{Loh:2013aa} also presented an empirical two-side CI, based
on (\ref{eq: arcones}).

In this paper, we utilize inequalities (\ref{eq: one side hoeffding})
and (\ref{eq: one side bernstein}) in order to construct empirical
one-sided CIs for U-statistics. We demonstrate that these bounds are
able to produce one-sided CIs that are analogous to the two-sided
bounds of \citet{Peel:2010aa}. Furthermore, the two-sided form of
our construction provides tighter bounds than those obtained by \citet{Peel:2010aa}.
Using our constructed CIs, we also demonstrate how one can obtain
empirical CIs for the mean from a recently-derived variance-dependent
improved Hoeffding-type concentration inequality of \citet{Bercu:2015aa}.
Our work can be seen as an extension and refinement of the results
of \citet{Peel:2010aa} and as an addition to the literature on empirical
variance-dependent bounds, as pioneered by \citet{Audibert:2009aa}
and \citet{Maurer:2009aa}. We make numerous comments regarding the
relationship between our work and previously obtained outcomes in
the final section of the article.

The paper proceeds as follows. In Section 2, we present the main results
of the paper. New empirical CIs are derived in Section 3. Concluding
remarks regarding our exposition and results are presented in Section
4.

\section{Main results}

\subsection{Technical preliminaries}

We begin the presentation of our main results by providing a pair
of lemmas that are used throughout the remainder of the paper.
\begin{lem}
[Union bound]\label{lem union bound}For random variables $X$, $Y$,
$Z$,
\[
\Pr\left(X>Z\right)\le\Pr\left(X>Y\right)+\Pr\left(Y>Z\right)\text{.}
\]
\end{lem}
\begin{proof}
For any events $\mathscr{A}$, $\mathscr{B}$, we have the union bound:
\[
\Pr\left(\mathscr{A}\cup\mathscr{B}\right)\le\Pr\left(\mathscr{A}\right)+\Pr\left(\mathscr{B}\right)\text{.}
\]
Consider that $\left(X>Z\right)$ is a subset of $\left(X>Y\right)\cup\left(Y>Z\right)$.
Thus
\[
\Pr\left(X>Z\right)\le\Pr\left(\left(X>Y\right)\cup\left(Y>Z\right)\right)\le\Pr\left(X>Y\right)+\Pr\left(Y>Z\right)\text{.}
\]
\end{proof}
\begin{lem}
[Inequality reversal]\label{lem inequality reversal}Suppose that
$X$ is a random variable, and let $A,B>0$ and $C,D\ge0$, such that
for every $\epsilon>0$,
\[
\Pr\left(X\ge\epsilon\right)\le A\exp\left(-\frac{B\epsilon^{2}}{C+D\epsilon}\right)\text{.}
\]
Then, with probability at least $1-\delta$, we have
\[
X\le\sqrt{\frac{C}{B}\log\frac{A}{\delta}}+\frac{D}{B}\log\frac{A}{\delta}\text{.}
\]
\end{lem}
\begin{proof}
Let $\Pr\left(X\ge\epsilon\right)=\delta$, then we have
\[
\delta\le A\exp\left(-\frac{B\epsilon^{2}}{C+D\epsilon}\right)\text{,}
\]
which we can solve for $\epsilon$, to get
\[
\epsilon\le\frac{1}{2B}\left(D\log\frac{A}{\delta}+\sqrt{D^{2}\log^{2}\frac{A}{\delta}+4BC\log\frac{A}{\delta}}\right)\text{.}
\]
We then apply the square root inequality $\sqrt{x+y}\le\sqrt{x}+\sqrt{y}$
in order to obtain the desired result.
\end{proof}
Although the two preceding results appear elsewhere in the literature,
we include them for the convenience of the reader. Lemma \ref{lem inequality reversal}
appears as Lemma 1 in \citet{Peel:2010aa}.

\subsection{\label{subsec:The-sample-variance}Concentration of the variance}

Let $X_{1},\dots,X_{n}\in\mathbb{R}$ be IID random variables and
note that the symmetric kernel
\[
h\left(x_{1},x_{2}\right)=\left(x_{1}-x_{2}\right)^{2}/2
\]
corresponds to the U-statistic
\[
U_{n}={n \choose 2}^{-1}\sum_{\bm{\kappa}\in\mathbb{K}_{2}}\frac{\left(X_{\kappa\left(1\right)}-X_{\kappa\left(2\right)}\right)^{2}}{2}=\left(n-1\right)^{-1}\sum_{i=1}^{n}\left(X_{i}-\bar{X}_{n}\right)^{2}=S_{n}^{2}\text{,}
\]
which is the unbiased estimator of the variance. That is $\theta=\mathbb{E}_{F}S_{n}^{2}=\mathbb{V}_{F}X$,
where $X$ arises from the same DGP as the sample $X_{1},\dots,X_{n}$.
Furthermore, assume that $X\in\left[0,1\right]$, so that $h\left(X_{1},X_{n}\right)\in\left[0,1/2\right]$.
Using (\ref{eq: one side hoeffding}), we can obtain our first result.
\begin{prop}
\label{prop Hoeffding bounds for the variance} If $X_{1},\dots,X_{n}\in\left[0,1\right]$
are IID, then for any $\epsilon>0$,
\begin{equation}
\Pr\left(S_{n}^{2}-\mathbb{V}_{F}X\ge\epsilon\right)\le\exp\left(-8\left\lfloor n/2\right\rfloor \epsilon^{2}\right)\text{.}\label{eq: var hoeffding 1}
\end{equation}
Alternatively, with probability at least $1-\delta$,
\begin{equation}
S_{n}^{2}-\mathbb{V}_{F}X\le\sqrt{\frac{1}{8\left\lfloor n/2\right\rfloor }\log\frac{1}{\delta}}\text{.}\label{eq: var hoeffding ci 1}
\end{equation}
Both (\ref{eq: var hoeffding 1}) and (\ref{eq: var hoeffding ci 1})
hold with $S_{n}^{2}-\mathbb{V}_{F}X$ replaced by $\mathbb{V}_{F}X-S_{n}^{2}$.
\end{prop}
\begin{proof}
Result (\ref{eq: var hoeffding 1}) is obtained by direct substitution
$S_{n}^{2}$ into (\ref{eq: one side hoeffding}), and (\ref{eq: var hoeffding ci 1})
arises via Lemma \ref{lem inequality reversal}.
\end{proof}
Assume the same hypothesis as Proposition \ref{prop Hoeffding bounds for the variance}.
We consider instead a bound for $S_{n}^{2}$ using (\ref{eq: one side bernstein}).
Since $h\left(X_{1},X_{2}\right)\le1/2$, it is true that
\begin{equation}
\mathbb{V}_{F}S_{n}^{2}=\mathbb{E}_{F}\left[S_{n}^{4}\right]-\mathbb{E}_{F}\left[S_{n}^{2}\right]\le\mathbb{E}_{F}S_{n}^{2}=\mathbb{V}_{F}X\text{,}\label{eq: variance less than mean}
\end{equation}
and thus, by direct substitution into (\ref{eq: one side bernstein}),
we obtain the bound
\[
\Pr\left(S_{n}^{2}-\mathbb{V}_{F}X\ge\epsilon\right)\le\exp\left(-\frac{\left\lfloor n/2\right\rfloor \epsilon^{2}}{2\mathbb{V}_{F}X+\left(2/3\right)\epsilon}\right)\text{,}
\]
since $c=\max\left\{ 0,1\right\} =1$.

Using Lemma \ref{lem inequality reversal}, we have, with probability
at least $1-\delta$,

\[
S_{n}^{2}-\mathbb{V}_{F}X\le\sqrt{\frac{2\mathbb{V}_{F}X}{\left\lfloor n/2\right\rfloor }\log\frac{1}{\delta}}+\frac{2}{3\left\lfloor n/2\right\rfloor }\log\frac{1}{\delta}\text{,}
\]
which we can complete the square to obtain
\begin{equation}
S_{n}^{2}\le\left[\sqrt{\mathbb{V}_{F}X}+\frac{1}{2}\sqrt{\frac{2}{\left\lfloor n/2\right\rfloor }\log\frac{1}{\delta}}\right]^{2}+\frac{1}{6\left\lfloor n/2\right\rfloor }\log\frac{1}{\delta}\text{.}\label{eq: intermediate 1}
\end{equation}
Taking the square root, and applying the square root inequality to
both sides of (\ref{eq: intermediate 1}) yields

\begin{align*}
\sqrt{S_{n}^{2}} & +\sqrt{\frac{1}{6\left\lfloor n/2\right\rfloor }\log\frac{1}{\delta}}\le\sqrt{\mathbb{V}_{F}X}+\sqrt{\frac{2}{4\left\lfloor n/2\right\rfloor }\log\frac{1}{\delta}}\\
 & \le\sqrt{\mathbb{V}_{F}X}+\left(\frac{\sqrt{2}}{2}+\frac{\sqrt{6}}{6}\right)\sqrt{\frac{1}{\left\lfloor n/2\right\rfloor }\log\frac{1}{\delta}}\text{.}
\end{align*}

From (\ref{eq: one side bernstein}), we also have, with probability
at least $1-\delta$,
\[
\mathbb{V}_{F}X-S_{n}^{2}\le\sqrt{\frac{2\mathbb{V}_{F}X}{\left\lfloor n/2\right\rfloor }\log\frac{1}{\delta}}+\frac{1}{3\left\lfloor n/2\right\rfloor }\log\frac{1}{\delta}\text{,}
\]
which rearranges to

\[
\mathbb{V}_{F}X-\sqrt{\frac{2\mathbb{V}_{F}X}{\left\lfloor n/2\right\rfloor }\log\frac{1}{\delta}}-\frac{1}{3\left\lfloor n/2\right\rfloor }\log\frac{1}{\delta}\le S_{n}^{2}
\]
and we can again complete the square to obtain

\[
\left[\sqrt{\mathbb{V}_{F}X}-\frac{1}{2}\sqrt{\frac{2}{\left\lfloor n/2\right\rfloor }\log\frac{1}{\delta}}\right]^{2}\le S_{n}^{2}+\frac{7}{6\left\lfloor n/2\right\rfloor }\log\frac{1}{\delta}\text{.}
\]
By the square root inequality, we have
\begin{align*}
\sqrt{\mathbb{V}_{F}X} & \le\sqrt{S_{n}^{2}}+\left(\frac{\sqrt{2}}{2}+\frac{\sqrt{42}}{6}\right)\sqrt{\frac{1}{\left\lfloor n/2\right\rfloor }\log\frac{1}{\delta}}\text{.}
\end{align*}
We therefore have the following empirical one-sided CIs.
\begin{prop}
\label{prop variance bernstein}If $X_{1},\dots,X_{n}\in\left[0,1\right]$
are IID, then with probability at least $1-\delta$,
\begin{equation}
\sqrt{S_{n}^{2}}\le\sqrt{\mathbb{V}_{F}X}+\left(\frac{\sqrt{2}}{2}+\frac{\sqrt{6}}{6}\right)\sqrt{\frac{1}{\left\lfloor n/2\right\rfloor }\log\frac{1}{\delta}}\text{ and }\label{eq: Prop 2 1}
\end{equation}
\begin{equation}
\sqrt{\mathbb{V}_{F}X}\le\sqrt{S_{n}^{2}}+\left(\frac{\sqrt{2}}{2}+\frac{\sqrt{42}}{6}\right)\sqrt{\frac{1}{\left\lfloor n/2\right\rfloor }\log\frac{1}{\delta}}\text{,}\label{eq: Prop 2 2}
\end{equation}
where $\sqrt{2}/2+\sqrt{6}/6\le1.116$ and $\sqrt{2}/2+\sqrt{42}/6\le1.788$.
\end{prop}

\subsection{Variance of a U-statistic}

Assume that $h\in\left[0,1\right]$. Following \citet{Peel:2010aa},
we introduce a new kernel of order $2m$:
\[
\eta\left(\bm{x}_{1},\dots,\bm{x}_{2m}\right)=\left[h\left(\bm{x}_{1},\dots,\bm{x}_{m}\right)-h\left(\bm{x}_{m+1},\dots,\bm{x}_{2m}\right)\right]^{2}/2\text{.}
\]

We note that $\eta$ can either be symmetric or otherwise. If $\eta$
is not symmetric, then we can define a symmetric version of $\eta$,
in the form
\[
\tilde{\eta}\left(\bm{x}_{1},\dots,\bm{x}_{2m}\right)=\frac{1}{\left(2m\right)!}\sum_{\bm{\pi}\in\Pi_{2m}}\eta\left(\bm{x}_{\pi\left(1\right)},\dots,\bm{x}_{\pi\left(2m\right)}\right)\text{,}
\]
where $\mathbb{E}_{F}\left[\eta\left(\bm{X}_{1},\dots,\bm{X}_{2m}\right)\right]=\mathbb{E}_{F}\left[\tilde{\eta}\left(\bm{X}_{1},\dots,\bm{X}_{2m}\right)\right]$
(cf. \citealp[Ch. 5]{Serfling:1980aa}). Furthermore, we can inspect
that
\begin{align*}
\mathbb{E}_{F}\eta & =\frac{1}{2}\mathbb{E}_{F}\left[h\left(\bm{X}_{1},\dots,\bm{X}_{m}\right)-h\left(\bm{X}_{m+1},\dots,\bm{X}_{2m}\right)\right]^{2}\\
 & =\frac{1}{2}\left[\mathbb{E}_{F}h^{2}-2\mathbb{E}_{F}h\mathbb{E}_{F}h+\mathbb{E}_{F}h^{2}\right]\\
 & =\mathbb{E}_{F}h^{2}-\left[\mathbb{E}_{F}h\right]^{2}=\mathbb{V}_{F}h=\sigma^{2}\text{,}
\end{align*}
and note also that $h\in\left[0,1\right]$ implies $\eta\in\left[0,1/2\right]$.

Consider the U-statistics
\[
W_{n}={n \choose 2m}^{-1}\sum_{\bm{\kappa}\in\mathbb{K}_{2m}}\eta\left(\bm{X}_{\kappa\left(1\right)},\dots,\bm{X}_{\kappa\left(2m\right)}\right)
\]
and
\[
\tilde{W}_{n}={n \choose 2m}^{-1}\sum_{\bm{\kappa}\in\mathbb{K}_{2m}}\tilde{\eta}\left(\bm{X}_{\kappa\left(1\right)},\dots,\bm{X}_{\kappa\left(2m\right)}\right)\text{.}
\]
Using $W_{n}$ and $\tilde{W}_{n}$, we seek to obtain $\left(1-\delta\right)\times100\%$
one-sided CIs for the comparison of the quantities $\mathbb{V}_{F}h$
and $W_{n}$.

Assume that $\eta$ is a symmetric kernel. As in Section \ref{subsec:The-sample-variance},
we may use (\ref{eq: one side hoeffding}) to obtain the following
result, analogous to Proposition \ref{prop Hoeffding bounds for the variance}.
\begin{prop}
\label{prop Hoeffding bounds for the U-stat var} If $\bm{X}_{1},\dots,\bm{X}_{n}\in\mathbb{X}$
are IID, $h\in\left[0,1\right]$, and $\eta$ is symmetric, then for
any $\epsilon>0$,
\begin{equation}
\Pr\left(W_{n}-\sigma^{2}\ge\epsilon\right)\le\exp\left(-8\left\lfloor n/\left(2m\right)\right\rfloor \epsilon^{2}\right)\text{.}\label{eq: var hoeffding U-stat}
\end{equation}
Alternatively, with probability at least $1-\delta$,
\begin{equation}
W_{n}-\sigma^{2}\le\sqrt{\frac{1}{8\left\lfloor n/\left(2m\right)\right\rfloor }\log\frac{1}{\delta}}\text{.}\label{eq: var hoeffding ci U-stat}
\end{equation}
Both (\ref{eq: var hoeffding U-stat}) and (\ref{eq: var hoeffding ci U-stat})
hold with $W_{n}-\sigma^{2}$ replaced by $\sigma^{2}-W_{n}$.
\end{prop}
Unfortunately, we cannot always assume that $\eta$ is symmetric.
However, by definition, $\tilde{\eta}$ is always symmetric and has
the same range as $\eta$. That is, if $h\in\left[0,1\right]$, then
$\tilde{\eta}\in\left[0,1/2\right]$. Thus, we have the following
Proposition.
\begin{prop}
\label{cor Hoeffding bounds for the U-stat asym} If $\bm{X}_{1},\dots,\bm{X}_{n}\in\mathbb{X}$
are IID and $h\in\left[0,1\right]$, then all of the conclusions from
Proposition \ref{prop Hoeffding bounds for the U-stat var} hold with
$W_{n}$ replaced by $\tilde{W}_{n}$.
\end{prop}
In order to obtain Bernstein bound analogs of Proposition \ref{prop Hoeffding bounds for the U-stat var}
and Proposition (\ref{cor Hoeffding bounds for the U-stat asym}),
we require the following versions of inequality \ref{eq: variance less than mean}:
\[
\mathbb{V}_{F}\eta=\mathbb{E}_{F}\left[\eta^{2}\right]-\left[\mathbb{E}_{F}\eta\right]^{2}\le\mathbb{E}_{F}\left[\eta^{2}\right]\le\mathbb{E}_{F}\eta=\sigma^{2}\text{ and }
\]
\[
\mathbb{V}_{F}\tilde{\eta}=\mathbb{E}_{F}\left[\tilde{\eta}^{2}\right]-\left[\mathbb{E}_{F}\tilde{\eta}\right]^{2}\le\mathbb{E}_{F}\left[\tilde{\eta}^{2}\right]\le\mathbb{E}_{F}\tilde{\eta}=\sigma^{2}\text{.}
\]
In the same manner in which Proposition \ref{prop variance bernstein}
was obtained, we may prove the following result.
\begin{prop}
\label{prop symmetric bernstein variance}If $\bm{X}_{1},\dots,\bm{X}_{n}\in\mathbb{X}$
are IID, $h\in\left[0,1\right]$, and $\eta$ is symmetric, then with
probability at least $1-\delta$,
\begin{equation}
\sqrt{W_{n}}\le\sqrt{\sigma^{2}}+\left(\frac{\sqrt{2}}{2}+\frac{\sqrt{6}}{6}\right)\sqrt{\frac{1}{\left\lfloor n/\left(2m\right)\right\rfloor }\log\frac{1}{\delta}}\text{ and }\label{eq: sym bern1}
\end{equation}
\begin{equation}
\sqrt{\sigma^{2}}\le\sqrt{W_{n}}+\left(\frac{\sqrt{2}}{2}+\frac{\sqrt{42}}{6}\right)\sqrt{\frac{1}{\left\lfloor n/\left(2m\right)\right\rfloor }\log\frac{1}{\delta}}\text{,}\label{eq: sym bern2}
\end{equation}
where $\sqrt{2}/2+\sqrt{6}/6\le1.116$ and $\sqrt{2}/2+\sqrt{42}/6\le1.788$.
More generally, (\ref{eq: sym bern1}) and (\ref{eq: sym bern2})
also hold with $W_{n}$ replaced by $\tilde{W}_{n}$.
\end{prop}

\subsection{Empirical confidence intervals}

Assume that $\bm{X}_{1},\dots,\bm{X}_{n}\in\mathbb{X}$ are IID and
$h\in\left[0,1\right]$. Via inequality (\ref{eq: one side bernstein})
and Lemma \ref{lem inequality reversal}, we have
\begin{equation}
U_{n}-\theta>\sqrt{\frac{2\sigma^{2}}{\left\lfloor n/m\right\rfloor }\log\frac{2}{\delta}}+\frac{4}{3\left\lfloor n/m\right\rfloor }\log\frac{2}{\delta}\text{,}\label{eq: delta over 2 union}
\end{equation}
with probability at most $\delta/2$. If $\eta$ is symmetric, then
Proposition \ref{prop Hoeffding bounds for the U-stat var} and Proposition
\ref{cor Hoeffding bounds for the U-stat asym} imply that
\begin{equation}
\sigma^{2}>W_{n}+\sqrt{\frac{1}{8\left\lfloor n/\left(2m\right)\right\rfloor }\log\frac{2}{\delta}}\text{,}\label{eq: hoeffding u stat delta on 2}
\end{equation}
with probability at most $\delta/2$. We may apply Lemma \ref{lem union bound}
along with the square root inequality in order to prove that
\begin{equation}
U_{n}-\theta\le\sqrt{\frac{2W_{n}}{\left\lfloor n/m\right\rfloor }\log\frac{2}{\delta}}+\sqrt{\frac{1}{\left\lfloor n/m\right\rfloor }\sqrt{\frac{1}{2\left\lfloor n/\left(2m\right)\right\rfloor }}\log^{3/2}\left(\frac{2}{\delta}\right)}+\frac{4}{3\left\lfloor n/m\right\rfloor }\log\frac{2}{\delta}\text{,}\label{eq: Hoeffding empirical bound}
\end{equation}
with probability at least $1-\delta$. Inequality (\ref{eq: Hoeffding empirical bound})
also holds with $U_{n}-\theta$ replaced by $\theta-U_{n}$. Furthermore,
(\ref{eq: Hoeffding empirical bound}) also holds with $W_{n}$ replaced
by $\tilde{W}_{n}$.

From Proposition \ref{prop Hoeffding bounds for the U-stat var},
we have
\[
\sqrt{\sigma^{2}}>\sqrt{W_{n}}+\left(\frac{\sqrt{2}}{2}+\frac{\sqrt{42}}{6}\right)\sqrt{\frac{1}{\left\lfloor n/\left(2m\right)\right\rfloor }\log\frac{2}{\delta}}\text{,}
\]
with probability at most $\delta/2$, for symmetric $\eta$. Using
Lemma \ref{lem union bound} in combination with (\ref{eq: delta over 2 union}),
we obtain the bound
\begin{equation}
U_{n}-\theta\le\sqrt{\frac{2W_{n}}{\left\lfloor n/m\right\rfloor }\log\frac{2}{\delta}}+\left(\frac{\sqrt{2}}{2}+\frac{\sqrt{42}}{6}\right)\sqrt{\frac{2}{\left\lfloor n/m\right\rfloor \left\lfloor n/\left(2m\right)\right\rfloor }}\log\frac{2}{\delta}+\frac{4}{3\left\lfloor n/m\right\rfloor }\log\frac{2}{\delta}\text{,}\label{eq: Bernstein empirical bound}
\end{equation}
with probability at least $1-\delta$. Generally, (\ref{eq: Bernstein empirical bound})
also holds when we replace $U_{n}-\theta$ by $\theta-U_{n}$ or $W_{n}$
by $\tilde{W}_{n}$, or both simultaneously. We summarize the results
of this section in the following theorem.
\begin{thm}
If $\bm{X}_{1},\dots,\bm{X}_{n}\in\mathbb{X}$ are IID, $h\in\left[0,1\right]$,
and $\eta$ is symmetric, then inequalities (\ref{eq: Hoeffding empirical bound})
and (\ref{eq: Bernstein empirical bound}) hold with probability at
least $1-\delta$. More generally, both (\ref{eq: Hoeffding empirical bound})
and (\ref{eq: Bernstein empirical bound}) hold with probability at
least $1-\delta$ when $U_{n}-\theta$ is replaced by $\theta-U_{n}$,
when $W_{n}$ is replaced by $\tilde{W}_{n}$, or when both quantities
are substituted, simultaneously.
\end{thm}

\section{Empirical confidence intervals based on an improved Hoeffding inequality}

For independent random variables $X_{1},\dots,X_{n}\in\left[a,b\right]$,
the inequalities
\begin{equation}
\Pr\left(\bar{X}_{n}-\mathbb{E}\bar{X}_{n}\ge\epsilon\right)\le\exp\left(-\frac{2n\epsilon^{2}}{\left(b-a\right)^{2}}\right)\text{ and }\label{eq: classic hoeffding 1side}
\end{equation}
\begin{equation}
\Pr\left(\left|\bar{X}_{n}-\mathbb{E}\bar{X}_{n}\right|\ge\epsilon\right)\le2\exp\left(-\frac{2n\epsilon^{2}}{\left(b-a\right)^{2}}\right)\label{eq: classic hoeffding 2side}
\end{equation}
were proved, for any $\epsilon>0$, in \citet{Hoeffding:1963aa}.
In \citet{Bercu:2015aa}, an interesting improvement to the Hoeffding
inequality of form (\ref{eq: classic hoeffding 1side}) was reported.
We present the IID expectation form of the inequality below, and note
that the more general summation form appears as Theorem 2.47 in \citet{Bercu:2015aa}.
\begin{thm}
[Bercu et al., 2015] \label{thm bercu}If $X_{1},\dots,X_{n}\in\left[a,b\right]$
are IID random variables, then
\begin{equation}
\Pr\left(\bar{X}_{n}-\mathbb{E}_{F}X\ge\epsilon\right)\le\exp\left(-\frac{3n\epsilon^{2}}{\left(b-a\right)^{2}+2\mathbb{V}_{F}X}\right)\text{.}\label{eq: bercu 2015}
\end{equation}
Furthermore, (\ref{eq: bercu 2015}) also holds when $\bar{X}_{n}-\mathbb{E}_{F}X$
is replaced by $\mathbb{E}_{F}X-\bar{X}_{n}$.
\end{thm}
Without loss of generality, suppose that $\left[a,b\right]=\left[0,1\right]$.
Then, we obtain
\[
\Pr\left(\bar{X}_{n}-\mathbb{E}_{F}X\ge\epsilon\right)\le\exp\left(-\frac{3n\epsilon^{2}}{1+2\mathbb{V}_{F}X}\right)
\]
and, with probability at least $1-\delta$,

\[
\bar{X}_{n}-\mathbb{E}_{F}X\le\sqrt{\frac{1+2\mathbb{V}_{F}X}{3n}\log\frac{1}{\delta}}\text{,}
\]
via Lemma \ref{lem inequality reversal}. Thus, with probability at
most $\delta/2$,
\begin{equation}
\bar{X}_{n}-\mathbb{E}_{F}X>\sqrt{\frac{1+2\mathbb{V}_{F}X}{3n}\log\frac{2}{\delta}}\text{.}\label{eq: improved hoeffding delta on 2}
\end{equation}
Similar to (\ref{eq: hoeffding u stat delta on 2}), we have
\begin{equation}
\mathbb{V}_{F}X>S_{n}^{2}+\sqrt{\frac{1}{8\left\lfloor n/2\right\rfloor }\log\frac{2}{\delta}}\label{eq: hoeffding var delta on 2}
\end{equation}
with probability at most $\delta/2$, via (\ref{eq: var hoeffding ci 1}).

Combining (\ref{eq: improved hoeffding delta on 2}) and (\ref{eq: hoeffding var delta on 2})
via Lemma \ref{lem union bound} and the square root inequality then
yields
\begin{equation}
\bar{X}_{n}-\mathbb{E}_{F}X\le\sqrt{\frac{1+2S_{n}^{2}}{3n}\log\frac{2}{\delta}}+\sqrt{\frac{1}{12n}\sqrt{\frac{8}{\left\lfloor n/2\right\rfloor }}\log^{3/2}\left(\frac{2}{\delta}\right)}\text{,}\label{eq: empirical improved hoeffding-hoeffding}
\end{equation}
with probability at least $1-\delta$. Given the symmetry of (\ref{eq: bercu 2015}),
we may also switch $\bar{X}_{n}-\mathbb{E}_{F}X$ with $\mathbb{E}_{F}X-\bar{X}_{n}$
in (\ref{eq: empirical improved hoeffding-hoeffding}).

Next, (\ref{eq: Prop 2 2}) implies that
\begin{equation}
\sqrt{\mathbb{V}_{F}X}>\sqrt{S_{n}^{2}}+\left(\frac{\sqrt{2}}{2}+\frac{\sqrt{42}}{6}\right)\sqrt{\frac{1}{\left\lfloor n/2\right\rfloor }\log\frac{2}{\delta}}\label{eq: bernstein var delta on 2}
\end{equation}
with probability at most $\delta/2$. Combining with (\ref{eq: bernstein var delta on 2})
via Lemma \ref{lem union bound} and the square root inequality then
yields

\begin{equation}
\bar{X}_{n}-\mathbb{E}_{F}X\le\sqrt{\frac{1}{3n}\log\frac{2}{\delta}}+\sqrt{\frac{2S_{n}^{2}}{3n}\log\frac{2}{\delta}}+\left(\frac{\sqrt{2}}{2}+\frac{\sqrt{42}}{6}\right)\sqrt{\frac{2}{3\left\lfloor n/2\right\rfloor n}}\log\frac{2}{\delta}\text{,}\label{eq: empirical improved hoeffding-bernstein}
\end{equation}
with probability at least $1-\delta$. Again, the inequality (\ref{eq: empirical improved hoeffding-bernstein})
holds if we switch $\bar{X}_{n}-\mathbb{E}_{F}X$ with $\mathbb{E}_{F}X-\bar{X}_{n}$.

Finally, note that we can obtain two-sided versions of (\ref{eq: empirical improved hoeffding-hoeffding})
and (\ref{eq: empirical improved hoeffding-bernstein}) by considering
the union of lower bounds on $\bar{X}_{n}-\mathbb{E}_{F}X$, $\mathbb{E}_{F}X-\bar{X}_{n}$,
and $\sqrt{\mathbb{V}_{F}X}$, simultaneously. We then obtain, with
probability at least $1-\delta$,
\[
\left|\bar{X}_{n}-\mathbb{E}_{F}X\right|\le\sqrt{\frac{1+2S_{n}^{2}}{3n}\log\frac{4}{\delta}}+\sqrt{\frac{1}{12n}\sqrt{\frac{8}{\left\lfloor n/2\right\rfloor }}\log^{3/2}\left(\frac{3}{\delta}\right)}
\]
and
\[
\left|\bar{X}_{n}-\mathbb{E}_{F}X\right|\le\sqrt{\frac{1}{3n}\log\frac{3}{\delta}}+\sqrt{\frac{2S_{n}^{2}}{3n}\log\frac{3}{\delta}}+\left(\frac{\sqrt{2}}{2}+\frac{\sqrt{42}}{6}\right)\sqrt{\frac{2}{3\left\lfloor n/2\right\rfloor n}}\log\frac{3}{\delta}\text{.}
\]

\section{Concluding remarks}
\begin{rem}
Inequalities (\ref{eq: one side bernstein}), (\ref{eq: two side bernstein}),
(\ref{eq: classic bernstein one-sided}), and (\ref{eq: classic bernstein two-sided})
are often stated with either the additional assumption that $\mathbb{E}_{F}h=0$
or that $\mathbb{E}\bar{X}_{n}=0$ (see, e.g., \citealp[Prop. 2.3]{Arcones:1993aa},
\citet{Arcones:1995aa}, and \citealp[Thm. 2.28]{Bercu:2015aa}).
As such, the constant $c$ is usually defined as $\max\left\{ \left|a\right|,\left|b\right|\right\} $.
However we define $c=2\max\left\{ \left|a\right|,\left|b\right|\right\} $,
since we allow for cases where $\mathbb{E}_{F}h\ne0$ or $\mathbb{E}\bar{X}_{n}\ne0$,
and due to the fact that $\left|X-\mathbb{E}X\right|\le\left|X\right|+\left|\mathbb{E}X\right|\le2\max\left\{ \left|a\right|,\left|b\right|\right\} $,
for any $X\in\left[a,b\right]$.
\end{rem}
\begin{rem}
To the best of our knowledge, Propositions \ref{prop Hoeffding bounds for the variance}
and \ref{prop Hoeffding bounds for the U-stat var} and Proposition
\ref{cor Hoeffding bounds for the U-stat asym} do not appear elsewhere
in the literature. Proposition \ref{prop variance bernstein} is a
special case of Proposition \ref{prop symmetric bernstein variance},
where we use the order one kernel $h=x$. Proposition \ref{prop variance bernstein}
can be viewed as a refinement of intermediate results from the proof
of \citet[Thm. 3]{Peel:2010aa}. Our refinements are as follows: we
correctly differentiated the cases where $\eta$is symmetric or where
we must replace it by $\tilde{\eta}$; we constructed our CIs using
the inequality (\ref{eq: one side bernstein}) instead of the two-sided
version (\ref{eq: two side bernstein}), and thus obtained better
constants; we obtained one-sided CIs for both $\sqrt{W_{n}}-\sqrt{\sigma^{2}}$
and $\sqrt{\sigma^{2}}-\sqrt{W_{n}}$ (and when $W_{n}$ is replaced
by $\tilde{W}_{n}$, whereas \citet{Peel:2010aa} only considered
the CI for $\sqrt{\sigma^{2}}-\sqrt{W_{n}}$; and lastly, we make
no assumptions on the divisibility of $n$, whereas \citet{Peel:2010aa}
assumes that $n$ is divisible by both 2 and $m$ in their presentation.
\end{rem}
\begin{rem}
We may compare Proposition \ref{prop variance bernstein} directly
to Theorem 10 of \citet{Maurer:2009aa}, which cannot be obtained
within the U-statistics framework. Under the same conditions as Proposition
\ref{prop variance bernstein}, \citet[Thm. 10]{Maurer:2009aa} states
that, with probability at least $1-\delta$,
\begin{equation}
\sqrt{S_{n}^{2}}\le\sqrt{\mathbb{V}_{F}X}+\sqrt{\frac{2}{n-1}\log\frac{1}{\delta}}\text{ and}\label{eq: Maurer 1}
\end{equation}
\begin{equation}
\sqrt{\mathbb{V}_{F}X}\le\sqrt{S_{n}^{2}}+\sqrt{\frac{2}{n-1}\log\frac{1}{\delta}}\text{.}\label{eq: Maurer 2}
\end{equation}
It is easy to see that the CIs of \citet[Thm. 10]{Maurer:2009aa}
achieve the same rates with respect to $n$ and $\delta$ as Proposition
\ref{prop variance bernstein}. Furthermore, for large $n$, (\ref{eq: Maurer 1})
and (\ref{eq: Maurer 2}) provide tighter bounds than (\ref{eq: Prop 2 1})
and (\ref{eq: Prop 2 2}), respectively, for almost all values of
$n$. We note that the only exception is when $n=2$ or $4$, where
(\ref{eq: Prop 2 1}) is tighter than (\ref{eq: Maurer 1}).
\end{rem}
\begin{rem}
By considering the lower bounds of $U_{n}-\theta$, $\theta-U_{n}$
and $\sqrt{\sigma^{2}}$, simultaneously, we may obtain the two-sided
version of (\ref{eq: Bernstein empirical bound}):
\begin{equation}
\left|U_{n}-\theta\right|\le\sqrt{\frac{2W_{n}}{\left\lfloor n/m\right\rfloor }\log\frac{3}{\delta}}+\left(\frac{\sqrt{2}}{2}+\frac{\sqrt{42}}{6}\right)\sqrt{\frac{2}{\left\lfloor n/m\right\rfloor \left\lfloor n/\left(2m\right)\right\rfloor }}\log\frac{3}{\delta}+\frac{4}{3\left\lfloor n/m\right\rfloor }\log\frac{3}{\delta}\text{,}\label{eq: twosided bernstein Ustat}
\end{equation}
with probability at least $1-\delta$. We may compare this directly
with \citet[Thm. 3]{Peel:2010aa}, which for symmetric $\eta$ and
under the assumption that $n$ is divisible by $2$ and $m$, implies
that
\[
\left|U_{n}-\theta\right|\le\sqrt{\frac{2mW_{n}}{n}\log\frac{4}{\delta}}+\frac{5m}{n}\log\frac{4}{\delta}\text{,}
\]
with probability at least $1-\delta$. Under the same assumptions
regarding the divisibility of $n$, we can write (\ref{eq: twosided bernstein Ustat})
as
\[
\left|U_{n}-\theta\right|\le\sqrt{\frac{2mW_{n}}{n}\log\frac{3}{\delta}}+\frac{\left(4+\sqrt{2}\left[3+\sqrt{21}\right]\right)m}{3n}\log\frac{3}{\delta}\text{,}
\]
where $\left(4+\sqrt{2}\left[3+\sqrt{21}\right]\right)/3\le4.908$.
Thus, (\ref{eq: twosided bernstein Ustat}) is tighter than the CI
of \citet[Thm. 3]{Peel:2010aa}.
\end{rem}
\begin{rem}
As noted in \citet[Sec. 2.5.4]{Bercu:2015aa}, (\ref{eq: bercu 2015})
is a strict improvement of (\ref{eq: classic hoeffding 1side}) since
$\left(b-a\right)^{2}\ge4\mathbb{V}_{F}X$. Furthermore, it is provable
that $\left(b-a\right)^{2}>4\mathbb{V}_{F}X$ for all random variables
$X\in\left[a,b\right]$, except for binary $X\in\left\{ a,b\right\} $,
where $\Pr\left(X=a\right)=\Pr\left(X=b\right)=1/2$.
\end{rem}
\begin{rem}
CIs (\ref{eq: empirical improved hoeffding-hoeffding}) and (\ref{eq: empirical improved hoeffding-bernstein})
may be compared directly to the empirical Bernstein-type inequalities
of \citet{Audibert:2009aa} and and \citet{Maurer:2009aa}. Under
the same conditions as those under which (\ref{eq: empirical improved hoeffding-hoeffding})
and (\ref{eq: empirical improved hoeffding-bernstein}) are established,
with probability at least $1-\delta$, we have
\begin{equation}
\bar{X}_{n}-\mathbb{E}_{F}X\le\sqrt{\frac{2\left(n-1\right)S_{n}^{2}}{n^{2}}\log\left(\frac{2}{\delta}\right)}+\frac{3}{n}\log\left(\frac{2}{\delta}\right)\text{ and}\label{eq: audibert bennett}
\end{equation}
\begin{equation}
\bar{X}_{n}-\mathbb{E}_{F}X\le\sqrt{\frac{2S_{n}^{2}}{n}\log\left(\frac{2}{\delta}\right)}+\frac{7}{3\left(n-1\right)}\log\left(\frac{2}{\delta}\right)\text{,}\label{eq: maurer bennett}
\end{equation}
via \citet[Thm. 1]{Audibert:2009aa} and \citet[Thm. 4]{Maurer:2009aa},
respectively. A visual comparison of the logarithms of the RHSs of
(\ref{eq: empirical improved hoeffding-hoeffding}), (\ref{eq: empirical improved hoeffding-bernstein}),
(\ref{eq: audibert bennett}), and (\ref{eq: maurer bennett}) is
provided in Figure \ref{fig: AudibertMaurer}. We compare the four
bounds for $S_{n}^{2}\in\left\{ 0.05,0.25\right\} $ and $\delta\in\left\{ 0.01,0.1\right\} $.
It is observable that (\ref{eq: empirical improved hoeffding-hoeffding})
was uniformly tighter than (\ref{eq: empirical improved hoeffding-bernstein}).
For smaller values of $S_{n}^{2}$, (\ref{eq: audibert bennett})
and (\ref{eq: maurer bennett}) were tighter than (\ref{eq: empirical improved hoeffding-hoeffding}),
for larger $n$. For larger $S_{n}^{2}$, (\ref{eq: empirical improved hoeffding-hoeffding})
was tighter than (\ref{eq: audibert bennett}) and (\ref{eq: maurer bennett})
over a middle range of $n$, however (\ref{eq: empirical improved hoeffding-bernstein})
remained uncompetitive. Changing $\delta$ did not tend to alter the
relative performance of the bounds. As noted by \citet{Maurer:2009aa},
(\ref{eq: maurer bennett}) has better constants than (\ref{eq: audibert bennett})
and thus provided tighter bounds for larger values of $n$. We finally
note that (\ref{eq: empirical improved hoeffding-bernstein}) can
be improved by using (\ref{eq: Maurer 2}) in the place of (\ref{eq: Prop 2 2})
in its derivation. This was not pursued because we wished for the
derivation to be self-contained within the U-statistics framework.

\begin{figure}
\begin{centering}
\includegraphics[width=16cm]{AudibertMaurerComparison}
\par\end{centering}
\caption{\label{fig: AudibertMaurer} Logarithms of the upper bounds of the
CIs (\ref{eq: empirical improved hoeffding-hoeffding}), (\ref{eq: empirical improved hoeffding-bernstein}),
(\ref{eq: audibert bennett}), and (\ref{eq: maurer bennett}), as
a function of $n$, for $S_{n}^{2}$ and $\delta$ set at various
levels. Each subplot A--D visualizes the relative performances of
the four bounds, for the values $S_{n}^{2}$ and $\delta$ that are
displayed in the title. The labels Improved Hoeffding 1 and 2 correspond
to (\ref{eq: empirical improved hoeffding-hoeffding}) and (\ref{eq: empirical improved hoeffding-bernstein}),
respectively, whereas Audibert and Maurer correspond to (\ref{eq: audibert bennett})
and (\ref{eq: maurer bennett}), respectively.}

\end{figure}
\end{rem}
\begin{rem}
For $k,n\in\mathbb{N}$, such that $k<n$, we have $\left\lfloor n/k\right\rfloor >\left(n-k+1\right)/k$.
Thus, we may remove the floor operator in each of the inequalities
where it appears by upper bounding its multiplicative inverse. For
instance,
\[
\sqrt{\frac{2mW_{n}}{n-m+1}\log\frac{3}{\delta}}+\left(\frac{\sqrt{2}}{2}+\frac{\sqrt{42}}{6}\right)\sqrt{\frac{4m^{2}}{\left(n-m+1\right)\left(n-2m+1\right)}}\log\frac{3}{\delta}+\frac{4m}{3\left(n-m+1\right)}\log\frac{3}{\delta}
\]
is an upper bound for the RHS of (\ref{eq: twosided bernstein Ustat}).
\end{rem}

\section*{Acknowledgements}

The author is funded by Australian Research Council grants DE170101134
and DP180101192.

\bibliographystyle{apalike2}