text:
\section{Introduction}
\label{sec:Intro}

LHCb has reported anomalies in the measured decay rates of the $B$ meson, which have been interpreted as hints of lepton flavour universality violation \cite{Aaij:2014pli,Aaij:2015esa}. The SM predicts equal rates for the processes $B\to K^{(*)} \mu^+\mu^-$ and $B\to K^{(*)} e^+e^-$, and it is customary to study the ratios of these branching ratios, defined as $R(K)$ and $R(K^*)$, since the dependencies on hadronic matrix elements (and associated uncertainties) cancel out \cite{Hiller:2003js}. The measurements of these hadronically clean observables deviate consistently (although perhaps with not enough statistical significance) from the SM prediction $R(K^{(*)})=1$ \cite{Bordone:2016gaq}. These hints are complemented by measurements of other observables that are more sensitive to hadronic physics. In particular, the differential branching fractions \cite{Wei:2009zv,Aaij:2014pli,Aaij:2015esa} and angular observables \cite{Khachatryan:2015isa,Abdesselam:2016llu,Lees:2015ymt,Sirunyan:2017dhj,Aaboud:2018krd,Wehle:2016yoi,Aaltonen:2011ja,Aaij:2015oid} associated to the processes $B\to \phi \mu^+\mu^-$ and $B\to K^{(*)} \mu^+\mu^-$ also deviate from the SM predictions. Interestingly, all the apparent anomalies involve the transition $b \rightarrow s\mu^{+}\mu^{-}$.

In order to account for these experimental results, one can modify the SM effective Hamiltonian, which involves penguin and box diagrams, by including one-loop contributions from new exotic particles. A full classification of the various particle combinations, considering different gauge representations, was presented in  Refs.~\cite{Gripaios:2015gra,Arnan:2016cpy}. Among the different models, some featured neutral scalar or fermions that, if stable, could play the role of dark matter (DM)\footnote{An alternative to this one-loop solution is to consider $ Z^\prime$~\cite{Buras:2013qja,Gauld:2013qja} or leptoquark~\cite{Bauer:2015knc,Angelescu:2018tyl} tree-level contributions, see e.g., Ref.~\cite{Capdevila:2017bsm} and references therein. The DM problem has been addressed in the framework of these constructions \cite{Vicente:2018xbv}, see e.g.,  Refs.~\cite{Sierra:2015fma,Belanger:2015nma,Altmannshofer:2016jzy,Celis:2016ayl,Cline:2017lvv,Ko:2017yrd,Ellis:2017nrp,Baek:2017sew,Fuyuto:2017sys,Cox:2017rgn,Falkowski:2018dsl,Darme:2018hqg,Singirala:2018mio,Baek:2018aru,Kamada:2018kmi} for the $Z^\prime$, and Refs.~\cite{Hati:2018fzc,Choi:2018stw,Cline:2017aed,Varzielas:2015sno} for the leptoquark models.}. The first possibility was investigated in Ref.~\cite{Kawamura:2017ecz}, where it was found that the large new couplings required to reproduce the correct DM relic abundance induce sizeable 1-loop contributions to DM-nucleon scattering, leading to very strong limits from direct detection experiments. In addition, as reported by \cite{Bhattacharya:2015xha}, the Higgs portal coupling typically dominates over other new physics effects. The second possibility was addressed in Ref.~\cite{Cline:2017qqu}, where the fermionic dark matter field was accompanied by one additional scalar and one additional coloured fermion.

In this work, we consider a modification of the model of Ref.~\cite{Cline:2017qqu}. Namely, we will also assume a fermionic dark matter particle, but with two extra scalar fields, one of which has a colour charge. On top of this, we include the latest SM theoretical prediction for the mass difference in $B_{s}-$mixing \cite{DiLuzio:2017fdq}, which differs from the experimental observation by $1.8\,\sigma$. In order to reduce this tension and provide an explanation for the $B$ anomalies, complex couplings are needed, leading to new CP-violation sources, a scenario that has not been studied in the context of one-loop models so far. We explore the parameter space of this model, taking into account all the flavour observables, DM constraints, and LHC collider signatures.

This paper is organised as follows. In Section~\ref{sec:model}, we introduce the details of the particle physics model, address the constraints from the observed DM relic abundance and $B_{s}-$mixing and discuss the implications on the model's parameter space. In Section~\ref{sec:lhc}, we investigate the possibility of observing this scenario at the LHC, for which we take into account dijet and dimuon searches. We also include a projection of the potential reach of the High Luminosity phase of the LHC.  Finally, in Section~\ref{sec:direct}, we compute the DM-nucleus scattering cross section and study current constraints and the future reach of direct DM detection experiments. The conclusions are presented in Section~\ref{sec:conclusions}.

\section{The model}
\label{sec:model}

In this article, we consider a model in which the DM particle is a Majorana fermion, $\chi$, with two extra scalar fields, $\phi_q$ and $\phi_l$, which couple to left-handed quarks and leptons, respectively\footnote{As we will comment in Section \ref{sec:direct}, the alternative construction with Dirac DM is ruled out mainly by experimental results from direct DM detection.}. The interactions between the new particles and the SM are described by the Lagrangian,
\begin{equation}
  \pazocal{L}_{\text{int}}^{\text{NP}} = \lambda_{Q_i}\bar{Q}_{i}\phi_{q}P_{R}\chi + \lambda_{L_i}\bar{L}_{i}\phi_{l}P_{R}\chi + \rm{h.c.} \, ,
  \label{eq:lagrangian}
\end{equation}
where $Q_i$ and $L_i$ denote the SM left-handed quark and lepton doublets of each generation, and $\lambda_{Q_i}$ and $\lambda_{L_i}$ are the corresponding new couplings. The quantum numbers for the new fields are summarised in Table~\ref{tab:quantum_numbers}. We impose a ${\cal Z}_2$ parity under which the SM fields are invariant, and which guarantees the stability of the DM candidate, as long as $m_{\phi_{q,l}} > m_{\chi}$. Upon rotation from the electroweak to the quark mass eigenbasis, the couplings $\lambda_{Q_i}$ are rotated in flavour space. Assuming that the electroweak and mass eigenbasis are aligned for the leptons and down-type quarks, the couplings to the up-type quarks are generated by the CKM rotation as follows:
\begin{equation}
  \lambda_{Q_i}\bar{Q}_{i} \rightarrow \lambda_{Q_j}(\bar{u}_{L,i}V_{ij}, \bar{d}_{L,j}) \, .
\label{eq:rotation}
\end{equation}
From now on, we will denote the couplings in the mass eigenbasis with the corresponding quark or lepton label. These couplings are, in general, complex.

\begin{table}[t!]
\begin{center}
\begin{tabular}{|c|c|c|c||c|}
\hline
& $SU(3)$ & $SU(2)_{L}$ & $U(1)_{Y}$ & ${\cal Z}_2$ \\
\hline
\hline $\phi_{q}$ & $3$ & $2$ & $1/6$ & $-1$ \\
\hline $\phi_{l}$ & $1$ & $2$ & $-1/2$ & $-1$ \\
\hline $\chi$ & $1$ & $1$ & $0$ & $-1$ \\
\hline
\end{tabular}
\caption{Quantum numbers of the new fields. We also indicate the charges under ${\cal Z}_2$.}
\label{tab:quantum_numbers}
\end{center}
\end{table}

\begin{figure}[!t]
    \begin{center}
	\includegraphics[scale=0.8]{plots-final/bsmumu_paper.pdf}
    \end{center}
	\caption{One-loop diagram contribution from the new particles to the $b \rightarrow s\mu^{+}\mu^{-}$ transitions.}
    \label{fig:box}
\end{figure}

This model induces new physics contribution to flavour observables at the one loop level. In particular, a new box diagram appears for the  $b \rightarrow s\mu^{+}\mu^{-}$ transition, as shown in Figure~\ref{fig:box}. These effects  can be described using an effective field theory approach, thus parameterising the new contributions as corrections to the corresponding Wilson coefficients of the effective Hamiltonian,
\begin{equation}
  \pazocal{H}_{\text{eff}}^{\mu^{+}\mu^{-}} = -\frac{4G_{F}}{\sqrt{2}}V_{tb}V_{ts}^{*}(C_{9}\pazocal{O}_{9} + C_{9}^{\prime}\pazocal{O}_{9}^{\prime} + C_{10}\pazocal{O}_{10} + C_{10}^{\prime}\pazocal{O}_{10}^{\prime}) + \text{h.c.} \, ,
  \label{eq:hamiltonian}
\end{equation}
where the effective operators $\pazocal{O}_{9}$, $\pazocal{O}_{9}^{\prime}$, $\pazocal{O}_{10}$, $\pazocal{O}_{10}^{\prime}$ are defined as:
\begin{eqnarray}
  \pazocal{O}_{9} &=& \frac{\alpha_{\text{em}}}{4\pi}(\bar{s}\gamma^{\mu}P_{L}b)(\bar{\mu}\gamma_{\mu}\mu) \, , \\
  \pazocal{O}_{9}^{\prime} &=& \frac{\alpha_{\text{em}}}{4\pi}(\bar{s}\gamma^{\mu}P_{R}b)(\bar{\mu}\gamma_{\mu}\mu) \, ,\\
  \pazocal{O}_{10} &=& \frac{\alpha_{\text{em}}}{4\pi}(\bar{s}\gamma^{\mu}P_{L}b)(\bar{\mu}\gamma_{\mu}\gamma_{5}\mu) \, , \\
  \pazocal{O}_{10}^{\prime} &=& \frac{\alpha_{\text{em}}}{4\pi}(\bar{s}\gamma^{\mu}P_{R}b)(\bar{\mu}\gamma_{\mu}\gamma_{5}\mu) \, .
\end{eqnarray}
The Wilson coefficients $C_{9}$, $C_{9}^{\prime}$, $C_{10}$, $C_{10}^{\prime}$ contain both the SM and new physics (NP) contributions,
\begin{eqnarray}
  C_{9} &=&C_{9}^{\text{SM}} + C_{9}^{\text{NP}} \, ,\nonumber \\
  C_{10} &=& C_{10}^{\text{SM}} + C_{10}^{\text{NP}} \, ,
  \label{eq:c9c10}
\end{eqnarray}
with the primed coefficients defined in an equivalent way.

Global fits \cite{Descotes-Genon:2015uva,Hurth:2016fbr,Capdevila:2017bsm,Altmannshofer:2017yso,DAmico:2017mtc,Hiller:2017bzc,Geng:2017svp,Ciuchini:2017mik,Alok:2017sui,Hurth:2017hxg} have been used to determine the new physics contribution to the Wilson coefficients in order to reproduce the observed experimental results. These fits favour $C_{9}^{\text{NP}} = -C_{10}^{\text{NP}}$, and suggest that no new physics is required for operators involving electrons or tau leptons. Motivated by these results, we assume negligible couplings to the first quark generation (i.e., $\lambda_{Q_1} = 0$) and to the first and third lepton generations (i.e., $\lambda_{e}=\lambda_{\tau}=0$). This provides an explanation for the $R_{K^{(*)}}$ anomalies, while relaxing the bounds from other searches.

Therefore, in total, we are left with six free parameters in this model, namely the masses of the three new particles ($m_{\chi}$, $m_{\phi_l}$, $m_{\phi_q}$), and the couplings to $b-$type quarks, $s-$type quarks, and leptons ($\lambda_b$, $\lambda_s$, $\lambda_\mu$).

It should be noted that the couplings $\lambda_{1}|\phi_{l}|^{2}|H|^{2}$ and $\lambda_{2}|\phi_{q}|^{2}|H|^{2}$ are allowed by gauge symmetry in the Lagrangian of Eq.~(\ref{eq:lagrangian}). However, they only lead to an overall shift to the masses of $\phi_{l}$ and $\phi_{q}$ after electroweak symmetry breaking since the couplings to the Higgs play no phenomenological role in the relevant range of $\phi_{l,q}$ masses. Likewise, the terms $\lambda_{3}|\phi_{l}H|^{2}$ and $\lambda_{4}|\phi_{q}H|^{2}$ are also allowed by gauge symmetry. They typically induce a small split in the masses of the neutral and charged components of the doublets $\phi_{l}$ and $\phi_{q}$ in the range of $\phi_{l,q}$ masses that survive the collider constraints. Finally, a term of the form $(\phi_{l}H)^{2}$ can lead to large contributions to neutrino masses at one loop, which forces the corresponding coupling to be extremely small \cite{Cline:2017qqu}. We will neglect these couplings in the following.

As mentioned in the Introduction, similar models have been discussed in the literature, featuring either scalar DM \cite{Kawamura:2017ecz,Chiang:2017zkh,Barman:2018jhz,Grinstein:2018fgb} or fermionic DM \cite{Cline:2017qqu}. Our model differs from that of Ref.~\cite{Cline:2017qqu} in that we have two extra scalar fields which couple to the lepton or quark sectors.

\subsection{Dark matter relic abundance}
\label{sec:relic}

In order for $\chi$ to be a viable DM candidate, it must reproduce the observed relic abundance, which can be inferred from Planck satellite data to be $\Omega h^2=0.1199\pm 0.0022$ \cite{Ade:2015xua}. The pair-annihilation proceeds through the two $t-$channel diagrams with $\phi_q$ and $\phi_l$, shown in Figure~\ref{fig:annihilation}.

\begin{figure}[!t]
    \begin{center}
    \includegraphics[scale=0.82]{plots-final/relic_phiq.pdf}
	\includegraphics[scale=0.82]{plots-final/relic_phil.pdf}
    \end{center}
	\caption{Tree-level contributions to the DM pair annihilation.}
    \label{fig:annihilation}
\end{figure}

