\section{Introduction}
The discovery of neutrino oscillations and at least two non-vanishing neutrino
masses, provides clear evidence for Physics Beyond the Standard Model (SM).
The simplest extension of the SM accommodating two non-vanishing neutrino masses
involves the addition of at least two right-handed neutrinos. The most general
gauge invariant Lagrangian includes a right-handed bare Majorana mass matrix $M$.
As a result, the scale of $M$ can be much larger than the
electroweak scale, which leads to an elegant explanation for the smallness of
neutrino masses, through the seesaw mechanism \cite{Minkowski:1977sc},
\cite{Yanagida:1979as}, \cite{Glashow}, \cite{GellMann:1980vs}, \cite{Mohapatra:1979ia}.
The seesaw mechanism necessarily implies violations from $3 \times 3$ unitarity of the
Pontecorvo-Maki-Nakagawa-Sakata  (PMNS) matrix, as well as Z-mediated lepton
flavour violating couplings. The introduction of these heavy right-handed neutrinos
can also have profound cosmological implications since they are a crucial component
of the Leptogenesis mechanism to create the observed Baryon Asymmetry of the Universe
(BAU) \cite{Fukugita:1986hr}. Leptogenesis is a very appealing scenario \cite{reviews1} but for
heavy neutrinos with masses many orders of magnitude higher than the electroweak scale
it is difficult or impossible to test it at low energies
\cite{delAguila:2008hw}, \cite{Deppisch:2015qwa}.
Furthermore, without a flavour model, CP violation at high energies, relevant for
Leptogenesis, cannot be related to CP violation at low energies \cite{Branco:2001pq},
\cite{Rebelo:2002wj}. This may be possible in the context of a flavour model, involving
symmetries which allow to establish a connection between high and low energies
\cite{Branco:2002kt},\cite{Frampton:2002qc}, \cite{Branco:2002xf}, \cite{Branco:2005jr},
\cite{Pascoli:2006ci}, \cite{Branco:2007nb}, \cite{Hagedorn:2016lva},
\cite{Fukugita:2016hzf}. \\

At this stage, it should be emphasised that within the seesaw type I framework,
the observed pattern of neutrino masses and mixing does not require that all the
heavy neutrino masses be much larger than the electroweak scale. In this paper,
we carefully examine the question of whether it is possible, within the seesaw
type I mechanism, to have  experimentally detectable violations of $3 \times 3$
unitarity, taking into account the present experimental constraints. In particular,
we address the following questions:

i) In the seesaw type I mechanism, is it possible to have significant deviations from
$3\times 3$ unitarity of the leptonic mixing matrix? By significant,
we mean deviations which are sufficiently small
to conform to all present stringent experimental constraints on these deviations,
but are sufficiently large to be detectable in the next round of experiments.
These experimental constraints arise from bounds on rare processes.

ii) In the case the scenario described in (i) can indeed be realised within the framework
of seesaw type I, what are the requirements on the pattern of heavy neutrino masses?

For definiteness, we will work in a framework where three right-handed neutrinos are added
to the spectrum of the SM.  Our analysis starts with the introduction of the unitary
$6 \times 6$ mixing matrix ${\cal V}$, characterising all the leptonic mixing. We write this
$6 \times 6$ mixing matrix in terms of four blocks of $3 \times 3$ matrices. Using unitarity
of ${\cal V}$, we show that the full matrix ${\cal V}$ can be expressed in terms of only
three blocks of $3 \times 3$ matrices.
Then we apply these results to the diagonalisation of the $6 \times 6$ neutrino mass
matrix, including both Dirac and Majorana mass terms.  Through the use of a
specially convenient exact parametrisation of the $6 \times 6$
leptonic unitary mixing matrix,
we evaluate deviations from $3 \times 3$ unitarity and derive the
maximum value of the lightest heavy neutrino mass which is
required in order to generate significant deviations from unitarity in the
framework of the seesaw type I mechanism.

The paper is organised as follows. In the next section, we review the seesaw mechanism,
define our notation
and introduce a specially convenient exact parametrisation of the $6 \times 6$
leptonic mixing matrix ${\cal V}$. We evaluate the size of the deviations of
$3 \times 3$ unitarity in the present framework and derive a constraint on the
magnitude of the mass of the heavy Majorana neutrinos, in order to have
significant deviations of unitarity. Numerical examples are given
in section 3 and some of the derivations of our results are put in an Appendix.
Finally we present our conclusions in the last section.

\section{Deviations from Unitarity in the Leptonic Sector}

\subsection{ Type I Seesaw mechanism}
In the context of the Type I seesaw mechanism, with only three right-handed neutrinos
added to the Lagrangian of the SM, the leptonic mass terms are given by:
\begin{eqnarray}
{\cal L}_m  &=& -[\overline{{\nu}_{L}^0} m \nu_{R}^0 +
\frac{1}{2} \nu_{R}^{0T} C \ M \nu_{R}^0+
\overline{l_L^0} m_l l_R^0] + h. c.  \nonumber \\
&=& - [\frac{1}{2}† n_{L}^{T} C \ {\cal M}^* n_L +
\overline{l_L^0} m_l l_R^0 ] + h. c.,
\label{lep}
\end{eqnarray}
There is no loss of generality in choosing a weak basis where $m_l$ is already
real and diagonal. The analysis that follows is performed in this basis.
The neutrino mass matrix $\cal M $ is a $6 \times 6$ matrix and has the form:
\begin{eqnarray}
{\cal M}= \left(\begin{array}{cc}
 0  & m \\
m^T & M \end{array}\right)
\label{calm}
\end{eqnarray}
This matrix is diagonalised by the unitary transformation
\begin{eqnarray}
{\cal V}^T {\cal M}^* {\cal V} = {\cal D} \label{dgm}  \qquad
\mbox{i.e.} \qquad {\cal V}^\dagger {\cal M} = {\cal D} {\cal V}^T
\label{vmvd}
\end{eqnarray}
where
\begin{equation}
{\cal D}=\left(\begin{array}{cc}
d & 0 \\
0 & D \end{array}\right)
\end{equation}
with $d = \mbox{diag}.(m_1, m_2, m_3)$ and $D = \mbox{diag}.(M_1, M_2, M_3)$ denoting
respectively, the light and the heavy Majorana neutrino masses.
The unitary $6 \times 6$ matrix ${\cal V} $ is often denoted in the literature as:
\begin{equation}
{\cal V} =\left(
\begin{array}{cc}
K & R \\
S & Z
\end{array}
\right)\
\label{unit}
\end{equation}
where $K$, $R$, $S$ and $Z$ are $3\times 3$ matrices.
For $K$ and $Z$ non singular, we may write
\begin{equation}
{\cal V} =\left(
\begin{array}{cc}
K & 0 \\
0 & Z
\end{array}
\right) \left(
\begin{array}{cc}
{1\>\!\!\!\mathrm{I}}\, & Y \\
-X & {1\>\!\!\!\mathrm{I}}\,
\end{array}
\right);\quad -X=Z^{-1}S;\quad Y=K^{-1}R  \label{u1}
\end{equation}
From the unitary relation ${\cal V} \ {\cal V} ^{\dagger}={1\>\!\!\!
\mathrm{I}}_{(6\times 6)}$, we promptly conclude that
\begin{equation}
Y=X^{\dagger }  \label{xy}
\end{equation}
The matrix ${\cal V} $ can thus be written:
\begin{eqnarray}
{\cal V} =\left(\begin{array}{cc}
K & KX^\dagger \\
-ZX & Z \end{array}\right) \label{eq15}
\end{eqnarray}
We have thus made clear that the unitary $6 \times 6$ matrix ${\cal V}$
can be expressed in terms of three
independent $3 \times 3$ matrices. From the unitarity
of ${\cal V} $, we obtain:
\begin{equation}
\begin{array}{c}
K\ \left( {1\>\!\!\!\mathrm{I}}\,+X^{\dagger }X\right) \ K^{\dagger }
={1\>\!\!\!\mathrm{I}} \\
\\
Z\ \left( {1\>\!\!\!\mathrm{I}}\,+X\ X^{\dagger }\right) \ Z^{\dagger }
={1\>\!\!\!\mathrm{I}}
\end{array}
\label{kt}
\end{equation}
showing that the matrix $X$ parametrizes the deviations from unitary
of the matrices $K$ and $Z$. More explicitly:
\begin{equation}
\begin{array}{c}
K\ K^{\dagger} = {1\>\!\!\!\mathrm{I}} - K\ X^{\dagger }X \ K^{\dagger}
\\
Z\ Z^{\dagger} = {1\>\!\!\!\mathrm{I}} - Z\ X X^{\dagger } \ Z^{\dagger}
\end{array}
\label{3lin}
\end{equation}

From Eq.~(\ref{dgm}) we derive:
\begin{eqnarray}
-X^\dagger Z^\dagger m^T = d K^T \label{13a} \\
K^\dagger m -X^\dagger Z^\dagger M = - d X^T Z^T  \label{13b} \\
Z^\dagger m^T = D X^* K^T  \label{13c} \\
XK^\dagger m + Z^\dagger M = D Z^T \label{13d}
\end{eqnarray}
replacing $Z^\dagger m^T$ from Eq.~(\ref{13c}) into Eq.~(\ref{13a}) we get
\begin{eqnarray}
d = - X^T \ D \ X
\label{dxdx}
\end{eqnarray}
which implies that:
\begin{eqnarray}
X = \pm i \sqrt{D^{-1}} O_c \sqrt{d}
\label{xc}
\end{eqnarray}
where $O_c$ is a complex orthogonal matrix, i.e., $O_c^T O_c = 1\>\!\!\!\mathrm{I}$,
or explicitly:
\begin{eqnarray}
|X_{ij}| = \left| (O_c)_{ij} \sqrt{\frac{m_j}{M_i}} \right|
\end{eqnarray}
It should be stressed that the parametrisation of the $6 \times 6$ unitary matrix $\cal{V}$
given by  Eq.~(\ref{eq15}) has the especial property of allowing to connect in a
straightforward and simple way the masses of the light and the heavy neutrinos through an
orthogonal complex matrix $O_c$, as can be seen from Eqs.~(\ref{dxdx}) and (\ref{xc}).
This is an important new result which plays a crucial r\^ ole in our analysis.

