\section{Introduction}

The purpose of this paper is to prove the following theorem:
\begin{thm}\label{main-thm}
Let $\Omega$ be a bounded domain in $\mathbb{R}^n$ with Lipschitz boundary. Let $V$ be the volume of $\Omega$ and $\delta$ be its packing constant. Let $\lambda_k$ be its $k^{th}$ Dirichlet eigenvalue. Then for all $k\in\mathbb{N}$,
$$ \lambda_k \geq C_n\bigg(\delta\frac{k}{V}\bigg)^{2/n}.$$
\end{thm}

Here is a brief survey of the relevant literature.

Weyl \cite{Weyl1911} proved the asymptotic estimate
\begin{equation}\label{weyl-law}
\lim_{k\to\infty}\frac{\lambda_k}{k^{2/n}} = C_n\bigg(\frac{1}{V}\bigg)^{2/n}
\end{equation}
as $k\to\infty$. Here $C_n = (2\pi)^2/(\mbox{volume of unit $n$-ball})^{2/n}.$

Polya \cite{Polya1961} considered domains which tile Euclidean space by reflection and translation. He used this tiling property and Weyl's law to prove that, for such domains,
\begin{equation}\label{polya}
\lambda_k \geq C_n\bigg(\frac{k}{V}\bigg)^{2/n}
\end{equation}
for every $k\in\mathbb{N}$. He conjectured that this holds for any planar domain. This conjecture is still open.

Urakawa \cite{Urakawa1984} proved a weaker version of Theorem 1, bounding Dirichlet eigenvalues in terms of the \emph{lattice packing constant} of a domain.

Li and Yau \cite{Li1983} proved for an arbitrary domain in $\mathbb{R}^n$ that
\begin{equation}\label{li-yau} \lambda_k \geq \frac{n}{n+2}C_n \bigg(\frac{k}{V}\bigg)^{2/n}
\end{equation}
for every $k\in\mathbb{N}$. They proved this as a corollary of an inequality about the Riesz mean of the eigenvalue sequence, $\frac{1}{k}\sum_{j=1}^k \lambda_j$. Kroger \cite{Kroger1994} proved the corresponding Li-Yau inequality for Neumann eigenvalues.

Recent improvements to inequality \eqref{li-yau} have proceeded by adding terms to the inequality involving Riesz means of eigenvalues and extending the inequality to more general settings. For details, we refer the reader to \cite{Laptev1997}, \cite{Melas2003}, \cite{Wei2010}, \cite{Geisinger2011}, \cite{Hatzinikitas2013}, \cite{Yolcu2013a}, \cite{Yolcu2013}, \cite{Yolcu2014}, and \cite{Kovarik2015}.

Following Polya's original argument and Urakawa's generalization, we use packings of $\Omega$. By doing so, the packing constant $\delta$ of $\Omega$ enters the inequality. In fact, we are able to replace $n/(    n+2)$ in Li-Yau's estimate with $\delta^{2/n}$. This replacement sacrifices universality, but strengthens the inequality for domains with high packing constant.

\section{Discussion}

Theorem 1 is a generalization of Polya's theorem, as the packing constant of a tiling domain is $1$. Theorem 1 also permits us to replace the factor $n/(n+2)$ in inequality \eqref{li-yau} with the factor $\delta^{2/n}$. In particular, if $\delta > [n/(n+2)]^{n/2}$, then the inequality in Theorem 1 is stronger than inequality \eqref{li-yau}. Such domains are not difficult to construct; for instance, Theorem 1 is stronger than inequality \eqref{li-yau} for any domain in dimension $n\geq 2$ which has a bounding parallelopiped with less than twice the volume of the domain.

General lower bounds for various classes of domain in all dimensions tend to be weak. For instance, a theorem of Minkowski-Hlawka guarantees that the packing constant for a convex, centrally symmetric domain in $\mathbb{R}^n$ is no less than $\zeta(n)/2^{n-1}$. Schmidt then proved that there is a constant $c$ such that a convex domain in $\mathbb{R}^n$ has $\delta \geq cn^{3/2}/4^n$. See, for instance, the discussion in \cite{Toth1993}.

However, dimensions $2$ and $3$ are better-studied. As an example, we reproduce a portion of a table from the survey \cite{Bezdek2013}, modifying the last row with information from section 8.4 of the same paper:

\begin{tabular}{|c|c|}
\hline
Body & Lower bound for packing density \\
\hline
Unit ball & $\frac{\pi}{\sqrt{18}} = 0.7408\ldots$ \\
Regular octahedron & $\frac{18}{19} = 0.9473\ldots$ \\
Cylinder over a plane domain $K$ & $\delta(K)$ \\
Doubled cone & $\pi\sqrt{6}/9 = 0.855\ldots$ \\
Tetrahedron & $ 0.856\ldots $ \\
\hline
\end{tabular}

Here a cylinder over a plane domain $K$ is the Minkowski sum of $K\times\{0\}$ with a line segment $s$ (which is assumed non-parallel to $K$). Observe that the packing constant of cylinders implies the three-dimensional case of Laptev's proof of Polya's conjecture in \cite{Laptev1997} for products of tiling domains and arbitrary domains.

In \cite{Torquato2009}, the authors provide a survey of known lattice packing constants for Platonic and Archimedean solids. All the Platonic and Archimedean solids have packing densities in excess of $0.5$.

Specializing to $n=2$, Kuperberg-Kuperberg \cite{Kuperberg1990}, later improved by Doheny \cite{Doheny1995}, found lower bounds for packing constants of convex planar domains:
\begin{thm}
If $\Omega$ is a convex planar domain, then its packing constant is at least $\sqrt{3}/2$.
\end{thm}
This gives the following corollary to Theorem 1.
\begin{cor}
Let $\Omega$ be a convex planar domain. Then for all $k$, its Dirichlet eigenvalues satisfy
$$ \lambda_k > 2\sqrt{3}\pi\frac{k}{V}. $$
\end{cor}
Note that this improves inequality \eqref{li-yau} by a factor of $\sqrt{3}$. (After the author uploaded a first draft of this preprint to the arxiv, Iosif Polterovich informed the author that this result is known to him and Olivier Mercier.)

We also note that in general, packing constants are greater than lattice packing constants. For instance, the regular tetrahedron has a low lattice packing constant and admits non-lattice packings with much, much higher density; see section 8 of \cite{Bezdek2013} for more information. In fact, many domains have high packing constants but relatively low lattice packing constants. This is the case even in $\mathbb{R}^2$; there are families of convex polygons (such as triangles) which tile the plane, but whose lattice packing constants are strictly less than one.

\section{Proof of Result}

Let $\Omega$ be an open, bounded domain in $\mathbb{R}^n$ with volume $V = |\Omega|$. Let $\Delta = -\sum\d_i^2$ denote the Laplace operator in $\Omega$. Denote by
$$0<\lambda_1<\lambda_2\leq\cdots$$
the spectrum of the Laplace operator with Dirichlet boundary conditions.

The proof proceeds in two steps. The first step applies Dirichlet domain monotonicity and uses Weyl's law to prove Theorem 1. The second step equates the limit $\lim_{\epsilon\to 0}N(\epsilon)\epsilon^n$ with the packing constant of $\Omega$.

Now we prove Theorem 1.

\begin{proof}Let $\sigma > 0$ be given. Set $G = [-\sigma/2,\sigma/2]^n$. Note that $|G| = \sigma^n$. Denote by $N_G$ and $N_\Omega$ the eigenvalue counting functions of $G$ and $\Omega$, resp.

Let $\mathcal{P}$ be a maximal packing of $\Omega$. (This exists by work of Groemer, \cite{Groemer1986}.) Let $N(\sigma)$ be the number of components of $\mathcal{P}$ contained within $\Omega$. By Dirichlet domain monotonicity, for every $x$,
$$ \sum_{A\in\mathcal{P}, A\subset G} N_A(x) \leq N_G(x). $$

Since $G$ tiles $\mathbb{R}^n$, by Polya's theorem, inequality \eqref{polya}, we have
$$ N_\Omega(x) \leq \frac{N_G(x)}{N(\sigma)} \leq C_n \frac{\sigma^n}{N(\sigma)} x^{n/2}. $$
This is true for every $\sigma > 0$. Letting $\sigma\to\infty$ and using the lemma proved below,
$$ N_\Omega(x) \leq C_n \frac{V}{\delta} x^{n/2} $$
for all $x$. Equivalently,
$$ \lambda_k(\Omega) \geq C_n \bigg( \delta \frac{k}{V} \bigg)^{2/n} $$
for all $k$ (where $C_n$ is a different constant depending only on dimension).

This completes the proof.
\end{proof}

\begin{lemma}
Let $\Omega$ be a bounded domain in $\mathbb{R}^n$. Let $\mathcal{P}$ be a maximum-density packing of $\Omega$. For $\sigma > 0$, let $N(\sigma)$ be the number of components of $\mathcal{P}$ entirely contained in $[-\sigma/2,\sigma/2]^n$. Then
$$ \lim_{\sigma\to\infty} \frac{\sigma^n}{N(\sigma)} = \frac{V}{\delta}. $$
\end{lemma}

\begin{proof}
For a more thorough summary of basic concepts in packing, we refer the reader to \cite{Toth1993}, section 2, and to \cite{Groemer1986}.

A packing $\mathcal{P}$ of $\mathbb{R}^n$ by $\Omega$ is a collection pairwise disjoint congruent copies of $\Omega$. If $G$ is a domain, define the inner density $d_{\operatorname{inn}}$ and outer density $d_{\operatorname{out}}$ with respect to $G$ to be
$$ d_{\operatorname{inn}}(\mathcal{P}|G) = \frac{1}{|G|}\sum_{A\in\mathcal{P},A\subset G}|A| $$
$$ d_{\operatorname{out}}(\mathcal{P}|G) = \frac{1}{|G|}\sum_{A\in\mathcal{P},A\subset G}|A| $$

We define the inner (resp. outer) densities of $\mathcal{P}$ with respect to the gauge $G$ as
$$ d_-(\mathcal{P},G,o) = \liminf_{\lambda\to\infty} d_{\operatorname{inn}}(\mathcal{P}|\lambda G) $$
$$ d_+(\mathcal{P},G,o) = \limsup_{\lambda\to\infty} d_{\operatorname{out}}(\mathcal{P}|\lambda G) $$
where $\lambda G$ is the image of $G$ under a homothety of scale $\lambda$ fixing $o$.

Call $(G,o)$ a gauge for the density. Then the packing density $\delta$ of $\Omega$ is defined to be the supremum of the outer densities $d_+(\mathcal{P},G,o)$ over all packings of $\Omega$ and all choices of gauge $(G,o)$.

According to a theorem of Groemer \cite{Groemer1986}, c.f. also section 2 of \cite{Toth1993}, for every compact domain $\Omega$ in $\mathbb{R}^n$, there exists a packing $\mathcal{P}$ by congruent copies of $\Omega$ such that
$$ d_+(\mathcal{P},G,o) = d_-(\mathcal{P},G,o) = \delta $$
for every gauge $(G,o)$.

Therefore we may choose a suitable gauge pair: $([-1/2,1/2]^n,0)$. Then
$$ \delta = \lim_{\lambda\to\infty} d_{\operatorname{inn}}(\mathcal{P},\lambda G) = \lim_{\sigma\to\infty} \frac{1}{\sigma^n V}\sum_{A\in\mathcal{P},A\subset\lambda\Omega}|A|. $$
In view of the fact that every $A$ is a copy of $\Omega$, we have
$$ \delta = \lim_{\sigma\to\infty} \sigma^{-n}\sum_{A\in\mathcal{P},A\subset\lambda\Omega} V = \lim_{\sigma\to\infty} \frac{V N(\sigma)}{\sigma^n}. $$

Therefore,
$$ \lim_{\sigma\to\infty} \frac{\sigma^n}{N(\sigma)} = \frac{V}{\delta}. $$
\end{proof}

\bibliographystyle{plain}


\section{Introduction}

Weakly Interacting Massive Particles (WIMPs) are one of the most well-motivated candidates for dark matter (DM) and have guided dedicated experimental searches in recent years.  In particular, direct-detection experiments have been optimized to detect neutral, weak-scale $\mathcal{O}(100)$~GeV particles that scatter off nuclei in targets.  Experiments such as LUX~\cite{Akerib:2013tjd}, XENON100~\cite{Aprile:2012nq}, and SuperCDMS~\cite{Agnese:2013jaa} are now so sensitive that they are beginning to probe the highly motivated Higgs-exchange regime.  In the next few years, the sensitivities will improve by several more orders of magnitude, closing in on the most-relevant regions of parameter space for WIMPs~\cite{Cheung:2012qy}.

However, as limits continue to tighten with no hints of signal detection, it is becoming increasingly worthwhile to consider loosening the assumptions on the WIMP paradigm; one possible direction is to consider weakly interacting sub-GeV DM.  Such light DM is motivated by several classes of models, including asymmetric~\cite{Kaplan:1991ah,Kaplan:2009ag,Lin:2011gj}, WIMPless~\cite{Feng:2008ya} and other scenarios~\cite{Boehm:2003ha,Boehm:2003hm,Borodatchenkova:2005ct,Fayet:2007ua,Hooper:2008im,Pospelov:2007mp}.