The stringent flavour constraints force the couplings to quarks to be much smaller than the couplings to leptons (muons and neutrinos), and the combination of flavour and collider bounds impose $m_{\phi_q}>m_{\phi_l}$, with coloured scalars generally above 1~TeV. Therefore DM annihilation into a $\mu^-\mu^+$ or $\nu_\mu \bar\nu_\mu$ pair is the dominant channel. The thermally-averaged annihilation cross section, $\langle\sigma v\rangle$, can be expressed as a plane wave expansion in terms of the dimensionless parameter $x=m_{\chi}/T$. For the case of a Majorana fermion, the zero-velocity term is helicity suppressed, and the leading contribution comes from the linear term in $1/x$ \cite{Chang:2014tea},
\begin{equation}
  \langle\sigma v\rangle = 2 \frac{|\lambda_{\mu}|^{4}m_{\chi}^{2}\left(m_{\phi_{l}}^{4} + m_{\chi}^{4}\right)}{16\pi \left(m_{\phi_{l}}^{2} + m_{\chi}^{2}\right)^{4}} \, \frac{1}{x}\ ,
  \label{eq:sigmav}
\end{equation}
where we have neglected the muon and the neutrino masses. In order to reproduce the correct relic abundance, we can now impose $\langle\sigma v\rangle = 2.2\times 10^{-26}$~cm$^3$~s$^{-1}$ (where $x\sim 20$ at freeze-out).

We will use this relation to fix $m_{\phi_l}$ as a function of the other parameters, thus effectively reducing by one the number of free parameters. Furthermore, due to the suppression of the velocity-independent term for $\langle\sigma v\rangle$, indirect detection bounds are not expected to constrain our model.

\subsection{$B_{s}-$mixing and other flavour constraints}
\label{sec:bsmixing}

This model introduces new couplings to the $s$ and $b$ quarks (and to the rest of the quarks by rotation of the CKM matrix). We must therefore incorporate constraints from $B$ meson physics.

The most relevant bounds are those that involve $b \rightarrow s\mu^{+}\mu^{-}$ transitions. The new physics contribution to the Wilson coefficient comes from box and photon-penguin diagrams \cite{Gripaios:2015gra,Arnan:2016cpy}, $C_{9}^{\rm NP}=C_{9}^{\text{box}}+C_{9}^{\gamma}$, with\footnote{We have neglected the $Z$-penguin contribution to $C_{9}^{\text{NP}}$, since it is suppressed by $(m_{b}/m_{Z})^{2}$ and is subdominant compared to the photon exchange.}
\begin{eqnarray}
  C_{9}^{\text{box}} &=& \frac{\sqrt{2}}{128\pi\alpha_{\text{em}}G_{F}m^{2}_{\psi}} \frac{\lambda_{s}\lambda_{b}^{*}}{V_{tb}V_{ts}^{*}} |\lambda_{\mu}|^{2} \left(F(x_{q},x_{l}) + 2G(x_{q},x_{l}) \right) \, ,\nonumber\\
  C_{9}^{\gamma} &=& \frac{\sqrt{2}}{8G_{F}m^{2}_{\psi}} \frac{\lambda_{s}\lambda_{b}^{*}}{V_{tb}V_{ts}^{*}} F_{9}(x_{q}) \, ,
\end{eqnarray}
where we have defined the dimensionless variables $x_{q} = m_{\phi_{q}}^{2} / m_{\chi}^{2}$ and $x_{l} = m_{\phi_{l}}^{2} / m_{\chi}^{2}$, and the loop functions are:
\begin{eqnarray}
  F(x,y) &=& \frac{1}{(1-x)(1-y)} + \frac{x^{2}\log x}{(1-x)^{2}(x-y)} + \frac{y^{2}\log y}{(1-y)^{2}(y-x)} \, ,\nonumber\\
  G(x,y) &=& \frac{1}{(1-x)(1-y)} + \frac{x\log x}{(1-x)^{2}(x-y)} + \frac{y\log y}{(1-y)^{2}(y-x)} \, ,\nonumber\\
  F_{9}(x) &=& \frac{-2x^{3} + 9x^{2} - 18x + 6\log x + 11}{36(x-1)^{4}} \, .
\label{eq:fg}
\end{eqnarray}
The term  $G(x_{q},x_{l})$ vanishes if $\chi$ is a Dirac particle.

In order to constrain the Wilson coefficients we use the first global fit that takes into account the possibility that $C_{9}$ and $C_{10}$ are complex \cite{Alok:2017jgr}. This is a natural scenario that arises when new CP-violation sources are introduced, and has not been studied in detail in the literature so far.

Likewise, the new physics contribution to $B_{s}-$mixing can be parameterised in terms of an effective Hamiltonian,
\begin{equation}\label{eq:3.2.1}
  \pazocal{H}_{\text{eff}}^{b\bar{s}} = C_{B\bar{B}}^{\text{NP}} \, (\bar{s}_{\alpha}\gamma^{\mu}P_{L}b_{\alpha})(\bar{s}_{\beta}\gamma_{\mu}P_{L}b_{\beta}) \, ,
\end{equation}
where $\alpha$ and $\beta$ are colour indices. The new physics contribution to the Wilson coefficient is given by
\begin{equation}\label{eq:wilsonCBB}
  C_{B\bar{B}}^{\text{NP}} = \frac{1}{128\pi^{2}m^{2}_{\psi}} (\lambda_{s}\lambda_{b}^{*})^{2} \left(F(x_{q},x_{q}) + 2G(x_{q},x_{q}) \right) \, ,
\end{equation}
where the loop functions $F$ and $G$ were already defined in Eq.~(\ref{eq:fg}).

In order to quantify the allowed magnitude of the Wilson coefficient $C_{B\bar{B}}^{\text{NP}}$, we follow the steps of \cite{DiLuzio:2017fdq} and introduce a complex parameter $\Delta$ in the following way:
\begin{equation}\label{eq:3.2.3}
  \frac{M_{12}^{\text{SM}} + M_{12}^{\text{NP}}}{M_{12}^{\text{SM}}} \equiv |\Delta|e^{i\phi_{\Delta}} \, ,
\end{equation}
where $M_{12}^{\text{SM}}$ and $M_{12}^{\text{NP}}$ describe the SM and new physics contributions to $B_{s}-$mixing, and their values are given by the corresponding box diagrams. The complex phase, $\phi_{\Delta}$, quantifies the CP-violating effects introduced by the imaginary parts of the new couplings. We find:
\begin{eqnarray}
  |\Delta| &=& \frac{\Delta M_{s}^{\text{exp}}}{\Delta M_{s}^{\text{SM}}} = \left|1+\frac{C_{B\bar{B}}^{\text{NP}}}{C_{B\bar{B}}^{\text{SM}}} \right| \, ,\nonumber\\
  \phi_{\Delta} &=& \text{Arg} \left(1+\frac{C_{B\bar{B}}^{\text{NP}}}{C_{B\bar{B}}^{\text{SM}}} \right) \, ,
  \label{eq:deltadef}
\end{eqnarray}
where $\Delta M_{s}$ is the mass difference of the mass eigenstates of the $B_s$ meson.

The parameter $|\Delta|$ can be constrained using the most precise experimental measurement of $\Delta M_{s}$\cite{Amhis:2016xyh} and the last update on its theoretical prediction \cite{DiLuzio:2017fdq}, which show a $1.8 \sigma$ difference,
\begin{eqnarray}
\Delta M_{s}^{\text{exp}} &=& (17.757 \pm 0.021) \, \text{ps}^{-1} \, ,\nonumber\\
\Delta M_{s}^{\text{SM}} &=& (20.01 \pm 1.25) \, \text{ps}^{-1} \, .
\label{eq:deltams}
\end{eqnarray}
The dominant uncertainties in the calculation of $\Delta M_{s}^{\text{SM}}$ come from lattice predictions for the non-perturbative bag parameter, $\mathcal{B}$, and decay constant, $f_{B_{s}}$, and to a lesser extent from the uncertainty in the values of CKM elements. Both of these errors have been considerably reduced since the last theory update for the mass difference \cite{Artuso:2015swg}. The last average given by the lattice community \cite{Aoki:2016frl} gives significantly more precise values for $\mathcal{B}$ and $f_{B_{s}}$.

From these values, one can infer $|\Delta| = 0.887 \pm 0.055$, and using the data provided in Ref.~\cite{DiLuzio:2017fdq} we obtain $C_{B\bar{B}}^{\text{SM}} = 4.897 \times 10^{-5} \; {\mathrm{TeV}}^{-2} \,$. Using Eq. (\ref{eq:deltadef}) we find that the Wilson coefficient has to satisfy
\begin{equation}
  \sqrt{\left(1 + \frac{\text{Re} \, C_{B\bar{B}}^{\text{NP}}}{C_{B\bar{B}}^{\text{SM}}} \right)^{2} +          \left(\frac{\text{Im} \, C_{B\bar{B}}^{\text{NP}}}{C_{B\bar{B}}^{\text{SM}}} \right)^{2}} \in [0.777, 0.998] \quad (2\sigma) \, .
  \label{eq:cbb}
\end{equation}
CP-violating effects are further constrained by the CP asymmetry of the golden mode $B_{s} \rightarrow J/\psi \, \phi$ \cite{Amhis:2016xyh},
\begin{equation}
A_{\text{CP}}^{\text{mix}}(B_{s} \rightarrow J/\psi \phi) = \sin(\phi_{\Delta} - 2\beta_{s}) \, = -0.021 \pm 0.031\ ,
\label{eq:acp}
\end{equation}
where $\beta_{s} = 0.01852 \pm 0.00032$\cite{Charles:2004jd}, and penguin contributions are neglected.
Using Eq.~(\ref{eq:deltadef}), this can be interpreted as an additional constraint on the real and imaginary parts of $C_{B\bar{B}}^{\text{NP}}$ (and in turn, on the real and imaginary parts of the couplings $\lambda_{s}\lambda_{b}^{*}$).

\begin{figure}[!t]
\begin{center}
	\includegraphics[scale=.4]{plots-final/Re_Im_Majorana_10_1500_4.pdf}
	\includegraphics[scale=.4]{plots-final/Re_Im_Majorana_50_1500_4.pdf}
	\includegraphics[scale=.4]{plots-final/Re_Im_Majorana_150_1500_4.pdf}\\ \vspace*{4ex}
	\includegraphics[scale=.4]{plots-final/Re_Im_Majorana_10_2500_4.pdf}
    \includegraphics[scale=.4]{plots-final/Re_Im_Majorana_50_2500_4.pdf}
    \includegraphics[scale=.4]{plots-final/Re_Im_Majorana_150_2500_4.pdf}
\end{center}
	\caption{
	 The dark (light) green area is the $1\sigma$ ($2\sigma$) allowed region by $b \rightarrow s \mu^+ \mu^- $  observables in the $(\operatorname{Re}(\lambda_s  \lambda_b^*), \operatorname{Im}(\lambda_s  \lambda_b^*) )$ plane. Dark (light) blue regions correspond to  $1\sigma \; (2\sigma)$  $B_s-$ mixing allowed regions. We take $\lambda_\mu = \sqrt{4\pi} $ and $m_{\phi_q} = 1.5 \, {\mathrm{TeV}}$ (top row),  $2.5 \, {\mathrm{TeV}}$ (bottom row). The  specific values of $m_\chi$, $m_{\phi_q}$ are given in the plot and $m_{\phi_l}$ is fixed to reproduce the measured DM relic abundance.
	}
	\label{fig:greensandblues1}
\end{figure}

In Figure~\ref{fig:greensandblues1}, the effect of all of these constraints on the real and imaginary parts of the couplings $\lambda_{s}\lambda_{b}^{*}$ for several benchmark points is shown. Regions that are allowed by $b \rightarrow s\mu^{+}\mu^{-}$ observables and $B_{s}-$mixing (given by Eqs.~(\ref{eq:cbb}) and (\ref{eq:acp})) are shaded in green and blue, respectively. For illustrative purposes, the figure shows the constraints for multiple values of DM and mediator masses, while keeping $\lambda_\mu=\sqrt{4\pi}$ fixed. We remind the reader that the mass of $m_{\phi_l}$ is fixed so as to reproduce the correct relic density using Eq.~(\ref{eq:sigmav}).

As we can observe, in order to simultaneously satisfy both types of constraints, complex couplings are needed ($\text{Im}(\lambda_{s}\lambda_{b}^{*})\ne 0$). Also, as the mass of the dark matter particle and the mediators increase, both areas are more difficult to reconcile. In practise, this leads to an upper bound on the masses of the exotic new particles. The precise limit depends on the choice of couplings, which we will discuss in Section~\ref{sec:lhc}.

Finally, the new physics couplings to the up-type quarks are generated via CKM rotation,
\begin{eqnarray}
\lambda_{u} &=& V_{us}\lambda_{s} + V_{ub}\lambda_{b} \, , \nonumber\\
\lambda_{c} &=& V_{cs}\lambda_{s} + V_{cb}\lambda_{b} \, .
\label{eq:upcharmCKM}
\end{eqnarray}
These couplings generate a new physics contribution to $D^{0}-$mixing, and the Wilson coefficient $C_{D\bar{D}}^{\text{NP}}$ is obtained replacing $\lambda_{s}$ and $\lambda_{b}^{*}$ in Eq.~(\ref{eq:wilsonCBB}) by $\lambda_{u}$ and $\lambda_{c}^{*}$, respectively.

In contrast to $B_{s}-$mixing, there is no precise theory determination for the mass difference in the $D^{0}$ system. Therefore, in order to constrain the new physics contribution to $C_{D\bar{D}}$ we use the measured value of the mass difference in $D^{0}-$mixing. The experimental bound on the mixing diagram is given by \cite{Bona:2017kam}
\begin{equation}\label{eq:M12exp}
  |M_{12}|_{D\bar{D}}^{\text{exp}} \in [0.6, 7.5] \times 10^{-3} \, \text{ps}^{-1} \quad (2\sigma) \, ,
