text:
\section{Introduction}

The standard thin accretion disc (\citealp{Shakura73})
has been extensively applied to a variety
of astrophysical systems, such as active galactic nuclei,
X-ray binaries, and cataclysmic variable stars.
The theory of the standard thin disc, which is based on the $\alpha$
stress prescription, predicts that the disc will
suffer both thermal and viscous instability when the pressure
is dominated by radiation pressure.
However, observations of X-ray binaries show little evidence
for such instabilities except for GRS 1915+105 (hereafter GRS 1915).
Since the stability is related to the description of stress,
the standard $\alpha$ stress prescription
was modified in some previous works to assure the disc to be stable
no matter whether gas or radiation pressure dominates.
For instance, \citet{Lightman74} introduced a so-called $\beta$
stress prescription, i.e., the stress is proportional to the gas pressure
rather than the sum of the gas and radiation pressure.
With the $\beta$ stress, the disc will remain stable
even for the case that radiation pressure dominates over gas pressure.
Furthermore, some other stress prescriptions were introduced by assuming
that the stress is relevant to both the total pressure and the gas pressure
(e.g., \citealp{Kato08}, p.~108).
As pointed out by \citet{Gierli04}, however,
the $\beta$ stress would be in conflict with observations as follows.
Observations show that the colour
temperature correction remains stable for a wide range of luminosity,
whereas discs with $\beta$ stress should have a colour temperature correction
that changes from $\sim 1.8$ to $\sim 2.7$ as the accretion rate increases.
Conversely, the standard $\alpha$ prescription may predict that
the correction should remain stable at $\sim 1.8$ for varying accretion rates.
Moreover, some recent simulations (e.g., \citealp{Hirose09})
showed that the stress is roughly proportional to the sum of gas and
radiation pressure. Thus, the standard $\alpha$ prescription seems to be
more acceptable than other prescriptions such as the $\beta$ stress.

The slim disc model was introduced by \citet{Abramowicz88},
in which it was predicted that
the thermal instability of the radiation-pressure-dominated thin
disc will trigger a limit cycle behaviour between the thin disc
and the slim disc. Such a limit cycle was confirmed by time-dependent
numerical calculations
(e.g., Szuszkiewicz \& Miller 2001; \citealp{Li07}).
On the other hand, \cite{Belloni97}
applied the limit cycle theory to interpret the observational
quasi-periodic variability of GRS 1915.
Recently, such limit cycle has also been applied to young radio
galaxies (e.g., \citealp{Czerny09,Wu09}).
Comparing the numerical results with the observational ones,
however, we can find the following inconsistency
between the theory and the observation. The numerical results of
Li et al.~(2007) showed that the duration of the outburst phase
$t_{\rm high}$ is less than 5 percent of the duration of the quiescent phase
$t_{\rm low}$, and the luminosity of the outburst phase $L_{\rm high}$
is around two orders of magnitude larger than that of the quiescent phase
$L_{\rm low}$; whereas observations of GRS 1915 showed that,
$t_{\rm high}$ is comparable to $t_{\rm low}$ and
$L_{\rm high}$ is only 3-20 times larger than $L_{\rm low}$
(e.g., \citealp{Belloni97, Wu10}).
Some efforts have been made to improve the theory in order to explain
observations either by some artifical viscosity prescription
(e.g., \citealp{Nayakshin00}) or by additional assumption of the
energy exchange between the disc and corona (Janiuk et al.~2002).

In the present paper, we will investigate this issue through
another way by taking the stress evolution process into account.
Our paper is organized as follows. In Section 2, we describe the
basic assumption for the stress evolution process. The time-dependent
equations and linear stability analysis are presented in Sections 3 and 4,
respectively. Conclusions and discussion are made in Section 5.

\section[]{Model}

The stress in MHD accretion discs can be expressed as
(e.g., Balbus \& Papaloizou 1999)
\begin{equation}
\tau_{r\varphi} =
<{\rho}{\delta\upsilon_r}{\delta\upsilon_{\varphi}}
- \frac{{{\delta}B_r}{{\delta}B_{\varphi}}}{4\pi}>_\varphi \ ,
\end{equation}
where ${\delta\upsilon_r}$(${{\delta}B_r}$) and
${\delta\upsilon_{\varphi}}$(${\delta}B_{\varphi}$)
are the fluctuating components of the velocity (the magnetic field)
in $r$ and $\varphi$ direction, respectively.
The angled bracket notation denotes the azimuthal average.

The above formula is a general form and is complicated for calculation.
In the study of accretion discs,
a more practical formula, i.e., the $\alpha$ stress prescription (\citealp{Shakura73}),
is widely adopted:
\begin{equation}
\tau_{r\varphi}=-{\alpha}p_{\rm tot} \ ,
\end{equation}
where $p_{\rm tot}$ is the total pressure, and $\alpha$ is regarded
as a constant less than unity.
With the above simple form, the $\alpha$ stress prescription has been
extensively applied to the research of accretion discs,
even in the disc stability analysis. But in the real system,
when the pressure varies, it will take a certain period
for the stress to adjust.
As mentioned in \S~1, the stability is sensitive to the description
of stress, thus the stress evolution process
may have essential effects on the stability.
In this paper, we will modify the standard $\alpha$ stress prescription
by taking the stress delay into account.