Since  $O_c$ is an orthogonal complex matrix, not all of its elements need to be small;
furthermore, not all the $M_i$ need to be much larger than the electroweak scale, in order
for the seesaw mechanism to lead to naturally suppressed neutrino masses.
These observations about the size of the elements of X are specially relevant
in view of the fact that some of the important physical implications of the seesaw model
depend crucially on $X$. In particular, the deviations from $3 \times 3$ unitarity are controlled
by $X$, as shown in Eq.~(\ref{3lin}).

Given the importance of the matrix $X$, one may ask whether it is possible to write the
$6 \times 6$ unitary matrix ${\cal V}$ in terms of $3 \times 3$ blocks, where only
$3 \times 3$ unitary matrices enter, together with the matrix $X$. In the Appendix,
we show that this is indeed possible, and that the matrix ${\cal V}$ can be written:
\begin{equation}
{\cal V}=\left(
\begin{array}{cc}
K & R \\
S & Z
\end{array}
\right) =\left(
\begin{array}{cc}
\Omega \left( \sqrt{{1\>\!\!\!\mathrm{I}}\,+X^{\dagger }X}\ \right)^{-1}
& \Omega \left( \sqrt{{1\>\!\!\!\mathrm{I}}
\,+X^{\dagger }X}\ \right)^{-1}X^{\dagger } \\
- \Sigma \left( \sqrt{{1\>\!\!\!\mathrm{I}}\,+X\ X^{\dagger }}\
\right)^{-1}X\  & \Sigma  \left( \sqrt{{1\>\!\!\!\mathrm{I}}
\,+X\ X^{\dagger }}\ \right)^{-1}
\end{array}
\right)  \label{v1}
\end{equation}
where $\Omega$ and $\Sigma$ are $3 \times 3$ unitary matrices given by:
\begin{equation}
\Omega = U_{K}U^{\dagger } \qquad \Sigma = W_{Z}\ W^{\dagger }
\end{equation}
and $U$, $W$ are the unitary matrices that diagonalise respectively
$X^\dagger  X$ and $X  X^\dagger$:
\begin{equation}
U^\dagger \ X^{\dagger }X \ U = d_{X}^{2}; \qquad
W^\dagger X X^\dagger \ W = d_{X}^{2}
\label{uwxx}
\end{equation}
It is also shown in the Appendix, that $U_K$ and $W_Z$ defined by:
\begin{equation}
U_{K} \equiv K\ U\ \sqrt{\left( {1\>\!\!\!\mathrm{I}}
\,+d_{X}^{2}\right) } \qquad
W_Z \equiv Z\ W\ \sqrt{\left( {1\>\!\!\!\mathrm{I}}
\,+d_{X}^{2}\right) }
\label{esta}
\end{equation}
are in fact unitary matrices.

As will be explained in the next section, this will allow us, in our analysis,
to trade the matrix $K$ by the combination $U_{K}U^{\dagger}$
which we identify as the best fit for $U_{PMNS}$ derived under the assumption of
unitarity, multiplied by the remaining factor that parametrises the deviations
from unitarity.

\subsection{On the size of deviations from unitarity}
In the framework of the type I seesaw, it is the block $K$ of
the matrix ${\cal V}$
that takes the r\^ ole  played by
$U_{PMNS}$ matrix at low energies in models with only Dirac-type neutrino masses.
Clearly, in this framework, $K$ is no longer a unitary
matrix. However, present neutrino experiments are putting stringent constraints
on the deviations from unitarity.
In our search for significant deviations from unitarity of $K$,
we must make appropriate choices for the matrix $X$ in order to comply with the
experimental bounds, while at the same time obtain deviations that are sizeable
enough to be detected  experimentally in the near future. It is our aim to show that,
contrary to common wisdom, we can achieve this result with at least one of the heavy
neutrinos with a mass at the TeV scale, without requiring unnaturally small Yukawa
couplings and still have light neutrino masses not exceeding one eV.

Deviations from unitarity \cite{Antusch:2006vwa}, \cite{FernandezMartinez:2007ms},\cite{Antusch:2014woa},
\cite{Fernandez-Martinez:2016lgt}, \cite{Blennow:2016jkn}
of $K$ have been parametrised as the product of an
Hermitian matrix by a unitary matrix \cite{Fernandez-Martinez:2016lgt}:
\begin{equation}
K = ({1\>\!\!\!\mathrm{I}} - \eta ) V
\label{khu}
\end{equation}
where $\eta $ is an Hermitian matrix with small entries.
In order to identify the different components of our matrix $K$, given in Eq.~(\ref{v1}),
with the parametrisation of Eq.~(\ref{khu}) we rewrite $K$ as:
\begin{equation}
K = U_{K}U^{\dagger }\left( \sqrt{{1\>\!\!\!\mathrm{I}}\,+X^{\dagger }X}\
\right) ^{-1} = \left[ U_{K}U^{\dagger }\left( \sqrt{{1\>\!\!\!\mathrm{I}}\,+X^{\dagger }X}\
\right) ^{-1} U U_{K}^{\dagger } \right]  U_{K}U^{\dagger }
\label{kz}
\end{equation}
inside the square brackets  we wrote the Hermitian matrix that we identify with
$({1\>\!\!\!\mathrm{I}} - \eta )$, and which will parametrise the deviations from unitarity.
The matrix $V \equiv U_{K}U^{\dagger }$ is a unitarity matrix which is identified with $U_{PMNS}$
obtained from the standard parametrisation \cite{Patrignani:2016xqp} for a unitary matrix.
One can also write:
\begin{equation}
\left[ U_{K}U^{\dagger }\left( \sqrt{{1\>\!\!\!\mathrm{I}}\,+X^{\dagger }X}\
\right) ^{-1} U U_{K}^{\dagger } \right] \equiv \left[ U_{K} \left(
\sqrt{{1\>\!\!\!\mathrm{I}}\,+ d^2_X }\
\right) ^{-1} U_{K}^{\dagger } \right]
\label{eq24}
\end{equation}
where $d^2_X$ is a $3 \times 3$ diagonal matrix, introduced in Eq.~(\ref{uwxx}).
Identifying the second expression of Eq.~(\ref{eq24}) to $({1\>\!\!\!\mathrm{I}} - \eta )$
we derive:
\begin{equation}
\label{etaa}
\eta = {1\>\!\!\!\mathrm{I}} - U_{K} \left( \sqrt{{1\>\!\!\!\mathrm{I}}\,+ d^2_X }\
\right) ^{-1} U_{K}^{\dagger } \approx \frac{1}{2} U_K \ d^2_X \  U_{K}^{\dagger }
\end{equation}
for small $d^2_X$ .
The matrix $U_{PMNS}$ is then fixed making use of the present best fit values obtained from
a global analysis based on the assumption of unitarity. As pointed out in
\cite{Fernandez-Martinez:2016lgt}, from the phenomenological point of
view it is very useful to parametrise $K$ with the unitary matrix on the right, due to the
fact that experimentally it is not possible to determine which physical light neutrino is
produced, and therefore, one must sum over the neutrino indices. As a result, most observables
depend on $K K^\dagger$ which depends on the following combination:
\begin{equation}
\left( K K^\dagger \right)_{\alpha \beta} = \delta_{\alpha \beta} -
2 \eta_{\alpha \beta} + \mathcal{O}(\eta^2_{\alpha \beta})
\end{equation}

The standard parametrisation  for $U_{PMNS}$ is given by \cite{Patrignani:2016xqp}:
\begin{eqnarray}
U_{PMNS}=\left(
\begin{array}{ccc}
c_{12} c_{13} & s_{12} c_{13} & s_{13} e^{-i \delta}  \\
-s_{12} c_{23} - c_{12} s_{23} s_{13}   e^{i \delta}
& \quad c_{12} c_{23}  - s_{12} s_{23}  s_{13} e^{i \delta} \quad
& s_{23} c_{13}  \\
s_{12} s_{23} - c_{12} c_{23} s_{13} e^{i \delta}
& -c_{12} s_{23} - s_{12} c_{23} s_{13} e^{i \delta}
& c_{23} c_{13}
\end{array}\right) \cdot  P
\label{std}
\end{eqnarray}
with P given by
\begin{equation}
P=\mathrm{diag} \ (1,e^{i\alpha_{21}}, e^{i\alpha_{31}})
\end{equation}
where $c_{ij} = \cos \theta_{ij}$, $s_{ij} = \sin \theta_{ij}$
and $\delta$ is a Dirac-type
CP violating phase, while $\alpha_{21}$, $\alpha_{31}$ denote Majorana
phases. Neutrino oscillation experiments are not sensitive to these
factorisable phases.

There are several groups performing global phenomenological fits on
$\theta_{12}$, $\theta_{23}$, $\theta_{13}$ and $\delta$, as well as on
neutrino mass differences \cite{deSalas:2017kay},
\cite{Capozzi:2017ipn}, \cite{Esteban:2016qun}.

