text:
\section{Introduction}
Masses of galactic bulges and their central black holes (BHs) follow a
tight relation in the local Universe \citep[e.g.][]{marc03,haer04} with
only 0.3~dex scatter -- strong evidence for a coupled formation and
evolution of galaxies and BHs. The source of this coupling is unclear, but
feedback mechanisms have been proposed involving the central potential
well depth regulating BH accretion, or more violent feedback from active
galactic nuclei (AGN) into their host galaxies
\citep[e.g.][]{hopk06c,some08,menc08}. While these scenarios potentially
provide ingredients for acquiring consensus with observations, all such
models include ad hoc assumptions and do not work from first
principles. Empirical constraints are urgently needed to investigate the
actual physical processes involved in the coupled evolution.

One strong constraint is the evolution of the
$M_\mathrm{BH}$--$M_\mathrm{*,bulge}$-relation over time.  While
circumstantial evidence grows that the value of
$M_\mathrm{BH}$/$M_\mathrm{*,bulge}$ was larger at earlier cosmic times
\citep{peng06a,peng06b,treu07,woo08,walt04,riec08a,riec08b,riec09},
studies are subject to biases \citep{laue07} and better statistics are
required to investigate where in $M_\mathrm{BH}$, or when in cosmic time, a
turnoff from the local $M_\mathrm{BH}$--$M_\mathrm{*,bulge}$-relation
occurs.

Broad-line AGN and their host galaxies are the only systems at higher
redshifts in which both the mass of the galaxy or its bulge as well as its
central BH can be estimated. Here we set constraints on an evolving
$M_\mathrm{BH}$--$M_\mathrm{*,bulge}$-relation by computing optical
color-based stellar masses (from HST/ACS and HST/NICMOS) and combine them
with virial $M_\mathrm{BH}$ (from Magellan/IMACS and zCOSMOS/VLT/VIMOS)
for 10 AGN in the redshift interval $1.06<z<1.92$, 3.2--5.5~Gyrs after the
Big Bang.

Throughout we use AB zero-points unless otherwise noted and a cosmology of
$H_0=70\,\mathrm{km\,s^{-1} Mpc^{-1}}$, $\Omega_M=0.3$, and
$\Omega_\Lambda=0.7$.

\section{Database}
In order to control selection effects on $M_\mathrm{*}$, $M_\mathrm{BH}$
or any special relation between these two, we require a transparent sample
definition.  Our selection of type-1 AGN is based on X-ray detection in
the XMM-COSMOS survey \citep{hasi07,capp09} and subsequent identification
of their optical counterparts \citep{brus07}. Classification as type-1 AGN
for this study uses both spectroscopic identification of
broad emission lines in the Magellan/IMACS \citep{trum09a} and zCOSMOS
\citep{lill07} surveys as well as photometric classification using the
long spectral baseline SED covered in COSMOS \citep[Capak et al., in
  prep.]{capa07,salv09,ilbe09}. As such, only objects with
high-confidence classification enter our sample.

