\section{Introduction}\label{sec1}

Since its discovery in 1998, the current accelerated expansion
 of our universe~\cite{r1} has been one of the most active
 fields in modern cosmology. As is well known, it could be due
 to an unknown energy component (dark energy) or a modification
 to general relativity (modified gravity)~\cite{r1,r2}. Before
 1998, it was commonly believed that in the future our universe
 will either expand forever or contract again into a final Big
 Crunch. However, the discovery of the current accelerated
 expansion of the universe significantly changed our view of
 the fate of the universe. In fact, many novel possibilities
 are under the active consideration in the community nowadays.

Today, there are many dark energy candidates in the market.
 Among them, the famous phantom dark energy~\cite{r3,r4} is
 very interesting. Its equation-of-state parameter (EoS) is
 smaller than $-1$. Although phantom dark energy is consistent
 with the current observational data~\cite{r1,r3}, it violates
 all the energy conditions. One of the consequences is that
 our universe will encounter a singularity at a finite time,
 namely the so-called Big Rip~\cite{r4}. At this singularity,
 the scale factor $a$, energy density and pressure of our
 universe are all divergent. In fact, besides the traditional
 Big Bang, Big Crunch, and the Big Rip, many novel
 singularities have been considered in the literature, such as
 Sudden singularities, Generalized sudden singularities,
 Quiescent singularities, Big Boost, Big Brake, Big Freeze,
 $w$ singularities, Inaccessible singularities, Directional
 singularities (see e.g.~\cite{r5,r6} and references therein
 for some brief reviews). These singularities arise at the
 price of violation of one or several energy conditions.
 As is well known, in~\cite{r6} (see also e.g.~\cite{r5,r7,r8})
 the future singularities have been classified into four types,
 namely
 \begin{itemize}
   \item Type I (Big Rip): $a\to\infty$, $\rho\to\infty$,
    $H\to\infty$, $|p|\to\infty$, when $t\to t_s<\infty\,$;
   \item Type II (Sudden singularity): $a\to a_s$, $\rho\to\rho_s$,
    $H\to H_s$, $|p|\to\infty$, $\dot{H}\to\infty$, when
     $t\to t_s<\infty\,$;
   \item Type III (Big Freeze): $a\to a_s$, $\rho\to\infty$,
    $H\to\infty$, $|p|\to\infty$, when $t\to t_s<\infty\,$;
   \item Type IV (Generalized sudden singularity): $a\to a_s$,
    $\rho\to\rho_s$, $H\to H_s$, $|p|\to p_s$,
     $\dot{H}\to\dot{H}_s$, and higher derivatives
     of $H$ diverge, when $t\to t_s<\infty\,$,
 \end{itemize}
 where $t_s$, $a_s$, $\rho_s$, $H_s$, $p_s$, $\dot{H}_s$ are
 all finite constants ($a_s\not=0$); $\rho$ and $p$ are energy
 density and pressure respectively; $H\equiv\dot{a}/a$ is the
 Hubble parameter, and a dot denotes a derivative with respect
 to the cosmic time $t$. In fact, the above four types can
 include almost all known future singularities.

Of course, singularities usually are not desirable in physics.
 Therefore, other possible fates of our universe are also
 considered in the literature, such as the cyclic/oscillatory
 cosmology. Recently, some interesting scenarios concerning the
 fate of the universe attracted much attention in the
 community, namely the so-called ``Little Rip''~\cite{r9} and
 ``Pseudo-Rip''~\cite{r10}. In~\cite{r9,r10,r11}, their authors
 showed that if the cosmic energy density will remain constant
 or monotonically increase in the future, then all the possible
 fates of our universe can be divided into four categories
 based on the time asymptotics of the Hubble parameter
 $H(t)$~\cite{r10}, namely
 \begin{itemize}
   \item Big Rip: $H(t)\to\infty$, when $t\to t_{rip}<\infty\,$;
   \item Little Rip: $H(t)\to\infty$, when $t\to\infty\,$;
   \item Cosmological Constant: $H(t)=const.\,$;
   \item Pseudo-Rip: $H(t)\to H_\infty <\infty$, when $t\to\infty\,$,
 \end{itemize}
 where $H_\infty$ is a constant. Obviously, the Big Rip singularity
 is not the only fate of our universe with the phantom-like
 dark energy. Both Little Rip and Pseudo-Rip are non-singular,
 and hence fall outside of the four categories in~\cite{r6}.
 Similar to the Big Rip, the Little Rip dissociates {\em all}
 bound structures, but the strength of dark energy is not
 enough to rip apart spacetime (unlike the Big Rip)~\cite{r9,r10}.
 On the other hand, the Pseudo-Rip dissociates the bound structures
 which are held together by a binding force at or below a
 particular threshold, and hence it is possible that only
 {\em some} bound structures are dissociated while the others
 are {\em not} dissociated (depending on the model
 parameters)~\cite{r10}. In fact, the Little Rip is an intermediate
 case between the cosmological constant and the Big Rip~\cite{r9},
 while the Pseudo-Rip is an intermediate case between the
 cosmological constant and the Little Rip~\cite{r10}.

It is worth noting that all the Big Rip, Little Rip and Pseudo-Rip
 arise from the assumption that the dark energy density
 $\rho(a)$ is monotonically increasing~\cite{r9,r10,r11}, i.e., the
 dark energy is phantom-like (its EoS $w<-1$). In the present work,
 we are interested to investigate what will happen if this
 assumption is broken. Obviously, in the case of the dark
 energy density $\rho(a)$ is monotonically decreasing (i.e.,
 the dark energy is quintessence-like with an EoS $w>-1$), no
 rip will happen and no bound structures will be dissociated.
 On the other hand, in the case of the dark energy density $\rho(a)$
 monotonically decreases (namely $w>-1$) in the first stage
 and then monotonically increases (namely $w<-1$) in the second
 stage (this is the case of the so-called ``quintom A'' dark energy
 in terminology of e.g.~\cite{r12} and references therein), the fate
 of our universe is the Big Rip, which is trivial in some sense. The
 third case is that the dark energy density $\rho(a)$ monotonically
 increases (namely $w<-1$) in the first stage and then
 monotonically decreases (namely $w>-1$) in the second stage
 (this is the case of the so-called ``quintom~B'' dark
 energy in terminology of e.g.~\cite{r12} and references
 therein). It can be expected that in the first stage some or
 all bound structures will be dissociated (similar to the case
 of Pseudo-Rip), but then the disintegration process will stop, and
 the already disintegrated structures have the possibility to
 be recombined in the second stage. We dub it ``Quasi-Rip'',
 which is the subject of the present work (here we temporarily
 do not consider the case of oscillatory quintom dark energy).
 Since the quintom-like dark energy~\cite{r13} (whose EoS can
 cross the so-called phantom divide $w=-1$) is slightly favored
 by the observational data~\cite{r1} (see e.g.~\cite{r12} for a
 comprehensive review), we note that the Quasi-Rip is well-motivated
 in fact.

This paper is organized as follow. In Sec.~\ref{sec2}, we
 discuss the disintegration of bound structures. In Sec.~\ref{sec3},
 we present an explicit model of Quasi-Rip. We constrain this model
 with the observational data, and then clearly show the Quasi-Rip in
 this model. In Sec.~\ref{sec4}, some concluding remarks are given.

\section{The disintegration of bound structures}\label{sec2}