The specific bounds vary slightly from group to group.
For definiteness we present in Table 1 the present bounds on neutrino
masses and leptonic mixing from \cite{deSalas:2017kay}. The quantities
$\Delta m^2_{ij}$ are defined by $(m^2_i - m^2_j)$.
\begin{center}
\begin{table}[h]
\caption{Neutrino oscillation parameter summary from \cite{deSalas:2017kay}.
For $\Delta m^2_{31}$,
$\sin^2 \protect\theta_{23}$ , $\sin^2 \protect\theta_{13}$, and $\protect
\delta$ the upper (lower) row corresponds to normal (inverted) neutrino mass
hierarchy.
${}^a$There is a local minimum in the second octant, at  $\sin^2 \theta_{23}= 0.596 $
with $\Delta \chi^2 = 2.08$ with respect to the global minimum.
${}^b$ There is a local minimum in the first octant, at  $\sin^2 \theta_{23}= 0.426 $
with $\Delta \chi^2 = 1.68$ with respect to the global minimum for IO.}
\label{reps}
\begin{tabular}{ccc}
\hline\hline
Parameter & Best fit & $1 \sigma $ range \\ \hline
$\Delta m^2_{21}$ $[10^{-5} eV^2 ] $ & 7.56 & 7.37 -- 7.75\\
$|\Delta m^2_{31}|$ $[10^{-3} eV^2 ] (NO)$ & 2.55 & 2.41 -- 2.59 \\
$|\Delta m^2_{31}|$ $[10^{-3} eV^2 ] (IO) $ & 2.49 & 2.45 -- 2.53 \\
$\sin^2 \theta_{12}$ & 0.321 & 0.305 -- 0.339 \\
$\sin^2 \theta_{23} (NO)$ & 0.430  & 0.412 -- $0.450^a$ \\
$\sin^2 \theta_{23} (IO)$ & 0.596 & 0.576 -- $0.614^b$ \\
$\sin^2 \theta_{13}$ (NO)& 0.02155 & 0.02080 --0.02245 \\
$\sin^2 \theta_{13}$ (IO) & 0.02140 & 0.02055 -- 0.02222 \\
$\delta$ (NO) & 1.40 $\pi $ & 1.20 --1.71 $\pi$ \\
$\delta$ (IO) & 1.44 $\pi$ & 1.70 --1.21 $\pi$ \\ \hline
\end{tabular}
\end{table}
\end{center}
In Ref.~\cite{Fernandez-Martinez:2016lgt} global constraints are derived on
the matrix $\eta$ through a fit of twenty eight observables including
the W boson mass, the effective mixing weak angle $\theta_W$, several
ratios of $Z$ fermionic decays, the invisible width of the $Z$, several
ratios of weak decays constraining EW universality, weak decays constraining
CKM unitarity and some radiative lepton flavour violating (LFV) processes.
The final result is translated into:
\begin{eqnarray}
\left| 2 \eta_{\alpha \beta} \right| \leq \left(
\begin{array}{ccc}
2.5 \times 10^{-3} & 2.4 \times 10^{-5}  & 2.7 \times 10^{-3}    \\
2.4 \times 10^{-5} & 4.0 \times 10^{-4}  & 1.2 \times 10^{-3}   \\
2.7 \times 10^{-3}  & 1.2 \times 10^{-3}  & 5.6 \times 10^{-3}
\end{array} \right)
\label{hen}
\end{eqnarray}
Ref.~\cite{Fernandez-Martinez:2016lgt} also compares these bounds with
those of previous studies \cite{Antusch:2014woa}, \cite{Antusch:2015mia},
pointing out that in general there is good agreement. Variations in the
scale of the masses of the heavy neutrinos lead to small effects and therefore,
do not significantly change our analysis.

\subsection{The elements of the neutrino Dirac mass matrix $m$ and
deviations from unitarity}

In this subsection, we show that there is a correlation among:
\begin{itemize}
\item The size of deviations from unitarity of the $3 \times 3$
leptonic mixing matrix.
\item The mass of the lightest heavy neutrino.
\end{itemize}
From Eq.~(\ref{13c}) we get
\begin{eqnarray}
m= K X^\dagger D {Z^*}^{-1}
\label{mmtu}
\end{eqnarray}
the experimental fact that $K$ is almost unitary implies that $Z$ is also
almost unitary. Therefore the Dirac mass matrix
$m$ is of the same order as $X$ times $D$. Notice
that the scale of $D$ may be of the order of the top quark mass, so that indeed
the Yukawa couplings need not be extremely small.

The elements of the neutrino Dirac mass
matrix $m$ are connected to the deviations from unitarity of the $3 \times 3$
leptonic mixing matrix.
From Eq.~(\ref{mmtu}) together with Eqs.~(\ref{esta}), (\ref{kz}) and     (\ref{eq24}), we obtain:
\begin{equation}
m = U_{K}\left( \sqrt{\left( {1\>\!\!\!\mathrm{I}}\,+d_{X}^{2}\right) }\ \right)
^{-1}d_{X}\ W^{\dagger }\ D\ W^{\ast }\left( \sqrt{\left( {1\>\!\!\!\mathrm{I
}}\,+d_{X}^{2}\right) }\ \right) W_{Z}^{T} \label{mss}
\end{equation}
where we have used $d_{X}=$ $W^{\dagger }\ X\ U$ from Eq. (\ref{uwxx}).
Thus, we find
\begin{equation}
Tr\left[ mm^{\dagger}\right] =
Tr\left[ \left( \sqrt{\left( {1\>\!\!\!\mathrm{I}}\,+d_{X}^{2}\right) }\
\right) ^{-1}d_{X}\ W^{\dagger }\ D\ W^{\ast }\ \left( {1\>\!\!\!\mathrm{I}}%
\,+d_{X}^{2}\right) \ W^{T}\ D\ W\ d_{X}\ \left( \sqrt{\left( {1\>\!\!\!
\mathrm{I}}\,+d_{X}^{2}\right) }\ \right) ^{-1}\right]
\label{mmt}
\end{equation}
As previously emphasised, deviations from $3 \times 3$ unitarity in $U_{PMNS}$
are controlled by the matrix $X$, as it is clear from Eq.~(\ref{3lin}). For
$X=0$, there are no deviations from unitarity. Small deviations from unitarity
correspond to $d_X$ small and, in that case, one has to a very good
approximation,
\begin{equation}
Tr\left[ mm^{\dagger}\right] =
Tr\left[ d_{X}\ W^{\dagger }D^{2}W\ d_{X}\right]
\end{equation}
which can be written as:
\begin{equation}
\begin{array}{l}
Tr\left[ mm^{\dagger}\right] =
d_{X_{1}}^{2}\left( M_{1}^{2}\left\vert W_{11}\right\vert
^{2}+M_{2}^{2}\left\vert W_{21}\right\vert ^{2}+M_{3}^{2}\left\vert
W_{31}\right\vert ^{2}\right) + \\
d_{X_{2}}^{2}\left( M_{1}^{2}\left\vert W_{12}\right\vert
^{2}+M_{2}^{2}\left\vert W_{22}\right\vert ^{2}+M_{3}^{2}\left\vert
W_{32}\right\vert ^{2}\right) + \\
d_{X_{3}}^{2}\left( M_{1}^{2}\left\vert W_{13}\right\vert
^{2}+M_{2}^{2}\left\vert W_{23}\right\vert ^{2}+M_{3}^{2}\left\vert
W_{33}\right\vert ^{2}\right)
\end{array}
\label{mmt1}
\end{equation}
It can be shown that, using the properties of orthogonal complex matrices, only one of the
$d_{X_i}$, corresponding to $d_{X_3}$, can have a significant value (e.g.
$d_{X_3} \approx 10^{-2}$),  while the other two are negligible.
Thus, we find in good approximation
\begin{equation}
Tr\left[ mm^{\dagger}\right] =
d_{X_{3}}^{2}\left( M_{1}^{2}\left\vert W_{13}\right\vert
^{2}+M_{2}^{2}\left\vert W_{23}\right\vert ^{2}+M_{3}^{2}\left\vert
W_{33}\right\vert ^{2}\right)  \label{mmt2}
\end{equation}
or using the unitary of $W$
\begin{equation}
Tr\left[ mm^{\dagger}\right] =
d_{X_{3}}^{2}M_{1}^{2}\left( 1+\left( \frac{M_{2}^{2}}{M_{1}^{2}}-1\right)
\left\vert W_{23}\right\vert ^{2}+\left( \frac{M_{3}^{2}}{M_{1}^{2}}
-1\right) \left\vert W_{33}\right\vert ^{2}\right)
\label{mmt3}
\end{equation}
which, with the choice $M_{3}\geq M_{2}\geq M_{1}$, leads to
\begin{equation}
d_{X_{3}}^{2}M_{1}^{2}\leq Tr\left[ mm^{\dagger }\right] =\sum_{i,j} \left\vert
m_{ij}\right\vert ^{2}   \label{dim}
\end{equation}
From Eq.~(\ref{dim}), it is clear that for significant values of $d_{X_3}$,
$M_1$ cannot be too large in order to avoid a too large value of
$Tr\left[ mm^{\dagger}\right] $, which in turn would imply that
at least one of the $\left\vert m_{ij}\right\vert ^{2}$ is too large.
This can be seen in both Fig.~1 and Fig.~3 where we plot $\frac{1}{2}d_{X_{3}}^{2}$ versus $M_1$.
Significant values of $d_{X_{3}}^{2}$  can only be obtained for
$M_1 \leq 1 - 2$ TeV.
\begin{figure}[h!]

\includegraphics[width=1.\linewidth]{aM1-dev.jpg}

\caption{Maximum deviations from unitarity as a function of $M_1$, generated under
the condition that Tr$( mm^\dagger)\leq m_t^2$ and $|\eta_{12}|\leq 10^{-4}$.  }
\label{fig:sub1}

