text:
\section{Introduction}
\label{SectionIntro}

The modern complexity theory of Convex Optimization
originated in the book of~\cite{nemirovskii1983problem},
where the first lower bounds were constructed for many different
classes of optimization problems and algorithms.
After this work, it became clear that
the best possible rate of convergence for a given family of methods
is fundamentally limited by a particular problem class.

For example, in the case of minimizing the convex functions with bounded subgradient,
by black-box local methods, the best rate for decreasing the functional
residual is $O(k^{-1/2})$,
assuming the problem dimension is sufficiently large.
This is the rate of the classical Subgradient Method (\cite{shor2012minimization}),
so the method is \textit{optimal}.
At the same time, for the differentiable convex
functions with Lipschitz continuous gradient, the best rate is $O(k^{-2})$,
which is much better. The optimal Fast Gradient Method
with this rate of convergence was developed by~\cite{nesterov1983method}.

During the last years, we started to see more and more examples of problem classes
that significantly vary the standard picture of complexity theory.
Thus, there were established
lower complexity bounds for the functions
with H\"older continuous gradient (w.r.t. different norms) by \cite{guzman2015lower},
for the functions with second- and high-order derivatives being Lipschitz continuous
by \cite{agarwal2018lower} and \cite{arjevani2019oracle},
for the relatively smooth functions by \cite{dragomir2021optimal}.
Different extensions of the complexity theory to the randomized algorithms
were presented by \cite{bubeck2019complexity},
\cite{diakonikolas2019lower}, and
\cite{garg2021near}.

Meanwhile, the framework of \textit{composite optimization},
that was proposed by \cite{nesterov2013gradient}, \cite{beck2009fast},
provided us with much more flexibility in formulation of problem parameters.
Indeed, in composite problems, we can add to the smooth part
of the objective additional terms such as indicator of a convex set,
or regularizers of a different kind. Consequently, smoothness characteristics of the objective
(e.g. the Lipschitz constant) need not necessary be consistent with
uniformly convex properties of the composite part.

In this work, we investigate the complexity of
the optimization problems with composite regularizer
that is given by a power of a norm.
While having a variety of examples for the problems of this type,
it appears that the lower complexity bounds are not really covered
in the literature.

Let us consider  convex minimization problems with the following structure,
\begin{equation} \label{MainProblem}
\begin{array}{rcl}
\min\limits_{x \in \mathbb{R}^n} \Bigl\{
F(x) & \stackrel{\mathrm{def}}{=} & f(x) + \frac{\sigma}{p}\|x\|^p
\Bigr\},
\end{array}
\end{equation}
with some parameters $p \geq 2$, $\sigma > 0$, and
the standard Euclidean norm:
$$
\begin{array}{rcl}
\|x\| & \stackrel{\mathrm{def}}{=} & \biggl(
\sum\limits_{i = 1}^n \bigl(  x^{(i)} \bigr)^2
\biggr)^{\! 1/2}.
\end{array}
$$
Function $f$ is the main source of complexity of solving~\eqref{MainProblem}.
We assume that its gradient is H\"older continuous,
for some $\nu \in [0, 1]$ and $H_{\nu} > 0$
\begin{equation} \label{HolderGrad}
\begin{array}{rcl}
\| \nabla f(x) - \nabla f(y) \| & \leq & H_{\nu}\|x - y\|^{\nu},
\qquad \forall x, y \in \mathbb{R}^n.
\end{array}
\end{equation}

When $p = 2$, we obtain the classical $\ell_2$-regularization by the squared Euclidean norm.
Finding the minimum of~\eqref{MainProblem} is equivalent to computing the value
of the proximal operator of $f$ (\cite{moreau1965proximity}).
We know that the problem is \textit{strongly convex}.
The rate of convergence becomes \textit{linear}
for the functions with Lipschitz continuous gradient $(\nu = 1)$,
and $O(k^{-1})$ for bounded variation of gradient $(\nu = 0)$ (see, e.g. \cite{nesterov2018lectures}).
Thus, the cost of solving the regularized problem is much cheaper than minimizing $f$ solely.
It happens that these observations work for $p > 2$ as well.

For $p = 3$, minimization problem~\eqref{MainProblem} appears to be a subproblem
for computing one step of the Newton Method with Cubic Regularization~(\cite{nesterov2006cubic}).
The smooth part then is a convex second-order Taylor polynomial, which is a quadratic function,
$$
\begin{array}{rcl}
f(x) & = & \frac{1}{2}\langle Ax, x \rangle - \langle b, x \rangle, \qquad A \succeq 0,
\end{array}
$$
and condition~\eqref{HolderGrad} is satisfied
with $\nu = 1$ and constant $H_1 = \lambda_{\max}(A)$.
In the following years, the Cubic Newton Method received a substantial interest
with the development of inexact, stochastic, and adaptive schemes
(see the papers of \cite{cartis2011adaptive1,grapiglia2017regularized,doikov2018randomized,cartis2018global,hanzely2020stochastic}).
In all of these methods, it is required to solve the subproblem with
our structure~\eqref{MainProblem}.
The Fast Gradient Methods with restarts that have the rate of convergence $O(k^{-6})$
for this problem were constructed by
\cite{roulet2017sharpness} and by \cite{nesterov2019inexact}
with applications to the Cubic Newton.
In our paper we justify that the latter rate is \textit{optimal}, i.e.
the best that can be achieved for the considered class of functions.
Note that in the methods of~\cite{grapiglia2017regularized,grapiglia2019accelerated},
it is also required to minimize the quadratic function
alongside the regularizer with \textit{arbitrary} $p \in [2, 3]$.

In a recent paper by~\cite{nesterov2022quartic},
there were considered new
second-order schemes based on regularization of \textit{degree} $p = 4$.
Then, the rate of the Fast Gradient Method for solving the corresponding subproblem is
$O(k^{-4})$. We prove that this rate is \textit{optimal}, by justifying the corresponding lower bound.

Finally, let us mention the framework of High-Order Proximal-Point Methods
(\cite{nesterov2021inexact})
that became influential for the development of \textit{super-fast} second-order schemes
(see also the works of \cite{kamzolov2020near,nesterov2021inexact2}).
By generalizing the construction of the classical Proximal-Point
Algorithm (\cite{rockafellar1976monotone})
onto arbitrary order $p \geq 2$, the methods assume
that the subproblem of the form~\eqref{MainProblem} can be solved
somehow at each iteration. Thus, it becomes an important question
to study the complexity of problems with arbitrary degree of regularization.

\smallskip

We are interested to analyse the worst-case behaviour
of the first-order algorithms on such problem classes.
Each optimization algorithm $\mathcal{A}$ can be associated
with a sequence of mappings
$$
\begin{array}{rcl}
\mathcal{A} & = & (A_0, A_1, A_2, \dots).
\end{array}
$$
At each iteration $k \geq 0$, mapping $A_k$ takes as input
\textit{oracle information} from the previous $k$ points,
and returns the next point $x_{k+1}$.
In other words, $x_1 := A_0(\varnothing)$, and
$$
\boxed{
\begin{array}{rcl}
x_{k + 1} & := & A_k( \mathcal{I}(x_1), \dots, \mathcal{I}(x_{k}) ), \qquad k \geq 1.
\end{array}
}
$$
In this paper, we consider so called \textit{first-order local oracle} of $f$, that is
$$
\begin{array}{rcl}
\mathcal{I}(x) & = & \{ f(x), \nabla f(x) \}.
\end{array}
$$
The parameters of the regularizer are directly available for the method.

When we provide a lower bound for the rate of convergence,
for a fixed method
we need to build an \textit{example of problem} from our class
such that the residual after the first $k$ iterations is bounded from below:
$$
\begin{array}{rcl}
F(x_k) - F^{*} & \geq & R_k,
\end{array}
$$
by some quantity $R_k$ that is called the \textit{risk}.
The inverse of it, i.e. the minimum number of iterations
required to have an $\varepsilon$-precision for the solution
is called the \textit{complexity}:
$$
\begin{array}{rcl}
C_{\varepsilon} & = & \min \bigl\{ k \; : \; R_k < \varepsilon  \bigr\}.
\end{array}
$$

In the following table, we list the new lower bounds that we prove in our work
alongside the already known ones
for the risk and complexity.
The initial distance to the solution is denoted by~$D$.
Numerical constants are hidden.

\begin{center}
	{
		\small
		\centering
		\renewcommand{\arraystretch}{1.5}
		\begin{tabular}{ | l | c| c | c | }
			\multicolumn{4}{ c  }{
				\begin{tabular}{c}
					{ Euclidean norms }
					\\[7pt]
				\end{tabular}
			}
			\\
			\hline
			\textbf{ \scriptsize Specification } &
			 \textbf{ \scriptsize Risk, }  $R_k$  &
			 \textbf{ \scriptsize Complexity, }  $C_{\varepsilon}$  &
			 \textbf{ \scriptsize Reference}
			\\
			\hline
			\begin{tabular}{@{}l@{}}
				\vspace*{10pt}
				$\;\, p > 1 + \nu$ \\[7pt]
			\end{tabular}
			 &
			$ \displaystyle \!\!\!\!\!
				\biggl( \, \frac{H_{\nu}}{\sigma^{\frac{1 + \nu}{p}} k^{\frac{1 + 3\nu}{2}} }
				\, \biggr)^{\!\!\frac{p}{p - 1 - \nu}} \!\!\!\!
			$
			&
			$
			\displaystyle \!\!
			 \biggl(
			 \frac{H_{\nu}}{\sigma^{\frac{1 + \nu}{p} } \varepsilon^{\frac{p - 1 - \nu}{p}} }
			 \biggr)^{\!\!\frac{2}{1 + 3\nu}} \!\!\!
			$
			&
			{ \scriptsize
			\textbf{Our result:}
			Theorem~\ref{TheoremMain} } \\
			\hline
			\begin{tabular}{@{}l@{}}
				\vspace*{10pt}
				\;{\scriptsize Lipschitz gradient}
				$(\nu = 1, p > 2)$ \\[7pt]
			\end{tabular}
			&
			$ \displaystyle \!\!\!\!\!
			\biggl( \, \frac{H_{1}}{\sigma^{2 / p} k^{2} }
			\, \biggr)^{\! \frac{p}{p - 2}} \!\!\!\!
			$
			&
			$ \displaystyle
			\sqrt{ \frac{H_1}{\sigma^{2/p} \varepsilon^{(p - 2) / p}} }
			$
			&
			{\scriptsize $\triangle$} \\
			\hline

			\begin{tabular}{@{}l@{}}
				\vspace*{10pt}
				\;{\scriptsize Bounded gradient}
				$(\nu = 0, p \geq 2)$ \\[7pt]
			\end{tabular}
			&
			$ \displaystyle \!\!\!\!\!
			\biggl( \, \frac{H_{0}}{\sigma^{1 / p} k^{1/2} }
			\, \biggr)^{\! \frac{p}{p - 1}} \!\!\!\!
			$
			&
			$ \displaystyle \!\!\!\!\!
			\biggl( \,  \frac{H_0}{\sigma^{1/p} \varepsilon^{(p - 1) / p}}   \, \biggr)^{\! 2}
			$
			&
			{ \scriptsize
			\cite{juditsky2014deterministic} } \\
			\hline
			\begin{tabular}{@{}l@{}}
				{ \scriptsize Strongly convex functions with }\\[-3pt]
				{ \scriptsize Lipschitz gradient } $(\nu = 1, p = 2)$
			\end{tabular} &
			$ \sigma D^2 \exp\Bigl( -k \sqrt{\frac{\sigma}{H_1}} \Bigr)$
			&
			$ \sqrt{\frac{H_1}{\sigma}} \log \Bigl( \frac{\sigma D^2}{\varepsilon} \Bigr)$
			&
			{ \scriptsize \cite{nesterov2018lectures} } \\
			\hline
			\begin{tabular}{@{}l@{}}
							\vspace*{5pt}
			 { \scriptsize No composite part } ($\sigma = 0$) \\[10pt]
			 \end{tabular}
			  &
			 $ \displaystyle  \frac{H_{\nu}D^{1 + \nu}}{k^{(3\nu + 1) / 2}} $
			  &
			$\displaystyle  \Bigl( \frac{H_\nu D^{1 + \nu}}{\varepsilon} \Bigr)^{\! \frac{2}{3 \nu + 1} }$
			& \!\!\!\!\!
			{ \scriptsize \cite{guzman2015lower}  }
			\!\!\!\!\!
			\\
			\hline
	\end{tabular}
}
\end{center}

Note that the complexity lower bounds of the same order for similar problem classes
appeared for the first time in the work of \cite{nemirovskii1985optimal} (see their equation (1.21)).
However, it seems that the corresponding proof was never published for a wide audience.

