\section{Introduction}

Gravitational waves (GWs) are fluctuations in the fabric of
spacetime produced by the motion of massive celestial objects. The
scientific curiosity and struggles to detect these waves by the
Earth based detectors lead to the invention of laser interferometer
detectors such as LIGO, VIRGO, GEO and LISA (Bassan 2014). The most
promising source for these detectors is merging the compact binaries
composed of neutron star-neutron star, neutron star-black hole and
black hole-black hole. These orbiting systems loose their energy in
the form of GWs which speed up their orbital motion and this process
ends up at the merging of orbiting objects. Recently, LIGO
scientific and Virgo collaborations (Abbott \textit{et al}. 2016)
detected these waves and provided two observational evidences (with
signals known as GW150914 and GW151226) for GWs each of which is the
result of a pair of colliding black holes.

Polarization of a wave gives information about the geometrical
orientation of oscillations. A common method to discuss polarization
modes (PMs) of GWs is the linearized theory consisting of metric
perturbations around Minkowski background. Newman and Penrose (1962)
introduced tetrad and spinor formalism in general relativity (GR) to
deal with radiation theory. Eardley \textit{et al}. (1973) used this
formalism for linearized gravity and showed that six Newman-Penrose
(NP) parameters for plane null waves represent six polarization
modes (amplitudes) of these GWs. They also introduced
Lorentz-invariant E(2) classification of plane null waves.

Hawking (1971) found an upper bound for the energy of gravitational
radiation emitted by the collision of two black holes. Wagoner
(1984) investigated gravitational radiation emitted by accreting
neutron stars. Culter and Flanagan (1994) explored the extent of
accuracy of the distance to source and masses as well as spin of two
bodies measured by the detectors LIGO and VIRGO from the
gravitational wave signal. Turner (1997) worked on GWs produced by
inflation and discussed the potential of cosmic microwave background
anisotropy as well as laser interferometers (LIGO, VIRGO, GEO and
LISA) for the detection of GWs. Langlois \textit{et al}. (2000)
studied the evolution of GWs for a brane embedded in
five-dimensional anti-de Sitter universe and showed that a discrete
normalizable massless graviton mode exists during slow roll
inflation.

Recent indications of accelerated expansion of the universe caused
by dark energy introduced much interest in cosmology. The mysteries
of dark energy and dark matter (invisible matter) leads to modified
theories of gravity obtained by either modifying matter part or
geometric part of the Einstein-Hilbert action. A direct
generalization of GR is the $f(R)$ theory in which the Ricci scalar
$R$ in the Einstein-Hilbert action is replaced by its generic
function $f(R)$. De Felice and Tsujikawa (2010) presented a
comprehensive study on various applications of f(R) theory to
cosmology and astrophysics. Starobinsky (1980) proposed the first
inflationary model in $f(R)$ gravity compatible with anisotropies of
cosmic microwave background radiation. Hu and Sawicki (2007)
proposed a class of $f(R)$ models without cosmological constant that
satisfy cosmological and solar system tests for small field limit of
the parameter space. Tsujikawa (2008) explored observational
consequences of $f(R)$ models that satisfy the local gravity
constraints. Bamba {\it et al.} (2010) introduced $f(R)$ model which
explains inflation and late cosmic expansion at the same time.

A lot of work has been done for PMs of GWs in $f(R)$ as well as in
other modified theories. Capozziello \textit{et al}. (2008)
investigated PMs of GWs in $f(R)$ gravity and concluded that for
every $f(R)$ model there is an extra mode than GR called massive
longitudinal mode. They also worked out the response function of GWs
with LISA. Alves \textit{et al}. (2009) discussed PMs of GWs for
particular $f(R)$ model concluding the same results. They showed
that five non-zero PMs exist for a specific form of quadratic
gravity. The topic of gravitational radiation for linearized $f(R)$
theory has also been discussed in literature (Berry and Gair 2011;
N$\ddot{a}$f and Jetzer 2011). Capozziello and Stabile (2015)
studied GWs in the context of general fourth order gravity and
discussed the states of polarization and helicity. Kausar \textit{et
al}. (2016) found that for any $f(R)$ model there are two extra
modes as compared to GR. Alves \textit{et al}. (2016) explored these
modes in $f(R,T)$ as well as $f(R,T^{\phi})$ theories concluding
that the earlier one reduces to $f(R)$ in vacuum while PMs for the
later depend on the expression of $T^{\phi}$.

Herrera {\it et al.} (2015a, 2015b) studied the presence of
gravitational radiation in GR for perfect as well as dissipative
dust fluid with axial symmetry using super-Poynting vector and
showed that both fluids do not produce gravitational radiation. We
have investigated that the axial dissipative dust acts as a source
of gravitational radiation in $f(R)$ theory (Sharif and Siddiqa
2017). This paper is devoted to find PMs for some viable dark energy
models of this gravity. The paper is organized as follows. In the
next section, we write down field equations and dark energy models
of $f(R)$ gravity. We then find PMs of GWs for three models in its
subsections. Finally, we conclude our results.

\section{Dark Energy Models in $f(R)$ Gravity}

The $f(R)$ gravity action is defined as
\begin{equation}\label{1}
S=\frac{1}{16\pi}\int\sqrt{-g}f(R)d^{4}x+S_{M},
\end{equation}
where $S_{M}=\int\sqrt{-g}L_{M}d^{4}x$ denotes the matter action and
$L_{M}$ represents the matter Lagrangian. To discuss PMs of GWs, one
needs to investigate the linearized far field vacuum field
equations. The vacuum field equations for the action (\ref{1}) are
given by
\begin{equation}\label{2}
F(R)R_{\beta\gamma}-\frac{1}{2}f(R)g_{\beta\gamma}-
\nabla_{\beta}\nabla_{\gamma}F(R)+g_{\beta\gamma}\Box F(R)=0,
\end{equation}
where $F=\frac{df}{dR}=f_{R}$ and
$\Box=\nabla^{\alpha}\nabla_{\alpha}$ is the D'Alembertian operator.
The trace of this equation is
\begin{equation}\label{3}
RF(R)-2f(R)+3\Box F(R)=0.
\end{equation}
We assume that waves are traveling in $z$-direction, i.e., each
quantity can be a function of $z$ and $t$.

Various models of $f(R)$ gravity have been proposed in literature
(Starobinsky 1980; Hu and Sawicki 2007; Tsujikawa 2008; Bamba
\textit{et al}. 2010) describing the phenomena of early inflation
and late cosmic expansion. The model proposed by Hu and Sawicki
(2007) is reduced to the model considered by Alves \textit{et al.}
(2009) in the weak field regime (i.e., when $R<<m^{2}$, $m$ stands
for mass). Thus Hu and Sawicki model which satisfies the
cosmological and solar system tests has been indeed examined for PMs
of GWs in the low curvature case or Minkowski background. Similarly,
the Starobinsky model having consistency with the temperature
anisotropies measured by CMBR (De Felice and Tsujikawa 2010) has
also been analyzed for PMs of GWs by Kausar \textit{et al}. (2016).

Amendola \textit{et al}. (2007) derived the conditions for
cosmological viability of some dark energy models in $f(R)$ gravity.
They divided $f(R)$ models into four classes according to the
existence of a matter dominated era and the final accelerated
expansion phase or geometrical properties of the $m(r)$ curves where
$m(r)=\frac{Rf_{,RR}}{f_{,R}}$. They concluded that models of class
I are not physical, class II models asymptotically approach to de
Sitter universe, class III contains models showing strongly phantom
era and models of class IV represent non-phantom acceleration
$(\omega>-1)$. They argued that only models belonging to class II
are observationally acceptable with the final outcome of
$\Lambda$CDM model. Here we consider these observationally
acceptable models having the similar geometry of $m(r)$ curves to
discuss the PMs of GW. There are four models among the considered
models that fall in class II while the model $R+\alpha R^{-n}$ has
already been discussed by (Alves \textit{et al}. 2009) so we discuss
the remaining three models in this paper.

\subsection{Polarization Modes for $f(R)=R+\xi R^{2}-\Lambda$}

We consider the model $f(R)=R+\xi R^{2}-\Lambda$, it is assumed that
$\xi$ (an arbitrary constant) and $\Lambda$ (cosmological constant)
have positive values (Amendola \textit{et al}. 2007). This model
corresponds to $\Lambda$CDM model in the limit $\xi\rightarrow0$ and
Starobinsky inflationary model for $\Lambda\rightarrow0$. In this
case, Eq.(\ref{3}) yields
\begin{equation}\label{4}
R(1+2\xi R)-2(R+\xi R^{2}-\Lambda)+3\Box (1+2\xi R)=0,
\end{equation}
which on simplification gives
\begin{equation}\label{5}
\Box R-\frac{1}{6\xi}R=-\frac{\Lambda}{3\xi}.
\end{equation}
For the sake of simplicity, we consider gravitational waves moving
in one direction, i.e., $z$-direction. Thus Eq.(5) can be
interpreted as a non-homogeneous two-dimensional wave equation or
Klein-Gordon equation and its solution can be found using different
methods like Fourier transform and Green's function etc. Here we
obtain its solution following the technique used to solve
Klein-Gordon and Sine-Gordon equations given in (Rajaraman 1986)
which is simple as compared to other methods. According to this
method, any static solution is a wave with zero velocity and for the
systems with Lorentz invariance, once a static solution is known,
moving solutions are trivially obtained by boosting, i.e.,
transforming to a moving coordinate frame. Since we are considering
the vacuum field equations and the background metric is Minkowski,
so we can apply Lorentz transformations to the Ricci scalar $R$ (a
Lorentz invariant quantity). Hence static solution of Eq.(\ref{5})
is obtained by solving
\begin{equation}\label{6}
\frac{d^{2}R}{dz^{2}}-\frac{1}{6\xi}R=-\frac{\Lambda}{3\xi},
\end{equation}
whose solution is
\begin{equation}\label{7}
R(z)=c_{1}e^{\sqrt{m}z}+c_{2}e^{-\sqrt{m}z}+2\Lambda,
\end{equation}
where $c_{1}$, $c_{2}$ are constants of integration and
$m=\frac{1}{6\xi}$.

Since our system is Lorentz invariant, the time dependent solution
is obtained from the static solution through Lorentz transformation
as
\begin{equation}\label{8}
R(z,t)=c_{1}e^{\sqrt{m}\frac{z-v
t}{\sqrt{1-v^{2}}}}+c_{2}e^{-\sqrt{m}\frac{z-v
t}{\sqrt{1-v^{2}}}}+2\Lambda,
\end{equation}
where $\sqrt{1-v^{2}}$ is the Lorentz factor and $v$ represents the
velocity of wave propagation. Also, Eq.(\ref{2}) can be rewritten as
\begin{equation}\label{9}
R_{\beta\gamma}=\frac{1}{F(R)}\left[\frac{1}{2}f(R)g_{\beta\gamma}+
\nabla_{\beta}\nabla_{\gamma}F(R)-g_{\beta\gamma}\Box F(R)\right].
\end{equation}
Replacing the values of $f(R)$ and $F(R)$, we obtain its linearized
form as
\begin{equation}\label{10}
R_{\beta\gamma}=\frac{1}{2}(R-\Lambda+2\xi\Lambda
R)g_{\beta\gamma}+2\xi \nabla_{\beta}\nabla_{\gamma}R-2\xi
g_{\beta\gamma}\Box R.
\end{equation}
The non-zero components of Ricci tensor are
\begin{eqnarray}\label{10a}
R_{tt}&=&\frac{1}{6(1-v^{2})}[3v^{2}(R-\Lambda)-(R+3\Lambda)]-\xi\Lambda
R,\\\label{10b} R_{xx}&=&\frac{1}{6}(R+\Lambda)+\xi\Lambda
R=R_{yy},\\\label{10c}
R_{zz}&=&\frac{1}{6(1-v^{2})}[3(R-\Lambda)-v^{2}(R+\Lambda)]+\xi\Lambda
R,\\\label{10d} R_{tz}&=&-\frac{v R}{3(1-v^{2})}.
\end{eqnarray}
With the help of Eqs.(\ref{A.5}) and (\ref{A.9}), we have
\begin{eqnarray}\label{10e}
\Psi_{2}=\frac{1}{12}R,\quad
\Psi_{3}=\frac{1}{2}R_{l\tilde{m}},\quad
\Phi_{22}=-\frac{1}{2}R_{ll}.
\end{eqnarray}