\end{figure}

In all our plots, we require $Tr\left[ mm^{\dagger }\right] \leq
m_{t}^{2}$. We consider the case of normal ordering
and vary over the values of light neutrinos masses
$m_{i}$, up to $m_3 = 0.5$ eV.  Concerning
the heavy Majorana masses $M_i$, we allow $M_3$ to reach values of the order of
$10^4 m_t$ and allow for all possible
forms of $O_{c}$. In Figs.~1,~2 we impose the condition that $\left\vert \eta _{12}\right\vert \leq 10^{-4}$, while in Figs.~3,~4, we chose
$\left\vert \eta _{12}\right\vert \leq 2\times 10^{-5}$.

In Figs.~2,~4, we also plot the $|\eta_{11}|$ deviations from unitarity.

\begin{figure}[h!]

\includegraphics[width=1.\linewidth]{aM1-eta11.jpg}

\caption{$|\eta_{11}|$ deviations from unitarity as a function of $M_1$, generated under the condition that Tr$(mm^\dagger)\leq m_t^2$ and $|\eta_{12}|\leq 10^{-4}$. }
\label{fig:sub2}

\end{figure}

\begin{figure}[h!]

\includegraphics[width=1.\linewidth]{bM1-dev.jpg}

\caption{Maximum deviations from unitarity as a function of $M_1$, generated under the condition
that Tr$(mm^\dagger)\leq m_t^2$ and $|\eta_{12}|\leq 2\times 10^{-5}$.  }
\label{fig:sub3}

\end{figure}

\begin{figure}[h!]

\includegraphics[width=1.\linewidth]{bM1-eta11.jpg}

\caption{$|\eta_{11}|$ deviations from unitarity as a function of $M_1$, generated under the condition that Tr$(mm^\dagger)\leq m_t^2$ and $|\eta_{12}|\leq 2\times 10^{-5}$. }
\label{fig:sub4}

\end{figure}

\section{Numerical Examples}

In this section, we present some illustrative results of our numerical analysis showing
that it is possible to obtain the observed pattern of neutrino masses and mixing without
requiring that all masses of the heavy Majorana neutrinos, $M_i$, be much larger
than the electroweak scale.
We have realistic examples even when all the three heavy neutrinos have masses below 2 TeV.
One might expect that lowering the scale of the heavy neutrino masses would
result in the need for extremely small Yukawa couplings for the Dirac mass terms,
thus defeating the rationale for the seesaw mechanism. However, this is not the case
and, to illustrate, we include for each example the corresponding moduli of the entries
of the neutrino Dirac mass matrices $m$ and the trace of the product $mm^\dagger$.  \\

\begin{table}[th]
\caption{Examples of deviations from unitarity, for the case of NH, expressed by the
Hermitian matrix $\eta$.
We consider three different hierarchies for the heavy Majorana neutrino masses $M_i$. For
each hierarchy we give two examples, where we vary the choice of the matrix $X$. The third
column contains the Dirac-type neutrino mass matrix expressed in GeV units. The light
neutrino  masses $m_i$ are given in eV.}
\begin{center}
\setlength{\tabcolsep}{0.5pc}
\resizebox{\textwidth}{!}{\begin{tabular}{|c|c|c|c|c|c|}
\hline
Neutrino Masses \\ $m_1=0.005050$\ , \  $m_2=0.01005$\ , \  $m_3=0.05075$& $|\eta|$ &$|m|$\\
\hline
$\begin{pmatrix}
	M_1= 3 \ m_t\\
	M_2= 6\ m_t\\\
	M_3=9\ m_t\
\end{pmatrix}$ &
$\begin{pmatrix}
	1.14\times 10^{-3} & 9.11\times 10^{-6}& 1.39\times 10^{-3} \\
	. & 7.27\times 10^{-8}& 1.11\times 10^{-5} \\
	.& . & 1.68\times 10^{-3}
\end{pmatrix}$
 &
$\begin{pmatrix}
	35.183 &  6 \times 10^{-5} &  7 \times 10^{-11} \\
	0.280 & 2 \times 10^{-4} & 3 \times 10^{-10} \\
	42.707 & 1 \times 10^{-4} & 3 \times 10^{-10}
\end{pmatrix}$ \ \mbox{tr} ($mm^\dagger$) = 0.102 $m^2_t$
\\
\hline

$\begin{pmatrix}
	M_1= 3 \ m_t\\
	M_2= 6\ m_t\\\
	M_3=9\ m_t\
\end{pmatrix}$ &
$\begin{pmatrix}
	5.6\times 10^{-7} & 9.87\times 10^{-6}& 3.28\times 10^{-5} \\
	. &1.74\times 10^{-4}& 5.78\times 10^{-4} \\
	.& . & 1.92\times 10^{-3}
\end{pmatrix}$ &
$\begin{pmatrix}
	0.779 &  4 \times 10^{-5} & 1 \times 10^{-10}  \\
	13.735 & 2 \times 10^{-4} & 4 \times 10^{-10} \\
	45.63&  9 \times 10^{-5} & 3 \times 10^{-10}
\end{pmatrix}$ \ \mbox{tr} ($mm^\dagger$) = 0.075 $m^2_t$

\\

\hline

$\begin{pmatrix}
	M_1=3\ m_t\ \\
	M_2=6\  m_t\\\
	M_3=1000\ m_t\
\end{pmatrix}$   &
$\begin{pmatrix}
	1.12\times 10^{-3} & 8.91\times 10^{-6}& 1.36\times 10^{-3} \\
	. & 7.10\times 10^{-8}& 1.08\times 10^{-5} \\
	.& . & 1.64\times 10^{-3}
\end{pmatrix}$
 & $\begin{pmatrix}
	34.779 &  2 \times 10^{-4} & 7 \times 10^{-11}  \\
	0.277 & 0.002 & 3 \times 10^{-10}  \\
	42.217& 0.002 & 3 \times 10^{-10}
\end{pmatrix}$ \ \mbox{tr} ($mm^\dagger$) =  0.099 $m^2_t$

\\

\hline

$\begin{pmatrix}
	M_1=3\ m_t\ \\
	M_2=6\  m_t\\\
	M_3=1000\ m_t\
\end{pmatrix}$   &
$\begin{pmatrix}
	5.6\times 10^{-7} & 9.87\times 10^{-6}& 3.28\times 10^{-5} \\
	. &1.74\times 10^{-4}& 5.78\times 10^{-4} \\
	.& . & 1.92\times 10^{-3}
\end{pmatrix}$
 & $\begin{pmatrix}
	0.779 &  5 \times 10^{-4}  & 1 \times 10^{-10} \\
	13.735 & 0.002 & 4 \times 10^{-10} \\
	45.643& 0.002 & 3 \times 10^{-10}
\end{pmatrix}$ \ \mbox{tr} ($mm^\dagger$) = 0.075 $m^2_t$

\\

\hline

$\begin{pmatrix}
	M_1=3\ m_t\ \\
	M_2=60\ m_t\\\
	M_3=1000\ m_t\
\end{pmatrix}$   &
$\begin{pmatrix}
	1.09\times 10^{-3} & 8.71\times 10^{-6}& 1.33\times 10^{-3} \\
	. & 6.95\times 10^{-8}& 1.06\times 10^{-5} \\
	.& . & 1.61\times 10^{-3}
\end{pmatrix}$
 & $\begin{pmatrix}
	108.753 &   2 \times 10^{-4}  & 7 \times 10^{-11}  \\
	0.867 & 0.002 & 3 \times 10^{-10}  \\
	132.014 & 0.002 & 3 \times 10^{-10}
\end{pmatrix}$ \ \mbox{tr} ($mm^\dagger$) = 0.9708 $m^2_t$

\\
\hline

$\begin{pmatrix}
	M_1=3\ m_t\ \\
	M_2=60\ m_t\\\
	M_3=1000\ m_t\
\end{pmatrix}$   &
$\begin{pmatrix}
	5.27\times 10^{-7} & 9.3\times 10^{-6}& 3.09\times 10^{-5} \\
	. &1.64\times 10^{-4}& 5.45\times 10^{-4} \\
	.& . & 1.81\times 10^{-3}
\end{pmatrix}$
 & $\begin{pmatrix}
	2.391 &  5 \times 10^{-4} & 1 \times 10^{-10}  \\
	42.159 & 0.002 & 4 \times 10^{-10} \\
	140.072& 0.002 & 3 \times 10^{-10}
\end{pmatrix}$ \ \mbox{tr} ($mm^\dagger$) = 0.710 $m^2_t$

\\

\hline
\end{tabular}}
\end{center}
\end{table}

In Table 2, we include examples with normal ordering of light neutrino masses, for
a particular fixed value of these three masses, common to all examples.
We consider three different hierarchies for the heavy Majorana neutrino masses $M_i$. For
each choice of heavy neutrino mass hierarchy we give two examples, where we vary
the matrix $O_c$. The other set of free parameters in our analysis are the
Majorana-type phases entering in the choice of matrix $\Omega$.

