\section{Introduction}

It is unequivocal that quantum key distribution (QKD)~\cite{QKDreview1,QKDreview2} needs to bridge the gap between theory and practice. In theory,
QKD offers perfect security. In practice,
however,
it does not, as most practical devices behave differently from the theoretical models assumed in the security proofs.
As a result, we face
implementation loopholes, or so-called side-channels, which may be used by adversaries without being detected,
as seen
in recent attacks against certain commercial QKD
systems~\cite{sidechannelrefs2a, sidechannelrefs1a, sidechannelrefs2b, sidechannelrefs2c, sidechannelrefs1b, sidechannelrefs1c, sidechannelrefs2d, sidechannelrefs2f, sidechannelrefs2e}.

There are two potential ways to
guarantee security in the
realisations of QKD. The first is to
develop mathematical models that perfectly match the behaviour of physical apparatuses,
and then incorporate this information into a new security proof.
While this is plausible in theory, unfortunately it is very hard to realise in practice, if not impossible.
The second alternative is to design new protocols and develop security proof techniques
that are
compatible with a wide class of device
imperfections.
This allows us to omit an
accurate characterisation of real apparatuses.
The most well-known example of such a solution is
(full) device-independent QKD (diQKD)~\cite{diqkdrefsa,diqkdrefsb,diqkdrefsc,diqkdrefsd,diqkdrefse}. Here, the legitimate users of the system
(typically called Alice and Bob) treat their devices as two quasi ``black boxes", i.e.,
they need to know which elements their boxes contain, but not
how they fully function~\cite{diqkdrefs2}.
The security of diQKD relies
on the violation of a Bell inequality~\cite{bell1,bell2}, which certifies the presence of quantum correlations.
Despite its beauty, however,
this approach is highly impractical because it
requires a loophole-free Bell test which at the moment is still unavailable~\cite{Pearle1970}.
Also, its secret key rate at practical distances is very limited~\cite{Gisin2010a,Gisin2010b}.

Very recently,
a novel approach has been
introduced, which is fully practical and feasible to implement.
This scheme is known as
measurement-device-independent QKD (mdiQKD)~\cite{mdiQKD} and
offers a clear avenue to bridge
the gap between theory and practice.
Its feasibility has been
promptly demonstrated both in laboratories and
via field-tests~\cite{Rubenok2012a,Rubenok2012b,Rubenok2012c,Rubenok2013a,Zhiyuan:2013}.
It successfully
removes
all (existing and yet to be discovered) detector side-channels~\cite{sidechannelrefs2a,sidechannelrefs2b,sidechannelrefs2c,sidechannelrefs2d,sidechannelrefs2e,sidechannelrefs2f},
which, arguably, is the most critical part of most QKD implementations.
Importantly,
in contrast to diQKD this solution does not require that Alice and Bob perform a loophole-free Bell test;
it is enough if they prove the presence of entanglement in a quantum state that is effectively distributed between them,
just like in standard QKD schemes~\cite{curty_prl2004}.
In addition, now Alice and Bob may treat the measurement apparatus
as a true ``black box", which may be fully controlled
by the adversary. A slight drawback is that
Alice and Bob need to characterise the quantum states
[e.g., the polarisation degrees of freedom of phase-randomised weak coherent pulses (WCPs)]
that they send through
the channel. But, as this process can be verified in a protected environment outside the influence of the
adversary, it is less likely to be a problem.
For completeness, the readers can refer to~\cite{charles2} where a characterisation of the prepared states is no longer required.

Nevertheless, so far the security of mdiQKD has only been proven
in the asymptotic regime~\cite{mdiQKD}, which assumes that Alice and Bob have access to an unlimited amount of resources,
or in the finite regime but only against particular types of attacks~\cite{finitea,finiteb}.
In summary, until now, a rigorous security proof of mdiQKD that takes full
account of the finite size effects~\cite{finite1,finite2,finite3} has appeared to be missing and, for this reason, the
feasibility of long-distance implementations of mdiQKD within a reasonable time-frame of signal transmission
has remained undemonstrated.

The main contributions of this work are twofold. First, in contrast to existing heuristic results on
mdiQKD, we provide, for the first time, a security proof in the finite-key regime
that is valid against general attacks, and
satisfies the composability definition~\cite{RennerThesis,MullerQuadeRenner}
of QKD. Second, we apply large deviation
theory, specifically a multiplicative form of the Chernoff bound~\cite{cher}, to perform the parameter estimation step.
The latter is crucial to demonstrate that a long-distance
implementation of
mdiQKD
(e.g., $150$ km of optical
fiber with $0.2$ dB/km) is feasible
within a reasonable time-frame. To obtain high secret key rates in this scenario,
it is common to use decoy state techniques~\cite{decoy1,decoy2,decoy3},
both for standard QKD protocols and mdiQKD.
Here a key challenge is to estimate
the transmittance
and the quantum bit error rate (QBER) of the
single-photon component of the signal at the presence of
high losses (e.g., $30$ dB). We show that such an estimation
problem can be solved
using the Chernoff bound, as it provides good bounds for the parameters above
even in the high-loss regime. We highlight that our results
can be applied to other QKD protocols
(e.g., the standard decoy state BB84 protocol~\cite{decoy1,decoy2,decoy3}) as well as to general
experiments in quantum information.

\section{Security Definition}

Prior to stating the protocol, let us quickly review the security framework~\cite{RennerThesis,MullerQuadeRenner}
 that we are considering here.
A general QKD protocol (executed by Alice and Bob) generates either a pair of bit strings $S_{\rm A}$ and $S_{\rm B}$,
or a symbol $\perp$ to indicate the abort of the protocol. In general, the string of Alice, $S_{\rm A}$, can be quantum mechanically
correlated with a quantum state that is held by the adversary. Mathematically, this situation is described by the classical-quantum state
\begin{eqnarray}
\rho_{\rm AE}=\sum_{s}\proj{s}\otimes \rho_{\rm E}^s, \nonumber
\end{eqnarray}
where $\{\ket{s}\}_s$ denotes an orthonormal basis for Alice's system, and the
subscript ${\rm E}$ indicates the system of the adversary.

Ideally, we say that a QKD protocol is secure if it satisfies two conditions, namely the correctness and the secrecy.
The correctness condition is met if $S_{\rm A}=S_{\rm B}$, i.e., Alice's and Bob's bit strings are identical.
The secrecy condition is met if $\rho_{\rm AE}=U_{\rm A} \otimes \rho_{\rm E}$, where $U_{\rm A}=\sum_s\frac{1}{|{\mathcal S}|}\proj{s}$
is the uniform mixture of all possible values of the bit string $S_{\rm A}$. That is, the system of the adversary is completely decoupled from
that of Alice.

Owing to the presence of errors, however, these two conditions can never be perfectly met.
For example, in the finite-key regime it is impossible to guarantee $S_{\rm A}=S_{\rm B}$ with certainty.
In practice, this implies that we need to allow for some minuscule errors. That is, we say that a QKD scheme is $\epsilon_{\textrm{cor}}$-correct
if $\Pr[S_{\rm A} \not= S_{\rm B}] \leq \epsilon_{\textrm{cor}}$, i.e.,
the probability that Alice's and Bob's bit strings are not identical is not greater than $\epsilon_{\textrm{cor}}$.
Similarly, we say that a protocol is $\epsilon_{\textrm{sec}}$-secret if
\begin{eqnarray}
\frac{1}{2}\|\rho_{\rm AE}-U_{\rm A}\otimes\rho_{\rm E}\|_1 \leq \epsilon_{\textrm{sec}}, \nonumber
\end{eqnarray}
where $\|\cdot\|_1$ denotes the trace norm.
That is, the state $\rho_{\rm AE}$ is $\epsilon_{\textrm{sec}}$-close to the ideal situation described by $U_{\rm A} \otimes \rho_{\rm E}$.
Thereby a QKD protocol is said to be $\epsilon$-secure if it is both $\epsilon_{\textrm{cor}}$-correct and $\epsilon_{\textrm{sec}}$-secret,
with $\epsilon_{\textrm{cor}}+\epsilon_{\textrm{sec}}\leq\epsilon$.

With this security definition we are able to guarantee that the security of the protocol holds even when combined with other protocols, i.e.,
the protocol is secure in the so-called universally composable framework~\cite{RennerThesis,MullerQuadeRenner}.

\section{Protocol definition}\label{sec_prot}

The setup is illustrated in Fig.~\ref{scheme}.
Alice and Bob use a laser source to generate quantum signals that are diagonal in the Fock basis.
Instances of such sources include attenuated laser diodes emitting phase-randomised WCPs, triggered spontaneous parametric down-conversion sources,
and practical single-photon sources.
Each pulse is prepared in a different BB84 state~\cite{bb84}, which is
selected, for example, uniformly at random from two mutually unbiased bases,
denoted as $\rm Z$ and $\rm X$.
The signals are then sent to an untrusted relay Charles, who is supposed to perform
a Bell state measurement that projects them into a Bell state.
Also, Alice and Bob apply decoy state techniques~\cite{decoy1,decoy2,decoy3}
to estimate the gain (i.e., the probability that the relay outputs a successful result)
and the QBER for various input photon-numbers.
\begin{figure}[h]
\begin{center}
 \includegraphics[scale=0.48]{Setup.eps}
 \end{center}
 \caption{A schematic diagram of mdiQKD.
Alice and Bob prepare quantum signals in different BB84 polarisation states~\cite{bb84}
with a polarisation modulator (Pol-M). Also, they use
an intensity modulator (Decoy-IM) to
generate decoy states.
The signals are sent to an untrusted relay Charles, who is supposed to perform
a Bell state measurement that projects the incoming signals into a Bell state.
See the main text for details.
\label{scheme}}
\end{figure}

Next, Charles announces whether or not his measurements are successful, including the Bell states obtained.
Alice and Bob keep the data that correspond to these instances and discard the rest. Also, they
post-select the events where they employ the same basis. Finally, either Alice or Bob
flips part of her/his bits to correctly correlate them with those of the other. See
Table~\ref{table_protocol} for a detailed description of the different steps of the protocol.
\begin{table}
\begin{description}
\item[1. State Preparation]
Alice and Bob repeat the first four steps of the protocol
for $i=1,\ldots,N$ till the conditions in the
Sifting step are met. 	For each $i$, Alice chooses an intensity $a\in\{a_{\rm s}, a_{\rm d_1},\ldots,a_{{\rm d}_n}\}$,
a basis $\alpha\in\{{\rm Z, X}\}$, and a random bit $r\in\{0,1\}$ with probability
$p_{a,\alpha}/2$. Here $a_{\rm s}$ ($a_{{\rm d}_j}$) is the intensity of the signal (decoy) states.
Next, she generates a
quantum signal (e.g., a phase-randomised WCP)
of intensity $a$ prepared in the basis state of $\alpha$ given
by $r$. Likewise, Bob does the same.

