# On the Performance of Jerk-Constrained Time-Optimal
Trajectory Planning for Industrial Manipulators

Jee-eun Lee1, Andrew Bylard2,∗, Robert Sun2, and Luis Sentis1,∗ This work was
supported by Dexterity, Inc.* Corresponding authors contributed equally.1 J.
Lee and L. Sentis are with The University of Texas at Austin, Austin, TX, USA,
{jelee<EMAIL_ADDRESS>A. Bylard, and R. Sun are with Dexterity, Inc.,
Redwood City, CA, USA, {andrew.bylard<EMAIL_ADDRESS>
###### Abstract

Jerk-constrained trajectories offer a wide range of advantages that
collectively improve the performance of robotic systems, including increased
energy efficiency, durability, and safety. In this paper, we present a novel
approach to jerk-constrained time-optimal trajectory planning (TOTP), which
follows a specified path while satisfying up to third-order constraints to
ensure safety and smooth motion. One significant challenge in jerk-constrained
TOTP is a non-convex formulation arising from the inclusion of third-order
constraints. Approximating inequality constraints can be particularly
challenging because the resulting solutions may violate the actual
constraints. We address this problem by leveraging convexity within the
proposed formulation to form conservative inequality constraints. We then
obtain the desired trajectories by solving an $\boldsymbol{n}$-dimensional
Sequential Linear Program (SLP) iteratively until convergence. Lastly, we
evaluate in a real robot the performance of trajectories generated with and
without jerk limits in terms of peak power, torque efficiency, and tracking
capability.

###### Index Terms:

Time-Optimal Trajectory Planning, Jerk Constraints, Smooth Trajectory
Generation, Industrial Robots

##  I INTRODUCTION

In automation, time-optimal motion planning plays a crucial role in maximizing
productivity while ensuring efficiency and safety. However, solving the
complete trajectory-planning problem is often difficult to perform in real
time [1, 2, 3, 4, 5]. To obtain satisfactory results using reasonable
computational resources, trajectory planning is often solved in two stages [6,
7, 8]. As shown in Fig. 1, during the initial stage, the path generator
creates a geometric path that incorporates task specifications and obstacle
avoidance. In the subsequent stage, a time-optimal trajectory is calculated
along the specified geometric path, enforcing constraints such as permissible
joint velocity and torque.

Figure 1: Two-Stage Motion Planning and Control Architecture: The path
generator creates a collision-free path, which is fed into the time-optimal
trajectory planner to compute the desired trajectory satisfying robot hardware
limits. Subsequently, the predictive proportional-integral (PPI) controller
generates high-frequency joint commands with feedback.

This paper focuses on solving the second-stage motion planning sub-problem,
also known as time-optimal trajectory planning (TOTP) or path tracking (TOPT)
[9, 10]. Although TOTP under first- and second-order constraints (such as
velocity, acceleration, and torque limits) has been thoroughly explored,
third-order constraints such as jerk or torque rate bounds are often
overlooked. This is primarily because their inclusion lacks clear guidance
(for instance, these bounds may not be specified in motor usage requirements),
and they introduce non-convex terms into the formulation, making problem-
solving much more challenging. However, without considering jerk limits,
solution trajectories may include large acceleration jumps, resulting in
various issues when applied to real robots [11].

In this paper, we propose a new approach for TOTP with practical third-order
constraints. By incorporating jerk limits into trajectory planning and control
strategies, robotic systems can achieve higher reliability and effectiveness.
Some well-known advantages include: 1) Enhanced path tracking: a smoother
trajectory allows the robot’s movements to more closely match the desired
trajectory, minimizing jerk-induced vibrations. 2) Reduced wear and tear:
Sudden changes in acceleration can subject mechanical components to higher
stress levels. By constraining jerk, the robot’s motor gears undergo gentler
transitions, extending component lifespans. 3) Energy efficiency: Jerk-limited
motions prevent unnecessary spikes in power consumption, leading to decreased
energy usage and reduced operational costs.

Finally, we perform real robot experiments to confirm the benefits of
including jerk limits. We analyze factors such as peak power and torque, which
affect the performance, durability, and maintainability of the robot. The
proposed motion planning and control framework is depicted in Fig. 1.

###  I-A Related Works

Traditionally, there are two widely adopted approaches for addressing time-
optimal trajectory planning (TOTP): numerical integration (NI), and convex
optimization (CO). NI-based approaches [9, 10, 12] solve the problem based on
the insight that one can obtain a time-optimal trajectory through bang-bang
control. However, this method requires identifying switch points where
acceleration changes, and it is prone to failure in certain cases. Conversely,
optimization provides a general and robust means to formulate and solve the
problem [13, 14]. However, it increases the problem’s complexity by at least
an order of magnitude, resulting in longer computation times. Other works
solve the problem using dynamic programming [15, 16] or by leveraging specific
conditions which may be too constraining, such as imposing constant first- or
second-order limits in the path parameter space [17, 18], which may not map to
desired (e.g., constant) limits in the joint space.

Recently, there have been efforts to simplify TOTP algorithms through methods
based in reachability analysis (RA) [19, 20]. These methods, instead of
solving a high-dimensional optimization problem (with $n$ optimization
variables, one for each timestep), can achieve the optimal solution by solving
a series of 2$n$ 2-D sub-problems. However, it’s worth noting that these
reachability-based algorithms are restricted to second-order constraints,
since they cannot handle situations where second-order constraints violate the
third-order constraints, known as third-order singularity [21].

Many works have attempted to incorporate third-order constraints into time-
optimal trajectory planning (TOTP). Many of them rely on the NI-based
approaches [21, 22, 23, 24]. In [22], however, this formulation is limited to
bounding the third-order derivative of the path parameter ($\dddot{s}$) rather
than directly limiting the robot jerk or torque rate, and the third-order
singularity issue is not addressed. Instead, one can carefully bridge the
maximum and minimum jerk profile to avoid failure from singularity [21].
However, this approach requires handling numerous exceptional cases, making
the algorithm overly complex when adding first and second-order constraints,
which are often more critical [24].

Our approach provides a more general formulation within a convex optimization
framework. A key challenge is enforcing nonlinear third-order inequality
constraints, as approximating them without careful consideration can lead to
solutions that violate real constraints. An alternative is to minimize both
jerk and time [25], but this requires careful tuning of the relative weight
for jerk and time in the cost function. Since the solutions of the
optimization problem are often very sensitive to the weight, the solution can
frequently sacrifice motion time or smoothness excessively depending on the
input path given, even under the same weight. In this paper, we formulate the
problem as an SLP using conservative linearization of the nonlinear
constraints at each iteration to ensure solution feasibility. The work in [26]
also used conservative approximation for third-order constraints, based on the
fact that torque obtained without third-order constraints is always bigger
than the one with third-order constraints. Thus, their approach is limited to
a one-time approximation, which can lead to suboptimal solutions, with
stricter limits often widening the optimality gap.

###  I-B Contribution

  * •

We formulate a unique TOTP with third-order constraints as a Sequential Linear
Program (SLP), exceeding the performance of existing methods.

  * •

We introduce a method for approximating general third-order constraints,
enabling iterative optimization to converge to the optimal solution.

  * •

We demonstrate that imposing jerk limits significantly reduces peak power,
enhances energy efficiency, and improves tracking performance.

###  I-C Organization

The paper is organized as follows. Section II revisits the concepts of TOTP
and the corresponding parameterization methods. Section III formulates the
jerk constraints for TOTP and outlines the overall optimization process with
the approximated jerk constraints. In Section IV, we discuss the numerical and
experimental results of the proposed jerk-bounded time-optimal trajectory
planning. Finally, in Section V we discuss conclusions and future work.

##  II PRELIMINARIES

###  II-A Problem Formulation

We formulate time-optimal trajectory planning (TOTP) as follows: Given a path
$\mathbf{q}(s)$ as a function of a scalar path coordinate $s\in[0,1]$, find a
monotonically increasing time scaling $s(t):[0,T]\rightarrow[0,1]$ that

  * •

satisfies an initial state $(s_{0},\dot{s}_{0})=(0,0)$ and a final state
$(s_{\textrm{end}},\dot{s}_{\textrm{end}})=(1,0)$,

  * •

minimizes the total travel time $T$ along the path,

  * •

enforces kinematics and dynamics constraints imposed by robot hardware
limitations.

###  II-B Time-Optimal Trajectory Planning (TOTP)

For the purpose of solving the TOTP problem to track the prescribed path
$\mathbf{q}(s)$ for a robot manipulator, we express the first- and second-
order time derivatives of $\mathbf{q}(s)$ as

$\displaystyle\dot{\mathbf{q}}(s)=\mathbf{q}^{\prime}(s)\dot{s},\quad\ddot{\mathbf{q}}(s)=\mathbf{q}^{\prime\prime}(s)\dot{s}^{2}+\mathbf{q}^{\prime}(s)\ddot{s}.$
(1)

We also transform the dynamics of the manipulator into

$\displaystyle\boldsymbol{\tau}$
$\displaystyle~{}=M(\mathbf{q})\ddot{\mathbf{q}}+C(\mathbf{q},\dot{\mathbf{q}})\dot{\mathbf{q}}+g(\mathbf{q})$
$\displaystyle~{}=M(\mathbf{q})(\mathbf{q}^{\prime\prime}\dot{s}^{2}+\mathbf{q}^{\prime}\ddot{s})+C(\mathbf{q},\mathbf{q}^{\prime})\mathbf{q}^{\prime}\dot{s}^{2}+g(\mathbf{q})$
$\displaystyle:=\mathbf{m}(s)\ddot{s}+\mathbf{c}(s)\dot{s}^{2}+\mathbf{g}(s).$
(2)

Note that from the linear nature of the Coriolis matrix, we have
$C(\mathbf{q},\mathbf{q}^{\prime}\dot{s})\mathbf{q}^{\prime}\dot{s}=C(\mathbf{q},\mathbf{q}^{\prime})\mathbf{q}^{\prime}\dot{s}^{2}$
above.

Furthermore, it is known that a second-order time differential equation can be
transformed into a first-order differential equation in $(\dot{s}^{2},s)$
based on the relation introduced in [10],

$\displaystyle\ddot{s}=\frac{d\dot{s}}{dt}=\frac{d\dot{s}}{ds}\frac{ds}{dt}=\dot{s}^{\prime}\dot{s}=\frac{1}{2}(\dot{s}^{2})^{\prime}.$
(3)

This is beneficial for solving TOTP problems, as it allows writing all the
equations in linear form by setting ($\dot{s}^{2},\ddot{s}$) as optimization
variables [19, 13].

###  II-C TOTP Constraints in the Discretized System

We now show that we can express all constraints up to second-order as linear
in $x=\dot{s}^{2}$ for the discretized system. By dividing the path interval
$[0,1]$ into $N$ segments, we have

$\displaystyle 0=:s_{0},s_{1},...,s_{N-1},s_{N}:=1.$

Then from the relation in (3), $\dot{s}_{k},\ddot{s}_{k}$ can be represented
as:

$\dot{s}_{k}^{2}:=x_{k},\quad\ddot{s}_{k}=\frac{x_{k+1}-x_{k}}{2\triangle_{k}}\quad\textrm{or}\quad\frac{x_{k}-x_{k-1}}{2\triangle_{k-1}},$
(4)

where $\triangle_{k}:=s_{k+1}-s_{k}$ and $\triangle_{k-1}:=s_{k}-s_{k-1}$.

For representational convenience, we simplify the notation of each quantity at
$s_{k}$ as $(\cdot)(s_{k})=(\cdot)_{k}$. Instead of introducing another
optimization variable $u=\ddot{s}$ like other works [19], we express all
quantities as linear in $x$, which is sufficient. Finally, we formulate the
following constraints:

  * •