Lower bounds for the cubic subproblem with a (possibly nonconvex) quadratic smooth part
were constructed by \cite{carmon2018analysis}. Comparing with that work,
we were interested in dependence of the complexity
on the \textit{free} parameter $\sigma$ in our sublinear rates.
Thus, we admit arbitrary choice of the regularization constant.

Besides, one significant novelty of our analysis is the use of the \textit{composite} formulation
for~\eqref{MainProblem} that provides a complete flexibility for the regularizer
and smoothness parameters.
Exploring this possibility, we also construct lower bounds for the regularization
by non-Euclidean norms (see the table below).
The employing of arbitrary norms as regularizers for the methods with Taylor's polynomials of different
order was considered in a recent paper by \cite{gratton2021adaptive}.
One step of their second-order algorithm requires to have an inexact solution
to the problem of our form with $p = 3$.

\begin{center}
	{
		\small
		\centering
		\renewcommand{\arraystretch}{1.5}
		\begin{tabular}{ | l | c| c | c | }
			\multicolumn{4}{ c  }{
				\begin{tabular}{c}
					{ Regularization by $\|\cdot\|_q$-norm, $q \geq 1$ }
					\\[7pt]
				\end{tabular}
			}
			\\
			\hline
			\textbf{ \scriptsize Specification } &
			\textbf{ \scriptsize Risk, }  $R_k$  &
			\textbf{ \scriptsize Complexity, }  $C_{\varepsilon}$  &
			\textbf{ \scriptsize Reference}
			\\
			\hline
			\begin{tabular}{@{}l@{}}
				\vspace*{10pt}
				$\;\, p > 1 + \nu$ \\[7pt]
			\end{tabular}
			&
			$ \displaystyle \!\!
			\biggl( \, \frac{H_{\nu}}{\sigma^{\frac{1 + \nu}{p}} k^{\frac{1 + \nu + \nu q}{q}} }
			\, \biggr)^{\!\!\frac{p}{p - 1 - \nu}} \!\!\!\!
			$
			&
			$
			\displaystyle \!\!
			\biggl(
			\frac{H_{\nu}}{\sigma^{\frac{1 + \nu}{p} } \varepsilon^{\frac{p - 1 - \nu}{p}} }
			\biggr)^{\!\!\frac{q}{1 + \nu + \nu q}} \!\!\!
			$
			&
			{ \scriptsize
				\textbf{Our result:}
				Theorem~\ref{TheoremGenMain} } \\
			\hline
			\begin{tabular}{@{}l@{}}
				\vspace*{10pt}
				\;{\scriptsize Lipschitz gradient}
				$(\nu = 1, p > 2)$ \\[7pt]
			\end{tabular}
			&
			$ \displaystyle \!\!\!\!\!
			\biggl( \, \frac{H_{1}}{\sigma^{2 / p} k^{ (2 + q) / q } }
			\, \biggr)^{\! \frac{p}{p - 2}} \!\!\!\!
			$
			&
			$ \displaystyle
			\biggl( \,
			\frac{H_1}{\sigma^{2 / p} \varepsilon^{(p - 2) / p}}
			\,\biggr)^{\! \frac{q}{2 + q}}
			$
			&
			{\scriptsize $\triangle$} \\
			\hline

			\begin{tabular}{@{}l@{}}
				\vspace*{10pt}
				\;{\scriptsize Bounded gradient}
				$(\nu = 0, p \geq 2)$ \\[7pt]
			\end{tabular}
			&
			$ \displaystyle \!\!\!\!\!
			\biggl( \, \frac{H_{0}}{\sigma^{1 / p} k^{1/q} }
			\, \biggr)^{\! \frac{p}{p - 1}} \!\!\!\!
			$
			&
			$ \displaystyle \!\!\!\!\!
			\biggl( \, \frac{H_{0}}{\sigma^{1 / p} \varepsilon^{(p - 1) / p} }
			\, \biggr)^{\! q} \!\!\!\!
			$
			&
			{\scriptsize $\triangle$} \\
			\hline
		\end{tabular}
	}
\end{center}

\bigskip

It appears that for $\ell_1$-norm we obtain the best lower bounds for the rate of convergence.
However, it is not clear how far these estimates from
the upper bounds that can be achieved by some optimization schemes.
At the same time, the rates for the Euclidean case
are tight. They can be reached by the Fast Gradient Methods with restarts
(see, e.g. the monograph of \cite{d2021acceleration}).

The rest of the paper is organized as follows.
In Section~\ref{SectionSmoothing} we review the basic properties of the local smoothing.
We use them in our analysis.
Section~\ref{SectionBound} contains the proof of our lower bounds for the Euclidean norm,
and in Section~\ref{SectionNonEuclidean} we discuss generalizations to non-Euclidean norms.
Section~\ref{SectionDiscussion} contains some final remarks.

\section{Local Smoothing}
\label{SectionSmoothing}

In this section, we summarise well known facts about the local smoothing,
that is the core component of the worst-case example.

Let $g: \mathbb{R}^n \to \mathbb{R}$ be a convex 1-Lipschitz functions, not-necessary differentiable.
Thus
$$
\begin{array}{rcl}
g(x) - g(y) & \leq & \|x - y\|, \qquad \forall x, y \in \mathbb{R}^n.
\end{array}
$$
For a fixed parameter $\mu > 0$, we denote the local smoothing of $g$ by
$$
\begin{array}{rcl}
\mathcal{S}_{\mu}[g](x) & \stackrel{\mathrm{def}}{=} &
\min\limits_{y \in \mathbb{R}^n} \Bigl\{  g(y) + \frac{\mu}{2}\|y - x\|^2  \Bigr\},
\qquad x \in \mathbb{R}^n.
\end{array}
$$
Let us review the basic properties of the function $G(x) := S_{\mu}[g](x)$.
The proofs can be found in the Appendix.

\begin{lemma} \label{LemmaLipschitz}
Function $G$ is convex and differentiable. It is $1$-Lipschitz,
and its gradient is $\mu$-Lipschitz continuous. Thus, for all $x, y \in \mathbb{R}^n$:
\begin{equation} \label{LipFuncSmooth}
\begin{array}{rcl}
G(x) - G(y) & \leq & \|x - y\|,
\end{array}
\end{equation}
and
\begin{equation} \label{LipGradSmooth}
\begin{array}{rcl}
\| \nabla G(x) - \nabla G(y) \| & \leq & \mu \|x - y \|.
\end{array}
\end{equation}
\end{lemma}

As a direct consequence of this lemma, we can conclude that $G$ has
a H\"older continuous gradient for \textit{any} $\nu \in [0, 1]$ with
constant
\begin{equation} \label{H_nu_constant}
\begin{array}{rcl}
H_{\nu}(G) & = & 2^{1 - \nu} \mu^{\nu}.
\end{array}
\end{equation}
Indeed, \eqref{LipFuncSmooth} implies that the gradients are bounded:
$\| \nabla G(x) \| \leq 1$, $\forall x \in \mathbb{R}^n$. Hence,
$$
\begin{array}{rcl}
\| \nabla G(x) - \nabla G(y) \|
& = &
\| \nabla G(x) - \nabla G(y) \|^{\nu} \cdot \| \nabla G(x) - \nabla G(y) \|^{1 - \nu} \\
\\
& \leq &
\mu^{\nu} \|x - y\|^{\nu} \cdot 2^{1 - \nu}
\;\; = \;\; H_{\nu}(G)\|x - y\|^{\nu}, \qquad \forall x, y \in \mathbb{R}^n.
\end{array}
$$

The following lemma shows that $G$ is also close to $g$.

\begin{lemma} \label{LemmaClose}
For any $x \in \mathbb{R}^n$, it holds
\begin{equation} \label{G_close_to_g}
\begin{array}{rcl}
g(x) & \geq & G(x) \;\; \geq \;\; g(x) - \frac{1}{2\mu}.
\end{array}
\end{equation}
\end{lemma}

We see that parameter $\mu$ provides a trade-off between the quality of approximation
and the level of smoothness for the new function.
The effect of varying this parameter is shown in Figure~\ref{Fig:Smoothing}.

\begin{figure}[h]
	\centering
	\includegraphics[width=0.50\textwidth ]{smoothing_1.pdf}
	\caption{ \small Local smoothing of a piecewise linear function.}
	\label{Fig:Smoothing}
\end{figure}

Finally, we state the following important property that $G$ depends on $g$
in \textit{a local way}:
the value and the gradient of $G$ at $x$ depend only on the restriction
$g|_{U_{x, \mu}}$ of $g$ onto the ball
$$
\begin{array}{rcl}
U_{x, \mu} & \stackrel{\mathrm{def}}{=} &
\Bigl\{
y \in \mathbb{R}^n \; : \;
\|y - x\| \leq \frac{2}{\mu}
\Bigr\}.
\end{array}
$$
More formally, the following proposition holds.
\begin{lemma} \label{LemmaLocal}
Let $z := \mathop{\rm argmin}\limits_{y \in \mathbb{R}^n}
\Bigl\{ \; g(y) + \frac{\mu}{2}\|y - x\|^2 \; \Bigr\}$
for some fixed $x \in \mathbb{R}^n$ and $\mu > 0$.
Then
$$
\begin{array}{rcl}
x - \frac{1}{\mu}\nabla G(x)
& \equiv & z \;\; = \;\;
\mathop{\rm argmin}\limits_{y \in {\rm int \,} U_{x, \mu}}
\Bigl\{
\;
g(y) + \frac{\mu}{2}\|y - x\|^2
\;
\Bigr\}.
\end{array}
$$
\end{lemma}

Therefore, in order to compute the oracle information
$\{ G(x), \nabla G(x) \}$ at some given point $x$,
it is enough to have an access only to the \textit{neighbourhood}
of function $g$ around this point.
This fact is crucial for proving the lower bound.

\section{Lower Complexity Bound}
\label{SectionBound}

Now, employing the construction introduced in \cite{guzman2015lower},
we are going to provide a resisting oracle strategy to establish the
lower complexity bound for problem~\eqref{MainProblem}.

Let us fix arbitrary
$\nu \in [0, 1]$ such that $\nu < p - 1$, and some constant $H_{\nu} > 0$.
We prove the following result.

\begin{theorem} \label{TheoremMain}
For any $T$-step algorithm with $T \leq n$,
there exists a convex differentiable function $f(\cdot)$
whose gradient is H\"older continuous of degree $\nu$ with constant $H_{\nu}$, such that
for problem~\eqref{MainProblem}, we have

\begin{equation} \label{MainResult}
\begin{array}{rcl}
F(x_T) - F^{*}
& \geq &
\bigl(  \frac{p - 1}{p} \bigr)^{\frac{(p - 1)(1 + \nu)}{p - 1 - \nu}}
 \bigl( \frac{1}{2} \bigr)^{ \frac{(2p - 1)(1 + \nu)}{p - 1 - \nu} } \cdot
\biggl(
\frac{ H_{\nu} }{  \sigma^{\frac{1 + \nu}{p} } T^{ \frac{1 + 3\nu}{2} }}
\biggr)^{\frac{p}{p - 1 - \nu}}.
\end{array}
\end{equation}
The size of the solution to that problem is bounded as
\begin{equation} \label{SolBound}
\begin{array}{rcl}
\|x^{*}\| & \leq &
\bigl(  3(p - 1) 2^{p - 3}  \bigr)^{\frac{1}{p}}
\cdot
\Bigl(\frac{1}{2} \bigl( \frac{p - 1}{4p} \bigr)^{\nu} \Bigr)^{\frac{1}{p - 1 - \nu}}
\biggl(
\frac{H_{\nu}}{\sigma T^{\frac{1 + 3\nu}{2}}}
\biggr)^{\frac{1}{p - 1 - \nu}}.
\end{array}
\end{equation}
\end{theorem}
\proof

\noindent
Let us fix some positive $\delta$,  $\mu$, and $\beta$.

Resisting oracle chooses a set of numbers $\xi_1, \dots, \xi_T \in \{-1, 1\}$
and a permutation $k \mapsto \alpha(k) \in \{1, 2, \dots, T\}$.

Having these parameters, we consider the family of convex functions, $x \in \mathbb{R}^n$:
$$
\begin{array}{rcl}
g_t(x) & = &
\max\limits_{1 \leq k \leq t} \Bigl[ \; \xi_k \langle e_{\alpha(k) }, x\rangle \, - \, (k - 1) \delta  \; \Bigr],
\qquad 1 \leq t \leq T,
\end{array}
$$
where $e_i$ is the standard basis vector in $\mathbb{R}^n$ (recall that we assume that $n \geq T$).
It is clear that all $g_t(\cdot)$ are Lipschitz continuous with constant $1$.