For $\sim$550 type-1 AGN selected this way we require the
following data, resulting in a random subsample: (1) coverage by ACS
F814W \citep{scov07b,koek07}, and (2) NICMOS
parallels\footnote{http://irsa.ipac.caltech.edu/data/COSMOS/images/nicmos/},
as well as (3) spectra for virial $M_\mathrm{BH}$-estimates from the MgII
broad emission line.

The limiting factors are the relative NIC3 coverage of 6.4\% of the ACS
area ($\sim$20 AGN) and the status of ongoing spectroscopic surveys. Black
hole masses have been calculated by \citet{trum09b} from IMACS spectra and
by \citet{merl09} using zCOSMOS/VIMOS. For this letter we have BH masses
available for 10 AGN, spanning the redshift range $1<z<2$. Two AGN were
observed by both instruments -- the $M_\mathrm{BH}$-estimates are
consistent within 0.2~dex in both cases. Sample information is listed in
Table~\ref{tab:all}, including BH masses and galaxy parameters derived
below.

\begin{deluxetable*}{crcccccccc}
\tabletypesize{\scriptsize}
\tablewidth{0pt}
\tablecaption{Sample summary\label{tab:all}}
\tablecolumns{10}
\tablehead{
 \colhead{}&\colhead{} &\colhead{}&\colhead{$M_\mathrm{BH}$\tablenotemark{b}}&\colhead{} &
 \colhead{$F814W_\mathrm{host}$}&\colhead{$F160W_\mathrm{host}$}
 &\colhead{S\'ersic $n$\tablenotemark{d}}&\colhead{$(B-V)_\mathrm{host}^\mathrm{rest}$}&
 \colhead{$\log(M_\mathrm{*,total})$\tablenotemark{e}}\\
\colhead{XMM-Newton name\tablenotemark{a}}
&\colhead{XID\tablenotemark{a}}&\colhead{$z$}&\colhead{($M_\odot$)}&
\colhead{Ref\tablenotemark{c}}&\colhead{(AB mag)}&\colhead{(AB mag)}&
\colhead{$F160W$}&\colhead{(Vega mag)}&\colhead{($M_\odot$)}
}
\startdata
XMMU J100118.5+022739&  14& 1.065&  8.52& 1   &20.70 &19.41&1.9 &0.38      & 11.18\\
XMMU J100046.8+020016&  59& 1.923&  8.72& 2   &--    &20.28&0.6 &0.15--0.75& 11.07--12.12\\
XMMU J095927.7+020010& 219& 1.248&  8.07& 2   &22.55 &20.61&2.1 &0.49      & 11.00\\
XMMU J100035.3+024303& 281& 1.177&  8.25& 1,2 &22.63 &20.25&3.5 &0.74      & 11.44\\
XMMU J095928.5+015934& 329& 1.166&  8.05& 1,2 &22.95 &21.03&1.7 &0.58      & 10.90\\
XMMU J100130.7+021147&2148& 1.526&  8.43& 1   &23.45 &21.15&1.5 &0.40      & 10.88\\
XMMU J100243.8+020502&2261& 1.260&  8.05& 1   &--    &21.38& -- &0.15--0.75& 10.25--11.09\\
XMMU J100226.9+015938&2637& 1.630&  8.35& 1   &24.12 &20.72&1.8 &0.73      & 11.64\\
XMMU J095903.2+022001&5049& 1.131&  8.40& 2   &22.66 &20.82&1.5 &0.57      & 10.93\\
XMMU J095908.1+024310&5230& 1.359&  8.22& 1   &--    &19.13& -- &0.15--0.75& 11.21--12.08\\
\enddata
\tablenotetext{a}{Original XMM-Newton source name and ID \citep{capp09}}
\tablenotetext{b}{Mean value where two measurements are available}
\tablenotetext{c}{Source for $M_\mathrm{BH}$: (1) Magellan/IMACS
  \citep{trum09b}; (2) zCOSMOS/VIMOS \citep{merl09};
  $M_\mathrm{BH}$-errors are quoted as 0.4~dex and 0.3~dex, respectively}
\tablenotetext{d}{From free-$n$ fit before fixing}
\tablenotetext{e}{Total uncertainty is $\pm$0.4~dex
  (Section~\ref{sec:uncertainty}).}

\end{deluxetable*}

\section{Host galaxy masses}

\subsection{Observed host galaxy photometry and colors}
We obtain information on the host galaxies using broad-band photometry
from the high-resolution HST ACS/WFC images in the F814W (=$I$)
filter with 0\farcs03/pixel sampling \citep{koek07}, and the NICMOS/NIC3
parallels in the F160W (=$H$) filter at 0\farcs101/pixel, both integrated
for one orbit.

In the NICMOS $H$-band the host galaxies of all AGN are clearly resolved,
visible already to the unaided eye. We hence extract the host galaxy flux
from the composite galaxy+AGN NICMOS image by modelling the
two-dimensional light distribution of each object using GALFIT
\citep{peng02,peng09} in Version~3.0 (C.~Y.\ Peng, priv.\ comm.). We
restrict our models to a point-source plus a single elliptical S\'ersic
(\citeyear{sers68}) profile. Previous simulations show that at our
resolution and depth it is unreliable to use many and/or complex galaxy
components \citep{simm08,sanc04a}. We carry out several passes of GALFIT,
first with free S\'ersic parameter\footnote{$n=1$ represents an
  exponential disk, $n=4$ a de Vaucouleurs spheroid.}  $n$, and
subsequently fixed $n = 1$, 2 or 4, depending on the best free-$n$-fit.

We require GALFIT to converge on a sensible solution, indicating that an
actual host galaxy is being described and not e.g.\ uncertainties in the
point spread function (PSF). This means we require a GALFIT
nucleus-to-host contrast of $<3.25$~mag, a half-light radius of $r>2$
pixel, and $0.5<n<8$.

XID~2261 and 5230 both show unrealistically compact host galaxy model
scale-lengths, indicating an unsuccessful GALFIT model. For these two
cases we simply subtract the best fitting single point-source model --
without a S\'ersic component -- from the original image, resulting in an
only slightly oversubtracted/underestimated host galaxy. We use aperture
photometry on the host galaxy in these cases, and the GALFIT galaxy model
magnitude for the other eight.

In total we resolve {\em all} 10 host galaxies in the sample -- no object
drops from the sample due to high nucleus-to-host contrast or other
reasons. Extracted host galaxy images are shown in
Figures~\ref{fig:allobj}+\ref{fig:allobj2}.
\medskip

In the ACS $I$-band data the contrast between AGN and host galaxy is less
favorable than in the $H$-band, as expected from the near-UV
SED. Therefore we take a two-step approach. First, we carry out ``peak
subtraction'' removal of the nuclear component for the $I$-band, by
scaling a PSF to the central 4 pixel aperture flux. From statistics on the
expected random residuals after subtraction of this scaled PSF for several
1000 stars \citep[][]{jahn04b}, we require a residual flux of $>$5\% for a
host galaxy to have a high probability of being resolved. This is the case
for 7/10 objects -- the hosts of XIDs~59, 2261, and 5230 remain unresolved
in the ACS image.

The seven resolved objects are again modeled using GALFIT, with S\'ersic
$n$ fixed from the $H$-band fit, thus minimizing S/N-dependent biases in
our extracted colors. As for the $H$-band we check the models for
successful convergence, which is the case for all seven resolved objects,
and again use the host model magnitude. The resulting host galaxy
photometry is also listed in Table~\ref{tab:all}.

\begin{figure*}
\begin{center}
\includegraphics[width=15cm,clip]{sampleobjects1.eps}
\end{center}
\caption{\label{fig:allobj}
Nucleus-removed host galaxies: Galaxy plus nucleus model fitted for 8/10
objects (HST/NIC3 F160W). Images are $7\arcsec\times7\arcsec$, some
objects lie near NICMOS tile edges.
}
\end{figure*}

\begin{figure}
\begin{center}
\includegraphics[width=7.5cm,clip]{sampleobjects2.eps}
\end{center}
\caption{\label{fig:allobj2}
As Figure~\ref{fig:allobj}: XID~2261 and 5230, only nucleus model fitted
and subtracted.
}
\end{figure}

\subsection{Host galaxy stellar masses}
Our general approach is to constrain the mass-to-light ratio (M/L) of each
host galaxy by a single optical color. We have successfully employed this
method before to study stellar populations of low-$z$ quasar host
galaxies \citep{jahn04a} or stellar ages and masses of quasar hosts at
$0.5<z<3$ \citep{sanc04a,jahn04b,kuhl05,schr08}. Here we compute stellar
masses from the rest-frame $V$-band luminosity in combination with the M/L
from the rest-frame $(B-V)$ color:

\begin{equation}
  L_\mathrm{V} = 10^{-0.4(V-4.82)}
\end{equation}

\begin{equation} \label{eqn:mass}
  M_\mathrm{*,total} = 10^{-0.952+1.710(B-V)}\times L_\mathrm{V},
\end{equation}

with $L$ and $M$ in solar units. This calibration is based on
template-fitted masses \citep[Chabrier IMF]{bruz03} and luminosities
derived for galaxies in the COSMOS field \citep{ilbe09b}. We convert
luminosities to Vega zeropoint and apply a $-$0.124~dex mass offset to
transform to the mass scale of the models by Charlot \& Bruzual (2007/2009
in prep.) that include contributions from TP-AGB stars. The linear
relation of eqn.~\ref{eqn:mass} is a fit to galaxies in redshift-
($1.0<z<1.6$) and color-range ($0.38<(B-V)<0.74$) of our 7 host galaxies
resolved in both HST bands. The scatter of the fit corresponds to an RMS
uncertainty in resulting stellar mass of $\pm$0.3~dex.

We convert our measured (F814W--F160W) colors to restframe
$(B-V)_\mathrm{host}$ by applying both $K$-correction and a color term,
individually for each object and its redshift. For this purpose we
identify for each galaxy the single stellar population model again from
Charlot \& Bruzual (Chabrier IMF, solar metallicity) with the closest
(F814W--F160W) at a given $z$. Since the interpolation intervals are rather
small and we are not using the interpolation SED to extract any
further information, this method is quite insensitive to the exact choice
of models, and errors from M/L-calibration dominate.

For the three host galaxies only detected in F160W we have to make
assumptions for the interpolation template to convert observed F160W to
$L_\mathrm{V}$, as well as the color (or M/L) in
equation~\ref{eqn:mass}. On the red/old side we assume a value of
$(B-V)=0.75$ (Vega ZP) corresponding to the red end of the red sequence at
$z\sim1.5$ \citep{krie08}, a value also consistent with the reddest values
of the rest of the sample. As a blue limit we use $(B-V)=0.15$,
corresponding to 3$\sigma$ from the mean value for COSMOS inactive
galaxies at $z\sim1.4$, $\log M_\mathrm{*,total}>10.7$ and age
$t<1$~Gyr. These assumptions should bracket the true M/L, and provide
robust limits on the host galaxies' stellar mass.

\subsection{Stellar mass uncertainties}\label{sec:uncertainty}
Uncertainties in stellar mass have several sources. In addition to the
$(B-V)$ calibration stated above, the strongest contributions come from
(1) GALFIT precision of extracting the host galaxy, (2) potential
influence from different spatial resolutions in the ACS and NIC3 images,
and (3) dependency of bandpass conversions on assumed templates.

(1) \citet{sanc04a} tested how reliably GALFIT can derive host
galaxy photometry in one-orbit ACS data. Given the typical host galaxy
magnitude of our sample in the F814W filter we conclude an uncertainty of
0.15~mag for F814W, and -- due to counteracting effects of lower
spatial resolution but more favorable galaxy-to-nucleus contrast -- also for
F160W.

(2) GALFIT's functionallity depends on the spatial resolution of an image
and the spatial difference between the AGN and galaxy component. If a
galaxy is compact and the PSF not well characterized, flux transfer
between components is possible, its amplitude potentially depending on
spatial resolution.
We test if the different resolutions of the ACS and NIC3 have a
significant influence, by rebinning a mock AGN+galaxy image resembling a
typical object to different spatial resolutions. The recovered host galaxy
photometry shows an rms variation of $\sim$0.15~mag, no systematic offset,
and a negligible trend with resolution. This rms scatter is consistent
with the GALFIT precision from (1) above and we conclude that resolution
effects are insignificant.

(3) The bandpass conversion and $K$-correction depend on the IMF and
metallicity of the assumed single population interpolation model. The
masses derived from Chabrier and Salpeter IMF generally differ by less than
0.1~dex. The metallicity of galaxies of the estimated masses at $1<z<2$ is
expected to range between solar and 3$\times$ solar
\citep[$8.8<12+\log(O/H)<9.2$,][]{trem04,maio08}. Changing from the solar
metallicity used, to the most metal-rich templates ($Z=0.05$), masses
change by $\la$0.05~dex. We conclude that our stellar masses should be
good to within $\sim$0.1~dex from the choice of bandpass conversion SEDs.
\medskip

In combination we find our error budget in stellar mass is 0.21~mag for
the $(B-V)$-color (0.15~mag from each band), corresponding to 0.27~dex
uncertainty in stellar mass -- dominating the three sources of uncertainty
above. Adding the uncertainty of the mass calibration from
equation~\ref{eqn:mass}, the total uncertainty in stellar mass is
$\pm$0.4~dex.

\section{Results}

While we can not estimate bulge masses directly, we find that the seven
objects with direct {\em total stellar mass} estimates fall directly onto
the $M_\mathrm{BH}$--$M_\mathrm{*,bulge}$-relation
(Fig.~\ref{fig:mbh_mstell}) of the local universe from \citet{haer04}. The
three objects with a bracketing range on stellar mass are also consistent
with the local relation. The objects have a maximum deviation of
$\sim$0.3~dex perpendicular to the $z=0$ relation. The seven non-limit AGN
in the $z$-range $1.06<z<1.65$ show a mean ratio
$M_\mathrm{BH}/M_\mathrm{*,total}=0.00178\pm0.0012$ and mean
$\log(M_\mathrm{BH})=8.31$. This is consistent with the value at $z=0$ of
$M_\mathrm{BH}/M_\mathrm{*,bulge}=0.00165$ at the same $M_\mathrm{BH}$,
and has exactly the same 0.3dex scatter.

\citet{merl09} compute stellar masses for a larger sample of COSMOS type-1
AGN, using an independent SED-decomposition method. For 18 galaxies where
stellar masses could be estimated with both methods (five objects are part
of this study, 13 have no BH mass estimates yet) their masses are smaller
by 0.1--0.2~dex. This agreement within the error bars reinforces our
conclusion that our mass estimates are robust. In total they find a mild
deviation from the $z=0$ relation, but most of their signal comes from
objects at $z>1.5$, which is not well covered by our study.

\begin{figure}
\centerline{\includegraphics[height=8cm,clip]{haeringrix_nicmos_lim_rev.eps}}
\caption{\label{fig:mbh_mstell}
$M_\mathrm{BH}$--$M_\mathrm{*,total}$-relation from COSMOS ACS+NIMOS:
  Shown are ten type-1 AGN in the redshift range $1.06<z<1.92$, seven with
  direct stellar mass estimates (red circles) and three with bracketing
  range (blue circles and lines). The points are overplotted over the
  local $M_\mathrm{BH}$--$M_\mathrm{*,bulge}$-relation by \citet{haer04}
  (black triangles) and its best fit (solid line, $\log M_\mathrm{BH} =
  8.2 + 1.12\times(\log M_\mathrm{*,bulge}-11)$). At a given mass there is
  no difference in the $M_\mathrm{BH}$/$M_\mathrm{*,total}$-ratio at
  $z\sim1.4$ and local $M_\mathrm{BH}$/$M_\mathrm{*,bulge}$-ratio, for the
  sampled BH mass range, $\log M_\mathrm{BH}\sim8.3$.
}
\end{figure}

\section{Discussion \& Conclusions}
\subsection{Completeness and systematics}
Are there systematic effects inherent in the dataset that prevent us from
detecting objects which deviate more strongly from this relation? X-ray
selection and subsequent broad-line AGN classification finds all AGN down
to a contrast of AGN/galaxy$\la$10\%, beyond which the stellar light
swamps any AGN signature. The general limit is set by the magnitude limit
of the spectral followup, modified by this contrast. At the observed
$M_\mathrm{BH}$ we expect generally high accretion rates
\citep[e.g.][]{zhen09}, thus this detection limit converts to a limit in
$M_\mathrm{BH}$ present in the data. The additional ACS and NIC3 coverage
are random processes. Being able to extract a host galaxy for every single
object of our sample means that aside from the $M_\mathrm{BH}$-limit we
have no completeness limitations with respect to galaxy mass. In summary,
we do not expect missing objects in the $M_\mathrm{BH}$-regime currently
populated.
\medskip

We also comment on the luminosity function or ``Lauer'' bias
\citep{laue07}: The $M_\mathrm{BH}$--$M_\mathrm{*,total}$-relation has a
scatter and the BH luminosity function drops rapidly towards higher
$L$. In combination, a flux-limited sample of AGN will select many more
massive BHs in lower mass galaxies than vice versa, biasing the sample
towards higher $M_\mathrm{BH}$.

Despite the local comparison sample not being selected by $M_\mathrm{BH}$
-- it consists of inactive galaxies -- we do not expect a large bias: The
\citet{merl09} and \citet{trum09b} BH masses follow the mass calibration
by \citet{onke04}. Their virial $M_\mathrm{BH}$ for active galaxies are
calibrated by a forced match to the $M_\mathrm{BH}$--$\sigma_*$-relation for
inactive galaxies. This might create ``wrong'' $M_\mathrm{BH}$ estimates,
but it compensates most of the expected offset $\Delta M_\mathrm{BH}$
between our AGN-sample and any (active or inactive) sample in the local
universe.  Together with only a small scatter in the local
$M_\mathrm{BH}$--$M_\mathrm{*,bulge}$ it is not surprising to find no
discernable signature of the LF-bias in our data.

\subsection{For Massive Galaxies: Relative Non-Evolution of BH Mass Versus Total Stellar Mass}

\citet{haer04} derive dynamical {\em bulge} masses -- we derive {\em
  total} stellar masses of the host galaxies. Interpreting the coincidence
of these two relations at $z=0$ and $z\sim1.4$ depends on how much stellar
mass in our sample is actually part of a disk component. The stellar
masses of $\log M_\mathrm{*,total}=11.1\pm0.3$ are already high in the
galaxy mass function and our sample galaxies are the likely progenitors of
giant ellipticals in the local universe. If they were already
bulge-dominated at the observed redshifts, nine billion years earlier, our
total masses revert to bulge masses and our observations mean zero
evolution (though consistent with evolution of factor $<$2.65 within
errors) in $M_\mathrm{BH}$--$M_\mathrm{*,bulge}$, at these
$M_\mathrm{BH}$.

\subsection{Evolution of BH Mass Versus Bulge Mass and Velocity Dispersion}
While the actual bulge-to-total mass-ratio $(B/T)_\mathrm{mass}$ of our
galaxies -- in fact of any galaxy at these redshifts -- is unknown, we
have circumstantial evidence that our galaxies could contain substantial
disk components: (a) Visual impression in
Figs.~\ref{fig:allobj}+\ref{fig:allobj2}, (b) some galaxies have S\'ersic
indices near $n=1$ (Table~\ref{tab:all}), (c) some best fitting SSP models
have ages down to 1~Gyr (oldest: 5~Gyr), (d) the $z=0.36$ host galaxies of
\citet{treu07}, with similar $M_\mathrm{BH}$, have
$(B/T)_\mathrm{mass}\sim1/3$, (e) more than 50\% of $\log
M_\mathrm{*,total}=11.1\pm0.3$ galaxies are strongly starforming at
$z\sim1$ \citep{noes07} and (f) of them $\la40\%$ are on the red sequence
\citep[at $z=1.5$,][]{tayl09}.

With a substantial disk component a different interpretation is possible:

\begin{itemize}
\item We know, independently of their exact evolutionary path, these
  galaxies have to end up on the local
  $M_\mathrm{BH}$--$M_\mathrm{*,bulge}$-relation.
 If they currently obey the same relation but with
  $M_\mathrm{*,total}$, then their masses in bulge+disk at $z=1.4$ can end
  up in their bulges 9~Gyrs later with no addition of mass to either the
  BH or the bulge from outside the galaxy or from star-formation.

\item What is required is a process redistributing disk stars to the
  bulge: major and minor galaxy merging \citep{hopk09} and disk
  instabilities \citep{parr09} together dominate bulge-creation at high
  masses.

\item In addition to the disk-to-bulge conversion process these galaxies can
  still grow. Merging of similar systems will just coadd BHs and stellar
  components separately, moving the system parallel along the local
  scaling relation. However total mass growth is limited since the
  observed evolution in space density of massive galaxies at $z<1.5$ is
  small \citep[e.g.,][and references therein]{ilbe09}.

\item Even wet mergers are allowed: The gas conversion efficiency is not
  very high -- in individual mergers \citep{crot06b} and also overall
  \citep[at $z\sim0.6$,][]{roba09} only 10\% of SF arise from merging.

\item Any subsequent (or ongoing) AGN phase will also less than double the
  BH mass, since high accretion states are rare and short
  \citep{hopk06c}. Any BH or stellar mass change below a factor of two
  will lie within the 0.3~dex scatter of the local relation.

\end{itemize}

In the case of a disk component in our galaxies at $z=1.4$, bulge masses are
smaller than the total stellar mass in Fig.~\ref{fig:mbh_mstell} and the
$M_\mathrm{BH}$--$M_\mathrm{*,bulge}$-relation will actually evolve, when
the galaxy-structures are changing over time through merging. They will
constantly move closer to the local relation -- consistent with
predictions from simulations of a merger-driven bulge evolution
\citep{crot06b}.

This has the implication that our result is consistent with the
non-evolution of the {\em bulge} mass relation at $z<1.7$ for massive
early type quasar host galaxies found by \citet{peng06a,peng06b} -- there
$M_\mathrm{*,bulge}=M_\mathrm{*,total}$. At the same time we agree with
the strong evolution in $M_\mathrm{BH}$--$M_\mathrm{*,bulge}$ and
$M_\mathrm{BH}$--$\sigma_\mathrm{*,bulge}$ claimed by \citet{treu07} and
\citet{woo08}: If we assume as a limit for our galaxies
$(B/T)_\mathrm{mass}\le1/3$, our results predict an evolution of
$M_\mathrm{BH}$/$M_\mathrm{*,bulge}\le (1+z)^{1.2}$, consistent with
$M_\mathrm{BH}$/$M_\mathrm{*,bulge}\le (1+z)^{1.5\pm1.0}$ found by
\citet{treu07} for Seyfert 1 galaxies of the same black hole mass at
$z=0.36$.

This path of converging co-evolution of bulge and BH is fully indepedent
of any interaction or feedback {\em between} bulge and BH -- the AGN can
even be switched off since $z=1.4$. What we are witnessing is the formation of
the bulge independent of the BH. Only before $z=1.4$ is a mechanism
required to connect total stellar mass and BH mass in massive galaxies. We
might see an indication of this in the large offset of the relation for
$z>4$ \citep{walt04,riec08a,riec08b,riec09}, but the Lauer-bias and the
different mass scale ($M_\mathrm{BH}>9.2$) complicate the picture. Since
stellar mass and BH mass buildup trace each other very well over cosmic
time with a factor similar to the local ratio of BH and stellar mass
\citep{zhen09}, the coupling mechanism can be very indirect and does not
need to be dominated by a strong version of AGN feedback.

\acknowledgments
KJ would like to thank E.~F.\ Bell, N.\ Neumayer, C.~Y.\ Peng, and
A.~van~der~Wel for very fruitful discussions.
KJ is supported through the Emmy Noether Programme of the German Science
Foundation (DFG).
DR acknowledges support from from NASA through a Hubble Fellowship.

{\it Facilities:} \facility{ESO VLT (VIMOS)}, \facility{HST (ACS, NICMOS)},
\facility{Magellan (IMACS)}, \facility{XMM-Newton}.


text:
\section{Introduction and motivation \label{sec:intro}}
The main recent developments on neutrino mixing\cite{Altarelli:2004za,Mohapatra:2005wg,Mohapatra:2006gs,Strumia:2006db,GonzalezGarcia:2007ib} are related to the relatively large value of $\theta_{13}$, which was measured by recent experiments\cite{Abe:2011sj,Adamson:2011qu,Abe:2012tg,An:2012eh,Ahn:2012nd},  and indications of significant deviation of $\theta_{23}$ from maximal mixing\cite{Altarelli:2012vq,Smirnov:2012ei}. In particular, non-maximal $\theta_{23}$ is strongly indicated by recent data, but the sign of $\theta_{23} - \pi/4$ is not yet determined. Results of global analysis are given in refs.\cite{Tortola:2012te,Fogli:2012ua,GonzalezGarcia:2012sz}.

To explain the observed pattern of neutrino mixing, the quark-lepton complementarity (QLC) has been widely investigated in the literature\cite{Raidal:2004iw,Minakata:2004xt,Frampton:2004vw,Kang:2005as,Xing:2005ur,Datta:2005ci,Antusch:2005ca,Minakata:2005rf,Jarlskog:2005jn,Everett:2005ku,Harada:2005km,Dighe:2006zk,Chauhan:2006im,Hochmuth:2006xn,Schmidt:2006rb,Plentinger:2006nb,Plentinger:2007px,Frampton:2008ep,Xing:2009eg,Zheng:2010kp,Barranco:2010we,Patel:2010hr,Shimizu:2010pg,Ahn:2011yj,Li:2011ag,Zhang:2012zh,Zhang:2012pv}. In particular, much attention has been paid to the class of models, $U_{\mbox{\scriptsize PMNS}} = V_{\mbox{\scriptsize CKM}}^\dagger V_M$, which can be obtained in grand unified theories (GUTs). The correlation matrix $V_M$ is simply defined by the product of the CKM and PMNS mixing matrices. In general $V_M$ is not determined by theories, because there is no relation between the Dirac and the Majorana mass operators~\cite{Harada:2005km}. In this class of models, the observed two large mixing angles $\theta_{12}$ and $\theta_{23}$ indicate that $V_M$ has two large mixing angles because all the mixing angles in $V_{\mbox{\scriptsize CKM}}$ are small. As the simplest possibility, for example, we can take $V_M$ being the bimaximal mixing matrix $V_{bm}$,
\begin{align}
	U_{\mbox{\scriptsize PMNS}} = V_{\mbox{\scriptsize CKM}}^\dagger V_{bm}.
	\label{eq:PMNS-CKM}
\end{align}

In this paper we only consider this minimal model. In ref.\cite{Harada:2005km}, expanding $V_{\mbox{\scriptsize CKM}}$ to ${\cal O}(\lambda^4)$ (where $\lambda \equiv \sin \theta_{\mbox{\scriptsize C}} \approx 0.2253$), it is found that
\begin{align}
	\sin^2 2\theta_{13} &= 2\lambda^2 + {\cal O}(\lambda^4) = 0.102 + {\cal O}(\lambda^4), \label{eq:13mixingLO} \\
	\sin^2 2\theta_{23} &= 1 + {\cal O}(\lambda^4), \label{eq:23mixingLO} \\
	\sin^2 2\theta_{12} &= 1 - 2\lambda^2 + {\cal O}(\lambda^4) = 0.898 + {\cal O}(\lambda^4). \label{eq:12mixingLO}
\end{align}
These correspond to $\theta_{13} \simeq \theta_{\mbox{\scriptsize C}}/\sqrt{2} \simeq 9^\circ$, $\theta_{23} \simeq 45^\circ$ and $\theta_{12} \simeq 36^\circ$, which are consistent with experimental results at leading order approximation. In particular, it is interesting that the predicted value of 1-3 mixing, $\theta_{13} \simeq 9^\circ$,  has been confirmed by recent experiments\cite{Abe:2011sj,Adamson:2011qu,Abe:2012tg,An:2012eh,Ahn:2012nd}.

However, the next-to-leading order corrections, the ${\cal O}(\lambda^4)$ terms in eqs.\eqref{eq:13mixingLO}--\eqref{eq:12mixingLO}, should be calculated, because:
\begin{itemize}
	\item the ${\cal O}(\lambda^4)$ corrections may not be enough small with respect to the recent experimental errors.
		Particularly, the Daya Bay Collaboration reported the precise value at $1\sigma$\cite{An:2012eh},
		\begin{align}
			\sin^2 2\theta_{13} = 0.092 \pm 0.016(\mbox{stat}) \pm 0.005(\mbox{syst})
			\label{eq:Daya-Bay}
		\end{align}
		where the magnitude of systematic errors is the same order of that of $\lambda^4$;
	\item the same value of 1-3 mixing, $\theta_{13} \simeq 9^\circ$,  can be obtained in various models with different schemes: flavor symmetries, texture, ans\"{a}tz etc\cite{Smirnov:2012ei}.
		To distinguish models, we need a more precise prediction;
	\item a deviation of $\theta_{23}$ from maximal mixing is the ${\cal O}(\lambda^4)$ correction in eq.\eqref{eq:23mixingLO}.
		To determine the magnitude of deviation, the ${\cal O}(\lambda^4)$ corrections should be calculated;
	\item in the earlier works, eq.\eqref{eq:PMNS-CKM} was analyzed numerically. Analytical formulae of higher-order terms are usually neglected; therefore, it is very unclear that which parameter is relevant at each order.
\end{itemize}

Motivated by these points, in this paper we perform analytical calculations of the ${\cal O}(\lambda^4)$ corrections. We show that the next-to-leading order corrections
\begin{itemize}
	\item reduce the value of $\sin^2 2\theta_{13}$ from 0.102 to 0.092, which is completely consistent with the result of Daya Bay Collaboration of eq.\eqref{eq:Daya-Bay}.
	\item lead to a deviation of $\theta_{23}$ from maximal mixing. The predicted value is $\sin^2 \theta_{23} = 0.446 + {\cal O}(\lambda^6)$, which corresponds to $\theta_{23} \simeq 41.9^\circ$.
	\item provide the value of $\sin^2 \theta_{12}$ of 0.335, which is very close to the predicted value of tri-bimaximal (TBM) mixing of 0.333.
\end{itemize}

\section{Neutrino mixing angles \label{sec:Mixingangles}}
\subsection{1-3 mixing $\theta_{13}$ \label{subsec:theta_13}}
The recent experimental results on $\theta_{13}$ from two accelerator experiments, T2K\cite{Abe:2011sj} and MINOS\cite{Adamson:2011qu}, and from three reactor experiments, Double-Chooz\cite{Abe:2012tg}, Daya Bay\cite{An:2012eh} and RENO\cite{Ahn:2012nd} were very important developments in neutrino physics.
The global fit value of $\theta_{13}$ is rather large, $\theta_{13} \simeq 9^\circ$ for the normal and inverted mass hierarchy, and $\theta_{13} = 0$ is now excluded at more than 10$\sigma$.

The 1-3 mixing angle $\theta_{13}$ is given by
\begin{align}
	\sin^2 2\theta_{13} = 4 |U_{e3}|^2\left( 1 - |U_{e3}|^2\right). \label{eq:sin22theta_13}
\end{align}
Eq.\eqref{eq:sin22theta_13} indicates that $\theta_{13}$ is determined only by $U_{e3}$.
In eq.\eqref{eq:sin22theta_13}, $\theta_{13}$ is a parameter defined by the standard parametrization, and the $U_{e3}$ is the matrix element of any parametrization of $U_{\mbox{\scriptsize PMNS}}$\cite{Ohlsson:2002rb}.

Using eq.\eqref{eq:Ue3^2} in the appendix, we find that
\begin{align}
	\sin^2 2\theta_{13}
	&= 2\lambda^2 - \left\{ 1 + 4A (1 - \overline{\rho})\right\} \lambda^4 + {\cal O}(\lambda^6), \nonumber \\
	&= \underbrace{0.1016}_{\lambda^2 \ \mbox{\scriptsize term}} - \underbrace{0.0098}_{\lambda^4 \ \mbox{\scriptsize term}} + {\cal O}(\lambda^6), \nonumber \\
	&= 0.092 + {\cal O}(\lambda^6), \label{eq:13mixingNLO}
\end{align}
where we have used the best-fit values of $\lambda$, $A$ and $\overline{\rho}$ of eq.\eqref{eq:Wolfenstein_values}. It was confirmed that varying values in the error range in eq.\eqref{eq:Wolfenstein_values} are almost negligible. Eq.\eqref{eq:13mixingNLO} corresponds to $\theta_{13} \simeq 8.8^\circ$ (an another solution of $\theta_{13} \simeq 81.2^\circ$ is excluded, because $\sin^2 \theta_{13} = |U_{e3}|^2 < 0.5$).

The key points of eq.\eqref{eq:13mixingNLO} are:
\begin{itemize}
	\item  $\sin^2 2\theta_{13}$ is determined only by one parameter $\lambda$ at leading order, and three parameters $\lambda$, $A$ and $\overline{\rho}$ are relevant at next-to-leading order. The  $\overline{\eta}$ is irrelevant up to ${\cal O}(\lambda^6)$;
	\item the next-to-leading order corrections reduce the value of $\sin^2 2\theta_{13}$ from 0.102 to 0.092, which is completely consistent with the result of Daya Bay Collaboration of eq.\eqref{eq:Daya-Bay}. Thus, the value of $\sin^2 2\theta_{13}$ becomes smaller by $9.8 \%$ due to ${\cal O}(\lambda^4)$ correction terms;
	\item the magnitude of the ${\cal O}(\lambda^4)$ correction terms is approximately twice larger than that of the systematic errors of Daya Bay Collaboration, and therefore the ${\cal O}(\lambda^4)$ corrections cannot be negligible.
\end{itemize}

In figure~\ref{fig:13mixing}, we summarize the theoretical and experimental values of $\sin^2 2\theta_{13}$.
The results from T2K\cite{Abe:2011sj} and MINOS\cite{Adamson:2011qu} are shown for $\sin^2 \theta_{23}=0.4$ and $\delta=\pi$. For the QLC predictions, eq.\eqref{eq:13mixingLO} is shown for $\sin^2 2\theta_{13} = 0.102 \pm 5 \lambda^4$ (leading order), and eq.\eqref{eq:13mixingNLO} is shown for $\sin^2 2\theta_{13} = 0.092 \pm 5\lambda^6$ (next-to-leading order), which correspond to $0.089 < \sin^2 2\theta_{13} < 0.114$ and $0.091 < \sin^2 2\theta_{13} < 0.093$, respectively. Figure~\ref{fig:13mixing} shows that the QLC prediction at next-to-leading order is consistent with the results of five experiments\cite{Abe:2011sj,Adamson:2011qu,Abe:2012tg,An:2012eh,Ahn:2012nd} and three global fits\cite{Tortola:2012te,Fogli:2012ua,GonzalezGarcia:2012sz} at the $1\sigma$ level.

\begin{figure}[tbp]
\centering
\hfill
\includegraphics[width=1\textwidth,origin=c]{13mixing.pdf}
\caption{\label{fig:13mixing} Determination of the 1-3 mixing. Shown are the results from T2K\cite{Abe:2011sj}, MINOS\cite{Adamson:2011qu}, Double Chooz\cite{Abe:2012tg}, Daya Bay\cite{An:2012eh}, RENO\cite{Ahn:2012nd} and global fits of Forero et al.\cite{Tortola:2012te}, Fogli et al.\cite{Fogli:2012ua} and Gonzalez-Garcia et al.\cite{GonzalezGarcia:2012sz} for the normal hierarchy case (in the inverse hierarchy case the values do not differ by much). The QLC predictions are shown at leading and next-to-leading order.}
\end{figure}

\subsection{2-3 mixing $\theta_{23}$ \label{subsec:theta_23}}

The recent global analysis indicates that there is a solid deviation of $\theta_{23}$ from maximal mixing\cite{Altarelli:2012vq,Smirnov:2012ei}. At present, the sign of $\theta_{23} - \pi/4$ is not yet determined\cite{Tortola:2012te, Fogli:2012ua, GonzalezGarcia:2012sz}.

Eq.\eqref{eq:23mixingLO} shows that the 2-3 mixing is maximal at leading order.
Eq.\eqref{eq:23mixingLO} also shows that a deviation from maximal mixing is not ${\cal O}(\lambda^2)$ but the ${\cal O}(\lambda^4)$ effect. Therefore, a deviation is not so large.

The 2-3 mixing angle $\theta_{23}$ is given by
\begin{align}
	\sin^2 2 \theta_{23} = \frac{4 |U_{\mu 3}|^2 |U_{\tau 3}|^2}{\left( 1 - |U_{e3}|^2\right)^2}. \label{eq:sin22theta_23}
\end{align}
Eq.\eqref{eq:sin22theta_23} shows that $\theta_{23}$ is determined by three matrix elements, $U_{\mu 3}$, $U_{\tau 3}$ and $U_{e3}$. Of these elements, two are independent because of the unitarity condition $|U_{e3}|^2 + |U_{\mu 3}|^2+|U_{\tau 3}|^2=1$.

From eqs.\eqref{eq:Umu3^2}, \eqref{eq:Utau3^2} and \eqref{eq:1/(1-Ue32)^2}, we find that
\begin{align}
	\sin^2 2\theta_{23} &
	= 1 - 4 \left(\frac{1}{4} +  A \right)^2 \lambda^4  + {\cal O}(\lambda^6), \nonumber \\
	&= 1 - \underbrace{0.0116}_{\lambda^4 \ \mbox{\scriptsize term}} + {\cal O}(\lambda^6), \nonumber \\
	&= 0.988 + {\cal O}(\lambda^6). \label{eq:23mixingNLO}
\end{align}
Eq.\eqref{eq:23mixingNLO} indicates that the value of $\sin^2 2\theta_{23}$ is determined only by two parameters, $\lambda$ and $A$ at next-to-leading order. The $\overline{\rho}$ and $\overline{\eta}$ are irrelevant up to ${\cal O}(\lambda^6)$.

To determine the sign of $\theta_{23} - \pi /4$, we calculate $\sin^2 \theta_{23} = |U_{\mu 3}|^2/(1-|U_{e3}|^2)$ .
From \eqref{eq:1/1-Ue32}, we find that
\begin{align}
	\sin^2 \theta_{23}
	&= \frac{1}{2} - \left(\frac{1}{4} +  A \right)\lambda^2
		- \frac{1}{2}\left(\frac{1}{4} + A\overline{\rho}\right)\lambda^4
		+{\cal O}(\lambda^6), \nonumber \\
	&= \frac{1}{2} - \underbrace{0.0539}_{\lambda^2 \ \mbox{\scriptsize term}}
		     - \underbrace{0.0005}_{\lambda^4 \ \mbox{\scriptsize term}}
		+ {\cal O}(\lambda^6), \nonumber \\
	&= 0.446 + {\cal O}(\lambda^6). \label{eq:23mixingNLO2}
\end{align}
Eq.\eqref{eq:23mixingNLO2} indicates $\theta_{23} < \pi/4$. It is easy to confirm that eqs.\eqref{eq:23mixingNLO} and \eqref{eq:23mixingNLO2} are consistent. The predicted value of $\theta_{23}$ is $\theta_{23} \simeq 41.9^\circ$.

In figure~\ref{fig:23mixing}, we summarize the values of $\sin^2 \theta_{23}$. The results of global analysis of Forero et al.\cite{Tortola:2012te} and Fogli et al.\cite{Fogli:2012ua} are shown for the normal (NH) and inverted (IH) hierarchy cases. For the QLC predictions, eq.\eqref{eq:23mixingLO} is shown for $1 - 5\lambda^4 < \sin^2 2\theta_{23} < 1$ (leading order), and eq.\eqref{eq:23mixingNLO} is shown for $\sin^2 2\theta_{23} = 0.988 \pm 5\lambda^6$ (next-to-leading order), which correspond to $0.443 < \sin^2 \theta_{23} < 0.500$ and $0.444 < \sin^2 \theta_{23} < 0.447$, respectively. Figure~\ref{fig:23mixing} shows that the QLC prediction at next-to-leading order is consistent with the global fits\cite{Tortola:2012te,Fogli:2012ua,GonzalezGarcia:2012sz} for $\theta_{23} < \pi /4$ at the 2$\sigma$ level.

\begin{figure}[tbp]
\centering
\hfill
\includegraphics[width=1\textwidth,origin=c]{23mixing.pdf}
\caption{\label{fig:23mixing} Determination of the 2-3 mixing. Shown are the results of global analysis from Forero et al.\cite{Tortola:2012te}, Fogli et al.\cite{Fogli:2012ua} and Gonzalez-Garcia et al.\cite{GonzalezGarcia:2012sz} for the normal and inverted hierarchy cases. The QLC predictions are shown at leading and next-to-leading order.}
\end{figure}

\subsection{1-2 mixing $\theta_{12}$ \label{subsec:theta_12}}
It has been known that a deviation of $\theta_{12}$ from maximal mixing is large\cite{Altarelli:2004za,Mohapatra:2005wg,Mohapatra:2006gs,Strumia:2006db,GonzalezGarcia:2007ib, Hall:2013yha}.
This can be naturally obtained in the model of eq.\eqref{eq:PMNS-CKM}: eq.\eqref{eq:12mixingLO} shows that a deviation of $\theta_{12}$ from maximal mixing is not ${\cal O}(\lambda^4)$ but the ${\cal O}(\lambda^2)$ effect, and therefore a large deviation is obtained.

The 1-2 mixing angle is given by
\begin{align}
	\sin^2 2\theta_{12} = \frac{4 |U_{e1}|^2 |U_{e2}|^2}{\left( 1 - |U_{e3}|^2 \right)^2}. \label{eq:sin22theta_12}
\end{align}
Using eqs.\eqref{eq:Ue1^2}, \eqref{eq:Ue2^2} and \eqref{eq:1/(1-Ue32)^2}, we find that
\begin{align}
	\sin^2 2\theta_{12}
	&= 1 - 2\lambda^2 - 4A(1 - \overline{\rho})\lambda^4 + {\cal O}(\lambda^6), \nonumber \\
	&= 1 - \underbrace{0.1016}_{\lambda^2 \ \mbox{\scriptsize term}} - \underbrace{0.0073}_{\lambda^4 \ \mbox{\scriptsize term}} + {\cal O}(\lambda^6), \nonumber \\
      &= 0.891 + {\cal O}(\lambda^6). \label{eq:12mixingNLO}
\end{align}
Eq.\eqref{eq:12mixingNLO} shows that the value of $\sin^2 2\theta_{12}$ is determined by three parameters, $\lambda$, $A$ and $\overline{\rho}$ to ${\cal O}(\lambda^6)$. The $\overline{\eta}$ is irrelevant up to ${\cal O}(\lambda^6)$.

Eq.\eqref{eq:12mixingNLO} corresponds to $\theta_{12} \simeq 35.4^\circ$(an another solution of $\theta_{12} \simeq 54.6^\circ$ is excluded, because $\sin^2 \theta_{12}=|U_{e2}|^2/(1-|U_{e3}|^2) < 0.5$. The $\sin^2 \theta_{12}$ written in terms of $\lambda$, $A$ and $\overline{\rho}$ is given in the appendix). It is interesting that the value of eq.\eqref{eq:12mixingNLO} is very close to the tri-bimaximal (TBM) value of $\sin^2 2\theta_{12}=0.889$\cite{Harrison:2002er,Harrison:2002kp}.

In figure~\ref{fig:12mixing}, we summarize the values of $\sin^2 \theta_{12}$. In addition to the results of global analysis\cite{Tortola:2012te,Fogli:2012ua,GonzalezGarcia:2012sz}, the predicted TBM value of 1/3\cite{Harrison:2002er,Harrison:2002kp}, two proposed golden ratio (GR) values of 0.276(GR1\cite{Datta:2003qg,Kajiyama:2007gx,Everett:2008et,Ding:2011cm}) and 0.345(GR2\cite{Rodejohann:2008ir,Adulpravitchai:2009bg,Ding:2011cm}) are shown. For the QLC predictions, eq.\eqref{eq:12mixingLO} is shown for $\sin^2 2\theta_{12} = 0.898 \pm 5\lambda^4$ (leading order), and eq.\eqref{eq:12mixingNLO} is shown for $\sin^2 2\theta_{12} = 0.891 \pm 5 \lambda^6$ (next-to-leading order), which correspond to $0.331 < \sin^2 \theta_{12} < 0.351$ and $0.334 < \sin^2 \theta_{12} < 0.336$, respectively. Figure~\ref{fig:12mixing} shows that the QLC prediction at next-to-leading order is consistent with the global fits at the 2$\sigma$ level.

\begin{figure}[tbp]
\centering
\hfill
\includegraphics[width=1\textwidth,origin=c]{12mixing.pdf}
\caption{\label{fig:12mixing} Determination of the 1-2 mixing. Shown are the results of global analysis\cite{Tortola:2012te,Fogli:2012ua,GonzalezGarcia:2012sz} and the predicted values of TBM\cite{Harrison:2002er,Harrison:2002kp}, GR1\cite{Datta:2003qg,Kajiyama:2007gx,Everett:2008et,Ding:2011cm} and GR2\cite{Rodejohann:2008ir,Adulpravitchai:2009bg,Ding:2011cm}. The QLC predictions are shown at leading and next-to-leading order.}
\end{figure}

\section{Summary and conclusions \label{sec:conclusions}}
In this paper we have presented analytical formulae of neutrino mixing angles at the next-to-leading order in the framework of eq.\eqref{eq:PMNS-CKM}. It has been shown that higher order corrections were important to explain the observed pattern of neutrino mixing. Some conclusions obtained in this paper are given below.
\begin{itemize}
	\item $\theta_{13}:$ $\sin^2 2\theta_{13}$ is determined by three parameters $\lambda$, $A$ and $\overline{\rho}$ up to ${\cal O}(\lambda^6)$. The $\overline{\eta}$ is irrelevant up to ${\cal O}(\lambda^6)$. The ${\cal O}(\lambda^4)$ corrections reduce the value of $\sin^2 2\theta_{13}$ from 0.102 to 0.092, which is consistent with five experiments and three global fits at the $1\sigma$ level. A summary is shown in figure~\ref{fig:13mixing}.
	\item $\theta_{23}:$ the ${\cal O}(\lambda^4)$ corrections lead to a deviation of $\theta_{23}$ from maximal mixing. A deviation of $\sin^2 2\theta_{23}$ from 1 is determined only by two parameters $\lambda$ and $A$ up to ${\cal O}(\lambda^6)$. The negative sign of $\theta_{23} - \pi/4$ is predicted, and the obtained value of $\sin^2 \theta_{12}$ is 0.446, which is consistent with three global fits at the $2\sigma$ level. A summary is shown in figure~\ref{fig:23mixing}.
	\item $\theta_{12}:$ $\sin^2 2\theta_{12}$ is determined by three parameters $\lambda$, $A$ and $\overline{\rho}$ up to ${\cal O}(\lambda^6)$. The $\overline{\eta}$ is irrelevant up to ${\cal O}(\lambda^6)$. The predicted value of $\sin^2 \theta_{12}$ at next-to-leading order is 0.335, which is very close to the TBM value of 0.333. These values are consistent with the global fits at the $2\sigma$ level. A summary is shown in figure~\ref{fig:12mixing}.
\end{itemize}

In this paper we have shown that eq.\eqref{eq:PMNS-CKM} is consistent with experimental results with high precision. If it will be confirmed by  experiments that some other corrections from eq.\eqref{eq:PMNS-CKM} are very small or negligible, it may indicate that there exists an unknown theoretical mechanism behind eq.\eqref{eq:PMNS-CKM}. Therefore, it is very interesting to test eq.\eqref{eq:PMNS-CKM} by near future experiments.

\section{Appendix: The matrix elements of the PMNS matrix \label{app:PMNSelements}}
In this appendix we present the PMNS matrix elements $U_{\alpha i}$$(\alpha = e,\mu,\tau, i=1,2,3)$. We use the Wolfenstein parametrization\cite{Wolfenstein:1983yz,Beringer:1900zz} of $V_{\mbox{\scriptsize CKM}}$.

The $V_{\mbox{\scriptsize CKM}}$ to ${\cal O}(\lambda^4)$ has been widely used in the literature, however, $V_{\mbox{\scriptsize CKM}}$ to ${\cal O}(\lambda^6)$  is necessary for calculations at next-to-leading order.
We can write $V_{\mbox{\scriptsize CKM}}$ to ${\cal O}(\lambda^6)$ in terms of the Wolfenstein parameters $\lambda$, $A$, $\overline{\rho}$ and $\overline{\eta}$,
\begin{align}
	 V_{\mbox{\scriptsize CKM}} &=
	\begin{pmatrix}
		V_{ud} & V_{us} & V_{ub} \\
		V_{cd} & V_{cs} & V_{cb} \\
		V_{td} & V_{ts} & V_{tb}
	\end{pmatrix},
\end{align}
\begin{align}
	V_{ud} & = 1-\frac{\lambda^2}{2} - \frac{\lambda^4}{8} + {\cal O}(\lambda^6), \\
	V_{us} & = \lambda +{\cal O}(\lambda^6), \\
	V_{ub} & = A\lambda^3\left(1+\frac{\lambda^2}{2}\right)(\overline{\rho} - i \overline{\eta}) +{\cal O}(\lambda^6), \\
	V_{cd} & = -\lambda +A^2 \lambda^5 \left( \frac{1}{2} - \overline{\rho} - i \overline{\eta} \right) +{\cal O}(\lambda^6), \\
	V_{cs} & = 1-\frac{\lambda^2}{2} - \frac{\lambda^4}{8}\left( 1 + 4A^2 \right) +{\cal O}(\lambda^6), \\
	V_{cb} & = A\lambda^2 +{\cal O}(\lambda^6), \\
	V_{td} & = A\lambda^3(1-\overline{\rho} - i \overline{\eta}) +{\cal O}(\lambda^6), \\
	V_{ts} & = -A\lambda^2 + A\lambda^4(\frac{1}{2}-\overline{\rho} - i\overline{\eta}) +{\cal O}(\lambda^6), \\
	V_{tb} & = 1 - \frac{A^2\lambda^4}{2} +{\cal O}(\lambda^6),
\end{align}
where the Wolfenstein parameters $\lambda$, $A$, $\overline{\rho}$ and $\overline{\eta}$ are defined by
\begin{align}
	s_{12} = \lambda, \quad s_{23} = A\lambda^2, \quad
	s_{13}e^{i\delta} = \frac{A\lambda^3(\overline{\rho} + i\overline{\eta})\sqrt{1-A^2\lambda^4}}{\sqrt{1-\lambda^2}[1-A^2\lambda^4(\overline{\rho}+i\overline{\eta})]}.
\end{align}
The quark mixing angles $s_{ij} = \sin \theta_{ij}$ and the Kobayashi-Maskawa CP phase $\delta$ are defined by the standard parametrization.
The CKM matrix written in terms of $\lambda$, $A$, $\overline{\rho}$ and $\overline{\eta}$ is unitary to all orders in $\lambda$.
The values of the Wolfenstein parameters are given by\cite{Beringer:1900zz}
\begin{align}
	\lambda = 0.22535 \pm 0 .00065, \quad
	A = 0.811^{+0.022}_{-0.012}, \quad
	\overline{\rho} = 0.131^{+0.026}_{-0.013}, \quad
	\overline{\eta} = 0.345^{+0.013}_{-0.014}. \label{eq:Wolfenstein_values}
\end{align}

The PMNS mixing matrix in the model of eq.\eqref{eq:PMNS-CKM} is
\begin{align}
	 U_{\mbox{\scriptsize PMNS}} =
	\begin{pmatrix}
		U_{e1} & U_{e2} & U_{e3} \\
		U_{\mu 1} & U_{\mu 2} & U_{\mu 3} \\
		U_{\tau 1} & U_{\tau 2} & U_{\tau 3}
	\end{pmatrix}
	=
	V_{\mbox{\scriptsize CKM}}^\dagger V_{bm}, \quad
	\label{eq:PMNS-bimaximal}
	V_{bm} =
	\begin{pmatrix}
		\frac{1}{\sqrt{2}} & \frac{1}{\sqrt{2}} & 0 \\
		-\frac{1}{2} & \frac{1}{2} & \frac{1}{\sqrt{2}} \\
		\frac{1}{2} & - \frac{1}{2} & \frac{1}{\sqrt{2}}
	\end{pmatrix}.
\end{align}
We summarize the PMNS matrix elements $U_{\alpha i} (\alpha = e, \mu, \tau, i = 1, 2, 3)$ and the squared $|U_{\alpha i}|^2$ to ${\cal O}(\lambda^6)$:
\begin{align}
	U_{e1} &=
			\frac{1}{\sqrt{2}}
			+ \frac{\lambda}{2}
			- \frac{\lambda^2}{2\sqrt{2}}
			+\frac{A\lambda^3}{2}(1 - \overline{\rho} + i \overline{\eta})
		       - \frac{\lambda^4}{8\sqrt{2}}
			-\frac{A^2 \lambda^5}{2} \left(\frac{1}{2} - \overline{\rho} + i \overline{\eta} \right)
			+ {\cal O}(\lambda^6),
			\label{eq:Ue1}\\
	U_{e2} &=
			\frac{1}{\sqrt{2}}
			- \frac{\lambda}{2} - \frac{\lambda^2}{2\sqrt{2}}
			-\frac{A\lambda^3}{2}(1 - \overline{\rho} + i \overline{\eta})
		       - \frac{\lambda^4}{8\sqrt{2}}
			+\frac{A^2 \lambda^5}{2} \left(\frac{1}{2} - \overline{\rho} + i \overline{\eta} \right)
			+ {\cal O}(\lambda^6),
			\label{eq:Ue2}\\
	U_{e3} &=
			- \frac{\lambda}{\sqrt{2}}
			+ \frac{A\lambda^3}{\sqrt{2}}(1-\overline{\rho}+i\overline{\eta})
		       +\frac{A^2 \lambda^5}{\sqrt{2}} \left(\frac{1}{2} - \overline{\rho} + i \overline{\eta} \right)
			+ {\cal O}(\lambda^6),
			\label{eq:Ue3} \\
	U_{\mu 1} &=
			-\frac{1}{2}
			+ \frac{\lambda}{\sqrt{2}}
			+ \frac{\lambda^2}{2}\left(\frac{1}{2} -  A \right)
		       +\frac{\lambda^4}{2}\left\{\frac{1}{8} + \frac{A^2}{2} + A\left(\frac{1}{2} - \overline{\rho} + i \overline{\eta}\right) \right\}
			+ {\cal O}(\lambda^6),
			\label{eq:Umu1}\\
	U_{\mu 2} &=
			\frac{1}{2}
			+ \frac{\lambda}{\sqrt{2}}
			- \frac{\lambda^2}{2}\left(\frac{1}{2} -  A \right)
		 	-\frac{\lambda^4}{2}\left\{\frac{1}{8} + \frac{A^2}{2} + A\left(\frac{1}{2} - \overline{\rho} + i \overline{\eta}\right) \right\}
			+ {\cal O}(\lambda^6),
			\label{eq:Umu2}\\
	U_{\mu 3} &=
			\frac{1}{\sqrt{2}}
			- \frac{\lambda^2}{\sqrt{2}}\left( \frac{1}{2} + A \right)
		       -\frac{\lambda^4}{\sqrt{2}}\left\{\frac{1}{8} + \frac{A^2}{2} - A\left( \frac{1}{2} - \overline{\rho} + i \overline{\eta} \right)\right\}
			+ {\cal O}(\lambda^6),
			\label{eq:Umu3}\\
	U_{\tau 1} &=
			\frac{1}{2}
			- \frac{A\lambda^2}{2}
			+ \frac{A\lambda^3}{\sqrt{2}} (\overline{\rho} + i\overline{\eta})
			- \frac{A^2 \lambda^4}{4}
		       +\frac{A\lambda^5}{2\sqrt{2}} (\overline{\rho} + i\overline{\eta} )
			+ {\cal O}(\lambda^6),
			\label{eq:Utau1}\\
	U_{\tau 2} &=
			-\frac{1}{2}
			+ \frac{A\lambda^2}{2}
			+ \frac{A\lambda^3}{\sqrt{2}} (\overline{\rho} + i\overline{\eta})
			+ \frac{A^2 \lambda^4}{4}
		       +\frac{A\lambda^5}{2\sqrt{2}} (\overline{\rho} + i\overline{\eta} )
			+ {\cal O}(\lambda^6),
			\label{eq:Utau2}\\
	U_{\tau 3} &=
			\frac{1}{\sqrt{2}}
			+ \frac{A\lambda^2}{\sqrt{2}}
			- \frac{A^2 \lambda^4}{2\sqrt{2}}
			+ {\cal O}(\lambda^6);
			\label{eq:Utau3}
\end{align}
\begin{align}
	|U_{e1}|^2 &=
			\frac{1}{2}
			+ \frac{\lambda}{\sqrt{2}}
			- \frac{\lambda^2}{4}
			- \frac{\lambda^3}{\sqrt{2}}\left\{\frac{1}{2} - A(1-\overline{\rho})\right\}
			+\frac{A\lambda^4}{2}(1 - \overline{\rho}) \nonumber \\
			&
			- \frac{\lambda^5}{8\sqrt{2}}(1+2A)\left\{1 + 2A(1 - 2\overline{\rho})\right\}
			+ {\cal O}(\lambda^6),
			\label{eq:Ue1^2}\\
	|U_{e2}|^2 &=
			\frac{1}{2}
			- \frac{\lambda}{\sqrt{2}}
			- \frac{\lambda^2}{4}
			+ \frac{\lambda^3}{\sqrt{2}}\left\{\frac{1}{2} - A(1-\overline{\rho})\right\}
			+\frac{A\lambda^4}{2}(1 - \overline{\rho}) \nonumber \\
			&
			+ \frac{\lambda^5}{8\sqrt{2}}(1+2A)\left\{1 + 2A(1 - 2\overline{\rho})\right\}
			+ {\cal O}(\lambda^6)
			\label{eq:Ue2^2}, \\
	|U_{e3}|^2 &= \frac{\lambda^2}{2} - A\lambda^4(1-\overline{\rho}) + {\cal O}(\lambda^6),	\label{eq:Ue3^2}
\end{align}
\begin{align}
	|U_{\mu 1}|^2 &=
			\frac{1}{4}
			- \frac{\lambda}{\sqrt{2}}
			+ \frac{\lambda^2}{2}\left(\frac{1}{2} + A \right)
			+\frac{\lambda^3}{\sqrt{2}}\left(\frac{1}{2} - A \right)
			- \frac{A\lambda^4}{2}\left(1 - \overline{\rho} \right) \nonumber \\
			&
			+\frac{\lambda^5}{\sqrt{2}}\left\{\frac{1}{8} + \frac{A^2}{2} + A\left(\frac{1}{2} - \overline{\rho}\right)\right\}
			+ {\cal O}(\lambda^6), \label{eq:Umu1^2}\\
	|U_{\mu 2}|^2 &=
			\frac{1}{4}
			+ \frac{\lambda}{\sqrt{2}}
			+ \frac{\lambda^2}{2}\left(\frac{1}{2} + A \right)
			-\frac{\lambda^3}{\sqrt{2}}\left(\frac{1}{2} - A \right)
			- \frac{A\lambda^4}{2}\left(1 - \overline{\rho} \right) \nonumber \\
			&
			-\frac{\lambda^5}{\sqrt{2}}\left\{\frac{1}{8} + \frac{A^2}{2} + A\left(\frac{1}{2} - \overline{\rho}\right)\right\}
			+ {\cal O}(\lambda^6), \label{eq:Umu2^2}\\
	|U_{\mu 3}|^2 &=
			\frac{1}{2}
			- \lambda^2\left( \frac{1}{2} + A \right)
			+ A\lambda^4\left( 1 - \overline{\rho} \right)
			+ {\cal O}(\lambda^6),
			\label{eq:Umu3^2}
\end{align}
\begin{align}
	|U_{\tau 1}|^2 &=
			\frac{1}{4}
			- \frac{A\lambda^2}{2}
			+ \frac{A\lambda^3 \overline{\rho}}{\sqrt{2}}
			+ \frac{A\lambda^5\overline{\rho}}{\sqrt{2}} \left( \frac{1}{2} - A \right)
			+ {\cal O}(\lambda^6), \label{eq:Utau1^2}\\
	|U_{\tau 2}|^2 &=
			\frac{1}{4}
			- \frac{A\lambda^2}{2}
			- \frac{A\lambda^3 \overline{\rho}}{\sqrt{2}}
			- \frac{A\lambda^5\overline{\rho}}{\sqrt{2}} \left( \frac{1}{2} - A \right)
			+ {\cal O}(\lambda^6), \label{eq:Utau2^2}\\
	|U_{\tau 3}|^2 &=
			\frac{1}{2}
			+ A\lambda^2
			+ {\cal O}(\lambda^6). \label{eq:Utau3^2}
\end{align}
For convenience, we present $1/(1-|U_{e3}|^2)$, $1/(1-|U_{e3}|^2)^2$ and $\sin^2 \theta_{12}$ in terms of $\lambda$, $A$ and $\overline{\rho}$,
\begin{align}
	\frac{1}{1-|U_{e3}|^2} &=
			1
			+ \frac{\lambda^2}{2}
			+ \left\{ \frac{1}{4} - A (1 - \overline{\rho}) \right\}\lambda^4
			+ {\cal O}(\lambda^6), \label{eq:1/1-Ue32} \\
	\frac{1}{\left(1-|U_{e3}|^2\right)^2} &=
			1
			+ \lambda^2
			+ \left\{ \frac{3}{4} - 2A (1 - \overline{\rho}) \right\}\lambda^4
			+ {\cal O}(\lambda^6), \label{eq:1/(1-Ue32)^2} \\
	\sin^2 \theta_{12} &= \frac{1}{2} - \frac{\lambda}{\sqrt{2}} - \frac{A}{\sqrt{2}}(1 - \overline{\rho})\lambda^3
	+ \frac{\lambda^5}{\sqrt{2}}\left\{\frac{1}{8} + A(1-\overline{\rho}) + A^2 \left(\frac{1}{2} - \overline{\rho}\right)\right\}
	+ {\cal O}(\lambda^6).
\end{align}

redpajama_set_name: RedPajamaArXiv