It is useful to introduce the concept of ``inertial force''
 when we discuss the disintegration of bound structures.
 According to~\cite{r10,r11,r23}, in a flat
 Friedmann-Robertson-Walker (FRW) universe dominated by dark energy,
 (motivated by the well-known Newton's second law) the inertial
 force $F_{inert}$ on a mass $m$ as seen by a gravitational
 source separated by a comoving distance $r_0$ is given by
 \be{eq1}
 F_{inert}\equiv mr_0\frac{\ddot{a}}{a}=mr_0\left(\dot{H}
 +H^2\right)=-mr_0\frac{4\pi G}{3}\left(\rho+3p\right)
 =mr_0\frac{4\pi G}{3}\left(2\rho+a\frac{d\rho}{da}\right),
 \end{equation}
 in which we have used the energy conservation equation
 $\dot{\rho}+3H(\rho+p)=0$. A bound structure dissociates when
 the inertial force $F_{inert}$ (dominated by dark energy) is
 equal to the force $F_{bound}$ holding together this bound
 structure~\cite{r10,r11}. For the gravitationally bound
 structure whose mass is $M$, it dissociates when
 \be{eq2}
 F_{inert}=mr_0\frac{\ddot{a}}{a}=F_{bound}=G\frac{Mm}{r_0^2}
 =mr_0\omega_0^2\,,
 \end{equation}
 namely
 \be{eq3}
 \frac{\ddot{a}}{a}=\omega_0^2=\frac{GM}{r_0^3}\,,
 \end{equation}
 or equivalently
 \be{eq4}
 f(a)\equiv\frac{1}{\rho_0}\left(2\rho+a\frac{d\rho}{da}\right)
 =\frac{1}{\rho_0}\left(2\rho+\frac{d\rho}{d\ln a}\right)
 =\frac{2H_0^{-2}}{\Omega_0}\frac{GM}{r_0^3}\,,
 \end{equation}
 where $\omega_0$ is the angular velocity; $\rho_0$ is the
 present density of dark energy;
 $\Omega_0\equiv (8\pi G\rho_0)/(3H_0^2)$ is the present
 fractional density of dark energy. Note that in this paper,
 the subscript ``0'' usually indicates the present value of
 corresponding quantity, and we have set $a_0=1$. In fact, one
 can check that the results Eqs.~(\ref{eq3}) and~(\ref{eq4})
 are coincident with the ones in the case of phantom dark
 energy whose EoS $w$ is a constant~\cite{r4} (we refer to
 e.g.~\cite{r14} for its detailed derivation; note that the
 result of~\cite{r4} is invalid for a non-constant $w$, while
 our results Eqs.~(\ref{eq3}) and~(\ref{eq4}) are still valid).

However, if the bound structure is massive enough to
 significantly affect the local spacetime metric, it is not
 accurate to express $F_{inert}$ in terms of FRW
 metric~\cite{r9,r10}. A more accurate method was presented
 in~\cite{r14}. Following~\cite{r14}, in the Newtonian limit,
 the interpolating metric of the local spacetime is given by
 \be{eq5}
 ds^2=\left[1-\frac{2GM}{a(t)\eta}\right]dt^2-a^2(t)\left[
 d\eta^2+\eta^2\left(d\theta^2+\sin^2\theta d\varphi^2\right)
 \right]\,,
 \end{equation}
 where $\eta$ is the comoving radial coordinate. The radial
 equation of motion for a test particle in the Newtonian limit
 reads
 \be{eq6}
 \ddot{r}=\frac{\ddot{a}}{a}r+\frac{L^2}{r^3}-\frac{GM}{r^2}\,,
 \end{equation}
 where $L=r^2\dot{\varphi}=r^2\omega=const.=r_0^2\omega_0$ is
 the constant angular momentum per unit mass. Noting that
 $\omega_0^2=GM/r_0^3\,$, Eq.~(\ref{eq6}) can be recast as
 \be{eq7}
 \ddot{r}=\frac{\ddot{a}}{a}r+\frac{GMr_0}{r^3}-\frac{GM}{r^2}
 =-\frac{dV_{\rm eff}}{dr}\,,
 \end{equation}
 where the time-dependent effective potential is given by
 \be{eq8}
 V_{\rm eff}=-\frac{1}{2}\frac{\ddot{a}}{a}r^2+
 \frac{GMr_0}{2r^2}-\frac{GM}{r}\,.
 \end{equation}
 The bound structure dissociates when the minimum of the
 time-dependent effective potential (including the centrifugal
 term) disappears. Note that the corresponding
 radius $r_{_{V_{min}}}$ at which the effective potential
 reaches its minimum is determined by
 \be{eq9}
 \frac{\ddot{a}}{a}r^4_{_{V_{min}}}-GMr_{_{V_{min}}}+GMr_0=0\,.
 \end{equation}
 One can find (with the help of e.g. Mathematica) that
 Eq.~(\ref{eq9}) has a real solution only for
 \be{eq10}
 \frac{\ddot{a}}{a}\leq\frac{27}{256}\frac{GM}{r_0^3}\,.
 \end{equation}
 Therefore, the minimum of the time-dependent effective
 potential disappears when
 \be{eq11}
 \frac{\ddot{a}}{a}=\frac{27}{256}\frac{GM}{r_0^3}\,,
 \end{equation}
 which is also the time when the bound structure dissociates.
 Similarly, we can also recast Eq.~(\ref{eq11}) as
 \be{eq12}
 f(a)\equiv\frac{1}{\rho_0}\left(2\rho+a\frac{d\rho}{da}\right)
 =\frac{1}{\rho_0}\left(2\rho+\frac{d\rho}{d\ln a}\right)
 =\frac{27}{128}\frac{H_0^{-2}}{\Omega_0}\frac{GM}{r_0^3}\,.
 \end{equation}
 As is shown in~\cite{r14}, the accurate disintegration time
 obtained here is usually earlier than the qualitative one
 obtained in~\cite{r4} (this point can also be seen by
 comparing Eqs.~(\ref{eq11}), (\ref{eq12}) with
 Eqs.~(\ref{eq3}), (\ref{eq4})). It is worth noting that most
 of the results in~\cite{r14} are valid for dark energy whose
 EoS $w$ is a constant, while our results obtained here are
 also valid for dynamical dark energy whose EoS $w$ is not a
 constant.

In the present work, we consider five bound structures, namely
 the Coma Cluster ($M=6\times 10^{48}\,{\rm g}$,
 $r_0=9\times 10^{24}\,{\rm cm}$), the Milky Way galaxy
 ($M=2\times 10^{45}\,{\rm g}$,
 $r_0=5\times 10^{22}\,{\rm cm}$), the Solar System
 ($M=2\times 10^{33}\,{\rm g}$,
 $r_0=7\times 10^{15}\,{\rm cm}$), the Earth
 ($M=6\times 10^{27}\,{\rm g}$,
 $r_0=6.4\times 10^8\,{\rm cm}$), and the Hydrogen atom. Note
 that the first three bound structures are suitable for
 Eqs.~(\ref{eq11}) and~(\ref{eq12}), while the Earth is
 suitable for Eqs.~(\ref{eq3}) and~(\ref{eq4}), since its
 surface material has no centrifugal term in the effective
 potential. Finally, the Hydrogen atom is also suitable for
 Eqs.~(\ref{eq3}) and~(\ref{eq4}), but in which the term $GM$
 should be replaced by a new term
 $q_e^2/(4\pi\varepsilon_0 m_e)$, since $F_{bound}$ is the
 electromagnetic force in this case. Note that
 $q_e=1.6\times 10^{-19}\,{\rm C}$, $m_e=0.511\,{\rm MeV}$
 ($1\,{\rm MeV}=1.7827\times 10^{-27}\,{\rm g}$),
 $1/(4\pi\varepsilon_0)=9\times 10^{18}\,{\rm g\,cm^3\,C^{-2}
 \,sec^{-2}}$, and $r_0=5.3\times 10^{-9}\,{\rm cm}$ for
 the Hydrogen atom. In Eqs.~(\ref{eq4}) and~(\ref{eq12}),
 $G=6.672\times 10^{-6}\,{\rm cm^3\,g^{-1}\,sec^{-2}}$,
 $H_0^{-1}=3.0856\times 10^{17}\,h^{-1}\,{\rm sec}$ (here $h$
 is the Hubble constant $H_0$ in units of
 $100\,{\rm km/sec/Mpc}$).  Note that $\Omega_0$ and $h$
 will be determined by the observational data (see below).

Since the current observational data are in the epoch $a<1$,
 we cannot ignore the contribution from pressureless matter
 in this stage, although it can be safely ignored in the above
 discussions when the universe is dominated by dark energy.
 From the Friedmann equation
 $3H^2=8\pi G\rho_{tot}=8\pi G(\rho+\rho_m)$ and
 $\rho_m=\rho_{m0}a^{-3}$, we have
 \be{eq13}
 E^2\equiv\frac{H^2}{H_0^2}=\Omega_0\frac{\rho}{\rho_0}+
 \Omega_{m0}a^{-3}=\Omega_0\frac{\rho}{\rho_0}
 +\left(1-\Omega_0\right)a^{-3}\,,
 \end{equation}
 where $\Omega_{m0}\equiv(8\pi G\rho_{m0})/(3H_0^2)$ is the
 present fractional density of pressureless matter. If
 $\rho(a)$ is given, we can use Eq.~(\ref{eq13}) to constrain
 the model parameters with the current observational data. On
 the other hand, from Eqs.~(\ref{eq12}) or~(\ref{eq4}), one can
 find the corresponding scale factor $a_\ast$ at which the
 bound structure dissociates. Then, noting $H=\dot{a}/a$, we
 can evaluate the disintegration time measuring from today
 ($a=1$), namely
 \be{eq14}
 t_\ast-t_0=H_0^{-1}\int_1^{a_\ast}\frac{da}{aE(a)}=H_0^{-1}
 \int_0^{\ln a_\ast}\frac{d\ln a}{E(\ln a)}\,,
 \end{equation}
 where $H_0^{-1}=9.7776\,h^{-1}\,{\rm Gyr}$.

\section{An explicit model of Quasi-Rip}\label{sec3}

Here, we consider an explicit model of Quasi-Rip. In some
 sense, constructing a suitable model is a smart task. There
 are various ways to this end. For example, one can specify
 the scale factor $a$ as a function of the cosmic time $t$
 (see e.g.~\cite{r5,r15,r16}), the pressure $p$ as a function
 of the energy density $\rho$ (see e.g.~\cite{r6,r17,r18}),
 the energy density $\rho$ as a function of the scale factor
 $a$ (see e.g.~\cite{r9,r10}), or the Hubble parameter $H$ as
 a function of the cosmic time $t$ (see e.g.~\cite{r11}). It
 is shown in~\cite{r9,r10,r11} that these ways are equivalent
 in fact. In our case, for convenience, we choose to specify
 the energy density $\rho$ as a function of the scale factor
 $a$, similar to the cases of Little Rip~\cite{r9} and
 Pseudo-Rip~\cite{r10}. As mentioned in Sec.~\ref{sec1}, here
 our task is to construct an explicit function $\rho(a)$,
 which monotonically increases (namely $w<-1$) in the first
 stage and then monotonically decreases (namely $w>-1$)
 in the second stage. A naive idea is to use a piecewise
 function, which is phantom-like in the first sector (namely
 $\rho(a)\propto a^{-3(1+w_1)}$ with a constant EoS $w_1<-1$)
 and is quintessence-like in the second sector (namely
 $\rho(a)\propto a^{-3(1+w_2)}$ with a constant EoS $w_2>-1$).
 Then, let us refine this naive idea with a smooth function
 $\rho(a)\propto a^{\mu(a)}$, in which $\mu(a<a_t)>0$ and
 $\mu(a>a_t)<0$, where $a_t$ is the transition point. Noting
 that it is more convenient to use $\ln a$ as the variable
 in Eqs.~(\ref{eq12}), (\ref{eq4}), and~(\ref{eq14}), as well
 as the fact $a=e^{\,\ln a}$, it is appropriate to write $\mu$
 as a function of $\ln a$ instead. Therefore, the simplest
 function $\rho(a)$ can be given by
 \be{eq15}
 \rho(a)=\rho_0 a^{\alpha-\beta\ln a}
 =\rho_0 e^{\,\ln a\;(\alpha-\beta\ln a)}\,,
 \end{equation}
 where $\alpha$ and $\beta$ are both constants. Let us have
 a closer observation. By requiring $d\rho/d\ln a=0$, we find
 that its extremum locates at $\ln a=\alpha/(2\beta)$. To
 ensure this extremum is a maximum, $d^2 \rho/d(\ln a)^2$
 should be negative here, therefore $\beta>0$ is required.

Substituting Eq.~(\ref{eq15}) into Eq.~(\ref{eq13}), we have
 \begin{eqnarray}
 &E^2=\displaystyle\frac{H^2}{H_0^2}&=\Omega_0\,a^{\alpha-\beta\ln a}
 +\left(1-\Omega_0\right)a^{-3}=\Omega_0\,e^{\,\ln a\;(\alpha
 -\beta\ln a)}+\left(1-\Omega_0\right)e^{-3\ln a}\nonumber\\
 & &=\Omega_0 (1+z)^{-\alpha-\beta\ln(1+z)}+
 \left(1-\Omega_0\right)(1+z)^3\,,\label{eq16}
 \end{eqnarray}
 where $z=1/a-1$ is the redshift. There are three free model
 parameters, namely $\Omega_0$, $\alpha$ and $\beta$. With
 this $E(z)$, following the methodology in
 e.g.~\cite{r19,r20,r21}, we can constrain this model with the
 latest Union2.1 Type~Ia supernovae (SNIa) dataset~\cite{r22}
 which consists of 580 SNIa. In fact, we find that the best fit
 has $\chi^2_{min}=562.224$, and the best-fit parameters are
 $\Omega_0=0.718859$, $\alpha=0.0294816$ and $\beta=0.0002525$,
 whereas the corresponding $h=0.698862$. In Fig.~\ref{fig1},
 we present the corresponding $68.3\%$ and $95.4\%$ confidence
 level~(C.L.) contours in the $\alpha-\beta$ parameter space.
 Although their best-fit values are small, $\alpha$~and $\beta$ can
 still be consistent with the observational data in the large
 $2\sigma$~C.L. region, i.e., $-1.2\, \lower .75ex \hbox{$\sim$} \llap{\raise .27ex \hbox{$<$}} \,\alpha\, \lower .75ex \hbox{$\sim$} \llap{\raise .27ex \hbox{$<$}} \,0.5$
 and $0\leq\beta\, \lower .75ex \hbox{$\sim$} \llap{\raise .27ex \hbox{$<$}} \,3.2$.

 \begin{center}
 \begin{figure}[tbhp]
 \centering
 \includegraphics[width=0.49\textwidth]{fitsn.eps}
 \caption{\label{fig1} The $68.3\%$ and $95.4\%$ confidence
 level contours in the $\alpha-\beta$ parameter space. The
 best-fit parameters are also indicated by a black solid
 point.}
 \end{figure}
 \end{center}

\vspace{-4mm}

Next, we show the Quasi-Rip in this model clearly. Substituting
 Eq.~(\ref{eq15}) into the reduced inertial force $f(a)$
 defined in Eqs.~(\ref{eq12}) and~(\ref{eq4}), we have
 \be{eq17}
 f(a)=\frac{\rho}{\rho_0}\left(2+\alpha-2\beta\ln a\right)
 =\left(2+\alpha-2\beta\ln a\right)a^{\alpha-\beta\ln a}
 =\left(2+\alpha-2\beta\ln a\right)
 e^{\,\ln a\;(\alpha-\beta\ln a)}\,.
 \end{equation}
 In Fig.~\ref{fig2}, we plot $\ln\,(\rho/\rho_0)$ and $\ln f$
 as functions of $\ln a$ for various model parameters. For
 convenience, we fix $\beta=0.0001$ which is at the same order
 of its best-fit value, while the various values of $\alpha$
 are taken from the $1\sigma$~C.L. region in Fig.~\ref{fig1}.
 From the left panel of Fig.~\ref{fig2}, it is easy to see that
 our $\rho(a)$ in Eq.~(\ref{eq15}) is the desirable one, which
 monotonically increases in the first stage and
 then monotonically decreases in the second stage. It is not
 surprising that the plot of $\ln f$ is very similar to the
 one of $\ln\,(\rho/\rho_0)$, since their difference is given
 by $\ln\,(2+\alpha-2\beta\ln a)$ which is a small quantity in
 fact (cf. Eq.~(\ref{eq17})). Depending on the model parameters
 $\alpha$ and $\beta$, the inertial force can successively
 disintegrate the Coma Cluster, the Milky Way galaxy, the Solar
 System, the Earth and the Hydrogen atom. In different cases,
 some bound structures will be disintegrated while the other
 bound structures will not. If none of the bound structures is
 disintegrated, the corresponding Quasi-Rip is a failed rip.
 From the right panel of Fig.~\ref{fig2}, we find that the
 innermost one is a failed rip, the outermost one is the
 strongest Quasi-Rip which can disintegrate all the five bound
 structures, and the others are moderate ones which can only
 disintegrate one or several of the five bound structures.
 From Eqs.~(\ref{eq12}) or~(\ref{eq4}), we can find the
 corresponding $\ln a_\ast$ at which the bound structure
 dissociates. Then, we can evaluate the disintegration time
 measuring from today $t_\ast-t_0$ by using Eq.~(\ref{eq14})
 with $E(\ln a)$ given in Eq.~(\ref{eq16}), and present the
 results in Tab.~\ref{tab1}. For a fixed $\beta$, a bound
 structure dissociates earlier with a larger $\alpha$, since
 the corresponding inertial force is stronger. In the case of
 $\alpha=0.27$, the difference between the disintegration time
 of Earth and Solar System is 4349 years, while the difference
 between the disintegration time of Hydrogen atom and Earth is
 only 15 years. Finally, it is easy to see that the most
 distinct feature of Quasi-Rip is that the inertial force
 monotonically decreases in the second stage. Eventually, it
 will become lower than all the thresholds to disintegrate
 the bound structure. Therefore, the already disintegrated
 structures have the possibility to be recombined in the second
 stage. This is the unique feature of Quasi-Rip different from
 Big Rip, Little Rip and Pseudo-Rip. Our universe has a chance
 to be rebuilt from the ashes after the terrible rip. This
 might be the last hope in the ``hopeless'' rip.

 \begin{center}
 \begin{figure}[tbhp]
 \centering
 \includegraphics[width=0.98\textwidth]{rip.eps}
 \caption{\label{fig2} $\ln\,(\rho/\rho_0)$ and $\ln f$ as
 functions of $\ln a$ for various model parameters. The curves
 from the innermost to the outermost are for $\alpha=0.04$,
 $\alpha=0.064$, $\alpha=0.1$, $\alpha=0.15$, $\alpha=0.2$
 and $\alpha=0.27$, respectively, while $\beta$ is fixed
 to be 0.0001. The dashed horizontal lines in the right
 panel indicate the corresponding thresholds to disintegrate
 the Coma Cluster, the Milky Way galaxy, the Solar System,
 the Earth and the Hydrogen atom (from bottom to top). Note
 that $\Omega_0$ and $h$ are taken to be their best-fit
 values.}
 \end{figure}
 \end{center}

\vspace{-7mm}

 \begin{table}[tbhp]
 \begin{center}
 \begin{tabular}{l|llllll}
 \hline\hline\\[-4.1mm]
 Bound structure~~ & ~~$\alpha=0.04$~~ & $\alpha=0.064$~~
 & $\alpha=0.1$~~~ & $\alpha=0.15$~~ & $\alpha=0.2$
 & $\alpha=0.27$ \\[0.5mm]
 \hline\\[-4mm]
 Coma Cluster & & 563.154 & 333.555 & 216.991
  & 161.28 & 118.715\\[1mm]
 Milky Way & & & 344.035 & 223.126 & 165.75 & 122.005\\[1mm]
 Solar System & & & & 223.281 & 165.857985136 & 122.083461362\\[1mm]
 Earth & & & & & 165.857991583~~~ & 122.083465711\\[1mm]
 Hydrogen atom & & & & & & 122.083465696~\\[1mm]
 \hline\hline
 \end{tabular}
 \end{center}
 \caption{\label{tab1} The disintegration time measuring from
 today $t_\ast-t_0$ (in units of Gyr) for models with various
 $\alpha$, while $\beta$ is fixed to be 0.0001, and $\Omega_0$
 and $h$ are taken to be their best-fit values.}
 \end{table}

\vspace{-4mm}

\section{Concluding remarks}\label{sec4}

The fate of our universe is an unceasing topic of cosmology
 and the human being. The discovery of the current accelerated
 expansion of the universe significantly changed our view of
 the fate of the universe. Recently, some interesting scenarios
 concerning the fate of the universe attracted much attention
 in the community, namely the so-called ``Little Rip'' and
 ``Pseudo-Rip''. It is worth noting that all the Big Rip,
 Little Rip and Pseudo-Rip arise from the assumption that the
 dark energy density $\rho(a)$ is monotonically increasing. In
 the present work, we are interested to investigate what will
 happen if this assumption is broken, and then propose a
 so-called ``Quasi-Rip'' scenario, which is driven by a type
 of quintom dark energy. In this work, we consider an explicit
 model of Quasi-Rip in detail. We show that Quasi-Rip has an
 unique feature different from Big Rip, Little Rip and
 Pseudo-Rip. Our universe has a chance to be rebuilt from the
 ashes after the terrible rip. This might be the last hope in
 the ``hopeless'' rip.

Some remarks are in order. Firstly, as is shown in
 Sec.~\ref{sec3}, our Quasi-Rip model is well consistent with
 the current observational data. However, even in the $1\sigma$
 C.L. region of $\alpha-\beta$ parameter space, the future behavior
 of our universe can be different enough (depending on the
 particular model parameters $\alpha$ and~$\beta$). In fact,
 as is well known, the current observational data can be
 consistent with all the phantom-like, quintessence-like and
 quintom-like dark energy models. Therefore, the current
 observational data cannot tightly tell what is the true fate
 of our universe. Most of the possibilities (including Big
 Rip, Little Rip, Pseudo-Rip, Quasi-Rip, de~Sitter expansion,
 other future singularities and so on) are still living.

Secondly, the explicit model of Quasi-Rip considered in the
 present work is the simplest case. One can construct other
 more complicated $\rho(a)$ to implement the Quasi-Rip. For
 example, one might construct an EoS $w(a)$ as a function
 of scale factor $a$, which is smaller than $-1$ when $a<a_t$
 and larger than $-1$ when $a>a_t$. Then the corresponding
 $\rho(a)$ can be found from the energy conservation equation
 $\dot{\rho}+3H\rho(1+w(a))=0$. Of course, other smart methods
 to construct the desirable $\rho(a)$ are awaiting us.

Thirdly, as mentioned in~\cite{r10}, in its second Pseudo-Rip
 model, the reduced inertial force $f(a)$ can also have a peak,
 similar to our Quasi-Rip model. However, we note that in the
 Pseudo-Rip model, after the peak, the reduced inertial force
 $f(a)\to const.$ which is still higher than the corresponding
 threshold to disintegrate the bound structure. Therefore, the
 already disintegrated structures have {\em no} possibility
 to be recombined in the Pseudo-Rip models. On the contrary, in
 the Quasi-Rip models, the reduced inertial force $f(a)$
 monotonically decreases in the second stage. Eventually, it
 will become lower than all the thresholds to disintegrate the bound
 structure. Therefore, the already disintegrated structures
 have the possibility to be recombined in the second stage.

Fourthly, it is well known that phantom is unstable at quantum
 level and hence the perturbations grow large. Noting that
 in the present work our discussions are at classical level
 instead, this problem could be set aside. In fact, the quantum
 stability of a phantom phase has been considered in~\cite{r24}. The
 authors of~\cite{r24} studied the perturbations in the
 quantum-corrected effective field equation at one- and
 two-loop order, and they found that the system is stable. On
 the other hand, it is claimed in~\cite{r25} that scalar
 perturbations can grow during a phantom phase if EoS $w<-5/3$.
 However, from Eq.~(\ref{eq15}) and Fig.~\ref{fig1}, it is easy
 to see that the corresponding $w$ is only slightly smaller
 than $-1$ in the phantom phase for our particular $\alpha$ and
 $\beta$ (n.b. Fig.~\ref{fig1}). So, $w>-5/3$ instead and hence our
 Quasi-Rip model can avoid the problem raised in~\cite{r25}.
 Further, in fact the dark energy considered in this work is
 not necessarily a scalar field. It could even be an
 effective dark energy from modified gravity, namely the
 so-called ``geometric dark energy''. So, it might avoid the
 corresponding problems in the phantom phase.

Finally, in the present work, we consider only the quintom
 dark energy which crosses the phantom divide $w=-1$ once.
 In fact, we can also consider the quintom dark energy which
 can cross the phantom divide for many times. The most
 attractive Quasi-Rip model might be the one driven by the
 oscillatory quintom dark energy. In this oscillatory Quasi-Rip
 model, our universe will be destroyed and then be rebuilt
 again and again.

\section*{ACKNOWLEDGEMENTS}
We thank the anonymous referee for quite useful comments and
 suggestions, which helped us to improve this work. We are
 grateful to Professors Rong-Gen~Cai and Shuang~Nan~Zhang
 for helpful discussions. We also thank Minzi~Feng for kind
 help and discussions. This work was supported in part by NSFC
 under Grants No.~11175016 and No.~10905005, as well as NCET
 under Grant No.~NCET-11-0790, and the Fundamental Research
 Fund of Beijing Institute of Technology.

\renewcommand{\baselinestretch}{1.1}


\section{Introduction}\label{intro}

Cosmological perturbations from primordial inflation
\cite{Starobinsky:1979ty,Mukhanov:1981xt,Starobinsky:1982ee} have a
crucial importance for fundamental theory because they are first
quantum gravitational effect ever detected \cite{Mukhanov:2007zz,
Woodard:2009ns}, because they provide strong evidence for primordial
inflation \cite{Mukhanov:1990me,Liddle:1993fq,Lidsey:1995np}, and
because they set the initial conditions for structure formation in
cosmology \cite{Dodelson:2003ft,Mukhanov:2005sc,Weinberg:2008zzc}.
Two sorts of perturbations are produced during single-scalar
inflation: a scalar perturbation characterized by the field
$\zeta(t,\vec{x})$, and a tensor perturbation characterized by the
transverse-traceless field $h_{ij}(t,\vec{x})$. The scalar signal
has been imaged in the anisotropies of the cosmic microwave
background \cite{Komatsu:2010fb,Keisler:2011aw}, and by measuring
the matter power spectrum with large scale structure surveys
\cite{Reid:2009xm}. The tensor signal has not been imaged so far
\cite{Komatsu:2010fb,Keisler:2011aw} but strenuous efforts are
underway to detect it through the polarization of the cosmic
microwave background \cite{PLANCK,EBEX,SPIDER,PIPER}.

The two perturbation fields are defined on the homogeneous,
isotropic and spatially flat geometry characterized by scale factor
$a(t)$, with Hubble parameter $H(t)$ and slow roll parameters
$\epsilon(t)$ and $\eta(t)$,
\begin{equation}
ds^2 = -dt^2 + a^2(t) d\vec{x} \!\cdot\! d\vec{x} \; , \; H(t)
\equiv \frac{\dot{a}}{a} \; , \; \epsilon(t) \equiv
-\frac{\dot{H}}{H^2} \; , \; \eta(t) \equiv \epsilon -
\frac{\dot{\epsilon}}{2 H \epsilon} \; . \label{geometry}
\end{equation}
Only the quadratic parts of their Lagrangians are relevant to
current measurements (see \cite{Kahya:2010xh} for an explanation of
the full formalism),
\begin{eqnarray}
\mathcal{L}^{(2)}_{\zeta} & = & \frac{\epsilon a^3}{8 \pi G}
\Biggl\{ \dot{\zeta}^2 \!-\! \frac{\partial_k \zeta \partial_k
\zeta}{a^2} \Biggr\} \; , \label{Lzeta} \\
\mathcal{L}^{(2)}_{h} & = & \frac{a^3}{64 \pi G} \Biggl\{
\dot{h}_{ij} \dot{h}_{ij} \!-\! \frac{\partial_k h_{ij} \partial_k
h_{ij}}{a^2} \Biggr\} \; . \label{Lh}
\end{eqnarray}
The reported results
\cite{Komatsu:2010fb,Keisler:2011aw,Reid:2009xm} for the scalar and
tensor power spectra are consistent with evaluating the following
2-point correlators long after the first horizon crossing time $t_k$
such that $k = H(t_k) \, a(t_k)$,
\begin{eqnarray}
\Delta^2_{\mathcal{R}}(t,k) & \equiv & \frac{k^3}{2 \pi^2} \int \!\!
d^3x \, e^{-i \vec{k} \cdot \vec{x}} \Bigl\langle \Omega \Bigl\vert
\zeta(t,\vec{x}) \zeta(t,\vec{0}) \Bigr\vert \Omega \Bigr\rangle \;
, \label{DeltaR} \\
\Delta^2_{h}(t,k) & \equiv & \frac{k^3}{2 \pi^2} \int \!\! d^3x \,
e^{-i \vec{k} \cdot \vec{x}} \Bigl\langle \Omega \Bigl\vert
h_{ij}(t,\vec{x}) h_{ij}(t,\vec{0}) \Bigr\vert \Omega \Bigr\rangle
\; . \label{Deltah}
\end{eqnarray}

At tree order the correlators (\ref{DeltaR}-\ref{Deltah}) can be
expressed in terms of the scalar mode function $v(t,k)$ and its
tensor cousin $u(t,k)$,
\begin{eqnarray}
\Delta^2_{\mathcal{R}}(t,k) & = & \frac{k^3}{2 \pi^2} \times 4\pi G
\times \vert v(t,k) \vert^2 + O(G^2) \; , \label{treezeta} \\
\Delta^2_{h}(t,k) & = & \frac{k^3}{2 \pi^2} \times 32 \pi G \times 2
\times \vert u(t,k) \vert^2 + O(G^2) \; . \label{treeh}
\end{eqnarray}
The relevant mode equations and normalization conditions are,
\begin{eqnarray}
\ddot{v} + \Bigl(3 H + \frac{\dot{\epsilon}}{\epsilon}\Bigr) \dot{v}
+ \frac{k^2}{a^2} v = 0 \quad & , & \quad W_v \equiv v \dot{v}^* -
\dot{v} v^* = \frac{i}{\epsilon a^3} \; , \label{veqns} \\
\ddot{u} + 3 H \dot{u} + \frac{k^2}{a^2} u = 0 \quad & , & \quad W_u
\equiv u \dot{u}^* - \dot{u} u^* = \frac{i}{a^3} \; . \label{ueqns}
\end{eqnarray}
The initial conditions for the mode functions (which correspond to
the choice of Bunch-Davies vacuum) derive from the WKB solutions in
the far ultraviolet (that is, for $k \gg H(t) a(t)$),
\begin{eqnarray}
v(t,k) & \longrightarrow & \frac{\exp[-ik \int^t dt'/a(t')]}{\sqrt{2
k a^2(t) \epsilon(t)}} \Biggl\{ 1 + O\Bigl(\frac{H a}{k}\Bigr) \Biggr\}
\; , \label{initialzeta} \\
u(t,k) & \longrightarrow & \frac{\exp[-ik \int^t dt'/a(t')]}{\sqrt{2
k a^2(t)}} \Biggl\{1 + O\Bigl(\frac{Ha}{k}\Bigr) \Biggr\} \; .
\label{initialh}
\end{eqnarray}
As is evident from these asymptotic forms, both mode functions
oscillate and fall off in the far ultraviolet. They become
approximately constant near the time of first horizon crossing
$t_k$. (One can infer the existence of constant solutions quite
generally when the $\frac{k^2}{a^2(t)}$ terms in equations
(\ref{veqns}-\ref{ueqns}) become irrelevant.) It is these constant
amplitudes which determine the crucial theoretical predictions for
the power spectra through relations (\ref{treezeta}-\ref{treeh}).

