\section{\label{intr}Introduction}

Almost two decades have passed from the first evidence for a late-time acceleration of the expansion of the Universe based on type Ia supernovae observations \cite{Perlmutter:1998np,Riess:1998cb}. Since then, more precise cosmological data have confirmed these first results and provided overwhelming evidence for such an acceleration \cite{Suzuki:2011hu,Parkinson:2012vd,Hinshaw:2012aka,Anderson:2012sa,Ade:2015xua}. Most cosmological observations are remarkably consistent with a six parameter spatially flat $\Lambda {\rm CDM}$ model where the late-time dynamics of the Universe is dominated by a cosmological constant $\Lambda$ accounting for nearly $70 \%$ of the total energy density of the Universe.

Despite its successes on the observational side, the $\Lambda {\rm CDM}$ model is faced with yet unsolved fundamental challenges, in particular regarding the small magnitude of $\Lambda$ and the coincidence of the era where it becomes dynamically relevant with the present epoch (see, e. g., \cite{Barreira:2011qi}). Hence, a dynamical scalar field violating the strong energy condition offers a better motivated alternative to the cosmological constant $\Lambda$ in the attempt to explain the current acceleration of the Universe \cite{Copeland:2006wr,Frieman:2008sn,Caldwell:2009ix,Li:2011sd,Bamba:2012cp}. On the other hand, it is reasonable to expect that such field could couple to other fields, possibly leading to measurable variations of nature's fundamental couplings \cite{Carroll:1998zi,Chiba:2001er,Wetterich:2002ic,Nunes:2003ff,Anchordoqui:2003ij,Copeland:2003cv,Parkinson,Doran:2004ek,Marra:2005yt,Avelino:2008dc,Avelino:2009fd,Dent:2008vd,Avelino:2011dh,Avelino:2014xsa}.

One of such couplings is the fine-structure constant whose dynamics over a wide redshift range is severely constrained using both cosmological, astrophysical and terrestrial data, as well as local laboratory experiments (see \cite{Uzan:2010pm} for a recent review). Despite a few positive claims for a detection of a spacetime variation of the fine-structure constant \cite{Webb:2000mn,Murphy:2003hw,King:2012id}, there is presently no unambiguous  evidence for such variation. Moreover, low redshifts laboratory experiments \cite{Rosenband:2008} (see also \cite{Luo:2011cf} for a recent review of atomic clock constraints on the variation of fundamental couplings) and the Oklo natural nuclear reactor \cite{Gould:2007au,Petrov:2005pu} provide stringent limits on the time variation of $\alpha$.

In \cite{Martins:2015ama,Martins:2015jta,Martins:2016oyv} it has been suggested that upper limits on the time variation of the fine-structure constant tightly constrain the coupling of a dark energy scalar field to the electromagnetic sector, and, indirectly, the violation of the weak equivalence principle. Furthermore, it has also been suggested that such limits, in combination with standard methods, could be used to improve the constraints on the equation of state of dark energy. In this paper we shall demonstrate that such constraints rely on the specification of priors consistent but not favoured by current data. We will  further show that, even in the case of an idealised model where the DE is described by a dynamical scalar field with a constant equation of state parameter $w$ linearly coupled to the electromagnetic sector through a dimensionless coupling $\zeta$, when such priors are relaxed only the combination $|\zeta| {\sqrt {|w+1|}}$ is tightly constrained by current upper limits on the time variation of the fine-structure constant.

Throughout this paper we shall use units with $\hbar=c=8 \pi G=1$ and a metric signature $(+,-,-,-)$.

\section{Dark energy and varying couplings \label{sec2}}