Now, we find the expressions of $\Psi_{3}$ and $\Phi_{22}$ using
Eq.(\ref{A.4}). For $\Psi_{3}$, it yields
\begin{eqnarray}\label{10f}
\Psi_{3}&=&\frac{1}{2}R_{l\tilde{m}}=\frac{1}{2}R_{\mu\nu}l^{\mu}\tilde{m}^{\nu},
\end{eqnarray}
which can also be written as
\begin{eqnarray}\label{10f1}
\Psi_{3}&=&\frac{1}{2}(R_{tt}l^{t}\tilde{m}^{t}+
R_{xx}l^{x}\tilde{m}^{x}+R_{yy}l^{y}\tilde{m}^{y}+
R_{zz}l^{z}\tilde{m}^{z}+R_{tz}l^{t}\tilde{m}^{z}).
\end{eqnarray}
From Eqs.(\ref{A.1}) and (\ref{A.2}), the component form of vectors
$k$, $l$, $m$ and $\tilde{m}$ can be written as
\begin{eqnarray}\label{10f2}
k^{\mu}&=&\frac{1}{\sqrt{2}}(1,0,0,1),\quad
l^{\mu}=\frac{1}{\sqrt{2}}(1,0,0,-1),\\\label{10f3}
m^{\mu}&=&\frac{1}{\sqrt{2}}(0,1,i,0),\quad
\tilde{m}^{\mu}=\frac{1}{\sqrt{2}}(0,1,-i,0).
\end{eqnarray}
Substituting all the required values in Eq.(\ref{10f1}), we obtain
\begin{eqnarray}\label{10f4}
\Psi_{3}=0.
\end{eqnarray}
Similarly, Eq.(\ref{A.4}) for $\Phi_{22}$ yields
\begin{eqnarray}\nonumber
\Phi_{22}&=&-\frac{1}{2}R_{ll}=-\frac{1}{2}R_{\mu\nu}l^{\mu}l^{\nu}=
-\frac{1}{2}(R_{tt}l^{t}l^{t}+2R_{tz}l^{t}l^{z}+R_{zz}l^{z}l^{z}).
\end{eqnarray}
Replacing the Ricci tensor components and components of $l^{\mu}$,
the above equation leads to
\begin{eqnarray}\label{10g}
\Phi_{22}=-\frac{R}{12}\left(\frac{1+v}{1-v}\right)+\frac{\Lambda(2v^{2}+3)}{12(1-v^{2})}.
\end{eqnarray}
Notice that $\Psi_{4}\neq0$ represents the tensor modes of GWs.
Since there is no expression of $\Psi_{4}$ in terms of Ricci tensor,
so it cannot be evaluated with the help of available values of Ricci
tensor and Ricci scalar (Alves \textit{et al}. 2016). It can be
observed that for $\Lambda$CDM model (when $\xi\rightarrow0$)
$\Psi_{2}$ and $\Phi_{22}$ remain non-zero.

The model, $f(R)=R+\xi R^{2}-\Lambda$, is always viable and reduces
to GR when both $\xi$ as well as $\Lambda$ approach to zero. In GR,
there are only two tensor modes of polarization associated with
$Re\Psi_{4}$ and $Im\Psi_{4}$, i.e., we have only $\Psi_{4}$
non-zero among six NP parameters. From Eq.(4), we have $R=0$ for
$\xi\rightarrow0$, $\Lambda\rightarrow0$, hence GR results are
retrieved.

\subsection{Polarization Modes for $f(R)=R^{p}(\ln \alpha R)^{q}$}

This model is observationally acceptable for $p=1$ and $q>0$
(Amendola \textit{et al}. 2007). Here we assume that $q=1$ such that
the model becomes $f(R)=R\ln \alpha R$. Substituting the values of
$f(R)$ and $F(R)$ in Eq.(\ref{3}), it gives
\begin{eqnarray}\label{15}
3\Box \ln \alpha R-R\ln\alpha R+R=0.
\end{eqnarray}
Assuming $\ln \alpha R=\phi$, this equation transforms to
\begin{eqnarray}\label{16}
\Box \phi&=&\frac{e^{\phi}}{3\alpha}(\phi-1),
\end{eqnarray}
which can also be written as
\begin{eqnarray}\label{16a}
\Box \phi&=&\frac{\partial U}{\partial \phi};\quad
U(\phi)=\frac{e^{\phi}}{3\alpha}(\phi-2).
\end{eqnarray}
First we seek for a static solution, i.e., consider $\phi=\phi(z)$
such that integration of Eq.(\ref{16a}) gives
\begin{eqnarray}\label{16b}
\frac{1}{2}\left(\frac{d\phi}{dz}\right)^{2}=U(\phi).
\end{eqnarray}
Substituting the value of $U(\phi)$ and then integrating, it follows
that
\begin{eqnarray}\label{16c}
\phi(z)=2\left[1+\text{InverseErf}\left[\frac{e z}{\sqrt{3\alpha
\pi}}+\frac{ec_{3}}{\sqrt{2\alpha \pi}}\right]^{2}\right],
\end{eqnarray}
where $c_{3}$ is a constant of integration, $e= 2.71828$ and Erf is
defined by
\begin{eqnarray}\label{16c}
Erf(z)=\frac{2}{\sqrt{\pi}}\int_{0}^{z}e^{-s^{2}}ds.
\end{eqnarray}
Using Lorentz transformation, we obtain time dependent solution
given by
\begin{eqnarray}\label{16d}
\phi(z,t)=2\left[1+\text{InverseErf}\left[
\frac{e(z-vt)}{\sqrt{1-v^{2}}\sqrt{3\alpha
\pi}}+\frac{ec_{3}}{\sqrt{2\alpha \pi}}\right]^{2}\right].
\end{eqnarray}
The expression for Ricci scalar is obtained as
\begin{eqnarray}\label{16e}
R(z,t)=\frac{1}{\alpha}\exp\left(2\left[1+\text{InverseErf}\left[
\frac{e(z-vt)}{\sqrt{1-v^{2}}\sqrt{3\alpha
\pi}}+\frac{ec_{3}}{\sqrt{2\alpha \pi}}\right]^{2}\right]\right).
\end{eqnarray}
The non-zero components of the Ricci tensor have the form
\begin{eqnarray}\nonumber
R_{tt}&=&-\frac{R}{6(v^{2}-1)}\left[\frac{(3v^{2}-1)\ln\alpha R
-2}{\ln\alpha R+1}\right],\\\nonumber
R_{xx}&=&\frac{R}{6}\left(\frac{\ln\alpha R+2}{\ln\alpha
R+1}\right)=R_{yy},\\\nonumber R_{tz}&=&\frac{Rv(\ln\alpha
R-1)}{3(v^{2}-1)(\ln\alpha R+1)},\\\nonumber
R_{zz}&=&\frac{R}{6(v^{2}-1)}\left[\frac{(v^{2}-3)\ln\alpha
R+2v^{2}}{\ln\alpha R+1}\right].
\end{eqnarray}
Finally, the NP parameters for this case are
\begin{eqnarray}\label{b1}
\Psi_{2}=\frac{1}{12}R,\quad \Psi_{3}=0,\quad
\Phi_{22}=-\frac{R}{12}\left(\frac{1+v}{1-v}\right)\frac{(\ln\alpha
R-1)}{\ln\alpha R+1}.
\end{eqnarray}
Here $\Psi_{4}$ is also a non-vanishing NP parameter as discussed in
the previous case.

\subsection{Polarization Modes for $f(R)=R^{p}e^{\frac{q}{R}}$}

This model is observationally acceptable for $p=1$, so we take
$f(R)=Re^{\frac{q}{R}}$ (Amendola \textit{et al}. 2007). This model
reduces to GR when $q=0$ and consequently gives no additional PMs.
Thus to find extra PMs, we consider $q\neq0$ in further
calculations. For this model, the trace equation (\ref{3}) becomes
\begin{eqnarray}\label{3a}
Re^{\frac{q}{R}}\left(1-\frac{q}{R}\right)-2Re^{\frac{q}{R}}+
3\Box\left(e^{\frac{q}{R}}\left(1-\frac{q}{R}\right)\right)=0.
\end{eqnarray}
In low curvature regime, we have $R<<q$ which reduces the above
equation to the following
\begin{eqnarray}\label{3b}
\Box\left(\frac{1}{R}e^{\frac{q}{R}}\right)+\frac{1}{3}e^{\frac{q}{R}}=0.
\end{eqnarray}
Replacing $\frac{1}{R}=u$ and $u=u(z)$ for static solution, we
obtain
\begin{eqnarray}\label{3c}
\frac{d^{2}}{dz^{2}}(ue^{qu})+\frac{1}{3}e^{qu}=0.
\end{eqnarray}
Solving the double derivative of the above equation, it becomes
\begin{eqnarray}\nonumber
(1+qu)\frac{d^{2}u}{dz^{2}}+q(qu+2)\left(\frac{du}{dz}\right)^{2}+\frac{1}{3}=0.
\end{eqnarray}
The is a non-homogeneous non-linear second order differential
equation and does not provide an exact analytic solution unless we
make some assumptions to simplify it. Since we are working in the
weak-field regime, so $R$ is very small. Assuming $q$ to be very
large, we have $qu=\frac{q}{R}$ (as $u=\frac{1}{R}$) to be very
large such that $(qu+1)\approx qu$ as well as $(qu+2)\approx qu$ and
the above equation reduces to
\begin{eqnarray}\nonumber
u\frac{d^{2}u}{dz^{2}}+qu\left(\frac{du}{dz}\right)^{2}+\frac{1}{3q}=0.
\end{eqnarray}
Here $\frac{1}{3q}\rightarrow 0$ as $q$ is very large, hence it
reduces to
\begin{eqnarray}\label{3d}
u\frac{d^{2}u}{dz^{2}}+qu\left(\frac{du}{dz}\right)^{2}=0
\end{eqnarray}
whose solution yields $(u=1/R)$
\begin{eqnarray}\label{3f}
R(z)&=&\left(\frac{1}{q}\ln[q(c_{4}z+c_5)]\right)^{-1},
\end{eqnarray}
where $c_{4}$ and $c_{5}$ are integration constants. The non-static
solution becomes
\begin{eqnarray}\label{3g}
R(z,t)&=&\left(\frac{1}{q}\ln\left[q\left(c_{4}\frac{(z-vt)}{\sqrt{1-v^{2}}}
+c_{5}\right)\right]\right)^{-1}.
\end{eqnarray}
The non-vanishing components of the Ricci tensor are
\begin{eqnarray}\nonumber
R_{tt}&=&\frac{q\left[q(z-vt)^{2}+2\ln\left[q\left(c_{4}\frac{(z-vt)}{\sqrt{1-v^{2}}}
+c_{5}\right)\right]\right]c_{4}^{2}}
{2(\ln\left[q\left(c_{4}\frac{(z-vt)}{\sqrt{1-v^{2}}}
+c_{5}\right)\right])^{2}((tv-z)c_{4}-\sqrt{1-v^{2}}c_{5})^{2}}\\\label{29a}
&+&\frac{q\left[2q\sqrt{1-v^{2}}(z-vt)c_{4}c_{5}-qc_{5}^{2}(v^{2}-1)\right]}
{2(\ln\left[q\left(c_{4}\frac{(z-vt)}{\sqrt{1-v^{2}}}
+c_{5}\right)\right])^{2}((tv-z)c_{4}-\sqrt{1-v^{2}}c_{5})^{2}},
\\\nonumber
R_{xx}=R_{yy}&=&\frac{-q\left[q(z-vt)^{2}-2(v^{2}-1)\ln\left[q\left(c_{4}\frac{(z-vt)}{\sqrt{1-v^{2}}}
+c_{5}\right)\right]\right]c_{4}^{2}}
{2(\ln\left[q\left(c_{4}\frac{(z-vt)}{\sqrt{1-v^{2}}}
+c_{5}\right)\right])^{2}((tv-z)c_{4}-\sqrt{1-v^{2}}c_{5})^{2}}\\\label{29b}
&-&\frac{q\left[2q\sqrt{1-v^{2}}(z-vt)c_{4}c_{5}-qc_{5}^{2}(v^{2}-1)\right]}
{2(\ln\left[q\left(c_{4}\frac{(z-vt)}{\sqrt{1-v^{2}}}
+c_{5}\right)\right])^{2}((tv-z)c_{4}-\sqrt{1-v^{2}}c_{5})^{2}},
\\\nonumber
R_{zz}&=&\frac{-q\left[q(z-vt)^{2}-2v^{2}\ln\left[q\left(c_{4}\frac{(z-vt)}{\sqrt{1-v^{2}}}
+c_{5}\right)\right]\right]c_{4}^{2}}
{2(\ln\left[q\left(c_{4}\frac{(z-vt)}{\sqrt{1-v^{2}}}
+c_{5}\right)\right])^{2}((tv-z)c_{4}-\sqrt{1-v^{2}}c_{5})^{2}}\\\label{29c}
&-&\frac{q\left[2q\sqrt{1-v^{2}}(z-vt)c_{4}c_{5}-qc_{5}^{2}(v^{2}-1)\right]}
{2(\ln\left[q\left(c_{4}\frac{(z-vt)}{\sqrt{1-v^{2}}}
+c_{5}\right)\right])^{2}((tv-z)c_{4}-\sqrt{1-v^{2}}c_{5})^{2}},\\\label{29d}
R_{tz}&=&-\frac{qvc_{4}^{2}}{(1-v^{2})\left(\frac{(z-vt)}{\sqrt{1-v^{2}}}c_{4}+c_{5}\right)}.
\end{eqnarray}
The corresponding NP parameters are
\begin{eqnarray}\label{b3}
\Psi_{2}&=&\frac{1}{12}R,\quad \Psi_{3}=0,\\\label{b4}
\Phi_{22}&=&-\frac{Rc_{4}^{2}(v-1)^{2}}{4((tv-z)c_{4}-\sqrt{1-v^{2}}c_{5})^{2}}.
\end{eqnarray}
$\Psi_{4}$ is also non-zero.

\section{Final Remarks}