Then, we are going to take as our final objective the function
$$
\boxed{
	\begin{array}{rcl}
	f(x) & := & \beta \mathcal{S}_{\mu}[ g_T ] (x)
	\end{array}
}
$$

Let us estimate the minimum $F^{*}$ and the size of the solution $\|x^{*}\|$.
Note that we have, for any $x \in \mathbb{R}^n$
\begin{equation} \label{FUpperBound}
\begin{array}{rcl}
F(x) & = &
\beta \mathcal{S}_{\mu}[g_T](x) + \frac{\sigma}{p}\|x\|^p
\;\; \overset{\eqref{G_close_to_g}}{\leq} \;\;
\beta g_T(x) + \frac{\sigma}{p}\|x\|^p \\
\\
& \leq &
\mathcal{H}(x)
\;\; \stackrel{\mathrm{def}}{=} \;\;
\max\limits_{1 \leq k \leq T} \beta \xi_k \langle e_{\alpha(k)}, x \rangle
+ \frac{\sigma}{p}\|x\|^p.
\end{array}
\end{equation}
At the same time, we can also bound $F$ by $\mathcal{H}$ from below,
as follows
\begin{equation} \label{FLowerBound}
\begin{array}{rcl}
F(x) & \overset{\eqref{G_close_to_g}}{\geq} &
\beta g_T(x) + \frac{\sigma}{p}\|x\|^p - \frac{\beta}{2\mu}
\;\; \geq \;\;
\mathcal{H}(x) - \frac{\beta}{2\mu} -  \beta (T - 1) \delta.
\end{array}
\end{equation}
We can compute the minimum of $\mathcal{H}(\cdot)$
by using the symmetry within the problem,
\begin{equation} \label{MinUpperBound}
\begin{array}{rcl}
F^{*} & \overset{\eqref{FUpperBound}}{\leq} &
\mathcal{H}^* \;\; = \;\; \min\limits_{x \in \mathbb{R}^n} \mathcal{H}(x)
\;\; = \;\;
\min\limits_{\gamma > 0} \Bigl\{  -\beta \gamma
+ \frac{\sigma}{p} T^{\frac{p}{2}} \gamma^{p}  \Bigr\} \\
\\
& = &
- \frac{p - 1}{p} \cdot \Bigl(
\frac{\beta^p}{\sigma T^{p / 2}}
\Bigr)^{\frac{1}{p - 1}}.
\end{array}
\end{equation}
From the uniform convexity of the regularizer (see, e.g. Lemma 2.5 in \cite{doikov2021minimizing})
we can bound the size of the solution,
\begin{equation} \label{SolSize}
\begin{array}{rcl}
\bigl( \frac{1}{2} \bigr)^{p - 2}  \frac{\sigma}{p}\|x^*\|^p
& \leq &
F(0) - F^{*}
\;\; \overset{\eqref{FUpperBound}, \eqref{FLowerBound} }{\leq} \;\;
\mathcal{H}(0) - \mathcal{H}^{*}
+ \frac{\beta}{2\mu} +  \beta (T - 1) \delta \\
\\
& = &
\frac{p - 1}{p} \cdot \Bigl(
\frac{\beta^p}{\sigma T^{p / 2}}
\Bigr)^{\frac{1}{p - 1}}
+ \frac{\beta}{2\mu} +  \beta (T - 1) \delta.
\end{array}
\end{equation}

Now, let us present a strategy for choosing $\xi_k$ and $\alpha(k)$.

\begin{itemize}
	\item At first step, the algorithm returns point $x_1$,
	which does not depend on the objective.
	Let us pick
	$$
	\begin{array}{rcl}
	\alpha(1) & \in & \mathop{\rm Argmax}\limits_{1 \leq k \leq T} |\langle e_k,  x_1 \rangle |.
	\end{array}
	$$
	In other words, $\alpha(1)$ is an index of
	a maximal element (in absolute value) among first $T$ coordinates of $x_1$.
	Then, we specify $\xi_1 \in \{-1, 1\}$ in a way that
	$$
	\begin{array}{rcl}
	\xi_1 \langle e_{\alpha(1)}, x_1 \rangle & = &
	|\langle e_{\alpha(1)}, x_1 \rangle|,
	\end{array}
	$$
	hence $\xi_1 =  {\rm sign \,}(\langle e_{\alpha(1)}, x_1 \rangle)$.
	So,
	$$
	\begin{array}{rcl}
	g_1(x) & = & \xi_1 \langle e_{\alpha(1)}, x \rangle.
	\end{array}
	$$

	\item At step $2 \leq t \leq T$, assume that we have built function $g_{t - 1}(x)$.
	Let $x_t$ be the point of the trajectory of $\mathcal{A}$ at iteration $t - 1$,
	applied to the current objective:
	$$
	\begin{array}{rcl}
	x_t & := & A_{t - 1}( \mathcal{I}_{t - 1}(x_1), \dots, \mathcal{I}_{t - 1}(x_{t - 1})),
	\end{array}
	$$
	where $\mathcal{I}_{t - 1}(x) \stackrel{\mathrm{def}}{=} \{ \beta \mathcal{S}_{\mu}[g_{t - 1}](x)  ,
	\beta \nabla \mathcal{S}_{\mu}[g_{t - 1}](x)   \}$.
	Let us choose as $\alpha(t)$ the index of a maximal element of $x_t$ (in absolute value) among
	first $T$ coordinates, except $\alpha(1), \dots, \alpha(t - 1)$. Thus,
	$$
	\begin{array}{rcl}
	\alpha(t) & \in &
	\mathop{\rm Argmax}\limits_{ \stackrel{1 \leq k \leq T}{k \notin \{\alpha(1), \dots, \alpha(t - 1) \}}}
	|\langle e_k, x_t \rangle|.
	\end{array}
	$$
	Then, we specify $\xi_t \in \{-1, 1\}$ such that
	$$
	\begin{array}{rcl}
	\xi_t \langle e_{\alpha(t)}, x_t \rangle & = & | \langle e_{\alpha(t)}, x_t \rangle |,
	\end{array}
	$$
	hence $\xi_t = {\rm sign \,}(\langle e_{\alpha(t)}, x_t \rangle)$, and thus we obtain the next $g_t(x)$.
\end{itemize}

We need to prove that for any $2 \leq t \leq T$, function $g_t(\cdot)$ coincides with $g_s(\cdot)$,
$1 \leq s < t$,
in the ball of radius $\frac{2}{\mu}$ around $x_s$:
\begin{equation} \label{LocalEquiv}
\begin{array}{rcl}
g_t(x) & = & g_s(x), \qquad x \in U_{x_s, \mu}
\;\; \stackrel{\mathrm{def}}{=} \;\;
\Bigl\{ x \in \mathbb{R}^n \; : \; \|x - x_s\| \leq \frac{2}{\mu}  \Bigr\},
\end{array}
\end{equation}
which means that $g_t(\cdot)$ is \textit{indistinguishable} from $g_s(\cdot)$ during the first $s$ steps,
and
hence by Lemma~\ref{LemmaLocal} their smoothings
$\mathcal{S}_{\mu}[g_t]$ and $\mathcal{S}_{\mu}[g_s]$ are also indistinguishable for the algorithm.

Indeed,
$$
\begin{array}{rcl}
g_t(x) & = &
\max\Bigl\{  g_s(x), \;
\max\limits_{s < k \leq t}\Bigl[\;  \xi_k\langle e_{\alpha(k)}, x \rangle \, - \, (k - 1) \delta    \;\Bigr]
\Bigr\}.
\end{array}
$$
By the definition of $g_s(x)$, we have
$$
\begin{array}{rcl}
g_s(x_s) & \geq & \Bigl[ \xi_k \langle e_{\alpha(k)}, x_s \rangle - (k - 1)\delta \Bigr] + \delta,
\qquad s < k \leq t.
\end{array}
$$
Hence, due to the Lipschitz continuity, it holds for all $x$ such that
$\|x - x_s\| \leq \frac{\delta}{2}$:
$$
\begin{array}{rcl}
g_s(x) & \geq & \Bigl[ \xi_k \langle e_{\alpha(k)}, x \rangle - (k - 1) \delta \Bigr],
\qquad s < k \leq t.
\end{array}
$$
Choosing
\begin{equation} \label{DeltaChoice}
\boxed{
	\begin{array}{rcl}
	\delta := \frac{4}{\mu}
	\end{array}
}
\end{equation}
we conclude that~\eqref{LocalEquiv} is true.

Thus, we established correctness for the construction of the resisting oracle. Namely, it holds,
for all $s \leq t \leq T$:
$$
\begin{array}{rcl}
 \mathcal{I}_{t}(x_s)
    & = &
 \mathcal{I}_T(x_s)
 \;\; \equiv \;\;
 \mathcal{I}(x_s),
\end{array}
$$
so the oracles $\mathcal{I}_{t}(\cdot)$ and $\mathcal{I}(\cdot)$
are identical along the trajectory of the method.

It remains to bound from below
the residual in the function value for the last iteration of the method.
From simple observations, we notice that
\begin{equation} \label{GTLastIter}
\begin{array}{rcl}
g_T(x_T) & \geq & | x_T^{(\alpha(T))} | - (T - 1) \delta
\;\; \geq \;\; - (T - 1) \delta
\;\; \overset{\eqref{DeltaChoice}}{=} \;\; -(T - 1) \frac{4}{\mu}.
\end{array}
\end{equation}
Hence, for the last function value, we have
\begin{equation} \label{FTLastIter}
\begin{array}{rcl}
F(x_T) & \geq & f(x_T)
\;\; = \;\;
\beta \mathcal{S}_{\mu}[g_T](x_T)
\;\; \overset{ \eqref{G_close_to_g}, \eqref{GTLastIter}}{\geq} \;\;
-(T - 1) \frac{4\beta}{\mu} - \frac{\beta}{2 \mu}
\;\; \geq \;\;
-\frac{4\beta T }{\mu}.
\end{array}
\end{equation}
Therefore, the residual can be bounded as follows,
\begin{equation} \label{ResBound}
\begin{array}{rcl}
F(x_T) - F^{*}
& \overset{\eqref{FTLastIter}, \eqref{MinUpperBound}}{\geq} &
\frac{p - 1}{p} \cdot \Bigl(
\frac{\beta^p}{\sigma T^{p / 2}}
\Bigr)^{\frac{1}{p - 1}}  - \frac{4\beta T }{\mu}
\;\; = \;\;
\frac{p - 1}{2p} \cdot \Bigl(
\frac{\beta^p}{\sigma T^{p / 2}}
\Bigr)^{\frac{1}{p - 1}},
\end{array}
\end{equation}
where the last equality holds by the following choice of the smoothing parameter:
\begin{equation} \label{MuChoice}
\boxed{
\begin{array}{rcl}
\mu & := &
\frac{8p}{p - 1} \Bigl(
\frac{\sigma T^{ (3p - 2) / 2 }}{\beta}
\Bigr)^{\frac{1}{p - 1}}
\end{array}
}
\end{equation}
The only parameter which remains to determine is $\beta$.
We know that $f$ must have H\"older continuous gradient with constant
$$
\begin{array}{rcl}
H_{\nu} & \overset{\eqref{H_nu_constant}}{=} &
\beta 2^{1 - \nu} \mu^{\nu}
\;\; \overset{\eqref{MuChoice}}{=} \;\;
2 \Bigl(  \frac{4p}{p - 1}  \Bigr)^{\nu}
\Bigl(  \sigma T^{(3p - 2) / 2}  \Bigr)^{\frac{\nu}{p - 1}}
\beta^{\frac{p - 1 - \nu}{p - 1}}.
\end{array}
$$
Therefore, we get the following value for the last free parameter:
\begin{equation} \label{BetaChoice}
\boxed{
\begin{array}{rcl}
\beta & := &
\Bigl(
\frac{1}{2} \bigl(  \frac{p - 1}{4p} \bigr)^{\nu} H_{\nu}
\Bigr)^{\frac{p - 1}{p - 1 - \nu}}
\Bigl(  \sigma T^{(3p - 2) / 2} \Bigr)^{- \frac{\nu}{p - 1 - \nu}}
\end{array}
}
\end{equation}
Substituting this value into~\eqref{ResBound},
we obtain the final lower bound on the convergence rate~\eqref{MainResult}.