One of the frustrating things about this formalism is the need to
employ approximation techniques, even for evaluating tree order
formulae such as (\ref{treezeta}-\ref{treeh}), because equations
(\ref{veqns}-\ref{ueqns}) for the mode functions cannot be solved
analytically for general scale factor $a(t)$
\cite{Tsamis:2002qk,Tsamis:2003zs}. Examples of such approximation
techniques include \cite{Stewart:1993bc,Wang:1997cw,Martin:1999wa}:
\begin{itemize}
\item{Assuming the slow roll parameter $\epsilon(t)$ is constant;}
\item{Matching the leading ultraviolet and infrared forms at $t = t_k$; and}
\item{Employing the full WKB solutions.}
\end{itemize}
When greater accuracy is needed, one must resort to numerical
evolution of (\ref{veqns}-\ref{ueqns}) from the known initial
conditions (\ref{initialzeta}-\ref{initialh}) until well past the
time of first horizon crossing \cite{Wang:1997cw}. Excellent
numerical solution techniques exist
\cite{Easther:2010qz,Mortonson:2010er,Easther:2011yq}, but a large
fraction of their power is wasted on reproducing the oscillations of
the mode functions, which contribute nothing to the power spectra.
This can be especially time-consuming for models in which there is
an extended phase of inflation, or when scanning properties of
classes of models.

