text:
\section{Introduction}

\noindent
The so-called November Revolution in 1974~\cite{Augustin:1974xw,Aubert:1974js}
triggered numerous theoretical developments in hadron spectroscopy.
Non-relativistic phenomenological models were justified due to the large mass of
the charm and bottom quarks inside mesons, and counterbalanced the lack of
development of Quantum Chromodynamics (QCD) for heavy quarkonium systems at that
time.
Among those early studies we highlight the works of
Eichten~\cite{Eichten:1978tg}, Godfrey~\cite{Godfrey:1985xj},
Stanley~\cite{Stanley:1980zm} or Bhanot~\cite{Bhanot:1978mj}, which employed an
unified and simple framework to phenomenologically study both light- and
heavy-meson spectroscopy.

Those aforementioned phenomenological models considered quarks as low-energy
spectators, which interact through
flavor-independent gluonic degrees of freedom. The short-distance interaction is
assumed to be of perturbative nature, dominated by a single t-channel gluon
exchange (that is, a Coulomb interaction proportional to $\alpha_s$, the strong
coupling constant at some scale). At long distances, non-perturbative effects are
expected to emerge, and are modelled with a linear confining interaction,
confirmed by lattice QCD calculations~\cite{Bali:2000vr}.
Hence, the simplest quark model for heavy quarkonium is the so-called {\it
Cornell potential}\,:
\begin{align}\label{eq:Cornellpot}
V_{\rm Cor}(r)=\sigma\,r\,-\,\frac{C_F\,\alpha_s^{\rm Cornell}}{r}\,.
\end{align}
Here $\alpha_s^{\rm Cornell}$ and $\sigma$ are purely phenomenological constants
of the model with, at first glance,
no direct connection with QCD parameters. $C_F$ is the first $SU(N_c)$ Casimir
[\,$C_F= (N_c^2-1)/(2N_c) =4/3$ for $N_c=3$\,].

Despite of its simplicity and ``ad hoc'' construction, such model is able to
reproduce the heavy quarkonium spectra, which suggests that a relation, even
vague, may exist between the Cornell model and the more theoretically-solid
non-relativistic QCD (NRQCD).
NRQCD~\cite{Lepage:1987gg} is obtained from QCD by integrating out the heavy
quark mass $m_Q$. A perturbative matching can be performed exploiting that
$m_Q\gg\Lambda_{\rm QCD}$, so NRQCD inherits all the light degrees of freedom from QCD.
However, NRQCD mixes soft and ultrasoft scales, complicating the power-counting.
Solutions to this problem appeared with the development of different
EFTs such as velocity NRQCD (vNRQCD)~\cite{Luke:1999kz} and
potential NRQCD (pNRQCD)~\cite{Pineda:1997bj,Brambilla:1999xf}, which describe
the interactions of a non-relativistic
system with ultrasoft gluons, systematically organizing the perturbative
expansions in $\alpha_s$ and the velocity of heavy quarks. Such EFTs only
include the relevant degrees of freedom for $Q\overline Q$ systems near
threshold, integrating out the rest of degrees of freedom.

In this work we explore the connection between the simplest realization of the
Cornell model against NRQCD.
To that end we compare the mass of the lowest-lying $Q\overline{Q}$ bound
states, observables that can be reliably predicted both in the theory and the
model, varying the quark mass and the strong coupling constant. We show that the
Cornell potential agrees for large values of $r$ with the QCD static potential
once the latter is expressed in terms of the MSR mass and improved with
all-order resummation of large renormalon-related logs via R-evolution~\cite{Hoang:2017suc}.
Our R-improved static potential also compares nicely with lattice QCD simulations
from Refs.~\cite{Bazavov:2014pvz,Bazavov:2017dsy}.

A more complete description of the methods and further results can be found in Ref.~\cite{Mateu:2018zym}.

\section{Cornell Potential and QCD Static Potential for $\mathbf{Q\overline Q}$}\label{sec:StaticPot}

\noindent
As mentioned above, our main goal is to compare the Cornell potential predictions with NRQCD,
so relations among QCD fundamental constants and Cornell model parameters can be obtained.

On the one hand, the Schr\"odinger equation for the Cornell potential of Eq.~\eqref{eq:Cornellpot} is
solved numerically using the Numerov algorithm~\cite{Numerov1}, from where we calculate the mass
of the bound states. However, such potential is not sufficient to predict the  $\Upsilon(1S)-\eta_b(1S)$
or  $\chi_{bJ}$ multiplet (with $J=\{0,1,2\}$) mass splittings. Further spin dependence is needed in the
Cornell static potential. For that reason, it is important to add $1/m_Q^2$ terms to take into account the
spin-spin, spin-orbit and tensor interactions, breaking the leading power degeneracy~\cite{DeRujula:1975qlm}.
Such contributions emerge from the s-channel gluon exchange and the leading relativistic
corrections of the t-channel gluon and confinement interactions~\cite{Godfrey:1985xj}\,:
\begin{align}
V_{\rm SS}^{\rm OGE}(r)&=\frac{8\alpha_s^{\rm Cornell}}{9m_Q^2\,r^2}(\vec S_1\cdot \vec S_2)\,\delta(r)\,,\label{eq:subleading}
&V_{\rm LS}^{\rm OGE}(r)&=\frac{2\alpha_s^{\rm Cornell}}{m_Q^2\,r^3}\,(\vec L\cdot\vec S)\,,\\
V_{\rm LS}^{\rm CON}(r)&=-\frac{\sigma}{2m_Q^2\,r}\,(\vec L\cdot \vec S)\,,
&V_{\rm T}^{\rm OGE}(r)&=\frac{\alpha_s^{\rm Cornell}}{3m_Q^2\,r^3}\,S_{12}\,,\nonumber
\end{align}
being $\vec S=\vec S_1+\vec S_2$ the total spin, $\vec L$ the relative orbital momentum and $S_{12}$ the
tensor operator of the $Q\overline Q$ bound state, defined as
$S_{12}=2\,(\vec S_1\cdot\hat r)(\vec S_2\cdot\hat r)-(\vec S_1\cdot\vec S_2)$, with $\hat r= {\vec r}/r$.
With no significant loss of precision, such contributions can be calculated
using first-order perturbation theory, given the large mass of the heavy quark.

Using the above potential we are able to fit the parameters to the low-lying experimental bottomonium spectrum, those with $n_p=\{1,2\}$, whose masses are taken from the PDG~\cite{Tanabashi:2018oca}. Such fit gives the values\,:
\begin{align}\label{eq:Cornell-Fit}
\{m_b^{\rm Cornell}\,,\,\sigma\,,\,\alpha_s^{\rm Cornell}\} & =
\{4.733 \pm 0.018\,{\rm GeV}\,,\, 0.207 \pm 0.011\,{\rm GeV}^2\,,\, 0.356 \pm 0.015\,\}
\end{align}
where the uncertainties correspond to the $68\,\%$ confidence level for each parameter. Given the extremely precise nature of the experimental masses, the uncertainties on the parameters that come out of the fit are penalized, increasing their values by the square root of the reduced $\chi^2$ at its minimum.

On the other hand, within NRQCD we can define the static potential as the color-neutral interaction
between two infinitely heavy color-triplet states. Such potential has been calculated up to $\mathcal{O}(\alpha_s^4)$,
at which point the potential becomes time-dependent and the static approximation breaks down.
Such feature derives in a dependence on the ``ultrasoft'' scale $\mu_{\rm us}$.

In position space and taking the charm quark as massless, the perturbative contribution to the static QCD potential can
be written as follows\,:
\begin{equation}\label{eq:static}
V^{(n_{\ell})}_{\rm QCD}(r,\mu) = -\,C_F\,\frac{\alpha^{(n_{\ell})}_s(\mu)}{r}\sum_{i=0}\sum_{j=0}^{i}
\biggl(\frac{\alpha^{(n_{\ell})}_s(\mu)}{4\pi}\biggr)^{\!\!i}\, a_{i,j}(n_{\ell})\,\log^{\,j}(r\,\mu\, e^{\gamma_E})-\,\frac{9\,C_F}{4\pi}\,\alpha_s^{(n_{\ell})}(\mu)^4\,\frac{1}{r}\log\,(\mu_{\rm us}\,r)\,,
\end{equation}
where the coefficients $a_{i,0}$ are known to four loops~\cite{Smirnov:2009fh,Anzai:2009tm}
and $a_{i,j>0}$ can be derived from the former requiring that $V_{\rm QCD}$ does not depend on $\mu$\,:
\begin{equation}
a_{i,j}(n_{\ell}) = \frac{2}{j}\sum_{k=j}^{i}\,k\,a_{k-1,j-1}(n_{\ell})\,\beta_{i-k}(n_{\ell})\,.
\end{equation}

