\section{Introduction}
We have experimental evidence for three generations of quarks and
leptons.  There are powerful constraints on the possibility of an
additional standard model like chiral generation of quarks and leptons
from LHC data.  These come mostly from the constraints on the
additional quarks.  At present we don't have a good understanding of
why there are three generations and so it is worth exploring the
possible experimental consequences of additional fermions. One
possibility is an additional vector-like generation, where a standard
model like generation is paired with one of opposite chirality. For
such fermions mass terms are allowed and so the masses of such
fermions could be much larger than the electroweak (EW) scale and
furthermore the masses of the quarks could be very different from
those of the leptons. For any value of the vector-like fermion masses
adding these particles is technically natural since there are no
quadratically divergent contributions to their masses. In this paper
we explore some of the experimental consequences of vector-like
leptons (assuming either the vector-like quarks are much heavier than
the leptons or don't exist).

Vector-like leptons have already been studied extensively in the
literature. They have been added to cancel anomalies in models that
extend the gauge group to include $U(1)$ lepton number and possibly
also baryon number (see, for example, \cite{FileviezPerez:2011pt}). We
considered the impact of vector-like leptons on the properties of the
Higgs boson under the assumption that almost all of their masses come
from their couplings to the Higgs boson~\cite{Ishiwata:2011hr}.  In
Ref.~\cite{Kannike:2011ng}, the contribution of additional leptons to
muon $g-2$ was studied. They considered several models, which consist
of extra leptons with different EW charges, and derived the Wilson
coefficient for $\bar{\mu} \sigma^{\mu\nu}\mu F_{\mu\nu}$. Here
$F_{\mu\nu}$ is field strength of the electromagnetic field. Bounds
from precision data were also discussed because the extra leptons
induce dimension six operators, such as $(H^\dagger D_\mu \tau^a
H)(\bar{L}\gamma^\mu \tau^a L)$, $(H^\dagger D_\mu H)
(\bar{L}\gamma^\mu L)$ and $(H^\dagger D_\mu H) (\bar{e}_R\gamma^\mu
e_R)$.  Here $H=(H^+,H^0)^T$ is the Higgs doublet, $L=(\nu_L,e_L)^T$ and
$e_R$ are the lepton doublets and singlets (flavor indices are
omitted), $\tau^a$ is a Pauli matrix and $D_\mu$ denotes the covariant
derivative. Constraints from precision EW measurements were also
discussed in
Refs.~\cite{delAguila:2008pw,Joglekar:2012vc,Kearney:2012zi}. In
Refs.~\cite{Joglekar:2012vc,Kearney:2012zi}, the impact of additional
vector-like matter on $h\rightarrow \tau\tau,~bb,~\gamma\gamma$ was
studied. The impact on $h\rightarrow \mu\mu$ and the renormalization
group equations for the Higgs quartic coupling, gauge couplings, and
production rate of the Higgs boson at the LHC, were also studied.
Ref.~\cite{Dermisek:2013gta} is similar.  In
Ref.~\cite{Harnik:2012pb}, flavor violating processes ($\tau
\rightarrow \mu \gamma$, $\tau \rightarrow e\gamma$, $\mu \rightarrow
e \gamma$, $\tau \rightarrow 3\mu$, $\tau \rightarrow 3e$, $\mu
\rightarrow 3e$), magnetic and electric dipole moments,
$\mu\rightarrow e$ conversion, and flavor violating Higgs decay were
studied from the effective theory point of view. They carefully
pointed out regions of parameter space which avoid a tuning of lepton
Yukawa couplings. See also Ref.~\cite{McKeen:2012av} which studied
Higgs to dilepton decay and CP and lepton flavor violating
process.

In this paper we focus on weakly coupled vector-like leptons in the
region of parameter space where they are heavy, {\it i.e.}, masses
significantly above one TeV. Although they change Higgs properties by
a very small amount and are not constrained by precision EW
measurement, we find that vector-like leptons with masses up to about
$100~{\rm TeV}$ can contribute at the present experimental limit to
lepton flavor violating muon decay and the electric dipole moment of
the electron without unnatural assignments for their Yukawa coupling
constants.  Even though our calculation of the amplitude for $\mu
\rightarrow e \gamma$ exists in the literature (see
Refs.~\cite{Kannike:2011ng,Dermisek:2013gta,delAguila:1982yu,delAguila:1982fs})
and its enhancement due to mass mixing of the vector-like fermions has
been noted previously, we do address a few new issues including the
tree-level decay rate for $\mu \rightarrow 3e$, the electric dipole
moment of the electron and the one loop correction to mass matrix for
the light leptons. Furthermore our focus on the multi-TeV regime is
different from most of the previous literature and this allows for
simplification of the formulas for dipole moments and the amplitudes
for flavor violating decays of charged leptons.

One goal of this paper is to explore whether the enhancement in the
amplitudes for the chirality flipping radiative processes $\mu
\rightarrow e \gamma$ and the electric dipole moment of the electron,
arising from heavy vector-like lepton mixing, results in sensitivity
to much higher mass scales for the new vector-like leptons. We find
that this is not the case if naturalness in the mass matrix for the
standard model charged leptons is imposed.\footnote{A similar
  situation occurs with lepton flavor violation induced at one loop by
  scalar leptoquarks, see Ref.~\cite{Arnold:2013cva} } We also stress
that for equal branching ratios $\mu \rightarrow 3e$ is sensitive to
higher masses for the vector-like leptons than $\mu \rightarrow e
\gamma$.

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

We add to the standard model a fourth generation of leptons, that have
the same quantum numbers under the gauge group as the standard model
leptons.  We call them $L_L^{\prime}$, $e_R^{\prime}$ and
$L_R^{\prime\prime}$, $e_L^{\prime\prime}$ where
$L_L^{\prime}=(\nu_L',e_L')^T$ and
$L_R^{\prime\prime}=(\nu_R'',e_R'')^T$. Here $L_L^{\prime}$,
$e_R^{\prime}$ is another copy of a standard model lepton generation,
while $L_R^{\prime\prime}$, $e_L^{\prime\prime}$ is the one with the
opposite chirality. Hence there is no gauge anomaly caused by these
extra leptons.
The new terms in the Lagrangian density that involve these additional
leptons are
\begin{eqnarray}
{\cal L}&=&-M_L\bar{L}'_L L''_R -M_E\bar{e}'_R e''_L
-h'_E\bar{L}'_L H e'_R \\
&& -h''_E\bar{L}''_R H e''_L- {\lambda}^i_E\bar{L}'_L H e^i_R
-{\lambda}^i_L\bar{L}^i_L H e'_R
+h.c. .\nonumber
\label{eq:L}
\end{eqnarray}
$L_L^i,e_R^i$ ($i=e,\mu,\tau$) are mostly the standard model leptons.
In addition there are the gauge invariant kinetic terms involving
these leptons. We have chosen a basis where terms such as
$\bar{L}_R''L_L^i$ and $\bar{e}_L''e_R^i$ are rotated away. In that
basis we write the Yukawa interactions of the ordinary generations as
\begin{equation}
\label{smyuk}
 {\cal L}_{Y}=-Y_e^{ij} {\bar L_L}^{i}H e_R^{j} + h.c..
\end{equation}
Introducing the four component fields, $e_1=(e'_L,e''_R)^T$,
$e_2=(e''_L,e'_R)^T$ and $\nu_{1}=(\nu'_L,\nu''_R)^T$, the
Lagrangian density becomes,
\begin{eqnarray}
&&{\cal L}=
  -M_L \bar{e}_{1L}e_{1R}  -M_L \bar{\nu}_{1L}\nu_{1R} -M_E\bar{e}_{2R}e_{2L}
\nonumber \\ &&
- h'_E (H^0\bar{e}_{1L} + H^+\bar{\nu}_{1L} ) e_{2R}
- h''_E(H^0\bar{e}_{1R} + H^+\bar{\nu}_{1R} ) e_{2L}
\nonumber \\ &&
- {\lambda}^i_E (H^0\bar{e}_{1L}  + H^+\bar{\nu}_{1L} ) e^i_{R}
  - {\lambda}^i_LH^0\bar{e}^i_{L} e_{2R}
+h.c..
\end{eqnarray}

We are primarily interested in the processes $\mu \rightarrow e
\gamma$, $\mu \rightarrow 3e$ and the electric dipole moment (EDM) of
the electron.  One reason is that the experimental bounds for the
processes are very stringent, {\it i.e.}, the branching fraction for
$\mu \rightarrow e \gamma$ and $\mu \rightarrow 3e$ should
satisfy~\cite{Adam:2011ch,Bellgardt:1987du}
\begin{eqnarray}
{\rm Br}(\mu \rightarrow e \gamma)  &<& 2.4 \times 10^{-12},
\label{Brmu2e_exp} \\
{\rm Br}(\mu \rightarrow 3e )  &<& 1.0 \times 10^{-12},
\label{Brmu23e_exp}
\end{eqnarray}
respectively, while the EDM of the electron is very suppressed
~\cite{Hudson:2011zz}
\begin{eqnarray}
|d_e| &<&1.05 \times 10^{-27}~ e~{\rm cm}.
\label{de_exp}
\end{eqnarray}
Another important reason is that the chirality flip in the radiative
processes ({\it i.e.}, $\mu \rightarrow e \gamma$ and electron EDM) is
not suppressed by the masses of the muon and electron respectively
since mass mixing of the heavy leptons in a loop diagram can cause the
chirality flip.  Regarding $\mu \rightarrow 3e$, we will see that this
process occurs at tree level, which enhances its importance compared
with the loop process $\mu \rightarrow e\gamma$. This is only
partially degraded by the phase space suppression associated with the
additional particle in the final state.

The standard model leptons are quite light and so if there are no
awkward cancellations in their mass matrices one expects hierarchical
dimensionless standard model Yukawa couplings $Y_e^{ij}$ in the range,
$10^{-2}-10^{-5}$. Naively one might expect that the same is true for
both the couplings $\lambda_E^i$ and $\lambda_L^{i}$ since the fields
$L_L^{\prime}$ and $e_R^{\prime}$ have the same quantum numbers as the
standard model fields. As will be shown later, the parts of the
amplitudes for $\mu \rightarrow e \gamma$ and electron EDM that are
not suppressed by the muon or electron mass are proportional to the
product of dimensionless couplings, $ \epsilon=\lambda_E \lambda_L
h_E^{\prime,\prime \prime}$. (Here the flavor indices of
$\lambda_{L,E}^i$ are omitted.) Hence, the above argument would give
rise to the expectation that $\epsilon$ is very small, ${\it i.e.}$,
in the range $10^{-4}-10^{-10}$. However this is not necessarily the
case. To make the argument quantitative we now show that it is
possible to find a global $U(1)$ symmetry that forces the standard
model Yukawa coupling constants $Y_e^{ij}$ to vanish but forces only
one of the couplings in the product $\epsilon$ to vanish.  Under this
symmetry the fields transform as: $L_L^i \rightarrow e^{i q_L \alpha}
L_L^i$, $e_R^i \rightarrow e^{i q_e \alpha} e_R^i$, $e_1 \rightarrow
e^{i\ q_1\alpha}e_1$ and $e_2 \rightarrow e^{i q_2\alpha}e_2$. We take
$q_e \ne q_L$ to forbid the standard model Yukawa interactions. Then
it is not possible to choose the charges $q_1$ and $q_2$ so all the
coupling constants $\lambda_E^{i}$, $\lambda_L^{i}$ and $h_E^{\prime,
  \prime \prime}$ are non-zero. But it is possible to choose values of
the charges so only one of these couplings is zero,
\begin{itemize}
\item[{\it i)}] $q_1=q_e$, $q_2=q_L$ forces $h_E^{\prime}=h_E^{\prime
    \prime}=0$ but allows non zero values for $\lambda_E^{i}$ and
  $\lambda_L^{i}$.
\item[{\it ii)}] $q_1=q_2=q_L$ forces $\lambda_E^{i}=0$ but allows non
  zero values for $\lambda_L^i$ and $h_E^{\prime, \prime \prime}$.
\item[{\it iii)}] $q_1=q_2=q_e$ forces $\lambda_L^{i}=0$ but allows non
    zero values for $\lambda_E^i$ and $h_E^{\prime, \prime \prime}$.
\end{itemize}
Given these results, it is perfectly reasonable to consider
the product of couplings $\epsilon$ in the range $10^{-2}-10^{-5}$ and
perhaps even larger if one takes a different point of view than we
have adopted.

We now transition to the mass eigenstate basis for the charged lepton
fields which we denote by a tilde. The $5\times5$ charged lepton mass
matrix $\cal M$ is defined by
\begin{eqnarray}
{\cal L}_{\rm mass} &=& -
\begin{pmatrix}
\bar{e}_1 & \bar{e}_2 & \bar{e}^i
\end{pmatrix}_L
\begin{pmatrix}
M_L   &  m'     & \mu^i_E \\
m''^* &  M_E    & 0_{1\times 3} \\
0_{3\times1}     & \mu^i_L & 0_{3\times 3}
\end{pmatrix}
\begin{pmatrix}
e_1 \\
e_2 \\
e^i
\end{pmatrix}_R
+h.c.
\nonumber \\
&\equiv& - \bar{E}_L {\cal M} E_R,
\end{eqnarray}
where $E=(e_1,e_2, e^i)^T$. Here we take $M_L$ and $M_E$ real and
define $m'=h'_E v$, $m''=h''_E v$, $\mu^i_E=\lambda^i_E v$,
$\mu^i_L=\lambda^i_L v$, where $v$ is the vacuum expectation value of
Higgs field. For our purposes we can set the standard model Yukawa
couplings $Y_e^{ij}$ to zero.  ${\cal M}$ is diagonalized  using
two unitary matrices as follows
\begin{eqnarray}
{\cal M}_{\rm diag} = V_L {\cal M} V_R^{\dagger},
\end{eqnarray}
where the mass eigenstates,
$\tilde{E}=(\tilde{e}_1,\tilde{e}_2,\tilde{e}^i)$, are given by
\begin{eqnarray}
\tilde{E}_R &=& V_R E_R ,
\label{ER}
\\
\tilde{E}_L &=& V_L E_L .
\label{EL}
\end{eqnarray}
Thus for $i=1,2,3$ the $\tilde{e}^i$ are the familiar $e,\mu,\tau$.
We derive $V_R$ and $V_L$ by treating off-diagonal elements in ${\cal
  M}$ as a perturbation.  More explicitly we assume that $v/M_{E,L}\ll