The bound~\eqref{SolBound} for the size of the solution follows immediately
from
$$
\begin{array}{rcl}
\|x^{*}\| & \overset{\eqref{SolSize}}{\leq} &
\biggl[ \frac{p 2^{p - 2}}{\sigma}
 \biggl(
 \frac{p - 1}{p} \cdot \Bigl(  \frac{\beta^p}{\sigma T^{p / 2}} \Bigr)^{\frac{1}{p - 1}}
 + \frac{\beta}{2\mu} + \beta(T - 1) \delta
 \biggr)
\biggr]^{1 / p} \\
\\
& \overset{\eqref{DeltaChoice}}{\leq} &
\biggl[ \frac{p 2^{p - 2}}{\sigma}
\biggl(
\frac{p - 1}{p} \cdot \Bigl(  \frac{\beta^p}{\sigma T^{p / 2}} \Bigr)^{\frac{1}{p - 1}}
+ \frac{4\beta T}{\mu}
\biggr)
\biggr]^{1 / p} \;\; \overset{\eqref{MuChoice}}{=} \;\;
\biggl[ \frac{3 (p - 1) 2^{p - 3}}{\sigma}
\cdot \Bigl(  \frac{\beta^p}{\sigma T^{p / 2}} \Bigr)^{\frac{1}{p - 1}}
\biggr]^{1 / p},
\end{array}
$$
by substituting the expression for $\beta$ \eqref{BetaChoice}.
\hfill $\Box$ \par \noindent \medskip

Let us consider the case $\nu = 1$. Then, Theorem \ref{TheoremMain}
states that for any $T$-step algorithm there is an instance of problem
with bounded solution,
\begin{equation} \label{ExplSolBound}
\begin{array}{rcl}
\| x^{*} \| & \overset{\eqref{SolBound}}{\leq} & O \Bigl(
\bigl( \frac{H_1}{\sigma T^2} \bigr)^p  \Bigr),
\end{array}
\end{equation}
such that the residual is bounded from below,
\begin{equation} \label{ExpLower}
\begin{array}{rcl}
F(x_T) - F^{*} & \geq &
O\Bigl(
\bigl(  \frac{H_1}{\sigma^{2 / p} T^2}  \bigr)^{\frac{p}{p - 2}}
\Bigr).
\end{array}
\end{equation}
At the same time, the standard rate of the composite Fast Gradient Method~\cite{nesterov2013gradient}
starting from $x_1 = 0$ is
$$
\begin{array}{rcl}
F(x_T) - F^{*} & \leq & O\Bigl( \frac{H_1 \|x^{*}\|^2}{T^2}  \Bigr)
\;\; \overset{\eqref{ExplSolBound}}{\leq} \;\;
O\Bigl( \bigl( \frac{H_1}{\sigma^{2/p} T^{2}} \bigr)^{\frac{p}{ p - 2}}  \Bigr),
\end{array}
$$
that matches the lower bound~\eqref{ExpLower}.

\newpage
\section{Regularization with Non-Euclidean Norms}
\label{SectionNonEuclidean}

In this section, we study a generalization
of problem~\eqref{MainProblem} to non-Euclidean norms.
For arbitrary $q \geq 1$, we denote by $\| \cdot \|_q$
the standard $\ell_q$-norm:
$$
\begin{array}{rcl}
\|x\|_q & \stackrel{\mathrm{def}}{=} & \biggl(
\sum\limits_{i = 1}^n
\big| x^{(i)} \big|^{q}  \biggr)^{1/q}.
\end{array}
$$

Now, let us consider the problem:
\begin{equation} \label{GenMainProblem}
\begin{array}{rcl}
\min\limits_{x \in \mathbb{R}^n} \Bigl\{
F(x) & \stackrel{\mathrm{def}}{=} & f(x) + \frac{\sigma}{p}\|x\|^p_q
\Bigr\},
\end{array}
\end{equation}
for some $p \geq 2$, $\sigma > 0$, and convex differentiable $f$
that has H\"older continuous gradient
(w.r.t. the standard Euclidean norm\footnote{Hence, we use
	different norms for the regularizer and for defining the smoothness properties of $f$.
	Generalizations to non-Euclidean norms for $f$ are also possible,
	see~(\cite{guzman2015lower}).}).
As before, we assume that $\nu < p - 1$.
We can prove the following result.

\begin{theorem} \label{TheoremGenMain}
For any $T$-step algorithm with $T \leq n$, there exists
a convex differentiable function $f(\cdot)$ whose gradient
is H\"older continuous of degree $\nu$ with constant $H_{\nu}$,
such that for problem~\eqref{GenMainProblem}, we have
\begin{equation} \label{GenProblemLowerBound}
\begin{array}{rcl}
F(x_T) - F^* & \geq &
\bigl( \frac{p - 1}{p} \bigr)^{\frac{(p - 1)(1 + \nu)}{p - 1 - \nu}}
\bigl( \frac{1}{2} \bigr)^{\frac{(2p - 1)(1 + \nu)}{p - 1 - \nu}}
\cdot
\biggl(
\frac{ H_{\nu} }{  \sigma^{\frac{1 + \nu}{p} } T^{ \frac{1 + \nu + \nu q}{q} }}
\biggr)^{\frac{p}{p - 1 - \nu}}.
\end{array}
\end{equation}
\end{theorem}
\proof
The proof is very similar to that one of Theorem~\ref{TheoremMain}.
We use the same strategy for the resisting oracle, and the
candidate for the smooth part is $f(x) := \beta \mathcal{S}_{\mu}[g_T](x)$,
with $\delta := \frac{4}{\mu}$ (see definition of $g_T(\cdot)$ in the previous theorem).

For the minimum, we have the following estimate,
$$
\begin{array}{rcl}
F^{*} & = &  \min\limits_{x \in \mathbb{R}^n} F(x)
\;\; \leq \;\;   \min\limits_{x \in \mathbb{R}^n}
\biggl\{ \;  \max\limits_{1 \leq k \leq T} \beta \xi_k \langle e_{\alpha(k)}, x \rangle
+ \frac{\sigma}{p}\|x\|_q^p
\; \biggr\} \\
\\
& = &
\min\limits_{\gamma > 0}
\Bigl\{   -\beta \gamma + \frac{\sigma}{p} T^{\frac{p}{q}} \gamma^p  \Bigr\}
\;\; = \;\;
- \frac{p - 1}{p} \cdot \Bigl(
\frac{\beta^p}{\sigma T^{p / q}}
\Bigr)^{\frac{1}{p - 1}}.
\end{array}
$$
Thus, we get the bound for the residual,
\begin{equation} \label{GenMainResBound}
\begin{array}{rcl}
F(x_T) - F^{*} & \geq &
\frac{p - 1}{p} \cdot  \Bigl(
\frac{\beta^p}{\sigma T^{p / q}}
\Bigr)^{\frac{1}{p - 1}} - \frac{4\beta T}{\mu}
\;\; = \;\;
\frac{p - 1}{2p} \cdot \Bigl( \frac{\beta^p}{\sigma T^{p / q}} \Bigr)^{ \frac{1}{p - 1}},
\end{array}
\end{equation}
where we made the following choice of the smoothing parameter to balance the two terms,
$$
\boxed{
\begin{array}{rcl}
\mu & := & \frac{8p}{p - 1} \Bigl(
\frac{\sigma T^{\frac{p + qp - q}{q}}}{\beta}
\Bigr)^{\frac{1}{p - 1}}
\end{array}
}
$$
Substituting this expression into the equation for the constant of H\"older continuity,
$H_{\nu} = \beta 2^{1 - \nu} \mu^{\nu}$, we get the value for the last unknown parameter,
$$
\boxed{
\begin{array}{rcl}
\beta & := &
\Bigl(
\frac{1}{2} \bigl(  \frac{p - 1}{4p}  \bigr)^{\nu} H_{\nu}
\Bigr)^{\frac{p - 1}{p - 1 - \nu}}
\Bigl( \sigma T^{\frac{p + qp - q}{q}} \Bigr)^{-\frac{\nu}{p - 1 - \nu}}
\end{array}
}
$$
Plugging it into~\eqref{GenMainResBound} completes the proof.
\hfill $\Box$ \par \noindent \medskip

\begin{figure}[h]
	\centering
	\includegraphics[width=0.28\textwidth ]{z1.pdf}
	\includegraphics[width=0.28\textwidth ]{z2.pdf}
	\includegraphics[width=0.28\textwidth ]{z3.pdf}
	\caption{ \small Level sets of $\|\cdot\|_1^p$ for different powers.}
	\label{Fig:Norms}
\end{figure}

\section{Discussion}
\label{SectionDiscussion}

In this paper, we proved new lower bounds on the
oracle complexity for minimizing regularized convex functions
by first-order methods.
As a particular case, we established that the
\textit{best possible} rate of convergence for minimizing a function with Lipschitz continuous gradient
 regularized by \textit{cube} of the Euclidean norm is $O(k^{-6})$,
and for the \textit{fourth power} of the norm we have $O(k^{-4})$.
We know that the Fast Gradient Method achieves these rates.
It remains to be an interesting open question: whether we can construct the lower bounds
with a quadratic function for the smooth part.

Another interesting observation that we discovered is a change in the rate of convergence
for the regularizers with different norms. Thus, by using the third power of $\ell_1$-norm,
we obtain the lower bound of order $O(k^{-9})$.
However, it is not clear how far this estimate from the upper bound
that can be reached by some optimization method.
We keep these questions for further investigation.
At the same time, taking a power of $\ell_1$-norm might preserve its sparsity properties (see Figure~\ref{Fig:Norms}),
which is desirable for applications.