Velocity limits (1st-order constraints) First, the joint velocity limit can be
represented as

$\displaystyle\dot{\mathbf{q}}_{k}\circ\dot{\mathbf{q}}_{k}=(\mathbf{q}_{k}^{\prime}\circ\mathbf{q}_{k}^{\prime})x_{k}$
$\displaystyle\leq$
$\displaystyle\dot{\mathbf{q}}_{k}^{max}\circ\dot{\mathbf{q}}_{k}^{max},$ (5)

where $\mathbf{a}\circ\mathbf{b}$ represents the Hadamard product (i.e., the
element-wise product) of two vectors.

  * •

Acceleration limits (2nd-order constraints) Similarly, from the parameterized
expression of acceleration

$\ddot{\mathbf{q}}_{k}=\mathbf{q}^{\prime\prime}_{k}\dot{s_{k}}^{2}+\mathbf{q}^{\prime}_{k}\ddot{s}_{k}=\mathbf{q}^{\prime\prime}_{k}x_{k}+\mathbf{q}^{\prime}_{k}\left(\frac{x_{k+1}-x_{k}}{2\triangle_{k}}\right),$

we can constrain two consecutive $x_{k}$ by

$\displaystyle\left|\frac{\mathbf{q}^{\prime}_{k}}{2\triangle_{k}}x_{k+1}+\left(\mathbf{q}^{\prime\prime}_{k}-\frac{\mathbf{q}^{\prime}_{k}}{2\triangle_{k}}\right)x_{k}\right|\leq\ddot{\mathbf{q}}_{k}^{max}$
$\displaystyle\textrm{or}\quad\left|\left(\mathbf{q}^{\prime\prime}_{k}+\frac{\mathbf{q}^{\prime}_{k}}{2\triangle_{k-1}}\right)x_{k}-\frac{\mathbf{q}^{\prime}_{k}}{2\triangle_{k-1}}x_{k-1}\right|\leq\ddot{\mathbf{q}}_{k}^{max}.$
(6)

  * •

Torque limits (2nd-order constraints) The parameterized dynamics can also be
rewritten with respect to $x$ by substituting (4) into (2):

$\displaystyle\left|\frac{\mathbf{m}_{k}}{2\triangle_{k}}x_{k+1}+\left(\mathbf{c}_{k}-\frac{\mathbf{m}_{k}}{2\triangle_{k}}\right)x_{k}+\mathbf{g}_{k}\right|$
$\displaystyle\leq\boldsymbol{\tau}_{k}^{max}$
$\displaystyle\textrm{or}~{}\left|\left(\mathbf{c}_{k}+\frac{\mathbf{m}_{k}}{2\triangle_{k-1}}\right)x_{k}-\frac{\mathbf{m}_{k}}{2\triangle_{k-1}}x_{k-1}+\mathbf{g}_{k}\right|$
$\displaystyle\leq\boldsymbol{\tau}_{k}^{max}.$ (7)

By stacking up all the constraints derived above, we can simply express first-
order constraints as $\boldsymbol{\omega}_{k}x_{k}\leq\boldsymbol{\nu}_{k}$,
where
$\boldsymbol{\omega}_{k}=\mathbf{q}_{k}^{\prime}\circ\mathbf{q}_{k}^{\prime}$
and
$\boldsymbol{\nu}_{k}=\dot{\mathbf{q}}_{k}^{max}\circ\dot{\mathbf{q}}_{k}^{max}$.
Similarly, second-order constraints can be represented in the form of
$\boldsymbol{\alpha}_{k}^{0}x_{k}+\boldsymbol{\alpha}_{k}^{1}x_{k+1}\leq\boldsymbol{\beta}_{k}$.

##  III Jerk-Constrained Time-Optimal Trajectory Planning

###  III-A Jerk Constraint Formulation

Suppose we want to enforce joint-level jerk bounds

$-\dddot{\mathbf{q}}_{\textrm{max}}\leq\dddot{\mathbf{q}}\leq\dddot{\mathbf{q}}_{\textrm{max}}.$
(8)

In this paper, we define jerk as change in acceleration between adjacent
trajectory points. While this formula is more forgiving than an instantaneous
rate, it is still effective in preventing abrupt changes in acceleration and
simplifies the formulation. Additional details can be found in Section V. Then
joint jerk given the discretized system can be formulated as:

$\displaystyle\dddot{\mathbf{q}}_{k}=\frac{\Delta\ddot{\mathbf{q}}}{\Delta
t}(s_{k})$
$\displaystyle=\frac{\left(\ddot{\mathbf{q}}_{k+1}-\ddot{\mathbf{q}}_{k}\right)}{\frac{1}{2}(\triangle
t_{k+1}+\triangle t_{k})}$
$\displaystyle=\frac{\mathbf{j}_{k}^{2}x_{k+2}+\mathbf{j}_{k}^{1}x_{k+1}+\mathbf{j}_{k}^{0}x_{k}}{h_{k}(\mathbf{x}_{k:k+2})},$
(9) where $\displaystyle\mathbf{x}_{k:k+2}$
$\displaystyle=[x_{k}~{}x_{k+1}~{}x_{k+2}]^{\top}$ $\displaystyle
h_{k}(\mathbf{x}_{k:k+2})$
$\displaystyle={\dfrac{\triangle_{k+1}}{\sqrt{x_{k+2}}+\sqrt{x_{k+1}}}+\dfrac{\triangle_{k}}{\sqrt{x_{k+1}}+\sqrt{x_{k}}}}$
$\displaystyle\mathbf{j}_{k}^{2}=\frac{\mathbf{q}^{\prime}_{k+1}}{2\triangle_{k+1}}$
$\displaystyle,\quad\mathbf{j}_{k}^{1}=\mathbf{q}^{\prime\prime}_{k+1}-\frac{\mathbf{q}^{\prime}_{k+1}}{2\triangle_{k+1}}-\frac{\mathbf{q}^{\prime}_{k}}{2\triangle_{k}},\mathbf{j}_{k}^{0}$
$\displaystyle=-\mathbf{q}_{k}^{\prime\prime}+\frac{\mathbf{q}_{k}^{\prime}}{2\triangle_{k}}.$

Unfortunately, unlike second- or lower-order constraints, a nonlinear term
such as $h_{k}(\mathbf{x}_{k:k+2})$ in Eqn. (9) always remains when expressing
jerk in the TOTP problem. We address this issue by linearizing the nonlinear
term $h_{k}(\mathbf{x}_{k:k+2})$ and leveraging its convexity to ensure the
solution always satisfies the original constraints, i.e.,

$\displaystyle h_{k}(\mathbf{x}_{k:k+2})$ $\displaystyle\geq
h_{k}(\bar{\mathbf{x}}_{k:k+2})+\nabla
h_{k}(\bar{\mathbf{x}}_{k:k+2})({\mathbf{x}}_{k:k+2}-\bar{\mathbf{x}}_{k:k+2})$
$\displaystyle:=\bar{h}_{k}+\bar{h}_{k}^{2}x_{k+2}+\bar{h}_{k}^{1}x_{k+1}+\bar{h}_{k}^{0}x_{k},$

where

$\displaystyle\bar{h}_{k}^{i}$ $\displaystyle=\frac{\partial
h_{k}(\mathbf{x}_{k:k+2})}{\partial
x_{k+i}}\bigg{|}_{\mathbf{x}=\mathbf{\bar{x}}},\quad i\in\\{0,1,2\\},$
$\displaystyle\bar{h}_{k}$
$\displaystyle=h_{k}(\bar{\mathbf{x}}_{k:k+2})-\bar{h}_{k}^{2}\bar{x}_{k+2}-\bar{h}_{k}^{1}\bar{x}_{k+1}-\bar{h}_{k}^{0}\bar{x}_{k}.$

We show the convexity of this function in the Appendix. From this, we can
reformulate the joint jerk constraints as

$\displaystyle\big{|}~{}\mathbf{j}_{k}^{2}x_{k+2}$
$\displaystyle+\mathbf{j}_{k}^{1}x_{k+1}+\mathbf{j}_{k}^{0}x_{k}~{}\big{|}$
$\displaystyle\leq\dddot{\mathbf{q}}_{\textrm{max}}\left(\bar{h}_{k}+\bar{h}_{k}^{2}x_{k+2}+\bar{h}_{k}^{1}x_{k+1}+\bar{h}_{k}^{0}x_{k}\right)$
$\displaystyle\leq\dddot{\mathbf{q}}_{\textrm{max}}h_{k}(\mathbf{x}_{k:k+2}),$

which finally leads to the linear form of 3rd-order constraints:

$\boldsymbol{\gamma}_{k}^{0}x_{k}+\boldsymbol{\gamma}_{k}^{1}x_{k+1}+\boldsymbol{\gamma}_{k}^{2}x_{k+2}\leq\boldsymbol{\eta}_{k}.$

We can now combine all first- to third-order constraints as a linear matrix
inequality:

1st order constraints: $\boldsymbol{\omega}_{k}x_{k}\leq\boldsymbol{\nu}_{k}$
2nd order constraints:
$\boldsymbol{\alpha}_{k}^{0}x_{k}+\boldsymbol{\alpha}_{k}^{1}x_{k+1}\leq\boldsymbol{\beta}_{k}$
(10) 3rd order constraints:
$\boldsymbol{\gamma}_{k}^{0}x_{k}+\boldsymbol{\gamma}_{k}^{1}x_{k+1}+\boldsymbol{\gamma}_{k}^{2}x_{k+2}\leq\boldsymbol{\eta}_{k}$

Figure 2: Test motion snapshots. A: Front place motion to move boxes from a
conveyor to the top of a virtual stack of packages. B: Pre-pick motion
preparing to pick a box from a conveyor: C. Bottom place motion to place boxes
onto the bottom of the truck. D: Post-place motion to return to the ready pose
after placing boxes. E and F show box unloading motions, which return the box
from a stack of packages to the conveyor with and without jerk constraints.
Without jerk constraints, the robot’s suction gripper is also more likely to
lose grasp of the box.

###  III-B Trajectory Optimization

The cost function for minimizing the total time required to follow the given
path can be expressed as

$f(\mathbf{x})=\sum_{i=0}^{N-1}\frac{s_{i+1}-s_{i}}{\sqrt{x_{i}}+\sqrt{x_{i+1}}},$
(11)

which can be linearized along the nominal trajectory $\bar{\mathbf{x}}$.

Now we formulate the trajectory optimization problem as a Sequential Linear
Program (SLP), where the optimization variables are
$\mathbf{x}=[x_{1},\cdots,x_{N-1}]^{\top}$ as follows:

$\displaystyle\min_{\mathbf{x}}\quad\mathbf{c}^{\top}\mathbf{x}$ (12)
$\displaystyle\textrm{subject to}\quad A\mathbf{x}\leq b,$

where