The last term in Eq.~\eqref{eq:static} depends on the ultrasoft factorization scale $\mu_{\rm us}$. Following e.g.\ Ref.~\cite{Tormo:2013tha}, we will take the following expression for the ultrasoft factorization scale\,:
\begin{align}
\mu_{\rm us}=\frac{N_c}{2}\,\mu\,\alpha^{(n_{\ell})}_s(\mu)\,,
\end{align}
that takes into account the power counting of pNRQCD~\cite{Pineda:1997bj,Brambilla:1999xf}.

If we plot the Cornell potential versus the QCD static potential (Fig.~\ref{fig:polebottom}) we appreciate a bad
perturbative behaviour of the latter, visible as a vertical shift of the potential in the region between $0.05\,$fm and $0.2\,$fm between different orders. Indeed, this is a clear manifestation of the $u=1/2$ renormalon of the QCD static potential, which is $r$-independent, depending only on the coefficients of the QCD beta function.
Furthermore, no addition of further orders brings the QCD static potential closer to the Cornell model.

It is important, thus, to cancel the static QCD potential renormalon, which exactly matches that of the pole mass except for a $-2$ factor, so that the static energy $E_{\rm stat}(r,\mu) = 2\,m_Q^{\rm pole} + V_{\rm QCD}(r,\mu)$ is renormalon free. A standard approach is to re-write the pole mass in terms of a short-distance scheme
to make the cancellation manifest. For the cancellation to take place one also needs to express the
perturbative series $\delta m_Q^{\rm SD}$ [\,Eq.~\eqref{eq:MSRdef}\,] that relates the pole mass with a short-distance mass $m_Q^{\rm SD}$ in terms of $\alpha_s(\mu)$.
When done so, powers of $\log(\mu/R)$\,\footnote{We denote $R$ as the infrared scale of the short-distance mass.} in $\delta m_Q^{\rm SD}$ will emerge  which
may become large if $\mu$ and $R$ are very different. Besides, from the static QCD potential definition in
Eq.~\eqref{eq:static}, $\mu$ should depend on $r$ such that $\log(r\,\mu\, e^{\gamma_E})\sim\mathcal{O}(1)$,
therefore mass schemes
with a fixed value of $R$, such as the $\overline{\mathrm{MS}}$, are disfavored.

In order to avoid the canonical $\mu=1/r$ to rapidly acquire values in the non-perturbative region,
we smoothly freeze it to $1\,$GeV once it reaches this value, using a transition function between $r=0.08\,$fm and $0.2\,$fm~\cite{Mateu:2018zym}. Fig.~\ref{fig:Profile} shows a graphical representation of the ``profile function'', as will be referred to.
This scale setting is similar to those implemented in the  renormalon-subtracted scheme at
Refs.~\cite{Pineda:2013lta,Peset:2018ria,Peset:2018jkf}, but with a smoother transition between the canonical and frozen regimes.

\begin{figure}[tbh!]\centering
\floatbox[{\capbeside\thisfloatsetup{capbesideposition={right,top},capbesidewidth=6cm}}]{figure}[\FBwidth]
{\caption{Dependence of the static potential renormalization scale in Eq.~\eqref{eq:static} with $r$. In solid
blue the profile function is shown, while dashed green and dashed red show the behavior
for small ($1/r$) and large (constant) $r$, respectively. The positions in which the piece-wise function changes its
functional form are signaled by vertical dashed black lines.}\label{fig:Profile}}
{\includegraphics[width=\linewidth]{figs/Profile}}
\end{figure}

Following Ref.~\cite{Mateu:2017hlz} we use the MSR mass~\cite{Hoang:2017suc} and choose $\mu = R$ to simultaneously minimize logs in the potential and in $\delta m_Q^{\rm MSR}$.
The $\mathrm{MSR}$ mass scheme was created as a natural extension of the $\overline{\mathrm{MS}}$ mass for renormalization scales
below the heavy quark mass. It can be directly defined from the $\overline{\mathrm{MS}}$-pole mass relation\,:

\begin{align}\label{eq:MSRdef}
\delta m_Q^\mathrm{MSR} \equiv
m_Q^{\rm pole}-m_Q^\mathrm{MSR}(R) & =R\sum_{n=1}^\infty \sum_{k=0}^n d_{n,k}(n_{\ell})\biggl(\frac{\alpha^{(n_{\ell})}_s(\mu)}{4\pi}\biggr)^{\!\!n}
\ln^k\Bigl(\frac{\mu}{R}\Bigr)\,,
\end{align}
where we see that, apart from the logarithmic dependence on the scale $\mu$, a logarithmic and linear dependence on an infrared scale R appears. The $d_{n,0}>$ are derived from the $\overline{\mathrm{MS}}$-pole mass relation, and depend on the specific method to change from a scheme with $(n_\ell+1)$ dynamical flavors to another with only $n_\ell$. On the one hand, the \emph{practical $\mathrm{MSR}$ mass} (MSRp) employs the threshold matching relations of the strong coupling to express $\alpha_s^{(n_\ell+1)}(\overline{m})$ in terms of $\alpha_s^{(n_\ell)}(\overline{m})$ and, on the other hand, the \emph{natural $\mathrm{MSR}$ mass} (MSRn) directly integrates out the heavy quark Q from the $\overline{\mathrm{MS}}$-pole relation, setting to zero all diagrams containing heavy quark loops.

The R dependence of the $\mathrm{MSR}$ mass is described by\,:
\begin{equation}\label{eq:Revol}
-\,\frac{{\rm d}}{{\rm d} R}m_Q^{\mathrm{MSR}}(R)=\gamma^R[\alpha_s^{(n_{\ell})}(R)]=\sum_{n=0}^\infty \gamma_n^R \left(\frac{\alpha_s^{(n_{\ell})}(R)}{4\pi}\right)^{\!\!n+1},
\end{equation}
where $\gamma_n^R(n_{\ell}) = d_{n+1,0}(n_{\ell})-2\sum_{j=0}^{n-1} (n-j)\,\beta_j(n_{\ell})\, d_{n-j,0}(n_{\ell})\,$ are the R-anomalous dimension coefficients~\cite{Hoang:2017suc}, and the renormalon cancels between the first term and the sum. The solution of the RGE in Eq.~\eqref{eq:Revol} sums up
powers of $\log(R_1/R_2)$ to all orders in perturbation theory\,:
\begin{align}\label{eq:MSRsum}
\Delta^{\rm MSR}(n_{\ell},R_1,R_2)\equiv m_Q^{\rm MSR}(R_2) - m_Q^{\rm MSR}(R_1) =
\!\int_{R_1}^{R_2}\!{\rm d} R\,\gamma_n^R[n_{\ell},\alpha_s^{(n_{\ell})}(R)]\,.
\end{align}
Besides, the MSR mass can be easily related to the $\overline{\mathrm{MS}}$ mass at the scale \mbox{$R = \overline{m}_Q(\overline{m}_Q)$}, from where it can then run down to any value of $R<\overline{m}_Q$.

Using the previous mass scheme we can define a short-distance potential, which is renormalon free and independent
of the heavy quark mass. It will depend, though, on a fixed scale where the renormalon is subtracted, denoted as $R_0$.
Manipulating the static energy formula we arrive to the convenient expression\,:
\begin{align}\label{eq:MSRstaticE}
E_{\rm stat}(r,\mu) \,&\,=\, 2\,m_Q^{\rm pole} + V_{\rm QCD}(r,\mu)
=\, 2\,m_Q^{\rm MSR}(R_0) + 2\,\delta m_Q^{\rm MSR}(R_0,\mu)
+ V_{\rm QCD}(r,\mu)\nonumber\\
\,&\,\equiv\,2\,m_Q^{\rm MSR}(R_0) + V^{\rm MSR}_{\rm QCD}(r,\mu,R_0)\,.
\end{align}
However, to avoid the appearance of large logs of $R_0/\mu$ in $\delta m_Q^{\rm MSR}(R_0,\mu)$ we can use \mbox{R-evolution}
to sum up large logs and express it in terms of $\delta m_Q^{\rm MSR}(R,\mu)$\,:
\begin{align}
\delta m_Q^{\rm MSR}(R_0) &= \delta m_Q^{\rm MSR}(R_0) + \delta m_Q^{\rm MSR}(R) - \delta m_Q^{\rm MSR}(R) = \delta m_Q^{\rm MSR}(R) + m_Q^{\rm MSR}(R) - m_Q^{\rm MSR}(R_0) \nonumber\\
& = \delta m_Q^{\rm MSR}(R) + \Delta^{\rm MSR}(R,R_0)\,.
\end{align}
Hence, the following R-improved expression for the MSR-scheme Static QCD Potential can be written as\,:
\begin{align}\label{eq:MSRStat}
V^{\rm MSR}_{\rm QCD}(r,\mu,R_0) =  V_{\rm QCD}(r,\mu) + 2\, \delta^{\rm MSR}(R,\mu) +2\, \Delta^{\rm MSR}(R,R_0)\,.
\end{align}
The R-improved static potential is similar to the Renormalon Subtracted scheme used in Ref.~\cite{Brambilla:2009bi},
based on \cite{Pineda:2001zq}, and to the analysis of Ref.~\cite{Sumino:2003yp} based on renormalon dominance.