In all examples given in Table 2 the orthogonal complex matrix is of the form:
\begin{equation}
O_{c}=\left(
\begin{array}{ccc}
0 & \sqrt{x^{2}+1} & i\ x  \\
0 & i\ x & -\sqrt{x^{2}+1}  \\
1 & 0 & 0
\end{array}
\right) \  \left(
\begin{array}{ccc}
\cos \tau & 0 & \sin \tau  \\
0 & 1 & 0  \\
-\sin \tau & 0 & \cos \tau
\end{array}
\right)
\label{OcO}
\end{equation}
with different choices of $x$ and $\tau$ for the different cases, respectively:
\begin{equation}
\begin{array}{cc}
x= 2.80 \times 10^{5}  & \ \tau = \frac{\pi}{2.9} \\
x= 2.48 \times 10^{5}  & \ \tau = \frac{\pi}{2.8} \\
x= 2.77 \times 10^{5}  & \ \tau = \frac{\pi}{2.9} \\
x= 2.48 \times 10^{5} & \ \tau = \frac{\pi}{2.8} \\
x= 3.27 \times 10^{5} & \ \tau = \frac{\pi}{2.9} \\
x= 2.88 \times 10^{5}  & \ \tau = \frac{\pi}{2.8}
\end{array}
\end{equation}

The modulus of $m$ is obtained from Eq.~(\ref{mmtu}) and requires
a choice for the $W_Z$ matrix. This choice has implications for the entries of the
blocks $S$ and $Z$ of the matrix $\cal{V}$. Since, at the moment, there are no direct
experimental constraints on these entries, providing guidance for this choice, we made
the simplest one, by fixing $W_Z$ to be equal to the identity. With a
different choice of $W_Z$  we could in principle homogenise the orders of magnitude of
the entries of $m$ so that all of the Yukawa couplings would be of the same
order of magnitude or close. \\

In Table 3 we include examples with inverted ordering of light neutrino masses, for
a particular fixed value of these three masses, common to all examples.  We consider
the same three different hierarchies for the heavy Majorana neutrino masses $M_i$, as
in the cases of Table 2. For each choice of heavy neutrino
mass hierarchy we give two examples, where we vary the matrix $O_c$. In the second,
fourth and sixth examples of Table 3, $O_c$ is of the form given by Eq.~(\ref{OcO}).
In the first, third and fifth examples $O_c$ is of the form:
\begin{equation}
O_{c}=\left(
\begin{array}{ccc}
\sqrt{x^{2}+1}& 0 & i\ x  \\
i\ x & 0 & -\sqrt{x^{2}+1}  \\
0 & 1 & 0
\end{array}
\right) \  \left(
\begin{array}{ccc}
\cos \tau  & \sin \tau & 0  \\
-\sin \tau &  \cos \tau & 0 \\
0 & 0 & 1
\end{array}
\right)
\label{1c1}
\end{equation}
The choices of $x$ and $\tau$ are respectively:
\begin{equation}
\begin{array}{cc}
x= 1.68 \times 10^{5} &  \ \tau = \frac{\pi}{10} \\
x= 9.65 \times 10^{4}  &  \ \tau = \frac{\pi}{3.9} \\
x= 1.68 \times 10^{5}  &  \ \tau = \frac{\pi}{10} \\
x= 9.65 \times 10^{4}  &  \ \tau = \frac{\pi}{3.9} \\
x=1.99 \times 10^{5}  &  \ \tau = \frac{\pi}{10} \\
x= 1.14 \times 10^{5}  & \ \tau = \frac{\pi}{3.9}
\end{array}
\end{equation}
\begin{table}[th]
\caption{Examples of deviations from unitarity, for the case of IO, expressed by the
Hermitian matrix $\eta$.
We consider three different hierarchies for the heavy Majorana neutrino masses $M_i$. For
each hierarchy we give two examples, where we vary the choice of the matrix $X$. The third
column contains the Dirac-type neutrino mass matrix expressed in GeV units. The light
neutrino masses $m_i$ are given in eV.}
\begin{center}
\setlength{\tabcolsep}{0.5pc}
\resizebox{\textwidth}{!}{\begin{tabular}{|c|c|c|c|c|c|}
\hline
Neutrino Masses \\ $m_1= 0.05064$\ , \  $m_2= 0.05138$\ , \  $m_3= 0.00864$& $|\eta|$ &$|m|$\\
\hline
$\begin{pmatrix}
	M_1= 3 \ m_t\\
	M_2= 6\ m_t\\\
	M_3=9\ m_t\
\end{pmatrix}$ &
$\begin{pmatrix}
	1.17\times 10^{-3} & 9.96\times 10^{-6}& 1.2\times 10^{-3} \\
	. & 8.47\times 10^{-8}& 1.02\times 10^{-5} \\
	.& . & 1.22\times 10^{-3}
\end{pmatrix}$
 &
$\begin{pmatrix}
	35.630 &  2 \times 10^{-4} &  5 \times 10^{-10} \\
	0.303 & 2 \times 10^{-4} & 2.5 \times 10^{-10} \\
	36.356 & 2 \times 10^{-4} & 3 \times 10^{-10}
\end{pmatrix}$ \ \mbox{tr} ($mm^\dagger$) = 0.086 $m^2_t$
\\
\hline

$\begin{pmatrix}
	M_1= 3 \ m_t\\
	M_2= 6\ m_t\\\
	M_3=9\ m_t\
\end{pmatrix}$ &
$\begin{pmatrix}
	1.32\times 10^{-7} & 5.34\times 10^{-6}& 1.08\times 10^{-5} \\
	. &2.15\times 10^{-4}& 4.35\times 10^{-4} \\
	.& . & 8.80\times 10^{-4}
\end{pmatrix}$ &
$\begin{pmatrix}
	0.379 &  2 \times 10^{-4} & 1 \times 10^{-9}  \\
	15.272 & 1 \times 10^{-4} & 3 \times 10^{-10} \\
	30.896 &  2 \times 10^{-5} & 2 \times 10^{-10}
\end{pmatrix}$ \ \mbox{tr} ($mm^\dagger$) = 0.039 $m^2_t$

\\

\hline

$\begin{pmatrix}
	M_1=3\ m_t\ \\
	M_2=6\  m_t\\\
	M_3=1000\ m_t\
\end{pmatrix}$   &
$\begin{pmatrix}
	1.17\times 10^{-3} & 9.96\times 10^{-6}& 1.2\times 10^{-3} \\
	. & 8.47\times 10^{-8}& 1.02\times 10^{-5} \\
	.& . & 1.22\times 10^{-3}
\end{pmatrix}$

 & $\begin{pmatrix}
	35.630 &  2 \times 10^{-3} & 5 \times 10^{-10}  \\
	0.303 & 0.002 & 3 \times 10^{-10}  \\
	36.356 & 0.002 & 3 \times 10^{-10}
\end{pmatrix}$ \ \mbox{tr} ($mm^\dagger$) = 0.086 $m^2_t$

\\

\hline

$\begin{pmatrix}
	M_1=3\ m_t\ \\
	M_2=6\  m_t\\\
	M_3=1000\ m_t\
\end{pmatrix}$   &
$\begin{pmatrix}
	1.32\times 10^{-7} & 5.34\times 10^{-6}& 1.08\times 10^{-5} \\
	. &2.15\times 10^{-4}& 4.35\times 10^{-4} \\
	.& . & 8.80\times 10^{-4}
\end{pmatrix}$
 & $\begin{pmatrix}
	0.379 &  0.002  & 1 \times 10^{-9} \\
	15.272 & 0.001 & 3 \times 10^{-10} \\
	30.896 & 0.0004 & 2 \times 10^{-10}
\end{pmatrix}$ \ \mbox{tr} ($mm^\dagger$) = 0.039 $m^2_t$

\\

\hline

$\begin{pmatrix}
	M_1=3\ m_t\ \\
	M_2=60\ m_t\\\
	M_3=1000\ m_t\
\end{pmatrix}$   &
$\begin{pmatrix}
	1.16\times 10^{-3} & 9.83\times 10^{-6}& 1.18\times 10^{-3} \\
	. & 8.36\times 10^{-8}& 1.0\times 10^{-5} \\
	.& . & 1.20\times 10^{-3}
\end{pmatrix}$
 & $\begin{pmatrix}
	111.946 &  0.002 & 5 \times 10^{-10}  \\
	0.952 & 0.002 & 3 \times 10^{-10}  \\
	114.227 & 0.002 & 3 \times 10^{-10}
\end{pmatrix}$ \ \mbox{tr} ($mm^\dagger$) = 0.849 $m^2_t$

\\
\hline

$\begin{pmatrix}
	M_1=3\ m_t\ \\
	M_2=60\ m_t\\\
	M_3=1000\ m_t\
\end{pmatrix}$   &
$\begin{pmatrix}
	1.29\times 10^{-7} & 5.19\times 10^{-6}& 1.05\times 10^{-5} \\
	. &2.09\times 10^{-4}& 4.23\times 10^{-4} \\
	.& . & 8.56\times 10^{-4}
\end{pmatrix}$
 & $\begin{pmatrix}
	1.181 &  0.002 & 1 \times 10^{-9}  \\
	47.622 & 0.001 & 3 \times 10^{-10} \\
	96.343 & 0.0003 & 3 \times 10^{-10}
\end{pmatrix}$ \ \mbox{tr} ($mm^\dagger$) = 0.383 $m^2_t$

\\

\hline
\end{tabular}}
\end{center}
\end{table}

In all our examples the matrix $X$
will have several entries of order at most $10^{-2}$.

\section{Conclusions}
We have studied the possibility of having significant deviations from $3\times 3$ unitarity of
the leptonic mixing matrix in the framework of type-I seesaw mechanism.
The analysis was done in the framework of an extension of the Standard Model where
three right-handed neutrinos are added to the spectrum. We have shown that the
$6 \times 6$
unitary leptonic mixing matrix ${\cal V}$ can be written in terms of two unitary $3\times 3$ matrices
and a matrix denoted X,  which controls the deviations from unitarity of the $3\times 3$
PMNS matrix. This parametrisation of the matrix ${\cal V}$,  played a crucial role in showing
that one may have significant deviations from $3\times 3$ unitarity while conforming to all present
data on neutrino masses and mixing, as well as respecting all stringent bounds on
deviations from $3\times 3$ unitarity of the PMNS matrix.