$\mathbf{c}=\frac{\partial
f}{\partial\mathbf{x}}\bigg{|}_{\mathbf{x}=\mathbf{\bar{x}}},\;A=\begin{bmatrix}A^{1}\\\
A^{2}\\\ A^{3}\end{bmatrix},\;b=\begin{bmatrix}b^{1}\\\ b^{2}\\\
b^{3}\end{bmatrix}$ $\displaystyle A^{1}$
$\displaystyle=\begin{bmatrix}\boldsymbol{\omega}_{1}&0&\cdots&0\\\
\vdots&&&\vdots\\\
0&\cdots&0&\boldsymbol{\omega}_{N-1}\end{bmatrix},\quad\;\;b^{1}=\begin{bmatrix}\boldsymbol{\nu}_{1}\\\
\vdots\\\ \boldsymbol{\nu}_{N-1}\end{bmatrix},$ $\displaystyle A^{2}$
$\displaystyle=\begin{bmatrix}\boldsymbol{\alpha}^{1}_{0}&0&\cdots&0\\\
\boldsymbol{\alpha}^{0}_{1}&\boldsymbol{\alpha}^{1}_{1}&\cdots&0\\\
\vdots&&&\vdots\\\
0&\cdots&\boldsymbol{\alpha}^{0}_{N-2}&\boldsymbol{\alpha}^{1}_{N-2}\\\
0&\cdots&\cdots&\boldsymbol{\alpha}^{0}_{N-1}\end{bmatrix},b^{2}=\begin{bmatrix}\boldsymbol{\beta}_{0}-\boldsymbol{\alpha}_{0}^{0}x_{0}\\\
\boldsymbol{\beta}_{1}\\\ \vdots\\\ \boldsymbol{\beta}_{N-2}\\\
\boldsymbol{\beta}_{N-1}-\boldsymbol{\alpha}_{N-1}^{1}x_{N}\end{bmatrix},$
$\displaystyle A^{3}$
$\displaystyle=\begin{bmatrix}\boldsymbol{\gamma}^{2}_{0}&0&0&0&\cdots&0\\\
\boldsymbol{\gamma}^{1}_{1}&\boldsymbol{\gamma}^{2}_{1}&0&0&\cdots&0\\\
\boldsymbol{\gamma}^{0}_{2}&\boldsymbol{\gamma}^{1}_{2}&\boldsymbol{\gamma}^{2}_{2}&0&\cdots&0\\\
\vdots&&\vdots&\ddots&\ddots&\vdots\\\
0&\cdots&0&\boldsymbol{\gamma}^{0}_{N-3}&\boldsymbol{\gamma}^{1}_{N-3}&\boldsymbol{\gamma}^{2}_{N-3}\\\
0&\cdots&0&0&\boldsymbol{\gamma}^{0}_{N-2}&\boldsymbol{\gamma}^{1}_{N-2}\end{bmatrix},$
(13) $\displaystyle b^{3}$
$\displaystyle=\begin{bmatrix}\boldsymbol{\eta}_{0}-\boldsymbol{\gamma}^{0}_{0}x_{0}-\boldsymbol{\gamma}^{1}_{0}x_{1}\\\
\boldsymbol{\eta}_{1}-\boldsymbol{\gamma}^{0}_{1}x_{1}\\\
\boldsymbol{\eta}_{2}\\\ \vdots\\\ \boldsymbol{\eta}_{N-3}\\\
\boldsymbol{\eta}_{N-2}-\boldsymbol{\gamma}^{2}_{N-2}x_{N}\end{bmatrix}.$ (14)

Note that due to the linearization, we need to iteratively update the nominal
trajectory to reformulate the coefficients for cost function and 3rd-order
constraints. The steps of the resulting algorithm are shown in Algorithm 1.

Algorithm 1 TOTP with 3rd-order constraints (TOTP3)

Nominal variables $\bar{\mathbf{x}}_{1:N}$

while True do

Update $\mathbf{c},A^{3},b^{3}$ $\triangleright$ Eqn. (III-B),(14)

$\mathbf{x}_{1:N}$ = solve LP $\triangleright$ Eqn. (12)

if $\|\mathbf{x}_{1:N}-\bar{\mathbf{x}}_{1:N}\|<\epsilon$ then break

end if

Update $\bar{\mathbf{x}}_{1:N}$

end while

Figure 3: Comparison of the velocity curve for the pre-pick motion (Fig. 2B)
computed from TOPP-RA [19] without jerk limits (black), [26] with pseudo jerk
limits (blue) and the proposed approach with jerk limits (red).

##  IV EXPERIMENT RESULTS

In this section, we provide experimental results to validate the effectiveness
of the proposed approach.

Figure 4: Algorithm performance comparison on real-robot for the front place
motion (Fig. 2A): A shows tracking performance of the algorithms with and
without jerk limits by comparing the desired, commanded, and actual joint
position and velocity. B shows the corresponding estimated torque and power.
Figure 5: Time-optimal trajectory generated for the front place motion (Fig.
2A) without jerk limits (blue) and with jerk limits (red).

###  IV-A Experiment Setup

For practical evaluation, we used an example motion from a real-world box-
loading application. We performed the test motions depicted and described in
Fig. 2 with and without the jerk limit. We used the output from TOPP-RA [19]
as the initial nominal trajectory for our algorithm and as a comparison target
for the optimal trajectory without the jerk limit. Finally, the jerk limits
were heuristically chosen to be $100\sim 1000\textnormal{
rad}/\textnormal{s}^{3}$, balancing the trade-off between motion duration and
overall performance metrics such as energy efficiency and tracking (note that
some applications may require more specific strict limitations on jerk).
Lastly, we used 7-DOF robotic system composed of a 6-DOF RS020N Kawasaki arm
mounted on an additional revolute joint for our experiment.

###  IV-B Comparison With TOPP-RA

####  IV-B1 Simulation Results

We compare the TOTP outcomes obtained by three different methods in the phase
plane. These methods include TOPP-RA [19] without jerk limits, TOTP3 [26] with
pseudo jerk limits, and our proposed approach with jerk limits. While the
original work in [26] focused on torque rate constraints, we adapted it for
our experiment by introducing jerk limits. We employed a similar methodology
to the one used for pseudo-torque rate to establish the pseudo-jerk limit.
Fig. 3 shows that the curve tends to have a smoother shape when jerk limits
are enforced. Note that since we formulate jerk limits in joint space and then
map them to the space of the path parameter $s$, the smoothness of the joint-
space curves along the path may vary at different path points. The result
shows that the proposed approach yields faster motion than the pseudo-torque
rate method under the same limits. To highlight the optimality gap, we applied
a very stringent jerk limit of $100\textnormal{ rad}/\textnormal{s}^{3}$.

Fig. 5 illustrates the computed velocity, acceleration, and jerk with the
corresponding limits for the significant axes of the robot. It demonstrates
that by enforcing a jerk limit, rapid changes in acceleration are mitigated,
leading to a smoother velocity profile.

####  IV-B2 Robot Implementation Results

While we can observe jerk limits producing a smoother trajectory in
simulation, we must still establish their tangible benefits on a real robot.
We performed a comparative test on the real Kawasaki robot system described
above, running the same path with and without jerk limits. The results are
shown in Fig. 4.

Fig. 4A shows that the robot deviated from its desired path during rapid
acceleration changes, leading to overshooting and increased fluctuations in
the velocity plots. However, applying jerk limits significantly smoothed the
motions, resulting in differences in the generated torque and power, as
illustrated in Fig. 4B. The jerk-limited motions prevented aggressive velocity
spikes, lowering peak power and torque usage, with significant reductions for
some joints.

We measured peak power and RMS torque and the results are presented in Table
I, showing that peak power was reduced by about 25%, and RMS torque was
reduced to half of its original value by limiting jerk.

TABLE I: RMS torque and peak power during front place motion | Axis
---|---
| 0 | 1 | 2 | 3 | 4 | 5 | 6
RMS Trq [Nm] | TOPP-RA | 3624 | 965 | 32 | 73 | 795 | 100 | 28
TOTP3 | 1504 | 377 | 14 | 77 | 267 | 35 | 6.4
Peak Pwr [$\textnormal{Nm}^{2}$/s] | TOPP-RA | 1119 | 113 | 57 | 40 | 102 | 18 | 8.5
TOTP3 | 612 | 67 | 27 | 38 | 34 | 6.6 | 1.2

####  IV-B3 Algorithm Efficiency

We ran experiments via a C++ implementation running on a 64-bit system with a
2.35 GHz AMD EPYC 7452 32-core processor. To show the algorithm’s efficiency,
we measured the computation time of TOTP3 and TOPP-RA by averaging 20 trials
of box-place motions. Each trial involved 20-30 path points beginning from
slightly different initial box lift poses. TOTP3 (w/ jerk limit) took
$7.533\pm 5.848$ ms while TOPP-RA (w/o jerk limit) took $0.273\pm 0.089$ ms to
obtain the solution. The TOPT3 algorithm’s runtime was longer than TOPP-RA due
to the additional complexity of enforcing third-order constraints. However, it
could still be executed within 10 ms, very reasonable for real-time
applications. We were able to accelerate our computation times by warm-
starting our algorithm with the non-jerk-limited TOPP-RA result and
recomputing the only changed portion at each iteration.

Figure 6: Motion duration for several motions computed with and without jerk
limit.

Lastly, we measured the motion time generated by each algorithm for several
different motions. Notably, imposing a jerk limit of $1000\textnormal{
rad}/\textnormal{s}^{3}$ increased overall motion duration by just 3-5%. For
the “Post-Place” motion, we observed lower motion time with jerk limit, but
this is just due to a known case of suboptimality in TOPP-RA, which guarantees
optimality only if path points exclude zero-inertia point (see [19] for
details), which requires denser path segmentation. On the other hand, the
proposed algorithm’s SLP formulation allows it to converge to the optimal
solution.

##  V DISCUSSION & CONCLUSIONS

In summary, we introduced a novel approach to solve the time-optimal
trajectory planning problem while enforcing third-order constraints. In
particular, we constrained the acceleration change between path points to
produce smoother and more trackable trajectories. Using the proposed
formulation, we removed the need for high-order dynamics variables which are
computationally difficult to handle. For instance, the jerk formulation in
[19] requires $\mathbf{q}^{\prime\prime\prime}$, and the torque rate
formulation in [26] requires $\mathbf{M}^{\prime}$ and $\mathbf{C}^{\prime}$.
These quantities can be challenging to compute from standard robot dynamics
libraries. Though this paper tested only jerk limits, the provided formulation
is versatile and can be extended to other useful third-order constraints such
as torque rate.

Furthermore, our study has shown a significant reduction in peak power and
torque efficiency from these jerk limits. This comes at the expense of some
increased computation time, though it remains well within real-time limits.
Additionally, this trade-off can be managed by computing the next step while
the robot is executing the preceding path and using warm-starting to recompute
solutions as needed. Additionally, although jerk limits result in slightly
longer motion durations, this can be compensated by reduced time required for
the robot to stabilize and come to a stop, due to reduced tracking fluctuation
and overshoot.

During this study, we also observed that the results were quite sensitive to
path parameter function discretization and splining. Splines can interpolate
the discretized solution to generate robot commands at the desired high
control frequency, but some spline types may induce jerk discontinuities
between spline segments. This can be mitigated by constraining acceleration
changes in transitions between spline segments, which can be formulated as
additional third-order constraints. Furthermore, it is crucial to ensure
sufficient discretization resolution to avoid high jerk between knot points
where jerk is not directly constrained. In the future, we should explore more
optimal and constrained parameterization methods for TOTP’s path parameter
function and refine the overall optimization process to accommodate
corresponding changes effectively.

## APPENDIX

###  V-A Proof of convexity

Proposition 1.
$f(x)=\dfrac{a}{\sqrt{x_{1}}+\sqrt{x_{2}}}+\dfrac{b}{\sqrt{x_{2}}+\sqrt{x_{3}}}$
is convex in $x=(x_{1},x_{2},x_{3})\in\mathbb{R}_{+}^{3}$.

Proof. To show a twice-differentiable function $f$ is convex on a convex set,
it is sufficient to show that its Hessian matrix of second partial derivatives
is positive semi-definite in the interior of the set [27]. The Hessian matrix
is