Generally, the growth rate of the stress should be related to
$\tau_{r\varphi}^{\rm exp}$ and $\tau_{r\varphi}$,
where $\tau_{r\varphi}^{\rm exp}$ is the expected
stress and $\tau_{r\varphi}$ is the current stress.
Thus we have the following general function for the evolution of stress:
\begin{equation}
\frac{{\partial}\tau_{r\varphi}}{{\partial}t}+\upsilon_r\frac{{\partial}
\tau_{r\varphi}}{{\partial}r}=f(\tau_{r\varphi}^{\rm exp}, \tau_{r\varphi})
\ ,
\end{equation}
where $\upsilon_r$ is the radial velocity.
In a quasi-stationary state, the function $f$ should match
the following conditions:
\begin{equation}
\left \{
\begin{array}{cc}
f(\tau_{r\varphi}^{\rm exp}, \tau_{r\varphi}) < 0;
& \tau_{r\varphi} > \tau_{r\varphi}^{\rm exp} \\
f(\tau_{r\varphi}^{\rm exp}, \tau_{r\varphi}) = 0;
& \ \ \tau_{r\varphi} = \tau_{r\varphi}^{\rm exp} \ . \\
f(\tau_{r\varphi}^{\rm exp}, \tau_{r\varphi}) > 0 ;
& \tau_{r\varphi} < \tau_{r\varphi}^{\rm exp}
\end{array} \right.
\end{equation}
In the present paper, we adopt a simple form as follows:
\begin{equation}
f(\tau_{r\varphi}^{\rm exp}, \tau_{r\varphi}) =
\frac{\tau_{r\varphi}^{\rm exp}-\tau_{r\varphi}}{t_{\rm {ps}}} \ ,
\end{equation}
where $t_{\rm ps}$ is the timescale of stress delay for varying pressure.
We further adopt $\tau_{r\varphi}^{\rm exp} = -\alpha{p_{\rm tot}}$
according to the $\alpha$ prescription, where $\alpha$ is a constant.
Combining equations (3) and (5),
the evolution equation of the stress tensor is
\begin{equation}
\frac{{\partial}\tau_{r\varphi}}{{\partial}t}
+ \upsilon_r\frac{{\partial}\tau_{r\varphi}}{{\partial}r}
= \frac{-{\alpha}p_{\rm {tot}}
- \tau_{r\varphi}}{t_{\rm {ps}}} \ .
\end{equation}

\section{Equations}

We first present the time-dependent equations of accretion discs
for stability analysis.
We consider an axisymmetric, Keplerian rotating disc
($\Omega=\Omega_{\rm K}$) under cylindrical coordinates
($r$, $\varphi$, $z$).
The basic vertically integrated equations are described as follows
(e.g., Kato et al.~2008):
\begin{equation}\label{MassE}
\frac{{\partial}\Sigma}{{\partial}t}
+ \frac{1}{r}\frac{{\partial}}{{\partial}r}(r{\Sigma}\upsilon_r) = 0 \ ,
\end{equation}
\begin{equation}\label{AzimuthalE}
{\Sigma}\upsilon_r\frac{1}{r}\frac{\partial}{{\partial}r}(r^2\Omega)
= \frac{1}{r^2}\frac{\partial}{{\partial}r}(r^2T_{r\varphi}) \ ,
\end{equation}
\begin{equation}\label{VerticalE}
\Omega_{\rm K}^2H^2 = \frac{\Pi}{\Sigma} \ ,
\end{equation}
\begin{equation}\label{EnergyE}
\frac{{\partial}E}{{\partial}t}
- (E+\Pi)\frac{{\partial}\ln\Sigma}{{\partial}t}
+ \Pi\frac{{\partial}{\ln}H}{{\partial}t}
= Q_{\rm {vis}}^+ - Q_{\rm {adv}}^- - Q_{\rm {rad}}^- \ ,
\end{equation}
\begin{equation}\label{StateE}
\Pi = \Pi_{\rm {gas}}+\Pi_{\rm {rad}} = \frac{k_{\rm B}}
{\mu m_{\rm H}}{\Sigma}T + \frac{2a}{3}T^4H \ ,
\end{equation}
where $\Sigma = 2 \rho H$ is the surface mass density,
$\Pi= 2 p_{\rm tot} H$, $T_{r\varphi} = 2\tau_{r\varphi} H$, and $E$
are the vertically integrated pressure, stress, and internal energy,
respectively. $T$ is the temperature, $H$ is the half-thickness
of the disc, $k_{\rm B}$ is Boltzmann constant, $\mu$ is the mean
molecular weight, $m_{\rm H}$ is the hydrogen mass, and $a$ is
radiation constant.
The quantities
$Q_{\rm vis}^+$, $Q_{\rm adv}^-$, and $Q_{\rm rad}^-$ represent the viscous
heating, advective cooling, and radiative cooling rates, respectively.
The expressions of $E$, $Q_{\rm vis}^+$, $Q_{\rm adv}^-$,
and $Q_{\rm rad}^-$ are the following:
\begin{displaymath}
E= \left [ 3(1-\beta)+\frac{\beta}{\gamma-1} \right ] \Pi \ ,
\end{displaymath}
\begin{displaymath}
Q_{\rm {vis}}^+ = T_{r\varphi}(r\frac{d\Omega}{dr}) \ ,
\end{displaymath}
\begin{displaymath}
Q_{\rm {adv}}^- = \upsilon_r\left [
-(E+\Pi)\frac{\partial{\ln\Sigma}}{{\partial}r}
+ \frac{\partial{E}}{{\partial}r}
+ \Pi\frac{\partial{\ln{H}}}{{\partial}r}\right ] \ ,
\end{displaymath}
\begin{displaymath}
Q_{\rm {rad}}^- = \frac{16acT^4}{3\bar{\kappa}{\Sigma}} \ ,
\end{displaymath}
where $\beta$ is defined as the ratio of the gas pressure
to the total pressure, i.e., $\beta \equiv \Pi_{\rm gas}/\Pi$,
$\gamma$ is the ratio of specific heating, and $\bar{\kappa}$ is
the opacity.

Since the timescale for establishment of vertical equilibrium is
the dynamic timescale, which is much shorter than the thermal timescale,
we modify the stress evolution process (Eq.~[6]) as
\begin{equation}\label{StressE}
\frac{{\partial}T_{r\varphi}}{{\partial}t}
+ \upsilon_r\frac{{\partial}T_{r\varphi}}{{\partial}r}
= \frac{-\alpha\Pi-T_{r\varphi}}{t_{\rm {ps}}} \ .
\end{equation}

\section{Stability analysis}

Now we have the set of six equations, i.e., Eqs.~(7)-(12) for linear
stability analysis, and the corresponding six physical quantities
are $\Sigma$, $\Pi$, $\upsilon_r$, $H$, $T$, and $T_{r\varphi}$.
With subscript 0 denoting the quasi-stationary solution, and subscript 1
denoting perturbation quantities, we define the following dimensionless
variables:
\begin{displaymath}
\sigma = \frac{\Sigma_1}{\Sigma_0},\
\xi = \frac{\Pi_1}{\Pi_0},\
u = \frac{\upsilon_{r, 1}}{r\Omega},\
h = \frac{H_1}{H_0},\
m = \frac{T_{1}}{T_0},\
l = \frac{T_{r\varphi, 1}}{-\alpha\Pi_0}.
\end{displaymath}
With the assumption that all the perturbation quantities are proportional
to $\exp(\omega t - ikr)$,
the perturbed equations corresponding
to Eqs.~(7)-(12) can be written as follows:
\begin{equation}\label{MassP}
\frac{\omega}{\Omega}\sigma-ikru=0 \ ,
\end{equation}
\begin{equation}\label{AzimuthalP}
\frac{\kappa^2}{2\Omega^2}u
= i\alpha kr(\frac{1}{r{\Omega}})^2\frac{\Pi_0}{\Sigma_0}l \ ,
\end{equation}
\begin{equation}\label{VerticalP}
2h=\xi-\sigma \ ,
\end{equation}
\begin{eqnarray}\label{EnergyP}
\frac{\omega}{\Omega}[3(1-\beta)+\frac{\beta}{\gamma-1}]\xi
+ \frac{\omega}{\Omega}{\beta}\frac{4-3\gamma}{\gamma-1}
\frac{1-\beta}{4-3\beta}(-3\xi+4\sigma-h)
\nonumber \\
- \frac{\omega}{\Omega}(4-3\beta+\frac{\beta}{\gamma-1})\sigma
+ \frac{\omega}{\Omega}h
= -\frac{3}{2}{\alpha}l - 6\alpha m + \frac{3}{2}\alpha \sigma \ ,
\end{eqnarray}
\begin{equation}\label{StateP}
\xi = \beta(\sigma+m)+(1-\beta)(4m+h) \ ,
\end{equation}
\begin{equation}\label{StressPS}
l = \frac{1}{{\omega}t_{\rm {ps}}+1}\xi
-({\Omega}t_{\rm {ps}})\frac{{\partial}{\ln}\Pi_0}
{{\partial}{\ln}r}\frac{1}{{\omega}t_{\rm {ps}}+1}u \ ,
\end{equation}
where ${\kappa}$ is the epicyclic frequency defined as
${\kappa^2} \equiv 2\Omega(2\Omega+r{d\Omega}/{dr})$.

Before numerical calculations, we would point out that,
for long-wavelength unstable modes,
the second term on the right hand side of Eq.~(18) is significantly
less than the first term, thus the equation can be reduced to
\begin{equation}
l = \frac{1}{{\omega}t_{\rm {ps}}+1}\xi \ .
\label{StressP}
\end{equation}
The reason for the above simplification is as follows.
Substituting equation (\ref{AzimuthalP}) into (\ref{StressPS}),
we will have
\begin{equation}\label{Reduction}
\left [\frac{\kappa^2}{2\Omega^2}\frac{({r{\Omega}})^2{\Sigma_0}}
{i\alpha kr{\Pi_0}}+({\Omega}t_{\rm {ps}})\frac{{\partial}{\ln}\Pi_0}
{{\partial}{\ln}r}\frac{1}{{\omega}t_{\rm {ps}}+1} \right ]
u
=  \frac{1}{{\omega}t_{\rm {ps}}+1}\xi,
\end{equation}
For standard thin discs,
${{\partial}{\ln}\Pi_0}/{{\partial}{\ln}r}$
and ${\kappa^2}/{(2\Omega^2)}$ are around unity.
Thus, for long-wavelength unstable modes, on the left hand side of Eq.~(20),
the ratio of the second term to the first term is
\begin{eqnarray}
\left|
\left[ ({\Omega}t_{\rm {ps}})\frac{{\partial}{\ln}\Pi_0}
{{\partial}{\ln}r}\frac{1}{{\omega}t_{\rm {ps}}+1} \right]/
\left [ \frac{\kappa^2}{2\Omega^2}\frac{({r{\Omega}})^2{\Sigma_0}}
{i\alpha kr{\Pi_0}} \right ] \right| \nonumber \\
\approx
\left| \frac{{\partial}{\ln}\Pi_0}{{\partial}{\ln}r}\right|
\left| \frac{1}{{\omega}t_{\rm {ps}}+1}\frac{1}{\kappa^2
/({2\Omega^2})}\right|
\left| 2\frac{t_{\rm {ps}}}{t_{\rm {th}}}\frac{H}{r} \right|
\left| kH \right|
{\ll}1 \ .
\end{eqnarray}
Thus, equation (\ref{Reduction}) can be simplified as
\begin{displaymath}
\left [ \frac{\kappa^2}{2\Omega^2}\frac{({r{\Omega}})^2{\Sigma_0}}
{i\alpha kr{\Pi_0}} \right ]
u
= \frac{1}{{\omega}t_{\rm {ps}}+1}\xi \ ,
\end{displaymath}
which is equivalent to Eq.~(19).

For non-trivial solutions of $\sigma$, $\xi$, $u$, $h$, $m$, and $l$,
we obtain the following dispersion relation from equations
(\ref{MassP})-(\ref{StateP}) and (\ref{StressP}),
\begin{eqnarray}\label{DissipationE}
c_1(\frac{\omega}{\Omega})^3
+ c_2(\frac{\omega}{\Omega})^2
+ c_3(\frac{\omega}{\Omega})
+ c_4 = 0 \ ,
\end{eqnarray}
where
\begin{displaymath}
c_1 = ({\Omega}t_{\rm {ps}})\left ( \frac{7-6\beta}{2}+\frac{\beta}{\gamma-1}
- \frac{7}{2}\beta\frac{4-3\gamma}{\gamma-1}\frac{1-\beta}{4-3\beta} \right )
\ ,
\end{displaymath}
\begin{eqnarray*}
c_2= \left ( \frac{7-6\beta}{2}+\frac{\beta}{\gamma-1}
- \frac{7}{2}\beta\frac{4-3\gamma}{\gamma-1}\frac{1-\beta}{4-3\beta} \right )\\
+ \frac{3\alpha(1+\beta)}{(4-3\beta)}({\Omega}t_{\rm {ps}})\ ,
\end{eqnarray*}
\begin{eqnarray*}
c_3 =2\alpha (kH)^2\left( \frac{\Omega}{{\kappa}} \right)^2\left (
\frac{9-6\beta}{2}+\frac{\beta}{\gamma-1}-
\right.\\
\left.\frac{9}{2}\beta\frac{4-3\gamma}{\gamma-1}\frac{1-\beta}{4-3\beta}  \right)
- \frac{3\alpha}{2}+\frac{3\alpha_(1+\beta)}{(4-3\beta)}
\ ,
\end{eqnarray*}
\begin{displaymath}
c_4 = 3\alpha^2 (kH)^2\left( \frac{\Omega}{{\kappa}} \right)^2
\frac{2+3\beta}{4-3\beta}
\ .
\end{displaymath}
When $t_{\rm ps} \to 0$, the above dispersion relation is reduced to
the form in the $\alpha$ stress scenario (e.g., Kato et al.~2008, p.~180).

\begin{figure}
\begin{center}
{\includegraphics[width=3.00in,height=2.70in]{f1.eps}}
\caption{Variation of the growth rate with the radial wavelength
of perturbation.
The thin and thick lines correspond to the standard $\alpha$ stress
prescription and the stress evolution process with $t_{\rm ps} = t_{\rm th}$,
respectively.
The solid, dashed, and dotted lines correspond to $\beta=0.01$, $0.1$,
and $0.4$, respectively.
\label{fig1}}
\end{center}
\end{figure}

\begin{figure}
\begin{center}
{\includegraphics[width=3.00in,height=2.70in]{f2.eps}}
\caption{Variation of the growth rate with the timescale of stress delay
in the long-wavelength limit ($kH \to 0$).
\label{fig2}}
\end{center}
\end{figure}

We numerically solve the dispersion equation (\ref{DissipationE})
to obtain the growth rates of unstable modes.
In our calculation, we fix $\alpha = 0.1$ and $\gamma = 5/3$.
There exist three modes for each given radial wavelength of perturbations.
Among the three modes, one is always stable corresponding to a negative
$\omega$, and is therefore neglected in the following analysis.
We study the dispersion relation for two cases: (i) with stress
evolution process ($t_{\rm ps} = t_{\rm th}$); (ii) with the standard
$\alpha$ stress ($t_{\rm ps}=0$).
The growth rates of the remaining two modes are shown in Figure~1.
In this figure, the thick and thin lines correspond to the stress
evolution process and the $\alpha$ stress prescription, respectively.
The solid, dashed, and dotted lines represent the solutions with
$\beta=0.01$, $0.1$, and $0.4$, respectively, where $\beta=0.4$ is
the well known critical value for
the thermal stability with the $\alpha$ prescription.
It is shown that the growth rate of thermally unstable modes (the upper branch)
in the stress evolution scenario is significantly smaller than that in
the $\alpha$ stress case. On the contrary,
the growth rates of the viscously unstable modes
(the lower branch) of the above two cases are similar.
Here, the thermal timescale $t_{\rm th}$ is calculated by
\begin{equation}
t_{\rm th} = \frac{E}{Q_{\rm vis}^+}
= \frac{2[3(1-\beta)(\gamma-1)-\beta]}{3(\gamma-1)}\frac{1}{\alpha\Omega} \ .
\end{equation}
Figure~1 indicates that, even though the stress evolution process does not
change the stability criterion, it may have
significant influence on the growth rate of the thermally unstable mode.

For further investigation, we calculate the solutions under
the long-wavelength limit, i.e., $kH \to 0$. In such case,
the non-zero roots of the dispersion equation (\ref{DissipationE}) are,
\begin{equation}\label{GrowthRate}
\omega = {\frac{1}{2A t_{\rm th}}}
\left [-(\frac{A t_{\rm th}}{t_{\rm ps}}+1)
{\pm} \sqrt{(\frac{At_{\rm th}}{t_{\rm ps}}+1)^2
+ 2\frac{A t_{\rm th}}{t_{\rm ps}}\frac{2-5\beta}{1+\beta}} \right ] \ ,
\end{equation}
where
\begin{displaymath}
A = \frac{\left ( \frac{7-6\beta}{2}+\frac{\beta}{\gamma-1}
- \frac{7}{2}\beta\frac{4-3\gamma}{\gamma-1}\frac{1-\beta}{4-3\beta} \right )}
{[3(1-\beta)+\frac{\beta}{\gamma-1}]}
\frac{4-3\beta}{2(1+\beta)} \ .
\end{displaymath}
The above equation shows the relationship between the growth rate
$\omega$ and the timescale of stress delay $t_{\rm ps}$.

For the case of $t_{\rm ps}/t_{\rm th} \gg 1$, Eq.~(\ref{GrowthRate})
is simplified as
\begin{equation}
\omega = {\frac{1}{2t_{\rm ps}}}
\frac{2-5\beta}{1+\beta} \sim {\frac{1}{t_{\rm ps}}} \ ,
\end{equation}
which means that if $t_{\rm ps}$ is sufficiently larger than $t_{\rm th}$,
the growing timescale of thermal instability is around $t_{\rm ps}$.
On the contrary, for the case of $t_{\rm ps}/t_{\rm th}{\ll}1$,
Eq.~(\ref{GrowthRate}) is reduced to
\begin{equation}
\omega = {\frac{1}{2At_{\rm {th}}}} \frac{2-5\beta}{1+\beta}
\sim {\frac{1}{t_{\rm th}}} \ .
\end{equation}
The above growth rate is exactly the same as that in the $\alpha$ stress case,
in which the growing timescale of thermal instability is around
$t_{\rm th}$.

Figure~2 shows the variation of the growth rate of
the thermally unstable mode in the long-wavelength limit with
$t_{\rm ps}$ for $\beta = 0.01$, 0.1, and 0.3.
The figure clearly illustrates that, for the case of
$t_{\rm ps} \ga t_{\rm th}$,
the stress evolution process will have essential effects on the thermal
instability since the growth rate is significantly decreased.

\section{Conclusions and discussion}

In this paper, the stress evolution process, which was ignored
in previous works, is taken into account in the linear stability analysis
of standard thin accretion discs.
We show the variation of the growth rate of thermally unstable modes
with the wavelength of perturbations for two cases: with and without stress
delay. We find that the growth rate with stress delay can be
apparently lower than that without stress delay.
We also make some analytical approximation for the growth rate in
long-wavelength limit, and present the relationship between the
specific growth rate and the timescale of stress delay.
In conclusion, the stress evolution process may have essential
influence on the thermal instability by significantly decreasing
the growth rate, in particular for the case in which
the timescale of stress delay is comparable to or even larger than
the thermal timescale.

In a real system the timescale of stress delay remains
unclear. We would argue that this timescale is possibly
comparable to the thermal timescale as follows.
In simulations with an initially weak toroidal or poloidal magnetic field,
the magnetic energy may first experience
an exponential growth during the first few orbits due
to the linear instability, and then will be followed by the nonlinear evolution.
Finally, a saturated quasi-steady state phase may form.
Totally, it will take $15 \sim 20$ orbits for the system to enter
a fully turbulent state (e.g., \citealp{Hawley96,Fromang07}).
This implies that the evolution timescale of the
magnetoturbulence is around $15 \sim 20$ times of the dynamical timescale
$t_{\rm dyn}$, which can be regarded as a duration around
the thermal timescale $t_{\rm th}$
(e.g., $t_{\rm th} \sim t_{\rm dyn}/\alpha$ and $\alpha = 0.1$).

As mentioned in the first section, when the theory of the limit cycle behaviour
between the standard thin disc and the slim disc is applied to the
observational quasi-periodic variability of GRS 1915, there exists
some conflict between the theory and the observation
on the duration ratio $t_{\rm high}/t_{\rm low}$ and
the luminosity ratio $L_{\rm high}/L_{\rm low}$.
In our opinion, the stress evolution process may improve the theory
to explain observations
due to the following reasons. When the flow suffers thermal instability,
the growth rate of $\dot M$ can be significantly decreased by
the stress delay.
We therefore can expect a relatively lower $\dot M$
for the outburst phase, thus a corresponding lower $L_{\rm high}/L_{\rm low}$.
As a consequence, for a certain fixed mass accretion rate supply at outer
boundary, if $\dot M$ at outburst phase drops, the duration of this phase
will therefore become longer, thus a moderate value of
$t_{\rm high}/t_{\rm low}$ should appear.
In other words, the influence of the stress delay
on the growth rate may result in moderate ratios of duration and luminosity
of the outburst phase to the quiescent one,
which will be more likely to explain the observational results than
previous calculations. The further investigation
on this issue may require time-dependent numerical calculations including
the stress evolution process.

Moreover, we would point out the possible application of the stress evolution
process to another geometrically thin disc,
namely the Shapiro-Lightman-Eardley disc (the SLE disc, Shapiro et al.~1976),
which was originally introduced as the inner region of a thin
disc to provide hard X-ray emission. However, it is known that the SLE disc
is viscously stable but suffers thermal instability.
If the stress evolution process is taken into
consideration in the SLE model, and the timescale of stress delay is
comparable to the thermal timescale, we can therefore expect the existence
of a quasi-stable SLE disc since the decreased growth rate of
thermal instability may not completely destroy the disc and the viscous
stability may help to suppress the thermal instability
in the viscous timescale.

\section*{Acknowledgments}

We thank Shoji Kato and Sheng-Ming Zheng for beneficial discussion.
This work was supported by the National Basic Research Program of China
under grant 2009CB824800, and the National Natural Science Foundation
of China under grants 10833002 and 11073015.


text:
\section{Introduction}
\label{sec:introduction}

Turbulence above sloping ocean boundaries has been the subject of several studies since the work of \citet{munk_abyssal_1966}, who hypothesised that enhanced turbulence close to the solid boundary plays an important role in maintaining the ocean stratified in density from the surface to the bottom.
The dynamics of the benthic boundary layer above a slope has been the subject of several studies \citep[e.g.][]{garrett_role_1990,thorpe_variability_1990,garrett_marginal_1991,van_haren_measurements_1994,alford_observations_2000,slinn_modeling_2003,gayen_boundary_2011,gayen_negative_2011}.
Turbulence close to sloping boundaries is enhanced by the breaking of internal waves \citep[see the review of][]{lamb_internal_2014} which, from a broader point of view, dissipates tidal energy \citep{munk_abyssal_1998} and, to a lesser extent, eddy kinetic energy \citep{zhai_significant_2010}.

Here, we consider turbulence above a deep-ocean sloping boundary from the particular perspective of scalar (temperature) statistics.
The interest in the scalar statistics embedded in a turbulent flow is manifold.
The scalar field is highly intermittent, and its statistical moments bear the marks of this intermittency.
In relation to this, the study of the statistics of a scalar is the basis for understanding the dynamics driving its mixing, and for reliably estimating and modelling mixing itself.
The statistics of scalar tracers in a turbulent flow is a widely studied topic in the laboratory as well as numerically and analytically \citep[see the reviews of][]{shraiman_scalar_2000,warhaft_passive_2000}; most often the focus is on grid-generated turbulence.

On the other hand, intermittency of scalar fields in strongly stratified turbulent flows, where the buoyancy force is expected to be dominant, has received limited attention \citep[e.g.][using numerical simulations]{brethouwer_passive_2008}.
While we are unaware of laboratory studies on this topic, \citet{seuront_multifractal_1999} discussed temperature fluctuations in the surface ocean boundary layer in terms of multifractal intermittency.
\citet{frola_experimental_2014} recently reported observations of the probability density function (pdf) of temperature increments in the atmospheric boundary layer, discussing it in terms of a bifractal intermittency model.

We must stress that observations in the ocean and atmosphere often include some temperature statistics, but these are virtually always limited to second order, i.e.~to spectra \citep[see, for instance,][and references therein]{riley_stratified_2008}.
The oceanographic literature usually considers observations of temperature at the dissipation scale with the aim of estimating turbulent diffusivity and dissipation rate by using the \citet{osborn_oceanic_1972} model, and only second-order statistics are presented.
An exception is \citet{thorpe_skewness_1991}, where a possible link between the skewness of temperature increments and convection in benthic boundary layers was discussed.
Possibly more often, ocean temperature (or density) observations are collected at lower resolution, comparable to the Ozmidov scale \citep[e.g.][]{alford_observations_2000}.

This work presents a detailed study of the temperature statistics above a deep-ocean sloping boundary, focusing on the intermittency properties of the scalar.
The results are discussed comparing them to laboratory experiments with both passive (shear turbulence) and active (convective turbulence) scalars.
By necessity, our observations provide an average view of fluid regions having different Reynolds number, different mean shear and different mean stratification.
While this is to some extent a limitation, it should be kept in mind that geophysical turbulence is intrinsically sporadic i.e.~intermittent also at scales larger than those of the inertial range \citep{rorai_turbulence_2014}.
This fact must be taken into account when studying the scalar statistics in a stratified flow, and is naturally included in real-ocean observations.
Furthermore, laboratory and numerical experiments mainly focus on grid-generated turbulence or on the evolution of a particular kind of instability, usually Kelvin--Helmholtz.
A unique forcing scale (grid length, most unstable wavenumber) can usually be identified.
In the present observations, on the other hand, there is no obvious separation between breaking waves (the main forcing) and turbulence.
Instabilities such as those produced in the laboratory are seldom recognised in the ocean, even if new observations may change this \citep[][and references therein]{van_haren_deep-ocean_2010,haren_stratified_2013,smyth_ocean_2012}.
For this reason, the relevance of the conclusions drawn from the study of Kelvin--Helmholtz instabilities for the real ocean is still not clear \citep{thorpe_kelvinhelmholtz_2012}.
In this context, the main aim of this work is to help bridge the gap between laboratory and geophysical flows, by analysing the observations in a way that is familiar to the fluid dynamics community (pdf's, moments, structure functions), and thus could serve as a reference.

The data used in this work have been collected using a moored array of thermistors deployed above the slopes of Seamount Josephine, in the north-eastern Atlantic Ocean.
Motions at this location are dominated by the tide.
The thermistors used, built in-house at the Royal Netherlands Institute for Sea Research (NIOZ), measure temperature variations at a relatively high temporal resolution ($1\s$) for long periods of time (several months).
While their temporal resolution is generally insufficient to resolve the complete inertial range of turbulence (nor to resolve the dissipation scale), they do capture part of the turbulence cascade.
In this work, we concentrate on time-averaged statistics, distinguishing between the upslope and the downslope phases of the tide, and distinguishing different parts of the moored array, i.e.~different depth ranges.
On the other hand, we do not consider time modulation of turbulence properties other than the tidal one, i.e.~we do not consider the modulation by inertial waves, mesoscale eddies or other low-frequency processes.
Time averages are a necessary starting point, enabling better-converged statistics, and thus enabling higher-order statistics to be reliably computed.

We will consider almost exclusively fluctuations and increments of temperature within the horizontal plane.
This is motivated by the fact that the resolution of our data is substantially better in the horizontal direction ($0.2\m$ by using Taylor's hypothesis; see section \ref{sec:velocity}) than in the vertical ($0.7\m$).
The difference in resolution is particularly relevant if the anisotropy of turbulence in a strongly stratified environment is taken into account, with length scales much smaller along the vertical direction (approximately parallel to the mean density gradient) than in the horizontal (approximately perpendicular to the mean density gradient) \cite{billant_self-similarity_2001}.

The paper is organised as follows.
The temperature data used throughout the paper are described in section \ref{sec:data}.
Velocity observations are briefly described in section \ref{sec:velocity}.
The statistical observables considered are concisely presented in section \ref{sec:stats}.
The results are presented in section \ref{sec:results}; discussion and conclusions follow in section \ref{sec:conclusions}.

\section{Methods}
\label{sec:Methods}

\subsection{Temperature data}
\label{sec:data}

The data used in this work are obtained from a mooring deployed between spring and late summer 2013 in the north-eastern Atlantic Ocean on the slopes of Seamount Josephine (see table \ref{tab:mooring} for details).
The mooring has 144 ``NIOZ4'' thermistors \citep[an evolved version of the ones described in][]{van_haren_nioz3:_2009} taped on a nylon-coated steel cable at intervals of $0.7\m$.
The thermistors sampled temperature at a rate of $1\Hz$ with a precision better than $1 \,\mathrm{mK}$.
The time series from each thermistor is approximately $10^7\s$ long, with the entire dataset containing approximately $10^9$ measurements.
The cable was attached to a $500\kg$ ballast weight at the bottom, and to an elliptical buoy at the top.
The high net buoyancy (approximately $300\kg$) together with the small size of the sensors, limiting the drag on the cable, guarantee that measurements are very close to Eulerian, as checked also in similar previous moorings \citep{van_haren_high-resolution_2009,van_haren_detailed_2012}.
The pressure and tilt sensors of the acoustic Doppler current profiler (ADCP) mounted in the top buoy, $210\m$ above the bottom, recorded movements smaller than $0.2\m$ in the vertical direction and $3\m$ in the horizontal one.
All the thermistors performed satisfactorily for the whole deployment, with a noise level of approximately $5 \times 10^{-5}\dC$.
Calibration is applied to the raw data and the drift in the response of the thermistor electronics, visible over periods longer than a few weeks, is compensated for.

At the mooring location, temperature and velocity variance are dominated by the semidiurnal tidal period (in particular by the $M_2$ lunar period of $12.42\hrs$), as shown by the velocity spectrum in figure \ref{fig:velocity} (left panel).
Tidal variations are visible in the temperature record as successive cooling and warming trends.
Throughout this paper, we consider various quantities computed separately for the cooling and warming tidal phases.
For simplicity, we will usually refer to the cooling phase as ``upslope'' tide and to the warming phase as ``downslope'' tide.
The two phases are considered separately because statistical properties vary markedly between the two.
The upslope and downslope phases are defined by the sign of the time derivative (negative and positive, respectively), of the vertically averaged temperature signal, band pass-filtered at the semidiurnal frequency.
This definition produces two subsets of the data that have virtually the same size.
Two snapshots are shown in figure \ref{fig:temperature}.
The two panels show approximately $30\mn$ of data each, during downslope (labelled ``down'') and upslope (labelled ``up'') tidal phases, respectively.
Figure \ref{fig:temperature} also shows the four vertical segments of the mooring (A, B, C, D) considered separately in the rest of the work, each approximately $25\m$ high and including 36 thermistors (table \ref{tab:segments}).
From figure \ref{fig:temperature}, we can see that, during both the upslope and the downslope phases, water is stratified at the mooring location.
This is confirmed by table \ref{tab:stratification}, reporting mean, first and ninety-ninth percentiles of the buoyancy frequency ($N$) in the four segments of the mooring and during the two tidal phases separately.
Stratification is computed from $25\s$ averages of temperature, with profiles vertically reordered to obtain a statically stable stratification\protect\footnote{This procedure, used to obtain an estimate of the background stratification, removes all instabilities from the vertical profiles. On the other hand, unstratified regions are not removed from profiles and would still produce vanishing values of the buoyancy frequency, if present. See \citet{winters_diascalar_1996} for a formal justification of the reordering procedure in the context of mixing estimation.}.
This observation is consistent with those by \citet{van_haren_detailed_2012}, which extended down to within $0.5\m$ from the solid boundary, differently from those discussed here, having the deepest thermistor $5\m$ above the bottom.
More details on the computation of buoyancy frequency using these data can be found in \citet{cimatoribus_comparison_2014}.
The mean buoyancy frequency is slightly higher during the upslope tidal phase, with the exception of the lowest segment, segment A.

\subsection{Velocity data and Taylor's hypothesis}
\label{sec:velocity}

A Teledyne/RDI 75 kHz ADCP was mounted in the top buoy of the mooring, $210\m$ above the bottom.
The ADCP was pointed downwards, resolving current and echo intensity observations every $5 \m$ vertically, between $15$ and $185\m$ above the bottom.
Horizontal vector velocity is measured, but only the magnitude within the depth range of temperature measurements is considered here.
We will refer to the magnitude of velocity simply as ``velocity'' in the following, unless stated otherwise.
The sampling interval was $600 \s$.

Despite the limited resolution, velocity data can be used to transform the measured temperature from the temporal to the spatial domain, according to the Taylor's hypothesis of frozen turbulence.
Here, we must take into account that the mean velocity is not constant, following the approach of \citet{pinton_correction_1994}.
The velocity field $v$ is first decomposed into a mean component $\mean{v}$ and a fluctuating one $v'$ using a low-pass Butterworth filter with stop band $1/\sigma=3000\s$.
On average, this can be considered the forcing frequency, based on the frequency spectra of temperature fluctuations (not shown).
Approximately at this value, the spectral slope decreases moving towards higher frequencies, approaching the value of $-5/3$ typical of the inertial range of turbulence in both homogeneous and stratified environments \citep{tennekes_first_1972,lindborg_vertical_2009}.
A similar transition is observed in the velocity spectrum in figure \ref{fig:velocity}, but in this case the turbulence range is mostly undetected and contaminated by noise (light blue line in the figure).
The pdf's of $v$, $\mean{v}$ and $v'$ are shown in figure \ref{fig:velocity} (right panel).
Mean velocity is higher during the upslope phase of the tide, while the pdf's of velocity fluctuations are similar for the upslope and downslope phases.

In order to transform the temperature time series of each thermistor into the spatial domain, the mean velocity is averaged, along the vertical coordinate, within each of the four segments of the mooring.
Since velocity is not constant in time, the spatial series thus obtained has a non-constant sampling rate.
Linear interpolation is performed to a constant time step of $0.2\m$, which reduces the size of the dataset by approximately a factor of 2.

The results of this transformation have to be evaluated with care, since we cannot assess the goodness of the Taylor's hypothesis, lacking information on the turbulent velocity.
This is particularly true closer to the bottom, where turbulence and shear are stronger.
Turbulence intensity ($\mean{v'^2}^{1/2}/\mean{v}$) computed from the ADCP data is much smaller than unity most of the time.
While this gives support for the use of Taylor's hypothesis, we do not emphasise this result since the ADCP resolution is insufficient (both in time and in the vertical) to detect most of the turbulence range, and noise contaminates the highest frequency range of ADCP measurements\protect\footnote{The presence of noise artificially increases the estimated turbulence intensity.}.
In the next subsections, the results are mostly based on the spatial (Taylor-transformed) data.
The same quantities reported here have been computed also on the original (time series) data, with qualitatively similar results.
Relatively short spatial separations are considered, which correspond in the time domain to time intervals no longer than approximately one hour.

\subsection{Statistical description}
\label{sec:stats}

A power spectrum can fully characterise a dataset only if the distribution of the data is approximately Gaussian.
This is usually not the case for temperature (and scalars in general) in a turbulent flow, nor is it the case here.
For this reason, our description of the dataset is extended to higher-order moments.
In doing so, we follow the example of the analyses of laboratory data available in the literature, mostly concerning passive scalars in grid-generated turbulence; see for example the works by \citet{thoroddsen_exponential_1992} and by \citet{mydlarski_passive_1998}, and  the reviews by \citet{warhaft_passive_2000} and \citet{shraiman_scalar_2000}.
We also make use of some of the concepts developed in \citet{zhou_plume_2002} for convectively generated turbulence.
We thus apply the now-classic methods of these works to the case of a strongly stratified ocean environment, with the aim of contributing to the understanding of the latter.
In fact, we find that the results from controlled environments are often comparable to those from our observational data.

We define here a common notation for the normalised moments of a distribution, defined for a generic variable $x$ at order $q$ by
\begin{equation}
  \muq{x}{q} = \frac{\mean{\left(x-\mean{x}\right)^q}}{\mean{\left(x-\mean{x}\right)^2}^{\frac{q}{2}}},
  \label{eq:norm_mom}
\end{equation}
with the average indicated by the angular brackets.
With this notation, the third-order moment of $x$, that is the skewness, is written as $\muq{x}{3}$ and the fourth-order moment, the flatness, as $\muq{x}{4}$, and both are defined by equation \ref{eq:norm_mom}.
These two moments have received particular attention in the literature \citep{warhaft_passive_2000}, as a way to characterise intermittency and anisotropy (skewness).

In section \ref{sec:fluctuations}, we consider the temperature fluctuations, $\theta$, using the full dataset in the time domain (not using Taylor's hypothesis).
Fluctuations are defined as the deviation from the linear trend computed during a single tidal phase ($6\hrs$ interval).
We compute the skewness and flatness of the fluctuations for each thermistor in order to assess the dependence on the height above the bottom.
Results from upslope and downslope tidal phases are averaged separately, after normalisation by the standard deviation of the single tidal phase.

Using the full dataset, skewness and flatness are determined by a combination of waves and turbulence that characterises a strongly stratified environment.
In particular, there is a strong influence from the waves at periods shorter than the tidal one, not removed by the linear detrending performed before computing the moments of the distribution.
In an attempt to separate turbulence from waves, the data can be filtered.
However, there is no obvious scale separation between waves and turbulence; this distinction may in fact be ill-defined.
The fluctuations seem to be a continuum of waves of different degree of nonlinearity at frequencies lower than the (local) buoyancy frequency and turbulent motions at frequencies both below and, mainly, above $N$ \citep[see e.g.][]{billant_self-similarity_2001}.
Keeping in mind that the filtering cannot separate, or define, these components, it still proves to be an instructive exercise, discussed in section \ref{sec:fluctuations}.
We use a Butterworth filter removing frequencies below $1/5600\s^{-1}$, a frequency approximately equal to the lowest values of the buoyancy frequency (see table \ref{tab:stratification}); qualitatively similar results are obtained using a higher stop-band frequency).

If the quantity $x$ in equation \ref{eq:norm_mom} is an increment, whose value depends on the distance between two measurements (in either time or space), then $\muq{x}{q}$ depends on the same distance too, and is usually called a structure function.
In particular, in section \ref{sec:sk_and_fl} we consider the temperature increments:
\begin{equation}
    \Dtr=\theta(r_0+r)-\theta(r_0),
    \label{eq:incrm}
\end{equation}
as a function of the separation distance $r$ (assuming homogeneity), usually along the horizontal direction (i.e.~using spatial data exploiting Taylor's hypothesis).
Temperature increments are particularly useful when dealing with non-stationary time series \citep{davis_multifractal_1994}.
Increments enable one to naturally study turbulence and waves together without the need for filtering.

Further insight into the turbulence dynamics can be gained by studying the behaviour of generalised structure functions of temperature increments \citep{warhaft_passive_2000}
\begin{equation}
    \gamma_q\equiv\gamma_q(r)=\mean{\left|\Dtr\right|^q}.
    \label{eq:gsf}
\end{equation}
This is done in section \ref{sec:saturation}.
Generalised structure functions are a standard tool to characterise the intermittency of the velocity and scalar fields in a turbulent flow \citep{frisch_turbulence_1996,shraiman_scalar_2000,warhaft_passive_2000}.
Here, we obviously focus on the intermittency of the scalar (temperature).
According to the non-intermittent theory of Kolmogorov--Obukhov--Corrsin, the temperature field is expected to be fully self-similar within the inertial range, resulting in a structure function $\gamma_q$ having a scaling behaviour of the form \citep{warhaft_passive_2000}
\begin{equation}
    \gamma_q \sim r^{\zeta(q)}\mathrm{, with\; } \zeta(q) = q/3,
    \label{eq:zeta}
\end{equation}
for every order $q$ \citep[see also the discussion of this in][]{frola_experimental_2014}.
Note that a self-similar behaviour is compatible with a Gaussian distribution, but does not imply it.
For $q=2$, the structure function embeds the same information as the spectrum, with the slope $\zeta(2)=2/3$ of the structure function being the Fourier transform of the $-5/3$ slope of the fluctuations in wavenumber space.
While for $q=2$ the prediction of the theory is usually close to the values measured in the laboratory, for $q>2$ there is an increasing discrepancy between the predicted $q/3$ and the observed $\zeta(q)$, the latter being systematically below the theoretical value, possibly reaching an asymptotic value $\zeta_\infty$, i.e.~saturating, for sufficiently high orders ($q\approx 10$).
This saturation is related to the presence of ramp--cliff structures in the scalar field, in connection with an accumulation of sharp gradients at small scales.
A value of $\zeta(q)$ lower than expected is the consequence of an excess of these large gradients at small scales (the cliff structures) and of a lack of large increments at large scales (the ramp structures).
Evidence for saturation of moments has been reported in both passive scalars (shear-dominated turbulence) and active scalars (convectively driven turbulence).
However, the estimated asymptotic value of $\zeta_\infty$ is different for the two cases.
In particular, values of $\zeta_\infty\approx 1.4$ are reported for passive scalars in grid turbulence \citep{warhaft_passive_2000,celani_fronts_2001}, while smaller values ($\zeta_\infty\approx 0.8$) are observed in convectively driven turbulence \citep{zhou_plume_2002,celani_scaling_2002}.

In section \ref{sec:plus_and_minus}, we will consider ``plus'' and ``minus'' increments, introduced by \citet{zhou_plume_2002} in the context of convectively driven scalar turbulence, and used here to better investigate the passive/active nature of temperature.
They are defined as $\mathrm{\Delta}_\tau\theta^{\pm}=\left( \left| \mathrm{\Delta}_\tau\theta \right|\pm \mathrm{\Delta}_\tau\theta \right)/2$, and simply separate positive from negative increments in time (subscript $\tau$), measured at a fixed position.
A plus (minus) increment is zero for all negative (positive) increments, otherwise it is equal to the absolute value of the increment.
A hot plume rising through a thermistor produces a steep increase in temperature with its leading cap, but a relatively gentle temperature decrease with its tail.
This asymmetry between positive and negative temperature increments should show up in the skewness of the plus and minus increment distributions, with higher skewness for the plus increments (note that plus and minus increments are positive definite).
The laboratory results of \citet{zhou_plume_2002} indeed show this asymmetry in the convective region of a Rayleigh--B\'enard cell, for increments at inertial and dissipation separation distances.
More generally, here we use plus and minus increments to study if ramp--cliff structures have a preferential orientation.
In our case, convective plumes, if present, are advected through the thermistor array according to Taylor's hypothesis, rather than vertically rising through the array.
Consequently, we have to consider spatial increments of temperature rather than time increments, and we do so along both the horizontal and along the vertical directions.
Plus and minus increments are thus defined here as $\Dtr^{\pm}=\left( \left| \Dtr \right|\pm \Dtr \right)/2$ for increments in the horizontal direction, and similarly as $\Dtz^{\pm}=\left( \left| \Dtz \right|\pm \Dtz \right)/2$ along the vertical.
Here, $\Dtz$ is defined by equation \ref{eq:incrm}, with $r$ substituted by $z$, the upward vertical coordinate.
Given the relatively small number of sensors within each segment and the relatively low vertical resolution, $\Dtz^\pm$ is computed only at the thermistor spacing, i.e.~$0.7\m$.

\section{Results}
\label{sec:results}

Velocities are higher during the upslope phase of the tide (see figure \ref{fig:velocity}, right panel, black triangles).
Deviations of velocity from the vertical average explain approximately 25\% of the variance.
Shear is computed at the $5\m$ vertical resolution of the ADCP, as the average of the magnitude of the vertical gradient of the horizontal (vector) velocity.
Shear is similar during the two tidal phases, approximately equal to $3.2\times10^{-3}\s^{-1}$ (as a mean of the instantaneous shear).
Mean velocity is virtually constant in the upper half of the mooring, and decreases by approximately 20\% in the lower half.
The temporal and spatial resolution of the ADCP is such that turbulent motions cannot be studied in any detail.
Velocity has a strong component along the local isobath during both tidal phases.

Figure \ref{fig:temperature} hints at some important differences between the two tidal phases.
During the downslope phase, the sharp interface above the cold near-bottom layer is lower and thinner.
Overturns at the sharp interface are also smaller than those observed during the upslope phase.
Both these elements suggest that turbulence is weaker during the downslope phase of the tide, as already observed by \citet{van_haren_detailed_2012} and also by \citet{moum_convectively_2004} in a shallower environment.
This fact is confirmed by calculations on this dataset using similar methods to \citet{van_haren_detailed_2012} (not shown).
A peculiar feature of both panels in figure \ref{fig:temperature} is that some of the turbulent structures, in particular in the upper half of the mooring, remind one of convective plumes.
Particularly striking for its size is the mushroom-like plume during the upslope phase (black arrow in the figure).

\subsection{Temperature horizontal wavenumber spectra}
\label{sec:spectra}

As a starting point in the statistical characterisation of the temperature field, horizontal wavenumber spectra of temperature fluctuations compensated by $k^{-5/3}$ are shown in figure \ref{fig:spectra}.
The figure shows spectra averaged during the two tidal phases (upslope, downslope) and in each mooring segment.

During the upslope tide (upper four spectra in figure \ref{fig:spectra}), an inertial range is clearly visible in all mooring segments as a horizontal slope, approximately in the $k^{-1}$ range from $5$ to $100\m$.
For segments A and B, there is also a hint of a local minimum at $k^{-1}\approx 100\m$, possibly the sign of a separation between the gravity-wave spectrum (lower wavenumbers, higher slope) and the turbulence band.
High-wavenumber roll-off of the spectrum is visible above $k^{-1}\approx5\m$, but the roll-off should be interpreted with care, since the $1\Hz$ sampling rate of the thermistors is insufficient to resolve the dissipation scales (of both momentum and temperature) in the most turbulent portions of the dataset.
Spectra from single tidal phases ($6\hrs$ slices of data), show that the high-wavenumber roll-off is not always detected.
In many cases, the spectrum is flat up to almost the highest wavenumber, where the spectrum often bends upwards.
This upward bend of the spectrum is due to the signal approaching the noise level of the measuring device in weaker stratification and less turbulent regions.

In the upper half of the mooring, the spectra do not reach the theoretical $-5/3$ slope, i.e.~they slope downwards to the right in figure \ref{fig:spectra}.
We suggest that this is the result of averaging together data regions having inertial ranges spanning different scales.
This effect is particularly strong during the downslope phase, being less turbulent than the upslope one.
We conclude that the thermistors resolve at least the lower-wavenumber part of the turbulence inertial range.
Turbulence is strongly intermittent, in particular during the downslope phase of the tide, as judged from the fact that the spectral slope is steeper than $-5/3$.

Similar conclusions could be drawn from the frequency spectra (avoiding the use of Taylor's hypothesis).
However, the spatial data have spectral slopes closer to the classical value of $-5/3$ during the upslope phase, and also better highlight the differences between the two tidal phases.

\subsection{Skewness and flatness of fluctuations}
\label{sec:fluctuations}

Figure \ref{fig:sk_ku_fluctuations} shows the skewness and flatness of temperature fluctuations, $\theta$, from the full dataset in the time domain, as a function of the height above the bottom and of the tidal phase.
The skewness is always negative and small (figure \ref{fig:sk_ku_fluctuations}, left panel).
The most negative skewness appears in the region above the near-bottom cold layer, most likely in connection with mixing of the cold layer itself; it decays towards zero further above.
The flatness (right panel in figure \ref{fig:sk_ku_fluctuations}) is close to the Gaussian value of 3, in particular in the upper half of the mooring and during the upslope phase of the tide.
The Gaussian value of the flatness is clearly exceeded only close to the bottom during the downslope phase.

As discussed in section \ref{sec:stats}, filtering the data is a naive attempt to separate turbulence from wave motions.
Skewness of the high-pass-filtered data (figure \ref{fig:sk_ku_fluctuations_filt}, left panel) is positive close to the bottom, while it is negative further above.
Deviations from zero are larger than in the unfiltered data, as expected for a turbulent signal where nonlinear dynamics dominate.
Moving towards the upper end of the mooring, skewness slowly approaches zero.
After filtering, skewness is similar during the two tidal phases, indicating that the differences observed in figure \ref{fig:sk_ku_fluctuations} can mostly be ascribed to internal wave motions.
The change of sign of skewness above the bottom is the signature of the overturns that mix the interface above the cold near-bottom layer.
Skewness has slightly larger negative values during the upslope phase.

The flatness of the filtered dataset (figure \ref{fig:sk_ku_fluctuations_filt}, right panel) is much larger than in the full data; it is also larger during the downslope phase.
Flatness increases towards the bottom, in particular in the lowest $10\m$ of the mooring.
A flatness larger that the Gaussian value of 3, with exponential tails of the pdf as the ones found in our data (not shown), is a common feature of turbulent fluctuations in the inertial and dissipation ranges of both active scalars \citep{ching_probabilities_1991} and passive scalars in the presence of a mean gradient \citep{shraiman_scalar_2000,warhaft_passive_2000}.

From the comparison of filtered and unfiltered data, we conclude that even if the internal waves dominating the total signal seem highly nonlinear and continuously breaking, as suggested by figure \ref{fig:temperature}, their statistics are not too far from Gaussian.
In the filtered data, on the other hand, turbulent fluctuations are strongly non-Gaussian.
Furthermore, the skewness and flatness of the wave motions are affected by the phase of the tide more than the skewness and flatness of the turbulent motions.
Both filtered and unfiltered data depend on the height above the bottom, in particular in segment A.

\subsection{Skewness and flatness of increments}
\label{sec:sk_and_fl}

We now consider the skewness and flatness of temperature increments, defined by equations \ref{eq:norm_mom} and \ref{eq:incrm}, plotted in figure \ref{fig:sk_fl_increments} as a function of the separation distance $r$.
Here, the values are presented as averages in the four segments.
Note that $\Dtr$ has a mean different from zero, since the tidal wave produces a trend in $\theta$; the mean increment is greater or smaller than zero depending on the tidal phase (downslope or upslope respectively).

Skewness has a rather complex dependence on $r$.
During both tidal phases, skewness increases from the bottom (A) to the top (D) of the mooring.
For $r$ smaller than  $5-10\m$, skewness is higher during the upslope phase compared to the downslope phase, being positive only for segments C and D during the upslope phase.
For separation distances larger than $10\m$, skewness changes in opposite directions for the two tidal phases, towards more positive values during the downslope phase and towards more negative values during the upslope phase.
This behaviour for larger scales is due to the opposite mean trend of temperature due to the tidal wave itself.

In the inertial range, a skewness of the temperature increments different from zero is considered the signature of ramp--cliff structures \citep[][]{warhaft_passive_2000}.
These ramp--cliff structures are most visible if increments are computed parallel to the mean scalar gradient (vertical), but they have a projection also on the horizontal plane.
Here we focus on horizontal increments, owing to the limitations of the vertical resolution.
The observed behaviour of skewness shows that the anisotropy of turbulence dynamics (as measured by skewness) depends on the phase of the tide and on the height above the bottom at the smallest scales resolved in the dataset.
If skewness is determined by small-scale convective plumes \citep[producing structures that can be assimilated to ramp--cliff ones,][]{zhou_plume_2002}, the observed differences in skewness could be caused by the dominance of either warm or cold convective plumes, depending on the tidal phases.
\citet{thorpe_skewness_1991} suggested that negative skewness is associated with convection, using observations from near-surface and near-bottom oceanic boundary layers.
Here, the clearest signs of convection are found in segments C and D during the upslope tide (as will be discussed in sections \ref{sec:saturation} and \ref{sec:plus_and_minus}), having positive skewness at small separations.
The difference between our results and those of \citet{thorpe_skewness_1991} could be caused by the prevalence here of warm rather than cold plumes, a hypothesis that is discussed in more detail below.
More generally, we suggest that the skewness of the temperature increments is not a good indicator of the presence or absence of convective activity.

Flatness decreases for increasing separation distances during both tidal phases, approaching without actually reaching the Gaussian value of 3 at the largest scales considered.
A broad range of power-law decay is evident, as often observed within the inertial range in the laboratory \citep[e.g.~figure 3 in][]{warhaft_passive_2000}.
Flatness is expected to saturate for $r$ comparable to the dissipation scale, but we barely see a hint of such a saturation at the smallest $r$, confirming that the dissipation scale is not resolved here.
Flatness is very similar throughout the mooring during the downslope phase, while it is higher away from the bottom during the upslope phase.

In general, skewness and flatness strongly depend on $r$.
Skewness in particular has a rich behaviour within the turbulence inertial range identified in the spectra in figure \ref{fig:spectra}.
However, variations in skewness are generally smooth, and overall these quantities do not provide any clear suggestion of a ``scale break'' between turbulence and wave dynamics.
In the next section, we will show that, on the other hand, higher-order moments give a clear indication of such a scale break.

\subsection{High-order moments}
\label{sec:saturation}

We extend the description of the statistics of the temperature increments to higher-order moments, using in particular the generalised structure function defined in \ref{eq:gsf}.
In order to measure the convergence of the computation of these higher-order moments, we consider the quantity
\begin{equation}
    \frac{\max{\left|\Dtr\right|^q}}{M\,\gamma_q(r)},
    \label{eq:conv}
\end{equation}
where $M$ is the number of data points used to compute the average $\gamma_q$.
This ratio is the relative contribution of the largest resolved increment (at scale $r$) to the generalised structure function.
A well-converged generalised structure function implies a small value of this ratio, since the increments most relevant for the order $q$ are then well resolved, i.e.~they are smaller than $\max{\left|\Dtr\right|}$.
In order to improve the convergence of $\gamma_q$ for higher $q$, we consider the statistics averaged in two vertical subsets of the mooring instead of four as done in the previous sections.
We consider the lower half (segments A and B combined together) and the upper half (segments C and D combined together).

Based on the ratio \ref{eq:conv}, $\gamma_q$ is well converged at scales larger than $1\m$ at all orders (ratio $\mathcal{O}(10^{-2})$), with the partial exception of the highest orders in the lower half of the mooring during the upslope phase (not shown).
For scales smaller than $1\m$, the convergence is worse in particular during the upslope phase, and the orders higher than $q\approx 6$ are unreliable.
For this reason, we will focus on scales larger than $1\m$.

Figure \ref{fig:Rm} shows, in log--log space, the values of $\gamma_q$ as a function of $r$, for $q$ ranging from 1 to 10, in the lower and upper halves of the mooring and during the two tidal phases separately.
The scaling exponents $\zeta_q$ (equation \ref{eq:zeta}) are estimated in the turbulence scaling range and are shown in figure \ref{fig:slopes} as a function of the order $q$.
Considering the lower half (A+B panels in figure \ref{fig:Rm}), $\gamma_2$ has a scaling exponent close to $2/3$ over almost two decades during the upslope phase.
A scaling range with exponent $1.4$ (dotted lines) up to $r\approx 100\m$ is approached for high $q$ in the case reported in the lower left panel of figure \ref{fig:Rm}.
Figure \ref{fig:slopes}, however, shows that the scaling exponents do not saturate at this value, but keep increasing with $q$.
During the downslope phase (figure \ref{fig:Rm}, lower right panel), $\gamma_2$ has a scaling exponent higher than $2/3$, and $\gamma_q$ has an exponent larger than $1.4$ at the highest orders computed.
A clear break of the scaling behaviour is visible at $r\approx 100\m$ for $q>5$, in particular during the downslope phase.
It is unlikely that this is a consequence of the violation of Taylor's hypothesis, because the scale of the break is smaller during the upslope phase, having larger mean velocity.
We would expect on the other hand the scale of the break to be proportional to the mean velocity, if it were caused by a breakdown of Taylor's hypothesis.
Higher-order moments thus identify two separate ranges, which could be interpreted as turbulence-dominated for $r$ smaller than approximately $100\m$ and wave-dominated for larger scales.

Considering now the upper half of the mooring (C+D panels in figure \ref{fig:Rm}), $\gamma_2$ has a scaling exponent greater than $2/3$, markedly so during the downslope phase.
The scaling behaviour at high $q$ extends up to scales not as large as in segments A+B ($30$--$50\m$).
Wave motions may thus reach smaller scales in the upper half of the mooring, and turbulent motions may be limited to a narrower range.
The value of $\zeta_q$ for the highest $q$'s is well below $1.4$ during the upslope phase (see figure \ref{fig:slopes}), while it is higher than $1.4$ during the downslope one.
Considering that $\gamma_2$ has a slope in excess of $2/3$, we hypothesise that also the slopes for higher $q$ are similarly biased with respect to controlled environment results.
If this is actually the case, the asymptotic exponent during the upslope phase would in fact approach the value of 0.8 measured for an active scalar in convective turbulence \citep{zhou_plume_2002}.

Considering the scaling exponents in figure \ref{fig:slopes} more generally, we see strong deviations from the non-intermittent $\zeta_q=q/3$ prediction, marked by the dashed line in the figure.
The strongest deviations are observed during the upslope tidal phase, in particular in the upper half of the mooring (C+D).
Saturation of the moments is clearly visible in this part of the dataset, already for $q=5$.
These higher-order moments thus show that the tidal phase actually has a strong impact on intermittency, and that the turbulence behaviour is far from universal in the dataset.

\subsection{Plus and minus increments}
\label{sec:plus_and_minus}

In the previous section, we suggested that the behaviour of temperature is similar to that of a passive scalar in the lower half of the mooring, while in the higher part the statistics deviate towards those of active scalars, in particular during the upslope phase.
To further investigate this issue, we consider the skewness of plus and minus increments, introduced in section \ref{sec:stats}.
The skewness of the distributions of $\Dtr^\pm$ and $\Dtz^\pm$ are shown in figure \ref{fig:sk_increments_pm}, as a function of the separation distance and for the four segments of the mooring separately.

During the downslope phase, the skewness of the vertical minus increments, $\muq{\Dtz^-}{3}$ (large, filled red triangles), is greater than $\muq{\Dtz^+}{3}$ (large, filled red squares) in all segments.
During the upslope phase of the tide this is also true, but the difference is much smaller in the lower segments that in  upper ones.
Considering the horizontal increments during the downslope phase of the tide, the skewness of the minus increments (small, red filled triangles) is larger at all scales and throughout the mooring.
In the upslope phase closer to the bottom (segments A and B), the skewness of plus and minus horizontal increments are the same, or close, with the exception of larger scales.
On the other hand, plus skewness (small, blue open squares) is larger in segments C and D at all scales.

The difference in skewness between plus and minus increments is of order 1, consistent with that reported in the laboratory by \citet{zhou_plume_2002}.
On the other hand, laboratory data show that $\muq{\mathrm{\Delta}_\tau\theta^+}{3}$ and $\muq{\mathrm{\Delta}_\tau\theta^-}{3}$ converge to the same value for large $\tau$.
This does not happen here, where an increase of the difference is observed for $r$ larger than $10$--$100\m$.
We should not conclude from this that plumes of all sizes are present, but rather that wave motions are also associated with asymmetric temperature gradients.
We also cannot rule out that contamination between the two tidal phases and the use of the Taylor's hypothesis may play a role at the largest scales.

These results are discussed with the help of the sketch in figure \ref{fig:plumes}.
We consider two ``plumes'', which we regard here generically as an asymmetric temperature anomaly, with a sharp leading front and a smoother tail.
In the upper part of the figure, the isotherms of two of these plumes are sketched in the spatial domain: (1) on the left a hot plume with an upper sharper front; and (2) on the right, a cold plume with a lower sharper front.
Two transects are taken through each of the two plumes, one vertical (green dashed arrow) and one horizontal (red continuous arrow).
The resulting temperature transects are sketched in the two panels below, with green dashed and red continuous lines for vertical and horizontal directions, respectively.
Consider the hot plume on the left.
Its leading front is sharper than its tail; the vertical transect thus shows a slower increase of the temperature and a faster decrease to the mean value.
This structure of the plume shows up as larger values of skewness for minus rather than plus vertical increments, up to separations comparable with the size of the plume.
The same reasoning applies to the horizontal direction and for a plume with an inclination to the mean flow given in the figure.
The result is a larger skewness for plus than for minus horizontal increments.
Note that a cold plume with opposite tilt would produce different results for the vertical increments, so the asymmetry is characterised unambiguously.
This is consistent with the estimates in segments C and D during the upslope tidal phase.
A similar reasoning can be applied to the cold plume in figure \ref{fig:plumes}, suggesting that in this case $\muq{\Dtz^+}{3}<\muq{\Dtz^-}{3}$ and $\muq{\Dtr^+}{3}<\muq{\Dtr^-}{3}$.

We can now interpret the results from figure \ref{fig:sk_increments_pm} as follows.
During the upslope phase, warm plume-like temperature anomalies are dominant in segments C and D, well above the bottom.
During the downslope phase, cold plume-like anomalies are instead dominant throughout the mooring.
The asymmetry of temperature anomalies has both a vertical and a horizontal component.
In other words, plumes also have a preferential horizontal orientation, pointing towards the direction of the mean flow, as shown in figure \ref{fig:plumes}.
On the other hand, in segments A and B during the upslope tide, temperature anomalies are more symmetric, in particular along the horizontal direction.
Overall, we do not have any clear indication of the size of the plumes.
We can only speculate that their typical size is at least of $\mathcal{O}(1\m)$, i.e.~at least comparable to the values of $r$ for which skewness departs from its saturated, small-scale value.
The large plume of figure \ref{fig:temperature} is thus an extreme case, and scanning through the temperature record confirms this.

\section{Discussion and conclusions}
\label{sec:conclusions}

We described the temperature statistics using observations from a moored thermistor array, deployed in the deep ocean above a sloping bottom.
The dynamics at the mooring location are strongly affected by the semidiurnal tidal wave, and we thus considered separately the statistics for two tidal phases, cooling (upslope) and warming (downslope).
The thermistor array does not reach the bottom, with the lowest sensor approximately $5\m$ above the bottom, and a bottom mixed layer, if present, is not visible in the data.
The results are one of the few studies on the high-order scalar statistics in stably, strongly stratified fluids.

In the dataset, temperature fluctuations and temperature increments are strongly non-Gaussian at (horizontal) scales smaller than $10$--$100\m$.
At larger scales, deviations from Gaussian are smaller, but there generally is still a scaling behaviour, with a higher value of the scaling exponent $\zeta_q$.
The surprising result here is that high-order moments show an abrupt transition between the two regimes (see in particular panels A+B up and down, C+D down in figure \ref{fig:Rm}).
The length scale of the transition is of order $100\m$, but the actual value is state-dependent, with significant variations depending on the tidal phase and on the height above the bottom.
The $100\m$ scale is equal to the buoyancy length obtained as the ratio of the velocity scale ($10^{-1}\m\s^{-1}$; see figure \ref{fig:velocity}) and the buoyancy frequency ($10^{-3}\s^{-1}$, see table \ref{tab:stratification}), and is approximately 100 times larger than the Ozmidov scale \citep[estimated from the Thorpe scale, averaged over the whole dataset, computed in][]{cimatoribus_comparison_2014}.
Since the layer thickness scales with the buoyancy length in a stratified turbulent environment \citep{billant_self-similarity_2001}, we hypothesise that motions at larger horizontal scales are less nonlinear, i.e.~more wave-like, as overturns are limited by the layer thickness.
In this view, we argue that this break in the scaling behaviour could be a characteristic feature of stratified turbulence.
The fact that a sharp break is visible only in the high-order moments, at least in this dataset, could result from the fact that the buoyancy length is relevant only for turbulent motions.
On the other hand, internal waves are scale-free, at least according to the dispersion relation from linear theory, and thus have scales both smaller and larger than the break.
Since waves are less effective than turbulent motions at producing sharp gradients, their signature on higher-order moments is smaller than for turbulent fluctuations, and thus smooth the transition between the two regimes only in the lower-order moments.
The possibility to define a (state-dependent) scale separation between waves and turbulence may have important practical applications, e.g.~for large eddy simulations \citep{frola_experimental_2014}.
It is also noteworthy that, assuming we can indeed identify large-scale motions with waves, deviations from Gaussian are small even if there is ample evidence of nonlinearity and breaking.

We also considered generalised structure functions and skewness of plus/minus increments with respect to the passive/active nature of temperature dynamics.
There is a clear indication of a passive scalar behaviour during the upslope phase close to the bottom; shear is strong enough to dominate the buoyancy forces.
On the other hand, in the upper half of the mooring, temperature statistics approach those of an active scalar during the upslope phase.
The difference in skewness of plus/minus increments is consistent with the presence of warm convective plumes.
The numerical simulations of \citet{slinn_modeling_2003} and \citet{gayen_boundary_2011,gayen_negative_2011} suggested that tidal motions above slopes can produce strong convection and mixing, in particular during the upslope tidal phase.
Our results provide support for their conclusions.

During the downslope phase, results are more ambiguous.
The difference in skewness of plus/minus increments indicates an asymmetry of temperature anomalies consistent with the presence of cold plumes throughout the mooring.
However, the scaling exponents of the generalised structure function exceed those of passive scalars (figure \ref{fig:slopes}).
The fact that in this case low-order moments exceed the $q/3$ prediction suggests, however, that this less turbulent tidal phase is strongly inhomogeneous, and a comparison with laboratory results is more difficult.
We hypothesise that the observed scaling behaviour is typical of strongly stratified, sporadically turbulent flows.
In this respect, it would be useful to compare our results with numerical simulations similar to those reported in \citet{rorai_turbulence_2014}, if scalar statistics can be computed from them.

In conclusion, we want to return to the initial, broader, issue of mixing in geophysical flows.
Overall, the indication from our results is that temperature statistics are highly dependent on the state of the flow, even if stratification is always present.
Differences in the statistics, in particular in relation to convective motions, are very likely to be connected to changes in mixing properties, in particular mixing efficiency.
A clearer understanding of mixing in a stably stratified flow is an important open question with implications in the fields of oceanography and climate.
We hope that the present results can help in this effort.\\

  The authors would like to thank the crew of RV~{\em Pelagia} for enabling the deployment and recovery of the instruments, and all the technicians, in particular Martin Laan, indispensable for the success of the measurements.
The authors would like to thank the referees for their useful and detailed comments on the previous versions of the manuscript.

redpajama_set_name: RedPajamaArXiv


text:
\section{Conclusion}
\vspace{-2mm}
We described a pipeline that combines deep features, explicitly
trained to discriminate between different poses with the poselets
framework of \cite{Bourdev09} to generate high quality
bounding box proposals for people. Using the RCNN person detector of
\cite{Ross}, we then score these proposals and refine the bounding box
coordinates. The system produces state-of-the-art performance, beating the RCNN
detector alone by 0.6--1.2\% on the highly challenging PASCAL
datasets. Surpassing the RCNN detector is significant,
given that it gave a $\sim$10\% AP gain over the previous state-of-the-art.
The results show the importance of maintaining the correct aspect
ratio when classifying regions and give further evidence for the power
of deep features over hand-crafted features such as HOG.
\section{Experiments}

We start by comparing the performance of our deep poselet features to
the HOG-based poselets, introduced in \cite{Bourdev09}.

\subsection{Evaluation of Feature Representation}
\label{sec:feat_eval}
For a given seed window we extracted up to 500 positive examples for
several poselets capturing profile face, legs, etc. The examples are
pulled from the poselets training set which includes keypoint
annotations. We then computed HOG features and deep features coming
from fc1 of our CNN. HOG features had 1476 dimensions whereas the deep
features had only 256 dimensions.  We randomized their indices and
split them into training and test sets using 25\% of the positive and
25\% of the negative examples as test and the rest as training and
trained a linear SVM as a binary classification problem (the poselet
vs the background). We computed the mean Average Precision on the test
set (which has 125 positive and 2500 negative examples). We measured
performance by training on the full training set as well as on smaller
portions of it. The results are shown in \fig{deep_vs_hog} for four
different poselets (profile face, hands in front of body, legs, head
and shoulders).

\begin{figure}[h!]
\begin{center}
\includegraphics[width=4.5in]{deep_vs_hog.pdf}
\end{center}
\vspace*{-0.3cm}
\caption{A comparison between HOG-based poselets (dashed) and our deep
  feature-based poselets (solid), for 4 different types of
  poselet. The deep features give superior performance, particularly
  when only a few training examples are available. }
\label{fig:deep_vs_hog}
\vspace*{-0.3cm}
\end{figure}

\subsubsection{Resilience to Jittering}

We also explored the performance of our learned representation to
small misalignments that will occur at test time. We use the frontal
face poselet and generate jittered versions which have random
rotations of +/-20 degrees, scalings of 0.7 to 1.3 and translations of
+/- 16 pixels (out of 64) in both x and y dimensions. The ranges are sampled
uniformly for each poselet example. \tab{jitter} shows the performance
of the HOG and learned representations as the number of training
examples is varied.

\begin{minipage}{\textwidth}
  \begin{minipage}[b]{0.35\textwidth}
    \centering
\includegraphics[width=1.5in]{jitter.pdf}
    \captionof{figure}{Some jittered examples.}
  \end{minipage}
  \hfill
\vspace{-3mm}
  \begin{minipage}[b]{0.60\textwidth}
    \centering
    \small
\vspace{-3mm}
    \begin{tabular}{| c | c | c | c |}
    \hline
Pos. & Neg. & HOG AP \% & Deep AP \% \\ \hline
375 & 7500 & 70.59 & 99.44 \\ \hline
187 & 3750 & 64.87 & 98.92 \\ \hline
93 & 1875 & 59.1 2 & 98.47 \\ \hline
46 & 937 & 55.79 &  97.13 \\ \hline
23 & 468 & 48.84 & 95.22  \\ \hline
11 & 234 & 45.6 & 92.94 \\ \hline
5 & 117 & 32.45 & 88.35 \\ \hline
2 & 58 & 29.83 & 82.58 \\
    \hline
\end{tabular}
\label{tab:jitter}
      \captionof{table}{Comparison between HOG-based and deep feature-based poselets
  for a jittered test set. }
    \end{minipage}
  \end{minipage}
\vspace{3mm}

In summary, our learned representation:
\begin{itemize}
\item Gives significantly higher accuracy than HOG.
\item Perform very impressively with few training examples. This is
  particularly important because some poselets are very rare and we
  cannot train HOG-based models on them.
\item Is more compact (1476 vs 256 dimensions) and therefore SVM evaluation is faster
\item Is more resilient to misalignment of the data. This is important
  because it implies we can run them at coarser step and the system
  could run much faster. However, they may be less precise than HOG in
  predicting the keypoints of the detected person.
\end{itemize}

Note that here we did not try fine-tuning the poselet network, instead
just training a linear SVM on the fixed final-layer
features. Furthermore, the learned features were also trained on a different
dataset to recognize different sets of poselets, albeit of very
similar patterns.

\subsection{PASCAL Evaluation}

We evaluate our algorithm of the PASCAL VOC datasets \cite{PASCAL}, which show
people in a wide range of natural settings. The large variation in
scale, occlusion and pose make the dataset highly challenging.
These datasets have also been used to evaluate many of the leading
approaches. In particular, the recent work of Girshick {\textit{et~al.~}}
\cite{Ross} which also make used of convolutional networks, has the
current leading performance on the person detection task with a score
of 58.0\% average precision. As shown in \tab{person_ap}, our
algorithm outperforms them by 0.6-1.2\%.

\begin{table}[h]
\begin{center}
    \begin{tabular}{| c | c | c | c |}
    \hline
Test set  &  HOG Poselets & RCNN & Deep Poselets        \\ \hline
VOC 2007  & 46.9 & 58.7  &  \textbf{59.3}     \\
VOC 2010  & 48.5 & 58.1  &  \textbf{59.3}     \\
VOC 2011  & - & 57.8  &  \textbf{58.7}      \\ \hline
    \end{tabular}
\end{center}
\caption{mAP of person detection compared with R-CNN.}
\label{tab:person_ap}
\end{table}

Note that we achieved this high performance even though our poselet model is significantly simplified. In~\tab{differences_with_hog} we show the differences of our method vs the original poselets of~\cite{BourdevPoseletsECCV10}. The original poselets have various enhancements -- multiple aspect ratios, contextual poselets, dynamic selection of the number of training examples per poselet, etc. Each of these choices adds about 1-2\% to the accuracy but at the expense of making the system more complicated and slower to train and evaluate. We have decided to disable these extensions favoring simplicity over higher performance. We used a fixed square aspect ratio and our patch resolution is on average 40\% smaller. For object level classification and regression we use the pre-trained CNN of~\cite{Ross}. In this work we use the CNN without fine-tuning it for our task.

\begin{table}[h]
\begin{center}
    \begin{tabular}{| c | c | c |}
    \hline
  &  HOG Poselets~\cite{BourdevPoseletsECCV10} & Deep Poselets        \\ \hline
Poselet classifier features & HOG & PDF \\
Pixels per patch & 4096 - 8192 & 3721 \\
Aspect ratios & 64x96 64x64 96x64 128x64 & 61x61 \\
Contextual poselets (Q-poselets) & enabled & not used \\
num training examples per poselet & dynamic 200 - 1000+ & 200 \\
Object-level classification features & poselet scores & FC7 layer of R-CNN~\cite{Ross} \\
Object-level regression features & consensus bounds & Pool5 layer of R-CNN~\cite{Ross}\\ \hline
    \end{tabular}
\end{center}
\caption{Detailed differences from the HOG poselets of~\cite{BourdevPoseletsECCV10}.}
\label{tab:differences_with_hog}
\vspace{-3mm}
\end{table}

\begin{figure}[h!]
 \begin{center}
\includegraphics[width=2.32in]{screen1.png}
\includegraphics[width=2.5in]{screen2.png}
\end{center}
 \vspace*{-0.3cm}
 \caption{Some detection examples. }
 \label{fig:ex}
 \vspace*{-0.3cm}
 \end{figure}

\section{Introduction}

Detecting humans in natural scenes is challenging due to the huge
variability in appearance, pose and occlusion. The recent state-of-the-art approaches of Girschik {\textit{et~al.~}} \cite{Ross} and Sermanet {\textit{et~al.~}} \cite{Sermanet13} differ mostly based on the algorithms they use to select candidate locations -- sliding windows or bottom-up segmentation.

The R-CNN
work of Girschik {\textit{et~al.~}} \cite{Ross} shows how localization can be
performed by first proposing possible object locations and then
presenting each in turn to a convolutional network classifier. Although this
approach is somewhat slow, it currently has the leading performance
Imagenet and PASCAL detection tasks.

The OverFeat system of Sermanet {\textit{et~al.~}} \cite{Sermanet13} also uses a
convolutional network, but scans it over the entire image, rather than
a small set of locations. The top part of the network directly
predicts the coordinates of the bounding box. The proposed bounding
boxes are agglomerated into distinct object detections, which yields very
good performance on the Imagenet detection task.

\begin{figure}[h!]
\begin{center}
\includegraphics[width=5in]{teaser.pdf}
\end{center}
\vspace*{0cm}
\caption{\textbf{Left:} The two state-of-the-art methods handle differently the problem of variable aspect ratios to fit them into the fixed-dimensional CNN. R-CNN's distorting the data makes the problem harder as the system has to learn everything under various degrees of distortions that don't occur in the natural world. OverFeat's solution doesn't distort the data but introduces sometimes lots of noisy background. Neither methods leverage the natural correspondences of parts.~\textbf{Right:} Our method doesn't distort the image and doesn't include unnecessary background. By factoring the pose and viewpoint from the appearance, poselets significantly simplify the variability of the positive class for each part and makes training part classifiers easier. Furthermore, poselets leverage the semantic correspondences in the input images, which which makes it easier to do high-level semantic reasoning. }
\vspace{-5mm}
\label{fig:regions}
\vspace*{0cm}
\end{figure}

One of the challenging aspects of person detection is that the visible portion of the person can have an arbitrary aspect ratio and it is hard to fit into the fixed dimensional input of the convolutional neural net. OverFeat and R-CNN take different approaches. In the case of OverFeat, the square scanning window approach can result in much background being included, and other times fitting only part of the person in the scanning window. To detect people, the holistic system of OverFeat has to learn to detect a huge number of patterns -- arbitrary parts of a person under arbitrary viewpoints, poses and appearances, which is particularly challenging because it has to learn to distinguish each from arbitrary non-person patterns. In the case of R-CNN the region proposals come from bottom-up segmentation and can have different aspect ratios. To fit them into the fixed size dimensions, the method resizes them, which sometimes results in very distorted patterns.  So the training examples for R-CNN are made more complex by including patterns that do not exist in the natural world, which makes it harder for the system to learn the correspondences across instances of the same object. Since convolutional nets are not aspect-ratio invariant, R-CNN would have to include plenty of examples with each possible aspect ratio to compensate for it. Furthermore, bottom-up methods are constrained by the ability of the segmentation system to detect suitable regions for the target object. If the segmentation fails to propose a suitable candidate region for a person, the CNN will not be given the opportunity to classify it.

We propose a third approach for detection -- a system that does not have to train with distorted aspect ratios or handle a large variability of positive examples within the same classifier.  Our method is inspired by the {\em Poselets} approach of Bourdev {\textit{et~al.~}} \cite{Bourdev09}. Poselets are parts trained to respond to parts of a person under specific viewpoint and pose. Examples of poselets are shown in \fig{teaser}. Instances of each poselet type are detected in the test image, each of which votes on the location of the person. Poselet detections with compatible votes are clustered together into a person hypothesis. By constraining the training set of each poselet to examples of the same part of the body under the same viewpoint, the variability of the positive class of each poselet classifier is significantly reduced, training is easier, and the aspect ratio need not change. These benefits allowed poselets to maintain state-of-the-art results on the PASCAL person detection tasks until very recently before being overtaken by convolutional neural net methods, which, as numerous recent works demonstrate, work much better than the HOG features traditionally used by poselets.

We propose a deep learning version of the poselet architecture. CNNs require lots of training examples to avoid overfitting but building a large set of training examples for poselets is expensive. \emph{Our first contribution is a bootstrapping scheme} where we use the traditional poselets to collect lots of examples, automatically provide weak labels for them, and use this much larger training set to train deep poselet classifiers. Pose combinations have a heavy-tailed distribution. Lots of patterns are rare and it is hard to collect enough training examples to train a deep classifier for them. Furthermore, different vision tasks (detection, pose estimation, attribute extraction, segmentation) may require different set of optimal poselets. Training a deep poselet for each of them is time consuming. \emph{Our second contribution is using the deep representation to provide a generic pose-discriminative compact feature vector} of only 256 dimensions. We show that this feature vector allows us to quickly, and with very few examples, train an effective poselet classifier that the original HOG-based method was unable to do. Our method achieves state-of-the-art results on the challenging PASCAL datasets for the person detection category.

\begin{figure}[h!]
\begin{center}
\includegraphics[width=1.7in]{poselets_ex1.pdf}
\includegraphics[width=1.7in]{poselets_ex2.pdf}
\includegraphics[width=1.7in]{poselets_ex3.pdf}
\end{center}
\vspace*{-0.3cm}
\caption{Examples of poselet detections: legs of a front-facing person, back-facing head-and-shoulders, and a semi-profile face.}
\label{fig:teaser}
\vspace*{-0.3cm}
\end{figure}

\subsection{Related work}

Notable early approaches to object detection focused on faces~\cite{Viola2001} and pedestrians~\cite{DalalTriggs}. A natural approach is to find individual parts of the person and then combine them into an overall detection. An early example is the constellation model of~\cite{Fergus03} whose parts location and appearance are jointly trained via EM. The implicit shape model of~\cite{leibe2004combined} employs a variable number of parts that use Hough voting to suggest the object location. A more recent popular example is the deformable parts
model by Felzenszwalb {\textit{et~al.~}} \cite{dpm}, which consists of a small set of
part detectors plus a relative geometry term. Different views are
handled as components in a mixture model. Given bounding boxes in
training, the part appearances can be learned using a latent SVM
framework. However, even with several mixture components, the number
of distinct parts in the model is relatively small, i.e.~ a few dozen. Other notable approaches include: Filder {\textit{et~al.~}} \cite{Fidler13}, Uijlings {\textit{et~al.~}} \cite{Uijlings}, Dollar {\textit{et~al.~}} \cite{dollar2010fastest}, and the Regionlets approach of Wang {\textit{et~al.~}} \cite{wang}.

Lecun {\textit{et~al.~}} introduced convnets \cite{lecun} for hand-writing recognition, but they have recently been shown to be very effective across a range of vision tasks -- classification~\cite{Krizhevsky12}, detection~\cite{Szegedy13,Sermanet13,sermanet-cvpr-13}, scene labeling~\cite{Farabet12}, face recognition~\cite{taigman14}, and attributes~\cite{panda} to name a few. Several works \cite{Ross,Josef,Matt}
have shown how deep networks can learn features that are far more
powerful than hand-crafted features such as HOG or SIFT.