\begin{figure*}[tbh!]
\center
{\label{fig:polebottom}\includegraphics[width=0.49\textwidth]{figs/Static-pole}~~}
{\label{fig:MSRnbottom}\includegraphics[width=0.47\textwidth]{figs/MSRn-stat-pot} }
\caption{\label{fig:Static} Comparison of the Cornell potential (dashed black line) with the static QCD potential at
N$^3$LO (red), N$^2$LO (green), NLO (cyan) and LO (blue) for bottomonium.
The left panel shows the static QCD potential in the pole scheme, and the right panel uses the MSR scheme with the
reference scales \mbox{$R_0=1\,$GeV}. The parameters of the Cornell potential are fixed to best-reproduce the $n_p=\{1,2\}$
experimental bottomonium states [\,Eq.~\eqref{eq:Cornell-Fit}\,], while the Static QCD potential uses the
world average value $\alpha_s(m_Z) = 0.1181$. The left plot uses the canonical scale $\mu = 1/r$, but the right plot uses
the profile function shown in Fig.~\ref{fig:Profile}. Results for charmonium are shown in Ref.~\cite{Mateu:2018zym}.}
\end{figure*}

The result for Eq.~\eqref{eq:MSRStat} is shown in Figs.~\ref{fig:MSRnbottom} for bottomonium for the value
\mbox{$R_0=1\,$GeV}.\,\footnote{Results for charmonium are similar and can be found in Ref.~\cite{Mateu:2018zym}.}
The choice of $R_0$ is set so it is identical to the value at which the renormalization scale freezes.
We can see how the static MSR potential converges nicely towards the Cornell model for moderate values of $r$. The
agreement for large distances improves with each new perturbative order addition.
For larger values of the radius $\log(r\mu)$ becomes large, as $\mu$ freezes, which makes perturbation theory unreliable. At small distances all orders agree very well due to the small value of $\alpha_s$,
but disagree with the Cornell model. In conclusion, the Cornell model and QCD agree for moderate values of
$r$, but disagree in the ultraviolet as the model does not incorporate logarithmic modifications due to the running
of $\alpha_s$. An interesting question is, then, if this difference in the UV can be absorbed in the definition of the quark mass. This confirms the claims of Refs.~\cite{Sumino:2003yp,Sumino:2004ht} where, using renormalon dominance arguments and in the framework of
the operator product expansion of pNRQCD, it is shown that perturbation theory alone should be capable of describing both the Coulomb and linear behavior of the static potential, and that non-perturbative corrections start at
$\mathcal{O}(\Lambda_\mathrm{QCD}^3r^2)$.

The aforementioned R-improved QCD static potential can also be compared to lattice simulations.
Our result for $\sigma$ in Eq.~\eqref{eq:Cornell-Fit} is in very nice agreement with the lattice determinations of Ref.~\cite{Koma:2006fw} ($\sigma = 0.2098 \pm 0.0009\,$GeV$^2$) and
Ref.~\cite{Kawanai:2013aca} ($\sigma = 0.206 \pm 0.010\,$GeV$^2$). The direct comparison of
$\alpha_s^{\rm Cornell}$ parameter and lattice analyses differs by a factor of roughly $2$.
However, such result should be taken with caution, as in the static potential at short distances loop corrections modify the short-distance $1/r$ behavior. Thus, this discrepancy should be of little concern.

Regarding the static potential, we use the results of Ref~\cite{Bazavov:2014pvz} with lattice
spacing \mbox{$a=0.04\,$fm}. These range between $0.039\,{\rm fm} \leq r \leq 0.84\,$fm, with an average relative precision
of $2.6\,\%$. This dataset is complemented with results from \cite{Bazavov:2017dsy} with a smaller lattice spacing of
$a = 0.025\,$fm, covering values of the radius as small as $0.024\,$fm. For the latter we only consider data with
\mbox{$r\le0.25\,$fm}, since uncertainties in Ref.~\cite{Bazavov:2017dsy} are larger for higher values of $r$.
Furthermore, in the range \mbox{$0.024\,{\rm fm}\leq r \leq0.25\,$fm}, Ref.~\cite{Bazavov:2017dsy}
has more density of points than \cite{Bazavov:2014pvz}.
The full dataset is shown in Fig.~\ref{fig:Lattice} as black dots with error bars.
The static potential is only dependent on $\alpha_s$ and an arbitrary additive constant, a vertical offset which can be related to the subtraction scale $R_0$ in our R-improved version. For that reason, we perform a two-parameter fit to the lattice data of our \mbox{R-improved} static QCD potential with the scale setting
shown in Fig.~\ref{fig:Profile}. For this simple comparison we use $\mu = R$ and no large ultrasoft logs resummation. The lattice data is assumed to be statistically independent, as the correlation matrix is currently unknown. The results of the fit are\,:
\begin{align}\label{eq:alphaSfit}
\alpha_s^{(n_f = 5)}(m_Z) = 0.1168\,,\qquad R_0 = 1.024\,{\rm GeV}.
\end{align}
Note that the remarkable agreement of the value of $R_0$ to that previously employed to compare with the Cornell model potential.
No fit uncertainties are shown as they would not reflect the actual accuracy reachable by this procedure.
The comparison using the above best-fit values for $\alpha_s$ and $R_0$ is shown as a red
line in Fig.~\ref{fig:Lattice}. The blue line shows the R-improved static QCD potential with the same
values for the parameters, but implementing a fully canonical profile $\mu = 1/r$. Our profile completely agrees with lattice QCD results up to distances of approximately
$1\,$fm, while the canonical profile behaves in an unphysical manner for $r\gtrsim 0.2\,$fm. Such findings can have a direct impact on the precise determination of $\alpha_s$ constant from fits to lattice simulations.

\begin{figure}[tbh!]\centering
\floatbox[{\capbeside\thisfloatsetup{capbesideposition={right,top},capbesidewidth=6cm}}]{figure}[\FBwidth]
{\caption{Comparison of the R-improved perturbative static potential with lattice QCD results for $n_f = 3$ dynamical flavors. We use $\mu=R$, with $\mu$ set to the profile shown in Fig.~\ref{fig:Profile}. We determine
$\alpha_s$ and $R_0$ fitting our pNRQCD theoretical expression to the lattice
simulations~\cite{Bazavov:2014pvz,Bazavov:2017dsy}.}\label{fig:Lattice}}
{\includegraphics[width=\linewidth]{figs/Lattice-fits}}
\end{figure}

\section{Bottomonium spectrum with a floating bottom quark mass}