$\displaystyle H$ $\displaystyle=\frac{\partial^{2}f}{\partial
x^{2}}=\begin{bmatrix}s_{11}&s_{12}&0\\\ s_{12}&s_{22}&s_{23}\\\
0&s_{23}&s_{33}\end{bmatrix},$

where

$\displaystyle s_{11}=$
$\displaystyle\frac{a}{2x_{1}(\sqrt{x_{1}}+\sqrt{x_{2}})^{2}}\left(\frac{1}{2\sqrt{x_{1}}}+\frac{1}{\sqrt{x_{1}}+\sqrt{x_{2}}}\right)$
$\displaystyle s_{22}=$
$\displaystyle\frac{a}{2x_{2}(\sqrt{x_{1}}+\sqrt{x_{2}})^{2}}\left(\frac{1}{\sqrt{x_{1}}+\sqrt{x_{2}}}+\frac{1}{2\sqrt{x_{2}}}\right)$
$\displaystyle+\frac{b}{2x_{2}(\sqrt{x_{2}}+\sqrt{x_{3}})^{2}}\left(\frac{1}{2\sqrt{x_{2}}}+\frac{1}{\sqrt{x_{2}}+\sqrt{x_{3}}}\right)$
$\displaystyle s_{33}=$
$\displaystyle\frac{b}{2x_{3}(\sqrt{x_{2}}+\sqrt{x_{3}})^{2}}\left(\frac{1}{\sqrt{x_{2}}+\sqrt{x_{3}}}+\frac{1}{2\sqrt{x_{3}}}\right)$
$\displaystyle s_{12}=$
$\displaystyle\frac{a}{2\sqrt{x_{1}}\sqrt{x_{2}}(\sqrt{x_{1}}+\sqrt{x_{2}})^{3}},\;s_{23}=\frac{b}{2\sqrt{x_{2}}\sqrt{x_{3}}(\sqrt{x_{2}}+\sqrt{x_{3}})^{3}}.$

As $H$ is a real symmetric matrix, $\lambda$s (eigenvalues of $H$) are real,
which can be obtained from the characteristic equation:

$\displaystyle\lambda^{3}-\text{Tr}(H)\lambda^{2}+\Bigg{(}\sum_{i\neq
j}s_{ii}s_{jj}-s_{12}^{2}-s_{23}^{2}\Bigg{)}\lambda-|H|=0.$

First, given the domain of $x_{i}>0$, it is clear that $\text{Tr}(H)>0$ and
$\sum_{i\neq j}s_{ii}s_{jj}-s_{12}^{2}-s_{23}^{2}>0$. The determinant of $H$
is also positive in the given domain:

$\displaystyle|H|=\dfrac{3ab(a(4\sqrt{x_{2}}\sqrt{x_{3}}+x_{2}+3x_{3})+b(4\sqrt{x_{1}}\sqrt{x_{2}}+3x_{1}+x_{2}))}{64\sqrt{x_{1}^{3}x_{2}^{3}x_{3}^{3}}(\sqrt{x_{1}}+\sqrt{x_{2}})^{4}(\sqrt{x_{2}}+\sqrt{x_{3}})^{4}}>0.$

When considering the relationship between the roots and coefficients of a
cubic equation, it is known that if all coefficients above are positive, then
all $\lambda$s satisfying the equation are also positive. Thus we have shown
that $H$ is positive semi-definite, which implies that $f$ is convex.

## ACKNOWLEDGMENTS

We thank Dexterity and the HCRL personnel for their support of this project.
Jee-eun Lee was a robotics intern for Dexterity, Inc. during the summer of
2023, and Luis Sentis was a consultant for Dexterity, Inc. during the summer
of 2022.

## References

  * [1] T. Chettibi, H. Lehtihet, M. Haddad, and S. Hanchi, “Minimum cost trajectory planning for industrial robots,” _European Journal of Mechanics-A/Solids_ , vol. 23, no. 4, pp. 703–715, 2004.
  * [2] M. Diehl, H. G. Bock, H. Diedam, and P.-B. Wieber, “Fast direct multiple shooting algorithms for optimal robot control,” _Fast Motions in Biomechanics and Robotics_ , pp. 65–93, 2006.
  * [3] J. Schulman, J. Ho, A. X. Lee, I. Awwal, H. Bradlow, and P. Abbeel, “Finding locally optimal, collision-free trajectories with sequential convex optimization.” in _Robotics: Science and Systems_ , vol. 9, no. 1. Berlin, Germany, 2013.
  * [4] Y. Zhao, H.-C. Lin, and M. Tomizuka, “Efficient trajectory optimization for robot motion planning,” in _International Conference on Control, Automation, Robotics and Vision (ICARCV)_. IEEE, 2018, pp. 260–265.
  * [5] X. Zhang, F. Xiao, X. Tong, J. Yun, Y. Liu, Y. Sun, B. Tao, J. Kong, M. Xu, and B. Chen, “Time optimal trajectory planning based on improved sparrow search algorithm,” _Frontiers in Bioengineering and Biotechnology_ , vol. 10, 2022.
  * [6] K. Shin and N. McKay, “Minimum-time control of robotic manipulators with geometric path constraints,” _IEEE Transactions on Automatic Control_ , vol. 30, no. 6, pp. 531–541, 1985.
  * [7] Z. Shiller and S. Dubowsky, “Robot path planning with obstacles, actuator, gripper, and payload constraints,” _The International Journal of Robotics Research_ , vol. 8, no. 6, pp. 3–18, 1989.
  * [8] ——, “On computing the global time-optimal motions of robotic manipulators in the presence of obstacles,” _IEEE Transactions on Robotics and Automation_ , vol. 7, no. 6, pp. 785–797, 1991.
  * [9] J. E. Bobrow, S. Dubowsky, and J. S. Gibson, “Time-optimal control of robotic manipulators along specified paths,” _The international journal of robotics research_ , vol. 4, no. 3, pp. 3–17, 1985.
  * [10] F. Pfeiffer and R. Johanni, “A concept for manipulator trajectory planning,” _IEEE Journal on Robotics and Automation_ , vol. 3, no. 2, pp. 115–123, 1987.
  * [11] D. Eager, A.-M. Pendrill, and N. Reistad, “Beyond velocity and acceleration: Jerk, snap and higher derivatives,” _European Journal of Physics_ , vol. 37, no. 6, 2016.
  * [12] Q.-C. Pham, “A general, fast, and robust implementation of the time-optimal path parameterization algorithm,” _IEEE Transactions on Robotics_ , vol. 30, no. 6, pp. 1533–1540, 2014.
  * [13] D. Verscheure, B. Demeulenaere, J. Swevers, J. De Schutter, and M. Diehl, “Time-optimal path tracking for robots: A convex optimization approach,” _IEEE Transactions on Automatic Control_ , vol. 54, no. 10, pp. 2318–2327, 2009.
  * [14] K. Hauser, “Fast interpolation and time-optimization with contact,” _The International Journal of Robotics Research_ , vol. 33, no. 9, pp. 1231–1250, 2014.
  * [15] K. Shin and N. McKay, “A dynamic programming approach to trajectory planning of robotic manipulators,” _IEEE Transactions on Automatic Control_ , vol. 31, no. 6, pp. 491–500, 1986.
  * [16] C. Singh, “Optimality conditions in multiobjective differentiable programming,” _Journal of Optimization Theory and Applications_ , vol. 53, pp. 115–123, 1987.
  * [17] R. Haschke, E. Weitnauer, and H. Ritter, “On-line planning of time-optimal, jerk-limited trajectories,” in _2008 IEEE/RSJ International Conference on Intelligent Robots and Systems_. IEEE, 2008, pp. 3248–3253.
  * [18] L. Berscheid and T. Kröger, “Jerk-limited real-time trajectory generation with arbitrary target states,” _arXiv preprint arXiv:2105.04830_ , 2021.
  * [19] H. Pham and Q.-C. Pham, “A new approach to time-optimal path parameterization based on reachability analysis,” _IEEE Transactions on Robotics_ , vol. 34, no. 3, pp. 645–659, 2018.
  * [20] L. Consolini, M. Locatelli, A. Minari, A. Nagy, and I. Vajk, “Optimal time-complexity speed planning for robot manipulators,” _IEEE Transactions on Robotics_ , vol. 35, no. 3, pp. 790–797, 2019.
  * [21] H. Pham and Q.-C. Pham, “On the structure of the time-optimal path parameterization problem with third-order constraints,” in _IEEE International Conference on Robotics and Automation (ICRA)_. IEEE, 2017, pp. 679–686.
  * [22] J. Ma, S. Gao, H. Yan, Q. Lv, and G. Hu, “A new approach to time-optimal trajectory planning with torque and jerk limits for robot,” _Robotics and Autonomous Systems_ , vol. 140, 2021.
  * [23] J. Mattmüller and D. Gisler, “Calculating a near time-optimal jerk-constrained trajectory along a specified smooth path,” _The International Journal of Advanced Manufacturing Technology_ , vol. 45, pp. 1007–1016, 2009.
  * [24] M. Wang, J. Xiao, S. Liu, and H. Liu, “A third-order constrained approximate time-optimal feedrate planning algorithm,” _IEEE Transactions on Robotics_ , vol. 38, no. 4, pp. 2295–2307, 2021.
  * [25] S. Lu, J. Zhao, L. Jiang, H. Liu, _et al._ , “Solving the time-jerk optimal trajectory planning problem of a robot using augmented lagrange constrained particle swarm optimization,” _Mathematical Problems in Engineering_ , vol. 2017, 2017.
  * [26] Q. Zhang, S. Li, and X. Gao, “Practical smooth minimum time trajectory planning for path following robotic manipulators,” in _American Control Conference_. IEEE, 2013, pp. 2778–2783.
  * [27] S. P. Boyd and L. Vandenberghe, _Convex optimization_. Cambridge University Press, 2004.


