#  PyRoboCOP: Python-based Robotic Control & Optimization Package for
Manipulation

Arvind U. Raghunathan1, Devesh K. Jha1and Diego Romeres1 1All authors are with
Mitsubishi Electric Research Laboratories (MERL), Cambridge, MA 02139. Email–
{raghunathan,jha<EMAIL_ADDRESS>
###### Abstract

PyRoboCOP is a Python-based package for control, optimization and estimation
of robotic systems described by nonlinear Differential Algebraic Equations
(DAEs). In particular, the package can handle systems with contacts that are
described by complementarity constraints and provides a general framework for
specifying obstacle avoidance constraints. The package performs direct
transcription of the DAEs into a set of nonlinear equations by performing
orthogonal collocation on finite elements. PyRoboCOP provides automatic
reformulation of the complementarity constraints that are tractable to NLP
solvers to perform optimization of robotic systems. The package is interfaced
with ADOL-C [1] for obtaining sparse derivatives by automatic differentiation
and IPOPT [2] for performing optimization. We evaluate PyRoboCOP on several
manipulation problems for control and estimation.

##  I Introduction

Most manipulation applications are characterized by presence of constrained
environments while dealing with challenging underlying phenomena like
unilateral contacts, frictional contacts, impact and deformation [3]. With
this understanding, we present a python-based robotic control and optimization
package (called PyRoboCOP) that allows solution to a large class of
mathematical programs with nonlinear and complementarity constraints. The
current paper and package only considers systems which can be represented by
DAEs. Integration with physics engines is left as a future work as that
requires additional development.

Contact-rich robotic manipulation tasks could be modeled as complementarity
systems. Obtaining a feasible, let alone an optimal trajectory, can be
challenging for such systems. An effective integration of the high-level
trajectory planning in configuration space with physics-based dynamics is
necessary in order to obtain optimal performance of such robotic systems. To
the best of author’s knowledge, none of the existing python-based open-source
optimization packages can provide support for trajectory optimization with
support for complementarity constraints that arise from contact-rich
manipulation and the easy specification of obstacle avoidance constraints.
Such optimization capability is, however, highly desirable to allow easy
solution to optimization problems for a large-class of contact-rich robotic
systems.

In this paper, we present PyRoboCOP , a lightweight but powerful Python-based
package for control and optimization of robotic systems. The formulation in
PyRoboCOP allows us to handle contact and collision avoidance in an unified
manner. PyRoboCOP uses ADOL-C [1] and IPOPT [2] at its backend for automatic
differentiation and optimization respectively. The main features of the
package are:

  * •

Contact modeling by complementarity constraints

  * •

Obstacle avoidance modeling by complementarity constraints

  * •

Automatic differentiation for sparse derivatives

  * •

Support for minimum time problems

  * •

Support for optimization over fixed mode sequence problems with unknown
sequence time horizons

  * •

Support for parameter estimation in linear complementarity systems.

The features described above should convince the reader that PyRoboCOP
addresses lots of important optimization problems for manipulation systems. By
bringing together ADOL-C [1] and IPOPT [2] we believe that PyRoboCOP would be
very useful for real-time model-based control of manipulation systems. Codes
and instructions for installing and using PyRoboCOP could be found in [4]
under a research only license.
Contributions. The main contributions of the paper are:

  1. 1.

We present a python-based package for optimization and control of robotic
systems with contact and collision constraints.

  2. 2.

We evaluate our proposed package, PyRoboCOP , over a range of different
manipulation systems for control and estimation.

For the sake of space, we do not report the formulation of collision avoidance
with complementarity constraints. The interest reader can refer to the arXiv
version of the work [5].

##  II Related Work

Our work is closely related to various optimization techniques proposed to
solve contact-implicit trajectory optimization. Some related examples could be
found in [6, 7, 8, 9, 10, 11, 12, 13]. In a more general setting, our work is
related to trajectory optimization in the presence of non-differentiable
constraints. These problems are common in systems with constraints like non-
penetrability [14], minimum distance (e.g., in collision avoidance) [15], or
in some cases robustness constraints [16].

Some of the existing open-source software for dynamic optimization are
Optimica [17], ACADO Toolkit [18], TACO [19], pyomo.dae [20], Drake [21] and
CasADi [22]. All of the cited software leverage automatic differentiation to
provide the interfaced NLP solvers with first and second-order derivatives.
However, these software do not provide convenient interfaces for handling
complementarity constraints and obstacle-avoidance which are key requirements
in robotic applications. More recently, some packages have been proposed to
perform contact-rich tasks in robotics [23]. However, the solver proposed in
[23] uses DDP-based [24] techniques which suffer from sub-optimality and
difficulty in constraint satisfaction. Compared to most of the other
techniques in open literature, PyRoboCOP implements different formulations for
handling complementarity constraints using NLP solvers robustly.

The optimization method presented in our work is most closely related to the
direct trajectory optimization method for contact-rich systems earlier
presented in [11, 12, 13]. The authors pose contact dynamics as a measure
differential inclusion and employ an augmented Lagrangian to solve the
resulting complementarity constrained optimization problem. [14] handle the
complementarity constraint by relaxing to an inequality and solving using an
active-set solver. In an analogous manner, other optimization packages like
CasADi and Pyomo can also be extended for the solution to trajectory
optimization in presence of complementarity constraints through a similar
reformulation of the complementarity constraints. We provide an adaptive
approach for relaxing the complementarity constraints. Further, we also
provide a novel formulation for trajectory optimization in the presence of
minimum distance constraints for collision avoidance. To the best of our
knowledge, there is no other existing open-source, python-based optimization
package that can handle constraints arising due to frictional contact
interaction and collision avoidance.

##  III Problem Description

PyRoboCOP solves the dynamic optimization problem

$\displaystyle\min\limits_{x,y,u,p}$
$\displaystyle\,\int\limits_{t_{0}}^{t_{f}}c(x(t),y(t),u(t),p)dt+\phi(x(t_{f}),p)$
(1a) s.t.
$\displaystyle\,f(\dot{x}(t),x(t),y(t),u(t),p)=0,\,x(t_{0})=\widehat{x}$ (1b)
$\displaystyle\,([y(t)]_{\sigma_{l,1}}-\nu_{l,1})([y(t)]_{\sigma_{l,2}}-\nu_{l,2})=0\,\forall\,l\in{\cal
L}$ (1c) $\displaystyle\,\underline{x}\leq
x(t)\leq\overline{x},\underline{y}\leq y(t)\leq\overline{y},\underline{u}\leq
u(t)\leq\overline{u}$ (1d)

where $x(t)\in\mathbb{R}^{n_{x}}$, $y(t)\in\mathbb{R}^{n_{y}}$,
$u(t)\in\mathbb{R}^{n_{u}}$, $\dot{x}(t)\in\mathbb{R}^{n_{x}}$,
$p\in\mathbb{R}^{n_{p}}$ are the differential, algebraic, control, time
derivative of differential variables and time-invariant parameters
respectively. The function
$\phi:\mathbb{R}^{n_{x}+n_{p}}\rightarrow\mathbb{R}$ represents Mayer-type
objective function [25] term and is not a function of the entire trajectory.
In addition, $\underline{x},\overline{x}$, $\underline{y},\overline{y}$,
$\underline{u},\overline{u}$ are the lower and upper bounds on the
differential, algebraic and control variables. The initial condition for the
differential variables is $\widehat{x}$. Constraints (1b)-(1c) are the
Differential Algebraic Equations (DAEs) modeling the dynamics of the system
with
$f:\mathbb{R}^{2n_{x}+n_{y}+n_{u}}\rightarrow\mathbb{R}^{n_{x}+n_{y}-n_{c}}$
with $n_{c}=|{\cal L}|$. Each $l\in{\cal L}$ defines a pair of indices
$\sigma_{l,1},\sigma_{l,2}\in\\{1,\ldots,n_{y}\\}$ that specifies the
complementarity constraint between the algebraic variables
$[y(t)]_{\sigma_{l,1}}$ and $[y(t)]_{\sigma_{l,2}}$. In (1c)
$\nu_{l,1},\nu_{l,2}$ correspond to either the lower or upper bounds on the
corresponding algebraic variables. For example, if they are set respectively
to the lower and upper bounds of corresponding algebraic variables then (1c)
in combination with the bounds (1d) model the complementarity constraint

$0\leq[y(t)-\underline{y}]_{\sigma_{l,1}}\perp[\overline{y}-y(t)]_{\sigma_{l,2}}\geq
0.$

The dynamic optimization problem in (1) is transcripted to a NonLinear Program
(NLP) using the Implicit Euler time-stepping scheme. The time interval
$[t_{0},t_{f}]$ is discretized into $N_{e}$ finite elements of width $h_{i}$
such that $\sum_{i=1}^{N_{e}}h_{i}=t_{f}-t_{0}$. Let
$t_{i}=t_{0}+\sum_{i^{\prime}\leq i}h_{i^{\prime}}$ denote the ending time of
the finite elements $i$. The NLP that results is

$\displaystyle\min$
$\displaystyle\,\sum\limits_{i=1}^{N_{e}}h_{i}c(x_{i},y_{i},u_{i-1})+\phi(x_{N_{e}},p)$
(2a) s.t.
$\displaystyle\,f(\dot{x}_{i},x_{i},y_{i},u_{i-1})=0,\,x_{0}=\widehat{x}$ (2b)
$\displaystyle\,([y_{i}]_{\sigma_{l,1}}-\nu_{l,1})([y_{i}]_{\sigma_{l,2}}-\nu_{l,2})=0\,\forall\,l\in{\cal
L}$ (2c) $\displaystyle\,\underline{x}\leq
x_{i}\leq\overline{x},\underline{y}\leq
y_{i}\leq\overline{y},\underline{u}\leq u_{i}\leq\overline{u}$ (2d)
$\displaystyle\,x_{i+1}=x_{i}+h_{i}\dot{x}_{i}$ (2e)

where the decision variables in (2) are $x_{i}$, $\dot{x}_{i}$, $y_{i}$ and
$u_{i}$. The subscript $i$ approximates the value of the corresponding
variable at time $t_{i}$. The constraints (2b)-(2c) are imposed for $i\in{\cal
N}_{e}$. The constraint in (2e) models Implicit Euler time-stepping scheme and
is imposed for $i\in{\cal N}_{e}\setminus\\{N_{e}\\}$. If complementarity
constraints are present then (2) is an instance of a Mathematical Program with
Complementarity Constraints (MPCC).

MPCCs are well known to fail the standard Constraint Qualification (CQ) such
as the Mangasarian Fromovitz CQ (MFCQ), see [26]. Hence, solution of MPCCs has
warranted careful handling of the complementarity constraints when used in
Interior Point Methods for NLP (IPM-NLP) using relaxation [27, 28] or penalty
formulations [29]. In the case of active set methods, the robust solution of
MPCC relies on special mechanism such as the elastic mode [30].

PyRoboCOP implements two possible relaxation schemes for complementarity
constraints

$\displaystyle\alpha_{l}([y_{i}]_{\sigma_{l,1}}-\nu_{l,1})([y_{i}]_{\sigma_{l,2}}-\nu_{l,2})\leq\delta\,\forall\,l\in{\cal
L}$ (3a) $\displaystyle\sum\limits_{l\in{\cal
L}}\alpha_{l}([y_{i}]_{\sigma_{l,1}}-\nu_{l,1})([y_{i}]_{\sigma_{l,2}}-\nu_{l,2})\leq\delta$
(3b)

where $\alpha_{l}=1$ if the involved bounds ($\nu_{l,1},\nu_{l,2}$) are either
both lower or both upper bounds. If one of the bounds is a lower bound and
other is an upper bound then $\alpha_{l}$ is set to $-1$. Note that the choice
of $\alpha_{l}$ ensures that the resulting product is nonnegative whenever
(2d) are satisfied. The first approach relaxes each complementarity constraint
by a positive parameter $\delta$ [27] while the second approach imposes the
relaxation on the summation of all the complementarity constraints over a
finite element $i$ [30].

Figure 1: Workflow in PyRoboCOP. The dynamics provided by the user to create a
MPCC which is then optimized using IPOPT and the gradients are evaluated using
automatic differentiation via ADOL-C.

In addition, we also have flexibility to keep the $\delta$ fixed to a constant
parameter through out the optimization or link this with the barrier parameter
in IPM-NLP [27, 28].

##  IV Software Description

Figure 1 provides a high-level summary of the flow of control in PyRoboCOP.
Detailed descriptions on the software API and classes are available in the
Software Description in [4]. A user provided class specifies the dynamic
optimization problem (1). This is also briefly described in Figure 1. The user
needs to provide the equality constraints for the dynamical system. These
constraints could include the dynamics information for the system, the bounds
on the system state and inputs, and information about complementarity
constraints, if any. Furthermore, a user needs to provide the objective
function, and also has the option to provide derivative information (note the
derivative information is optional). PyRoboCOP expects the user provided class
to implement the following methods in order to formulate an MPCC (or NLP)
(also shown in Figure 1).

  * •

get_info: Returns information on (1) including $n_{d}$, $n_{a}$, $n_{u}$,
$|{\cal L}|$, $n_{p}$, $N_{e}$, $h_{i}$.

  * •

bounds: Returns the lower and upper bounds on the variables $x(t)$,
$\dot{x}(t)$, $y(t)$, $u(t)$ at a time instant $t$.

  * •

initialcondition: Returns the initial conditions for the variables $x(t_{0})$,
i.e. values of the differential variables at initial time instant $t_{0}$.

  * •

initialpoint: Returns the initial guess for the variables $x(t)$,
$\dot{x}(t)$, $y(t)$, $u(t)$ at a time instant $t$. This initial guess is
passed to the NLP solver.

  * •

objective: Implements method to evaluate and return $c(x(t),y(t),u(t),p)$ at a
time instant $t$.

  * •

constraint: Implements method to evaluate and return
$(f(x(t),\dot{x}(t),y(t),u(t),p)$ at a time instant $t$.

We provide a description of optional methods that are expected if certain
specified conditions are satisfied.

  * •

bounds_finaltime: Returns the bounds on the variables $x(t_{f})$ at the final
time. This method allows to specify a final time condition on a subset or all
of the differential variables.

  * •

bounds_params: Returns information on lower and upper bounds on the parameters
$p$. This method must be implemented if $n_{p}>0$.

  * •

initialpoint_params: Returns the initial guess for the parameters $p$. This
method must be implemented if $n_{p}>0$. This initial guess is passed to the
NLP solver.

  * •

get_complementarity_info: Returns information on the complementarity
constraints in (1) i.e. ${\cal L}$ and also information on whether the lower
or upper bound is involved in the complementarity constraint. This method must
be implemented if ${\cal L}\neq\emptyset$.

  * •

objective_mayer: Implements method to evaluate and return $\phi(x(t_{f}),p)$.

  * •

get_objects_info: Returns the information on number of objects $n_{O}$, flags
to indicate if these obstacles are static or dynamic and the number of
vertices $n_{vi}$ for the polytope bounding the objects.

  * •

get_object_vertices: Implements and returns the matrix
$V_{i}(x(t),y(t))\in\mathbb{R}^{3\times n_{vi}}$ representing the vertices of
the polytope bounding the objects. This method is called only when
get_objects_info is implemented and $n_{O}>0$.

PyRoboCOP is interfaced with ADOL-C [1] to compute derivatives (see the
Backend block in Figure 1). Note that the ADOL-C can also provide the sparsity
pattern of the constraint Jacobian and Hessian of the Lagrangian. As mentioned
earlier, the exploitation of sparsity in computations of the NLP is critical
to solve large problems. To provide derivatives PyRoboCOP used ADOL-C to set
up tapes [1] for evaluating: (i) the objective (2a), (ii) constraints
including the DAE (2b) and a reformulation of (2c), and (iii) the Hessian of
the Lagrangian of the NLP (2). The set-up of the tape is done prior to passing
control to the NLP solver. The advantage of this approach is that evaluation
of (2a), (2b) are now C-function calls instead of Python-function calls. This
considerably reduced the time spent in function evaluations for the NLP
solver. As shown in Figure 1, PyRoboCOP uses IPOPT as the optimization solver.
The choice of the solver agrees with our desire of having an open source
software package, for this reason IPOPT is chosen over other competitors like
SNOPT [31] and filterSQP [32].

The user has flexibility in specifying how the complementarity constraints are
solved. The choices are: (i) (3a) with $\delta$ fixed, (ii) (3b) with $\delta$
fixed, (iii) (3a) with $\delta$ set equal to the interior point barrier
parameter, (iv) (3b) with $\delta$ set equal to the interior point barrier
parameter, and (v) the objective function is appended with complementarity
terms as $\sum_{i=1}^{N_{e}}\sum\limits_{l\in{\cal
L}}\alpha_{l}([y_{i}]_{\sigma_{l,1}}-\nu_{l,1})([y_{i}]_{\sigma_{l,2}}-\nu_{l,2})$.
The convergence behavior of formulations can be quite different and we provide
these implementations so the user can choose one that works best for the
problem at hand. As empirical guidelines, we recommend options (iii) and (iv)
since they enforce the complementarity constraints in a gradual manner as the
algorithm converges to a solution.

##  V Numerical results

In this section, we test PyRoboCOP in several robotic simulations providing
solutions to trajectory optimization problems including several systems with
complementarity constraints. To foster reproducibility, the source code for
each of the following examples is available in [4].

(a) Trajectory for $\mathbf{x_{g}}=(0,0.5,\pi)$

(b) Optimal Controls

(c) Trajectory for $\mathbf{x_{g}}=(0,0,\pi)$

(d) Optimal Controls

Figure 2: Optimal pushing sequences and control inputs obtained by solving the
MPCC for two different goal conditions. The switching sequence between
sticking and slipping contact formation could be visualized by the trajectory
of $\dot{p}_{y}$. The pusher maintains a sticking contact with the slider when
$\dot{p}_{y}=0$. For clarity, we show only few frames of the trajectory in
Fig. 2(c). Figure 3: A schematic of cartpole with softwalls and planar pusher-
slider system we study in this paper.

###  V-A Planar Pushing

In this section, we show some results for planar pushing. We briefly describe
the dynamic model here .For more detailed description of the pushing model,
readers are referred to [33] and [34]. A schematic for a pusher-slider system
is shown in Figure 3. The pusher interacts with the slider by exerting forces
in the normal and tangential direction denoted by $f_{\overrightarrow{n}}$,
$f_{\overrightarrow{t}}$ (as shown in Figure 3) as well as a torque $\tau$
about the center of the mass of the object. Assuming quasi-static interaction,
the limit surface [35] defines an invertible relationship between applied
wrench $\mathbf{w}$ and the twist of the slider $\mathbf{t}$. The applied
wrench $\mathbf{w}$ causes the object to move in a perpendicular direction to
the limit surface $\mathbf{H(w)}$. Consequently, the object twist in body
frame is given by $\mathbf{t}=\nabla\mathbf{H(w)}$, where the applied wrench
$\mathbf{w}=[f_{\overrightarrow{n}},f_{\overrightarrow{t}},\tau]$ could be
written as
$\mathbf{w}=\mathbf{J}^{T}(\overrightarrow{n}f_{\overrightarrow{n}}+\overrightarrow{t}f_{\overrightarrow{t}})$.
For the contact configuration shown in Figure 3, the normal and tangential
unit vectors are given by $\overrightarrow{n}=[1\quad 0]^{T}$ and
$\overrightarrow{t}=[0\quad 1]^{T}$.

The equations of motion of the pusher-slider system are

$\mathbf{\dot{x}}=\mathbf{f(\mathbf{x},\mathbf{u})}=\begin{bmatrix}\mathbf{Rt}\quad\dot{p}_{y}\end{bmatrix}^{\top}$
(4)

where $\mathbf{R}$ is the rotation matrix. Since the wrench applied on the
system depends of the point of contact of pusher and slider, the state of the
system is given by $\mathbf{x}=[x\quad y\quad\theta\quad p_{y}]^{T}$ and the
input is given by $\mathbf{u}=[f_{\overrightarrow{n}}\quad
f_{\overrightarrow{t}}\quad\dot{p}_{y}]^{T}$. The elements of the input vector
must follow the laws of coulomb friction which can be expressed as
complementarity conditions as follows:

$\displaystyle
0\leq\dot{p}_{y+}\perp(\mu_{p}f_{\overrightarrow{n}}(t)-f_{\overrightarrow{t}}(t))\geq
0$ $\displaystyle
0\leq\dot{p}_{y-}\perp(\mu_{p}f_{\overrightarrow{n}}(t)+f_{\overrightarrow{t}}(t))\geq
0$ (5)

where $\dot{p}_{y}=\dot{p}_{y+}-\dot{p}_{y-}$ and the $\mu_{p}$ is the
coefficient of friction between pusher and slider. The complementarity
conditions in Eq. (5) mean that both $\dot{p}_{y+}$ and $\dot{p}_{y-}$ are
non-negative and only one of them is non-zero at any time instant and
$\dot{p}_{y}$ is non-zero only at the boundary of friction-cone.

Two pushing trajectories with different goal configurations from the same
initial state are shown in Figure 2. In both these examples, the initial pose
of the slider is $\mathbf{x}_{\texttt{init}}=(0,0,0)$ and the desired goal
pose of the slider is $\mathbf{x_{g}}=(0,0.5,\pi)$ and $(0,0,\pi)$. The
initial point of contact between the pusher and the slider is $p_{y}=0$. For
all these examples, the maximum normal force is set to $0.5$ N and the
coefficient of friction is $\mu_{p}=0.3$. The corresponding control trajectory
shows the sequence of forces $f_{n}$ and $f_{t}$ used by the slider to obtain
the desired trajectory. The plot of $\dot{p}_{y}$ shows the sequence of
sticking and slipping contact as found by PyRoboCOP and thus this also decides
the contact point between the pusher and the slider. Note that the pusher
maintains sticking contact with slider whenever $\dot{p}_{y}=0$, and slipping
contact otherwise.

###  V-B Assembly of Belt Drive Unit

(a) Real Setup of the Belt Drive Unit system.

(b) Visualization of two subtasks decomposition, $S_{1}$ and $S_{2}$. $P_{1}$
and $P_{2}$ are the pulleys.

(c) Optimal trajectory to assemble the belt.

Figure 4: BDU system is shown in (a). In (b) the blue lines represent the belt
gripped by a robot at keypoint $K_{1}$, and $K_{2}$ is the lower keypoint. In
figure (c) the orange points represent the trajectory of the upper keypoint,
$K_{1}$, and the blue points represent the lower keypoint, $K_{2}$, which
together represent the model of the belt. The green and red points are the
starting and the final points, respectively. The belt approaches the first
pulley (not shown), then there is a downwards movement to hook the pulley with
the lower keypoint during $S_{1}$. $K_{1}$ is now hooked onto the pulley and
will not move. Then, the higher keypoint, $K_{2}$, moves toward the second
pulley (not shown) stretching the belt and wraps around the pulley to conclude
$S_{2}$.

(a) Pushing sequence with $\mathbf{x}_{\texttt{init}}$ $(0,0,0)$ and
$\mathbf{x_{g}}=(0.5,0.5,0)$

(b) Optimal Controls obtained for Example 5(a)

(c) Pushing sequence $\mathbf{x}_{\texttt{init}}$ $(0,0,0)$ and
$\mathbf{x_{g}}=(0.45,-0.1,\frac{3\pi}{2})$

(d) Optimal Controls obtained for Example 5(c).

Figure 5: Planar pushing in the presence of obstacles. Our proposed
formulation in PyRoboCOP allows us to solve the collision avoidance. The
trajectory of $\dot{p}_{y}$ shows the slipping contact sequence between the
pusher and the slider. Pusher maintains a sticking contact when
$\dot{p}_{y}=0$.

An example of a complex manipulation problem that involves contacts, elastic
objects and collision avoidance is provided by the Belt Drive Unit system.
This assembly challenge was presented as one of the most challenging
competition in the World Robot Summit
2018111https://worldrobotsummit.org/en/about/ [36]. The real world system is
represented in Figure 4(a) where the objective of the manipulation problem is
to wrap the belt, held by a robotic manipulator around the two pulleys. The
elastic belt is modeled through a 3D keypoint representation. The hybrid
behavior of the model generated by the contacts between the belt and the
pulleys and the elastic properties of the belt is captured by the
complementarity constraints.

The full manipulation task has been divided into two subtasks as shown in
Figure 4(b). The goals of the first and second subtask are to wrap the belt
around the first and second pulley, respectively. We formulate two trajectory
optimization problems one for each of the two subtasks as a MPCC of the form
in (1).

Details on the modeling assumptions, the division into the two subtasks, the
exact formulation including the explanation of the dynamics and
complementarity constraints can be found in our previous paper [37]. In Figure
4(c) we report successful trajectories computed by PyRoboCOP to assemble the
belt drive unit combining the optimal trajectories obtained in the two
subtasks. The optimal trajectory was implemented on the real system with a
tracking controller, see [37] for further details.

###  V-C Planar Pushing With Obstacles

In this section, we show the solution to some planar pushing scenarios in the
presence of obstacles and show that our proposed method can handle
complementarity constraints as well as obstacle avoidance constraints
simultaneously. We demonstrate our approach on two different pushing scenarios
with same initial condition for the slider but different location of the
obstacles and different goal states for the slider. In particular, the initial
state of the slider in both these examples was set to
$\mathbf{x}_{\texttt{init}}=(0,0,0)$ and the goal state for the two conditions
was specified as $\mathbf{x}_{g}=(0.5,0.5,0)$ and $(-0.1,-0.1,3\pi/2)$. We add
the obstacles next to the goal state so that PyRoboCOP has to find completely
different solution compared to the case when there are no obstacles. The
initial point of contact between the pusher and the slider is $p_{y}=0$. The
optimal pushing sequence to reach the goal states for the slider are shown in
Figures 5(a) and 5(c). To provide more insight about the solution, we also
provide the plot of the input sequences in Figures 5(b) and 5(d). The slipping
contact sequence between the slider and the pusher is seen in the plot of
$\dot{p}_{y}$. Sticking contact occurs when $\dot{p}_{y}=0$. We show that the
proposed solver can optimize for the desired sequence of contact modes in
order to reach the target state. For the example in Figure 5(a), the objective
is a function of target state and control inputs. For the example in Figure
5(c), the Mayer objective function is used.

###  V-D Optimization with Mode Enumeration

(a) Optimal pushing sequence.

(b) Optimal control inputs.

Figure 6: Optimal pushing sequence and control inputs obtained by optimizing
mode sequence with PyRoboCOP.

We show our optimization approach over fixed mode sequences using the quasi-
static pushing model presented in Section V-A while considering sticking
contact at the $4$ faces of the slider (see Figure 3). In particular, we use
the dynamics model and the problem described in [38] to show solutions
obtained by PyRoboCOP in the case the mode sequence is pre-specified and the
resulting problem is feasible.

The contact model in this case can be obtained from the model described in
Section V-A, Eq 4 with $\dot{p}_{y}=0$. The modes appear based on which face
the pusher contacts with the slider, and thus we have four different modes
that could be used during any interaction. For a given mode, state-space of
the pusher-slider system is then $3$ dimensional while the input is only $2$
dimensional. The optimization process ensures continuity of dynamics and
selection of final time for each mode in a trajectory. The initial state of
the slider is $\mathbf{x}_{\texttt{init}}=(0,0,0)$ and the goal state of the
slider is $\mathbf{x_{g}}=(0,0,\pi)$. The two modes we use for this example
are pushing from the left face followed by pushing from the top face of the
slider. The trajectory obtained by PyRoboCOP is shown in Figure 6(a). The
inputs used in different modes is shown in Figure 6(b). The objective function
used is the Mayer objective function, i.e., the minimum-time problem. The time
spent in mode $1$ is $12.36$[s] and in mode $2$ is $0.56$[s].

###  V-E System Identification For Complementarity Systems

The system identification problem for systems with complementarity constraints
is a particular case of the optimization problem (2) and therefore can be
solved in PyRoboCOP. The objective is to identify the physical parameters of a
system given a set of collected data. As a case of study, we consider the
cart-pole with softwalls depicted in Figure 3. The interactions of the pole
with the soft walls is modeled as complementarity constrains. The dynamical
equations and more details on the system can be found in [39].

We formalize the parameter estimation problem as MPCC (2) where the cost
function is the normalized Root Mean Square Error (nRMSE) between the observed
trajectory and the trajectory obtained from the estimation procedure. The
parameters we aim to identify are the mass of the pole, $m_{p}$, and the
spring constants of the two walls $k_{1}$ and $k_{2}$. In PyRoboCOP these
parameters are implemented as time independent parameters $p$. We validated
the method with a Monte Carlo simulation on 4 different sets of parameters
$p^{i}=[m_{p}^{i},k_{1}^{i},k_{2}^{i}]$ with $i=\\{1,\ldots,4\\}$ sampled
independently from a uniform distribution, each of which was tested with
different levels of independent Gaussian noise added to the trajectories
collected. The trajectories are generated with an input sequence that is
computed as a sum of sinusoids. Each MC simulation has 50 random realization
of the noise. The results are shown in Figure 7 where on the x-axis we have
the cart-pole system defined with one of the parameter set $p^{i}$ and the
standard deviation of the noise for each system is in order
$[0.0001,0.001,0.01,0.05]$.

Figure 7: Distributions of the estimation errors for 4 different cart-pole
with softwalls with increasing noise level.

We can observe how PyRoboCOP is able to identify both the parameters in the
dynamics equations as well as in the complementarity constraints with the
lower levels of noise. As seen in the figure, we observe that with higher
amounts of noise, the estimation method starts diverging.

##  VI Concluding Remarks

This paper presented PyRoboCOP , a Python-based optimization package for
model-based control of robotic systems. We showed that PyRoboCOP can be used
to solve trajectory optimization problems of a number of dynamical systems in
different configurations in presence of contact and collision avoidance
constraints. A description of the functions that a potential user needs to
implement in order to solve their control or optimization problem has been
provided. More details are available in the software package [4].

Strengths of PyRoboCOP PyRoboCOP can handle systems with contact as well as
collision constraints with a novel complementarity formulation. PyRoboCOP also
allows automatic differentiation by using ADOL-C. To the best of our
knowledge, PyRoboCOP is the only Python-based, open-source software that
allows handling of contact & collision constraints and automatic
differentiation for control and optimization. Unlike most of the competing
optimization solvers which are available in Python, PyRoboCOP allows users to
provide dynamics information in Python through a simple script, refer to the
software description in [4].

Limitations of PyRoboCOP Since PyRoboCOP uses IPOPT as the solver for the
resulting MPCC problems, it borrows limitations of IPOPT. In particular, one
of the main limitations is that only local solutions can be found.
Furthermore, good initialization to find even the local solutions might be
required. Finally, infeasibility of the underlying optimization problem
provided by the user cannot be detected. Another possible limitation is given
by interfacing PyRoboCOP with ADOL-C. While, as described above, this is one
of the strengths of PyRoboCOP it also carries some limitations as we have to
rely on an external code to do automatic differentiation, while other software
like CasADi have built-in source code transformation into C and can handle the
differentiation internally with faster performance.

## References

  * [1] A. Griewank, D. Juedes, and J. Utke, “Algorithm 755: Adol-c: A package for the automatic differentiation of algorithms written in c/c++,” _ACM Trans. Math. Softw._ , vol. 22, no. 2, p. 131–167, Jun. 1996.
  * [2] A. Wächter and L. T. Biegler, “On the implementation of an interior-point filter line-search algorithm for large-scale nonlinear programming,” _Math. Program._ , vol. 106, p. 25–57, 2006.
  * [3] M. T. Mason, “Toward robotic manipulation,” _Annual Review of Control, Robotics, and Autonomous Systems_ , vol. 1, pp. 1–28, 2018.
  * [4] A. U. Raghunathan, D. K. Jha, and D. Romeres. (2021) Pyrobocop. [Online]. Available: https://www.merl.com/research/license/PyRoboCOP
  * [5] ——, “Pyrobocop : Python-based robotic control & optimization package for manipulation and collision avoidance,” 2021.
  * [6] Z. Manchester and S. Kuindersma, “Variational contact-implicit trajectory optimization,” in _Robotics Research_. Springer, 2020, pp. 985–1000.
  * [7] A. Patel, S. L. Shield, S. Kazi, A. M. Johnson, and L. T. Biegler, “Contact-implicit trajectory optimization using orthogonal collocation,” _IEEE Robotics and Automation Letters_ , vol. 4, no. 2, pp. 2242–2249, 2019\.
  * [8] T. Erez and E. Todorov, “Trajectory optimization for domains with contacts using inverse dynamics,” in _2012 IEEE/RSJ International Conference on Intelligent Robots and Systems_. IEEE, 2012, pp. 4914–4919.
  * [9] I. Mordatch, Z. Popović, and E. Todorov, “Contact-invariant optimization for hand manipulation,” in _Proceedings of the ACM SIGGRAPH/Eurographics symposium on computer animation_ , 2012, pp. 137–144.
  * [10] I. Mordatch, E. Todorov, and Z. Popović, “Discovery of complex behaviors through contact-invariant optimization,” _ACM Transactions on Graphics (TOG)_ , vol. 31, no. 4, pp. 1–8, 2012.
  * [11] K. Yunt and C. Glocker, “Trajectory optimization of mechanical hybrid systems using sumt,” in _9th IEEE International Workshop on Advanced Motion Control_ , 2005, p. 665–671.
  * [12] ——, “A combined continuation and penalty method for the determination of optimal hybrid mechanical trajectories,” in _IUTAM Symposium on Dynamics and Control of Nonlinear Systems with Uncertainty_. Netherlands: Springer, 2007, p. 187–196.
  * [13] K. Yunt, “An augmented lagrangian-based shooting method for the optimal trajectory generation of switching lagrangian systems,” _Dynamics of Continuous, Discrete and Impulsive Systems Series B: Applications & Algorithms_, vol. 18, p. 615–645, 2011.
  * [14] M. Posa, C. Cantu, and R. Tedrake, “A direct method for trajectory optimization of rigid bodies through contact,” _The International Journal of Robotics Research_ , vol. 33, no. 1, pp. 69–81, 2014.
  * [15] X. Zhang, A. Liniger, A. Sakai, and F. Borrelli, “Autonomous parking using optimization-based collision avoidance,” in _2018 IEEE Conference on Decision and Control (CDC)_. IEEE, 2018, pp. 4327–4332.
  * [16] P. Kolaric, D. K. Jha, A. U. Raghunathan, F. L. Lewis, M. Benosman, D. Romeres, and D. Nikovski, “Local policy optimization for trajectory-centric reinforcement learning,” in _2020 IEEE International Conference on Robotics and Automation (ICRA)_. IEEE, 2020, pp. 5094–5100.
  * [17] J. Åkesson, K.-E. Årzén, M. Gäfvert, T. Bergdahl, and H. Tummescheit, “Modeling and optimization with optimica and jmodelica.org—languages and tools for solving large-scale dynamic optimization problems,” _Computers & Chemical Engineering_, vol. 34, no. 11, pp. 1737–1749, 2010\. [Online]. Available: https://www.sciencedirect.com/science/article/pii/S009813540900283X
  * [18] B. Houska, H. Ferreau, and M. Diehl, “Acado toolkit – an open source framework for automatic control and dynamic optimization,” _Optimal Control Applications and Methods_ , vol. 32, no. 3, p. 298–312, 2011.
  * [19] S. Leyffer and C. Kirches, “Taco - a toolkit for ampl control optimization,” _Mathematical Programming Computation_ , p. 1–39, 2013.
  * [20] B. Nicholson, J. Siirola, J. Watson, Z. V.M., and L. Biegler, “pyomo.dae: a modeling and automatic discretization framework for optimization with differential and algebraic equations,” _Mathematical Programming Computation_ , vol. 10, p. 187–223, 2018.
  * [21] R. Tedrake and the Drake Development Team, “Drake: Model-based design and verification for robotics,” 2019. [Online]. Available: https://drake.mit.edu
  * [22] J. A. E. Andersson, J. Gillis, G. Horn, J. B. Rawlings, and M. Diehl, “CasADi – A software framework for nonlinear optimization and optimal control,” _Mathematical Programming Computation_ , vol. 11, no. 1, pp. 1–36, 2019.
  * [23] C. Mastalli, R. Budhiraja, W. Merkt, G. Saurel, B. Hammoud, M. Naveau, J. Carpentier, L. Righetti, S. Vijayakumar, and N. Mansard, “Crocoddyl: An efficient and versatile framework for multi-contact optimal control,” in _2020 IEEE International Conference on Robotics and Automation (ICRA)_ , 2020, pp. 2536–2542.
  * [24] D. M. Murray and S. J. Yakowitz, “Constrained differential dynamic programming and its application to multireservoir control,” _Water Resources Research_ , vol. 15, no. 5, pp. 1017–1027, 1979.
  * [25] A. E. Bryson and Y.-C. Ho, _Applied Optimal Control: Optimization, Estimation, and Control_ , 1st ed. Talor & Francis, 1975.
  * [26] Z.-Q. Luo, J.-S. Pang, and D. Ralph, _Mathematical programs with equilibrium constraints_. Cambridge University Press, 1996.
  * [27] A. Raghunathan and L. Biegler, “An interior point method for mathematical programs with complementarity constraints (mpccs),” _SIAM J. Optimization_ , vol. 15, no. 3, pp. 720–750, 2005.
  * [28] A. De Miguel, M. Friedlander, F. Nogales, and S. Scholtes, “A two-sided relaxation scheme for mathematical programs with equilibrium constraints,” _SIAM J Optimization_ , vol. 16, no. 2, p. 587—609, 2005.
  * [29] S. Leyffer, G. López-Calva, and J. Nocedal, “Interior methods for mathematical programs with complementarity constraints,” _SIAM J. Optimization_ , no. 1, p. 52–77, 2006.
  * [30] M. Anitescu, P. Tseng, and S. Wright, “Elastic-mode algorithms for mathematical programs with equilibrium constraints: global convergence and stationarity properties,” _Mathematical Programming_ , vol. 110, no. 2, pp. 337–371, 2007.
  * [31] P. E. Gill, W. Murray, and M. A. Saunders, “Snopt: An sqp algorithm for large-scale constrained optimization,” _SIAM review_ , vol. 47, no. 1, pp. 99–131, 2005.
  * [32] R. Fletcher and S. Leyffer, “Nonlinear programming without a penalty function,” _Mathematical programming_ , vol. 91, no. 2, pp. 239–269, 2002\.
  * [33] F. R. Hogan, E. R. Grau, and A. Rodriguez, “Reactive planar manipulation with convex hybrid mpc,” 2018.
  * [34] M. Bauza, F. R. Hogan, and A. Rodriguez, “A data-efficient approach to precise and controlled pushing,” in _Proceedings of The 2nd Conference on Robot Learning_ , ser. Proceedings of Machine Learning Research, A. Billard, A. Dragan, J. Peters, and J. Morimoto, Eds., vol. 87. PMLR, 29–31 Oct 2018, pp. 336–345. [Online]. Available: http://proceedings.mlr.press/v87/bauza18a.html
  * [35] S. Goyal, A. Ruina, and J. Papadopoulos, “Planar sliding with dry friction part 1. limit surface and moment function,” _Wear_ , vol. 143, no. 2, pp. 307–330, 1991. [Online]. Available: https://www.sciencedirect.com/science/article/pii/0043164891901043
  * [36] F. von Drigalski, C. Schlette, M. Rudorfer, N. Correll, J. Triyonoputro, W. Wan, T. Tsuji, and T. Watanabe, “Robots assembling machines: Learning from the world robot summit 2018 assembly challenge,” 2019.
  * [37] S. Jin, D. Romeres, A. Ragunathan, D. K. Jha, and M. Tomizuka, “Trajectory optimization for manipulation of deformable objects: Assembly of belt drive units,” in _2021 IEEE International Conference on Robotics and Automation (ICRA)_ , 2021. [Online]. Available: https://arxiv.org/abs/2106.00898
  * [38] N. Doshi, F. R. Hogan, and A. Rodriguez, “Hybrid differential dynamic programming for planar manipulation primitives,” 2020.
  * [39] A. Aydinoglu, V. M. Preciado, and M. Posa, “Stabilization of complementarity systems via contact-aware controllers,” _arXiv preprint arXiv:2008.02104_ , 2020.


# Thermal Noise-Limited Laser Stabilization to an 8 mL Volume Fabry-Pérot
Reference Cavity with Microfabricated Mirrors

Charles A. McLemore  Department of Physics, University of Colorado Boulder,
Boulder, CO 80309, USA  National Institute of Standards and Technology,
Boulder, CO 80305, USA<EMAIL_ADDRESS>Naijun Jin  Department
of Applied Physics, Yale University, New Haven, CT 06520, USA  Megan L.
Kelleher  Department of Physics, University of Colorado Boulder, Boulder, CO
80309, USA  National Institute of Standards and Technology, Boulder, CO 80305,
USA  James P. Hendrie  Department of Physics, University of Colorado Boulder,
Boulder, CO 80309, USA  National Institute of Standards and Technology,
Boulder, CO 80305, USA  David Mason  Department of Applied Physics, Yale
University, New Haven, CT 06520, USA  Yizhi Luo  Department of Applied
Physics, Yale University, New Haven, CT 06520, USA  Dahyeon Lee  Department of
Physics, University of Colorado Boulder, Boulder, CO 80309, USA  National
Institute of Standards and Technology, Boulder, CO 80305, USA  Peter Rakich
Department of Applied Physics, Yale University, New Haven, CT 06520, USA
Scott A. Diddams  Department of Physics, University of Colorado Boulder,
Boulder, CO 80309, USA  National Institute of Standards and Technology,
Boulder, CO 80305, USA  Department of Electrical, Computer and Energy
Engineering, University of Colorado Boulder, Boulder, CO 80309, USA  Franklyn
Quinlan  Department of Physics, University of Colorado Boulder, Boulder, CO
80309, USA  National Institute of Standards and Technology, Boulder, CO 80305,
USA<EMAIL_ADDRESS>
###### Abstract

Lasers stabilized to vacuum-gap Fabry-Pérot optical reference cavities display
extraordinarily low noise and high stability, with linewidths much less than 1
Hz. These lasers can expand into new applications and ubiquitous use with the
development of compact, portable cavities that are manufacturable at scale.
Here we demonstrate an 8 mL volume Fabry-Pérot cavity constructed with mirrors
that are fabricated lithographically with finesse near 1 million. A laser
locked to the cavity exhibits phase noise limited by the cavity thermal noise
for offset frequencies ranging from 1 Hz to $\approx$ 1 kHz, with a fractional
frequency stability of 7$\times$10-15 at 1 second. Furthermore, the use of
microfabricated mirrors allows us to expand the design space of centimeter-
scale cavities, and we explore the noise implications of pushing towards
cavity volumes of 2 mL or less.

##  I Introduction

Low-noise, highly stable Fabry-Pérot reference cavities are integral to
ultrastable laser systems used in optical clocks Ludlow  _et al._ (2015);
Oelker  _et al._ (2019); McGrew  _et al._ (2018); Takano  _et al._ (2016),
optically derived low-noise microwaves Fortier  _et al._ (2011); Xie  _et al._
(2017), and gravitational wave detection Willke  _et al._ (2008), with
applications in fundamental physics Huntemann  _et al._ (2014); Delva  _et
al._ (2017); Derevianko and Pospelov (2014), coherent radar Ghelfi  _et al._
(2014); Scheer and Kurtz (1993), and the redefinition of the SI second Riehle
(2015). Through careful design with ultra-low expansion materials, Fabry-Pérot
reference cavities are engineered to maintain unparalleled length stability
between their end mirrors. Locking a laser to a reference cavity transfers
this length stability to the laser frequency, with state-of-the-art systems
demonstrating laser linewidths below 10 mHz and fractional frequency
instabilities below 10-16 Matei  _et al._ (2017). This level of performance is
achieved by exploiting large cavity mode volumes and/or cryogenic operation,
relegating their use to well-controlled laboratory environments. However,
there are a growing number of applications that demand compact, portable
systems capable of low-noise operation in diverse and unpredictable
environments, ranging from ground-based geodesy Grotti  _et al._ (2018);
Takamoto  _et al._ (2020) and earthquake detection Marra  _et al._ (2018) to
space-based tests of fundamental physics Sanjuan  _et al._ (2019). For these
applications, significant scientific and technical impact can be achieved with
laser instability at least an order of magnitude better than what is
achievable either with an unlocked laser or with a system locked to a
microwave oscillator; in other words, with laser fractional instability below
10-14 for timescales up to $\approx$ 1 second. Moreover, repeatable and
scalable cavity manufacture can enable widespread deployment of ultrastable
laser systems, further expanding their use through integration with other
chip-scale photonic devices.

Approaches to address the need for portable ultrastable lasers can be
separated into two broad categories: dielectric resonators and vacuum-gap
Fabry-Pérot cavities, each of which has its own advantages. While dielectric
resonators offer compact and manufacturable device structures, their frequency
stability is fundamentally limited by thermo-refractive noise Lim  _et al._
(2017), as the optical mode is confined within material at finite temperature,
and they typically have greater susceptibility to long-term thermal drift.
Thus, a significant degree of active temperature control is often necessary to
maintain frequency stability near the thermo-refractive noise limit. In
contrast, vacuum-gap Fabry-Pérot cavities have an optical mode that only
interacts with the mirror surfaces, reducing the influence of stochastic
thermal fluctuations inherent to all materials, and they can be constructed
with low thermal expansion materials that greatly reduce the required
temperature stability. Fig. 1a shows a comparison of the fractional frequency
instability of some of the most stable solid-state resonators Lee  _et al._
(2013); Zhang  _et al._ (2019); Stern  _et al._ (2020); Loh  _et al._ (2020);
Alnis  _et al._ (2011) relative to what has been achieved with compact vacuum-
gap Fabry-Pérot cavities Davila-Rodriguez  _et al._ (2017); Leibrandt,
Bergquist, and Rosenband (2013). The improved performance of Fabry-Pérot
cavities comes at the expense of larger cavity volume, small-batch lap and
polish manufacturing, and individual hand assembly. Thus, the ideal portable
ultrastable laser would inherit the frequency stability of bulk Fabry-Pérot
reference cavities but with compactness and manufacturability associated with
microresonator designs.

Figure 1: (a) Comparison of existing compact frequency reference technologies.
FP: Fabry-Pérot cavity, WGM: whispering-gallery-mode resonator. Bold text in
legend indicates which results relied on vacuum (V) and/or active temperature
control (T). (b) Cross-section of the finalized cavity geometry with an
overall volume of 8 mL. (c) Photograph of the reference cavity. Inset: surface
profile data of the three micromirrors fabricated on a single substrate. Color
scale on inset ranges from 0 $\mu$m (blue) to 0.75 $\mu$m (red).

To this end, here we present a near 1 million finesse Fabry-Pérot reference
cavity constructed with lithographically fabricated micromirrors. The
fabrication process allows for highly parallel manufacture of mirrors on a
single substrate with user-defined radius of curvature. The constructed
cavity, shown in schematically in Fig. 1b and pictured in Fig. 1c, has a
volume of only 8 mL. We confirm the cavity thermal noise is not adversely
affected by the manufacturing process, with the phase noise of a laser locked
to the cavity operating at the predicted thermal noise limit out to $\approx$
1 kHz offset. The corresponding fractional frequency instability is close to
the thermal limit of 5$\times$10-15 and is better than 10-14 up to several
seconds of averaging. We use finite element simulation to further explore the
design space, predicting thermal noise and cavity thermal expansion as a
function of mirror substrate thickness. This demonstration of a high-
performance optical reference cavity using microfabricated mirrors represents
an essential step towards the integration of ultrastable frequency references
in numerous field-based experiments.

##  II Noise Modeling & Cavity Design

_Mirror design & fabrication._ Several compact Fabry-Pérot cavity designs have
been proposed and demonstrated to specifically address the need for
portability while meeting high stability requirements Davila-Rodriguez  _et
al._ (2017); Leibrandt, Bergquist, and Rosenband (2013); Didier  _et al._
(2016); Webster and Gill (2011). These cavities exploit a high degree of
mechanical symmetry to reduce acceleration and holding force sensitivity and
use large radius of curvature (ROC) mirrors to reduce the cavity thermal noise
limit. Importantly, these cavities use standard mirrors that are fabricated by
lapping and polishing a curved surface on a glass substrate. To permit bonding
of this substrate to a planar cavity spacer, it is necessary to polish a flat
annulus on the outer rim of the mirror. Both of these features (i.e., curved
mirror and bonding annulus) we instead create using lithographic techniques,
as we describe here.

Central to our cavity is a curved mirror formed using a novel lithographic
technique that utilizes a carefully engineered photoresist reflow and etching
to form mm-scale diameter micromirrors with customizable ROC Kharel  _et al._
(2018); Jin  _et al._ . The process consists of patterning a disk of
photoresist onto a super-polished glass substrate, which is then exposed to a
solvent vapor reflow in a custom chamber. Gradually, the solvent vapor is
absorbed into the photoresist, which begins to flow as surface tension
reshapes the disk. The reflow process can be halted at a prescribed time by
baking out the solvent from the photoresist, resulting in a small dimple shape
on the top of the photoresist. This dimple approximates a parabolic surface
near its center and serves as a template for a concave mirror. Fine tuning of
the reflow duration along with the initial dimensions of the photoresist disk
allow for a wide range of mirror curvatures to be achieved with this
technique; parabolic surfaces have been produced with radii of curvature
ranging from less than 1 mm to over 1 m. Following the reflow step, a reactive
ion etch is used to transfer the parabolic micromirror shape into the glass
substrate below with angstrom-level surface roughness. Finally, a highly
reflective dielectric mirror coating is applied via ion beam sputtering to
complete the micromirror. Full details of this process are available in Ref.
Jin  _et al._ , . Since the photoresist is applied lithographically and the
whole substrate is etched at once, this technique allows for a high degree of
parallel fabrication. Furthermore, the bonding annulus can be lithographically
defined as well, leading to a high degree of repeatability in the cavity
assembly.

Since we seek to demonstrate the viability of these microfabricated mirrors
for use in ultrastable frequency references, we maximized the micromirror ROC
to $\approx$ 1 m to correspondingly increase the optical spot size. This
averages thermal fluctuations over a greater area, thereby lowering the
thermal noise floor Levin (1998). Furthermore, to evaluate the repeatability
and parallel manufacture capabilities of the micromirror fabrication
technique, we included 3 individual large ROC micromirrors on a single mirror
substrate, each capable of forming an optical cavity when paired with one
optical flat. To optimize the remainder of the cavity design, we performed
finite element simulations to assess both the thermal noise and temperature
sensitivity of various geometries. With our microfabricated mirrors, the
design space includes mirror substrate thicknesses much less than that
available with standard polished mirrors, as well as substrates with smaller
cross-sectional diameter.

_Noise modeling._ The most crucial considerations in the cavity design are
those that affect the thermal noise floor, which must be below $1.3/f^{3}$
dBc/Hz (where $f$ is the offset frequency) to reach fractional frequency
instability of below 10-14. Thermal noise in Fabry-Pérot reference cavities is
rooted in thermally driven stochastic fluctuations in the cavity mirrors and
spacer. By individually considering how thermal fluctuations affect the length
of the optical axis through different mechanisms, we can calculate the thermal
noise floor as the sum of these contributions. In our analysis, we consider
Brownian noise arising from mechanical damping within the materials Levin
(1998); Kessler, Legero, and Sterr (2012), thermo-elastic noise resulting from
the coupling of thermal fluctuations to a finite coefficient of thermal
expansion (CTE) in the cavity materials Liu and Thorne (2000); Cerdonio  _et
al._ (2001), and thermo-refractive noise due to a temperature-dependent index
of refraction within the mirror coatings Evans  _et al._ (2008). While
analytical models exist to calculate these noise sources, the models often
depend on simplifying assumptions that typically hold for bulk Fabry-Pérot
cavities. However, as we design more compact cavities, some of these
assumptions break down, such as assuming the mirror substrates are large
compared to the characteristic length scale of thermal diffusion, and assuming
the noise from the mirror substrates and spacer can be treated independently.
Following Ref. Kessler, Legero, and Sterr, 2012, we used finite element
analysis (COMSOL) to construct a model of our cavity and apply a force to the
mirrors with a Gaussian profile matching that of our optical mode. After
computing the resulting deformation of the model cavity, we leverage the
fluctuation-dissipation theorem Callen and Welton (1951) to extract the power
spectral density of mirror surface fluctuations. Summing the contributions
from Brownian, thermo-elastic, and thermo-refractive noise sources, we arrive
at the total expected thermal noise of a given cavity design. As pointed out
in Ref. Kessler, Legero, and Sterr, 2012, the Brownian noise from a small
cavity spacer can be grossly underestimated when using simplified analytical
models. We additionally find that, for the thin mirror substrates available
with microfabrication, the Brownian noise of the coating and mirror substrate
are affected by the contact area between the mirror and spacer. A larger
contact area tends to restrict the drumhead-like motion, driven by Brownian
noise, that otherwise increases the mirror surface displacement.

_Cavity thermal expansion._ In addition to calculating the thermal noise floor
of the cavity, we also utilized finite element simulations to determine the
temperature sensitivity of various cavity designs. While stochastic thermal
noise dominates the stability of a Fabry-Pérot reference cavity at short
timescales, temperature drift of the optical cavity length typically dominates
over long timescales. To maintain sufficient fractional frequency stability
over several seconds, it was therefore necessary to reduce the magnitude of
the effective cavity CTE through careful design of the cavity geometry and
materials. We focused our CTE simulations on a cavity structure that includes
our test substrate with 3 micromirrors. This required a total cavity diameter
of 25.4 mm, with a 10 mm diameter bore hole in the center, while the mirror
substrate thickness and material were left as free parameters. We then
simulated the change in length of the optical axis (sampled with a Gaussian
profile matching the optical mode) as the temperature of the cavity was swept.
The slope of this length change with respect to temperature is proportional to
the cavity’s effective CTE Legero, Kessler, and Sterr (2010), which we used to
evaluate the relative temperature sensitivity.

Using noise and temperature sensitivity simulations in conjunction, we
considered several design trade-offs to arrive at a compact, low-noise design.
For example, while using fused silica as the mirror substrate material
provides a lower thermal noise floor (due to the material’s high mechanical
quality-factor $Q_{M}$ and correspondingly low internal damping), it also
increases the magnitude of the effective cavity CTE when compared to other
possible substrate materials like ultra-low expansion (ULE) glass. To
counteract this increased cavity CTE, we explored thinner mirror substrates,
which improve the CTE, although making the substrates too thin results in a
significant increase to the thermal noise (Fig. 2). Furthermore, we found
that, for a given cavity outer diameter, increasing the bond area reduces the
effective cavity CTE. It is worth noting that in future designs, the inclusion
of only a single micromirror will allow for a large reduction in cavity and
bore hole diameter, and therefore total volume, without sacrificing noise
performance. For example, reducing the cavity diameter by half to 12.7 mm
would result in a cavity volume of only 2 mL while maintaining a comparable
level of thermal noise.

Figure 2: (a) Simulated trade-off between fractional frequency instability and
cavity temperature sensitivity for various potential mirror substrate
thicknesses. The CTE zero-crossing of the 10 mm ULE spacer used in the
simulation is 45 ∘C, whereas the zero-crossing of the ULE used in the
experiment is specified at 45 ∘C $\pm$ 5 ∘C.

The ability to lithographically tailor the shape of the substrate surface
allows us to overcome some of the limitations of conventional fabrication
techniques. Considering the importance of the optical contact bond area on
both the noise and cavity CTE, it is worth noting the level of control our
microfabrication technique provides over this parameter. Two factors set the
upper bound on the bond area: the spacer bore hole diameter and the inner
diameter of the mirror substrate annulus. In our experience, standard curved
mirrors created through lapping and polishing display significant variability
in the mirror contact annulus, leading to unpredictability in CTE and thermal
noise when forming bonded cavities. Additionally, large ROC mirrors have a
central recess of only a few micrometers, with smaller mirror diameters
yielding smaller recess for a given ROC. This puts a lower limit on the
diameter of the mirror, since otherwise the annulus polishing will damage the
pristine surface at the mirror’s center. In our microfabrication procedure, an
annular area around the edge of the substrate is protected with a thick
coating of photoresist, blocking the reactive ion etch. This leaves an intact
bonding annulus of precisely controlled area. Furthermore, this annulus
protection scheme guarantees that the surface quality of the bonding area
remains sufficient for optical contact bonding. As we push towards building
more compact cavities with thinner mirror substrates, larger ROC, and smaller
overall diameter, the precise control over bond area that this technique
allows will be essential for repeatably achieving target performance. Combined
with the numerical simulation methods discussed above, these tools will help
realize a new regime of compact cavities that push beyond traditional
geometries where analytical solutions apply.

Implementing this annulus protection scheme, we finalized the cavity design
with a 10 mm long ULE spacer (CTE zero-crossing: 45 ∘C $\pm$ 5 ∘C) to minimize
thermal sensitivity while still maintaining sub-10-14 fractional frequency
instability. The minimum substrate thickness we could readily obtain with
suitably low roughness (for high finesse) and high surface figure (for optical
contact bonding) was 3 mm. This choice of substrate thickness also provides a
good compromise between noise and CTE. The completed cavity design was
assembled by optical contact bonding the fused silica mirror substrates (with
dielectric coating reflectivity >99.999%) to the ULE spacer, resulting in an
overall volume of 8 mL, effective CTE of 1.5$\times$10-7 K-1 at room
temperature, thermal phase noise floor of $-4.4/f^{3}$ dBc/Hz, and fractional
frequency instability floor of $5~{}\times~{}10^{-15}$.

##  III Experimental Results & Discussion

_Cavity finesse._ For a compact cavity of short length, high finesse is
particularly important to achieve a narrow cavity resonance, which enables
greater rejection of electronic noise sources in the laser frequency lock
Black (2001) and improves rejection of frequency instabilities due to residual
amplitude modulation Zhang  _et al._ (2014). The finesse is dependent on the
mirror transmission and excess loss due to scattering, absorption, and beam
clipping Hood, Kimble, and Ye (2001). To maintain high finesse, it is critical
that the micromirror etching process does not degrade the original angstrom-
level roughness of the super-polished substrate, as this could lead to excess
scattering losses. Furthermore, while scattering and absorption are largely
intrinsic to a given mirror, determined by the roughness of the surface and
quality of the dielectric coating, significant additional beam clipping loss
can be introduced when rigidly bonding a cavity. Due to the mm-scale diameter
of these micromirrors, even minor errors in parallelism (greater than 0.1
mrad) between the end faces of the spacer can result in significant finesse
reduction due to beam clipping. Small temperature gradients in the reflow
chamber while forming the micromirror shapes could also introduce asymmetries
and mirror tilt, which, once bonded to a rigid spacer, would shift the
position of the cavity’s optical mode transversely, incurring clipping losses
as the optical mode overlaps the edge of the micromirror.

In light of these considerations, we ensured that the bonded cavity maintained
consistently high finesse by performing cavity ringdown (CRD) measurements on
all three micromirrors Anderson, Frisch, and Masser (1984). By measuring all
three, we confirmed that the micromirrors exhibit pristine surface quality and
minimal tilt asymmetries with a high degree of repeatability. To perform the
CRD measurements, we swept the frequency of a commercial fiber laser through a
TEM00 mode of the cavity while sending the light through an acousto-optic
modulator (AOM). When a photodetector monitoring the transmission through the
cavity witnessed a spike in intensity, it triggered a high-speed switch, which
cut off the radio-frequency (RF) power to the AOM, blocking the light incident
on the cavity. We then monitored the exponential decay of resonant light out
of the cavity using the transmission detector. With the exponential decay time
constant $\tau$ (equivalent to the photon lifetime) and the cavity length $L$,
we then determined the finesse as $F=\pi\tau c/L$, with $c$ the speed of light
in vacuum. Fig. 3a shows CRD data for all three micromirrors in this cavity,
along with exponential fits, which are used to extract the photon lifetime in
each case. In order to assess the repeatability of the measurement, we
performed twenty CRD measurements on each micromirror. The results are
presented in Fig. 3b, with the average finesse of micromirrors 1, 2, and 3
being 980 000, 900 000, and 920 000, respectively. We separately measured the
coating transmission to be $\approx$ 1.9 ppm Jin  _et al._ , indicating single
mirror loss of the contacted cavity of only $\approx$ 1.3 ppm, $\approx$ 1.6
ppm, and $\approx$ 1.5 ppm for mirrors 1, 2, and 3, respectively. In addition
to the mirrors used to form this cavity, over a dozen other mirrors were
tested without bonding to spacers (allowing for optimal alignment), yielding a
mean finesse value near 1 million and below 1 ppm scattering loss, with the
highest finesse mirrors reaching 1.3 million Jin  _et al._ . These results
confirm that the etching process does not introduce substantial scattering
loss and rigid bonding of the fabricated substrates in a cavity introduces
significantly less than 1 ppm of additional clipping loss per mirror. With
finesse values consistently near 1 million, all three micromirrors are
suitable for locking a laser at the thermal noise limit.

Figure 3: (a) Cavity ringdown measurements for each of the 3 micromirrors
(offset for clarity), including exponential fits used to extract the photon
lifetime $\tau$. (b) Histogram showing the repeatability of the CRD
measurement for each micromirror, with mean values of 900 000, 920 000, and
980 000.

_Phase noise & frequency stability._ Having established the high finesse of
all three micromirrors in the cavity, we chose micromirror 3 for the remaining
measurements since it has the largest ROC at 1.1 m, leading to the lowest
thermal noise limit, as discussed in Section II. However, this thermal noise
limit in part depends on the high mechanical $Q_{M}$ of the underlying fused
silica used as a mirror substrate. While high $Q$ phonon resonances in quartz
have been demonstrated with a similar surface etch process Kharel  _et al._
(2018), and the finesse measurements confirmed that the micromirror
fabrication technique maintains the original surface quality of the substrate,
it was also necessary to measure the fundamental thermal noise of the mirrors
to confirm that the etch process used in the micromirror fabrication did not
significantly degrade the $Q_{M}$ of the substrates. For our cavity, the phase
noise sensitivity to changes in $Q_{M}$ is limited by the coating noise,
estimated to be about 10 dB above the mirror substrate noise. Thus, a tenfold
reduction in $Q_{M}$ would increase the total noise approximately 3 dB.
Coatings with lower thermal noise Cole  _et al._ (2013) or a cavity with a
larger spot size would reduce the separation between coating and substrate
noise and increase the sensitivity to small changes in $Q_{M}$ of the
substrate.

Using the TEM00 mode of the cavity formed by this mirror, we locked the output
frequency of a fiber laser to the cavity using the Pound-Drever-Hall (PDH)
method Drever  _et al._ (1983), depicted in Fig. 4a. The output from a
commercial fiber laser is sent through an AOM and an electro-optic modulator
(EOM), which applies sidebands at 8 MHz. The light is then coupled to the
cavity, which is housed in a rigid holding structure with a heat shield, held
in a vacuum enclosure at 10-5 Pa with active temperature stabilization applied
to the outside of the vacuum enclosure. For further environmental isolation,
the entire apparatus is contained in an acoustic damping box with active
vibration cancellation. The reflected signal from the cavity is photodetected
and demodulated with a mixer to extract an error signal, which is filtered and
amplified to feed back to both the laser and the AOM. Fast feedback is applied
to the AOM via frequency modulation implemented with a voltage-controlled
oscillator (VCO), while low-pass filtered slow feedback is applied to a
piezoelectric tuning port on the laser. To measure the stability of the locked
laser, stabilized light is split off after the AOM and divided into two
channels for comparison with two optical frequency combs stabilized to
independent ultrastable optical references Nakamura  _et al._ (2020).
Heterodyne beat notes from each arm are digitally sampled with software-
defined radio (SDR) at 2 MSa/s and cross-correlated to remove noise not common
to both channels, leaving just the noise of the stabilized light Davila-
Rodriguez  _et al._ (2017); Sherman and Jördens (2016). As Fig. 4b shows, the
phase noise of the stabilized light follows the simulated thermal noise limit
of the cavity out to around 1 kHz offset, confirming that the micromirror
fabrication process does not substantially alter the high mechanical $Q_{M}$
of the fused silica substrate and demonstrating the utility of these
micromirrors for low-noise laser stabilization. The noise at higher offset
frequencies is dominated by residual noise of the fiber laser. Finally, to
document the longer-term stability of the cavity, the Allan deviation is
extracted from the phase record of one of the two heterodyne beat notes
sampled with SDR (Fig. 4c). At 1 second, the Allan deviation is $7\times
10^{-15}$. Following removal of a 1.3 Hz/s linear drift, the fractional
frequency instability is shown to be 10-14 or better out to 7 seconds of
averaging, after which nonlinear drift dominates.

Figure 4: (a) Experimental setup for PDH locking and phase noise measurement.
PDH, Pound-Drever-Hall; VCO, voltage-controlled oscillator; AOM, acousto-optic
modulator; EOM, electro-optic modulator; SDR, software-defined radio. (b)
Phase noise measurement of the stabilized light. Roll-off above 400 kHz is a
measurement artifact. (c) Fractional frequency instability with and without
removal of 1.3 Hz/s linear drift. The effective bandwidth of this measurement
is 2 kHz.

##  IV Conclusion

These results demonstrate ultrastable laser performance from a compact 8 mL
cavity using lithographically fabricated mirrors. Finite element noise
modeling indicates that the volume can be reduced to 2 mL without sacrificing
the noise performance by simply reducing the cavity’s cross-sectional
diameter. With consistently high finesse across all three micromirrors, we
show that this scalable microfabrication technique maintains very high surface
quality, with minimal excess losses due to absorption, scattering, or clipping
from poor alignment tolerances or micromirror shape. The finesse values
obtained for all three micromirrors are consistent with total excess loss of
less than 2 ppm per mirror. Furthermore, the thermal noise limited performance
out to $\approx$ 1 kHz offset confirms that the mirror fabrication technique
does not significantly degrade the fused silica substrate’s naturally high
mechanical $Q_{M}$. Finally, the fractional frequency instability reaches
$6\times 10^{-15}$ between 0.1 s and 1 s, and remains below 10-14 out to
several seconds of averaging, making this cavity useful for portable optical
clock applications.

While recent years have seen impressive efforts to harness some of the
precision of state-of-the-art reference cavities into more portable systems,
these efforts have largely relied on individually polished bulk optics. In
order for a new generation of stable, portable frequency references to
proliferate to a wide range of field applications, it will be necessary to
make many of them, which requires a path towards scalable manufacturing. We
envision the possibility of fabricating many cavities in parallel: by
patterning a grid of micromirrors onto a single substrate, bonding this with a
spacer possessing a matching grid of bore holes, and capping the other end
with a flat mirror substrate, we could simultaneously assemble a large number
of reference cavities. Then, using foundry-scale dicing techniques, the
reference cavities could be divided into individual units. Here, we take a
first step towards the mass manufacture of portable, ultrastable Fabry-Pérot
reference cavities by demonstrating better than 10-14 instability and
thermally limited performance over several decades using micromirrors
fabricated in a highly repeatable lithographic process.

###### Acknowledgements.

We thank Chun-Chia Chen and Amit Agrawal for helpful comments on the
manuscript. We also thank Andrew Ludlow and the NIST Yb clock team for
providing stable reference light. This work was funded by the DARPA A-PhI
program and NIST. This work is a contribution of an agency of the U.S.
government and not subject to copyright in the USA. Commercial equipment is
identified for scientific clarity only and does not represent an endorsement
by NIST.

## Data Availability Statement

The data that support the findings of this study are available from the
corresponding author upon reasonable request.

## References

  * Ludlow  _et al._ (2015) A. D. Ludlow, M. M. Boyd, J. Ye, E. Peik, and P. O. Schmidt, “Optical atomic clocks,” Reviews of Modern Physics 87, 637 (2015).
  * Oelker  _et al._ (2019) E. Oelker, R. Hutson, C. Kennedy, L. Sonderhouse, T. Bothwell, A. Goban, D. Kedar, C. Sanner, J. Robinson, G. Marti, _et al._ , “Demonstration of 4.8$\times 10^{-17}$ stability at 1 s for two independent optical clocks,” Nature Photonics 13, 714–719 (2019).
  * McGrew  _et al._ (2018) W. McGrew, X. Zhang, R. Fasano, S. Schäffer, K. Beloy, D. Nicolodi, R. Brown, N. Hinkley, G. Milani, M. Schioppo, _et al._ , “Atomic clock performance enabling geodesy below the centimetre level,” Nature 564, 87–90 (2018).
  * Takano  _et al._ (2016) T. Takano, M. Takamoto, I. Ushijima, N. Ohmae, T. Akatsuka, A. Yamaguchi, Y. Kuroishi, H. Munekane, B. Miyahara, and H. Katori, “Geopotential measurements with synchronously linked optical lattice clocks,” Nature Photonics 10, 662–666 (2016).
  * Fortier  _et al._ (2011) T. M. Fortier, M. S. Kirchner, F. Quinlan, J. Taylor, J. Bergquist, T. Rosenband, N. Lemke, A. Ludlow, Y. Jiang, C. Oates, _et al._ , “Generation of ultrastable microwaves via optical frequency division,” Nature Photonics 5, 425–429 (2011).
  * Xie  _et al._ (2017) X. Xie, R. Bouchand, D. Nicolodi, M. Giunta, W. Hänsel, M. Lezius, A. Joshi, S. Datta, C. Alexandre, M. Lours, _et al._ , “Photonic microwave signals with zeptosecond-level absolute timing noise,” Nature Photonics 11, 44–47 (2017).
  * Willke  _et al._ (2008) B. Willke, K. Danzmann, M. Frede, P. King, D. Kracht, P. Kwee, O. Puncken, R. Savage, B. Schulz, F. Seifert, _et al._ , “Stabilized lasers for advanced gravitational wave detectors,” Classical and Quantum Gravity 25, 114040 (2008).
  * Huntemann  _et al._ (2014) N. Huntemann, B. Lipphardt, C. Tamm, V. Gerginov, S. Weyers, and E. Peik, “Improved limit on a temporal variation of $m_{p}/m_{e}$ from comparisons of Yb+ and Cs atomic clocks,” Physical Review Letters 113, 210802 (2014).
  * Delva  _et al._ (2017) P. Delva, J. Lodewyck, S. Bilicki, E. Bookjans, G. Vallet, R. Le Targat, P.-E. Pottie, C. Guerlin, F. Meynadier, C. Le Poncin-Lafitte, _et al._ , “Test of special relativity using a fiber network of optical clocks,” Physical Review Letters 118, 221102 (2017).
  * Derevianko and Pospelov (2014) A. Derevianko and M. Pospelov, “Hunting for topological dark matter with atomic clocks,” Nature Physics 10, 933–936 (2014).
  * Ghelfi  _et al._ (2014) P. Ghelfi, F. Laghezza, F. Scotti, G. Serafino, A. Capria, S. Pinna, D. Onori, C. Porzi, M. Scaffardi, A. Malacarne, _et al._ , “A fully photonics-based coherent radar system,” Nature 507, 341–345 (2014).
  * Scheer and Kurtz (1993) J. A. Scheer and J. L. Kurtz,  _Coherent Radar Performance Estimation_ (Artech House, 1993).
  * Riehle (2015) F. Riehle, “Towards a redefinition of the second based on optical atomic clocks,” Comptes Rendus Physique 16, 506–515 (2015).
  * Matei  _et al._ (2017) D. Matei, T. Legero, S. Häfner, C. Grebing, R. Weyrich, W. Zhang, L. Sonderhouse, J. Robinson, J. Ye, F. Riehle, _et al._ , “1.5 $\mu$m lasers with sub-10 mHz linewidth,” Physical Review Letters 118, 263202 (2017).
  * Grotti  _et al._ (2018) J. Grotti, S. Koller, S. Vogt, S. Häfner, U. Sterr, C. Lisdat, H. Denker, C. Voigt, L. Timmen, A. Rolland, _et al._ , “Geodesy and metrology with a transportable optical clock,” Nature Physics 14, 437–441 (2018).
  * Takamoto  _et al._ (2020) M. Takamoto, I. Ushijima, N. Ohmae, T. Yahagi, K. Kokado, H. Shinkai, and H. Katori, “Test of general relativity by a pair of transportable optical lattice clocks,” Nature Photonics 14, 411–415 (2020).
  * Marra  _et al._ (2018) G. Marra, C. Clivati, R. Luckett, A. Tampellini, J. Kronjäger, L. Wright, A. Mura, F. Levi, S. Robinson, A. Xuereb, _et al._ , “Ultrastable laser interferometry for earthquake detection with terrestrial and submarine cables,” Science 361, 486–490 (2018).
  * Sanjuan  _et al._ (2019) J. Sanjuan, K. Abich, M. Gohlke, A. Resch, T. Schuldt, T. Wegehaupt, G. P. Barwood, P. Gill, and C. Braxmaier, “Long-term stable optical cavity for special relativity tests in space,” Optics Express 27, 36206–36220 (2019).
  * Lim  _et al._ (2017) J. Lim, A. A. Savchenkov, E. Dale, W. Liang, D. Eliyahu, V. Ilchenko, A. B. Matsko, L. Maleki, and C. W. Wong, “Chasing the thermodynamical noise limit in whispering-gallery-mode resonators for ultrastable laser frequency stabilization,” Nature Communications 8, 1–7 (2017).
  * Lee  _et al._ (2013) H. Lee, M.-G. Suh, T. Chen, J. Li, S. A. Diddams, and K. J. Vahala, “Spiral resonators for on-chip laser frequency stabilization,” Nature Communications 4, 1–6 (2013).
  * Zhang  _et al._ (2019) W. Zhang, F. Baynes, S. A. Diddams, and S. B. Papp, “Microrod optical frequency reference in the ambient environment,” Physical Review Applied 12, 024010 (2019).
  * Stern  _et al._ (2020) L. Stern, W. Zhang, L. Chang, J. Guo, C. Xiang, M. A. Tran, D. Huang, J. D. Peters, D. Kinghorn, J. E. Bowers, _et al._ , “Ultra-precise optical-frequency stabilization with heterogeneous III–V/Si lasers,” Optics Letters 45, 5275–5278 (2020).
  * Loh  _et al._ (2020) W. Loh, J. Stuart, D. Reens, C. D. Bruzewicz, D. Braje, J. Chiaverini, P. W. Juodawlkis, J. M. Sage, and R. McConnell, “Operation of an optical atomic clock with a Brillouin laser subsystem,” Nature 588, 244–249 (2020).
  * Alnis  _et al._ (2011) J. Alnis, A. Schliesser, C. Y. Wang, J. Hofer, T. J. Kippenberg, and T. Hänsch, “Thermal-noise-limited crystalline whispering-gallery-mode resonator for laser stabilization,” Physical Review A 84, 011804 (2011).
  * Davila-Rodriguez  _et al._ (2017) J. Davila-Rodriguez, F. N. Baynes, A. D. Ludlow, T. M. Fortier, H. Leopardi, S. A. Diddams, and F. Quinlan, “Compact, thermal-noise-limited reference cavity for ultra-low-noise microwave generation,” Optics Letters 42, 1277–1280 (2017).
  * Leibrandt, Bergquist, and Rosenband (2013) D. R. Leibrandt, J. C. Bergquist, and T. Rosenband, “Cavity-stabilized laser with acceleration sensitivity below 10${}^{-12}~{}g^{-1}$,” Physical Review A 87, 023829 (2013).
  * Didier  _et al._ (2016) A. Didier, J. Millo, C. Lacroûte, M. Ouisse, J. Delporte, V. Giordano, E. Rubiola, and Y. Kersalé, “Design of an ultra-compact reference ule cavity,” in  _Journal of Physics: Conference Series_ , Vol. 723 (IOP Publishing, 2016) p. 012029.
  * Webster and Gill (2011) S. Webster and P. Gill, “Force-insensitive optical cavity,” Optics Letters 36, 3572–3574 (2011).
  * Kharel  _et al._ (2018) P. Kharel, Y. Chu, M. Power, W. H. Renninger, R. J. Schoelkopf, and P. T. Rakich, “Ultra-high-Q phononic resonators on-chip at cryogenic temperatures,” APL Photonics 3, 066101 (2018).
  * (30) N. Jin, C. A. McLemore, D. Mason, J. P. Hendrie, Y. Luo, M. Kelleher, P. Kharel, F. Quinlan, S. A. Diddams, and P. Rakich, “Scalable fabrication of ultrahigh-finesse micromirrors with variable geometry,” Submitted for publication.
  * Levin (1998) Y. Levin, “Internal thermal noise in the LIGO test masses: A direct approach,” Physical Review D 57, 659 (1998).
  * Kessler, Legero, and Sterr (2012) T. Kessler, T. Legero, and U. Sterr, “Thermal noise in optical cavities revisited,” JOSA B 29, 178–184 (2012).
  * Liu and Thorne (2000) Y. T. Liu and K. S. Thorne, “Thermoelastic noise and homogeneous thermal noise in finite sized gravitational-wave test masses,” Physical Review D 62, 122002 (2000).
  * Cerdonio  _et al._ (2001) M. Cerdonio, L. Conti, A. Heidmann, and M. Pinard, “Thermoelastic effects at low temperatures and quantum limits in displacement measurements,” Physical Review D 63, 082003 (2001).
  * Evans  _et al._ (2008) M. Evans, S. Ballmer, M. Fejer, P. Fritschel, G. Harry, and G. Ogin, “Thermo-optic noise in coated mirrors for high-precision optical measurements,” Physical Review D 78, 102003 (2008).
  * Callen and Welton (1951) H. B. Callen and T. A. Welton, “Irreversibility and generalized noise,” Physical Review 83, 34 (1951).
  * Legero, Kessler, and Sterr (2010) T. Legero, T. Kessler, and U. Sterr, “Tuning the thermal expansion properties of optical reference cavities with fused silica mirrors,” JOSA B 27, 914–919 (2010).
  * Black (2001) E. D. Black, “An introduction to Pound–Drever–Hall laser frequency stabilization,” American Journal of Physics 69, 79–87 (2001).
  * Zhang  _et al._ (2014) W. Zhang, M. Martin, C. Benko, J. Hall, J. Ye, C. Hagemann, T. Legero, U. Sterr, F. Riehle, G. Cole, _et al._ , “Reduction of residual amplitude modulation to 1$\times 10^{-6}$ for frequency modulation and laser stabilization,” Optics Letters 39, 1980–1983 (2014).
  * Hood, Kimble, and Ye (2001) C. J. Hood, H. Kimble, and J. Ye, “Characterization of high-finesse mirrors: Loss, phase shifts, and mode structure in an optical cavity,” Physical Review A 64, 033804 (2001).
  * Anderson, Frisch, and Masser (1984) D. Z. Anderson, J. C. Frisch, and C. S. Masser, “Mirror reflectometer based on optical cavity decay time,” Applied Optics 23, 1238–1245 (1984).
  * Cole  _et al._ (2013) G. D. Cole, W. Zhang, M. J. Martin, J. Ye, and M. Aspelmeyer, “Tenfold reduction of Brownian noise in high-reflectivity optical coatings,” Nature Photonics 7, 644–650 (2013).
  * Drever  _et al._ (1983) R. Drever, J. L. Hall, F. Kowalski, J. Hough, G. Ford, A. Munley, and H. Ward, “Laser phase and frequency stabilization using an optical resonator,” Applied Physics B 31, 97–105 (1983).
  * Nakamura  _et al._ (2020) T. Nakamura, J. Davila-Rodriguez, H. Leopardi, J. A. Sherman, T. M. Fortier, X. Xie, J. C. Campbell, W. F. McGrew, X. Zhang, Y. S. Hassan, _et al._ , “Coherent optical clock down-conversion for microwave frequencies with 10-18 instability,” Science 368, 889–892 (2020).
  * Sherman and Jördens (2016) J. A. Sherman and R. Jördens, “Oscillator metrology with software defined radio,” Review of Scientific Instruments 87, 054711 (2016).