1$ and that $|M_L-M_E|\gg v$ and perform perturbation theory in the
small quantity $v/M_{E,L}$. At the leading non-trivial order of
perturbation theory,
\begin{eqnarray}
V_R &=&
\begin{pmatrix}
1      & -p_R  & q_R^i \\
p_R^*  &  1    & 0_{1 \times 3} \\
-q_R^{* i}&  0_{3 \times 1}    & 1_{3 \times 3}
\end{pmatrix},
\label{VR}
\\
V_L &=&
\begin{pmatrix}
1      & - p_L  & 0_{1 \times 3} \\
p_L^*  &  1     & q_L^{* i}\\
0_{3 \times 1}      &  -q_L^i  & 1_{3 \times 3}
\end{pmatrix},
\label{VL}
\end{eqnarray}
where
\begin{eqnarray}
&&
p_R = \frac{M_Lm'+M_Em''}{M_E^2-M_L^2},\ \ \ q^i_R = \frac{\mu^i_E}{M_L},
\\ &&
p_L = \frac{M_Lm''+M_Em'}{M_E^2-M_L^2},\ \ \ q^i_L = \frac{\mu^i_L}{M_E}.
\end{eqnarray}

In the neutrino sector, on the other hand, the heavy neutrino $\nu_1$
and linear combinations of ordinary neutrinos are mass eigenstates.

We will perform calculations of loop-level processes in $R_{\xi}$
gauge with $\xi=1$ where not only the contribution from the Higgs
boson $h$ but also the fictitous scalar $\phi^0$ and $\phi^+$ ({\it
  i.e.}, the neutral and charged Goldstone boson) must be included.
Using Eqs.~(\ref{ER})-(\ref{VL}) we find that their Yukawa couplings
in the mass eigenstate basis are (to the order we require)
\begin{eqnarray}
&&
{\cal L}_{\rm Yukawa} =
\nonumber \\ &&
 \frac{h}{\sqrt{2}}
\Bigl[
-\lambda^i_E \bar{\tilde{e}}_{1L} \tilde{e}^i_R
+(h''^*_E q_R^i-\lambda^i_Ep^*_L) \bar{\tilde{e}}_{2L} \tilde{e}^i_R
\nonumber \\
&& +(h''_Eq^{i*}_L+\lambda^{i*}_Lp_R) \bar{\tilde{e}}_{1R} \tilde{e}^i_L
-\lambda^{i*}_L \bar{\tilde{e}}_{2R} \tilde{e}^i_L
\Bigr]
\nonumber \\ &&
+  i\frac{\phi^0}{\sqrt{2}}
\Bigl[
-\lambda^i_E \bar{\tilde{e}}_{1L} \tilde{e}^i_R
-(h''^*_E q_R^i+\lambda^i_Ep^*_L) \bar{\tilde{e}}_{2L} \tilde{e}^i_R\nonumber \\
&&+(h''_Eq^{i*}_L-\lambda^{i*}_Lp_R) \bar{\tilde{e}}_{1R} \tilde{e}^i_L
+\lambda^{i*}_L \bar{\tilde{e}}_{2R} \tilde{e}^i_L
\Bigr]
\nonumber \\ &&
+\phi^+\Bigl[
-\lambda^i_E\bar{\nu}_{1L}\tilde{e}^i_R
+ h''_E q_L^{i*}\bar{\nu}_{1R}\tilde{e}^i_L
\Bigr]
+h.c. .
\end{eqnarray}
The flavor changing couplings of the $Z$ boson to the light charged
leptons is important for the calculation of $\mu \rightarrow 3e$. The
relevant terms are
\begin{eqnarray}
&&
{\cal L}_Z =
-\frac{g_Z}{2} Z_\mu
\left[
q^{i*}_Rq^j_R \bar{\tilde{e}}^i_R \gamma^\mu\tilde{e}_R^j
-q^{i}_Lq^{j*}_L \bar{\tilde{e}}^i_L \gamma^\mu\tilde{e}_L^j
\right]
\nonumber \\ &&
+g_Z Z_\mu \left[
\sin^2 \theta_W \bar{\tilde{e}}_R^i\gamma^\mu \tilde{e}_R^i+
\Bigl(-\frac{1}{2}+\sin^2 \theta_W \Bigr)
\bar{\tilde{e}}_L^i\gamma^\mu \tilde{e}_L^i
\right],
\nonumber \\
\end{eqnarray}
where $g_Z=\sqrt{g_1^2+g_2^2}$ ($g_1$ and $g_2$ are gauge coupling
constants of $U(1)_Y$ and $SU(2)_L$, respectively) and $ \theta_W$ is
the Weinberg angle.  Integrating out the heavy leptons induces flavor
mixing in the $Z$ couplings to the light standard model leptons which
causes $\mu \rightarrow 3e$ at tree level.

\section{Results}
The effective Lagrangian density that describes the impact of
integrating out the heavy leptons on the light charged lepton dipole
moments and on radiative transitions between light charged lepton
flavors is,
\begin{eqnarray}
{\cal L}_{\rm eff} =
C_{ij} \bar{\tilde{e}}^i_L \sigma^{\mu\nu} \tilde{e}^j_R F_{\mu\nu}
+C_{ji}^* \bar{\tilde{e}}^i_R \sigma^{\mu\nu} \tilde{e}^j_L F_{\mu\nu},
\end{eqnarray}
where $\sigma^{\mu\nu}=i[\gamma^\mu,\gamma^\nu]/2$.  The coefficients
$C_{ij}$ are derived by computing the amplitude for $\tilde{e}^i_R
\rightarrow \tilde{e}^j_L \gamma$ process
(Fig.~\ref{fig:ampMuEGam}). We find that,
\begin{eqnarray}
C_{ij} = \frac{e}{64\pi^2}
\left[
\sum_{B,I} c_{ij}^{BI} \frac{K_e(M_I,m_B)}{M_I}
-c_{ij}^{\phi^+1}\frac{K_\nu(M_1,M_W)}{M_1}
\right],
\nonumber \\
\end{eqnarray}
where $e$ is electric charge, $B=h,\phi^0$ and $I=1,2$ ({\it
  i.e.}, the index labeling heavy lepton type) and recall that
$M_1\simeq M_L$ and $M_2\simeq M_E$.

\begin{figure}[t]
  \begin{center}
    \includegraphics[scale=0.55]{ampMuEGam.eps}
  \end{center}
  \caption{Diagram for the $\tilde{e}^i_R \rightarrow \tilde{e}^j_L
    \gamma$. Here arrows show the chirality of fermions and chirality
    flip is indicated by a cross. }
  \label{fig:ampMuEGam}
\end{figure}

The coefficients $c_{ij}^{BI}$ are,
\begin{eqnarray}
c_{ij}^{h1}&=& \frac{1}{2}(h''^*_Eq^i_L+\lambda^i_Lp_R^*) \lambda^j_E,
\label{eq:ch1} \\
c_{ij}^{\phi^01}&=& \frac{1}{2}(h''^*_Eq^i_L-\lambda^i_Lp_R^*) \lambda^j_E,
\label{eq:cphi1}
\end{eqnarray}
for $h$-$\tilde{e}_1$ and $\phi^0$-$\tilde{e}_1$ loop diagrams, and
\begin{eqnarray}
c_{ij}^{h2}&=& \frac{1}{2}\lambda^i_L(h''^*_Eq^j_R-\lambda^j_Ep_L^*) ,
\label{eq:ch2} \\
c_{ij}^{\phi^02}&=& \frac{1}{2}\lambda^i_L(h''^*_Eq^j_R+\lambda^j_Ep_L^*) ,
\label{eq:cphi2}
\end{eqnarray}
for $h$-$\tilde{e}_2$ and $\phi^0$-$\tilde{e}_2$ loop diagrams, and
\begin{eqnarray}
c_{ij}^{\phi^+1} = h''^*_Eq^i_L \lambda^j_E,
\end{eqnarray}
for $\phi^+$-$\nu_1$ loop diagrams.  $K_e(M_I,m_B)$ and
$K_\nu(M_1,M_W)$ are loop function given by
\begin{eqnarray}
K_e(M_I,m_B)&=&\frac{1}{(1-r)^3}
\Bigl[
(1-r)(1-3r)-2r^2\log r
\Bigr],
\nonumber \\ \\
K_\nu(M_1,M_W)&=& \frac{1}{(1-r)^3}
\Bigl[
1-r^2+2r\log r
\Bigr],
\end{eqnarray}
where $r=m_B^2/M_I^2$ and $r=M_W^2/M_1^2$ respectively. $K_{e,\nu}
\simeq 1$ when $r\ll 1$.  Since we computed in $R_{\xi}$ gauge with
$\xi=1$ the appropriate mass for the fictitious scalars are
$m_{\phi^0}=M_Z$ and $m_{\phi^+}=M_W$. In this gauge the contribution
from the $Z$ and $W$ bosons are negligible. For numerical results, we
will take $m_h=125~{\rm GeV}$~\cite{Aad:2012tfa,Chatrchyan:2012ufa}.
In the region where $M_{L,E}$ are multi-TeV, the Wilson coefficients
$C_{ij}$ are given by the very simple expression,
\begin{eqnarray}
C_{ij} \simeq
\frac{e\lambda^i_L\lambda^j_Em''^*}{64\pi^2M_EM_L}.
\label{Cij}
\end{eqnarray}
Here we have neglected contributions to the Wilson coefficients
$C_{ij}$ that are suppressed by factors of the light lepton masses.

The process $\tilde{e}^i \rightarrow
\tilde{e}^j\tilde{e}^k\bar{\tilde{e}}^k$ is induced at tree level
(Fig.~\ref{fig:ampMuEEE}). The decay rate of this process is
given by (assuming $m_i\gg m_j,m_k$),
\begin{eqnarray}
\label{etoeee1}
&&
\Gamma(\tilde{e}^i \rightarrow
\tilde{e}^j\tilde{e}^j\bar{\tilde{e}}^j) =
\frac{m_i^5}{1536\pi^3}
\nonumber \\ &&
\times \left[
\Bigl|\frac{\lambda_E^i\lambda_E^j}{M_L^2}\Bigr|^2 (\kappa^2_L+2\kappa^2_R)
+\Bigl|\frac{\lambda_L^i\lambda_L^j}{M_E^2}\Bigr|^2 (2\kappa^2_L+\kappa^2_R)
\right],
\nonumber \\
\end{eqnarray}
for $k=j$, where $\kappa_L=-1/2+\sin^2 \theta_W$,
$\kappa_R=\sin^2 \theta_W$, and
\begin{eqnarray}
\label{etoeee2}
\Gamma(\tilde{e}^i \rightarrow
\tilde{e}^j\tilde{e}^k\bar{\tilde{e}}^k) =
\frac{(\kappa^2_L+\kappa^2_R)m_i^5}{1536\pi^3}
\left[
\Bigl|\frac{\lambda_E^i\lambda_E^j}{M_L^2}\Bigr|^2
+\Bigl|\frac{\lambda_L^i\lambda_L^j}{M_E^2}\Bigr|^2
\right],
\nonumber \\
\end{eqnarray}
for $j\neq k$.

The radiative processes which we are interested in are obtained from
the Wilson coefficients $C_{ij}$. The decay rate for
$\tilde{e}^i\rightarrow \tilde{e}^j \gamma$ is given by
\begin{eqnarray}
\Gamma(\tilde{e}^i\rightarrow \tilde{e}^j \gamma)
= \frac{1}{4\pi}(|C_{ij}|^2+|C_{ji}|^2) m_i^3,
\label{mu2e}
\end{eqnarray}
where $m_i$ is the mass of $\tilde{e}^i$ and we have assumed $m_i \gg
m_j$.

EDMs $d_i$ of the charged leptons, on the other hand, are coefficients
of the effective Lagrangian
\begin{eqnarray}
{\cal L}_{\rm EDM} = - \frac{i}{2} d_i
\bar{\tilde{e}}^i \sigma^{\mu\nu}\gamma_5\tilde{e}^i F_{\mu\nu},
\end{eqnarray}
and so $d_i = -2{\rm Im}(C_{ii})$.  For later discussion we also give
anomalous magnetic dipole moments (MDM) $a_i$ of the charged leptons,
which are similarly given by
\begin{eqnarray}
{\cal L}_{\rm MDM} =  \frac{e}{4m_{i}} a_i
\bar{\tilde{e}}^i \sigma^{\mu\nu}\tilde{e}^i F_{\mu\nu},
\end{eqnarray}
implying that  $ a_i =  (4m_i/e){\rm Re}(C_{ii})$.

\begin{figure}[t]
  \begin{center}
   \includegraphics[scale=0.55]{ampMuEEE.eps}
  \end{center}
  \caption{Diagram of $\tilde{e}^i \rightarrow
    \tilde{e}^j\tilde{e}^k\bar{\tilde{e}}^k$ process. When $k\neq j$,
    only left diagram is relevant, meanwhile both diagrams contribute
    when $k=j$.}
  \label{fig:ampMuEEE}
\end{figure}

To summarize,  we find that
\begin{eqnarray}
\label{mainresultsa}
\Gamma(\tilde{e}^i\rightarrow \tilde{e}^j \gamma) &=&
\frac{\alpha v^2 m_i^3|h^{\prime \prime}_E|^2
(|\lambda_L^i\lambda_E^j|^2+|\lambda_L^j\lambda_E^i|^2)}
{4096\pi^4M_E^2M_L^2},
\nonumber \\
 \\
d_i &=&
-\frac{e v}{32\pi^2}
\frac{{\rm Im}(h_E^{\prime \prime *}\lambda^i_L\lambda^i_E)}{M_E M_L},
\\  \label{mainresultse}
a_i &=&
\frac{m_i v}{16\pi^2}
\frac{{\rm Re}(h_E^{\prime \prime *}\lambda^i_L\lambda^i_E)}{M_E M_L}.
\label{mainresultsf}
\end{eqnarray}
Eqs.~(\ref{etoeee1}), (\ref{etoeee2}) and
(\ref{mainresultsa})-(\ref{mainresultsf}) are the principal analytic
results of this paper. They are valid in the limit that $v/M_{L,E}\ll
1$ and $m_{i}/v\ll 1$

Plugging in the values of the standard model parameters and focusing
on the quantities ${\rm Br}(\mu \rightarrow e \gamma$), ${\rm Br}(\mu
\rightarrow 3e$), $|d_e|$ and $|a_{\mu}|$ the above implies that,
\begin{eqnarray}
\label{target1}
{\rm Br}(\mu \rightarrow e \gamma)&\simeq&
2.3 \times 10^{-12}\left|{h_E^{\prime \prime} \over 10^{-4}}\right|^2
\left( {10^2 ~{\rm TeV^2} \over M_EM_L} \right)^2 \nonumber \\
&&\times (| \lambda_L^{\mu}\lambda_E^e|^2+| \lambda_L^e\lambda_E^{\mu}|^2),
\end{eqnarray}

\begin{eqnarray}
\label{target2}
{\rm Br}(\mu \rightarrow 3e)&\simeq&
1.7 \times 10^{-8}
\nonumber \\ &&  \! \! \! \! \! \! \! \! \! \! \! \! \! \! \! \!
\! \! \! \! \! \! \! \! \! \! \! \! \! \! \! \! \! \! \! \!
\times
\left[
\Bigl| \lambda^e_E\lambda^\mu_E{10^2 ~{\rm TeV^2} \over M_L^2} \Bigr|^2
+ 1.1\times
\Bigl| \lambda^e_L\lambda^\mu_L{10^2 ~{\rm TeV^2} \over M_E^2} \Bigr|^2
 \right],
\nonumber \\
\end{eqnarray}
\begin{equation}
\label{target3}
|d_e| \simeq 1.1\times 10^{-26} e~{\rm cm}
\left|{\rm Im}\left(
{h_E^{\prime \prime *}\lambda^e_L\lambda^e_E \over 10^{-4}}
\right)\right|\left( {10^2 ~{\rm TeV^2} \over M_EM_L} \right),
\end{equation}
and
\begin{equation}
\label{target4}
|a_{\mu}|\simeq  1.2\times 10^{-13}
\left|{\rm Re}\left(
{h_E^{\prime \prime *}\lambda^{\mu}_L\lambda^{\mu}_E \over 10^{-4}}
\right) \right|\left( {10^2 ~{\rm TeV^2} \over M_EM_L} \right).
\end{equation}
Eqs. (\ref{target1}), (\ref{target2}), (\ref{target3}) and
(\ref{target4}) are compared with the measured values given in
Eqs.~(\ref{Brmu2e_exp}), (\ref{Brmu23e_exp}) and (\ref{de_exp}) and the
deviation of the observed muon $g-2$ ~\cite{Bennett:2006fi} from its
standard model prediction~\cite{Davier:2010nc},
\begin{eqnarray}
\Delta a_\mu^{\rm exp} = a_\mu^{\rm exp}-a_\mu^{\rm SM} &\simeq&
(2.8\pm0.8)\times 10^{-9}.
\end{eqnarray}
Note that while it is expected that a 10 TeV vector-like lepton can
give an observable contribution to $\mu \rightarrow e \gamma$, $\mu
\rightarrow 3e$ and the electron EDM, such a heavy vector-like lepton
doesn't give a large enough contribution to the anomalous magnetic
moment to explain the discrepancy between standard model value and its
observed value.

Before going on to the numerical results, we comment on the
perturbative corrections to light lepton mass matrix that come from
integrating out the heavy leptons.  There is a tree-level correction
to light lepton mass matrix elements, which is induced by mixing with
heavy leptons. It is given by
\begin{eqnarray}
|\Delta m_{ij}^{\rm tree}| = |m'' q^i_Rq^j_L|.
\label{mijTree}
\end{eqnarray}
In addition to this correction, there is also a logarithmically
divergent one loop correction to the light lepton mass matrix that
arises from one loop self-energy diagrams analogous to the Feynman
diagrams that contribute to the radiative transitions and dipole
moments but without the photon. Using a momentum cutoff $\Lambda$, it
is
\begin{eqnarray}
&&\Delta m_{ij}^{\rm loop}
= -\frac{1}{16\pi^2}\sum_{B,I}c_{ij}^{BI}
M_I {\rm log}\left({\Lambda^2 \over M_I^2} \right) \nonumber \\
&&\simeq -\frac{\lambda_L^i\lambda_E^j m^{\prime \prime *}}{16\pi^2}
\left[{2M_L \over M_E}{\rm log}\left({\Lambda^2 \over  M_L^2}\right)
+{M_E \over M_L}
   {\rm log}\left({\Lambda^2 \over M_E ^2}\right)\right],
\nonumber \\
\label{mijLoop}
\end{eqnarray}
assuming $\Lambda \gg M_I\gg m_h,m_Z$. This quantum correction is more
important than the tree-level one when $M_{L,E}$ is greater than about
a few hundred GeV. In the region of heavy vector-like leptons masses
we are focusing on the tree level contribution from mixing with the
heavy leptons can be neglected and so in order to avoid awkward
cancellations between $\Delta m_{ij}=\Delta m_{ij}^{\rm tree}+ \Delta
m_{ij}^{\rm loop} \simeq \Delta m_{ij}^{\rm loop}$ and the standard
model tree level light lepton mass matrix elements from
Eq.~(\ref{smyuk}) (which are proportional to $Y_e^{ij}$), we adopt the
naturalness criteria~\cite{Cheng:1987rs}:
\begin{eqnarray}
|\Delta m_{ij}\Delta m_{ji}|\lesssim m_i m_j.
\label{finetune}
\end{eqnarray}
A more precise analysis of naturalness constraints would involve
solving renormalization group equations for the evolution of the
coupling constants in the theory

It is important to remember that the region of parameter space
excluded by Eq.~(\ref{finetune}) is not ruled out and can be
consistent with experiment. We view the region of parameter space that
does not satisfy Eq.~(\ref{finetune}) as disfavored by the naturalness
criteria that the one loop correction to the light lepton mass matrix
calculated with a large momentum space cutoff not overwhelm the
contribution which is coming from Eq.~(\ref{smyuk}).

We have computed the parts of the amplitudes for $\mu \rightarrow e
\gamma$ and the electric dipole moment of the electron that are not
suppressed by the light lepton masses.  The radiative correction to
the light lepton mass matrix in Eq.~(\ref{mijLoop}) is proportional to
the same combination of couplings constants that occur in the
amplitudes for these radiative processes.  Hence the naturalness
condition in Eq.~(\ref{finetune}) restricts the size of this part of
their amplitudes. There are situations where the radiative correction
to the charged lepton mass matrix is negligible and the rate for $\mu
\rightarrow 3e$ not suppressed, for example, if $h_E^{\prime \prime}$
is extremely small but $\lambda_{E,L}$ order unity. In this case the
contributions to the amplitudes for the radiative processes
proportional to the light lepton masses, that we have not calculated,
dominate. However, in that situation the rate for $\mu \rightarrow 3e$
is sensitive to higher mass scales for the vector-like leptons than
the rate for $\mu \rightarrow e \gamma$ since ${\rm Br}(\mu
\rightarrow e \gamma) \sim (\alpha/4\pi){\rm Br}(\mu \rightarrow 3
e)$.

\begin{widetext}

\begin{figure}[t]
  \begin{center}
    \includegraphics[scale=0.6]{MuEGam1.eps} \! \! \! \! \! \!
    \includegraphics[scale=0.6]{MuEEE1.eps}
    \includegraphics[scale=0.6]{de1.eps} \! \! \! \! \! \!
    \includegraphics[scale=0.6]{amu1.eps}
  \end{center}
  \caption{Contours of fixed ${\rm Br}(\mu \rightarrow e \gamma)$ (top
    left), ${\rm Br}(\mu \rightarrow 3e)$ (top right), electron EDM
    (bottom left) and muon MDM (bottom right) in the $M_L$-$|\lambda|$
    plane for the parameters discussed in the text.  Lines are ${\rm
      Br}(\mu \rightarrow e \gamma)={\rm Br}(\mu \rightarrow
    3e)=10^{-10},10^{-12},10^{-14},10^{-16}$,
    $|d_e|=10^{-26},10^{-28},10^{-30}~e~{\rm cm}$, and
    $|a_\mu|=10^{-12},10^{-14},10^{-16}$ from top to
    bottom. The dark-shaded region is excluded by experiments, while
    light-shaded (blue) indicates fine-tuned parameter space.  In the
    top left, the prospective MEG sensitivity is also plotted. }
  \label{fig:case1}
\end{figure}
\end{widetext}

\section{Numerical Results}

First we consider case {\it i)} discussed in Sec.~\ref{sec:model}. We
take $h'_E=h''_E=h=10^{-3}$, $\lambda_E^i=\lambda_L^i=\lambda$ and
$M_E=1.5M_L$.  Fig.~\ref{fig:case1} shows contours of fixed ${\rm
  Br}(\mu \rightarrow e \gamma)$, ${\rm Br}(\mu \rightarrow 3e)$,
electron EDM and muon MDM in the $M_L$-$|\lambda|$ plane. The
dark-shaded region is excluded by experiments.  We also shade the
fine-tuned parameter space where Eq.~(\ref{finetune}) is not
satisfied, taking $\Lambda = 10^{15}~{\rm GeV}$. Note that, $\mu
\rightarrow 3e$ gives a more stringent constraint on the parameter
space. This is because ${\rm Br}(\mu \rightarrow 3e)$ is not
suppressed by $h$.  The current experimental bound constrains
$|h\lambda^2|\lesssim 10^{-5}$ when $M_{L}\sim 10~{\rm TeV}$. On the
other hand, the region of parameter space where $|h\lambda^2|\gtrsim
10^{-5}$ is fine-tuned.  For the electron EDM, we take $\delta=\pi/4$
where ${\rm arg}(\lambda^2)=\delta$.  The region of the plane excluded
by experiment is almost the same as in the $\mu \rightarrow 3e$
case. Note that the rate for $\mu \rightarrow 3e$ does not depend on
$h$ and if we take $h$ much smaller then the fine-tuned region
disappears.

There are also experimental constraints on flavor changing $\tau$
lepton decays. Since the weak decay width of a charged lepton of type
$i$ is proportional to $m_i^5$, the branching fractions for $\tau
\rightarrow \mu \gamma$ and $\tau \rightarrow e \gamma$ satisfy ${\rm
  Br}(\tau \rightarrow \mu \gamma)\simeq {\rm Br}(\tau \rightarrow e
\gamma)\simeq (m_{\mu}/m_{\tau})^2{\rm Br}(\mu \rightarrow e \gamma)$
in  case {\it i)}. Current experimental bounds for these process
are by ${\rm Br}(\tau \rightarrow \mu \gamma) < 4.4\times 10^{-8}$ and
${\rm Br}(\tau \rightarrow e \gamma) < 3.3\times 10^{-8}$
~\cite{Aubert:2009ag}. Thus these bounds are less stringent than for
muon decays.