\noindent
In this section we will describe the procedure to fit the Cornell potential with NRQCD using synthetic
spectra with a floating bottomonium mass. This complete bottomonium spectrum, up to $n_p = 2$ for
arbitrary \mbox{$\overline{m}_b\equiv\overline{m}_b(\overline{m}_b)$} is constructed using NRQCD up to N$^3$LO~\cite{Kiyo:2014uca}.
In the pole mass scheme~\cite{Penin:2002zv,Beneke:2005hg,Kiyo:2014uca}, the energy of a non-relativistic $Q\overline{Q}$
bound state, univoquely determined by the $(n_p,j,\ell,s)$ quantum numbers and with $n_{\ell}$ massless active flavors, can be written
as\,:\,\footnote{See Ref.~\cite{Mateu:2018zym} for more details.}
\begin{align}\label{eq:EXpole}
& M^{(n_{\ell})}_{n_p,j,\ell,s}(\mu) = 2\,m_Q^{\rm pole}\!
\Bigg[1-\frac{C_F^2\,\alpha^{(n_{\ell})}_s(\mu)^2}{8n_p^2}\sum_{i=0}^{\infty}
\bigg(\frac{\alpha^{(n_{\ell})}_s(\mu)}{4\pi}\bigg)^{\!\!i}\,
\varepsilon^{i+1}P_i(L_{n_{\ell}})\Bigg],
\end{align}
with
\begin{align}\label{eq:NRlog}
& L_{n_{\ell}}=\log\!\bigg(\frac{n_p\mu}{C_F\alpha_s^{(n_{\ell})}(\mu)m_Q^{\rm pole}}\bigg)+ H_{n_p+\ell}\,,\qquad
P_i(L) = \sum_{j=0}^i\,c_{i,j}\,L^j\,,
\end{align}
where $H_{n}$ is the harmonic number. In Eq.~\eqref{eq:EXpole} $\varepsilon$ acts as a bookkeeping
parameter that properly implements the so called $\Upsilon$-expansion~\cite{Hoang:1998hm}.
The $c_{i,0}$ coefficients have been computed up to $i=3$~\cite{Brambilla:2001qk,Kiyo:2013aea},
while the $c_{i,j>0}$ coefficients can be directly obtained from the latter $c_{i,0}$ imposing $\mu$
independence of the quarkonia mass. We denote the sum in Eq.~\eqref{eq:NRlog} truncated to
$\mathcal{O}(\varepsilon^{n+1})$ as the N$^n$LO result. This formula does not include the resummation of
large ultrasoft logarithms.\,\footnote{See Ref.~\cite{Peset:2018jkf} for a recent summation to N$^3$LL precision for P-wave states.}

Due to the already discussed $u=1/2$ renormalon in the QCD static potential, it is essential to express the static energy in Eq.~\eqref{eq:EXpole} in terms of a short-distance mass, exactly as in the static energy.
Following the results of Sec.~\ref{sec:StaticPot} and the analysis in Ref.~\cite{Mateu:2017hlz} we will employ the MSR mass~\cite{Hoang:2017suc} in our analysis.

The final goal of this study is to obtain relations among the Cornell model constants (especially the Cornell mass) and the QCD fundamental parameters $\alpha_s$ and $m_b$. To achieve this, we will perform a calibration scanning over these two parameters.
Specifically, we will create templates of bottomonium spectra in reasonable ranges of values for the bottom mass and the strong coupling constant, and from there we will obtain the dependence of the Cornell model parameters in terms of $m_b$ and $\alpha_s$. For the former we will vary the $\overline{\mathrm{MS}}$ mass between $4$ and $8\,$GeV in steps of $500\,$MeV.
For the strong coupling constant, we consider $\alpha_s^{(n_\ell=5)}(m_Z)$ between $0.114$ and $0.12$ for $\overline{m}_b = 4.2\,$GeV.

We generate QCD predictions for the $n_p=\{1,2\}$ bottomonium states varying the two renormalization scales $\mu$ and $R$ in a correlated way (see Ref~\cite{Mateu:2018zym} for details). The ranges are selected so the argument of the logarithm in Eq.~\eqref{eq:NRlog} in the MSR scheme ranges between $1/2$ and $2$. Indeed, this range in general depends on the value of $m_b$. Hence, we set a lower limit dependent on the bottom mass:
\begin{align}\label{eq:mui}
\mu_1^{\rm min} = 0.638\,{\rm GeV} + 0.209\,\overline{m}_b\,\,,\,\,\,\,
\mu_2^{\rm min} = 0.510\,{\rm GeV} + 0.120\,\overline{m}_b\,,
\end{align}
where the subindex $1,\,2$ denotes the principal quantum number of the heavy quarkonium state.
The upper limit is fixed to $4\,$GeV, as no dependence of $m_b$ is found for $\mu$ or $R$.

Given the synthetic bottomonium spectrum for a specific $m_b$ and $\alpha_s$, we will fit the Cornell model parameters via a simple statistical regression analysis. The generated QCD pseudo-data at LO, NLO, N$^2$LO and
N$^3$LO can be thought of as a set of highly correlated experimental measurements, which makes a traditional $\chi^2$
fit with a non-diagonal covariance matrix impossible, due to the d'Agostini bias. However, the theoretical
covariance matrix is the only uncertainty source we have, so it is impossible to write down a $\chi^2$ function.
We follow the strategy of Ref.~\cite{Mateu:2017hlz}, where the QCD renormalization scales were varied in a correlated way in terms of two dimensionless variables \mbox{$\mu_n = \mu_n^{\rm min} + x\,(4\,{\rm GeV} - \mu_n^{\rm min})$},
\mbox{$R_n = \mu_n^{\rm min} + y\,(4\,{\rm GeV} - \mu_n^{\rm min})$}, with \mbox{$n = \{1, 2\}$, $0\le \{x,y\}\le1$} and
$\mu_i$ defined in Eq.~\eqref{eq:mui}. The employed regression $\chi^2$ function is, then\,:
\begin{align}\label{eq:chi2Reg}
\chi^2(x,y,\overline{m}_b)=\frac{\sum_i\big(M_i^{\rm QCD}(x,y,\overline{m}_b)-M_i^{\rm Cornell}\big)^{2}}
{\frac{\chi^2_{\rm min}(x,y)}{\rm d.o.f.}}\,,
\end{align}
where, actually, $\chi^2_{\rm min}(x,y)$ is known only after the minimization is carried out, but allows the $\chi^2$
function to be dimensionless. For a given value of $\overline{m}_b$, all possible
values of $\{x,y\}$ are scanned, and the parameters of the Cornell model are determined for each pair. Scanning over all values of $\overline{m}_b$ and $\alpha_s$, functions of the Cornell model
parameters in terms of these QCD fundamental quantities are obtained.

\section{Calibration of Cornell model}\label{sec:results}
\noindent
First we focus on the dependence of the Cornell mass parameter with the short-distance QCD mass.
For convenience, we have produced our bottomonium spectra as a function of $\overline{\mathrm{MS}}$ mass, although this scheme should be thought as a coupling constant rather than a kinematic mass.
Following the analysis of Secs.~\ref{sec:StaticPot}, we use the MSR mass with small values of $R$, which is a kinetic mass free from renormalon ambiguities. When we calibrate the Cornell mass versus
\mbox{$m_b^{\rm MSR}(R={1\rm GeV})$} we observe a linear relation, with a slope very close to (and compatible with) unity ($0.995\pm 0.026$) and a constant term compatible with zero ($0.05 \pm 0.19$).
Furthermore, this correspondance is achieved for all considered values of $\alpha_s$ and different orders~\cite{Mateu:2018zym}.
In Fig.~\ref{fig:Money} we show the linear relation at N$^3$LO for the world average value of $\alpha_s$, even though this correspondance is achieved for all considered values of $\alpha_s$ and different perturbative orders.

\begin{figure}[tbh!]\centering
\floatbox[{\capbeside\thisfloatsetup{capbesideposition={right,top},capbesidewidth=6cm}}]{figure}[\FBwidth]
{\caption{Upper plot\,: Dependence of the Cornell mass parameter with the MSR bottom mass at
$R=1\,$GeV at N$^3$LO, using the world average value for the strong coupling constant. \newline
Lower plot\,: Difference between the
Cornell and MSR masses as a function of $m_b^{\rm MSR}(R=1\,{\rm  GeV})$. \newline
The error bars represent the quadratic sum of fit and perturbative
uncertainties. The blue line and band in the upper plot shows a linear fit to the points, where the individual uncertainties are taken as uncorrelated. In the lower plot, those correspond to the weigthed average of the differences and the regular average of the uncertainties, respectively.}
\label{fig:Money}}
{\includegraphics[width=0.92\linewidth]{figs/Money-Plot}}
\end{figure}
The most relevant result of this work is, however, the difference between the Cornell model mass parameter and the MSR mass\,:
\begin{align}\label{eq:main-result}
m_b^{\rm Cornell} = m_b^{\rm MSR}(R = 1\,{\rm GeV}) + [\,0.023 \pm 0.086\,{\rm GeV}\,]\,,
\end{align}
calculated as the weigthed average of the individual
$m_b^{\rm Cornell} - m_b^{\rm MSR}(1\,{\rm GeV})$ values, taking the uncertainty as the regular average of individual uncertainties.

