\section{Introduction}
Understanding the vertical disk structure of spiral galaxies is crucial for
learning about how galaxies form and evolve over time
\citep{Samland03,Moster10}, as well as producing an accurate estimate for the
baryon content of disks \citep{Bershady10a,Bershady10b}. In the Milky Way,
with survey measurements of many thousands of individual stars, disk structure
can be examined in great detail \citep[e.g.][]{Bovy12}. When viewed as if it
were an external galaxy, the Milky Way's light distribution can be fit with just three main disk components:
the thick disk, containing old stars (scaleheight h$_{z}\sim1$ kpc); the
thin disk, the most domininant light-weighted component (h$_{z}\sim300$ pc);
and a super-thin disk, consisting mostly of young stars \citep[h$_z\sim100$ pc,][]{vanDokkum94}. The latter component goes by many names in the
literature; it is frequently referred to as the `young' or `star-forming'
disk. We prefer the term `super-thin' as it carries no assumption about the age of the disk.

The existence of both thin and thick disks in other spiral galaxies, while an
uncertainty for many years \citep{vanDokkum94}, has now been confirmed
\citep{Dalcanton02,Comeron11}. However, the ubiquitous presence of dust near
the midplane of these galaxies \citep[especially for massive ones most similar
to the Milky Way,][]{Dalcanton04} prevents easy study of any super-thin
components.

Here we present initial results from our ongoing work to map edge-on disk
structure, using a combination of
high-resolution near-infrared (NIR) imaging and advanced Monte Carlo radiation
transfer models to remove the effects of dust and probe vertical surface
brightness profiles all the way down to the midplane. Our method is
illustrated by its application to NGC 891, but is broadly applicable to all
edge-on galaxies.

\articlefigure[scale=0.45]{schechtman-rooka_fig1_color.eps}{fig:truncation}{Vertical
  profiles of NGC 891 at different radii. Data is shown as black dots, while
  dashed dark gray lines indicate the best fitting three-component disk model with
  exponential radial and vertical profiles. Inner regions deviating from this
  idealized model are marked by arrows indicating the additional structure
  required to fit the observed surface-brightness profiles. The light gray solid
  lines show our best fitting model including the inner disk truncation, bar,
  and nuclear disk. }
\section{Data and Model}
The following is a brief overview of our data and our radition transfer
modeling. Full details of our methodology, both for data reduction and producing our
attenuation correction, can be found in \citet{Schechtman-Rook13}.
\subsection{NIR Data}
Data were taken on the WIYN High-resolution InfraRed Camera
\citep[WHIRC;][]{Meixner10} on the WIYN 3.5 m telescope in 2011 October. NGC
891 is significantly larger in radius than the WHIRC field of view, so three
pointings were required to image out to just over $\pm$10 kpc in radius. Sky
images were interspersed with data exposures in order to compute an accurate
background estimate for such an extended object. The
images were processed using in-house reduction software, which performs all
aspects of data reduction from initial trimming of the images to creating the
final mosaic.
\subsection{Attenuation Correction}
To compute the attenuation correction we first model NGC 891 using the
radiation transfer software HYPERION \citep{Robitaille11}; a model was created
which matched the integrated spectral energy distribution of NGC 891, and
images of this model with and without dust were compared to produce a
pixel-by-pixel attenuation map. The attenuation in a given pixel was
well-predicted by the model's K$_{\mathrm{s}}$-4.5$\mu$m color.
\articlefigure[scale=0.53]{schechtman-rooka_fig2_color_compressed.eps}{fig:fastrotators}{JHK$_{\mathrm{s}}$
  false-color images of the isolated fast
  rotators (other than NGC 891) in our sample. The dashed ellipse
  indicates the inner truncation in visible in our images of NGC 4565.}

In order to use this color to correct our WHIRC data, we obtained an archival
{\it Spitzer} IRAC image of NGC 891. While this image is at significantly
lower resolution compared to our NIR data, the intrinsic smoothness of the luminosity
profile at 4.5$\mu$m allows us to preserve the resolution of our WHIRC data.

\section{Results}
\subsection{NGC 891}
Observed and attenuation-corrected JHK$_{\mathrm{s}}$ false-color images are
given in\\ \citet{Schechtman-Rook13}. We fit models with up to three discrete disk components, producing
best-fitting models using a Levenberg-Marquardt non-linear least-squares
fitter. For these initial models we ignored the central 3 kpc of NGC 891 to
avoid bulge/bar contamination. We find that three disks are necessary in order
to reproduce NGC 891's intrinsic K$_{\mathrm{s}}$-band surface-brightness at all
heights. These three components have very scale-heights to the Milky Way, with
the super-thin disk scale-height $h_{z}\approx$60 pc.

We then investigated the ability of our models to fit the data within 3 kpc,
and found that models with just disks and R$^{1/4}$-law bulges were
  insufficient. In all, we had to add three new features to our best-fitting 3
  disk model: an inner disk truncation for all three disk components at R$\approx$3.1 kpc; a bar, which was
  well-approximated by an exponential disk with $h_{z}\approx$400 pc; and
  a nuclear extension of the super-thin disk, which has a scale-length of only
  250 pc. An R$^{1/4}$-law bulge is not necessary to reproduce the light
  profile of NGC 891, and actually results in significantly worse fits if used
  in place of the nuclear super-thin disk. The light from the bar is roughly
  equivalent to the light lost due to the inner disk truncation, implying that
  the bar is composed of stars swept up from NGC 891's disks.

\articlefigure[scale=0.53]{schechtman-rooka_fig3_color_compressed.eps}{fig:slowrotators}{JHK$_{\mathrm{s}}$
  false-color images of the two slow-rotators (V$<$90 km sec$^{-1}$) in our
  sample. These galaxies show little evidence for dust-lane attenuation but
  have extensive star-clusters throughout the vertical range of the disk.}
\subsection{Other Galaxies}
In addition to NGC 891, we have a small sample of nearby edge-on galaxies
suitable for this analysis. We generally selected systems with properties
similar to the Milky Way and NGC 891 (Figure \ref{fig:fastrotators}); however
we have also included some low-mass spirals (Figure \ref{fig:slowrotators}) as well as a galaxy undergoing a
minor interaction (bottom panel of Figure \ref{fig:fastrotators}). In both
figures the images are not attenuation corrected.

Of special note is NGC~4565 (second from bottom panel of Figure \ref{fig:fastrotators}),
which has a clearly visible inner disk truncation. Interestingly this also
appears to be the situation in the Milky Way \citep[at least in dust
emission, see][]{Robitaille12}. While three galaxies is too small a sample to draw
definitive conclusions, it appears inner disk truncations may be a common
feature of massive spiral galaxies.

\acknowledgements This research was supported by NSF AST-1009471.


\section{Introduction}
\label{sec:intro}

As of today, about 250 extrasolar planets have been detected through
both radial velocity and transit measurements.  The combination of
these two detection methods makes it possible to determine the angle
between the sky projection of the stellar spin axis and the orbit
normal (that we will hereafter call the {\em projected inclination
  angle}) using the Rossiter--McLaughlin effect.  Among the 70 planets
for which this angle has been obtained, about 30 show significant
misalignment, with a few even displaying retrograde orbits.  The
distance of most of these planets to their host star is below 0.1~au,
and most of these objects have a mass larger than half that of
Jupiter.  The sample is therefore highly biased towards massive and
close planets.  However, these detections indicate that misalignement
is not rare.

If planets on short orbits have formed in and migrated through a disc,
they would be in a plane misaligned with that of the star's equator
only if the disc itself were misaligned, or if interactions had taken
place that would put the planets on misaligned orbits.  Such
interactions, occuring after dissipation of the disc, include secular
perturbation by a distant stellar or planetary companion that produces
Kozai cycles with tidal friction (Fabrycky \& Tremaine 2007; Wu,
Murray \& Ramsahai 2007; Naoz et al. 2011), planet--planet scattering
(Chatterjee et al.  2008) or secular chaos (Wu \& Lithwick 2011).
Processes leading to the misalignment of the disc itself include
accretion  (at later times) of material with angular momentum
misaligned with that of the star (Bate, Lodato \& Pringle 2010), close
encounter with a star while the disc accretes from an extended
envelope (Thies et al.  2011), or precession induced by the
interaction with a stellar companion in a binary system (Batygin
2012).  Planetary orbits inclined with respect to the stellar
equatorial plane would also result if the stellar spin axis were
tilted due to interaction with the disc (Foucart \& Lai 2011; Lai,
Foucart \& Lin 2011).

The interactions studied so far, which would misalign an orbit
initially in the stellar equatorial plane, are presumed to take place
after the disc has dissipated.  Note, however, that if the orbital
plane were inclined relative to the disc while it is still present,
the orbit would not necessarily (re)align as a result of the
frictional force exerted by the disc on the planet when it crosses it.
This is because the secular gravitational interaction between an
inclined planet and the disc's outer parts leads to Kozai cycles that
maintain either the inclination or the eccentricity at large values,
thus delaying alignment with the disc and circularization of the orbit
(Terquem \& Ajmia 2010; Teyssandier, Terquem \& Papaloizou 2013).
This Kozai effect could not be seen in previous simulations of planets
on inclined orbits in which the disc was truncated at rather small
radii (e.g., Marzari \& Nelson 2009), but was observed in the SPH
simulations of Xiang-Gruess \& Papaloizou (2013) which used a more
extended disc and higher orbital inclinations.

So far, the mechanisms that have been put forward to produce a
misalignment of a planetary orbit relative to a disc rely on the
interaction of the planet with at least one stellar or planetary
companion.  Here, we investigate the interaction of the planet with
the disc itself, and show that a warping of the disc's outer parts may
lead to a significant inclination of the planet's orbit relative to
the disc's inner parts.

\subsection{Warping of protoplanetary discs}

The outer parts of protostellar discs may become warped as a result of
tidal interaction with a binary star (Terquem \& Bertout 1993, 1996;
Papaloizou \& Terquem 1995).  Most stars are in multiple systems
(e.g., K\"ohler 2008) and the distribution of the separation of
pre--main sequence binaries has a peak around 30 au (Leinert et al.
1993).  It is therefore expected that tidal interactions between a
circumstellar disc and the companion are common.  If the disc and the
orbit are not coplanar, such interactions lead to a warping of the
disc's outer parts.  Evidence for warped discs have been seen, e.g.,
in the young triple system WL~20 (Barsony 2002) and in the Bok globule
CB~26 (Launhardt \& Sargent 2001).  Also, in the young binary system
HK~Tau, the disc around one of the stars is observed not to be
coplanar with the orbital plane (Stapelfeldt et al. 1998), which
implies that it may be warped.

Warping may also occur if the disc accrete at later times material
with angular momentum misaligned with that of the star (Bate, Lodato
\& Pringle 2010).  This is expected as stars form within clusters in
which gravitational interactions are important.  It was suggested by
Tremaine (1991) that the obliquities of the outer planets in our solar
system could have been produced by a change in the direction of the
angular momentum of the planetary system after it had formed.  The
proposed causes for this change were a time-dependant angular momentum
of the material falling into the disc, possibly due to inhomogeneities
in the collapsing core, or gravitational torques due to mass
concentrations within the core.  Whether asymmetric infall might
actually result in a warped disc or not is still a matter of debate.
The accretion of material with misaligned angular momentum onto the
disc outer parts does produce an inclination of these outer parts
with respect to the inner regions, but such an inclination may be
smoothed out by viscous diffusion through the disc.  Whether the final
configuration is planar or not depends on whether the viscous
diffusion timescale is smaller or larger than the precession timescale
(see discussion below).  In the simulations by Bate, Lodato \& Pringle
(2010), the disc would settle into a planar configuration for the
range of parameters used, but whether it is always the case or not is
still an open question.

It has been proposed that, in a stellar binary system in which a disc
is misaligned with respect to the orbital plane, planets forming into
this disc would themselves be inclined relative to the central star's
equatorial plane (Batygin 2012).  This is because the disc is tilted
due to the precession of the disc's angular momentum vector around the
orbital angular momentum vector.  However, warping was neglected in
this study, under the argument that a self--gravitating disc similar
to the one considered would only precess as a rigid body if it
remained unwarped (see also Batygin et al. 2011).

There have been extensive studies of the dynamics of both galatic and
protostellar discs that have shown that, for a disc to precess as a
rigid body in a non spherically symmetric potential, it actually has
to bend to alter the precession frequency at each radius so that the
rate is everywhere the same (Sparke \& Casertano 1988; Hofner \&
Sparke 1994).  The disc therefore settles into a discrete bending mode
(representing a warp) that is referred to as the modified tilt mode
because, in the limit that the external potential is spherically
symmetric, it reduces to the trivial rigid--tilt mode.  Differential
precession can be smoothed out by gravitational torques in a
self--gravitating disc (Toomre 1983), and also by radial pressure
forces and viscous diffusion (Papaloizou \& Terquem 1995, Larwood et
al. 1996).  These processes are efficient if the timescale on which
the different parts of the disc communicate with each other is shorter
than the precession timescale.  In a protostellar disc in which
pressure forces are likely to dominate over self-gravitating forces,
rigid body precession is achieved if the sound crossing time through
the disc is much smaller than the precession period, so that bending
waves propagate through the disc sufficiently fast for the different
parts of the disc to adjust their precession rate to a constant value.
Communication becomes diffusive in a disc in which turbulent stresses
dominate.

In general, in a disc dominated by pressure forces, rigid body
precession of the whole disc is achieved with only a small degree of
warping (Larwood et al. 1996, Larwood 1997).  In that case, as
envisionned by Batygin (2012), the relative inclination of the disc
with the central star's equatorial plane is due mainly to the
precession, not to the warp.  However, it was found in the numerical
simulations of Larwood et al. (1996) that the disc could maintain
itself as a long--lived and coherent structure even when a strong warp
develops and the inner and outer parts of the disc precess at
different rates.  In that case, there is a coupling between the two
parts of the disc due to pressure gradients, and the rigidly
precessing outer parts drag the inner parts behind it.

Planets forming in the disc's inner parts would therefore interact
with the strongly misaligned outer parts, and this is the type of
situation we investigate in this paper.

\subsection{Setup}

We consider a planet that has formed in the inner parts of a disc
which outer parts have been warped according, possibly, to one of the
processes described in the previous section.

Note that, if the disc's warp were due to a companion on a noncoplanar
orbit in a binary system, gravitational interaction between this
companion and the planet would affect the orbital evolution of the
planet.  This is not taken into account in this paper, where we focuss
on the interaction between the warped disc and the planet
independently of the cause of the warp.  The dynamics of a planet in a
noncoplanar binary system, taking into account disc warping and the
presence of the companion, will be the subject of a further paper.

We consider the stage of evolution after an inner cavity has been
opened in the disc, and we assume that the planet's orbit lies in this
cavity.  Inner holes are commonly observed in the so--called
transitional discs (Andrews et al.  2013).  They are usually assumed
to be opened either by X--ray photoevaporation (Owen et al. 2011) or
by tidal interaction with low mass object(s).  A giant planet is
expected to clear--out a gap, so that the parts of the disc interior
to the orbit are no longer replenished efficiently from the outer disc
and accrete onto the central star on a viscous timescale if the disc
is active.  Note that, in this scenario, the planet may be expected to
migrate in as a result of tidal interaction with the outer parts of
the disc.  However, migration would be slow if the mass of gas in the
vicinity of the planet were significantly smaller than the mass of the
planet itself.  Both planet's migration and accretion of the disc
material would also be very inefficient if the disc were passive in
the region where the planet is located.  In that case, the planet
would move in a gap but there would still be material interior to its
orbit.

In the calculations presented in this paper, we assume for simplicity
that there is no mass interior to the planet's orbit, but our results
would hardly be changed if only a gap were cleared, with still
material orbiting in the region inside the orbit.  This is because the
gravitational interaction with the disc is dominated by the parts of
the disc which are outside the orbit (Terquem \& Ajmia 2010).  The
important point here is that we assume there is no friction between
the disc and the planet, so that the planetary orbit evolves only
under the gravitational potential from the disc.

\subsection{Plan of the paper}

In section~\ref{sec:interaction}, we develop a secular analysis of the
interaction between a planet and a warped disc.  We first write the
gravitational potential exterted by the disc on the planet and
describe the disc model used in this study in
section~\ref{sec:gravitation}.  In section~\ref{sec:expansion}, we
expand the potential in the case where the radius of the inner edge of
the disc is large compared to the planet's semimajor axis, and
calculate the secular terms as a function of the orbital elements to
fourth order in the eccentricity in section~\ref{sec:secular}.
Lagrange's planetary equations, valid to second order in the
eccentricity, are then written in section~\ref{sec:Lagrange}.  We
solve these equations in section~\ref{sec:solutions}, both for a flat
and a warped disc, to obtain the evolution of the orbital elements.
In the case of a flat disc, we recover the results of Terquem \& Ajmia
(2010), who showed that perturbation by an outer flat disc inclined
relative to the planetary orbit may lead to Kozai cycles.  In the case
of a warped disc, we find that, in general, the interaction bewteen
the planet and the disc leads to the precession of the planetary
orbit.  The orbit therefore becomes inclined relative to the disc's
inner parts, with no increase of its eccentricity.  We find that there
is an instability when the inclination approaches $90^{\circ}$.  In
that case, the eccentricity starts to increase and both the
inclination and the eccentricity have very large amplitude variations,
with the orbit becoming retrograde.
In section~\ref{sec:simulations}, we perform simulations that enable
us to follow the evolution of the orbit when the eccentricity gets
large.  Finally, in section~\ref{sec:discussion}, we summarize and
discuss our results.

\section{Secular analysis}
\label{sec:interaction}

\subsection{Gravitational potential and disc model}
\label{sec:gravitation}

We consider a planet of mass $M_p$ orbiting around a star of mass
$M_{\star}$ which is itself surrounded by a disc of mass $M_d$.  We denote
by $(x,y,z)$ a Cartesian coordinate system centred on
the star and $(r, \theta, z)$ the associated cylindrical polar
coordinates.
The warped disc has its central parts tangent to the $(x,y)$--plane
and its elevation above this plane is taken to be:
\begin{eqnarray}
z(r, \theta) = \left( \frac{r}{R_o} \right)^p H_o \sin \theta,
\label{zwarp}
\end{eqnarray}

\noindent where $R_o$ is the disc's outer radius, $H_o$ is the
maximum elevation at $r=R_o$, and $p>1$ is an integer which controls
the warp's curvature ($p=1$ would correspond to a disc rigidly tilted
relative to the $(x,y)$--plane).  We introduce $\alpha$, the inclination of
the disc's outer edge with respect to the $(x,y)$--plane, so that: \begin{eqnarray}
H_o= R_o \tan \alpha,
\label{Ho}
\end{eqnarray}

\noindent and let $R_i$ be the disc's inner edge.  Initially, the orbit of
the planet is in the $(x,y)$--plane and has a semimajor axis $a<R_i$.
We suppose that the angular momentum of the disc is large compared to
that of the planet's orbit so that the effect of the planet on the
disc is negligible.

Note that if the disc as a whole is tilted with respect to the
stellar equatorial plane, the $(x,y)$--plane does not coincide
with the equatorial plane of the star.  This would be the case if what
caused the warp also induced a precession of the disc (e.g.,
fly--by of a star).

The gravitational potential exerted by the disc at the location
of the planet is:

\begin{eqnarray} \Phi_d = -G \int_{R_i}^{R_o} \Sigma(r) r dr \int_{0}^{2\pi}
\frac{d\theta}{\sqrt{r^2+r_p^2 - 2rr_p \cos \left( \theta - \theta_p
    \right) + \left[ z_p - \left( \frac{r}{R_o} \right)^p R_o \tan \alpha
    \sin \theta \right]^2}},
\label{potential}  \end{eqnarray}

\noindent where the subscript $p$ refers to the planet, $G$ is the
gravitational constant and $\Sigma$
is the mass density in the disc.  We assume:

\begin{eqnarray}
\Sigma(r)= \Sigma_0 \left( \frac{r}{R_o} \right)^{-q},
\label{sigma}
\end{eqnarray}

\noindent where $q$ is a positive number.  In the numerical
simulations described below, we will set $q=1/2$.  The coefficient
$\Sigma_0$ can be expressed in term of the disc mass $M_d$:

\begin{eqnarray} \Sigma_0 = \frac{(-q+2) M_d}{ 2 \pi R_o^2 \left[ 1 - \left(R_i/R_o
    \right)^{-q+2} \right]}.  \end{eqnarray}

The total potential at the location of the planet has four
contributions: one from the disc, which is given by
equation~(\ref{potential}), one from the central star, an indirect
contribution since the center of the coodinate system
is accelerated, and another indirect contribution since
that, in general, the $(x,y)$--plane rotates as a result of the disc's
precession.  However, these last three terms do not yield any secular
contribution and will be ignored in this section.

\subsection{Expansion of the gravitational potential}
\label{sec:expansion}

We define $u \equiv r_p/r$, $v \equiv z_p/r$ and $m \equiv \left(
  r/R_o \right)^{p-1} \tan \alpha$.  The potential in
equation~(\ref{potential}) can then be written in the form:

\begin{eqnarray}
\Phi_d = -G \int_{R_i}^{R_o} \Sigma(r) dr {\cal I} \left( r,u,v \right) ,
\label{potential2}
\end{eqnarray}

\noindent with:

\begin{eqnarray} {\cal I} \left(r,u,v \right) = \int_0^{2 \pi} \frac{ d
  \theta}{\sqrt{1+u^2-2u \cos \left( \theta - \theta_p \right) + \left( v -
    m \sin \theta \right)^2}} .
\label{Idefinition}
\end{eqnarray}

We assume that $R_i \gg r_p$ and $R_i \gg z_p$ so that $u$ and $v$ are
small compared to unity everywhere in the disc.  To second--order in
$u$ and $v$, ${\cal I} \left( r,u,v \right)$ can then be expanded as:

\begin{eqnarray} {\cal I} \left( r,u,v \right) = {\cal I} \left( r,0,0 \right) +
\frac{u^2}{2} \frac{\partial^2 {\cal I}}{\partial u^2} \left( r,0,0 \right)
 + \frac{v^2}{2} \frac{\partial^2 {\cal I}}{\partial v^2} \left( r,0,0
\right)  + u v \frac{\partial^2 {\cal I}}{\partial u \partial
  v} \left( r,0,0 \right) ,
\label{Iexpansion}
\end{eqnarray}

\noindent where the omitted first--order terms vanish.  It is shown in the
appendix~\ref{appendixA} that this expansion yields the following
expression for the potential~(\ref{potential2}):

\begin{eqnarray} \Phi_d = -\frac{G \Sigma_0}{R_o} \left( {\cal I}_1 R_o^2 + {\cal I}_2
    r_p^2  + {\cal I}_3 r_p^2 \cos 2 \theta_p
   + {\cal I}_4 z_p^2  + {\cal I}_5 r_p z_p \sin
    \theta_p  \right) ,
\label{potential3}
\end{eqnarray}

\noindent where ${\cal I}_1, {\cal I}_2, {\cal I}_3, {\cal I}_4$ and
${\cal I}_5$ are dimensionless integrals over $r$ which depend only on
disc parameters.  Their expression in term of elliptic
integrals is given in appendix~\ref{appendixA}.

\subsection{Secular terms}
\label{sec:secular}

We wish to calculate the long--term effect of the interaction between
the disc and the planet, obtained by averaging the
potential~(\ref{potential3}) over the mean longitude of the planet. In
order to perform this averaging, we first express the planet's polar
coordinates $r_p, \theta_p$ and $z_p$ in terms of the planet's orbital
elements.  We denote by $a$, $e$, $I$, $\omega$ and $\Omega$ the
semimajor axis, eccentricity, inclination, argument of pericentre and
longitude of ascending node, respectively, of the orbit, defined with
respect to the reference frame $(x,y,z)$.  We then have (Murray \&
Dermott 1999):

\begin{eqnarray}
r_p \cos \theta_p & = & R \left[ \cos \Omega \cos \left( \omega + f \right)
- \sin \Omega \sin \left( \omega + f \right) \cos I \right] ,
\label{rpcos} \\
r_p \sin \theta_p & = & R \left[ \sin \Omega \cos \left( \omega + f \right)
+ \cos \Omega \sin \left( \omega + f \right) \cos I \right] ,
\label{rpsin}
\\
z_p & = & R \sin \left( \omega + f \right) \sin I,
\label{zp}
\end{eqnarray}

\noindent where $R$ is the distance between the planet and the star
center and $f$ is the true anomaly of the planet (polar angle in the
orbital plane relative to the pericentre).  To perform a time
averaging of the potential~(\ref{potential3}), we further express $R$,
$\sin f$ and $\cos f$ in terms of the mean anomaly $M$, which is $2
\pi$--periodic and a linear function of time.  Expansions of $R$, $\sin
f$ and $\cos f$ to fourth order in $e$ are given in
appendix~\ref{appendixB}.  We insert these series expansion in
equations~(\ref{rpcos}), (\ref{rpsin}) and~(\ref{zp}) above, and
expand $r_p^2$, $r_p^2 \cos 2 \theta_p $, $z_p^2$, $r_p
z_p \sin \theta_p$, and therefore the potential~(\ref{potential3}), in
terms of $M$ to fourth order in $e$.  Finally we perform a time
average to obtain the secular potential $\left< \Phi_d \right>$.
Time--averaging eliminates from the expansion the terms which are odd
in $e$.  As it happens, the coefficients of $e^4$ also vanish.  The
series expansions we obtain are therefore valid up to fifth order in $e$.
The time--averages of $r_p^2$, $r_p^2 \cos 2 \theta_p $,
$z_p^2$, $r_p z_p \sin \theta_p$, which are needed to compute $\left<
  \Phi_d \right>$ from equation~(\ref{potential3}), are given below:

\begin{eqnarray}
  \frac{ \left< r_p^2 \right>}{a^2} & = & \frac{1}{4} \left( 3 + \cos 2I
  \right) + \frac{e^2}{8} \left[ 3 \left( 3 + \cos 2I \right) + 10
    \cos 2 \omega \sin^2 I \right] + {\cal O}(e^6),
\label{rp2}
\end{eqnarray}

\begin{eqnarray}
  \frac{ \left< r_p^2 \cos 2 \theta_p \right>}{a^2} = & &
  \frac{1}{2} \cos 2 \Omega \sin^2 I \nonumber \\ & + &
  \frac{e^2}{8} \left[ 5 \left( 3 + \cos 2I \right)
    \cos 2 \omega \cos 2 \Omega + 6
    \cos 2 \Omega \sin^2 I - 20 \cos I \sin 2 \omega \sin 2 \Omega \right]
\nonumber \\
  & + & {\cal O}(e^6),
\label{rp2cos}
\end{eqnarray}

\begin{eqnarray}
\frac{ \left< z_p^2 \right>}{a^2} & = & \frac{1}{2} \sin^2 I +
\frac{e^2}{4} \left( 3 - 5 \cos 2 \omega \right) \sin^2 I +  {\cal O}(e^6),
\label{zp2}
\end{eqnarray}

\begin{eqnarray}
  \frac{ \left< r_p z_p \sin \theta_p \right>}{a^2} = & & \frac{1}{4}
\sin 2I \cos \Omega  \nonumber \\
& + & \frac{e^2}{4} \sin I \left( 3 \cos I \cos \Omega -
5 \cos I \cos 2 \omega \cos \Omega + 5 \sin 2 \omega \sin \Omega \right)
\nonumber \\
& + & {\cal O}(e^6).
\label{rpzpsin}
\end{eqnarray}

\subsection{Lagrange's planetary equations}
\label{sec:Lagrange}

In secular theory, the semimajor axis $a$ is constant, and the
time variation of $e$, $\Omega$, $\omega$ and $I$ are given by the
Lagrange's planetary equations (e.g., Roy 1978):

\begin{eqnarray}
\frac{de}{dt} & = & - \frac{ \sqrt{1-e^2}}{n a^2 e}
\frac{\partial {\cal R}}{ \partial \omega} ,
\label{dedt}
\\
\frac{d \Omega}{dt} & = & \frac{1}{n a^2 \sqrt{1-e^2} \sin I}
\frac{\partial {\cal R}}{ \partial I} ,
\label{dOmegadt}
\\
\frac{d \omega}{dt} & = & \frac{1}{n a^2} \left(
\frac{\sqrt{1-e^2}}{e} \frac{\partial {\cal R}}{ \partial e}
- \frac{\cot I}{\sqrt{1-e^2}} \frac{\partial {\cal R}}{ \partial I} \right),
\label{dwdt}
\\
\frac{d I}{dt} & = & \frac{ 1}{n a^2 \sqrt{1-e^2}} \left(
\cot I \frac{\partial {\cal R}}{ \partial \omega} -
\frac{1}{\sin I} \frac{\partial {\cal R}}{ \partial \Omega} \right),
\label{dIdt}
\end{eqnarray}

\noindent where ${\cal R} = - \left< \Phi_d \right>$ is the disturbing
function and $n=\sqrt{G(M_{\star} + M_p)/a^3}$ is the planet's mean motion.
From these equations and the expression~(\ref{potential3}) of
$\Phi_d$, a natural timescale emerges:

\begin{eqnarray}
T = n^{-1} \frac{M_{\star} + M_p}{ \Sigma_0 a^2} \frac{R_o}{a}.
\label{T}
\end{eqnarray}

\noindent Note that, when the evolution of the orbital elements is
periodic, in general $T$ is not the period of these oscillations (see
section~\ref{sec:solutions}).  However, equation~(\ref{T}) is
important as it gives the scaling of the evolutionary timescale with
the different parameters.  We see in particular that $T \propto
a^{-3/2} M_d^{-1}$.

Having expanded ${\cal R}$ to fourth order in $e$, we obtain $(\partial
{\cal R} / \partial e ) /e$ to second order in $e$.  We insert the
time--averaged expressions of $r_p^2$, $r_p^2 \cos 2 \theta_p$,
$z_p^2$ and $r_p z_p \sin \theta_p$ given by
equations~(\ref{rp2})--(\ref{rpzpsin}) into the
potential~(\ref{potential3}).  We then perform the derivatives of
${\cal R}= - \left< \Phi_d \right>$ with respect to $e$, $\Omega$,
$\omega$ and $I$ to obtain the dimensionless form of the Lagrangian
equations of motion~(\ref{dedt})--(\ref{dIdt}) to second order in $e$:

\begin{eqnarray}
  T \frac{de}{dt} = e  & & \left\{ \left( {\cal I}_2 - {\cal I}_4 \right)
    \frac{5}{2} \sin 2 \omega \sin^2 I \right. \nonumber \\
  & & \left. + \frac{5 {\cal I}_3}{4} \left[ \left( 3+ \cos 2I \right)
      \sin 2 \omega
      \cos 2 \Omega + 4 \cos I \cos 2 \omega \sin 2 \Omega \right] \right.
  \nonumber \\
  & & \left.  -
    \frac{5 {\cal I}_5}{2} \sin I \left( \cos I \sin 2 \omega \cos \Omega
      + \cos 2 \omega \sin \Omega \right)   \right\} ,
\label{Tdedt}
\end{eqnarray}

\begin{eqnarray}
  T \frac{d \Omega}{dt} = & & \left( - {\cal I}_2+{\cal I}_4 \right)
  \frac{\cos I}{2} \left( 2 + 4 e^2 -5 e^2
    \cos 2 \omega \right) \nonumber \\ & + & {\cal I}_3
  \left[ \frac{1}{2} \cos I
    \left( 2 + 4 e^2 - 5 e^2 \cos 2 \omega \right) \cos 2 \Omega
    + \frac{5}{2} e^2 \sin 2 \omega \sin 2 \Omega \right]
  \nonumber \\
  & + & \frac{{\cal I}_5}{4 \sin I} \left[ \left( 2 + 4 e^2 - 5 e^2
\cos 2 \omega \right)
    \cos 2 I \cos \Omega  + 5 e^2 \cos I \sin 2 \omega \sin \Omega \right],
\label{TdOmegadt}
\end{eqnarray}

\begin{eqnarray}
  T \frac{d \omega}{dt} = & &
  \left( {\cal I}_2 - {\cal I}_4 \right) \frac{1}{8} \left[
    5 \left( 2 -3e^2 \right) \cos 2 \omega +
    10 \left( 2 + e^2 \right)
    \cos 2 I \sin^2 \omega \right] \nonumber \\
  & + &  \frac{ {\cal I}_2}{8} \left( 22-e^2 \right)+ \frac{ {\cal I}_4}{8}
  \left( 2-11e^2 \right)
  \nonumber \\
  & + & {\cal I}_3 \left\{ \frac{1}{8} \cos 2 \Omega \left[2 -11e^2
     +5 \left( 6 - e^2 \right)  \cos 2 \omega
      - 10 \left( 2 + e^2 \right) \cos 2I \sin^2 \omega \right]  \right.
  \nonumber \\
  & & \left. \hspace{9.cm} -  5
    \cos I \sin 2 \omega \sin 2 \Omega \right\} \nonumber \\
  & + & \frac{{\cal I}_5}{8 \sin I} \left\{
    \cos \Omega \cos I
    \left[ \left( -2 + e^2  \right) \left( -3 + 5 \cos 2 \omega \right)
      -10 \left( 2+ e^2  \right) \cos 2I \sin^2 \omega \right]
  \right.
  \nonumber \\
  &  & \left. \hspace{6.8cm} + 10 \left( 1 - e^2 - \cos 2I \right)
    \sin 2 \omega \sin \Omega
  \right\}
\label{Tdwdt}
\end{eqnarray}

\begin{eqnarray}
  T \frac{d I}{dt} = & & \left( -{\cal I}_2 + {\cal I}_4 \right) \frac{5}{4}
e^2 \sin 2I  \sin 2 \omega \nonumber \\
& + & {\cal I}_3 \frac{ \sin I}{2} \left[
5 e^2 \cos I \cos 2 \Omega \sin 2 \omega
+ \left( 2 + 4 e^2 + 5 e^2 \cos 2 \omega  \right) \sin 2 \Omega \right]
\nonumber \\
& + & \frac{{\cal I}_5 }{4} \left[ 5 e^2 \cos 2I \cos \Omega \sin 2 \omega
+ \cos I \left( 2 + 4e^2 + 5 e^2 \cos 2 \omega \right) \sin \Omega \right]
\label{TdIdt}
\end{eqnarray}

The perturbing function depends on the third Delaunay variable
$\Omega$.  Thus the Delaunay momentum, which is $\propto \sqrt{1-e^2}
\cos I$, is generally not constant.  This is in contrast to the
classical Kozai case (e.g., Innanen et al. 1997 and see
section~\ref{sec:Kozai} below).

\section{Solution of the secular Lagrangian equations: orbital evolution}
\label{sec:solutions}

We solve the set of coupled differential
equations~(\ref{Tdedt})--(\ref{TdIdt}) to obtain the evolution of the
orbital elements when the planet is subject to the secular
gravitational perturbation by the disc.  We first focus on the case of
a flat disc and then investigate the case of a finite warp.

\subsection{Case of a flat disc}
\label{sec:Kozai}

The flat disc case is obtained by taking $\tan \alpha=0$ in
equation~(\ref{Ho}).  Therefore, $m=0$ in
equations~(\ref{F1})--(\ref{F4}).  Since $F_1(0)=F_2(0)=F_3(0)= \pi
/2$ and $F_4(0)=0$, equations~(\ref{I2})--(\ref{I5}) give:

\begin{eqnarray} {\cal I}_2 = \frac{\pi}{2} \int_{\rho_i}^1 \frac{ \Sigma
  (\rho)}{\Sigma_0} \rho^{-2} d \rho = \frac{\pi}{2(q+1)} \left(
  \rho_i^{-q-1} - 1 \right),
\label{I2F}
\end{eqnarray}

\noindent ${\cal I}_3 = {\cal I}_5 = 0$ and ${\cal I}_4 = -2 {\cal
  I}_2$.  In this limit, the potential given by
equation~(\ref{potential3}) reduces to the following expression:

\begin{eqnarray} \left< \Phi_d \right> = - \frac{G \Sigma_0}{R_o} \left[ {\cal I}_1
  R_o^2 + 2 {\cal I}_2 \left( \frac{1}{2} \left< r_p^2 \right> -
    \left< z_p^2 \right> \right) \right], \end{eqnarray}

\noindent which is of the same form as in the classical Kozai case
(Kozai 1962, Lidov 1962), in which an inner body is perturbed by a
distant companion on an inclined orbit.  When only the secular terms
are considered, the Kozai case is equivalent to averaging the mass of
the outer companion over its orbit.  As pointed out by Terquem \&
Ajmia (2010), this is analogous to having the inner body perturbed by
a disc, as long as most of the disc's mass is beyond the orbit of the
body.  We therefore expect that the perturbing potential will have the
same form in both cases.  Note that here $\left< \Phi_d \right>$ does
not depend on the third Delaunay variable $\Omega$, and therefore
$\sqrt{1-e^2} \cos I$ is constant.

With the above values of the ${\cal I}_i$, the Lagragian equations of
motion~(\ref{Tdedt})--(\ref{TdIdt}) reduce to:

\begin{eqnarray}
T \frac{de}{dt} & = &   \frac{15}{2} {\cal I}_2 e \sin 2 \omega \sin^2 I ,
\label{TdedtF}
\\
T \frac{d \Omega}{dt} & = & -  \frac{3}{2} {\cal I}_2 \cos I \left(
2 + 4 e^2 - 5 e^2 \cos 2 \omega \right),
\label{TdOmdtF}
\\
T \frac{d \omega}{dt} & = & \frac{3}{8} {\cal I}_2
\left[ 6 + 7 e^2 + 5 \left( 2 - 3 e^2 \right) \cos 2 \omega +
   10 \left( 2 + e^2 \right) \cos 2I \sin^2 \omega  \right],
\label{TdwdtF}
\\
T \frac{dI}{dt} & = &  -  \frac{15}{4} {\cal I}_2  e^2 \sin 2 \omega \sin 2I .
\label{TdIdtF}
\end{eqnarray}

\noindent To second order in $e$, these
equations are the same as equations~(5) from Innanen et al. (1997),
which were derived in the classical Kozai case, provided that we correct an
error in their expression for $d \omega / d \tau$, in which the factor
$\sqrt{1-e^2}$ should be in the denominator, not in the numerator.

Figure~\ref{fig1} shows $e$ and $I$ versus time, obtained by solving
equations~(\ref{TdedtF})--(\ref{TdIdtF}) with the following initial
values: $e_0=10^{-3}$, $\Omega_0=\omega_0=0$ and $I_0=55^{\circ}$.  (In
this section, we only aim at showing that the equations we use reduce
to the classical Kozai case in the appropriate limit, and thefore we
do not address the issue of how a planet could start on an inclined
orbit.)  We have fixed $q=1/2$ in the expression~(\ref{sigma}) of
$\Sigma$, $R_i=60$~au, $R_o=80$~au, $M_d =10^{-2}$~M$_{\odot}$,
$a=10$~au and $M_p =10^{-3}$~M$_{\odot}$.  For these parameters,
equation~(\ref{I2F}) gives ${\cal I}_2=0.565$ and the
timescale~(\ref{T}) is $T=3.8 \times 10^5$~years for
$M_{\star}=1$~M$_\odot$.  As expected, since $I_0$ is larger than the
critical angle $I_c= \arccos \sqrt{3/5} =39^{\circ}.23$, $e$ and $I$
undergo antiphase oscillations (Kozai 1962, Innanen et al. 1997).
Note that $T$ is not the period of these oscillations.  The
inclination angle varies between $I_c$ and $I_0$.  In theory, as
$\sqrt{1-e^2} \cos I$ is constant, the maximum eccentricity should be
$\sqrt{1-5 \cos^2 I_0/3} = 0.67$.  We see in figure~\ref{fig1} that
the maximum eccentricity is actually a bit larger.  This is because
when $e$ reaches the largest values, the analysis above, which is
valid only to second order in $e$, breaks down and $\sqrt{1-e^2} \cos
I$ is not constant (it varies by about 14\% is this particular case).

\begin{figure}
\begin{center}
\includegraphics[scale=0.7]{CT_figure1.ps}
\end{center}
\caption{Eccentricity $e$ (solid line) and inclination angle $I$ (in
  degrees, dotted line) versus time $t$ in Myr for $\alpha=0$ (flat
  disc), $a=10$~au, $M_p =10^{-3}$~M$_{\odot}$, $M_d
  =10^{-2}$~M$_{\odot}$, $R_i = 60$~au, $R_o = 80$~au, $e_0=10^{-3}$
  and $I_0=55^{\circ}$.  This illustrates the classical Kozai case, in
  which $e$ and $I$ undergo antiphase oscillations.  }
\label{fig1}
\end{figure}

As noted, $p=1$ in equation~(\ref{zwarp}) corresponds to a disc that
is rigidly tilted with respect to the $(x,y)$--plane.  In this case,
the term proportional to $\sin \theta$ in the denominator of the
integral over $\theta$ in the potential~(\ref{potential}) can be
eliminated by the appropriate choice of coordinates, and this reduces
to the Kozai case.  Here however, we have defined the orbital elements
with respect to the $(x,y)$--plane and {\em not} to the plane in which
the disc lies when $p=1$, so we cannot simply recover Kozai's results
when $p=1$.

\subsection{Case of a warped disc}

To investigate the warped disc case ($\alpha \ne 0$), we solve
equations~(\ref{Tdedt})--(\ref{TdIdt}).  Here again we fix $q=1/2$ in
the expression~(\ref{sigma}) of $\Sigma$ and we adopt $p=2$ or 3 in
equation~(\ref{zwarp}).  In principle, the analysis in this section is
valid only for $a \ll R_i$.  However, as we will see below, some of
the results in this section do not depend on $a$, and as it happens, any
value of $a < R_i$ could be considered.

When the disc reduces to a narrow ring ($R_o$ only slightly larger
than $R_i$), we expect the orbital angular momentum vector ${\bf
  L}_{\rm orb}$ to precess around the total angular momentum vector.
Since the ring's angular momentum is large compared to that of the
orbit, this means that ${\bf L}_{\rm orb}$ precesses around the ring's
angular momentum vector ${\bf L}_{\rm ring}$.  As the ring is inclined
by an angle $\alpha$ with respect to the $(x,y)$--plane, the
inclination $I$ of the orbit varies periodically between its initial
value $I_0$ and $I_{\rm max}=2 \alpha-I_0$ during the precession,
independently of $a$ and $M_p$.

This is illustrated in figure~\ref{fig2}, which shows $I$, $\omega$
and $\Omega$ versus time for the initial values $e_0=10^{-3}$,
$I_0=10^{\circ}$ and $\Omega_0=\omega_0=0$, and for $a=10$~au, $M_p
=10^{-3}$~M$_{\odot}$, $\alpha=30^{\circ}$, $p=3$, $M_d
=10^{-2}$~M$_{\odot}$, $R_i = 60$~au and $R_o = 61$~au.  We have
chosen a finite value of $I_0$ because, as equations~(\ref{TdOmegadt})
and~(\ref{Tdwdt}) show, $d\Omega/dt$ and $d\omega/dt$ become infinite
when $I \rightarrow 0$.  However, we have checked that exactly the
same variations of $\Omega$ and $I$ are obtained for $I_0=1^{\circ}$,
for example.  In this latter case, only the variations of $\omega$ are a bit
more irregular around the times when $I$ is minimum.  As we had
anticipated, $I_{\rm max}=50^{\circ}= 2 \alpha-I_0$.  Also, as
expected for a precessing motion, $\omega$ circulates uniformly.  As
can be seen from equation~(\ref{Tdedt}), this implies that changes in
the eccentricity due to the torque exerted by the disc average to zero
over a precession period, so that no secular change of the
eccentricity occurs.

Note that exactly the same curves are obtained for different values of
$a$, $M_p$ or $M_d$.  Only the precession period $T_{\rm prec} \propto
a^{-3/2} M_d^{-1}$ changes when these parameters are varied.

\begin{figure}
\begin{center}
\includegraphics[scale=0.7]{CT_figure2.ps}
\end{center}
\caption{Inclination angle $I$ (solid line), longitude of ascending
  node $\Omega$ (dotted line) and argument of pericentre $\omega$
  (dashed line) versus time $t$ in Myr for $\alpha=30^{\circ}$,
  $p=3$, $a=10$~au, $M_p =10^{-3}$~M$_{\odot}$, $M_d
  =10^{-2}$~M$_{\odot}$, $R_i = 60$~au, $R_o = 61$~au, $e_0=10^{-3}$
  and $I_0=10^{\circ}$.  The orbital angular momentum vector precesses
  around the ring's angular momentum vector, which causes $I$ to
  oscillate bewteen 10 and $50^{\circ}$ and $\omega$ to circulate.}
\label{fig2}
\end{figure}

In the case displayed in figure~\ref{fig2}, the ring is narrow and
essentially all of its mass is in a plane inclined by $\alpha$ with
respect to the $(x,y)$--plane.  If we now increase the ratio
$R_o/R_i$, while keeping $M_d$ fixed, more of the disc's mass will be
at smaller inclination with respect to the $(x,y)$--plane, and
therefore the angle between ${\bf L}_{\rm ring}$ and ${\bf L}_{\rm
  orb}$ will be smaller.  The maximum value reached by $I$ during the
precession will accordingly be reduced.  For the same parameters as
above but with $R_o=80$ or 100~au, we obtain $I_{\rm max}=36^{\circ}$ or
$27^{\circ}$, respectively.  If we keep $R_0=61$~au but take
$R_i=50$~au, $I_{\rm max}$ is reduced to $40^{\circ}$.  These values
of $I_{\rm max}$ would be higher for $p=2$, as in that case the inner
parts of the disc would have a higher inclination (see table below).

From equations~(\ref{zwarp}) and~(\ref{Ho}), we see that the disc can
be modelled as a collection of rings of width $dr$ located at a
distance $r$ from the central star and inclined with respect to the
$(x,y)$--plane by an angle:

\begin{eqnarray} \beta(r)=\arctan \left[ \left( \frac{r}{R_o} \right)^{p-1} \tan
  \alpha \right]. \end{eqnarray}

\noindent The planet will be more affected by rings which are more
massive and closer.  We therefore define an average inclination angle
$\overline{\beta}$ over the whole disc by weighting $\beta(r)$ by the
mass of the ring and a function $f(r)$ that decreases with $r$ and
such that $f(R_i)=1$:

\begin{eqnarray}
\overline{\beta} = \frac{1}{M_d} \int_{R_i}^{R_o} 2 \pi r dr \Sigma(r) f(r)
\beta(r).
\end{eqnarray}

\noindent We would expect the orbital angular momentum vector to
precess around the vector inclined by this angle with respect to the
vertical axis, and therefore $I_{\rm max}= 2 \overline{\beta} - I_0$.
We have found that this is the case to within a good approximation if
we take $f(r)=\sqrt{R_i/r}$, which we will adopt thereafter.  For the
parameters used above, with $R_i=60$~au and $R_o=80$ or 100~au, we obtain
$\overline{\beta}=22^{\circ}$ or $18^{\circ}$, respectively, and
therefore $2 \overline{\beta} - I_0= 34^{\circ}$ or $26^{\circ}$,
reasonably close to the values of $I_{\rm max}$ which were found to be
$36^{\circ}$ and $27^{\circ}$.

Table~\ref{table} gives $I_{\rm max}$ and $\overline{\beta}$ for
different values of $\alpha$, $R_i$ and $R_o$ and for $p=2$ or 3.  In
this table, we have considered parameters for which the interaction
between the disc and the planet leads to a precession of the orbit.
As noted above, $I_{\rm max}$ does not depend on $a$, and therefore
the results listed here are valid for any $a < R_i$, not just $a \ll
R_i$.  This has been checked by doing numerical simulations using the
exact disc potential, i.e. without assuming that $a$ is very small
compared to $R_i$ (see section~\ref{sec:simulations}).

We use equation~(\ref{T}) to write the precession period $T_{\rm
  prec}$ as a function of $a$ and $M_d$ and define $\overline{T}_{\rm
  prec}$ such that:

\begin{eqnarray}
 T_{\rm prec} = \overline{T}_{\rm prec} \left(
\frac{ 1 \; {\rm au}}{a} \right)^{3/2}
\frac{ 10 \; {\rm M}_{\rm J}}{M_d} .
\end{eqnarray}

\noindent
$\overline{T}_{\rm prec}$ is given in table~\ref{table} for the cases
considered and is equal to the precession period for $a=1$~au and
$M_d=10$~$M_{\rm J}$.  Note that this precession period is divided by about
10 if we take $a=5$~au.  Neither $I_{\rm max}$ nor $T_{\rm prec}$
depend on $M_p$.  For the initial inclination angle $I_0=1^{\circ}$
used here, we expect $I_{\rm max} = 2 \overline{\beta}$.  We note that
this is approximately the case, except when $R_o/R_i$ is large and
$p=3$, which corresponds to a disc with very low inclinations in the
inner parts.

\begin{table}
  \caption{Maximum inclination and precession period
in the regime where the orbit precesses for $p=3$ and $p=2$.}
\label{table}
\begin{tabular}{ccccccccccc}
  \cline{5-7}
\cline{9-11}
  & & & & \multicolumn{3}{c}{p=3} & & \multicolumn{3}{c}{p=2} \\
  \hline
  $\alpha$ &  $R_i$  & $R_o$ & &  $\overline{\beta}$ & $I_{\rm max}$ &
  $\overline{T}_{\rm prec}$ & &  $\overline{\beta}$ & $I_{\rm max}$ &
  $\overline{T}_{\rm prec}$ \\
  (degrees) & (au) & (au)  & & (degrees) & (degrees) & (Myr) & &
  (degrees) & (degrees) & (Myr) \\
  \hline
  30  & 60  & 61  & & 29.5 & 59 & 36  & & 29.5 & 59 & 36  \\
  --- & --- & 80  & &  22  & 45 & 50  & & 25   & 52 & 52  \\
  --- & --- & 100 & &  18  & 36 & 68  & & 21.5 & 46 & 70  \\
  --- & --- & 150 & & 12.5 & 23 & 120 & & 16.5 & 36 & 124  \\
  --- & --- & 200 & &  10  & 16 & 180 & & 14   & 30 & 187  \\
  30  & 15  & 100 & & 6.5  & 6  & 7.8 & & 9.5  & 18 & 8.0  \\
  50  & --- & --- & & 12   & 12 & 8.1 & & 17.5 & 35 & 8.6  \\
  70  & --- & --- & & 21   & 24 & 8.8 & & 28.5 & 68 & 11  \\
  \hline
\end{tabular} \\
\medskip
The initial conditions are $I_0=1^{\circ}$, $e_0=10^{-3}$ and
$\Omega_0=\omega_0=0$.  $I_{\rm max}$ is independant of $M_d$, $M_p$ and
$a$ as long as $a<R_i$.  The precession period is $T_{\rm prec}=
\overline{T}_{\rm prec} \times (1 \; {\rm au} / a)^{3/2}
(10 \; {\rm M}_{\rm J} / M_d)$.
\end{table}

When $I_{\rm max}$ gets larger than $90^{\circ}$, the orbit becomes
retrograde (more precisely, it would become retrograde relative to the
stellar spin if stellar rotation were taken into account).  In the
case of the narrow ring and for $I_0=10^{\circ}$, this happens for
$\alpha > 50^{\circ}$.  For $\alpha=49^{\circ}$, $I_{\rm
  max}=88^{\circ}$ and $e$ stays small (smaller than $3 \times
10^{-3}$).  In that case, $\omega$ still circulates, although not as
smoothly as for smaller values of $\alpha$.  For $\alpha=55^{\circ}$,
$I_{\rm max}=100^{\circ}$ and $e$ grows to 0.23.  This is illustrated
in figure~\ref{fig3}, which shows $I$, $\omega$, $\Omega$ and $e$ for
the same parameters as figure~\ref{fig2} except for
$\alpha=55^{\circ}$.  At the beginning of the evolution, $\omega$
increases .  But as $I$ approaches $90^{\circ}$, $\omega$ no longer
circulates, so that changes in $e$ no longer average to zero over time
and a net increase of the eccentricity is observed.

\begin{figure}
\begin{center}
\includegraphics[scale=0.7]{CT_figure3.ps}
\end{center}
\caption{Inclination angle $I$ (solid line, upper plot), longitude of
  ascending node $\Omega$ (dotted line, upper plot), argument of
  pericentre $\omega$ (dashed line, upper plot) and eccentricity $e$
  (lower plot) versus time $t$ in Myr for the same parameters as
  figure~\ref{fig2} except for $\alpha=55^{\circ}$.  When $I$
  approaches $90^{\circ}$, $\omega$ no longer circulates and $e$
  grows.  }
\label{fig3}
\end{figure}

As $\alpha$ is increased, $e$ becomes too large for the second order
analysis of this section to apply.

\section{Numerical simulations}
\label{sec:simulations}

The analysis presented in the above section is valid only for small
values of $e$.  In principle, it is also valid only for $a \ll R_i$.
However, as pointed out in the previous section, when the planetary
orbit precesses, the variations of $I$ do not depend on $a$ and
therefore any value of $a$ may be considered in the analysis.  To
study the evolution of the planet's orbit in a more general case, we
use the $N$--body code described in Papaloizou~\& Terquem~(2001) in
which we have added the gravitational force exerted by the disc onto
the planet.

The equation of motion for the planet is:
\begin{equation}
{d^2 {\bf r} \over dt^2} = -{GM_\star{\bf r} \over |{\bf r}|^3} -
\mbox{\boldmath $\nabla$} \Phi_d - {GM_p{\bf r} \over |{\bf r}|^3} + {\bf
\Gamma}_{t,r} \; , \label{emot}
\end{equation}

\noindent where ${\bf r}$ is the position vector of the planet and
$\Phi_d$ is the gravitational potential of the disc given by
equation~(\ref{potential}).  Note that here we use the exact disc
potential, so that the results of this section are valid in the
general case $a<R_i$, not just $a \ll R_i$.  The third term on the
right--hand side is the acceleration of the coordinate system based on
the central star.  This neglects the net force of the disc on the
central star, which would of course require a calculation of the disc
response to the planet's perturbation.  Tides raised by the star in
the planet and relativistic effects are included through ${\bf
  \Gamma}_{t,r}$, but they are unimportant here as the planet does not
approach the star closely.  Equation~(\ref{emot}) is integrated using
the Bulirsch--Stoer method and the integrals involved in
$\mbox{\boldmath $\nabla$} \Phi_d$ are calculated with the Romberg
method (Press et al.  1993).

The planet is set on a circular orbit in the $(x,y)$--plane at a
distance $a$ from the star.  As here there is no frictional force
acting on the planet, the orbital energy is conserved and the
semimajor axis $a$ is constant throughout the evolution of the system.
In this paper, we suppose that the planet's orbit is inside the disc
inner cavity ($a<R_i$).  If the orbit were still embedded in the disc,
the frictional force exerted by the disc on the planet would damp the
planet's inclination and eccentricity, and this would have to be taken
into account.

In the simulations reported here, we have taken $q=1/2$ in the
expression~(\ref{sigma}) of $\Sigma$, $p=3$ in equation~(\ref{zwarp})
and $M_{\star}=1$~M$_{\odot}$, $M_p =10^{-3}$~M$_{\odot}$ and $M_d
=10^{-2}$~M$_{\odot}$.

Figure~(\ref{fig4}) shows the evolution of $e$ and $I$ for
$\alpha=70^{\circ}$, $a=10$~au, $R_i = 60$~au and $R_o = 80$~au.  Note
that exactly the same evolution is obtained for different values of
$M_p$.  The orbit becomes retrograde and the eccentricity is pumped up
to very high values.  We have stopped the calculation when, due to the
very high eccentricity, the code failed to conserve the energy to
within a specified accuracy .  As $e$ approaches unity, we expect
tidal interaction with the central star to become important and
dissipate orbital energy so that $e$ decreases.  Ultimately, if the
disc does not disappear, the orbit would become circularized on a
short orbit.  We have verified that we obtain a similar evolution by
solving equations~(\ref{Tdedt})--(\ref{TdIdt}) in the early stages
when $e \ll 1$.

\begin{figure}
\begin{center}
\includegraphics[scale=0.7]{CT_figure4.ps}
\end{center}
\caption{Eccentricity $e$ (solid line) and inclination angle $I$ (in
  degrees, dotted line) versus time $t$ in Myr for
  $\alpha=70^{\circ}$, $a=10$~au, $M_p =10^{-3}$~M$_{\odot}$, $M_d
  =10^{-2}$~M$_{\odot}$, $R_i = 60$~au and $R_o = 80$~au. The orbit
  becomes retrograde and the eccentricity is increased to very high
  values.  }
\label{fig4}
\end{figure}

We observe exactly the same variations of $I$ and $e$ if we take
$a=20$~au instead of $a=10$~au.  The evolution is about three
times faster for this wider orbit, which is consistent with an
evolutionary timescale propotional to $a^{-3/2}$.  Similar
behaviour is observed for $p=2$ (instead of $p=3$ used here), also on
a shorter timescale.

Similar evolution is observed for other values of $\alpha$.  The case
displayed in figure~\ref{fig4} is typical, illustrating
what happens when the inclination of the orbit passes through
$90^{\circ}$: the eccentricity is pumped up to high values and $I$
varies over a broad range of values, including retrograde
inclinations.

\section{Discussion}
\label{sec:discussion}

In this paper, we have studied the secular interaction between a
planet and a warped disc.  We have considered the stage of evolution
after the planet has cleared--out a gap.  Therefore, there is no
friction between the disc and the planet, and the planetary orbit
evolves only under the effect of the gravitational potential from the
disc.  In general, we find that the interaction between the planet and
the disc leads to the precession of the orbit.  The orbital plane
therefore becomes inclined relative to the plane which is tangent to
the disc's inner parts, which we will hereafter call the ``reference
plane''.  If we define some average inclination of the disc with
respect to the reference plane, then the inclination of the orbital
plane is twice this averaged angle.  The precession of the orbit does
not modify its eccentricity: if the planet started on a circular
orbit, the eccentricity stays small.

There is an instability when the inclination of the orbital plane
relative to the reference plane approaches 90$^{\circ}$.  The
eccentricity of the planet is then built up and both the inclination
and the eccentricity develop very large amplitude variations.  That
happens when the disc is very severely warped, or if there is a
significant amount of mass in a ring inclined by a least 45$^{\circ}$
relative to the initial orbital plane.

Note that the inclination of the orbital plane is measured relative to
the disc's inner region.  The disc as a whole will in general precess
around the central star's equatorial plane, and this will have to be
taken into account to compute the inclination of the orbital plane
relative to this equatorial plane.

The process studied in this paper shows that the orbit of a planet
that has formed in a warped disc could develop a significant
inclination when there is no frictional damping.

Although the inclination of the orbital plane does not depend on the
planet's semimajor axis $a$, the timescale for generating a
significant inclination does increase as $a$ decreases.  For a disc of
$\sim 10$~$M_{\rm J}$, we typically find that for this timescale to be at
most of a few Myr, $a$ has to be larger than 5~au.  We therefore
expect only planets beyond the snow line to become significantly
misaligned through this process.  If the warp is such that the orbital
plane precesses, the eccentricity is not built up and the planet would
be left on its original semimajor axis as the disc dissipates.
However, if large eccentricities are generated, the pericentre would
become small enough for tidal interaction with the star to circularize
the orbit.  The planet would then end up on a short orbit, as are the
misaligned planets observed so far.

The process described in this paper therefore leads to two distinct
populations of misaligned planets: one with objects beyond the snow
line and at most moderate eccentricities, and another with objects on
short and circular orbits.  Of course, the presence of other planets in
the system would modify this picture as planet--planet interactions
would then change the orbital elements.  So far, only planets on short
orbits have been detected on inclined orbits.  Measurements of the
projected inclination angle for more distant objects would help to put
constraints on the process discussed here.

\section*{Acknowledgements}

I would like to thank J. Silk for suggesting that the interaction
between a planet and a warped disc may drive the inclination of the
planetary orbit.  It is a pleasure to thank S. Balbus for a critical
reading of the manuscript and useful comments.


\section{Introduction}  In this work we consider finite, undirected, simple graphs $G=(V,E)$ with $n$ vertices and $q$ edges.
The \emph{neighbourhood} of a vertex $v\in V(G)$ is the set $N_{G}(v)$ of all the vertices adjacent to $v$ in $G$.
For a set $X\subseteq V(G),$ \emph{the open neighbourhood} $N_{G}(X)$ is defined to be $\bigcup_{v\in X}N_{G}(v)$ and
\emph{the closed neighbourhood} $N_{G}[X]=N_{G}(X)\cup X.$

The \emph{degree} $d_{G}(v)$ of a vertex $v\in V(G)$ is the number of
edges incident to $v;$ $d_{G}(v)=|N_{G}(v)|.$  If $d_G(v)=n-1,$ then $v$ is a \emph{universal vertex} of $G$ and we call $v$ a \emph{semi-universal vertex} if every vertex of $G$ is either a neighbour of $v$ or has a common neighbour with $v.$

If $d_G(v)=1,$ then $v$ is an \emph{end-vertex} of $G.$ We denote by $\Omega(G)$ the set of end vertices of $G$. The neighbour of an end vertex is called a \emph{support vertex.} Support vertex which has more than one end-vertex as a neighbour is a \emph{strong support} vertex. For two sets $X, Y\subset V(G),$ we denote by $E(X,Y)$  the set of edges $uv\in E(G)$ such that $u\in X$ and $ v\in Y.$

For $X \subseteq  V (G)$  and $ x \in X,$  the set
$PN(x,X) =  N_G[x]- N_G[X -\{x\}]$
is called the \emph{private neighborhood} of $x$ with respect to $X$.
An element $u$ of $PN(x, X)$ is called a \emph{private neighbour of x relative to X} or a $X$-private neighbour of $x.$ An  $X$-private neighbour of $x$ is either  $x$ itself, in wich case $x$ is an isolate vertex of $G[X]$, or is a neighbour of $x$ in $G$ which is not adjacent to any vertex of $X$. This latter type will be called an  $X$-\emph{external private neighbour}  of $x$.

A subset $D \subseteq V(G)$  is a \emph{dominating set} if for every vertex $u \in V(G)-D$, there exists $v \in D$ such that $u$ and $v$ are adjacent.  The minimum cardinality of a dominating set in $G$ is  the \emph{domination number}  of $G$ and is denoted by $\gamma(G)$.

It is well-known (Chapter 1,\cite{b4}) that a dominating set $D$ is minimal if and only if for any vertex $v \in  D$ there exists a private neighbourd of $v$ with respect $D.$  This fact was our motivation for the next definition.

\begin{df}\label{SDS}
A subset $D$ of $V(G)$ is called a super dominating set if for every vertex $v \in V(G)-D$ there exists an external private neighbour of $v$ with respect to $V(G)-D.$
\end{df}

The minimum cardinality of a super dominating set in $G$ is called
the \emph{super domination number} of $G$ and is denoted by
$\gamma_{sp}(G).$ A $\gamma_{sp}(G)-$set is a super dominating set
of $G$ with cardinality $\gamma_{sp}(G) .$
\begin{ex}
We can consider a super dominating set $D$ as a set of workers and $V(G)-D$ as a set of managers. Two workers (resp. manager) can be related or not. A worker has category A if it is adjacent to one and only one manager. By definition of superdominating set, each manager in $V(G)-D$ is adjacent at least one worker with category A.   \\

In the same way, we can consider a super dominating set $D$ as a set of students and $V(G)-D$ as a set of professors. Two students (resp. professors) can be related or not. A student is a Ph.D students if it has relation with one and only one professor. By the definition of superdomination set, each professor in  $V(G)-D$ has at least one Ph. D students.   \\

\end{ex}

The undefined terms in this paper may be found in \cite{b1}, \cite{b4}.

\begin{ex}
For $n\geq 2,$ the super domination number of the complete graph $K_n$ is $\gamma_{sp}(K_{n})=n-1$ and also for the star $K_{1,n-1}$ is $\gamma_{sp}(K_{1,n-1})=n-1$.\\

For a complete bipartite graph $K_{m,n}$ with $min\{ m,n\} \geq 2,$ the super domination number is $\gamma_{sp}(K_{m,n})=m+n-2 $.

\end{ex}

\begin{obsn} A set $D \subseteq  V(G)$ is a super dominating set of $G$ if and only if   for every $v \in V(G)-D$ there exists $u \in N_G(v) \cap D$ such that $N_G(u) \subseteq D \cup \{v\}$.
\end{obsn}

\begin{df}\label{SDPS}
A super dominating set $D$ is a perfect set if $E(D, V(G)-D)$ is a perfect matching.
\end{df}

\section{Preliminary results}

In this section we give some  observations about the bounds on the super domination number.

\begin{obsn}\label{LowerB}
For any graph $G$, $\gamma_{sp}(G)\geq \frac{n}{2}$ since
for  any $\gamma_{sp}(G)-$set $D,$ we have that $|D|\geq |V(G)-D|$.
\end{obsn}

The next observation follows immediately from the Definition $\ref{SDS}$.

\begin{obsn}\label{TrivialB}
For any graph $G$, the super domination number $\gamma_{sp}(G)=1$ if and only if $G\cong K_1$ or $G\cong K_2$ and $\gamma_{sp}(G)=n$ if and only if $G=\overline{K_n}.$
\end{obsn}

From the above  observations we have that for  any connected graph $G$ is $$\frac{n}{2}\leq \gamma_{sp}(G) \leq n-1.$$

The following theorem gives us a characterization of the connected graphs with $\gamma_{sp}(G)=\frac{n}{2}.$

\begin{thm}\label{twr1} For any connected graph $G$,  $\gamma_{sp}(G)=\frac{n}{2}$ if and only if every minimum super dominating set is a perfect set.

\end{thm}

\begin{proof}   Let $D$ be a $\gamma_{sp}-$set of $G$. If $D$ is a perfect set, then $|D|=|V(G)-D|$ and $\gamma_{sp}(G)=\frac{n}{2}.$ \\

Conversely, suppose $\gamma_{sp}(G)=\frac{n}{2}.$ Let $D$ be a minimum super dominating set of $G.$ If there exists a vertex $z\in D$ such that $|N(z)\cap (V(G)-D)|>1,$ then  from  Definition \ref{SDS} we obtain $|D|>|V(G)-D|,$ a contradiction. If there exists a vertex $z\in D$ such that $|N(z)\cap (V(G)-D)|=0$, then, since $|D|=|V(G)-D|,$ there is a vertex $w\in D$ such that $w$ has more than one neighbour in $V(G)-D,$ a contradiction. Thus every vertex of $D$ has exactly one neighbour in $V(G)-D.$

Suppose there is a vertex $a\in V(G)-D$ such that $|N(a)\cap D|\geq 2.$ Then either $|D|>\frac{n}{2}$ or there exists a vertex from $D$ which has more than one neighbour in $V(G)-D,$ again a contradiction. Thus every vertex of $V(G)-D$ has exactly one neighbour in $D.$ By
Definition $\ref{SDPS},$  $D$ is a perfect set.
\end{proof}

The next lemma gives us a relation between the super domination number of a connected graph $G$ and its diameter. Recall that the diameter of a graph $G$, denoted by $diam (G)$,  is defined to be the maximum distance between any two vertices $x,y\in V(G)$.

\begin{lem}
If $G$ is a connected graph with $diam (G)\geq 3,$ then $\gamma_{sp}(G) \leq n-2.$
\end{lem}

\begin{proof}
Suppose $diam(G)=k$, $k\geq 3$, $d_G(x,y)=k$ and $\alpha =(x,w_1,\cdots, w_{k-1},y)$ is any $xy-path$ with minimum length. Then $V(G)-\{ x,y\} $ is a  super dominating set of $G$ and $\gamma_{sp}(G) \leq n-2.$

\end{proof}

\begin{cor}\label{Diam2}
By  Observation $\ref{TrivialB}$ and the above Lemma, for any connected graph $G,$ if $\gamma_{sp}(G) =n-1,$ then $diam(G)\leq 2.$

\end{cor}

In particular, if $\gamma_{sp}(G) =n-1$ and $G$ is a connected graph, then for any $x\in V(G)$, $x$ is a universal vertex or $x$ is a semi-universal vertex.\\

The converse of Corollary $\ref{Diam2}$ is not true.  Consider a graph $G$ with
$V(G)=\{x_1, x_2, ..., x_k, y_1, y_2,..., y_k, u\}$ and  $E(G)=\{ x_1y_1, x_2y_2,
..., x_ky_k\} \cup \{ux_i, uy_i\}$ for  $1\leq i\leq k$ and $k\geq 2$. Then $diam(G)=2$, but
$\{x_1, x_2, ....,x_k, u\}$ is a super dominating set of $G$, so $\gamma_{sp}(G)\leq k+1<2k.$\\

The next theorem describes the super domination number of  a cycle with $n$ vertices.

\begin{thm}
For a cycle $C_n$ is

$$\gamma_{sp}(C_{n})=\left\{ \begin{array}{cl}
\left\lceil \frac{n}{2}\right\rceil   & \textrm{if $n\equiv 0,3\ (mod 4)$} \\
                                       &         \\
\left\lceil \frac{n+1}{2}\right\rceil  & \textrm{otherwise.}
\end{array} \right.$$

\end{thm}

\begin{proof} Let $V(G)=\{ u_1,u_2,...,u_n\}$ be the set of vertices of $C_4$. We consider the following four cases.

\begin{enumerate}
    \item Let $n \equiv 0(mod 4)$. Then $n=4k$ for some positive integer $k.$ Consider the set $D=\{ u_{2},u_{3},u_{6},u_{7},\cdots,u_{4k-2},u_{4k-1}\}$. Clearly, $D$ is a super dominating set with $|D|=\frac{n}{2}$. By  Observation $\ref{LowerB}$, $D$ is minimum and $\gamma_{sp}(C_{n})=\left\lceil \frac{n}{2}\right\rceil .$

    \item Let $n \equiv 1 (mod 4)$. Then $n=4k+1$ for some positive integer $k.$ If $n=5,$ then it is easy to check that $\gamma_{sp}(C_5)=3=\frac{6}{2}=\frac{n+1}{2}.$ Let $n>5$ and consider $D=\{ u_{2},u_{3},u_{6},u_{7},\cdots,u_{4k-2},u_{4k-1}, u_{4k+1}\}$. Then $D$ is a  super dominating set with $|D|=\frac{n+1}{2}$. By  Observation  $\ref{LowerB}$, $D$ is minimum and $\gamma_{sp}(C_{n})=\left\lceil \frac{n+1}{2}\right\rceil $.

    \item Let $n \equiv 2 (mod 4)$. Then $n=4k+2$ for some positive integer $k.$ Consider $D=\{ u_{2},u_{3},u_{6},u_{7},\cdots,u_{4k-1},u_{4k+2}\}$ . Then $D$ is a  super dominating set with $|D|=2k+2$. By Def \ref{SDS}, it is not possible to have a super dominating set with $2k+1$ elements. Therefore, $D$ is minimum and $\gamma_{sp}(C_{n})=\left\lceil \frac{n+1}{2}\right\rceil $.

\item Let $n \equiv 3 (mod 4)$. Then $n=4k+3 $ for some positive integer $k.$ Consider $D=\{u_{2},u_{3},u_{6},u_{7},\cdots,u_{4k+2},u_{4k+3}\}$. Clearly, $D$
is a super dominating set with $|D|=\frac{n+1}{2}$. By Observation $\ref{LowerB}$, $D$ is minimum and $\gamma_{sp}(C_{n})=\left\lceil \frac{n}{2}\right\rceil .$
\end{enumerate}
\end{proof}

\begin{cor}
For a path $P_n$ with $n\geq 3$, $\gamma_{sp}(P_{n})=\left\lceil \frac{n}{2}\right\rceil.$
\end{cor}

Now we present a lemma showing the relation between the super domination number of $G$ and the number of its edges.

\begin{lem}\label{t1} For any connected graph $G$ with $ n > 1$, $\gamma_{sp}(G) \leq 2q - n + 1$. If $\gamma_{sp}(G)=2q-n+1$, then $G$ is a tree.
\end{lem}

\begin{proof} For any connected graph $G$, $\gamma_{sp}(G) \leq n-1 = 2(n-1)-n+1 \leq 2q-n+1$.  If $\gamma_{sp}(G)=2q-n+1$, then $q = n-1$ and $G$ is a tree.
\end{proof}

\begin{thm} For any graph $G$, $\gamma_{sp}(G) \geq n - \frac{1}{2}-\sqrt{\frac{2n^2-2n-4q+1}{4}}$ and the bound is sharp.
\end{thm}

\begin{proof} Let $D$ be a $\gamma_{sp}-$set of $G$.  Since $D$ is a super dominating set, for every $u \in V(G)-D$ there exists $v \in N(u) \cap D$ such that $N(v) \subseteq D \cup \{u\}$.  For every $u \in V(G)-D$, we can find an element $v \in D$ adjacent only to $u$.  Hence, $v$ is not adjacent to $n-\gamma_{sp}-1$ elements of $V(G)-D$.  Since there are $n-\gamma_{sp}$ elements in $V(G)-D$, we can find $n-\gamma_{sp}$ vertices in $D$ such that each of  $n-\gamma_{sp}$ vertices in $D$ is not adjacent to $n-\gamma_{sp}-1$ elements in $V(G)-D$. Therefore
\begin{center}
$q \leq \frac{n(n-1)}{2} - (n-\gamma_{sp}(G))(n-\gamma_{sp}(G)-1)$
\end{center}
and after some calculations we obtain
\begin{center}
 $\gamma_{sp}(G) \geq n - \frac{1}{2}-\sqrt{\frac{2n^2-2n-4q+1}{4}}$.
\end{center}
The bound is attained for the graph $C_{4}$.
\end{proof}

We finish this section with the following lemma that gives us a relation between the super domination number of a graph $G$ and its complement $\overline{G}.$ In $1956$ the original paper \cite{NG} by Nordhaus and Gaddum appeared. In it they gave
sharp bounds on the sum and product of the chromatic numbers of a graph and its
complement. Since then such results have been given for several parameters; see
for example \cite{NG2}. Here we have similar inequalities for the super domination number.

\begin{thm} For any graph $G,$ $n\leq \gamma_{sp}(G)+\gamma_{sp}(\overline{G})\leq 2n-1.$ The equality of the upper bound holds if and only if  $\{G,\overline{G}\}\cong K_n.$
\end{thm}

\begin{proof} The lower bound follows from Observation $\ref{LowerB}$.  Since $G$ and $\overline{G}$ can not be simultaneously  isomorphic to $\overline{K_n},$ by  Observation $\ref{TrivialB}$ we have the upper bound.

If $G\cong K_n,$ then $\overline{G}\cong \overline{K_n}$ and  $\gamma_{sp}(G)+\gamma_{sp}(\overline{G})= 2n-1.$ Assume now   $\gamma_{sp}(G)+\gamma_{sp}(\overline{G})= 2n-1.$ Then, without loss of generality, $\gamma_{sp}(G)=n$ and $\gamma_{sp}(\overline{G})=n-1.$ From the Observation $\ref{TrivialB}$, $G\cong \overline{K_n}$ and $\overline{G} \cong K_n .$

The lower bound is also attained, for example if $G\cong P_4.$
\end{proof}

\section{Super domination in trees}

In this section we are interested in giving a characterization of  trees which attain the lower and the upper bounds in terms of the number of vertices. The following Theorem gives us an upper bound for  trees.

\begin{thm}\label{twr2} If $T$ is a tree with at least three vertices, then $\frac{n}{2}\leq\gamma_{sp}(T)\leq n-s,$ where $s$ is the number of support vertices in $T.$
\end{thm}

\begin{proof}
The lower bound follows from Observation \ref{LowerB}. Let $Supp$ be the set of supports in $T,$ $|Supp|=s.$ Let $S_1\subseteq \Omega(T),$ $S_1=\{v_1,\ldots,v_{|Supp|}\}$ such that every $v_i, 1\leq i \leq |Supp|$ is adjacent to a different vertex from $Supp.$ Then $V-S_1$ is a super dominating set of $T$ and $\gamma_{sp}(T)\leq n-s.$
\end{proof}

Now we characterize the extremal trees of Theorem \ref{twr2}. We begin with the following observation.

\begin{obsn}\label{drzewa} Let $T$ be a tree such that $V(T)=2m$ and $\gamma_{sp}=m$ for some $m\geq 1.$ If there are two trees $T_1, T_2$ with $V(T_1)\cup V(T_2)=V(T), V(T_1)\cap V(T_2)=\emptyset$ and $E(T_1), E(T_2)\subseteq E(T)$ such that $\gamma_{sp}(T)=\gamma_{sp}(T_1)+\gamma_{sp}(T_2),$ then $\gamma_{sp}(T_1)=\frac{|V(T_1)|}{2}$ and $\gamma_{sp}(T_2)=\frac{|V(T_2)|}{2}.$
\end{obsn}
\begin{proof} Without loss of generality, suppose $\gamma_{sp}(T_1)>\frac{|V(T_1)|}{2}.$ Then $\gamma_{sp}(T)=\gamma_{sp}(T_1)+\gamma_{sp}(T_2)>\frac{|V(T_1)|}{2}+\frac{|V(T_2)|}{2}>m.$
\end{proof}

Let $\mathcal{R}$ be the family of trees $T$ that can be obtained from a sequence $T_1,\ldots,T_j$ ($j\geq 1$) of trees such that:\\

1. The tree $T_1=P_2=(a_1,b_1).$\\

2. If $j\geq 2,$ define the tree $T_j$ such that $V(T_j)=V(T_{j-1})\cup \{a_j,b_j\}$ and $E(T_j)=E(T_{j-1})\cup \{a_jb_j\}\cup \{e\},$ where $e=a_ia_j$ or $e=b_ib_j$ for some $1\leq i\leq j-1.$\\

If $T\in \mathcal{R}$ with $|V(T)|=2m,$ denote by $A=\{a_1,\ldots, a_m\}$ and $B=\{b_1,\ldots,b_m\}.$

Observe that for all $1\leq i\leq m,$ $N_T(a_i)\cap B=\{b_i\}$ and $N_T(b_i)\cap A=\{a_i\}.$

It is clear that $A$ and $B$ are minimum super dominating sets of $T.$ Moreover, by construction of $T,$ $T$ does not have strong support vertices.

We  say that a vertex $v\in A$ has a status $a$ ($sta(v)=a$) and a vertex $v\in B$ has a status $b$ ($sta(v)=b$).

\begin{thm} Let $T$ be a tree with $|V(T)|=n\geq 2.$ Then $\gamma_{sp}(T)=\frac{n}{2}$ if and only if $T\in \mathcal{R}.$
\end{thm}
\begin{proof} If $T\in \mathcal{R},$ then by definition of the family $\mathcal{R},$  $\gamma_{sp}(T)=\frac{n}{2}.$
Let $T$ be a tree with  $\gamma_{sp}(T)=\frac{n}{2};$ then $n=2m$ for $m\geq 1$. We show that $T\in \mathcal{R}$ by induction on $m.$  If $m=1,$ then $T=P_2=(a_1,b_1)\in \mathcal{R}.$ If $m=2$ and $\gamma_{sp}(T)=2,$ then $T=P_4=(a_1, b_1, b_2, a_2)\in \mathcal{R}.$ Assume that $|V(T)|=2m\geq 6$ and the Theorem holds for any $T'$ with $|V(T')|=2m'$ with $m'<m.$

Let $P=(x_1,\ldots,x_p)$ be a longest path of $T$ and let $D$ be a minimum super dominating set of $T.$ By definition of $D,$ $|\{x_1,x_2\}\cap D|\geq 1.$ We consider the following cases. \\

$Case$ $1.$ The vertices $x_1, x_2\in D.$ Let $T'=T-\{x_2x_3\}=T_1\cup T_2,$ where $x_2\in T_1,$ $x_3\in T_2.$ Since $P$ is a longest path of $T,$ $T_1$ is a star and $\gamma_{sp}(T_1)=|V(T_1)|-1.$ \\

$a)$ Suppose $x_3\in D.$

 Then $d_T(x_2)>2$ (in the other case $D-\{x_1\}$ is a super dominating set of $T,$ a contradiction).  We have  $\gamma_{sp}(T_1)=|V(T_1)|-1>\frac{|V(T_1)|}{2}.$ Since  $\gamma_{sp}(T)=\gamma_{sp}(T_1)+\gamma_{sp}(T_2),$ we have a contradiction with Observation \ref{drzewa}.

$b)$ Suppose $x_3\notin D.$

Observe that $\Omega(T)\cap N_T(x_2)\subseteq D.$ If $N_T(x_3)\cap
D\neq \{x_2\},$ then $D-\{x_2\}$ is a super dominating set of $T.$
 Therefore $x_2$ is
the only neighbour of $x_3$ belonging to $D$ and $x_4\notin D.$

Let $T'=T-\{x_3x_4\}=T_1\cup T_2,$ where $x_3\in T_1,$ $x_4\in T_2.$ Since $x_3, x_4\notin D,$ $\gamma_{sp}(T)=\gamma_{sp}(T_1)+\gamma_{sp}(T_2).$

If $T_1$ is a star, then we have the same contradiction as above.
Thus $T_1$ is not a star. Since $N_T(x_3)\cap D=\{x_2\},$ there
exists $v\in \Omega(T)$ such that $d_T(v,x_3)=2$ and $v\in D.$
Similarly to the case above,
$\gamma_{sp}(T_1)>\frac{|V(T_1)|}{2},$ a contradiction.

$Case$ $2.$ Suppose $x_2\notin D.$
 If $d_T(x_2)>2,$ then $N_T(x_2)\cap \Omega(T)\subseteq D.$ Let $T'=T-{x_2x_3}=T_1\cup T_2,$ where $x_2\in T_1, x_3\in T_2.$  Then $T_1$ is a star with three or more vertices and $\gamma_{sp}(T_1)=|V(T_1)|-1>\frac{|V(T_1)|}{2}.$

Therefore $d_T(x_2)=2.$ Let $T^{''}=T-\{x_1, x_2\}.$ We have $|V(T^{''})|=n-2=2(m-1)$ and $\gamma_{sp}(T^{''})=\gamma_{sp}(T)-1=m-1.$ By the induction hypothesis, $T^{''}\in \mathcal{R}.$ Then $T=T^{''}\cup \{x_1, x_2, x_1x_2, x_2x_3\}\in \mathcal{R}$ with $sta(x_2)=sta(x_3)$ and $sta(x_1)\neq sta(x_2).$

$Case$ $3.$ Suppose $x_1\notin D.$ Then $x_2\in D.$ If
$d_T(x_2)>2,$ then there exists $z\in N_T(x_2)\cap \Omega(T)$ such that $z\neq x_1$ and $z\in D.$ Let $T^{'''}=T-zx_2=T_1\cup T_2, T_1=\{z\}.$ Then we have $\gamma_{sp}(T)=\gamma_{sp}(T_1)+\gamma_{sp}(T_2)\geq 1+\frac{n(T_2)}{2}=\frac{n+1}{2}>\frac{n}{2},$ a contradiction.

 Therefore
$d_T(x_2)=2$ and the proof $T\in \mathcal{R}$ is the same as in
the above case.
\end{proof}

Now we are in position to characterize all trees for which the upper bound is attained.
Let  $\mathcal{S}$ be the family of trees $T$ that can be obtained from a sequence $T_1,\ldots,T_j$ ($j\geq 1$) of trees such that $T_1=P_3=(a_1,b_0,b_1)$ and if $j\geq 2,$ define the tree $T_j$ in the following three ways:\\

$1.$ $V(T_j)=V(T_{j-1})\cup \{a_l, b_k\}$ and $E(T_j)=E(T_{j-1})\cup \{a_lb_k\}\cup \{b_ib_k\}$ for some $b_i\in V(T_{j-1});$ or \\

$2.$ $V(T_j)=V(T_{j-1})\cup \{a_l\}$ and $E(T_j)=E(T_{j-1})\cup \{a_lb_i\}$ for some non-support vertex  $b_i\in V(T_{j-1});$ or \\

$3.$  $V(T_j)=V(T_{j-1})\cup \{ b_k\}$ and $E(T_j)=E(T_{j-1})\cup \{b_kb_i\}$ for some support vertex  $b_i\in V(T_{j-1}).$ \\

If $T\in \mathcal{S},$  denote by $A=\{a_1,\ldots, a_l\}$ and $B=\{b_1,\ldots,b_k\},$ where $l+k=n.$

Similarly like in the family $\mathcal{R},$ we will say that a vertex $v\in A$ has a status $a$ ($sta(v)=a$) and a vertex $v\in B$ has a status $b$ ($sta(v)=b$).

Observe that for all $1\leq i, |N_T(a_i)\cap B|=1$ and $|N_T(b_i)\cap A|\leq 1.$

By definition of the family $\mathcal{S},$ it is clear that the
set of vertices $B$ is a minimum super dominating set of $T$ and
$|A|=|Sp|,$ where $Sp$ is the set of supports in $T$. Therefore
$|B|=|V-A|=|V-Sp|=n-s.$ Moreover, if $v$ is a non-end vertex  of
$T,$ then $sta(v)=b.$

\begin{thm} If $T$ is a tree with at least three vertices and $s$ is the number of supports in $T,$ then $\gamma_{sp}(T)=n-s$ if and only if $T\in \mathcal{S}.$
\end{thm}
\begin{proof} By definition of the family $\mathcal{S},$ if $T\in \mathcal{S},$ then $\gamma_{sp}(T)=n-s.$

Assume $T$ is a tree with $|V(T)|=n$ and $\gamma_{sp}(T)=n-s.$ We
show that $T\in \mathcal{S}$ by induction on $n.$

If $n=3,$ then $T=P_3=(a_1, b_0, b_1)$ and $T\in \mathcal{S}.$ If $n=4,$ then $T=P_4$ or $T=K_{1,3};$ in both cases $T\in \mathcal{S}.$

Assume $n\geq 4$ and the Theorem holds for every $m<n.$

Let $P=(v_1,\ldots,v_p)$ be a longest path in $T$ such that $d_T(v_2)$ is as big as possible. Let $D$ be a minimum super dominating set of $T.$ We consider the following cases.\\

$Case$ $1.$ Assume $d_T(v_2)>2.$ Let $T'=T-\{v_1\}.$ Then
$|V(T')|=n-1$ and $s(T')=s(T)=s.$ If $v_1\in D,$ then $D-\{v_1\}$
is a minimum  super dominating set of $T'$ and
$\gamma_{sp}(T')=\gamma_{sp}(T)-1=(n-1)-s.$ By the induction
hypothesis, $T'\in \mathcal{S}.$ Since $v_2$ is a non-end vertex
of $T',$ we have $sta(v_2)=b.$ Putting $sta(v_1)=b$ we have that
$T\in \mathcal{S}.$

Now, if $v_1\notin D,$ then $v_2\in D$ and since $d_T(v_2)>2,$ there exists a vertex $x\in N_T(x_2)\cap \Omega(T)$ such that $x\in D.$ Define $T'=T-\{x\}$ and similarly like in the above case we obtain $T\in \mathcal{S}.$\\

$Case$ $2.$ Suppose $d_T(v_2)=2.$ Here we have three cases.\\

$1.$ The set $\{v_1, v_2\}\subseteq D.$ Since $D$ is minimum, $v_3\notin D.$ If $v_3$ is a support vertex, then all of its end neighbours belong to $D$ and $D-\{v_2\}$ is a super dominating set of $T,$ a contradiction. If $v_3$ is not a support vertex, $d_T(v_3,\Omega(T))=2$ and also  $D-\{v_2\}$ is a super dominating set of $T,$ a contradiction.\\

$2.$ The vertex $v_1\in D$ and $v_2\notin D.$ Let $T'=T-\{v_1,v_2\}.$ We have $\gamma_{sp}(T')=\gamma_{sp}(T)-1,$ $|V(T')|=|V(T)|-2, s(T')=s-1.$ Then $\gamma_{sp}(T')=|V(T')|-s(T').$ From the induction hypothesis, $T'\in \mathcal{S}.$ We can obtain $T$ from $T',$ putting $sta(v_1)=a, sta(v_2)=b$ and thus $T\in \mathcal{S}.$\\

$3.$ The vertex $v_1\notin D$ and $v_2\in D.$ By definition of $D,$ $v_3\in D.$ Similarly to the previous case we can consider $T'=T-\{v_1, v_2\}$ and prove that $T\in \mathcal{S}.$
\end{proof}