For reference, we have given the contours of fixed muon MDM. Clearly
the impact of vector-like leptons on the muon MDM is negligible in the
parameter region we have focused on. In order to explain the muon $g-2$
anomaly by vector-like leptons, they should be lighter than 1 TeV, for
example, $M_{L,E}\sim O(100~{\rm GeV})$ with $|h\lambda^2 \cos
\delta|\sim 10^{-3}$.  Note that having vector-like leptons with
masses around the weak scale doesn't necessarily lead to a conflict
with the limits on the rates for $\mu \rightarrow e \gamma$, $\mu
\rightarrow 3e$ and the electron electric dipole moment.  For example,
it could be the case that $|\lambda_{E,L}^e|$ is much smaller than
$|\lambda_{E,L}^\mu|$.

\begin{widetext}

\begin{figure}[t]
  \begin{center}
   \includegraphics[scale=0.6]{MuEGam2.eps} \! \! \! \! \! \!
   \includegraphics[scale=0.6]{MuEEE2.eps}
   \includegraphics[scale=0.6]{de2.eps}\! \! \! \! \! \!
    \includegraphics[scale=0.6]{amu2.eps}
  \end{center}
  \caption{ Same as Fig.~\ref{fig:case1} but with different choice
    of parameters.}
  \label{fig:case2}
\end{figure}

\end{widetext}

Next we consider the case {\it ii)}. (Parametrically case {\it iii)}
is identical to case {\it ii)}.) We take $h'_E=h''_E=h=1$,
$\lambda^e_L=\lambda^\mu_L=\lambda$, $\lambda_E^i=m_i/v$ and
$M_E=1.5M_L$.  Fig.~\ref{fig:case2} shows the results. (For the
electron EDM, we take $\delta={\rm arg}(\lambda)=\pi/4$.) In this case
there is no fine-tuned region (taking $\Lambda=10^{15}~{\rm
  GeV}$). Since the electron EDM is suppressed by $|\lambda_E^e|$, the