For DM candidates below the GeV mass range, elastic nuclear recoil energies fall below current detection thresholds. In this range, the possibility of inelastic processes such as DM-electron scattering leading to ionization become important, because the total energy available to scattering
is still appreciable.  Direct detection of electron-scattering events was studied in~\cite{Bernabei:2007gr,Dedes:2009bk, Kopp:2009et}.  Recently, there has been a renewed focus on the application of electron scattering experiments to the search for sub-GeV DM~\cite{Essig:2011nj,Graham:2012su}, and the first limits have been set using XENON10 data~\cite{Essig:2012yx}.

Two aspects of DM-electron scattering affect the phenomenology of such signals.  First, the scattering event is inelastic.  Inelastic scattering events have also been explored in the context of nuclear excitations~\cite{Ellis:1988nb}, as well as DM scattering to an excited state~\cite{TuckerSmith:2001hy}.  For an electron scattering event to occur, the minimum velocity that the DM must have to excite the electron depends on the bound-state energy of the electron, which of course depends on the detector target.  As a result, the experiments are not as sensitive to lower-velocity DM, which does not have the requisite minimum energy needed to excite the electron.

Secondly, the detection rate depends on an ionization form factor, which describes the likelihood that a given momentum transfer results in a particular electron recoil energy.  This form factor can be challenging to calculate, as it depends on the wave function of the scattered electron.  This ionization form factor is target-dependent and shapes the energy dependence of the event rate.

A significant challenge that will be faced by these experiments is the presence of background.  Identifying unique features of the DM signal that distinguish it from potential sources of background is therefore of the utmost importance.   An annually modulating signal, which arises due to the motion of the Earth around the Sun~\cite{Drukier:1986tm, Lee:2013wza}, is the most prominent example and may be a potentially crucial test of the DM origin of a potential signal.  Characteristics of the modulating signal, such as its phase and fractional amplitude, are expected to depend sensitively on physical parameters such as the assumed DM velocity distribution and the form of the coupling between the DM and target particles. Interactions between the target particle and its surroundings (for example, a target electron in a binding potential) can also have a significant effect.

As a result of the inelasticity of the signal, combined with the ionization form factor, it is non-trivial to extend our intuition of annual modulation from the nuclear-scattering case to the electron-scattering case.  To address this, this paper presents the first detailed study of the time-dependent rate in electron-scattering scenarios.  We complete a detailed study of the expected modulation fraction, higher-harmonic modes, and expected phase, for both atomic and semiconductor targets.  We find that electron-scattering events typically have large modulation fractions, and that the modulation phase may be affected by gravitational focusing  and local substructure.

We begin in Sec.~\ref{sec:rate} by introducing a semi-analytic approach to calculating the electron-scattering event rate for both atomic and semiconductor targets, building on previous methods.  Our approach for semiconductor targets should be more accurate than previous analytic approximations to the event rate, such as those in~\cite{Graham:2012su}, while at the same time being more tractable than the full numerical calculations presented in, for example,~\cite{Essig:2011nj}. Sec.~\ref{sec:mod} then applies these techniques to study the time-dependent characteristics of the signal. Sec.~\ref{sec: Earth effects} explores the effects of DM interactions inside the Earth.  We point out that in certain models, DM-nucleus scattering cross sections can be much larger than the DM-electron cross section.  Even though DM-nucleus scattering is not observable directly in the lab for these scenarios because of the low thresholds necessary, DM-nucleus scattering inside of the Earth can modify the lab-frame DM phase-space distribution for large enough cross sections.  We conclude in Sec.~\ref{sec:conc}.

\section{Calculating the Event Rate}
\label{sec:rate}

The kinematics of the inelastic process whereby DM ionizes an atomic electron is more complicated than that of DM-nuclear elastic scattering because the bound electron does not carry a fixed momentum.  As a result, the scattering process may take place with any momentum transfer ${\bf q}$ between the initial and final DM state.  However, when $q = |{\bf q}|$ deviates too far above the inverse Bohr radius, $a_0^{-1} \approx 3.7$ keV, the scattering rate receives a strong wave-function suppression, arising from the fact that it is unlikely for the atomic electron to be found with such a high momentum.

The relevant momentum transfers are significantly smaller than the nuclear masses we consider, which means that the nuclear recoil energy does not significantly contribute to energy conservation.  As a result, the energy conservation equation reads
\es{}{
({\bf p_\chi} + {\bf q})^2 = {\bf p_\chi}^2 - 2 m_\chi (E_{\text{er}} + E_b^i)\,,
}
where $E_{\text{er}}$ is the electron recoil energy, $E_b^i$ is the negative binding energy of the bound initial state (labeled by the index $i$), $m_\chi$ is the DM mass, and ${\bf p_\chi}$ is the initial DM momentum.  For a fixed $q$, the lowest DM speed $v_\text{min}$ that could induce an electron recoil $E_{\text{er}}$ is found by taking ${\bf q}$ to be antiparallel to ${\bf p_\chi}$:
\es{}{
v_\text{min} = {q \over 2 m_\chi} + {E_{\text{er}} + E_b^i \over q} \,.
}

