\section{Introduction}
The gravitational collapse for a spherically symmetric stellar object has been analyzed in detail in the standard General Relativity (GR) theory \cite{Weinberg}. By assuming the metric of the space-time to be spherically symmetric and that the collapsing fluid is pressureless, the found metric interior to the object turns to be Robertson-Walker type with a parameter playing the role of spatial curvature and proportional to initial density. The time lapse and the size of the object are given by  a cycloid parametric equation with an angle parameter $\psi$.

In spite of the fact that GR has been one of the most successful theories of the twentieth century, it does not give a satisfactory explanation to some of the latest cosmological and astrophysical observations with usual matter sources. On the one hand, the baryonic matter content has to be supplemented by a dark matter (DM) component. Although there are many possible origins for this additional component \cite{DM}, DM is usually assumed to be in the form of thermal relics that naturally freeze-out with the right abundance in many extensions of the standard model of particles \cite{WIMPs}. Future experiments will be able to discriminate among the large number of candidates and models, such as direct and indirect detection designed explicitly for their search \cite{isearches}, or even at high energy colliders, where they could be produced \cite{Coll}.

On the other hand, a dark energy contribution needs to be considered to provide cosmological acceleration \cite{DE}. Alternatively,
a modified gravitational interaction can potentially generate the present expansion of the Universe without invoking the presence of any exotic dark energy \cite{otras}. In this context, fourth order theories of gravity defined by adding to the action a function of the scalar curvature:
\begin{equation}
S_G=\frac{1}{16\pi G}\int \text{d}^4x \sqrt{\vert g \vert}\left(R+f(R)\right)\,,
\label{Modified action}
\end{equation}
has been extensively studied during the last years \cite{R2a,Rm1,varia,R2b}.