A more effective technique was developed recently for the tensor
power spectrum \cite{Romania:2011ez} in order to work out the
gravitational wave signal from a novel model of inflation in which
the Hubble parameter oscillates for a brief time
\cite{Tsamis:2009ja,Romania:2010zq}. The key to the new technique is
to convert relations (\ref{ueqns}) into an equation for the
norm-squared of the mode function, $M(t,k) \equiv \vert u(t,k)
\vert^2$. This is the quantity that actually enters the tensor power
spectrum (\ref{treeh}), and it has a much more sedate evolution than
the mode function. Evolving $M(t,k)$ avoids the wasted effort of
numerically simulating the oscillations of $u(t,k)$. Further, an
excellent asymptotic expansion can be derived for $M(t,k)$ which
must be valid until only a few e-foldings before first horizon
crossing, no matter how long inflation persists
\cite{Romania:2011ez}. It is therefore only necessary to numerically
evolve $M(t,k)$ for a few e-foldings, from just before $t_k$ until
it becomes constant.

The purpose of this paper is to extend this more effective technique
to the scalar mode functions using a trick \cite{Tsamis:2003zs} for
converting solutions to $u(t,k)$ into solutions for $v(t,k)$. To
further simplify the formalism, we express our results in terms of
the tree order power spectra (\ref{treezeta}-\ref{treeh}) directly,
without ever mentioning the mode functions. Section \ref{eqn}
derives the key equations for $\Delta^2_{\mathcal{R}}(t,k)$ and
$\Delta^2_{h}(t,k)$, and section \ref{expand} gives their asymptotic
expansions. Our discussion comprises section \ref{discuss}.

\section{Differential Equations for $\Delta^2_{\mathcal{R}}$ and
$\Delta^2_{h}$}\label{eqn}

We begin by reviewing how one derives an equation for $M(t,k) \equiv
\vert u(t,k) \vert^2$. The first step is to write out the first and
second time derivatives,
\begin{eqnarray}
\dot{M}(t,k) & = & \dot{u} u^* + u \dot{u}^* \; , \label{Mdot} \\
\ddot{M}(t,k) & = & \ddot{u} u^* + 2 \dot{u} \dot{u}^* + u
\ddot{u}^* \; . \label{Mddot}
\end{eqnarray}
One next eliminates $\ddot{u}$ and $\ddot{u}^*$ in (\ref{Mddot}) by
using (\ref{ueqns}), and then recognizing factors of $M$ and
$\dot{M}$,
\begin{equation}
\ddot{M}(t,k) = -3 H \dot{M} - \frac{2 k^2}{a^2} M + 2 \dot{u}
\dot{u} \; . \label{Mstep2}
\end{equation}
The final factor involving $\dot{u} \dot{u}^*$ can be expressed in
terms of $M$ by subtracting the square of the Wronskian $W_u$ in
expression (\ref{ueqns}) from the square of $\dot{M}$ in expression
(\ref{Mdot}),
\begin{equation}
\dot{M}^2 - W^2_u = 2 M \times 2 \dot{u} \dot{u}^* = \dot{M}^2 +
\frac1{a^6} \qquad \Longrightarrow \qquad 2 \dot{u} \dot{u}^* =
\frac{\dot{M}^2}{2 M} + \frac1{2 M a^6} \; . \label{dotuu*}
\end{equation}
Substituting (\ref{dotuu*}) in relation (\ref{Mstep2}) gives the
desired equation for $M(t,k)$,
\begin{equation}
\ddot{M} + 3 H \dot{M} + \frac{2 k^2}{a^2} M = \frac{ \dot{M}^2}{2
M} + \frac1{2 M a^6} \; . \label{Meqn}
\end{equation}

Deriving an equation for $N(t,k) \equiv \vert v(t,k)\vert^2$ entails
the same first step,
\begin{eqnarray}
\dot{N}(t,k) & = & \dot{v} v^* + v \dot{v}^* \; , \label{Ndot} \\
\ddot{N}(t,k) & = & \ddot{v} v^* + 2 \dot{v} \dot{v}^* + v
\ddot{v}^* \; . \label{Nddot}
\end{eqnarray}
We next employ expression (\ref{veqns}) to eliminate $\ddot{v}$ and
$\ddot{v}^*$ in favor of $N$ and $\dot{N}$,
\begin{equation}
\ddot{N} = -\Bigl(3 H + \frac{\dot{\epsilon}}{\epsilon}\Bigr)
\dot{N} - \frac{2 k^2}{a^2} N + 2 \dot{v} \dot{v}^* \; .
\label{Nstep2}
\end{equation}
Subtracting the squares of the $\dot{N}$ from (\ref{Ndot}) and the
Wronskian $W_v$ in equation (\ref{veqns}) gives $\dot{v} \dot{v}^*$,
\begin{equation}
\dot{N}^2 - W^2_v = 2 N \times 2 \dot{v} \dot{v}^* = \dot{N}^2 +
\frac1{\epsilon^2 a^6} \qquad \Longrightarrow \qquad 2 \dot{v}
\dot{v}^* = \frac{\dot{N}^2}{2 N} + \frac1{2 N \epsilon^2 a^6} \; .
\label{dotvv*}
\end{equation}
And the final equation for $N(t,k)$ comes from substituting
(\ref{dotvv*}) in (\ref{Nstep2}),
\begin{equation}
\ddot{N} + \Bigl(3 H + \frac{\dot{\epsilon}}{\epsilon}\Bigr) \dot{N}
+ \frac{2 k^2}{a^2} N = \frac{\dot{N}^2}{2 N} + \frac1{2 N
\epsilon^2 a^6} \; . \label{Neqn}
\end{equation}

Because equations (\ref{Meqn}) and (\ref{Neqn}) are nonlinear, they
incorporate the Wronskian normalization as well as the evolution
equations for the mode functions. Converting them to equations for
$\Delta^2_{\mathcal{R}}(t,k)$ and $\Delta^2_{h}(t,k)$ requires only
changing the coefficients of the final terms,
\begin{eqnarray}
\ddot{\Delta}^2_{\mathcal{R}} + \Bigl(3 H +
\frac{\dot{\epsilon}}{\epsilon} \Bigr) \dot{\Delta}^2_{\mathcal{R}}
+ \frac{2 k^2}{a^2} \Delta^2_{\mathcal{R}} =
\frac{(\dot{\Delta}^2_{\mathcal{R}})^2}{\Delta^2_{\mathcal{R}}} +
\frac{2 G^2 k^6}{\pi^2 \epsilon^2 a^6}
\frac1{\Delta^2_{\mathcal{R}}} \; , \label{scalareqn} \\
\ddot{\Delta}^2_{h} + 3 H \dot{\Delta}^2_{h} + \frac{2 k^2}{a^2}
\Delta^2_{h} = \frac{(\dot{\Delta}^2_{h})^2}{\Delta^2_{h}} +
\frac{2^9 G^2 k^6}{\pi^2 a^6} \frac1{\Delta^2_{h}} \; .
\label{tensoreqn}
\end{eqnarray}

\section{Asymptotic Expansions for $\Delta^2_{\mathcal{R}}$ and
$\Delta^2_{h}$} \label{expand}

Equations (\ref{scalareqn})-\ref{tensoreqn}) require initial values
for the power spectra and their first time derivatives. Of course
these derive from the ultraviolet limits
(\ref{initialzeta}-\ref{initialh}), which imply the zeroth order
terms for an expansion in powers of $(H a/k)^2$,
\begin{eqnarray}
\Delta^2_{\mathcal{R}}(t,k) & \longrightarrow & \frac{G k^2}{\pi
\epsilon a^2} \Biggl\{ 1 + O\Bigl( \frac{H^2 a^2}{k^2}\Bigr)
\Biggr\} \; , \label{zeta0term} \\
\Delta^2_{h}(t,k) & \longrightarrow & \frac{16 G k^2}{\pi a^2}
\Biggl\{ 1 + O\Bigl( \frac{H^2 a^2}{k^2} \Bigr) \Biggr\} \; .
\label{h0term}
\end{eqnarray}
We shall develop the expansion for $\Delta^2_{h}(t,k)$ to second
order, and then employ a trick for converting tensor mode functions
to scalar ones \cite{Tsamis:2003zs} to obtain the corresponding
expansion of $\Delta^2_{\mathcal{R}}(t,k)$. These expansions are so
accurate that there seems little point to numerically evolving for
anything except the last few e-foldings before first horizon
crossing.

From the asymptotic limit (\ref{h0term}) one can see that the
``big'' terms in the evolution equation (\ref{tensoreqn}) are those
with no derivatives. It is best to segregate these terms,
\begin{equation}
\Delta^2_{h} - \Bigl( \frac{16 \pi G k^2}{\pi a^2} \Bigr)^2
\frac1{\Delta^2_{h}} = \frac{a^2}{2 k^2} \Biggl[-\ddot{\Delta}^2_{h}
- 3 H \dot{\Delta}^2_{h} + \frac{ (\dot{\Delta}^2_{h})^2}{2
\Delta^2_{h}} \Biggr] \; . \label{perteqn}
\end{equation}
The desired expansion is now easy to read off \cite{Romania:2011ez},
\begin{eqnarray}
\lefteqn{ \Delta^2_{h}(t,k) = \frac{16 G k^2}{\pi a^2} \Biggl\{ 1 +
\Bigl[1 - \frac12 \epsilon\Bigr] \Bigl( \frac{H a}{k}\Bigr)^2 }
\nonumber \\
& & \hspace{.5cm} + \Biggl[ \frac94 \epsilon \!-\! \frac{21}8
\epsilon^2 \!+\! \frac34 \epsilon^3 \!+\! \Bigl( \frac74 \!-\!
\frac34 \epsilon\Bigr) \frac{\dot{\epsilon}}{H} \!+\! \frac18
\frac{\ddot{\epsilon}}{H^2} \Biggr] \Bigl( \frac{H a}{k}\Bigr)^4 +
O\Bigl( \frac{H^6 a^6}{k^6} \Bigr) \Biggr\} \; . \qquad
\label{tensorexp}
\end{eqnarray}
Note that this expansion is in powers of $(\frac{H a}{k})^2$, rather
than the expansion in powers of $\frac{H a}{k}$ one gets for the
mode function. The better convergence for small $\frac{H a}{k}$ is
one more indication of the superiority of our method.

The trick for converting $\Delta^2_{h}(t,k)$ into
$\Delta^2_{\mathcal{R}}(t,k)$ involves simultaneously changing the
scale factor and the meaning of time \cite{Tsamis:2003zs},
\begin{eqnarray}
a(t) & \longrightarrow & \sqrt{\epsilon(t)} \times a(t) \; ,
\label{newa} \\
\frac{\partial}{\partial t} & \longrightarrow &
\frac1{\sqrt{\epsilon(t)}} \times \frac{\partial}{\partial t} \; .
\label{newdt}
\end{eqnarray}
One can see that these two replacements convert the $u(t,k)$ mode
relations (\ref{ueqns}) into the relations (\ref{veqns}) for
$v(t,k)$. Replacements (\ref{newa}-\ref{newdt}) also convert the
$\Delta^2_{h}(t,k)$ equation (\ref{tensoreqn}) into expression
(\ref{scalareqn}) for $\Delta^2_{\mathcal{R}}(t,k)$. We could have
used them to derive that equation had the straightforward derivation
not been so trivial.