Finally, we show the dependence of the remaining Cornell parameters, $\alpha_s^{\rm Cornell}$
and  $\sigma$, with the QCD quantities $\alpha_s$ and $m_b^{\rm MSR}$.
The variation of $\alpha_s^{\rm Cornell}$ with the latter QCD parameters is shown in the upper plots of Fig.~\ref{fig:alpha-sigma} as blue dots. The red solid line corresponds to the QCD strong coupling constant evaluated at the non-relativistic scale $\mu_{\rm NR}$, chosen so the argument of the log in Eq.~\eqref{eq:NRlog} is of the order of one. As our fit employs $n_p=\{1,2\}$ states, the average value $\bar n_p=1.5$ is taken.
\begin{figure*}[tbh!]
\center
{\includegraphics[width=0.46\textwidth]{figs/alpha-vs-mass}~~~~}
{\includegraphics[width=0.47\textwidth]{figs/alpha-vs-alpha} }\vspace*{-0.3cm}
{\includegraphics[width=0.46\textwidth]{figs/sigma-vs-mass}~~~~}
{\includegraphics[width=0.47\textwidth]{figs/sigma-vs-alpha} }
\caption{Dependence of the Cornell parameters $\alpha_s^{\rm Cornell}$ (upper two panels)
and $\sigma$ (lower two panels) with the MSR bottom mass (leftmost two panels) and the QCD coupling constant at the
Z-pole (rightmost two panels). The two upper plots also show, with a solid red line, the strong coupling constant evaluated at a characteristic non-relativistic scale.}
\label{fig:alpha-sigma}
\end{figure*}

The scale $\mu_{\rm NR}$ is obtained by choosing $R=\mu$ and solving numerically the equation\,:
\begin{align}\label{eq:muNR}
C_F\alpha_s^{(n_f=4)}(\mu_{\rm NR})\,m_b^{\rm MSR}(\mu_{\rm NR}) = \frac{3}{2}\,\mu_{\rm NR}\,.
\end{align}
This scale is determined for different values of $\alpha_s^{(n_f=5)}(m_Z)$ and the bottom mass, finding a
remarkable agreement between $\alpha_s^{\rm Cornell}$ and $\alpha_s^{(n_f=4)}(\mu_{\rm NR})$ within perturbative uncertainties. Besides the similar order of magnitude, its dependence on $m_b$ and $\alpha_s(m_Z)$ follows the same trend. This simple analysis disfavors quark models which consider a QCD-like running for $\alpha_s^{\rm Cornell}$, with a scale given by the reduced mass of the $Q\overline Q$ pair.

The lower plots in Fig.~\ref{fig:alpha-sigma} represents the dependence of $\sigma$ with the bottom mass and the strong coupling constant. It shows a remarkable independence on the bottom mass, as expected for a static potential. Taking the average of all dots we obtain \mbox{$\sigma^{\rm Cal.} = 0.176 \pm 0.088\,$GeV$^2$}, which compares well with that obtained from the fit to the bottomonium experimental data in
Eq.~\eqref{eq:Cornell-Fit} \mbox{$\sigma^{\rm fit} = 0.207 \pm 0.011\,{\rm GeV}^2$}. Regarding the dependence on  $\alpha_s^{(n_f=5)}(m_Z)$, a small positive correlation is envisaged, although given the uncertainties no strong conclusions can be drawn. Anyway, some dependence of $\sigma$ with $\alpha_s$ is expected since, as argued in Refs.~\cite{Sumino:2003yp,Sumino:2004ht} and confirmed in Sec.~\ref{sec:StaticPot} of this work, the linear rising term in the static potential is of perturbative nature.

\bibliographystyle{JHEP}


text:
\section{Motivation}

Scale invariance is a global spacetime symmetry under the rescaling of coordinates and fields according to their canonical dimensions:
\begin{equation}
x^{\mu}\to tx^{\mu}~,~~S(x)\to tS(tx)~,~~V_{\mu}(x)\to tV_{\mu}(tx)~,~~F(x)\to t^{3/2}F(tx)~,
\label{1}
\end{equation}
where $S(x),~V_{\mu}(x),~F(x)$ denote scalar, vector and fermion fields, respectively\footnote{Relativistic scale-invariant theories are typically accompanied by a symmetry under special conformal transformations. Together with the $ISO(1,3)$ group of relativistic invariance, scale (dilatation) and special conformal transformations form the 15-parametric conformal group $SO(2,4)$.  In particle physics, however, we are interested in (spontaneously) broken conformal invariance, $SO(2,4)\to ISO(1,3)$, and the only physical manifestation of this breaking is represented by a (pseudo)Goldstone scalar particle, the dilaton, associated with the scale invariance. Therefore, in what follows we will be interested only in scale transformations. }.  In phenomenologically relevant theories, however, scale invariance is an anomalous symmetry, that is, it is an exact symmetry of the classical action, $\sim{\cal O}(\hbar^0)$, and is broken by quantum corrections $\sim {\cal O}(\hbar)$.  An obvious question about scale invariance arises then: Does it matter if it is not actually a symmetry of a full quantum theory?

Firstly, we have a precedent in nature. It is quite remarkable, that most of the mass of the visible matter in the universe carried by protons and neutrons is essentially generated within the scale-invariant QCD through the mechanism of dimensional transmutation. Therefore, it is conceivable to think that the masses of all elementary particles   are also generated in a similar manner. In fact, in technicolour models of electroweak symmetry breaking and alike theories this idea is realised in a non-perturbative fashion. The simplest, QCD-like technicolour models predict (super)heavy composite Higgs boson, which is incompatible the LHC discovery of a Higgs-like particle with mass $m_h\approx 125.5$ GeV announced at this conference \cite{:2012gk}. Dimensional transmutation, however, is a generic phenomenon in clasically scale-invariant theories and is realised also in weekly-coupled, perturbative theories \cite{Coleman:1973jx}. In this talk I will be discussing such perturbative Coleman-Weinberg-type models.

Secondly, scale-invariance may be a low-energy remnant symmetry within a fully conformally-invariant fundamental theory. The best known example of such a theory is string theory, for which conformal invariance, alongside with supersymmetry, is an underlying symmetry. In string theory-motivated phenomenological models, it is more common to assume that conformal invariance is broken at a very high energy scale and supersymmetry survives down to  $\sim$TeV-scale in order to stabilize the electroweak scale against radiative corrections. It seems, however, that the announced mass for the Higgs-like particle is somewhat higher than it is predicted by the simplest supersymmetric extensions of the Standard Model. On the other hand, the radiative stability of electroweak scale can be ensured also by classical scale invariance \cite{bardeen}. Therefore, it is conceivable to consider a version of string theory where supersymmetry is broken at a high energy scale (perhaps not admissible at LHC), while the classical scale invariance survives in the low-energy limit.

\section{Scale invariance and naturalness}