\acks{This project has received funding from the European Research Council (ERC)
	under the European Union's Horizon 2020 research and innovation programme (grant agreement No. 788368).}


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

Since the first planet discoveries in the 1990s \citep{Wolszczan1992,Mayor1995}, more than 1400 planets have been discovered, according to the Exoplanet Archive \citep{Wright2011}.  \textit{Kepler} was launched in 2009 and obtained precise photometric measurements for $\sim160,000$ stars with nearly continuous coverage for four years.  More than 3500 candidates have been discovered via the planet transit method \citep{Borucki2011,Batalha2013,Burke2014}, with nearly 1,000 of them now confirmed \citep[e.g.,][]{Lissauer2014,Rowe2014}. The transit planet search (TPS) algorithm \citep{Jenkins2002,Jenkins2010} has been used to search \textit{Kepler} stellar light curves for the characteristic dips in flux indicative of a planet.  Those which meet a certain significance threshold are placed on the Threshold Crossing Event (TCE) list \citep{Tenenbaum2013,Tenenbaum2014}. The TCEs are further examined and can be upgraded to \textit{Kepler Object of Interest} (KOI) candidate status or downgraded as false positives (FPs), which may be erroneous signals or real, astrophysical phenomena like eclipsing binaries \citep{Prsa2011,Slawson2011}, either on the target star or in the background.

The photometric transit technique can determine the radius of a planet, but generally not the mass and hence does not immediately indicate if a transit signal is due to a planet or a binary star system.  In certain cases, these transiting planet candidates can be confirmed as planets through a statistical elimination of astrophysical false positives, such as eclipsing binaries \citep{Fressin2012}.  Validation by multiplicity has been used on hundreds of candidates in multiple planet candidate systems \citep{Lissauer2014,Rowe2014}.  This method uses the fact that there are very few false positives in multiple transiting planet candidate systems due to the rare occurrence of two independent, unlikely events (either multiple false positives or a false positive and a true planet) both occurring for the same star.  However, neither of these methods can determine the planetary masses or orbital parameters, which is necessary for further characterization of the planets.  To date, the only way to determine masses purely from a photometric light curve is through the detection and inversion of TTVs, deviations from a perfectly linear ephemeris due to gravitational interactions of neighboring planets. (see Section~\ref{sec:ttv}).  In this paper, strong TTVs of the planets orbiting PH3 allow for the determination of planetary masses and thus their confirmation as planets.

\subsection{Planet Hunters}

The Planet Hunters (PH) project\footnote{http://www.planethunters.org/} is one of the many projects within the Zooniverse\footnote{https://www.zooniverse.org/} citizen science program \citep{Lintott2008,Lintott2011,Fortson2012}.  The goal of this project is to assist in the analysis of \textit{Kepler} lightcurves by complementing the TPS algorithmic search with human eyes.  PH breaks \textit{Kepler} lightcurves into 30 day increments and asks users to mark transit-like signals.  Additionally, the Talk\footnote{http://talk.planethunters.org} discussion page allows users to interact with the science team and each other, discuss lightcurves, and even collect certain categories of potential phenomena, such as eclipsing binaries,  heartbeat stars, variable stars, microlensing, and circumbinary planets.  Since December, 2010, more than 22 million \textit{Kepler} 30-day increment lightcurves have been examined by more than 290,000 public volunteers, amounting to 180 years of 40-hour work weeks.

PH has detected more than 60 new planet candidates, including a large number in their host star's habitable zone (HZ)	 \citep{Fischer2012,Lintott2013,Wang2013a,Schmitt2014}, and two newly confirmed planets \citep{Schwamb2013,Wang2013a}.  The first known \textit{Kepler} seven planet candidate \citep{Cabrera2014} was also independently discovered by PH and orbits in a compact, solar system-like arrangement \citep{Schmitt2014}.  PH's main contribution to the population of \textit{Kepler} candidates has been in long period candidates.  This is due to the fact that PH users detect one transit at a time, making PH sensitive to even one and two transit systems, whereas computer algorithms typically require three or more transits and build up signal with an increasing number of transits, thus decreasing sensitivity to longer periods.  However, computer algorithms are more robust at identifying smaller planets with transits that are hard to detect by eye.  As such, most PH candidates are approximately Neptune-sized or larger, which is where PH approaches completeness for short periods \citep{Schwamb2012}.

\subsection{KOI-1353}

KIC 7303287 was first found to have a planet candidate (KOI-1353.01) in \citet{Borucki2011}, an $18.6\pm5.3R_{\oplus}$ planet in a 125.87 day orbit.  Later, a 34.54 day planet candidate (KOI-1353.02) with a radius of $3.8\pm1.1R_{\oplus}$ was also discovered \citep{Batalha2013}.  PH users then noticed an additional 66 day planet candidate, which was later detected as a 330 day TCE with a loosely constrained radius of $4.23\pm3.98$ \citep{Tenenbaum2014}, an alias of the true period.  PH volunteers also discovered its five hour amplitude TTVs in the transit data and brought it to the attention of the science team via the Talk page (see Section~\ref{sec:agol}).  Closer examination of KOI-1353.01 showed that outer planet also exhibits TTVs, but with an amplitude of order minutes.  In addition to a sinusoidal component with period equal to the super period, we also can identify a chopping signal in the residuals of the middle planet, allowing us to break the mass-eccentricity degeneracy.  As described below in Section~\ref{sec:agol}, the inner planet is also confirmed by having a mass upper limit well within the planetary regime.  Therefore, KOI-1353 is confirmed as PH3~b, c, and d.  While in preparation of this paper, \citet{Rowe2014} confirmed the inner and outer planets, PH3~b and d, via ``validation by multiplicity'' based on the statistical assessment that multiple transiting signals are unlikely to be false positives.  They were given the names Kepler-289 b and c.  See Table~\ref{tab:names} for cross-matching the names of the planets.   Furthermore, PH3 d's TTVs were detected at a significant level in \citet{Mazeh2013}, but they were only able to fit the signal to a parabola rather than a full sine curve.  As such, they do not measure the TTV period or its amplitude.  PH3~c's TTVs were undetected because \citet{Mazeh2013} only searched through KOIs. This system has period ratios of 1.91 (both the outer/middle and middle/inner period ratios), which is somewhat far from the 1:2:4 mean motion resonance (MMR).

\begin{deluxetable}{lccc}
\tablewidth{0pt}
\tablecaption{Stellar and planetary designations. \label{tab:names}}
\tablehead{
\colhead{Star} &
\colhead{Inner planet} &
\colhead{Middle planet} &
\colhead{Outer planet} }
\startdata

PH3          & PH3~b         & PH3~c        & PH3~d         \\
Kepler-289   & Kepler-289~b  & \nodata      & Kepler-289~c  \\
KOI-1353     & KOI-1353.02   & \nodata$^*$  & KOI-1353.01   \\

\enddata

\tablecomments{$^*$ Identified as a TCE with a period five times the true 66 day period.}

\end{deluxetable}

\subsection{Transit Timing Variations and Mean Motion Resonances}
\label{sec:ttv}

Since the transit technique provides a planet's radius, a mass measurement can both confirm the planetary nature of the candidate as well as measure the planetary density.  In some systems with multiple planets, TTVs can determine or provide constraints on the planet masses \citep{Miralda2002,Agol2005,Holman2005}.   TTVs of transiting  planets can even be used to indicate the presence of non-transiting planets \citep{Ballard2011,Nesvorny2012}, co-orbital (Trojan) planets \citep{Ford2007}, or exomoons \citep{Kipping2009}.  TTVs were first used to confirm a pair of planet candidates orbiting Kepler-9 \citep{Holman2010}. Since then, numerous other \textit{Kepler} candidates have been confirmed using TTVs \citep[e.g.,][]{Lissauer2011,Steffen2012,Xie2014}.

One potential disadvantage of measuring planet masses with TTVs is that, in some circumstances, the masses become degenerate with the eccentricity \citep{Lithwick2012}.  This occurs when: 1) the planet system is close to, but not within, a first-order mean motion resonance; and 2) the measurement errors are significant enough to only detect the sinusoidal TTV caused by the nearby resonance.   In this case, the amplitude of the TTV depends both on the planet masses and on the free eccentricity, $z_{\rm{free}}$, which is the component of eccentricity that is not driven  by the resonant forcing.  If the $z_{\rm{free}}$ is known, the mass is also known.  \citet{Lithwick2012} examines the degeneracy closely for systems near MMR. In a $z_{\rm{free}}=0$ system, the TTVs of pairs of interacting planets will be anti-correlated.  If the two sets of TTVs are not anti-correlated, this implies that significant free eccentricities exist.  However, the absence of a phase shift does not necessarily mean that free eccentricities are zero \citep{Lithwick2012}.  This could instead be an unlikely moment in time when free eccentricities do exist, but their phases cancel out.  If $z_{\rm{free}}$ is small, then one can calculate an upper limit on the mass by assuming $z_{\rm{free}}=0$ \citep{Lithwick2012}.  Determining $z_{\rm{free}}$ allows one to calculate the true mass.

The mass-$z_{\rm{free}}$ degeneracy can also be broken with more precise data.  An $O-C$ TTV signal for a pair of planets near MMR results in a periodic signal (sinusoidal for cases of low $z_{\rm{free}}$).  This period is equal to:

\begin{center}
\begin{equation}
\label{eq:ttv}
P_{\rm{TTV}}=\frac{1}{j/P_j - (j-1)/P_{j-1}}
\end{equation}
\end{center}

\noindent where the outer and inner planets are in a near $j:j-1$ resonance, respectively, and $j$ is a positive integer \citep{Agol2005,Lithwick2012}.  This period is called the super-period or TTV period.  After one fits the transit midpoints with a linear period and the TTV period, a combination of high quality data and large amplitude TTVs can show the proportionally smaller residual chopping signal \citep{Agol2005,Fabrycky2012b}.  This chopping signal is best seen for the inner planet of a near-resonant system and is not seen as strongly in the outer planet, as the position of the inner planet at the time of the outer planet's transit changes very slowly when near a $j:j-1$ MMR.  The chopping signal scales with the mass of the outer planet \citep{Holman2010} and with the eccentricity \citep{Nesvorny2014}.  This additional constraint allows one to break the mass-$z_{\rm{free}}$ degeneracy to calculate the mass of each planet.

\section{Orbital and Stellar Characterization}

\subsection{Orbital fit}

Using quarters 1-16 of the \textit{Kepler} data, we extracted and flattened each transit using the IDL \texttt{AutoKep} program \citep{Gazak2012}.  For the high signal-to-noise transits of the outer planet, we used short cadence data where available.  We then used a new, modified version of the IDL program \texttt{TAP} \citep{Carter2009,Gazak2012,Eastman2013} to fit for the orbital parameters of each planet:  impact parameter ($b$), duration ($T$), the ratio of planet radius to stellar radius ($R_{\rm{p}}/R_*$), the midtransit times, linear limb darkening, quadratic limb darkening \citep{Kipping2013}, and white and red noise.  The ratio of semi-major axis to the radius of the star ($a/R_*$) and the inclination ($i$) can be derived from these parameters.  For the purpose of transit fitting, circular orbits were assumed, which we found to be a good assumption (see Section~\ref{sec:agol}).  See Figures~\ref{fig:intransit}, \ref{fig:midtransit}, and \ref{fig:outtransit} for the transit fits of PH3~b, c, and d, respectively.

The linear and quadratic limb darkening is very poorly constrained for PH3~b and c when analyzing each planet individually. Therefore, we used the posterior distribution for the linear and quadratic limb darkening from the high signal-to-noise transits (depth $\approx12500$~ppm) of PH3~d's fit as a prior for the inner and middle transit fits.  Starspot anomalies are seen in a few transits of PH3~d.  We masked out obvious spots, but we did not account for smaller scale starspot anomalies. This causes a slight residual in the light curve of PH3~d, primarily due to the ninth transit (see Figure~\ref{fig:outtransit}). For PH3~b, \texttt{TAP} was unable to simultaneously fit for both midtransit times and the orbital properties.  Therefore, we fit the phase-folded model first for the orbital parameters.  Then, holding the orbital parameters fixed, we fit for the midtransit times.

The parameter $a/R_*$ from the \texttt{TAP} fit is poorly constrained for PH3~b and PH3~c, and their best-fit values are inconsistent with PH3~d's best fit.  For example, PH3~c's $a/R_*$ is greater than PH3~d's, but it is on an interior orbit (see Table~\ref{tab:planets}).  Therefore, we revise the $a/R_*$ for each planet using Kepler's Third Law with the planet's period and the stellar parameters $M_*$ and $R_*$ (derived from stellar density; see Section~\ref{sec:stellarfit}) to get $(a/R_*)_{\rm{rev}}$.  The two values agree within errors, the unrevised $a/R_*$ approximately one $\sigma$ higher than $(a/R_*)_{\rm{rev}}$.  We use the revised $a/R_*$ to calculate the planet's incident flux ($S$) and semi-major axis ($a$).  The best fit transit light curves for PH3~b, c, and d are shown in Figures~\ref{fig:intransit}, \ref{fig:midtransit}, and \ref{fig:outtransit}, respectively.  The orbital and planetary properties are shown in Table~\ref{tab:planets}.

\begin{deluxetable*}{lccc}
\tablewidth{0pt}
\tablecaption{Orbital and planetary parameters of PH3 b, c, and d. \label{tab:planets}}
\tablehead{
\colhead{} &
\colhead{PH3~b} &
\colhead{PH3~c} &
\colhead{PH3~d} }

\startdata

$P$ (days)                          & $34.5450\pm0.0005$            & $66.0634\pm0.0114$            & $125.8518\pm0.0076$              \\
$T_0$ (JD$-$2454000)                & $965.6404\pm0.0040$           & $975.6436\pm0.0068$           & $1069.6528\pm0.0077$             \\
$M_{\rm{P}}$ ($M_\oplus$)           & $7.3\pm6.8$                   & $4.0\pm0.9$                   & $132\pm17$                       \\
$R_{\rm{P}}$ ($R_\oplus$)           & $2.15\pm0.10$                 & $2.68\pm0.17$                 & $11.59\pm0.19$                   \\
$\rho$ ($g/cm^3$)                   & $4.1\pm3.9$                   & $1.2\pm0.3$                   & $0.47\pm0.06$                    \\
$a_{\rm{rev}}$ ($AU$)               & $0.21\pm0.01$                 & $0.33\pm0.02$                 & $0.51\pm0.03$                    \\
$S$ ($S_\oplus$)                    & $24.8\pm4.4$                  & $10.7\pm1.8$                  & $4.4\pm0.8$                      \\
$t_{\rm{dur}}$ (hours)              & $3.178^{+0.055}_{-0.053}$     & $3.557^{+0.072}_{-0.065}$     & $8.067^{+0.026}_{-0.024}$        \\
$a/R_*$                             & $71.1^{+10}_{-20}$            & $117.8^{+21}_{-42}$           & $109.5\pm1.2$                    \\
$(a/R_*)_{\rm{rev}}$                & $45.9\pm0.5$                  & $70.7\pm0.7$                  & $108.6\pm1.1$                    \\
$b$                                 & $0.04^{+0.66}_{-0.68}$        & $0.05^{+0.70}_{-0.76}$        & $0.394^{+0.026}_{-0.029}$        \\
$i$ (deg)                           & $89.59^{+0.30}_{-0.48}$       & $89.73^{+0.20}_{-0.38}$       & $89.794^{+0.017}_{-0.016}$       \\
$R_{\rm{P}}/R_*$                    & $0.0197^{+0.0011}_{-0.0006}$  & $0.0246^{+0.0022}_{-0.0009}$  & $0.10620^{+0.00049}_{-0.00050}$  \\
$M_{\rm{P}}/M_*$ ($\times10^{-5}$)  & $2.0\pm1.9$                   & $1.1\pm0.2$                   & $36.43\pm4.66$                   \\
$e\cos(\omega)$                     & $-0.0215\pm0.0255$            & $-0.0035\pm0.0022$            & $0.0032\pm0.0066$                \\
$e\sin(\omega)$                     & $-0.0113\pm0.0239$            & $-0.0108\pm0.0122$            & $0.0033\pm0.0086$                \\