\end{equation}
whereas the new physics contribution to $D^{0}-$mixing is described by
\begin{equation}\label{eq:M12lattice}
  |M_{12}|_{D\bar{D}} = \frac{|C_{D\bar{D}}|}{2M_{D^{0}}}\langle{D^{0}}| \pazocal{O} |{\bar{D}^{0}}\rangle \, ,
\end{equation}
where $\pazocal{O}$ is a combination of operators containing all possible SM and new physics contributions to $D^{0}-$mixing. Using the last results from \cite{Bazavov:2017weg} we get the following bound on the Wilson coefficient:
\begin{equation}\label{eq:CDD}
  |C_{D\bar{D}}^{\text{exp}}| \le 5.695 \times 10^{-8} \; {\mathrm{TeV}}^{-2} \quad (2\sigma) \, .
\end{equation}

Although this model induces new physics contributions to other flavour observables (such as $b \rightarrow s\gamma$, $b \rightarrow s\nu\bar{\nu}$ and effective $Z\mu^{+}\mu^{-}$ and $Zq_{i}q_{j}$ couplings), their size is very small and does not produce significant deviations from current experimental searches.

\subsection{Benchmark scenarios}
\label{sec:benchmark}

All the new physics contributions to the observables described above depend on five independent parameters: the three masses of the new particles, $m_{\chi}$, $m_{\phi_{q}}$ and $m_{\phi_{l}}$, the product of the couplings $\lambda_{s}\lambda_{b}^{*}$ and the absolute value of the coupling $|\lambda_{\mu}|$.

The three masses only enter the Wilson coefficients through the factor $m_{\chi}^{-2}$ and the dimensionless loop functions. In addition, all the Wilson coefficients are proportional to $\lambda_{s}\lambda_{b}^{*}$ or $|\lambda_{\mu}|^{2}$ or both. In order to constrain our model, we consider two scenarios by fixing the value of $|\lambda_{\mu}|$. Then we scan over the mass parameters $m_{\chi}$ and $m_{\phi_{q}}$, with $m_{\phi_{l}}$ fixed by the requirement of reproducing the correct relic abundance, and check all the flavour observables described in Section \ref*{sec:bsmixing}. In this way, for any combination of masses and a fixed value of $|\lambda_{\mu}|$ we get a set of allowed values for $\lambda_{s}\lambda_{b}^{*}$.
We consider two hierarchies between $|\lambda_{s}|$ and $|\lambda_{b}|$ that lead to different constraints from $D^{0}-$mixing, and, ensuring that $\text{Im}(\lambda_s\lambda_b^*)\ne0$, we define the following benchmark scenarios:

\begin{minipage}{0.48\linewidth}
    \begin{itemize}
        \item[(A1)] $|\lambda_{\mu}| = 2$, with $\lambda_{b} = \lambda_{s}^*$;

        \item[(A2)] $|\lambda_{\mu}| = 2$, with $\lambda_{b} = 4\lambda_{s}^*$;
    \end{itemize}
\end{minipage}
\begin{minipage}{0.48\linewidth}
    \begin{itemize}
        \item[(B1)]  $|\lambda_{\mu}| = \sqrt{4\pi}$, with $\lambda_{b} = \lambda_{s}^*$;

        \item[(B2)]  $|\lambda_{\mu}| = \sqrt{4\pi}$, with $\lambda_{b} = 4\lambda_{s}^*$,
    \end{itemize}
\end{minipage}

\noindent
where $|\lambda_{\mu}| = \sqrt{4\pi}$ is the perturbative limit. After establishing a hierarchy between $|\lambda_{s}|$ and $|\lambda_{b}|$, we calculate their maximum and minimum allowed values from the corresponding maximum and minimum allowed values of $\lambda_{s}\lambda_{b}^{*}$. Scenarios with $|\lambda_{s}| > |\lambda_{b}|$ are excluded by $D^{0}-$mixing constraints. Likewise, as we will see in Section~\ref{sec:lhc}, smaller values of $\lambda_\mu$ are constrained by LHC bounds.

\section{LHC constraints and prospects for high-luminosity}
\label{sec:lhc}

In this section, we study the experimental signatures that this model would produce at the LHC. DM search strategies in both ATLAS and CMS involve analysing final states containing jets and leptons produced in association with a DM particle, identified from missing transverse energy. In this model, direct production of the coloured and leptonic scalar doublets $\phi_q$ and $\phi_l$, respectively, typically leads to such final states.

Let us first consider production processes that involve the coloured scalar, $\phi_q$. In this case, our model could lead to visible signals in final states with both $\text{monojet} \, / \, \text{dijet} + E_{\rm T} \hspace{-1.2em}/\;\;$ signatures. When the new physics coupling $\lambda_q$ is smaller than the strong interaction coupling, $\alpha_{\text{QCD}}$, pure QCD processes constitute the main contribution to the cross section \cite{An:2013xka}. In this model, this implies that QCD diagrams dominate over those with new physics couplings. As a consequence, monojet searches for this model are less effective than dijet searches and we will concentrate on the latter. The $\text{dijet} + E_{\rm T} \hspace{-1.2em}/\;\;$ processes are shown in Figure~\ref{fig:prod-colour}, where diagrams (a) correspond to the QCD contributions, and diagrams (b) and (c) involve new physics couplings. The main production channel is the pair production of the coloured scalar particles, that subsequently decays into a DM particle and a quark,
\begin{equation}
pp \to \phi_q\phi_q^* \, / \, \phi_q\phi_q \, / \, \phi_q^*\phi_q^* \to q q + E_{\rm T} \hspace{-1.2em}/\;\; \; .
\end{equation}
In addition, the scalar doublet $\phi_q$ has the same quantum numbers as squarks in supersymmetric (SUSY) models. Therefore, the kinematics in its production and decay in diagrams (a) of Figure~\ref{fig:prod-colour} mimic those of squarks in SUSY models with decoupled gluinos. As a consequence, limits from ATLAS and CMS squark searches can be used to constrain the model.

One can also consider the pair production of the leptonic scalar, $\phi_l$. In this case, the production process is mediated by $W$ or $Z$ bosons and involves the electroweak coupling, as shown in Figure~\ref{fig:prod-ew}. The decays of $\phi_l$ lead to clean final states with one or two leptons and missing energy. Although flavour constraints require $\lambda_\mu \gg \lambda_q$, the cross section of this process is smaller than the production of the coloured mediator for similar mediator masses. However, since $m_{\phi_l}$ is fixed for every value of $m_{\chi}$ to reproduce the correct relic abundance, there are regions of the parameter space where both searches are complementary. We will here consider the  process
\begin{equation}
pp \to \phi_l\phi_l^* \to \mu \mu / \mu \nu + E_{\rm T} \hspace{-1.2em}/\;\; \; ,
\end{equation}
where the dimuon channel leads to the strongest constraints. As in the previous case, we can exploit the analogy between $\phi_l$ and sleptons to use the limits from slepton searches to constrain this model.

\begin{figure}[t!]
\centering
\renewcommand{\thesubfigure}{a1}
\begin{subfigure}[t]{.32\textwidth}
  \centering
  \includegraphics[width=\linewidth]{plots-final/production1.pdf}
  \caption{}
\end{subfigure}
\hspace{-10pt}
\renewcommand{\thesubfigure}{a2}
\begin{subfigure}[t]{.32\textwidth}
  \centering
  \includegraphics[width=\linewidth]{plots-final/production2.pdf}
  \caption{}
\end{subfigure}
\hspace{-10pt}
\renewcommand{\thesubfigure}{a3}
\begin{subfigure}[t]{.32\textwidth}
  \centering
  \includegraphics[width=\linewidth]{plots-final/production3.pdf}
  \caption{}
\end{subfigure}

\vspace{0pt}
\renewcommand{\thesubfigure}{a4}
\begin{subfigure}[t]{0.32\textwidth}
  \centering
  \includegraphics[width=\linewidth]{plots-final/production4.pdf}
  \caption{}
\end{subfigure}
\hspace{-10pt}
\renewcommand{\thesubfigure}{b}
\begin{subfigure}[t]{0.32\textwidth}
  \centering
  \includegraphics[width=\linewidth]{plots-final/production5.pdf}
  \caption{}
\end{subfigure}
\hspace{-10pt}
\renewcommand{\thesubfigure}{c1}
\begin{subfigure}[t]{0.32\textwidth}
  \centering
  \includegraphics[width=\linewidth]{plots-final/production6.pdf}
  \caption{}
\end{subfigure}

\vspace{0pt}
\renewcommand{\thesubfigure}{c2}
\begin{subfigure}[t]{0.32\textwidth}
  \centering
  \includegraphics[width=\linewidth]{plots-final/production7.pdf}
  \caption{}
\end{subfigure}
\hspace{-10pt}
\renewcommand{\thesubfigure}{c3}
\begin{subfigure}[t]{0.32\textwidth}
  \centering
  \includegraphics[width=\linewidth]{plots-final/production8.pdf}
  \caption{}
\end{subfigure}
\hspace{-10pt}
\renewcommand{\thesubfigure}{c4}
\begin{subfigure}[t]{0.32\textwidth}
  \centering
  \includegraphics[width=\linewidth]{plots-final/production9.pdf}
  \caption{}
\end{subfigure}
\caption{Diagrams for the pair production of the coloured scalar mediator, $\phi_q$,  leading to $\text{dijet} + E_{\rm T} \hspace{-1.2em}/\;\;$ signatures in the final state.  Diagrams (a1)--(a4) are generated by purely QCD interactions, and diagrams (b), (c1)--(c4) are generated by DM t-channel exchange.}\label{fig:prod-colour}
\end{figure}

\begin{figure}[t!]
\centering
\renewcommand{\thesubfigure}{d1}
\begin{subfigure}[t]{.32\textwidth}
  \centering
  \includegraphics[width=\linewidth]{plots-final/productionew2.pdf}
  \caption{}
\end{subfigure}
\hspace{0pt}
\renewcommand{\thesubfigure}{d2}
\begin{subfigure}[t]{.32\textwidth}
  \centering
  \includegraphics[width=\linewidth]{plots-final/productionew1.pdf}
  \caption{}
\end{subfigure}
\caption{Diagrams for the pair production of the leptonic scalar mediator, $\phi_l$, leading to $\mu \mu / \mu \nu + E_{\rm T} \hspace{-1.2em}/\;\;$ signatures in the final state.}\label{fig:prod-ew}
\end{figure}

\subsection{Simulation details}

We have implemented this model in \texttt{Feynrules~2.3}\cite{Alloul:2013bka}. The calculation of the matrix elements and the event generation is done using \texttt{MadGraph5\_aMC@NLO~2.6.3}\cite{Alwall:2014hca}. Production and decay of the new particles are considered independently using the narrow width approximation, as implemented in \texttt{MadSpin}\cite{Artoisenet:2012st}, which further accounts for spin correlations in decay chains\footnote{The narrow width approximation is not valid in benchmark points B1 and B2, for which we have taken interference effects into account. }. We then use \texttt{Pythia 8.235}\cite{Sjostrand:2014zea} to shower the parton-level events and we pass the output to \texttt{CheckMATE 2.0.26} \cite{Dercks:2016npn}, which compares the expected signal with supersymmetric searches at the LHC and derives an exclusion limit. As we have explained above, we can apply squark and slepton searches to constrain the coloured and leptonic mediator, respectively.

In order to describe initial and final state radiation and reproduce the correct jet structure precisely, we consider leading order (LO) production with parton shower matching and multijet merging when needed. The LO multijet merging techniques describe how parton shower emissions can be combined with full matrix element calculations to achieve a better accuracy in the description of the radiation spectrum. Using this technique, every jet is classified according to its $p_T$ and then compared to a hardness scale $Q_{\text cut}$. In this way, emissions above the hardness scale $Q_{\text cut}$ are described at LO accuracy using the corresponding matrix element calculation for an extra hard, wide-angle QCD emission in the final state, while emissions below this scale are defined as soft or collinear jets and the all-orders resummation description from the parton shower is preserved. Note that even though $\pazocal{O}(\alpha_s)$ corrections are included using this procedure, the calculation remains formally $\text{LO}+\text{LL}$ accurate after parton shower due to missing virtual corrections.

After hadronization, the showered events and the production cross sections are passed to \texttt{CheckMATE}. Each model point is tested against all the implemented experimental analyses to determine the optimal signal region. For this signal region, \texttt{CheckMATE} compares the simulated signal with the actual experimental observation and determines whether the model point is excluded at the $90\%$ confidence level.

\subsection{Results}

Constraints from LHC searches for the four benchmark points defined in Section~\ref{sec:benchmark} are presented in Figure~\ref{fig:lhc} on the $(m_{\chi},\,m_{\phi_q})$ plane, for all the points that satisfy the flavour constraints of Section~\ref{sec:bsmixing} and that reproduce the correct DM relic abundance. This figure shows the complementarity between the experimental limits obtained from the $pp \to jj + E_{\rm T} \hspace{-1.2em}/\;\;$ and $pp \to \mu\mu + E_{\rm T} \hspace{-1.2em}/\;\;$ searches. The experimental results used in our analysis are summarised in Table~\ref{tab:lhc_searches}. The colour code represents the average value of the coupling $|\lambda_b|$ in the region allowed by flavour constraints, defined as $|\lambda_b|_{\rm mean} = (|\lambda_b|_{\rm max} + |\lambda_b|_{\rm min}) / 2$, where $|\lambda_b|_{\rm max}$ and $|\lambda_b|_{\rm min}$ are the maximum and minimum allowed values respectively. The variation of our results when choosing either the minimum or maximum value for $|\lambda_b|$ has been checked and is insignificant.