* [12] P. Drineas, R. Kannan, and M. W. Mahoney.  Fast Monte Carlo Algorithms for Matrices I: Approximating Matrix Multiplication.  SIAM Journal on Computing, 36(1):132–157, 2006.
  * [13] P. Drineas, M. Magdon-Ismail, M. W. Mahoney, and D.P. Woodruff.  Fast approximation of matrix coherence and statistical leverage.  Journal of Machine Learning Research, 13(Dec):3475–3506, 2012.
  * [14] P. Drineas and M. W. Mahoney.  RandNLA: Randomized numerical linear algebra.  Communications of the ACM, 59:80–90, 2016.
  * [15] P. Drineas and M. W. Mahoney.  Lectures on randomized numerical linear algebra.  In M. W. Mahoney, J. C. Duchi, and A. C. Gilbert, editors, The Mathematics of Data, IAS/Park City Mathematics Series, pages 1–48. AMS/IAS/SIAM, 2018.
  * [16] P. Drineas, M. W. Mahoney, and S. Muthukrishnan.  Sampling algorithms for $\ell_{2}$ regression and applications.  In Proceedings of the 17th Annual ACM-SIAM Symposium on Discrete Algorithms, pages 1127–1136, 2006.
  * [17] P. Drineas, M. W. Mahoney, and S. Muthukrishnan.  Relative-error CUR matrix decompositions.  SIAM Journal on Matrix Analysis and Applications, 30:844–881, 2008\.
  * [18] P. Drineas, M. W. Mahoney, S. Muthukrishnan, and T. Sarlos.  Faster least squares approximation.  Numerische Mathematik, 117(2):219–249, 2011.
  * [19] J. Durbin.  Efficient Estimation Of Parameters In Moving-Average Models.  Biometrika, 46(3-4):306–316, 1959.
  * [20] A. Eshragh, S. Alizamir, P. Howley, and E. Stojanovski.  Modeling the dynamics of the COVID-19 population in Australia: A probabilistic analysis.  PLoS ONE, 15(10):e0240153, 2020.
  * [21] A. Eshragh, B. Ganim, T. Perkins, and K. Bandara.  The importance of environmental factors in forecasting Australian power demand.  Environmental Modeling & Assessment, 27:1–11, 2022.
  * [22] A. Eshragh, G. Livingston Jr, T. McCann, and L. Yerbury.  Rollage: Efficient rolling average algorithm to estimate ARMA models for big time series data.  arXiv:2103.09175, 2023.
  * [23] A. Eshragh, Pietro, O.D., and M. Saunders.  Toeplitz least squares problems, fast algorithms and big data.  arXiv:2112.12994, 2021.
  * [24] A. Eshragh, F. Roosta, A. Nazari, and M. W. Mahoney.  LSAR: Efficient leverage score sampling algorithm for the analysis of big time series data.  Journal of Machine Learning Research, 23:1–36, 2022.
  * [25] R. Estrin, D. Orban, and M. Saunders.  LSLQ: An iterative method for linear least-squares with an error minimization property.  SIAM Journal on Matrix Analysis and Applications, 40(1):254–275, 2019.
  * [26] B. Fahimnia, J. Sarkis, and A. Eshragh.  A tradeoff model for green supply chain planning: A leanness-versus-greenness analysis.  Omega, 54:173–190, 2015.
  * [27] D. Fong and M. Saunders.  LSMR: An iterative algorithm for sparse least-squares problems.  SIAM Journal on Scientific Computing, 33(5):2950–2971, 2011.
  * [28] G. Golub and C. Loan.  Matrix Computations.  Johns Hopkins University Press, 1983.
  * [29] J. Hamilton.  Time Series Analysis.  Princeton University Press, Princeton, New Jersey, 1994.
  * [30] E. Hannan and J. Rissanen.  Recursive Estimation of Mixed Autoregressive-Moving Average Order.  Biometrika, 69(1):81–94, 1982.
  * [31] G. Kutyniok, M. Saunders, S. Wright, and O. Yilmaz.  Sparse representations, numerical linear algebra, and optimization.  2014\.
  * [32] G. Li, M.and Miller and R. Peng.  Iterative row sampling.  In 2013 IEEE 54th Annual Symposium on Foundations of Computer Science, pages 127–136. IEEE, 2013.
  * [33] E. Liberty, F. Woolfe, P. Martinsson, V. Rokhlin, and M. Tygert.  Randomized algorithms for the low-rank approximation of matrices.  Proceedings of the National Academy of Sciences, 104(51):20167–20172, 2007.
  * [34] Y. Liu and S. Zhang.  Fast quantum algorithms for least squares regression and statistic leverage scores.  Theoretical Computer Science, 657:38–47, 2017.
  * [35] M. W. Mahoney.  Randomized algorithms for matrices and data.  Foundations and Trends® in Machine Learning, 2011\.
  * [36] M. W. Mahoney and P. Drineas.  CUR matrix decompositions for improved data analysis.  Proc. Natl. Acad. Sci. USA, 106:697–702, 2009.
  * [37] X. Meng and M. W. Mahoney.  Low-distortion subspace embeddings in input-sparsity time and applications to robust linear regression.  In Proceedings of the 45th Annual ACM Symposium on Theory of Computing, pages 91–100, 2013.
  * [38] X. Meng, M.A. Saunders, and M. W. Mahoney.  LSRN: A parallel iterative solver for strongly over-or underdetermined systems.  SIAM Journal on Scientific Computing, 36(2):C95–C118, 2014.
  * [39] R. Murray, J. Demmel, M. W. Mahoney, N. B. Erichson, M. Melnichenko, O. A. Malik, L. Grigori, P. Luszczek, M. Dereziński, M. E. Lopes, T. Liang, H. Luo, and J. Dongarra.  Randomized Numerical Linear Algebra – a perspective on the field with an eye to software.  Technical Report Preprint: arXiv:2302.11474v2, 2023.
  * [40] C. Musco and C. Musco.  Recursive sampling for the Nystrom method.  Advances in neural information processing systems, 30, 2017.
  * [41] J. Nelson and H. Nguyên.  Osnap: Faster numerical linear algebra algorithms via sparser subspace embeddings.  In 2013 ieee 54th annual symposium on foundations of computer science, pages 117–126. IEEE, 2013.
  * [42] C. Paige and M. Saunders.  LSQR: An Algorithm for Sparse Linear Equations and Sparse Least Squares.  ACM Transactions on Mathematical Software, 8(1):43–71, 1982.
  * [43] P. Rousseeuw and M. Hubert.  Robust statistics for outlier detection.  WIREs Data Mining and Knowledge Discovery, 1(1):73–79, 2011.
  * [44] A. Rudi, D. Calandriello, L. Carratino, and L. Rosasco.  On fast leverage score sampling and optimal learning.  Advances in Neural Information Processing Systems, 31, 2018.
  * [45] T. Sarlós.  Improved approximation algorithms for large matrices via random projections.  In Proceedings of the 47th Annual IEEE Symposium on Foundations of Computer Science, pages 143–152, 2006.
  * [46] V. Sharan, P. Gopalan, and U. Wieder.  Efficient anomaly detection via matrix sketching.  Advances in neural information processing systems, 31, 2018.
  * [47] R. Shumway and D. Stoffer.  Time Series Analysis and Its Applications: With R Examples.  Springer Texts in Statistics. Springer International Publishing, 4 edition, 2017.
  * [48] A. Sobczyk and E. Gallopoulos.  Estimating leverage scores via rank revealing methods and randomization.  SIAM Journal on Matrix Analysis and Applications, 42(3):1199–1228, 2021.
  * [49] J. Sybrandt, M. Shtutman, and I. Safro.  Moliere: Automatic biomedical hypothesis generation system.  In Proceedings of the 23rd ACM SIGKDD international conference on knowledge discovery and data mining, pages 1633–1642, 2017.
  * [50] D. Woodruff.  Sketching as a tool for numerical linear algebra.  Foundations and Trends® in Theoretical Computer Science, 2014.
  * [51] J. Yang, X. Meng, and M. W. Mahoney.  Implementing randomized matrix algorithms in parallel and distributed environments.  In Proceedings of the IEEE, pages 58–92, 2016.
  * [52] Q. Zuo and H. Xiang.  A quantum-inspired algorithm for approximating statistical leverage scores.  arXiv preprint arXiv:2111.08915, 2021.


# Jupiter and Saturn as Spectral Analogs for Extrasolar Gas Giants and Brown
Dwarfs

Daniel J. Coulter Department of Physics, University of Idaho, Moscow, ID
83844, USA  Jason W. Barnes Department of Physics, University of Idaho,
Moscow, ID 83844, USA  Jonathan J. Fortney Department of Astronomy &
Astrophysics, University of California Santa Cruz, Santa Cruz, CA 95064, USA

###### Abstract

With the advent of direct imaging spectroscopy, the number of spectra from
brown dwarfs and extrasolar gas giants is growing rapidly. Many brown dwarfs
and extrasolar gas giants exhibit spectroscopic and photometric variability,
which is likely the result of weather patterns. However, for the foreseeable
future, point-source observations will be the only viable method to extract
brown dwarf and exoplanet spectra. Models have been able to reproduce the
observed variability, but ground truth observations are required to verify
their results. To that end, we provide visual and near-infrared spectra of
Jupiter and Saturn obtained from the _Cassini_ VIMS instrument. We disk-
integrate the VIMS spectral cubes to simulate the spectra of Jupiter and
Saturn as if they were directly imaged exoplanets or brown dwarfs. We present
six empirical disk-integrated spectra for both Jupiter and Saturn with phase
coverage of $1.7^{\circ}$ to $133.5^{\circ}$ and $39.6^{\circ}$ to
$110.2^{\circ}$, respectively. To understand the constituents of these disk-
integrated spectra, we also provide end member (single feature) spectra for
permutations of illumination and cloud density, as well as for Saturn’s rings.
In tandem, these disk-integrated and end member spectra provide the ground
truth needed to analyze point source spectra from extrasolar gas giants and
brown dwarfs. Lastly, we discuss the impact that icy rings, such as Saturn’s,
have on disk-integrated spectra and consider the feasibility of inferring the
presence of rings from direct imaging spectra.

Jupiter, Saturn, Extrasolar gaseous planets, Brown Dwarfs, Direct Imaging,
Spectrophotometry

††journal: Astrophysical Journal Supplement

##  1 Introduction

###  1.1 Brown Dwarfs and Gas Giants

Brown dwarfs are substellar objects that are not massive enough to sustain the
fusion of hydrogen into helium but are massive enough to sustain deuterium
fusion ($13-75$ $\mathrm{M_{J}}$). Their temperatures range from over $2000$ K
for the hottest L-class brown dwarfs down to $\sim 250$ K for the coldest
Y-class brown dwarf discovered to date (Luhman, 2014). Due to their relatively
low temperatures, brown dwarfs have atmospheres that allow for the
condensation of molecules into clouds. Hot L-class brown dwarfs ($1300-2200$
K) have atmospheric spectra characterized by strong water vapor opacity,
carbon dioxide, metallic hydrides, alkali metals, and optically think clouds
thought to be composed of silicates and iron (Kirkpatrick et al., 1999; Marley
& Robinson, 2014). For T-class brown dwarfs ($600-1300$ K), these signatures
are dominated by prominent methane and water vapor absorption bands
(Kirkpatrick et al., 1999; Burgasser et al., 2002). One model for the L to T
spectral class transition is marked by a shift from a largely uniform cloud
layer to a mix of cloudy and cloudless regions, which eventually results in
atmospheres potentially devoid of clouds (Burgasser et al., 2002; Marley et
al., 2010). However, this specific model appears inconsistent with recent
observations (e,g, Apai et al., 2013; Lew et al., 2020), although by some
mechanism the T dwarfs are indeed mostly devoid of clouds. The spectra of cool
Y-class brown dwarfs (less than $600$ K) show even more pronounced methane and
water absorption bands and the coolest members even show evidence of ammonia
(Cushing et al., 2011). Extrasolar gas giant planets cover this same
temperature range, and are theorized to show similar molecular and cloud
signatures (Sudarsky et al., 2000, 2003; Fortney et al., 2020).

Understanding how clouds form and evolve is a primary objective of brown dwarf
and extrasolar gas giant research (see Question B3 in the white paper by Apai
et al. (2017)). Clouds are of first-order importance in these atmospheres, as
they dramatically impact radiation transport and the atmospheric depth that
one probes in spectral observations (Marley et al., 1999). The energy radiated
from a brown dwarf’s atmosphere taps into the energy reservoir from past and
present gravitational contraction. The emitted thermal radiation is mediated
by the absorption due to molecules and cloud opacity in brown dwarf
atmospheres. Clouds act as a barrier to this radiation, with patchy clouds
yielding non-uniform emission across the disk. Over the past several years,
time-domain studies of the varied thermal emission of brown dwarfs have
expanded the matured (e.g. Metchev et al., 2015; Zhou et al., 2016; Yang et
al., 2016; Cushing et al., 2016; Apai et al., 2017), solidifying the view
brown dwarfs have non-homogeneous photospheres (Tan & Showman, 2019, 2021).