Here we shall assume that the late-time acceleration of the Universe is due to a dynamical dark energy scalar field $\phi$ non-minimally coupled to the electromagnetic field. A particularly interesting class of dark energy models may be defined by the action
\begin{equation}
S=\int d^4x\,\sqrt{-g}\mathcal\,{\cal L}\,, \label{S1}
\end{equation}
where $g$ is the determinant of the metric tensor,
\begin{eqnarray}
{\cal L} &=& {\cal L}_\phi + {\cal L}_{\phi F} + {\cal L}_{\rm other}\,, \label{S2}\\
{\cal L}_\phi&=&\pm X-V(\phi)\,, \label{S2a}\\
X&=&\pm \frac{1}{2}\partial^\mu \phi \partial_\mu \phi\,, \label{S3}
\end{eqnarray}
$V(\phi)$ is the scalar field potential,
\begin{equation}
{\cal L}_{\phi F}= -\frac{1}{4} B_F (\phi)F_{\mu \nu}F^{\mu \nu}\,, \label{S4}
\end{equation}
$B_F(\phi)$ is a gauge kinetic function, $F_{\mu \nu}$ are the components of
the electromagnetic field tensor, and ${\cal L}_{\rm other}$ is the Lagrangian density of the other fields. In this class of models the fine-structure constant is given by
\begin{equation}
\alpha(\phi)=\frac{\alpha_0}{B_F(\phi)}\,,
\end{equation}
where the subscript `$0$' denotes the present time ($B_F(\phi_0)=1$ today).

\subsection{Time variation of the fine-structure constant}

In the family of models described by Eqs. (\ref{S1}-\ref{S4}) the evolution of $\phi$ induced solely by its coupling to electromagnetically interacting matter is so small, given weak equivalence principle constraints (see  \cite{Uzan:2010pm} and references there in), that the resulting time variation of $\alpha$ can be neglected. Hence, we shall assume that the dynamics of $\phi$ is fully driven by the scalar field potential $V(\phi)$ (and damped by the expansion). On the other hand, since the sound speed of the scalar field $\phi$ equals the speed of light, the spatial variations of the scalar field $\phi$ are small and their contribution to the variation of $\alpha$ may also be neglected in this context.

Consider a flat homogeneous and isotropic Friedmann-Robertson-Walker universe whose dynamics obeys the Friedmann equation given by
\begin{equation}
H^2 = \frac{\rho}{3}=\frac{\rho_\phi+\rho_{[\phi F+{\rm other}]}}{3}\,, \label{fried}
\end{equation}
where $\rho_\phi$ is the dark energy density (associated with ${\cal L}_{\phi}$), $\rho_{[\phi F+{\rm other}]}$ is the energy density associated with the remaining lagrangian components (${\cal L}_{\phi F}$ and $ {\cal L}_{\rm other}$), $H={\dot R}/R$ is the Hubble parameter, $R$ is the scale factor, and a dot represents a derivative with respect to the  physical time, $t$.

Taking into account that the energy density and pressure associated with the scalar field $\phi$ are given respectively by
\begin{equation}
\rho_\phi=\pm {\dot \phi}^2/2+V(\phi) \,, \quad  p_\phi=\pm {\dot \phi}^2/2-V(\phi)\,, \label{rhop}
\end{equation}
one obtains
\begin{equation}
w \equiv \frac{p_\phi}{\rho_\phi}= -1 \pm  \frac{{\phi'}^2H^2}{\rho_\phi}=-1 \pm \frac{{\phi'}^2}{3\Omega_\phi}\,, \label{wphi}
\end{equation}
where $\Omega_\phi=\rho_\phi/\rho$ and a prime represents a derivative with respect to $\ln R$ (${\dot \phi}=\phi' H$).

If the gauge kinetic function $B_F$ is a linear function of $\phi$ with $|B_F(\phi)-1| \ll 1$ then
\begin{equation}
\frac{\dot \alpha}{\alpha}=\zeta {\dot \phi}\,, \label{linear}
\end{equation}
where $\zeta$ is a constant. In the following we shall only consider the solution with $\dot \phi > 0$, so that the sign of $\dot \alpha$ is the same as that of $\zeta$.  Note, however, that this assumption may be relaxed since both  $ {\cal L}_{\phi}$ and ${\cal L}_{\phi F}$ are invariant under the transformation $\phi \to -\phi$, $V(\phi) \to V(-\phi)$, $\zeta \to -\zeta$.

Eqs. (\ref{wphi}) and (\ref{linear}) imply that
\begin{equation}
\frac{1}{H}\frac{\dot \alpha}{\alpha} =\zeta \sqrt{3\Omega_\phi(z)|1+w|}\,. \label{alphadot}
\end{equation}

\subsubsection{Constant $w$ models}

The time evolution of the energy density associated with the dark energy scalar field $\phi$ obeys the equation
\begin{equation}
{\dot \rho_\phi}+3H(1+w)=0\,. \label{adiabatic}
\end{equation}
In this paper, for simplicity, we shall consider a constant equation of state parameter $w$ smaller than unity (see \cite{Avelino:2009ze,Avelino:2011ey} for a discussion of constant $w$ models) --- relaxing this assumption would only strengthen our conclusions. Then, from Eqs. (\ref{rhop}), (\ref{wphi}) and (\ref{adiabatic}), it is possible to show that
\begin{equation}
{\dot \phi}=\left(2V\frac{|1+w|}{1-w}\right)^{1/2}\,.
\end{equation}
with
\begin{equation}
V=V_0 R^{-3(1+w)}\,,
\end{equation}
where the scale factor at the present time is normalised to unity ($R_0=1$).
Together with Eq. (\ref{linear}) this implies that
\begin{equation}
\frac{\dot \alpha}{\alpha}=\zeta {\dot \phi} \propto (1+z)^{3(1+w)/2}\,, \label{alphazev}
\end{equation}
has a very slow evolution evolution with the redshift $z$ ($1+z=1/R$) for values of $w$ close to $-1$ consistent with the current observations.

\subsection{Constraints on the time variation of the fine-structure constant}

In this paper, for simplicity, we shall focus on the atomic clock constraint on the present variation of $\alpha$ with time of Rosenband {\it et al.}  \cite{Rosenband:2008}
\begin{equation}
\left. \frac{\dot\alpha}{\alpha} \right |_{z=0} =(-1.6\pm2.3)\times10^{-17}\,{\rm yr}^{-1}\,. \label{Rosen}
\end{equation}
which is currently the strongest laboratory constraint on $\alpha$ alone. Eq. (\ref{alphazev}) implies that, for a constant $w$, this local laboratory constraint on the value of $\alpha$ is significantly more constraining than most other astrophysical and cosmological constraints (in particular than the astrophysical constraints on the variation of $\alpha$ considered in \cite{Martins:2015ama,Martins:2015jta,Martins:2016oyv}).

The constraint given in Eq. (\ref{Rosen}) may be rewritten in a dimensionless form as
\begin{equation}
\left. \frac{1}{H_0}\frac{\dot\alpha}{\alpha} \right |_{z=0} =(-2.3\pm3.3)\times10^{-7}\,,
\end{equation}
taking into account that $H_0=100 \, h \, {\rm km \, s^{-1} \, Mpc^{-1}}$, with $h=0.678$ \cite{Ade:2015xua}, and neglecting the small uncertainty on the current value of the Hubble parameter $H_0$. Assuming a value of $\Omega_{\phi 0} = 0.692$ consistent with the Planck 2015 results \cite{Ade:2015xua} and neglecting the relatively small error bar associated with $\Omega_{\phi 0}$, one finally finds (using Eq. (\ref{alphadot}) evaluated at the present time) that
\begin{equation}
x = \zeta {\sqrt {|1+w|}}=(-1.6\pm 2.3)\times10^{-7}\,. \label{limitx}
\end{equation}

Eq. (\ref{limitx}) implies that a nontrivial lower limit on the value of $|1+w|$ ($|1+w| \ge |1+w|_{\rm min}>0$), if it existed, could be used to obtain an upper bound on the value of $|\zeta|$,
\begin{equation}
|\zeta| \le \frac{3.9}{\sqrt {|1+w|_{\rm min}}} \times10^{-7}\,.
\end{equation}
Analogously, Eq. (\ref{limitx}) combined with a nontrivial lower limit on $|\zeta|$ would imply an upper bound on the value of $|1+w|$ associated to constraints on the time evolution of $\alpha$. Nevertheless, there is currently no unambiguous observational evidence favouring a nontrivial lower limit on either $|w+1|$ or $|\zeta|$. As we shall demonstrate in the following section this precludes the use of constraints on the time variation of $\alpha$ to set realistic upper bounds on $|\zeta|$ or $|w+1|$.

Alternatively, as recognised in \cite{Dvali:2001dd,Chiba:2001er}, a nontrivial lower bound on the value of $|x|$ together with an upper bound on the value of $|w+1|$ could be used to obtain a nontrivial lower bound on $|\zeta|$. However, presently there is also no unambiguous nontrivial lower bound on the value of $|x|$ and, consequently, no non-trivial lower bound on $|\zeta|$ from varying $\alpha$ constraints.

Although in \cite{Martins:2015ama,Martins:2015jta,Martins:2016oyv} cosmological constraints on $|w+1|$ from Type Ia supernova data were considered, here we shall simplify the analysis and account for standard constraints on the value of $|w+1|$ by incorporating them in the $|w+1|$ prior (except if stated otherwise, we shall conservatively assume that $|w+1| \le 0.1$ \cite{Ade:2015xua}).

\section{Role of priors \label{sec3}}

We label the random variables associated with the parameters $\theta={\sqrt {|1+w|}} \ge 0$,  $\zeta$, and $x=\theta \zeta$ by $\Theta$, $Z$, and $X=\Theta Z$, respectively. We shall investigate the impact of the prior on the random variable $\Theta$ on the estimation of $|Z|$, assuming, for simplicity, that $\Theta$ and $X$ are independent random variables. This means that the prior on $\Theta$ is the same as the posterior, since, in this case, the probability density function for the variable $\Theta$ is not altered by a measurement of $X$. This assumption allows us to derive analytical expressions for the probability density function of $|Z|$ for various priors of $\Theta$, but does not otherwise affect our main results.

Given that $\Theta$ and $X$ are assumed to be independent, the cumulative distribution function of the random variable $|Z|=|X|/\Theta$ is given by
\begin{eqnarray}
F_{|Z|}(|\zeta|) &\equiv& P(|Z| \le |\zeta|)= \nonumber \\
& &\int_0^\infty \left( \int_0^{\theta |\zeta|}  f_{|X|} (|x|)  f_{\Theta} (\theta)  d |x| \right) d \theta\,, \label{cdis}
\end{eqnarray}
so that the corresponding probability density function is
\begin{equation}
f_{|Z|}(|\zeta|)=\frac{F_{|Z|}(|\zeta|)}{d|\zeta|} = \int_0^\infty \theta f_{|X|} (\theta|\zeta|)  f_{\Theta} (\theta) d \theta\,. \label{fZ}
\end{equation}

We shall consider a probability density function for the random variable $|X|$,
\begin{eqnarray}
f_{|X|}(|x|)&=&\frac{1}{\sqrt {2 \pi} \sigma} \left(\exp \left[-\frac{(|x|-\mu)^2}{2\sigma^2}\right]+\right. \nonumber \\
& &\left. \exp \left[-\frac{(|x|+\mu)^2}{2\sigma^2}\right]\right)\,, \label{fX}
\end{eqnarray}
consistent with the atomic clock constraint on the present variation of $\alpha$ with time of Rosenband {\it et al.}  \cite{Rosenband:2008} discussed in the previous section ($\mu=-1.6 \times 10^{-7}$ and $\sigma=2.3 \times 10^{-7}$).

\begin{figure}
\includegraphics[width=3.4in]{Fig1.pdf}
\caption{The red solid line represents the probability density function $f_{|Z|}(|\zeta|)$ given by Eqs. (\ref{fu1}) and  (\ref{fu2}) with $\mu=-1.6 \times 10^{-7}$, $\sigma=2.3 \times 10^{-7}$ and $b \equiv {\sqrt {|w+1|}_{\rm max}}=10^{-1/2}$. The vertical dotted line is defined by $\zeta_b \equiv \sigma/b$. $f_{|Z|}$ is nearly constant for $|\zeta| < \zeta_b$ and decays roughly proportionally to $\zeta^{-2}$ for $|\zeta| > \zeta_b$.}
\label{Fig1}
\end{figure}

\subsection{Uniform prior}

Let us start by assuming that the probability density function of the variable $\Theta$ is uniform in the interval $[0,b]$ and vanishes outside it, so that
\begin{equation}
f_{\Theta}(\theta)=\frac{1}{b}\,, \qquad 0 \le \theta \le b\,,
\end{equation}
where $b \equiv {\sqrt {|w+1|}_{\rm max}}$.

In this case, the probability density function of the variable $|Z|$ may be computed analytically, using Eqs. (\ref{cdis}) and (\ref{fZ}), and it is given by
\begin{equation}
f_{|Z|}(|\zeta|)=f^+_{|Z|}(|\zeta|)+f^-_{|Z|}(|\zeta|)\,, \label{fu1}
\end{equation}
where
\begin{eqnarray}
f^\pm_{|Z|}(|\zeta|)&=&\frac{1}{{\sqrt {8 \pi}} b |\zeta|^2} \left(2\sigma\left[\exp\left(-\frac{\mu^2}{2\sigma^2}\right)-\right. \right. \nonumber \\
& &\left. \exp\left(-\frac{(\mu \pm b|\zeta|)^2}{2\sigma^2}\right)\right]+  \nonumber \\
& &\left. {\sqrt {2 \pi}} \mu \left( {\rm erf}\left[\frac{\mu}{{\sqrt 2}\sigma}\right]-{\rm erf}\left[\frac{\mu \pm b|\zeta|}{{\sqrt 2}\sigma}\right] \right)\right)\,, \label{fu2}
\end{eqnarray}
and the error function is defined by
\begin{equation}
{\rm erf}(y) \equiv \frac{2}{\sqrt \pi} \int_0^y e^{-u^2} du\,.
\end{equation}

\begin{figure}
\includegraphics[width=3.4in]{Fig2.pdf}
\caption{The value of $\zeta_*$ such that $F_{|Z|}(\zeta_*) \equiv P(|Z| \le \zeta_*) =\chi$ (calculated using Eqs. (\ref{fu1}) and  (\ref{fu2}) with $\mu=-1.6 \times 10^{-7}$, $\sigma=2.3 \times 10^{-7}$) as a function of $b \equiv {\sqrt {|w+1|}_{\rm max}}$, for $\chi=0.99$ (red solid line), $\chi=0.95$ (blue dashed line) and $\chi=0.68$ (black dot-dashed line).}
\label{Fig2}
\end{figure}

Fig. 1 shows $f_{|Z|}$, given by Eqs. (\ref{fu1}) and  (\ref{fu2}) with $\mu=-1.6 \times 10^{-7}$, $\sigma=2.3 \times 10^{-7}$ and $b \equiv {\sqrt {|w+1|}_{\rm max}}=10^{-1/2}$, as a function of $|\zeta|$ (red solid line). The vertical dotted line is defined by $\zeta_b \equiv \sigma/b$. The probability density function $f_{|Z|}$ is nearly constant for $|\zeta| < \zeta_b$ and decays roughly  proportionally to $\zeta^{-2}$ for $|\zeta| > \zeta_b$. Fig. 2 shows the value of $\zeta_*$ such that $F_{|Z|}(\zeta_*) \equiv P(|Z|  \le  \zeta_*)=\chi$ (calculated using Eqs. (\ref{fu1}) and  (\ref{fu2}) with $\mu=-1.6 \times 10^{-7}$, $\sigma=2.3 \times 10^{-7}$) as a function of $b \equiv {\sqrt {|w+1|}_{\rm max}}$, for $\chi=0.99$ (red solid line), $\chi=0.95$ (blue dashed line) and $\chi=0.68$ (black dot-dashed line).  Fig. 2 shows that a weaker prior on $|\Theta|$ leads to stronger constraints on $|Z|$. This is directly associated with the choice of a uniform prior in the interval $[0,b]$ which disfavours very small values of $|\zeta|$, specially if $b$ is large. Also note that, due to the heavy tail of $f_{|Z|}$, the constraints on $|Z|$ degrade very rapidly as one increases the confidence level $\chi$. For $b \equiv {\sqrt {|w+1|}_{\rm max}}=10^{-1/2}$ and $\chi=0.68$ one obtains $\zeta_*=2.3 \times 10^{-6}$, which is in reasonable agreement with the results obtained in \cite{Martins:2015ama,Martins:2015jta,Martins:2016oyv} considering a flat prior for the equation of state parameter of the dark energy.

\subsection{Logarithmic prior}

Let us now consider the case of a uniform probability density function of the variable $\ln \Theta$ for $\theta$ in the interval $[a,b]$, with $0 < a < b$. The corresponding probability density function of the variable $\Theta$ is
\begin{equation}
f_{\Theta}(\theta)=\left[\ln \frac{b}{a}\right]^{-1} \frac{1}{\theta}\,, \qquad a \le \theta \le b\,,
\end{equation}
and it is equal to zero outside the interval $[a,b]$.
In this case, the probability density function of the variable $|Z|$ may also be computed analytically, using Eqs. (\ref{cdis}) and (\ref{fZ}), and it is given by
\begin{eqnarray}
& &f_{|Z|}(|\zeta|)= \left[\ln \frac{b}{a}\right]^{-1} \frac{1}{2|\zeta|} \left({\rm erf}\left[\frac{\mu+b|\zeta|}{{\sqrt 2}\sigma}\right]-\right. \nonumber \\
& &\left. {\rm erf}\left[\frac{\mu-b|\zeta|}{{\sqrt 2}\sigma}\right]-{\rm erf}\left[\frac{\mu+a|\zeta|}{{\sqrt 2}\sigma}\right]+{\rm erf}\left[\frac{\mu-a|\zeta|}{{\sqrt 2}\sigma}\right]\right)\,. \label{fl}
\end{eqnarray}

\begin{figure}
\includegraphics[width=3.4in]{Fig3.pdf}
\caption{The red solid line represents the probability density function $f_{|Z|}(|\zeta|)$ given by Eq. (\ref{fl}) with $\mu=-1.6 \times 10^{-7}$, $\sigma=2.3 \times 10^{-7}$, $a \equiv {\sqrt {|w+1|}_{\rm min}}=10^{-3}$ and $b \equiv {\sqrt {|w+1|}_{\rm max}}=10^{-1/2}$. The vertical dotted lines are defined by $\zeta_b \equiv \sigma/b$ and $\zeta_a \equiv \sigma/a$ (left and right dotted lines, respectively). $f_{|Z|}(|\zeta|)$  is nearly constant for $|\zeta| < \zeta_b$ and decays roughly proportionally to $\zeta^{-1}$ for $\zeta_b < |\zeta| < \zeta_a$ (and much faster than that for $|\zeta| > \zeta_a$).}
\label{Fig3}
\end{figure}

Fig. 3 shows $f_{|Z|}$, given by Eq. (\ref{fl}) with $\mu=-1.6 \times 10^{-7}$, $\sigma=2.3 \times 10^{-7}$, $a \equiv {\sqrt {|w+1|}_{\rm min}}=10^{-3}$ and $b \equiv{\sqrt {|w+1|}_{\rm max}}=10^{-1/2}$, as a function of $|\zeta|$ (red solid line). The vertical dotted lines are defined by $\zeta_b \equiv \sigma/b$ and $\zeta_a \equiv \sigma/a$ (left and right dotted lines, respectively). The probability density function is nearly constant for $|\zeta| < \zeta_b$ and decays roughly proportionally to $\zeta^{-1}$ for $\zeta_ b < |\zeta| <  \zeta_a$ (and much faster than that for $|\zeta| > \zeta_a$). Fig. 4 shows the value of $\zeta_*$ such that $F_{|Z|}(\zeta_*) \equiv P(|Z| \le |\zeta_*|) =\chi$ (calculated using Eq. (\ref{fl}) with $\mu=-1.6 \times 10^{-7}$, $\sigma=2.3 \times 10^{-7}$ and $b \equiv {\sqrt {|w+1|}_{\rm max}}=10^{-1/2}$), as a function of $a$, for $\chi=0.99$ (red solid line), $\chi=0.95$ (blue dashed line) and $\chi=0.68$ (black dot-dashed line).  Fig. 2 shows that the lower the value of $a$, the weaker the constraints on $|Z|$ become. For $a \equiv {\sqrt {|w+1|}_{\rm min}} = 10^{-2}$ and $\chi=0.68$ one obtains $\zeta_*=6.0 \times 10^{-6}$, in reasonable agreement with the results obtained in \cite{Martins:2016oyv} using a logarithmic prior for the equation of state parameter of the dark energy. However, in the $a \equiv {\sqrt {|w+1|}_{\rm min}} \to 0$ limit the logarithmic prior favours values of $|w+1|$ extremely close to zero, which is the reason why, in this limit, the constraints on $|Z|$ become extremely weak. Again, note the rapid degradation of the constraints on $|Z|$ with the increase of the confidence level $\chi$.

\begin{figure}
\includegraphics[width=3.4in]{Fig4.pdf}
\caption{The value of $\zeta_*$ such that $F_{|Z|}(\zeta_*) \equiv P(|Z| \le \zeta_*) =\chi$ (calculated using Eq. (\ref{fl}) with $\mu=-1.6 \times 10^{-7}$, $\sigma=2.3 \times 10^{-7}$ and $b \equiv {\sqrt {|w+1|}_{\rm max}}=10^{-1/2}$), as a function of $a \equiv {\sqrt {|w+1|}_{\rm min}}$, for $\chi=0.99$ (red solid line), $\chi=0.95$ (blue dashed line) and $\chi=0.68$ (black dot-dashed line).}
\label{Fig4}
\end{figure}

\begin{figure}
\includegraphics[width=3.4in]{Fig5.pdf}
\caption{The value of $\zeta_*$ such that $F_{|Z|}(\zeta_*) \equiv P(|Z| \le \zeta_*) =\chi$ (calculated using Eq. (\ref{fZpower}) with $\sigma=2.3 \times 10^{-7}$ and $b \equiv {\sqrt {|w+1|}_{\rm max}}=10^{-1/2}$), as a function of $\beta$, for $\chi=0.99$ (red solid line), $\chi=0.95$ (blue dashed line) and $\chi=0.68$ (black dot-dashed line).}
\label{Fig5}
\end{figure}

\subsection{Power law prior}

A more general class of probability density functions, which includes uniform probability density functions for the variable $\Theta$ as a particular sub-class and uniform probability density functions for the variable $\ln \Theta$ as a special limit for $\beta \to -1$, is given by
\begin{equation}
f_{\Theta}(\theta)=\frac{\beta+1}{b} \left(\frac{\theta}{b}\right)^\beta \,, \qquad \theta \le b\,,
\end{equation}
with $f_{\Theta}(\theta)=0$ for $\theta > b$ (here, $\beta > -1$). Taking the probability density function for the variable $X$ given in Eq. (\ref{fX}), but now assuming $\mu=0$, $f_{|Z|}(|\zeta|)$ may be computed analytically using Eqs. (\ref{cdis}) and (\ref{fZ}). The result is
\begin{eqnarray}
f_{|Z|}(|\zeta|)&=& \frac{b(\beta+1)}{\sqrt {2 \pi} \sigma} \left(\frac{b^2 \zeta^2}{2 \sigma^2}\right)^{-\beta/2-1} \times \nonumber \\
& &\left( \Gamma\left[\frac{\beta}{2}+1\right]-\Gamma\left[\frac{\beta}{2}+1,\frac{b^2 \zeta^2}{2 \sigma^2}\right]\right)\,, \label{fZpower}
\end{eqnarray}
where
\begin{equation}
\Gamma(y) \equiv \int_0^\infty u^{y-1} e^{-u} du
\end{equation}
is the Gamma function and $\Gamma(y,z)$ is the upper incomplete Gamma function defined by
\begin{equation}
\Gamma(y,z) \equiv \int_z^\infty u^{y-1} e^{-u} du\,.
\end{equation}
The value of $\mu$ implied by Eq. (\ref{limitx}) is less than one sigma away from zero. Consequently, the error committed in assuming that $\mu=0$ is relatively small, thus justifying the use of this approximation in order to obtain the analytical result for $f_{|Z|}(|\zeta|)$ given by Eq. (\ref{fZpower}).

Fig. 5 shows the value of $\zeta_*$ such that $F_{|Z|}(\zeta_*) \equiv P(|Z| \le \zeta_*) =\chi$ (calculated using Eq. (\ref{fZpower}) with $\sigma=2.3 \times 10^{-7}$ and $b \equiv {\sqrt {|w+1|}_{\rm max}}=10^{-1/2}$), as a function of $\beta$, for $\chi=0.99$ (red solid line), $\chi=0.95$ (blue dashed line) and $\chi=0.68$ (black dot-dashed line). It shows that for values of $\beta$ sufficiently close to $-1$ the constraints on $|Z|$ become extremely weak. Notice that the the $68\%$, $95\%$ and $99\%$ constraints on $|Z|$, may span several orders of magnitude, in particular for $\beta$ close to $-1$.

\section{\label{conc} Conclusions}

In this paper we critically assessed recent claims suggesting that upper limits on the time variation of $\alpha$ could be used to tightly constrain the dynamics a dark energy scalar field (in particular, its coupling to the electromagnetic sector). We have shown that such constraints rely on assumptions which are consistent but not favoured by current data. This situation could be improved if there was i) a nontrivial lower bound on the value of $|x|$ or ii) a nontrivial lower bound on the value of $|w+1|$ or $|\zeta|$. Although i) may in principle be accomplished with a new generation of high-resolution ultra-stable spectrographs, such as ESPRESSO and ELT-HIRES, and ii) may, in principle, be achieved respectively by forthcoming missions to map the geometry of the Universe, such as Euclid, or to test the equivalence principle, such as MICROSCOPE or STEP, there is a priori no guarantee that these missions will make a detection rather than significantly improving current bounds. In the later case the analysis reported in the present paper will remain pertinent, despite the improved constraints. As demonstrated in \cite{Avelino:2008dc,Avelino:2011dh,Avelino:2014xsa}, even in the more optimistic case of a significant detection, the dependence of dark energy constraints from the time variation of $\alpha$ on crucial assumptions, including i) that general relativity provides an accurate description of gravity on cosmological scales ii) that dark energy may be described by a dynamical scalar field obeying Eqs. (\ref{S1}-\ref{S4})  iii) that the coupling between the dark energy scalar field and the electromagnetic field is linear, should not neglected.

\begin{acknowledgments}

The author thanks Margarida Cunha for enlightening discussions. This work was supported by Funda{\c c}\~ao para a Ci\^encia e a Tecnologia (FCT) through the Investigador FCT contract of reference IF/00863/2012 and POPH/FSE (EC) by FEDER funding through the program Programa Operacional de Factores de Competitividade - COMPETE. Funding of this work was also provided by the FCT Grant No. UID/FIS/04434/2013.
\end{acknowledgments}


\section{Introduction

Magnetic fields generated at the base of the Sun's convective zone are thought to form toroidal flux tubes that become buoyant and rise to the surface \citep{par55, cha97, cha05}.  Sunspots often appear where the flux loops break the surface.  On average, bipolar sunspots show leading spots to be closer to the Equator than following spots.  \citet{hal19} first published observations of this phenomenon, now known as Joy's law, after statistical analysis showed that the mean tilt angle of bipolar sunspots increased with latitude. Joy's law has traditionally been interpreted as the Coriolis force operating on divergent plasma at the apex of a rising magnetic flux tube.  Rough calculations of the Coriolis effect on a rising flux tube by \citet{dsi93} measured a deflection in the tilt angle from an E-W orientation over time in terms of the rotational frequency of the Sun and emergence latitude.  The subsurface pitch angle of the toroidal field has also been proposed as a cause of tilt angle prior to the rise of flux tubes through the convection zone \citep{bab61}.

Numerical simulations of toroidal flux loops \citep{fan94} show a dependence of the tilt angle on the emerging latitude and field strength (B) where B $>$ 20 kG. Negative tilt angles start to occur when B $<$ 20 kG, where the weak field strength and low flux ($10^{20}$ Mx) host a converging flow at the apex in contrast to the standard divergent flow model \citep{web13}.  The thin flux tube approximation has been used to study rising magnetic loops in the convection zone and to explore the tilt angles and latitudes of emergence \citep{spr81, mor86}.  Studies show that the strength of the toroidal magnetic field should be around 30 - 100 kG in addition to indicating that the Coriolis force could explain the tilt angle described by Joy's law \citep{dsi93, cal95}. See \citet{fan09} for a review of these topics, including models of toroidal flux tubes rising with and without the influence of convective zone turbulence. \citet{web13} also compared flux tube simulations and found that tilt has a dependence on magnetic field strength in the flux tubes but no dependence on the total flux in the tube.

Thin flux tube models that include convection and radiative diffusion predict shorter rise times through the convective zone \citep{web15}, producing tilt angles consistent with observed active regions. These simulations are able to reproduce tilt angles consistent with Joy's law without the need for anchored footpoints in the overshoot region.  However, tilt angle scatter for magnetic field strengths $\le 40$ kG are higher than observations.

\citet{fis95} found that tilt was proportional to latitude and magnetic flux in the tube and inversely related to the magnetic field strength in the initial toroidal flux tube at the base of the convection zone. This version of Joy's law includes flux and initial field strength because aerodynamic drag balances the Coriolis force for rapidly emerging regions and magnetic tension balances the Coriolis force for slowly emerging regions. However, \citet{kos08} studied MDI magnetograms and found no evidence of a Joy's law dependence on magnetic flux or for a relaxation of the tilt angle toward zero after emergence.  Introducing magnetic field line twist stabilizes the cohesion of the rising tubes but also affects the tilt angle \citep{fan08}.  Since sunspot area and sunspot flux are highly correlated, area is used as a proxy for flux in order to study the relationship between the tilt angle and flux. \citet{jia14} binned Kodaikanal and Mt. Wilson Observatory tilt angle data according to sunspot size and found a weak correlation between mean tilt angle and sunspot group size, while the standard deviations significantly decrease with sunspot group size.

The amount of scatter in the tilt angles of bipolar regions provides information about the magnetic structures that produce sunspots.   \citet{wan89} analyzed 2710 bipolar magnetic regions (BMRs) and found no noticeable dependence of tilt angle on flux but higher deviations from the mean tilt angle for weaker BMRs. \citet{lon02} explained the departures of tilt angle from Joy's law during emergence as caused by upper convective zone turbulence.  Tilt angle scatter about Joy's law introduced by convection increases as flux decreases in thin flux tube models by \citet{web13}. \citet{ill15} showed the significant scatter in bipolar regions with areas less than 300 MSH (including ephemeral regions without sunspot activity).  A transition occurs between 300 and 400 MSH where the distribution of larger regions becomes dominated by sunspot activity with substantially less tilt angle scatter and more tilt angles that follow Joy's law.

 It has been proposed that the magnetic field could become so weak below the surface that convective zone turbulence dominates the flux tube.  \citet{fan94} suggested a mechanism of dynamic disconnection where a submerged portion of an emerged flux tube collapses after achieving hydrostatic equilibrium, disconnecting the tilt angle from the influence of the initial toroidal field pitch angle.  \citet{sch99} proposed subsurface reconnection of the untethered legs of the flux tube at depth to explain surface activity during decay.  \citet{lon02} rule out dynamic disconnection at shallow depths in their model of rising flux tubes.  \citet{sch05} modify the dynamic disconnection model based on strong, buoyancy-driven upflow and radiative cooling in a rising flux loop prior to emergence, finding disconnection depths around 5 Mm and disconnection times less than 3 days.

As a flux loop ends its emergence and the footpoints stop separating, as observed in the photosphere as the centroid locations of the two polarities, Coriolis forces should end and the tilt angles should relax toward zero due to magnetic tension restoring field lines to an E-W orientation.  If the field is significantly frozen into the plasma and the differential rotation force is not strong enough, then the region maintains the tilt angle established prior to the end of emergence. The higher cadence of recently collected sunspot data presents a more complete picture of the tilt angle and footpoint separation over the lifetime of active regions.  A limited study of six active regions by \citet{pev03} using Solar and Heliospheric Observatory (SOHO)-MDI data reported separation distances around 25 Mm but generally increasing over time.  \citet{fan09} discussed the post-emergence evolution of subsurface fields (section 8.3) and observed that the photospheric portion of the footpoints stop separating at 100 Mm which cannot be explained in the $\Omega$-loop model.

In order to better understand the complexities of Joy's law, it is helpful to keep the following in mind.

\begin{enumerate}
\item  The pitch angle of the toroidal field beneath the surface may be a cause of the tilt angle prior to the rise of flux tubes through the convection zone.

\item Coriolis forces act on flows from the expanding plasma in the apex of the flux tubes rising through the bulk of the convection zone.  Coriolis forces increase with latitude and conversely should decrease near the Equator.

\item The high scatter in the tilt angle is attributed to the interaction of a rising flux tube with convection.  After emergence, subsurface convection should no longer impart scatter in tilt over the lifetime of the region.

\item A disconnection of the flux tube from the source field would cause the tilt angle to no longer be affected by the intial pitch angle but instead relax to the angle held by the legs at the disconnection depth.

\end{enumerate}

\newpage
\section{Data

Images taken by the Helioseismic and Magnetic Imager on board NASA's Solar Dynamics Observatory were used to calculate bipolar sunspot tilt angles for HMI - Debrecen Data (HMIDD)\footnote{$http://fenyi.solarobs.unideb.hu/ESA/HMIDD.html$} from 2010 April 30 to the present.  The calculation methods employed were an extension of those used by \citet{gyo11} on the SOHO/MDI - Debrecen Data (SDD). After correcting for limb darkening and flat field effects, Sunspot Automatic Measurement (SAM) software determined the penumbra borders from the first contour having a local maximum in an averaged gradient along contour (AGAC) with the umbra border contour having the global maximum in the AGAC.  Umbral area is defined by the number of pixels within the umbral border and reported in the HMIDD as millionths of solar hemispheres (MSH).  Within the umbra border, the centroid of the pixels weighted by intensity determined the umbral latitude and longitude of the spots.

Line-of-sight magnetic field information and umbral area measurements were used to calculate the mean latitude and longitude of the leading and following sunspot groups.  The HMIDD data used the polarity of sunspots only to separate the following and leading groups and do not indicate whether Hale's polarity rule is observed or not. After grouping sunspots by polarity, only longitude determines which group is considered the leading group.\footnote{$ftp://fenyi.solarobs.unideb.hu/pub/SDO/additional/tilt\_angle/Readme.txt$}  The HMIDD calculation of the tilt angle included the area-weighted latitude, area-weighted longitude, and latitude of the centroid of the entire bipolar region \citep[Eq. 1]{bar15}.  The latitude and longitude of the leading and following spot groups were determined by averaging the positions of all of the individual spots (weighted by area) within their respective group. We calculate the separation in degrees from the latitude and longitude of the leading and following spot groups and convert to Mm by equating 1\degr\ in separation with 12.13 Mm on the solar surface.

HMIDD data report tilt angles as positive in either hemisphere if the leading spot group is closer to the Equator than the following group.  Joy's law would be observed in the tilt angle as a function of the unsigned latitude.  However, we emphasize that near-Equator measurements of the tilt are incomplete since bipolar regions are assigned a hemisphere by latitude without regard to polarity \citep{mcc14}. At the time of publication, the HMIDD data only contained the beginning of Solar Cycle 24 when sunspot activity occurs at higher latitudes.

\section{Bipolar Sunspot Behavior During Emergence and Decay

Trends in tilt angle are difficult to observe in individual active regions and become more apparent in larger samples.  We identify 1151 NOAA regions in the HMIDD data and determine that 1111 regions had at least one umbral tilt angle reported. We use the umbral calculations of the latitude and longitude for the leading and following sunspot groups to determine footpoint separation, which we report as the distance between the centroids of the opposite polarities in the photosphere.  NOAA active regions often contain new emergence activity after previous umbral activity has stopped.  We exclude new activity in a particular NOAA region if the umbral activity was not reported for more than 24 hr. To minimize the foreshortening distortion of active regions observed near the limb, we limit the data to observations taken within 0.7 solar radii from the center of the Sun's disk. The hourly cadence of HMIDD data allow for binning of the tilt angle, total umbral area, and footpoint separation over 8 hr intervals for each active region.

Individual active region information is recorded at various stages of development and decay depending on where the activity occurs on the solar surface in relation to the observation sight lines.  It would not be as useful to observe tilt angle behavior over the lifetime of an active region unless we calibrate the data to an active region characteristic that is observable in each region.  The onset of emergence would be an ideal reference point, but requiring regions to emerge on disk significantly limits the number of viable regions for study and emphasizes the emergence period over decay. The observational data of an active region that emerges on disk are more likely to exclude decay information, especially for longer-lived regions, as the regions rotate out of sight.

The umbral area bin with the maximum value establishes the end of emergence and the beginning of sunspot decay for that region, creating a suitable common reference point across all of the active regions and placing equal emphasis on emergence and decay observations.   Data which start or end with peak umbral area for a region are excluded as these regions most likely began their decay prior to appearing on disk or did not complete their emergence period before vanishing off disk due to solar rotation away from the observational line of sight.  According to this restriction and all previously stated parameters, the number of viable regions to date available for study is limited to 235.

\citet{mun15} reconcile the area and flux distributions of the photospheric magnetic structures from multiple sunspot and active region databases into a composite of Weibull and log-normal distributions for flux below $10^{21}$ Mx and above $10^{22}$ Mx, respectively. They suggest that two separate mechanisms are ``giving rise to visible structures on the photosphere: one directly connected to the global component of the dynamo (and the generation of bipolar active regions), and the other with the small-scale component of the dynamo (and the fragmentation of magnetic structures due to their interaction with turbulent convection)'' \citet[p.18]{mun15}. Mu\~noz-Jaramillo et al. show that a shift in the HMI sunspot data from a Weibull distribution ($<10^{21}$ Mx) to a log-normal distribution ($>10^{22}$ Mx) occurs around 90 MSH in the umbral area. We use 90 MSH to separate tilt angle and footpoint separation data by peak umbral area into two data sets. A substantial number of regions fall below this threshhold, so we use 45 MSH to distinguish between small and midsize regions less than 90 MSH in peak umbral area.

For each active region with a maximum umbral area ($UA_{max}$) of less than 45 MSH, we the sort tilt angle into 8 hr bins and find the median of each bin.  We do the same for the mean footpoint separation and mean umbral area. The median serves as a better measure of the center for low-sampled degree measurements that might include positive and negative values, however, the mean is preferred for non-negative measures of the separation and area. The time at which the peak umbral area is observed is noted such that all data before (after) that time are considered as emergence (decay).  The tilt, separation, and area values corresponding to the time of maximum umbral area are plotted at the t=0 point along the x-axis, see Figure \ref{hmi9x}.  We repeat the process to create two more data sets of midsize (45 $\le$ $UA_{max}$ $<$ 90 MSH) and large ($UA_{max}$ $\ge$ 90 MSH) regions.  Error bars are overplotted as the standard error of the mean. We excluded bins more than 72 hr (96 hr) from time zero for small and midsize (large) regions due to low sampling sizes at these times.

\begin{figure}  [!ht]
\centerline{\includegraphics[width=0.99\textwidth,clip=]{hmi9x}}
\caption{Umbral tilt angle, footpoint separation, and total umbral area binned in 8 hr intervals are plotted in reference to the maximum umbral area ($UA_{max}$). Regions are separated into small ($UA_{max} <$ 45 MSH), midsize (45 $\le UA_{max} <$ 90 MSH), and large regions ($UA_{max} \ge$ 90 MSH). (a) Median tilt angle of small, (b) midsize, (c) and large regions. (d) Mean footpoint separation of small, (e) midsize, (f) and large regions. (g) Mean total umbral area for small, (h) midsize, (i) and large regions. Standard error of the mean overplotted as error bars.} \label{hmi9x}
\end{figure}

At the peak in umbral area, tilt angles are least noisy and typical for early cycle (higher latitude) activity.  Midsize regions have the lowest median tilt angle (5.9\degr) compared to small (7.6\degr) and large regions (9.3\degr), see Table \ref{areabin}.  Given the concentration of all of the regions at higher latitudes (Table \ref{areabin}), it is not useful to normalize the tilt angle for latitude until an entire solar cycle is observed. The dip in small region tilt angles around -48 hr (also visible in the separation and area plots) is likely an indication of regions that are smaller than average for this group beginning their emergence with negative tilt angles.  For two days before and after the peak area, the median tilt angles of midsize regions remain fairly constant, whereas the large regions show tilt angles occurring near zero earlier in the emergence period before steadily increasing over 4-5 days well into decay.     Higher variations are observed for all of the area classifications closer to the onset of emergence and toward the end of the decay period.  It is not unexpected to see more variability in the median tilt angles of small regions, especially at the beginning of their lifetimes when the influence of external factors such as convective turbulence is greater. Larger organized regions are more resistant to convective turbulence as seen in the more stabilized tilt angles.

Tilt angle increases in larger regions during emergence and decay but remains relatively steady for the small and midsize regions. We cite this as evidence of supergranular convective flows influencing the formation and evolution of smaller regions.  It may be that the evolution of larger regions is better explained by the $\Omega$-loop model, whereas smaller regions begin to form bipoles primarily from supergranular convection, although this does not preclude smaller regions from evolving beyond the supergranule model into larger regions.  This remains a topic of interest beyond the scope of our study.

\begin{table} [!ht]
\begin{center}
 \caption{Tilt Angle, Area, and Latitude at Maximum Umbral Area ($UA_{max}$)}\label{areabin}
 \begin{tabular}{l l l l}
 \tableline\tableline
 & Small & Midsize & Large\\
 & ($UA_{max} < 45$ MSH) & (45 $\le UA_{max} <$ 90 MSH) & ($UA_{max} \ge$ 90 MSH)\\
 \tableline

Median Tilt & 7.6\degr & 5.9\degr & 9.3\degr \\
Mean Area & 18 MSH & 52 MSH & 107 MSH \\
Median Latitude & 16.7\degr & 16.1\degr & 15.2\degr \\
n (viable) & 149 & 60 & 26\\
N (original) & 618 & 261 & 232\\
 \tableline
 \end{tabular}
 \end{center}
\end{table}

All of the regions demonstrate a 3-4 day period of increase in footpoint separation with the onset of this period varying by region size.  Approximately 3 days before the umbral area peaks and lasting a day into decay, large regions increase in separation from around 45 to 75 Mm.  Midsize regions show a 30 Mm increase over 4 days as well, although the separation values are shifted downward slightly (40-70 Mm) and forward in time by about a day.  Small regions maintain 35-40 Mm in separation throughout the observed portion of the emergence period, then begin a 3 day period of increase in separation at the onset of decay, peaking near 70 Mm. At -16 hr, small regions organize around 35 Mm in separation with almost no variation.  This coincides with the beginning and end of two trends in separation for these regions: the onset of a steady increase in separation after a period of relatively constant mean separation values.

The footpoint separation for small regions during the observed portion of the emergence phase remains relatively constant before beginning a period of separation at the start of decay, whereas the larger regions begin their separation prior to decay.  We suggest that these smaller regions have not yet accumulated enough flux to overcome the influences of supergranular convection before decay begins.  The divergent flow at the top of a supergranule cell pushes the flux to the cell boundaries as described by \citet{sch68} while the magnetic structure of the bipole begins to dictate the size and shape of the supergranule cell.  Smaller regions cannot begin their separation phase until the supergranular cell that aided in its formation begins to dissipate, typically after 1-2 days \citep{hir08}. We suggest that larger regions are less influenced by supergranular convection during the emergence phase and simply separate beyond the typical size of a supergranule cell as a result.

The mean peak umbral areas for small regions ($\approx18$ MSH) and midsize regions ($\approx 52$ MSH) skew significantly lower in their respective area bins (Table \ref{areabin}).  \citet{mun15} determined that smaller sunspot regions ($<90$ MSH) display a Weibull distribution that also skews toward lower umbral areas. An empirical distribution of HMI data (Figure \ref{MUN15fig10g}) shows the transition at 90 MSH from the Weibull distribution of smaller regions to a log-normal distribution for larger regions.

Including BMRs in the discussion extends the observation of tilt angle behavior to smaller regions that may or may not include sunspots.  It should be noted that BMR areas are reported to be up to 44 times larger than sunspot areas in the same active region \citep{cha11}.  Figure \ref{ILLfig1} shows a BMR distribution of the area and tilt angle from MDI data at higher latitudes ($| \theta| \ge 10 \degr $) where the color intensity indicates the number of bipoles relative to the total number of bipoles with the same areas. Two distinct distributions are visible due to the substantial amount of tilt angle scatter for BMRs less than 300 MSH as well as the trend toward negative tilt angles for smaller regions.

\begin{figure}[!htb]
    \centering
    \begin{minipage}{0.43\textwidth}
        \centering
        \includegraphics[width=0.99\linewidth, height=0.3\textheight]{MUN15fig10g}
        \caption{Empirical distribution of HMI sunspot areas (red), with Weibull (dashed blue line) and log-normal distributions (dotted yellow line) fitted to the darker red shade. Note the transition in distributions near 90 MSH. Reproduced with permission from \citet{mun15}.}
        \label{MUN15fig10g}
    \end{minipage}
\hfill
    \begin{minipage}{.53\textwidth}
        \centering
        \includegraphics[width=1.\linewidth, height=0.33\textheight]{ILLfig1}
        \caption{Distribution of BMR area and tilt angle from MDI data (1998-2007, latitudinal zone $| \theta| \ge 10 \degr $).  Color intensity indicates the number of bipoles relative to the total number of bipoles with the same areas.  The red line corresponds to the linear least-square fit.  Reproduced with permission from \citet{ill15}.}
        \label{ILLfig1}
    \end{minipage}
\end{figure}

\FloatBarrier

\section{Very Small Bipolar Sunspot Regions

We have already seen evidence of smaller regions emerging with negative tilt angles (Figure \ref{hmi9x}a), particularly regions with shorter emergence periods and areas which are lower than average.  We limit our observations to very small peak umbral areas ($< 9$ MSH) in the HMIDD data.  Since these regions do not last as long, we are more likely to observe the entire life of the region within 0.7 solar radii from the center of the Sun's disk.  The behavior of very small active regions gives insight into how all active regions behave when they are first forming with smaller areas.  Due to the shorter lifetimes of these regions, we can apply the same binning techniques previously described but at a higher cadence of 2 hr.

The median tilt angle, mean separation, and mean area are binned every 2 hr and calibrated to the time of peak umbral area according to previously described methods.  The standard error of the mean is overplotted as error bars for the tilt, separation, and area in Figure \ref{hmi3x}.  Given the previous parameters established for viable active regions as well as the size restriction to less than 9 MSH, our sample is small (n = 12) but still worthy of inclusion in the discussion.  Median tilt angles are anti-Joy (negative) at first, increasing to expected values in the decay period for activity at a median latitude of 16.1\degr.  Footpoint separation increases from 20 Mm to about 35 Mm in a 14 hr span ($\approx 298$ m s$^{-1}$), which is more rapidly than larger regions during the previously observed 3-4 periods of steady increase ($\approx 87$ m s$^{-1}$). The umbral area averages between 3 and 5 MSH throughout the lifetimes of these regions.

\begin{figure}  [!ht]
\centerline{\includegraphics[width=0.4\textwidth,clip=]{hmi3x}}
\caption{(a) Median umbral tilt angle, (b) mean footpoint separation, and (c) total umbral area binned in 2 hr intervals are plotted in reference to the peak umbral area for regions of less than 9 MSH peak umbral area. Standard error of the mean overplotted as error bars.} \label{hmi3x}
\end{figure}

Using MDI magnetograms, \citet{tla13} noted the anti-Joy behavior (negative tilt angles) of BMRs with areas less than 300 MSH at higher latitudes that was not observed in larger regions. \citet{cha11} report an average facular-to-sunspot ratio of approximately 44 to 1.  As a rough conversion of total umbral area to BMR area, 5 MSH in the umbral area equates to 220 MSH in the BMR area, which is less than the 300 MSH threshhold set by \citet{tla13}.  This permits a comparison of tilt angle results for our respective observations of smaller regions.  It appears that anti-Joy tilt angles occur in smaller regions during the emergence period and then rotate toward positive values as the decay period begins. Coriolis forces from the divergent plasma flow in the apex of a rising flux tube deflect E-W oriented bipoles toward the positive mean tilt angles typically reported. In contrast, any existing convergent flows would deflect toward negative tilt angles. The apex of a flux tube nearing emergence in the upper convective zone with weak toroidal field strength (15-30 kG) and low flux ($10^{20}$ Mx) hosts a converging flow \citep{fan94, web13}, producing tilt away from expected Joy's law angles.

\FloatBarrier

\section{Discussion

When the peak umbral area is used to sort sunspot regions by size and to define the emergence and decay periods for each region, evidence of two size distributions emerges in the tilt angle behavior. After a period of near-zero tilt angles during emergence, a consistent increase in the large region ($UA_{max}$ $\ge 90$ MSH) tilt angle begins 48 hr before peak umbral area and lasts several days into the decay period. Regions smaller than 90 MSH show more consistency in tilt during this same time frame. It may be that larger regions eventually accumulate enough flux to form a new dynamic, with Coriolis forces acting on plasma draining from the apex of a more organized flux tube.  At some point during the accumulation of flux as larger regions emerge, Coriolis force induced tilt overtakes the influence of toroidal fields oriented in the E-W direction such that the increase in tilt angle persists into the decay period.

We attribute the negative tilt angles in small regions ($UA_{max} < 45$ MSH) at -48 hr as evidence of new activity emerging with weaker initial field strengths \citep{web13} and consistent with other observational studies of smaller active regions \citep{tla13,ill15}. Our preliminary observations of very small regions ($UA_{max} < 9$ MSH) at a higher cadence indicate the presence of negative tilt angles during emergence. As the current solar cycle progresses, further study of small regions in the HMIDD data will be of interest.

A sustained period of increase in footpoint separation lasting several days occurs for all regions, but the onset of that period varies with region size. Small regions ($UA_{max}< 45$ MSH) do not increase in separation until the end of the emergence period, whereas larger regions begin several days earlier. Separation distances during the observed portion of the emergence period for small regions are consistently near 35-40 Mm, with almost no variation in the binned values near the peak in umbral area.  We attribute this behavior to the influence of supergranular convective cells as suggested by \citet{sch68} and discussed below. These small regions then steadily increase in separation during the first three days of decay, separating to at least 70 Mm. Midsize ($45 \le UA_{max}< 90$ MSH) and large regions ($UA_{max}\ge 90$ MSH) were observed to steadily increase in separation from 40 to 70 Mm and 45 to 75 Mm, respectively, from the start of the observed emergence period and lasting into decay. It should be noted that our observations of bipolar sunspot separation distances early in the emergence period are approximately 40 Mm, which are larger than the $\approx25$ Mm distances reported by \citet{pev03} for BMRs. This can be attributed in part to magnetic bipoles forming before the umbral intensity signature necessary for measuring bipolar sunspot separation distances is reached.

Dynamo theories of the solar magnetic cycle use the pitch angle of subsurface toroidal fields as a source of the initial tilt in active regions prior to emergence \citep{bab61, lei69}, however, pitch angle is not sufficient to explain all tilt angle behavior on the surface. \citet{sch68} suggests that the Coriolis forces from a divergent supergranular flow influences the tilt angle. He states that ``the horizontal Coriolis force needs a day of unimpeded uniform motion to rotate the resulting displacement by 6 degrees. Such motion can be found in a supergranulation cell which seems to last for 1 day (cf. Simon and Leighton, 1964). The same timescale holds for the appearance of a new active center" \citet[pp. 96-97]{sch68}. In modeling magnetic activity in the convective zone, \citet{wei71} considered Coriolis-induced cyclonic motion formed by the divergent flow at the top of convection cells, also noting the concentration of magnetic flux around the perimeters of the convection cells. The supergranule divergent flow and the resultant anti-cyclonic motion push magnetic flux to the boundaries of the cell and build up tilt angles as long as the flow persists. The footpoint separation sizes of very small regions or regions early in the emergence period are reasonably comparable to the observed distribution of supergranule diameter sizes \citep{hir08}. The separation velocities of the bipolar portion of very small regions are consistent with observed supergranular divergent flows of around 300 m s$^{-1}$ \citep{rie10, lan15}. It is likely that many active regions begin their emergence period with small umbral areas and that supergranular divergent flows influence the initial tilt angle and separation distances in these just forming regions.

We cite the distinct variations in tilt angle and footpoint separation behaviors by region size as confirmation that two distributions of sunspot sizes exist. In order for a rising flux tube model to adequately describe all of the magnetic activity as it appears on the surface, some artificial assumptions are necessary. As more detailed observational data of magnetic regions become available, questions arise about the flux tube model. \citet{get15} discuss difficulties in comparing rising-tube model to observations and how a convective mechanism using in situ amplification and structuring of magnetic fields by convection avoids these difficulties. The rising flux tube model may be sufficient to describe the formation of larger bipolar sunspot regions, whereas the initial tilt angle and footpoint separation of smaller regions are heavily influenced by supergranular convection.  Further research into the cause of tilt angles in either sunspot size distribution is worth pursuing, especially the probability that a portion of the larger region distribution may be the result of smaller regions accumulating enough flux to expand beyond the influence of supergranular convection.

\FloatBarrier


\section{Application Examples}\label{sec:apps}

\section{Conclusion}
In this paper, we presented a method for realtime hand tracking using detection-guided optimization.
Our method is robust and tracks the hand at 50~fps without using a GPU.
We contribute to the tracking literature by proposing a novel representation of the input data and hand model using a mixture of Gaussians.
This representation allows us to formulate pose estimation as an optimization problem and efficiently optimize it using analytic gradient.
We also showed how additional evidence from part detection can be incorporated into our tracking framework to increase robustness.
We evaluated our method on a publicly available dataset and compared with other state-of-the-art methods.
An important direction for future work is the tracking of multiple hands interacting with each other or with objects.
We believe that the strong analytic formulation offered by our method can help solve this.

\parahead{Acknowledgments}
This research was funded by the ERC Starting Grant projects CapReal (335545) and COMPUTED (637991), and the Academy of Finland.
We would like to thank Christian Richardt.

\section{Late Fusion}\label{sec:fusion}
The goal of optimization is to find the pose $\Theta$ such that $-\mathcal{E}(\Theta)$ is minimized.
Our energies---both with and without detection---are well suited for gradient based optimization because we can derive the analytic gradient with
respect to the DOFs $\Theta$.
For efficiency, we adopt the fast gradient-based optimizer with adaptive step length proposed by~\cite{stoll_fast_2011}.

To improve robustness, especially with changes in direction and global rotation, we use multiple \emph{pose particles} for optimizing each frame.
Multiple particles improve the chances of a good initialization for optimization.
Each particle $P_i$ is initialized using the pose parameters from two previous
time steps $\Theta^{t-1}$ and $\Theta^{t-2}$ with different extrapolation factors $\alpha_{ij}$ for each DOF $j$.
This is given as
$P_i = \theta^{t-1}_j + \alpha_{ij} \, \theta^{t-2}_j, \, \forall j$.
We sample $\alpha_{ij}$ from a normal distribution with mean fixed at the initial value of $\theta_j$.
All but one of these particles is optimized using the depth-only energy.
Finally, the pose particle which converges with the best energy value is chosen
as the winning pose. In all our experiments, we found that 2--3 particles
were sufficient to obtain more robust results. Increasing the particles
had a negative effect and caused jitter in the final pose.
Each particle used 10--30 iterations per frame.
We justify the choice of these parameters in Section~\ref{sec:results}.

\section{User Specific Hand Modeling}\label{sec:handmodeling}
Accounting for the fact that there are large variations in anthropometric dimensions, our pose optimization method works best with a customized hand model for each user.
Our method does not necessitate laser scans, manual tuning of the model, nor semi-automatic bone model optimization as used by existing methods~\cite{stoll_fast_2011,sridhar_interactive_2013}.

We observed in our experiments that the primary variations in hand dimensions are finger thickness, hand length and width.
We developed a simple strategy where a default hand model is scaled using three parameters: hand length, width, and variance of Gaussians.
To find the scaling parameters for a user, we perform a greedy search over a fixed range for each scaling parameter.
At each point on this parameter grid we evaluate the energy function value from Equation~\ref{eqn:objective}.
The parameters that obtain the best energy are selected as the model scaling parameters.
This method is fast and takes less than a second to find a user-specific hand model.
Figure~\ref{fig:modelfit} shows some qualitative results from our model fitting strategy for different users.
\begin{figure}
  \centering
  \subfigure{\includegraphics[height=0.15\textwidth]{content/images/results/modelfit/user1_spheres.png}}
  \subfigure{\includegraphics[height=0.15\textwidth]{content/images/results/modelfit/user2_spheres.png}}
  \subfigure{\includegraphics[height=0.15\textwidth]{content/images/results/modelfit/user3_spheres.png}}
  \subfigure{\includegraphics[height=0.15\textwidth]{content/images/results/modelfit/user4_spheres.png}}
  \caption{Automatic fitting of user specific hand model for 4 subjects, one of whom is wearing a thick glove to simulate variability in hand dimension. The red spheres denote 3D Gaussians.}
  \label{fig:modelfit}
\end{figure}

\section{Input and Model Representation}\label{sec:input}
In the past, representations such as spheres or cylinders have been used to represent the hand model~\cite{oikonomidis_efficient_2011,qian_realtime_2014}.
Similarly, downsampled images~\cite{wang_6d_2011,wang_real-time_2009} or silhouettes~\cite{ballan_motion_2012} have been used as representations of input data.
However, such representations make pose optimization energies discontinuous and difficult to optimize.
Our novel representation of depth and 3D model data uses a mixture of \emph{weighted} Gaussian functions to
represent both depth data and the hand shape.
We were inspired by~\cite{stoll_fast_2011} who use multiple 2D RGB images and \cite{kurmankhojayev_monocular_2013} who use depth data.
Both methods rely on a uniformly weighted Gaussian mixture, and a 2D or 3D error metric for pose estimation.
However, we make important modifications that allows representing 3D depth data using a 2.5D formulation since data from depth sensors contains
information only about the \emph{camera-facing} parts of the scene.
Thus, we enable pose estimation based on alignment to a single depth image using a 2.5D error metric.

An instance of the input depth or the hand model can be represented as a mixture of Gaussian functions
\begin{align}
  \mathcal{\mathcal{C}}(\mathbf{x}) = \sum_{i=1}^n w_i \, \mathcal{G}_i(\mathbf{x}; \sigma, \boldsymbol{\mu}),
  \label{eqn:gmm}
\end{align}
where $\mathcal{G}_i(.)$ denotes a unnormalized Gaussian function with isotropic variance, $\sigma^2$, in all dimensions of $\mathbf{x} \in \mathcal{R}^n$, and mean $\boldsymbol{\mu}$.
The Gaussian mixture representation has many advantages.
First, it enables a mathematically smooth pose estimation energy which is analytically differentiable.
Second, only a few Gaussians are needed for representing the input depth and the hand model, an implicit data reduction which makes optimization extremely fast.
Finally, it provides a natural way to compute collisions using an analytically differentiable energy.
We show later that collisions are important for pose estimation (Section~\ref{sec:tracking}).
To aid visualization we henceforth represent each Gaussian in the mixture as a sphere ($\mathbf{x} \in R^3$) or circle ($\mathbf{x} \in R^2$) with a radius of $1 \, \sigma$.
However, Gaussians have infinite support ($\mathcal{C}(\mathbf{x}) > 0 $ everywhere) and can produce long range attractive or repulsive \emph{force} during pose optimization.

\subsection{Depth Data Representation}\label{sec:depth_rep}
The depth camera outputs depth maps, \emph{i.e.\ } each pixel has an associated depth value.
Depth maps contain only the camera-facing parts of the scene and information about occluded parts is unavailable.
We therefore only represent the camera-facing pixels using Gaussian mixtures, which are computed in real-time.

First, we decompose the depth image into regions of homogeneous depth using a quadtree.
The quadtree recursively decomposes depth image regions further, until the depth difference between
the furthest and nearest point in a region is below a threshold $\epsilon_c$ ($\epsilon_c = 20$~mm in all experiments).
To each quad in the tree, we fit a Gaussian function with $\boldsymbol{\mu}$ set to the center of the quad, and $\sigma = {a}/\sqrt{2}$, where $a$ is the side length of the quad.
We also set each Gaussian function to have unit weight $w_i$ since we consider all input data to be equally important.
This leads us to an analytic representation of the \emph{camera-facing surface} of the input depth, $\mathcal{C}_I(\mathbf{x}) = \sum_{q=1}^n \mathcal{G}_q(\mathbf{x})$,
where $\mathbf{x} \in \mathcal{R}^2$ and $n$ is the number of leaves in the quadtree.
Additionally, each quad has an associated depth value, $d_q$, which is the mean of all depth pixels within the quad.
Figure~\ref{fig:pipeline} illustrates the process of converting input depth to a Gaussian mixture.

\subsection{Hand Model}
We model the volumetric extent of the hand analytically using a mixture of 3D Gaussian functions, \mbox{$\mathcal{C}_h(\mathbf{x}) = \sum_{h=1}^m w_h \, \mathcal{G}_h(\mathbf{x})$}
where $\mathbf{x} \in \mathcal{R}^3$ and $m$ is the number of Gaussians.
We assume that the best fitting model has Gaussians whose isosurface at $1 \, \sigma$ coincides with the surface of the hand.
In Section~\ref{sec:handmodeling} we present a fully automatic procedure to fit such a hand model to a user.
Additionally, $\mathcal{C}_h$, is attached to a \emph{parametric}, kinematic skeleton similar to that of~\cite{simo_serra_kinematic_2011}, \emph{i.e.\ } each 3D Gaussian is attached to
a bone which determines its mean position in 3D.
We use $|\Theta| = 26$ skeletal pose parameters in twist representation, including 3 translational DOFs, 3 global rotations, and 20 joint angles.

\parahead{Model Surface Representation}
$\mathcal{C}_I$ is a representation of the \emph{camera-facing surface} while $\mathcal{C}_h$ represents the full volumetric extent of the hand.
In order to create an equivalent representation of the hand model that approximates the camera-facing parts, which we later use in pose optimization (Section~\ref{sec:tracking}).
For each model Gaussian in $\mathcal{C}_h$, we create a new projected Gaussian such that the projected hand model has the form $\mathcal{C}_p = \sum_{p=1}^m w_p \, \mathcal{G}_p(\mathbf{x})$
where $\mathbf{x} \in R^2$ and $w_p = w_h \, \forall \, h$.
$\mathcal{C}_p$ is a representation of the hand model as seen from the perspective of the depth camera and is defined over the depth image domain.
The parameters of each Gaussian $\mathcal{G}_p$ are set to be $(\boldsymbol{\mu}_p, \sigma_p)$,
where $\boldsymbol{\mu}_p = \mathbf{K} \, [ \,\mathbf{I} \, |\, \mathbf{0} \,] \, \boldsymbol{\mu}_h$.
Like~\cite{stoll_fast_2011} we approximate the perspective projection with a scaled orthographic projection, yielding 2D Gaussians
with $\sigma_p = \sigma_h \, f / \left[\boldsymbol{\mu}_p\right]_z$.
Here $f$ is the focal length of the camera, and $\left[\boldsymbol{\mu}_p\right]_z$ denotes the $z$-coordinate of the Gaussian mean.

\section{Introduction}
There is increasing interest in using markerless hand tracking in human-computer interaction, for instance when
interacting with 3D applications, augmented reality, smart watches, and for gestural input~\cite{kim_digits:_2012,lee_spacetop:_2013,wang_6d_2011}.
However, flexible, realtime markerless tracking of hands presents several unique challenges.
First, natural hand movement involves simultaneous control of several ($\geq$ 25) degrees-of-freedom (DOFs), fast motions with rapid changes in direction, and self-occlusions.
Tracking fast and complex \emph{finger articulations} combined with global motion of the hand at \emph{high framerates} is critical but remains a challenging problem.
Second, many methods use dense camera setups~\cite{oikonomidis_full_2011,sridhar_interactive_2013} or GPU acceleration~\cite{oikonomidis_efficient_2011}, \emph{i.e.\ } have high \emph{setup costs} which limits deployment.
Finally, applications of hand tracking demand tracking across many camera-to-scene configurations including desktop, egocentric and wearable settings.

This paper presents a novel method for hand tracking with a single depth camera that aims to address these challenges.
Our method is extremely fast (nearly equalling the capture rate of the camera), reliable, and supports varying close-range camera-to-hand arrangements
including desktop, and moving egocentric (camera mounted to the head).

The main novelty in our work is a new \textbf{detection-guided optimization} strategy that
combines the benefits of two common strands in hand tracking research---model-based generative tracking and discriminative hand pose detection---into a unified framework that yields high efficiency and robust performance and minimizes their mutual failures (see Figure~\ref{fig:pipeline}).
The first contribution in this strategy is a novel, efficient representation of both the input depth and the hand model shape as a mixture of Gaussian functions.
While previous work used primitive shapes like cylinders~\cite{oikonomidis_efficient_2011,oikonomidis_full_2011} or spheres~\cite{qian_realtime_2014} to represent
the hand model, we use Gaussian mixtures for both the depth data and the model.
This compact, mathematically smooth representation allows us to formulate pose estimation as a 2.5D generative optimization problem in depth.
We define a new \textbf{depth-only} energy, that optimizes for the similarity of the input depth with the hand model.
It uses additional prior and data terms to avoid finger collisions and preserve the smoothness of reconstructed motions.
Importantly, since the energy is smooth, we can obtain analytic gradients and perform rapid optimization.
While pose tracking on this energy alone could run in excess of 120~fps using gradient-based local optimization, this often results in a wrong local pose optimum.

The second contribution in our strategy is thus to incorporate evidence from trained randomized decision forests that label depth pixels into predefined parts of the hand.
Unlike previous purely detection-based approaches~\cite{Fanello_2014,shotton_real-time_2011}, we use the part labels as additional constraints in an augmented version of the aforementioned depth-only energy, henceforth termed \textbf{detection-guided} energy.
The part labels include discriminative detection evidence into generative pose estimation.
This enables the tracker to better recover from erroneous local pose optima and prevents temporal jitter common to detection-only approaches.
The precondition for recovery is reliability of the part labels.
However, even with large training sets it is hard to obtain perfect part classification (per-pixel accuracy is usually around 60\%).
Thus, pose estimation based on this additional discriminative evidence is also not sufficient.

Our third contribution therefore, is a new \textbf{late fusion approach} that combines particle-based multi-hypothesis optimization with an efficient local gradient-based optimizer.
Previous work has used particle-based optimizers, but they tend to be computationally expensive~\cite{oikonomidis_efficient_2011,oikonomidis_full_2011}.
Our approach is fast because we combine the speed of local gradient-based optimization with the robustness of particle-based approaches.
At each time step of depth video, a set of initial pose hypotheses (particles) is generated, from which a subsequent local optimization is started.
Some of these local optimizers use the depth-only pose energy, some others use the detection-guided energy.
In a final late fusion step the best pose is chosen based on the pose fitting energy.

Our approach results in a temporally stable and efficient tracker that estimates full articulated joint angles of
even rapid and complex hand motions at previously unseen frame rates in excess of 50~fps, even with a CPU implementation.
Our tracker is resilient to erroneous local convergence by resorting to the detection-guided solution when labels can be trusted, and it
is not misled by erroneous detections as it can then switch to the depth-only tracking result.

We show these improvements with (1) qualitative experiments, (2) extensive evaluation on public datasets, and (3) comparisons
with other state-of-the-art methods.

\section{Overview}

\section{Related Work}
In this brief review we focus on previous approaches to markerless hand tracking from depth images.
First, we briefly discuss marker-based and multi-camera techniques.
Gloves fitted with retro-reflective markers or color patches were used to estimate the
kinematic skeleton using inverse kinematics~\cite{sturman_survey_1994,wang_real-time_2009,zimmerman_hand_1986}.
Research on \emph{markerless} tracking was made popular in the early 2000s (\emph{e.g.\ }~\cite{athitsos_estimating_2003,wu_capturing_2001}).
Some recent solutions assume a multi-camera setup with offline
processing~\cite{ballan_motion_2012,oikonomidis_full_2011,wang_video-based_2013},
while others track at interactive rates~\cite{sridhar_interactive_2013,wang_6d_2011} of up to 30~fps~\cite{sridhar2014real}.
However, calibrated multi-camera setups make these methods difficult to adopt for practical applications.
The recent introduction of consumer depth sensors has resulted in a number of methods that require only a single depth camera.
Some commercial solutions exist, such as the Leap Motion.
Although Leap Motion is fast, the approach uses strong priors and fails with complex self-occlusions and non-standard motions (we show an example in Section~\ref{sec:results}).

The main approaches to real-time hand tracking can be divided into two classes: (1) generative and (2) discriminative methods.\footnote{There are algorithmic
parallels to full-body tracking~\cite{baak_data-driven_2011,hutchison_real-time_2012,kurmankhojayev_monocular_2013,shotton_real-time_2011}.}
First, \cite{hamer_tracking_2009} proposed a method to track a hand manipulating an object that takes 6.2~s/frame.
\cite{oikonomidis_efficient_2011} proposed a model-based method that made use of particle-swarm optimization.
This method requires GPU acceleration to achieve 15~fps and uses skin color segmentation which is sensitive to lighting.
They showed an extension to interacting hands, although only offline~\cite{Oikonomidis:CVPR2012,oikonomidis_evolutionary_2014}.
\cite{melax_dynamics_2013} proposed a tracking method directly in depth by efficient parallel physics simulations.
While this method is fast, finger articulations are often incorrectly tracked, as we demonstrate later.
Recent real-time surface tracking methods from depth~\cite{Zollhofer:2014} were applied to hands, but are limited to simple motions with no occlusions.

Second, decision forests were used with great success for full body tracking~\cite{girshick_efficient_2011,shotton_real-time_2011} and later adopted to hand tracking with varying success.
\cite{keskin_real_2011} proposed a method for recognizing finger spelling in depth data using classification forests.
\cite{Fanello_2014,tang_latent_2014,tang_real-time_2013,xu_efficient_2013} also proposed methods based on variants of random forests.
Tompson \emph{et al.\ }~\cite{Tompson_TOG_2014} track hand motion from depth at $\leq$ 25~fps using feature detections from a
convolutional network and further pose refinement through inverse kinematics.
However, a common problem with these approaches is jitter due to missing temporal information at each time step.
We provide a direct comparison with one recent method~\cite{tang_latent_2014} to demonstrate this.
Moreover, most methods estimate joint positions with temporally varying bone lengths, limiting applicability

\cite{sridhar_interactive_2013} proposed combining discriminative and generative hand pose estimation.
Their approach detected only fingertips, which could easily be occluded or misdetected.
Offline tracking in RGBD using a combination of discriminative and generative pose estimation was shown in~\cite{Gall_GCPR_2014}.
\cite{qian_realtime_2014} proposed a method based on optimization in combination with discriminative fingertip detection, achieving 25~fps.
However, tracking would be hard with this method when one or more of the fingertips are occluded.

In this paper we present a method that combines decision forests and pose estimation in a unified optimization framework.
To our knowledge, ours is the first method to track rapid articulations at 50~fps using a single depth camera and yet achieve state-of-the-art accuracy.

\section{Results and Evaluation}\label{sec:results}
We provide quantitative and qualitative evidence for performance with fast motions and finger articulations.
Evaluation of hand tracking algorithms is challenging because ground truth data is difficult to obtain.
Marker-based motion capture is often problematic due to self-occlusions.
Many methods have therefore resorted to evaluation on synthetic data~\cite{oikonomidis_efficient_2011,oikonomidis_full_2011} which, however, is not representative of real hand motions.
There are also no established benchmark datasets with accepted error metrics, and only a few implementations have been made public.

We use the dataset from \cite{sridhar_interactive_2013} which consists of seven challenging sequences (abduction--adduction, finger counting, finger waving, flexion--extension, pinching, random motions, grasping)
that are further split into slow and fast parts.
The fingertips are annotated manually in the depth data thus making it possible to compare with the multi-view approaches of \cite{sridhar_interactive_2013} and \cite{sridhar2014real}.
Additionally, we also compare with the discriminative method of \cite{tang_latent_2014} on 3 sequences.
We also motivate the need for our fusion strategy, parameter selection in optimization, and analyze the effects of different components of our objective function.
We also provide details about our framerate and qualitative evidence of improvements over \cite{melax_dynamics_2013} and the Leap Motion.
Please see the supplementary material for more results.

\parahead{Error Metrics}
Our evaluations concern the average fingertip localization error which correlates well with overall pose accuracy.
For each sequence, we compute Euclidean error of the 5 fingertip positions averaged over all frames.
Additionally, we use a second error metric~\cite{qian_realtime_2014} which is the percentage of frames that have an error of less than $x$~mm where $x \in \{15, 20, 25, 30\}$.
This is a stricter measure that highlights reliability.

\subsection{Quantitative Evaluation}
\parahead{Accuracy}
Figure~\ref{fig:dexter1} shows our average error compared with that of \cite{sridhar_interactive_2013}, \cite{sridhar2014real}, and \cite{tang_latent_2014}.
Our method produces the lowest average error of $\textbf{19.6}$~mm while using only a single depth camera.
The multi-view approaches of \cite{sridhar2014real} and \cite{sridhar_interactive_2013} have errors of $\textbf{24.1}$~mm and $\textbf{31.8}$~mm respectively.
The detection-based discriminative method of \cite{tang_latent_2014} has an error of $\textbf{42.4}$~mm (3 sequences only) highlighting the need for using temporal information.
We observe that our method does particularly well for motions that involve articulation of fingers such as \texttt{flexex1}.
Our worst performance was on the \texttt{random} sequence involving fast global hand rotation.
\begin{figure}
  \centering
  \includegraphics[width=\columnwidth]{content/results/accuracy.png}
  \caption{Average error over the 7 sequences in Dexter 1 and comparison with the \textbf{multi-view} methods of \cite{sridhar_interactive_2013} and \cite{sridhar2014real}, and the \textbf{detection-based} method of \cite{tang_latent_2014}. Our method achieves the lowest error on 5 of the 7 sequences and the best average error ($\textbf{19.6}$~mm).}
  \label{fig:dexter1}
\end{figure}

\parahead{Error Frequency}
Table 1 confirms the trend that our method performs well for finger articulations.
In 6 out of 7 sequences, our method results in tracking errors of less than 30~mm in 85\% of the frames.
A closer examination shows that these sequences contain complex finger articulations.
\begin{table}
  \centering
  \begin{tabular}{|c||c|c|c|c|}
    \hline
    \tabhead{Error $<$ (mm)} & \tabhead{adbadd} & \tabhead{fingercount} & \tabhead{fingerwave}\\
    \hline
    15 & 56.6 & 50.0 & 56.2\\
    20 & 70.6 & 66.5 & 71.2 \\
    25 & 76.2 & 77.7 & 78.3 \\
    30 & 84.9 & 85.8 & 85.0 \\
    \hline
  \end{tabular}

  \begin{tabular}{|>{\centering\arraybackslash}m{1.57cm}||c|c|c|c|}

    \hline
    \tabhead{Error $<$ (mm)} & \tabhead{flexex1} & \tabhead{pinch} & \tabhead{random} & \tabhead{tigergrasp}\\
    \hline
    15 &  53.7 & 56.7 & 19.1 & 62.9\\
    20 &  68.1 & 83.9 & 40.7 & 80.6 \\
    25 &  76.7 & 93.1 & 59.0 & 87.3 \\
    30 &  85.5 & 97.4 & 70.6 & 91.8\\
    \hline
  \end{tabular}
  \caption{Percentage of total frames in a sequence that have an error of less $x$~mm.}
\end{table}

\parahead{Robustness}
We measure robustness as the ability of a tracker to recover from tracking failures.
To demonstrate how our late fusion strategy and the different terms in the energy help achieve this, we show the frame-wise error
over the \texttt{flexex1} sequence (Figure~\ref{fig:flexex}).
Using the depth-only energy with all terms except $E_{sim}$ disabled (2 particles) results
in catastrophic tracking failure as shown by the accumulating error.
Adding the other terms, especially the collision penalty ($E_{col}$) term,  improves accuracy but results are still unsatisfactory.
The results from the late fusion approach show large gains in accuracy.
The errors also remain more uniform which results in temporally stable tracking with less jitter.
\begin{figure}
  \centering
   \subfigure{\includegraphics[height=0.15\textwidth]{content/results/samples.png}}
   \subfigure{\includegraphics[height=0.15\textwidth]{content/results/iter.png}}
   \caption{Effect of varying the number of particles and iterations during optimization. We found that increasing the number of particles resulted in diminishing returns.}
   \label{fig:iter}
\end{figure}

\parahead{Number of Particles and Iterations}
Figure~\ref{fig:iter} shows the effect of varying the number of particles and iterations during optimization.
As the number of particles increased we noticed very little increase in accuracy. In fact, the best accuracy was with 2 particles which we use throughout.
We noticed a reduction in error when using more number of iterations per particle but at the expense of runtime.
We therefore fixed the number of iterations to be 10.
\begin{figure*}
  \centering
  \includegraphics[width=\textwidth]{content/results/robustness.png}
  \vspace{-4mm}
  \caption{Plot of the error for the depth-only tracking and late fusion approach. Each approach was run with only the similarity term $E_{sim}$ and with all terms. Notice the catastrophic tracking failure with the depth-only energy. The late fusion strategy is robust and prevents error accumulation. The collision penalty term also results in large accuracy gains. Best viewed in color.}
  \label{fig:flexex}
  \vspace{-0.5cm}
\end{figure*}
\begin{figure*}[!htb]
  \centering
   \subfigure{\includegraphics[height=0.15\textwidth]{content/images/results/pose1.png}}
   \subfigure{\includegraphics[height=0.15\textwidth]{content/images/results/texas.png}}
   \subfigure{\includegraphics[height=0.15\textwidth]{content/images/results/pose2.png}}
   \subfigure{\includegraphics[height=0.15\textwidth]{content/images/results/pose3.png}}
   \subfigure{\includegraphics[height=0.15\textwidth]{content/images/results/pinch.png}}
   \subfigure{\includegraphics[height=0.15\textwidth]{content/images/results/pose6.png}}
   \subfigure{\includegraphics[height=0.15\textwidth]{content/images/results/cross.png}}
   \subfigure{\includegraphics[height=0.15\textwidth]{content/images/results/vsign.png}}
   \subfigure{\includegraphics[height=0.15\textwidth]{content/images/results/pose01.png}}
   \subfigure{\includegraphics[height=0.15\textwidth]{content/images/results/pose03.png}}
\framebox{
   \subfigure{\includegraphics[height=0.15\textwidth]{content/results/Intel2.png}}
   \subfigure{\includegraphics[height=0.15\textwidth]{content/results/Intel2_comp.png}}
}
\framebox{
   \subfigure{\includegraphics[height=0.15\textwidth]{content/results/Intel1.png}}
   \subfigure{\includegraphics[height=0.15\textwidth]{content/results/Intel1_comp.png}}
}
\framebox{
   \subfigure{\includegraphics[height=0.15\textwidth]{content/results/leap03.png}}
   \subfigure{\includegraphics[height=0.15\textwidth]{content/results/leap01.png}}
   \subfigure{\includegraphics[height=0.15\textwidth]{content/results/leap_comp03.png}}
   \subfigure{\includegraphics[height=0.15\textwidth]{content/results/leap_comp06.png}}
}
\vspace{2mm}
  \caption{Qualitative results from our tracking approach (top row and four leftmost in the second row). The highlighted boxes show comparison with \cite{melax_dynamics_2013} and the Leap Motion both of which produce a finger sliding effect. Our method tracks the pinch faithfully.}
  \label{fig:quals}
\end{figure*}

\parahead{Tracking Speed}
We tested the tracking speed of different variants of our method on a 3.6~GHz Intel Processor with 16~GB of RAM.
Our method was parallelized using OpenMP but no GPU was used.
All tests were done with the Intel Senz3D depth camera with a depth resolution of 320 $\times$ 240 and capture rate of 60~fps.
The decision forest when loaded in memory used 1~GB because the trees were stored as full trees.
This can be avoided by loading only nodes that are valid.
The depth-only energy when used with 1 particle, and 10 iterations per particle ran at 120~fps.
When 2 particles were used, the speed came down to 60~fps.
The late fusion approach, when used with 2 particles (10 iterations per particle), achieved a framerate of \textbf{50}~fps.
Image acquisition, part labeling, preprocessing, and creating the Gaussian mixture representation took 2~ms.
The optimization took between 18 and 20~ms.

\subsection{Qualitative Results}
We present several qualitative results from realtime sequences in Figure~\ref{fig:quals}.
The examples show motions with a wide range of finger articulations involving
abduction, adduction, flexion, extension, and considerable occlusions.
They also include common gestures such as the \emph{v-sign} and pointing.
In the boxes, we also show comparison with \cite{melax_dynamics_2013} and the Leap Motion on similar poses.
We observe a finger sliding effect in both these methods.
Pinching is an important gesture, but the Leap Motion produces sliding fingertips
which makes it hard to detect pinching gestures from the tracked skeleton.
Our method reproduces pinching faithfully as is evident from the skeleton overlaid on the depth image.
Occasional tracking failures occur with large global rotations but the detection-guided energy eventually reinitializes to the correct pose.
Please see the supplementary materials for more results on different camera arrangements, failure cases, and tracking with different subjects.

\section{Hand Pose Optimization}\label{sec:tracking}
In this section we describe our new formulation of pose estimation as an optimization problem using the Gaussian mixture representation of 2.5D depth data (See Figure~\ref{fig:pipeline}).
Our algorithm uses two variants of a model-to-image similarity energy, one that is
only based on depth data (Section~\ref{sec:depth_based_pose}), and another that is guided by decision forests-based part detection (Section~\ref{sec:detection_guided_pose}).
Pose estimates obtained with each energy are used by a late fusion approach to find the final pose estimate (Section~\ref{sec:fusion}).
Input to pose optimization at each time step of depth video is the 2.5D mixture of Gaussians representation of a depth image  $\mathcal{C}_I$.
The latter is computed after median filtering the depth (to remove flying pixels in time-of-flight data), and for a constrained working volume in
depth between 150~mm and 600~mm from the camera. The 3D Gaussian mixture of the hand model is denoted by $\mathcal{C}_h$
and its projected version is denoted by $\mathcal{C}_p$.

\subsection{Depth-Only Pose Optimization}\label{sec:depth_based_pose}
Our goal is to optimize for the skeleton pose parameters $\Theta$ that best explain the input data and are anatomically plausible.
We frame an energy that satisfies our goal while being mathematically smooth and differentiable. These properties make the energy ideal for fast optimization.

\subsection{Objective Function}\label{sec:objective}
Our new energy has the following general form:
\begin{align}
  \mathcal{E}(\Theta) &= E_{sim} - w_c \, E_{col} \nonumber\\&- w_l \, E_{lim} - w_s \, E_{smo},
  \label{eqn:objective}
\end{align}
where $E_{sim}$ is a measure of 2.5D similarity between $\mathcal{C}_I$ and $\mathcal{C}_p$,
$E_{col}$ is a penalty for collisions between Gaussians in $\mathcal{C}_h$,
$E_{lim}$ enforces a soft constraint on the skeleton joint limits,
$E_{smo}$ enforces smoothness in the tracked motion.
In all our experiments, we used fixed weighting factors chosen by searching for the best accuracy over the dataset: $w_c = 1.0$, $w_l = 0.2$, and $w_s = 1.0$.
Before describing each of the terms in detail we first introduce a measure of similarity between two Gaussian mixtures which is the basis for many of the terms in the objective.

\parahead{Gaussian Similarity Measure}
We define a similarity measure between any two pairs of Gaussian mixtures $\mathcal{C}_a$ and $\mathcal{C}_b$ as,
\begin{align}
  E(\mathcal{C}_a,\mathcal{C}_b)  \label{eqn:gausssim}
  & = \sum_{p \in \mathcal{C}_a} \sum_{q \in \mathcal{C}_b} D_{pq},\\
  \text{where}, D_{pq} & = \, w_p \, w_q \, \int_\Omega \mathcal{G}_p(\mathbf{x}) \, \mathcal{G}_q(\mathbf{x}) \, \mathrm{d}\mathbf{x},\label{eqn:dpq}
\end{align}
$\Omega$ denotes the domain of integration of $\mathbf{x}$.
This Gaussian similarity measure has a high value if the spatial support of the two Gaussian mixtures aligns well.
It bears resemblance to the Bhattacharyya Coefficient~\cite{bhattacharyya_measure_1946} used to measure the similarity of probability distributions while being computationally less expensive.

\parahead{Depth Similarity Term ($E_{sim}$)}
The 2.5D depth similarity term measures the quality of overlap between the projected model Gaussian mixture $\mathcal{C}_p$ and the image Gaussian mixture $\mathcal{C}_I$.
Additionally, this measure also incorporates the depth information available for each Gaussian in the mixture.
Figure~\ref{fig:esim} explains this term intuitively.
Two Gaussians that are close (in 2D pixel distance) in the depth image obtain a high value if their depth values are also close.
On the other hand, the same Gaussians obtain a low value if their depths are too far apart.
Formally, this term is defined as,
\begin{figure}
  \centering
  \includegraphics[height=8cm]{content/images/figures/esim.png}
  \caption{\textbf{Depth Similarity Term}: Consider the similarity value ($E_{sim}$) for a cylindrical shape represented by 3 Gaussians ($x \in \mathcal{R}^3$). The top figure shows a case where the value of $E_{sim}$ is high since the image overlap is high and the depth difference $\Delta_{pq}$ is low.
The bottom figure shows a case where the image overlap is moderate but $\Delta > 2 \, \sigma_h$ thus making $E_{sim} = 0$.}
  \label{fig:esim}
\end{figure}
\begin{align}
  E_{sim}(\mathcal{C}_p,\mathcal{C}_I)
  & = \frac{1}{E(\mathcal{C}_I, \mathcal{C}_I)} \, \sum_{p \in \mathcal{C}_p} \sum_{q \in \mathcal{C}_I} \Delta(p, q) \, D_{pq}
\end{align}
where $D_{pq}$ is as defined in Equation~\ref{eqn:dpq} and the \emph{depth similarity factor} is
\begin{align}
  \Delta(p, q) =
  \begin{cases}
    0, & \mbox{if } |d_p - d_q| \ge 2 \, \sigma_h\\ \label{eqn:dsf}
    1 - \frac{|d_p - d_q|}{2 \, \sigma_h}, & \mbox{if } |d_p - d_q| < 2 \, \sigma_h
  \end{cases}.
\end{align}
Here, $d_p$ and $d_q$ are the depth values associated with each Gaussian in $\mathcal{C}_p$ and $C_q$ respectively,
and $\sigma_h$ is the standard deviation of the \emph{unprojected} model Gaussian $\mathcal{G}_h$.
The surface depth value of each Gaussian in $C_p$ is computed as $d_p = \left[\boldsymbol{\mu}_h\right]_z - \sigma_h$.
The factor ${E(\mathcal{C}_I, \mathcal{C}_I)}$ is the similarity measure from equation~\ref{eqn:gausssim} of the depth image with itself and serves to normalize the similarity term.
The $\Delta$ factor has a support $[0, 1]$ thus ensuring the similarity between a projected model Gaussian and an image Gaussian is $0$ if they lie too far apart in depth.

\parahead{Collision Penalty Term ($E_{col}$)}
The fingers of a hand are capable of fast motions and often come in close proximity with one another causing aliasing of corresponding depth pixels in the input.
Including a penalty for collisions avoids fingers \emph{sticking} with one another and Gaussian interpenetration.
The 3D Gaussian mixture representation of the hand model ($\mathcal{C}_h$) offers an efficient way to penalize collisions because they implicitly act as collision proxies.
We define the penalty for collisions as,
\begin{align}
  E_{col}(\Theta)
  & = \frac{1}{E(\mathcal{C}_h, \mathcal{C}_h)} \, \sum_{p \in \mathcal{C}_h} \sum_{\substack{q \in \mathcal{C}_h\\ q > p}} D_{pq},
\end{align}
where $E(\mathcal{C}_h, \mathcal{C}_h)$ is the similarity measure from equation~\ref{eqn:gausssim} for the hand model and serves to normalize the collision term.
The collision term penalizes model Gaussians that collide with others but not if they collide with themselves.
As we show in the results, the collision term has a large impact on tracking performance.

\parahead{Joint Limit Penalty Term ($E_{lim}$)}
We add a penalty for poses that exceed predefined joint angle limits.
This forces biomechanically plausible poses to be preferred over other poses.
The joint limit penalty is given as,
\begin{align}
  E_{lim}(\Theta) = \sum_{\theta_j \in \Theta}
  \begin{cases}
    0, & \mbox{if } \theta_j^l \le \theta_j \le \theta_j^h\\
    ||\theta_j^l - \theta_j||^2, & \mbox{if } \theta_j < \theta_j^l\\
    ||\theta_j - \theta_j^h||^2, & \mbox{if } \theta_j > \theta_j^h\\
  \end{cases}
\end{align}
where $\theta_j^l$ and $\theta_j^h$ are the lower and higher limits of the parameter $\theta_j$ which is defined based on anatomical studies of the hand~\cite{simo_serra_kinematic_2011}.
The result is a tracked skeleton that looks biomechanically plausible.

\parahead{Smoothness Penalty Term ($E_{smo}$)}
During frame-by-frame pose optimization, noise is introduced which manifests as jitter in tracking.
To prevent this we penalize fast motions by adding a penalty as done by~\cite{stoll_fast_2011}. This term is given as,
\begin{align}
  E_{smo}(\Theta) = \sum_{j = 0}^{|\Theta|-1} \left( 0.5 \, \left( \Theta_j^{t-2} + \Theta_j^t\right) - \Theta_j^{t-1}\right)^2
\end{align}
where, $\Theta^t$ denotes the pose at time $t$. This term acts as a regularizer and prevents jitter in the tracked pose.

\subsection{Detection-Guided Pose Optimization}\label{sec:detection_guided_pose}
To increase chances of recovery when the estimated pose is at a wrong local pose optima, we use a second pose optimization energy
that includes evidence from hand part detection.
In particular we use pixel labels computed with a trained random forest~\cite{criminisi2013decision}.
Decision forests have been used before for 3D pose and joint position detection~\cite{keskin_real_2011,tang_latent_2014,tang_real-time_2013,xu_efficient_2013}.
We are interested in part labels and therefore follow an approach similar to \cite{shotton_real-time_2011} and \cite{keskin_real_2011}.
The evidence from the part labels is incorporated in our tracking.

We use 12 part labels for the hand (see Figure~\ref{fig:pipeline}) and found this to be an ideal trade-off between classification accuracy and sufficient evidence for detection-guided optimization.
We adopt the same depth features as \cite{shotton_real-time_2011}.
We use 50,000 labeled training images spanning the hand pose space.
As opposed to previous work~\cite{keskin_real_2011} that use  synthetic data, we use
\textbf{real} hand motions with part labels which were obtained using the depth-only version of our method and tracking motions slowly without causing tracking failure.
During training we trained 3 trees, each with a maximum depth of 22.
For each training image, we sampled 2000 random, foreground pixels, and evaluated 4000 candidate threshold-feature response pairs.

During quadtree clustering of the depth (Section~\ref{sec:depth_rep})
each quad is endowed with a part label, $l_q$ which is the label with the highest number of votes among all pixels in the quad.
We can now tightly integrate the part labels in the optimization by defining a pose fitting energy
identical to Equation~\ref{eqn:objective} with one exception: the depth similarity factor from Equation~\ref{eqn:dsf}
is replaced by the following \emph{label similarity factor}.
\begin{align}
  \Delta_l(p, q) =
  \begin{cases}
    0, & \mbox{if } l_p \neq l_q \mbox{ or } |d_p - d_q| \ge 2 \, R_i\\
    1 - \frac{|d_p - d_q|}{2 \, R_i}, & \mbox{if } l_p = l_q\nonumber
  \end{cases},
\end{align}
where $l_p$ and $l_q$ are the part labels
, $d_p$ and $d_q$ are the depth values
, $R_i$ refers to the radius of influence which is set to 200~mm in all our experiments.
Intuitively, $\Delta_l$ has a value of zero if the labels are different
and a value of one if the two Gaussians have identical labels and are perfectly aligned in 2.5D.
The labels $l_p$ are obtained from preassigned labels of each Gaussian in the hand model.