Observations suggest that our universe is facing an accelerated
expansion phase due to a mysterious factor of dark energy. Moreover,
direct observation of GWs opens up a new window of research. It
would be worthwhile to discuss combine effect of both dark energy
and GWs. In this paper, we have found PMs of GWs in the context of
$f(R)$ dark energy models. For each of the three models, we have
first obtained a static solution of differential equation in $R$ and
then applied Lorentz transformation to obtain a time dependent
solution. It is observed that due to Lorentz transformation, a
factor of $\left(\frac{1+v}{1-v}\right)$ appears in the value of
$\Phi_{22}$ mode in first and second case but it has negligible
effect because the speed of GW is comparable with the speed of
light. It can be seen that the expressions of $\Psi_{2}$ and
$\Phi_{22}$ in all cases are directly proportional to $R$ implying
that increase in $R$ enhances these modes or amplitudes. The mode
$\Phi_{22}$ for the first model (\ref{10g}) depends directly on the
model parameter $\Lambda$, for the second model (\ref{b1}), it
depends on $\ln\alpha$ (as $\ln\alpha R=\ln\alpha+\ln R$) while for
the third model (\ref{b3}), this depends on constants $c_{4}$ and
$c_{5}$.

In each case, we have found four non-zero PMs of GWs $\Psi_{2}$
(longitudinal scalar mode), $\Psi_{4}$ ($+$,$\times$ tensorial
modes) and $\Phi_{22}$ (breathing scalar mode) which is in agreement
with the results of (Kausar \textit{et al}. 2016). We have
non-vanishing $\Psi_{2}$ for each model implying that GWs for $f(R)$
dark energy models correspond to class $II_{6}$ (as mentioned in
Table \textbf{1}). This is the only observer dependent mode,
remaining modes are all observer independent (Eardley \textit{et
al}. 1973). These expressions of NP parameters representing the
amplitudes of GWs are significant due to the presence of dark energy
dominated era.

Here we elaborate the PMs of GW for some modified theories. The six
non-zero PMs are found only for the quadratic gravity with
Lagrangian density $\mathcal{L}=R+\alpha R^{2}+\gamma
R_{\mu\nu}R^{\mu\nu}$ in (Alves \textit{et al}. 2009). For $F(T)$
theory (where T is torsion scalar in teleparallelism), there are no
extra PMs from GR as shown in (Bamba \textit{et al}. 2013) and in
$f(R,T^{\phi})$ theory the number of PMs of GW depend on the
functional form of $f(R,T^{\phi})$ (Alves \textit{et al}. 2016). On
the other hand, the PMs of GW in scalar-tensor theory (Kausar 2017)
and massive Brans-Dicke theory (Sathyaprakash and Schutz 2009) are
same as in $f(R)$ theory.

The LIGO instruments in Livingston and Hanford have similar
orientations and the possibility of extra PMs than GR cannot be
excluded. Moreover, with two detectors having no electromagnetic and
neutrino counterpart, a large uncertainty is expected about the
source of event and consequently in the speed of GW. Thus the
possibility that speed of GW is less than the speed of light cannot
be excluded. Hence from this perspective, the modified theories of
gravity cannot be ruled out. The improvements to automated pipelines
and analysis techniques for the detection of future GW events are
continuously made for accurate measurements. Recently, two more
events of GWs, GW170104 (Abbott \textit{et al}. 2017a) and GW170817
(Abbott \textit{et al}. 2017b) have been detected by the advanced
interferometers. The event GW170104 is consistent with merging black
holes of masses $31M_{\odot}$ and $19M_{\odot}$ in GR while the
second one GW170817 is consistent with the binary neutron star
inspiral having masses in the range $1.17M_{\odot}-1.60M_{\odot}$.
The signal GW170817, has the association with GRB170817A detected by
Fermi-GBM and provides the first direct evidence of a link between
these mergers and short $\gamma$-ray bursts. It is expected that
future GW observations made by a network of the Earth based
interferometers could actually measure the polarization of GWs and
thus constrain $f(R)$ deviations from GR.

\vspace{1.0cm}

{\bf Acknowledgments}

\vspace{0.25cm}

We would like to thank the Higher Education Commission, Islamabad,
Pakistan for its financial support through the {\it Indigenous Ph.D.
5000 Fellowship Program Phase-II, Batch-III}. We are also grateful
to the anonymous referee for his constructive comments.

\renewcommand{\theequation}{A\arabic{equation}}
\setcounter{equation}{0}
\section*{Appendix A}

In this appendix, we first briefly describe the Newman-Penrose
formalism (Newman and Penrose 1962) to discuss gravitational waves
and then PMs as well as classification of null waves is developed
(Eardley \textit{et al}. 1973).

Newman and Penrose developed a new technique in GR with the help of
tetrad formalism and applied this to resolve the issue of outgoing
gravitational radiation. They defined the following relations
between the Cartesian $(\hat{t},\hat{x},\hat{y},\hat{z})$ and
null-tetrads $(k,l,m,\tilde{m})$
\begin{eqnarray}\label{A.1}
k&=&\frac{1}{\sqrt{2}}(\hat{t}+\hat{z}),\quad
l=\frac{1}{\sqrt{2}}(\hat{t}-\hat{z}),\\\label{A.2}
m&=&\frac{1}{\sqrt{2}}(\hat{x}+i\hat{y}),\quad
\tilde{m}=\frac{1}{\sqrt{2}}(\hat{x}-i\hat{y}),
\end{eqnarray}
which satisfy the relations
\begin{eqnarray}\label{A.3}
-k.l=m.\tilde{m}=1,\quad k.m=k.\tilde{m}=l.m=l.\tilde{m}=0.
\end{eqnarray}
Any tensor can be transformed from Cartesian to null basis by the
formula (Alves \textit{et al}. 2009)
\begin{eqnarray}\label{A.4}
S_{abc...}=S_{\alpha\beta\gamma...}a^{\alpha}b^{\beta}c^{\gamma...},
\end{eqnarray}
where $(a,b,c,...)$ vary over the set
$\left\{k,l,m,\tilde{m}\right\}$ and $(\alpha,\beta,\gamma,...)$
vary over the set $\left\{t,x,y,z\right\}$. In (Newman and Penrose
1962), the irreducible parts of the Riemann tensor, also called the
NP parameters, are defined by ten $\Psi$'s, nine $\Phi$'s and a term
$\Lambda$ (these are all algebraically independent). Eardley et al.
(1973) showed that for plane null waves (due to differential and
symmetry properties of the Riemann tensor) these NP quantities are
reduced to the set
$\left\{\Psi_{2},\Psi_{3},\Psi_{4},\Phi_{22}\right\}$. This set
consists of six NP parameters or PMs because $\Psi_{3}$ and
$\Psi_{4}$ are complex and thus represent two independent modes.
They also give formulas of these NP quantities in terms of
null-tetrad components of the Riemann tensor as
\begin{eqnarray}\label{A.5}
\Psi_{2}=-\frac{1}{6}R_{lklk},\quad
\Psi_{3}=-\frac{1}{2}R_{lkl\tilde{m}},\quad
\Psi_{4}=-R_{l\tilde{m}l\tilde{m}},\quad
\Phi_{22}=-R_{lml\tilde{m}}.
\end{eqnarray}
Following are some helpful relations of null-tetrad components of
the Riemann and Ricci tensors
\begin{eqnarray}\label{A.9}
R_{lk}=R_{lklk},\quad R_{ll}=2R_{lml\tilde{m}},\quad
R_{lm}=R_{lklm},~R_{l\tilde{m}}=R_{lkl\tilde{m}},~ R=-2R_{lk}.
\end{eqnarray}
The classification of weak plane null waves (Eardley \textit{et al}.
1973) obtained for standard observer (i.e., each observer sees the
waves traveling in $z$-direction and each observer measures the same
frequency) is given in Table \textbf{1}.\\\\\\\\

Table \textbf{1}: The E(2) Classes of Weak Plane Null Waves
\begin{table}[bht]
\centering
\begin{small}
\begin{tabular}{|c|c|}
\hline Classes&Condition for NP Parameters\\
\hline $II_{6}$&$\Psi_{2}\neq0$\\
\hline $III_{5}$&$\Psi_{2}=0$ and $\Psi_{3}\neq0$\\
\hline $N_{3}$&$\Psi_{2}=0=\Psi_{3}$, $\Psi_{4}\neq0$ and $\Phi_{22}\neq0$ \\
\hline $N_{2}$&$\Psi_{2}=0=\Psi_{3}=\Phi_{22}$ and $\Psi_{4}\neq0$\\
\hline $O_{1}$&$\Psi_{2}=0=\Psi_{3}=\Psi_{4}$ and $\Phi_{22}\neq0$\\
\hline $O_{0}$&$\Psi_{2}=0=\Psi_{3}=\Phi_{22}=\Psi_{4}$\\
\hline
\end{tabular}
\end{small}
\end{table}


\section{Introduction}

Clever observing strategies and techniques, together with technological progress, are moving comparative exoplanetology towards increasingly Earth-like planets. The coolest stars, in particular, offer clear observational advantages: compared to FGK stars, and everything else being equal, planets around M dwarfs have larger reflex motions, deeper transits (for well-aligned systems), and more favourable star-planet contrast ratios. This has long motivated our radial velocity (RV) search for planets around M dwarfs, which started with the discovery of the first planet around such a star \citep[GJ876b;][]{1998A&A...338L..67D}\footnote{also detected by \citet{Marcy:1998jv}}. This now amounts to almost 40 detections, which include a few Earth-mass planets and a few super-Earths located in the habitable zones of their host \citep[e.g.][]{2017A&A...602A..88A}. M dwarfs have also been the focus of several other planet searches with spectacular discoveries, including Proxima Cen b \citep{2016Natur.536..437A}, TRAPPIST-1 planets \citep{2017Natur.542..456G}, and LHS1140b \citep{2017Natur.544..333D}.

Considering their number and their well-characterised selection function, these detections provide us with statistical insights into planet formation \citep{2013A&A...549A.109B}. At the same time, many of these individual detections, and all the more so when the planetary properties such as liquid water might exist on their surface, call for follow-up studies to characterise their atmosphere and constrain their structure, composition, and chemistry.

For the subset of planets that transit, transmission and occultation spectroscopy are the characterisation methods of choice. James Webb Space Telescope (JWST) transmission spectroscopy of a few dozen coadded transits of the TRAPPIST-1 planets $b$, $c,$ and $d$, for instance, is expected to have sufficient sensitivity to detect O$_3$ in putative Earth-like atmospheres for these planets \citep{2016MNRAS.461L..92B}. This makes these planets strong candidates for a biomarker detection within the next few years, but one should remember that the TRAPPIST-1 star emits intense extreme ultraviolet (XUV) radiation and frequently produces powerful stellar flares, which together might have sterilised, if not completely stripped out, the atmospheres of at least its closer-in planets \citep{2017A&A...599L...3B, 2017ApJ...841..124V}. With stellar activity factored in, the quiet M dwarf LHS1140 and its temperate super-Earth become an appealing alternative. Both TRAPPIST-1 and LHS1140 have been given top priority for JWST Guaranteed Time observing.

Planets that do not transit are generally more difficult to characterise, but can be found closer to our Sun. This translates into both increased brightness and wider angular separation and the closest non-transiting exo-Earths might thus be amenable to characterisation. The maximum angular separation between Proxima Cen and its $b$ planet, for instance, is 37 milli-arcsec and can be resolved at visible wavelength by an 8 m class telescope. To match the daunting $10^{-7}$ planet-to-star contrast ratio, \citet{2017A&A...599A..16L} have proposed to couple the SPHERE extreme adaptive optics system and the ESPRESSO high-resolution spectrograph, which combines the contrast enhancements that one can achieve with high-resolution spectroscopy and high-contrast imaging \citep{2014Natur.509...63S, 2015A&A...576A..59S}. Under slightly optimistic assumptions, \citet{2017A&A...599A..16L} have concluded that a few dozen observing nights at the VLT would detect O$_2$, H$_2$O, and possibly CH$_4$, which
like TRAPPIST-1 represents a historic opportunity to detect biomarkers in the near future. Like TRAPPIST-1, however, Proxima Cen flares strongly and often,
which likewise challenges the habitability of its planet \citep{2016ApJ...829L..31D}.

In that context, we report the detection of a planet orbiting a 2\sfrac{1}{2} times more distant but much quieter M dwarf, Ross~128. The planet is only slightly more massive than our Earth, is temperate, and orbits a very nearby, slowly rotating, quiet M dwarf.  We discuss the properties of the star in Sect.~\ref{sect:star}, present the data in Sect.~\ref{sect:data}, and use archive photometry in Sect.~\ref{sect:rot} to determine the stellar rotation period. In Sect.~\ref{sect:rv}, we analyse the RVs and demonstrate the presence of a planet and an additional periodicity likely caused by stellar activity. The final model parameters are derived from a Markov Chain Monte Carlo (MCMC)\ algorithm with Gaussian processes in Sect.~\ref{sect:model}. In Sect.~\ref{sect:concl}, we conclude that, although K2 photometry excludes transit, the low stellar activity and moderate distance from Earth make Ross 128 $b$ a good target for biomarker searches with forthcoming telescopes.

\section{\label{sect:star}Star}