It is likely that all brown dwarfs are variable to some degree. The most
variable appear to be those at the L-to-T transition, where cloud opacity is
changing dramatically, suggesting a close connection between spatially
heterogeneous cloud opacity and variability. A recent review of time-domain
observations is provided by Biller (2017), and impressive overview of models
and observations is found in Zhang (2020). In a manner analogous to the brown
dwarfs, Jupiter and Saturn also radiate thermal energy due to leftover heat
from their formation (Aumann et al., 1969). These planets also have strongly
non-uniform surfaces, marked by regions of thick and thin clouds. These
atmospheres have been studied in increasing detail over the past several
decades, in particular via space missions. The Cassini Mission flyby of
Jupiter in 2000 and orbit around Saturn during 2004-2017 provided fresh data
from new viewing geometries that allow us a unique probe of the reflected and
thermal flux from these non-homoegenous atmospheres. In this paper we study
these solar system spectra in the context of exoplanets and brown dwarfs to
draw connections between these close-by and more distant objects.

###  1.2 Direct Imaging Spectroscopy

Progress in exoplanet atmosphere characterization via high-contrast imaging
has been challenging due to technological limitations, but refined techniques
and the next era of space telescopes show promise of significant progress in
the field. Recent reviews on direct imaging include Traub & Oppenheimer (2010)
and Biller & Bonnefoy (2018). There is a strong interest in the direct imaging
and subsequent characterization of all varieties of exoplanets. While brown
dwarfs can typically be observed in isolation, without the bright glow of a
parent star, this is not true for exoplanets. The James Webb Space Telescope
(JWST) and Nancy Grace Roman Space Telescope (formerly WFIRST) will be fitted
with coronagraphs, which will enable direct imaging of exoplanets near bright
stars. In addition, the Astro2020 decadal survey recommended a large
$\sim$6-meter class space telescope with the driving science case to image and
characterize exoplanets Earth-exoplanets with habitable zones of nearby
Sunlike stars (National Academies of Sciences, Engineering, and Medicine,
2021). With the influx of direct imaging data expected in the coming decades,
it is imperative that the exoplanet and brown dwarf communities are able to
interpret this data using tools and modeling frameworks trained with ground
truth observations from the solar system.

Trying to assess the reflection and emission from solar system planets in the
exoplanet context has expanded as well. Recently, spectral analysis studies
have been conducted for Earth (Livengood et al., 2011; Jiang et al., 2018; Gu
et al., 2021) and Neptune (Simon et al., 2016), providing solar system analogs
for terrestrial planets and ice giants, respectively. Similar studies have
been done for Jupiter (Gelino & Marley, 2000; Ge et al., 2019), which along
with Saturn, are our only spatially resolvable proxies for brown dwarfs and
extrasolar gas giants. These studies present photometric variability of an
unresolved Jupiter, with the most significant variability in the infrared, the
result of patchy clouds. Our study builds upon these results in four key ways.
First, we include spectra of Saturn in addition to Jupiter, which provides an
additional analog to extrasolar gas giants and brown dwarfs. Second, previous
studies have very low spectral resolution, two bands in Gelino & Marley (2000)
and twelve bands in Ge et al. (2019), so high resolution spectral analogs are
badly needed. We provide disk-integrated spectra of Jupiter and Saturn with
330 spectral channels across the visible and near-infrared. Third, we present
end member spectra of Jupiter and Saturn (see Section 2.4) in addition to
disk-integrated spectra. And finally, we make our spectra publicly available
for use in training atmosphere models.

Table 1: VIMS spectral cubes selected for this study. All data presented here were retrieved from the NASA Planetary Data System (PDS) archive, except the resolution values marked with an asterisk, which were not available in PDS and we calculated ourselves. Color images created from each cube can be found in Figure 1. Target | Cube Name | Acquisition Date | VIS, IR Exposure Duration (s) | Resolution (km/pixel) | Phase Angle (${}^{\circ})$
---|---|---|---|---|---
Jupiter | V1355313318_3 | 2000 Dec 12 | 640, 20 | 9660 | 1.66
| V1353677179_2 | 2000 Nov 23 | 640, 20 | 17,478* | 15.19
| V1357119162_1 | 2001 Jan 2 | 640, 20 | 5011 | 76.96
| V1357335218_1 | 2001 Jan 4 | 640, 20 | 5455 | 89.26
| V1357767306_1 | 2001 Jan 9 | 640, 20 | 6872 | 107.27
| V1359577908_1 | 2001 Jan 30 | 640, 20 | 15,536* | 133.52
Saturn | V1564773481_1 | 2007 Aug 2 | 1280, 20 | 1862 | 39.64
| V1565414411_1 | 2007 Aug 10 | 1280, 20 | 2031 | 49.61
| V1565931915_1 | 2007 Aug 16 | 1280, 20 | 1883 | 57.76
| V1469259344_1 | 2004 Jul 23 | 1280, 120 | 3304 | 91.90
| V1518276136_1 | 2006 Feb 10 | 2560, 40 | 1925 | 101.07
| V1517685790_1 | 2006 Feb 3 | 2560, 40 | 1990 | 110.20

In this paper we present visual and near infrared disk-integrated spectra of
Jupiter and Saturn at various phase angles, defined hereafter as the sun-
planet-_Cassini_ angle. In Section 2 we discuss the process of selecting and
reducing the spectral cubes used in this study. In Sections 3 and 4 we present
the disk-integrated and end member spectra for Jupiter and Saturn,
respectively. The disk-integrated spectra are proxies for point-source spectra
of extrasolar gas giants obtained via direct imaging. The end member spectra
for different combinations of illumination and cloudiness will be useful when
attempting to deconstruct point-source spectra into constituent atmospheric
features. In Section 4, we also discuss the ability to discern the presence of
icy rings from point source spectra. Lastly, in Section 5, we conclude by
reiterating how this study improves upon previous work in the field and
comment on potential future work in exo-ring detection using the results
presented here.

Figure 1: Color images of all twelve cubes presented in this study with green
mapped to 2.0 $\mu m$, red mapped to 5.0 $\mu m$, and blue mapped to 2.4 $\mu
m$. In this color scheme, green represents reflected light, red represents
thermal emission, and blue indicates the presence of water. The images of
Saturn were created from the cubes before the ring removal algorithm was
implemented. Note the lack of blue in all images except for V1469259344_1, in
which Saturn’s rings, primarily composed of water ice, are clearly visible in
blue.

##  2 Observations & Methods

###  2.1 Cassini VIMS Observations

Spectral data comes from the _Cassini_ Visual and Infrared Mapping
Spectrometer (VIMS) instrument (Brown et al., 2004), which acquired individual
point spectra with 352 spectral channels between 0.35 and 5.2 $\mu m$. VIMS is
composed of two separate optical and detector assemblies – the visual half, or
VIMS-V, and the infrared half VIMS-IR – that employ a common pointing and
spacecraft interface such that they operate as if they were a single
instrument. Using a targeting mirror, VIMS was able to build spectral image
cubes up to 64x64 pixels in size, providing adequate spatial resolution for
this study. _Cassini_ VIMS obtained over 14,000 cubes during its flyby of
Jupiter and over 550,000 cubes of Saturn’s disk during its 13 years in orbit.
We compiled cubes suitable for this study based on three criteria:

  1. (i)

The image showed the entire disk of the target planet in order to obtain
complete disk-integrated spectra

  2. (ii)

There were no pixels on the target planet’s disk that were saturated due to
overexposure

  3. (iii)

There were no objects or artifacts that would affect the spectra, such as a
transiting moon or severe cosmic ray spikes

Under these constraints, the number of suitable cubes for both Jupiter and
Saturn dropped to only a few dozen each, owing to the frequent overexposure of
some wavelengths and the nature of _Cassini_ ’s orbit, which rarely allowed
for full disk imaging of Saturn. These remaining cubes were clustered in
groups of images taken in quick succession at similar phase angles, with many
phase angles absent. Thus, full phase angle coverage is not possible without
relaxing the selection criteria. In the end, six cubes each were selected for
Jupiter and Saturn so that they provided the most complete phase angle
coverage possible for each planet and for the two combined. Color images of
the cubes selected for this study are shown in Figure 1 and important
observational data can be found in Table 1.

###  2.2 Data Reduction

We retrieved raw VIMS cubes from the NASA Planetary Data System (PDS) archive
and reduced them using the pipeline described in Barnes et al. (2013). The
processing proceeds as follows:

  1. 1.

On-board the spacecraft, for each line VIMS observes it measures a background
level. Raw cubes as downloaded from the spacecraft have the average value of
this background measurement already subtracted from them. However since the
on-board backgrounds have proven to be noisy, we first add the mean background
back in to each pixel for the visual half of VIMS and the full measured
background back in for the infrared half. Then for each cube we create a
linear fit to the instrument-measured background observations acquired at the
end of each row in the cube as a function of time and subtract the linearly
interpolated background from each pixel. This method has proven to be much
more robust against cosmic ray spikes and the slowly varying background level
as compared to the on-board subtraction.

  2. 2.

We begin the automatic pipeline process by marking those pixels that have
saturated.

  3. 3.

Next, the pipeline employs a despiking algorithm to identify pixels whose flux
is dominated by cosmic ray hits. We assign the offending pixels an
interpolated value based on a 3-D polynomial fit of the twenty-six nearest
neighbor pixels. Results of the automated despiker compare favorably to
manually despiked cubes for extended objects like Jupiter and Saturn, although
the despike algorithm does incorrectly flag specular reflections off of Titan
lakes as cosmic ray hits (Barnes et al., 2013).

  4. 4.

We divide through by a spatial flatfield that was taken on the ground before
launch at each wavelength. Lacking an instrumental mechanism for generating an
in-flight flatfield, we tried using hundreds of uniformly illuminated science
observation cubes to generate an empirical flatfield, but found the result
inferior to the ground-based calibration.

  5. 5.

We convert to specific energy ($I_{\lambda}$) using the instrument’s spectro-
radiometric response function $R(\lambda)$ in photons per analog-to-digital
data number ($DN$) (Brown et al., 2004):

$I_{\lambda}=\frac{DN}{\tau}\times
R(\lambda)\times\frac{hc}{\lambda}\times\frac{1}{A\Omega\delta\lambda}$ (1)

where $\tau$ is the exposure time in seconds, $h$ is Planck’s constant, $c$ is
the speed of light, $\lambda$ is the wavelength of this particular VIMS
channel, $\delta\lambda$ is the size of the wavelength bin for this VIMS
channel, $A$ is the area of the VIMS mirror, and $\Omega$ is the solid angle
subtended by a pixel in steradians.

  6. 6.

Using SPICE (Acton, 1999) data provided by the NASA Jet Propulsion Laboratory
(JPL), we calculate the latitude, longitude, x- and y-direction resolution,
phase angle, incidence angle, emission angle, and north azimuth of every point
individually. These data are stored in the cube backplane. Since VIMS is a
spot-scanner (Brown et al., 2004), we are able to compensate for spacecraft
pointing changes during an exposure. The SPICE geometry kernel files also
contain knowledge of the spacecraft position relative to the target planet,
which we ultimately use to convert observed specific energy into the net
spectral flux of the planet.

Some aspects of the VIMS team data reduction pipeline were hardcoded for
reducing cubes obtained during _Cassini_ ’s time orbiting Saturn. Therefore,
for the background subtraction in the visible wavelengths we used a separate
but equivalent procedure when reducing Jupiter observations. To calibrate for
visible wavelengths, we calculated the median value of the background (non-
disk) pixels for each spectral band and then subtracted that value from all
pixels in that band.

###  2.3 Disk-Integrated Spectra