\section{Gravitational collapse}
In order to study the gravitational collapse inside this modified gravity theories, so-called $f(R)$ theories, we can make use of the spherically symmetric metric. The collapsing object is considered to be approximated pressureless and homogeneous. Therefore, we can search a separable solution for this metric as follows \cite{Weinberg}:
\begin{eqnarray}
&&\text{d}s^2 =\text{d}t^2-U(r,t)\text{d}r^2-V(r,t)(\text{d}\theta^2+\text{sin}^2\theta \text{d}\phi^2),\\
&&U(r,t)=A^2(t)h(r),\quad V(r,t)=A^2(t)r^2.
\label{Intervalo colapso, U y V}
\end{eqnarray}
Components $tt$, $rr$ and $\theta\theta$ for the modified tensorial equations may be written respectively in terms of the functions $A(t)$ and $h(r)$. Once we have calculated $h(r)$, the resulting metric is formally the same as the one obtained in GR. We consider $\rho\,(t)=\rho\,(t=0)/A(t)^3$ (given by the energy motion equation for dust matter). Furthermore, provided that the fluid is assumed to be at rest for $t=0$, initial conditions $\dot{A}(t=0)=0$ and $A(t=0)=1$ hold. We define $R(t=0)\equiv R_0$ and $\rho\,(t=0)\equiv \rho_0$ to simplify notation. Taking into account the previous conditions, the evolution equation for $A$ is found to be \cite{Cembranos:2012fd}:
\begin{eqnarray}
\dot{A}^2&=&-\frac{1}{6(1-f_R^2(R_0))}\big[\vphantom{\ddot{f}_R(R_0)}\,8\pi G\rho_0\left(2-f_R(R_0)\right)
- f(R_0)(1+f_R(R_0))
\\
&&
-3\ddot{f}_R(R_0)f_R(R_0)\,\big]
+\frac{1}{(1-f_R^2)}\,\frac{8\pi G}{3}\rho_0\,A^{-1}
-\frac{1}{6(1-f_R^2)}\big[8\pi G\rho_0\,A^{-1}f_R
\nonumber\\
&&
+3A^2\ddot{f}_{R}f_R-3A\dot{A}\dot{f}_{R}(1-f_R)+A^2f(R)(1+f_R)\big]\,.\nonumber
\label{Simplified Eq for dotA}
\end{eqnarray}
Let us remind at this stage that the zeroth order solution of GR is given by the parametric equations of a cycloid \cite{Weinberg}:
\begin{equation}
t=\frac{\psi+\text{sin}\,\psi}{2\sqrt{k}},\;\;\; A_G=\frac{1}{2}(1+\text{cos}\,\psi).
\label{Eqs cicloide}
\end{equation}
Expression \eqref{Eqs cicloide} clearly implies that a sphere with initial density $\rho_0$ and negligible pressure will collapse from rest to a state of infinite proper energy density in a finite time that we will denote $T_G$. This time is obtained for the first value of $\psi$ such as $A_G=0$, i.e. for $\psi=\pi$. In order to study the modification to the gravitational collapse in $f(R)$ theories, we have expanded $A$ around $A_G$ and $f(R)$ around the scalar curvature in GR ($R=R_G$) \cite{Cembranos:2012fd},i.e. $A=A_G+g(\psi)$, and $f(R)\simeq f(R_G)+f'(R_G)(R-R_G)$.
The presence of a function $f(R)$ in the gravitational Lagrangian represents a correction of first order with respect to the usual Einstein-Hilbert action. By substituting the above series expansions in expression \eqref{Simplified Eq for dotA} until first order in $\varepsilon$,  we find that equation (\ref{Simplified Eq for dotA}) becomes \cite{Cembranos:2012fd}:
\begin{eqnarray}
\text{tg}\left(\frac{\psi}{2}\right)g'\,&=&\,-\frac{1}{2}\text{cos}^{-2}\left(\frac{\psi}{2}\right) g
+\frac{1}{12k}\text{cos}^{2}\left(\frac{\psi}{2}\right)\left(f(R_{G0})+3 k f_{R}(R_{G0})\right)
\\
&-&\frac{1}{4}f_R(R_G)-\frac{1}{12k}\text{cos}^{6}\left(\frac{\psi}{2}\right) f(R_G)+\frac{1}{4\sqrt{k}}\,\text{sin}^{}\left(\frac{\psi}{2}\right)\text{cos}^{3}\left(\frac{\psi}{2}\right)\dot{f}_{R}(R_G)
\,. \nonumber
\label{Eq g' with generic Rs}
\end{eqnarray}
By studying this equation for particular $f(R)$ models, it is possible to constraint an important region of the parameter space of
this gravitational theories (See Fig. \ref{validity}). This validity tests are interesting since the higher contraction that these models generally present in the perturvative regime, can alleviate some tension of the standard $\Lambda$CDM model with observations of structures at high redshift \cite{clusters}.
\begin{figure}[tp]
	\centering
		\includegraphics[width=0.48\textwidth]{Validityregimemodel1.eps}\;\;\;\;
		\includegraphics[width=0.48\textwidth]{Validityregimemodel2.eps}
		\caption{\footnotesize{
Validity of the perturbative regime for $f(R)= \epsilon R^2$ \cite{R2a,R2b} (left panel) and $f(R)= \epsilon R^{-1}$ \cite{Rm1} (right panel). As it is indicated in the figures, in the first case, it is the high density region, the one that exhibits a non linear regime; whereas it is the low density region for the second model (on the right). The right area on the left panel is excluded by the condition $\varepsilon\leq 2.3\times10^{22}\,\text{GeV}^{-2}$ \cite{R2b}. Finally, the density marking the beginning of structure formation (SF) $\rho_{\text{SF}}\simeq1.5 \times 10^{-38}$ GeV$^4$, and the dark energy (DE) density $\rho_{\text{DE}}\simeq2.8 \times 10^{-47}$ GeV$^4$, have also been plotted for reference \cite{Cembranos:2012fd}.}}
\label{validity}
\end{figure}

\begin{theacknowledgments}
This work has been supported by MICINN (Spain) project numbers FIS 2008-01323, FIS2011-23000, FPA 2008-00592, FPA2011-27853-01, Consolider-Ingenio MULTIDARK CSD2009-00064 (Spain) and URC (South Africa).
\end{theacknowledgments}

\bibliographystyle{aipproc}


\section{Introduction and the main result}\label{sec_introduction}

For any strictly convex flux function $A\in C^2(\mathbb{R})$, we consider the family of
 viscous scalar conservation laws in one space dimension:
\begin{equation}\label{scalar_cl}\begin{cases}
&\partial_tU+\partial_xA(U)=\varepsilon\partial_{xx}^2U\quad \mbox{ for } t>0,x\in\mathbb{R},\\
 &U(0,x)= U_0(x)  \mbox{ for } x\in\mathbb{R},
\end{cases}
\end{equation}
for any $\varepsilon>0$ and  $U_0\in L^\infty$.
Global unique solutions to (\ref{scalar_cl}) have been constructed by
Hopf \cite{MR0047234}  and
 Ole{\u\i}nik \cite{MR0094541}. The inviscid case, $\varepsilon=0$, is covered by the theory of
Kru{\v{z}}kov \cite{Kruzkov}. Kuznetsov   showed in \cite{Kuznetsov} that, for fixed initial
 values $U_0$, the solutions of (\ref{scalar_cl})
 converges, when $\varepsilon$ goes to zero,  to the solution of
 the inviscid scalar conservation law (the equation  (\ref{scalar_cl}) with $\varepsilon=0$).
 He showed also that, for general initial data, the optimal rate of convergence in $L^1$
 is $\sqrt{\varepsilon}$.  In the case of  the convergence to a shock, however, the rate is
 better. This is, usually,  linked to the formation of layers.
  \vskip0.3cm
 In this paper we consider the asymptotic limit for general initial values.
We are particularly interested in the cases where the initial values carry
too much entropy  for the structure of the layer to be preserved asymptotically.
The shocks solutions of the inviscid case ($\varepsilon=0$) can be described as follows.
Consider two constants $C_L>C_R$, and the associated function defined by
 \begin{equation}\label{def_S}
 S_0(x)=\begin{cases}
&C_L \mbox{ if } x<0,\\
&C_R \mbox{ if } x\geq 0.
\end{cases}
\end{equation}
Then, the Rankine-Hugoniot conditions ensures that the function
\begin{equation}\label{def_sigma}
S_0(x-\sigma t) \qquad \mbox{with }\quad\sigma:=\frac{A(C_L)-A(C_R)}{C_L-C_R},
\end{equation}
is a solution to the inviscid equation (\ref{scalar_cl}) with $\varepsilon=0$.
The condition $C_L>C_R$ implies that they verify the entropy conditions, that is:
$$
\partial_t \eta(U)+\partial_x G(U)\leq 0,\qquad t>0, \ x\in \mathbb{R},
$$
for any convex functions $\eta$, and
$G'=\eta'A'.$\\

Our main result is the following.
\begin{thm}\label{main_thmL2}
Let $C_L>C_R$  and $U_0\in L^\infty(\mathbb{R})\cap BV_{loc}(\mathbb{R}) $ be such that
$$
(U_0-S_0)\in L^2(\mathbb{R})\qquad\mathrm{ and}\qquad
  (\frac{d}{dx}U_0)_+\in L^2(\mathbb{R}).
 $$
 Then, for any $T>0$, there exist $\varepsilon_0>0$  and $C^*>0$ such that  the followings
 hold true.
 \begin{itemize}
 \item[I.]
For any  $U$ solution to
 \eqref{scalar_cl} with $0<\varepsilon\leq \varepsilon_0$,
there exists a curve $X\in L^\infty(0,T)$ such that $X(0)=0$ and  for any $0<t<T$:
\begin{equation}\begin{split}\label{main_estiamteL2}
 &\|U(t)-S(t)\|^2_{L^2(\mathbb{R})}\leq
\|U_0-S_0\|^2_{L^2(\mathbb{R})}
 +C^*  \varepsilon\log({1}/{\varepsilon}),
\end{split}\end{equation}
 where $S(t,x):=S_0(x-X(t))$, and $S_0$ is defined by (\ref{def_S}).

\item[II.] Moreover, this curve satisfies
\begin{equation}\label{estimate_curve_lishitzL2}
|\dot{X}(t)|\leq C^* \quad\mbox{ and}
\end{equation}
\begin{equation}\label{estimate_curve_w.r.t._L^2_differL2}
|X(t)-\sigma t|^2\leq C^* t^{2/3}
\left(\|U_0-S_0\|^2_{L^2(\mathbb{R})}
+\varepsilon\log({1}/{\varepsilon})\right).
 \end{equation}

\item[III.] The constant $\varepsilon_0$ depends only on
$\|(\frac{d}{dx}U_0)_+\|_{L^2}$, $C_L$, $C_R$, $\|U_0\|_{L^\infty}$ and
the flux function $A$, while $C^*$ depends only on the same set as well as T.
  \end{itemize}
 \end{thm}

\begin{rem}
Note that our estimates do not depend on any local  $BV$ norms of
$U_0$. The assumption $U_0\in BV_{loc}$
ensures that  $\frac{d}{dx}U_0$ is a Radon
measure. Hence,  $(\frac{d}{dx}U_0)_+$ is also a Radon measure, and the
condition $(\frac{d}{dx}U_0)_+\in L^2$ makes sense.\\
\end{rem}

This result shows a rate of convergence  slightly worse than $\varepsilon$ (to the
$\log$), for the inviscid limit to a shock, measured via  the $L^2$ norm (squared).
In the case of the limit to a regular solution of the inviscid case,
the rate of convergence is $\sqrt{\varepsilon}$ (see \cite{vass:handbook}, for instance).
We also refer to
\cite{Kuznetsov}, Goodman and Xin \cite{MR1188982},
Bressan, Liu and Yang \cite{MR1723032},
Lewicka \cite{MR2091511},
 Bressan and Yang \cite{MR2053759},
Christoforou and Trivisa \cite{MR2861664}.

\vskip0.3cm
An easy layer study shows that $\varepsilon$ is the optimal rate for shocks with  special
initial data. Indeed, one can construct an associated steady viscous
 layer (see for example Il{$'$}in and Ole{\u\i}nik \cite{Oleinik}) $S_1$
solution to
\begin{equation}\label{layer}
\begin{cases}
&\displaystyle{A(S_1)-A(C_L)-\sigma(S_1-C_L)=S_1',\qquad x\in\mathbb{R},}\\
&\displaystyle{\lim_{x\to-\infty}S_1=C_L,\qquad \lim_{x\to+\infty}S_1=C_R.}
\end{cases}
\end{equation}
It is easy to show that $S_1((x-\sigma t)/\varepsilon)$ is a solution to (\ref{scalar_cl})
 with initial data
$S_1(x/\varepsilon)$. In this case, the rate of convergence is of order $\varepsilon$ since:
$$
\int_\mathbb{R}|S_1((x-\sigma t)/\varepsilon)-S_0(x-\sigma t)|^2\,dx\leq \varepsilon \int_{\mathbb{R}}|S_1(x)-S_0(x)|^2\,dx=C\varepsilon.
$$
This layer study can be extended to the case of small initial perturbation where:
$$
\int_\mathbb{R}|U_0(x)-S_0(x)|^p\,dx\leq C\varepsilon,
$$
for a $1\leq p<\infty$.
In this case, for a  solution $U$ to \eqref{scalar_cl}, we can consider
$$
V(t,x)=U(\varepsilon t,\varepsilon x),
$$
and study the asymptotic for large time.
The function $V$ is a solution to the equation
\begin{equation*}
\begin{cases}&\partial_t V+\partial_x A(V) -\partial^2_{xx}V=0,\\
&V(0,x)=U(0,x/\varepsilon). \end{cases}
\end{equation*}
The convergence to $S_1$, up to a (constant) drift,
 in this setting,  has been extensively studied
(see for instance  \cite{Oleinik}, Freist{\"u}hler and Serre \cite{Freis_Serre},
 Kenig and Merle
\cite{Kenig_Merle}).
In this situation of small perturbation of the initial shock, those results show
that the convergence with rate
$\varepsilon$ for the system (\ref{scalar_cl}) is due
to the asymptotic limit in large time  of the layer
function $U(\cdot/\varepsilon)$ to $S_1((\cdot-\sigma t)/\varepsilon)$.

 \vskip0.3cm

 This layer study, however, collapses when
 $$
 \int_\mathbb{R} |U_0(x)-S_0(x)|^2\,dx\gg \varepsilon.
 $$
In this situation, there is too much entropy for the asymptotic limit
of the layer structure to be true. The physical layer may be destroyed.
Theorem \ref{main_thmL2} shows that,
nevertheless, the sharp convergence (up to the $\log$) still holds.
\vskip0.3cm
Taking a limit as $\varepsilon$ goes to $0$ in Theorem  \ref{main_thmL2}, we recover
 the $L^2$ stability  of
shocks (up to a drift) first showed by Leger in \cite{nick}.  Note that the stability
result has to be up to a drift which depends on the solution itself (and may be not unique).
This feature is also true for our result.
The drift cannot be taken constant, as in the case of the layer problem.

\vskip0.3cm

Our result  is based on the relative entropy method first used by Dafermos and DiPerna
 to show $L^2$ stability and uniqueness of Lipschitzian solutions to conservation laws
\cite{Dafermos4,Dafermos1,DiPerna}. They showed, in particular, that if $\overline{U}$
 is a Lipschitzian solution of a suitable conservation law on a lapse of time $[0,T]$,
then for any bounded weak entropic solution $U$ it holds:
\begin{equation}\label{entropy_decrease}
\int_{\mathbb{R}}|U(t)-\overline{U}(t)|^2\,dx  \leq C \int_{\mathbb{R}}|U(0)-\overline{U}(0)|^2\,dx,
\end{equation}
for a constant $C$ depending on $\overline{U}$ and $T$.\\

The relative entropy method is also an  important tool in the study of asymptotic
 limits. The main idea is that convergence holds thanks to the strong stability of the
solutions of the limit equations. Roughly speaking, if we have good consistency of $\varepsilon$
 models, with respect to the limit one, then non linearities are driven by the strong
 stability of the solution of the limit equation.
Applications of the relative entropy method in this context began with the work of Yau
 \cite{Yau} and have been studied by many others. For incompressible limits, see Bardos,
 Golse, Levermore \cite{Bardos_Levermore_Golse1,Bardos_Levermore_Golse2}, Lions and
 Masmoudi \cite{Lions_Masmoudi}, Saint Raymond et al.
\cite{SaintRaymond1,SaintRaymond2,SaintRaymond3,SaintRaymond4}. For  compressible models,
see Tzavaras \cite{Tzavaras_theory} in the context of relaxation
and \cite{BV,BTV,MV} in the context of hydrodynamical limits. However, in all
those cases, the method works as long as the limit solution is Lipschitz.  This is
due to the fact that strong stability as
(\ref{entropy_decrease}) is not true when $\overline{U}$ has a discontinuity. It has
been proven in  \cite{nick, nick_vasseur}, however, that some shocks are  strongly
stable up to a shift (see also related works from Chen and Frid  \cite{Frid2,Frid1}
and Chen, Frid
and Li \cite{Chen1}).  This article is the first extension  of
those results of stability, to the study of  asymptotic limits to a shock.
This is a part of the program initiated in \cite{vass:handbook}.
\vskip0.3cm
The result can be extended to any entropy in the following way.
Fix any strictly convex function $\eta\in C^2$ as an entropy.
We define the associated  relative entropy functional $\eta(\cdot|\cdot)$ as
$$
\eta(x|y):=\eta(x)-\eta(y)-\eta^\prime(y)(x-y).
$$
We then have the following extension.
\begin{thm}\label{main_thm} Consider a strictly convex entropy functional $\eta\in C^2(\mathbb{R})$.
Let $C_L>C_R$  and $U_0\in L^\infty(\mathbb{R})\cap BV_{loc}(\mathbb{R}) $ be such that
$$
(U_0-S_0)\in L^2(\mathbb{R})\qquad\mathrm{ and}\qquad
  (\frac{d}{dx}U_0)_+\in L^2(\mathbb{R}).
 $$
 Then, for every $T>0$,  there exist $\varepsilon_0>0$  and $C^*>0$ such that  the followings
hold true.
 \begin{itemize}
 \item[I.]
For any  $U$ solution to
 \eqref{scalar_cl} with $0<\varepsilon\leq \varepsilon_0$,
there exists a curve $X\in L^\infty(0,T)$ such that $X(0)=0$, and  for any $0<t<T$,
and for any $\alpha$ verifying $\varepsilon\leq \alpha\leq\varepsilon_0$, we have:
\begin{equation}\begin{split}\label{main_estiamte}
 &\int_{\{|x-X(t)|\geq C^*\alpha\}}\eta(U(t,x)|S(t,x))\,dx\leq
\int_\mathbb{R}\eta(U_0(x)|S_0(x))\,dx
 +C^* e^{-\alpha/\varepsilon},
\end{split}\end{equation}
 where $S(t,x):=S_0(x-X(t))$, and $S_0$ is defined by (\ref{def_S}).

\item[II.] Moreover, this curve satisfies
\begin{equation}\label{estimate_curve_lishitz}
|\dot{X}(t)|\leq C^* \quad\mbox{ and}
\end{equation}
\begin{equation}\label{estimate_curve_w.r.t._L^2_differ}
|X(t)-\sigma t|^2\leq C^* t^{2/3}
\left(\int_\mathbb{R} \eta(U_0(x)|S_0(x))\,dx
+\varepsilon\log(1/\varepsilon)\right).
 \end{equation}

  \item[III.] The constant $\varepsilon_0$ depends only on
$\|(\frac{d}{dx}U_0)_+\|_{L^2}$, $C_L$, $C_R$, $\|U_0\|_{L^\infty}$,
the flux function $A$ and the entropy functional $\eta$, while $C^*$ depends only on
the same set as well as T.
  \end{itemize}
 \end{thm}
Theorem \ref{main_thmL2} is a direct application of Theorem \ref{main_thm} with
$\eta(x):=x^2$, and $\alpha=\varepsilon\log(1/\varepsilon)$. Indeed, in this case we have
$\eta(x|y)=(x-y)^2$, and
$$
\int_{\{|x-X(t)|\leq C^*\alpha\}}\eta(U(t,x)|S(t,x))\,dx\leq
C |\{|x-X(t)|\leq C^*\alpha\}|\leq C C^* \alpha.
$$
\vskip0.3cm
For the rest of the paper, we will assume that the
initial value $U_0$ lies not only $BV_{loc}$ but also $C^1$. It allows us to work with
smooth solutions $U\in C^1([0,T]\times \mathbb{R})$.  The general
 $BV_{loc}$ case can be obtained by a density argument.
\vskip0.3cm
The idea of the proof is to study the evolution of the relative entropy of the solution with respect to the shock, outside of a small region centered at $X(t)$ (this small region corresponds to the layer localization):
\begin{equation}\label{def_quantity}
\int_{-\infty}^{X(t)-\delta\varepsilon}\eta(U(t,x)|C_L)\,dx+\int_{X(t)+\delta\varepsilon}^{\infty}\eta(U(t,x)|C_R)\,dx.
\end{equation}
The change in time involves two effects. One is due to the hyperbolic part of the equation, and the second involves the parabolic part (or order $\varepsilon$).
In \cite{nick}, it was shown that,
for the hyperbolic case $\varepsilon=0$, with zero layer width $\delta=0$, this quantity is non-increasing when we choose wisely the drift $X(t)$. When considering  the viscous term, the layer with width $\delta\varepsilon$ is introduced to mollify its effect. The idea is then that the stability induced by the hyperbolic part is enough to compensate the parabolic one, provided that we consider a layer fat enough. For technical considerations, we will fix $\delta=\log(1/\varepsilon)$. The drift $X(t)$ is still chosen with respect to the hyperbolic part of the equation in a similar way
 as in \cite{nick}. The stability effect is preserved, despite the non zero layer width, thanks to a monotonicity property
induced in the layer  by the additional assumption $(\frac{d}{dx}U_0)_+\in L^2(\mathbb{R})$.

\section{Evolution of the relative entropy}

For $\delta>0$, we consider a Lipschitz nondecreasing
function $\phi_\delta$ to localize the layer, verifying
\begin{equation*}
\phi_\delta(x)=\begin{cases}&0\qquad \mathrm{if} \ \ x\leq 0,\\
&1\qquad \mathrm{if} \ \ x\geq \delta.
\end{cases}\end{equation*}
To get the optimal result, we will later fix a special function (see (\ref{pd})). \\

For any fixed $\delta>0$ and $X\in C^1([0,T])$, we are interesting in the evolution  of
\begin{equation}\label{def_H}
H(t)=\int_{-\infty}^{\infty}\phi_\delta^2(|x-X(t)|/\varepsilon)\eta(U(t,x)|S(t,x))\,dx,
\end{equation}
where the function $S$ depends on $X$ as in the statement of the theorems. A special value of $\delta$ (depending on $\varepsilon$),  and of the function $X$ will be chosen later.
Note that $H(t)$ is controlling the quantity (\ref{def_quantity}).
From now on, we drop the subindex $\delta$ for $\phi_\delta$ for a simple presentation.
\vskip0.3cm
Let us denote $F(\cdot,\cdot)$ the flux of the relative entropy $\eta(\cdot|\cdot)$  defined by
\begin{equation}\begin{split}
  &F(x,y):=G(x)-G(y)-\eta^\prime(y)(A(x)-A(y)).
 \end{split}\end{equation}
Note that \begin{eqnarray*}
 &&\partial_U\eta(U|C)=\eta^\prime(U)-\eta^\prime(C),\\
 &&\partial_U F(U,C)=G'(U)-\eta'(C)A'(U)=(\eta^\prime(U)-\eta^\prime(C))A'(U).
 \end{eqnarray*}
 So, for any solution  $U$ of \eqref{scalar_cl} and any constant $C$, multiplying \eqref{scalar_cl} by $\eta^\prime(U)-\eta^\prime(C)$, we get
\begin{equation}\label{relative_entropy_scalar_cl}
\partial_t\eta(U|C)+\partial_xF(U,C)=\varepsilon(\eta^\prime(U)
-\eta^\prime(C))\partial_{xx}^2
\end{equation}
We have the following lemma.
\begin{lem}\label{baby_lemma}
The function $H$, defined in (\ref{def_H}), satisfies the following on $(0,T)$
\begin{equation*}
\begin{split}
&\qquad\qquad\qquad \frac{dH}{dt}(t)=\\
&{\int_{X(t)-\delta\varepsilon}^{X(t)}
\frac{2}{\varepsilon}
\phi\Big(\frac{-x+X(t)}{\varepsilon}\Big)
\phi^\prime\Big(\frac{-x+X(t)}{\varepsilon}\Big)
\Big[\dot{X}(t)\eta(U(t,x)|C_L)-F(U(t,x),C_L)
\Big]dx}\\
& + {\varepsilon \int_{-\infty}^{X(t)}
\Big[\phi\Big(\frac{-x+X(t)}{\varepsilon}\Big)\Big]^2
\partial^2_{xx}U(t,x) (\eta^\prime(U(t,x))-\eta^\prime(C_L))dx}\\
\end{split}
\end{equation*}
\begin{equation*}
\begin{split}
&{-\int_{X(t)}^{X(t)+\delta\varepsilon}
\frac{2}{\varepsilon}
\phi\Big(\frac{x-X(t)}{\varepsilon}\Big)
\phi^\prime\Big(\frac{x-X(t)}{\varepsilon}\Big)
\Big[\dot{X}(t)\eta(U(t,x)|C_R)-F(U(t,x),C_R)
\Big]dx}\\
& +{\varepsilon \int_{X(t)}^\infty
\Big[\phi\Big(\frac{x-X(t)}{\varepsilon}\Big)\Big]^2
\partial^2_{xx}U(t,x) (\eta^\prime(U(t,x))-\eta^\prime(C_R))dx}\\
&:={(L)_{\mbox{\it{Hyp}}}}+{(L)_{\mbox{\it{Dif}}}}+{(R)_{\mbox{\it{Hyp}}}}+{(R)_{\mbox{\it{Dif}}}}.
\end{split}
\end{equation*}
\end{lem}
\begin{proof}
First we split the above integral into the two parts:
  \begin{equation*}\begin{split}
  &H(t)=\int_{-\infty}^{\infty}
\Big[\phi\Big(\frac{|x-X(t)|}{\varepsilon}\Big)\Big]^2\eta(U|S)dx\\
&=\int_{-\infty}^{\infty}\Big(
\Big[\phi\Big(\frac{-x+X(t)}{\varepsilon}\Big)\Big]^2+\Big[\phi\Big(\frac{x-X(t)}{\varepsilon}\Big)\Big]^2\Big)\eta(U|S)dx\\
&={\int_{-\infty}^{\infty}
\Big[\phi\Big(\frac{-x+X(t)}{\varepsilon}\Big)\Big]^2\eta(U|C_L)dx}
+{\int^{\infty}_{-\infty}
\Big[\phi\Big(\frac{x-X(t)}{\varepsilon}\Big)\Big]^2\eta(U|C_R)dx}\\
&:=H^L+H^R.
\end{split}\end{equation*}

To compute $\frac{d}{dt}(H^{L})$,  we put $C=C_L$ in \eqref{relative_entropy_scalar_cl}, multiply by
$\Big[\phi\Big(\frac{-x+X(t)}{\varepsilon}\Big)\Big]^2$,
 and integrate in $x$. Then   we have
\begin{equation*}\begin{split}
\frac{d}{dt}(H^{L})
&=\int_{-\infty}^{\infty}
\partial_t
\Big(\Big[\phi\Big(\frac{-x+X(t)}{\varepsilon}\Big)\Big]^2\Big) \eta(U|C_L)dx\\
&\quad +\int_{-\infty}^{\infty}
\partial_x\Big(\Big[\phi\Big(\frac{-x+X(t)}{\varepsilon}\Big)\Big]^2\Big)
F(U,C_L)dx\\
&\quad +\varepsilon\int_{-\infty}^{\infty}
\Big[\phi\Big(\frac{-x+X(t)}{\varepsilon}\Big)\Big]^2
\partial^2_{xx}U (\eta^\prime(U)-\eta^\prime(C_L))dx
\end{split}\end{equation*}
\begin{equation*}\begin{split}
&={\int_{X(t)-\delta\varepsilon}^{X(t)}
(\frac{2}{\varepsilon})
\cdot
\phi\Big(\frac{-x+X(t)}{\varepsilon}\Big)
\phi^\prime\Big(\frac{-x+X(t)}{\varepsilon}\Big)
\Big[\dot{X}(t)\eta(U|C_L)-F(U,C_L)
\Big]dx}\\
&\quad +{\varepsilon\cdot\int_{-\infty}^{X(t)}
\Big[\phi\Big(\frac{-x+X(t)}{\varepsilon}\Big)\Big]^2
\partial^2_{xx}U (\eta^\prime(U)-\eta^\prime(C_L))dx}
=(L)_{\mbox{\it{Hyp}}}+(L)_{\mbox{\it{Dif}}}.
\end{split}\end{equation*}
We get the result for $\frac{d}{dt}(H^{R})={(R)_{\mbox{\it{Hyp}}}}+{(R)_{\mbox{\it{Dif}}}} $
in the same way.
\end{proof}

\section{Study of the contribution of the hyperbolic structure}

In this section, we show that by choosing a special drift function $X(\cdot)$, the hyperbolic effects become
nonpositive. This will be used in section \ref{section_parabolic} to control
the parabolic effects.\\

Following \cite{nick},  we define the normalized relative entropy flux
$f(\cdot,\cdot)$ by
$$f(x,y):=\frac{F(x,y)}{\eta(x|y)}.$$ We have the following properties.
 \begin{lem}\label{lem_basic_propperties} For any $L>0$, there exists  a constant
 $\Lambda>0$,  such that for any $x,y$ with $|x|,|y|\leq L$, we have
 \begin{equation}\begin{split}\label{properties_eta}.\\
&1/\Lambda \leq\eta^{\prime\prime}(x)\leq \Lambda,\\
&\frac{1}{2\Lambda}(x-y)^2
\leq\eta(x|y)\leq \frac{1}{2}{\Lambda}(x-y)^2,\\
   &|F(x,y)|\leq \Lambda(x-y)^2,\\
 &0\leq (\partial_x f)(x,y)\leq \Lambda, \\
&1/\Lambda\leq(\partial_y f)(x,y).
 \end{split}\end{equation}
 \end{lem}
The proof of this lemma can be found in   \cite{nick}.
\vskip0.3cm
We now define the shift function $X$. It is the solution of the following O.D.E.
\begin{equation}
\begin{cases}\label{ode}
&\dot{X}(t)=f\Big(U(t,X(t)),\frac{C_L+C_R}{2}\Big)\\
&X(0)=0\end{cases}. \end{equation}
Note that for any $\varepsilon>0$, $U\in C^1([0,T]\times\mathbb{R})$ (since $U_0\in C^1(\mathbb{R})$). The existence and uniqueness of $X$ comes from the Cauchy-Lipschitz theorem.
\vskip0.3cm
First, $X$ is Lipschitz, since we have from Lemma \ref{lem_basic_propperties}
\begin{equation}\begin{split}
|\dot{X}(t)|&
\leq\frac{\Big|F\Big(U(t,X(t)),\frac{C_L+C_R}{2}\Big)\Big|}
{\eta\Big(U(t,X(t))\Big |\frac{C_L+C_R}{2}\Big)}
\leq
2\Lambda^2
\end{split}\end{equation} where we used the fact
$\|U(t)\|_{L^\infty}\leq\|U_0\|_{L^\infty}$ for $t>0$. It proves \eqref{estimate_curve_lishitz}.

\vskip0.3cm

Note that from the definition of $X$, if the the velocity $U$ were constant in the layer (that is $U(t,x)\sim U(t,X(t))$ for $-\delta\varepsilon\leq x-X(t)\leq \delta\varepsilon$), then, from the
 last property of Lemma \ref{lem_basic_propperties}, we would have directly that
$$
(L)_{\mbox{\it{Hyp}}}+(R)_{\mbox{\it{Hyp}}}\leq -\frac{C_L-C_R}{\Lambda} (\eta(U(t,X(t))|C_L)+\eta(U(t,X(t))|C_R)).
$$
However, this is too much to hope, since the layer characterize
the region where the function $U(t,\cdot)$ is expected
to drop from $C_L$ to $C_R$. We still can show that the
hyperbolic terms are negative, provided that the behavior
of $U$ in the layer is not too much oscillatory
(the values can drop, but not much bounce back). This last property of $U$ is proved in the following lemma which can be seen as a weak version of the Ole{\u\i}nik's principle
.
\begin{lem}\label{lem_decreasing_L_2_positive_derivative}
$\|(\partial_xU(t))_+\|_{L^2(\mathbb{R})}\leq\|(\frac{d}{dx}U_0)_+\|_{L^2(\mathbb{R})}$ for any $t>0$.
\end{lem}
\begin{proof}
We differentiate \eqref{scalar_cl} w.r.t. $x$, multiply $(\partial_xU)_+$ and
integrate in $x$ to get
\begin{equation*}\begin{split}
&0=\int(\partial_xU)_+ \Big[\partial_t\partial_xU+
A^{\prime\prime}(U) |\partial_xU|^2+A^{\prime}(U) \partial^2_{xx}
U-\varepsilon \partial_{xxx}^3U\Big]dx\\
&=\int\Big[\frac{1}{2}\partial_t([(\partial_xU)_+]^2)
+A^{\prime\prime}(U) (\partial_xU)_+^3\\
&\quad\quad\quad\quad\quad\quad\quad
+A^{\prime}(U) \partial_{x}\Big(\frac{[(\partial_xU)_+]^2}{2}\Big)+\varepsilon
|\partial_x((\partial_xU)_+)|^2
\Big]dx.
 \end{split}\end{equation*} Then, we use the integration by parts to get
 \begin{equation*}\begin{split}
&=\int\Big[\frac{1}{2}\partial_t([(\partial_xU)_+]^2)
+\frac{1}{2}A^{\prime\prime}(U) (\partial_xU)_+^3
+\varepsilon
|\partial_x((\partial_xU)_+)|^2
\Big]dx\\
&\geq\frac{1}{2}\frac{d}{dt}\int[(\partial_xU)_+]^2
dx.\\
 \end{split}\end{equation*}

\end{proof}

We now prove the main proposition of this section.

\begin{prop}\label{prop_Hyp}
Let $(L)_{\mbox{\it{Hyp}}}$ and $(R)_{\mbox{\it{Hyp}}}$ be such as in Lemma \ref{baby_lemma}.
There exists a constant $\theta>0$ such that, for any $\varepsilon, \delta$ satisfying
$$
\varepsilon\delta\leq {\theta},
$$
we have
\begin{eqnarray*}
&&\qquad\qquad (L)_{\mbox{\it{Hyp}}}+(R)_{\mbox{\it{Hyp}}}\\
&&\leq -\frac{\theta}{\varepsilon}\int_{X(t)-\delta\varepsilon}^{X(t)+\delta\varepsilon}\phi\left(\frac{|x-X(t)|}{\varepsilon}\right)\phi'\left(\frac{|x-X(t)|}{\varepsilon}\right)(U(t,x)-S(t,x))^2\,dx.
\end{eqnarray*}

\end{prop}
\begin{proof}
We use the definition of $X(t)$ to get
\begin{equation*}\begin{split}
&(L)_{\mbox{\it{Hyp}}}
=\int_{X(t)-\delta\varepsilon}^{X(t)}
(\frac{2}{\varepsilon})
\cdot
\phi\Big(\frac{-x+X(t)}{\varepsilon}\Big)\cdot
\phi^\prime\Big(\frac{-x+X(t)}{\varepsilon}\Big)
\cdot\eta(U|C_L)
\cdot h(t,x)dx
\end{split}\end{equation*} where $h(t,x):=\Big[f\Big(U(t,X(t)),\frac{C_L+C_R}{2}\Big)-f(U(t,x),C_L)
\Big]$.\\

In order to make the function $h(t,x)$ strictly negative over the
domain of the above integral, we use the condition
$(\frac{d}{dx}U_0)_+\in{L^2(\mathbb{R})}$. Indeed, we
observe that, for any $x\in[X(t)-\delta\varepsilon,X(t)]$,
\begin{equation*}\begin{split}
U(t,X(t))-U(t,x)&=\int_x^{X(t)}(\partial_xU)(t,y)dy
\leq\int_x^{X(t)}(\partial_xU)_+(t,y)dy\\
&\leq\|(\partial_xU(t))_+\|_{L^2(\mathbb{R})}\sqrt{|X(t)-x|}
\leq\|(\frac{d}{dx}U_0)_+\|_{L^2(\mathbb{R})}\sqrt{\delta\varepsilon},
\end{split}\end{equation*} where we used that $\|(\partial_xU(t))_+\|_{L^2}$ is not increasing (see Lemma \ref{lem_decreasing_L_2_positive_derivative}).

We can rewrite the function $h$ as
\begin{equation*}\begin{split}
h(t,x)&=f\Big(U(t,X(t)),\frac{C_L+C_R}{2}\Big)
-f\Big(U(t,x),\frac{C_L+C_R}{2}\Big)\\
&\quad\quad\quad\quad\quad
+f\Big(U(t,x),\frac{C_L+C_R}{2}\Big)
-f(U(t,x),C_L).\\
\end{split}\end{equation*}
 Since $f$ is increasing with respect to the first variable, we have
\begin{equation*}\begin{split}
&h(t,x)\leq f\Big(U(t,x)+\|(\frac{d}{dx}U_0)_+\|_{L^2(\mathbb{R})}\sqrt{\delta\varepsilon},\frac{C_L+C_R}{2}\Big)
-f\Big(U(t,x),\frac{C_L+C_R}{2}\Big)\\
&\quad\quad\quad\quad\quad+f\Big(U(t,x),\frac{C_L+C_R}{2}\Big)
-f(U(t,x),C_L).\\
\end{split}\end{equation*}
Then, thanks to Lemma \ref{lem_basic_propperties},
  we get
\begin{equation*}\begin{split}
&h(t,x)\leq \Lambda\|(\frac{d}{dx}U_0)_+\|_{L^2(\mathbb{R})}
\sqrt{\delta\varepsilon}
-\frac{{C_L-C_R}}{2\Lambda}\leq -\theta<0\\
\end{split}\end{equation*}
for $\sqrt{\delta\varepsilon}$ small enough.\\

\noindent Since $\phi(\cdot),\phi^\prime(\cdot)$ and $\eta(\cdot|\cdot)\geq0$,  we get
 \begin{equation*}\begin{split}
&(L)_{\mbox{\it{Hyp}}}\leq -\theta\int_{X(t)-\delta\varepsilon}^{X(t) }
 \frac{2}{\varepsilon}
\phi\Big(\frac{-x+X(t)}{\varepsilon}\Big)
\phi^\prime\Big(\frac{-x+X(t)}{\varepsilon}\Big)
\eta(U|C_L)  dx.\\
\end{split}\end{equation*}
Then, from Lemma \ref{lem_basic_propperties}, we have (changing the constant $\theta$ if necessary)
 \begin{equation*}\begin{split}
&(L)_{\mbox{\it{Hyp}}}\leq-\theta \int_{X(t)-\delta\varepsilon}^{X(t) }
( \frac{2}{\varepsilon})
\phi\Big(\frac{-x+X(t)}{\varepsilon}\Big)
\phi^\prime\Big(\frac{-x+X(t)}{\varepsilon}\Big)
(U-C_L)^2
 dx.\\
\end{split}\end{equation*}

In a similar way, we obtain the following estimate on $(II)_{Hyp}$.
 \begin{equation*}\begin{split}
&(R)_{\mbox{\it{Hyp}}}\leq-\theta \int_{X(t)}^{X(t)+\delta\varepsilon }
( \frac{2}{\varepsilon})
\phi\Big(\frac{x-X(t)}{\varepsilon}\Big)
\phi^\prime\Big(\frac{x-X(t)}{\varepsilon}\Big)
(U-C_R)^2
 dx.\\
\end{split}\end{equation*}
Combining the two last inequalities gives the desired result.
\end{proof}

\section{Control of the parabolic terms}\label{section_parabolic}

For any $\delta\geq 1/\theta$, we now fix the function $\phi_\delta$  in the following explicit way.
\begin{equation}\label{pd}
\phi_\delta(x)=
\begin{cases}
& \theta e^{1-\theta\delta}x,\qquad \mathrm{for} \ \ x\in [0,1/\theta),\\
&e^{\theta(x-\delta)},\qquad \mathrm{for} \ \ x\in [1/\theta, \delta].
\end{cases}
\end{equation}
It can be proved that this functional minimize
$$
\int_0^\delta (\phi'(x))^2 \chi_{\{\phi'>\theta\phi\}}\,dx
$$
among all increasing Lipschitz functions with $\phi(0)=0$, $\phi(\delta)=1$ and
$\delta\geq 1/\theta$.
We will only use the straightforward computation:
\begin{equation}\label{computation_pd}
\int_0^\delta (\phi_\delta'(x))^2 \chi_{\{\phi_\delta'>\theta\phi_\delta\}}\,dx=C_\theta \cdot e^{-2\theta \delta}.
\end{equation}

This section is dedicated to the proof of the following proposition.
\begin{prop}\label{prop_diffusion}
There exists constants $\theta, C>0$ such that for any $\varepsilon,\delta$ verifying
$$
\frac{1}{\theta}\leq \delta\qquad\mbox{and}\qquad \varepsilon\delta\leq {\theta},
$$
we have
$$
\frac{dH(t)}{dt}\leq Ce^{-\theta\delta}.
$$
\end{prop}
\begin{proof}
First, we estimate the term  $(L)_{\mbox{\it{Dif}}}$. Integrating  by parts,  we obtain
\begin{equation*}\begin{split}
(L)_{\mbox{\it{Dif}}}&=\int_{-\infty}^{X(t)}
 2\phi\Big(\frac{-x+X(t)}{\varepsilon}\Big)
\phi^\prime\Big(\frac{-x+X(t)}{\varepsilon}\Big)
\partial_{x}U  (\eta^\prime(U)-\eta^\prime(C_L))
dx\\
&\quad -2\varepsilon \int_{-\infty}^{X(t)}
\Big[\phi\Big(\frac{-x+X(t)}{\varepsilon}\Big)\Big]^2
\eta^{\prime\prime}(U)
|\partial_{x}U|^2dx\\
\end{split}\end{equation*}
 Then, using H\"{o}lder's inequality and  Lemma \ref{lem_basic_propperties}, we get
 \begin{equation*}\begin{split}
&(L)_{\mbox{\it{Dif}}}\leq
\frac{2\varepsilon}{\Lambda} \int_{-\infty}^{X(t)}
\Big[\phi\Big(\frac{-x+X(t)}{\varepsilon}\Big)\Big]^2
|\partial_{x}U|^2dx\\
&\quad+\frac{\Lambda}{8\varepsilon}\int_{\infty}^{X(t)}
\Big[2 \phi^\prime\Big(\frac{-x+X(t)}{\varepsilon}\Big)(\eta^\prime(U)-\eta^\prime(C_L))\Big]^2dx\\
&\quad-\frac{2\varepsilon}{\Lambda} \int_{-\infty}^{X(t)}
\Big[\phi\Big(\frac{-x+X(t)}{\varepsilon}\Big)\Big]^2
|\partial_{x}U|^2dx\\
&\leq\frac{C}{\varepsilon}\int_{X(t)-\delta\varepsilon}^{X(t)}
\Big[ \phi^\prime\Big(\frac{-x+X(t)}{\varepsilon}\Big)\Big]^2
|U-C_L|^2 \,dx.
\end{split}\end{equation*}
In the same way, we obtain the following estimate for $(R)_{\mbox{\it{Dif}}}$.
$$
(R)_{\mbox{\it{Dif}}}\leq \frac{C}{\varepsilon}\int_{X(t)}^{X(t)+\delta\varepsilon}
\Big[ \phi^\prime\Big(\frac{x-X(t)}{\varepsilon}\Big)\Big]^2
|U-C_R|^2 \,dx.
$$
Combining the two last inequalities, we find
\begin{equation}\label{estimate_diffusion}
(L)_{\mbox{\it{Dif}}}+(R)_{\mbox{\it{Dif}}}\leq \frac{C}{\varepsilon}\int_{X(t)-\delta\varepsilon}^{X(t)+\delta\varepsilon}
\Big[ \phi^\prime\Big(\frac{|x-X(t)|}{\varepsilon}\Big)\Big]^2
|U(t,x)-S(t,x)|^2 \,dx.
\end{equation}

\vskip0.3cm
\noindent Using Lemma \ref{baby_lemma}, Proposition \ref{prop_Hyp}, and (\ref{estimate_diffusion}), we find
\begin{equation}\begin{split}\label{eq_keep_charac}
&\frac{dH(t)}{dt}\leq
 \frac{1}{\varepsilon}\int_{X(t)-\delta\varepsilon}^{X(t)+\delta\varepsilon}
\Big[\phi^\prime (C\phi'-\theta\phi)\Big]\left(\frac{|x-X(t)|}{\varepsilon}\right)|U(t,x)-S(t,x)|^2\,dx.
 \end{split}\end{equation}
 Using that $U-S$ is a bounded function, and doing the change of
variables $z=(x-X(t))/\varepsilon$, we find:
\begin{equation*}\begin{split}
 \frac{dH(t)}{dt} &
 \leq  \frac{C}{\varepsilon}\int_{X(t)-\delta\varepsilon}^{X(t)+\delta\varepsilon}
\Big[(\phi^\prime)^2 \chi_{\{C\phi'-\theta\phi>0\}}\Big] \left(\frac{|x-X(t)|}{\varepsilon}\right)|U(t,x)-S(t,x)|^2\,dx\\
&\leq \frac{C\|U(t)-S(t)\|^2_{L^\infty}}{\varepsilon}\int_{X(t)-\delta\varepsilon}^{X(t)+\delta\varepsilon}
\Big[(\phi^\prime)^2 \chi_{\{C\phi'-\theta\phi>0\}}\Big] \left(\frac{|x-X(t)|}{\varepsilon}\right) \,dx\\
&\leq C\int_{0}^{\delta}
(\phi^\prime)^2(z) \chi_{\{C\phi'-\theta\phi>0\}}(z)\,dz.
 \end{split}\end{equation*}
 Changing the constant $\theta$ if needed, and using (\ref{computation_pd}), gives the desired result.
 \end{proof}

 \section{Proof of Theorem \ref{main_thm}}
 This section is devoted to the proof of Theorem \ref{main_thm}.
  Integrating the estimate of Proposition \ref{prop_diffusion} between $0$
  and $t\in(0,T)$ gives the result of (I).
  Indeed,  for any $\varepsilon,\delta$ with
$
\frac{1}{\theta}\leq \delta\mbox{ and } \varepsilon\delta\leq {\theta},
$
 where
  $\theta$ is the constant from Proposition \ref{prop_diffusion},
  we have
  \begin{equation*}\begin{split}
  \int_{\{|x-X(t)|\geq \delta\varepsilon\}}\eta(U(t,x)|S(t,x))\,dx
  &\leq H(t)
  \leq H(0)  +\int_0^t \frac{d}{dt}H(s)\,ds\\
  &\leq
\int_\mathbb{R}\eta(U_0|S_0)\,dx
 +CT e^{-\theta\delta}
 \end{split}\end{equation*}
  By taking $\varepsilon_0:=\theta^2$,
     we have  for any $\varepsilon\leq\alpha\leq\varepsilon_0$,
  \begin{equation*}\begin{split}
  \int_{\{|x-X(t)|\geq \alpha/\theta\}}\eta(U(t,x)|S(t,x))\,dx
   &\leq
\int_\mathbb{R}\eta(U_0|S_0)\,dx
 +CT e^{-\alpha/\varepsilon}.
 \end{split}\end{equation*} It proves our main estimate \eqref{main_estiamte}
 by taking $C^*$ large enough.\\%=\max\{1/\theta, CT\}$.\\

Observe that
\begin{equation*}\begin{split}
 \int_{\mathbb{R}} \eta(U|S)\,dx&=
 \int_{\{|x-X(t)|\geq C^*\alpha\}}\eta(U|S)\,dx
 +\int_{\{|x-X(t)|< C^*\alpha\}}\eta(U|S)\,dx\\
\end{split}\end{equation*}
and the second term is bounded by $CC^*\alpha$.
  Thus, by taking $\alpha=\varepsilon\log(1/\varepsilon)$, we obtain for any $t\in(0,T)$,
 \begin{equation}\label{estimate_temp}
  \int_{\mathbb{R}} \eta(U|S)\,dx\leq \int_{\mathbb{R}} \eta(U_0|S_0)\,dx+C^*\varepsilon\log(1/\varepsilon)
 \end{equation} for any $\varepsilon\leq\varepsilon_0$
(changing $\varepsilon_0$ and $C^*$ if needed).

 \vskip0.3cm
  It only remains to prove \eqref{estimate_curve_w.r.t._L^2_differ}.
We define first
$\psi(x):=\begin{cases}
&0 \mbox{ if } |x|>2,\\
&1 \mbox{ if } |x|\leq1\\
&2-|x| \mbox{ if } 1<|x|\leq2
\end{cases}$. Let $s\in(0,t)$ and $R>0$. We multiply $\Psi_R(s,x):=\psi(\frac{x-X(s)}{R})$
to the  equation  \eqref{scalar_cl} and integrate in $x$ to get
\begin{equation*}\begin{split}
0=&-\frac{d}{ds}\int \Psi_R\cdot U dx+\int\partial_x(\Psi_R) A(U) dx
+\int\partial_t(\Psi_R) U dx+\varepsilon\int\Psi_R
\cdot\partial^2_{xx}Udx\\
&=-
\underbrace{\frac{d}{ds}\int \psi(\frac{x-X(s)}{R})\cdot U(s,x) dx}_{(I)}\\
&\quad\quad+\underbrace{\frac{1}{R}\int\psi^\prime(\frac{x-X(s)}{R})\cdot \Big(
A(U(s,x))-\dot{X}(s)U(s,x) \Big)dx}_{(II)}\\
&\quad\quad-\underbrace{\varepsilon\frac{1}{R}\int\psi^\prime(\frac{x-X(s)}{R})
\cdot\partial_{x}U(s,x)dx}_{(III)}
 \end{split}\end{equation*}
By using the above observation,
 we have
\begin{equation*}\begin{split}
&(\sigma-\dot{X}(s))=\frac{1}{C_L-C_R}\Big(A(C_L)-A(C_R)-(C_L-C_R)\dot{X}(s)\Big)\\
&=\frac{1}{C_L-C_R}\Big(A(C_L)-A(C_R)-(C_L-C_R)\dot{X}(s)
-(II)+(I)+(III)\Big).\\
 \end{split}\end{equation*} Then we integrate the above equation in time on $[0,t]$
to get:
\begin{equation}\begin{split}\label{estimate_sum}
|\sigma t-{X}(t)|&\leq
C\Big(t\cdot\max_{s\in(0,t)}\underbrace{\Big|A(C_L)-A(C_R)-(C_L-C_R)\dot{X}(s)
-(II)\Big|}_{(II^\prime)}\\
&\quad+\Big|\int_0^t(I)ds\Big|+t\cdot\max_{s\in(0,t)}\Big|(III)\Big|\Big).\\
 \end{split}\end{equation}

We observe
\begin{equation*}\begin{split}
{(II^\prime)}\leq&\underbrace{\Big|A(C_L)-A(C_R)
-{\frac{1}{R}\int\psi^\prime(\frac{x-X(s)}{R})\cdot
A(U)dx}\Big|}_{(II^\prime_1)}\\
&+\underbrace{\Big|-(C_L-C_R)\dot{X}(s)+\frac{1}{R}\int\psi^\prime(\frac{x-X(s)}{R})\cdot \Big(
\dot{X}(s)U(s,x) \Big)dx\Big|}_{(II^\prime_2)}.\\
 \end{split}\end{equation*}
 For  $(II^\prime_1)$, we compute
\begin{equation*}\begin{split}
(II^\prime_1)
&=\Big|A(C_L)
-{\frac{1}{R}\int_{-2R+X(s)}^{-R+X(s)}
A(U)dx}
-A(C_R)
+{\frac{1}{R}\int_{R+X(s)}^{2R+X(s)}
A(U)dx}\Big|\\
&\leq
{\frac{1}{R}\Big[\int_{-2R+X(s)}^{-R+X(s)}
|A(C_L)-A(U)|dx}
+{\int_{R+X(s)}^{2R+X(s)}
|A(U)-A(C_R)|dx}\Big].
\end{split}\end{equation*} We use  $|A(y)-A(z)|\leq
C
|y-z|$ for $|y|,|z|\leq {M_1}$ to get
\begin{equation*}\begin{split}
\leq
{\frac{
C
}{R}\int_{-2R+X(s)}^{2R+X(s)}
|U-S|dx}.
\end{split}\end{equation*} We use the H\"{o}lder's inequality
and Lemma
 \ref{lem_basic_propperties} to get
\begin{equation*}\begin{split}
(II^\prime_1)^2&
\leq\frac{C}{{R}}\cdot
{\int_{\mathbb{R}}\eta(U(s)|S(s))dx}.
 \end{split}\end{equation*}
Likewise, for the second term $(II^\prime_2)$, we have
\begin{equation*}\begin{split}
(II^\prime_2)&=|\dot{X}(s)|\cdot\Big|-(C_L-C_R)
+\frac{1}{R}\int\psi^\prime(\frac{x-X(s)}{R})\cdot
U(s,x) dx\Big|\\
&\leq
{\frac{
C
}{R}\int_{-2R+X(s)}^{2R+X(s)}
|U-S|dx}
\leq\frac{C}{\sqrt{R}}\cdot
\|U(s)-S(s)\|_{L^2(\mathbb{R})}
 \end{split}\end{equation*} where we used $|\dot{X}(s)|\leq
C$.
Thus we get
\begin{equation}\begin{split}\label{estimate_(II)}
(II^\prime)^2\leq\frac{C}{{R}}\cdot
{\int_{\mathbb{R}}\eta(U(s)|S(s))dx}.
 \end{split}\end{equation}
On the other hand, we compute
\begin{equation*}\begin{split}
&\Big|\int_0^t(I)ds\Big|
=\Big|\int \psi(\frac{x-X(t)}{R})\cdot U(t,x) dx-
\int \psi(\frac{x}{R})\cdot U_0(x) dx\Big|\\
&=\Big|\int \psi(\frac{x-X(t)}{R})\cdot \Big(U(t,x)-S(t,x)\Big) dx+
\int \psi(\frac{x-X(t)}{R})\cdot S(t,x) dx\\
&-\int \psi(\frac{x}{R})\cdot S_0(x) dx
-\int \psi(\frac{x}{R})\cdot \Big(U_0(x)-S_0(x) \Big)dx\Big|.
 \end{split}\end{equation*} Note that
$\int \psi(\frac{x-X(t)}{R})\cdot S(t,x) dx=\int \psi(\frac{x}{R})\cdot S_0(x) dx.$ Thus, we have
\begin{equation*}\begin{split}
&\leq\Big|\int \psi(\frac{x-X(t)}{R})\cdot \Big(U(t,x)-S(t,x)\Big) dx\Big|
+\Big|\int \psi(\frac{x}{R})\cdot \Big(U_0(x)-S_0(x) \Big)dx\Big|.
 \end{split}\end{equation*} We use the H\"{o}lder's inequality and Lemma
 \ref{lem_basic_propperties}
 to get
\begin{equation}\begin{split}\label{estimate_(I)}
\Big|\int_0^t(I)ds\Big|^2&\leq C{R}\Big({\int_{\mathbb{R}}\eta(U(t)|S(t))dx}+
{\int_{\mathbb{R}}\eta(U_0|S_0)dx}\Big).
 \end{split}\end{equation}
Also, we have
\begin{equation}\begin{split}\label{estimate_(III)}
&\Big|(III)\Big|=
\frac{\varepsilon}{R}\Big|\int\psi^\prime(\frac{x-X(s)}{R})
\cdot\partial_{x}U(s,x)dx\Big|\\
&=
\frac{\varepsilon}{R}\Big|\int_{-2R+X(s)}^{-R+X(s)}
\partial_{x}U(s,x)dx-
\int_{R+X(s)}^{2R+X(s)}
\partial_{x}U(s,x)dx\Big|\\
&\leq \frac{\varepsilon}{R}\cdot 4\cdot\|U(s)\|_{L^\infty}
\leq \frac{C\cdot\varepsilon}{R}.
 \end{split}\end{equation}

Finally, by using \eqref{estimate_temp}, we combine
\eqref{estimate_(II)}, \eqref{estimate_(I)} and  \eqref{estimate_(III)} with
\eqref{estimate_sum}
to get, for any  $t\in(0,T)$,
\begin{equation*}\begin{split}
|\sigma t-{X}(t)|^2&\leq C\Big(
\frac{ t^2}{{R}}+
{R}
\Big)\cdot\Big(
{\int_{\mathbb{R}}\eta(U_0|S_0)dx}+{
\varepsilon\log(\frac{1}{\varepsilon})}\Big)
+\frac{C\cdot\varepsilon^2\cdot t^2}{R^2}.\\
 \end{split}\end{equation*} Since the above estimate holds for any $0<R<\infty$,
 the estimate
 \eqref{estimate_curve_w.r.t._L^2_differ} follows once we take  $R:=t^{2/3}$
(changing $C^*$ if needed).

\begin{ack}
The second author was partially funded by the NSF.
\end{ack}

\bibliographystyle{plain}


\section{Introduction}

The study of transiting exoplanets is rapidly advancing our understanding of planets beyond our own Solar system. Planets' bulk densities are obtained via measurement of the radius and mass using transits and radial velocity measurements; this is the first step in understanding their structure and composition. However, to understand planetary systems more fully and explore their diversity we need spectroscopic measurements of their atmospheres; luckily, transiting planets allow such measurements without requiring the star and planet to be spatially resolved.

Transmission spectroscopy is a measurement of the effective size of the planet as a function of wavelength during primary transit. Due to wavelength dependent opacities in the atmosphere, the planet appears larger at wavelengths where the atmosphere absorbs or scatters light. We can therefore probe for the presence of atomic and molecular species, as well as clouds or hazes \citep{Seager_2000,Brown_2001}. Until recently, transmission spectroscopy has only been feasible using space-based telescopes, which have been tremendously successful, particularly at optical wavelengths \citep[e.g.][]{Charbonneau_2002,Pont_2008,Sing_2008,Sing_2011,Huitson_2012}. Measurements in the NIR have proved more controversial, with different groups reporting conflicting conclusions from the same dataset \citep[e.g.][]{Swain_2008, Gibson_2011}; however, advances in data analysis techniques are starting to resolve the issue \citep[e.g.][]{Gibson_2012,Waldmann_2012}, along with the availability of more stable NIR cameras such as WFC3 \citep[e.g.][]{Berta_2012, Gibson_2012b}.

Most recently, the use of a multi-object spectrograph (MOS) to perform differential spectro-photometry has started to show that ground based observations can also play an important role. This technique was pioneered by \citet{Bean_2010}, using the VLT to produce a transmission spectrum of the super-Earth GJ-1214b. Such observations have benefitted not only from advances in observing strategy, but also the availability of new bright targets with nearby comparison stars. Similar observations have since been made in the near-infrared \citep{Bean_2011} and with long-slit spectrographs \citep{Sing_2012}. Here we report Gemini observations of WASP-29 using this MOS technique.

WASP-29b \citep{Hellier_2010} is a Saturn-size planet with a mass and radius of $0.24\pm0.02\,M_J$ and $0.79\pm0.05R_J$. It orbits a K4 dwarf with a period of $\sim$3.9\,days. Given its relatively low equilibrium temperature of $\sim$980\,K, it does not have a particularly large scale height compared to the hottest and lowest-density hot Jupiters. It was in fact chosen as a test case for the Gemini GMOS instrument, to test its stability before more favourable targets are observed. WASP-29 has a nearby, bright comparison star, making it ideal for differential photometry or spectroscopy.

Despite its recent successes, MOS differential spectroscopy, like other techniques, is expected to suffer from correlated noise originating from poorly understood instrumental effects. Although the differential nature of the measurements is intended to reduce this problem, systematic effects are not yet fully understood for ground based differential photometry, and the extra complexity in spectroscopic work may introduce more systematic effects.
In this paper, we explore the use of several methods to analyse transit light curves to obtain useful measurements in the presence of significant correlated noise. These methods include the wavelet method of \citet{Carter_2009}, and the Gaussian process (GP) model of \citet{Gibson_2012} adapted to model time-correlated noise.

This paper is structured as follows: Sect.~\ref{sect:observations} presents our observations and data reduction procedure, in Sect.~\ref{sect:analysis} we present our transit analysis methods, and finally in Sects.~\ref{sect:results} and \ref{sect:discussion} we present our results and conclusions.

\section{GMOS Observations and Data Reduction}
\label{sect:observations}

A transit of WASP-29 was observed using the 8-m Gemini-South telescope with the Gemini Multi-Object Spectrograph (GMOS) on October 19 2011. Data were taken as part of program GS-2011B-Q-13 (PI. Gibson). GMOS has an imaging field-of-view of 5.5$\times$5.5 arcmin squared, and consists of three 2048$\times$4608 pixel CCDs arranged side by side with small gaps in-between. We used GMOS in multi-object mode to observe the target plus two comparison stars simultaneously and continuously for $\sim$5.1 hours, covering the 2.66 hour transit plus 1.5 hours prior to ingress and 1.0 hours after egress. Conditions were photometric for the duration of the observations.

Observations used the R400 grism + OG515 filter with a central wavelength of 725\,nm. The dispersion is 0.07\,nm\ per (unbinned) pixel, giving wavelength coverage from $\sim$515--940\,nm. Exposure times were 30.5 seconds and 2$\times$2 binning was used. The full frame readout of the GMOS chips in 2$\times$2 binning is $\sim$55 seconds (in the recommended `slow' readout). In order to reduce the overheads we read out only three regions of interest (ROI) on the chip containing the three stellar spectra, resulting in a readout time of $\sim$22 seconds, and therefore a cadence of $\sim$53 seconds. This allowed 348 exposures over the 5.1 hours of observation. To minimise slit losses we created a mask with slits of 30 arcsec length and 10 arcsec width for the three stars, giving seeing limited (therefore variable) resolution ranging from R$\sim$430--860 at 700\,nm.
Calibrations were taken before and after the science observations, and consisted of flat fields and arc lamp exposures. Further arcs were taken with a calibration mask. This was almost identical to the science mask but with narrow slits of 1 arcsec to enable more precise wavelength calibration.

The data were reduced using the standard GMOS pipeline contained in the Gemini {\sc IRAF}\footnote{{\sc IRAF} is distributed by the National Optical Astronomy Observatory, which is operated by the Association of Universities for Research in Astronomy (AURA) under cooperative agreement with the National Science Foundation}/{\sc PyRAF}\footnote{{\sc PyRAF} is a product of the Space Telescope Science Institute, which is operated by AURA for NASA} package. First the ROI images were processed to be in the standard GMOS format. Basic reductions included bias subtraction and wavelength calibration. Fringing is particularly large at the red end of the spectra\footnote{Fringing is not so problematic using the GMOS-North detectors, where future observations are planned.} ($>$750\,nm). We tried correcting for this using the flat fields, however this proved particularly problematic given that our spectral flat fields were taken with slit widths much larger than the PSF of the star.
Flat fields taken with the calibration mask were corrupted. Given this, we decided not to apply flat-fielding. In differential spectrophotometry, flat fielding is not necessary providing the sensitivity ratios between the target and comparison stars in each wavelength channel remain constant. Of course, this is often not the case if there is movement of the spectra on the CCDs, seeing variations, etc. However, correcting all the data using the same flat field cannot account for such effects anyway (although it might mitigate against them), given spectral flat-fields are PSF (at least when using wide slits) and wavelength dependent. Simple flat-fielding is also unable to correct for severe fringing at the level required for transmission spectroscopy.

Wavelength calibrated spectra for the 3 stars were extracted using the {\sc gsextract} routine, with an aperture of 8 pixels (varying the aperture by a few pixels had little effect) after sky subtraction. A few pixel columns (the spatial direction) showing significant temporal variation were masked from the extraction.
Examples of extracted spectra of WASP-29 and the two comparison stars are shown in Fig.~\ref{fig:eg_spec}, showing uncorrected fringing effects at the red end. The spectra were divided up into several spectral regions, defining independent wavelength channels. To extract the light curves for each wavelength channel, the spectra were summed in each of the spectral regions to produce a flux time-series of the target and comparison stars for each channel. The target flux was then divided through by the sum of the comparison stars' flux to produce multiple light curves at each wavelength channel, which we hereafter refer to as the `spectral' light curves. We experimented with different numbers of wavelength channels, which affect the resolution and signal to noise, and finally settled on 30 across the whole spectral range, which are marked by the vertical dashed lines in Fig.~\ref{fig:eg_spec}. We divided each of the three chips into 10 channels each of the same pixel widths, ignoring the gaps between detectors in the images.

Correlated noise is present to some degree in all the light curves, but is particularly bad at the red end due to the fringing (and possibly the strong O$_2$ telluric feature at $\sim$7590\,\AA). We therefore decided to analyse only 15 spectral light curves at the blue end and discard the remaining 15 for the remainder of this analysis. We also experimented with the two comparison stars. Given that the second comparison star is significantly fainter, it was excluded from the analysis and we only used the brighter one. Finally, we also produced a `white' light curve, by summing up the flux over the first 10 wavelength channels (the other 5 contained larger systematics), prior to dividing through the target flux by the single comparison star's flux. The white light curve is shown in Fig.~\ref{fig:GP_white_lcv}, and the spectral light curves are shown in Fig.~\ref{fig:all_lightcurves}.

We calculated the theoretical photon noise for the white light curve and each spectral bin, taking into account the read noise and sky contribution (although negligible for high signal-to-noise data). Typical integrated electron counts per exposure for each spectral bin range from $\sim$3.7$\times10^6$ to 2.1$\times10^7$ for the target star, and $\sim$7.9$\times10^6$ to 2.6$\times10^7$ for the comparison star, and vary slightly throughout the night with airmass. For the white light curve the typical integrated counts per exposure were $\sim$1.3$\times10^8$ and 2.0$\times10^8$ electrons for the target and comparison star, respectively. This results in (time-averaged) theoretical precision on the relative flux per exposure ranging from 2.9$\times10^{-4}$ to 6.2$\times10^{-4}$ for the spectral light curves, and $\sim$1.1$\times10^{-4}$ for the white light curve.

We also extracted auxiliary data from our target and comparison star's spectra to investigate the cause of the systematics present in the light curves. The relative shift in the dispersion axis was measured by cross correlating each star's spectra with the first in the time-series. The relative shift in the cross-dispersion axis and the width of the spectral trace were found by fitting a Gaussian function to each column of the spectrum, and finding the average values per exposure. Over the course of the observations, the shift in the dispersion and cross-dispersion axes were $\sim$1.0 and 0.5 (binned) pixels respectively, and were the same for both stars (within the measurement error of $\sim0.1$ pixels), showing that no significant rotation of the field occurred. The FWHM in the cross dispersion direction ranged from 6--12 (binned) pixels, resulting in varying spectral resolution. We found no obvious correlations between these measurements and the systematics in the light curves, although the seeing variations are likely the most significant contributor given the relatively large changes.

\begin{figure}
\includegraphics[width=85mm]{ExSpec.pdf}
\caption{Example spectra extracted from a single exposure. The black line is WASP-29, and the blue and green lines are the two comparison stars. The vertical dashed lines mark the extraction regions, with the shaded regions marking the gaps in the target spectrum between the detectors. Only the first 15 channels were used in the final analysis, due to significant fringing at the red end.}
\label{fig:eg_spec}
\end{figure}

\begin{figure}
\includegraphics[width=85mm]{RawFit_wlc.pdf}
\includegraphics[width=85mm]{GPFit_wlc.pdf}
\caption{The `white' light curve produced from the first 10 wavelength channels of the GMOS spectra. Top: best fit model using a white noise analysis and the residuals. Clearly, significant systematics are present and must be accounted for in the fitting process. Bottom: best fit using a Gaussian process model with a Mat\'ern 3/2 covariance kernel (see text). The red line represents the best fit model, with the grey regions representing the 1 and 2\,$\sigma$ limits of the instrument model plus white noise times the transit model. The dashed red line is the projection of the systematics model without the transit, along with the 1 and 2\,$\sigma$ limits (now excluding the white noise term).}
\label{fig:GP_white_lcv}
\end{figure}

\section{Light Curve Analysis}
\label{sect:analysis}

The white light curve and the 15 spectral light curves were modelled in a variety of ways to account for the instrumental systematics. We first analysed them using a simple white noise model to inspect the residuals and establish the nature of the systematics. In all cases the transit model was constructed using the analytic equations of \citet{mandel_agol_2002}, and is similar to that described in \citet{Gibson_2008}. We assumed a circular orbit, and fixed the period ($P$) to 3.922727 days as given in \citet{Hellier_2010}. The remaining parameters of the transit model were the central transit time ($T_\text{C}$), the system scale ($a/R_\star$), the planet-to-star radius ratio ($\rho=R_\text{p}/R_\star$), the impact parameter ($b$), the two quadratic limb darkening parameters ($c_1$,$c_2$), and two parameters of a linear baseline model of time ($f_\text{oot}$, $T_\text{grad}$). We hereafter denote the transit model as $T(\bmath t,\bphi)$, where $\bmath t$ is a vector of time, $\bphi$ is the vector of transit parameters, and $\bmath f$ as the vector of flux measurements.

Unless otherwise stated, all of these parameters except $P$ were fitted to the white light curve, and $\rho$, $c_1$, $c_2$, $f_\text{oot}$ and $T_\text{grad}$ were fitted to the spectral light curves with the remaining parameters fixed to the {\it final} white light curve values (see Sect.~\ref{sect:results_wlc}). This is because we are interested in finding the relative values for $\rho$ as a function of wavelength to produce our transmission spectrum, and we can condition on the values of $T_\text{C}$, $a/R_\star$ and $b$, that change the inferred values for $\rho$ in the same way.

Given the lack of correlations found between the systematics and the auxiliary information extracted (see Sect. \ref{sect:observations}), we model the systematics as time-correlated noise. The only differences in the following analyses are the noise models used to account for instrumental systematics. These are described in turn.

\begin{figure*}
\includegraphics[width=170mm]{CombinedLCV_withres.pdf}
\caption{Spectral light curves with the central wavelengths marked, with a linear trend in time removed. Left panel: light curves fitted with the simple white noise model.  Middle panel: residuals from the white noise fits. Right panel: same light curves as left panel with their best fit Gaussian process model in red. The grey shading represents the 1 and 2\,$\sigma$ limits of the GP model (including white noise).}
\label{fig:all_lightcurves}
\end{figure*}

\subsection{White noise analysis}
\label{sect:whiteanalysis}

We first analysed all of the light curves using a simple white noise model, with likelihood given by
\begin{equation}
\label{eq:white_noise_likelihood}
p(\bmath f| \bmath t,\bphi) = \mathcal{N} \left (T(\bmath t,\bphi) , \sigma_\text{w}^2\mathbfss{I} \right),
\end{equation}
where $\mathcal{N}(\bmu,\mathbf\Sigma)$ is the multivariate normal distribution with mean $\bmu$ and covariance matrix $\mathbf\Sigma$, $\sigma_\text{w}$ is the uncertainty of each data point and \mathbfss{I} is the identity matrix. In other words we have a diagonal covariance matrix with all diagonal terms equal to $\sigma_\text{w}^2$, representing a standard i.i.d. noise model as is commonly used to model transits.

We then multiply the likelihood by the prior on the transit and noise parameters, $p(\bphi,\sigma_\text{w}^2)$, to produce the joint posterior probability distribution, and use a Monte-Carlo Markov-Chain (MCMC) to explore the posterior distribution and produce marginal probabilities for each of the model parameters. In practice we do not explicitly state priors for most parameters, implying uniform, improper priors. The exceptions are for the limb darkening parameters and impact parameter\footnote{strictly speaking we should apply the same prior to $\rho$ and $a/R_\star$, but this would have no effect on the inference.}, where we restrict the parameter to be positive using a step function of the form;
\[
p(x) = \left\{
\begin{array}{l}
0,~\mathrm{if}~x < 0 \\
1,~\mathrm{if}~x \geq 0 \\
\end{array}
\right.,
\]
specifying another improper prior. We also restrict the sum of the two limb darkening parameters $c_1+c_2 \leq 1$ in a similar way, to ensure that the brightness of the stellar surface is positive. Four MCMC chains were run for each light curve, of length 100\,000. We excluded the first 10\% of each chain, and verified convergence by checking the Gelman \& Rubin (GR) statistic \citep{GelmanRubin_1992}. The light curves along with their best fit models are shown in Figs.~\ref{fig:GP_white_lcv} and \ref{fig:all_lightcurves} for the white and spectral light curves, respectively.

Clearly, the white noise model is incapable of accounting for the correlated noise in the light curves, as seen in the residuals. We therefore analyse the residuals of each of the best fit models in an effort to understand the form of the systematics. We first use the time-averaging method to obtain a simple estimate of the red noise \citep{Pont_2006}, following the procedure of \citet{Winn_2008}, where the residuals are averaged into bins of width $N$, and the RMS is calculated as a function of N. The noise should drop by $1/\sqrt{N}$ if it is uncorrelated in time\footnote{actually, $\frac{1}{\sqrt{N}}\sqrt{\frac{M}{M+1}}$, where $M$ is the number of bins}. See \citet{Gibson_2009} for a more detailed description of this procedure. Fig.~\ref{fig:residuals_analysis} shows an example of this for one of the light curves, clearly showing that there is time-correlated noise in the light curves. As a first attempt to account for correlated noise, we calculate the factor $\beta$, which is the ratio of the RMS vs $N$ plot to the theoretical noise in the white case. We chose the maximum value for this, and then scale the noise parameter, $\sigma_\text{w}$, by this value, fix it, and re-fit the light curves using artificially inflated error bars to account for the systematics using the same MCMC procedure as before. $\beta$ ranged from $\sim$2.4 to 3.8. We will hereafter refer to this as the `white noise plus $\beta$' model.

This method is a useful way to estimate the additional uncertainties expected in the presence of systematic noise. However, it does not allow the form of the correlations to be modelled and therefore cannot produce more accurate parameter estimates (this is discussed in \citealt{Carter_2009}), and therefore we consider more sophisticated models in the following sections.

\begin{figure}
\includegraphics[width=85mm]{Residuals_13.pdf}
\caption{Residuals analysis of one of the spectral light curves (6895\,\AA).
Top panel: Residuals from the best fit white noise model. The dashed lines represent the white noise RMS.
Second panel: Plot of the RMS as a function of bin size $N$. The dashed line represents the theoretical curve expected for white noise, and the ratio of the two curves gives the factor $\beta$, used to re-scale the uncertainties for the simple rescaling systematics model.
Third panel: Autocorrelation function of the residuals. The solid green line marks the best fit Mat\'ern 3/2 kernel from the GP model, and the dashed red line zero correlation.
Bottom panel: Power spectral density of the residuals. The dashed line marks a $1/f$ envelope.
}
\label{fig:residuals_analysis}
\end{figure}

\subsection{Gaussian process analysis}
\label{sect:GPanalysis}

We use a Gaussian process to model the time-correlated noise, similar to that described in \citet{Gibson_2012} and \citet{Gibson_2012b}, allowing us to model instrumental systematics as a stochastic rather than a deterministic process. This avoids the need to specify a parametric form for the systematics, which is often impossible to do, whilst also allowing for a much more flexible model. Furthermore, GPs are intrinsically Bayesian, thus avoiding the possibility of over-fitting systematics. The combination of a very flexible model and principled Bayesian inference effectively allows one to marginalise out any ignorance about the form of the systematics model and account for it whilst inferring transit parameters. This is extremely challenging to achieve using parametric models as it requires calculation of the Bayesian evidence (and therefore proper, usually informative, priors on the model parameters), and perhaps even marginalisation over many possible instrument models.

Here we briefly describe the GP model, and refer to \citet{Gibson_2012} and references therein for further details. A GP is a collection of random variables, any finite subset of which have a joint Gaussian distribution. Therefore we can write our GP as
\begin{equation}
\label{eq:gp}
p(\bmath f| \bmath t,\bphi,\btheta) = \mathcal{N} \left (T(\bmath t,\bphi) , \mathbf{\Sigma} (\bmath t,\btheta) \right).
\end{equation}
The only difference to Eq.~\ref{eq:white_noise_likelihood} is that we now consider the off diagonal elements of the covariance matrix. As well as specifying a {\it mean function}, in this case the transit function, with a GP model we must also specify a {\it kernel function} which populates the elements of the covariance matrix and has hyperparameters\footnote{A Gaussian process, $\mathcal{GP}(\bmu,\mathbf\Sigma)$, is fully specified by its mean $\bmu$ and covariance $\mathbf\Sigma$. Parameters of both the mean function and kernel are known as hyperparameters.} $\btheta$, written as:
\[
\mathbf\Sigma_{nm} = k({t}_n,{t}_m | \btheta).
\]
We will discuss the choice of the kernel function later.

The above GP uses the kernel to model the {\it residuals} from the light curve model. Alternatively, we can model the light curve as the transit light curve {\it times} a GP, if we wish our systematics model to be multiplicative rather than additive, i.e.:
\[
\bmath f = T(\bmath t,\bphi, \btheta) \times \mathcal{GP}(\mathbf 1, \mathbf{\Sigma} (\bmath t,\btheta)).
\]
In this case the joint probability distribution can be written for $\bmath f / T(\bmath t,\bphi)$ as
\[
p(\bmath f / T(\bmath t,\bphi, \btheta) | \bmath t,\bphi) = \mathcal{N} \left (\mathbf 1 , \mathbf{\Sigma} (\bmath t,\btheta) \right).
\]
In practice it makes little difference which model we choose, given that transit light curves are shallow, and we could also combine multiplicative and additive GPs. For the remainder of this paper we use the latter, multiplicative model.

We tested several different types of kernels to model the time-correlated noise, including the squared exponential, Mat\'ern and rational quadratic (see \citealt{Rasmussen_Williams} for a detailed discussion of kernels, which is beyond the scope of this paper). In a fully Bayesian analysis, we could calculate the Bayesian evidence for each kernel, and use it to choose the best kernel, or alternatively even marginalise over them. However, not only is this computationally prohibitive, but it would also require us to specify proper (therefore informative) priors on the hyperparameters, which would make the evidence somewhat subjective. We therefore selected a kernel based on analysis of the white noise model residuals, and by running tests on simulated light curves.

In the end we decided to use a Mat\'ern 3/2 kernel function, given by:
\begin{equation}
\label{eq:kernel}
k({t}_n, {t}_m | \btheta) = \xi^2 \left( 1+{\sqrt{3}\eta\Delta t} \right) \exp \left( -{\sqrt{3}\eta\Delta t}\right) + \delta_{nm}\sigma_\text{w}^2,
\end{equation}
where $\xi$ is a hyperparameter that specifies the maximum covariance, $\Delta t = |t_n-t_m|$ is the time difference, $\eta$ is the inverse characteristic length scale, and $\delta$ is the Kronecker delta. This kernel can be seen as a rougher version of the commonly used squared exponential kernel (i.e. that used in \citealt{Gibson_2012,Gibson_2012b}). One major motivation for this kernel is that is gives the best match to the auto-correlation function (ACF) for most of the residuals from the white noise model. This provides an estimate of how the data points are correlated with one another as a function of time lag, and an example is shown in Fig.~\ref{fig:residuals_analysis}. The green line in the ACF function marks the best fit covariance kernel (although in practice we marginalise over the kernel parameters). We also ran a series of tests on simulated light curves, described in the Appendix, which further validate the use of the Mat\'ern kernel. We stress that selecting a kernel, whilst in some ways analogous to selecting a parametric model, allows a much more flexible model than any parametric form and is also intrinsically Bayesian. In addition, we ran much of the same analysis using the squared exponential kernel. This gave similar results, and we therefore conclude that the choice of kernel is not critical for this particular dataset.

Analogously to the white noise case, we can specify priors for all the hyperparameters of the model, and multiply by the marginal likelihood\footnote{called the marginal likelihood because in a GP we have marginalised over all the possible functions for each set of hyperparameters} to produce the posterior joint probability distribution. This can then be optimised with respect to the hyperparameters using a Nelder-Mead simplex algorithm, or alternatively the marginal parameter distributions for each hyperparameter can be obtained by exploring the posterior distribution with an MCMC in just the same way as for the white noise model. The same priors were applied to the limb darkening parameters and $b$, and we also specified hyperpriors for the hyperparameters $\xi$ and $\eta$. These took the form of Gamma distributions with shape parameter unity, given by
\[
p(x) = \left\{
\begin{array}{ll}
0, & ~\mathrm{if}~x < 0 \\
\dfrac{1}{l} \exp\left(-x/l\right), & ~\mathrm{if}~x \geq 0 \\
\end{array}
\right.,
\]
where $l$ is the length scale of the hyperprior. We set the length scales for $\xi$ and $\eta$ as $10^{-3}$ and $200$, respectively. These were not specified to influence the results of the inference, but rather to ease convergence of the MCMC chains (when both parameters are small they are unconstrained by the likelihood). Indeed, we checked that the length scales of the hyperpriors did not affect the inferred transit parameters.

Whilst GPs are rather simple in theory, each evaluation of the marginal likelihood requires inversion of the covariance matrix, which makes a full marginalisation tedious as it requires $\mathcal{O}(n^3)$ computations, and limits full GP analyses to relatively small datasets. Therefore in addition to the full GP marginalisation, we also use a technique called Maximum Likelihood type-II (ML-II), where the hyperparameters of the covariance kernel are fixed to their maximum-posterior values, and we marginalise over the remaining transit parameters. Once the covariance hyperparameters are fixed this negates the need to invert the covariance matrix. This is a valid approximation when the posterior is sharply peaked with respect to the covariance hyperparameters, and is particularly useful when running many tests on the data, although for our final results we always use fully marginalised GPs (i.e. we marginalised over the covariance hyperparameters as well).

We ran four chains for all light curves of length 100\,000 and 50\,000 for the ML-II and full marginalisations, respectively. Each chain of length 50\,000 for the full marginalisation took about 17.5 minutes to compute using a single core on a standard desktop, compared to about 2 minutes and 1 minute for 100\,000 length chains with the ML-II method and white noise model, repectively. The run times do not scale as badly as one might expect with $\mathcal{O}(n^3)$ complexity, because the computation of the light curve model rather than the likelihood dominates for the simpler noise models. We tested for convergence in the same way as before. Due to degeneracy in the linear baseline and the GP model, we decided to fix the $T_\text{grad}$ parameter for the full GP marginalisation. This requires shorter chains to reach convergence, and doesn't effect the results as $\rho$ did not significantly correlate with $T_\text{grad}$. We verified this with longer chains for a subset of the light curves. In a few cases, $\xi$ and/or $\eta$ did not fully converge (the GR statistic was a few percent from unity), but only when one or both parameters were consistent with zero; however, in all cases $\rho$ converged and therefore the transmission spectrum is not affected. The best fit GP models to the light curves are shown in Figs.~\ref{fig:GP_white_lcv} and \ref{fig:all_lightcurves}. To illustrate correlations in the parameters and convergence of the final marginal distributions, 1D and 2D marginal distributions are shown in Fig.~\ref{fig:wlc_correlations} for the white light curve and one spectral light curve.

\begin{figure*}
\includegraphics[width=170mm]{GPFit_both_correlations.pdf}
\caption{1D and 2D marginal distributions from the posterior probability distribution for the Gaussian process noise model. The lower left shows the probability distributions for the white light curve, and the upper right for one of the spectral light curves. The black lines in the 2D distributions mark the 1 and 2\,$\sigma$ limits, respectively. Distributions for the four separate MCMC chains are shown in the 1D histograms.}
\label{fig:wlc_correlations}
\end{figure*}

\subsection{Wavelet analysis}
\label{sect:waveletanalysis}

As we have used a GP to fit for time-correlated noise only, it is worth exploring other methods in the literature to account for time-correlated noise. One such method is the wavelet method introduced by \citet{Carter_2009}. It is valid when the power spectral density (PSD) of the noise takes the form:
\[
\mathrm{PSD} \propto \frac{1}{f^\gamma},
\]
where $f$ is the frequency and $\gamma$ is the spectral index. The method is based on taking a wavelet transform of the residuals from the best fit model, and the likelihood is computed from the wavelet coefficients. The central idea of this method is that it diagonalises the covariance matrix of Eq.~\ref{eq:gp} when specified in the wavelet domain. Like GPs it models the systematics as a stochastic process; however, it has the significant advantage that a burdensome matrix inversion is not required for each evaluation of the likelihood and is therefore much faster. The extra computation required as compared to the white noise analysis is a fast wavelet transform, which is of order $\mathcal{O}(n)$, and therefore does not significantly extend the computation time.

In order to justify the wavelet method we analysed the PSD of all the residuals from the white noise fit as recommended by \citet{Carter_2009}. An example of this is shown in Fig.~\ref{fig:residuals_analysis}. In most cases the noise appeared consistent with $1/f^\gamma$, by which we mean that low frequency components dominate, but the shape itself is hard to determine (the PSD should appear as noise within a $1/f^\gamma$ envelope). This supports the use of the wavelet based likelihood for our light curves, or at least suggests the noise properties are nearly $1/f^\gamma$.

Our likelihood therefore took the form of Eq.~32 from \citet{Carter_2009}, and we multiplied by priors similarly to Sect.~\ref{sect:whiteanalysis} to produce a posterior distribution. In addition to the white noise parameter $\sigma_\text{w}$, the wavelet likelihood has a red noise parameter $\sigma_\text{r}$ specifying the amplitude of the red noise component, and the spectral index $\gamma$. We fixed $\gamma$ to 1 (as in \citealt{Carter_2009}), and optimised and explored the joint posterior distribution with respect to the free transit parameters, plus $\sigma_\text{w}$ and $\sigma_\text{r}$. Chains lengths and run times were approximately the same as for the white noise analysis.

\subsection{Removal of common mode systematics}
\label{sect:commonmode}

Finally, we experimented with a simple method to remove common signals observed in the spectral light curves, to see if we could increase the precision in our transmission spectrum.
We were motivated to do this as we obtain a reduced $\chi^2$ significantly lower than one for our transmission spectrum using all of the noise models (see Sect.~\ref{sect:trans_spectrum}). This is extremely unlikely to happen by chance, and is probably the result of similar systematic signals in the spectral light curves. Noise models will take all signals into account when calculating the uncertainties in $\rho$ (as they should); however, common signals will increase the uncertainties for each point in the transmission spectrum. As we are trying to find the relative change in $\rho$ with transmission spectroscopy (hence why we condition on fixed values of $T_\text{C}$, $a/R_\star$ and $b$), we tried to remove a common signal in all the light curves.

This common signal is evident in Fig.~\ref{fig:all_lightcurves}, where many of the residuals appear to have the same shape. This signal is also very similar to the GP model fitted to the white light curve (red dashed line, Fig.~\ref{fig:GP_white_lcv}). Taking this signal into account with the GP (or any other method) will increase the uncertainties in the calculated values of $\rho$, but this signal should not affect the relative values for $\rho$. We therefore divided through each spectral light curve by the GP systematics model fitted to the white light curve prior to fitting each light curve with the methods described above. This will remove the common signal and allow the noise models to produce more independent data points. Of course the underlying physical signal may not be identical for all the light curves, but the noise models described above can also take into account any excess signal added or not removed by the procedure in the same way that they account for `normal' systematics. In theory we could model this in more principled ways, e.g. by modelling all light curves simultaneously with a common signal plus independent signals, or perhaps by using the white light curve instrument model as an input for the systematics model, but we choose not to pursue them here and follow this simple procedure. Each fit for the spectral light curves as described in Sects.~\ref{sect:whiteanalysis}, \ref{sect:GPanalysis} and \ref{sect:waveletanalysis} was repeated by first dividing the light curve by the GP noise model for the white light curve. This procedure appeared to remove much of the common signal, and allowed for a more precise determination of the transmission spectrum, as described in the following section.

\section{Results and Conclusions}
\label{sect:results}

\subsection{White light curve analysis}
\label{sect:results_wlc}

The results from the white light curve fit are given in Tab.~\ref{tab:wlc_results} for three of the noise models used: white noise plus $\beta$ rescale, wavelets, and the Gaussian process after marginalising over all the transit parameters and covariance hyperparameters. The derived parameters are consistent with those reported in \citet{Hellier_2010} and \citet{Dragomir_2011}. In general, the three noise models gave consistent results. The simple white noise model produced the smallest uncertainties (which we do not reproduce here), whereas the white noise plus $\beta$ model gave the largest uncertainties. The wavelet and GP models both gave something in between. This is expected given that they both try to account for the form of the systematics model when inferring transit parameters, rather than just scaling the uncertainties to account for it. The fact that the wavelet and GP methods are not only consistent, but produce similar uncertainties is a strong validation of both techniques. We briefly discuss the relative merits of both methods in Sect.~\ref{sect:discussion}, and for the remainder of this paper we adopt the GP results, given that we cannot verify the $1/f$ nature of the noise, and they provide more conservative uncertainties. This maybe indicates that we are taking into account a larger range of possible systematic signals. The white noise values fitted for the GP and wavelet methods are $\sim$3.8 and 2.9 times the theoretical photon noise. The difference is perhaps due to the wavelet method absorbing some of the white noise into the systematic component (see e.g. Figs. 1 and 4 of \citealt{Carter_2009}, where the $1/f$ noise appears to contain a white component). Further sources of noise are accounted for by the time-correlated component of these models, given by the max covariance and red noise parameters quoted in Tab.~\ref{tab:wlc_results}.

Using the distributions from the MCMC chains, we calculate further system parameters of WASP-29b. These are given in Tab.~\ref{tab:wlc_results_derived}. Where the distributions were not available, we generated normal distributed values from the values given in the literature to propagate the uncertainties in the calculations properly. Again these results are consistent with the values reported in \citet{Hellier_2010} and \citet{Dragomir_2011} within the uncertainties, further confirming WASP-29b's status as a Saturn-like exoplanet.

\begin{table*}
\caption{Parameters from the MCMC fits of the white light curves, given for three different noise models: white noise plus $\beta$, wavelets, and the fully marginalised Mat\'ern 3/2 Gaussian process.}
\label{tab:wlc_results}
\begin{tabular}{llll}
\hline
\noalign{\smallskip}
~ & \multicolumn{3}{c}{Noise model}\\
Parameter & White Noise $+\beta$ & Wavelet & GP\\
\hline
\noalign{\smallskip}
Central transit time, $T_\text{C}$ (HJD$_{UTC}$) & $2455853.72515^{+0.00028}_{-0.00028}$ & $2455853.72469^{+0.00015}_{-0.00015}$ & $2455853.72442^{+0.00016}_{-0.00017}$ \\\noalign{\smallskip}
Period, $P$ (days) & $3.922727$ (fixed) & $3.922727$ (fixed) & $3.922727$ (fixed) \\\noalign{\smallskip}
System scale, $a/R_\star$ & $11.95^{+0.34}_{-0.67}$ & $12.36^{+0.12}_{-0.21}$ & $12.36^{+0.13}_{-0.22}$ \\\noalign{\smallskip}
Planet-star radius ratio, $R_\text{p}/R_\star$ & $0.0998^{+0.0026}_{-0.0017}$ & $0.0984^{+0.0010}_{-0.0009}$ & $0.0982^{+0.0015}_{-0.0015}$ \\\noalign{\smallskip}
Impact parameter, $b$ & $0.25^{+0.16}_{-0.16}$ & $0.14^{+0.10}_{-0.08}$ & $0.14^{+0.11}_{-0.09}$ \\\noalign{\smallskip}
Linear limb darkening parameter, $c_1$ & $0.721^{+0.076}_{-0.091}$ & $0.745^{+0.039}_{-0.047}$ & $0.698^{+0.059}_{-0.089}$ \\\noalign{\smallskip}
Quad limb darkening parameter, $c_2$ & $0.154^{+0.135}_{-0.094}$ & $0.090^{+0.095}_{-0.051}$ & $0.142^{+0.142}_{-0.084}$ \\\noalign{\smallskip}
Out-of-transit flux, $f_\text{oot}$ & $0.99969^{+0.00009}_{-0.00009}$ & $0.99963^{+0.00008}_{-0.00008}$ & $0.99965^{+0.00031}_{-0.00033}$ \\\noalign{\smallskip}
Time gradient, $T_\text{grad}$ & $-0.00014^{+0.00005}_{-0.00005}$ & $-0.00008^{+0.00004}_{-0.00003}$ & $-0.00008$ (fixed) \\\noalign{\smallskip}
GP max covariance, $\xi$ & $-$ & $-$ & $0.00058^{+0.00033}_{-0.00018}$ \\\noalign{\smallskip}
GP inverse length scale, $\eta$ & $-$ & $-$ & $25.7^{+14.9}_{-9.5}$ \\\noalign{\smallskip}
White noise, $\sigma_\text{w}$ & $0.00119$ (fixed) & $0.000324^{+0.000015}_{-0.000014}$ & $0.000426^{+0.000017}_{-0.000016}$ \\\noalign{\smallskip}
Red noise, $\sigma_\text{r}$ & $-$ & $0.00251^{+0.00033}_{-0.00030}$ & $-$ \\\noalign{\smallskip}\noalign{\smallskip}
\hline
\end{tabular}
\end{table*}

\begin{table}
\caption{WASP-29 parameters derived from the MCMC posterior distribution of the GP fits.}
\label{tab:wlc_results_derived}
\begin{tabular}{lll}
\hline
Parameter & Value & Unit \\
\hline
\noalign{\smallskip}
Transit epoch, $T_0$ & $2455830.18811^{+0.00016}_{-0.00016}$ & HJD$_{UTC}$ \\\noalign{\smallskip}
Period, $P$ & $3.9227186^{+0.0000068}_{-0.0000068}$ & days \\\noalign{\smallskip}
Transit duration, $T_{14}$ & $0.11036^{+0.00071}_{-0.00063}$ & days \\\noalign{\smallskip}
Inclination, $i$ & $89.17^{+0.50}_{-0.56}$ & deg \\\noalign{\smallskip}
Semi-major axis, $a$ & $0.04565^{+0.00060}_{-0.00062}$ & AU \\\noalign{\smallskip}
Stellar radius$^\alpha$, $R_\star$ & $0.808\pm0.044$ & $R_\odot$ \\\noalign{\smallskip}
Planet mass$^\alpha$, $M_\text{p}$ & $0.244\pm0.020$ & $M_J$ \\\noalign{\smallskip}
Planet radius, $R_\text{p}$ & $0.776^{+0.043}_{-0.043}$ & $R_J$ \\\noalign{\smallskip}
Planet density, $\rho_\text{p}$ & $0.53^{+0.11}_{-0.09}$ & $\rho_J$ \\\noalign{\smallskip}
Surface gravity, $\log g_\text{p}$ & $3.00^{+0.06}_{-0.06}$ & [cgs] \\\noalign{\smallskip}
Equilibrium temp, $T_\text{p}$ & $970^{+32}_{-31}$ & K \\\noalign{\smallskip}\hline
\noalign{\smallskip}
\multicolumn{3}{l}{{\footnotesize\,$^\alpha$\,Adopted from \citet{Hellier_2010}.}} \\
\noalign{\smallskip}
\end{tabular}
\end{table}

\subsection{Transit Ephemeris}

A new ephemeris was calculated for WASP-29 using the transit time derived for the white Gemini transit light curve, plus the ephemeris reported in \citet{Hellier_2010} and the transit time from \citet{Dragomir_2011}. These were converted to HJD$_{UTC}$ format, and a straight line of the form
\[
T_\text{C}(E) = T_\text{C}(0) + PE
\]
was fitted to the three transit times. The zero-point epoch was set equal to as near to the centre of mass of the three points as possible, weigthed as $1/\sigma_{T_\text{C}}^2$. This was to minimise the covariance between the transit epoch and the period, which was verified after the fit. The chosen $E=0$ transit was 6 periods prior to the Gemini transit, giving epochs of -130, -98 and 6 for the three transits. The new ephemeris is reported in Tab.~\ref{tab:wlc_results_derived}.

\subsection{Transmission spectrum}
\label{sect:trans_spectrum}

The transmission spectra produced via the noise models discussed in Sect.~\ref{sect:analysis} are shown in Fig.~\ref{fig:transpec}. These are prior to removal of the common mode systematic.
The horizontal dashed lines represent the weighted average and plus and minus three scale heights, calculated to be $\sim$360\,km (one scale height of 360\,km corresponds to $\sim1.3\times10^{-4}$ in transit depth). In general the transmission spectra are all in broad agreement, and are remarkably flat, showing a featureless spectrum at the few parts in $10^{-4}$ level. Uncertainties were smallest for the simple white noise model (not shown) and largest for the full GP marginalisation. However, the dispersion of the points (i.e. the scatter around the average) was smallest for the full GP marginalisation. This implies that the GP model is doing a particularly good job at determining the correct value for $\rho$, but for some reason overestimates the uncertainty (it is highly unlikely a draw from random noise would lead to such small dispersion). In fact, all of the noise models presented give reduced $\chi^2$ significantly smaller than 1. We propose that this is due to a common mode systematic that the noise models take into account in a similar way for each wavelength channel, and led to the correction for this as discussed in Sect.~\ref{sect:commonmode}.

Fig.~\ref{fig:transpec_subsig} shows the transmission spectra produced after the common mode systematic correction. The spectra produced by all noise models were again consistent, with the full GP model giving uncertainties typically $30-40\%$ larger than the other models. We therefore choose to adopt these as our final uncertainties, for the reasons discussed in Sect.~\ref{sect:results_wlc}. These results are given in Tab.~\ref{tab:trans_spec_single}. The transmission spectrum is still consistent with a flat model, but we are able to place stronger constraints on the atmosphere of WASP-29b using the common mode correction. For the common mode corrected light curves, the white noise values fitted for the GP model ranges from 1.51 to 2.12 times the theoretical noise, and from 1.11 to 1.50 for the wavelet model. Similarly to the white light curve, further sources of noise are accounted for by the systematic component of these models. The lower ratios between the actual and theoretical white noise for the spectral light curves as compared to the white light curve perhaps indicate that atmospheric transmission corrections using the comparison star are best done in narrow wavelength ranges.

Fig.~\ref{fig:transpec_single} shows the transmission spectrum produced with the full GP model, now with several model transmission spectra of WASP-29b overplotted. These forward models were produced using the {\sc nemesis} retrieval tool \citep{Irwin_2008},  a radiative transfer code originally developed to investigate the atmospheres of Solar system planets, and recently adapted for exoplanet transmission spectra \citep[][Barstow et al. {\it submitted}]{Lee_2012}. The grey line shows a model containing a purely H$_2$ and He atmosphere. The green, red and blue lines are with 100 ppmv H$_2$O and 1, 5 and 10 ppmv of Na and K added, respectively. Gas absorption line data is from \citet{Rothman_2010} and \citet{Kupka_2000} for the H$_2$O and the alkali metals, respectively. Given the precision of the transmission spectrum, we do not attempt a detailed retrieval here, rather the models are plotted for reference to show the scale of potential features.

Using our data we can only realistically rule out cloud-free atmospheres with significant amounts of Na, given the lack of a pressure broadened feature.
However, a Na rich atmosphere with thick clouds or Rayleigh scattering haze is not ruled out.
Of course, another explanation is that elemental Na is simply not present. This is likely for atmospheres of $\sim$1000\,K or cooler \citep[e.g.][]{Burrows_2000}, where atomic Na can be lost in compounds such as disodium monosulphide (Na$_2$S) or ansite (NaAlSi$_3$O$_8$), depending on the presence of other species in the atmosphere. Similar observations at shorter wavelengths could distinguish between a flat featureless spectrum and one dominated by a Rayleigh scattering haze, such as the prototypical HD 189733b \citep{Pont_2008,Sing_2012}, and higher resolution (at similar S/N) is required to rule out the presence of a Na core if there are clouds or haze present. For comparison the HD 189733b transmission spectrum varies by about 2 scale heights over this spectral range due to a Rayleigh scattering haze. Despite the lack of constraints we can place on the atmosphere, we can rule out WASP-29b having a similar atmosphere to the other prototypical hot Jupiter HD 209458b, given the lack of a strong pressure broadened Na feature. As stated, Na is likely to form compounds below $\sim$1000\,K. Given its prominence in HD 209458b's transmission spectrum, this could indicate another significant transition in classes of hot Jupiter atmospheres.

\begin{table}
\caption{Transmission spectrum of WASP-29b using the full GP marginalisation and after removal of the common mode systematic.}
\label{tab:trans_spec_single}
\begin{tabular}{ccc}
\hline
\noalign{\smallskip}
\noalign{\smallskip}
Wavelength & $\rho$ & $\Delta\rho$\\
\hline
5211.1\AA & 0.0961 & 0.0027\\
5348.0\AA & 0.0979 & 0.0016\\
5485.0\AA & 0.0971 & 0.0018\\
5621.9\AA & 0.0970 & 0.0016\\
5758.8\AA & 0.0973 & 0.0013\\
5895.8\AA & 0.0974 & 0.0013\\
6032.7\AA & 0.0976 & 0.0013\\
6169.6\AA & 0.0981 & 0.0013\\
6306.6\AA & 0.0970 & 0.0014\\
6444.9\AA & 0.0969 & 0.0014\\
6616.0\AA & 0.0973 & 0.0016\\
6755.7\AA & 0.0975 & 0.0032\\
6895.4\AA & 0.0969 & 0.0019\\
7035.0\AA & 0.0970 & 0.0032\\
7174.7\AA & 0.0970 & 0.0018\\
\noalign{\smallskip}
\hline
\end{tabular}
\end{table}

\begin{figure}
\includegraphics[width=85mm]{transpec_plot.pdf}
\caption{Transmission spectra of WASP-29b produced by the various noise models prior to removal of the common mode systematic. The horizontal dashed lines represent the weighted average and plus and minus three scale heights, calculated to be $\sim$360\,km (one scale height of 360\,km corresponds to $\sim1.3\times10^{-4}$ in transit depth). In all cases the $\chi^2$ of a flat model is considerably lower than 1, indicating that the uncertainty in the relative planet-to-star radius ratio might be overestimated due to common mode systematics.}
\label{fig:transpec}
\end{figure}

\begin{figure}
\includegraphics[width=85mm]{transpec_plot_subsig.pdf}
\caption{Same as Fig.~\ref{fig:transpec} with the common mode systematic removed from all light curves prior to fitting.}
\label{fig:transpec_subsig}
\end{figure}

\begin{figure*}
\includegraphics[width=142mm]{transpec_plot_single.pdf}
\caption{Transmission spectrum of WASP-29b using the full Mat\'ern 3/2 GP noise model, and after removal of the common mode systematic. The horizontal dashed lines represent the weighted average and plus and minus three scale heights, calculated to be $\sim$360\,km (one scale height of 360\,km corresponds to $\sim1\times10^{-4}$ in transit depth).
The grey line shows a model containing a purely H$_2$ and He atmosphere. The green, red and blue lines are with 100 ppmv H$_2$O and 1, 5 and 10 ppmv of Na and K added, respectively.
These are not fitted to the data, but simply over plotted for reference. The models are plotted with a resolution of 50\AA; the points do not significantly change when plotted at the resolution of the spectrum.}
\label{fig:transpec_single}
\end{figure*}

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

We have presented Gemini GMOS observations of the transmission spectrum of WASP-29b, using the technique of differential spectro-photometry. Using a single comparison star, we reached precision on the transit depth of $\sim1\times10^{-4}$ showing that GMOS can provide precision spectrophotometry at the level needed to probe the atmospheres of extrasolar planets.

Using the `white' light curve, we refined the system parameters and ephemeris for WASP-29, finding them to be consistent with previous studies \citep{Hellier_2010,Dragomir_2011}.
Despite picking WASP-29 as a test case for GMOS transmission spectroscopy, the precision attained allows us to rule out a Na rich, cloud/haze free atmosphere, given the lack of a pressure broadened Na feature. This indicates that Na is not present in the atmospheres of cooler `hot' Jupiters, or that clouds and/or hazes play an important role and mask the pressure broadened alkali metal signatures in WASP-29b's upper atmosphere. The former explanation is perhaps more likely, although a spectrum covering a larger wavelength region is required to confirm this, and it is of course possible that both are true. A higher resolution spectra (at similar S/N) is required to rule out the presence of a narrow Na core if clouds or hazes dominate. We note that this represents the first transmission spectrum of a hot-Saturn planet.

We have also presented a detailed analysis and comparison using various types of noise models to account for the GMOS systematics. Rather than search for correlations with observational parameters such as seeing, airmass, etc., we decided to focus on blind methods to account for the systematics, i.e. with no additional input parameters used to model the systematics. The methods tried include a white noise model, a simple rescaling of the photometric uncertainties, the wavelet method of \citet{Carter_2009}, and the Gaussian process model of \citet{Gibson_2012} applied to time correlated noise. The more sophisticated methods gave similar uncertainties, verifying their usefulness for analysis of time-correlated systematics.

In general, we restate some of the conclusions of \citet{Carter_2009}, that any model taking into account time-correlated noise is better than ignoring it, and that an analysis of the residuals using ACFs and PSDs are especially useful in guiding the choice of noise model required.
The wavelet and GP models give similar results and uncertainties for the white light curve, and within about $40\%$ for the uncertainties in the transmission spectrum. Despite GPs giving a slightly more conservative estimate of the uncertainties for the GMOS data, the wavelet method is perhaps preferred in general for the analysis of time-correlated noise when the PSD of the residuals follow a $1/f^\gamma$ distribution, given its much faster execution time (although we note that this is hard to verify for individual datasets, and care must be taken for low significance results). The GP method is more general, and can be applied to almost any noise model given a suitable choice of kernel (potentially even non-stationary noise), non-regularly spaced data, and can incorporate arbitrary numbers of input vectors into the stochastic function. This allows physical systematics models to be folded into the stochastic part of the GP, therefore allowing principled Bayesian inference of the instrument model and negating the need to specify the instrument model in closed form \citep{Gibson_2012,Gibson_2012b}. However, this added functionality comes at a significant runtime cost, and restricts the use of GPs up to datasets of $\sim$1000 points (at least using full marginalisation over the hyperparmeters with MCMC methods). Investigations into sparse GP models may allow their application to larger datasets \citep[e.g.][]{quinonero2005unifying}. We finally note that given the difficulty in dealing with systematic noise, the use of multiple, complimentary techniques is desirable where possible, although perhaps the only truly robust way to confirm results is to repeat measurements.

\section*{Acknowledgments}

Based on observations obtained at the Gemini Observatory, which is operated by the Association of Universities for Research in Astronomy (AURA) under a cooperative agreement with the NSF on behalf of the Gemini partnership: the National Science Foundation (United States), the Science and Technology Facilities Council (United Kingdom), the National Research Council (Canada), CONICYT (Chile), the Australian Research Council (Australia), CNPq (Brazil) and CONICET (Argentina). We are extremely grateful for the support provided by the Gemini staff. We also thank the referee for comments that improved the paper.
N. P. G and S. A. acknowledge support from STFC grant ST/G002266/2. J. K. B. acknowledges the support of the John Fell Oxford University Press (OUP) Research Fund for this research. L.N.F. was supported by a Glasstone Fellowship at the University of Oxford.