experimental bound becomes weaker compared to the previous case.  The
current bounds from the limits on the $\mu \rightarrow e \gamma$ and
$\mu \rightarrow 3e$ branching ratios, on the other hand, give a more
severe constraint on the parameter space. For $\mu \rightarrow e
\gamma$, this is because the branching ratio is only suppressed by
$|\lambda_E^\mu|$.  The numerical result indicates that the
$M_{E,L}\sim 10-100~{\rm TeV}$ region may be probed in the future MEG
experiment, which has prospective sensitivity of ${\rm Br}(\mu
\rightarrow e \gamma)\simeq 5.0\times 10^{-13}$. (For $\tau$ radiative
decay, we note that ${\rm Br}(\tau \rightarrow \mu \gamma)={\rm
  Br}(\mu \rightarrow e \gamma)$ in this case.)  The process $\mu
\rightarrow 3e$, on the other hand, is not suppressed by any of the
$|\lambda^i_E|$ (see Eq.~(\ref{target2})), which gives rise to large
branching fraction. It is especially interesting that the $M_{L,E}\sim
100~{\rm TeV}$ region with $|\lambda|\sim O(1)$ is already excluded by
current experiments, which means vector-like leptons with masses up to
100 TeV can have an observable impact on experiments.

\section{Concluding Remarks}
We have considered the phenomenological impact of vector-like leptons
that are much heavier than the weak scale. We discussed expectations
for their couplings based on the observed masses of the charged
leptons. With couplings that are reasonable from this perspective we
found that vector-like leptons with masses up to $100~ {\rm TeV}$ can
have an impact on the electric dipole moment of the electron and the
rates for $\mu \rightarrow 3e$ and $\mu \rightarrow e \gamma$ that is
detectable in the next generation of experiments. When the vector-like
leptons are this heavy their impact on the anomalous magnetic moment
of the muon and Higgs properties is very small.

For the chirality flipping radiative processes $\mu \rightarrow e
\gamma$ and the electric dipole moment of the electron we have
focused on the parts of their amplitudes that are not suppressed by
the light lepton masses. One of the goals of this work has been to
explore whether the absence of light lepton mass suppression results
in sensitivity to much higher masses for the vector-like leptons. We
found that one-loop radiative corrections to the light lepton mass
matrix are proportional to the same combination of couplings and
demanding that the one loop contribution not overwhelm the
contribution from Yukawa couplings restricts the reach in vector-like
lepton mass. Finally we noted that, for similar branching ratios, the
three body flavor changing decay $\mu \rightarrow 3e$ has sensitivity
to greater vector-like lepton masses than the radiative two body decay
$\mu \rightarrow e \gamma$.

\subsection*{Acknowledgment}
This paper is funded by the Gordon and Betty Moore Foundation through
Grant $\# 776$ to the Caltech Moore Center for Theoretical Cosmology
and Physics. The work of the authors was supported also in part by the
U.S. Department of Energy under contract No. DE-FG02-92ER40701. The
research of MBW was supported in part by Perimeter Institute for
Theoretical Physics. Research at Perimeter Institute is supported by
the Government of Canada and by the Province of Ontario through the
Ministry of Economic Development \& Innovation.


\section{Introduction}
This paper is a contribution to the theory of integral variational functionals,
depending on submanifolds of a given manifold $X$. The theory is based on geometric notions,
such as the bundles of (skew-symmetric) multivectors, and Grassmann fibrations.
Conceptually, it extends local parametric integrals of Finsler-Kawaguchi and areal geometries
(see e.g. Chern, Chen, Lam~\cite{ChernChenLam}, Davies~\cite{Davies}, Kawaguchi~\cite{AK5},
and Tamassy~\cite{Tamassy4})
to global functionals, depending on (global) submanifolds.
In Section 2 we summarize integration theory of differential forms along submanifolds.
Section 3 is devoted to vector bundles of $k$-vectors; we show how mappings of Euclidean spaces into manifolds
({\it parametrisations}) can be lifted to the bundles of $k$-vectors.
In Section 4 we introduce, using the Pl\"ucker embedding,
underlying spaces for parameter-invariant variational problems, the Grassmann fibrations.
In Section 5 we show that any $k$-form on the Grassmann fibration defines an integral variational functional,
depending on $k$-dimensional submanifolds. An example is the {\it Hilbert form}, a well-known first-order construction
in Finsler geometry and its generalisations (Chern, Chen, Lam \cite{ChernChenLam},
Crampin, Saunders~\cite{CrampinSaunders}).

It should be pointed out that the theory can be further generalised.
To this purpose one should consider higher-order Grassmann fibrations,
endowed with Lagrangians, satisfying relevant homogeneity conditions
(Zermelo conditions, see e.g. Saunders~\cite{Saunders3}, and Urban and Krupka~\cite{Krupka-Urban}).

\section{Integration over submanifolds}
Let $X$ be an $n$-dimensional manifold, $S$ a  subset of $X$, ${x_0} \in S$ a point.
A chart $(U,\varphi )$,$\varphi  = ({x^i})$, at ${x_0}$ is a {\it submanifold chart} for $S$,
if there exists a non-negative integer $k \leqslant n$ such that
$\varphi (U \cap S) = \{ x \in U|{x^{k + 1}}(x) = {c_1},{x^{k + 2}}(x) = {c_2}, \ldots ,{x^n}(x) = {c_{n - k}}\} $.
If such a chart exists, we say that $S$ is a {\it submanifold} of $X$ at the point ${x_0}$;
$k$ is the {\it dimension} of $S$ at ${x_0}$. If such a submanifold chart exists at every point of $X$,
we say $S$ is a {\it submanifold} of $X$ and call $k$ the {\it dimension} of $S$.

Denote by $({t^1},{t^2}, \ldots ,{t^n})$ the canonical coordinates on the Euclidean space
${R^n}$, and $R_{( - )}^n = \{ {t_0} \in {R^n}|{t^n}({t_0}) \leqslant 0\} $, $\partial R_{( - )}^n
= \{ {t_0} \in R_{( - )}^n|{t^n}({t_0}) = 0\} $. $R_{( - )}^n$ is the {\it halfspace} of ${R^n}$,
$\partial R_{( - )}^n$ is the {\it boundary} of $R_{( - )}^n$. Let $\Omega $ be a non-void subset of $X$,
and ${x_0} \in \Omega $ a point. A chart $(U,\varphi )$ at ${x_0}$ is said to be {\it adapted} to $\Omega $,
if the set $\varphi (U \cap \Omega )$ is an open set in $R_{( - )}^n$. $\Omega $ is a {\it piece} of $X$,
if it is compact and each point $x \in \Omega $ admits a chart, adapted to $\Omega $.