The count rate for DM-induced electron ionization events is proportional to the average over the DM velocity distribution of the ionization cross section times the DM speed, $\langle \sigma^i_\text{ion} v \rangle$.
In Ref.~\cite{Essig:2011nj} (see also~\cite{Kopp:2009et}), it was shown that
\begin{equation}
\frac{d\langle\sigma_{\text{ion}}^{{i}}v\rangle}{d\ln E_{\text{er}}}=\frac{\bar{\sigma}_{e}}{8\mu^2_{e\chi}}\int dq\,q|f_{\text{ion}}^{{i}}(k',q)|^{2}|F_\text{DM}(q)|^{2}\eta(v_\text{min},t) \, ,
\label{eq:cross_section}
\end{equation}
where $\mu_{e\chi}$ is the reduced mass of the DM-electron system, and $\eta$  is the mean inverse speed.  The normalized cross section $\bar \sigma_e$ and the DM form factor $F_\text{DM}(q)$ may be calculated from the relevant matrix element for DM-free-electron scattering.  The function $|f_{\text{ion}}^{{i}}(k',q)|^{2}$ is the wave-function suppression factor to ionize an electron in the bound state labeled by $i$ to a final state with momentum $k'$, through a momentum transfer $q$.  We will discuss $|f_{\text{ion}}^{{i}}(k',q)|^{2}$ more later in this section.  However, for now, note that if the final state is a plane wave, then $k'=\sqrt{2m_{e}E_{\text{er}}}$, where $m_e$ is the mass of the electron.

The differential scattering rate involves a sum over the differential cross sections for all possible initial electron states, accounting for any degeneracies in the states:
\begin{equation}
\frac{dR}{d\ln{E_{\text{er}}}}= N_{T}\frac{\rho_{\chi}}{m_{\chi}} F(k') \sum_{i}\frac{d\langle\sigma_{ion}^{{i}}v\rangle}{d\ln E_{\text{er}}} \, ,
\label{eqRate}
\end{equation}
where $N_T$ is the number of target nuclei and \mbox{$\rho_\chi \approx 0.4$~GeV/cm$^{3}$} is the local DM density~\cite{Catena:2009mf,Pato:2010yq,Bovy:2012tw}.  As in the case of nuclear beta decay, the wave function of the scattered electron is distorted by the presence of the nearby atom, requiring that the rate be corrected by the Fermi factor, $F(k')$.
In the non-relativistic limit,
\begin{equation}
F(k')= \frac{2\pi\nu}{1-e^{-2\pi\nu}} \, ,
\end{equation}
where $\nu=Z_{\text{eff}}\,(\alpha m_{e}/k')$ and $\alpha$ is the fine-structure constant. $Z_\text{eff}$ is the effective charge that is felt by the scattered electron. Although this is expected to be somewhat larger than unity due to the imperfect shielding of the escaping electron by the remaining electrons, we set $Z_\text{eff} = 1$ throughout. As was shown in~\cite{Essig:2011nj}, this is expected to be a good approximation for outer-shell electrons, while inner-shell electrons may require somewhat higher values of $Z_\text{eff}$. Our choice of $Z_\text{eff} = 1$ is conservative, since larger values are expected to further enhance the rate.

The differential scattering rate depends on the convolution of the atomic physics factor, $|f_{\text{ion}}^{{i}}(k',q)|^{2}$, the particle physics term $\bar{\sigma}_e\,|F_\text{DM}(q)|^{2}$, and the astrophysical input $\eta(v_\text{min},t)$.  With this factorization in mind, we begin by summarizing the astrophysical input.
The mean inverse speed
\begin{equation}
{\eta(v_\text{min},t)\equiv \int_{v_\text{min}}^{\infty}\frac{ f_\oplus \left( \V{v}, t\right)}{v}\,d^3v,}
\end{equation}
depends on the Earth-frame velocity distribution of the DM, $f_\oplus \left( \V{v}, t\right)$, which acquires a time dependence as the Earth orbits the Sun.
In the Galactic frame, and asymptotically far away from the Sun's gravitational potential, we take the velocity distribution $f_\infty (\V{v})$ to be that of the Standard Halo Model (SHM):
\es{SHM}{
 f_{\infty} (\V{v}) = \left\{ \begin{array}{ll}
{1 \over N_{\text{esc}} } \left( {1 \over \pi v_0^2 } \right)^{3/2} e^{- \V{v}^2 / v_0^2 } \qquad &|\V{v}| < v_\text{esc} \\
0 \, \qquad &\text{else} \,,
\end{array}
\right.
}
where $N_{\text{esc}}$ is a normalization factor, and we take \mbox{$v_{0}\approx220$}~km/s~\cite{Kerr:1986hz} and the escape velocity \mbox{$v_\text{esc} \approx 550$}~km/s~\cite{Piffl:2013mla}.
\begin{figure}[tb]
   \centering
   \includegraphics[width=3.5in]{Figures/FormFactor.pdf}
   \caption{The integration range in the rate computation includes values of the momentum transfer $q$ for which $v_\text{min}$ dips below the Earth-frame escape velocity ($v_\text{min} \lesssim 770$ km/s, dashed green line).
At the same time, the ionization form factor $|f_{ion}|^{2}$ (scale on the right) strongly suppresses larger values of $q$. Form factors are illustrated for 100 MeV DM scattering off a $5p$ shell Xenon electron and a $3d$ shell Germanium electron, both with $E_{\text{er}} = 15$ eV.}
   \vspace{-0.15in}
   \label{fig:qPlot}
\end{figure}

To a first approximation, the velocity distribution at the Earth's location may be found simply by applying a Galilean transformation to $f_\infty (\V{v})$ to transform from the Galactic frame to the lab frame, so that
\es{eq:Liouville}{
 f_\oplus (\V{v},t) \approx f_{\infty} \left( \V{v}_\odot +  \V{v} +  \V{V}_\oplus(t) \right) \,,
}
where \mbox{$\V{v_\odot} = (11, 232, 7)$} km/s~\cite{Schoenrich:2009bx} is the velocity of the Sun in Galactic coordinates and $\V{V}_\oplus(t)$ is the time-dependent velocity of the Earth in the Solar frame.  Equation~(\ref{eq:Liouville}) is corrected by the fact that the trajectories of slow-moving DM are deflected in the Sun's gravitational potential~\cite{Alenazi:2006wu,Lee:2013wza}.  This phenomenon, known as gravitational focusing, has important implications for annual modulation, and will be discussed in more detail in Sec.~\ref{sec:mod}.

The particle physics input for the differential rate consists of the normalized cross section $\bar \sigma_e$ and the form factor $F_\text{DM}(q)$.  Following~\cite{Essig:2011nj}, we define
\es{sigma_bar_e}{
\bar \sigma_e \equiv {\mu_{e\chi}^2 \over 16 \pi m_\chi^2 m_e^2} \left. \overline{ | {\cal M}_{e\chi} (q) |^2} \right|_{q^2 = \alpha^2 m_e^2 } \,
}
in terms of the squared and spin-averaged momentum-space matrix element $\overline{| {\cal M}_{e\chi} (q) |^2}$, so that $\bar \sigma_e$ is equal to the non-relativistic scattering rate between DM and free electrons at momentum transfer $\alpha m_e$.  The form factor $F_\text{DM}(q)$ then captures the $q$-dependence of the matrix element:
\es{FDMq}{
|F_\text{DM}(q)|^2 = \overline{ | {\cal M}_{e\chi} (q) |} /  \overline{ | {\cal M}_{e\chi} (\alpha m_e) |}  \,.
}

For most of this paper, we will frame our conclusions model-independently in terms of $\bar \sigma_e$ and $F_\text{DM}(q)$.  However, it is useful to keep certain models in mind that give rise to specific form factors.  One example is the class of dark-photon mediated models (see, for example,~\cite{Essig:2011nj,Lin:2011gj}) in which the DM is charged under a hidden $U(1)$ gauge group that has a small mixing with the Standard Model photon.  There are two interesting limits in this model: when the dark-photon mass, $m_A$, is much greater than the momentum transfer, and when it is much less.

To be concrete, let's consider the case where the DM is a Dirac fermion $\chi$.
Then, in the case where $m_A$ is much greater than the momentum transfer $m_\chi v \sim 10^{-3} m_\chi$, the effective Lagrangian for the DM-Standard Model fermion interactions may be written as
\es{Leff1}{
{\cal L_\text{eff}} = {\epsilon g_D \over m_A^2} \sum_{i}  Q_i (\bar \chi \gamma_\mu \chi)(\bar \psi_i \gamma^\mu \psi_i) \,,
}
where $Q_i$ is the electromagnetic charge of the Standard Model fermion $\psi_i$, $g_D$ is the charge of $\chi$ under the dark photon, and $\epsilon$ is the small mixing parameter between the dark and visible photons.  In this case, $F_\text{DM}(q)  = 1$ and $\bar \sigma_e \approx 16 \pi \alpha \alpha_D \epsilon^2 \mu_{e\chi}^2 / m_A^4$, where $\alpha_D = g_D^2 / 4 \pi$.
The second interesting limit is when $m_A \ll m_\chi v$.  In this case,  $\bar \sigma_e \approx 16 \pi \alpha \alpha_D \epsilon^2 \mu_{e\chi}^2 /  (m_e \alpha)^4$ and $F_\text{DM}(q) = \alpha^{2} m_{e}^{2}/ q^2$.
Motivated by these models, we will consider DM form factors $F_\text{DM}(q)$ = 1 and $\alpha^{2} m_{e}^{2}/ q^2$ throughout the paper; these form factors are generic for heavy and light mediator models, respectively.

Now, we turn our attention to the ionization form factor $|f_{\text{ion}}^{{i}}(k',q)|^{2}$.  The form factors are calculated differently in atomic and semiconductor targets.  We begin by reviewing the simpler case of atomic targets, and then present an improved semi-analytic approximation for the case of DM-electron scattering in semiconductor targets.

\subsection{Atomic Target}
\label{Xenon}
\begin{figure}[t]
   \centering
   \includegraphics[width=3.3in]{Figures/Rate.pdf}
   \caption{Differential event rate for DM form factors $F_\text{DM}(q)=1$ (\textbf{top}) and $\alpha^2 m_e^2/q^2$ (\textbf{bottom}), shown for Germanium (red) and Xenon (blue) targets assuming a cross section \mbox{$\bar\sigma_{e}=10^{-37}$~cm$^{2}$} and the SHM.  The solid(dashed) lines correspond to example masses $m_{\chi}=10(100)$~MeV.  The bands in the Germanium curves come from varying the valence-band initial wave functions between $4s$ and $4p$ states. The Germanium lines include the effect of inner-shell $3d$ electrons.}
      \vspace{-0.15in}
   \label{fig:a0}
\end{figure}

The main challenge in calculating the DM-electron scattering rate is determining the ionization form factor, $|f_{\text{ion}}^{i}(k',q)|^{2}$, which depends on the bound-state wave function of the electron, as well as its final-state wave function.  It is most easily determined for the case of an atomic target, so we briefly review this calculation before proceeding to the more challenging case of a semiconductor target.

Assuming plane-wave final states for the scattered electron, so that $k' = \sqrt{2 m_e E_\text{er}}$, the ionization form factor for spherically symmetric full shells with quantum numbers $(n\,\ell)$ reduces to
\begin{equation}
|f_{\text{ion}}^{n\ell}(k',q)|^{2}=\frac{(2\ell+1)k'^{2}}{4\pi^3 q}\int dk\,k|\chi_{n\ell}(k)|^{2} \, ,
\label{atomicFF}
\end{equation}
where $\chi_{n\ell}$ is the radial part of the bound-state wave function and the integral runs over allowed values of $k$ between $|k'\pm q|$ (see, for example,~\cite{Essig:2011nj}).

As an example, Fig.~\ref{fig:qPlot} shows the ionization form factor for the $5p$ states in Xenon (dashed red) and the $3d$ core-electron states in Germanium (dotted red),
determined using the Roothaan-Hartree-Fock (RHF) ground-state wave functions and fixing $E_\text{er} = 5$~eV.  The radial RHF wave functions are described as a linear combination of Slater-type orbitals and take the form
\begin{equation}
R_{n\ell}(r) =\sum_{j}C_{j\ell n}r^{n_{j\ell}-1} e^{-\zeta_{j\ell} r} \, ,
\end{equation}
where the coefficients are tabulated in~\cite{Bunge1993113}. Notice that the form factors fall steeply with momentum recoil;
the ionization form factors strongly bias the scattering towards low-momentum recoil.  In addition, the form factors do not necessarily fall monotonically and thus shape the differential scattering rate.

Also in Fig.~\ref{fig:qPlot}, we show $v_\text{min}$ (solid orange) as a function of $q$ for an example with 100 MeV DM.  When $v_\text{min} \gtrsim 770$~km/s,
no DM is moving fast enough in the Earth's frame to induce ionization, to a good approximation;
for the Xenon example in Fig.~\ref{fig:qPlot},
the allowed momentum transfer is constrained to be $5~\text{keV}\lesssim q \lesssim 500~$keV.

Fig.~\ref{fig:a0} shows the differential scattering rate for form factors $F_{\text{DM}}=1$ (top panel) and $F_{\text{DM}}=\alpha^{2}m_{e}^{2}/q^{2}$ (bottom panel) for $m_{\chi}=10$ and 100 MeV assuming a Xenon target.  Only the three outermost orbitals ($5p, 5s$ and $4d$), with respective binding energies $\sim$12, 26 and 76~eV, were used to calculate the rate.  We have verified that the contributions from more tightly bound electrons are negligible.

\subsection{Semiconductor Target}

Next, we consider the case of DM scattering off electrons in a semiconductor target, exciting them above the band gap. Semiconductor materials provide an ideal target to study DM-electron scattering because their band structure allows for electron ionization energies of $\mathcal{O}$(1)~eV compared to noble gas targets with binding energies of $\mathcal{O}$(10)~eV.  For example, any interaction depositing energy above the band gap of $\sim$0.67 eV results in ionization of electron-hole pairs to the conduction band in Germanium.

Several current detector technologies have the potential to take advantage of semiconductor targets to achieve sensitivity to sub-GeV DM.  In these experiments, the scattering signal is  amplified by drifting ionized electrons to induce detectable phonons. The CDMSlite mode of operation of the SuperCDMS experiment~\cite{Agnese:2013jaa}, for example, relies on voltage-assisted amplification of the ionization energy deposited by particle interactions in order to achieve an ionization threshold of 170 eV, which makes it sensitive to sub-GeV DM. Reduction in background levels
planned for SuperCDMS would further increase the sensitivity~\cite{Sander:2012nia}. The potential to use fully-depleted CCDs to achieve thresholds of 40 eV has been demonstrated by DAMIC~\cite{Barreto:2011zu}.
Given the progress that has been made with the CDMSlite detectors, we will focus mainly on Germanium targets.

Calculating the ionization form factor for a semiconductor target carries with it particular challenges, as the electrons are described by Bloch wave functions in a periodic lattice.  The ionization form factor can be determined using special packages, such as \texttt{Quantum Espresso}~\cite{QE-2009}, as was done in~\cite{Essig:2011nj}.  This is computationally intensive, however, and so an analytic method is also useful for obtaining estimates of experimental sensitivity.  One approach, presented in~\cite{Graham:2012su}, is to model each lattice site in a Germanium crystal by a hydrogen atom with a variable binding energy and calculate the appropriately binding-energy-averaged scattering rate off the $1s$ electron.
The simplifying assumption of a hydrogenic wave function allows one to obtain analytic expressions for the scattering cross section.

It is important to understand how well the analytic results in~\cite{Graham:2012su} describe the real ionization form factor in Germanium.  Towards that end, we present an alternate semi-analytic approach here to calculating the ionization form factor, which is related to that in~\cite{Graham:2012su} but relies on using the RHF wave functions for the electrons in Germanium instead of hydrogenic wave functions.  We find significant differences between the results of the two approaches.

Another difference to keep in mind when considering semiconductor targets lies in the experimental method of detection of a signal. While for atomic targets the final state involves an electron-ion pair, for semiconductor targets it involves creation of electron-hole pairs.
These final-state charge carriers are drifted using an applied electric field, generating Luke-Neganov phonons.  The energy of the phonons is detected~\cite{Agnese:2013jaa}, giving a direct measure of the number of electron-hole pairs created by the DM scattering.

The number of electron-hole pairs is a function of the total energy $E_d$ deposited into the material by the scattering DM, which is simply related to the electron recoil energy $E_\text{er}$ and the binding energy $E_b$: \mbox{$E_{d}=E_\text{er}+E_{b}$}.   The average energy deposited in order to create an electron-hole pair for Germanium is $\sim$2.9~eV above the band gap.  Thus, we may define the effective number of electrons in the conduction band to be \mbox{$n_{e} = 1 +  (E_d - 0.67 \, \, \text{eV} ) / (2.9 \, \, \text{eV}) $}, taking into account that the initial scattering event promotes one electron from the valence band to the conduction band. We will present results both in terms of $E_d$ and $n_e$.

The electronic states in a semiconductor lattice are described by Bloch wave functions, $\Psi_{\mathbf k}(\mathbf r)$, which may be expressed using Wannier functions:
\begin{equation}
\Psi_{\mathbf k}(\mathbf r) = \sum_{N}e^{{i\mathbf k\cdot\mathbf R_{N}}}\phi(\mathbf r-\mathbf R_{N}) \, ,
\label{bloch_wannier}
\end{equation}
where $\phi(\mathbf r)$ is a Wannier function localized at the site $\mathbf R_{N}$, $\mathbf k$ are the wavevectors in the first Brillouin zone (BZ) consistent with the lattice periodicity, and $N$ is the number of lattice sites.  In the tight-binding approximation, the electrons at a given lattice site are assumed to have limited interactions with the neighboring atoms.  In this case, an atom at a given lattice site is effectively isolated, and the Wannier functions are simply the free atomic orbitals.  Therefore, the Bloch wave function for a given band is the sum over all lattice sites of the associated atomic orbital.

For our purposes, the expression for the Bloch wave function simplifies even further.  The DM-electron interaction is localized to a single lattice site so long as the momentum transfer is
\begin{equation}
q \gtrsim (\text{Ge lattice constant})^{-1} \sim~0.4 \text{ keV} \, .
\end{equation}
In this case, the sum over lattice sites in (\ref{bloch_wannier}) disappears and the Bloch wave function is simply the free atomic orbital at the scattering site.

For large enough momentum transfers, the wave function of the scattered electron can be approximated as a plane wave.  Therefore, the total scattering cross section is obtained by considering the transition of an electron from a localized initial-state atomic wave function--with a ${\bf k}$-dependent binding energy--to a final-state wave function with plane-wave solution, at some energy $E_\text{er}$ above the conduction band minimum.  The atomic scattering can be calculated using the same prescription as that described in Sec.~\ref{Xenon}, with the appropriate RHF wave functions for Germanium.

We are interested in the directionally-averaged rate, for which the variability of the initial bound-state energy $E_{b}$ with ${\bf k}$ may be captured by the valence-band density of states $\rho(E_{b})$ (see, for example,~\cite{Graham:2012su});  the total differential event rate is then obtained by integrating over all binding energies, weighted by the density of states:
\begin{equation}
\frac{dR}{d\ln{E_{\text{er}}}}\approx N_{T}\frac{\rho_{\chi}}{m_{\chi}} F(k') \int dE_b\,\rho(E_b)\frac{d\langle\sigma_{\text{ion}}v\rangle}{d\ln E_{\text{er}}} \, .
\label{eq:dR}
\end{equation}
The isotropic valence-band density of states (Fig.~\ref{fig:dosGe}) is computed using the \texttt{GPAW} package~\cite{Mortensen:2005ab}, a density-functional theory code based on the projector-augmented wave method.  For Germanium, the density of states is peaked at bound-state energies of $\sim$4, 8, and 12~eV. These peaks correspond to predominantly $p$-like (III and IV, red), an admixture of $s$- and $p$-like (II, green) and predominantly $s$-like (I, blue) states in the band structure~\cite{Chadi75tight-bindingcalculations}.

\begin{figure}[t]
   \centering
   \includegraphics[width=3.3in]{Figures/GeDoS.pdf}
   \caption{Band structure of Germanium (\textbf{left}) and the resulting density of states (\textbf{right}) used in the cross-section calculation. Shown is the valence band associated with predominantly $p$-like (III and IV, red), a combination of $s$- and $p$-like (II, green) and predominantly $s$-like (I, blue) states.  The $k$-vectors in the band diagram correspond to a chosen set of high-symmetry points in the first Brillouin zone, with $\Gamma$ being the BZ center. The reference level for the binding energy is taken to be the bottom of the minimum-energy conduction band.}
   \label{fig:dosGe}
      \vspace{-0.15in}
\end{figure}

It would seem that---depending on the binding energy---we should take a different combination of $s$- and $p$-like atomic wave functions when calculating the expression for $d\langle\sigma_{\text{ion}}v\rangle / d\ln E_{\text{er}}$ that enters into~\eqref{eq:dR}.  However, we find in practice that taking either the $s$- or $p$-like wave functions independently leads to very similar results for the scattering rate as can be seen from Fig.~\ref{fig:a0}.  As such, we will work with either pure $s$- or $p$-like wave functions for simplicity and estimate our error in this approximation by the difference between the results computed from the respective wave functions.

In addition to the valence electrons, we also consider scattering off the inner-core electrons in Germanium. In general, these electrons do not participate in the bonding process which determines the crystalline structure, with the inter-atomic bonds in Germanium composed of the outer $4s$ and $4p$ states~\cite{Ferry:2013sc}. The least tightly bound set of core electrons reside in the $3d$ shell and can then be accurately described using the appropriate atomic RHF wave function and a shifted binding energy of 30 eV~\cite{Cardona:1978ab}. Including these states can affect the rate, while corrections from subsequent shells are found to be negligible.

We emphasize here that there are many detailed properties of the semiconductor band structure that our calculation ignores.  For example, the energy bands vary along different crystal axes, leading to anisotropy in the density of states.  In addition, depending on the momentum transfer, its effective mass might vary from that of a free electron.  For example, near the minimum of a band, where the curvature is large, the effective mass will be less than $m_e$.  In~\cite{Graham:2012su}, the final-state effective mass was taken to be $m_{*}=0.56\,m_{e}$, corresponding to the density-of-states effective mass at the bottom of the conduction band, for all electron recoil energies $E_{\text{er}}$. However, if the energy of the electron is sufficiently above the band gap, this effective mass does not accurately describe the final-state dispersion relation.  We assume that the final-state electron is sufficiently above the band gap such that it is well-modeled by a free-electron plane-wave solution, with free-electron mass.  We expect our approximations to break down when the energy of the scattered electron is very close to the conduction band minimum.

Returning to Fig.~\ref{fig:a0}, we can compare the predicted differential scattering rate for a 10 and 100~MeV DM scattering off Germanium with that for Xenon.  For the 10~MeV case, the rate is larger in Germanium for all electron recoil energies; this is a direct effect of the lower binding energy in the semiconductor.  The benefit of the semiconductor target is still present, but less pronounced, for the case of the heavier 100~MeV candidate.  However, the improvement in Germanium's sensitivity over Xenon's is quite dramatic for the form-factor suppressed case in both mass examples.
\begin{figure}[t]
   \centering
   \includegraphics[width=3.0in]{Figures/Comparison.pdf}
   \caption{The binned event rates for $m_{\chi} = 100$ MeV DM obtained using the hydrogenic approach (dashed line) and the Germanium $4s/p$ wave functions (solid band). The latter also includes the effect of inner-shell $3d$ electrons. The differential rates are shown for the cases of a heavy mediator with $F_{\text{DM}}=1$ (\textbf{top}) and a light mediator with $F_{\text{DM}}=\alpha^{2}m_{e}^{2}/q^{2}$ (\textbf{bottom}) assuming $\bar\sigma_{e}=10^{-37}$ cm$^{2}$.  The energy-bin width is that of a single effective electron.}
   \label{fig:comparison}
      \vspace{-0.15in}
\end{figure}

Next, we compare the differential scattering rate calculated using (\ref{eq:dR}) against that obtained by treating each lattice site in a Germanium crystal as a free hydrogenic wave function~\cite{Graham:2012su}.
By making the assumption that the DM scatters off a $1s$ hydrogenic state ($\psi_{1s}^\text{H}$), it is possible to obtain an analytic solution for the scattering cross section.  From~\cite{Graham:2012su}, the full analytic expression for the $F_\text{DM} = 1$ case is, for example,
\begin{equation}
\frac{d\sigma}{dE_{\text{er}}}\approx\frac{8\,\bar\sigma_{e}}{3 \pi\alpha\mu_{e\chi}^{2}}\frac{k'F(k')}{v^{2}\,(1+a_{0}^{2}\,q_{\text{min}}^{2})^{3}} \, ,
\end{equation}
with \mbox{$q_{\text{min}} = m_{\chi}v-\sqrt{m_{\chi}^{2}v^{2}-2m_{\chi}(E_{\text{er}} + E_{b})}$} the minimum kinematically allowed momentum exchange. The associated binned scattering rate for a 100~MeV DM is shown in Fig.~\ref{fig:comparison} (dashed blue).
The solid blue band in Fig.~\ref{fig:comparison} shows the range obtained using initial-state $4s$ and $4p$ RHF wave functions for Germanium ($\psi_{4s/4p + 3d}^{\text{Ge}}$).   The $3d$ inner-shell electrons are also included in the Ge calculation and account for the increase in rate at $E_d \gtrsim 30$~eV.  There are clear differences between the two calculations.  Namely, the hydrogenic approach underestimates the number of high-energy scattering events.  Heuristically, this difference can be attributed to outer-shell electrons being on average further away from the nucleus and hence less tightly bound. Using the appropriate $4s/p$ wave functions, as well as the inner-shell electrons, enables us to better model the tail of the differential rate at these energies.

Fig.~\ref{fig:sensitivity} shows the implications of underestimating the high-energy scattering rate on the projected 95\% C.L. sensitivities.  For a low-threshold ($n_e = 1$) experiment, the effect is minimal, as would be expected from the fact that the differential rate for the $\psi_{1s}^\text{H}$ and $\psi_{4s/4p+3d}^{\text{Ge}}$ cases are roughly comparable at low deposited energies.  However, as the threshold energy of the experiment increases, one becomes more sensitive to large $E_\text{d}$, and the differences between the two approximations become much more apparent.  Note that Fig.~\ref{fig:sensitivity} assumes 1 kg$\cdot$year of exposure.
The shaded region in Fig.~\ref{fig:sensitivity} corresponds to the XENON10 excluded region from the analysis in~\cite{Essig:2012yx}, which was performed with $\sim$15 kg-days of exposure on a Xenon target.

\section{Annual Modulation}
\label{sec:mod}

This section explores the annual modulation of the DM-electron scattering rate.  The time-dependence enters the rate via the Earth-frame DM phase-space distribution $\rho_\chi  f_\oplus(\V{v},t)$, which is determined not only by the velocity of the Earth with respect to the Galactic frame, but also by the position of the Earth in the gravitational potential of the Sun. This latter phenomenon is referred to as gravitational focusing (GF) and is especially important for slower-moving DM particles that linger in the Sun's potential~\cite{Alenazi:2006wu,Lee:2013wza}.
\begin{figure}[tb]
   \centering
   \includegraphics[width=3.2in]{Figures/Sensitivity.pdf}
   \caption{Sensitivities expected at 95\% confidence level (corresponding to $\sim$3.6 expected events~\cite{Essig:2011nj})
   assuming 1 kg$\cdot$year of exposure at a Germanium low-threshold experiment.   These results were obtained using the Germanium RHF $4s/4p$+$3d$ wave functions (solid bands, this work) and the hydrogenic wave functions (dashed lines) for different detector threshold energies, assuming no background.  The thresholds are designated by the number of effective electrons, $n_e$ (1$e^-$, 5$e^-$, 15$e^-$).  The results are shown for the cases of a heavy mediator with $F_{\text{DM}}=1$ (\textbf{top}) and a light mediator with $F_{\text{DM}}=\alpha^{2}m_{e}^{2}/q^{2}$ (\textbf{bottom}).  The XENON10 excluded region at 90\% CL is shown in shaded blue~\cite{Essig:2012yx}.  The hydrogenic approach underestimates the sensitivity at high thresholds.}
   \label{fig:sensitivity}
      \vspace{-0.15in}
\end{figure}

To properly account for GF,~(\ref{eq:Liouville}) must be corrected to include the fact that the trajectories of slow-moving DM are deflected in the Sun's gravitational potential.  In this case, Liouville's theorem requires that
\begin{equation} \label{eq:Liouville1}
\rho_\chi  f_\oplus(\V{v},t) = \rho_\infty f_{\infty} \left( \V{v}_\odot + \V{v}_{\infty} \left[ \V{v} +  \V{V}_\oplus(t) \right]\right) \, ,
\end{equation}
where $\rho_\infty$ is the DM density asymptotically far away from the Sun's potential well.  In addition, $\V{v_\infty}[\V{v_{\text{S}}}]$ is the velocity that a particle must have at asymptotic infinity in order to have a Solar-frame velocity $\V{v}_\text{S}$ when it reaches the Earth; it is given by
\begin{equation}  \label{eq:LRL}
\V{v_\infty}[\V{v_{\text{S}}}] = {v_{\infty}^2\V{v}_\text{S} + v_{\infty}(G M_{\odot} / r_{\text{E}}) \V{\hat r}_{\text{E}} - v_{\infty} \V{v}_{\text{S}} (\V{v}_{\text{S}} \cdot \V{ \hat r}_{\text{E}}) \over v_{\infty}^2 + (G  M_{\odot} / r_{\text{E}})-v_{\infty} (\V{v}_{\text{S}} \cdot {\bf \hat r_{\text{E}}}) }\,,
\end{equation}
where $\V{r}_{\text{E}}$ is the position of the Earth in the Solar frame ($\V{\hat r}_{\text{E}}$ is the unit vector and $r_\text{E}$ is the distance between the Earth and the Sun).  Energy conservation requires that \mbox{$v_{\infty}^2 = v^2 - 2  G  M_\odot / r_{\text{E}}$}, with $G$ the gravitational constant and $M_\odot$ the mass of the Sun.

Gravitational focusing can have a significant effect on the phase of the modulation, as was shown in~\cite{Lee:2013wza} for the case of nuclear scattering, and so we include it in the following analysis, which explores the modulation amplitude and phase of a DM-electron scattering signal in some detail.

\begin{figure*}[tb]
\leavevmode
\begin{center}$
\begin{array}{cc}
   \includegraphics[width=.95\textwidth]{Figures/Harmonics.pdf}
 \end{array}$
\end{center}
\vspace{-.50cm}
   \caption{Magnitude of the mode coefficients $A_{1}$ and $A_{2}$ relative to the unmodulated rate $A_{0}$, assuming a Xenon target, for DM masses $m_{\chi}=10$ MeV (\textbf{left}) and 100 MeV (\textbf{right}).  A momentum-independent DM form factor $F_\text{DM}(q) = 1$ is also assumed.}
 	\vspace{-0.15in}

   \label{fig:harmonics}
\end{figure*}

\begin{figure*}[tb]
\leavevmode
\begin{center}$
\begin{array}{cc}
   \includegraphics[width=.97\textwidth]{Figures/Exposure.pdf}
 \end{array}$
\end{center}
\vspace{-.50cm}   \caption{Exposure $E$ needed to observe $A_{1}$ (annual modulation) and $A_{2}$ at 95\% significance relative to that needed to observe the unmodulated rate $A_{0}$, assuming a Xenon target and a momentum-independent DM form factor, for DM masses $m_{\chi}=10$~MeV (\textbf{left}) and 100~MeV (\textbf{right}) and various energy thresholds, labeled in terms of $n_e$.}
	\vspace{-0.15in}

   \label{fig:exposure}
\end{figure*}

\subsection{Modulation Amplitude}

The amplitude of the modulation of the DM-electron scattering rate can be quantified by decomposing the rate into Fourier modes as
\begin{equation}
 \frac{d R}{d E_{\text{er}}}  = A_0 + \sum_{n = 1}^\infty\big[ A_n \cos n\, \omega (t - t_n)\big]\,\,,
\end{equation}
where $\omega$ is the angular frequency of the Earth's orbit.  Furthermore, each $A_n$ is a function of the recoil energy specifying the amplitude of the $n$th mode, which has a maximum at time $t_n$.

The leading mode $A_1$ is conventionally referred to as the annual modulation.  In the case of DM-electron scattering, annual-modulation fractions of ${|A_1/A_0| \sim \mathcal{O}(10\%)}$ are expected and are slightly larger than the $\sim$2--5\% fractions expected in the case of DM-nuclear scattering~\cite{Essig:2011nj,Freese:2012xd}. This results from the strong enhancement of low-$q$ events by the ionization form factor; as demonstrated in Fig.~\ref{fig:qPlot}, scattering events in the energy range of interest are primarily induced by DM from the tail of the velocity distribution (\emph{i.e.}, at large $v_\text{min}$), where the corresponding unmodulated rate $A_0$ becomes relatively small.
The modes beyond annual modulation can also provide valuable information about the dark sector~\cite{Gelmini:2000dm,Freese:2003na,Savage:2006qr,Fornengo:2003fm,Green:2003yh,Alves:2010pt,Chang:2011eb,Lee:2013xxa}. \
However, because their amplitudes are generally suppressed as $|A_n / A_0| \sim (V_\oplus/v_\odot)^n$, detection of these modes typically requires large exposures.

To illustrate these points, we plot in Fig.~\ref{fig:harmonics} the energy dependence of $|A_{1}/A_0|$ and $|A_{2}/A_0|$ for a Xenon target.  We assume the SHM velocity distribution and a momentum-independent DM form factor, and we consider DM masses $m_{\chi}=10$ MeV and 100 MeV.  The amplitude ratios are relatively flat and featureless over most of the relevant energy range, but do increase at high recoil energies that probe the extreme tail of the velocity distribution.  Fig.~\ref{fig:exposure} shows the exposure $E$ (detector mass times measurement period) required to detect the first two modes at 95\% significance relative to that needed to observe the unmodulated rate $A_0$, for various values of the energy threshold, computed as in~\cite{Lee:2013xxa}.  The required exposure grows exponentially with $n$, making detection of the higher modes increasingly difficult.

\subsection{Modulation Phase}

\begin{figure*}[t]
\leavevmode
\begin{center}$
\begin{array}{cc}
   \includegraphics[width=.95\textwidth]{Figures/Phase.pdf}
 \end{array}$
\end{center}
\vspace{-.50cm}
\caption{Time of year when the event rate is maximized as a function of the total deposited energy in a Germanium target with DM mass of $m_{\chi}=100$ MeV for $F_\text{DM} = \alpha^2 m_e^2/q^2$ (left) and  $F_\text{DM} = 1$ (right).  The solid black curve assumes the SHM velocity distribution, while the blue curve also includes an additional dark-disk component with a relative density of $\rho_\text{DD}/\rho_\text{SHM} \sim 0.25$. For reference, the dashed line indicates the maximum predicted for the SHM when GF is neglected.}
\vspace{-0.15in}
\label{fig:phase}
\end{figure*}

We now turn our attention to the modulation phase, which gives the time at which the event rate is maximized as a function of recoil energy.  Using \eqref{eq:Liouville1} and \eqref{eq:LRL}, the energy-dependent phase can be calculated for a given Galactic-frame velocity distribution $f_{\infty} (\V{v})$.

We first consider the SHM velocity distribution and review the behavior of the phase of the corresponding DM-nuclear scattering rate.   In this case, the rate is maximized on $\sim$June 1 at high $v_\text{min} \gtrsim 200$ km/s (and hence, at the corresponding high recoil energies).  When GF is neglected, the date of maximal rate abruptly shifts to $\sim$December 1 for $v_\text{min} \lesssim 200$ km/s.  However, when GF is properly accounted for, the date of maximum instead shifts gradually with decreasing $v_\text{min}$, asymptoting to $\sim$December 25 at low $v_\text{min}$~\cite{Lee:2013wza}.  Thus, including GF leads to a shift of nearly a month in the predicted phase at low recoil energies.

The effect of GF is more subtle in the case of DM-electron scattering.  This is because GF more strongly affects slow-moving DM and hence becomes important at ${v_\text{min}\lesssim 200 \text{ km/s}}$; however, as illustrated in Fig.~\ref{fig:qPlot}, events with momentum transfer $q$ corresponding to such $v_\text{min}$ are suppressed by the ionization form factor relative to those corresponding to high $v_\text{min}$.  Despite this, the phase shift (with GF) is similar to the nuclear-scattering case for $F_\text{DM} = 1$, as shown in Fig.~\ref{fig:phase}.  For $F_\text{DM} = \alpha^2 m_e^2/q^2$, the phase shift from June 1 is less enhanced, for the same deposited energy, due to the additional momentum suppression.

If we relax the assumption of the SHM velocity distribution and consider velocity substructures that increase the proportion of slow-moving DM, the effect of GF on the modulation phase may be enhanced.  For example, simulations find that DM subhalos may be disrupted by the stellar disk of their host galaxy, subsequently merging to form a dark disk that corotates with the stellar disk~\cite{Read:2008fh, Read:2009iv, Purcell:2009yp,Bruch:2008rx}.  These simulations suggest that such a disk might exist in the Milky Way and contribute to the local DM density at the level of $\rho_\text{DD} / \rho_\text{SHM} \sim 0.5$--2, where $\rho_\text{SHM}$ is the non-rotating component of the local DM density.  Other studies suggest that observations may constrain $\rho_\text{DD} / \rho_\text{SHM}$ below this range, however there are significant systematic uncertainties in this result~\cite{Read:2014qva}.

We shall consider a fiducial case of $\rho_\text{DD} / \rho_\text{SHM} \sim 0.25$.  The dark disk can then be modeled by an additional truncated Maxwellian component, which is boosted in velocity space so that it corotates with the stellar disk and added to the SHM velocity distribution in the correct proportion.
Taking typical values observed in simulations, we assume a dispersion $v_{0}=70$~km/s and a corotation lag speed $v_\textrm{lag} = 50$~km/s.  The presence of such a corotating component increases the number of slow-moving DM particles in the Solar rest frame, which results in the somewhat larger phase shifts shown in Fig.~\ref{fig:phase}.  It is interesting to compare these results with those predicted for the phase shift in a DM-nucleus scattering experiment in the presence of a dark disk~\cite{Lee:2013wza}.  In that case, it was found that GF can lead to phase shifts at low $v_\text{min}$ of order a month or more.  The same is true here for the DM-electron scattering scenario.

In this section, we have discussed the modulation of the DM-electron scattering rate resulting from the Earth's motion and GF.  We have examined the energy dependence of the modulation amplitude and phase, highlighting the differences from the case of DM-nuclear scattering that arise due to the different scattering kinematics and the ionization form factor.  We next explore the potential implications of DM interactions inside the Earth on the DM phase-space distribution in the lab frame.

\section{DM interactions inside the Earth}
\label{sec: Earth effects}
   \begin{figure*}[htb]
	\leavevmode
	\begin{center}
	\includegraphics[width=.98\textwidth]{Figures/EarthEffects.pdf}

	\end{center}
	\vspace{-.50cm}
	\caption{(\textbf{Left}) The $\bar \sigma_e$-$m_\chi$ parameter space in the dark-photon mediated DM model with DM form factor $F_\text{DM}(q)=1$.  For cross sections above the black curve, we estimate that the DM will scatter at least once off of nuclei with a large angle while traversing the Earth's interior.  (\textbf{Right}) The same as the left panel, but taking DM form factor $F_\text{DM}(q) = (\alpha m_e / q)^2$.  The XENON10 excluded regions~\cite{Essig:2012yx} are shown in blue and the parameter space for the dark-photon mediated models~\cite{Essig:2012yx} is indicted by orange.  }
	\vspace{-0.15in}
	\label{Fig: Fq}
\end{figure*}

The DM phase-space distribution near the surface of the Earth can be distorted if the DM scatters with nuclei while traversing the Earth's interior.  In the lab frame, these distortions may even acquire time dependence as the Earth rotates, leading to a daily modulation.  A daily modulation induced by scattering in the Earth's interior has been discussed before in the context of nuclear recoils by~\cite{Cline:2012is,Collar:1993ss,Foot:2003iv,Hasenbalg:1997hs,Kouvaris:2014lpa,Mack:2007xj,Sigurdson:2004zp,Zaharijas:2004jv}.

A key feature that experiments searching for DM-induced ionization signals should keep in mind is that, even though the DM may be so light that the DM nuclear recoils are undetectable in the lab, the nuclear recoil cross section may be significantly larger than the ionization cross section.  This opens the possibility that the light DM may have nuclear recoils inside the Earth before being detected in the lab by electron ionization. If the nuclear recoil cross section is strong enough, this effect can significantly alter the lab-frame DM phase-space distribution.

The dark-photon model discussed in Sec.~\ref{sec:rate} yields DM-nucleon scattering at tree level, in addition to DM-electron scattering.  We now estimate the cross sections where Earth effects are expected to become important in this scenario, for the mass range $m_e < m_\chi < m_N$, where $m_N$ is the mass of the nucleus. Note, however, that this model merely serves as an illustration, and  the results also apply to other DM models with similar low-energy physics.

In the case where the dark-photon mass $m_A$ is much greater than the momentum transfer $m_\chi v \sim 10^{-3} m_\chi$, the differential cross section for nuclear recoils is
\es{dsigdEr}{
{d \sigma_{N \chi} \over d E_{\text{nr}}} \approx {8 \pi Z^2 \alpha \alpha_D \epsilon^2   m_N \over m_A^4 v^2} \,,
}
where $Z$ is the atomic number of the nucleus and $E_{\text{nr}}$ is the nuclear recoil energy. The atomic number enters instead of the mass number because the DM only interacts with the protons.

The recoil energy is related to the scattering angle $\theta$ in the center-of-mass of the DM-nucleus system by
\es{}{
E_{\text{nr}} = {2 \mu_{N \chi}^2 \over m_N} v^2 \sin^2(\theta/2) \,,
}
where $\mu_{N \chi}$  is the reduced mass of the DM-nucleus system.  The total cross section is obtained by integrating~\eqref{dsigdEr} over all allowed recoil energies.
  It follows that the total nuclear-scattering cross section obeys the relation
\es{echi}{
{\bar \sigma_{e} \over \sigma_{N\chi}} \approx \left({ \mu_{e\chi} \over Z \, \mu_{N \chi}} \right)^2 \,,
}
where ${\bar \sigma_{e}}$ is the electron-scattering cross section for the heavy mediator case.  Importantly, $\sigma_{N\chi}$ is much greater than $\bar \sigma_e$ for typical elements inside the Earth and for DM masses between the electron and nuclear mass scales.

Now consider a DM particle traveling through the Earth.  There is some probability that this DM particle will have a single scatter with nuclear matter.  Because the differential cross section is isotropic in the center-of-mass frame---and the center-of-mass frame almost coincides with the lab frame---there is a high chance that after the scattering event, the DM will recoil at a large angle.
To estimate the number of large-angle scattering events, we model the Earth simply as a uniform sphere comprised of approximately 32\% Fe, 30\% O, 15\% Si, and 14\% Mg, ignoring all other trace elements.
Given that the Earth has a density $\sim$$5.5$ g/cm$^3$ and radius $\sim$$6\times10^3$~km, the DM-electron scattering cross section must satisfy
\es{IE1}{
{\bar \sigma_e } \gtrsim & \,\, 4 \times10^{-39} \left({ 100 \, \, \text{MeV} \over m_\chi} \right)^2 \left( {\mu_{e\chi} \over 0.5 \, \, \text{MeV}} \right)^2 \, \, \text{cm}^2\, , \\
}
in the limit $m_\chi \ll m_N$, in order for there to be at least one scattering event as the DM traverses the Earth's interior.  This bound is illustrated in the left panel of Fig.~\ref{Fig: Fq}; Earth effects are important for $\bar \sigma_e$ above the black curve.  The XENON10 excluded region (shaded blue) and the allowed parameter space for the dark-photon mediated model (shaded orange) are also shown~\cite{Essig:2012yx}.  As Fig.~\ref{Fig: Fq} illustrates, there is a clear region of allowed parameter space where Earth effects should be important.\footnote{In specific scenarios, collider constraints from \emph{e.g.} monojet searches may also be constraining~\cite{Khachatryan:2014rra, Aad:2015zva}.}

Note that~\eqref{IE1} is specific to a dark-photon model where the DM has a tree-level coupling to quarks.  If the only Standard Model fermion $\psi$ appearing in~\eqref{Leff1} is the electron, then the nucleon-DM coupling arises at one loop~\cite{Kopp:2009et}.  In this case, the right-hand side of~\eqref{IE1} should be divided by a factor $\mathcal{O}(\alpha^{2})$, which significantly reduces the allowed parameter space for Earth effects in Fig.~\ref{Fig: Fq}.

Next, we consider the limit where the momentum transfer $q$ is much greater than the dark-photon mass $m_A$.  In this case, the differential cross section for DM-nucleon scattering is
\es{diff1}{
{d \sigma_{N \chi}  \over d E_{\text{nr}}} \approx {2 \pi Z^2 \alpha \alpha_D \epsilon^2 \over v^2 m_N E_{\text{nr}}^2} \, .
}
Importantly, the total cross section is IR divergent, which results from the fact that there is a new long-range force.  However, since the electrons are also charged under the dark photon, we expect the force to be screened over a distance $\sim$$a_0$, just like in ordinary electromagnetism.  We take this into account by requiring\footnote{In these DM models, the DM-electron ionization scattering events are also screened for momentum transfers \mbox{$q \lesssim a_0^{-1} \approx 4$}~keV.  This may affect the detection prospects for these models relative to what is shown in Fig.~\ref{fig:sensitivity}. }
\es{Enrscreen}{
E_{\text{nr}} > E_{\text{nr}}^\text{screen} \approx {1  \over 2 m_N a_0^2} \,.
}
  Additionally, we require the DM to be deflected at an angle greater than $\theta_\text{min} \sim 1$, which amounts to imposing a lower bound
  \es{Enrhard}{
  E_{\text{nr}} > E_{\text{nr}}^\text{hard} = 2 {\mu_{N \chi}^2 \over m_N} v^2 \sin^2(\theta_\text{min} / 2)
  }
   on the recoil energy.  The hard-scattering lower energy bound is greater than the screening bound for \mbox{$m_\chi \gtrsim 1 / (2 \sin(\theta_\text{min} / 2) a_0 v)$}.

   In the large $m_A$ scenario, there is no screening bound $E_\text{nr}^\text{screen}$ because the dark force is short range.  Similarly, in that scenario the bound $E_\text{nr}^\text{hard}$ is not necessary because the differential cross section~\eqref{dsigdEr} is independent of $E_\text{nr}$.  However, in the light mediator case the differential cross section~\eqref{diff1} rises steeply for low-angle scatters.  These low-angle scattering events do not significantly modify the trajectories of the DM, and so we exclude them by imposing the bound $E_\text{nr}^\text{hard}$ (see~\cite{Agashe:2014yua} for related discussions).

When $E_{\text{nr}}^\text{hard} > E_{\text{nr}}^\text{screen}$, we find that
\es{}{
{\bar \sigma_e \over \sigma_{N \chi}} = \tan^2\left(\theta_\text{min} / 2 \right) {16 v^4 \mu_{N\chi}^2 \, \mu_{e\chi}^2 \over Z^2 \,(m_e \alpha)^4} \, ,
}
while when $E_{\text{nr}}^\text{hard} < E_{\text{nr}}^\text{screen}$,
\es{}{
{\bar \sigma_e \over \sigma_{N \chi}} = {16 \, v^4 \mu_{N\chi}^2 \, \mu_{e\chi}^2 \over Z^2 \,(m_e \alpha)^4}  \left( {1 \over 4 \mu_{N \chi}^2 a_0^2 v^2 - 1} \right)\,.
}
Requiring that there be at least one hard scatter (\mbox{$\theta_\text{min} \sim \pi/4$}) while the DM traverses the full length of the Earth's diameter  then gives the lower bound (black line) shown in the right panel of Fig.~\ref{Fig: Fq}.

The estimates in this section are meant to roughly approximate the DM-electron scattering cross sections for which nuclear scattering in the Earth's interior becomes relevant.  As seen in Fig.~\ref{Fig: Fq}, these effects may be important over a wide range of the allowed parameter space.  The results presented here, however, require that the DM have tree-level couplings with the quarks; if this interaction is loop-suppressed, then the prospect of observing Earth effects is less optimistic.  Additionally, the current XENON10 exclusion limits are contained within the region where these Earth effects are relevant; this may mean that it is necessary to reinterpret these experimental constraints in light of the modified phase-space distribution for these particular sub-GeV DM models.  Along those lines, we have not attempted to quantify the observable consequences of the DM interactions inside the Earth.  We suspect, however, that the phenomenology could be manifested by a daily modulation of the rate.  We leave such investigations to future work.

\section{Conclusion}

For sub-GeV DM, the main avenue for discovery in a direct-detection experiment is the detection of DM scattering off target electrons.  The properties of such a signal are affected by the inelastic nature of the electron excitation, as well as the electron ionization form factor.  We studied signals in both atomic (\emph{e.g.}, Xenon) and semiconductor (\emph{e.g.}, Germanium) targets.  In particular, we presented a new semi-analytic approach to calculate the scattering rates for semiconductor targets; this approach makes it tractable to estimate an experiment's sensitivity without relying on numerical packages to obtain detailed modeling of the semiconductor's properties.  We argue that most of the detailed band-structure physics should have little effect on the DM-electron scattering predictions.

The annual-modulation fraction for DM-electron scattering is found to be ${\sim\mathcal{O}(10\%)}$ over a large range of candidate masses, which is somewhat higher than that expected from nuclear recoils. Observation of higher harmonic modes ($n \gtrsim 2$) is challenging, typically requiring $\sim10^{2}$--$10^{4}$ times more exposure to observe after the annual modulation. The phase of annual modulation is shown to be affected by gravitational focusing due to the Sun, similar to the case of nuclear scattering.

We showed that DM-nuclear interactions inside the Earth can cause sub-GeV DM candidates to scatter before reaching the detector, possibly leading to a daily modulation and directional dependence of the count rate over a significant portion of the motivated parameter space. Further study is required to accurately quantify this effect, which we leave to future work.\\

\section*{Acknowledgements}

\noindent We are especially grateful for the many enlightening discussions we have shared with Rouven Essig and Tien-Tien Yu.  We also thank Ilya Belopolski, David Huse, Christopher Malbon, Jeremy Mardon, Jesse Thaler, Tomer Volansky, and David Wallace for helpful discussions.  B.R.S was supported by a Pappalardo Fellowship in Physics at MIT.  ML is supported by the U.S. Department of Energy under grant Contract Number DE-SC0007968.

\label{sec:conc}

\newpage
\onecolumngrid
\vspace{0.3in}
\twocolumngrid
\def\bibsection{}
\bibliographystyle{utphys}


\section{Introduction}
We call a length metric space, $(Y,d,p)$, a Ricci limit space if it is a Gromov-Hausdorff limit of a sequence of pointed complete Riemannian n-manifolds, $(M_i^n,p_i)$, with Ricci curvature $\mathrm{Ric}_{M_i^n}\geq -(n-1)$. In this note, we always assume that $Y$ is not a point. In papers \cite{B1,B2,B3,B4}, many important results about such limit spaces have been proven by Cheeger and Colding. By Gromov's compactness theorem, for any point $y\in Y$, given any sequence $\epsilon_j \rightarrow 0$, passing to a subsequence, $(Y,\epsilon_j^{-1}d, y)$ converges to a length space, $(T_y, 0_y)$, denoted by $(Y,\epsilon_j^{-1}d, y)\rightarrow (T_y, 0_y)$, and $(T_y, 0_y)$ is called a tangent cone at $y$ associated to $\epsilon_j\rightarrow 0.$ A tangent cone at $y$ may depend on a choice of $\epsilon_j\rightarrow 0$, even in a non-collapsing Ricci limit space $Y$ i.e., the volume of unit balls at $p_i$, $\mathrm{vol}(B_1(p_i))\geq v>0$ for all large $i$. Examples of non-collapsing Ricci limit spaces were constructed in \cite{B2},\cite{B6} whose tangent cones at a point in $Y$ are not isometric or even not homeomorphic. In a collapsing sequence, $\mathrm{vol}(B_1(p_i))\rightarrow 0$, tangent cones at a point in $Y$ may have different Hausdorff dimensions (\cite{B2}). If all tangent cones at a point, $y\in Y$, are isometric to $\mathbb{R}^k$, for some integer $k$, $y$ is called a k-regular point. The set of k-regular points in $Y$ is denoted by $\mathcal{R}_k$. A point in $Y$ is called singular if it is not a regular point. The set of singular points in $Y$ is denoted by $\mathcal{S}$. By Cheeger-Colding (\cite{B2}) and Colding-Naber (\cite{B5}), there exists an integer $k$ such that $\mathcal{R}_k$ is dense in $Y$ and has a
full renormalized measure (Lemma 2.2(2)). However, it is still open whether $\mathcal{R}_l\neq \emptyset,$ for some $l\neq k.$

The main result in this note asserts that if $\mathcal{R}_1 \neq \emptyset$, then $\bigcup_{l=2}^{n-1}\mathcal{R}_l$ is empty.
Furthermore, $Y$ is a one dimensional topological manifold with or without boundary.

\begin{Thm}
Let $Y$ be a Ricci limit space. If $\mathcal{R}_1 \neq \emptyset$ , then $Y$ is a one dimensional topological manifold
with or without boundary.
\end{Thm}

\begin{Cor} Let $Y$ be a Ricci limit space of Hausdorff-dimension $k\ge 2$. Then $\mathcal{R}_1=\emptyset$.
\end{Cor}
Let $\nu$ denote a renormalized measure (see section 1 of \cite{B2} for the definition of renormalized measure).
A geodesic in $Y$ is called a limit geodesic if it is the limit of geodesics in $M_i$.
\begin{Cor} Let $Y$ be a Ricci limit space such that the renomalized measure $\nu(Y\setminus \mathcal{R}_2)=\emptyset$. Then, for any $x, y \in Y$,
and any limit geodesic $\gamma$ connecting $x$ and $y$, interior points of $\gamma$ are all in $\mathcal{R}_2$ or
all in the singular set $\mathcal{S}$.
\end{Cor}

 Theorem 1.1 improves Theorem 1.1 in \cite{Honda}, where Honda proved that if $1\leq \mathrm{dim}_H(Y)<2$,
  then $Y$ is a one dimensional topological manifold. Note that $1\leq \mathrm{dim}_H(Y)<2$ is equivalent to that $\bigcup_{l=2}^{n-1}\mathcal{R}_l = \emptyset$ (\cite{B3}, \cite{Honda}), both imply that $\mathcal{R}_1 \neq \emptyset$. We noticed that in \cite{Honda}, Honda claimed (page 3 in \cite{Honda}) that in a paper under preparation, he can proved Theorem 1.1.

 Our proof uses the following splitting theorem of Cheeger-Colding (\cite{B1}),
 \begin{Thm}(\cite{B1})
 Let $(Y,p)$ be a Gromov-Hausdorff limit of a sequence, $\{(M_i^n, p_i)\},$ with $\mathrm{Ric}_{M_i^n}\geq -(n-1)\delta_i, \delta_i\rightarrow 0.$ If $Y$ contains a line, then $Y$ splits isometrically with a $\Bbb R^1$-facor i.e., $Y=\mathbb{R}^1\times X$ for some length metric space $X$.
 \end{Thm}

Using Theorem 1.2, we show that for $x\in \mathcal{R}_1$, there exist $\epsilon = \epsilon(x)>0, r(\epsilon)>0,$ such that if $r\leq r(\epsilon)$, $y^{\pm}\in A_{r-2\epsilon r, r+2\epsilon r}(x)=B_{r+2\epsilon r}(x)\setminus B_{r-2\epsilon r}(x),$  with $d(y^{+},x)+d(y^{-},x)-d(y^{-},y^{+})\leq 7\epsilon r,$ then $x$ lies in every geodesic connecting $y^{-}$ and $y^{+}$ (Lemma 2.1). Together with some applications of the sharp H\"older continuity on tangent cones obtained by Colding-Naber (Lemma 2.2)(\cite{B5}) and a similar argument to \cite{Honda}, theorem 1.1 can be proved.
\begin{Rem}
In \cite{B9}, Kitabeppu and Lakzian proved a similar result of Theorem 1.1 for a more general metric spaces (including Ricci limit spaces) via a complicated argument, because
the continuity of tangent cones along the interior of geodesic may not hold.
Our proof is simple but may not extend to their case. Nevertheless, Lemma 2.1 still holds in the setting of \cite{B9}.
\end{Rem}

The author would like to thank professor Shouhei Honda for his interests in this paper and for his useful suggestion that leads to Corollary 1.1.2. The author is grateful to Professors S. Honda and Tapio Rajala for bring her attention to \cite{B9}.

\section{Local structure of 1-regular points}
Let $x\in \mathcal{R}_1$. By definition, for every $\epsilon > 0$, there is $r(\epsilon)>0$ such that
for any $r\leq r(\epsilon)$,
$$d_{GH}(B_r(x),B_r^1(0))\leq \epsilon r,$$
where $B_r^1(0)\subset \mathbb{R}^1.$ Let $e_1$ be the standard basis in $\mathbb{R}^1$. Hence for every $\epsilon > 0$,
$r\leq r(\epsilon)$, there are $x_r^{+}, x_r^{-}\in B_r(x)$ such that $d(x_r^{\pm},\pm re_1)\leq \epsilon r$.
Consequently,
$$|d(x_r^{\pm},x)-r|\leq \epsilon r,$$
$$d(x_r^{+},x)+d(x_r^{-},x)-d(x_r^{+},x_r^{-})\leq 3\epsilon r.$$
\begin{Lem} Let $(Y ,d)$, $x\in \mathcal{R}_1, x_r^{+}, x_r^{-}$ be as in the above. Then there exists $\epsilon_0 > 0$
such that for any $\epsilon \leq \epsilon_0$, $r\leq r(\epsilon)$,  and any $$(y^{+},y^{-})\in
B_{\epsilon r}(x_r^{+})\times B_{\epsilon r}(x_r^{-}),$$ any geodesic $\gamma$ from $y_1$ to $y_2$ passing through $x$.
\end{Lem}
\begin{proof}
Arguing by contradiction, assume that there are $ \epsilon_i \rightarrow 0, r_i \rightarrow 0$, $(y_i^{+}, y_i^{-})\in B_{\epsilon_i r_i}(x_{r_i}^{+})\times B_{\epsilon_i r_i}(x_{r_i}^{-})$, and a geodesic, $\gamma_i$, from $y_i^{+}$ to $y_i^{-}$ such that $x\notin \gamma_i$.

By the triangle inequality,
$$|d(y_i^{+}, x)-r_i|\leq 2\epsilon_i r_i, |d(y_i^{-}, x)-r_i|\leq 2\epsilon_i r_i,$$
$$d(y_i^{+}, x)+d(y_i^{-},x) - d(y_i^{+}, y_i^{-})\leq 7\epsilon_i r_i.$$
Passing to a subsequence, we may assume that
$$(Y, r_i^{-1}d,x)\rightarrow (\mathbb{R},0),$$
and $y_i^{\pm}\rightarrow \pm 1.$

By assumption, $s_i = d(x,\gamma_i)> 0$. From the above there is $\gamma_i(t_i)$ such that
$d(x,\gamma_i)=d(x,\gamma_i(t_i))$ and
$$s_i\leq \min\{d(x, y_i^{+}),d(x, y_i^{-})\}.$$
Thus $s_i\rightarrow 0,$ as $i\rightarrow \infty$.
Without loss of generality, we can assume that
$$(Y, s_i^{-1}d, x)\rightarrow (T_x, 0_x).$$
Observe that $s_i^{-1}\gamma_i$ converges to a line in $T_x$ which has distance $1$ from $0_x$.
By Theorem 1.2 (\cite{B1}), $T_x$ is isometric to $\mathbb{R}\times W$, and because
$0_x$ is not on the line, $W$ is not a point; a contradiction to $x\in \mathcal{R}_1$.
\end{proof}
To conclude Theorem $1.1$ from Lemma 2.1, we will also need the following results of Colding-Naber (\cite{B5}).

\begin{Lem}(\cite{B5})
 If $Y$ is a Ricci limit space, then
 \begin{enumerate}
\item[(1)]
for $\nu \times \nu$ almost every pair $(a_1, a_2)\in A_1 \times A_2$, where $A_1$
and $A_2$ are subset of $Y$ that contained in a bounded ball, there exists a limit
geodesic from $a_1$ to $a_2$ whose interior lies in some $\mathcal{R}_l$,
$l$ is an integer;
  \item[(2)]
  There is an integer $k$, such that $\nu(Y\setminus \mathcal{R}_k)=0$.

   \end{enumerate}
\end{Lem}
  Let $\mathcal{WE}_k$ be the set of points, $y\in Y$, such that there exists a tangent cone at $y$ which is isometric to $ \mathbb{R}^k \times W$ and $W$ is a length metric space. Let $\underline{\mathcal{WE}}_k = \mathcal{WE}_k \setminus \mathcal{R}_k$.
\begin{Rem}$ $
\begin{enumerate}
\item[(1)]
 As an essentially direct consequence of Cheeger-Colding's work in \cite{B2,B4}, Honda (\cite{Honda}) pointed out that if $x \in \underline{\mathcal{WE}}_k$, then for all $r>0$, $\nu(B_r(x)\cap (\bigcup_{l\geq k+1}\mathcal{R}_l))>0$. Consequently, if $\nu(\cup_{l\geq k+1}\mathcal{R}_l)=0$, then $\underline{\mathcal{WE}}_k = \emptyset,$ especially, $\mathcal{R}_l = \emptyset,$ for all $l\geq k+1.$
 \item[(2)]
 By Lemma 2.2 (2), almost every pair of points in $Y$ can be connected by a limit geodesic whose
 interior is in $\mathcal{R}_k.$
\end{enumerate}
\end{Rem}
\begin{proof}[Proof of theorem 1.1]
 We complete the proof in two steps.

 First, we claim that for $x\in \mathcal{R}_1$, there exists $\tau=\tau(x)>0,$ such that $B_{\tau}(x)$ is isometric to $(-\tau,\tau).$

 Let $x\in \mathcal{R}_1, x_r^{+}, x_r^{-}, \epsilon_0, r(\epsilon)$ be chosen as in Lemma 2.1. Then for any $$(y^{+},y^{-})\in B_{\epsilon r}(x_r^{+})\times B_{\epsilon r}(x_r^{-}),$$ where $r\leq r(\epsilon)$ and $\epsilon\leq \epsilon_0\leq \frac{1}{4},$ any geodesic, $\gamma$, from $y^{+}$ to $y^{-}$, we have that $x\in \gamma$. By Lemma 2.2(1), we can assume $(y^{+},y^{-})\in B_{\epsilon r}(x_r^{+})\times B_{\epsilon r}(x_r^{-})$ and a limit geodesic $\gamma$
 from $y^{+}$ to $y^{-}$ whose interior is contained in $\mathcal{R}_l$ for some integer $l$. Since $x\in \gamma\cap \mathcal{R}_1$,
 we have that $l=1$ i.e. the interior of $\gamma$ is contained in $\mathcal{R}_1$.

 Next, we shall show that for $\tau\leq \frac{r(\epsilon_0)}{4}$, $B_\tau(x)$ is isometric to $\gamma|_{(-\tau,\tau)}$, where $\gamma$ is
 in above with a reparametrization such that $x=\gamma(0)$. If $B_\tau(x)$ is not isometric to $\gamma|_{(-\tau,\tau)}$, then there is
 a point $y \in B_{\tau}(x)\setminus \gamma$. We may assume $z\in \gamma$ such that $d(y,z)= d(y,\gamma)$. Observe that
 for every small $\delta>0$, there are $z^{\pm}\in \gamma$ such that $d(z^-,z^+)=d(z^-,z)+d(z,z^+)$ and $d(z^\pm,z)=\delta$,
 and $w$ in a geodesic from $z$ to $y$ such that $d(w,z)=\delta$ (note that $d(z^{\pm},w)\geq \delta$).
 Consequently, the tangent cone at $z$ associate to $\delta\to 0$ is isometric to $\mathbb{R}\times W$ and
 $W$ is not a point; a contradiction to $z\in \mathcal{R}_1$ (compare to the proof of Theorem 4.3, \cite{Honda}).

 A by-product of the above is that geodesics in $Y$ do not branch. By Lemma 2.2(2) and Remark 2.1(1),
 $\underline{\mathcal{WE}}_1 =\emptyset$. Consequently, any interior point $y$ in a geodesic has a
 a neighborhood isometric to $(-\tau,\tau)$ (no branching at $y$). This allow us to uniquely extend
 any geodesic in $Y$ to a maximal geodesic.

 Secondly, we will show that $Y\setminus \gamma = \emptyset$ i.e., $Y$ is a one dimensional topological manifold.
 Arguing by contradiction, assume $y\in Y\setminus \gamma$ and $z\in \gamma$ such that $d(y,z)=d(y,\gamma)$.
 If $z$ is an interior point of $\gamma$ and $\sigma$ is a geodesic from $z$ to $y$, then
 $z\in \mathcal{R}_1$ and there exists $\delta>0$, such that $B_{\delta}(z)$ is isometric to $(-\delta, \delta)$. Thus $\sigma(\delta)\in \gamma$ and $d(y,\sigma(\delta))\leq d(y,z)-\delta$, a contradiction to that $d(y,z)=d(y,\gamma)$. If $z$ is one end of $\gamma$, then
 a similar discussion yields that for any $w$ lies in the interior of $\gamma$, a geodesic from $w$ to $y$ is
 the union of the piece of $\gamma$ from $w$ to $z$ and a geodesic from $z$ to $y$ i.e., $\gamma$ extends through $z$,
 a contradiction.
\end{proof}

\begin{proof}[Proof of Corollary 1.1.2]
Let $x, y\in Y$, and $\gamma$ be a limit geodesic connecting $x$ and $y$. Let $p$ be an interior point of $\gamma$, then $p\in \mathcal{WE}_1$.
Since $\nu(Y\setminus \mathcal{R}_2)= 0$, by Remark 2.2(1) and Theorem 1.1, for any $i\neq 2$, $\mathcal{R}_i=\emptyset$. If $p \in \mathcal{S}$, there is a tangent cone, $T_p = \mathbb{R}^1\times W$, at $p$ satisfies that $\mathrm{dim}_H(W)\geq 1$ (otherwise $p \in \mathcal{R}_1$ \cite{Honda}). Thus by Proposition 3.37 in \cite{B8} and Theorem 1.1, for any renormalized measure $\nu_{\infty}$ in $T_p$, $\nu_{\infty}(T_p\setminus \mathcal{R}_2)= 0$. By splitting theorem (Theorem 1.2), using an argument as the proof of Theorem 1.1 (see also Theorem 4.3 \cite{Honda} ), we have that $T_p$ is isometric to the half plan or $\mathbb{R}^2$. Then by the sharp H\"older continuity on tangent cones obtained by Colding-Naber (\cite{B5}), we easily conclude the proof.

\end{proof}


\section{Introduction}

The HERA collider, which operated at DESY from 1992 to 2007, has
produced many important physics results, first of all perhaps the most
precise determinations to date of the quark and gluon densities in the
proton from single experiments (H1, ZEUS) \cite{Adloff:2000qk,%
Chekanov:2001qu} and their combined data sets \cite{Aaron:2009aa}. These
data, taken in deep-inelastic electron-proton scattering, are
complemented by a wealth of data from photoproduction at low virtuality
$Q^2$ of the exchanged photon, in particular on jet production
\cite{Klasen:1994bj}, giving access
also to the distributions of partons in the photon \cite{Albino:2002ck}
and to measurements
of the strong coupling constant \cite{Klasen:2014}.

Using the data set of the HERA-II phase of the HERA collider
from 2003-2007 with an integrated
luminosity of $\unit{351}{pb^{-1}}$, the H1 collaboration have recently published
final measurements of inclusive jet, dijet and three-jet production in deep-inelastic
scattering (DIS)
\cite{Andreev:2014wwa,Britzger:2013}
and used them to determine the strong coupling constant (at the mass $M_Z$ of
the $Z$-boson) to be
\begin{equation}
 \alpha_s(M_Z)=0.1185\pm0.0016 ({\rm exp.})\pm0.0040\, ({\rm th.}),
 \label{eq:1}
\end{equation}
taking into account absolute double-differential inclusive jet, dijet and three-jet cross section
data as functions of $Q^2$ and the jet transverse momentum $p_T$. A more precise value was
obtained from normalized jet cross sections, yielding
\begin{equation}
 \alpha_s(M_Z)=0.1165\pm0.0008  ({\rm exp.})\pm0.0038\, ({\rm th.})\, .
 \label{eq:2}
\end{equation}
Unsatisfactorily, only the value obtained by unrenormalized results is in agreement with the
current world average
of $\alpha_s(M_Z)=0.1185\pm0.0006$ \cite{Agashe:2014kda}.
The latter uses only observables that
are known to next-to-next-to-leading order (NNLO) of perturbative QCD, while
the analysis of the H1 collaboration was done in next-to-leading order (NLO) accuracy.
The lack of knowledge of higher-order contributions becomes manifest in a bigger
theoretical uncertainty due to scale variation in Eq.\ (\ref{eq:1}).
The absolute double-differential cross section measurement of inclusive jets
alone led H1 to a value of the strong coupling constant of
\begin{equation}
 \alpha_s(M_Z)=0.1174\pm0.0022 ({\rm exp.})\pm0.0050\, ({\rm th.}).
 \label{eq:3}
\end{equation}

In this Letter, we compute the inclusive jet production DIS cross
section for the first time including NNLO contributions, obtained
from a unified threshold resummation formalism \cite{Kidonakis:2003tx},
and extract an approximate NNLO (aNNLO) value for the strong coupling constant.
Our calculations are based
on previous work on inclusive jet production in deep-inelastic scattering up to NLO
\cite{Potter:1997}. They reduce, as we will see,
the theoretical prediction and increase the central fit value of the strong
coupling constant, improving
the description of the final HERA data in particular at high photon virtuality $Q^2$.

\section{NNLO contributions to jet production in DIS}
\label{sec:2}

The QCD factorization theorem allows to write the differential cross section for inclusive
jet production in neutral-current DIS with high momentum transfer $Q^2=-q^2$
as a convolution of the partonic cross section ${\rm d}\sigma_{\gamma a}$
with the parton densities in the proton $f_{a/P}$ and the flux of photons in electrons
$f_{\gamma/e}$ as
\begin{align}
{\rm d}\sigma=&\sum_a\int{\rm d} y\,f_{\gamma/e}(y) \nonumber \\
&\times\int{\rm d} x_P \,f_{a/P}(x_P,\mu_F){\rm d}\sigma_{\gamma a}(\alpha_s,\mu_R,\mu_F) \, ,
\end{align}
where we define $y=(p\cdot q)/(p\cdot k)$ with $p$ and $k$ the momenta of the
incoming proton and electron, respectively, and $q$ the momentum of the exchanged photon.
In deep-inelastic scattering the highly off-shell photon has no time to decay,
so resolved photon contributions can safely be neglected.

From a unified threshold resummation formalism a master formula
can be obtained that permits to compute soft and virtual corrections
to arbitrary partonic hard scattering cross sections \cite{Kidonakis:2003tx}.
At NLO it reads
\begin{equation}
 {\rm d}\sigma_{ab}={\rm d}\sigma_{ab}^B{\alpha_s(\mu)\over\pi}
 \left[  c_3D_1(z)+c_2D_0(z)+c_1\delta(1-z)\right]  \, ,
\end{equation}
where for just one color-charged parton in the initial state we only need the formula
for simple color flow. The functions
\begin{eqnarray}
 D_l(z)&=&\left[ {\ln^l(1-z)\over 1-z}\right] _+
\end{eqnarray}
with decreasing $l$ are the leading and subleading logarithms at
partonic threshold ($z\to1$) in pair-invariant-mass kinematics.
The NNLO master formula is given in the reference cited above,
as are the general formul\ae\ for the
coefficients $c_i$.

We state here the coefficients for the two partonic processes that contribute to jet
production in DIS. For $\gamma^*q\rightarrow qg$, where $\gamma^*$ represents the off-shell
photon, $g$ a gluon and $q$ a quark or an anti-quark, we find
\begin{equation}
 c_3=C_F-N_C \, ,
\end{equation}
\begin{align}
 c_2 =& \,2C_F \ln{\left(\frac{-u}{M^2}\right)}+N_C\ln{\left(\frac{t}{u}\right)} \nonumber \\
      &\qquad\quad-C_F\ln{\left(\frac{\mu_F^2}{M^2}\right)}-\frac{3}{4}C_F-\frac{\beta_0}{4}
\end{align}
and $c_1=c_1^\mu+T_1$ with
\begin{equation}
 c_1^\mu=-\frac{3}{4}C_F\ln{\left(\frac{\mu_F^2}{M^2}\right)}
 +\frac{\beta_0}{4}\ln{\left(\frac{\mu_R^2}{M^2}\right)} \, .
\end{equation}
For the second process $\gamma^* g \rightarrow q \bar{q}$ we find
\begin{equation}
 c_3=2\left( N_C-C_F\right) \, ,
\end{equation}
\begin{equation}
 c_2 = N_C \ln{\left(\frac{tu}{M^4}\right)}-N_C\ln{\left(\frac{\mu_F^2}{M^2}\right)}
 -\frac{3}{2}C_F \, ,
\end{equation}
and
\begin{equation}
 c_1^\mu = \frac{\beta_0}{4}\left[ \ln{\left(\frac{\mu_R^2}{M^2}\right)}
 -\ln{\left(\frac{\mu_F^2}{M^2}\right)} \right] \, .
\end{equation}
These coefficients agree with those found in photoproduction for massless jets for the direct
part
\cite{Klasen:2014}. For massive jets, additional logarithms depending on the jet
radius $R$ appear (e.g.\ in $c_2$) \cite{deFlorian:2007fv}, which are however
irrelevant in the case $R=1$ as in the H1 analysis considered here \cite{Andreev:2014wwa}.
The calculation is analogous to the case of single-jet production in hadron-hadron
collisions \cite{Kidonakis:2000gi}.
The above coefficients further depend on the
QCD color factors $C_F=4/3$ and $N_C=3$, on the one-loop
$\beta$-function $\beta_0=(11N_C-2n_f)/3$ with $n_f$ the number of active quark flavors,
the Mandelstam variables $t$ and $u$, the renormalization and factorization
scales $\mu_R$ and $\mu_F$ and the fixed large invariant scale $M^2$, that in DIS
is equal to $Q^2$, whereas in photoproduction it was the Mandelstam variable $s$. The part
$T_1$ of $c_1$ does not contain any dependence either on the renormalization or the
factorization scale. It includes the NLO virtual corrections and is not predicted by
the threshold resummation formalism. If available, it can be read off from a full NLO
calculation. For the case of transverse photon polarization it could be found in reference
\cite{Potter:1997}. For longitudinal polarization we took the formula directly from the
source code of the corresponding program JetViP \cite{Klasen:1998,Potter:2000}, that
calculates inclusive jet production in DIS to NLO accuracy. Some two-loop quantities
appearing in the NNLO master formula were not given explicitly in \cite{Kidonakis:2003tx}
and could not been found in the respective literature. As in our previous work on jet
photoproduction \cite{Klasen:2014}, they were therefore neglected.

\section{Comparisons to H1 data}
\label{sec:3}

The NNLO contributions have been implemented in the code JetViP for inclusive jet and dijet
production in DIS, where the convolution over $z$ was already included for NLO initial-state
corrections on the proton side. At NLO, we use of course our complete calculation and not
only the logarithmically enhanced terms described above. As a numerical check, we have
repeated the NLO analysis of inclusive single-differential jet production of the H1 collaboration,
performed with NLOJet++ \cite{Nagy:2001xb} and
presented in reference \cite{Andreev:2014wwa,Britzger:2013}, and found excellent agreement,
confirming previous successful comparisons of different NLO programs for jet production
in DIS \cite{Potter:1999gg,Nagy:2001xb}.

The measurement took place during the HERA-II running period with an integrated luminosity of
$\unit{351}{pb^{-1}}$. The beam energies were $\unit{27.6}{GeV}$ for electrons or positrons and
$\unit{920}{GeV}$ for protons, which gives a center-of-mass energy of $\unit{319}{GeV}$. The
leptonic phase space was given by $\unit{150}{GeV^2}<Q^2<\unit{15\,000}{GeV^2}$ and $0.2<y<0.7$.
The jet phase space was restricted to the rapidity interval $-1.0<\eta_{lab}<2.5$, where
$\eta_{lab}$ is the pseudorapidity of a jet in the HERA lab frame. The cross section was
measured differentially in the jet transverse momentum $p_T$ and the virtuality $Q^2$.
Jets were reconstructed using the $k_T$-clustering algorithm \cite{Catani:1993hr} in the
Breit frame, where exlusively electroweak processes can be ruled out by demanding a minimum
of jet transverse momentum (here $p_T>\unit{7}{GeV}$). In inclusive jet production in DIS,
an almost identical fit result for $\alpha_s(M_Z)$ was obtained with the anti-$k_T$ algorithm.
The jet radius was $R=1$. The perturbative scales were chosen to be
\begin{equation}\label{eq:scales}
 \mu_R^2=(Q^2+p_T^2)/2 \quad \text{and} \quad \mu_F^2=Q^2\, .
\end{equation}
The perturbative hard-scattering functions were convoluted with the MSTW2008 set of parton
distribution functions in the proton with different fixed $\alpha_s(M_Z)$ values
\cite{Martin:2009}. This PDF set especially offers the possibility to determine a best-fit
$\alpha_s(M_Z)$. The number of active flavors was $n_f=5$, since sea contributions of a heavy
top quark inside the proton can safely be neglected. Following the H1 analysis, we choose the
PDF member with $\alpha_s(M_Z)=0.118$ in all plots shown in this Letter, which together with
the scale choice from Eq.\ (\ref{eq:scales}) defines our central fit.

\begin{figure}
\centering
\includegraphics[width=\columnwidth]{unnormiert}
\caption{\label{fig:nnloscv}Single-differential inclusive jet cross sections
  as a function of photon virtuality $Q^2$ (top) and jet transverse momentum $p_T$
 (bottom) in NLO (green dashed lines) and aNNLO (blue full lines) with the corresponding
 scale uncertainty bands, obtained by varying $\mu_R$ and $\mu_F$ simultaneously by a
 factor of two up and down, compared to the final H1 data (red points, color online).}
\end{figure}

In Fig.\ \ref{fig:nnloscv}, we compare our NLO (green dashed lines) and aNNLO (blue full lines)
results to the experimental data of the H1 collaboration (red points). The uncertainty bands are
obtained by varying both scales simultaneously
about the central scales up and down by a factor of two as it was also done by the H1
collaboration. We have verified that our NLO results and uncertainty bands agree very well with
those of the H1 analysis. Comparing our calculations to the measurements, we see that the data
points lie in the error bands of both the NLO and the aNNLO calculations.
For the $Q^2$ distribution (top), the scale uncertainty is of similar size at NLO and aNNLO
at low $Q^2$ and is considerably reduced in aNNLO at higher $Q^2$ as expected.
We checked that this tendency continues and becomes more pronounced for even higher $Q^2$,
where unluckily we have no data to compare with.
For the $p_T$ distribution (bottom), the scale uncertainty is not reduced from NLO to aNNLO,
since analytically no logarithms of ratios of $\mu_R$ or $\mu_F$ over $p_T$ appear
(see above). While in our previous analysis of jet photoproduction the jet transverse momentum
$p_T$ ranged from 17 to 71 GeV \cite{Klasen:2014}, its range is restricted in DIS and this
analysis to lower values of 7 to 50 GeV, while the photon virtuality reaches values up to
$(122.5$ GeV)$^2$.

\begin{figure}
\centering
\includegraphics[width=\columnwidth]{normiert}
\caption{\label{fig:nnloscvn}Same as Fig.\ \ref{fig:nnloscv}, but normalized
 to the central NLO predictions.}
\end{figure}

In Fig.\ \ref{fig:nnloscvn}, we show the same comparison normalized to the central NLO result.
We also depict here the central aNNLO result, which is always smaller than the NLO central result
by approximately 6\%.
The central $Q^2$ distribution and even more the central $p_T$ distribution agree better with
the H1 data at aNNLO than at NLO, as they then lie right within the experimental uncertainties.
While the central NLO results overestimate the measured cross sections, indicating that the
value of the strong coupling constant $\alpha_s(M_Z)$ used at this order is too large, the
central aNNLO results underestimate the data and require a slightly larger value of $\alpha_s$.
The data are thus clearly sensitive to the strong coupling constant and can be used for an
extraction not only at NLO, but also at aNNLO. However, we do not expect a significant
reduction of the theoretical error from the scale uncertainty, in particular from the $p_T$
distribution.

\section{Determination of {$\boldmath{\alpha_s}$}}
\label{sec:4}

To determine the strong coupling constant from these comparisons, the
theoretical calculations have to be performed with a set of parton
densities in the proton obtained from global fits assuming different
values of $\alpha_s(M_Z)$. For our analysis at aNNLO, we employ the
latest fits of the MSTW group, which have been obtained
with NNLO running of the coupling, evolution of the parton densities, deep-inelastic
scattering and vector-boson production matrix elements \cite{Martin:2009}.
22 different MSTW2008 NNLO members were used, which correspond to values of
$\alpha_s(M_Z)=0.107$ to 0.127. To compare the aNNLO best-fit $\alpha_s$ to a
corresponding NLO $\alpha_s$ based on our full NLO calculation, we carry out the same
approach with the NLO MSTW2008 parton distribution functions. These are available for the
range of $\alpha_s(M_Z)=0.110$ to 0.130.

The strong coupling constant $\alpha_s$ was determined by comparing the theoretical predictions
at NLO and aNNLO to the experimental measurements by H1 and then finding the minimum value of the
reduced $\chi^2$. We present here our results for the single-differential $p_T$-distribution, but
have verified that fitting the $Q^2$-distribution yields very similar results. At NLO we find
\begin{equation}
 \alpha_s^{\rm NLO}(M_Z)=0.115\pm0.002({\rm exp.})\pm0.005({\rm th.}),
\end{equation}
where the central value is slightly lower than the one obtained by H1 from the unnormalized
double-differential inclusive jet cross section (cf.\ Eq.\ (\ref{eq:3})), but where the
total experimental error and the theoretical error obtained from a simultaneous variation of
the renormalization and factorization scales agree very well. At aNNLO, $\alpha_s$ gets shifted
upwards, since we demonstrated above that the aNNLO contributions reduce the differential cross
sections compared to NLO. We find
\begin{equation}
 \alpha_s^{\rm aNNLO}(M_Z)=0.122\pm0.002({\rm exp.})\pm0.013({\rm th.}),
\end{equation}
where the central value is now slightly above the world average, the experimental error
is of course unchanged, and the theoretical error is slightly larger, reflecting the
observation made above that the aNNLO calculation is not yet sufficiently stabilized by threshold
logarithms at these values of $p_T$ and $Q^2$. The numerical situation would only
improve at higher values of $Q^2$, where unfortunately no experimental data are
available.

\section{Conclusions}
\label{sec:5}

In conclusion, we have presented here a first calculation of inclusive-jet production in
neutral-current deep-inelastic scattering up to NNLO of perturbative QCD.
Leading and subleading logarithmic contributions were extracted from
a unified threshold resummation formalism for virtual photon-parton scattering processes
and shown to agree with those appearing in our full NLO calculations.
The aNNLO contributions implemented in our NLO program improve the description of the
final H1 data on inclusive-jet production in the $Q^2$ distribution and even more in the
$p_T$ distribution, when the world average value of $\alpha_s$ is used and for central
scale choices. The scale uncertainties are reduced only at the highest values of $Q^2$,
where threshold corrections are most important. An aNNLO fit of these data with the
MSTW2008 set of parton densities resulted in a new determination of the strong coupling
constant at the mass of the $Z$-boson that increased the central fit value from below
to above the current world average, but did not reduce the theoretical error.

We thank W.\ Vogelsang for useful discussions.
This work has been supported by the BMBF Theorie-Verbund ``Begleitende
theoretische Untersuchungen zu den Experimenten an den Gro\ss{}ger\"aten der
Teilchenphysik.''