Ross~128 entered the literature as the 128th entry in the \citet{1926AJ.....36..124R} catalogue of high-proper motion stars, and has since acquired denominations including Proxima Virginis, FY~Virginis, GJ~447, HIP~57548, and LHS~315. The spectral type of this object is M4 and, owing to its proximity, it is one of the brightest representatives of this subclass (V$^{\rm mag}$=11.15, J$^{\rm mag}$=6.51, H$^{\rm mag}$=5.95, K$^{\rm mag}$=5.65). With a distance of just 3.4 parsec ($\pi=295.80\pm0.54$~mas; GAIA \citeyear{Collaboration:2016gd}), Ross 128 is the closest star in the Virgo constellation ($\alpha$=11$^{\rm h}$47$^{\rm m}$44.4$^{\rm s}$, $\delta$=$+$00$^{\rm o}$48${\rm '}$16.4${\rm ''}$; Epoch=2000). Including brown dwarfs, it is the 13$^{th}$ closest (sub-)stellar system to the Sun. Ross 128 is moving towards us and will actually become our closest neighbour in just 71,000 years from now \citep[$D_{ca}=1.9$~pc;][]{2001A&A...379..634G}.

 \citet{2015ApJ...804...64M} have derived its effective temperature $T_{\rm eff}=3192\pm60$~K, mass $M_\star=0.168\pm0.017$, radius $R_\star=0.1967\pm0.0077,$ and metallicity $[Fe/H]=-0.02\pm0.08$. Accordingly, its luminosity is $L_\star=0.00362\pm0.00039~L_\odot$. In \citet{2017A&A...600A..13A}, we measured a low \ion{Ca}{ii} emission level $\log (R'_{HK})=-5.573\pm0.082$. The calibration between $\log (R'_{HK})$ and the stellar rotation period $P_{rot}$ in the same paper converts this low calcium-line emission to an estimated rotation period of approximately 100 days, which is indicative of an age of the order of a few Gyr \citep{Newton:2016ea}.

\begin{table}
\caption{
\label{table:stellar}
Observed and inferred stellar parameters for Ross 128.
}
\begin{tabular}{l@{}lcc}
\hline
Spectral type$^{(1)}$        &                  & M4                    \\
Epoch$^{(2)}$              & & 2000\\
Right ascension, $\alpha^{(2)}$ && 11$^{\rm h}$47$^{\rm m}$44.3974$^{\rm s}$\\
Declination, $\delta^{(2)}$ && $+$00$^{\rm o}$48${\rm '}$16.395${\rm ''}$\\
Parallaxe, $\pi$$^{(2)}$       & [mas]        & $295.80\pm0.54$ \\
Distance, d$^{(2)}$                &[pc]           & $3.3806 \pm 0.0064  $    \\

\multicolumn{3}{l}{Stellar photometry}\\
~~~~V$^{(3)}$                            &[mag]         & $11.15$ \\
~~~~J$^{(4)}$                             & [mag]        & $6.505\pm0.023$ \\
~~~~H$^{(4)}$                           & [mag]        & $5.945\pm0.024$ \\
~~~~K$^{(4)}$                           &  [mag]        & $5.654 \pm 0.024$     \\
Effective temperature, T$_{\rm eff}^{(5)}$                 & [K]    & 3192$\pm$60\\
Mass, $M_\star^{(5)}$       & [M$_\odot$]             & 0.168$\pm$0.017\\
Radius, $R_\star^{(5)}$       & [R$_\odot$]             & 0.1967$\pm$0.0077\\
Metallicity, [M/H]$^{(5)}$                 &     & $-0.02\pm0.08$\\
Luminosity, $L_\star$       & [$\mathrm{L_\odot}$]    &  $0.00362\pm0.00039$\\
$\log(R'_{HK})^{(6)}$  && $-5.573\pm0.082$\\
Rotation period, $P_{Rot}^{(6,7,8)}$&  [days] & 101, 121, 123\\
Age, $\tau^{(9)}$ & [Gyr] & $\gtrsim$5 & \\
\hline
\end{tabular}\\
(1): \citet{Henry:2002cf}; (2): GAIA \citeyear{Collaboration:2016gd}; (3) \citet{Landolt:1992hl}; (4): \citet{Cutri:2003tp}; (5): \citet{2015ApJ...804...64M}; (6) \citet{2017A&A...600A..13A}; (7): this work using ASAS photometry (see Sect.~\ref{sect:rot}); (8): this work from RV (see Sect.~\ref{sect:rv}) (9) \citet{Newton:2016ea} given $P_{Rot}$;
\end{table}

\section{\label{sect:data}Data}
From July 26, 2005 (BJD=2453578.46) to April 26, 2016 (BJD=2457504.7), we collected 157 observations with the HARPS spectrograph \citep{2003Msngr.114...20M, 2004A&A...423..385P}. Exposure times were fixed to 900 sec. We discarded the 158th measurement that appears in the ESO archives, which is a just a 5 second exposure (March 23, 2015; BJD=2456740.68).  We used the high-resolution mode (R=115'000), with the scientific fibre illuminated by the target and calibration fibre either unused or illuminated by the sky. The data reduction followed the same steps as in all our recent papers. Spectral extraction and calibration relied on the on-line pipeline \citep{2007A&A...468.1115L}, which also gives an initial guess for the RV. An offline processing then refines the RV measurements and their uncertainties \citep[e.g.][]{2015A&A...575A.119A, 2017A&A...602A..88A}. The line spread function changed significantly when the May 2015 upgrade of HARPS replaced its fibre link with octogonal fibres. In this work, we treat the pre- and post-upgrade data as independent time series, which appear in the figures  in red and blue, respectively. Table ~\ref{tab:rv} (only available electronically) gives the RV time series in the barycentric reference frame. Before proceeding to the next section, however, we removed the small but significant secular acceleration ($dRV/dt = 0.14 ~m/s/yr$), which we computed using the distance and proper motion of Ross 128 \citep[$\mu_\alpha=0.60526''/yr$, $\mu_\delta=-1.21926''/yr$;][]{vanLeeuwen:2007dc} and Eq. 2 of \citet{Zechmeister:2009kw}.

To complement our HARPS observations, we used archive photometry from both ASAS and K2. The All Sky Automated Survey \citep[ASAS;][]{1997AcA....47..467P} observed Ross 128 for over nine~years. We retrieved its V-band photometry extracted through the smallest ASAS aperture, ASAS MAG 0. The K2 mission \citep{Howell:2014ju} observed Ross 128 for 82 days in its Campaign 1. We retrieved  the K2 light curves detrended with the EVEREST \citep{Luger:2016fj} and POLAR \citep{Barros:2016hj} pipelines from the Mikulski Archive for Space Telescopes (MAST)\footnote{https://archive.stsci.edu/k2/}.

\begin{figure}
\includegraphics[width=\linewidth]{Fig2a.pdf}\\
\caption{\label{fig2a}Photometry of Ross 128. {\it Top:}  periodogram of the ASAS V-band photometry. {\it Middle :} ASAS photometry phase-folded to P=121.2 days. Grey filled circles are median values in 0.1-phase bins. {\it Bottom :} K2 photometry extracted with the Everest pipeline \citep{Luger:2016fj} as a function of time. We only show corrected photometry FCOR with the highest quality flag.}
\end{figure}

\section{\label{sect:rot}Stellar rotation}
Since inhomogeneities such as spots, plages, or inhibition of the convection at the surface of a rotating star can induce apparent Doppler shifts, prior knowledge of the stellar rotation helps eliminate false positive planets. The low $\log (R'_{HK})$ of Ross 128 already indicates that its rotation period is long, ${\sim}100$ days. Here, we used ASAS and K2 photometry to refine its value.

We only retained the last seven years of the more than nine years of ASAS photometry, since Ross 128 was sampled infrequently prior to BJD$=$2452500. We subtracted the median value of each observing season, clipped out all $4\sigma$ outliers, and computed the generalised Lomb-Scargle periodogram \citep[GLS; ][]{2009A&A...496..577Z}. As seen in Fig.~\ref{fig2a} (top panel), the GLS has obvious power excess for periods around 121 days with power $p_{\rm max}=0.08$. We evaluated the power threshold for a given false alarm probability (FAP) on virtual data sets generated by bootstrap with replacement. The 1\% FAP threshold is $p_{1\%}=0.07$ and the 121~days periodic signal is therefore significant. The phased photometry (middle panel of Fig.~\ref{fig2a}) shows a $\sim1\%$ semi-amplitude.

The K2 photometry has orders of magnitude better precision than the ground-based measurements and provides quasi-continuous observations during 80 days, but does not cover a full stellar rotation. The $\sim$0.4\% trend of this photometry over 80~days is compatible with ${\sim}1\%$ variations on a 121~day period (bottom panel of Fig.~\ref{fig2a}).

\section{\label{sect:rv}Evidence for an orbiting planet and additional stellar activity}
The raw pre- and post-upgrade RV time series (Fig.~\ref{fig1a}, top panel) have r.m.s. dispersions of 2.1 and 3.0~m/s, respectively, i.e. well in excess of the $\sim1.2$ m/s expected from the photon noise on the individual measurements, and a constant model has a Bayesian information criterion BIC=618. The GLS periodogram shows a prominent power excess around period of 9.9 days and several other significant peaks (Fig.~\ref{fig1a}, middle panel). The maximum power $p_{\rm max}=0.28$ is well in excess of the 1\% FAP threshold $p_{1\%}=0.17$ and the detection of a periodic signal is thus highly significant. For a visual sanity check, we phase the RVs to a 9.9 day period (Fig 1, bottom panel) and see that the signal is well sampled at every phase.

The 9.9~day period is comfortably away from the 121~day stellar rotation period (Sect.~\ref{sect:rot}) and its first few harmonics, which by itself already lends considerable confidence to its interpretation as a planet detection. A Levenberg-Marqardt adjustment of a Keplerian model has r.m.s. residuals of 1.9 and 2.6 m/s for the two time series and an overall BIC=429. The planet's orbital period is $P_b=9.9$~ days,  RV semi-amplitude is $K_1=1.7$~m/s, and eccentricity is compatible with zero.

The residuals from the 9.9~day Keplerian model thus remain well in excess of the dispersion expected from pure photon noise, and we searched for periodicities in those residuals using both GLS and Keplerian-GLS \citep[KGLS; ][]{2009A&A...496..577Z} periodograms.  Whereas the GLS measures the power of a sine fit at each period, the KGLS does so for Keplerian signals, therefore exploring periodicities for a wider range of functional shapes. The GLS (Fig.~\ref{fig3a}, top panel) has its maximum power, $p_{max}=0.23$, at a $P=51.8$~day period, and multiple other peaks above the 1\% FAP threshold. The $\sim$52 days peak remains significant in the KGLS (Fig.~\ref{fig3a}, middle panel), but the most powerful peak ($p_{max}=0.28$) is now at $P=123$~days. This period is very close to the $\sim$~121 day stellar rotation period inferred from the ASAS photometry (Sect. 4), and is certainly compatible with it after accounting for the effect of differential rotation. The residuals phase-folded with a $P=123$~day period (bottom panel of Fig.~\ref{fig3a}) suggest a coherent signal at that period with an approximate symmetry around phase 0.5. This approximate symmetry predicts excess power in the 123/2~day second harmonic of the rotation period, as is indeed observed, and the $\sim$52~day peak of the GLS additionally matches a 1~year alias of 123/2 days. The power excess in the residuals of the Keplerian fit is therefore entirely consistent with two-spotted stellar activity modulated by stellar rotation with a $\sim$120 day period. Since the spot configuration is likely to have evolved over the $\sim$11~years of the HARPS measurements, the next section models the effect of this configuration using Gaussian process regression rather than a deterministic physical model.

\begin{figure}
\includegraphics[width=\linewidth]{Fig1a.pdf}\\
\caption{\label{fig1a}HARPS radial velocities: as a function of time ({\it top}), periodogram ({\it middle}) and phase-folded to P=9.9 days ({\it bottom}). The red and blue points represent pre- and post-upgrade measurements, respectively. }
\end{figure}

\begin{figure}
\includegraphics[width=\linewidth]{Fig3a.pdf}\\
\caption{\label{fig3a}Radial velocity residuals after subtraction of the best-fit Keplerian and drift. {\it Top:} GLS periodogram is shown. {\it Middle}: KGLS periodogram is shown. {\it Bottom:} Residuals phase-folded to a $P=123$~days period are shown. The vertical red dashed line indicates the 121~day rotation period inferred from the ASAS photometry. The red and blue points represent pre- and post-upgrade measurements, respectively.}
\end{figure}

\begin{figure}
\includegraphics[width=\linewidth]{Transit.pdf}\\
\caption{\label{fig4a}Transit search in the K2 photometry.}
\end{figure}

\section{\label{sect:model}Modelling}

\newcommand{\epsilon}{\epsilon}
\newcommand{\mathcal{P}}{\mathcal{P}}

Our model of the HARPS RVs consists of a single Keplerian function representing the effect of the planetary companion, and we explore the effect of including an additional distant body modelled as a linear velocity drift. Since the RV time series contains additional signals with frequencies close to the rotational rate of the star, its harmonics, and aliases, we modelled the error term as a multivariate Gaussian distribution with a covariance matrix produced by an appropriate kernel function. This includes the effect of correlation between the data points into the model.

For the kernel function, we chose a quasi-periodic kernel,
$$
k_\mathrm{QP}(t_i, t_j) = A^2 \exp\left(-\frac{(t_i - t_j)^2}{2\tau^2} - \frac{2}{\epsilon}\sin^2{\left(\frac{\pi (t_i - t_j)}{\mathcal{P}}\right)}\right)\;,
$$
which is known to represent adequately the covariance produced by active regions rotating in and out of view \citep[e.g.][]{Haywood:2014hs, 2015MNRAS.452.2269R}. This kernel function has four hyperparameters, corresponding to the amplitude of the covariance term ($A$),  rotational period of the star ($\mathcal{P}$), covariance decay time ($\tau$), and shape parameter ($\epsilon$). To test the robustness of our results with respect to the choice of kernel function, we also explored models employing the simpler squared-exponential kernel,
$$
k_\mathrm{SE}(t_i, t_j) = A^2 \exp\left(-\frac{1}{2}\frac{(t_i - t_j)^2}{\tau^2}\right)\;,
$$
with only two hyperparameters, $A$ and $\tau$, and no periodic term. In addition, an extra white noise component was added to the model by adding the following term to each of the previous kernels:
$$
k_\mathrm{WN}(t_i, t_j) = \delta_{ij} \left[\sigma^2_i + S_i \sigma^2_J + S^+_i \left(\sigma^+_J\right)^2\right]\;,
$$
where $\delta_{ij}$ is the Kronecker delta function, $\sigma_i$ is the internal incertitude of the data point taken at time $t_i$; $\sigma_J$ and $\sigma^+_J$ are the width of the additional noise component for the pre-, and post-upgrade data, respectively; and $S_i$ is an indicator variable, whose value is one if observation $i$ is taken before the HARPS fibre upgrade and zero otherwise, and vice versa for $S^+_i$.

In summary, four models were tested and were constructed by combining the two variants for the data model, i.e. a single Keplerian (k1) or a Keplerian plus a linear velocity drift (k1d1), and the two options for the noise term model, i.e the squared-exponential kernel (sek) and quasi-periodic kernel (qpk). For the Bayesian inference of the model parameters, we set the priors listed in Table~\ref{table.priors}. The pre- and post-upgrade velocities were treated independently with a different extra white noise amplitude for each, and an offset between these velocities.

The model parameters were sampled using the MCMC algorithm described in \citet{Goodman+10} and implemented by \citet{2013ascl.soft03002F}. The initial positions of 300 walkers were randomly drawn from the prior distribution. The algorithm was run for 40'000 iterations, and the walkers were evolved to different posterior maxima. The separate maxima were identified by clustering the samples in parameter space and the marginal likelihood of each mode was estimated using the importance sampling estimator described by \citet{2013arXiv1311.0674P}. In all cases, a single mode exhibited overwhelming evidence with respect to all other secondary maxima. The walkers in the secondary maxima were then replaced by new walkers initiated in the main maximum and the algorithm was run until no further evolution of the samplers was seen. This step took between 15'000 and 40'000 iterations, depending on the model. Then 100'000 additional steps were run, on which the final inference was performed.

Results are reported in Table~\ref{table.resultsallmodels} for each tested model. The inferred results on most planet parameters are independent of the choice of model (see e.g. the marginal posterior of the velocity amplitude in Fig.~\ref{fig.posteriors}). The most notable exception is the orbital period, which exhibits a bimodal distribution, with modes centred on 9.86 and 9.88 days, in which the difference relative weight of the modes depends on whether the model includes a linear drift or not (Fig.~\ref{fig.posteriors}). The other parameters that change slightly with the inclusion of a linear drift are the velocity zero-point, the mean longitude at epoch, the amplitude of the covariance, $A$, and the offset between pre- and post-upgrade velocities. The evolution timescale hyperparameter $\tau$ marginal distribution varies significantly between models with different kernel functions.

The relative merits of each model was studied by estimating the marginal likelihood of each model using the importance sampling estimator of \citet{2013arXiv1311.0674P}. This is a biased estimator, so we explored the evolution of the estimation for each model as the size of the sample increased (Fig.~\ref{fig.evidences}). After around 5000 samples, the estimator seems to have converged. All models are approximately equally good at explaining the data with a slight preference for the squared-exponential kernel. The final inference on the model parameters was carried out by combining the samples from the four models weighted by their posterior probability; this probability was, in turn, computed assuming all four tested models form an exhaustive set, i.e. their probabilities add up to the value one. The results are listed in table~\ref{table.results} and the MAP velocity model is presented in Fig.~\ref{fig.phasefolded}.

\begin{figure*}
\center
\includegraphics[width=0.9\columnwidth]{Gl447_hist_planet1_period_merged_ylabel}
\includegraphics[width=0.9\columnwidth]{Gl447_hist_planet1_k1_merged_ylabel}
\caption{Marginal posterior distribution of the orbital period (\emph{left}) and RV semi-amplitude (\emph{right}) for the four tested models and the weighted average of all four. The period marginal posterior distribution has a more marked bimodality for models including a non-zero acceleration term.}
\label{fig.posteriors}
\end{figure*}

\begin{figure}
\center
\includegraphics[width=0.9\columnwidth]{Gl447_marginallikelihood_vs_samplesize_inset_bigfont}
\caption{Evolution of the marginal likelihood estimation for the four tested models with sample size. The known bias of the importance sampling estimator used \citep{2013arXiv1311.0674P} is evident, but seems to become negligible for sample sizes larger than around 2000. The inset shows an enlargement of the region with sample sizes between 1000 and 10'000, where a slight preference for the squared exponential models is seen. A small random noise was added in the x-direction to facilitate viewing.}
\label{fig.evidences}
\end{figure}

\begin{figure}
\center
\includegraphics[width=0.9\columnwidth]{Gl447_map_phasefolded}
\caption{Radial velocity induced by the planet on its host star, phase-folded to the maximum-a posteriori (MAP) orbital period. The effect of the secular acceleration and the Gaussian process prediction were subtracted from the data (using the MAP parameter values). The orange and blue points correspond to the pre- and post-upgrade velocities, respectively. The inferred MAP additional white noise term ($\sigma_J$ and $\sigma_J^+$) were added in quadrature to the velocity uncertainties. The grey shaded region extends between the 5th and 95th percentile of the model at each orbital phase, computed over 10,000 samples of the merged posterior distribution.}
\label{fig.phasefolded}
\end{figure}

\input{prior_tables.tex}

\begin{landscape}

\begin{table}
\caption{\label{table.resultsallmodels}Orbital elements inferred from each tested model. Values correspond to the posterior sample mean, with errors being the standard deviation of the MCMC samples. In the second line we report the 95\% highest density interval (HDI). }

\begin{tabular}{llcccc}
\hline
\hline\noalign{\smallskip}
&& \multicolumn{4}{c}{\bf Models}\\
&& k1 sek & k1d1 sek & k1 qpk & k1d1 qpk\\
\hline\noalign{\smallskip}
\multicolumn{6}{l}{\bf Noise model parameters}\\
\hline\noalign{\smallskip}
$\sigma_J$ & [m/s] & $0.14 \pm 0.13$ & $0.13 \pm 0.13$ & $0.85 \pm 0.13$ & $0.04 \pm 0.13$  \\
&& $[0.00, 0.43]$& $[0.00, 0.43]$& $[0.00, 0.42]$& $[0.00, 0.43]$ \\
\noalign{\smallskip}
$\sigma^+_J$ & [m/s] & $1.07 \pm 0.35$ & $0.43 \pm 0.35$ & $0.15 \pm 0.35$ & $0.14 \pm 0.35$  \\
&& $[0.00, 1.17]$& $[0.00, 1.12]$& $[0.00, 1.15]$& $[0.00, 1.12]$ \\
\noalign{\smallskip}
A & [m/s] & $1.76 \pm 0.34$ & $1.78 \pm 0.31$ & $1.76 \pm 0.36$ & $1.90 \pm 0.33$  \\
&& $[1.38, 2.73]$& $[1.19, 2.46]$& $[1.33, 2.78]$& $[1.22, 2.53]$ \\
\noalign{\smallskip}
$\log_{10}\tau$ & [d] & $1.155 \pm 0.092$ & $1.144 \pm 0.097$ & $1.46 \pm 0.63$ & $1.57 \pm 0.62$  \\
&& $[0.957, 1.315]$& $[0.904, 1.295]$& $[1.03, 2.31]$ $\cup$ $[2.36, 3.00]$& $[0.98, 2.34]$ $\cup$ $[2.44, 3.00]$ \\
\noalign{\smallskip}
$\mathcal{P}$ & [d] & --& --& $110.0 \pm 223.2$ & $254.7 \pm 222.6$  \\
&& & & $[60.9, 840.2]$& $[41.0, 840.2]$ \\
\noalign{\smallskip}
$\epsilon$ &  & --& --& $1.4 \pm 2.8$ & $3.3 \pm 2.8$  \\
&& & & $[0.7, 1.5]$ $\cup$ $[1.8, 10.0]$& $[0.9, 1.6]$ $\cup$ $[1.8, 10.0]$ \\
\noalign{\smallskip}
\hline\noalign{\smallskip}
\multicolumn{6}{l}{\bf Planet parameters}\\
\hline\noalign{\smallskip}
Orbital period, $P$ & [d] & $9.8596 \pm 0.0056$ & $9.8634 \pm 0.0078$ & $9.8650 \pm 0.0054$ & $9.8824 \pm 0.0077$  \\
&& $[9.8544, 9.8702]$ $\cup$ $[9.8759, 9.8830]$& $[9.8562, 9.8825]$& $[9.8542, 9.8700]$ $\cup$ $[9.8761, 9.8822]$& $[9.8558, 9.8825]$ \\
\noalign{\smallskip}
RV amplitude, $K$ & [m/s] & $1.34 \pm 0.18$ & $1.51 \pm 0.18$ & $1.28 \pm 0.18$ & $1.31 \pm 0.18$  \\
&& $[1.05, 1.78]$& $[1.05, 1.77]$& $[1.01, 1.75]$& $[1.02, 1.76]$ \\
\noalign{\smallskip}
Orbital eccentricity, $e$ & [0..1[ & $0.036 \pm 0.092$ & $0.030 \pm 0.090$ & $0.065 \pm 0.090$ & $0.044 \pm 0.089$  \\
&& $[0.000, 0.299]$& $[0.000, 0.295]$& $[0.000, 0.291]$& $[0.000, 0.290]$ \\
\noalign{\smallskip}
Mean longitude, $\lambda_0$ & [deg] & $68.0 \pm 7.7$ & $74.9 \pm 8.5$ & $73.6 \pm 7.8$ & $76.2 \pm 8.6$  \\
&& $[62.9, 93.6]$& $[62.5, 96.9]$& $[61.8, 92.8]$& $[62.3, 97.3]$ \\
\noalign{\smallskip}
\noalign{\smallskip}
$\sqrt{e} \cos{\omega}$ &  & $-0.08 \pm 0.21$ & $0.10 \pm 0.21$ & $0.25 \pm 0.21$ & $0.11 \pm 0.21$  \\
&& $[-0.20, 0.58]$& $[-0.21, 0.56]$& $[-0.23, 0.55]$& $[-0.23, 0.55]$ \\
\noalign{\smallskip}
$\sqrt{e} \sin{\omega}$ &  & $0.17 \pm 0.19$ & $-0.14 \pm 0.19$ & $0.01 \pm 0.19$ & $0.18 \pm 0.19$  \\
&& $[-0.41, 0.33]$& $[-0.39, 0.34]$& $[-0.40, 0.33]$& $[-0.38, 0.35]$ \\
\noalign{\smallskip}
Minimum mass, $M_p \sin i$ & [$\mathrm{M}_\oplus$] & $1.35 \pm 0.20$ & $1.49 \pm 0.20$ & $1.31 \pm 0.20$ & $1.30 \pm 0.20$  \\
&& $[0.98, 1.80]$& $[0.98, 1.81]$& $[0.99, 1.83]$& $[1.00, 1.82]$ \\
\noalign{\smallskip}
Semi-major axis, $a$ & [AU] & $0.0493 \pm 0.0017$ & $0.0487 \pm 0.0017$ & $0.0497 \pm 0.0017$ & $0.0492 \pm 0.0017$  \\
&& $[0.0460, 0.0529]$& $[0.0462, 0.0532]$& $[0.0462, 0.0528]$& $[0.0459, 0.0529]$ \\
\noalign{\smallskip}
\hline\noalign{\smallskip}
\multicolumn{6}{l}{\bf Planet parameters}\\
\hline\noalign{\smallskip}
Systemic velocity, $\gamma_0$ & [km/s] & $-30.89022 \pm 0.00048$ & $-30.88978 \pm 0.00049$ & $-30.89032 \pm 0.00053$ & $-30.89056 \pm 0.00053$  \\
&& $[-30.89175, -30.88976]$& $[-30.89125, -30.88914]$& $[-30.89189, -30.88970]$& $[-30.89135, -30.88913]$ \\
\noalign{\smallskip}
Linear drift, $\gamma_1$ & [m/s/year] & --& $0.48 \pm 0.15$ & --& $0.47 \pm 0.15$  \\
&& & $[0.06, 0.68]$& & $[0.06, 0.68]$ \\
\noalign{\smallskip}
Velocity offset, $\delta_{12}$ & [m/s] & $-3.6 \pm 1.1$ & $-3.0 \pm 1.2$ & $-3.5 \pm 1.2$ & $-3.9 \pm 1.2$  \\
&& $[-4.1, 0.5]$& $[-5.3, -0.4]$& $[-4.2, 0.8]$& $[-5.6, -0.6]$ \\
\noalign{\smallskip}
\hline
\end{tabular}

\end{table}

\end{landscape}

\input{paramtable.tex}

\section{\label{sect:concl}Discussion}

The $m \sin i = 1.35~m_\oplus$ Ross~128~b planet orbits Ross~128 with a 9.86~day period, and at 0.049~AU is $\sim20$ times closer to its star than the Earth is to the Sun. Since the star is $\sim280$ times less luminous than the Sun, Ross 128~$b$ receives just 1.38 times more energy than our Earth. For assumed albedos of 0.100, 0.367, or 0.750, its equilibrium temperature would thus be 294, 269, or 213 K. Using theoretically motivated albedos, the \citet{2017ApJ...845....5K} criteria place the planet firmly outside the habitable zone, while \citet{2013ApJ...765..131K}, \citet{2014ApJ...787L...2Y}, and \citet{2016ApJ...819...84K} find it outside, inside and just at the inner edge of the habitable zone. The precise location of the inner edge is therefore still uncertain, as it depends on subtle cloud-albedo feedbacks and on fine details in complex GCM models. The habitable zone most likely will not be firmly constrained until liquid water is detected (or inferred) at the surface of many planets. Meanwhile, it is probably preferable to refer to Ross 128~$b$ as a temperate planet rather than as a habitable zone planet.

A planet just 3.4~parsecs away either having liquid water or just shy of having some makes an extremely appealing characterisation target. From the occurrence rate of temperate planets measured by Kepler, \citet{2015ApJ...807...45D} estimated that the closest habitable zone planet that transits its star is approximately 11~parsecs away. Yet, a stroke of luck could certainly align a closer temperate planet to undergo transits from our position in space, and all RV detections are therefore worth following up with photometry. As for Ross~128~$b$, existing K2 photometry readily answers whether it transits or not. We phase-folded the de-trended, low-frequency filtered, POLAR K2 photometry \citep{Barros:2016hj} to the ephemeris computed in the previous section (Fig.~\ref{fig4a}).  Ross~128~$b$ unfortunately does not transit, with central transits of any planet bigger than 0.19 R$_\oplus$ excluded at least at the 99\% confidence level. Non-grazing transits of a more realistic 0.5- or a 1.0-R$_\oplus$ planet are excluded with very high confidence.

Transit spectroscopy being excluded, we turn to the potential of measuring phase curves. \citet{2017AJ....154...77S} estimated that five days of JWST observations could detect the putative atmosphere of Proxima~Cen~b \citep[see also][]{2016ApJ...832L..12K}. Ross 128~b is not as favorable however, since its host star is 1.4 times larger and, at near- or mid-infrared wavelengths, $3-4$ times fainter than Proxima Centauri. Similar JWST observations for Ross 128~b are thus likely to be prohibitively expensive.

The best odds of characterising Ross~128~b are most likely through combining the contrast improvements achieved with high-angular resolution and with high-spectral dispersion \citep{2002ApJ...578..543S}. \citet{2015A&A...576A..59S} investigated the potential of this strategy for rocky planets around our nearest neighbours and found that a putative temperate exo-Earth orbiting Proxima Cen could be detected in just 10 hours on the European ELT (E-ELT). A year later, \citep{2016Natur.536..437A} detected an actual planet with very similar properties using RV measurements, and the technique was immediately contemplated to characterise that planet. \citet{2016arXiv160903082L} proposed to upgrade the SPHERE adaptive optics system of the VLT and inject light from the location of the planet into the ESPRESSO high-resolution spectrograph to detect the planet in few tens of nights, and possibly detect its atmospheric O$_2$ with 60 nights of observations. Ross~128~b again is not quite as favorable as Prox~Cen~b, since it cannot be resolved by a 10~m-class telescope. Its 15~mas angular separation, however, will be resolved by the 39~m E-ELT at optical wavelengths ($>3\lambda/D$ in the O$_2$ bands) and its expected contrast is similar to that of Prox~Cen~b, owing to their similar radii and semi-major axes. The two host stars have similar optical apparent magnitudes, leading to similar planetary apparent magnitudes. A realistic investment of E-ELT resources can therefore most likely detect Ross~128~b with high-angular resolution plus high-dispersion spectroscopy, although not as easily as Prox~Cen~b.

On the flip side, Ross~128 is one of the quietest stars to host a temperate exo-Earth. \citet{2017ApJ...834...85N} measured an H$_\alpha$ equivalent width EW~H$_\alpha=-0.068~\AA$ which makes Ross~128 one the most quiescent M dwarfs. They classified stars as active when EW H$_\alpha < -1 \AA$ and, for comparison, measured EW H$_\alpha = -4.709 \AA$ for Proxima Cen b. Stellar activity is probably the highest concern regarding the emergence of life, and even the survival of an atmosphere, on planets orbiting M dwarfs. Restricting the target list to quiet stars would disqualify Proxima~Cen~b and leave Ross 128 b as the best temperate planet known to date. This will certainly make this new temperate exo-Earth a top target for characterisation with the ELTs.

\begin{acknowledgements}
We wish to thank the referee for his remarks that led to an improved manuscript. XB, JMA, and AW acknowledge funding from the European Research Council under the ERC Grant Agreement n. 337591-ExTrA.
N.C.S. acknowledges the support from Funda\c{c}\~ao para a Ci\^encia e a Tecnologia (FCT) through national funds and by FEDER through COMPETE2020 by grants UID/FIS/04434/2013\&POCI-01-0145-FEDER-007672 and PTDC/FIS-AST/1526/2014\&POCI-01-0145-FEDER-016886. N.C.S. was also supported by FCT through Investigador FCT contract reference IF/00169/2012/CP0150/CT0002.
\end{acknowledgements}

\bibliographystyle{aa}


\section{Introduction}\label{intro}
Observations of  type Ia supernova
\citep{Perlmutter:1998np,Riess:2004nr,Astier:2005qq,Jha:2006fm}, cosmic microwave background (CMB) radiation
\citep{Spergel:2006hy,Ade:2015rim} and baryon acoustic oscillation (BAO)
	\citep{Seo:2005ys,Blake:2011en} indicate that around $25\%$ of matter content of the universe is made of dark matter (DM). There are many scenarios to explain the nature of DM and one of the models proposes
DM might be composed totally or partially by the primordial black holes (PBHs) \citep{Blais:2002nd,Khlopov:2008qy,Frampton:2010sw}.

There are several mechanisms to explain the formation of PBHs including sharp peaks in density fluctuations \citep{GarciaBellido:1996qt}, phase transitions \citep{Jedamzik:1999am}, resonant reheating \citep{Suyama:2004mz}, tachyonic preheating \citep{Suyama:2006sr} and curvaton scenarios \citep{Kohri:2012yw,Kawasaki:2012wr,Bugaev:2013vba}. PBHs smaller than about $10^{12}$kg should have already evaporated through the Hawking radiation \citep{2010RAA....10..495K,2010PhRvD..81j4019C}. However, the massive PBHs, are unaffected by the Hawking radiation might have various cosmological consequences, such as seeds for supermassive black holes \citep{2002PhRvD..66f3505B}, generation of the large-scale structures \citep{Afshordi:2003zb} and affects on the thermal and ionization history of the universe \citep{2008ApJ...680..829R}.

The observations for searching the Massive Astrophysical Compact Halo Objects (MACHOs) by gravitational microlensing and femtolensing of $\gamma$-ray bursts excluded PBHs in the mass range of $\sim10^{-7}M\odot$ -$1 M\odot$ and $10^{14}$-$10^{17}$kg \citep{Tisserand:2006zx,Nemiroff:2001bp}.  However, assuming an extended mass function for the compact objects, it seems that various observational data  along with the dynamical constraints are consistent with PBHs as the dark matter candidate within the mass range of $10 M_\odot<M<10^3 M_\odot$ and/or $10^{17}~\text{kg} <M<10^{21}~\text{kg}$ \citep{carr,green}.  The later range for the PBHs is convenient to be written in terms of lunar
 mass, roughly as $10^{-6} <\bar{M}<10^{-2}$ (assuming the lunar mass of $M_{\rm m} \sim 7\times 10^{23}$ kg) where $\bar{M}={M}/{M_{\rm m}}$.

Gravitational lensing provides an exceptional tool for investigating the astrophysical phenomena including indirect detection of the compact objects \citep{rahvar:2015}.
The light deflection in gravitational lensing depends on the configuration of the lens mass distribution and might produce several images. The term of gravitational microlensing is used when the images from the lensing cannot be resolved by the conventional telescopes. In this case, the result of lensing is the magnification of light receiving from the source star. Taking into account the relative motion of the lens, source and the observer results in a bell shape light curve, so-called Paczynski light curve \citep{pac86}. In recent years, microlensing has been used for discovering extra-solar planets \citep{Cottle:1991zza,Gould:1992aj,Gould:2008zu,batista,muraki}, investigating the properties of a distant source stars \citep{Afonso:2001gh,Gould:2001bg,Abe:2003xb,Fields:2003zx,sajadian} and studying the structure of the Milky Way galaxy \citep{moniez}.
Moreover, in the cosmological scales, the gravitational microlensing provides a useful method for studying compact objects.

The quasar microlensing in the cosmological scales uses the caustic crossing of an ensemble of lenses \cite{2001ASPC..237..185W} for studying the distribution of stars and micro-halos around the galaxies \cite{2007A&A...475..453Z,rahvar2014}. The possibility of detection of PBHs have been studied through observation of Quasars in X-ray \citep{pbh2017} and they found no-candidate in the mass range of   $0.05 M_\odot<M<0.45 M_\odot$. Also, some microlensing observation of Quasars has been done in the
survey mode where the aim was the detection of the caustic crossing of the lenses in the halo of the strong lensed galaxies \citep{meld,Giannini}. All these observations have been done in the optical and shorter wavelengths. Here, we propose extending observations to the infrared and millimeter wavelengths where effects of wave optics is important in the light curve of sub-lunar mass lenses. The advantage of wave optics is that we can obtain more information about the parameters of lenses compared to the geometric microlensing. From the observational point of view, Spitzer  space-based telescope and the Atacama Large Millimeter/Submillimeter
Array (ALMA)  are ideal tools for studying the light curve of quasars in our desired wavelengths  \citep{2015MNRAS.452...88K, Venemans:2015hyr,2016ApJ...824..132N} .

In section (\ref{geometric}) we introduce the gravitational microlensing in the geometric optics regime. In section (\ref{WO}), we introduce the wave optics feature of gravitational microlensing and calculate the diffraction pattern from the scattering of the electromagnetic wave from a PBH on the observer plane. In section (\ref{obs}), we study the observational features of the diffraction of light from a PBH as a lens as well as we calculate the optical depth for the observation of this event. Conclusions are given in (\ref{sec:conclude}).

\section{Geometric microlensing}
\label{geometric}
The standard gravitational lensing formalism uses the geometric optics for the limit of
$\lambda \ll R_{\rm sch}$, where $\lambda$ is the wavelength of the light and $R_{\rm sch}$ is the Schwarzschild radius of the lens. However, when $ \lambda \approx R_{\rm sch}$, the wave optics features of lensing such as interference of the light from different images produces the interference pattern on the observer plane. The relative
motion of the observer with respect to the fringes results in a time variation in the intensity of light where the time-scale and the amplitude of these fringes provide more information compared to that of geometric optics.
The diffractive gravitational lensing has been studied for a system with a galaxy as a lens and a point radio source \cite{ohanian:1983}. Following this work, the caustic-crossing of quasars in the wave optics regime has been investigated, where they put a limit on the size of quasars \cite{Jaroszynski:1995cd}. The wave optics aspect of microlensing with a sub-stellar mass lens, like a free-floating planet in the galaxy, have been studied in a series of papers \citep{Heyl:2009av,Heyl:2011b,Heyl:2010hm}. Also, \citet{Mehrabi:2012dy} investigated the wave optics features for a binary lensing system near the caustic lines.

The gravity of a lens deflects light ray from a distant source and this deflection produces multiple images from a single source.
It is more convenient to write the lens equation in terms of angular scales \citep{Schneider1985}:
\begin{equation}\label{eq:lens}
\boldsymbol{\beta}=\boldsymbol{\theta}-\boldsymbol{\alpha}(\boldsymbol{\theta})\;,
\end{equation}
where $\boldsymbol{\beta}$ and $\boldsymbol{\theta}$ are the source and the image positions and $\boldsymbol{\alpha}(\boldsymbol{\theta})$ is the deflection angle. Notice that all angles are normalized to the Einstein-angle
\begin{equation}\label{eq:theta-e}
\theta_{E}=\sqrt{\frac{2R_{\rm sch}}{D_s}\frac{1-x}{x}},
\end{equation}
where $x={D_l}/{D_s}$ is the ratio of the comoving distance of the lens to the comoving distance of the source
 and $D(z)$ in $\Lambda$CDM model is given by
      \begin{equation}\label{eq:dis-cos}
      D(z)=\frac{c}{H_0}\int_0^z \frac{dz}{\sqrt{\Omega_m(1+z)^3+\Omega_{\Lambda}}}\;.
      \end{equation}
In the geometric optics limit, after
 solving equation (\ref{eq:lens}), one can find the corresponding map between the source position to the image positions. The Jacobian of transformation in the equation of $\theta = \theta(\beta)$, provides the ratio of areas in the image space to the source space which is equivalent to the magnification in the gravitational microlensing.

 In what follows, we concern the wave optics regime of the gravitational microlensing by considering the interferometry of the light rays. It is convenient to use the Fermat potential for the light ray which is proportional to the time delay between a given trajectory and a straight path. The Fermat potential for a generic lens is given by
    \begin{equation}\label{eq:ferma}
    \phi(\boldsymbol{\theta},\boldsymbol{\beta})=\frac{1}{2}(\boldsymbol{\theta}-\boldsymbol{\beta})^2-\psi(\boldsymbol{\theta})\;,
    \end{equation}
    where $\psi(\boldsymbol{\theta})$ is the gravitational potential in 2D and is defined on the lens plane as
    \begin{equation}\label{eq:psi}
    \psi(\boldsymbol{\theta})=\frac{1}{\pi}\int\frac{\Sigma(\boldsymbol{\theta'})}{\Sigma_{\rm cr}}\ln\vert\boldsymbol{\theta}-\boldsymbol{\theta}^{\prime}\vert d^2\boldsymbol{\theta}^{\prime},\;
    \end{equation}
    and $\Sigma(\boldsymbol{\theta})$ is the surface mass density of the lens and the critical mass density is given by
 \begin{equation}
 \nonumber
    \Sigma_{\rm cr}=\frac{c^2}{4\pi GD_{\rm l}(1-x)}.
    \end{equation}
    The lens equation gives rise from the Fermat principle,
   $ \nabla_{\boldsymbol{\theta}}\phi(\boldsymbol{\theta},\boldsymbol{\beta})=0,$
  and in terms of the Fermat potential, the deflection angle in equation (\ref{eq:lens}) is given by
     $\boldsymbol{\alpha(\boldsymbol{\theta})}=\boldsymbol{\nabla}\psi(\boldsymbol{\theta})$.

 \section{Diffractive microlensing}
 \label{WO}
In wave optics limit where the time delay between trajectories from a source to
observer is less than a period of light, the light rays can be considered temporally coherent and
the result is the production of the interference pattern on the observer plane. Under condition where the source and deflector are far from the observer, we can use the Huygens-Fresnel
principle for analyzing gravitational lensing. Then, every point on the
lens plane can be taken as a secondary source and the amplitude of the electromagnetic wave at each point on the observer plane is
the superposition of the light from various sources on the lens plane. This analysis can be done for a point source, however, for an extended realistic source the amplification is calculated by the superposition of the infinitesimal
incoherent sources. Finally, multiplying the superposition
of the electromagnetic wave by its complex conjugate results in
the magnification on the observer plane. The magnification for a point
source \citep{Falco:1992} is given by
\begin{equation}\label{eq:mag-wave}
\mu(\boldsymbol{\beta};k)=\frac{f^2}{4\pi^2}\vert\int e^{if\phi(\boldsymbol{\theta},\boldsymbol{\beta})}d^{2}\theta\vert^2,\;
\end{equation}
where $f=2kR_{\rm sch}$ and $k$ is the wave-number.
$\phi(\boldsymbol{\theta},\boldsymbol{\beta})$ is the
Fermat potential for a single lens and is given by:
\begin{equation}\label{eq:ferma:point}
\phi(\boldsymbol{\theta},\boldsymbol{\beta})=\frac{1}{2}(\boldsymbol{\theta}-\boldsymbol{\beta})^2-\ln\vert\boldsymbol{\theta}\vert.\;
\end{equation}

Substituting equation~(\ref{eq:ferma:point}) in equation~(\ref{eq:mag-wave}), the  magnification is given
as follows
\begin{equation}\label{eq:mag-poi-main}
\mu(\beta;f)=\frac{\pi \frac{f}{2}}{\sinh(\pi \frac{f}{2})}e^{\pi \frac{f}{2}}\vert 1F1(1-i\frac{f}{2},1,i\frac{f\beta^2}{2})\vert^2, \;
\end{equation}
where $1F1(a,b,x)$ is the confluent hypergeometric function \citep{Falco:1992}. Fig.~(\ref{fig:mag-point}) presents the magnification in terms of $\beta$ for different values of $f$.
\begin{figure}
	\centering
	\includegraphics[width=.45\textwidth]{f}
	\caption{Magnification of a point source for varies values of $f$. The solid red, dot-dashed blue and dashed green lines represent the magnification for $f=5, 10$ and $f=15$, respectively. }
	\label{fig:mag-point}
\end{figure}
The diffraction pattern is observable when $\lambda\approx R_{\rm sch}$ (i.e. $f \simeq \mathcal{O}(1)$).
By increasing $f$ (smaller wavelength or massive lens), the fringes shrink and the diffraction pattern converge to the geometric optics magnification.

To simplify equation (\ref{eq:mag-poi-main}), we expand the  Fermat potential in equation~(\ref{eq:ferma:point})
around the critical point of ($\boldsymbol{\theta}=1$, $\boldsymbol{\beta}=0$) where according to Fig. (\ref{fig:mag-point}), the light curve has peak around it, as follows:
\begin{equation}\label{eq:fer-expand}
\phi(\boldsymbol{\theta},\boldsymbol{\beta})=\theta^2 -2\theta -\theta\beta\cos\gamma,\;
\end{equation}
where  polar coordinate $(\theta,\gamma)$ is used on the lens plane.
Then, from equation~(\ref{eq:mag-wave}) the magnification simplifies to
\begin{equation}\label{eq:mag-appr}
\mu(\beta;k)=\pi f J_0^2(f\beta),
\end{equation}
where $J_0$ is the Bessel function of the first kind.  We note that the relative difference of magnification  from Eq. (\ref{eq:mag-appr}) compare to the exact equation is less than $1\%$ for sources with $\boldsymbol{\beta}<0.5$ and this difference decreases rapidly when the source moves toward the lens position (i.e. $\beta \rightarrow 0$).
In practice, it is possible to observe the light curve in two different wavelengths say $\lambda_1$ and $\lambda_2$. In this case, the relative magnification is given by
  \begin{equation}\label{eq:mag-rel}
  \frac{\mu_1}{\mu_2}=\frac{\lambda_2}{\lambda_1}\frac{J_0^2(f_1\beta)}{J_0^2(f_2\beta)},
  \end{equation}
  where close to the maximum magnification, using the series of $J_0(x)\approx 1 - \frac{x^2}{4}+\frac{x^4}{64}$, equation (\ref{eq:mag-rel}) simplifies as
     \begin{equation}\label{eq:mag-rel2}
     \frac{\mu_1}{\mu_2}=\frac{\lambda_2}{\lambda_1}\frac{1-\frac{1}{2}(f_1\beta)^2+\frac{3}{32}(f_1\beta)^4}{1-\frac{1}{2}(\frac{\lambda_2}{\lambda_1})^2(f_1\beta)^2+\frac{3}{32}(\frac{\lambda_2}{\lambda_1})^4 (f_1\beta)^4}.
     \end{equation}
 In this case the right-hand side of this equation is a function of $f_1\beta(t)$ where from the measurement of $\mu_1$ and $\mu_2$  as a function of time, we can extract $f_1\beta(t)$. This parameter depends on
 the lensing parameters as follows:
   \begin{equation}\label{eq:fy}
     f_1\beta(t)= f_1 \left(u_0^2 + (\frac{t}{t_E})^2\right)^{\frac{1}{2}},
     \end{equation}
 where $u_0$ is the minimum impact parameter
 and $t_E$ is the Einstein-crossing time. From the observation of a microlensing event in the regime of geometric optics (i.e. $\lambda\ll R_{\rm sch}$), we can extract $u_0$ and $t_E$. On the other hand, knowing the left-hand side of equation (\ref{eq:fy}) from the wave optics and right-hand side from the geometric optics at $t=0$, we  determine directly $f_1$ or mass of the lens. We note that unlike to the geometric optics, the mass of lens determine from this method is independent of the distance of lens and source as well as their relative velocities with respect to the observer.

In reality, quasars as the source in the lensing have finite sizes and this effect should be taken into account in the wave-optics calculation. For a given source, the total magnification is calculated by integrating over all individual elements on a source where these elements are independent and incoherent. Then the magnification of an extended source \citep{Mehrabi:2012dy} is given as
\begin{equation}\label{eq:mag-finite}
\mu(\boldsymbol{\beta},\rho;k)=\int_{s<\rho}\frac{I_{\rm w}(\boldsymbol{\beta})\mu(\boldsymbol{\beta};k)d^2 s}{I_{\rm w}(\boldsymbol{\beta})d^2 s},\;
\end{equation}
where $I_{\rm w}(\boldsymbol{\beta})$ is the surface brightness of the source (that might depends on the wavelength) and  $\rho={\theta_s}/{\theta_E}$ is the angular size of source normalized to the Einstein angle.  Here the integration in equation (\ref{eq:mag-finite}) is taken over the source area (i.e. $s<\rho$).

In Fig.~(\ref{fig:mag-finite}), we depict the magnification in terms of $\beta$ for three sources with different sizes in the wave-optics regime.
For the small sources, the magnification resembles a point source in the wave-optics regime and by increasing the source size, the fringes are smeared out and the magnification looks like the geometric optics.
The distance between the fringes for a point-lens is $\Delta \beta={2\pi}/{f}$ \citep{Falco:1992} and for a typical extended source with $\rho>\Delta \beta$, fringes smear out due to integration over a highly oscillating function. Hence, the fringes are observable only for  the extended source that satisfies the condition of $\rho<\Delta \beta$. Summarizing this part,  in the wave optics regime fringes for an extended source can be produced under the following condition of
  \begin{equation}\label{eq:sour-size1}
  \theta_{\rm s} <\frac{1}{2}\frac{\lambda}{R_{\rm sch}}\theta_{E},\;
  \end{equation}
  \begin{figure}
 	\centering
 		\includegraphics[width=.45\textwidth]{finite}
\caption{Magnification of a uniform luminous source as a function of $\beta$ for $f=10$. The solid red, dot-dashed blue and dashed green lines show the magnification for $\rho=0.01$, $\rho = 0.08$ and $\rho=0.15$, respectively. By increasing the size of source, the incoherent light of the source results in magnification pattern converge to the geometric optics limit profile.}
\label{fig:mag-finite}
 \end{figure}
where taking into account the redshift of deflector at $z_d$, this condition can be written as
\begin{equation}\label{eq:sour-size2}
	\theta_{E}>2\theta_s\left(1+z_d\right){R_{\rm sch}}/{\lambda_{obs}},\;
\end{equation}
  where $\lambda_{obs}$ is the wavelength of observation.
Also, since $\theta_E\propto\sqrt{M}$ and $R_{\rm sch}\propto M$, the detection of fringes is in favour of small mass PBHs.

 Now, let us assume the lens mass to be in the range of $10^{-6}\lesssim \bar{M}\lesssim 10^{-2}$.
Then we rewrite the wave optics parameter of microlensing $f$, as follows
\begin{equation}\label{eq:f-pbh}
f=4\pi(1+z_{\rm d})\frac{R_{\rm sch}}{\lambda_{obs}}=4\pi(1+z_{\rm d})(\frac{\lambda_{obs}}{0.1\text{mm}})^{-1}{\bar{M}}.
\end{equation}
In the case of strong lensing of a quasar by a galaxy, it is more likely that PBH resides in the halo of the lensed galaxy, which allows us to measure the redshift of lens and from
equation (\ref{eq:f-pbh}) directly obtain the mass of PBH. We note that we had also another method
of mass measurement from equation (\ref{eq:mag-rel2}),  if we use at least two different wavelengths for the observation.

There are other observables in the geometric optics that can be used for breaking the degeneracy between the lens parameters. Let us take the finite size effect in the geometric optics which smoothes the peak of a light curve \citep{2011AJ....141..105P,Paris:2016xdm,2017PKAS...32..305S}.
Knowing the physical size of a quasar as a source from the astrophysical informations, from the finite-size effect, we can extract the projected Einstein radius on the source plane (i.e. $R_E^{(s)}=\theta_{E}D_s$) as
 \begin{equation}
 \label{eq1}
 R_E^{(s)} = 1.65 \text{A.U.}(\frac{D_{\rm s}}{6\text{Gpc}})^{1/2}(\frac{1-x}{x})^{1/2}{\bar M}^{1/2},
 \end{equation}
 where in this case we assume a source at redshift $z\approx 3$. We can combine equation (\ref{eq1}) with the parameter $f = 2kR_{\rm sch}$ from the wave optics observations to obtain the distance of the lens as well as its mass.

\section{Observational prospect and optical depth}
\label{obs}
From observational point of view, the cadence between the data points in the light curve of a quasar should be small enough to measure the oscillations due to the diffraction pattern in the light curve. In order to estimate the time scale between the fringes, we use
 $\Delta \beta$ where in terms of $f$ in equation (\ref{eq:f-pbh}) is given by
\begin{equation}\label{eq:mag-osc}
\Delta \beta = \frac{2\pi}{f}=\frac12\frac{\lambda _{obs}}{0.1\text{mm}}\frac{1}{1+z_d}\frac{1}{\bar M}.\;
\end{equation}
 By multiplying $\Delta \beta$ to the Einstein crossing time of lens, $t_{E}$, the time-scale for the transit of fringes can be obtained. The Einstein crossing time for typical parameters of a lens with lunar mass at the cosmological distances is
 \begin{equation}\label{eq:te-geo}
 t_{E} = 5.7{\text d}~(\frac{D_{\rm s}}{6\text{Gpc}})^{1/2}(\frac{1-x}{x})^{1/2}{\bar M}^{1/2}(\frac{500}{v_t}),
 \end{equation}
where $v_t$ is the relative transverse velocity of the lens-source-observer, which is $\sim 1000$ km/s for a rich cluster and  $\sim 200$ km/s in the galactic scales. Hence, the time-scale for transit of fringes (i.e. $\Delta t = t_E \Delta\beta $)  is given by

\begin{equation}\label{eq:dely-time}
\Delta t =  2.9\text{d}~ (\frac{\lambda_{obs}}{0.1\text{mm}})\frac{1}{1+z_d}(\frac{D_{\rm s}}{6\text{Gpc}})^{1/2}(\frac{1-x}{x})^{1/2}(\frac{500}{v_t}){\bar M}^{-1/2}.
\end{equation}

 For a lunar mass PBH located at the distance $z_d\sim1$ and wavelength of $\lambda_{obs}=100\mu m$, we have $\Delta\beta\sim 10^{-1}$ and the time scale of fringe-transit  is of the order of $\sim 1.5$ days. According to (\ref{eq:dely-time}), the transit time is proportional to the $\lambda$ and decreases for shorter wavelengths. For PBHs in the mass range of $10^{-6} <\bar{M}<10^{-2}$ the time scale of fringe-transit is within the range of $15 \text{d}<\Delta t< 4 \text{yr}$.

One of the important technical issues in the observations of quasars is the filtering of intrinsic variabilities compare to the diffraction signals. For a
quasar with the variability time scale shorter or in the same order of fringe transit time-scale, it is difficult to filter out the background signals.  Some of the quasars with very rapid variabilities have been detected in the timescales shorter than hours to minutes, so-called micro-variability \citep{Whiting:2001vj}.
One solution is to survey those quasars with the low variabilities.
The other possibility is to study the quasars in the strong lensing systems and remove any intrinsic variations in the light curve by shifting the light curves according to the time delay between the images \citep{2001ASPC..237..185W}. This method in recent years is used for detecting microlensing signals in the geometric optics regime \citep{Ricci2011,Ricci,Giannini}.

 \begin{figure}
  	\centering
  	\includegraphics[width=.48\textwidth]{optical_depth}
  	\caption{The wave optics optical depth as a function of PBH mass for sources with the size of 50 AU (solid blue line), 100AU (dashed red line) and 150AU (dot-dashed green line). Here we consider concordance $\Lambda$CDM with $\Omega_m=0.3$, $h=0.7$ and use $\lambda_{\rm obs}=100 \mu m$. While for a lunar mass PBH the optical depth is negligible, for the smaller PBH the optical depth is larger. In addition, the optical depth increases by decreasing the size of source.}
  	\label{fig:opt-dep}
  \end{figure}

    \begin{figure}
  	\centering
  	\includegraphics[width=.48\textwidth]{rate}
  	\caption{The number of expected events per year as a function of PBH mass. Here the wave length of the observation is adapted $\lambda_{\rm obs}=100 \mu m$ and the rate is plotted for three different size of sources as 50 AU (solid blue line), 100AU (dashed red line) and 150AU (dot-dashed green line).  We note that $\tau_{\text{w}}$ is larger for the smaller PBHs as $z_{max}$ is getting larger and the optical depth is larger for the smaller mass from Fig. (\ref{fig:opt-dep}). On the other hand $\Delta t$ also is getting larger to the smaller PBHs. So the ratio of these two terms in (\ref{eq:rate-geo}) is a function of mass and  has a peak for the rate of number of events as depicted in this figure.

	}
  	\label{fig:rate}
  \end{figure}

 In order to estimate the number of detectable events, we calculate the microlensing optical depth for detection of PBH. The optical depth  \citep{rahvar:2015,2007A&A...475..453Z} is defined by
  \begin{equation}\label{eq:opt-dep1}
  \tau=\int\pi R_E^2n(M,z)c\frac{dt}{dz}dz,\;
  \end{equation}
  where $n(M,z)$ is the number density of PBHs and it follows the spatial clustering of the cold dark matter. For the lower redshifts, the optical depth is a function of the direction of line of sight, depending on the cosmic density perturbations that cross the line of sight. However, for quasars at the higher redshifts, we can take almost a uniform number density of PBHs that is proportional to the dark matter density of the universe.  This assumption has been carefully investigated in \citep{Zackrisson}, where by considering 6 different models for halos and sub-halos, for quasars at higher redshifts (i.e. $z>0.25$) the optical depth from the clustered and uniform distribution of lenses converge. Here, we assume a uniform distribution for the density of PBHs in the optical depth calculation.

  For a uniform distribution of PBH, we define a new optical depth, so-called the wave optics optical depth by $\tau_{\text w}$ where in $\Lambda$CDM, it is given by
  \begin{equation}\label{eq:opt-dep2}
    \tau_{\text{w}}=\frac{3}{2}\frac{D_H}{D_s}\Omega_{pbh}\int_0^{z_{max}}\frac{(1+z)^2(D_s-D_l)D_ldz}{\sqrt{\Omega_m(1+z)^3+\Omega_{\Lambda}}},\;
    \end{equation}
  where $\Omega_{pbh}$ is the density parameter of PBHs and $D_H$ is the present horizon size of universe.  The difference between this definition with the conventional optical depth is that in this equation, $z_{max}$ is not assigned to the position of the quasar while that is the largest distance for a lens that satisfies the detection of the wave optics condition ($\rho<\Delta y$).  Moreover, in the geometric optics $\tau$ is independent of the mass function of the lenses and it depends on the overall mass density of the lenses. However, for the wave optics regime, the optical depth depends
on the mass of lens as well as $z_{max}$. In Fig.(\ref{fig:opt-dep}) we plot the optical depth in unite of $\Omega_{pbh}$ for three different values of the source sizes. In this plot, we consider the concordance cosmology model of $\Omega_{m}=0.3~,~h=0.7$ and put the source at the redshift $z=3$. As it is expected, the small mass PBH and small size sources are in favor of wave optics microlensing detection.

 For a lunar mass PBH, the optical depth is very small, however it grows rapidly by decreasing the mass of lenses. For the case of $\bar{M}\ll 1$, from equation (\ref{eq:dely-time}) and equation (\ref{eq:te-geo}), $t_E \ll \Delta t$ and we take $\Delta t$ as the corresponding time-scale for the microlensing events in the wave optics regime instead of $t_E$.

 Now we define the rate of events in the regime of wave optics microlensing as
 \begin{equation}\label{eq:rate-geo}
 \Gamma_{\text w}=\frac{2}{\pi}\frac{\tau_{\text{w}}}{\Delta t},\;
 \end{equation}
where both $\tau_{\text{w}}$ and $\Delta t $ depend on the mass of PBH. Assuming Dirac-Delta function for the mass function of PBHs, in Fig. (\ref{fig:rate}) the rate of events per year is depicted as a function of $\bar{M}$. In equation (\ref{eq:rate-geo}), the optical depth increases for the smaller masses. On the other hand $\Delta t$ also increase for the smaller masses with the factor of ${\bar M}^{-1/2}$, so the ratio of these two terms in  (\ref{eq:rate-geo}) results in a peak as depicted in Fig. (\ref{fig:rate}).
For a given quasar, the number of detectable events is $N_{obs}=\Gamma T_{obs}$, where $T_{obs}$ is the duration of observation. Let us take a quasar with the size of $50$ AU, then from Fig. (\ref{fig:rate}), we expect to detect $N_{obs} \simeq 0.9\Omega_{\rm pbh}/{\text yr}$ for PBHs with $\bar{M}\sim 10^{-3}$ and $N_{obs} = 0.3\Omega_{\rm pbh}/{\text yr}$ for the sources with the radius $r_s=100AU$. Now if hundred percent of the dark matter is made of PBHs  (i.e. $\Omega_{\rm pbh} = 0.3$), all with the mass of $\bar{M}\sim 10^{-3}$, we expect to detect $0.27$ and $0.09$ event per year, respectively. Also for this mass, equation (\ref{eq:dely-time}) provides $\Delta t \simeq 46$~d and a cadence in the observation of light curve with one day can reveal the oscillation mode of diffraction pattern.

Let us define the contribution of PBH on the density of dark matter as $f_{pbh} =  \Omega_{pbh}/\Omega_{m}$.
Then we also define the parameter of $df_{pbh}/d\bar{M}$ which provide the fraction of dark matter in form of PBHs within the range of $(\bar{M}, \bar{M} + d\bar{M})$. This function can be measured by a long-term survey of quasars with cadence of order of $\sim$ days.
 Fig.(\ref{fig:sampl-data}) demonstrates the simulation of data points for a microlensing event with the parameters of $t_{E}=10$~hr, $\rho=0.1$, $u_0=0.05$ and $\bar{M}=0.007$.  Here we adapt the wavelength of $\lambda=2\mu m$ which results in the transit time scale of fringes of order of a few hours.  We assumed a photometric signal to noise ratio of $S/N=50$ and recover the parameters of light curve, using the Markov chain Monte Carlo method.
The best values of parameters with 1-$\sigma$ uncertainty for the light curve in Fig. (\ref{fig:sampl-data}) are given in Tab.(\ref{tab:res}).
The maximum magnification for this light curve is around $\sim 17$ and increases rapidly by decreasing $u_0$ and $\rho$. We note that the ratio of light in the anti-nodes to the nodes  (where $A_{node}\rightarrow 0$) of the interference pattern in Fig. (\ref{fig:sampl-data}) is a large number,  much larger than the intrinsic variations of a quasar which is about  $\sim 50\%$ \citep{Soldi:2008ev}. So, once we have enough photometric accuracy and the cadence shorter than the interference crossing time-scale, we can detect our desired signals.

 \begin{figure}
 	\centering
 	\includegraphics[width=.48\textwidth]{light_curve}
 	\caption{Simulated data points and the best fit model for a typical event. Here we use $t_E=10$ hr, $\rho=0.1$, $u_0=0.05$, $\bar{M}=7\times 10^{-3}$ and $\lambda=2\mu m$ to simulate data points.
		}
 	\label{fig:sampl-data}
 \end{figure}

\begin{table}
 	\centering
\begin{tabular} { l  c}
	Parameter &  68\% limits\\
	\hline
	{\boldmath$t_E            $} & $10.035\pm 0.036           $\\

	{\boldmath$u_0            $} & $0.0557\pm 0.0033          $\\

	{\boldmath$\rho           $} & $0.09995\pm 0.00070        $\\

	{\boldmath$\bar{M}              $} & $0.007027\pm 0.000025      $\\
	\hline
\end{tabular}
  	\caption{The best fit parameters with 1-$\sigma$ confidence level recovered from the fitting to the light curve in Fig. (\ref{fig:sampl-data}), using the Markov chain Monte Carlo method.}
  	\label{tab:res}
  \end{table}

\section{Conclusion}\label{sec:conclude}

Summarizing this work, we have proposed a new method for the microlensing observation of quasars from the far infrared to the millimeter wavelengths. For the small mass lenses where the Schwarzschild radius of the lens is of the order of the wavelength of observation, the gravitational lens can produce distortions on the wavefront of the light in the order of one wavelength. Since the lens and the source are far enough from the observer, this situation is similar to the Huygens-Fresnel approximation and the result is the diffraction pattern from the phase-shifted electromagnetic wave on the lens plane. A relative motion of the observer through the diffraction pattern on the observer plane produces a modulation in the light curve of the quasar. One of the problems with this wave-optics microlensing observation would be the intrinsic variations of quasars that might be mixed with our desired signals. In order to solve this problem, we suggested the observation of quasars with the multiple images from the strong lensing.  The advantage of using these quasars is that by shifting the time delay between the images, we can remove the intrinsic variations of the quasar and extract our desired signals.

We suggested the observation of wave-optics microlensing in two different wavelengths. This technic enable us to measure directly the mass of a lens. Also with single wavelength observation, from the measuring the redshift of the strong lensing galaxy and the redshift of the source, we can break the degeneracy between the lens parameters and extract the mass of lenses. One of the possible candidates for the dark matter is the PBHs in the mass range smaller than the lunar mass. In this work, we proposed the observation of quasars with suitable cadence and photometric accuracy to observe the transit of the fringes in the diffraction pattern by the observer. The optical depth and the rate of events depend on the fraction of dark matter made of PBHs as well as the mass of the PBHs. Assuming the mass of PBHs in the order of $10^{-3}$ lunar mass and hundred percent of dark matter is made of PBHs, we obtained the rate of event detection per year for a given quasar within the range of $\sim 0.1$ to $\sim 0.3$. The wave-optics quasar microlensing can put a constraint on the fraction of dark matter made of PBHs as well as the mass function of PBHs. A long term survey of quasars by the infrared telescopes such as Spitzer space-based telescope or millimeter and submillimeter wavelength ground-based telescopes such as ALMA was suggested for this project.

 \acknowledgments
We thank David Spergel for the useful discussions. Also we thank anonymous
referee for his/her useful comments and suggestions. S. Rahvar was supported by Sharif
University of Technology's Office of Vice President for Research
under Grant No. G950214.

 \bibliographystyle{apsrev4-1}