There are several theoretical (quantum gravity, unification of forces, electric charge quantization, strong CP problem, etc. ) motivations as well as observational evidences (neutrino masses, dark matter) for new physics that completes the Standard Model at high energies.   Extended theories typically involve different mass scales and the question of whether the electroweak scale is natural arises. A technical aspect of this naturalness problem lies in the fact that, considering the Standard Model as an effective theory valid below a cut-off energy scale $\Lambda$, perturbative corrections to the Higgs boson mass, $m_h$, (and, in fact, to the masses of other Standard Model particles) are of the order of $\Lambda$, i.e., $m_h\sim \Lambda$. Since particle mass is related to spacetime symmetries (an eigenvalue of the quadratic Casimir operator of the Poincar\'e group), an obvious way to solve the problem is to extend relativistic invariance by incorporating a new symmetry that forbids linear dependence of scalar masses on ultraviolet energy scales. The most discussed symmetry that serves this purpose is supersymmetry, another is scale invariance.

To demonstrate that light scalars are natural (in the above technical sense) within scale-invariant theories, let us consider a simple theory of a single scalar field $S(x)$ which is described the following generating functional:
\begin{equation}
Z_{\Lambda}\left[J_S\right]=\int DS \exp\left\lbrace i\int d^4x  \left({\cal L}_{\Lambda}+J_S S\right) \right \rbrace~,
\label{2}
\end{equation}
where ${\cal L}_{\Lambda}$ is an effective Wilsonian Lagrangian
\begin{equation}
{\cal L}_{\Lambda}=\frac{1}{2}\partial_{\mu}S\partial^{\mu}S-\frac{1}{2}m^2(\Lambda)S^2-\frac{\lambda(\Lambda)}{4}S^4+...
\label{3}
\end{equation}
 where ... denote infinite series of possible terms of mass dimension higher than 4, which are irrelevant at low-energies. In this effective theory scale invariance is badly broken both by the mass term and cut-off $\Lambda$. The effective Wilsonian theory is finite and no further regularization is required. One-loop quantum correction to the $\Lambda$-dependent bare mass in (\ref{3}) can be easily computed:
 \begin{equation}
 m_{\rm R}^{2}(\mu)=m^2(\Lambda)+ \frac{3\lambda(\Lambda)}{16\pi^2}\left[\Lambda^2-m^2(\Lambda)\ln\left( \Lambda^2/\mu^2\right) \right]~.
\label{4}
 \end{equation}
Thus, in order to have light scalar, $m_{\rm R}<<\Lambda$, an unnatural fine-tuning between the bare mass and the cut-off scale must be assumed. This is the above-mentioned naturalness (mass hierarchy) problem.

Suppose now that the effective theory described by (\ref{2}) is embedded in an underlying theory that contains heavy fields/field modes $H(x)$, which is classically scale invariant.
That is, the generating functional of this `fundamental' theory,
 \begin{equation}
Z\left[J_S, J_H\right]=\int \left[DSDH\right] \exp\left\lbrace i\int d^4x  \left({\cal L}[S,H]+J_S S+J_HH\right) \right \rbrace~,
\label{5}
\end{equation}
is not scale-invariant due to the non-invariance of the functional measure $\left[DSDH\right]$, while the action $\int d^4x {\cal}L[S,H]$ maintains the scale invariance, i.e.:
\begin{equation}
{\cal L}[tS,tH]=t^4{\cal L}[S,H]~.
\label{6}
\end{equation}
Since the bare Wilsonian effective action $\int d^4x {\cal L}_{\Lambda}[S]$ in (\ref{2}) results from integrating out $H$-fields in (\ref{5}), i.e., $\exp\left\lbrace i \int d^4x {\cal L}_{\Lambda}[S] \right\rbrace=\int DH  \exp\left\lbrace i\int d^4x  \left({\cal L}[S,H]\right) \right \rbrace $, the condition
\begin{equation}
m^2_{\rm R}(\Lambda)\equiv m^2(\Lambda)+\frac{3\lambda(\Lambda)}{16\pi^2}\Lambda^2=0~,
\label{7}
\end{equation}
 is a natural renormalization condition that is forced  due to the classical scale invariance of the underlying theory, see eq. (\ref{6}). Thus we are left only with a mild logarithmic dependence of the scalar mass in (\ref{4}) on the cut-off scale. Other related discussion on the absence of quadratic divergences can be found in \cite{bardeen}, \cite{Foot:2007as}, \cite{Meissner:2007xv}, \cite{Aoki:2012xs}.

\section{Scale-invariant models}

It is well-known that the scale-invariant version of the Standard Model is not phenomenologically viable since it predicts a very light ($m_h\lesssim 10$ GeV) Higgs boson and, in addition, requires light top quark $m_t\lesssim 40$ GeV. It has been realised for quite sometime  that phenomenologically successful scale-invariant models require extension of the bosonic sector of the Standard Model \cite{Hempfling:1996ht}. There are different physical motivations to do so, e.g., incorporation of neutrino masses, scalar dark matter, new gauge bosons, etc. In ref. \cite{Foot:2010et}  we have observed that incorporation of small cosmological constant within the scale-invariant models necessarily implies that the mass of the dilaton is generated at two-loop level. For the phenomenologically most interesting models with scale-invariance spontaneously broken at TeV-scale, this means that mass of the dilaton can be $\lesssim 10$ GeV. Some realistic electroweak scale-invariant models along these lines have been discussed in \cite{Foot:2011et}.

The minimal scale-invariant extension of the Standard Model is given by the addition to the Standard Model field content of a singlet scalar field. Demanding cancellation of the cosmological constant, one obtains the following prediction for the dilaton and Higgs mass, respectively
\begin{eqnarray}
m_{\rm dil}\approx 7-10~{\rm GeV}~, \\
 m_{h}=12^{1/4}m_t\approx 300~{\rm GeV}~.
\label{8}
\end{eqnarray}
If the announced LHC resonance is indeed a Higgs boson, the above minimal model is clearly excluded.

An interesting class of electroweak scale-invariant theories are those which incorporate type-II see-saw mechanism for neutrino mass generation. These models contain an extra electroweak triplet scalar field $\Delta$. The mass of this particle is predicted to be \cite{Foot:2011et}:
\begin{equation}
m_{\Delta}=\left(2m_t^4-m_h^4/6\right)^{1/4}\approx 190~{\rm GeV}~,
\label{9}
\end{equation}
for $m_h\approx 125$ GeV. Scale-invariant models with hidden/mirror sector dark matter have also been discussed in \cite{Foot:2011et}.

\section{Conclusion}

Particle physics models with classical scale invariance are attractive in many respects.  In scale-invariant theories all mass scale have a purely  quantum-mechanical origin.  Very simple models with classical scale invariance are capable of resolving the hierarchy problem without introducing supersymmetry and/or other exotics. Some of the predictions of these theories are  testable at LHC and/or future linear collider.

Phenomenologically the most interesting electroweak scale-invariant models with vanishing cosmological constant generically predict a light dilaton. Realistic scale-invariant models require extended bosonic sector with masses correlated with the masses of the Higgs boson and the top quark. Hopefully, some of the features of scale-invariant theories described in this talk will be observed at LHC.

\vspace{1cm}

I am grateful to Robert Foot, Kristian McDonald and Ray Volkas for collaboration. The work was partially supported by the Australian Research Council.

\vspace{1cm}

redpajama_set_name: RedPajamaArXiv


text:
\section{Introduction}
Evidence is mounting for an unseen planet in the Oort cloud \citep{Trujillo_2014, Fuente_Marcos_2014, Fuente_Marcos_2015, Batygin_2016}.  This so-called Planet Nine (a.k.a.\ Planet X, Bowie) betrays its presence by the anisotropic orbits of Sedna-like objects. Astronomers are refining the predictions of Planet Nine's location based on orbital dynamics \citep{Fuente_Marcos_2016}, constraints on gravitational perturbations \citep{Fienga_2016}, and constraints on reflected light \citep{Brown_2016}.

Any object with internal power exceeding incident starlight emits more energy at thermal wavelengths than it reflects.  This is likely the case for the purported Planet Nine, a $\gtrsim$10$M_\oplus$ planet at $\sim$700~AU \citep[][]{Batygin_2016}.  Even if a planet has negligible internal heat, thermal emission exceeds reflected light if its Bond albedo is less than 50\%, as for all Solar System planets other than Venus.  These simple facts motivate the search for Planet Nine's blackbody emission.

In this Letter, we propose that cosmology experiments at mm wavelengths are a viable way to search for---and characterize---Planet Nine.  In Section~2 we estimate the temperature and thermal radiation of Planet Nine and in Section~3 we estimate its motion through the sky. In Section~4 we show that the point-source sensitivity and angular resolution of many mm cosmology experiments are sufficient to detect and identify Planet Nine.  We consider false positives in Section~5, and discuss our findings in Section~6.

\section{Thermal Flux from Planet Nine}
\cite{Batygin_2016} suggest that Planet Nine is Neptune-mass and at a distance of approximately 700~AU.  We use the energy budgets of Neptune and Uranus to place constraints on the internal heat of Planet Nine.

The effective temperature of a planet, $T_{\rm eff}$, is the sum of internal power and absorbed solar power: $T_{\rm eff}^4 = T_{\rm iso}^4 + T_{\rm rad}^4$, where $T_{\rm iso}$ is the temperature the planet would have in isolation, and $T_{\rm rad}= (1-A)^{\frac{1}{4}} T_{\rm equ}$ is the temperature of the planet due solely to stellar radiation ($A$ is Bond albedo and $T_{\rm equ}$ is equilibrium temperature).

We may therefore solve for the temperature a planet would have in the absence of an external energy source:
\begin{equation}
T_{\rm iso} = \left(\frac{\zeta-1}{\zeta} \right)^{\frac{1}{4}} T_{\rm eff},
\end{equation}
where $\zeta$ is the ratio of total energy output to stellar energy input: $\zeta \equiv (T_{\rm eff}/T_{\rm rad})^4$.

Table 1.3 of \cite{Taylor_2010} lists the values of $\zeta$ for Uranus and Neptune: 1.06 and 2.52, respectively.  In other words, Neptune has greater internal heat than Uranus.  \cite{Hildebrand_1985} measured the effective temperatures of Uranus and Neptune: $T_{\rm eff}= 58$~K and 60~K, respectively.
This allows us to estimate the isolation temperature of Uranus and Neptune to be $T_{\rm iso} = 28$~K and 53~K, respectively.   We adopt this as a plausible range for Planet Nine.

We next estimate the solar heating at 700~AU.  We find that Planet Nine would have a temperature of $\sim$10~K, in the absence of internal heat and assuming that Planet Nine has the same Bond albedo as Uranus and Neptune (0.30 and 0.29, respectively).  In any case, even in the zero-albedo limit, the energy budget of Planet Nine is dominated by internal heat.  The exact albedo is therefore irrelevant.

Since $T_{\rm iso}>T_{\rm rad}$ for Planet Nine (and therefore $T_{\rm iso}^4\gg T_{\rm rad}^4$), we expect that its effective temperature will be approximately given by its internal temperature: $T_{\rm eff} \gtrsim T_{\rm iso}$.  We therefore expect Planet Nine to have an effective temperature of $T_{\rm eff} = 28$--53~K, corresponding to a blackbody peak of 50--100~micron. Anyone wishing to detect the thermal emission of this planet had better look in the far-IR and longward.

\cite{Hildebrand_1985} presented far-infrared fluxes of the Solar System giant planets. Their Figure~8 shows the flux of Uranus and Neptune from 30 to 1000~micron.  They found that Uranus had a flux density of $F_\nu= \{200, 70\}$~Jy and Neptune has a flux density of $F_\nu= \{60, 25\}$~Jy at $\lambda= \{500, 1000\}$~micron.

Scaling from Neptune, we obtain a 1~mm flux density of:
\begin{equation}
F_{\rm 300GHz} = 30~{\rm mJy} \left(\frac{T_{\rm eff}}{40~{\rm K}}\right) \left(\frac{R_p}{R_{\rm N}}\right)^2 \left(\frac{d}{700~{\rm AU}}\right) ^{-2},
\end{equation}
where $R_p$ is the radius of Planet Nine, $R_{\rm N}$ is the radius of Neptune, and $d$ is Planet Nine's current heliocentric ($\sim$geocentric) distance. This flux is consistent with the planetary evolution models of \cite{Linder_2016} and \cite{Ginzburg_2016}.

Most cosmology experiments designed to map the cosmic microwave background (CMB) have channels at/near 1~mm (300~GHz), but they are typically optimized for 2~mm radiation (150~GHz).  Since $F_\nu \propto \lambda^{-2}$ on the Rayleigh-Jeans tail, we can extrapolate Planet Nine's flux density to these longer wavelengths:
\begin{equation}
F_{\rm 150GHz} = 8~{\rm mJy} \left(\frac{T_{\rm eff}}{40~{\rm K}}\right) \left(\frac{R_p}{R_{\rm N}}\right)^2 \left(\frac{d}{700~{\rm AU}}\right) ^{-2}.
\end{equation}

\section{Proper Motion of Planet Nine}
Solar System objects move in the sky relative to the Earth due to both their own and Earth's orbital motion.

Neglecting Earth's orbital motion, the proper motion of a Solar System body on a circular orbit is $\mu_{\rm orb} = 2\pi/P$, where $P$ is the object's orbital period.  Using Kepler's Third Law, we get
\begin{equation}
\mu_{\rm orb} = 0.3~{\rm arcsec/day}\left(\frac{d}{700~{\rm AU}} \right)^{-3/2}
\end{equation}
or, in the units of Kuiper Belt aficionados,
\begin{equation}
\mu_{\rm orb} = 10^{-2}~{\rm arcsec/hr}\left(\frac{d}{700~{\rm AU}} \right)^{-3/2}.
\end{equation}
Given the eccentricity range predicted by \cite{Batygin_2016}, the instantaneous orbital motion of planet nine could be a factor of 2--3 greater than this estimate.

Neglecting an object's orbital motion, its annual parallax as seen from Earth is $p = (d/{\rm pc})^{-1}$~arcseconds.  For Planet Nine, this gives a parallax of
\begin{equation}
p = 5~{\rm arcmin}  \left( \frac{d}{700~{\rm AU}}\right)^{-1}.
\end{equation}

The amplitude of the parallax motion seen from Earth is roughly 2$p$ per 6~months, or
\begin{equation}
\mu_{\rm par} = 3.4~{\rm arcsec/day} \left( \frac{d}{700~{\rm AU}}\right)^{-1}
\end{equation}
or
\begin{equation}
\mu_{\rm par} = 10^{-1}~{\rm arcsec/hr}\left(\frac{d}{700~{\rm AU}} \right)^{-1}.
\end{equation}
The instantaneous motion can be somewhat greater than this at opposition. Note, however, that mm telescopes exhibit no preference for opposition since they can operate during the day.  We therefore expect Planet Nine to move a few arcseconds per day, mostly due to parallax.

\section{Detecting Planet Nine with CMB Experiments}
In order to find Planet Nine, a mm telescope needs sufficient point-source sensitivity to detect the object, and sufficient angular resolution to see its parallactic motion. The sensitivity of mm telescopes depends not only on their collecting area, but also the number of detectors. The angular resolution of mm telescopes, on the other hand, is simply related to their diameter, since they are nearly diffraction limited in most cases.  In other words, a small telescope with cutting-edge detectors may be sensitive to the mm flux of Planet Nine, but it is unlikely to resolve its sky motion. We now briefly discuss current and planned CMB experiments (also summarized in Figure~\ref{money_plot}).

Figure~4 of \citep{Fich_2014} suggests that the South Pole Telescope (SPT) and CCAT (with first generation instruments) have a 5$\sigma$ point-source detection limit of 1~mJy at 1~mm for a fiducial 1~hour scan of a 1 square degree field.  At that scan rate, a 1~year cosmological experiment could survey 365$\times$24 = 8760 square degrees, or 20\% of the sky.  This is large enough to put useful constraints on the location and luminosity of Planet Nine.

Alternatively, the same instruments could survey 80\% of the sky in one year down to 2~mJy at 1~mm (in practice SPT can only access half the sky because it is located at the South Pole). The optimal strategy to search for Planet Nine would be to map the maximum area of the sky every few months in order to resolve its parallactic motion. For example, SPT could scan the entire Southern sky down to 4~mJy at 1~mm in under 3 months.

If Planet Nine has a 1~mm flux of 30~mJy, then it should be robustly detected by SPT or CCAT within an hour using the fiducial scan strategy of 1 square degree per hour, or any of the other scan strategies described above. SCUBA-2 \citep{Thompson_2007} and JCMT \citep{Moore_2015}, however, are unlikely to be sensitive to Planet Nine unless it is significantly brighter than we predict.

SPT has an angular resolution of one arcminute \citep{Ruhl_2004}, so the parallax of Planet Nine is resolvable, and its motion is detectable with a few months baseline.  CCAT will have an angular resolution of ten arcseconds \citep{Fich_2014}, so the motion of Planet Nine would be easily detectable within a couple weeks. Planck has a resolution of 5~arcmin \citep{Cremonese_2002}, and can therefore just resolve the parallactic motion of Planet Nine.  The mission's point-source sensitivity, however, is insufficient to robustly detect Planet Nine unless it is hotter and/or closer than expected.

Advanced ACTPol \citep{Henderson_2015} will observe 40\% of the sky and will have both the point source sensitivity and angular resolution to detect and identify Planet Nine.  Finally, CMB-S4 \citep{Abazajian_2013} will observe half of the sky with exquisite sensitivity, but it may only have an angular resolution of 3~arcmin, which would barely allow it to resolve Planet Nine's parallax. As a general rule, mm telescopes are nearly diffraction limited, so a 5~m telescope has $\lesssim$1~arcmin resolution at 1~mm; this is the smallest telescope that could convincingly resolve the parallax of Planet Nine if it is at 700~AU. The 12~arcmin beam of BICEP \citep{Ade_2015}, for example, is unable to resolve the parallax motion of Planet Nine.

\begin{figure*}[htb]
\begin{center}
      \includegraphics[width=120mm]{money_plot_2.pdf}
    \caption{The search space for Planet Nine at mm wavelengths.  There is no other object with both parallax and mm flux within an order of magnitude of Planet Nine.  Stars are mostly fainter and all move much slower.  There are roughly 4000 asteroids at least as bright as Planet Nine, and many more that are fainter, but they all move orders of magnitude faster.  The largest Trans-Neptunian Objects, Pluto and Eris, are most similar to Planet Nine, but still move $>$10$\times$ faster.  The blue lines show the angular resolution and point-source sensitivity of current and planned CMB experiments.  The sensitivities are approximate because of the unknown scan strategy of future experiments, and because the reported numbers are for stationary sources: a moving source is less susceptible to confusion noise, but will require difference imaging.  Only experiments that can detect Planet Nine's mm flux and resolve its parallax motion can hope to discover it---this rules out Planck, unless Planet Nine is closer and brighter than we have assumed here.  Of the cosmology experiments that could identify the planet, those with greater sky coverage are better bets.  This tips the scales in favor of Advanced ACTPol or CMB-S4. \label{money_plot}}
    \end{center}
\end{figure*}

\section{False Positives}
Planet Nine will not be the only moving mm source. Here we consider false positives, some of which may be worthy targets in their own right.  Researchers have previously studied the use of CMB experiments to study Solar System objects, and the non-Gaussian noise that unresolved bodies contribute to mm maps \citep{Babich_2007, Babich_2009, Maris_2009}.

\subsection{Stars}
Even the stars with the greatest proper motions move orders of magnitude slower than Planet Nine: $\lesssim$0.03~arcsec/day. Stars are also radio quiet and hence much fainter than Planet Nine: there are only a handful of stars in the 10--100~mJy range at 1~mm \citep{Muller_2014}.

\subsection{Asteroids}
\cite{Redman_1995} reported the mm flux densities of asteroids at opposition. They find that the cumulative distribution of asteroids is roughly a power law down to 30~mJy, and they estimate that there are $4000$ asteroids at least this bright after accounting for incompleteness (a factor of 2).  By comparison with \cite{Terai_2013}, this roughly corresponds to all asteroids greater than 20~km in diameter.

Most asteroids at any given moment are not at opposition. Since the nightsides of asteroids are cooler, they are less bright than the daysides. But mm-radiation originates from relatively deep in the asteroid, where the diurnal cycle is attenuated, so asteroids should have muted phase variations at these wavelengths. Asteroids at opposition are roughly a factor of two closer to Earth, however, and hence a factor of four brighter than asteroids that are 60$^\circ$ from the Sun.  Applying this factor to the steep power law in size and mm brightness, we expect 10$\times$ fewer asteroids at sky locations 60$^\circ$ from the Sun than at opposition.

We treat the asteroid belt as being 10$^\circ$ wide and centered on the ecliptic \citep[a rough approximation of][]{Minton_2010}.  The asteroid belt therefore occupies approximately $360\times10= 3600$ square degrees on the sky.  At the SPT angular resolution of 1~arcmin, the asteroid belt therefore spans $60\times60\times3600 \approx 1.3\times10^7$ resolution elements.  We therefore expect an $F_\nu > 30$~mJy asteroid per every few thousand resolution elements, so we are far from the confusion-limit.

For the most distant Main Belt asteroids \citep[4~AU;][]{Minton_2010}, we estimate proper motions of $\mu_{\rm orb} \approx 12$~arcmin/day and $\mu_{\rm par} \approx 10$~arcmin/day.  In other words, although there will be thousands of asteroids comparably bright to Planet Nine at mm wavelengths, they will move orders of magnitude faster, and their motion in the sky will be roughly equal parts proper motion and parallax.

\subsection{Trans-Neptunian Objects}
Like Planet Nine, Trans-Neptunian Objects (TNOs) move slowly, mostly due to parallax.  Fortunately, only the very brightest of these objects would have comparable mm flux to Planet Nine: Pluto is 10--20~mJy \citep{Greaves_2015}, while Eris is about 1~mJy \citep{Bertoldi_2006}. But even here, TNOs betray themselves by their relatively large parallax.  The only way for an object to have Planet Nine's parallax is for it to be at Planet Nine's distance from the Sun.  And at those distances, only a bona fide planet could have a 1~mm flux density of 30~mJy.  In other words, the only false positive for Planet Nine is \emph{another} planet, which would be an equally exciting discovery.

\section{Discussion}
In order to detect Planet Nine, one needs to construct a mm map of the sky every few months and search for a $\sim$30~mJy source with annual parallax of $\sim$5~arcmin (Figure~\ref{money_plot}).  This is within reach of many current and planned cosmology experiments, provided they map large swaths of the sky and return to the same regions every few months.

This is a mode of operation that is easily accessible---but not guaranteed---for cosmology surveys. Some experiments will tile smaller fields and then move on, coming back once in a while or never.   Others survey large swaths over and over---this strategy is ideally suited for piggybacking Solar System science. Even in the worst case scenario of one and done, it takes a month or two to get to survey depth on a few hundred square degrees, during which the parallax motion of Planet Nine might betray its presence in a datacube.

Difference imaging at mm wavelengths will reveal stationary flickering AGN \citep[e.g.,][]{Hobbs_1977} and rapidly moving asteroids. Critically, asteroids move orders of magnitude faster on the sky than Planet Nine.  If Planet Nine is as radio bright as we predict, then there are sufficiently few foreground asteroids that they can be spatially resolved and tracked with hourly observations. Note that CMB experiments need not ``discover'' these asteroids: any object as mm-bright as 30~mJy is already in the Minor Planet Center (MPC) archive.\footnote{www.minorplanetcenter.net/iau/lists/ArchiveStatistics.html}

Current and planned cosmology experiments could observe thousands to millions of asteroids at mm wavelengths, constraining their vertical and horizontal temperature gradients, shapes, and albedo markings.  These experiments could also discover Planet Nine and measure its effective temperature.  Using Uranus and Neptune to bracket plausible internal heat fluxes, we expect Planet Nine to have an effective temperature of 30--50~K. Although we have adopted 40~K as our fiducial value throughout this Letter, our conclusions are essentially unchanged for any temperature in that range.  More importantly, estimating the actual effective temperature of Planet Nine by measuring its mm flux would help inform interior models of Uranus, Neptune, and possibly sub-Neptune exoplanets.

\subsection{What if Planet Nine is Fainter?}
If our assumptions about Planet Nine are wildly off then it could be much fainter. Lack of any internal heat source would lead to an effective temperature of $\sim$10~K if it has a similar Bond albedo to Uranus and Neptune, while a rocky planet of the same mass would be smaller.  Lastly, if Planet Nine lies at 1000 instead of 700~AU, then its flux would be halved.  If we conservatively adopt all of these scenarios, then the 1~mm flux density of Planet Nine could be as low as 1~mJy. This pessimistic scenario does not change our basic conclusion: Planet Nine---and its parallax motion---would still be detectable with many mm telescopes. Nonetheless, the search would be more challenging.

If Planet Nine is only a few mJy then there are many more false positives.  The brightness distribution of asteroids is steep: extrapolating the power law of \cite{Redman_1995} down to flux densities of 1~mJy yields $10^6$ asteroids.   About $7\times10^{5}$ of these asteroids are already in the MPC archives, of which roughly $5\times10^5$ have sufficiently precise ephemerides to predict in which arcmin$\times$arcmin resolution element they lie \citep{Tedesco_2002, Bottke_2005}. The Large Synoptic Sky Survey (LSST) will find and track all of these---and many more---within half a decade of operations \citep{Jones_2016}. Nonetheless, near-term CMB experiments with sensitivities of 1~mJy have the potential to significantly increase the number of known asteroids.

If we adopt the same 10 degree wide asteroid belt as above, we find that next-generation mm telescopes with 10~arcsec resolution will not be confusion limited, even for a fainter Planet Nine. Current telescopes such as SPT, however, are at the confusion limit near the ecliptic plane: 1 asteroid for every 10 resolution elements \citep{Franceschini_1982}.    In this more daunting regime, there are two approaches for dealing with foreground asteroids: 1) treat them as noise that can be averaged over, or 2) identify and remove asteroids from the map.  We consider each strategy in turn.

If $\sim$20~days of mm observations are combined into a single mm-map, the asteroids will have moved a few degrees over the course of the observations.  They will therefore appear as streaks in the map.  Subtracting two such maps with millions of streaks each will therefore produce a proverbial haystack of streaks (half negative and half positive) with millions of apparent point sources at the intersections of streaks (since no two asteroids have exactly the same mm flux).  This haystack is further complicated by the fact that the flux from individual asteroids varies in time due to rotational effects (non-spherical shape and non-uniform albedo).  It is not obvious how to identify a faint slowly-moving point-source in this mess.

We instead advocate to search for Planet Nine in the RA-DEC-time data cube. Asteroids will move from one resolution element to the next in approximately 2~hours.  By stacking images taken within an hour of each other, and comparing this hourly average to the following hourly average, asteroids can be identified.  Planet Nine will not move appreciably on this timescale, but will show up when comparing images taken weeks--months apart.  This seems like a surmountable problem of matched filtering that should be studied in more detail.

If Planet Nine is smaller and/or more distant than we have assumed, then it will be harder to detect with mm telescopes, but doubly so at optical and near-IR wavelengths, given the $1/d^4$ dependence of reflected sunlight.  A search at mm wavelengths is only hopeless if Planet Nine has no internal heat \emph{and} has an albedo much greater than 50\%; such an ice-ball would be easier to discover via an optical search. Fortunately, cosmologists have designed optical all-sky surveys that may be up to the task \citep{LSST_2009}.

\acknowledgements
We thank Jonathan Fortney, \'Etienne Artigau, Tracy Webb, Matt Dobbs, and Alan Fitzsimmons for timely and constructive feedback.