Disk integration of Jupiter’s spectral cubes was straightforward. For both
spectral radiance and apparent reflectance (I/F), we simply coadded all pixels
on the disk and divided by the number of pixels coadded. However, Saturn was
more challenging, since information on the geometry of the rings is not
present within the cube backplanes. Therefore, we manually masked out the
pixels containing the rings and their shadows from the image. We replaced each
missing pixel by swapping it with the pixel with the nearest incidence angle
value. This process was done separately for VIMS-V and VIMS-IR since they had
slightly different viewing geometries. Once the rings were synthetically
removed from the images, we disk-integrated them as we did for Jupiter.
Integrated images of Saturn with rings included were also produced and the
comparison to the “ringless” Saturn is discussed in Section 4. This process
provides spatially resolved examples of a planet with and without rings since
some planets and brown dwarfs will host impressive ring systems while others
will not.

Figure 2: Visual representation of the four end members outlined in Section
2.4. Incident solar light (white) can either be reflected by dense clouds
(green) or cloudless regions (orange). Thermal emission (red) are either
directly emitted or scattered by cloud layers (blue).

###  2.4 End Member Spectra

We obtained end member, or single feature, spectra by selecting an individual
pixel of the target feature and determining its radiance and apparent
reflectance (I/F) values across all VIMS spectral channels, once again
accounting for the difference in viewing geometry between VIMS-V and VIMS-IR.
For each end member, we chose three such pixels and averaged them to reduce
any noise present in a single pixel. Four end members are common to Jupiter
and Saturn:

  1. 1.

Dayside and Cloudy: Regions that are on the illuminated area of the disk and
have thick clouds that both reflect incident light and block thermal emission
from the planet (upper left in Figure 2). On Jupiter, these are the cloud-
dense zones on the illuminated side of the disk (bright green in Figure 1). On
Saturn, this is the illuminated part of the cloudy region straddling its
equator (yellow in Figure 1).

  2. 2.

Dayside and Clear: Regions that are on the illuminated area of the disk and
have few clouds to block thermal emission (lower left in Figure 2). On
Jupiter, these are the cloud-free belts on the illuminated side of the disk
(red on the illuminated side in Figure 1). On Saturn, this is the darker
region just above the bright equatorial zone (reddish-brown in Figure 1).

  3. 3.

Nightside and Cloudy: Regions that are on the dark side of the disk and have
thick clouds that block thermal emission from the planet (upper right in
Figure 2). On Jupiter, these are the cloud-dense zones on the dark side of the
disk (black in Figure 1). On Saturn, this is the dark part of the cloudy
region straddling its equator (red in Figure 1).

  4. 4.

Nightside and Clear: Regions that are on the dark side of the disk and have
fewer clouds to block thermal emission (lower right in Figure 2). On Jupiter,
these are the belts on the dark side of the disk (red on the illuminated side
in Figure 1), which have thinner clouds than the zones. On Saturn, this is the
darker region just above the equatorial zone (dark red in Figure 1).

For Saturn, we produced two additional end member spectra for the A Ring and B
Ring. All end member spectra from Jupiter and Saturn were obtained from cubes
V1357335218_1 and V1469259344_1, respectively.

Figure 3: Signal to noise ratio as a function of wavelength for both Jupiter
and Saturn. As expected, the signal to noise ratio is very low in the major
absorption bands. The data used to generate this figure can be found in the
supplementary information.

###  2.5 Signal to Noise Ratio Calculation

Uncertainties in the spectra are dominated by astrophysical aspects (i.e.
angular position, small atmospheric variations) rather than instrumental ones.
Therefore, viewing Jupiter or Saturn at the same phase angle, but at different
times (thus different angular positions) will result in slightly different
spectra. To account for these variations, we averaged the spectra of four
nearby cubes for Jupiter and five such cubes for Saturn. The phase variations
in these sets are $0.43^{\circ}$ for Jupiter and $0.68^{\circ}$ for Saturn. We
calculated the means and standard deviations for both sets and divided the
mean spectra by the standard deviation spectra (taken here to be the noise) to
get a signal to noise ratio for both planets. Our results are shown in Figure
3.

##  3 Jupiter

###  3.1 Results

Figure 4: _Top_ : Disk-integrated spectra of Jupiter in units of flux for six
phase angles ranging from $1.7^{\circ}$ to $133.5^{\circ}$. _Bottom_ : Same as
above, but in units of I/F. Prominent absorption bands are labelled.
Absorption bands in the visual and infrared wavelengths are taken from
Vdovichenko et al. (2021) and Baines et al. (2005), respectively. The source
data for this figure can be found in the supplementary information.

Disk-integrated spectra for six phase angles of Jupiter are shown in Figure 4.
The visual and near infrared (below 3.5 $\mu m$) spectrum of Jupiter is
dominated by reflected solar light with absorption bands primarily from
methane ($CH_{4}$), phosphine ($PH_{3}$), and ammonia ($NH_{3}$). Of note is
the ammonia absorption band at 3.0 $\mu m$, which is much more pronounced in
Jupiter, leading to a single peak structure in the spectrum as opposed to the
dual peak structure seen in Saturn (see Figure 6). Beyond 3.5 $\mu m$ the
influence of solar light wanes dramatically and the spectrum is now primarily
thermal emission. Therefore, the phase-dependent variations are no longer
present in this regime. Despite the thermal radiation being minimal in
comparison to reflected light in flux, it dominates the I/F spectrum (note the
axis changes in both spectra of Figure 4).

Figure 5: Jupiter end member spectra obtained from cube V1357335218_1 (see
Figure 1). The data used to generate this figure can be found in the
supplementary information.

Jupiter end member I/F spectra are shown in Figure 5. In the visual and near-
infrared, both the dayside and nightside spectra are structurally similar,
however their cloudy and clear members diverge dramatically in the mid-
infrared. This dichotomy is the result of the extreme banding present in
Jupiter’s atmosphere, which can be seen visually in Figure 1. Jupiter’s zones
are particularly cloud-dense while the clouds are relatively sparse in their
counterpart, the belts. The result is a significant percentage of thermal
emission in zones are blocked, while thermal energy can more easily escape
from belts. Beyond 4.5 $\mu m$, thermal emission dominates, so there should be
no difference between the clear spectra (dashed lines in Figure 5) in this
wavelength regime.

###  3.2 Discussion

The intensity variations observed in Figure 4 are the combination of phase-
dependent illumination with smaller effects from Jupiter’s scattering phase
function. As the primary purpose of this study is to provide the exoplanet and
brown dwarf communities with disk-integrated spectra as proxies for direct
imaging observations, the data was not altered to account for any of these
effects. A full analysis of the effects of phase-dependent scattering in
Jupiter’s atmosphere is beyond the scope of this study. However, studies by
Mayorga et al. (2016) and Heng & Li (2021) provide a robust analysis of phase-
dependent effects in Jupiter’s atmosphere, specifically in the context of
comparison to direct imaging data.

Over the past several years, atmosphere models that aim to assess the
temperature structure, composition, and reflection/emission spectra of
directly imaged planets, as well as “retrieval” codes that aim to turn
reflection spectra into atmospheric constraints, have been developed in
anticipation of upcoming spectra (Lupu et al., 2016; Nayak et al., 2017;
MacDonald et al., 2018; Damiano et al., 2020; Damiano & Hu, 2020). The disk-
integrated spectra presented here will be vitally important for training such
models with ground truth spectra from spatially resolvable solar system gas
giants. Direct imaging requires that the target planet be further from the
host star than other detection methods. Therefore, direct imaging spectroscopy
is biased towards planets with modest incident fluxes, perhaps $\sim$1-5 AU,
making the cool gas giants in our solar system ideal analogs. Planets as cool
as Jupiter $T_{\rm eff}=125$ K are dominated by gaseous methane and ammonia
clouds (Sudarsky et al., 2000), so abundances of methane and other molecules
can be determined by the depth of their absorption bands in reflection
spectra. Therefore, future observing campaigns should have wavelength coverage
and resolution to observe these absorption bands (e.g. the large methane bands
centered on 2.3 $\mu m$ and 3.5 $\mu m$ and the ammonia band at 3.1 $\mu m$).

Since it emits thermal radiation, Jupiter also acts as a brown dwarf analog
when viewed in the mid-infrared. Specifically, Jupiter is compositionally
similar to ultra-cool Y-class brown dwarfs, so its utility as a brown dwarf
analog increases as we probe cooler and cooler brown dwarfs such as WISE
J035000.32-565830.2 (Leggett et al., 2016), WISE J085510.83-071442.5 (Luhman,
2014; Skemer et al., 2016; Morley et al., 2018), and others (Leggett et al.,
2015). As shown in Figure 4, this thermal emission is most prominent beyond
4.5 $\mu m$, where cool brown dwarfs are also known to be quite bright. The
intensity of the thermal emission is highly dependent upon cloud density,
which is apparent in Figure 5. Jupiter’s emission near 5 $\mu$m emerges almost
exclusively from “holes” in the clouds, where cloud opacity is low. Y-class
brown dwarf atmospheres may exhibit the same banded structure as Jupiter
(Zhang & Showman, 2014), making it an excellent test case for future brown
dwarf observations. In addition, the end member spectra presented here can aid
in constraining the relative amounts of cloudy and cloudless regions in their
atmospheres.

Figure 6: _Top_ : Disk-integrated spectra of Saturn in units of flux for four
phase angles ranging from $39.6^{\circ}$ to $110.2^{\circ}$. Due to the nature
of Cassini’s orbit around Saturn, full disk, full phase spectra were
impossible to obtain. _Bottom_ : Same as above, but in units of I/F. Once
again, prominent absorption bands, courtesy of Vdovichenko et al. (2021) and
Baines et al. (2005), are labelled. The data used to generate this figure, as
well as for two additional phase angles, can be found in the supplementary
information.

##  4 Saturn

###  4.1 Results

Figure 7: _Top_ : Saturn end member spectra obtained from cube V1469259344_1
(see Figure 1). _Bottom_ : Same as above but the y-scale has been adjusted to
show detail. The features in the nightside spectra in the visible wavelengths
are evidence of ringshine from Saturn’s rings. The data used to generate this
figure can be found in the supplementary information.

Disk-integrated spectra for four phase angles of Saturn are shown in Figure 6.
As with Jupiter, the observed intensity variations are the combination of
phase-dependent illumination and scattering phase function. Studies by Tomasko
et al. (1980) and Pérez-Hoyos et al. (2016) provide analysis of Saturn’s
atmospheric phase function. Saturn’s visual and near-infrared spectrum is
dominated by reflected light and many of the same methane, phosphine, and
ammonia absorption bands as Jupiter, while the mid-infrared is dominated by
thermal radiation. Owing to its much hazier atmosphere and cooler temperature,
we do not see the same dramatic spike beyond 4.5 $\mu m$ that we do in
Jupiter’s spectrum. This difference could potentially indicate the severity of
banding due to haziness of brown dwarfs and gas giants from point-source
spectra alone. Saturn end member spectra are shown in Figure 7. Below 4.5 $\mu
m$, the end members resemble those of Jupiter, but beyond 4.5 $\mu m$, we do
not see the same wild divergence between cloudy and clear spectra. This lack
of divergence in the infrared is once again the result of Saturn’s hazy
atmosphere.

Figure 8: Spectral comparison of two ring geometries: edge-on rings from cube
V1518276136_1 (green) and rings at an inclination of approximately
$20^{\circ}$ from cube V1469259344_1 (magenta). Spectra with the rings and
their shadows removed are shown for comparison.

Composed almost entirely of water ice, Saturn’s rings are compositionally
different than Saturn itself. Water has few absorption bands in the visual and
near-infrared, giving the A Ring and B Ring end member spectra a relatively
flat appearance in Figure 7. Therefore, the rings represent a minority share
of the flux outside absorption bands, but are the contributing factor within
them. This result illustrates the distinctive effect rings have on disk-
integrated spectra and is expanded upon in the following section.