\item[2. Distribution] Alice and Bob send their states to Charles via the quantum channel.

\item[3. Measurement]
If Charles is honest, he measures the signals received with a Bell state measurement.
In any case, he informs Alice and Bob (via a public channel) of whether or not his measurement was successful. If successful,
he reveals the Bell state obtained.

\item[4. Sifting]
If Charles reports a successful result, Alice and Bob broadcast (via an authenticated channel)
their intensity and basis settings.
For each Bell state $k$, we define two groups of sets:
${\mathcal Z}^{a,b}_k$ and ${\mathcal X}^{a,b}_k$. The first (second) one identifies
signals
 where Charles declared the Bell state $k$ and Alice and Bob selected the
intensities $a$ and $b$ and the basis $\rm Z$
($\rm X$).
The protocol repeats these steps until
$|{\mathcal Z}^{a,b}_k|\geq N^{a,b}_k$ and $|{\mathcal X}^{a,b}_k|\geq M^{a,b}_k$
$\forall a, b, k$.
Next, say Bob flips part of his bits to correctly correlate them with those of Alice
(see Table~\ref{table1}).
Afterwards, they execute the last steps of the protocol for each $k$.

\item[5. Parameter Estimation]
Alice and Bob use $n_k$ random
bits from
${\mathcal Z}^{a_{\rm s},b_{\rm s}}_k$
to form the code bit strings
$Z_k$ and $Z'_k$, respectively. The
remaining
$R_k$ bits from ${\mathcal Z}^{a_{\rm s},b_{\rm s}}_k$ are used to compute
the error rate
$E^{a_{\rm s},b_{\rm s}}_k=\frac{1}{R_k}\sum_l r_l\oplus r'_l$,
where $r'_l$ are Bob's bits.
If $E^{a_{\rm s},b_{\rm s}}_k>E_{\rm tol}$, Alice and Bob assign an empty string to
$S_k$ and abort steps $6$ and $7$ for this
$k$.
The protocol only aborts if $E^{a_{\rm s},b_{\rm s}}_k>E_{\rm tol}$ $\forall k$.
If $E^{a_{\rm s},b_{\rm s}}_k\leq{}E_{\rm tol}$, Alice and Bob
use ${\mathcal Z}^{a,b}_k$ and ${\mathcal X}^{a,b}_k$ to
estimate $n_{k,0}$, $n_{k,1}$ and $e_{k,1}$. The parameter $n_{k,0}$ ($n_{k,1}$)
is a lower bound for the number of bits
in $Z_k$ where Alice (Alice and Bob) sent a vacuum (single-photon) state.
$e_{k,1}$ is
an upper bound for the single-photon phase error rate.
If $e_{k,1}>e_{{\rm tol}}$,
an empty string is assigned to $S_k$ and
steps $6$ and $7$ are aborted for this $k$, and
the protocol only aborts if $e_{k,1}>e_{{\rm tol}}$ $\forall k$.

\item[6. Error Correction] For those $k$ that passed the parameters estimation step,
Bob obtains an estimate ${\hat Z}_k$
of $Z_k$ using
an information reconciliation scheme. For this, Alice sends him
${\rm leak}_{{\rm EC},k}$ bits of error correction data.
Next, Alice computes a hash of $Z_k$ of length $\lceil\log_2(4/\epsilon_{\rm cor})\rceil$ using
a random universal$_2$
hash function, which she sends to Bob together with the hash~\cite{RennerThesis}.
If ${\rm hash}({\hat Z}_k)\neq {\rm hash}(Z_k)$, Alice and Bob assign an empty string to $S_k$ and abort step $7$ for this $k$.
The protocol only aborts if
${\rm hash}({\hat Z}_k)\neq {\rm hash}(Z_k)$
$\forall k$.

\item[7. Privacy Amplification] If $k$ passed the error correction step,
Alice and Bob apply a random universal$_2$ hash function
to $Z_k$ and ${\hat Z}_k$ to extract two shorter strings of length $\ell_k$~\cite{RennerThesis}.
Alice obtains $S_k$ and Bob
${\hat S}_k$.
The concatenation of
$S_k$ (${\hat S}_k$) form the secret key
$S_{\rm A}$ ($S_{\rm B}$).

\end{description}
\caption{Protocol Definition.}
\label{table_protocol}
\end{table}
\begin{table}[htb]
  \centering
  \begin{tabular}{ccccc}
  \quad & \multicolumn{4}{c}{Bell state reported by Charles}\\
    \hline\hline
  Alice \& Bob & $\quad \ket{\psi^-}$ \quad & \quad $\ket{\psi^+}$ \quad & \quad $\ket{\phi^-}$ \quad & \ \quad $\ket{\phi^+}$ \ \quad \\
  \hline
  {\rm {\rm Z} basis} & \quad {\rm Bit flip}  \quad & \quad {\rm Bit flip} \quad  & \quad - \quad  & \ \quad\quad - \ \quad\quad  \\
  {\rm {\rm X} basis} & \quad {\rm Bit flip} \quad & \quad - \quad & \quad {\rm Bit flip} \quad &\ \quad\quad - \ \quad\quad \\
  \hline\hline
\end{tabular}
\caption{Post-processing of data in the sifting step. To guarantee that their bit strings are correctly correlated,
say Bob applies a bit flip to part of his data, depending on the Bell state reported by Charles
and the basis setting selected.
  }\label{table1}
\end{table}

Since Charles' measurement is basically used to post-select entanglement
between Alice and Bob, the security of
mdiQKD can be proven by using the idea of time reversal. Indeed, mdiQKD builds on the earlier proposals of time-reversed EPR protocols
by Biham et al.~\cite{biham} and Inamori~\cite{inamori}, and combine them with the decoy state technique.
The end result is the best of both worlds---high performance and high security. We note on passing that the idea of
time reversal has also been previously used in other quantum information protocols including one-way quantum computation.

\section{Security analysis}\label{sec_proof}

We now present one main result of our paper. It states that the protocol introduced above is both
$\epsilon_{\textrm{cor}}$-correct and $\epsilon_{\textrm{sec}}$-secret, given that the length $\ell$ of the secret key
$S_{\rm A}$
is selected appropriately for a given set of observed values. See
Table~\ref{table_protocol}
for the definition of the different parameters that we consider in this section.

The correctness of the protocol is guaranteed by its error correction step, where,
for each possible Bell state $k$,
Alice sends
a hash of $Z_k$ to Bob, who compares it with the hash of ${\hat Z}_k$. If both
hash values are equal, the protocol gives $S_k={\hat S}_k$ except with error
probability $\epsilon_{\textrm{cor}}/4$. If ${\rm hash}({\hat Z}_k)\neq {\rm hash}(Z_k)$, it outputs the empty string
(i.e., the protocol is trivially correct). Moreover, if the protocol aborts, the result is
$\perp$, i.e.,
it is also correct. This guarantees that $S_{\rm A}=S_{\rm B}$ except with
error probability less or equal than $\epsilon_{\textrm{cor}}$.
Alternatively to this method, Alice and Bob may also guarantee the
correctness of the protocol by exploiting
properties of the error correcting code employed \cite{lut}.