We have presented specific examples where the above deviations from unitarity are
sufficiently large to have the potential for being observed at the next round
of experiments. An important feature of our analysis is the fact that the mass of
the lightest heavy sterile neutrino is, in principle, within experimental reach.
This can be achieved without unnaturally small neutrino Yukawa couplings, thus
defying the conventional wisdom that heavy neutrino masses many orders of magnitude
above the electroweak scale are needed for this purpose.

\section*{Appendix}

We can write the Hermitian matrices $\left( {1\>\!\!\!\mathrm{I}}
\,+X^{\dagger }X\right) $ and $\left( {1\>\!\!\!\mathrm{I}}\,+X\ X^{\dagger
}\right) $ as,
\begin{equation}
\begin{array}{c}
{1\>\!\!\!\mathrm{I}}\,+X^{\dagger }X=U\ \left( {1\>\!\!\!\mathrm{I}}
\,+d_{X}^{2}\right) \ U^{\dagger } \\
\\
{1\>\!\!\!\mathrm{I}}\,+X\ X^{\dagger }=W\ \left( {1\>\!\!\!\mathrm{I}}
\,+d_{X}^{2}\right) \ W^{\dagger }
\end{array}
\label{uw}
\end{equation}
where $U$ and $W$ are, by definition the unitary matrices
which diagonalise the Hermitian matrices $({1\>\!\!\!\mathrm{I}}\,+X^{\dagger}X)$
and $({1\>\!\!\!\mathrm{I}}\,+X\ X^{\dagger})$ respectively and
$d_{X}^{2}$ is a diagonal matrix. Inserting Eq.~(\ref{uw}) into Eq.~(\ref{kt})
we obtain:
\begin{equation}
\begin{array}{c}
K\ U\ \left( {1\>\!\!\!\mathrm{I}}\,+d_{X}^{2}\right) \ U^{\dagger }\
K^{\dagger }=K\ U\ \sqrt{\left( {1\>\!\!\!\mathrm{I}}\,+d_{X}^{2}\right) }\
\cdot \sqrt{\left( {1\>\!\!\!\mathrm{I}}\,+d_{X}^{2}\right) }U^{\dagger }\
K^{\dagger }={1\>\!\!\!\mathrm{I}} \\
\\
Z\  W\ \left( {1\>\!\!\!\mathrm{I}}\,+d_{X}^{2}\right) \ W^{\dagger }\
Z^{\dagger }=Z\ W\ \sqrt{\left( {1\>\!\!\!\mathrm{I}}\,+d_{X}^{2}\right) }\
\cdot \sqrt{\left( {1\>\!\!\!\mathrm{I}}\,+d_{X}^{2}\right) }\ W^{\dagger }\
Z^{\dagger }={1\>\!\!\!\mathrm{I}}
\end{array}
\end{equation}
We therefore conclude that $K\ U\ \sqrt{\left( {1\>\!\!\!\mathrm{I}}
\,+d_{X}^{2}\right) }=U_{K}$ and $Z\ W\ \sqrt{\left( {1\>\!\!\!\mathrm{I}}
\,+d_{X}^{2}\right) }=W_{Z}$ are unitary matrices, and thus also that
\begin{equation}
\begin{array}{l}
K=U_{K}\ \left( \sqrt{\left( {1\>\!\!\!\mathrm{I}}\,+d_{X}^{2}\right) }\
\right) ^{-1}U^{\dagger }=U_{K}U^{\dagger }\left( U\left( \sqrt{\left( {%
1\>\!\!\!\mathrm{I}}\,+d_{X}^{2}\right) }\ \right) ^{-1}U^{\dagger }\right)
\\
\\
Z=W_{Z}\ \left( \sqrt{\left( {1\>\!\!\!\mathrm{I}}\,+d_{X}^{2}\right) }\
\right) ^{-1}W^{\dagger }=W_{Z}\ W^{\dagger }\left( W\left( \sqrt{\left(
{1\>\!\!\!\mathrm{I}}\,+d_{X}^{2}\right) }\ \right) ^{-1}W^{\dagger }\right)
\end{array}
\label{kt1}
\end{equation}
or using the diagonalisation of the Hermitian matrices $\left( {1\>\!\!\!
\mathrm{I}}\,+X^{\dagger }X\right) $ and $\left( {1\>\!\!\!\mathrm{I}}\,+X\
X^{\dagger }\right) $ in Eq.(\ref{uw}), we write this, as
\begin{equation}
\begin{array}{c}
K=U_{K}\ \left( \sqrt{\left( {1\>\!\!\!\mathrm{I}}\,+d_{X}^{2}\right) }\
\right) ^{-1}U^{\dagger }=U_{K}U^{\dagger }\left( \sqrt{{1\>\!\!\!\mathrm{I}}
\,+X^{\dagger }X}\ \right) ^{-1} \\
\\
Z=W_{Z}\ \left( \sqrt{\left( {1\>\!\!\!\mathrm{I}}\,+d_{X}^{2}\right) }\
\right) ^{-1}W^{\dagger }=W_{Z}\ W^{\dagger }\left( \sqrt{{1\>\!\!\!\mathrm{I
}}\,+X\ X^{\dagger }}\ \right) ^{-1}\
\end{array}
\label{kt2}
\end{equation}
which lead to Eq.~(\ref{v1}).

\section*{Acknowledgments}
This work was partially supported by Funda\c c\~ ao para a Ci\^ encia e a
Tecnologia (FCT, Portugal) through the projects CERN/FIS-NUC/0010/2015,
CFTP-FCT Unit 777 \\
(UID/FIS/00777/2013) which are partially funded through
POCTI (FEDER), COMPETE, QREN and EU.
N.R.A has received funding/support from the European Union's
Horizon 2020 research and innovation programme under the Marie
Sklodowska-Curie grant agreement No 674896.
P.M.F.P. has a BIC-type Fellowship under project CERN/FIS-NUC/0010/2015.


\section{Introduction}

Understanding the physical parameters of dust aggregates is important in planetary science.
Specifically, the thermal conductivity of dust aggregates is key for determining the thermal evolution of planetary bodies, influencing the thermal evolution pathways of both rocky and icy planetesimals \citep[e.g.,][]{Henke+2012,Sirono2017}.
The thermal evolution and activity of cometary nuclei also depend on the thermal conductivity of icy aggregates \citep[e.g.,][]{Haruyama+1993,Guilbert-Lepoutre+2011}.

Dust aggregate thermal conductivity depends on many parameters, and many previous experimental studies have researched the thermal conductivity of dust aggregates with filling factors above $10^{-1}$.
The thermal conductivity of porous aggregates in vacuum is given by two terms: the thermal conductivity through the solid network $k_{\rm sol}$ and the thermal conductivity owing to radiative transfer $k_{\rm rad}$.
\citet{Krause+2011} showed that the thermal conductivity through the solid network $k_{\rm sol}$ is exponentially dependent on the filling factor of dust aggregates $\phi$ for $0.15 < \phi < 0.54$, and concluded that the coordination number of monomer grains $C$ influences the efficiency of heat flux within the aggregates.
\citet{Sakatani+2016} revealed that $k_{\rm sol}$ is also dependent on the contact radius between monomers $r_{\rm c}$.
The thermal conductivity owing to radiative transfer $k_{\rm rad}$ is affected by the temperature of dust aggregates $T$ and the mean free path of photons $l_{\rm p}$ \citep[e.g.,][]{Schotte1960,Merrill1969}.
Moreover, $l_{\rm p}$ depends on $R$ and $\phi$ when we apply the geometrical optics approximation for the evaluation of $l_{\rm p}$ \citep[e.g.,][]{Skorov+2011,Gundlach+2012}.

There are also several theoretical studies on the thermal conductivity of dust aggregates \citep[e.g.,][]{Chan+1973,Sirono2014,Sakatani+2017}.
However, no previous research has been conducted on the thermal conductivity of porous aggregates with filling factors of less than $10^{-1}$, although \citet{Kataoka+2013b} and \citet{Arakawa+2016b} revealed that the collisional growth of dust aggregates leads to planetesimal formation via highly porous aggregates with filling factors of much less than $10^{-1}$.
Therefore, the purpose of this study is to investigate the thermal conductivity and thermal evolution of fluffy dust aggregates in protoplanetary disks.

In this letter, we calculate thermal conductivity through the solid network $k_{\rm sol}$ for highly porous aggregates with filling factors in the range of $10^{-2}$ to $10^{-1}$.
We use the snapshot data of \citet{Kataoka+2013a} for calculation of $k_{\rm sol}$.
We then validate our results through a comparison with the experimental data of \citet{Krause+2011}.
We also derive the thermal conductivity owing to radiative transfer $k_{\rm rad}$ for porous aggregates of submicron-sized monomers.
Our results show that the thermal conductivity of highly porous aggregates is significantly lower than previously assumed.

\section{Method}

\subsection{Arrangement of monomer grains}

The arrangement of monomer grains depends on the coagulation history of the aggregates.
During initial dust aggregate coagulation in protoplanetary disks, both experimental \citep[e.g.,][]{Wurm+1998} and theoretical \citep[e.g.,][]{Kempf+1999} studies have shown that hit-and-stick collisions lead to the formation of fractal aggregates with a fractal dimension $D \sim 2$, which is called ballistic cluster-cluster aggregation \citep[BCCA;][]{Meakin1991}.
Furthermore, \citet{Kataoka+2013a} performed three-dimensional numerical simulations of static compression of BCCA aggregates constituted from 16384 spherical grains using a periodic boundary condition.
In this study, we use snapshots of the compressed BCCA aggregates calculated by \citet{Kataoka+2013a}.