Let $\eta $ be a $k$-form on $X$. Our aim now will be to introduce an integral of $\eta $ on a piece of a
$k$-dimensional submanifold $S$ ({\it $k$-piece} of a $X$). Express $\eta $ in a submanifold chart
$(U,\varphi )$, $\varphi  = ({x^i})$, as
$\eta  = {\eta _{{i_1}{i_2} \ldots {i_k}}}d{x^{{i_1}}}
\wedge d{x^{{i_2}}} \wedge  \ldots  \wedge d{x^{{i_k}}}.$
Then restricting $\eta $ to $S$
we get from the equations $x^{k + 1} = 0$, ${x^{k + 2}} = 0$, $\cdots$ , ${x^n} = 0$
\begin{align}
\eta  = fd{x^1} \wedge d{x^2} \wedge  \ldots  \wedge d{x^k},
\end{align}
where we write $f = f({x^{{i_1}}},{x^{{i_2}}}, \ldots ,{x^{{i_k}}})$
for the component of $\eta $ restricted to $S$.
From now on we suppose that $S$ is {\it orientable}, and is endowed with an orientation
${\operatorname{Or} _S}X$; only submanifold charts on $X$ belonging to ${\operatorname{Or} _S}X$ are used.
The integral of $\eta $ on a compact set $\Omega  \subset S$ is defined in a standard way.
There exist a finite family $\{ ({U_1},{\varphi _1}),({U_2},{\varphi _2}), \ldots ,({U_N},{\varphi _N})\} $
of submanifold charts on $X$, such that the family $\{ {U_1} \cap S,{U_2} \cap S, \ldots ,{U_N} \cap S\} $
covers $\Omega $. Let $\{ {\chi _1},{\chi _2}, \ldots ,{\chi _N}\} $ be a partition of unity,
subordinate to this covering. Then
\begin{align}
\int_\Omega  \eta   = \sum\limits_{j = 1}^N {\int_{\operatorname{supp} {\chi _j} \cap \Omega } {{\chi _j}\eta } } .
\end{align}
The following basic properties of the integral are needed in the calculus of variations.
\begin{lemma} {\bf (Transformation of integration domain)}
Let $X$ and $Y$ be two smooth n-dimensional oriented manifolds,
$\alpha :X \to Y$ an orientation-preserving diffeomorphism.
Then for any compact set $\Omega  \subset S$ and any continuous differential $n$-form on $Y$
\begin{align}
\int_\Omega  \eta   = \int_{{\alpha ^{ - 1}}(\Omega )} {\alpha *\eta } .
\end{align}
\end{lemma}
\begin{lemma} {\bf (Leibniz rule)}
Let $X$ be an oriented n-dimensional manifold, ${\eta _t}$ a family of $n$-forms on $X$,
differentiable on a real parameter $t$, $\Omega  \subset S$ a compact set.
Then the function
\begin{align}
I \ni t \to \int_\Omega  {{\eta _t}}  \in R
\end{align}
is differentiable, and
\begin{align}
\frac{d}{{dt}}\int_\Omega  {{\eta _t}}  = \int_\Omega  {\frac{{d{\eta _t}}}{{dt}}} .
\end{align}
\end{lemma}
\begin{lemma} {\bf (The Stokes' formula) }
Let $X$ be an $n$-dimensional manifold, $S$ a $k$-dimensional oriented submanifold of $X$,
$\eta $ a $(k - 1)$-form on $X$.
Let $\Omega $ be a piece of $S$ with boundary $\partial \Omega $, endowed with induced orientation. Then
\begin{align}
\int_{\partial \Omega } \eta   = \int_\Omega  {d\eta } .
\end{align}
\end{lemma}

\section{Bundles of $k$-vectors}
Let $X$ be an $n$-dimensional manifold,
${\Lambda ^k}{T_x}X$ the $k$-exterior product of the tangent space
${T_x}X$, $x \in X$ a point. We put
\begin{align}
{\Lambda ^k}TX = \bigcup\limits_{x \in X} {{\Lambda ^k}{T_x}X} .
\end{align}
This set has a natural vector bundle structure over $X$, with type fibre ${\Lambda ^k}{R^n}$.
We denote by ${\tau ^k}$ the vector bundle projection of ${\Lambda ^k}TX$.

Let $X$ (resp. $Y$) be a smooth manifold of dimension $n$ (resp. $m$),
and let $f:X \to Y$ be a differentiable mapping.
Choose a point $x \in X$ and a $k$-vector $\Xi  \in {\Lambda ^k}{T_x}X$.
Then choose a chart $(U,\varphi )$, $\varphi  = ({x^i})$, at $x$ and a chart $(V,\psi )$,
$\psi  = ({y^\sigma })$, at $f(x) \in Y$ such that $f(U) \subset V$. Expressing $\Xi $
in components and setting
\begin{align}
&{\Lambda ^k}{T_x}f \cdot \Xi  = \frac{1}{{{{(k!)}^2}}}{\left( {\frac{{\partial {y^{{\sigma _1}}}
f{\varphi ^{ - 1}}}}{{\partial {x^{{i_1}}}}}} \right)_{\varphi (x)}}{\left( {\frac{{\partial {y^{{\sigma _2}}}
f{\varphi ^{ - 1}}}}{{\partial {x^{{i_2}}}}}} \right)_{\varphi (x)}} \ldots {\left(
{\frac{{\partial {y^{{\sigma _k}}}f{\varphi ^{ - 1}}}}{{\partial {x^{{i_k}}}}}} \right)_{\varphi (x)}} \hfill \nonumber \\
&\quad  \cdot {\Xi ^{{i_1}{i_2} \ldots {i_k}}}{\left( {\frac{\partial }{{\partial {y^{{\sigma _1}}}}}}
\right)_{f(x)}} \wedge {\left( {\frac{\partial }{{\partial {y^{{\sigma _2}}}}}} \right)_{f(x)}} \wedge
 \ldots  \wedge {\left( {\frac{\partial }{{\partial {y^{{\sigma _k}}}}}} \right)_{f(x)}},
\end{align}
we get a $k$-vector ${\Lambda ^k}{T_x}f \cdot \Xi  \in {\Lambda ^k}{T_{f(x)}}Y$,
and a vector bundle homomorphism ${\Lambda ^k}Tf:{\Lambda ^k}TX \to {\Lambda ^k}TY$ over  $f$ (the {\it lift} of $f$).

It is easily seen that differentiable mappings of a Euclidean space into a manifold can be canonically lifted
to the bundles of $k$-vectors. To this purpose we use the {\it canonical $k$-vector field} on ${R^n}$
\begin{align}
{R^n} \ni t \to \theta (t) = \frac{1}{{k!}}{\varepsilon ^{{i_1}{i_2} \ldots {i_k}}}
{\left( {\frac{\partial }{{\partial {t^{{i_1}}}}}} \right)_t} \wedge
{\left( {\frac{\partial }{{\partial {t^{{i_2}}}}}} \right)_t} \wedge
\ldots  \wedge {\left( {\frac{\partial }{{\partial {t^{{i_k}}}}}} \right)_t} \in {\Lambda ^k}T{R^n}.
\end{align}
Identifying ${\Lambda ^k}T{R^n}$ with ${R^n} \times {\Lambda ^k}{R^n}$,
the canonical section becomes the mapping $t \to (t,{\varepsilon ^{{i_1}{i_2} \ldots {i_k}}})$.

Consider a differentiable mapping $f:U \to Y$, where U is an open subset of ${R^n}$.
For any point $t \in U$, ${\Lambda ^k}{T_t}f \cdot \theta (t)$ is an element of the vector space
${\Lambda ^k}{T_{f(t)}}Y$. We get the {\it canonical lift} ${\Lambda ^k}f$ of $f$ to ${\Lambda ^k}TY$, defined by
\begin{align}
{\Lambda ^k}f = {\Lambda ^k}Tf \cdot \theta .
\end{align}
The canonical lift of the {\it parametrisation} $U \ni t \to ({\psi ^{ - 1}} \circ {\iota _{k,m}})(t) \in V \cap S$
is expressed in a chart $(V,\psi )$, $\psi  = ({y^\sigma })$, as
\begin{align}
&{\Lambda ^k}({\psi ^{ - 1}} \circ {\iota _{k,m}})(t) \hfill \nonumber \\
&\quad  = {\left( {\frac{\partial }{{\partial {y^1}}}} \right)_{{\psi ^{ - 1}}
\circ {\iota _{k,m}}(t)}} \wedge {\left( {\frac{\partial }{{\partial {y^2}}}} \right)_{{\psi ^{ - 1}}
\circ {\iota _{k,m}}(t)}} \wedge  \ldots  \wedge {\left( {\frac{\partial }{{\partial {y^k}}}}
\right)_{{\psi ^{ - 1}} \circ {\iota _{k,m}}(t)}}. \label{form5}
\end{align}
Formula (\ref{form5}) also defines the mapping
$V \ni y \to ({\Lambda ^k}\psi )(y) \in {({\tau ^k})^{ - 1}}(V)$ by
\begin{align}
{\Lambda ^k}\psi  = {\Lambda ^k}({\psi ^{ - 1}} \circ {\iota _{k,m}}) \circ {\operatorname{pr} _{m,k}}\psi ,
\end{align}
the {\it canonical section along} $S$, associated with $(V,\psi )$. ${\Lambda ^k}\psi $ is expressed by
\begin{align}
&({y^1},{y^2}, \ldots ,{y^k},{y^{k + 1}},{y^{k + 2}}, \ldots ,{y^m}) \to {\Lambda ^k}({\psi ^{ - 1}}
\circ {\iota _{k,m}})({y^1},{y^2}, \ldots ,{y^k}) \hfill \nonumber \\
&\quad  = (({y^1},{y^2}, \ldots ,{y^k},0,0, \ldots ,0),(1,0,0, \ldots ,0)).
\end{align}
Writing in the multi-index notation $({({\tau ^r})^{ - 1}}(V),\Phi )$,
$\Phi  = ({\dot y^I})$, and setting ${I_0} = (12 \ldots k)$,
we get the image of this mapping as a subset of
${({\tau ^r})^{ - 1}}(V)$, defined by the equations
${y^{k + 1}} = 0$, ${y^{k + 2}} = 0$, $ \ldots $, ${y^m} = 0$,
${\dot y^{{I_0}}} = 1$, ${\dot y^I} = 0$, $I \ne {I_0}$.
\begin{lemma}
Let $(V,\psi )$, $\psi  = ({y^\sigma })$, and $(\bar V,\bar \psi )$,
$\bar \psi  = ({\bar y^\sigma })$, be two charts on $Y$, adapted to $S$, such that $V \cap \bar V \ne \O $.
\begin{enumerate}
\item The canonical sections along $S$ satisfy
\begin{align}
{\Lambda ^k}\bar \psi  = \det {\left( {\frac{{\partial {y^i}}}{{\partial {{\bar y}^j}}}}
\right)_{\psi (y)}}_{\bar \psi (y)}{\Lambda ^k}\psi .
\end{align}

\item The differential forms ${d y^\sigma }$ and $d{y^I}$ satisfy
$({\Lambda ^k}\psi )^* d{y^i} = d{y^i}$, $1 \leqslant i \leqslant k$,
$({\Lambda ^k}\psi )^* d{y^\nu } = 0$, $k + 1 \leqslant \nu  \leqslant m$,
$({\Lambda ^k}\psi )^* d{\dot y^I} = 0$.
In particular, on the set $V \cap \bar V$,
\begin{align}
& ({\Lambda ^k}\bar \psi )^*d{{\bar y}^1} \wedge d{{\bar y}^2} \wedge  \ldots  \wedge d{{\bar y}^k} \hfill \nonumber \\
&  \quad  = \det {\left( {\frac{{\partial {{\bar y}^i}}}{{\partial {y^j}}}} \right)_{\psi (y)}}({\Lambda ^k}\psi )*d{y^1} \wedge d{y^2} \wedge  \ldots  \wedge d{y^k}.
\end{align}
\end{enumerate}
\end{lemma}

\section{Grassmann fibrations}
Consider the vector bundle ${\Lambda ^k}TY$ and the subset $\Lambda _0^kTY \subset {\Lambda ^k}TY$,
consisted of {\it non-zero} $k$-vectors. We have an equivalence relation on $\Lambda _0^kTY$
``{\it ${\Xi _1}$ is equivalent with ${\Xi _2}$, if there exists a real number $\lambda  > 0$
such that ${\Xi _1} = \lambda {\Xi _2}$}''.
The quotient set has the structure of a fibration over $Y$, called the {\it Grassmann fibration} of degree $k$,
and denoted by ${G^k}Y$.

To describe the structure of the set ${G^k}Y$, we proceed similarly as in the case of classical projective spaces.
If in a chart $(V,\psi )$, $\psi  = ({y^\sigma })$,
\begin{align}
{\Xi _i} = \frac{1}{{k!}}\Xi _i^{{\sigma _1}{\sigma _2} \ldots {\sigma _k}}
{\left( {\frac{\partial }{{\partial {y^{{\sigma _1}}}}}} \right)_y} \wedge
{\left( {\frac{\partial }{{\partial {y^{{\sigma _2}}}}}} \right)_y} \wedge
\ldots  \wedge {\left( {\frac{\partial }{{\partial {y^{{\sigma _k}}}}}} \right)_y},\quad i = 1,2,
\end{align}
are two nonzero $k$-vectors, then ${\Xi _1}$ is equivalent with ${\Xi _2}$ if and only if in this chart,
$\Xi _1^{{\sigma _1}{\sigma _2} \ldots {\sigma _k}} = \lambda \Xi _2^{{\sigma _1}{\sigma _2}
\ldots {\sigma _k}}$ for some $\lambda  > 0$ and all
${\sigma _1},{\sigma _2}, \ldots ,{\sigma _k}$.
We denote ${V^{{\nu _1}{\nu _2} \ldots {\nu _k}}} = \{ \Xi
\in {({\tau ^k})^{ - 1}}(V) \,|\, {\Xi ^{{\nu _1}{\nu _2} \ldots {\nu _k}}} > 0\} $.
Then a $k$-vector belonging to the set ${V^{{\nu _1}{\nu _2} \ldots {\nu _k}}} \subset \Lambda _0^k TY$
can be expressed by
\begin{align}
&\Xi  = {\Xi ^{{\nu _1}{\nu _2} \ldots {\nu _k}}}{\left( {\frac{\partial }{{\partial {y^{{\nu _1}}}}}}
\right)_y} \wedge {\left( {\frac{\partial }{{\partial {y^{{\nu _2}}}}}} \right)_y} \wedge  \ldots  \wedge
{\left( {\frac{\partial }{{\partial {y^{{\nu _k}}}}}} \right)_y} \hfill \nonumber \\
& \quad  + \frac{1}{{k!}}\sum\limits_{({\tau _1}{\tau _2} \ldots {\tau _k})
\ne ({\nu _1}{\nu _2} \ldots {\nu _k})} {{\Xi ^{{\tau _1}{\tau _2} \ldots {\tau _k}}}}
{\left( {\frac{\partial }{{\partial {y^{{\tau _1}}}}}} \right)_y} \wedge
{\left( {\frac{\partial }{{\partial {y^{{\tau _2}}}}}} \right)_y} \wedge
\ldots  \wedge {\left( {\frac{\partial }{{\partial {y^{{\tau _k}}}}}} \right)_y}
\end{align}
(no summation through ${\nu _1},{\nu _2}, \ldots ,{\nu _k}$).
Denoting by $\operatorname{sgn} {\Xi ^{{\nu _1}{\nu _2} \ldots {\nu _k}}}$
the sign of the component ${\Xi ^{{\nu _1}{\nu _2} \ldots {\nu _k}}}$, we can write
${\Xi ^{{\nu _1}{\nu _2} \ldots {\nu _k}}} = \operatorname{sgn} {\Xi ^{{\nu _1}{\nu _2} \ldots {\nu _k}}} \cdot |
{\Xi ^{{\nu _1}{\nu _2} \ldots {\nu _k}}}|$ and
\begin{align}
& \Xi  = \operatorname{sgn} {\Xi ^{{\nu _1}{\nu _2} \ldots {\nu _k}}} \cdot |{\Xi ^{{\nu _1}{\nu _2}
\ldots {\nu _k}}}|{\left( {\frac{\partial }{{\partial {y^{{\nu _1}}}}}} \right)_y}
\wedge {\left( {\frac{\partial }{{\partial {y^{{\nu _2}}}}}} \right)_y} \wedge  \ldots  \wedge
{\left( {\frac{\partial }{{\partial {y^{{\nu _k}}}}}} \right)_y} \hfill \nonumber \\
& \quad  + \frac{{|{\Xi ^{{\nu _1}{\nu _2} \ldots {\nu _k}}}|}}{{k!}}{\sum {\frac{{{\Xi ^{{\tau _1}{\tau _2}
\ldots {\tau _k}}}}}{{|{\Xi ^{{\nu _1}{\nu _2} \ldots {\nu _k}}}|}}{{\left(
{\frac{\partial }{{\partial {y^{{\tau _1}}}}}} \right)}_y}
\wedge {{\left( {\frac{\partial }{{\partial {y^{{\tau _2}}}}}} \right)}_y} \wedge  \ldots  \wedge
\left( {\frac{\partial }{{\partial {y^{{\tau _k}}}}}} \right)} _y},
\end{align}
with the summation through
$({\tau _1}{\tau _2} \ldots {\tau _k}) \ne ({\nu _1}{\nu _2} \ldots {\nu _k})$. But
$\operatorname{sgn} {\Xi ^{{\nu _1}{\nu _2} \ldots {\nu _k}}} = 1$,
so we see the class of $\Xi $ can be represented as
\begin{align}
& [\Xi ] = {\left( {\frac{\partial }{{\partial {y^{{\nu _1}}}}}} \right)_y}
\wedge {\left( {\frac{\partial }{{\partial {y^{{\nu _2}}}}}} \right)_y} \wedge
\ldots  \wedge {\left( {\frac{\partial }{{\partial {y^{{\nu _k}}}}}} \right)_y} \hfill \nonumber \\
&  \quad  + \frac{1}{{k!}}{\sum {\frac{{{\Xi ^{{\tau _1}{\tau _2} \ldots {\tau _k}}}}}{{{\Xi ^{{\nu _1}{\nu _2}
\ldots {\nu _k}}}}}{{\left( {\frac{\partial }{{\partial {y^{{\tau _1}}}}}} \right)}_y}
\wedge {{\left( {\frac{\partial }{{\partial {y^{{\tau _2}}}}}} \right)}_y} \wedge  \ldots  \wedge
\left( {\frac{\partial }{{\partial {y^{{\tau _k}}}}}} \right)} _y}.
\end{align}
We set for any $\Xi  \in {V^{{\nu _1}{\nu _2} \ldots {\nu _k}}}$
\begin{align}
&  {w^\sigma }(\Xi ) = {y^\sigma }(\Xi ),\quad {w^{{\nu _1}{\nu _2} \ldots {\nu _k}}}(\Xi )
= {{\dot y}^{{\nu _1}{\nu _2} \ldots {\nu _k}}}(\Xi ), \hfill \nonumber \\
&  {w^{{\sigma _1}{\sigma _2} \ldots {\sigma _k}}}(\Xi ) = \frac{{{{\dot y}^{{\sigma _1}{\sigma _2}
\ldots {\sigma _k}}}(\Xi )}}{{{{\dot y}^{{\nu _1}{\nu _2} \ldots {\nu _k}}}(\Xi )}},\quad ({\sigma _1}{\sigma _2}
\ldots {\sigma _k}) \ne ({\nu _1}{\nu _2} \ldots {\nu _k}). \label{form5-2}
\end{align}
Then the pair $({V^{{\nu _1}{\nu _2} \ldots {\nu _k}}},
{\Psi ^{{\nu _1}{\nu _2} \ldots {\nu _k}}})$,
${\Psi ^{{\nu _1}{\nu _2} \ldots {\nu _k}}} = ({w^\sigma }, {w^{{\nu _1}{\nu _2} \ldots {\nu _k}}}, {w^{{\sigma _1}{\sigma _2}
\ldots {\sigma _k}}})$, where the indices satisfy
$({\sigma _1}{\sigma _2} \ldots {\sigma _k}) \ne ({\nu _1}{\nu _2} \ldots {\nu _k})$,
is a chart on $\Lambda _0^kTY$; we call this chart $({\nu _1}{\nu _2} \ldots {\nu _k})$-{\it associated}
with $(V,\psi )$. The pair $({V^{{\nu _1}{\nu _2} \ldots {\nu _k}}}, {W^{{\nu _1}{\nu _2} \ldots {\nu _k}}})$,
${W^{{\nu _1}{\nu _2} \ldots {\nu _k}}} = ({w^\sigma },{w^{{\sigma _1}{\sigma _2} \ldots {\sigma _k}}})$,
$({\sigma _1}{\sigma _2} \ldots {\sigma _k}) \ne ({\nu _1}{\nu _2} \ldots {\nu _k})$,
is a fibred chart on ${G^k}Y$.
Writting formulas (\ref{form5-2}) in a different way, we have the transformation equations
\begin{align}
{w^\sigma } = {y^\sigma },\quad {w^{{\nu _1}{\nu _2} \ldots {\nu _k}}} = {\dot y^{{\nu _1}{\nu _2}
\ldots {\nu _k}}},\quad {w^{{\sigma _1}{\sigma _2} \ldots {\sigma _k}}}
= \frac{{{{\dot y}^{{\sigma _1}{\sigma _2} \ldots {\sigma _k}}}}}{{{{\dot y}^{{\nu _1}{\nu _2}
\ldots {\nu _k}}}}}.
\end{align}
The projection ${\kappa ^k}:{\Lambda ^k}TY \to {G^k}Y$ of ${\Lambda ^k}TY$ onto
$G^k Y$ is the Cartesian projection
$({w^\sigma },{w^{{\nu _1}{\nu _2} \ldots {\nu _k}}}, {w^{{\sigma _1}{\sigma _2} \ldots {\sigma _k}}})
\to ({w^\sigma },{w^{{\sigma _1}{\sigma _2} \ldots {\sigma _k}}})$.
Combining ${\Lambda ^k}({\psi ^{ - 1}}{\iota _{k,m}})$ and
${\kappa ^k}$ we get the canonical lift of ${\psi ^{ - 1}}{\iota _{k,m}}$ to the Grassmann fibration,
\begin{align}
{G^k}({\psi ^{ - 1}}{\iota _{k,m}}) = {\kappa ^k} \circ {\Lambda ^k}({\psi ^{ - 1}}{\iota _{k,m}}).
\end{align}
\begin{lemma}
Let $(V,\psi )$, $\psi  = ({y^\sigma })$, and $(\bar V,\bar \psi )$, $\bar \psi  = ({\bar y^\sigma })$,
be two rectangle charts, adapted to $S$ at a point $y \in Y$.
Suppose that $(V,\psi )$ and $(\bar V,\bar \psi )$ are consistently oriented. Then
\begin{align}
{G_k}({\bar \psi ^{ - 1}}{\iota _{k,m}}) = {G_k}({\psi ^{ - 1}}{\iota _{k,m}}).
\end{align}
\end{lemma}
We set
\begin{align}
{G^k}S = \{ [\Xi ] \in {G^k}Y|[\Xi ] = {G^k}({\psi ^{ - 1}}{\iota _{k,m}})({\operatorname{pr} _{m,k}}\psi (y)),y \in S\} .
\end{align}
We associate to a given chart $(V,\psi )$, $\psi  = ({y^\sigma })$, the induced chart
$({({\tau ^k})^{ - 1}}(V),\Phi )$, $\Phi  = ({y^\sigma },{\dot y^{{\sigma _1}{\sigma _2} \ldots {\sigma _k}}})$,
on ${\Lambda ^k}TY$; the associated charts on the Grassmann fibration ${G^k}Y$ are $(V_0^{{\nu _1}{\nu _2}
\ldots {\nu _k}},{W^{{\nu _1}{\nu _2} \ldots {\nu _k}}})$,
${W^{{\nu _1}{\nu _2} \ldots {\nu _k}}} = ({w^\sigma }, {w^{{\sigma _1}{\sigma _2} \ldots {\sigma _k}}})$,
$({\sigma _1}{\sigma _2} \ldots {\sigma _k}) \ne ({\nu _1}{\nu _2} \ldots {\nu _k})$.
Then it is easily seen that each of the charts $(V_0^{{\nu _1}{\nu _2} \ldots {\nu _k}}, {W^{{\nu _1}{\nu _2}
\ldots {\nu _k}}})$ is adapted to the submanifold ${G^k}S$.
\begin{theorem} \label{th1}
Suppose $S$ is oriented. Then the subset ${G^k}S$ of the Grassmann fibration ${G^k}Y$ is a $k$-dimensional
oriented submanifold, diffeomorphic with $S$.
\end{theorem}

Theorem \ref{th1} allows us to integrate over $k$-dimensional submanifolds of $Y$ directly on the
Grassmann fibration ${G^k}Y$.

\section{Variational functionals depending on submanifolds}
As before, we denote by ${G^k}S$ (resp. ${G^k}\Omega $) the canonical lift of a $k$-dimensional submanifold
$S \subset Y$ (resp. $k$-piece $\Omega  \subset Y$) to the Grassmann fibration ${G^k}Y$.
Denote by ${\Gamma ^k}Y$ the set of all $k$-pieces $\Omega $ of the manifold $Y$.

Let $\eta $ be a $k$-form on the Grassmann fibration ${G^k}Y$.
The form $\eta $ defines a {\it variational functional}
\begin{align}
{\Gamma ^k}Y \ni \Omega  \to {\eta _\Omega }(S) = \int_{{G^k}\Omega } \eta   \in R. \label{vari}
\end{align}
We roughly describe in this paper this construction for $k = 1$,
representing variational functionals of {\it Finsler geometry} in terms of differential forms
(cf. Urban and Krupka~\cite{Krupka-Urban2}).
Consider the tangent bundle ${\Lambda ^1}TY = TY$, a chart $(V,\psi)$,
$\psi  = ({y^\sigma })$, on $Y$, and the associated chart ${({\tau ^1})^{ - 1}}(V)$,
$\Psi  = ({y^\sigma },{\dot y^\sigma })$, on $TY$. A function $F:TY \to R$
is said to be a {\it Finsler fundamental function}, if it satisfies
\begin{align}
F(\lambda \xi ) = \lambda F(\xi )
\end{align}
for all tangent vectors $\xi$ and every positive
$\lambda  \in R$ ({\it homogeneity condition}). The same can be stated in coordinates, requiring that
\begin{align}
F({y^\nu },\lambda {\dot y^\nu }) = \lambda F({y^\nu },{\dot y^\nu }).
\end{align}
\begin{theorem} \label{th2}
\begin{enumerate}
\item For any function $F:TY \to R$, the chart expressions
\begin{align}
\eta  = \frac{{\partial F}}{{\partial {{\dot y}^\nu }}}d{y^\nu }
\end{align}
define a global $1$-form on $TY$.
\item If $F$ satisfies the homogeneity condition, then $\eta $ is projectable on the Grassmann fibration
${G^1}TY$.
\item If $F$ satisfies the homogeneity condition, then for any curve $\zeta :I \to Y$
\begin{align}
({\Lambda ^1}\zeta )*\eta  = (F \circ {\Lambda ^1}\zeta )dt. \label{form5-3}
\end{align}
\end{enumerate}
\end{theorem}
The form $\eta $ (\ref{form5-3}) is known as the {\it Hilbert form} (Chern, Chen and Lam~\cite{ChernChenLam},
Crampin and Saunders~\cite{CrampinSaunders}).
Theorem \ref{th2} (2) characterizes its basic property when F is positive homogeneous:
 namely, in this case the Hilbert form is defined on the {\it Grassmann fibration} ${G^1}TY$.
 One can also easily verify that $\eta$ is a special case of the {\it Lepage-Cartan form}.
 This fact completely determines the behaviour of the variational functional (\ref{vari})
 under variations of submanifolds, extremal submanifolds, and their invariance properties.

\section*{Acknowledgement}
Erico Tanaka acknowledges the support by the JSPS Institutional Program for Young Researcher Overseas Visits,
Palacky University (PrF-2011-022) and Yukawa Institute Computer Facility.
The second author acknowledges support of National Science Foundation of China (grant No. 109320020).
He is thankful to the School of Mathematics, Beijing Institute of Technology,
for kind hospitality and collaboration during his stay in China.
He also acknowledges support from grant 201/09/0981 of the Czech Science Foundation and from the
IRSES project GEOMECH (project no. 246981) within the 7th European Community Framework Programme.


\section{Introduction}

Observations of the large-scale structure in the universe, most notably
of the cosmic microwave background (CMB) anisotropies and the galaxy distribution, are a
powerful probe of both the expansion history and the energy content
of the universe. The precision of these observations is such that in some cases
cosmological constraints on particle physics models that could plausibly account for some of the
energy content are more stringent than limits obtained from direct laboratory measurements. One notable example is the
cosmological bound on the neutrino mass which is currently close to
an order of magnitude more restrictive than existing laboratory
bounds \cite{Lesgourgues:2006nd,Hannestad:2010kz,Wong:2011ip,Lesgourgues:2012uu}. Given this
ability of cosmology to probe neutrino physics, it is of significant
interest to investigate how well other hot dark matter (HDM) models with
similar behaviours can be tested.

One example of particular interest is the case of hot dark matter
axions. In reference~\cite{Hannestad:2003ye} some of the present authors found that
cosmological data do in fact  provide an upper bound on the mass of
hadronic axions comparable to that for  neutrinos. The analysis was
subsequently updated on several occasions to include progressively more recent data
releases~\cite{Hannestad:2005df,Hannestad:2007dd,Hannestad:2008js,Hannestad:2010yi},
again with similar conclusions.   Other authors also found comparable
results~\cite{Melchiorri:2007cd}. Here, we provide another update incorporating
 the most recent measurement of the CMB temperature anisotropies by the Planck mission~\cite{Ade:2013xsa,Ade:2013lta}, as well as
the matter power spectrum extracted from
the seventh data release of the Sloan Digital Sky Survey (SDSS-DR7)~\cite{Reid:2009xm},
and the present-day Hubble parameter obtained from the most recent recalibration of
the Hubble Space Telescope (HST) Key Project sample as part of the Carnegie Hubble Program~\cite{Freedman:2012ny}.

The exercise of extracting information about low-mass particles from precision
cosmological data has received a new twist in the  past few
years because combinations of different data sets are not necessarily
consistent with the minimal (``vanilla'') $\Lambda$CDM cosmology.
In particular, many pre-Planck observations prefer at the $2 \sigma$ level a ``dark radiation'' component equivalent to an additional
$\sim 0.5$ massless neutrino species (see, e.g., \cite{Hou:2012xq,Archidiacono:2013lva,Calabrese:2013jyk,Abazajian:2012ys}
and references therein).
The Planck data alone and in combination with baryon acoustic oscillations
do not support the presence of a dark radiation component.  However, the preferred value of the Hubble parameter $H_0$ is across the board significantly lower than that from direct measurements in our local neighbourhood~\cite{Freedman:2012ny}; if these measurements were to be brought in line with one another, then the presence of a dark radiation equivalent to $\sim 0.6 \pm 0.25$ could achieve the required result through its degeneracy with $H_0$~\cite{Ade:2013lta}.  However, hovering around the $2 \sigma$ level the case for dark radiation, pre-Planck or post-Planck, is certainly not strong.   Indeed, it could be as much a signature for new physics as it is a signature for poorly controlled systematics in some cosmological data sets.  Nonetheless, since we are dealing here with low-mass particles which do contribute to the primordial radiation density,
some additional care is required when dealing with local~$H_0$ measurements.

The rest of the paper is organised as follows. In section~\ref{sec:axions} we briefly discuss the hot dark matter axion scenario.  We describe the cosmological parameter space to be analysed in section~\ref{sec:model}, while section~\ref{sec:data} summarises the data sets.  Section~\ref{sec:results} contains our results, notably,  bounds on the axion mass $m_a$ and the neutrino mass sum $ \Sigma \, m_\nu$ from various data combinations, as well as a brief discussion of the impact of hot dark matter on the CMB observables and the corresponding degeneracies with $H_0$.  We conclude in section~\ref{sec:conclusions}, followed by an appendix~\ref{sec:appendix}  in which we discuss the degeneracies of $m_a$ and $\Sigma \, m_\nu$ in more detail.

\section{Hot dark matter axions}                    \label{sec:axions}

The Peccei--Quinn solution of the CP problem of strong interactions
predicts the existence of axions, low-mass pseudoscalars that are
very similar to neutral pions, except that their mass and interaction
strengths are suppressed relative to the pion case
by a factor of order $f_\pi/f_a$, where
$f_\pi\approx 93$~MeV is the pion decay constant, and $f_a$ is a large
energy scale known as the axion decay constant or Peccei--Quinn
scale~\cite{Peccei:2006as,Kim:2008hd,Kawasaki:2013ae}. Specifically,
 the axion mass~is
\begin{equation}\label{eq:axmass}
 m_a=\frac{z^{1/2}}{1+z}\,\frac{f_\pi}{f_a} m_\pi
 =\frac{6.0~{\rm eV}}{f_a/10^6~{\rm GeV}}\,,
\end{equation}
where $z=m_u/m_d$ is the mass ratio of the up and down quarks. We
shall follow the previous axion literature and assume the canonical
value $z=0.56$. We note that variation in the range
0.38--0.58 is possible~\cite{Beringer:1900zz}, although this uncertainty has no
strong impact on our discussion. Henceforth we shall mostly use $m_a$
as our primary phenomenological axion parameter, instead of the more
fundamental parameter~$f_a$.

Axions with masses in the $10~\mu{\rm eV}$ region are well motivated
as a cold dark matter (CDM) candidate, albeit with considerable uncertainty
as regards the precise mass value~\cite{Kawasaki:2013ae,Sikivie:2006ni,Wantz:2009it,Hiramatsu:2012gg}.
Such CDM axions would have been produced non-thermally by the
re-alignment mechanism and, depending on the cosmological scenario,
by the decay of axion strings and domain walls. In addition, a sizeable hot
axion population can be produced by thermal processes if the Peccei--Quinn scale
is sufficiently low~\cite{Chang:1993gm,Turner:1986tb,Masso:2002np,Graf:2010tv}.  For $f_a\alt10^8$~GeV, or, equivalently, $m_a\agt60~{\rm meV}$,
axions attain thermal equilibrium at the QCD phase transition or later, and contribute to the cosmic radiation density and subsequently
to the cosmic hot dark matter  along with massive neutrinos.

In principle, the parameter region $f_a\mbox{\;\raisebox{.3ex 10^9$~GeV ($m_a\agt10~{\rm meV}$) has already been excluded
by the SN~1987A neutrino burst
duration~\cite{Mayle:1987as,Raffelt:1987yt,Turner:1987by,Raffelt:2006cw}
based on the axion--nucleon interaction. Furthermore, axions in the same parameter region
would also contribute strongly to neutron-star
cooling~\cite{Umeda:1997da,Keller:2012yr} and engender a cosmic diffuse
supernova axion background (DSAB)~\cite{Raffelt:2011ft}. However,
the sparse data sample, our poor understanding of the nuclear medium
in the supernova interior, and simple prudence suggest that one
should not base far-reaching conclusions about the existence of
axions in this parameter range on a single argument or experiment
alone.  For those axion models in which the axion couples nonvanishingly
to the electron (DFSZ-type axion models~\cite{Zhitnitsky:1980tq,Dine:1981rt}),
numerous other stellar energy loss limits based on the axion--electron coupling
also exclude axions with masses exceeding some 10~meV~\cite{Raffelt:1985nj,Raffelt:1994ry,Isern:1992,Corsico:2012ki,%
Corsico:2012sh,Isern:2008nt,Isern:2008fs,Barth:2013sma}.  We therefore
do not consider these (well-constrained) models, but rather focus on hadronic models (KSVZ-type models
\cite{Kim:1979if,Shifman:1979if}) in which the axion does not directly
couple to ordinary quarks and leptons, and in so doing bypasses constraints from stellar evolution arguments.
Hadronic axion models are also simple to analyse in that all axion properties depend on $f_a$ alone.

If axions do not couple to charged leptons, the main thermalisation
process in the post-QCD epoch is $a+\pi\leftrightarrow\pi+\pi$
\cite{Chang:1993gm}. The axion--pion interaction is given by the
Lagrangian~\cite{Chang:1993gm}
\begin{equation}\label{eq:axionpionlagrangian}
{\cal L}_{a\pi}=\frac{C_{a\pi}}{f_\pi f_a}\,
\left(\pi^0\pi^+\partial_\mu\pi^- +\pi^0\pi^-\partial_\mu\pi^+
-2\pi^+\pi^-\partial_\mu\pi^0\right)
\partial_\mu a\,,
\end{equation}
where, in hadronic axion models, the coupling constant is
$C_{a\pi}=(1-z)/[3\,(1+z)]$. Based on this interaction, the axion
decoupling temperature $T_{a,{\rm dec}}$ was calculated earlier in~\cite{Hannestad:2005df},
from which follows the present-day axion number density,
\begin{equation}
n_a=\frac{g_{*S}(T_{\rm today})}{g_{*S}(T_{a,{\rm dec}})} \times \frac{n_\gamma}{2},
\end{equation}
where $g_{*S}(T)$ denotes the effective entropy degrees of freedom, and $n_\gamma$ is the present-day photon number density.
The top panel of figure~\ref{fig:nama} shows $n_a$ as a function of the axion mass~$m_a$.
Before thermal axions become nonrelativistic by cosmological redshift, they
contribute to the total radiation density an amount
\begin{equation}
\rho_a = \frac{2}{3}\(\frac{3}{2}\frac{n_a}{n_\nu}\)^\frac{4}{3}\, \frac{7}{8}\(\frac{4}{11}\)^{\frac{4}{3}} \rho_\gamma
=
\Delta N_{\rm eff} \ \frac{7}{8}\(\frac{4}{11}\)^{\frac{4}{3}} \rho_\gamma
\label{eq:neff}
\end{equation}
equivalent to $\Delta N_{\rm eff}$ additional light neutrino species,
where $n_\nu=113~{\rm cm}^{-3}$ is the present-day neutrino number
density in one flavour, and $\rho_\gamma$ is the photon energy density.
We stress that
equation~(\ref{eq:neff}) is strictly valid only when the axion temperature far exceeds $m_a$.
For order 0.1--1~eV-mass axions this corresponds to a time well before photon decoupling;
the phenomenology of such axions in the context of the CMB anisotropies and the matter power spectrum is closely linked to the axion mass scale, and
{\it cannot} be reparameterised in terms of
$\Delta N_{\rm eff}$ {\it massless} neutrino species.
Indeed, as we shall show in section~\ref{sec:results} and appendix~\ref{sec:appendix}, the degeneracy of $m_a$ especially in the direction of the Hubble parameter is vastly different from that of $\Delta N_{\rm eff}$  massless neutrinos.

\begin{figure}
\centering
\begin{tabular}{c}
\includegraphics[scale=0.4]{na_vs_ma.eps}\\
\includegraphics[scale=0.4]{oa_vs_ma_loglog.eps}
\end{tabular}
\caption{{\it Top}: Axion number density $n_a$ as a function of the axion mass $m_a$.
{\it Bottom}: Present-day axion energy density $\omega_a$ as a function of $m_a$.}
\label{fig:nama}
\end{figure}

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

We perform the analysis in the context of a
$\Lambda$CDM model extended to include hot dark matter neutrinos and axions, and
vary a total of eight parameters:
\begin{equation}\label{eq:model}
{\bm \theta} = \{\omega_{\rm cdm},\omega_{\rm b},\theta_{\rm s},\tau,
\ln(10^{10}A_{s}),n_{s},
\Sigma\,m_\nu,m_a\}.
\end{equation}
Here, $\omega_{\rm cdm} \equiv \Omega_{\rm cdm} h^2$ and $\omega_{\rm b}
\equiv \Omega_{\rm b} h^2$ are the present-day physical CDM and baryon densities
respectively, $\theta_{\rm s}$ the angular the sound horizon, $\tau$  the optical
depth to reionisation, and $\ln(10^{10}A_{s})$ and $n_s$ denote respectively the amplitude and spectral index
of the initial scalar fluctuations.  We use the  neutrino mass sum $\Sigma\,m_\nu$ and the axion mass $m_a$ as fit parameters, noting that these can be easily converted to their corresponding present-day energy density via
\begin{eqnarray}
\omega_\nu &=& \frac{\Sigma \, m_\nu}{94 \ {\rm eV}}, \label{eq:omeganu} \\
\omega_a & = & \frac{2}{3}\(\frac{m_a}{94 \ {\rm eV}}\)\(\frac{T_a}{T_\nu}\)^3, \label{eq:omegaa}
\end{eqnarray}
where $T_\nu$ and $T_a$ are the neutrino and the axion present-day temperature respectively,
and $\omega_a$ is shown in the bottom panel of figure~\ref{fig:nama} as a function of $m_a$.
Table~\ref{tab:priors} shows the priors adopted in the analysis.
We stress that $\Sigma \, m_\nu$ is always let free to vary in our analysis, because neutrino masses are an experimentally established fact, not some novel untested physics that can be discarded at will.

\begin{table}[t]
\caption{Priors for the cosmological fit parameters considered in this work. All priors are uniform (top hat) in the given intervals.}
\label{tab:priors}
\begin{center}
\begin{tabular}{lc}
\hline
 Parameter & Prior\\
\hline
$\omega_{\rm b}$ & $0.005 \to 0.1$\\
$\omega_{\rm cdm}$ & $0.001 \to 0.99$\\
$\theta_{\rm s}$ & $0.5 \to 10$\\
$\tau$ & $0.01 \to 0.8$\\
$n_{s}$ & $0.9 \to 1.1$\\
$\ln{(10^{10} A_{s})}$ & $2.7 \to 4$\\
$\Sigma\,m_\nu$ [eV] &  $0 \to 7$\\
$m_a$ [eV] &  $0 \to 7$\\
\hline
\end{tabular}
\end{center}
\end{table}

For convenience we assume the neutrino mass sum $\Sigma \, m_\nu$ to be shared equally among three standard model neutrinos.  Note however that implementing
any other mass spectrum will return essentially the same conclusions,
since current observations are primarily sensitive to $\Sigma \, m_\nu$, not to the exact mass splitting.  In fact, discerning the neutrino mass spectrum---and thereby distinguishing between the normal and the inverted neutrino mass hierarchies---is unlikely to be possible even with the extremely precise measurements
from the future {\sc Euclid} mission~\cite{Hamann:2012fe}.

\section{Data and analysis}\label{sec:data}

We consider three types of measurements: temperature and polarisation power spectra of the CMB anisotropies, the large-scale matter power spectrum, and direct measurements of the local Hubble expansion rate.  These are discussed in more detail below.  To these data sets we apply a Bayesian statistical inference analysis using the publicly available Markov Chain Monte Carlo
parameter estimation package {\sc CosmoMC}~\cite{Lewis:2002ah} coupled to the CAMB~\cite{Lewis:1999bs} Boltzmann solver modified to accommodate two hot dark matter components.  With the exception of the local Hubble parameter measurement, the likelihood routines and the associated window functions are supplied by the experimental collaborations.

\subsection{CMB anisotropies}

Our primary data set is the recent measurement of the CMB temperature (TT) power spectrum by the Planck mission~\cite{Planck:2013kta}, which we implement into our
likelihood analysis following the procedure reported in~\cite{Ade:2013lta}.
This data is supplemented by measurements of the CMB polarisation from the WMAP nine-year data release~\cite{Bennett:2012fp}, in the form of
an autocorrelation (EE) power spectrum at $2<\ell<32$ and a cross-correlation (TE) with the Planck temperature measurements in the same multipole range.  We denote this supplement ``WP''.

\subsection{The matter power spectrum}

We use the matter (halo) power spectrum (HPS) determined from the
luminous red galaxy sample of the seventh data release of the SDSS-DR7~\cite{Reid:2009xm}. The full data set
consists of 45 data points, covering wavenumbers from $k_{\rm min} =
0.02\ h {\rm Mpc}^{-1}$ to $k_{\rm max} = 0.2\ h {\rm Mpc}^{-1}$. We
fit this data set following the procedure of
reference~\cite{Reid:2009xm}, using an adequately smeared power
spectrum to model nonlinear mode-coupling constructed according to the
method described in reference~\cite{Hamann:2010pw}.

\subsection{Local Hubble parameter measurements}

We further impose a constraint on the present-day Hubble parameter based on measurements of the universal expansion rate in our local neighbourhood.  The most
recent value released as part of the Carnegie Hubble Program~\cite{Freedman:2012ny} is
\begin{equation}
H_0 = [74.3 \pm 1.5 \, {\rm (stat)}\,\pm 2.1 \, {\rm (syst)}]~{\rm km}~{\rm s}^{-1}~{\rm Mpc}^{-1},
\end{equation}
which uses observations from the Spitzer Space Telescope to calibrate the Cepheid distance scale in the HST Key Project sample.
Following the analysis of reference~\cite{Ade:2013lta}, we add the statistical and systematic errors in quadrature,
assuming both uncertainties to be Gaussian-distributed.

\section{Results and discussions}\label{sec:results}

The main results of our inference analysis are presented in table~\ref{tab:results} and figure~\ref{fig:mamnu}, which show, respectively, the
one-dimensional marginal limits on $\Sigma \, m_\nu$ and $m_a$ derived from various data combinations, and the corresponding two-dimensional marginal probability density contours in
the $(\Sigma \, m_\nu,m_a)$-subspace.

\begin{table}
\center
{
\caption{One-dimensional marginal 95\% bounds on the neutrino mass sum and the axion mass from
various data combinations.\label{tab:results}} \medskip
\centering
\begin{tabular}{lcc}
\hline
Data set  & $\Sigma\,m_\nu$ [eV] & $m_a$ [eV]\\
\hline
Planck+WP  & $ <0.84$ & $<1.01$    \\
Planck+WP+HPS   & 0.007--0.48  & $<0.86$    \\
Planck+WP+$H_0$ & $<0.25$ & $<0.78$ \\
Planck+WP+HPS+$H_0$ & $<0.27$  & $<0.67$ \\
\hline
\end{tabular}}
\end{table}

\subsection{Axion mass limit from CMB alone}\label{sec:mafromcmb}

Most strikingly, while it was previously not possible to obtain an upper limit on the axion mass from the
CMB anisotropies alone, we now find that the combination of Planck+WP is able to constrain $m_a$
quite stringently to $m_a<1.01$~eV~(95\% C.L.).  This is in contrast with the CMB-only upper limit on the neutrino mass sum $\Sigma \, m_\nu$, which has been hovering around 1~eV since the WMAP five-year data release: here, we find  $\Sigma\,m_\nu<0.78$~eV (95\% C.L.) after marginalising over the unknown axion mass,  somewhat stronger than the pre-Planck value of 1.19~eV \cite{Hannestad:2010yi}.

\begin{figure}[t]
\centering
\includegraphics[scale=0.45]{msc_mnu.eps}
\caption{Two-dimensional marginal 68\% and 95\% probability density contours in the $(\Sigma\,m_\nu,m_a)$-plane derived from Planck temperature and WMAP polarisation (Planck+WP)
measurements (blue shading) and in combination with the halo power spectrum (red shading).}
\label{fig:mamnu}
\end{figure}

The fundamental and irreducible effect of 0.1--1~eV-mass hot dark matter neutrinos and axions on the CMB temperature anisotropies are similar in that both particle species
 transit from a ultrarelativistic state to a nonrelativistic one in close proximity to the photon decoupling era, as well as to the epoch of matter--radiation equality.  This transition impacts strongly on the evolution of the radiation-to-matter energy density ratio, which in turn affects the evolution of the spacetime metric perturbations and ultimately the photon temperature fluctuations.  To quantify the precise effect of this transition on the CMB temperature power spectrum naturally necessitates solving the relevant Boltzmann equations.  However, a general rule of thumb is that, assuming a fixed present-day total matter density,
those perturbation modes that enter the horizon  before the hot dark matter particle becomes fully nonrelativistic will suffer more strongly from decay of the metric perturbations in a mixed CDM+HDM cosmology than in a pure CDM scenario.  This is translated firstly into an additional ``uplift'' of the primary CMB temperature fluctuations at last scattering on the same length scales.  Secondly, the same decaying metric perturbations also alter the so-called early integrated Sachs--Wolfe (ISW) effect, which contributes to the CMB temperature anisotropies primarily around the first acoustic peak.

Following from this reasoning, what distinguishes between a hot dark matter neutrino population and an axion one of the same present-day (nonrelativistic) energy density is essentially two-fold.
Firstly, the axion population is always colder than the neutrino population by virtue of having decoupled from the cosmic plasma at an earlier time.  Being a scalar degree of freedom the axion must also have a larger particle mass than the neutrino in order to make up the same present-day energy density [see equations~(\ref{eq:omeganu})and~(\ref{eq:omegaa})].  Together these push the axion's relativistic-to-nonrelativistic transition to a higher redshift, so that its irreducible features on the CMB anisotropies occur on smaller scales (or higher multipoles).
Secondly, while the axion is a new source of energy density {\it in addition} to the energy content of the $\Lambda$CDM model, in endowing neutrinos with masses we are effectively ``stealing'' from the model's radiation content.  This effect combined with the axion's relative coldness and heaviness means that a hot axion population will always cause less deviation to the radiation-to-matter energy density ratio, and therefore leave a weaker signature on the CMB than a neutrino population of the same present-day energy density.  Figure~\ref{fig:cls} shows the CMB temperature angular power spectra for a neutrino and an axion HDM scenario, juxtaposed against their $\Lambda$CDM counterpart.

\begin{figure}[t]
\centering
\includegraphics[scale=0.4]{clsmatlab}
\includegraphics[scale=0.4]{clsadjustedmatlab}
\caption{CMB temperature angular power spectra for three different cosmological
models, all with the same total dark matter density $\omega_{\rm dm} =
\omega_{\rm cdm} + \omega_{\rm \nu} + \omega_{\rm a} = 0.112$.  The blue/dot-dash line corresponds to the reference $\Lambda$CDM model ($\omega_{\rm cdm}=0.112$), the green/dashed line an axion HDM scenario with $m_a=2.4$~eV ($\omega_{\rm cdm}=0.099$,  $\omega_a=0.013$), and the red/solid line the case of three degenerate neutrinos with masses
summing to $\Sigma \, m_\nu=1.2$~eV  ($\omega_{\rm cdm}=0.099$, $\omega_\nu=0.013$).
In the top panel we hold all other cosmological parameters fixed, while in the bottom panel we further adjust $H_0$
in order to match the peak positions---the corresponding $H_0$ values are, in units of ${\rm km} \ {\rm s}^{-1} \ {\rm Mpc}^{-1}$,
70, 72, and 64.5 for the $\Lambda$CDM, the axion and the neutrino model respectively.
}
\label{fig:cls}
\end{figure}

Thus what has enabled Planck to constrain the axion mass for the first time from CMB data alone can be traced to its significantly better measurement of the higher multipole moments from the third acoustic peak onwards (recall that WMAP is only cosmic-variance limited up to $\ell  \sim 500$, roughly the second acoustic peak), which rules out much of the axion parameter space in the $>1$~eV region manifesting on these scales.  In contrast, the neutrino mass could already be constrained to order 1~eV by CMB data alone
since the WMAP three-year data release \cite{Spergel:2006hy}
via the early ISW effect on the third-to-first acoustic peak height ratio; new, high multipole measurements from Planck offer no significant improvement because there are no new unique signatures from sub-eV-mass neutrinos on these scales.  See also discussions in appendix~\ref{sec:appendix}.

\subsection[Correlation with $H_0$]{Correlation with \boldmath{$H_0$}}\label{sec:correlation}

Combining CMB measurements with the halo power spectrum from SDSS DR-7 tightens the axion mass bound to
$m_a < 0.86$~eV~(95\%~C.L.), comparable to the pre-Planck value of 0.82~eV~\cite{Hannestad:2010yi}.
The neutrino mass bound likewise improves to $\Sigma \, m_\nu < 0.48$~eV~(95\%~C.L.).  Note that formally the combination of Planck+WP+HPS yields also
a  95\% lower limit of $\Sigma \, m_\nu> 0.007$~eV.  However, because neutrino oscillation experiments already limit the neutrino mass sum to $> 0.06$~eV, this cosmological lower limit is of little phenomenological consequence.

An interesting effect arises when we consider also local measurements of the Hubble parameter $H_0$;
while the neutrino mass bound continues to tighten to $\Sigma \, m_\nu < 0.25$~eV~(95\% C.L.), inclusion of the local $H_0$ value leads only to a marginal 10\% improvement
on the axion mass bound.  This effect can be understood from an inspection of
figures~\ref{fig:mnuH0} and~\ref{fig:maH0}, which show, respectively, the two-dimensional marginal density contours in
the $(\Sigma \, m_\nu,H_0)$- and $(m_a,H_0)$-subspaces.  Clearly, a strong anti-correlation exists between $\Sigma \, m_\nu$ and $H_0$, so that adding to the analysis the local $H_0$ value---which is higher than the Planck-inferred value---tends to tighten, perhaps artificially, the neutrino mass bound.  In contrast,
no such correlation exists between $m_a$ and $H_0$, making the inference of $m_a$ independent of the value of~$H_0$.  Thus we conclude that the axion hot dark matter bound is robust against the
unresolved discrepancy between the Planck-inferred  $H_0$ value and that measured in our local neighbourhood.
Figure~\ref{fig:mamnuH0} shows the two-dimensional marginal probability density contours in the $(\Sigma \, m_\nu, m_a)$-subspace with the inclusion of local $H_0$ data.

\begin{figure}[t]
\centering
\includegraphics[scale=0.45]{mnu_H0.eps}
\caption{Same as figure~\ref{fig:mamnu}, but in the $(\Sigma\,m_\nu,H_0)$-plane.}
\label{fig:mnuH0}
\end{figure}

\begin{figure}[t]
\centering
\includegraphics[scale=0.45]{msc_H0.eps}
\caption{Same as figure~\ref{fig:mamnu}, but in the $(m_a,H_0)$-plane.}
\label{fig:maH0}
\end{figure}

A deeper question now is, what causes $\Sigma \, m_\nu$ and $m_a$ to exhibit different correlations with~$H_0$ in the first place.  We defer a more detailed discussion of the parameter degeneracies in the CMB temperature anisotropies to appendix~\ref{sec:appendix}, but note here that at the most basic level, any degeneracy (or lack thereof) of $\sum m_\nu$ and $m_a$ with $H_0$ can be traced to the parameter dependence of the angular sound horizon $\theta_{\rm s}  \equiv r_{\rm s} (z_\star)/D_{\rm A} (z_\star)$,
where
\begin{eqnarray}
r_{\rm s} (z_\star) &\equiv& \int^\infty_{z_\star} {\rm d} z \ \frac{c_{\rm s} (z)}{H(z)} , \label{eq:rs}\\
 D_{\rm A}(z_\star) & \equiv&\int^{z_\star}_0 \frac{{\rm d}z}{H(z)} \label{eq:da}
\end{eqnarray}
are the sound horizon at decoupling and the angular diameter distance to the last scattering surface respectively, $c_{\rm s}(z)$ is the sound speed in the tightly-coupled photon--baryon fluid, and $z_\star \simeq 1100$ is the photon decoupling redshift.  The angular sound horizon determines the location of the acoustic peaks in the CMB angular power spectrum, and is arguably the most well-measured CMB observable.

\begin{figure}[t]
\centering
\includegraphics[scale=0.45]{msc_mnu_addHSTprior.eps}
\caption{Same as figure~\ref{fig:mamnu}, but including the $H_0$ value determined from local measurements.}
\label{fig:mamnuH0}
\end{figure}

Clearly, the presence of the Hubble expansion rate $H(z)$ in the expressions~(\ref{eq:rs}) and~(\ref{eq:da}) indicates {\it a priori} a three-way degeneracy between the neutrino (or axion) energy density (and hence mass), the CDM energy density $\omega_{\rm dm}$,%
\footnote{The baryon density in principle also affects $H(z)$.  However, $\omega_{\rm b}$ can be pinned down by the CMB odd-to-even peak ratios.  We therefore consider it to be a well-measured parameter.}
and the present-day Hubble parameter~$H_0$,%
\footnote{In a flat spatial geometry, the dependence of $H(z)$ on the vacuum energy $\Omega_\Lambda$ can be reparameterised in terms of $H_0$.}
meaning that for any choice of $\Sigma \, m_\nu$ or $m_a$, we can tweak $\omega_{\rm dm}$ and/or $H_0$ to maintain the same $\theta_{\rm s}$ value.
However, the amount of tweaking required for the axion scenario is generally less than that needed by massive neutrinos.  This is because, as discussed in section~\ref{sec:mafromcmb}, an axion population generally causes less deviation to the evolution of the radiation and matter energy densities than a neutrino population of the same present-day energy density.  Indeed, assuming a fixed present-day total matter density, figure~\ref{fig:cls} shows that while the presence of a $\Sigma \, m_\nu=1.2$~eV hot neutrino population causes a noticeable shift in the peak locations that can be restored by lowering the $H_0$ value to $64.5~{\rm km} \ {\rm s}^{-1} \ {\rm Mpc}^{-1}$ (compared to
the reference $\Lambda$CDM value of $70~{\rm km} \ {\rm s}^{-1} \ {\rm Mpc}^{-1}$),
the shift due to a $m_a=2.4$~eV axion population is of a much lesser extent and requires only a small {\it positive} tuning of $H_0$ to $72~{\rm km} \ {\rm s}^{-1} \ {\rm Mpc}^{-1}$.
As we shall demonstrate in appendix~\ref{sec:appendix}, the same negative correlation with $H_0$ remains true for the neutrino case even after accounting for uncertainties in the total matter density, while for the axion scenario the small tuning of $H_0$ may be negative or positive depending on the choice of $\omega_{\rm dm}$ so that overall no correlation exists between $m_a$ and $H_0$.

\section{Summary and conclusions}\label{sec:conclusions}

We have re-examined cosmological bounds on thermal axions based on
the recent CMB temperature anisotropy measurement provided by the Planck mission, as well as other types of cosmological
observations.
For the first time CMB data alone provide a
restrictive limit on the axion mass of $m_a<1.01$~eV (95\% C.L.), which
improves to 0.86~eV with the inclusion of the SDSS matter power spectrum, and to
0.78~eV when combined with local measurements of~$H_0$.
The changes brought about by combining CMB data with other observations are obviously fairly minor, suggesting that
the systematic uncertainties in our limits are small.  In particular,
the unsettled discrepancy between the Planck-inferred $H_0$ value and that measured directly in our local neighbourhood
does not strongly affect the axion limit. This is to be contrasted with the corresponding limit on the neutrino mass sum, which sees a
dramatic reduction from $\Sigma \, m_\nu < 0.84$~eV (95\% C.L.) from CMB alone to $0.25$~eV when combined with local $H_0$ measurements.

Note that when inferring the axion mass bound, we always marginalise over the unknown
neutrino masses which contribute an unavoidable but unknown HDM fraction to the
universe; in this sense, HDM contributions by neutrinos and axions are not
alternatives, since the former must always be present.  However,  unless neutrinoless double-beta decay eventually shows a
signal, laboratory experiments are unlikely to provide neutrino mass information on the cosmologically relevant
scale in the foreseeable future.  Thus marginalisation over $\Sigma \, m_\nu$ is arguably the only consistent approach to cosmological parameter estimation, and
this situation is unlikely to change any time soon.

In comparison with the pre-Planck state of affairs, although we now have an axion mass bound
from post-Planck CMB data alone, the overall limit has not moved significantly away from order~1~eV.
This situation is of interest for solar axion searches
by the Tokyo axion helioscope~\cite{Moriyama:1998kd,Inoue:2002qy,Inoue:2008zp} and by the CAST
experiment at
CERN~\cite{Zioutas:2004hi,Andriamonje:2007ew,Arik:2008mq,Arik:2011rx}.
The axion--photon conversion efficiency in a dipole magnet directed
towards the sun strongly degrades if the photon--axion mass difference
is too large. This mass-mismatch can be overcome by filling the
conversion pipe with a buffer gas at variable pressure that provides
the photons with an adjustable refractive mass. The final CAST
search range using $^3$He as buffer gas has reached a search mass of
1.17~eV, which is hard to push up any further~\cite{Arik:2013}. In this
sense, cosmological and helioscope axion mass limits are still largely
complementary.

In the future, further improvements in HDM bounds on the cosmological front
will likely come from the next generation of large-scale galaxy and cluster surveys, notably the ESA
{\sc Euclid} mission~\cite{Laureijs:2011mu}.  {\sc Euclid} has the potential to probe neutrino masses
with a sensitivity of $\sigma(\Sigma \, m_\nu) \simeq 0.01$~eV, sufficient to measure at $5\sigma$ the minimum mass sum of $\Sigma \, m_\nu \simeq 0.05$~eV
established by neutrino flavour oscillation experiments~\cite{Hamann:2012fe}.  A similar sensitivity to HDM axions should likely be possible.

\section*{Acknowledgements}

We thank Jan Hamann for help with the initial setup of the Planck
likelihood code, and the anonymous referee for constructive suggestions that improve the clarity of the manuscript.
We acknowledge use of computing resources from the
Danish Center for Scientific Computing (DCSC), and
partial support by the Deutsche Forschungsgemeinschaft through grant
No.~EXC~153 and by the European Union through the Initial Training
Network ``Invisibles,'' grant No.\ PITN-GA-2011-289442.