If the length $\ell_k$ of each bit string $S_k$, which forms the secret key $S_{\rm A}$, satisfies
\begin{equation}\label{klength}
\ell_k \leq n_{k,0}+n_{k,1}\left[1-h\left(e_{k,1}\right)\right]-{\rm leak}_{{\rm EC},k}-\log_2\frac{8}{\epsilon_{\rm cor}}
-2\log_2\frac{2}{\varepsilon'_k\hat{\varepsilon}_k}-2\log_2\frac{1}{2\varepsilon_{k,\tn{PA}}},
\end{equation}
the protocol is $\epsilon_{\rm sec}$-secret, with
$\epsilon_{\rm sec}=\sum_k \epsilon_{k,{\rm sec}}$ and $\epsilon_{k,{\rm sec}}=
2(\varepsilon'_k+2\varepsilon_{k,e}+\hat{\varepsilon}_k)+\varepsilon_{k,\tn{b}}+\varepsilon_{k,0}
+\varepsilon_{k,1}+\varepsilon_{k,\tn{PA}}$. In Eq.~(\ref{klength}), $h(x)=-x\log_2{(x)}-(1-x)\log_2{(1-x)}$ is the binary Shannon entropy, and
the parameters $\varepsilon_{k,0} $, $\varepsilon_{k,1}$, and $\varepsilon_{k,e}$ quantify, respectively,
the probability that the estimation of the terms $n_{k,0}$,  $n_{k,1}$ and $e_{k,1}$ is
incorrect. A sketch of the proof of Eq.~(\ref{klength}) can be found in Appendix~\ref{secre}.
There it is also explained the meaning of all the epsilons contained in the term
$\epsilon_{k,{\rm sec}}$, which we omit here for simplicity.
In the asymptotic limit of very large data blocks,
the terms reducing the length of $S_{\rm A}$
due to statistical fluctuations may be neglected, and thus $\ell$
  satisfies $\ell\leq\sum_k\max{\{n_{k,0}+n_{k,1}\left[1-h\left(e_{k,1}\right)\right]-{\rm leak}_{{\rm EC},k},0\}}$,
as previously obtained in \cite{mdiQKD}. That is,
$n_{k,0}$ and $n_{k,1}$ provide a positive contribution
to the secret key rate, while $n_{k,1}h\left(e_{k,1}\right)$ and
${\rm leak}_{{\rm EC},k}$ reduce it.
The term $n_{k,1}h\left(e_{k,1}\right)$ corresponds to the information
removed from $Z_k$ in the privacy amplification step of the protocol, while ${\rm leak}_{{\rm EC},k}$ is the
information revealed by Alice in the error correction step.

The second main contribution of this work is an estimation method to obtain the
relevant parameters $n_{k,0}$,  $n_{k,1}$ and $e_{k,1}$ needed to evaluate the key rate formula above,
when Alice and Bob send Charles a finite number $N$ of signals and use a finite number of decoy states.
We solve this problem using techniques
in large deviation theory.
More specifically, we employ
the Chernoff bound~\cite{cher}.
It is important to note that standard techniques such as Azuma's
inequality~\cite{azuma} do not give very good bounds here. This is because this result does not consider the properties of the
a priori distribution. Therefore, it is far from optimal for situations such as high loss or a highly bias coin flip,
which are relevant in long-distance QKD. In contrast,
the Chernoff bound takes advantage of the property of
the distribution and provides good bounds
even in a high-loss regime.

More precisely,
we show that the estimation of $n_{k,0}$,
$n_{k,1}$ and $e_{k,1}$
can be
formulated as
a linear program,
which
 can be solved efficiently
in polynomial time and gives the exact optimum even for large dimensions \cite{lpr}.
Importantly, this general method is valid for any finite number of decoy states used by Alice and Bob, and
for any photon-number distribution of their signals. Also, for the typical
scenario where
Alice and Bob send phase-randomised WCPs together with two decoy states each, we
solve analytically the linear program, and obtain analytical
expressions for the parameters above, which can be used directly in
an experiment. A sketch of the estimation method is given in Appendix~\ref{sketch}. For a detailed
analysis of both estimation techniques we refer to the Appendices~\ref{ap_anal} and~\ref{ap_numer}.

\section{Discussion}

In this section we analyse the behaviour of the secret key rate
provided in Eq.~(\ref{klength}).
In our simulation,
we consider that Alice and Bob encode their bits in the
polarisation degrees of freedom of phase-randomised WCPs.
Also, we assume
that Charles uses the linear optics quantum relay
illustrated in Fig.~\ref{fig_relay},
which is able to
identify
two of the four Bell states.
With this setup, a successful Bell state measurement
corresponds to the observation of precisely two detectors
(associated to orthogonal polarisations) being triggered.
Note, however, that the results presented in this paper can be
applied to other
types of coding schemes like, for instance,
phase or time-bin coding~\cite{QKDreview1,QKDreview2},
and to
any quantum operation that Charles may perform,
as they solely depend on the measurement results that he announces.
\begin{figure}[t]
\centerline{\includegraphics[width=0.24\textwidth]{Setup_Charles.eps}}
\vspace{-0.2cm}
\caption{A schematic diagram of Charles' measurement device.
The signals from Alice and Bob interfere at a $50:50$ beam splitter (BS),
which has on each end a polarising beam splitter (PBS) that projects
the incoming photons into either horizontal ($H$) or
vertical ($V$) polarisation states. A click in the single-photon
detectors $D_{1H}$ and
$D_{2V}$, or in $D_{1V}$ and $D_{2H}$, indicates a projection
into the Bell state $\ket{\psi^{-}}=1/\sqrt{2}(\ket{HV}-\ket{VH})$,
while a click in $D_{1H}$ and $D_{1V}$, or in
$D_{2H}$ and $D_{2V}$, implies a projection into the Bell state
$\ket{\psi^{+}}=1/\sqrt{2}(\ket{HV}+\ket{VH})$.
}\label{fig_relay}
\end{figure}

We use experimental parameters from~\cite{param}.
But, whereas~\cite{param} considers a free-space channel,
we assume
a fiber-based channel with a loss of $0.2$ dB/km.
The detection efficiency of the
relay (i.e., the transmittance of its optical components together with the efficiency of its
detectors)
is $14.5\%$, and the
background count rate is $6.02\times{}10^{-6}$.
Moreover, we use a rather generic channel model
that includes an intrinsic error rate which
simulates the misalignment and instability of the optical system.
This is done by placing a unitary rotation in both input arms of the $50:50$ beam splitter,
and another unitary rotation in one of its output arms \cite{feihu}.
In addition,
we
fix the security bound to
$\epsilon=10^{-10}$.
\begin{figure}[t]
\centerline{\includegraphics[width=0.44\textwidth]{Rate_vs_distance_wasymptotic.eps}}
\vspace{-0.2cm}
\caption{Expected key rate as function of the distance. Secret key rate $\ell/N$ in logarithmic scale for the
protocol introduced in section~\ref{sec_prot} with phase-randomised WCPs as a
function of the distance. The solid lines correspond to different values for the total number of signals $N$
sent by Alice and Bob. The overall misalignment
in the channel is $1.5\%$,
and the security bound $\epsilon=10^{-10}$.
For simulation purposes we consider the following experimental
parameters~\cite{param}: the loss coefficient of the channel is
$0.2$ dB/km, the detection efficiency of the relay is $14.5\%$,
and the
background count rate is $6.02\times{}10^{-6}$.
Our results show clearly that even with a realistic finite size of data,
say $N=10^{12}$ to $10^{14}$, it is possible to
achieve secure mdiQKD at long distances.
In comparison, the dotted line represents a lower bound on the
secret key rate for the asymptotic case where Alice and Bob send
Charles
infinite signals and use an infinite number of decoy settings.
}\label{fig1}
\end{figure}
\begin{figure}[t]
\centerline{\includegraphics[width=0.44\textwidth]{figc.eps}}
\vspace{-0.2cm}
\caption{Expected key rate as function of the block size. The plot shows the secret key rate $\ell/N$ in logarithmic scale
as a function of the total number of signals $N$
sent by Alice and Bob in the limit of zero distance. The security bound $\epsilon=10^{-10}$.
The solid lines correspond to different values for the intrinsic error
rate due to the misalignment and instability of the optical system.
The horizontal dotted lines show the asymptotic rates.
The experimental parameters are the ones described in the
caption of Fig.~\ref{fig1}.
Our results show that, even for a finite size of signals sent
by Alice and Bob, mdiQKD is robust to intrinsic errors due to basis
misalignment and instability of the optical system.}\label{fig2}
\end{figure}

The results are shown in Figs.~\ref{fig1}-\ref{fig2} for the situation where Alice and Bob use two decoy states each.
In this scenario,
we obtain
the
parameters $n_{k,0}$, $n_{k,1}$ and $e_{k,1}$ using the analytical estimation procedure introduced above
(see Appendix~\ref{ap_anal} for more details).
The first figure illustrates the secret key rate (per pulse) $\ell/N$ as a function of the distance
between Alice and Bob
for different values of the total number of signals $N$ sent.
We fix $\epsilon_{\rm cor}=10^{-15}$; this corresponds to a realistic
hash tag size in practice \cite{RennerThesis}.
Also, we fix
the intensity of the weakest
decoy states to $a_{\rm d_2}=b_{\rm d_2}=5\times10^{-4}$, since, in practice, is
difficult to generate a vacuum state due to imperfect extinction. This value for $a_{\rm d_2}$ and $b_{\rm d_2}$ can be easily achieved with
a standard intensity modulator.
Moreover, for simplicity, we assume an
error correction leakage that is a fixed fraction of the sifted key length $n_k$, i.e.,
${\rm leak}_{{\rm EC},k}=n_k\zeta{}h(E_k^{a_{\rm s},b_{\rm s}})$,
with $\zeta=1.16$ and
where
$h(\cdot)$
is again the binary Shannon entropy \cite{finite1}.
In a realistic scenario, however, the value of $\zeta$ typically depends on the
value of $n_k$, and when $n_k< 10^5$ the parameter
$\zeta$ may be bigger than $1.16$.
For a given distance, we optimise numerically $\ell/N$ over all the
free parameters of the protocol. This includes the intensities
$a_{\rm s}, a_{\rm d_1}, b_{\rm s}$ and $b_{\rm d_1}$,
the probability
distributions $p_{a,\alpha}$ and $p_{b,\beta}$
in the state preparation step, the parameters
$N_k^{a,b}$ and $M_k^{a,b}$ in the sifting step, the term
$n_k$ in the parameter estimation step,
and the different epsilons contained in $\epsilon_{\rm sec}$.
Our simulation result shows clearly that mdiQKD is feasible
with current technology and does not
require high efficiency detectors for its implementation.
If Alice and Bob use laser diodes operating at $1$ GHz repetition rate,
and each of them sends $N=10^{13}$ signals,
we find, for instance, that they can distribute a
$1$ Mb secret key over a $75$ km fiber link in less than $3$ hours.
This scenario corresponds to the
red line shown in Fig.~\ref{fig1}.
Notice that, at telecom wavelengths, standard InGaAs detectors have
modest detection efficiency of about $15\%$. Since mdiQKD requires
two-fold coincidence rather than single-detection events, as is the case in the
standard decoy state protocol, the key rate of mdiQKD is
lower than that of the standard decoy state scheme. However,
with high-efficiency detectors such as silicon detectors~\cite{det1}
in $800$nm or high-efficiency SSPDs~\cite{det2}, the key rate of
mdiQKD can be made comparable to that of the standard decoy state protocol.

The second figure illustrates $\ell/N$ as a function of $N$ for different
values of the misalignment in the limit of zero distance.
For comparison, this figure also includes the asymptotic secret key rate
when Alice and Bob send an infinite number of signals and
use an infinite number of decoy states \cite{mdiQKD}.
Our results show that significant secret key rates are already possible with
$10^{11}$ signals, given that the error rate is not too large.

\section{Conclusion}

We have proved the security of
mdiQKD in the finite-key regime against general attacks.
This is the only known
fully practical QKD protocol that offers an avenue to bridge the
gap between theory and practice in QKD implementations.
Importantly, our results clearly demonstrate that even with practical
signals [e.g., phase-randomised weak coherent pulses (WCPs)] and
a finite size of data (say $10^{12}$ to $10^{14}$ signals) it is possible to
perform secure mdiQKD over long distances (up to about $150$ km).

To achieve high secret key
rates in such high-loss regime, it is typical
both for standard QKD schemes and
mdiQKD
to use decoy state techniques. A main challenge in this scenario is to obtain tight bounds
for the gain and quantum bit error rate (QBER) of the single-photon components
sent by Alice and Bob. We have shown that this estimation problem can be successfully
solved using techniques in large deviation theory, more precisely,
the Chernoff bound. This result takes advantage of the property of the
distribution, and thus provides good bounds for the relevant
parameters even in the presence of high losses, as is the case in QKD realisations.

Using the Chernoff bound we have rewritten the problem of
estimating the gain and QBER of the single-photon signals
as a linear program,
which
 can be solved efficiently
in polynomial time. This general method is valid for any finite number of decoy states, and
for any photon-number distribution of the signals. It can be used,
for instance, with
laser diodes emitting phase-randomised WCPs, triggered spontaneous parametric
down-conversion sources,
and practical single-photon sources.
Also, for the common
scenario where
Alice and Bob send phase-randomised WCPs together with two decoy states each, we
have obtained tight analytical bounds for the quantities above. These results apply
to different types of coding schemes like, for example,
polarisation, phase or time-bin coding.

\section{Acknowledgements}

We thank Xiongfeng Ma and Johan L\"ofberg for valuable comments
and stimulating discussions, and
Lina M. Eriksson
for comments on the writing and presentation of the paper.
F. Xu. thanks the Paul Biringer Graduate Scholarship for financial support.
We
acknowledge support from the European Regional Development Fund (ERDF),
the Galician Regional Government (projects CN2012/279 and CN 2012/260,
``Consolidation of Research Units: AtlantTIC"),
NSERC, the CRC program,
the
National Centre of Competence in Research QSIT,
the Swiss NanoTera project QCRYPT, the FP7 Marie-Curie IAAP
QCERT project, and CHIST-ERA project DIQIP.
K.T. acknowledges support from the project ``Secure photonic network technology''
as part of ``The project UQCC'' by the National
Institute of Information and Communications Technology
(NICT) of Japan, and from the Japan Society for
the Promotion of Science (JSPS) through its Funding
Program for World-Leading Innovative R\&D on Science
and Technology (FIRST Program).


\section{Introduction}

Neutron stars can be used as probes to gain insight into the structure of
matter at high densities ($10^{15} g/cm^{3}$ at their core). Knowledge on mass, radius, and any other global
property of the neutron star can give us important information on the state
and composition of matter at those extreme physical conditions. But, although the
mass
can be determined with quite good precision, a direct measurement of the
radius is
currently very difficult using astrophysical observations.

Nevertheless, the detection of gravitational waves originated in the neutron
stars could be used to obtain estimations, or at least constrains, on the
value of the gravitational radius of the star. With knowledge on the mass and
radius, we would be able to have a better understanding of the behavior of
matter at high density and pressure. So the study of these
gravitational waves are crucial in the understanding of the physics of matter at
extreme densities.

Neutron stars present a layer structure which can be essentially described by
two very different regions:
a core (where exotic matter could be found) enveloped by a \textit{crust}
(with a solid
crystalline
structure similar to a metal).
The relativistic fluid that compose the neutron star can be described as
a perfect fluid, for which the equation of state is needed. Because not so
much information is available beyond the
nuclear densities, the EOS is very model dependent, since different
particle populations may appear at those energies
\cite{haensel2007neutron,glendenning2000compact}. Another interesting feature
is that along the core of the star
and specially in the core-\textit{crust} interface,
first order phase transitions could be found in
realistic equations of state. These phase transitions result in small
discontinuities in
the energy density of the star matter
\cite{haensel2007neutron,heiselberg2000}.

Very interesting information can be obtained from the recent measurement of
the mass of PSR J1614-2230, of $1.97 M_\odot$, which impose constraints on
equations of
state for exotic matter that could be found on neutron stars
\cite{Lattimer_Prakash}. Several
equations of state that are able to produce neutron star configurations with
$2 M_\odot$ and exotic matter in their cores have been proposed
\cite{Bednarek,Sedrakian,Weissen1,Weissen2}. If other global parameters could
be measured, more important constraints could be imposed to the equation of
state.

If an internal or external event perturbate a neutron star, it may oscillate
non radially, emitting gravitational waves to the space
\cite{GW_Living_Review2009}. Large-scale interferometric gravitational wave detectors LIGO, GEO, TAMA and VIRGO, have reached the
original design sensitivity, and the detection of the first gravitational waves could
happen in the next years \cite{GW_Living_Review2011}. These detections will provide useful information about the composition of the astrophysical objects that generate them, like neutron stars  \cite{GW_Living_Review2009}.

The gravitational waves emitted by a neutron star present
dominant frequencies which can be studied using the quasi-normal mode formalism
\cite{Kokkotas_Schmidt1999,Nollert1999,Rezzolla2003ua}. These eigen-frequencies
are given by a
complex number. The real part gives us the oscillation frequency of the mode,
while the
imaginary part gives us the inverse of the damping time. The quasi-normal mode
spectrum is quite dependent on the properties of the star, i.e. the equation
of state.

Hence, the detection of gravitational radiation from neutron stars combined
with a complete theoretical study of the possible spectrums that can be
obtained with different equations of
state can be very
useful in the
determination of the matter behavior beyond nuclear matter
\cite{Benhar_Berti_Ferrari_1999,Kokkotas2001,Benhar_Ferrari_Gualtieri_2004,Benhar2005,Ferrari2007,Chatterjee2009,Wen2009}.

The formalism necessary to study quasi-normal modes were first developed for black holes by Regge and Wheeler \cite{ReggeI} and by Zerilli
\cite{PhysRevLett.24.737}. Quasi-normal modes can be differentiated into polar
and axial modes. The equation
describing the quasi-normal mode perturbation of the Schwarzschild metric is
essentially a Schrodinger like equation: the Regge-Wheeler equation for axial
perturbations and the Zerilli equation for polar ones. For black holes, both
types of modes are space-time modes. The formalism was
studied in the context of neutron stars first by Thorne
\cite{ThorneI,ThorneII,ThorneIII,ThorneIV,ThorneV}, Lindblom
\cite{Lindblom1983,DetweillerLindblom1985},
and then reformulated by Chandrasekhar and Ferrari
\cite{Chandrasekhar08021991,Chandrasekhar09091991,Chandrasekhar08081991},
Ipser and Price \cite{ipser1990}, and
Kojima \cite{Kojima1992}. In neutron stars, axial modes are purely
space-time modes of oscillation (w-modes), while polar modes can be coupled
to fluid oscillations (fundamental f-modes, pressure p-modes, and also a
branch of spatial w-modes). In this paper we will consider only polar modes of
oscillation. The corresponding study of axial modes can be in found \cite{blazquez2013}.

Although the perturbation equations for neutron stars can be simplified into a
Regge-Wheeler equation or a Zerilli equation for the vacuum part of the
problem, it is complicated to obtain the quasi-normal modes spectrum of
neutron stars. The equations must be solved numerically and no
analytical
solution is known for physically acceptable configurations. The main difficulty is found on the diverging and oscillatory nature of the quasi-normal modes. These functions are
not well handled numerically.

Several methods have been developed to deal with these
and other difficulties. For a
complete review on the methods see the review by Kokkotas and Schmidt
\cite{Kokkotas_Schmidt1999}.

The axial part of the spectrum has been extensively studied for simplified constant density models and polytropes (\cite{Kokkotas1991,Andersson1995,Kokkotas1994}). Important results for the polar spectrum have been obtained using several methods, for example the continued fraction method, which have been used to study realistic equations of state \cite{Benhar_Berti_Ferrari_1999,Benhar_Ferrari_Gualtieri_2004,Wen2009}. More
recently, Samuelsson et al \cite{Samuelsson2007} used a complex-radius
approach for a constant
density configuration to study axial quasi-normal modes.

In order to use the results for future observations, Andersson and Kokkotas,
first for polytropes \cite{Anderssonprl1996} and later
for some realistic
equations of state \cite{Andersson1998}, proposed some empirical
relations between the modes and the mass and the radius of the neutron star
that could be used together with observations to
constrain the equation of state. In later works Benhar, Ferrari and Gualtieri
\cite{Benhar_Ferrari_Gualtieri_2004,benhar2007,Ferrari2007} reexamined those
relations for other realistic equations of state. Recently the authors studied
these relations for axial modes, using new equations of state
satisfying the 2 $M_{\odot}$ condition \cite{blazquez2013}.

In this work we study similar relations for the polar modes that could be used
to estimate global parameters of the star, constraining the equation of state
of the core of the neutron star. We present a new approach
to calculate quasi-normal modes of realistic neutron
stars. Essentially it is the implementation of the method used
previously for axial modes in \cite{blazquez2013} to polar modes.
We make use of several well-known
techniques, like the use of the phase
for the exterior solution and the use of a complexified coordinate to deal with
the divergence of the outgoing wave. We also introduce some new techniques not
used before in this context: freedom in the angle of the exterior complex path
of integration, use of
Colsys package to integrate all the system of equations at once
with proper boundary and junction conditions and possibility of implementation
of phase transition discontinuities. These new techniques  allow us to  enhance
precision,  to obtain more modes in shorter times, and also to
study several realistic equations of state,
comparing results for different
compositions.

In Section II we present a brief review on the quasi-normal mode formalism in
order to fix notation.
In section III we present the
numerical method, which has been implemented in Fortran based double-precision
programs. In Section IV we present the numerical results used from the
application of the method to the study of realistic equations of
state. We finish in Section V by presenting a summary of the main points of
the paper.

\section{Overview of the Formalism}
In this section we will give a brief review of the formalism used to describe
polar quasi-normal modes and the method used to calculate them. We will use
geometrized units $(c=1,G=1)$.

We consider polar perturbations on the static spherical space-time
describing the star: $ds^{2}=e^{2\nu}dt^{2}-e^{2\lambda}dr^{2} - r^{2}(d\theta^{2}
+ \sin^{2}\theta d\varphi^{2})$. The matter is considered as a perfect fluid with barotropic
equation of state. As it can be demonstrated, only polar modes do couple to
the matter perturbations of the star \cite{ThorneI}-\cite{Kojima1992}, i.e.,
first order perturbations of the energy density and pressure. Axial modes,
which were studied in \cite{blazquez2013}, are only coupled to the Lagrangian
displacement of the matter. So the spectrum
of polar modes is much more dependent of the matter
content of the star than the axial modes.

Up to first order in perturbation theory, and once the gauge is fixed
(Regge-Wheeler gauge \cite{Nollert1999}) , the
number of functions needed to describe the polar perturbations is eight, which
are, following Lindblom notation \cite{Lindblom1983,DetweillerLindblom1985}:
$(H^{0}_{lm},H^{1}_{lm},H^{2}_{lm},K_{lm},W_{lm},X_{lm},\Pi_{lm},E_{lm})$. Once
the Einstein equation is solved, together with the conservation of energy and
barion number, we obtain four first order differential equations for
$(H^{1}_{lm},K_{lm},W_{lm},X_{lm})$. There is also four algebraic relations for
$(H^{0}_{lm},H^{2}_{lm},\Pi_{lm},E_{lm})$.

The angular dependence of the functions is known from the tensor expansion and
we will only consider in this work the $l=2$ case. We can
extract the time dependence from the functions by making a Laplace
transformation resulting in equations explicitly dependent on the eigen-value
$\omega$.

Inside of the star we will need to
solve the well known zero order system of equations for hydrostatic
equilibrium of an spherical star. We need to specify the equation of state for
the matter. Once we have an static
solution we can solve the system of four differential equations for the
perturbations.

Outside the star (no matter), it can be seen that the system of equations can
be rewritten and it is reduced to a single second order differential
equation (Zerilli equation):
\begin{equation}
\frac{d^{2}Z_{lm}}{dy^{2}}+\left[\omega^{2} - V(r) \right]Z_{lm}=0, \label{eq_Z}
\end{equation}
where $y$ is the tortoise coordinate,
\begin{equation}
y=\int_{0}^{r} e^{\lambda - \nu}dr,
\end{equation}
and the eigen-frequency of the polar mode is a complex number $\omega =
\omega_{\Re} + i\omega_{\Im}$.
The potential
can be written as:
\begin{equation}
V(r) =  2(r-2M)\tfrac{n^2(n+1)r^3+3Mn^2r^2+9M^2nr+9M^3}{r^{4}(nr+3M)^2},
\end{equation}
where $n=(l+2)(l-1)/2$ and $l=2$ in this work.

\section{Numerical Method}

Note that in general the perturbation functions are complex
functions. So numerically we will have to integrate eight real first order
differential equations for the perturbation functions inside the star. These
also translates into a set of two real second order differential equations
outside the star.

The perturbation has to satisfy  a set of boundary conditions that
can be obtained from the following two requirements
\cite{Chandrasekhar08021991}: i) the perturbation must be regular at the
center of the star, and ii) the resulting quasi-normal mode must be a pure
outgoing wave.

 In general a quasi-normal mode will be
a composition of incoming and outgoing waves, i.e.
\begin{equation}
\lim_{r_{*}\to\infty}Z^{in} \sim e^{i\omega r_{*}},
\lim_{r_{*}\to\infty}Z^{out} \sim e^{-i\omega r_{*}}. \label{asymp_out}
\end{equation}
Note that, while the real part of $\omega$ determines the oscillation frequency
of the wave, the imaginary part of the eigen-value determines the asymptotic
behavior of the quasi-normal mode.
Outgoing quasi-normal modes are divergent
at radial infinity, while ingoing ones tends exponentially to zero as the
radius grows. We should impose the solution to behave as a purely outgoing quasi-normal mode
at a very
distant point. This has the inconvenience that any small contamination of the outgoing signal gives rise to an important incoming wave component.

We have adapted a method based on Exterior Complex Scaling method, that we
will describe briefly in the following paragraphs. It is an adaptation to
polar modes of the method previously employed for axial quasi-normal modes \cite{blazquez2013}.

Outside the star, and following
\cite{Andersson1992}-\cite{Samuelsson2007}, we study the phase function
(logarithmic derivative of the Z function),
which does not oscillate towards asymptotic infinity.
Although the phase allows us to
eliminate the oscillation, we can not distinguish between the mixed
incoming and outgoing waves and the purely outgoing ones.
To do so we can study an analytical extension of the phase function, rotating the radial coordinate into the
complex plane. This technique of
complexification of the integration variable is called Exterior Complex
Scaling \cite{aguilar_combes_1971,balslev_combes_1971,simon1972}.

In our method, the angle of the
integration path is treated as a free parameter of the solution.
This angle
can be adjusted appropriately to enhance precision and reduce
integration time.

The use of the phase and Exterior Complex Scaling with variable angle allows
us to compactify, so
we can impose the outgoing
quasi-normal mode behavior as a
boundary condition at infinity without using any cutoff for the radial
coordinate.

The interior part is integrated satisfying the regularity
condition at the origin. We integrate together
the zero order and the Zerilli equation. In this way, we can automatically
generate more points in the zero order functions if we require a more precise
result in our quasi-normal modes.

Because we are interested in realistic equations of state, the implementation
of the EOS is a fundamental step of the procedure. We use two
possibilities. The first one, a piece-wise polytrope
approximation, done by Read et al \cite{Jocelyn2009} for 34 different equations of state, where
the EOS is approximated by a polytrope in different
density-pressure intervals. For densities below the nuclear density, the SLy
equation of state is considered.
The second kind of implementation for the equation of state, more
general, is a piecewise monotone
cubic Hermite interpolation satisfying local thermodynamic conditions
\cite{Fritsch-Carlson}.

We implement this method into
different Fortran programs and routines, making use of the Colsys package
\cite{colsys1979} to solve numerically the
differential equations. The advantage of this package is that it allows the
utilization of quite flexible multi-boundary conditions
and
an adaptative mesh that increases precision

The full solution is generated using two independent solutions of the perturbation equations
for the same static configuration and
value of $\omega$.
The proper matching of a combination of these two interior solutions with
the exterior phase will give us an outgoing quasi-normal
mode. After a careful study of the matching conditions
it can be demonstrated that the junction can be written in
terms of a determinant of a $2\times2$ matrix.
The determinant is zero only if
the
junction conditions are satisfied. That is, a combination of the interior
solutions can be matched to the exterior phase, an outgoing wave. The
determinant of this matrix, once
the static configuration is fixed, only depends on the values
$(\omega_{\Re},\omega_{\Im})$. So if for a fixed configuration, a value of the
eigen-frequency makes the determinant null, that couple of values
$(\omega_{\Re},\omega_{\Im})$ corresponds to a quasi-normal mode of the static
configuration. Then, the behavior of the perturbation functions can be
analyzed to deduce if the quasi-normal mode is a w-mode, a p-mode or a f-mode.

In figures \ref{example} we show a particular solution for the interior of a
neutron star of
$1.4 M_{\odot}$, equation of state GNH3. They correspond to the f-mode and the fundamental p-mode.
\begin{figure}
\includegraphics[angle=-90,width=0.4\textwidth]{fig1a.eps}
\includegraphics[angle=-90,width=0.4\textwidth]{fig1b.eps}
\caption{Typical form of the functions integrated inside of the star for
  f-modes and p-modes.}
\label{example}
\end{figure}

\section{Results}
In this section we present our results for the polar quasi-normal modes of
neutron stars with realistic equations of state. We consider a wide range of
equations of state in order to
study the influence of exotic matter on the polar
quasi-normal modes. For a similar study performed
on the axial component of the spectra see \cite{blazquez2013}.

We start presenting the equations of state used in the study:

Using the parametrization presented by Read et all \cite{Jocelyn2009}, we have
considered SLy, GNH3, H4, ALF2, ALF4.

After the recent measurement of the $1.97 M_\odot$ for the pulsar PSR
J164-2230, several new equations of state have been proposed. These EOS take into account the presence of exotic matter in the core of the neutron star, and they have a maximum
mass stable configuration over $1.97 M_\odot$. We have
considered the following ones:   two equations of state presented by
Weissenborn et al with hyperons in \cite{Weissen1},  we call them  WCS1 and
WCS2;
three EOS presented by  Weissenborn et al with quark matter in
\cite{Weissen2}, we call them WSPHS1, WSPHS2,  WSPHS3; four equations of state
presented by L. Bonanno and A. Sedrakian in \cite{Sedrakian}; we call them BS1,
BS2, BS3, BS4; and one EOS presented by Bednarek et al in \cite{Bednarek}, we
call it BHZBM.
\begin{figure}
\includegraphics[angle=-90,width=0.45\textwidth]{fig2.eps}
\caption{Pressure versus density in logarithmic scale for the 15
  equations of state considered, in the high density region.}
\label{eos_full}
\end{figure}

In Figure \ref{eos_full} we plot all the 15 equations of state we have
studied in order to give an idea of the range considered.

We follow the notation previously used in \cite{blazquez2013}. We discuss their characteristics in the following paragraph.

For plain $npe\mu$ nuclear matter we use:
\begin{itemize}
\item SLy \cite{Haensel2001_SLy} with $npe\mu$ using a potential-method to obtain the EOS.
\end{itemize}
For mixed hyperon-nuclear matter we use:
\begin{itemize}
\item GNH3 \cite{Glendenning1985}, a relativistic mean-field theory EOS
  containing hyperons.
\item H4 \cite{lackey_2006}, a variant of the GNH3 equation of state.
\item WCS1 and WCS2 \cite{Weissen1}, two equations of state with
  hyperon matter, using ``model
  $\sigma\omega\rho\phi$'', and considering ideal mixing, SU(6) quark model, and
  the symmetric-antisymmetric couplings ratio $\alpha_v=1$ and $\alpha_v=0.2$
  respectively.
\item BHZBM \cite{Bednarek}, a non-linear relativistic mean field model
  involving baryon octet coupled to meson fields.
\end{itemize}
For Hybrid stars we use:
\begin{itemize}
\item ALF2 and ALF4 \cite{Alford2005}, two hybrid EOS with mixed APR nuclear
  matter and color-flavor-locked quark matter.
\item WSPHS3 \cite{Weissen2}, a hybrid star calculated using bag
  model, mixed with NL3 RMF hadronic EOS. The parameters employed are:
  $B_{eff}^{1/4}=140 MeV$, $a_4=0.5$, and a Gibbs phase transition.
\end{itemize}
For Hybrid stars with hyperons and quark color-superconductivity we use:
\begin{itemize}
\item BS1, BS2, BS3, and BS4, \cite{Sedrakian}, four equations of state
  calculated using a combination of phenomenological relativistic
  hyper-nuclear density functional and an effective NJL model of quantum
  chromodynamics. The parameters considered are: vector coupling $G_V/G_S=0.6$
  and quark-hadron transition density $\rho_{tr}/\rho_{0}$ equal to
  $2$, $3$, $3.5$ and $4$ respectively, where $\rho_0$ is the density of
  nuclear saturation.
\end{itemize}
For quark stars we use:
\begin{itemize}
\item WSPHS1 and WSPHS2 \cite{Weissen2}. The first equation of
  state is for unpaired quark matter, and we have considered the parameters
  $B_{eff}^{1/4}=123.7 MeV$,$a_4=0.53$. The second equation of state considers
  quark matter in the CFL phase (paired). The parameters considered are
  $B_{eff}^{1/4}=130.5 MeV$,$a_4=0.66$,$\Delta=50 MeV$
\end{itemize}

In Figure \ref{all_frequencies} we plot the frequencies of all the f-modes,
p-modes an wI0-modes obtained for the equations of state considered, in order
to present the range in which each modes are found. In the next subsection we
will explain the features of each class of modes depending on the equation of
state composition: hyperon matter or quark matter.

\begin{figure}
\includegraphics[angle=-90,width=0.45\textwidth]{fig3.eps}
\caption{Frequency vs compactness for all the modes studied in this work.}
\label{all_frequencies}
\end{figure}

\subsection{Hyperon matter}

\textbf{Fundamental wI mode}:
At the maximum mass configuration of each EOS (around $2 M_{\odot}$),
all the equations of state yields frequencies near $8.5 Khz$. This is found for stars containing hyperon matter or only nuclear matter. The only exception is the interesting case of the WCS2 (maximum mass of $2.4 M_{\odot}$) with the smallest frequency ($7.5 Khz$). In general, for one solar mass configurations, the frequency of stars containing hyperon matter is much smaller ($11 Khz$) than the frequencies for plain nuclear matter stars ($14 Khz$).

In order to use future observations of gravitational waves to estimate  the
mass and  the gravitational radius of the neutron star, as well as
 to discriminate between different families of equations of state, we obtain
 empirical relations between the frequency and damping time of
quasi-normal modes and the compactness of the star, following
\cite{Benhar_Berti_Ferrari_1999}, \cite{Kokkotas2001}, \cite{Anderssonprl1996} and \cite{Andersson1998}. In Figure
\ref{eos_wI0_real_scaled_hyp} we present  the
frequency of the fundamental mode scaled to the radius of each
configuration. The range of the plot is between $1 M_\odot$ and the maximum
mass configuration, so we consider only stable stars. It is interesting to note that even for the
softest equations of state with hyperon matter, the scaling relation is
very well satisfied. Since the w modes do not couple to matter oscillations, it is not a surprise that this scaled relations are universal.

A linear fit can be made for each equation of state. We fit to the following
phenomenological relation:
\begin{equation}
\omega(Khz) = \frac{1}{R(Km)}\left(A\frac{M}{R} + B\right).
\label{empiricalfrec}
\end{equation}
In Table \ref{tab:wI0fit_hyp} we present the fitting parameters $A$,$B$ for each
one of the equations of state studied.
For all the EOS very similar results are obtained.
The empirical parameters are compatible
with the empirical relation  obtained by Benhar et al \cite{Benhar1999} for
six equations of state.

A similar study can be done to the damping time of the wI0 modes. In Figure
\ref{eos_wI0_imaginary_scaled_hyp} we present the
damping time of the fundamental mode scaled to the mass of each
configuration.
In this case, the results can be fitted to a empirical quadratic  relation on $M/R$, as follows:
\begin{equation}
\frac{1}{\tau(\mu s)} = \frac{1}{M(M_{\odot})}\left[a\left(\frac{M}{R}\right)^{2} + b\frac{M}{R} + c\right].
\label{empiricaltau}
\end{equation}
In Table \ref{tab:wI0fit_hyp} we present the fitting parameters a,b and c. Note
that they are quite similar for all the equations of state, and also in accordance
with the
results obtained in \cite{Benhar1999}.

In this relations, as well as in the rest of empirical relations studied in
the paper, we have used $G=c=1$, so both mass and radius must be given in
units of length (Km), unless otherwise specified.
 \begin{figure}
 \includegraphics[angle=-90,width=0.45\textwidth]{fig4.eps}
 \caption{Scaled frequency of the fundamental wI mode vs M/R for hyperon matter EOS. The phenomenological relation is quite independent of the matter composition, as expected from spatial modes}
 \label{eos_wI0_real_scaled_hyp}
 \end{figure}
 \begin{figure}
 \includegraphics[angle=-90,width=0.45\textwidth]{fig5.eps}
 \caption{Scaled damping time of the fundamental wI mode vs M/R for hyperon matter EOS.}
 \label{eos_wI0_imaginary_scaled_hyp}
 \end{figure}

For a similar study of the axial wI modes using these equations of state see \cite{blazquez2013}.

\textbf{f mode}:
Now we will present the results for the f-modes of these configurations.
For these modes, the frequency of the SLy stars are always above the hyperon stars considered.
The Sly configurations range from $1.6
kHz$ for $1 M_{\odot}$ to $2.5
kHz$ for the maximum mass configuration. Hyperon stars range from $1.3
kHz$ for $1 M_{\odot}$ to $2.2
kHz$, with again the exception of WCS2, that reaches only $2.1
kHz$.

In this case one can also study some empirical relations. It is interesting to
study the dependence of the frequency with the square root of the mean density.
The empirical relations considered for the f modes are the following. For the frequency:
\begin{equation}
\omega(Khz) = U + V\sqrt{\frac{M}{R^{3}}}.
\label{empiricalfrec_f}
\end{equation}
For the damping time:
\begin{equation}
\tau(s) = \frac{R^{4}}{cM^{3}}\left[u + v\frac{M}{R}\right]^{-1}.
\label{empiricaltau_f}
\end{equation}
In figures \ref{eos_f_real_scaled} and \ref{eos_f_imaginary_scaled} we plot
these relations. In Figure \ref{eos_f_real_scaled} it can be seen that
(\ref{empiricalfrec_f}) is quite different
depending on the EOS. Again WCS2 present a very
different
behavior close to the maximum mass configurations.
For the damping time, relation
(\ref{empiricaltau_f}) is very well satisfied for every EOS along all
the range.

In Table \ref{tab:ffit_hyp} we present the fitting parameters U and V for the frequency, and u and v for the damping time.

 \begin{figure}
 \includegraphics[angle=-90,width=0.45\textwidth]{fig6.eps}
 \caption{Scaled frequency of the fundamental f mode vs M/R for hyperon matter. Although the frequency is linear with the square root of the mean density, each equation of state has its own particular phenomenological parameters, and so the empirical relation is less useful.}
 \label{eos_f_real_scaled}
 \end{figure}
 \begin{figure}
 \includegraphics[angle=-90,width=0.45\textwidth]{fig7.eps}
 \caption{Scaled damping time of the fundamental f mode vs M/R for hyperon matter. In this case, the scaled damping time is quite independent of the equation of state.}
 \label{eos_f_imaginary_scaled}
 \end{figure}
 \begin{figure}
 \includegraphics[angle=-90,width=0.45\textwidth]{fig8.eps}
 \caption{Scaled frequency of the fundamental p mode vs M/R. At high compactness the scaled frequency is more sensible to the presence of hyperon matter, with the exception of WCS2.}
 \label{eos_p_real_scaled}
 \end{figure}

\textbf{Fundamental p mode:}
The frequency of the fundamental p mode is the most sensitive to the EOS of the configuration. For SLy EOS, the frequency ranges from $5.7 kHz$ for $1 M_\odot$ to  $7.8 kHz$ for the maximum mass configuration. In the presence of hyperon matter, the range of frequencies is much shorter, variating from $4.6 kHz$ for $1 M_\odot$ to  $5.5 kHz$ for the maximum mass configuration.

In principle, an interesting empirical
relation now could be
 \begin{equation}
 \omega(Khz) = \frac{1}{M(Km)}\left(K\frac{M}{R} + K_{0}\right).
 \label{empiricalfrec_p}
 \end{equation}
In figure \ref{eos_p_real_scaled}, we present the frequency scaled to the mass, versus the compactness. There is an important resemblance between the Sly EOS and the WCS2, which present a similar behavior for high compact stars. The relation is quite different with respect the rest of hyperon matter EOS. This effect
could be used to constraint the value of $\alpha_v$ for compact enough stars.

In Table \ref{tab:pfit_hyp} we present the fitting parameters K and K0 for the frequency.

\subsection{Quark matter}

\textbf{Fundamental wI mode:} The frequencies for the hybrid EOS ALF4 gives almost the same results as the SLy EOS. The four BS1-4 EOS are essentially similar, with a quite low frequency: $8.8 kHz$ for the maximum mass configuration ($2 M_\odot$). Nevertheless is the WSPHS3 the EOS with lowest frequency: $8 kHz$ for the maximum mass configuration ($2.2 M_\odot$). For pure quark stars WSPHS1-2, even lower frequencies ($7.2 kHz$) are obtained at the maximum mass configuration ($2.5 M_\odot$).

We study the empirical relations (\ref{empiricalfrec}) and (\ref{empiricaltau})
for these equations of state. We present them in figures
\ref{eos_wI0_real_scaled_quarks} and \ref{eos_wI0_imaginary_scaled_quarks} for
frequency and damping time respectively. It can be seen that the empirical
relation is again very well satisfied for all the hybrid equations of state. For the
pure quark stars, a similar relation with different parameters is found. The
fits can be found in table \ref{tab:wI0fit_quark}.

Note that, although the empirical relation is also satisfied for pure quark
stars, it is slightly different. This is due to the different layer structure
of these pure quark stars. These stars are essentially naked quark cores,
where the pressure quickly drops to zero in the outer layers, although the
density is almost constant.

The study of the axial wI modes using these equations of state can be found in
\cite{blazquez2013}.

\textbf{f mode:}
The lowest frequency is around $1.2
kHz$ for $1 M_{\odot}$ stars composed only of quark matter (WSPHS1-2). The
higher frequencies are reached at $2.2 kHz$ for hybrid stars beyond $2
M_{\odot}$. For ALF2-4 frequencies, they approach those of pure nuclear
matter. The configurations with EOS WSPH3 and BS1-4  have always lower
frequencies than the rest of hybrid stars, but are larger than the pure quark
stars. In the case of the damping time
there is
essentially no difference between the EOS ($0.3 s$ for maximum mass configurations), except for WSPHS1-2, which have
slightly larger damping times for the most compact ones ($0.4 s$ for maximum mass configurations).

We study the dependence of the frequency with square root of the mean density as we did before for hyperon matter.
The empirical relations to consider are again (\ref{empiricalfrec_f}) and (\ref{empiricaltau_f}).
In Figures \ref{eos_f_real_scaled_quarks} and \ref{eos_f_imaginary_scaled_quarks} we plot
these relations. In Figure \ref{eos_f_real_scaled_quarks} it can be seen that
(\ref{empiricalfrec_f}) is quite different
depending on the EOS. ALF4 is almost identical to SLy, while the other hybrid equations are in-between these ones and the pure quark EOS. Regarding the damping time (Figure \ref{eos_f_imaginary_scaled_quarks}), relation
(\ref{empiricaltau_f}) is very well satisfied for every hybrid EOS considered along all
the range. In table \ref{tab:ffit_quark} we present the fits to (\ref{empiricalfrec_f})
and (\ref{empiricaltau_f}) of these results.

\textbf{Fundamental p mode:}
At high compactness, ALF4 has higher frequencies (near $7 khz$) than the rest of hybrid
EOS (for example BS1-4 and WSPHS1-2 have similar frequencies around $5.2 kHz$). The highest frequencies are reached for pure quark stars at $1 M_\odot$, with $10 kHz$.

In Figure \ref{eos_p_real_scaled_quarks} we present the empirical relations
(\ref{empiricalfrec_p}). It can be seen that all hybrid equations lay in the
same range although the linear relation is not very well satisfied for higher
compactness. In the case of pure quark stars, the situation is much worse. In table \ref{tab:pfit_quark} we present the fits to (\ref{empiricalfrec_p}) for each equation of state.

 \begin{figure}
 \includegraphics[angle=-90,width=0.45\textwidth]{fig9.eps}
 \caption{Scaled frequency of the fundamental wI mode vs M/R for EOS containing quarks. Note that the pure quark configurations of EOS WSPHS1-2 present a different set of phenomenological parameters.}
 \label{eos_wI0_real_scaled_quarks}
 \end{figure}
 \begin{figure}
 \includegraphics[angle=-90,width=0.45\textwidth]{fig10.eps}
 \caption{Scaled damping time of the fundamental wI mode vs M/R for EOS containing quarks. Hybrid and quark stars has different parameters in the empirical relations.}
 \label{eos_wI0_imaginary_scaled_quarks}
 \end{figure}

 \begin{figure}
 \includegraphics[angle=-90,width=0.45\textwidth]{fig11.eps}
 \caption{Scaled frequency of the fundamental f mode vs M/R for EOS containing quarks. In this case, the empirical relation is dependent of the equation of state considered, and so is less useful.}
 \label{eos_f_real_scaled_quarks}
 \end{figure}
 \begin{figure}
 \includegraphics[angle=-90,width=0.45\textwidth]{fig12.eps}
 \caption{Scaled damping time of the fundamental f mode vs M/R for hybrid and quark stars. The empirical relation is almost insensitive of the EOS, except for the pure quark stars.}
 \label{eos_f_imaginary_scaled_quarks}
 \end{figure}
 \begin{figure}
 \includegraphics[angle=-90,width=0.45\textwidth]{fig13.eps}
 \caption{Scaled frequency of the fundamental p mode vs M/R. Once again a general tendency is seen, but the particular phenomenological relation of each EOS is different from the rest. As expected the p mode is very sensitive to the equation of state.}
 \label{eos_p_real_scaled_quarks}
 \end{figure}

\squeezetable
\begin{table*}
\footnotesize
\caption{Fitting for the wI0 modes. Parameters $A$ and $B$ correspond to the
  linear
  empirical relation for the frequency (\ref{empiricalfrec}). Parameters $a$,$b$
  and $c$ correspond to the quadratic empirical relation for the damping time
  (\ref{empiricaltau}).}
\begin{tabular*}{\linewidth}{@{\extracolsep{\fill}} c c c c c c c }
    \hline \hline
wI0 & SLy & GNH3 & H4 & WCS1 & WCS2 & BHZBM
 \\
  \hline
$A$         &$-365.3\pm8.2$&$-371.0\pm8.6$&$-364.4\pm1.6$&$-374.4\pm8.2$&$-365.2\pm7.2$&$-384.2\pm8.2$\\
$B$         &$195.3\pm1.8$&$195.9\pm1.6$&$194.9\pm1.6$&$196.7\pm1.6$&$195.3\pm1.5$&$198.0\pm1.8$\\
$\chi^{2}$&2.523&1.547&1.479&1.06&1.84&2.03\\
  \hline
$a$          &$-155\pm31$&$-339\pm31$&$-330\pm35$&$-401\pm16$&$-217\pm30$&$-253\pm41$\\
$b$          &$-11\pm13$&$49\pm11$&$48\pm13$&$70\pm6$&$13\pm12$&$21\pm17$\\
$c$          &$19.7\pm1.3$&$15.1\pm1.0$ &$15.1\pm1.1$&$13.34\pm0.51$&$17.6\pm1.2$&$17.2\pm1.6$\\
$\chi^{2}$ &$0.10$&$0.041$&$0.0308$&$0.005$&$0.0827$&$0.0968$\\
  \hline \hline
\end{tabular*}
\label{tab:wI0fit_hyp}
\end{table*}
\begin{table*}
\footnotesize
\caption{Fitting for the f modes. Parameters $U$ and $V$ correspond to the
  linear
  empirical relation for the frequency (\ref{empiricalfrec_f}). Parameters $u$
  and $v$ correspond to the linear empirical relation for the damping time
  (\ref{empiricaltau_f}).}
\begin{tabular*}{\linewidth}{@{\extracolsep{\fill}} c c c c c c c }
    \hline \hline
f & SLy & GNH3 & H4 & WCS1 & WCS2 & BHZBM
 \\
  \hline
$U$         &$0.7041\pm0.0063$&$0.338\pm0.018$&$0.337\pm0.055$&$0.257\pm0.055$&$0.5330\pm0.0083$&$0.328\pm0.026$\\
$V$         &$32.44\pm0.15$&$43.07\pm0.55$&$41.5\pm1.1$&$44.5\pm1.6$&$34.52\pm0.24$&$41.93\pm0.68$\\
$\chi^{2}$&$1.58\cdot10^{-5}$&$1.82\cdot10^{-4}$&$3.75\cdot10^{-4}$&$7.41\cdot10^{-4}$&$2.77\cdot10^{-5}$&$3.14\cdot10^{-4}$\\
  \hline
$u$          &$0.065\pm0.005$&$0.079\pm0.004$&$0.071\pm0.006$&$0.070\pm0.004$&$0.067\pm0.006$&$0.063\pm0.005$\\
$v$          &$-0.221\pm0.023$&$-0.295\pm0.021$&$-0.255\pm0.030$&$-0.251\pm0.022$&$-0.231\pm0.021$&$-0.213\pm0.021$\\
$\chi^{2}$ &$1.98\cdot10^{-5}$&$1.32\cdot10^{-5}$&$1.72\cdot10^{-5}$&$7.18\cdot10^{-6}$&$2.42\cdot10^{-5}$&$1.33\cdot10^{-5}$\\
  \hline \hline
\end{tabular*}
\label{tab:ffit_hyp}
\end{table*}
\begin{table*}
\footnotesize
\caption{Fitting for the p modes. Parameters $K$ and $K_0$ correspond to the
  linear
  empirical relation for the frequency (\ref{empiricalfrec_p}).}
\begin{tabular*}{\linewidth}{@{\extracolsep{\fill}} c c c c c c c }
    \hline \hline
p & SLy & GNH3 & H4 & WCS1 & WCS2 & BHZBM
 \\
  \hline
$K$         &$74.9\pm4.5$&$62.2\pm3.3$&$61.8\pm3.3$&$59.8\pm4.8$&$66.8\pm2.7$&$54.4\pm4.5$\\
$K_0$         &$-0.68\pm1.0$&$0.368\pm0.634$&$0.74\pm0.64$&$1.54\pm0.90$&$0.616\pm0.553$&$2.33\pm0.96$\\
$\chi^{2}$&$0.65$&$0.218$&$0.249$&$0.401$&$0.260$&$0.590$\\
  \hline \hline
\end{tabular*}
\label{tab:pfit_hyp}
\end{table*}
\begin{table*}
\caption{Fitting for the wI0 modes. Parameters $A$ and $B$ correspond to the
  linear
  empirical relation for the frequency (\ref{empiricalfrec}). Parameters $a$,$b$
  and $c$ correspond to the quadratic empirical relation for the damping time
  (\ref{empiricaltau}).}
\begin{tabular*}{\linewidth}{@{\extracolsep{\fill}} c c c c c c c c c c}
    \hline \hline
wI0 & ALF2 & ALF4 & WSPHS1 & WSPHS2 & WSPHS3 & BS1 & BS2 & BS3 & BS4
 \\
  \hline
$A$         &$-365.3\pm11.9$&$-368.9\pm6.7$&$-351.9\pm12.3$&$-365.5\pm10.3$&$-342.1\pm6.7$&$-351.4\pm4.7$&$-350.0\pm0.30$&$349.6\pm2.7$&$353.2\pm3.2$\\
$B$         &$197.1\pm2.3$&$193.3\pm1.4$&$199.0\pm2.5$&$202.0\pm2.2$&$191.8\pm1.3$&$194.2\pm0.8$&$194.0\pm0.5$&$193.9\pm0.5$&$194.4\pm0.6$\\
$\chi^{2}$&2.227&1.164&3.326&2.039&0.687&0.251&0.11&0.086&0.149\\
  \hline
$a$          &$-347\pm18$&$-281\pm19$&$-333.9\pm8.7$&$-332\pm24$&$-339\pm21$&$-339\pm21$&$-473\pm21$&$-492\pm29$&$-432\pm34$\\
$b$          &$56.1\pm6.9$&$35.9\pm7.4$&$52.6\pm3.2$&$53.1\pm9.6$&$56.3\pm7.6$&$92.3\pm6.6$&$99.0\pm9.2$&$80\pm11$&$79.9\pm8.0$\\
$c$          &$14.66\pm0.63$&$15.93\pm0.67$&$16.10\pm0.27$&$16.03\pm0.92$&$13.97\pm0.65$&$11.98\pm0.51$&$11.46\pm0.71$&$12.79\pm0.85$&$12.82\pm0.63$\\
$\chi^{2}$ &$0.0079$&$0.013$&$0.0037$&$0.0204$&$0.0133$&$0.0049$&$0.0096$&$0.0148$&$0.0102$\\
  \hline \hline
\end{tabular*}
\label{tab:wI0fit_quark}
\end{table*}
\begin{table*}
\caption{Fitting for the f modes. Parameters $U$ and $V$ correspond to the
  linear empirical relation for the frequency (\ref{empiricalfrec_f}). Parameters $u$
  and $v$ correspond to the linear empirical relation for the damping time
  (\ref{empiricaltau_f}).}
\begin{tabular*}{\linewidth}{@{\extracolsep{\fill}} c c c c c c c c c c}
    \hline \hline
f & ALF2 & ALF4 & WSPHS1 & WSPHS2 & WSPHS3 & BS1 & BS2 & BS3 & BS4 \\
  \hline
$U$         &$0.38\pm0.06$&$0.79\pm0.04$&$-0.24\pm0.07$&$-0.28\pm0.07$&$0.71\pm0.06$&$0.37\pm0.04$&$0.36\pm0.04$&$0.38\pm0.04$&$0.36\pm0.03$\\
$V$         &$38.7\pm1.6$&$30.0\pm1.1$&$48.9\pm2.1$&$49.4\pm2.1$&$29.2\pm1.9$&$98.3\pm1.2$&$38.6\pm1.3$&$37.9\pm1.1$&$38.8\pm1.1$\\
$\chi^{2}$&$6.03\cdot10^{-4}$&$4.66\cdot10^{-4}$&$2.15\cdot10^{-4}$&$1.67\cdot10^{-4}$&$5.58\cdot10^{-4}$&$1.69\cdot10^{-4}$&$2.04\cdot10^{-4}$&$1.58\cdot10^{-4}$&$1.77\cdot10^{-5}$\\
  \hline
$10u$         &$0.59\pm0.04$&$0.70\pm0.05$&$0.32\pm0.01$&$0.31\pm0.02$&$0.79\pm0.08$&$0.72\pm0.04$&$0.73\pm0.04$&$0.73\pm0.04$&$0.71\pm0.04$\\
$v$         &$0.21\pm0.02$&$0.25\pm0.02$&$0.10\pm0.01$&$0.09\pm0.01$&$0.30\pm0.04$&$0.28\pm0.03$&$0.29\pm0.02$&$0.28\pm0.02$&$0.27\pm0.02$\\
$\chi^{2}$&$6.91\cdot10^{-6}$&$1.45\cdot10^{-5}$&$7.74\cdot10^{-7}$&$6.69\cdot10^{-7}$&$2.48\cdot10^{-5}$&$7.55\cdot10^{-6}$&$6.58\cdot10^{-6}$&$6.84\cdot10^{-6}$&$7.03\cdot10^{-6}$\\
  \hline \hline
\end{tabular*}
\label{tab:ffit_quark}
\end{table*}
\begin{table*}
\caption{Fitting for the p modes. Parameters $K$ and $K_0$ correspond to the
  linear
  empirical relation for the frequency (\ref{empiricalfrec_p}).}
\begin{tabular*}{\linewidth}{@{\extracolsep{\fill}} c c c c c c c c c c}
    \hline \hline
p & ALF2 & ALF4 & WSPHS1 & WSPHS2 & WSPHS3 & BS1 & BS2 & BS3 & BS4
 \\
  \hline
$K$         &$54.5\pm8.9$&$79.9\pm5.0$&$36.1\pm5.3$&$22.8\pm10.1$&$92.1\pm3.0$&$64.1\pm2.4$&$65.0\pm2.2$&$65.1\pm2.0$&$63.7\pm2.2$\\
$K_0$         &$3.86\pm1.73$&$-1.0\pm1.1$&$10.9\pm1.1$&$13.8\pm2.2$&$-3.56\pm0.60$&$1.5\pm0.4$&$1.38\pm0.40$&$1.36\pm0.35$&$1.58\pm0.41$\\
$\chi^{2}$ &$1.24$&$0.556$&$0.322$&$1.66$&$0.153$&$0.067$&$0.059$&$0.047$&$0.072$\\
  \hline \hline
\end{tabular*}
\label{tab:pfit_quark}
\end{table*}

 \begin{figure}
 \includegraphics[angle=-90,width=0.45\textwidth]{fig14.eps}
 \caption{All fundamental wI modes normalized to the square root of the central pressure (equation \ref{empiricalpc}). This scaling make all the wI0 modes lie in the same curve, and could be used to estimate the central pressure as explained in the text.}
 \label{all_wI0}
 \end{figure}
 \begin{figure}
 \includegraphics[angle=-90,width=0.45\textwidth]{fig15.eps}
 \caption{All fundamental f modes scaled as explained in the text. The real
   part of the eigen-value is scaled to the inverse of the radius. The
   imaginary part to the inverse of the mass. Since here we choose $c=G=1$,
   $\omega_R/R$ and $\omega_I/M$ have units of $cm^{-2}$. This universal relation can be very useful to estimate the radius or mass of a neutron star. Note that all the equations of state, with the exception of the pure quark stars, are found to satisfy the empirical relation \ref{empiricalf}.}
 \label{all_f}
 \end{figure}

\subsection{Universal Phenomenological Relations}

In order to use wave detections coming from neutron stars to estimate global
properties like radius or mass, we would like to have empirical relations as
much independent from the matter content of the configuration as possible. The
empirical relations for the w modes are quite universal, and hence are useful
for asteroseismology. On the other hand, the scaled relations for the f mode are still EOS dependent, specially for the frequency. Here we propose some new phenomenological relations that we think could be useful for asteroseismology.

First we present a new empirical relation for the wI
fundamental mode. A similar empirical relation for the axial wI modes were first
studied in \cite{blazquez2013}.

As it can be seen in figure \ref{all_wI0}, if the real and the imaginary part
are scaled in the following way:
\begin{eqnarray}
\bar\omega_{R} &=&2 \pi  \frac{1}{\sqrt{p_c(cm^{-2})}}\frac{10^3}{c} \omega(Khz),\nonumber\\
\bar\omega_{I}&=& \frac{1}{\sqrt{p_c(cm^{-2})}}\frac{10^6}{c} \frac{1}{\tau(\mu
  s)} ,\label{empiricalpc}
\end{eqnarray}
that is, normalized in units of the central pressure, taking $c=G=1$, then we obtain the following relation between real and imaginary part of the eigen-value for all the EOS (except WCSHS1-2):
\begin{eqnarray}
\bar\omega_{I}=(6.146\pm0.039)10^{-3}\bar\omega_{R}^2 + \nonumber\\(-5.57\pm0.18)10^{-7}\bar\omega_{R}^4 \label{frequencypcwI0},
\end{eqnarray}
with $\chi^{2}=0.0214$.
The relation is very well satisfied independently of the equation of state, in
particular for high density stars. This could be used to approximate the
central pressure of the star.

Note that although the empirical relation between $\bar \omega_{R}$ and $\bar
\omega_{I}$ is quite independent of the EOS, the parametrization of the curve
in terms of the central pressure
is EOS dependent. So if the frequency $\omega(Khz)$ and
the damping time $\tau(\mu s)$ are known, we can parametrize a line defining $\bar \omega_{R}$ and $\bar \omega_{I}$, with parameter $p_c$, using the observed frequency and damping time. The observed values will give us the slope of the line. The
crossing point of this line with the empirical relation presented in Figure \ref{all_wI0} give us an estimation of the central pressure
$p_c$ independent of the EOS. Now, we can check which EOS is compatible with this $p_c$, i.e., which one
have the measured wI0 mode near the crossing point for the estimated central
pressure. Hence, this method could be used to constrain the equation of state. Note that if
mass and radius are already measured, we would have another filter to impose to the EOS.

Concerning the f mode, an analogous relation can be obtained. In this case, by rescaling the real part and the imaginary part
in units of the radius and mass respectively. The empirical relation obtained is EOS independent (excluding again WSPHS1-2). It is:
\begin{eqnarray}
(\omega_{I}/M)=(-5.16\pm0.24)10^{-3} + \nonumber\\(164.63\pm0.83)10^{-5}(\omega_{R}/R) + \nonumber\\(-3.14\pm0.12)10^{-10}(\omega_{R}/R)^4,
\label{empiricalf}
\end{eqnarray}
with $\chi^{2}=1.42\cdot10^{-6}$. Note that in the formula we have chosen
$c=G=1$, so that $\omega_{R}$, $\omega_{I}$, $R$ and $M$ are all given in
units of length ($cm$). It is an almost linear relation, although
the $(\omega_{R}/R)^4$ plays an important role at low compactness
configurations. This relation is very well satisfied for
every configuration independently of the equation of state, as it can be seen in \ref{all_f}. For pure quark
stars the relation is similar but with different parameters, due to the
different density-pressure distribution at lower densities. But if only
nuclear, hyperon matter and hybrid stars are considered, the empirical
relation is
universal for these families of equations of state. If the mass of the star
emitting the gravitational wave is known by other measurement, this relation
could be
used to obtain the radius of
the star emitting the gravitational wave. Following a procedure similar to the
one explained previously for the wI modes, a detection of frequency a damping
time of a neutron star emission, whose mass is already known by other
measurements, would draw a horizontal line in \ref{all_f}.
The intersection of this line with the empirical relation  (\ref{empiricalf}) would provide
a estimation of the radius of the star.

Combining the two empirical relations (\ref{frequencypcwI0}) and (\ref{empiricalf}), the particular equation of state of the
source of the
detected emission could be constrained.

\section{Conclusions}

In this paper we have considered the polar quasi-normal modes for realistic
neutron stars. In particular, we have obtained results for the fundamental wI
mode, the f mode, and the fundamental p mode. The study has been realized for
15 realistic equations of state, all of them
satisfying the $2 M_{\odot}$ condition, and various compositions: plain nuclear matter,
hyperon matter, hybrid and pure quark matter. The numerical procedure
developed to obtain the quasi-normal modes is based on Exterior Complex
Scaling. A similar study for the axial part of the spectrum can be found in \cite{blazquez2013}.

We have considered empirical relations between the scaled frequencies and damping times of the modes with the compactness and mean density. We have studied in which cases the obtained relations are more independent of the equation of state, and hence, more useful for asteroseismology. We have found that the frequencies of the fundamental p mode and f mode are quite sensitive to the matter composition of the star, so these empirical relations are less useful.

New phenomenological relations have been studied, for the fundamental wI mode and f mode, between the real part and the imaginary part of the fundamental modes. We have found universal empirical relations, independent of the particular equation of state. These relations can be useful in neutron star asteroseismology, allowing to constrain the equation of state, and so providing insight into the state of matter at high densities.

\begin{acknowledgements}
We would like to thank I. Bednarek for kindly providing
us with the EOS BHZBM, I. Sagert for EOS WSPHS1-3, A. Sedrakian for EOS BS1-4,
and S. Weissenborn for EOS WCS1-2.
We thank Dr. Daniela Doneva for valuable comments on our work, and Prof. Gabriel A. Galindo for his help concerning the Exterior Complex Scaling method.
This work was supported by  the Spanish Ministerio de Ciencia e Innovacion, research project FIS2011-28013. J.L.B was supported by the Spanish Universidad Complutense de Madrid.
\end{acknowledgements}