Regarding the $pp \to jj + E_{\rm T} \hspace{-1.2em}/\;\;$ search, the limits in every scenario show that for the lightest DM mass, coloured mediators with masses below $\sim$1~TeV are excluded.  Even though heavier DM produces larger amounts of missing energy in final states, the cross section decreases rapidly with the $m_{\chi}$, leading to similar exclusion limits. It is interesting to note that exclusion limits are slightly stronger for the scenarios with $|\lambda_{b,t}| > |\lambda_{s,c}|$, where mediators with masses below $\sim$1.1~TeV are excluded. The reason for this is that final states with either top or bottom quarks are more sensitive to some experimental searches. The most stringent experimental search involves final states with at least two ($b \bar b$ production) or four ($t \bar t$ production) jets or exactly two leptons and missing energy \cite{Aaboud:2017rzf}. In particular, the most sensitive signal region is optimised to detect events featuring a DM particle produced in association with a $t \bar t$ pair, which decays fully hadronically.

\begin{figure}[t!]
\centering
  \includegraphics[width=0.49\linewidth]{plots-final/testcollider_mpsi_gmu2_x1_uniform.png}
  \includegraphics[width=0.49\linewidth]{plots-final/testcollider_mpsi_gmu4_x1_uniform.png}\\[2ex]
  \includegraphics[width=0.49\linewidth]{plots-final/testcollider_mpsi_gmu2_x4_uniform.png}
  \includegraphics[width=0.49\linewidth]{plots-final/testcollider_mpsi_gmu4_x4_uniform.png}
\caption{LHC limits from the $pp \to jj + E_{\rm T} \hspace{-1.2em}/\;\;$ (green) and $pp \to \mu\mu + E_{\rm T} \hspace{-1.2em}/\;\;$ (blue) searches. On the left (right), results for the scenario with $|\lambda_\mu| = 2$ ($|\lambda_\mu| = \sqrt{4\pi}$) are presented. The figures in the upper panel correspond to $\lambda_b = \lambda_s^*$, while the bottom panel shows limits for $\lambda_b = 4\lambda_s^*$. The colour code represents the average value of the coupling $|\lambda_b|$ in the region allowed by flavour constraints, as defined in the text. Solid lines represent the current exclusion limits, whereas dashed ones correspond to the projected reach of the LHC High luminosity phase. }\label{fig:lhc}
\end{figure}

\vspace{10pt}
\begin{table}[t!]
\begin{center}
\begin{tabular}{|c|l|l|c|}
\hline
$\sqrt{s}$ & Search & Final state & $\pazocal{L}$ [fb$^{-1}$] \\
\hline
\hline
13 TeV & 1710.11412 \cite{Aaboud:2017rzf} & $t\bar t \, / \, b\bar b + E_{\rm T} \hspace{-1.2em}/\;\;$ & 36.1 \\
 & ATLAS-CONF-2017-039 \cite{ATLAS:2017uun} & $2l \, / \, 3l + E_{\rm T} \hspace{-1.2em}/\;\;$ & 36.1 \\
\hline
\hline
14 TeV & atlas\_phys\_pub\_013\_011 \cite{ATL-PHYS-PUB-2013-011} & $t\bar t + E_{\rm T} \hspace{-1.2em}/\;\;$ & 3000 \\
& atlas\_phys\_2014\_010\_hl\_3l \cite{ATL-PHYS-PUB-2014-010} & $2l \, / \, 3l + E_{\rm T} \hspace{-1.2em}/\;\;$ & 3000 \\
\hline
\end{tabular}
\caption{List of experimental searches sensitive to our model, where $l$ denotes electron and muon. The third column describes the final state targeted by the analysis and the last column displays the total integrated luminosity. }
\label{tab:lhc_searches}
\end{center}
\end{table}
\vspace{-15pt}

Regarding the $pp \to \mu\mu + E_{\rm T} \hspace{-1.2em}/\;\;$ search, the limits show that models with dark matter masses below approximately $ 30 \; {\mathrm{GeV}}$ are ruled out for $|\lambda_\mu| = 2$, with the exclusion limit going down to $\sim 13 \; {\mathrm{GeV}}$ for $|\lambda_\mu| = \sqrt{4\pi}$. This corresponds to mediator masses below $360 \; {\mathrm{GeV}}$ for $|\lambda_\mu| = 2$ and $410 \; {\mathrm{GeV}}$ for $|\lambda_\mu| = \sqrt{4\pi}$. The $pp \to \mu\mu + E_{\rm T} \hspace{-1.2em}/\;\;$ cross section mainly depends on $m_{\phi_l}$, so the limits on $m_\chi$  can be understood through its relation with  $m_{\phi_l}$ given by the DM relic condition  (\ref{eq:sigmav}) for a particular value of $\lambda_\mu$. The most stringent search involves final states with $2l + 0 j$, $2l$ and at least $2$ jets, or $3l$ and missing energy \cite{ATLAS:2017uun}. In particular, the most sensitive signal region is characterised by $2l + 0j$ and a dilepton invariant mass $m_{ll} > 300 \; {\mathrm{GeV}}$, and it is optimised to target slepton pair production.

The most remarkable result is that LHC limits completely exclude the scenario with $|\lambda_\mu| = 2$ and $\lambda_b = \lambda_s^*$, as well as a sizeable region of the scenario with $\lambda_b = 4\lambda_s^*$ for the same $|\lambda_\mu|$. These constraints become weaker for larger values of $|\lambda_\mu|$ and, for the scenarios with $|\lambda_\mu| = \sqrt{4\pi}$, most of the parameter space is allowed. It is crucial to note that the limits coming from final states with jets and leptons are complementary to each other. While the former exclude regions of the parameter space with large $m_\chi$ and small $m_{\phi_q}$, the latter rule out models with very heavy mediator masses $m_{\phi_q}$ and light dark matter. Importantly, these limits are also complementary to the ones coming from direct detection, where dark matter masses below $12 \; {\mathrm{GeV}}$ lie below the neutrino floor. Therefore, it is fundamental to consider both approaches to explore the model.

It is worth mentioning that the small couplings required by flavour constraints lead to decay widths slightly below the QCD scale for $m_{\phi_q} \lesssim 370 \; {\mathrm{GeV}}$. Strictly speaking, this means that the computation of the decay width cannot be handled perturbatively and that the new particle $\phi_q$ may hadronize into bound states with SM quarks, analogous to R-hadrons \cite{FARRAR1978575}, before decaying. However, the typical width involved is $\Gamma_{\phi_q} \sim \pazocal{O}(10^{-2})-\pazocal{O}(10^{-3}) \; {\mathrm{GeV}}$, which means lifetimes of the order $\tau \sim 10^{-22} \; \text{s}$, so any potential bound state would decay promptly in the detector. This region of the parameter space is excluded by ATLAS and CMS R-hadron searches \cite{Aad:2013gva,CMS-PAS-EXO-16-004}.

We have also studied the limits that could be obtained with $3000 \; \text{fb}^{-1}$ of 14~TeV data once the LHC High Luminosity phase \cite{Apollinari:2017cqg} is completed. As we can observe in the plots, the main gain would come from the leptonic channels, which would allow to test a considerable amount of the model's parameter space. In particular, scenarios with $|\lambda_\mu| < 2$ would be completely excluded. The experimental searches giving the strongest exclusion limits target the same final states and are shown in the low panel of Table \ref{tab:lhc_searches}.

\section{Direct DM detection prospects}
\label{sec:direct}

Finally, in this section we discuss whether our model is expected to produce an observable response in direct detection experiments. We have calculated this response, by matching the model parameters to effective DM-nucleon interaction terms,
\begin{equation}
{\cal L}_{\text{int}}=\sum_{N} \sum_i c_i^{N}\mathcal{O}_i \overline{\chi} \chi \overline{N} N\ ,
\label{eq:eft}
\end{equation}
where $N$ is the corresponding nucleon, and $\mathcal{O}_i$ is the set of non-relativistic operators \cite{Fan:2010gt,Fitzpatrick:2012ix}. The values for the coefficients $c_i^N$ can be derived as the non-relativistic limit of the original interaction Lagrangian, and the differential rate can be computed using the corresponding nuclear form factors from Refs.~\cite{Fitzpatrick:2012ix,Anand:2013yka}, and for a given choice of the DM halo properties. We  have adopted the so-called standard halo model~\cite{Drukier:1986tm} with local DM density $\rho_{\chi}=0.4\ \textrm{GeV/cm}^3$, a central velocity of $v_0=220\ \textrm{km s}^{-1}$, and a escape speed of $v_{\textrm{esc}}=544\ \textrm{km s}^{-1}$ to calculate the number of expected recoils in a specific experiment.

The leading tree-level DM-quark interactions are given by scalar ($\overline{\chi}\chi\overline{\psi}\psi$) and vector ($\overline{\chi}\gamma^{\mu}\chi\overline{\psi}\gamma_{\mu}\psi$) type interactions. The latter is the leading contribution to $\mathcal{O}_1$ for Dirac DM \cite{Ibarra:2015fqa}, but it vanishes in the case of Majorana DM. For scalar type interactions Majorana DM does not in general vanish, but with our models chiral structure, it does.  With sub-dominant couplings to the first generation of quarks, and given that $m_{\phi_q}>m_{\phi_l}$, one-loop contributions to the DM-nucleon scattering cross sections will generally be larger than the tree level process. The loop contributions for a generic fermionic DM that involve the exchange of a photon can be classified as electric and magnetic dipoles ($\overline{\chi}i\sigma^{\mu\nu}\gamma^5\chi F_{\mu\nu}$ and $\overline{\chi}\sigma^{\mu\nu}\chi F_{\mu\nu}$, respectively), anapole ($\overline{\chi}\gamma^{\mu}\gamma^5\chi\partial^{\nu}F_{\mu\nu}$), and charge radius ($\overline{\chi}\gamma^{\mu}\chi\partial^{\nu}F_{\mu\nu}$). However, in the particular case of Majorana DM considered in this work, the magnetic dipole and charge radius effective couplings are forbidden by charge conjugation symmetry. Thus, the dominant one-loop interaction is the anapole moment \cite{Ho:2012bg}. When taking the non-relativistic limit,  the anapole moment gives contributions to the $\mathcal{O}_8$ and $\mathcal{O}_9$ operators \cite{Gresham:2014vja,DelNobile:2018dfg}, which are velocity and momentum dependent. In terms of the fundamental parameters of the model, the corresponding couplings read
\begin{eqnarray}
    c_8&=& 2 e \mathcal{A} Q_N\nonumber \, ,\\
    c_9&=& - e \mathcal{A} g_N \, ,
\end{eqnarray}
where $e$ is the electron charge, $Q_N$ is the nucleon charge, and $g_N$ are the nucleon g-factors ($g_p=5.59$ and $g_n=3.83$). The  effective coupling to the anapole interaction term, $\mathcal{A}$, reads \cite{Kopp:2014tsa}
\begin{eqnarray}
    \mathcal{A} = -\frac{e\,|\lambda_\mu|^2}{96\pi^2m_{\chi}^2}
    \left[
    \frac32 \log\frac{\mu}{\epsilon}
    - \frac{1+3\mu-3\epsilon}{\sqrt{(\mu-1-\epsilon)^2-4\epsilon}}
    \arctanh{\left(\frac{\sqrt{(\mu-1-\epsilon)^2-4\epsilon}}{\mu-1+\epsilon}\right)}
    \right]\ ,
\end{eqnarray}
with $\mu\equiv m_{\phi_l}^{2}/m_{\chi}^2$ and $\epsilon\equiv m_l^2/m_{\chi}^2$. The nuclear responses to the $\mathcal{O}_8$ and $\mathcal{O}_9$ operators are markedly weaker than that of $\mathcal{O}_1$, which implies that, in general, the scattering cross section is very small and beyond current experimental limits. Furthermore, because our DM particle interacts with the quark sector, it is not a priori clear that the spin-independent $\mathcal{O}_1$ and spin-dependent $\mathcal{O}_4$ arising from the so-called twist-2 operator \cite{Hisano:2010ct,Hisano:2015bma,jubbthesis} and the axial vector operator respectively are still negligible.

Given the range of DM masses that we consider in this study, the main constraint is due to Xenon1T results \cite{Aprile:2018dbl}, which we simulate using the prescription outlined in appendix A of Ref.~\cite{Kavanagh:2018xeh}, achieving good agreement. As we can see in Figure~\ref{fig:anapole}, the theoretical predictions for this model are beyond the reach of current experimental searches. We also show the reach of future direct detection experiments. The LZ detector, will employ 5.6 tons of liquid xenon with 1000 days exposure as outlined in \cite{Akerib:2015cja,Mount:2017qzi}. The DarkSide-20k experiment\cite{Aalseth:2017fik}, is an argon detector which will employ $20$ tons of fiducial mass for a duration of 10 years. We have assumed that the DarkSide collaboration will be able to achieve a threshold energy of $5$ keV, a reasonable assumption considering the results from DarkSide-50 \cite{Agnes:2018ves}. For reference we have also calculated the neutrino floor for anapole interactions in the $(\mathcal{A},\, m_{\chi})$ plane. We have used the prescription described in Ref.\,\cite{Billard:2013qya} and the expected neutrino fluxes from Refs.~\cite{Gaisser:2002jj,Battistoni:2005pd,Horiuchi:2008jz,Serenelli:2011py,2013arXiv1301.0365H,Gelmini:2018ogy}. It is clear that our model favourably lays in a region of parameter space that would be probed by a generation of experiments with multi-ton targets, that can probe near or even slightly beyond the neutrino floor. Spectral analysis with the neutrino background compounded with annual modulation data, could provide complete discrimination between model and the anapole moment which is both velocity and momentum dependent.

\begin{figure}[t!]
    \centering
    \hspace*{-0.5cm}\includegraphics[width=0.53\textwidth]{plots-final/anapole_plot_noLHC.png}\hspace*{-0.3cm}
    \includegraphics[width=0.53\textwidth]{plots-final/anapole_plot_LHC.png}
    \caption{Theoretical predictions for the anapole coupling, $\mathcal{A}$, as a function of the DM mass, $m_{\chi}$ for the four benchmark points: A1 (red points), A2 (green), B1 (orange), and B2 (blue). For comparison, we show the current exclusion line by Xenon1T \cite{Aprile:2018dbl} and the predicted reach of LZ \cite{Akerib:2015cja,Mount:2017qzi} and DarkSide-20k \cite{Aalseth:2017fik}. The shaded area represents the neutrino floor. The plot on the right-hand side incorporates LHC constraints, explained in more detail in Section\,\ref{sec:lhc}.}
    \label{fig:anapole}