\subsection{Temperature structure of the dust aggregate}

To calculate the thermal conductivity through the solid network of an aggregate $k_{\rm sol}$, we have to determine the temperature of each grain in a cubic periodic boundary.
We calculate the temperature of each grain using the method of \citet{Sirono2014}.
Here, we consider one-dimensional heat flow from the lower boundary plane to the upper boundary plane.
There are three choices regarding the pair of lower and upper planes, and we calculate $k_{\rm sol}$ from three directions.
Then, we average these values for each snapshot.

We define $R$ as the monomer radius and $L^{3}$ as the volume of each cubic space.
The location of the $i$-th grain $(x_{i}, y_{i}, z_{i})$ satisfies $|x_{i}| < L/2$, $|y_{i}| < L/2$, and $|z_{i}| < L/2$ for $i = 1, 2, ..., N$, where $N = 16384$ is the number of grains in the periodic boundary.
A sketch of a dust aggregate in a cubic periodic boundary is shown in Fig.\ \ref{fig1}.
Here, we assume that heat flow occurs along the z-direction.
The grains located in $- L/2 < z_{i} < - (L/2 - R)$ are on the lower boundary (\#1 in Fig.\ \ref{fig1}), and the grains located in $+ (L/2 - R) < z_{i} < + L/2$ are on the upper boundary (\#40 in Fig.\ \ref{fig1}).
If the $i$-th grain is located on the lower (upper) boundary, we add a new grain on the upper (lower) boundary.
The location of the new grain is $(x_{i}, y_{i}, z_{i} + L)$ if the $i$-th grain is located on the lower boundary (\#X in Fig.\ \ref{fig1}) and $(x_{i}, y_{i}, z_{i} - L)$ if the $i$-th grain is located on the upper boundary (\#Y in Fig.\ \ref{fig1}).
We set the temperature of grains located on the lower (\#1 and \#Y in Fig \ref{fig1}) and upper (\#40 and \#X in Fig \ref{fig1}) boundary as $T_{0} + {\Delta T}/2$ and $T_{0} - {\Delta T}/2$, respectively.

\begin{figure}[h]
\centering
\includegraphics[width=0.9\columnwidth]{fig1.eps}
\caption{
Sketch of a dust aggregate in a cubic periodic boundary.
The temperature of grains located on the lower (\#1 and \#Y) and upper (\#40 and \#X) boundary is set to $T_{0} + {\Delta T}/2$ and $T_{0} - {\Delta T}/2$, respectively.
The temperature of each grain is calculated by solving Eq.\ (\ref{eq1}) simultaneously for each grain.
}
\label{fig1}
\end{figure}

Heat flows through the monomer-monomer contacts, and for steady state conditions, the equation of heat balance at the $i$-th grain is given by
\begin{equation}
\sum_{j} F_{i, j} = 0,
\label{eq1}
\end{equation}
where $F_{i, j}$ is the heat flow from the $j$-th grain to the $i$-th grain, given by
\begin{equation}
F_{i, j} = H_{\rm c} {(T_{j} - T_{i})},
\end{equation}
where $H_{\rm c}$ is the heat conductance at the contact of two grains, and $T_{i}$ and $T_{j}$ are the temperatures of the $i$-th and $j$-th grains, respectively.
We consider the contacts not only inside the periodic boundary but also on the side boundaries (e.g., the contacts between \#9 and \#10 and between \#21 and \#22 in Fig.\ \ref{fig1}).
The heat conductance at the contact of two grains $H_{\rm c}$ is \citep{Cooper+1969}
\begin{equation}
H_{\rm c} = 2 k_{\rm mat} r_{\rm c},
\end{equation}
where $k_{\rm mat}$ is the material thermal conductivity and $r_{\rm c}$ is the contact radius of monomer grains.
The contact radius $r_{\rm c}$ depends on the monomer radius $R$ and the material parameters \citep{Johnson+1971}.
The heat conductance within a grain $H_{\rm g}$ is also given by \citep{Sakatani+2017}
\begin{equation}
H_{\rm g} = {\left( \frac{4 \pi}{3} \right)}^{1/3} k_{\rm mat} R.
\end{equation}
However, we neglect the effect of $H_{\rm g}$ because $H_{\rm g}$ is sufficiently larger than $H_{\rm c}$ for (sub)micron-sized grains.
Therefore, the temperature structure of the aggregate in the cubic periodic boundary can be calculated by solving Eq.\ (\ref{eq1}) simultaneously for all $N$ grains except lower and upper boundary grains.

\subsection{Thermal conductivity through the solid network}

Once the temperature structure is obtained, we can calculate the total heat flow at the upper boundary $\sum_{\rm upper} F_{i, j}$, where we take the sum of contacts between the upper boundary $i$-th grain and internal $j$-th grain (for the case of Fig.\ \ref{fig1}, $\sum_{\rm upper} F_{i, j} = F_{X, 27} + F_{40, 39}$).
The total heat flow at the upper boundary $\sum_{\rm upper} F_{i, j}$ can be rewritten using the thermal conductivity through the solid network $k_{\rm sol}$ as
\begin{equation}
\sum_{\rm upper} F_{i, j} = k_{\rm sol} \frac{\Delta T}{L} L^{2}.
\end{equation}
In this study, we discuss $k_{\rm sol}$ as a function of the filling factor $\phi$, and rewrite $L$ using $\phi$ as
\begin{equation}
L = {\left( \frac{4 \pi N}{3 \phi} \right)}^{1/3} R.
\end{equation}
Therefore, we obtain $k_{\rm sol}$ as a function of $\phi$ as follows:
\begin{eqnarray}
k_{\rm sol} &=& \frac{1}{L {\Delta T}} {\sum_{\rm upper} F_{i, j}}, \nonumber \\
&=& 2 k_{\rm mat} \frac{r_{\rm c}}{R} \cdot {\left( \frac{3 \phi}{4 \pi N} \right)}^{1/3} {\sum_{\rm upper} \frac{T_{j} - T_{i}}{{\Delta T}}}, \nonumber \\
&\equiv& 2 k_{\rm mat} \frac{r_{\rm c}}{R} \cdot f {(\phi)},
\label{eq7}
\end{eqnarray}
where $f {(\phi)}$ is a dimensionless function of $\phi$.
Note that the total heat flow at the lower boundary $- \sum_{\rm lower} F_{i, j}$ is clearly equal to $\sum_{\rm upper} F_{i, j}$ considering the heat balance.

\section{Results}

Here, we present the dimensionless function $f {(\phi)}$ for nine snapshots from three runs and three densities (Table \ref{table1}).
We calculated $f {(\phi)}$ in three directions for each snapshot and took the arithmetic mean values.
Note that the compressed BCCA aggregates might be isotropic if the number of monomer grains is sufficiently large; thus, we only discuss the mean values of $f {(\phi)}$.

\begin{table}[h]
\caption{Numerical calculation results.}
\label{table1}
\centering
\begin{tabular}{c c c}
\hline\hline
run & $\phi$ & $f (\phi)$ \\
\hline\hline
    $\Box$ & $1.01 \times 10^{-2}$ & $1.44 \times 10^{-4}$ \\
    $\Box$ & $3.02 \times 10^{-2}$ & $1.12 \times 10^{-3}$ \\
    $\Box$ & $9.91 \times 10^{-2}$ & $1.04 \times 10^{-2}$ \\
\hline
    $\triangle$ & $1.00 \times 10^{-2}$ & $3.86 \times 10^{-5}$ \\
    $\triangle$ & $3.00 \times 10^{-2}$ & $5.24 \times 10^{-4}$ \\
    $\triangle$ & $1.00 \times 10^{-1}$ & $6.24 \times 10^{-3}$ \\
\hline
    \scalebox{0.8}{$\bigcirc$} & $9.98 \times 10^{-3}$ & $4.14 \times 10^{-5}$ \\
    \scalebox{0.8}{$\bigcirc$} & $2.99 \times 10^{-2}$ & $5.15 \times 10^{-4}$ \\
    \scalebox{0.8}{$\bigcirc$} & $9.99 \times 10^{-2}$ & $9.05 \times 10^{-3}$ \\
\hline
\end{tabular}
\end{table}

Figure \ref{fig2} shows the dimensionless function $f {(\phi)}$ as a function of the filling factor $\phi$.
The best-fit line given by the least-squares method (green dashed line) is
\begin{equation}
f {(\phi)} = 1.18 \phi^{2.14}.
\end{equation}
Hereafter, we utilize the following more simple relationship between $f {(\phi)}$ and $\phi$ (magenta solid line)
\begin{equation}
f {(\phi)} = \phi^{2}.
\label{eq9}
\end{equation}
\citet{Sakatani+2017} predicted that $f {(\phi)} = {(2 / \pi^{2})} C \phi$, where $C$ is the coordination number.
For highly porous aggregates, the coordination number $C$ is approximately two, and the filling factor dependence on $C$ is weak.
Hence, $f {(\phi)}$ would be proportional to $\phi$ in the model of \citet{Sakatani+2017}; however, in reality, $f {(\phi)}$ is approximately proportional to $\phi^{2}$.

\begin{figure}[h]
\centering
\includegraphics[width=0.9\columnwidth]{fig2.eps}
\caption{
Fitting of the dimensionless function of thermal conductivity $f {(\phi)}$ as a function of the filling factor $\phi$.
The green dashed line is the best-fit line and the magenta solid line represents the simple function $f {(\phi)} = \phi^{2}$.
}
\label{fig2}
\end{figure}

In the context of the thermal conductivity of colloidal nanofluid and nanocomposites, \citet{Evans+2008} revealed that thermal conductivity is strongly affected by the fraction of linear chains that contribute to heat flow in the aggregates.
The contributing grains are called backbone grains, and non-contributing grains are called dead-end grains \citep[\#8, \#32, and \#33 in Fig.\ \ref{fig1};][]{Shih+1990}.
We will discuss the effects of different fractions of backbone and dead-end grains on the thermal conductivity in future research.

By comparing our model to the experimental data of \citet[][]{Krause+2011}, we can confirm the validity of our model (Fig.\ \ref{fig3}).
The magenta line in Fig.\ \ref{fig3} represents the calculated thermal conductivity from Eqs.\ (\ref{eq7}) and (\ref{eq9}), the blue curve is the exponential fitting of experimental data, $k_{\rm sol} = 1.4 e^{7.91 {(\phi - 1)}}\ {\rm W}\ {\rm m}^{-1}\ {\rm K}^{-1}$ \citep{Krause+2011}, and the black dashed line is a model commonly used to study the thermal evolution of planetary bodies, that is, $k_{\rm sol} = \phi k_{\rm mat}$ \citep[e.g.,][]{Sirono2017}.
Both experimental (crosses) and numerical (squares, triangles, and circles) data are plotted.

\begin{figure}[h]
\centering
\includegraphics[width=0.9\columnwidth]{fig3.eps}
\caption{
Experimental (crosses) and numerical (squares, triangles, and circles) data of thermal conductivity through the solid network $k_{\rm sol}$.
Our model (magenta line) was compared to the experimental fitting data of \citet[][blue curve]{Krause+2011}.
}
\label{fig3}
\end{figure}

When we consider the dust aggregates of (sub)micron-sized monomers, the contact radius between monomers $r_{\rm c}$ is given by \citep{Johnson+1971}
\begin{equation}
r_{\rm c} = {\left( \frac{9 \pi \gamma {(1 - \nu^{2})}}{2 Y R} \right)}^{1/3} R,
\label{eq10}
\end{equation}
where $\gamma = 25\ {\rm mJ}\ {\rm m}^{-2}$, $\nu = 0.17$, and $Y = 54\ {\rm GPa}$ are the surface energy, Poisson's ratio, and Young's modulus of ${\rm Si}{\rm O}_{2}$ grains, respectively \citep{Wada+2007}.
We set $R = 0.75\ \mu{\rm m}$ and $k_{\rm mat} = 1.4\ {\rm W}\ {\rm m}^{-1}\ {\rm K}^{-1}$ to the same values as \citet{Krause+2011}.
Figure \ref{fig3} clearly shows that our empirical model is applicable to the $k_{\rm sol}$ of porous aggregates with filling factors of $\phi \sim 0.1$.
Moreover, our model is applicable not only for $\phi \lesssim 0.1$ but also for the range $0.1 \lesssim \phi \lesssim 0.5$.

Note that most of the experimental data of $k_{\rm sol}$ are fitted using exponential functions of $\phi$, which pass the material thermal conductivities \citep[e.g.,][]{Krause+2011,Henke+2013}.
However, when we consider the thermal conductivity through the solid network, i.e., thermal conductivity limited by the necks between two monomers, $k_{\rm sol}$ must be given by $k_{\rm sol} \sim {(r_{\rm c} / R)} k_{\rm mat}$, even for dense dust aggregates whose filling factors are close to unity \citep[e.g.,][]{Chan+1973}.

\section{Discussion}

Finally, we evaluate the total thermal conductivity of porous icy aggregates under vacuum conditions.
Thermal conductivity through the solid network $k_{\rm sol}$ is given by
\begin{equation}
k_{\rm sol} = 2 k_{\rm mat} \frac{r_{\rm c}}{R} \phi^{2},
\end{equation}
and thermal conductivity owing to radiative transfer $k_{\rm rad}$ is given by \citep{Merrill1969}
\begin{equation}
k_{\rm rad} = \frac{16}{3} \sigma_{\rm SB} T^{3} l_{\rm p},
\end{equation}
where $\sigma_{\rm SB} = 5.67 \times 10^{-8}\ {\rm W}\ {\rm m}^{-2}\ {\rm K}^{-4}$ is the Stefan-Boltzmann constant.
We calculate the mean free path of photons in fluffy aggregates of submicron-sized grains $l_{\rm p}$ as follows:
\begin{equation}
l_{\rm p} = \frac{1}{{(\kappa_{\rm abs} + \kappa_{\rm sca})} \rho_{\rm mat} \phi},
\end{equation}
where $\kappa_{\rm abs}$ and $\kappa_{\rm sca}$ are the absorption and scattering mass opacities of monomers, respectively, and $\rho_{\rm mat} = 1.68 \times 10^{3}\ {\rm kg}\ {\rm m}^{-3}$ is the material density.
Here, the composition of icy dust aggregates is consistent with \citet{Pollack+1994}.
The total mass opacity of submicron-sized monomers $\kappa_{\rm abs} + \kappa_{\rm sca}$ is hardly dependent on the wavelength of the thermal radiation $\lambda = 2.9 \times 10^{-3}\ {(T / {\rm K})}^{-1}\ {\rm m}$ for $10^{-6}\ {\rm m} < \lambda < 10^{-4}\ {\rm m}$, and $\kappa_{\rm abs} + \kappa_{\rm sca}$ is on the order of $10^{2}\ {\rm m}^{2}\ {\rm kg}^{-1}$ \citep[e.g.,][]{Kataoka+2014}.
Then we set $l_{\rm p} = 10^{-5} \phi^{-1}\ {\rm m}$ in this study.
Note that, for the case of fluffy aggregates of submicron-sized monomers, the wavelength $\lambda$ is larger than the monomer radius $R$, even if the temperature is on the order of $10^{3}\ {\rm K}$.
Hence, we cannot apply the geometrical optics approximation for the evaluation of $l_{\rm p}$.

Figure \ref{fig4} shows the $k_{\rm sol}$ of crystalline and amorphous icy aggregates, $k_{\rm sol, cr}$ and $k_{\rm sol, am}$, and the thermal conductivity owing to radiative transfer $k_{\rm rad}$ for aggregates composed of icy monomers with a radius of $R = 0.1\ \mu{\rm m}$ and temperature of $T = 40\ {\rm K}$.
We set $\gamma = 100\ {\rm mJ}\ {\rm m}^{-2}$, $\nu = 0.25$, and $Y = 7\ {\rm GPa}$ for icy grains \citep{Wada+2007}.
The material thermal conductivities of crystalline and amorphous grains, $k_{\rm mat, cr}$ and $k_{\rm mat, am}$, are given by $k_{\rm mat, cr} = 5.67 \times 10^{2}\ {(T / {\rm K})}^{-1}\ {\rm W}\ {\rm m}^{-1}\ {\rm K}^{-1}$ and $k_{\rm mat, am} = 7.1 \times 10^{-8}\ {(T / {\rm K})}\ {\rm W}\ {\rm m}^{-1}\ {\rm K}^{-1}$, respectively \citep{Haruyama+1993}.

\begin{figure}[h]
\centering
\includegraphics[width=0.9\columnwidth]{fig4.eps}
\caption{
Comparison between $k_{\rm sol}$ (magenta solid line for crystalline icy aggregates and magenta dashed line for amorphous) and $k_{\rm rad}$ (blue dashed line).
The monomer radius is $R = 0.1\ \mu{\rm m}$ and the temperature is $T = 40\ {\rm K}$.
}
\label{fig4}
\end{figure}

For the case of $\phi < 4 \times 10^{-3}$, the thermal conductivity owing to radiative transfer $k_{\rm rad}$ is larger than the thermal conductivity through the solid network of crystalline icy aggregates $k_{\rm sol, cr}$, even if the temperature is sufficiently low ($T = 40\ {\rm K}$).
If the total thermal conductivity $k_{\rm sol} + k_{\rm rad}$ does not change substantially when crystallization occurs, then the internal temperature of icy planetesimals could still increase after crystallization, which might cause runaway crystallization due to latent heat.
In addition, when the temperature of an icy planetesimal increases, sintering can proceed inside the icy aggregate before monomer grains evaporate or melt \citep[e.g.,][]{Sirono2017}.
We cannot evaluate the contact radius $r_{\rm c}$ from Eq.\ (\ref{eq10}) when aggregates are sintered, and $k_{\rm sol}$ increases linearly as a consequence of the increase of $r_{\rm c}$.
Sintering might also affect the mechanical strength of aggregates \citep[e.g.,][]{Omura+2017} and the critical velocity for collisional growth \citep[e.g.,][]{Sirono+2017}.
Therefore, the growth pathways of icy planetesimals might be altered by sintering of icy aggregates.
Not only icy planetesimals but also rocky aggregates could experience sintering before growing into dm-sized bodies, which might explain the retainment of chondrules inside fluffy aggregates \citep{Arakawa2017}.
Note also that the total thermal conductivity might be controlled by the thermal conductivity due to gas diffusion for the case of fluffy aggregates in high gas density environments \citep[e.g., the innermost region of protoplanetary disks and/or planetary surfaces;][]{Piqueux+2009}.

In conclusion, we have revealed the filling factor dependence of the thermal conductivity of porous aggregates.
We showed that the thermal conductivity of highly porous aggregates is significantly lower than previously assumed.
In future work, we will reexamine the growth pathways of planetesimals in protoplanetary disks, and combine this with a density and thermal evolution analysis.

\begin{acknowledgements}
We thank the referee Gerhard Wurm for constructive comments.
This work is supported by JSPS KAKENHI Grant (15K05266).
S.A.\ is supported by the Grant-in-Aid for JSPS Research Fellow (17J06861).
\end{acknowledgements}

\bibliographystyle{aa}