To convert the asymptotic expansion (\ref{tensorexp}) into an
expansion for $\Delta^2_{\mathcal{R}}(t,k)$ we first note that the
replacements (\ref{newa}-\ref{newdt}) imply the following
geometrical replacements,
\begin{eqnarray}
H & \longrightarrow & \frac{H}{\sqrt{\epsilon}} \times \Bigl[ 1
\!+\! \epsilon \!-\! \eta \Bigr] \equiv \frac{H}{\sqrt{\epsilon}}
\times D \; , \label{newH} \\
\epsilon & \longrightarrow & \frac{ (1 \!+\! \epsilon \!-\! \eta) (2
\epsilon \!-\! \eta) \!-\! (\frac{\dot{\epsilon} - \dot{\eta}}{H})}{
(1 \!+\! \epsilon \!-\! \eta)^2 } \equiv \frac{N}{D^2} \; , \label{neweps} \\
\frac{\dot{\epsilon}}{H} & \longrightarrow & \frac{D \dot{N} \!-\! 2
N \dot{D}}{H D^4} \; , \label{newdoteps} \\
\frac{\ddot{\epsilon}}{H^2} & \longrightarrow & \frac{D^2 \ddot{N}
\!-\! 4 D \dot{N} \dot{D} \!-\! 2 D N \ddot{D} \!+\! 6 N
\dot{D}^2}{H^2 D^6} - \frac{ (\epsilon \!-\! \eta) (D \dot{N} \!-\!
2 N \dot{D})}{H D^5} \; . \label{newddoteps}
\end{eqnarray}
After many tedious manipulations the result is,
\begin{eqnarray}
\lefteqn{ \Delta^2_{\mathcal{R}}(t,k) = \frac{G k^2}{\pi \epsilon
a^2} \Biggl\{ 1 + \Biggl[ 1 \!-\! \frac12 (1 \!+\! \epsilon \!-\!
\eta) (2 \epsilon \!-\! \eta) - \Bigl( \frac{\dot{\epsilon} \!-\!
\dot{\eta}}{2 H} \Bigr) \Biggr] \Bigl( \frac{H a}{k} \Bigr)^2 }
\nonumber \\
& & \hspace{-.7cm} + \Biggl[ \frac38 (1 \!+\! \epsilon \!-\! \eta)
(3 \!-\! \epsilon \!-\! \eta) (2 \!-\! \eta) (2 \epsilon \!-\! \eta)
\!+\! \Bigl[ \frac54 \!+\! \frac92 \epsilon \!-\! \frac{17}4 \eta
\!-\! \frac78 \epsilon^2 \!-\! \frac32 \epsilon \eta \!+\!
\frac{13}8 \eta^2 \Bigr] \frac{\dot{\epsilon}}{H} \nonumber \\
& & \hspace{-.7cm} + \Bigl[ \frac12 \!-\! \frac{21}8 \epsilon \!+\!
\frac{13}8 \eta \!+\! \epsilon^2 \!+\! \frac12 \epsilon \eta \!-\!
\frac34 \eta^2 \Bigr] \frac{\dot{\eta}}{H} \!+\! \frac74 \Bigl(
\frac{\dot{\epsilon}}{H}\Bigr)^2 \!-\! \frac{25}8
\Bigl(\frac{\dot{\epsilon}}{H} \Bigr) \Bigl(
\frac{\dot{\eta}}{H}\Bigr) \!+\! \frac{11}8 \Bigl(
\frac{\dot{\eta}}{H} \Bigr)^2 \nonumber \\
& & \hspace{-.7cm} - \Bigl[ \frac32 \!+\! \frac38 \epsilon \!-\!
\frac34 \eta \Bigr] \frac{\ddot{\epsilon}}{H^2} \!+\! \Bigl[
\frac{13}8 \!+\! \frac{\epsilon}2 \!-\! \frac78 \eta\Bigr]
\frac{\ddot{\eta}}{H^2} \!-\! \Bigl( \frac{\dddot{\epsilon} \!-\!
\dddot{\eta}}{8 H^3} \Bigr) \! \Biggr] \! \Bigl( \frac{Ha}{k}
\Bigr)^4 \!\!\!+\! O\Bigl( \frac{H^6 a^6}{k^6} \Bigr) \! \Biggr\} .
\label{scalarexp} \quad
\end{eqnarray}
Unless $\epsilon(t)$ or $\eta(t)$ is changing rapidly on the Hubble
scale, the fractional error in (\ref{scalarexp}) should be of order
$(\frac{Ha}{k})^6$. Had we used this formula to estimate
$\Delta^2_{\mathcal{R}}(t,k)$ for $t$ as close as one e-folding
before first horizon crossing, the fractional error would only be
about $2.5 \times 10^{-3}$; for two e-foldings that would fall to
about $6.1 \times 10^{-6}$; and the fractional error would only be
about $1.5 \times 10^{-8}$ at three e-foldings before horizon
crossing. This is comparable to what loop corrections might give
\cite{Miao:2012xc}.

\section{Discussion}\label{discuss}

We have derived nonlinear, second order differential equations
(\ref{scalareqn}) and (\ref{tensoreqn}) for the correlators
$\Delta^2_{\mathcal{R}}(t,k)$ and $\Delta^2_{h}(t,k)$ whose late
time limits give the scalar and tensor power spectra at tree order.
Numerically evolving these correlators is more economical than
evolving the mode functions (\ref{veqns}-\ref{ueqns}) because one no
longer has to keep track of the phases which drop out of the power
spectra. This economy shows up in the asymptotic expansions
(\ref{scalarexp}) and (\ref{tensorexp}) pertinent to the ultraviolet
regime of $k \gg H(t) \, a(t)$. Our expansions for the correlators
are in powers of the small parameter $(\frac{Ha}{k})^2$, whereas the
analogous expansions for the mode functions are only in powers of
$\frac{Ha}{k}$. In fact our asymptotic expansions
(\ref{scalarexp}-\ref{tensorexp}) should be so accurate that they
could be used to provide the initial value data for
$\Delta^2_{\mathcal{R}}(t,k)$ and $\Delta^2_{h}(t,k)$ only a few
e-foldings before the time $t_k$ of first horizon crossing.