\end{figure}

\begin{figure}[t!]
    \centering
    \hspace*{-0.3cm}\includegraphics[width=0.53\textwidth]{plots-final/c1_noLHC4.png}\hspace*{-0.5cm}
    \includegraphics[width=0.53\textwidth]{plots-final/c1_LHC_4.png}
    \caption{The same as in Figure~\ref{fig:anapole}, but for the spin-independent coupling, $c_1$, that originates from the twist-2 coupling.
    }
    \label{fig:twist2}
\end{figure}

For completeness, we have also calculated the effect on the total scattering cross section from aforementioned twist-2 operator and spin-dependent interaction. The former contribute to the spin-independent scattering cross section (operator ${\cal{O}}_1$) and can be sizeable if the new coupling to quarks is large or the colour mediator is very light. We have explicitly checked that once LHC constraints are included in the parameter space of the model, these terms are always subdominant to the anapole term discussed above. We represent in Figure~\ref{fig:twist2} the theoretical predictions for $c_1$ as a function of the DM mass from this contribution. For the spin-dependent interaction, we found that the predicted rate for our sampled parameter space is always sub-dominant.

Had we chosen to work with a Dirac fermion, the dipole and charge radius contributions should have been added. As it has been pointed out in Ref.~\cite{Ibarra:2015fqa}, the fairly large coupling to muons that is required to explain the flavour anomalies leads to effective DM couplings that are orders of magnitude higher than those coming from the tree level contribution, the most important being the charge-radius interaction. This we have checked, and in fact above $m_{\chi}\sim 10$ GeV, all our parameter points are excluded by Xenon1T. Below $m_{\chi}\sim 10$, the model is excluded by both LHC constraints and indirect detection bounds. Unlike in the Majorana case, the s-wave contribution to the thermal cross section $\langle\sigma v\rangle$ is no longer helicity suppressed and hence excluded \cite{Leane:2018kjk}.

Our results suggest that future multi-ton direct detection experiments, such as DarkSide \cite{Aalseth:2017fik}, would be able to probe this model in the mass range $m_{\chi}\sim 10-60$~GeV. It is very interesting to point out that many of the points in this DM mass range feature very heavy $\phi_q$ and therefore would be beyond the reach of collider searches. In a sense, future direct DM detection and the LHC complement each other to probe a large part of the model's parameter space.

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

In this article, we have studied a particle physics model that addresses the hints of lepton flavour universality violation observed by LHCb in $b \rightarrow s\mu^{+}\mu^{-}$ transitions, and that provides a solution to the dark matter problem. The scenario that we have analysed incorporates two new scalar fields and a Majorana fermion that provide one-loop contributions to $B$ meson decays.

The Majorana fermion is stable and can reproduce the observed DM relic abundance. We have studied the effect of new physics in flavour observables, for which $B_{s}-$mixing and $b \rightarrow s\mu^{+}\mu^{-}$ processes provide the most important constraints. In order to find an explanation for the $B$ anomalies and to reduce the $1.8\,\sigma$ tension between the predicted and measured mass difference in $B_{s}-$mixing, complex couplings are needed.  We have used results from the first global fit that takes into account this possibility. The combination of flavour bounds and constraints on the DM relic abundance leads to upper limits on the masses of the exotic states, and in general points towards a rather light DM candidate (with a mass $m_{\chi}\lesssim 200$~GeV).

We have studied the signatures that this model would produce at the LHC. The dominant processes are the pair production of the coloured and leptonic scalars. For the former, the strongest exclusion limits are given by $\text{dijet} + E_{\rm T} \hspace{-1.2em}/\;\;$ searches. For the latter, the final states are very clean, containing $1$ or $2$ leptons and missing energy. Both searches are complementary and exclude different regions of the parameter space, setting lower bounds on DM and mediator masses. The high-luminosity phase improves bounds coming from both searches, with dilepton being the most pronounced. The collider constraints are weakened when the $\lambda_{\mu}$ parameter is pushed towards the perturbative limit.

Finally, we have investigated how DM direct detection experiments constrain this model. Given the range of DM masses that we consider in this study, the main constraint is due to Xenon1T results. The small new couplings required by flavour constraints means that one-loop contributions to the DM-nucleon scattering cross section are generally larger than the tree level process. In particular, the dominant loop induced interaction is the anapole moment. We have shown that this model is not excluded by current data and could be probed by the next generation of experiments with multi-ton targets in the mass range $m_{\chi}\sim 10-60$~GeV.

\vspace*{2ex}

\noindent
{\Large{\bf Acknowledgements}}