###  4.2 Discussion

Saturn, with its impressive ring system, provides an interesting test case for
our ability to infer the presence of icy rings from point-source spectra
alone. During its time in the Saturn system, _Cassini_ observed Saturn from a
range of orbital inclinations, from edge-on rings up to a $20^{\circ}$
inclination. A comparison of edge-on and inclined ring spectra is shown in
Figure 8. With edge-on rings, the spectrum shows little variation from its
ringless counterpart. The slight decrease in the visual and near-infrared is
the consequence of the rings casting shadows on the disk, resulting in less
reflected light. In contrast, the inclined ring spectrum varies significantly
from its ringless counterpart in the visual and near-infrared. The most
significant difference is the drastic percentage-wise increase within the
major absorption bands. This result shows that the rings prevent the planet
from appearing dark in the absorption bands, so the presence of water or ice
can be readily observed in directly imaged brown dwarfs and exoplanets.
Periodic temporal variations in spectra that oscillate between those shown in
Figure 8 would provide strong evidence for the presence of a significant ring
system. Arnold & Schneider (2004) and Dyudina et al. (2005) showed that rings
around exoplanets could be inferred from photometric variations, while this
study shows that they can be found from spectral analysis. Used in tandem,
these techniques could be a powerful tool for discovering exoplanet rings.

In addition, the nightside end member spectra have minor features that
resemble those in the daytime spectra (see the lower pane in Figure 7). We
conclude that these features are evidence of the nightside being illuminated
by ringshine. Such features are not present in the nightside end member
spectra of Jupiter, further supporting this conclusion. While this phenomenon
cannot be directly observed in exoplanets or brown dwarfs, it should be
accounted for when using these spectra for comparison to models.

##  5 Conclusions

Extrasolar gas giants and brown dwarfs are known to exhibit photometric and
spectroscopic variability, with the most significant variations primarily in
the infrared, the result of patchy clouds. Modeling the composition and
evolution of these clouds is a primary research objective in the brown dwarf
and exoplanet communities, but spatially resolvable spectra are needed to
verify current and future models. Such data is surprisingly sparse, so this
study aims to fill that data gap and builds upon previous work in four ways:

  1. (i)

We present spectra of Saturn in addition to Jupiter, providing an additional
analog object

  2. (ii)

Our work has much finer spectral resolution than previous studies

  3. (iii)

We present end member spectra in addition to disk-integrated spectra

  4. (iv)

Our data are publicly available for use in future studies, specifically for
training atmosphere models

The ground truth spectral data of Jupiter and Saturn provided in this study
will act as spatially resolved comparisons to distant directly imaged brown
dwarfs and extrasolar gas giants. This repository of empirical data will help
verify current and future atmospheric models, as well as provide a baseline
for novel instruments such as the Planet as Exoplanet Analog Spectrograph
(PEAS) described in Martin et al. (2020). Disk-integrated spectra are
analogous to current direct imaging spectra and those we will obtain from
future missions, while the end member spectra provide spatial context.

Another interesting result is the potential to infer the presence of icy rings
solely from temporal variations in point source spectra of extrasolar gas
giants. As the purpose of this study is to provide much needed empirical disk-
integrated spectra of the solar system gas giants and not to assess exo-ring
detection methods, the utility of this result was not explored further.
However, with the potential this technique has, especially when used alongside
photometric exo-ring detection methods, further investigation is certainly
merited.

D. C., J. B., and J. F. are all supported by the NASA _Cassini_ Data Analysis
Program, grant number 80NSSC19K0897.

## References

  * Acton (1999) Acton, C. H. 1999, Lunar and Planetary Science Conference
  * Apai et al. (2013) Apai, D., Radigan, J., Buenzli, E., et al. 2013, The Astrophysical Journal, 768, 121
  * Apai et al. (2017) Apai, D., Cowan, N., Kopparapu, R., et al. 2017, Exploring Other Worlds: Science Questions for Future Direct Imaging Missions (EXOPAG SAG15 Report)
  * Arnold & Schneider (2004) Arnold, L., & Schneider, J. 2004, Astronomy & Astrophysics, 420, 1153
  * Aumann et al. (1969) Aumann, H., Gillespie Jr, C., & Low, F. 1969, The Astrophysical Journal, 157, L69
  * Baines et al. (2005) Baines, K. H., Drossart, P., Momary, T. W., et al. 2005, Earth, Moon, and Planets, 96, 119
  * Barnes et al. (2013) Barnes, J. W., Clark, R. N., Sotin, C., et al. 2013, ApJ, 777, 161
  * Biller (2017) Biller, B. 2017, Astronomical Review, 13, 1
  * Biller & Bonnefoy (2018) Biller, B. A., & Bonnefoy, M. 2018, arXiv preprint arXiv:1807.05136
  * Brown et al. (2004) Brown, R. H., Baines, K. H., Bellucci, G., et al. 2004, Space Science Reviews, 115, 111
  * Burgasser et al. (2002) Burgasser, A. J., Marley, M. S., Ackerman, A. S., et al. 2002, The Astrophysical Journal, 571, L151
  * Burgasser et al. (2002) Burgasser, A. J., Kirkpatrick, J. D., Brown, M. E., et al. 2002, ApJ, 564, 421
  * Cushing et al. (2011) Cushing, M. C., Kirkpatrick, J. D., Gelino, C. R., et al. 2011, ApJ, 743, 50
  * Cushing et al. (2016) Cushing, M. C., Hardegree-Ullman, K. K., Trucks, J. L., et al. 2016, The Astrophysical Journal, 823, 152
  * Damiano & Hu (2020) Damiano, M., & Hu, R. 2020, The Astronomical Journal, 159, 175
  * Damiano et al. (2020) Damiano, M., Hu, R., & Hildebrandt, S. R. 2020, The Astronomical Journal, 160, 206
  * Dyudina et al. (2005) Dyudina, U. A., Sackett, P. D., Bayliss, D. D., et al. 2005, The Astrophysical Journal, 618, 973
  * Fortney et al. (2020) Fortney, J. J., Visscher, C., Marley, M. S., et al. 2020, The Astronomical Journal, 160, 288
  * Ge et al. (2019) Ge, H., Zhang, X., Fletcher, L. N., et al. 2019, The Astronomical Journal, 157, 89
  * Gelino & Marley (2000) Gelino, C., & Marley, M. 2000, in From Giant Planets to Cool Stars, Vol. 212, 322
  * Gu et al. (2021) Gu, L., Fan, S., Li, J., et al. 2021, The Astronomical Journal, 161, 122
  * Heng & Li (2021) Heng, K., & Li, L. 2021, The Astrophysical Journal Letters, 909, L20
  * Jiang et al. (2018) Jiang, J. H., Zhai, A. J., Herman, J., et al. 2018, The Astronomical Journal, 156, 26
  * Kirkpatrick et al. (1999) Kirkpatrick, J. D., Reid, I. N., Liebert, J., et al. 1999, ApJ, 519, 802
  * Leggett et al. (2015) Leggett, S., Morley, C. V., Marley, M., & Saumon, D. 2015, The Astrophysical Journal, 799, 37
  * Leggett et al. (2016) Leggett, S. K., Tremblin, P., Saumon, D., et al. 2016, The Astrophysical Journal, 824, 2
  * Lew et al. (2020) Lew, B. W., Apai, D., Marley, M., et al. 2020, The Astrophysical Journal, 903, 15
  * Livengood et al. (2011) Livengood, T. A., Deming, L. D., A’hearn, M. F., et al. 2011, Astrobiology, 11, 907
  * Luhman (2014) Luhman, K. L. 2014, ApJ, 786, L18
  * Lupu et al. (2016) Lupu, R. E., Marley, M. S., Lewis, N., et al. 2016, The Astronomical Journal, 152, 217
  * MacDonald et al. (2018) MacDonald, R. J., Marley, M. S., Fortney, J. J., & Lewis, N. K. 2018, The Astrophysical Journal, 858, 69
  * Marley et al. (1999) Marley, M. S., Gelino, C., Stephens, D., Lunine, J. I., & Freedman, R. 1999, The Astrophysical Journal, 513, 879
  * Marley & Robinson (2014) Marley, M. S., & Robinson, T. D. 2014, arXiv preprint arXiv:1410.6512
  * Marley et al. (2010) Marley, M. S., Saumon, D., & Goldblatt, C. 2010, ApJ, 723, L117
  * Martin et al. (2020) Martin, E. C., Skemer, A. J., Radovan, M. V., et al. 2020, 11447, 135
  * Mayorga et al. (2016) Mayorga, L., Jackiewicz, J., Rages, K., et al. 2016, The Astronomical Journal, 152, 209
  * Metchev et al. (2015) Metchev, S. A., Heinze, A., Apai, D., et al. 2015, The Astrophysical Journal, 799, 154
  * Morley et al. (2018) Morley, C. V., Skemer, A. J., Allers, K. N., et al. 2018, The Astrophysical Journal, 858, 97
  * National Academies of Sciences, Engineering, and Medicine (2021) National Academies of Sciences, Engineering, and Medicine. 2021, Pathways to Discovery in Astronomy and Astrophysics for the 2020s (The National Academies Press)
  * Nayak et al. (2017) Nayak, M., Lupu, R., Marley, M. S., et al. 2017, Publications of the Astronomical Society of the Pacific, 129, 034401
  * Pérez-Hoyos et al. (2016) Pérez-Hoyos, S., Sanz-Requena, J. F., Sánchez-Lavega, A., Irwin, P. G., & Smith, A. 2016, Icarus, 277, 1
  * Simon et al. (2016) Simon, A. A., Rowe, J. F., Gaulme, P., et al. 2016, The Astrophysical Journal, 817, 162
  * Skemer et al. (2016) Skemer, A. J., Morley, C. V., Allers, K. N., et al. 2016, The Astrophysical Journal Letters, 826, L17
  * Sudarsky et al. (2003) Sudarsky, D., Burrows, A., & Hubeny, I. 2003, The Astrophysical Journal, 588, 1121
  * Sudarsky et al. (2000) Sudarsky, D., Burrows, A., & Pinto, P. 2000, The Astrophysical Journal, 538, 885
  * Tan & Showman (2019) Tan, X., & Showman, A. P. 2019, The Astrophysical Journal, 874, 111
  * Tan & Showman (2021) —. 2021, Monthly Notices of the Royal Astronomical Society, 502, 678
  * Tomasko et al. (1980) Tomasko, M., McMillan, R., Doose, L., Castillo, N., & Dilley, J. 1980, Journal of Geophysical Research: Space Physics, 85, 5891
  * Traub & Oppenheimer (2010) Traub, W. A., & Oppenheimer, B. R. 2010, Exoplanets, 111
  * Vdovichenko et al. (2021) Vdovichenko, V. D., Karimov, A. M., Kirienko, G. A., et al. 2021, Solar System Research, 55, 35
  * Yang et al. (2016) Yang, H., Apai, D., Marley, M. S., et al. 2016, The Astrophysical Journal, 826, 8
  * Zhang (2020) Zhang, X. 2020, Research in Astronomy and Astrophysics, 20, 099
  * Zhang & Showman (2014) Zhang, X., & Showman, A. P. 2014, The Astrophysical Journal Letters, 788, L6
  * Zhou et al. (2016) Zhou, Y., Apai, D., Schneider, G. H., Marley, M. S., & Showman, A. P. 2016, The Astrophysical Journal, 818, 176