\enddata

\tablecomments{Best-fit parameters for the orbital and planetary properties.  The period is the mean period given over the length of observations.  The best fit $a/R_*$ from \texttt{TAP} is poorly constrained for PH3~b and PH3~c and are obviously inconsistent with PH3~d.  Therefore, we revise it using the Newton's version of Kepler's Third Law with $M_*$, $R_*$, and $P$ to get $(a/R_*)_{\rm{rev}}$.}

\end{deluxetable*}

\begin{figure*}[tbp]
	\centering
		\includegraphics[width=1.00\textwidth]{fig1.eps}
	\caption{The top panel shows the inner planet's (PH3~b's) phase-folded transit light curve (black data points) with the model overplotted in red.  The residuals are shown in the bottom panel.}
	\label{fig:intransit}
\end{figure*}

\begin{figure*}[tbp]
	\centering
		\includegraphics[width=1.00\textwidth]{fig2.eps}
	\caption{The top panel shows the middle planet's (PH3~c's) phase-folded transit light curve (black data points) with the model overplotted in red.  The residuals are shown in the bottom panel.  The extra scatter in the residuals are due to starspot anomalies.}
	\label{fig:midtransit}
\end{figure*}

\begin{figure*}[tbp]
	\centering
		\includegraphics[width=1.00\textwidth]{fig3.eps}
	\caption{The top panel shows the outer planet's (PH3~d's) phase-folded transit light curve (black data points) with the model overplotted in red.  The residuals are shown in the bottom panel.  There exists a slight residual due to starspots, primarily small-scale starspot crossings. Small spots in the ninth transit produce the majority of this effect, and due to its high density of short cadence data points, the residual is more readily noticeable.}
	\label{fig:outtransit}
\end{figure*}

\subsection{Stellar fit}
\label{sec:stellarfit}

We derived the stellar density, $\rho_*$, from the Markov chain orbital fit analysis of the outer planet, converting the observed transit characteristics into the density of the star at each link in the chain to derive the posterior distribution of the stellar density \citep{Seager2003}.  In computing the density, we calculate the sky velocity of the planet during transit; this in turn depends on the eccentricity, $e$, and longitude of periastron, $\omega$, or equivalently the eccentricity vector, $e = (e\cos{\omega},e\sin{\omega})$. Each component of the eccentricity vector of the planet was constrained by the standard deviation of its posterior distribution derived from the transit-timing analysis, which was found to have a negligible effect on stellar density (see Section \ref{sec:agol}).  This resulted in a stellar density of $\rho_* = 1.577 \pm 0.066$ g cm$^{-3}$.

We obtained a spectrum of KIC~7303287 with the HIRES instrument on the Keck I telescope \citep{Vogt1994} and performed Spectroscopy Made Easy (SME) on the data \citep{Valenti1996,Valenti2005}.  To further constrain the stellar mass ($M_*$), radius ($R_*$), and age ($t_*$), we fit the Padova PARSEC (v1.1) isochrone models \citep{Bressan2012} to the observed properties of the PH3 host star: 1) the effective temperature ($T_{\rm{eff}}$) and metallicity ($\rm{[Fe/H]}$) derived from the SME analysis; 2) the stellar density, $\rho_*$ (rather than $\log{g}$, as $\log{g}$ is poorly constrained and weakly dependent on the mass of the star on the main sequence);  3) the spectral energy distribution derived from measured stellar magnitudes in three available databases:  SDSS g, r, i, and z \citep{Abazajian2009}, 2MASS J, H, and Ks \cite{Cutri2003}, and WISE W1 and W2 \citep{Cutri2012}. For each band we used the reported uncertainty, $\sigma_i$, but also accounted for systematic uncertainty in the calibration of the flux by adding in quadrature a magnitude uncertainty, $\sigma_0$, that we allowed to float, and multiplied the likelihood by $\Pi_i(\sigma_i^2+\sigma_0^2)^{-1/2}$ \citep{Ford2006}.  We assumed a uniform prior on the extinction, $A_V$, and chose an $R_V=3.1$ Milky Way extinction law.  We also allowed the distance, $D$, and $t_*$ to vary.

There are six stellar model parameters: ($M_*$, $\rm{[Fe/H]}$, $t_*$, $D$, $A_V$, $\sigma_0$). We computed a Markov chain using the affine-invariant population Markov chain method \citep{Goodman2010,Foreman-Mackey2013} with eighteen chains (three times the number of free parameters). At each point in the chain, the likelihood was computed from the Gaussian probability of the agreement with the nine observed magnitudes, the stellar density, the effective temperature, and the metallicity.  To compute these properties, the stellar isochrones were linearly (or log-linearly) interpolated from a grid of stellar models. The stellar isochrones assumed scaled solar abundances ($\rm{[}\alpha\rm{/Fe]}=0$) and sampled metallicity in intervals of 0.1 dex, while we sampled age in intervals of 0.05 dex and mass in intervals that vary with age and metallicity.  Table~\ref{tab:star} gives the results of this analysis, which show that this is a young, solar-type star.  These stellar parameters are consistent with those derived by \citet{Rowe2014} and are moderately consistent with \citet{Huber2014}, which are also shown for comparison.   The small uncertainties on the stellar parameters are due to the precise determination of the density from the light curve and transit timing analysis and from the precise temperature and metallicity measurement from SME.  However,  these uncertainties do not account for possible systematic errors in the analysis.  Hence, we re-ran our isochrone analysis using the Dartmouth isochrones \citep{Dotter2008} and found results that were consistent within $1\sigma$ in the mean and yielded standard deviations of similar magnitude.  Also, see \citet{Torres2012} for a discussion of the systematic error and biases of SME and a comparison to stellar parameter classification (SPC) technique and MOOG.