We would like to thank Alejandro Ibarra for useful discussions. DC and AC are grateful for the support from the Science and Technology Facilities Council (STFC). DC acknowledges the partial support of the Centro de Excelencia Severo Ochoa Program through the IFT-UAM/CSIC Associate programme.
PMR and JMM acknowledge support from the Spanish Research Agency (Agencia Estatal de Investigaci\'on) through the contract FPA2016-78022-P and IFT Centro de  Excelencia Severo Ochoa under grant SEV-2016-0597.

\bibliographystyle{JHEP-cerdeno}


text:
\section{Introduction}

There are numerous applications in physics where an action
\eq{
I_{\rm bulk}[\phi]=\int_{\MM}\nts\nts d^nx \,\LL_{\rm bulk}(\phi,\nabla\phi)
}{eq:lp1}
has to be supplemented by boundary terms
\eq{
I_{\rm tot}[\phi]=I_{\rm bulk}[\phi]+\int_{\dM}\nts\nts\nts d^{n-1}x\,\LL_{\rm boundary}(\phi,\nabla_\bot\phi,\nabla_\parallel\phi)\,.
}{eq:lp2}
Here $\nabla_\bot$ and $\nabla_\parallel$ denote the normal and parallel components of the derivative with respect to the boundary $\dM$.
The simplest example is quantum mechanics, where
\eq{
I_{\rm bulk}[q,p] = \int_{t_0}^{t_1}\nts\nts dt\,\left[-q\dot{p}-H(q,p)\right]
}{eq:lp3}
has to be supplemented by a  boundary term
\eq{
I_{\rm tot}[q,p]=I_{\rm bulk}[q,p]+qp|_{t_0}^{t_1}
}{eq:lp4}
if Dirichlet boundary conditions are imposed on the coordinate, $\delta q|_{t_i}=0$. In addition to this ``Gibbons-Hawking-York'' boundary term one can add another boundary term
\eq{
\Gamma[q,p] = I_{\rm tot}[q,p] - \FF(t,q)|_{t_0}^{t_1}
}{eq:}
which depends only on quantities held fixed at the boundary. This seems to be a superfluous addition, as it does neither change the equations of motion nor the variational principle (as opposed to the ``Gibbons-Hawking-York'' boundary term), but in some applications such a term is crucial and determined almost uniquely from consistency requirements: symmetries and accessibility of the classical approximation. One such application is the Euclidean path integral for black holes (BHs), which provides a convenient shortcut to BH thermodynamics.

\section{Hamilton-Jacobi counterterm in two-dimensional gravity}\label{se:2}

The essential features and difficulties arise already in low dimensions. For transparency we focus on two-dimensional (2D) models. The bulk action for 2D dilaton gravity \cite{Grumiller:2002nm},
\begin{equation}\label{eq:bulkaction}
  I_{\rm bulk}[g,X] = - \frac{1}{16\pi G_2}\,\int_{\MM} \nts \nts d^{\,2}x \,\sqrt{g}\, \left[ X\,R - U(X)\,
		\left(\nabla X\right)^2 - 2 \, V(X) \raisebox{12pt}{~}\right]
\end{equation}
has to be supplemented by a Gibbons-Hawking-York boundary term
\eq{
I_{\rm tot}[g,X] = I_{\rm bulk}[g,X] -  \frac{1}{8\pi G_2}\, \int_{\dM} \bns dx \, \sqrt{\gamma}\,X\,K ~,
}{eq:GHYaction}
if Dirichlet boundary conditions are imposed on the dilaton field $X$ and the induced metric at the boundary $\gamma$. The meaning of all symbols is standard and our notation is consistent with \cite{Grumiller:2007ju}. In addition one could add another boundary term
\eq{
\Gamma[g,X]=I_{\rm tot}[g,X] - \underbrace{\frac{1}{8\pi G_2}\, \int_{\dM} \bns dx \, \sqrt{\gamma}\,\FF(X,\nabla_\parallel X,\gamma,\nabla_\parallel\gamma)}_{I_{\rm CT}[\gamma,X]}\,.
}{eq:angelinajolie}
We demonstrate now why such a term is needed and show that it is determined essentially uniquely from consistency requirements: symmetries and accessibility of the classical approximation.

\subsection{Symmetries}
\label{sec:symmetries}

Diffeomorphism covariance along the boundary requires that $\FF$ in \eqref{eq:angelinajolie} transforms as a scalar. Since there are no scalar invariants constructed from $\gamma$ in one dimension, $\mathcal F$ can be reduced to
$\FF=\FF(X,\nabla_\parallel X)$.
Another simplification arises if we restrict ourselves to isosurfaces of the dilaton field, which is sufficient for our purposes. Then $X$ is constant along the boundary, $\nabla_\parallel X=0$, so that we are left with a function
\eq{
\FF=\FF(X)\,.
}{eq:lp6}
Symmetry requirements have reduced the dependence on four variables in \eqref{eq:angelinajolie} to a dependence on only one variable, $X$.

\subsection{Accessibility of the classical approximation}
\label{sec:accessibility}
While symmetries help to reduce the ambiguities in $\FF$ they do not explain why such a term is needed in the first place. To this end we consider the Euclidean path integral,
\begin{equation}\label{PartitionFunction}
  \ZZ = \int \CurlD g \CurlD X \, \exp\Big(-\frac{1}{\hbar}\,I[g,X] \Big) ~.
\end{equation}
The path integral is evaluated by imposing boundary conditions on the fields and then performing the weighted sum over all relevant space-times $(\MM,g)$ and dilaton configurations $X$. In the classical limit it is dominated by contributions from stationary points of the action. This can be verified by expanding it around a classical solution
\begin{equation}\label{Saddle}
  I[g_{cl} + \delta g, X_{cl} + \delta X] = I[g_{cl},X_{cl}] + \de I[g_{cl}, X_{cl}; \delta g, \delta X] + \frac12 \de^2I[g_{cl}, X_{cl}; \delta g, \delta X] + \ldots
\end{equation}
where $\de I$ and $\de^2 I$ are the linear and quadratic terms in the Taylor expansion. The saddle point approximation of the path integral
\begin{equation}\label{ApproxPF}
 \ZZ \sim \exp\Big(-\frac{1}{\hbar}\,I[g_{cl},X_{cl}]\Big) \, \int \CurlD \delta g \CurlD \delta X \, \exp\Big(-\frac{1}{2\hbar}\,\de^2I[g_{cl}, X_{cl}; \delta g,\delta X] \Big)
\end{equation}
is defined if:
\begin{enumerate}
\item The on-shell action is bounded from below, $I[g_{cl},X_{cl}]>-\infty$.
\item The first variation vanishes on-shell, $\de I[g_{cl}, X_{cl}; \delta g, \delta X] =0$ for all variations $\delta g$ and $\delta X$ preserving the boundary conditions.
\item The second variation has the correct sign for convergence of the Gaussian in \eqref{ApproxPF}.
\end{enumerate}
The last condition actually means consistency of the semi-classical approximation, and we shall not discuss it here. Instead, we focus on the first two conditions. Both are violated for typical BH solutions of \eqref{eq:GHYaction} if the boundary is located in the asymptotic region $X\to\infty$:
\begin{enumerate}
\item The on-shell action behaves as $I_{\rm tot}[g_{cl},X_{cl}\to\infty]=2M/T-S-w(X_{cl}\to\infty)/T$, and $\lim_{X\to\infty} w(X)\to\infty$ for most models of interest.\footnote{$T$ is the Hawking temperature and $S$ the Bekenstein-Hawking entropy. Both are determined from the mass $M$ and the functions $Q(X)$ and $w(X)$ defined in the Appendix.}
\item The first variation of the action receives a boundary contribution \begin{equation}\label{ActionVariationIntro}
   \left.\delta I_{\rm tot} \raisebox{12pt}{\,}\right|_{\rm on-shell} \sim \int_{\dM} \bns dx \, \sqrt{\gamma}\,
	\left[ \pi^{ab} \, \delta \gamma_{ab} + \pi_{X} \, \delta X \, \right] \neq 0
\end{equation}
because the product $\pi^{ab} \, \delta \gamma_{ab}$ is non-vanishing: the variation of the induced metric $\gamma_{ab}$ does not fall off sufficiently fast to compensate for the divergence of the momenta $\pi^{ab}$.
\end{enumerate}
We emphasize that an ad-hoc subtraction $I_{\rm ren}[g_{cl},X_{cl}\to\infty]:=I_{\rm tot}[g_{cl},X_{cl}\to\infty]+w(X_{cl}\to\infty)/T$ is inconsistent: while it leads to a finite on-shell action it does not address the second problem. Both can be solved by choosing $\FF(X)$ adequately.

Since the on-shell action solves the Hamilton-Jacobi equation one can expect cancellations if also $\FF$ is a solution to the Hamilton-Jacobi equation. Thus, our guiding principle is to demand that $\FF$ be a solution of the Hamilton-Jacobi equation (see the next Section for details).
This method was applied first to the Witten BH and to type 0A string theory \cite{Davis:2004xi} and later generalized to generic 2D dilaton gravity \cite{Grumiller:2007ju}. The result is
\begin{equation}\label{LCT}
  \FF(X) =  - \sqrt{\left(w(X) + c\right)e^{-Q(X)}}
\end{equation}
where $c$ is an integration constant. It can be absorbed into a redefinition of $w$ (cf.~Appendix) and reflects the freedom to choose the ground state of the system. Thus, without loss of generality we can set it to zero and finally obtain a consistent action \cite{Grumiller:2007ju}
\eq{
\begin{split}
  \Gamma[g,X] &= - \frac{1}{16\pi G_2}\,\int_{\MM} \nts \nts d^{\,2}x \,\sqrt{g}\, \left[ X\,R - U(X)\,
		\left(\nabla X\right)^2 - 2 \, V(X) \right] \\
&\quad - \frac{1}{8\pi G_2}\, \int_{\dM} \bns dx \, \sqrt{\gamma}\,X\,K
+ \frac{1}{8\pi G_2}\, \int_{\dM} \bns dx \,\sqrt{\gamma}  \,
	\sqrt{w(X) \, e^{-Q(X)}} \,.
\end{split}
}{ActionConclusion}
The classical approximation is now well-defined because
\begin{enumerate}
\item $\Gamma[g_{cl},X_{cl}\to\infty]=M/T-S$ is finite.
\item $\delta\Gamma[g_{cl},X_{cl};\delta g,\delta X]=0$ for all $\delta g$ and $\delta X$ preserving the boundary conditions.
\end{enumerate}
Moreover, as opposed to \eqref{eq:GHYaction} the action \eqref{ActionConclusion} is consistent with the first law of thermodynamics. Perhaps the most remarkable property of the counterterm in \eqref{ActionConclusion} is its universality: while usually different subtraction methods are employed depending on whether spacetime is asymptotically flat, $AdS$ or neither of both, our result is not sensitive to the asymptotics. This universality does not appear to exist in higher dimensions or even in 2D if `standard' subtraction methods are used \cite{Regge:1974zd}.

\section{Cartan formulation}
In many applications a first order formulation in terms of Cartan variables is advantageous \cite{Grumiller:2002nm}. Therefore we now derive the Hamilton-Jacobi counterterm in this formulation. The corresponding action (we set $8\pi G_2=1$)
\begin{equation}
\label{FOaction}
 \begin{split}
 I^{\rm FO}_{\rm tot}[X, Y^a, e_a, \omega] &= - \int_{\mathcal M} \left[ Y^a D e_a + Xd\omega + \epsilon\,(\frac12 U(X)Y^aY_a+V(X)) \right] \\ &\quad + \int_{\partial \mathcal M} \left[ X\omega_{\|} + \frac{i}{2} Xd\ln \frac{e_{\|}}{\bar{e}_{\|}}\right]\ ,
\end{split}
\end{equation}
contains the Cartan variables $\omega$ and $e_a$, as well as the scalar fields $X$ and $Y^a$ (we use a complexified dyad, $\bar{e}=e^\ast$, cf.~\cite{Bergamin:2004pn} for the details of our notation). As \eqref{FOaction} is classically equivalent to \eqref{eq:GHYaction} (cf.~the Appendix of Ref.~\cite{Bergamin:2004pn}) it will suffer from the same problems as described in Section \ref{sec:accessibility}. We follow the same strategy as in the second order formulation \cite{Davis:2004xi,Grumiller:2007ju} to find the corresponding Hamilton-Jacobi counterterm $I_{\rm CT}^{\rm FO}$, which, using the arguments in Section \ref{sec:symmetries}, can be reduced to
\begin{equation}
\label{FOCT}
\Gamma^{\rm FO}[X, Y^a, e_a, \omega] =  I^{\rm FO}_{\rm tot}[X, Y^a, e_a, \omega] - \underbrace{\int_{\dM} \bns dx^{\|} \, \sqrt{2e_{\|} \bar e_{\|}}\,\FF(X)}_{I_{\rm CT}^{\rm FO}[X,e_{\|}\bar{e}_{\|}]}\, .
\end{equation}
The variation of the action produces the equations of motion plus the boundary term\footnote{Contributions emerging from the logarithm in the boundary term are dropped as we assume that the boundary is an isosurface of the dilaton and that there is no boundary of the boundary. The coordinate along the boundary, $x^{\|}$, can be thought of as Euclidean time $\tau$.}
\begin{equation}
 \int_{\partial \mathcal M} \bns d x^{\|} \,\left[ Y \delta \bar{e}_{\|} + \bar{Y} \delta \bar e_{\|} - \delta X \omega_{\|}\right]\, .
\end{equation}
To cancel it we assign Dirichlet boundary conditions to $X$, $e_{\|}$ and $\bar{e}_{\|}$. As in \cite{Davis:2004xi,Grumiller:2007ju} it is possible to write the momenta which are not fixed at the boundary,
\begin{align}
\label{OnShellVar}
 \omega_{\|} &= - \left.\frac{\delta I^{\rm FO}_{\rm tot}}{\delta X}\right|_{\rm on-shell} \,, & Y &= \left.\frac{\delta I^{\rm FO}_{\rm tot}}{\delta \bar{e}_{\|}}\right|_{\rm on-shell} \,,  & \bar{Y} &= \left.\frac{\delta I^{\rm FO}_{\rm tot}}{\delta e_{\|}} \right|_{\rm on-shell} \,,
\end{align}
as variations of the on-shell action. The Hamilton constraint
\eq{
- \omega_{\|} \frac{Y\bar{e}_{\|}+\bar{Y}e_{\|}}{2e_{\|}\bar{e}_{\|}} + U(X) Y\bar{Y}+ V(X) = 0\,,
}{eq:hc}
follows from a standard constraint analysis of the first order action \eqref{FOaction}. By construction the Hamilton-Jacobi counterterm must be a solution of this constraint. Replacing in \eqref{OnShellVar} the on-shell action $I_{\rm tot}^{\rm FO}$ by the counterterm $I_{\rm CT}^{\rm FO}$, plugging this into the Hamilton constraint \eqref{eq:hc} and exploiting that $I_{\rm CT}^{\rm FO}$ depends solely on the combination $e_{\|}\bar{e}_{\|}$ establishes
\eq{
\frac{\delta I_{\rm CT}^{\rm FO}}{\delta X} \frac{\delta I_{\rm CT}^{\rm FO}}{\delta(e_{\|}\bar{e}_{\|})} + U(X) e_{\|}\bar{e}_{\|} \left(\frac{\delta I_{\rm CT}^{\rm FO}}{\delta(e_{\|}\bar{e}_{\|})}\right)^2 + V(X) = 0\,.
}{eq:FDEFOWTF}
This functional differential equation for the counterterm by virtue of the Ansatz \eqref{FOCT} simplifies to
\eq{
\frac{d}{dX} \FF^2(X) + U(X) \FF^2(X) + 2V(X) = 0\,.
}{eq:ODEFFO}
The solution of this first order ordinary differential equation is given by\footnote{There is a sign ambiguity since \eqref{eq:ODEFFO} yields only $\FF^2(X)$. The sign choice in \eqref{eq:FFO} gives an action with a consistent classical limit.}
\eq{
\FF(X) =  - \sqrt{\left(w(X) + c\right)e^{-Q(X)}}\,.
}{eq:FFO}
This coincides with \eqref{LCT} and thus we conclude that the Hamilton-Jacobi counterterms in second and first order formalisms are identical, as might have been anticipated on general grounds. Setting again $c=0$, the consistent first order action is
\begin{equation}
\label{GammaFO}
 \begin{split}
 & \Gamma^{\rm FO}[X, Y^a, e_a, \omega] = - \int_{\mathcal M} \left[ Y^a D e_a + Xd\omega + \epsilon\,(\frac12 U(X)Y^aY_a+V(X)) \right] \\ &\quad + \int_{\partial \mathcal M} \left[ X\omega_{\|} + \frac{i}{2} Xd\ln \frac{e_{\|}}{\bar{e}_{\|}}\right] + \int_{\dM} \bns dx^{\|} \, \sqrt{2e_{\|} \bar e_{\|}}\,\sqrt{w(X)e^{-Q(X)}}\, .
\end{split}
\end{equation}

\section{Black hole thermodynamics and further applications}

An immediate consequence of our result \eqref{ActionConclusion} is the Helmholtz free energy \cite{Grumiller:2007ju}
\begin{equation}\label{FreeEnergy}
  F_c(T_c,X_c) = \frac{1}{8\pi G_2}\,\sqrt{w(X_c) e^{-Q(X_c)}}\left(1-\sqrt{1-\frac{2M}{w(X_c)}}\right)  - \frac{X_h}{4 G_2} \, T_c \,,
\end{equation}
which is related to the on-shell action in the usual way, $F_c= T_c \, \Gamma_c$. Here $X_c$ denotes the value of the dilaton field at the location of a cavity wall in contact with a thermal reservoir, while $X_h$ denotes the value of the dilaton at the BH horizon. The local temperature $T_c$ is related to the Hawking temperature $T$ by the standard Tolman factor, $T_c=T/\sqrt{\xi(X)}$. All other quantities are defined in the Appendix. The entropy,
\eq{
S=-\left.\frac{\partial F_c}{\partial T_c}\right|_{X_c} = \frac{X_h}{4 G_2} = \frac{A_h}{4 G_{\rm eff}}\,,
}{eq:S}
is in agreement with the Bekenstein-Hawking result. Here $A_h=1$ because we are in 2D, and $G_{\rm eff}=G_2/X_h$. For dimensionally reduced models \eqref{eq:S} can be interpreted also from a higher-dimensional perspective: $A_h\propto X_h$ and $G_{\rm eff}\propto G_2$, with the same proportionality constants. The result \eqref{eq:S} is well-known and was obtained by various methods \cite{Gibbons:1992rh}.
However, the free energy \eqref{FreeEnergy} contains a lot of additional information and allows a quasi-local treatment of BH thermodynamics (where applicable in agreement with \cite{Brown:1991gb}),
including stability considerations. For an extensive study of thermodynamical properties and more Refs.~we refer to \cite{Grumiller:2007ju}.

The class of BHs described by the action \eqref{ActionConclusion} or \eqref{GammaFO} is surprisingly rich (cf.~e.g.~table 1 in \cite{Grumiller:2006rc}), and includes spherically symmetric BHs (like Schwarzschild or Schwarzschild-$AdS$) in any dimension, spinning BHs in three dimensions \cite{Banados:1992wn}
and string BHs in two dimensions \cite{Elitzur:1991cb,
Dijkgraaf:1992ba}. As an example we consider now the exact string BH \cite{Dijkgraaf:1992ba}, and review some of its properties. Its target-space action \cite{Grumiller:2005sq} is given by \eqref{ActionConclusion} with $8\pi G_2=1$ and the potentials
 \eq{
U(X) = -\frac{\rho}{\rho^2+2(1+\sqrt{1+\rho^2})}\,,\qquad V(X) = -2 b^2 \,\rho\,.
}{eq:EBH1}
Here the canonical dilaton $X$ is related to a new field $\rho$ by
\eq{
X=\rho+{\rm arcsinh}\,{\rho} \,,
}{eq:EBH2}
and the parameter $b$ is related to the level $k$ and $\alpha'$ by
$\alpha'\,b^{2} = 1/(k-2)$.
In order for the background following from \eqref{ActionConclusion}, \eqref{eq:EBH1} and \eqref{eq:EBH2} to be a solution of string theory it must satisfy the condition $D - 26 + 6\,\alpha'\,b^2 = 0$. Because the target space here is two-dimensional, $D=2$, requiring the correct central charge fixes the level at the critical value $k_{\rm crit} = 9/4$. Following \cite{Kazakov:2001pj}, we vary $k$ by allowing for additional matter fields that contribute to the total central charge, so that $k \in [2,\infty)$ is possible. The Witten BH arises in the limit $k\to\infty$. Since it is not possible to place an abrupt cut-off on the space-time fields in string theory we have to consider the limit $X_c\to\infty$ in the Helmholtz free energy \eqref{FreeEnergy},
\eq{
F^{ESBH}=-b\,\sqrt{1-\frac{2}{k}}\,{\rm arcsinh}\,{\sqrt{k(k-2)}}
}{eq:EBH10}
and its thermodynamical descendants. It is straightforward to show that \eqref{eq:EBH10} leads to a positive specific heat for any $k\in(2,\infty)$, and that it vanishes in the limit $k\to 2$ in accordance with the third law.

We comment now briefly on the inclusion of semi-classical corrections from fluctuations of massless matter fields on a given BH background. We therefore add to the classical action \eqref{eq:bulkaction} the Polyakov action,
\eq{
I_{\rm bulk}^{\rm semi} = I_{\rm bulk} + c \int_{\MM}\bns d^2x\sqrt{g}\left[\psi R + \frac12 (\nabla\psi)^2\right]\,,
}{eq:Pol}
where we have introduced an auxiliary field $\psi$ which fulfills the on-shell relation $\square\psi=R$, and the constant $c$ depends on the number and type of massless matter fields.  Obviously, the addition of \eqref{eq:Pol} requires a reconsideration of boundary issues. One possibility is to demand that $\psi$ is a function of $X$ \cite{Zaslavsky:1998ca}.
Then the action \eqref{eq:Pol} reduces to a standard dilaton gravity action (we set $8\pi G_2=1$)
\eq{
I_{\rm bulk}^{\rm semi} = -\frac{1}{2} \int_{\MM}\bns d^2x\sqrt{g}\left[\hat{X}R-(U(X)+c (\psi'(X))^2)(\nabla X)^2-2V(X)\right]\,,
}{eq:redefaction}
upon introducing a redefined dilaton $\hat{X}=X-2\,c\,\psi(X)$. Therefore, as long as the assumption $\psi=\psi(X)$ is meaningful, the discussion of boundary terms in Section \ref{se:2} is still valid.  We note in this context that the large $X$ expansion of the exact string BH \eqref{eq:EBH1}, \eqref{eq:EBH2} can be interpreted as a semi-classical correction to the Witten BH, with $\rho$ playing the role of the unperturbed dilaton $X$ and $-\ln{(2\rho)/(2\,c)}$ playing the role of the auxiliary field $\psi$. This is consistent with the fact that the conformal factor $\psi$ scales logarithmically with the dilaton and concurs with semi-classical corrections \cite{Grumiller:2003mc} to the specific heat of the Witten BH, which also show qualitative agreement with the specific heat of the exact string BH. Another, more general, possibility is to treat $\psi$ as an independent field. In that case boundary issues have to be reconsidered. We expect  them to be relevant whenever the boundary term
\eq{
\int_{\dM}\bns dx\sqrt{\gamma} \,\psi n^\mu\partial_\mu\psi
}{eq:scbt}
does not fall off sufficiently fast at the asymptotic boundary ($n^\mu$ is the outward pointing unit normal). Since $\psi$ typically scales logarithmically with $X$ this happens for all models where $w(X)$ grows linearly or faster than $X$. Interestingly, the Witten BH is precisely the limiting case where this issue is of relevance. More recently semi-classical corrections were considered in the context of large $AdS$ BHs \cite{Hemming:2007yq}. There the issue is complicated because the matter fields couple non-minimally to the dilaton. Since the results of \cite{Hemming:2007yq} agree with ours in the large $X$ limit only,\footnote{The counterterm, (4.5) in \cite{Hemming:2007yq}, for finite values of $r$ differs from \eqref{LCT}, which yields $\FF(r^2)\sim r\sqrt{1+r^2/\ell^2}=1/\ell\,(r^2 + \ell^2/2 + \dots)$.} it would be interesting to analyze the counterterms using the Hamilton-Jacobi method, possibly by adapting the strategy described and applied in \cite{deBoer:1999xf}.
This would also allow to reconsider the path integral quantization of 2D dilaton gravity with matter in the presence of boundaries \cite{Kummer:1998zs}
and to clarify the role of the Hamilton-Jacobi counterterm for observables beyond thermodynamical ones.

For further applications and an outlook to future research we refer to the discussion in Section 7 of \cite{Grumiller:2007ju}.

\section*{Acknowledgments}

The content of this proceedings contribution was presented by one of us at the conference QFEXT07 in Leipzig, and DG would like to thank Michael Bordag as well as Boris, Irina, Mischa and Sascha Dobruskin for the kind hospitality.

The work of DG is supported in part by funds provided by the U.S. Department of Energy (DoE) under the cooperative research agreement DEFG02-05ER41360. DG has been supported by the Marie Curie Fellowship MC-OIF 021421 of the European Commission under the Sixth EU Framework Programme for Research and Technological Development (FP6).
The research of RM was supported by DoE through grant DE-FG02-91ER
40688 Task A (Brown University) and Perimeter Institute for
Theoretical Physics. Research at Perimeter Institute is supported
through the government of Canada through Industry Canada and by the
province of Ontario through the Ministry of Research and Innovation.

redpajama_set_name: RedPajamaArXiv


text:
\section{Introduction}
\label{sec:intro}

Observations of the Galactic Center (GC) going back as far as 20 years
\citep[see][ for a review]{gen10} reveal three facts: (1) An isotropic cusp of
young O/B and Wolf-Rayet (WR) stars, starting at a distance of $30\arcsec$ from
SgrA$^*$ and extending inward to about $1\arcsec$ ($1\arcsec\simeq0.04$ pc);
(2) a mildly thick stellar {\em disk}, of about $100$ WR and O-type stars,
spanning from an inner radius of $1\arcsec$ to an outer radius of about
$10\arcsec$; and (3) populating the innermost region, within $1\arcsec$ from
SgrA$^*$, are a population of B stars, commonly referred to as the ``S-stars'',
but {\em no} WR/O stars.  A single star-formation (SF)
episode may have explained the formation of disk and cusp stars
\citep[][]{lu13}, the S-stars, however, cannot have been born in this scenario,
because the violent environmental conditions within $1\arcsec$ do not allow
in-situ SF.  One way of populating that region is by dynamical friction, but
the associated timescale is too long. For this reason, the problem has the
reputation of the ``paradox of youth'' \citep{morris93,ghez03}.

This issue has led to the idea that S-stars could have formed at larger radii
and brought in later by an efficient dynamical mechanism. One possibility is
tidal separation of binaries \citep{Hills1991,gould03,gin06}.  A binary, formed
at larger radius, can be set in such an orbit that at periapsis it will be
tidally separated by the central massive black hole (MBH), leaving one star,
which could be a B star, bound to the MBH at a typical radius of $\la
1\arcsec$.  However, the captured stars would have very high eccentricities,
typically about $0.93-0.99$ (see the original work of \citealt{Hills1991} and
\citealt{ama12} for a review).  It would require some $20-50$ Myr for them to
achieve the observed near-thermal distribution
\citep[e.g.][]{perets09,antonini13,zhang13}, with the aid of a very dense cusp
of segregated old stars which is in contradiction with current observations
\citep{buc09,do09}.  Even in the presence of a dense cusp, there are two
additional issues: (1) the same process would work for WR/O stars, and we do
{\em not} see them within $\la 1\arcsec$ \citep{Tal11} and (2) the oldest O/WR
stars are $\la 10$ Myr, so there has to be at least two SF episodes, since
S-stars inside $1\arcsec$ need $\ge 20$ Myr to thermalize.

Since the stellar disk initially must be gaseous, it has been proposed as
another possibility that B stars migrated in it towards the center
\citep{levin07,griv10}. However, (1) this cannot explain the eccentricities of
the S-stars because the migrating stars will remain on near-circular orbits
\citep{perets09,madigan11,antonini13}, and (2) WR/O stars would have migrated
towards the center due to the same mecahnism,  but we do not observe them
there.  Once SF is over (no more gas), the stars in the disk, including the
WR/O stars, would secularly torque each other and drift away from
nearly-circular orbits to rather eccentric ones \citep{madigan09}, and hence at
periapsis they would populate the central $1\arcsec$, but, still, we do not see
WR/O-stars there either.

In this Letter we show that, provided the disk was heavier and more extended in
the past \citep{nay07,wardle08,Bon08,hobbs09,alig11,map12}, it created a rapid
evolving region (RER) inside $1\,\arcsec$, where the angular momenta of
stars rapidly redistribute because of a Kozai-Lidov-like resonance. This RER
can explain both the eccentricities of S-stars  and the absence of WR/O-stars
because the latter are tidally disrupted.

\section{Disk-driven evolution}
\label{sec:model}

\subsection{Timescales}

To understand the effect of the disk, we first analyze the
torque exerted by a wire of mass $\delta m$ and radius $R$ on a background star
of semi-major axis $a$ \citep{iva05,subr05,lockmann08}. \cite{chen11} showed
that the timescale for the wire to change the angular momentum of the star by a
full cycle, i.e. to vary the eccentricity $e$ of the star from its minimum
value to the maximum and back, is

\begin{align}
 T_{K} &=  \left\{
  \begin{array}{lll}
  \frac{2}{3\pi}\frac{M_\bullet}{\delta m}\left(\frac{a}{R}\right)^{-3}P(a),            &{~\rm Kozai-Lidov,~}   a\le R/2, \\
   \frac{16\sqrt{2}}{3\pi}\frac{M_\bullet}{\delta m}\left(\frac{a}{R}\right)^{1/2}P(a), &{~\rm Non-determ.,~}   a> R/2,
        \end{array} \right.  \label{eqn:TK}
\end{align}

\noindent
where $M_\bullet=4\times10^6~M_\odot$ is the mass of the MBH, and
$P(a)=2\pi\left({a^3}/{GM_\bullet}\right)^{3/2}\simeq1.4\times10^3~\left({a}/{[0.1~{\rm
pc]}}\right)^{3/2}~{\rm yr}$ is the orbital period of the star. The reason for
$R/2$ is a requirement for having all orbits within the radius of the wire,
including the most eccentric ones, $R_{\rm apo}=a\left(1+e\right) \sim 2\,a$,
with $R_{\rm apo}$ the apocenter distance.  Equation~(\ref{eqn:TK}) is a
generalization of the secular Kozai-Lidov (KL) timescale \citep[see][ and
references therehin]{nao13}: (1) In the regime $a\le R/2$ we recover this
well-known secular phenomenon but (2) when $a\ga R/2$, i.e. when stellar orbits
{\em cross} a sphere with a radius of the wire, it provides good approximation
to the non-deterministic (ND), but not necessarily chaotic, evolution of the
stellar orbit.

Admitting that an extended disk is a superposition of wires, one can derive the
corresponding timescale $T'_K$ for the sum of torques to change the orbital
elements of a star in a full cycle \citep{chang09}:

\begin{equation}
1/T'_K=\int_{R_{\rm in}}^{R_{\rm out}}d(1/T_K),\label{eqn:TKprime}
\end{equation}

\noindent
where $R_{\rm in}$ and $R_{\rm out}$ denote the inner and outer radii of the
disk, $d(1/T_K)\propto \delta m=2\pi\Sigma_d(R)RdR$, and $\Sigma_d(R)$ is the
surface density of the disk.  During $T'_K$, when secular evolution
predominates, a star typically oscillates a full cycle between the maximum and
minimum eccentricities, which are predetermined by three orbital parameters,
namely eccentricity, position angle of periapsis ($\omega$), and inclination
angle relative to the disk ($\theta$). At any intermediate stage of that cycle,
the {``instantaneous''} evolution timescale, defined as $t_K(l)\equiv
l/|\dot{l}|$, can be derived from

\begin{equation}
t_K(l)\simeq lT'_K(a)\label{eqn:tK}
\end{equation}

\noindent
\citep[e.g.][]{chang09,chen11}, where $l\equiv\sqrt{1-e^2}$ is the
dimensionless angular momentum and the dot denotes the time derivative.  The
linear dependence on $l$ reflects the coherence of the disk torque during
$t_K(l)$.

The MBH and cusp stars affect the KL-like evolution by perturbing the
orbital parameters $(e,\,\omega,\,\theta)$. We must distinguish two regimes:
(1) At {\em high} $e$, $\omega$ is significantly perturbed, because of the
induced relativistic (GR) precession rate,

\begin{equation}
\dot{\omega}_{\rm GR}=3(GM_\bullet)^{3/2}/\left(l^2c^2a^{5/2}\right),
\end{equation}

\noindent
with $c$ the speed of light. It may even exceed the KL precession rate,

\begin{equation}
\dot{\omega}_K\simeq2\pi/\left(T'_K\,l\right)
\end{equation}

\noindent
\citep[e.g.][]{chang09}. When this happens, the disk coherence is broken and
the KL cycle quenched, hence it defines a boundary to the region in phase space
where the evolution is driven by the disk.  In some loose sense, this boundary
is analogous to the Schwarzschild barrier in galactic nuclei
\citep{mer11,brem14}.  (2) At {\em low} $e$, the perturbation on $\omega$
originates from the total stellar mass $M_*(a)$ enclosed by the orbit. The
Newtonian precession rate

\begin{equation}\label{eqn:omegaM}
\dot{\omega}_M\simeq2\pi lM_*(a)/[M_\bullet P(a)]
\end{equation}

\noindent
may exceed $\dot{\omega}_K$ in this regime, which imposes a second boundary
\citep{chen11}.

Outside these boundaries, evolution of angular-momentum will
be determined by either two-body relaxation, with a characteristic timescale
of

\begin{equation}\label{eqn:t2b}
t_{\rm 2b}(l)\equiv|l/\dot{l}|\simeq l^2\left(M_\bullet/m_*\right)^2P(a)/\left(N\ln\Lambda\right)
\end{equation}

\noindent
\citep[e.g.][]{koc11}, or (scalar) resonant relaxation \citep[RR, ][]{rau96},
on a timescale of

\begin{align}\label{eqn:tRR}
t_{\rm RR,s}(l)\equiv\left|\frac{l}{\dot{l}}\right|\simeq
\frac{l^2}{1-l^2}\left(\frac{M_\bullet}{m_*}\right)^2\frac{P^2(a)}{Nt_\omega}
\end{align}

\noindent
\citep[][who studied the dependence on $e$]{gur07}.  In
Equations~(\ref{eqn:t2b}) and (\ref{eqn:tRR}), $m_*$ denotes the average mass
of one star,  $N=M_*(a)/m_*$ is the number of stars enclosed by the stellar
orbit, $\ln\Lambda=\ln(M_\bullet/m_*)$ is the Coulomb logarithm, and
$t_\omega=2\pi/|\dot{\omega}_M-\dot{\omega}_{\rm GR}-\dot{\omega}_K|$ is the
joint precession timescale combining Newtonian, GR, and KL precessions
\citep{chen13}.

Between the two boundaries is the RER: Any star in it cycles between the
maximum and minimum eccentricities predetermined by $(e,\,\omega,\,\theta)$.
Moreover, the two extrema are evolving. The corresponding timescale is given by
vectorial RR \citep{rau96},  which changes $\theta$ on a timescale of

\begin{align}\label{eqn:tvRR}
t_{\rm RR,v}\equiv\left|\frac{1}{\dot{\theta}}\right|
\simeq \frac{0.3}{(0.5+e^2)^{2}}\frac{M_\bullet}{m_*} \frac{P(a)}{\sqrt{N}}
\end{align}

\noindent
\citep{gur07,eilon09}. Inside RER, $t_{\rm RR,v}$ is longer than the Newtonian and GR
precession timescales, so vectorial RR does not impact the boundaries.  Its
role is to characterize the required time for a star to explore in a
random-walk-fashion the range of maxima and minima in eccentricities fenced in
by the boundaries of the RER.

\subsection{A receding disk}
\label{sec:recedingdisk}

The boundaries of the RER are changing because the properties of the disk have
changed during the past ($1-10$) Myr. We can distinguish two stages in the
evolution of the disk.

\begin{itemize}

\item[(1)] An early phase in which  the disk was mostly gaseous and its
inner edge reached the innermost stable circular orbit (ISCO) at about
$6GM_\bullet/c^2\simeq10^{-6}~{\rm pc}$ \citep{nay05a,levin07}. This disk contained
at least $10^4~M_\odot$ of gas, to trigger fragmentation and star formation
\citep{nay05a}, and it could have been as massive as $(3-10)\times10^4~M_\odot$
according to recent simulations \citep{nay07,Bon08,hobbs09,map12}.  We will
adopt a disk mass of $M_d=3\times10^4~M_\odot$ for this phase.  It is worth
noting that stars formed in the outer disk may migrate inward
\citep{levin07,griv10}, so the disk inside $R=0.04$ pc could contain both gas
and stars.

\item[(2)] Today, after some $(1-10)$ Myr,  the central $0.04$ pc of the
disk is no longer present, because the gas is consumed by either star formation
\citep{nay05a} or black-hole accretion \citep{alex12}, and the stars have had
time to be scattered out of the disk plane due to vectorial RR
\citep{hop06,koc11}.  For this reason, we say the inner edge of the disk has
{\it receded} from the ISCO to the current location of $R_{\rm
in}\simeq1\arcsec\simeq 0.04~{\rm pc}$ \citep[e.g.][]{paumard06}, while the
outer edge is still the same, at $R_{\rm out}\simeq12\arcsec\simeq 0.5~{\rm
pc}$. The present mass of the disk is $M_d=10^4~M_\odot$
\citep{paumard06,bartko10}.

\end{itemize}

In both situations, we modeled the disk surface density as a power-law of
$\Sigma_d(R)\propto R^{-1.4}$ \citep{bartko10},  which leads to a mass of
$6\times10^3~M_\odot$ at $10^{-6}~{\rm pc}<R<0.04~{\rm pc}$ in the early phase.
To derive $M_*(a)$ and $N$ in Equations~(\ref{eqn:omegaM})-(\ref{eqn:tvRR}), we
adopted the broken-power-law model from observations \citep{gen10}, whose
density slop is $\gamma=1.3$ for the inner $0.25$ pc. We assumed an average
stellar mass of $m_*\simeq10~M_\odot$ \citep[also see][ for discussions]{koc11}.
In this model, we have
$t_{\rm RR,v}\simeq1.5\times10^6~(0.5+e^2)^{-2}(a/1\arcsec)^{0.65}~{\rm yr}$ for stars
in the central arcsec of the Galaxy.

\section{Sculpting the Galactic Center}

\subsection{Rapid Evolving Region}
\label{sec:rer}

In Figure~\ref{fig.1} we display the boundaries of the RER. The left panel
corresponds to $R_{\rm in}=10^{-6}~{\rm pc}$ and the right one to $R_{\rm
in}=0.04~{\rm pc}$.  In this $(1-e)$ -- $a$ plane, at any location,
we can estimate the instantaneous evolution timescale as:

\begin{equation}\label{eqn:thori}
\left|\frac{1-e}{\dot{e}}\right|=\frac{e(1-e)}{l^2}t_K(l)\simeq \frac{e(1-e)}{l}T'_K(a).
\end{equation}

We can then identify the lines with constant evolution timescales, i.e. the
contours. We call them ``isochrones'', and depict them as blue dotted curves.
Outside the RER the isochrones are shown in grey and determined by either
two-body scattering, $e\,\left(1-e\right)\,t_{\rm 2b}\,(l)/l^2$,  or RR
process, $e\,\left(1-e\right)\,t_{\rm RR,s}(l)/l^2$, whichever timescale is shorter.

\begin{figure*}
\plotone{fig1.eps}
\caption
   {Mapping evolution timescales in the $a$ -- $(1-e)$ plane.  The thick grey
line on the left-hand-side corresponds to the last stable orbit (LSO) around
SgrA$^*$.  The thick solid black curve is the result of equating the KL
precession rate, $\dot{\omega}_K$, to its relativistic equivalent,
$\dot{\omega}_{\rm GR}$.  Above that curve, dynamical evolution is determined
by KL effect up to the next solid black curve, which comes from equating
$\dot{\omega}_K$ to $\dot{\omega}_M$, the precession rate induced by the
enclosed stellar mass.  The dashed black parallel lines crossing the figures
from the top to the bottom indicate the typical tidal-disruption radii for B,
O, and WR stars.  The blue dotted isochrones, fenced in the region where the
evolution is governed by the KL mechanism, are associated with the logarithms
of the KL timescales given by Equation~(\ref{eqn:thori}).  The grey dotted
isochrones are associated with the logarithms of the two-body-relaxation or RR
timescales, whichever is shorter.  The small orange triangles at the top-right
corners depict the loci of the red giants in the GC. In the left panel, the two
grey boxes depict the expected birth places of S-stars in the binary-separation
and migration-in-disk models \citep[also see][]{antonini13}.  In the right
panel, the dots correspond to: S-stars not associated with the young stellar
disk \citep[small-blue][]{gil09}, the infalling G2 object \citep[also called DSO, see
][ for a different interpretation of its nature]{eckart13} measured at different
times or different wavelengths \citep[small-red][]{gil13rev}, S2/S0-2
\citep[big-blue, the brightest S-star,][]{ghez03,eis03}, and S102/S0-102
\citep[small-cyan,][]{mey12}, the S-star with the shortest period known.}
\label{fig.1} \end{figure*}

Two striking conclusions from a first look at this figure are  (1) stars in the
RER evolve on very short timescales, of the order of $10^{3-5.5}$ yrs, to
complete a full KL cycle.  As discussed previously, after a time of $t_{\rm RR,v}$,
any star at $a<1\arcsec$ would have {\em fully} explored the angular-momentum
range within the RER. (2) As the disk recedes, the boundaries come closer and
the RER shrinks.  Any star that finds itself out of the RER will be ``frozen''
from the point of view of another star which still is in it: The timescales
outside RER are long.

Since short evolution timescale leads to low probability of stellar
distribution, today (right panel of Figure~\ref{fig.1}), the absence of S-stars
within the RER boundaries may be a plausible observational corroboration that
the RER does exit at the GC. At present, the only measured object within the
RER boundaries is G2 \citep[red dots, ][]{gil13rev}.  From its nearby
isochrones, we see that G2 must have been formed less than $10^{5.5}$ years
ago.

\subsection{A close thermalization of the S-stars}\label{sec:sstar}

The two more successful scenarios of depositing B stars close to the GC, i.e.
binary separation and disk migration, place these stars well within the RER
(left panel of Figure~\ref{fig.1}).  These stars are able to sufficiently
mix in angular-momentum space, on a timescale of $t_{\rm RR,v}\simeq0.7$ Myr in the
binary-separation scenario and of $t_{\rm RR,v}\simeq6$ Myr in the disk-migration one.
The latter mechanism (disk migration) requires much longer time to reach the
superthermal distribution in eccentricities becasue of the $e$ dependence of
Equation~(\ref{eqn:tvRR}). We note that these timescales are at least $10$
times shorter than those from the earlier models, which neglected the RER.

The fully-mixed eccentricities do not necessarily have a thermal distribution
(Brem et al., in preparation).  Following the argument that longer evolution
timescale correlates with higher probability distribution, we will have $dN/de
\propto dt/de$, and substituting Equation~(\ref{eqn:thori}) for $dt/de$, we can
derive $dN/de \propto dt/de\propto e/l$.  This distribution function is steeper
than a thermal one, $dN/de\propto e$. The steepness stems from the linear
dependence of the evolution timescale ($l/\dot{l}$) on the orbital angular
momentum $l$,  whereas in the case of two-body relaxation and RR, the evolution
time scales with $1-e^2$.

Figure~\ref{fig.2} compares various cumulative probability distribution
functions (CPDFs) for $e$, derived from two theoretical models --a thermal one
and our RER model-- as well as from observations.  It is clear that compared to
the thermal distribution, the RER one is in better agreement with the
observations.

\begin{figure}
\plotone{fig2.eps}
\caption
   {Cumulative probability distribution functions of the eccentricities of
S-stars, derived from observations (cyan and purple), thermal distribution
($e^2$, black-dashed), our RER model ($1-\sqrt{1-e^2}$, orange-solid), and best
fit to observational data (orange-dotted).} \label{fig.2}
\end{figure}

\subsection{Depleting WR/O stars}\label{sec:wrostar}

For a star of mass $m_*$ and radius $r_*$, if its eccentricity becomes so high
that the orbital pericenter distance from SgrA$^*$ becomes smaller than the
tidal radius

\begin{equation}
r_t\simeq r_*\left(\frac{M_\bullet}{m_*}\right)^{1/3}
\simeq4\times10^{-6}~{\rm pc}\left(\frac{r_*}{R_\odot}\right)
\left(\frac{m_*}{M_\odot}\right)^{-1/3},
\end{equation}

\noindent
it will be tidally disrupted \citep{ree88}.  In Figure~\ref{fig.1} we
show typical $r_t$ associated with B, O, and WR stars. Once stars cross it from
the right to the left, they are lost. This drain is enhanced in the RER by the
shorter and shorter timescales, as the stars progressively move to the left.

To calculate $r_t$, we assume $m_*=(7,\,25,\,60)~M_\odot$ respectively for the
three types of stars \citep{zin07}.  Main-sequence stars less massive than
$7M_\odot$ are below the current detection limit of observations.
Correspondingly, we have adopted $r_*=4~R_\odot$ for main-sequence B stars, and
$r_*=40~R_\odot$ and $80~R_\odot$ respectively for O- and WR stars. B stars can be
envisaged as main-sequence stars, but O- and WR stars are more massive, and
shorter-lived. We hence adopt larger radii for them, 3--4 times larger than
typical radii on the main-sequence, since they have evolved off the main
sequence \citep{paumard06,bartko10}.

We can see in the left panel of Figure~\ref{fig.1} that any WR star in a stripe
defined between $0.15\arcsec\la a\la0.8\arcsec$ and any O star in
$0.2\arcsec\la a\la0.8\arcsec$ will be tidally disrupted, because it will have
explored all the $(1-e)$ space in $\sim 10^6$ yrs.  Similarly, for B stars, the
corresponding stripe is delimited by the narrower zone $0.5\arcsec\la
a\la0.8\arcsec$. In fact, this predicted gap (for B stars) does occur in the
current distribution of S-stars (right panel of Figure~\ref{fig.1}).  If we
had assumed a disk mass of $M_d>3\times10^4~M_\odot$, this gap would have
broadened to incorporat the region where $a<0.5\arcsec$, and it would
contradict current observations.  Therefore, an upper limit to the disk mass
can be derived, approximately $3\times10^4~M_\odot$.

By looking at the left panel again, we realize that only WR/O stars with
$a>0.8\arcsec$ and low $e$ can survive, because they are always outside of the
RER and cannot drift quickly enough to higher $e$.  Indeed, WR/O stars have
been discovered only at $a\ga1\arcsec$ but not inside.  In principle, our
model cannot deplete WR/O stars at $a<0.1\arcsec$, because at such small $a$
the RER does not reach the tidal radii. Observations did not find any WR/O star
there, maybe because the extrapolation of the disk density profile
$\Sigma(R)\propto R^{-1.4}$ results in  $<1$ WR/O star at $R<0.1\arcsec$.

\section{Discussions}\label{sec:dis}

In this Letter we have presented a picture that explains the distribution of
the eccentricities of S-stars and the absence of more massive stars within
$1\arcsec$ of SgrA$^*$.  Our {\em sole} hypothesis is that around $(1-10)$ Myr
ago, the disk had extended down to $R\ll0.04$ pc.  We find that the torque
exerted by the disk creates a region at the GC in which the dynamical evolution
is significantly accelerated as compared to other regions, by a factor ranging
from 10 to 100 times, and we call it the ``rapid evolving region''.

Our scenario agrees with current observations about the nonexistence of an old
segregated cusp in the GC \citep{buc09,do09}, contrary to other works, which
crucially rely on the cusp to thermalize the S-stars
\citep{perets09,madigan11,antonini13,zhang13}.  Because the time that is
needed to randomize angular momentum is now shortened to $(0.7-6)$ Myr, our
model is able to accommodate various possibilities for the formation of
S-stars, while other models rely heavily on when S-stars were brought to the
GC \citep{perets09,antonini13}.

Moreover, our RER scenario unifies two observational facts that have been
thought until now to be disconnected: We successfully populate the observed
range of $e$ for B stars and we can duplicate the observed discontinuity of
WR/O stars above and below $1\arcsec$. Both of them will be established in as
short as $(0.7-6)$ Myr, so we can even unify the origin of all the young
stellar populations in the GC to {\em only one single SF episode}.  This
unification does pose a problem for earlier models: If all B stars formed
simultaneously with WR/O stars, since this must be less than 6 Myr ago
(because WR/O stars cannot be older),  two-body relaxation and RR will fail
to explain the distribution of $e$.

At this stage, it is crucial to theoretically understand the dynamical
response of the old stellar population to the RER, and test it against the
observations of dimmer (than B-type), older stars.  If they match, it would be
a robust evidence that the RER has indeed played a role in sculpting the GC.

\acknowledgments

This work has been supported by the Transregio 7 ``Gravitational Wave
Astronomy'' financed by the Deutsche Forschungsgemeinschaft DFG (German
Research Foundation). We thank Bence Kocsis, Meng Su, and Scott Tremaine
for discussions, and the Kavli Institute for Theoretical
Physics where one part of this work has been completed. This research was
supported in part by the National Science Foundation under Grant No.  NSF
PHY11-25915.  PAS is indebted with Sonia P{\'e}rez for conversations and
extraordinary support.