We have employed ``Hubble parametrization'', in which the scale
factor $a(t)$ is assumed known and results are expressed in terms of
background geometrical quantities $H(t)$, $\epsilon(t)$ and
$\eta(t)$ defined in expression (\ref{geometry}). This is especially
useful for the tensor power spectrum because gravitational waves
depend only on the background geometry, no matter what caused it.
Many people prefer ``potential parametrization'', in which the
scalar potential $V(\phi)$ is assumed known (with the field starting
at rest from some fiducial value) and results are expressed as
derivatives of the potential. Approximate conversion formulae are,
\begin{equation}
H^2 \simeq \frac83 \pi G V \quad , \quad \epsilon \simeq \frac1{16
\pi G} \Bigl( \frac{V'}{V} \Bigr)^2 \quad , \quad \eta \simeq
\frac1{16 \pi G} \Biggl[ 2\frac{V''}{V} \!-\! \Bigl( \frac{V'}{V}
\Bigr)^2 \Biggr] \; .
\end{equation}

Finally, we note that the same technology can be applied to derive
an equation for the power spectrum of any field whose plane wave
mode functions obey the Mukhanov equation \cite{Tsamis:2003zs}. Of
course the associated asymptotic expansion is also easy to develop.
Both the equation and its solution in the regime of $ k \gg H(t) \,
a(t)$ can be obtained from the analogous tensor relations,
(\ref{tensoreqn}) and (\ref{tensorexp}), by changing the scale
factor and the time as we did in equations (\ref{newa}-\ref{newdt}).

\vskip .5cm

\centerline{\bf Acknowledgements}

We are grateful to W. Kinney and M. Sasaki for conversations on this
subject. This work was partially supported by European Union program
Thalis ESF/NSRF 2007-2013, by European Union Grant
FP-7-REGPOT-2008-1-CreteHEPCosmo-228644, by NSF grant PHY-1205591,
and by the Institute for Fundamental Theory at the University of
Florida.


\section{Introduction}

Spatial capture-recapture (SCR) models are relatively new methods for
inference about population density from capture-recapture data using
auxiliary information about individual capture locations
\citep{efford:2004,borchers_efford:2008, royle_young:2008}.  SCR
models posit that $N$ individuals are located within a region denoted
$\mathcal{S}$. Each individual has a home range or activity area
within which movement occurs during some well-defined time interval,
and the center of the animal's activity has Cartesian coordinates
$\mathbf{s}_i$ for individuals $i=1,\ldots,N$. The population is sampled using $J$
traps with coordinates ${\bf x}_{j}$ for $j=1,\ldots,J$, and encounter
probability is expressed as a function of the distance between trap
location (${\bf x}_{j}$), and individual activity center (${\bf
  s}_{i}$).
While SCR models are a relatively recent innovation, their use is already becoming widespread
\citep{efford_etal:2009ecol, gardner_etal:2010jwm,
  gardner_etal:2010ecol,kery_etal:2010, gopalaswamy_etal:2012,
  foster_harmsen:2012} because they resolve critical problems with
ordinary non-spatial capture-recapture methods such as ill-defined
area sampled and heterogeneity in encounter probability due to the
juxtaposition of individuals with traps
\citep{borchers:2011}.
Furthermore, unlike traditional capture-recapture methods, SCR models
allow for inference about the processes determining
spatial variation in population density.

Despite the increasing popularity of SCR models, every
application of them has been based on encounter probability
models, such as the bivariate normal distribution, that imply
symmetric and stationary (invariant to translation) models for home range.
While such simple models might be necessitated in practice by sparse
data,
home range size and shape are often not well represented by stationary
distributions because animals select
resources
that
are unevenly distributed in space. Therefore
more complex models are needed to relate the capture
process with the way in which individuals utilize
space.

In this paper, we extend SCR capture probability models to accommodate
models of space usage or resource selection, by extending them to
include one or more explicit landscape covariates, which the
investigator believes might affect how individual animals use space
within their home range (this is what \citep{johnson:1980} called {\it
  third-order} selection). We do this in a way that is entirely
consistent with the manner in which parameters of classical resource
selection functions (RSF) \citep{manly_etal:2002} or utilization
distributions (UD) \citep{worton:1989, fieberg:2005, fieberg:2007} are
estimated from animal telemetry data.  In fact, we argue that SCR
models and RSF/UD models estimated from telemetry are based on the
same basic underlying model of space usage. The important distinctions
between SCR and RSF studies are that (1) resource selection studies do
not result in estimates of population density and (2) in SCR studies,
encounter of individuals is imperfect (i.e., ``$p<1$'') whereas, with
RSF data obtained by telemetry, encounter is perfect.  With respect to
the latter point, we can think of the RSF and SCR studies as being
exactly equivalent either if we have a dense array of trapping
devices, or if our telemetry apparatus
samples time or space imperfectly.
A key concept that we must confront in order to unify and integrate
SCR and RSF data is that we need to formulate both models in terms of
a common latent variable so that we can make them consistent with
respect to some underlying space utilization process. As we will
explain, this latent variable is the number of times that an
individual uses a particular region
of the landscape over some period of time.

\begin{comment}
We begin by describing a model for space usage that exists independent
of how we obtain the data (whether by telemetry or by
capture-recapture). The model posits that individuals use space in
proportion to one or more measurable covariates, say $z({\bf x})$,
defined for every pixel
${\bf x}$ in the landscape.  For example, $z({\bf x})$ could be the
habitat type or elevation at pixel $\textbf{x}$. Then we introduce
observation models consistent with both standard capture-recapture
methods and telemetry studies.  This allows us to define a general
likelihood function that is the product of the two components, if they
are both simultaneously available.  In the absence of RSF data from
telemetry, the model reduces to an ordinary SCR model but with a
spatial covariate $z({\bf x})$ affecting encounter probability.
Following \citet{borchers_efford:2008}, we provide a formal
likelihood-based analysis of the model, and carry-out a simulation
study for various sample sizes of telemetered individuals to evaluate
the improvement from adding telemetry information to SCR studies and
{\it vice versa}.
\end{comment}

The modeling framework we develop here simultaneously resolves three
important
problems: (1) it generalizes all existing capture probability models
for SCR data to accommodate realistic patterns of space usage that
result in asymmetric and irregular home
ranges;
(2) it allows estimation of RSF parameters directly from SCR data,
i.e., {\it absent} telemetry data; and (3) it provides the basis for
integrating telemetry data directly into SCR models to improve
estimates of model parameters, including density.  Our model greatly
expands the applied relevance of SCR methods for conservation and
management, and for addressing applied and theoretical questions
related to animal space usage and resource selection.

\section{Spatial Capture-Recapture}

A number of distinct observation models have been proposed for spatial
capture-recapture studies \citep{borchers_efford:2008,
  royle_etal:2009ecol, efford_etal:2009ecol}, including Poisson,
multinomial, and binomial observation models.
Here we focus on the binomial model in which we suppose that the $J$
traps
are operated for $K$ periods (e.g., nights), and the observations are
individual- and trap-specific counts $y_{ij}$, which are binomial with
sample size $K$ and capture probabilities $p_{ij}$ which depend on
trap locations ${\bf x}_{j}$ and individual activity centers ${\bf
  s}_{i}$ as described subsequently.  The vector of trap-specific
counts for individual $i$, ${\bf y}_{i} = (y_{i1},\ldots,y_{iJ})$ is
its {\it encounter history}.  A standard encounter probability model
\citep{borchers_efford:2008} is the Gaussian model in which
\begin{equation}
\log(p_{ij})= \alpha_{0} + \alpha_{1} d_{ij}^{2}
\label{eq.encounter}
\end{equation}
or, equivalently,
$p_{ij} = \lambda_{0} \exp(-d_{ij}^{2} /(2\sigma^{2}) )$,
where
$d_{ij}$ is the Euclidean
distance between points ${\bf s}_{i}$ and ${\bf x}_{j}$,
$d_{ij} = \|\textbf{s}_i - \textbf{x}_j\| = \sqrt{(s_{i1}-x_{j1})^2 +
  (s_{i2}-x_{j2})^2}$,
and
 $\alpha_{0} = log(\lambda_{0})$ and $\alpha_{1} = -1/(2\sigma^2)$.
Alternative detection models are used, but all are
functions of Euclidean distance and so we do not consider them
further here.
\begin{comment}
XXX Isn't all this from our ED ms?? It doesn't seem to belong here XXXX
XXXX yea -- paraphrased -- but I think we have to say how estimation
proceeds
XXXXX and also we have to say that the model implies symmetric home
ranges etc..
\end{comment}

The primary motivation behind our work is that, in all previous
applications of SCR models, simple encounter probability models based
only on Euclidean distance have been used, with estimation based on
standard likelihood or Bayesian methods. These methods regard the
activity center for each individual $i$, ${\bf s}_{i}$, as latent
variables and remove them from the likelihood either under a model of
``uniformity'' in which ${\bf s} \sim \mbox{Unif}({\cal S})$ where
${\cal S}$ is a spatial region (the ``state-space'' of ${\bf s}$), or
a model in which covariates might affect the spatial distribution of
individuals \citep{borchers_efford:2008}. The state-space ${\cal S}$
defines the potential values for any activity center ${\bf s}$, e.g.,
a polygon defining available habitat or range of the species under
study.

A critical problem with standard SCR models is that the encounter
probability model based on Euclidean distance metric is unaffected by
habitat or landscape structure, and it implies that the space used by
individuals is stationary and symmetric, which may be unreasonable in
many applications.  For example, if the common detection model based
on a bivariate normal probability distribution function is used, then
the implied space usage by {\it all} individuals, no matter their
location in space or local habitat conditions, is symmetric with
circular contours of usage intensity.  Subsequently we provide an
extension of this class of SCR models that accommodates asymmetric,
irregular and spatially heterogeneous models of space usage.  Thus,
``where'' an individual lives on the landscape, and the state of the
surrounding landscape, will determine the character of its usage of
space. In particular, we suggest encounter probability models that
imply irregular, asymmetric and non-stationary home ranges of
individuals and that are sensitive to the local landscape being used
by an individual.

\section{Basic Model of Space Usage}
\label{rsf.sec.rsfmodel}

We develop the model here in terms of a discrete landscape purely for
computational expediency. This formulation will accommodate the vast
majority of actual data sets, as almost all habitat or landscape
structure data comes to us in the form of raster data.  Let ${\bf
  x}_{1},\ldots,{\bf x}_{nG}$ identify the center coordinates of a set
of $nG$ pixels that define a landscape.  In SCR studies, a {\it
  subset} of the coordinates ${\bf x}$ will correspond to trap
locations where we might observe individuals whereas, in telemetry
studies, animals are observable (by telemetry fixes) at potentially
{\it all} coordinates.

Let $z({\bf x})$ denote a covariate measured (or defined) for every
pixel ${\bf x}$. For clarity, we develop the basic ideas here in terms of a
single covariate but, in practice, investigators typically have more
than 1 covariate, which poses no additional problems.
We suppose that a population of individuals wanders around space in
some manner related to the covariate $z({\bf x})$, and their locations
accumulate in pixels by some omnipotent accounting mechanism. We will
define ``use of ${\bf x}$'' to be the event that an individual animal
appeared in some pixel ${\bf x}$.
This is equivalently stated in the
literature in terms of individual having {\it selected} ${\bf  x}$.
 As a biological matter,
use is the outcome of individuals moving around their home range \citep{hooten_etal:2010},
i.e., where an individual is at any point in time is the result of
some movement process. However, to understand space usage, it is not
necessary to entertain explicit models of movement, just to observe
the outcomes, and so we don't elaborate further on what could be
sensible or useful models of movement.
\begin{comment}
, but we imagine existing methods
of hierarchical or state-space
models are suitable for this purpose \citep{jonsen_etal:2005,
  forester_etal:2007, patterson_etal:2008, hooten_etal:2010, mcclintock_etal:2012}.
\end{comment}

Suppose that an individual  is monitored over some period of time
and a fixed number, say $R$, of use observations are recorded.
Let $n({\bf x})$ be the use frequency of pixel ${\bf x}$ for that individual.
i.e., the number of times that individual used pixel ${\bf x}$
during some period of time.
We assume the
following probability distribution for the $nG \times 1$ vector of use
frequencies:
\[
{\bf n} \sim \mbox{Multinom}(R, {\bm \pi})
\]
where ${\bm \pi}$ is the $nG \times 1$ vector of use probabilities
with elements (for each pixel):
\[
 \pi({\bf x}) = \frac{ \exp( \alpha_{2} z({\bf x}) ) }
   { \sum_{x}    \exp(\alpha_{2} z({\bf x})) }
\]
This is the standard RSF model \citep{manly_etal:2002} used to model
telemetry data.
The parameter $\alpha_2$ is the effect of the
landscape covariate $z({\bf x})$ on the relative probability of
use. Thus, if $\alpha_2$ is positive, the relative probability of use
increases as the value of the covariate increases.
In practice, we don't get to
observe $\{ n({\bf x}) \}$ for all individuals but, instead, only for a small
subset say $i=1,2,\ldots, N_{tel}$,  which we capture and install telemetry devices on.
For the telemetered individuals, we assume they behave according
to the same RSF model as the population as a whole, which might be
justified if individuals are randomly sampled from the population.

We extend this model slightly to make it more realistic spatially and
also consistent with standard SCR models. Let
${\bf s}$ denote the centroid of an individuals home range and let
$d({\bf s},{\bf x}) = ||{\bf x} - {\bf s}||$ be the distance from the home
range center ${\bf s}$ of some individual to pixel ${\bf x}$, and let
$n({\bf x},{\bf s})$ denote the use frequency of pixel ${\bf x}$ for an individual with
activity center ${\bf s}$.
We modify the space usage model to accommodate that space
use will be concentrated around an individual's home range center \citep{johnson_etal:2008,forester_etal:2009}:
\begin{equation}
 \pi({\bf x}|{\bf s})  = \frac{ \exp( -\alpha_{1} d({\bf x},{\bf s})^{2} +\alpha_{2} z({\bf x}) )}
{ \sum_{x} \exp(-\alpha_{1} d({\bf x},{\bf s})^{2} +\alpha_{2} z({\bf x}))}
\label{rsf.eq.rsf}
\end{equation}
where $\alpha_1=1/(2\sigma^2)$ describes the rate at which encounter
probability declines as a function of distance, $d({\bf x},{\bf
  s})$. From ordinary telemetry data, it would be possible to estimate
parameters $\alpha_{1}$, $\alpha_{2}$ and also the activity centers
${\bf s}$ using standard likelihood methods based on the multinomial
likelihood \citep{johnson_etal:2008}.

Note that Eq.~\ref{rsf.eq.rsf} resembles standard encounter models
used in spatial capture-recapture but with an additional covariate
$z({\bf x})$.  The main difference between this observation model and
the standard SCR model is that the model here includes the normalizing
constant $\sum_{x} \exp(-\alpha_{1} d({\bf x},{\bf s})^{2} +\alpha_{2}
z({\bf x}))$, which ensures that the use distribution is a proper
probability density function. Thus we are able to characterize the
probability of encounter in terms of both distance from activity
center and space use.
Note that, under this model for space usage or resource selection, if
there are {\it no} covariates, or if $\alpha_{2} = 0$, then the
probabilities $\pi({\bf x}|{\bf s})$ are directly proportional to the
SCR model for encounter probability.  For example, setting $\alpha_{2}
= 0$, then this implies probability of use for pixel ${\bf x}$ is:
\[
p({\bf x}|{\bf s}) \propto  \exp( -\alpha_{1} d({\bf x},{\bf s})^{2}).
\]
Therefore,
for whatever model we choose for
$p({\bf x},{\bf s})$ in an ordinary SCR model, we can modify the distance
component in the RSF function in Eq. \ref{rsf.eq.rsf} accordingly to
be consistent with that model, by choosing $\pi({\bf x}|{\bf s})$
according to
\[
\pi({\bf x}|{\bf s}) \propto \exp( log(p({\bf x}|{\bf s})) + \alpha_{2} z({\bf x}))
\]

As an illustration of space usage patterns under this model, we
simulated a covariate that represents variation in habitat structure
(Fig. \ref{rsf.fig.habitat}) such as might correspond to habitat
quality.
This was simulated by using a
simple kriging interpolator  of spatial noise.
Space usage patterns for
 8 individuals in this landscape are shown in Fig. \ref{rsf.fig.homeranges},
simulated with $\alpha_{1} = 1/(2\sigma^2)$ with $\sigma = 2$ and the
coefficient on $z({\bf x})$ set to $\alpha_{2} = 1$.
These space usage densities -- ``home ranges'' -- exhibit clear
non-stationarity in response to the structure of the underlying
covariate, and they are distinctly asymmetrical.  We note that if
$\alpha_{2}$ were set to 0, the 8 home ranges shown here would
resemble bivariate normal kernels with $\sigma = 2$.  Another
interesting thing to note is that the activity centers are not
typically located in the pixel of highest use or even the centroid of
usage. That is, the observed ``average'' location is not an unbiased
estimator of ${\bf s}$ under the model in Eq. \ref{rsf.eq.rsf}.

\subsection{Poisson use model}

A natural way to motivate this specific model of space usage is to
assume that individuals make a sequence of random resource selection
decisions so that the outcomes $n({\bf x})$ (for all ${\bf
  x}$) are marginally {\it
  independent} Poisson random variables:
\[
 n({\bf x})|{\bf s} \sim \mbox{Poisson}( \lambda({\bf x}|{\bf s}))
\]
where
\[
 \log(\lambda({\bf x}|{\bf s})) = a_{0} -\alpha_{1} d({\bf x},{\bf s})^{2} +  \alpha_{2} z({\bf x})
\]
 In this case, the number of visits to any particular cell is affected
by the covariate $z({\bf x})$ but has a baseline rate ($\exp(a_{0})$)
related to the amount of movement occurring over some time interval.
This is an equivalent model to the multinomial
model given previously in the sense that, if we condition on the total
sample size $R = \sum_{x} n({\bf x})$, then the vector of use
frequencies $\{ n({\bf x}) \}$ for individual with activity center ${\bf s}$,
has a multinomial distribution with probabilities
\[
 \pi({\bf x}|{\bf s}) = \frac{\lambda({\bf x}|{\bf s})}{\sum_{x}
   \lambda({\bf x}|{\bf s})}
\]
which is the same as Eq. \ref{rsf.eq.rsf} because $a_{0}$ cancels from
the numerator and denominator of the multinomial cell probabilities
and thus this parameter is not relevant to understanding space usage.
Note that if use frequencies are summarized over
$i=1,2,\ldots,N_{tel}$ individuals for each pixel, then a standard
Poisson regression model for the resulting ``quadrat counts'' is
reasonable. This corresponds to ``Design 1'' in
\citet{manly_etal:2002}.

\subsection{Random Thinning}

Suppose our sampling is imperfect so that we only observe a smaller
number of telemetry fixes than actual use frequency, $n({\bf x})$.
We express this ``thinning'' (or sampling) by assuming the observed number of
uses is a binomial random variable based on a sample of size $n({\bf
  x})$:
\[
 m({\bf x}) \sim \mbox{Bin}(n({\bf x}), \phi_{0}).
\]

 Then, the marginal distribution of the new random variable $m$ is
 also Poisson
but with mean
\[
 log(\lambda({\bf x}|{\bf s})) = log(\phi_{0}) + a_{0} -\alpha_{1}
 d({\bf x}|{\bf s})^{2} +  \alpha_{2} z({\bf x}).
\]
 Thus, the space-usage model (RSF) for the
thinned counts $m$ is the same as the space-usage model for the
original variables $n$.  This is because if we remove $n$
from the conditional
 model by summing over its possible values, then the vector of thinned
 use frequencies
${\bf m}$ (i.e., for all pixels) is {\it also}  multinomial with cell probabilities
\[
\pi({\bf x}|{\bf s}) = \frac{\lambda({\bf x}|{\bf s})}{\sum_{x}  \lambda({\bf x}|{\bf s})}
\]
and so the constants $a_{0}$ and $\phi_{0}$
cancel from both the numerator and
denominator. Thus, the underlying RSF model applies to the true
unobserved count frequencies ${\bf n}$ and also those produced
by a random thinning or sampling process, ${\bf m}$.

In summary, if we conduct a telemetry study of $i=1,2\ldots,N_{tel}$
individuals, the observed data are
the $nG \times 1$  vectors of use frequencies ${\bf m}_{i}$ for each individual.
 We declare these data to be
``resource-selection data'' which are typical of the type used to
estimate resource-selection functions (RSFs) \citep{manly_etal:2002}.
In fact, the situation we have described here in which
we obtain a random sample
of use locations and a complete census of available locations is
referred to as ``Design 2'' by \citep{manly_etal:2002}.

\begin{comment}
Sometimes in RSF modeling activities we might have
continuous covariates and so the denominator in Eq. \ref{rsf.eq.rsf}
involves an integration over a distribution for the covariate which is
the conditional intensity of observed point locations in a point
process model. However, in a discrete landscape, entertaining pdfs for
the covariates isn't necessary \citep{royle_etal:2012mee} when we
recognize that the denominator should be the expectation over {\it
  space} and not the pdf of some covariate.
\end{comment}

\subsection{Resource Selection in SCR Models}

The key to combing RSF data with SCR data is to work with this
underlying resource utilization process and formulate SCR models in
terms of that process.  Imagine that we have a sampling device, such
as a camera trap, in {\it every} pixel. If the device operates
continually then it is no different from a telemetry instrument.  If
it operates intermittently {\it or} does not expose the entire area of
each pixel then a reasonable model for this imperfect observation is
the ``thinned'' binomial model given above, where $\phi_{0}$
represents the sampling effectiveness of the device.  For data that
arise from SCR studies, the frequency of use for each pixel {\it where}
  a trap is located serves as an intermediate latent variable that we
don't observe. From a design standpoint, the main difference between
SCR studies and telemetry is that, for SCR data, we do not have
sampling devices in all locations (pixels) in the landscape. Rather,
the data are only recorded at a subsample of them, the trap locations,
which we identify by the specific coordinates ${\bf x}_{1},\ldots,{\bf
  x}_{J}$.

So we imagine
that the hypothetical perfect data from a camera trapping study are
the counts $m({\bf x})$ only at the specific trap locations
${\bf x}_{j}$, and for all individuals in the population
$i=1,2,\ldots,N$
where $N> N_{tel}$. We denote the individual- and
trap-specific counts by $m_{ij}$ for individual with activity center
${\bf s}_{i}$ and trap location ${\bf x}_{j}$.
 In practice, many (perhaps most)
of the $m_{ij} \equiv m({\bf x}_{j},{\bf s}_{i})$ frequencies will be 0,
corresponding to individuals not captured in certain traps.
We then construct our SCR encounter
probability model based on the
view that these frequencies $m_{ij}$ are {\it latent} variables. In particular,
under the SCR model with binary observations,
 we observe a random variable
$y_{ij} = 1$  if the individual $i$ visited the pixel
containing trap $j$ and was detected.
We imagine that $y_{ij}$ is related to the latent variable $m_{ij}$ being the
event $m_{ij}>0$, as follows:
\[
 y_{ij} \sim \mbox{Bern}(p_{ij})
\]
where
\[
 p_{ij} = \Pr(m_{ij}>0) =  1-\exp(- \lambda({\bf x}_{j}|{\bf s}_{i}))
\]
This is the complementary log-log link relating
$p_{ij}$ to $\log(\lambda_{ij})$, setting
$\lambda_{ij} \equiv \lambda({\bf x}_{j}|{\bf s}_{i})$:
\[
 cloglog(p_{ij}) = log(\lambda_{ij})
\]
where
\[
 \log(\lambda_{ij} ) = \log(\phi_{0}) + a_{0} -\alpha_{1} d({\bf
   x}_{j},{\bf s}_{i})^{2} +  \alpha_{2} z({\bf x}_{j}).
\]
and we collect the constants so that $\alpha_{0} = log(\phi_{0}) +
a_{0}$ is the
 baseline encounter rate which includes
the constant intensity of use by the individual and also the baseline
rate of detection, conditional on use.

\section{The Joint RSF/SCR Likelihood}

To construct the likelihood for SCR data when we have auxiliary
covariates on space usage {\it or} direct information on space usage
from telemetry data, we regard the two samples (SCR and RSF) as
independent of one another. In practice, this might not always be the
case but (1) the telemetry data often come from a previous study;
(2) Or, the individuals are not the same, or cannot be reconciled, even if telemetry study occurs simultaneously;
(3) In cases where we {\it
  can} match some individuals between the two samples, regarding them as
independent should only entail a minor
loss of efficiency
because we are disregarding more precise information on a small number
of activity centers. Moreover, we believe, it is unlikely in practice
to expect the two samples to be completely reconcilable and that the
independence formulation is the most generally realistic.

Regarding the two data sets as being independent, our approach here
is to form the likelihood for each set of observations as a function
of the same underlying parameters and then combine them. In
particular, let ${\cal L}_{scr}(\alpha_{0}, \alpha_{1}, \alpha_{2}, N;
{\bf y}_{scr})$
be the likelihood for the SCR data in terms of the basic encounter
probability parameters and the total (unknown) population size $N$,
and let ${\cal L}_{rsf}(\alpha_{1},\alpha_{2}; {\bf m}_{rsf})$ be the
likelihood for the RSF data based on telemetry which, because the
sample size of such individuals is fixed, does not depend on $N$.
Assuming independence of the two datasets, the
joint likelihood is the product of these two pieces:
\[
{\cal L}_{rsf+scr}(\alpha_{0},\alpha_{1},\alpha_{2},N; {\bf y}_{scr},{\bf
  m}_{rsf})  = {\cal L}_{scr} \times {\cal L}_{rsf}
\]
In what follows, we provide a formulation of each likelihood
component. An  {\bf R} function for obtaining the MLEs of
model parameters is given in Appendix 1.

We adopt the notation $f(\cdot)$ to indicate the probability
distribution of whatever observable quantity is in question. e.g., $f(u)$ is
the marginal distribution of $u$ and $f(u|v)$ is the conditional
distribution of $u$ given $v$, etc. We use $g(\cdot)$ to represent the
probability distribution of latent variables.
The observation model for the SCR data for individual $i$ and trap $j$,
from sampling over $K$ encounter periods, is:
\begin{equation}
f(y_{ij}| {\bf s}_{i}) = \mbox{Bin}(K, p_{ij}({\bm \alpha}))
\label{rsf.mle.eq.cond-on-s}
\end{equation}
where
\[
p_{ij} \equiv  p(d({\bf x}_{j},{\bf s}_{i}), z({\bf x}_{j}); {\bm \alpha})
 = 1-\exp(- \lambda_{ij} )
\]
and
\[
 \lambda_{ij} = \lambda_{0} \exp(- \alpha_{1} d_{ij}^{2} + \alpha_{2}
 z({\bf x}_{j}) )
\]
We emphasize that this is conditional on the latent variables ${\bf
  s}_{i}$ (which appear in the distances $d_{ij}$). For these latent variables we
adopt the standard assumption of uniformity, ${\bf s}_{i} \sim
\mbox{Unif}({\cal S})$ for each individual $i=1,2,\ldots,N$
\citep{royle_young:2008} where ${\cal S}$ is the state-space of the
random variable ${\bf s}$.

The joint distribution of the data for
individual $i$, conditional on ${\bf s}_{i}$, is the product of $J$
binomial terms (i.e., the contributions from each of $J$ traps):
\[
  f({\bf y}_{i} | {\bf s}_{i} , {\bm \alpha}) =
  \prod_{j=1}^{J} \mbox{Bin}(K, p_{ij} ( {\bm \alpha})).
\]
The marginal likelihood \citep{borchers_efford:2008} is
computed by removing ${\bf s}_{i}$, by integration, from the
conditional-on-${\bf s}$ likelihood and regarding the {\it marginal}
distribution of the data as the likelihood. That is, we compute:
\[
  f({\bf y}_{i}|{\bm \alpha}) =
\int_{{\cal S}}  f( {\bf y}_{i} |{\bf s}_{i},{\bm \alpha}) g({\bf s}_{i}) d{\bf s}_{i}
\]
 where, under the uniformity assumption, we have
$g({\bf s}) = 1/||{\cal S}||$.
The joint likelihood for all $N$ individuals,
is the product of $N$ such terms:
\[
{\cal L}_{scr}({\bm \alpha} | {\bf y}_{1},{\bf y}_{2},\ldots, {\bf y}_{N}) = \prod_{i=1}^{N}
f({\bf y}_{i}|{\bm \alpha})
\]
In practice, we don't know $N$ and so we can't just compute the SCR
likelihood in this manner. Instead, we compute the contributions of
the $n$ observed individuals directly as given above, but then we have
to compute the likelihood contribution for the ``all 0'' encounter
history, i.e., that corresponding to unobserved individuals.  The
mechanics of computing that are the same as for an ordinary observed
encounter history, requiring that we integrate a binomial probability
of ${\bf 0}$ over the state-space ${\cal S}$:
\[
\pi_{0} = \Pr({\bf y} = {\bf 0}) = \int_{{\cal S}}
  f({\bf 0} | {\bf s} , {\bm \alpha})  d{\bf s}.
\]
We
then have to deal with the issue that $n$ itself is a random variable,
and that leads to the combinatorial term in front of the likelihood
which involves the total population size $N$. This produces the
conditional-on-$N$ or ``binomial form'' of the likelihood
\citep{borchers_efford:2008,royle:2009}:
\[
\frac{N!}{n! (N-n)!}
\left\{ \prod_{i=1}^{n} f({\bf y}_{i}|{\bm \alpha}) \right\}
\pi_{0}^{N-n}
\]

For the RSF data from the sample of individuals with telemetry devices
we adopt the same basic strategy of describing the
conditional-on-${\bf s}$ likelihood and then computing the marginal
likelihood by averaging over possible values of ${\bf s}$.
We have ${\bf m}_{i}$, the $nG \times 1$ vector of pixel counts for individual $i$,
where these counts are derived from a telemetry study or similar. We
index these elements as $m_{ig}$ for individual $i$ and grid cell $g$,
noting that our index $j$ is reserved only for trap locations, which
are a subset of the $nG$ coordinates ${\bf x}_{1},\ldots, {\bf x}_{nG}$.
The conditional-on-${\bf s}_{i}$ distribution of the telemetry data
from individual $i$ is, omitting the multinomial combinatorial term which does not
depend on parameters,
\[
 f({\bf m}_{i}  | {\bf s}_{i}, {\bm \alpha} ) \propto
\prod_{g=1}^{nG}  \pi({\bf x}_{g}| {\bf s}_{i})^{m_{ig}}
\]
where
\[
 \pi({\bf x}_{g}|{\bf s}_{i})  = \frac{ \exp( -\alpha_{1} d_{ig}^{2} +\alpha_{2} z({\bf x}_{g}) ) }
{ \sum_{g} \exp(-\alpha_{1} d_{ig}^{2} +\alpha_{2} z({\bf x}_{g}))}
\]
The marginal distribution is
\[
f({\bf m}_{i}|{\bm \alpha}) =    \int_{{\cal S}}  f({\bf m}_{i} |{\bf s}_{i},{\bm \alpha}) g({\bf s}_{i}) d{\bf s}_{i}
\]
and therefore the likelihood for the RSF data is
\[
{\cal L}_{rsf}({\bm \alpha} | {\bf m}_{1},{\bf m}_{2},\ldots, {\bf m}_{Ntel}) = \prod_{i=1}^{Ntel}
f({\bf m}_{i}|{\bm \alpha}).
\]

A key technical aspect of computing these likelihoods is the
evaluation of the 2-dimensional integral over the state-space ${\cal
  S}$, which we approximate (Appendix 1) by a
summation over a fine mesh of points.
We note also that the binomial form of the likelihood here is
expressed in terms of the parameter
$N$, the population size for the landscape defined by ${\cal
  S}$. Given ${\cal S}$, density is computed as $D({\cal S}) =
N/\mbox{area}({\cal S})$. In our simulation study below we report $N$ as the
two are equivalent summaries of the data once ${\cal S}$ is
defined. \citet{borchers_efford:2008} develop a likelihood based on a
further level of marginalization, in which $N$ is removed from the
likelihood by averaging over a  Poisson prior for $N$.

\section{Simulation Analysis}

We carried-out a simulation study using the landscape shown in
Fig. \ref{rsf.fig.habitat}, and based on populations of size $N=100$ and $N=200$
individuals with activity centers distributed uniformly over the
landscape.  This covariate was simulated by generating a field
of spatially correlated noise to emulate a typical patchy habitat
covariate relevant to habitat quality for a species.  We subjected individuals
to sampling over $K=10$ sampling periods, using a $7 \times 7$ array
of trapping devices located on the the integer coordinates $(u*5,v*5)$
for $u,v = 1,2,3,4,5,6,7$. The SCR encounter model was of the form
\[
\mbox{ cloglog}(p_{ij}) = \alpha_{0}  -\frac{1}{2\sigma^{2}} d_{ij}^{2} + \alpha_{2}  z({\bf x}_{j})
\]
with $\alpha_{0} = -2$,  $\sigma =2$ and $\alpha_{2} = 1$.
In the absence of the covariate $z$, this corresponds
to a RSF that is bivariate normal with standard
deviation 2.
These settings yielded an average of about $n=61$ individuals captured for
the $N=100$ case and about $n=123$ for the $N=200$ case. The latter case
represents what we believe is an extremely large sample size based on
our own experience and thus it should serve to gauge the large sample
bias of the likelihood estimator.

In addition to simulating data from this capture-recapture study, we
simulated 2, 4, 8, 12, 16 telemetered individuals
 to assess the
improvement in precision as sample size increases.  For all cases we
observed 20 telemetry fixes {\it per} individual, assuming individuals
were using space according to a RSF model
with the same parameters as those generating the
SCR data. We simulated 500 data sets for each scenario and, for each
data set, we fit 3 models: (i) the SCR only model, in which the
telemetry data were not used; (ii) the integrated SCR/RSF model which
combined all of the data  for jointly estimating model parameters; and
(iii) the RSF only model which just used the telemetry data alone (and
therefore $\alpha_{0}$ and $N$ are not  estimable parameters).
The focus of the simulations was to address the following basic questions:
 (1) how much does
the root mean-squared error (RMSE) of $\hat{N}$ improve as we add or increase the number of
telemetered individuals?  (2) How well does the SCR model do at
estimating the parameter of the RSF with {\it no} telemetry data?  (3)
How much does the precision of the RSF parameter improve if we add SCR
data to the telemetry data?

Results for $N=100$, $N=200$ and $N_{tel}=(2,4,8,12,16)$ are presented
in Table \ref{tab.results1}. We note that the first row of each batch
(labeled ``\mbox{\tt SCR only}'') represent the same estimator and data
configuration. These replicate runs of the SCR-only situation give us
an idea of the inherent MC error in these simulations, which is
roughly about 0.25 and 0.89 on the $N$ scale for the $N=100$ and
$N=200$ cases, respectively.  The mean $N$ for the SCR-only
estimator across all 5 simulations for $N=100$ was
$\mbox{mean}(\hat{N}) = 99.418$, an empirical bias of $0.6\%$. For
$N=200$, the estimated $N$ across all 5 simulations (5 levels of
$N_{tel}$) was $\mbox{mean}(\hat{N}) = 199.712$, an empirical bias of
about $0.15\%$, within the MC error of the true value of $N=200$.  The
results suggests a very small bias of $< 1\%$ in the MLE of $N$ for
both the  SCR-only and combined SCR/RSF estimators.  In practice, we
expect a small amount of bias in MLEs as likelihood theory only
guarantees asymptotic unbiasedness.

In terms of RMSE for estimating $N$, we see that (Table
\ref{tab.results1}), generally, there is about a 5\% reduction in RMSE
when we have at least 2 telemetered individuals. And, although there
is a lot of MC error in the RMSE quantities, it might be as much as a
10\% reduction as the sample size of captured individuals increases
under the higher $N=200$ setting. This incremental improvement in RMSE
of $\hat{N}$
makes sense because, while the
telemetry provides considerable information about the structural
parameters of the model, it provides no information about mean $p$,
i.e. $\alpha_{0}$, which comes only from the SCR data. Thus estimating
$N$ benefits only slightly from the addition of telemetry data.

The MLE of the RSF parameter $\alpha_{2}$ exhibits negligible or no
bias under {\it both} the SCR only and SCR/RSF estimators. It is
well-estimated from SCR data alone and even better than RSF data alone
(in terms of RMSE) until we have more than 200 or so telemetry
observations.  The biggest improvement from the use of telemetry data
comes in estimating the parameter $\sigma$. We see that $\hat{\sigma}$
is effectively unbiased, and there is a very large improvement in RMSE
of $\hat{\sigma}$, perhaps as much as 50-60\% in some cases, when the
telemetry data are used in the combined estimator (that really doesn't
translate much into improvements in estimating $N$ as we saw
previously).  Improvement due to adding telemetry data diminishes as
the expected sample sizes increases, and so telemetry data does less
to improve the precision of $\hat{\sigma}$ and $\hat{\alpha}_{2}$ for
$N=200$ than for $N=100$. This is because the SCR data along are
informative about both of those parameters.

The results as they concern likelihood estimation of $N$ suggest that
there is not a substantial benefit to having telemetry
data. Estimators ``SCR only'' and ``SCR/RSF'' both appear
approximately unbiased for $N=100$ and $N=200$, and for any sample
size of telemetered individuals. The RMSE is only 5-10\% improved with
the addition of telemetry information.  However, we find that there is
substantial bias in $\hat{N}$ if we use the {\it misspecified} model
that contains no resource selection component. That is if we leave the
covariate $z({\bf x})$ out of the model and incorrectly fit a model
with symmetric and spatially constant encounter model, we see about
20\% bias in the estimates of $N$ in a limited simulation study that
we carried-out (Tab. \ref{tab.bias}). As such, accounting for resource
selection is important, even though, when accounted for, telemetry
data only improves the estimator incrementally.
In addition, we find that the importance of telemetry data is
relatively more important for smaller sample sizes. We carried-out one
simulation study for the $N=100$ case but with lower average encounter
probabilities, setting $\alpha_{0}=-3$. This produces
relatively smaller data sets with  $E[n] = 37$. The results are shown
in Tab. \ref{tab.lowp}. There are some important features evident from
this table. First, as a result of the small samples, the MLE of $N$ is
biased for both SCR only and SCR/RSF estimators although less biased
for the SCR/RSF estimator than for SCR only. The persistent bias in
$\hat{N}$ for both models results from the information about
$\alpha_{0}$ coming only from SCR data, and that estimator itself is
intrinsically biased in small samples.  Conversely, the estimator of
$\alpha_{2}$, the RSF parameter, appears unbiased for all 3 estimators
(SCR only, SCR/RSF and RSF only), as does the estimator of $\sigma$.
We see relatively larger improvements in RMSE (compared with Tab. \ref{tab.results1}) of $\hat{N}$, and those improvements
increase substantially as $N_{tel}$ increases.

\section{Discussion}

How animals use space is a fundamental interest to ecologists, and
important in the conservation and management of many species.
Normally this is done
by telemetry and models referred to as resource selection functions
\citep{manly_etal:2002}.  Conversely, spatial capture-recapture models
have grown in popularity over the last several years
\citep{efford:2004,borchers_efford:2008, royle:2008,
  efford_etal:2009ecol,royle_etal:2009ecol, gardner_etal:2010ecol,
  gardner_etal:2010jwm, kery_etal:2010,
  sollmann_etal:2011,mollet_etal:2012,gopalaswamy_etal:2012}. These,
and indeed, most,
 development and applications of SCR models have focused on density
estimation, not understanding space usage.  However, it is intuitive that space
usage should affect encounter probability and thus it should be highly relevant
to density estimation in SCR applications. Despite this, a description
of the
relationship between encounter
 probability and space usage has not
been developed in the literature on spatial
capture-recapture models.  Essentially all
published applications of SCR models to date have been based on
simplistic encounter
probability models that are symmetric
and do not vary across space. One exception is
\citet{royle_etal:2012ecol} who developed SCR models that use
ecological distance metrics (``least-cost path'')
instead of normal Euclidean distance. Here
we developed an SCR model in terms of a basic underlying model of
space or resource use, that is consistent with existing views of
resource selection functions (RSFs) \citep{manly_etal:2002}.

In developing the SCR model in terms of an underlying model of space
usage, we achieve a number of enormously useful extensions of existing
SCR and RSF methods:
(1) We have shown how to integrate classical RSF data from
telemetry with spatial capture-recapture data based on individual
encounter histories obtained by classical arrays of encounter devices
or traps. This leads to an improvement in our ability to estimate
density, and also an improvement in our ability to estimate parameters
of the RSF function.  Thus, the combined model is both an extension of
standard SCR models and also and extension of standard RSF models. As
many animal population studies have auxiliary telemetry information,
the ability to incorporate such information into SCR studies has
enormous applicability and immediate benefits in many studies.
While adding RSF data to SCR data may increase precision of the MLE of $N$
only incrementally, the effect can be more substantial in sparse data sets
and, generally, RSF produces
relatively huge gains in precision in the MLE of $\sigma$.
(2) We have shown that one can estimate RSF model parameters
directly from SCR data {\it alone}.
While further exploration of this point
is necessary,
it does establish clearly that SCR
models {\it are} explicit models of space usage. Because
capture-recapture studies are, arguably, more widespread than
telemetry studies alone, this greatly broadens the utility and
importance of data from those studies.
(3) It is also now clear
that one of the important parameters of SCR models, that controlling
``home range radius'', can be directly estimated from telemetry data
alone.
The combined RSF+SCR model does yield large improvements in estimation
of $\sigma$. As a practical matter, this suggests we could estimate
$\sigma$ entirely from data extrinsic to the SCR study which might
provide great freedom in the design of SCR studies. For example, traps
could be spaced far enough apart to generate relatively few (even no) spatial recaptures,
but dramatically increase the coverage of the population, i.e., the
observed sample size of captured individuals relative to $N$.
(4) Finally, we found that an
ordinary SCR model with symmetric encounter probability model produces
extremely biased estimates of $N$ when the population of individuals
does exhibit resource selection.  As such, it is important to account
for space usage when important covariates are known to influence
space usage patterns.

Use of telemetry data in capture-recapture studies has been suggested
previously. For example, \citet{white_shenk:2001} and
\citet{ivan:2012} suggested using telemetry data to estimate the
quantity ``probability that an individual is exposed to sampling'' but
their estimator requires that individuals are sampled in proportion to
this unknown quantity, which seems impossible to achieve in many
studies. In addition, they do not directly integrate the telemetry
data with the capture-recapture model so that common parameters are
jointly estimated. In fact, they don't acknowledge shared
parameters of the two models.  \citet{sollmann_etal:inprep} did
recognize this, and used some telemetry data to estimate directly the
parameter $\sigma$ from the bivariate normal SCR model in order to
improve estimates of density. This was an important conceptual
development in the sense that it recognized the relationship between
SCR models and models of space usage, but their model did not include
an explicit resource selection component, and they did not implement a
joint estimation framework.

We developed a formal analysis framework here based on marginal
likelihood \citep{borchers_efford:2008}.
In principle, Bayesian analysis does not pose any unique challenges
for this new class of models although we expect some loss of
computational efficiency due to the increased number of times the
components of the likelihood would need to be evaluated.  We imagine
that some problems would benefit from a Bayesian formulation,
however. For example, using an open population model that allows for
recruitment and survival over time \citep{gardner_etal:2010ecol} is
convenient to develop in the {\bf BUGS} language and incorporating
information on unmarked individuals has been done using Bayesian
formulations of SCR models \citep{chandler_royle:2012,
  sollmann_etal:inprep} but, so far, not likelihood methods.

In our formulation of the joint likelihood for RSF and SCR data, we
assumed the data from  capture-recapture and telemetry studies were
independent of one another. This implies that whether or not an
individual enters into one of the data sets has no effect on whether
it enters into the other data set. We cannot foresee situations in
which violation of this assumption should be problematic or invalidate
the estimator under the independence assumption.  In some cases it
might so happen that some individuals appear in {\it both} the RSF and
SCR data sets. In this case, ignoring that information should entail
only an incremental decrease in precision because a slight bit of
information about an individuals activity center is
disregarded. Heuristically, an SCR observation (encounter in a trap)
is like one additional telemetry observation, and so the
misspecification (independence)
regards the
two pieces of information as having separate activity centers.
 Our model pretends that we don't know anything
about the telemetered individuals in terms of their encounter history
in traps.  In principle it shouldn't be difficult to admit a formal
reconciliation of individuals between the two lists. In that case, we
just combine the two conditional likelihoods before we integrate ${\bf
  s}$ from the conditional likelihood. This would be almost trivial to
do if {\it all} individuals were reconcilable (or none as in the case
we have covered here) but, in general , we think you will always have
an intermediate case -- i.e., either none will be or at most a subset
of telemetered individuals will be known. More likely you have variations of ``well, that
guy looks telemetered but we don't know which guy it is....hmmm'' and
that case, basically a type of marking uncertainty or
misclassification, is clearly more difficult to deal with.

We conclude that the key benefit of our combined SCR/RSF model is its
ability integrate realistic patterns of space usage directly into SCR
models and avoid extreme bias in estimating $N$ and, secondarily, we
are able to obtain RSF information from SCR alone.  Therefore, our new
class of integrated SCR/RSF models allows investigators to model how
the landscape and habitat influence movement and space usage of
individuals around their home range, using non-invasively collected
capture-recapture data or capture-recapture data augmented with
telemetry data.  This should improve our ability to understand, and
study, aspects of space usage and it might, ultimately, aid in
addressing conservation-related problems such as reserve or corridor
design. And, it should greatly expand the relevance and utility of
spatial capture-recapture beyond simply its use for density
estimation.

\section*{Acknowledgments}