PH3 shows strong activity clearly attributable to starspots (see Figure~\ref{fig:star}).  The long, deep duration of PH3~d shows several strong starspot anomalies, with the starspots primarily appearing at phases of about 1.2 hours before and 3.5 hours after the transit midpoint (see Figure~\ref{fig:outtransit}).  This may allow for a future investigation of its spin-orbit obliquity.  Starspots may also cause small, minute-scale TTVs themselves \citep{Mazeh2013}.  However, the orbital period of the outer planet is long compared to the rotational period, making this investigation difficult and outside the scope of our study.  The starspots allow for a determination of the rotational period, $t_{*,\rm{rot}}=8.648\pm0.0009$ days \citep{McQuillan2013}.   Gyrochronology allows one to calculate the age of the star ($t_{*,\rm{gyro}}$) using corrected $B-V$ colors and the stellar rotation period \citep{Barnes2007}.  Using the rotational period, the $B$ and $V$ magnitudes from the MAST\footnote{http://archive.stsci.edu/index.html} catalog ($B=14.816$, $V=14.255$), the earlier derived $A_V=0.06$, and the Milky Way extinction law ($R_V=3.1$), we find $t_{*,\rm{gyro}}=1.0\pm0.3$ Gyr, which is somewhat larger but consistent with the age derived from the Padova isochrones.

\begin{figure*}[tbp]
	\centering
		\includegraphics[width=1.00\textwidth]{fig4.eps}
	\caption{The normalized light curve of PH3 during quarter 7 as a function of the reduced barycentric Julian day (BJD).  The strong variability is caused by starspots and can be used to calculate its rotation period.}
	\label{fig:star}
\end{figure*}

\begin{deluxetable*}{lccc}
\tablewidth{0pt}
\tablecaption{Stellar parameters of PH3 (KIC 7303287). \label{tab:star}}
\tablehead{
\colhead{KIC 7303287}        &
\colhead{This paper}         &
\colhead{\citet{Huber2014}}  &
\colhead{\citet{Rowe2014}}   \\ }

\startdata

$M_*$ ($M_\odot$)        & $1.08\pm0.02$    & $1.16^{+0.31}_{-0.17}$   & $1.04^{\dagger}$  \\
$R_*$ ($R_\odot$)        & $1.00\pm0.02$    & $1.60^{+0.83}_{-0.48}$   & $1.16\pm0.22$     \\
$T_{\rm{eff}}$ (K)       & $5990\pm38$      & $6279^{+171}_{-215}$     & $5930\pm107$      \\
$\log{g}$ (cgs)          & $4.47\pm0.01$    & $4.10^{+0.25}_{-0.27}$   & $4.33\pm0.15$     \\
$\rho_*$ ($g/cm^3$)      & $1.58\pm0.07$    & $0.40^{+0.61}_{-0.26}$   & $0.94\pm0.42$     \\
$\rm{[Fe/H]}$            & $0.05\pm0.04$    & $-0.08^{+0.24}_{-0.30}$  & $-0.06\pm0.10$    \\
$L_*$ ($L_\odot$)        & $1.15\pm0.06$    & $3.58^{\dagger}$         & $1.50^{\dagger}$  \\
$t_*$ (Gyr)              & $0.65\pm0.44$    & \nodata                  & \nodata           \\
$t_{*,\rm{gyro}}$ (Gyr)  & $1.0\pm0.3$      & \nodata                  & \nodata           \\
$D$ (pc)                 & $700\pm14$       & \nodata                  & \nodata           \\
$A_V$                    & $0.06\pm0.01$    & \nodata                  & \nodata           \\

\enddata

\tablecomments{For comparison, we also show the stellar parameters from \citet{Huber2014} and \citet{Rowe2014}.  We note that \citet{Rowe2014} used a different Keck HIRES spectrum than us as their observational input for stellar parameters. \\
$^{\dagger}$ Parameter was not explicitly given, so calculated from other parameters as appropriate for the sake of comparison.  Due to the unknown posterior distributions of these values, we do not propagate its error. }

\end{deluxetable*}

\section{Transit Timing Mass Determinations}
\label{sec:agol}

Interactions between transiting planets can result in TTVs, which can provide constraints on the orbital parameters and masses of the transiting planets \citep{Agol2005,Holman2005}.  We modeled the times of transit with TTVFast, a newly written code for computing transit times of multi-planet transiting systems \citep{Deck2014}.  A symplectic integrator was used to compute the positions of the planets as a function of time \citep{Wisdom1991,Wisdom2006}, while Keplerian interpolation was used to find the times of transit \citep{Nesvorny2013}.   The three planets were assumed to be on plane-parallel orbits,  and the likelihood was computed using a $\chi^2$ fit to the times of transit, with the error bars of each planet inflated by a constant factor, $f_i$, which was added to the model. For each planet, five parameters were used to describe the orbit: the time of the first transit, $t_i$, the initial orbital period, $P_i$, the eccentricity vector, ($e_i\cos{\omega_i}$,$e_i\sin{\omega_i}$), and the mass ratio of the planet to the star, $\mu_i$.  With $i=1, 2, 3$, there are a total of 18 model parameters.  A prior on the likelihood of $f_i^{-N_i/2}$, with $N_i$ being the number of transits of planet $i$, was added to penalize large values of $f_i$;  this resulted in a median $\chi^2_{\rm{red}}$ for each planet of order unity \citep{Ford2006}.  A prior of $1/e_i$ was added for each planet to prevent small eccentricities from being disfavored \citep{Ford2006}.  One transit for each PH3~b and PH3~c were excluded as outliers with their contribution to $\chi^2\gg1.0$; their midtransit times were approximately $2456209$ and $2455900$ JD, respectively.

An affine-invariant ensemble  MCMC code was used to compute the posterior distribution on the parameters \citep{Goodman2010,Foreman-Mackey2013}.  An ensemble of 31 chains was run for $10^6$ generations, reaching a high degree of convergence; the first 60,000 generations were discarded when computing the posterior parameters.  Figure~\ref{fig:ttvs} shows the resulting transit-timing variations of the three planets with the 1$\sigma$ confidence intervals computed from the chain;  the plotted errors are inflated by $f=(1.73, 1.79, 1.35)$ for the inner to outer planets, respectively.  There is evidence of starspot crossings in the transit light curves, which we attempted to address by masking these regions in the light curve modeling.  However, this approach cannot take into account the effect of smaller spots crossed by the planets during a transit that lead to variation in the depth of transit that is commensurate with the photometric uncertainties.  Such spot crossings do not affect the out-of-transit data, and hence are not accounted for in the red-noise model used by TAP.  Spot crossings cause the transit shape to be asymmetric, and hence can skew the fit for the mid-transit time if not accounted for.  Hence, we include $f_i$ to account for additional systematic uncertainties on the errors in the times of transit that are not accounted for in our light curve fitting. Figure~\ref{fig:lava} shows a ``lava plot'' (or ``river plot''), demonstrating the wavy nature of TTVs.  The results show that the phase of the 1325 day TTV super period is not \textit{exactly} anti-correlated, requiring a small free eccentricity in the system, but this has a neglible contribution to $\rho_*$.  Table~\ref{tab:planets} gives the best-fit parameters for the planets' orbital parameters and masses.  All three planet masses are well within the planetary regime, confirming them as true planets.

The outer two planets in KOI-1353 are dynamically interacting, and this effect is large since they are close to a $2:1$ period ratio.  Their TTVs appear as anti-correlated sinusoidal variations with a period of $P_{TTV} = 1/(2/P_3-1/P_2) = 1325\pm 5$ days, with $P_2=66.0634\pm 0.0114$ days, and $P_3=125.8518\pm 0.0076$ days (see Equation~\ref{eq:ttv} and Figure~\ref{fig:ttvs}).  These two planets are near commensurate, but not in resonance.  The observed chopping signal breaks the mass-eccentricity degeneracy for this pair of planets.  Since the middle planet has the largest amplitude TTV, the mass of the outer planet is constrained more precisely than the mass of the middle planet.

\begin{figure*}[tbp]
	\centering
		\includegraphics[width=0.90\textwidth]{fig5.eps}
	\caption{The TTVs of each planet (top: PH3~b, middle: PH3~c, bottom: PH3~d) and the residuals from the fit.  Notice the sinusoidal chopping signal PH3~c.  This breaks the mass-eccentricity degeneracy allowing for a unique determination of mass. The blue shaded region is the one sigma best-fit region.}
	\label{fig:ttvs}
\end{figure*}

\begin{figure*}[tbp]
	\centering
		\includegraphics[width=0.90\textwidth]{fig6.eps}
	\caption{The TTVs of each planet (left: PH3~b, middle: PH3~c, right: PH3~d) in ``lava plot'' (or ``river plot'') format. The colored pixels represent the relative flux, increasing from black to red to white.  Each row of the plot shows a segment of the light curve around each transit, starting from the first transit on the bottom to the top transit at the top.  The ingress and egress of each transit is marked with blue lines, and data gaps are represented by grey bars.  Transit timing variations manifest themselves as a dark path winding back and forth in the lava and are readily seen in PH3~c.}
	\label{fig:lava}
\end{figure*}

\section{Long-term Stability}

We tested the long-term stability of the system by numerically integrating a random set of 1000 initial conditions which produced TTVs consistent with the \textit{Kepler} data. These 1000 initial conditions and corresponding planetary masses were drawn from the posterior distribution resulting from fitting the \textit{Kepler} data using Markov Chain Monte Carlo and hence form a statistically representative set of solutions.  We integrated these solutions for $10^8$ orbits of the inner planet, or $\sim 10^7$ years, using a Wisdom-Holman mapping \citep{Wisdom1991} in combination with a third-order symplectic corrector \citep{Wisdom2006}. We employed a time step of one day which resulted in a maximum fractional energy error of $\sim \times10^{-10}$.

We then determined the maximum full oscillation amplitude for the semi-major axis and eccentricity of each of the planets. Larger variations in semi-major axis compared to the average value are a sign of instability, whereas large eccentricity oscillations may result simply from large free eccentricities and, in this case, do not indicate unstable behavior. During the $10^8$ orbits of the inner planet, we only found deviations of less than a percent in the semi-major axes. The median eccentricity oscillations of the inner and outermost planet during these integrations were on the order of $30\%$ of the mean value, but the middle planet had large eccentricity oscillations of full amplitude $\sim150\%$ the mean value. Out of all the initial conditions, the maximum eccentricity oscillations were on the order of $200\%$ of the mean value.  These results are summarized in Table \ref{tab:kat}.

\begin{deluxetable*}{lccc}
\tablewidth{0pt}
\tablecaption{Results for semi-major axis and eccentricity oscillations about the mean ($\overline{a}$ and $\overline{e}$, respectively) after $10^8$ orbits of PH3~b \label{tab:kat}}
\tablehead{
\colhead{Planet}                          &
\colhead{Max [$(\Delta a)/\bar{a}] (\%)$}     &
\colhead{Median [$(\Delta e)/\bar{e}] (\%)$}  &
\colhead{Max [$(\Delta e)/\bar{e}] (\%)$}     \\ }

\startdata

PH3~b  & 0.0443259  & 34.0028  & 219.367  \\
PH3~c  & 0.302436   & 143.627  & 212.864  \\
PH3~d  & 0.0144348  & 27.6211  & 202.432  \\

\enddata

\end{deluxetable*}

The integrations indicated that the orbits were quiescent. We followed a subset of 100 of them for 10 times longer, for $10^9$ orbits of the inner planet or $10^8$ years. We find no significant deviations from the oscillation amplitudes reported in Table \ref{tab:kat}.

Although there is no analytic criterion for the stability of three planet systems, one can gain some insight by studying each pair of planets individually. Both the inner pair and outer pair satisfy the two-planet Hill criterion for stability \citep[e.g.,][]{Gladman1993} as well as the heuristic three-planet stability criterion used by \citep{Fabrycky2012a} for other \textit{Kepler} systems. Furthermore, neither pair is near a low-order mean motion resonance.  Although integrations of a length shorter than the age of the system cannot prove the system's stability, our tests strongly suggest that these orbits are long-lived.

\section{Laplace resonance}

With such low eccentricities, masses, and period ratios of 1.91, neither pair of planets is close enough to the 2:1 mean motion resonance for it to substantially affect the dynamics of the system. It is interesting to also establish how close the system is to a Laplace (or three-body) resonance. Note, however, that typically, three-body resonances are only important if one or both of the pairs is itself very close to a two-body resonance.

The equation of motion for the three-body resonance angle of the form $\phi = p \lambda_1 -(p+q) \lambda_2 +q \lambda_3$, where $p$ and $q$ are integers and $\lambda_i$ is the mean longitude of the planet $i$, resembles at lowest order the angle of a simple pendulum \citep{Aksnes1988,Quillen2011}. The Hamiltonian for this is written as:

\begin{equation}
H = \frac{1}{2} A_{p,q} p_\phi^2 +\epsilon_{p,q} \cos{\phi}
\end{equation}

\noindent where the coefficients $A_{p,q}$ and $\epsilon_{p,q}$ depend on which $(p,q)$ resonance is being considered. Therefore, one can approximate the full width of the resonance region (where the angle $\phi$ oscillates) simply as $\Delta p_\phi \sim 4 \sqrt{\epsilon_{p,q}/A_{p,q}}$. This can then be converted into a width in terms of the semi-major axis of one of the planets.

We follow the work of \citet{Quillen2011}, which focuses on the three-body resonances of circular orbits not near any two-body resonances and which shows that for equal mass planets the width of the resonance in terms of the semi-major axis of the inner planet $a_1$ roughly scales as:

\begin{equation}
\label{eq:a1}
\frac{\Delta a_1}{a_1} \sim \frac{m}{m_\star} \frac{|\log{\delta}|}{\delta} \exp{[-(p+q)\frac{\delta}{2}]}
\end{equation}

\noindent with $m$ the mass of one of the planets and for semi-major axis ratios of order $\alpha \sim 1-\delta$. In other words, for close pairs of planets the width of the resonance depends linearly on the mass of the planet, relative to the mass of the star, and is significantly larger for closer pairs of planets (smaller $\delta$). Moreover, unless $\delta$ is very small, the widths of the resonances with large values of $p$, $q$ will be exponentially small. This already suggests that the Laplace resonance will not be important for PH3.

We follow the exact formulas given in \citet{Quillen2011}, which do not assume equal mass planets, for $\epsilon_{p,q}$ and $A_{p,q}$ to determine how far away the system PH3 is from each $(p,q)$ three-body resonance. Given the orbital periods of the outer two planets and the masses of all of them, we determine the semi-major axis $a_{1,r}$ that the innermost planet would need such that the system would be in exact resonance (satisfying $\dot{\phi} = 0$), and we also determine the width of the resonance in terms of the semi-major axis of the innermost planet $\Delta a_1$.  We then computed the number $(a_1 -a_{1,r})/\Delta a_1$, the dimensionless number of how many resonance widths in semi-major axis of the inner planet the system was away from exact resonance, for each $(p,q)$ pair, with $0<p$, $q<30$.

As expected from Equation~\ref{eq:a1}, the resonances with larger values of $(p,q)$ are entirely negligible. Their widths are too small for the system to be near any of them. The PH3 system is closest to the three body resonances of the form $(p,q) = i(1, 2)$, with $i$ being an integer greater or equal to unity. However, the system is already $\sim 60$ widths away even from the $(1, 2)$ resonance and is orders of magnitude further from the majority of the resonances looked at.  From this, we conclude that three-body resonances in this regime (not near a two-body resonance and with circular orbits) are unimportant for this system.

\section{Planet Composition}

To constrain the masses of the interior planets' gas envelopes, we consider a scenario wherein they formed by core-nucleated accretion inside the snow-line and consist of an Earth-like composition rocky core (32\% by mass iron and 68\% by mass silicate) surrounded by an H/He envelope. For a given point in planet mass-radius-incident flux parameter space, planet interior structure models are used to constrain the distribution of each planet's mass between the H/He envelope and the heavy element interior. We model each planet's interior structure following an approach similar to \citet{Rogers2011}, but with updated opacities from \citet{Valencia2013}. An envelope metallicity of $30\times$ solar and a Bond albedo of 0.2 are assumed. We sample the posterior distribution on the envelope mass fraction by 1) drawing $10^5$ samples from the planet mass-radius-flux posterior distribution returned by the orbital, stellar, and TTV MCMC fits, and 2) computing the envelope mass fraction for each sample from the planet interior structure models.

\subsection{PH3~b}

With the large relative uncertainty in its measured mass, the bulk compositions that are consistent with the measured properties of PH3~b encompass a wide range of possibilities, including rocky planet, water-planet, and H/He-enshrouded rocky core scenarios. PH3~b could be a rocky planet with its transit radius defined by a rocky surface; integrating the posterior probability distribution on PH3b's mass, radius, and incident flux, we find a 16\% posterior probability that PH~3b is more dense than a pure silicate sphere. Most of the posterior probability on PH3~b's properties, however, falls in a low-density regime in which the planet must have a volatile envelope comprised of some a combination of astrophysical ices (dominated by water), hydrogen, and helium.  At the low-mass extreme, there is a 12\% posterior probability that PH3~b requires a hydrogen-dominated envelope and is less dense than a pure H$_2$O sphere, but scenarios where the planet consists of a mixture of rock and water are viable at intermediate masses (and notably at masses near the best fit). In the scenario where PH3~b has a rocky Earth-like composition heavy element interior surrounded by H/He, the H/He mass fraction of PH3~b is determined to be $0.24^{+0.28}_{-0.22}\%$ (see Figure~\ref{fig:bc}) or $\Delta R_{\rm{env}}/R_{\rm{p}} = 22\pm14\%$ ($\Delta R_{\rm{env}} = 0.47^{+0.32}_{-0.30}~R_{\oplus}$).

\subsection{PH3~c}

We find that PH3~c is $2.1^{+0.8}_{-0.3}\%$ by mass H/He (see Figure~\ref{fig:bc}), assuming the planet has a rocky core with no ices.  This corresponds to a radial extent of the H/He envelope of $\Delta R_{\rm{env}}/R_{\rm{p}} = 49^{+5}_{-4}\%$, or $\Delta R_{\rm{env}} = 1.31^{+0.32}_{-0.15}~R_{\oplus}$.  PH3~c must have a hydrogen-dominated envelope of light gases; there is less than $\sim~0.001\%$ posterior probability that PH3~c is equal density to or more dense than a pure H$_2$O body.

\begin{figure*}[tbp]
	\centering
		\includegraphics[width=1.00\textwidth]{fig7.eps}
	\caption{The posterior probability density distribution for PH3~b (top) and c (bottom), with darker blue representing a relatively higher probability, as a function of planet mass, $M_{\rm{P}}$ and the ratio of envelope mass to planetary mass, $M_{\rm{env}}/M_{\rm{P}}$. For both planets, we assume the cores are composed of heavy elements with no ices. PH3~c requires a significant envelope contribution to the total mass, while PH3~c's mass uncertainty leads to a broad range of possible compositions.}
	\label{fig:bc}
\end{figure*}

\citet{Lopez2013} also explore the mass-radius relationship as a function of incident stellar flux, age, and composition.  They find that, at fixed age and flux, the radius acts as a proxy for the composition of Neptune-like planets, such as PH3~c ($R=2.68\pm0.17R_\oplus$).  PH3~c conveniently falls almost directly on top of one of their provided grid points, (mass, incident flux, age) $=$ ($3.6M_\oplus$, $10S_\oplus$, 1 Gyr), whereas our derived values are ($4.0\pm0.9M_\oplus$, $10.7\pm1.8S_\oplus$, $0.65\pm0.44$ Gyr (model fit) or $1.0\pm0.3$ Gyr (gyrochronology)).  \citet{Lopez2013} calculates the radius for two metallicities:  solar and $50\times$ solar.  PH3 has a metallicity of $\rm{[Fe/H]}=0.05\pm0.04$, although the atmospheres of Neptunes may be significantly enhanced in metal \citep{Fortney2013,Morley2013}.  A quadratic interpolation of their grid at ($3.6M_\oplus$, $10S_\oplus$, 1 Gyr) finds that the H/He mass fraction is 3.2\% for solar metallicity and 2.5\% for the enhanced metallicity model, which is consistent with the models we have computed.

The results that PH3~c has a couple percent by mass H/He (and that PH3~b has no more than a few tenths of a percent by mass H/He) are robust. However, errors in the assumed planet energy budget, heavy element interior composition, envelope metallicity, and albedo could lead to small systematic shifts in the quantitative values of the H/He mass fractions quoted.

\subsection{PH3~d}

PH3~d lies in the Jovian planet regime; its composition is dominated by hydrogen and helium. Interpolating among model grids from 1~Gyr-old Jovian planets from \citet{Fortney2007}, we compute the planet core mass for $10^5$ draws from the mass-radius-flux posterior distribution of PH3~d. In this way, we estimate a heavy-element core mass of $14\pm4~M_{\oplus}$ (see Figure~\ref{fig:d}), corresponding to an envelope mass fraction of $89^{+0.03}_{-0.02}\%$. We note that \citet{Fortney2007} consider fully differentiated planet with cores that are 50\% ice and 50\% rock by mass. If the heavy elements are distributed through the planet envelope or are made up of a different mixture of ice and rock, the inferred planet core mass will be affected.

\begin{figure*}[tbp]
	\centering
		\includegraphics[width=1.00\textwidth]{fig8.eps}
	\caption{The posterior probability density distribution of PH3~d's core mass as a function of its total mass.}
	\label{fig:d}
\end{figure*}

\section{Discussion}

PH3~c avoided proper detection by other transit search routines, very likely due to the bias against detecting planets with large amplitude TTVs \citep{Garcia-Melendo2011}.  This may be the reason why the \textit{Kepler} pipeline misidentified this planet as a TCE with a period five times too large.  The Quasiperiodic Automated Search (QATS) algorithm \citep{Carter2013}, designed for detecting planets with TTVs, originally missed PH3~c as well due to PH3's strong stellar variability.  However, an upgraded version of QATS can now successfully identify this planet with improved detrending (E. Kruse 2014, private communication).

It is curious that both planet pairs (outer/middle and middle/inner) have very near the same period ratio ($\sim1.91$), with $|1-(P_{\rm{out}}/P_{\rm{mid}})/(P_{\rm{mid}}/P_{\rm{in}})| = 0.385\%$.  There are 127 stars with $3+$ confirmed planets, in which there are 186 unique sets of three consecutive planets\footnote{http://exoplanetarchive.ipac.caltech.edu, last accessed June 19, 2014} \citep{Wright2011}.  Only four have the same ratio of period ratios to within 0.385\%.  The Laplace resonant GJ 876 is one of these, as is the Laplace resonant candidate system HR 8799, although it is a directly imaged planetary system with large error bars in period.  The other two are Kepler-207 and Kepler-229.  In the 174 KOIs with $3+$ confirmed candidates, there are 242 unique sets of three consecutive planets. Six of these have the same period ratio to within 0.385\%:  KOIs 665, 757, 869, 1151, 1358, and 2693.

For the half of the distribution of the relative difference in the pair-wise period ratios that is not mathematically truncated at 1.0 (see the red histogram Figure~\ref{fig:perrat}), a log-normal fit can be reasonably applied to the distribution.  A variable with a log-normal distribution indicates that the variable is the product of many independent, random variables.  However, in the blue population, which includes PH3, a log-normal does not fit the allowed region of parameter space ($<1$).  There is an excess of planet triplets near zero.  This implies that the two populations are affected by different processes.  However, any further analysis of this is beyond the scope of this paper and is left for future studies.

\begin{figure*}[tbp]
	\centering
		\includegraphics[bb=160 0 410 700,scale=0.9]{fig9.pdf}
	\caption{The difference of one minus the ratio of period ratios for consecutive sets of three planets for the combined set of both \textit{Kepler} candidates and confirmed planets.  The blue histogram represents the three-planet systems where $P_{\rm{out}}/P_{\rm{mid}} < P_{\rm{mid}}/P_{\rm{in}}$, and the red histogram represents the planets where $P_{\rm{out}}/P_{\rm{mid}} > P_{\rm{mid}}/P_{\rm{in}}$, where the value of the red population is multiplied by $-1$ in order to compare to the blue population.  The purple represents overlapping red and blue populations, while the black histogram represents PH3, which is also included in the blue population.  The blue population drops off at high values of the abscissa due to the mathematical truncation at 1.0 as $(P_{\rm{out}}/P_{\rm{mid}})/(P_{\rm{mid}}/P_{\rm{in}})$ approaches 0. The two lines represent best fits for a log-normal distribution. The log-normal fit fails for the blue population as there exists an excess of planet triplets near zero.}
	\label{fig:perrat}
\end{figure*}

As of now, there is only one confirmed exoplanet system in a Laplace resonance ($\sim$ 4:2:1 MMR):  GJ 876 \citep{Rivera2010,Marti2013}.  GJ 876 b has a period of 61.12 days \citep{Marcy1998,Delfosse1998}, GJ 876 c has a period of 30.09 days \citep{Marcy2001}, and GJ 876 e orbits every 124.26 days \citep{Rivera2010}, while the period ratios of outer/middle and middle/inner of these three planets are the same to three digits, 2.03.  The Laplace resonance of GJ 876 was further explored by \citet{Marti2013}, who concluded that GJ 876's Laplace resonant solutions are stable, but surrounded by extremely unstable regions.  HD 8799 has three confirmed planets \citep{Marois2008} and also has regions of parameter space that are stable to the 4:2:1 MMR \citep{Reidemeister2009}.

In Figure~\ref{fig:mr}, we show the mass-radius diagram of PH3's planets compared to other confirmed planets (those with masses, radii, and errors from \url{http://exoplanets.org/} \citep{Wright2011}), two recent notable results discussed below, and mass-radius-composition models provided by \citet{Zeng2013}.

\begin{figure*}[tbp]
	\centering
		\includegraphics[bb=160 00 410 700,scale=0.9]{fig10.pdf}
	\caption{The mass-radius diagram for all planets with measured masses, radii, and error bars from \url{http://exoplanets.org/} \citep{Wright2011} in black circles.  Blue circles represent PH3 candidates, green squares Kepler-51 \citep{Masuda2014}, and red diamonds Kepler-79 \citep{Jontof-Hutter2014}.  The colored lines represent the masses and radii of planets in different two-component models of planetary compositions from \citet{Zeng2013} as indicated on the right.  The 100\% water-ice model does not include any steam atmosphere that would occur for warm/hot water-ice planets.}
	\label{fig:mr}
\end{figure*}

With a density of 1.2 g/cm$^3$, PH3~c joins a growing population of low-mass ($M<10M_\oplus$), low-density planets that require significant H/He gaseous atmospheres, e.g., Kepler-11 \citep[$\rho=1.7$, 1.28, 0.66, 0.58, and 0.69 g/cm$^3$ for Kepler-11~b, c, d, e, and f, respectively,][]{Lissauer2013}, Kepler-87c \citep[$\rho=0.152$ g/cm$^3$,][]{Ofir2014}, GJ 1214~b \citep[$\rho=1.87$ g/cm$^3$,][]{Charbonneau2009}, and Kepler-36c \citep[$\rho=0.89$ g/cm$^3$,][]{Carter2012}.  In the Kepler-79 system \citep{Xie2014}, \citet{Jontof-Hutter2014} calculate more precise masses of its four planets via TTVs and finds all four with low densities.  The least dense, Kepler-79d, has $\rho=0.09$ g/cm$^3$ with a mass of just $6.0M_\oplus$. Kepler-51 is even more extreme with two confirmed planets \citep{Steffen2012} both with densities of 0.03 g/cm$^3$ and masses of $2.1M_\oplus$ and $4.0M_\oplus$ \citep{Masuda2014}.  The third planet in the system confirmed by \citet{Masuda2014} has $\rho=0.05$ g/cm$^3$ and $M=7.6M_\oplus$.

The dominant method for discovering these low-mass, low-density planets is via TTV analysis.  However, TTV discoveries may be biased to low density results since, at constant mass, planets with larger radius and thus lower density will be detected more readily and with more precise midtransit times \citep{Jontof-Hutter2014}. The likelihood of transit detection approximately increases with $R^2_{\rm{P}}$.  Therefore, although more massive planets would create larger, more detectable TTV signals, the balance is likely in favor of a bias towards discovering low density planets.

As noted in \citet{Masuda2014}, a number of low-density planets are found near MMR.  However, TTV detections are best-suited for planets near MMR as they provide the clearest and largest amplitude signals, so this feature may be a selection effect.  More discoveries of low-mass, low-density planets and future theoretical work examining this new population of low-mass, low density planets will be needed to address whether this is a truly a selection effect or a signature of planet formation or migration.

\noindent\textit{Acknowledgements}

The data presented in this paper are the result of the efforts of the PH volunteers, without whom this work would not have been possible. Their contributions are individually acknowledged at \url{http://www.planethunters.org/authors}. The authors thank the  PH volunteers who participated in identifying and analyzing the candidates presented in this paper.

E.A. acknowledges funding by NSF Career grant AST 0645416; NASA Astrobiology Institutes Virtual Planetary Laboratory, supported by NASA under cooperative agreement NNH05ZDA001C; and NASA Origins of Solar Systems grant 12-OSS12-0011. K.M.D. acknowledges the support of an NSF Graduate Research Fellowship.  DF acknowledges funding support for PlanetHunters.org from Yale University and support from the NASA Supplemental Outreach Award, 10-OUTRCH.210-0001 and the NASA ADAP12-0172. LAR gratefully acknowledges support provided by NASA through Hubble Fellowship grant \#HF-51313 awarded by the Space Telescope Science Institute, which is operated by the Association of Universities for Research in Astronomy, Inc., for NASA, under contract NAS 5-26555.  KS gratefully acknowledges support from Swiss National Science Foundation Grant PP00P2\_138979/1. MES is supported in part by an Academia Sinica postdoctoral fellowship.

This research has made use of the Exoplanet Orbit Database and the Exoplanet Data Explorer at exoplanets.org. The Zooniverse is supported by The Leverhulme Trust and by the Alfred P. Sloan foundation. PH is supported in part by NASA JPL's PlanetQuest program. The Talk system used by PH was built during work supported by the NSF under Grant No. DRL-0941610. We gratefully acknowledge the dedication and achievements of \textit{Kepler} science team and all those who contributed to the success of the mission. We acknowledge use of public release data served by the NASA/IPAC/NExScI Star and Exoplanet Database, which is operated by the Jet Propulsion Laboratory, California Institute of Technology, under contract with the National Aeronautics and Space Administration. This research has made use of NASA's Astrophysics Data System Bibliographic Services. This paper includes data collected by the \textit{Kepler} spacecraft, and we gratefully acknowledge the entire \textit{Kepler} mission team's efforts in obtaining and providing the light curves used in this analysis. Funding for the \textit{Kepler} mission is provided by the NASA Science Mission directorate. Support for MAST for non-HST data is provided by the NASA Office of Space Science via grant NNX13AC07G and by other grants and contracts.  Some of the data presented in this paper were obtained from the Mikulski Archive for Space Telescopes (MAST). STScI is operated by the Association of Universities for Research in Astronomy, Inc., under NASA contract NAS5-26555. Support for MAST for non-HST data is provided by the NASA Office of Space Science via grant NNX13AC07G and by other grants and contracts.  The data presented herein were partly obtained at the W. M. Keck Observatory, which is op-erated as a scientific partnership among the California Institute of Technology, the University of California, and the National Aeronautics and Space Administration. The Observatory was made possible by the generous financial support of the W. M. Keck Foundation.

redpajama_set_name: RedPajamaArXiv