# Fluid surface self-propulsion via confined Hocking radiation fields

Steven W. Tarr  School of Physics, Georgia Institute of Technology, 837 State
Street, Atlanta, Georgia 30332, USA  Joseph S. Brunner  School of Physics,
Georgia Institute of Technology, 837 State Street, Atlanta, Georgia 30332, USA
Daniel Soto  School of Physics, Georgia Institute of Technology, 837 State
Street, Atlanta, Georgia 30332, USA  Daniel I. Goldman  School of Physics,
Georgia Institute of Technology, 837 State Street, Atlanta, Georgia 30332, USA
<EMAIL_ADDRESS>
(September 3, 2024)

###### Abstract

Diverse physical systems, from maritime to quantum mechanical, experience
forces mediated by asymmetries in surrounding fields. Here we discover a new
phenomenon in which a vertically oscillating, floating robot can be attracted
to or repelled from a boundary via a complex interference of generated and
reflected gravity-capillary (GC) waves, whose dynamics were studied by Hocking
in the 1980s. Force measurements on the robot reveal that attraction increases
as oscillation frequency increases or particle-boundary separation decreases.
Reconstruction of Hocking wave dynamics rationalizes a field-based asymmetry
in wave amplitudes near boundaries that drives robot motion.

Gravity-capillary waves, Boundary effect, Field-mediated locomotion,
Robophysics

††preprint: APS/123-QED

Locomotion via internal driving sufficient to overcome dissipation is well-
studied in diverse systems purcell1977life ; maladen2009undulatory ;
kantsler2013ciliary ; chong2023self ; diaz2022water ; roh2019honeybees ;
lee2019milli ; rhee2022surferbot ; yuan2021wave ; karsai2022real . Inertial
self-propulsion requires asymmetric momentum generation and is typically
produced using changes in shape or mass distribution in the desired direction
of motion. For example, boats can use propellers to expel fluid
gerr1989propeller  while swimmers can propel via traveling waves of body
bending purcell1977life ; maladen2009undulatory ; kantsler2013ciliary ;
diaz2022water ; chong2023self . In both cases, Newton’s Third Law ensures
locomotion.

At fluid surfaces, movement is typically accompanied by wave generation
burghelea2001onset ; lee2019milli ; rhee2022surferbot ; wilcox1979sex ;
roh2019honeybees ; yuan2021wave ; diaz2022water ; benham2022gunwale . Such
wave fields are an important potential source of locomotor benefits, including
asymmetric self-propulsion in animals and robots roh2019honeybees ;
lee2019milli ; rhee2022surferbot , wave-drag reduction of surface-swimming
collectives yuan2021wave , and animal signaling wilcox1979sex . The physics of
coupled waves and sources has been carefully studied and exploited in the
transport of a fluid droplet; interactions with its own wave field
spontaneously produce asymmetric locomotion despite no imposed directional
bias oza2013trajectory ; harris2015pilot ; harris2018interaction ;
saenz2021emergent .

Figure 1: Wave-generating robot boat. (A) Photo of boat generating 17.1 Hz
waves. (B) Schematic of the eccentric motor vibrating the boat to generate
waves; propellers shown in (A) are not used in this study and thus omitted in
(B). (C) Diagram of the tank wherein all experiments were performed. A backlit
checkerboard enables Fast Checkerboard Demodulation (FCD) for spatiotemporal
surface reconstruction wildeman2018real . (Inset) FCD determines fluid surface
height using the instantaneous distortion of a checkerboard by surface
perturbations. (D1-2) Time series of repulsion from (17.1 Hz) and attraction
toward (33.5 Hz) wall, respectively. (E) Perpendicular hull-wall distance for
repulsive and attractive trials at 17.1 Hz.

An interesting question is whether an entity that produces a spatially
symmetric wave field will move on its own. The answer, of course, must be no
for a sufficiently free system. However, in physically confined systems,
emergent wave field asymmetries can produce nonzero net radiation forces on
boundaries; such forces are observed across scales, from the quantum
mechanical vacuum casimir1948influence ; klimchitskaya2009casimir ;
woods2016materials  to fluids boersma1996maritime ; lee2017fluctuation ;
denardo2009water ; parra2014casimir ; ray2014casimir ; ni2015tunable ;
harris2015pilot ; harris2018interaction ; saenz2021emergent . In quantum
mechanics, the Casimir effect demonstrates that nearby neutral plates confine
and modify zero-point-energy wave fields such that they attract one another
casimir1948influence ; klimchitskaya2009casimir ; woods2016materials . In
driven fluid systems with fluctuating surface waves, boundaries generate an
analogous downsampling of wave modes called the “maritime Casimir effect”
boersma1996maritime ; lee2017fluctuation . The downsampled wave modes reduce
the radiation pressure between objects at the fluid surface and can be
observed as reduced amplitude waves boersma1996maritime ; denardo2009water ;
lee2017fluctuation ; parra2014casimir .

In this Letter, we discover that locomotion can be induced in a free-floating,
oscillating robophysical system that does not directly generate an asymmetric
wave field. Symmetrically propagating waves undergo a complex interference
when reflected at a boundary, breaking symmetry and generating propulsive
radiation forces. In lab experiments, we probe these dynamics with a custom-
developed robot and map radiation forces as we vary both oscillation frequency
and confinement distance. We find that confinement on one side leads to a
reduction in wave field amplitude. The dependence of the consequent radiation
force on oscillation frequency can be quantitatively explained by theory for
gravity-capillary (GC) waves developed by Hocking in the 1980s
hocking1987reflection ; hocking1987damping ; hocking1988capillary .

Figure 2: Wave-generating boat experiences attraction and repulsion near
boundaries. (A-B) $\ddot{d}^{\perp}$ versus $\omega$ and $d^{\perp}_{0}$. Red
dotted lines denote the noise floor determined by behavior far from
boundaries. Simultaneous dependence on both parameters is shown in (C), where
each box corresponds to the average of 5 trials.

Apparatus & fundamental behaviors – The robotic boat (total mass $m=368$ g)
consists of a circularly symmetric hull of radius $R_{B}=6$ cm, a custom
circuit board, two fan motors (uxcell Coreless Micro Motor 412), and an
eccentric motor (Vybronics Inc. Cylindrical Vibration Motor VJQ24-35K270B).
The boat’s hull was 3D-printed in Polymaker PolyLite™ PLA and waterproofed
with marine epoxy. All electronics and batteries were mounted onboard, and
additional weights were added such that a free-floating boat at rest is level
to within $1\degree$. We mounted the eccentric motor beneath the electronics;
when enabled, the motor vibrates the boat with power-dependent frequency
$\omega$ primarily along the fore-aft axis (roll) with minimal vertical motion
or induced surface currents. Beyond $\omega=20$ Hz, the vibration tends toward
roll amplitude $0.15\degree\pm 0.02\degree$, pitch (left-right axis) amplitude
$0.05\degree\pm 0.01\degree$, and vertical oscillation amplitude $0.09\pm
0.02$ mm (see SI). The result is a left-right and fore-aft symmetric, radially
emanated ($\overrightharp{k}=k\hat{r}$), monochromatic wave train of
wavelength $\lambda(\omega)$ traveling along the fluid surface (Fig. 1A-B,
Movie S1).

Due to the symmetries of the emitted waves, a boat placed far from boundaries
experiences no net radiation force $F_{W}$. Upon breaking symmetry by
approaching a boundary, $F_{W}$ becomes nonzero, and the boat self-propels
(Figs. 1D1-2, Movie S2). We observe both repulsive and attractive behaviors
(Fig. 1E), with repulsion occurring more weakly such that it is often
indistinguishable from noise.

To probe these dynamics, we placed the boat near a planar boundary extending
from the floor above the water (61 cm long, 30 cm tall, vertical to within
$1\degree$), varied both $\omega$ and initial hull-boundary distance
$d^{\perp}_{0}$, and allowed the boat to move freely in response to $F_{W}$.
Though we were unable to prescribe wave amplitude $A$ independently from
$\omega$, we expect it to affect $F_{W}$ in accord with established theory on
the energy of surface waves lamb1924hydrodynamics ; longuet1964radiation . We
chose a wall with length $\ell\gg R_{B},\lambda$ such that we may treat our
system as quasi-1D and study the boat’s motion along the axis normal to the
wall. Any observed parallel motion had no clear bias. For all experiments, we
programmed a motor controller to ramp the eccentric motor up linearly to the
target frequency over $10$ s to minimize transients.

We recorded images of trials with a Logitech C920 webcam at 30 FPS and tracked
the boat’s lateral motion with color-thresholding code in MATLAB. We extracted
the boat’s perpendicular acceleration $\ddot{d}^{\perp}$ by fitting a
quadratic equation to the position-time data prior to any drag-induced
inflection point. Below a threshold $d^{\perp}_{T}(\omega)$, we observed an
increasingly attractive force with increasing $\omega$ and decreasing
$d^{\perp}_{0}$ (Figs. 2A-C). Near $d^{\perp}_{T}$, a lightly repulsive
$F_{W}$ emerged. Above $d^{\perp}_{T}$, the boat was considered to be “far
from boundaries”; the wave field symmetry was restored and the boat
experienced a near-zero $F_{W}$.

Direct measurement of wave force – Having observed
$\mathcal{O}(\ddot{d}^{\perp})\leq 10^{2}$ $\mu$m/$s^{2}$ across all tested
initial conditions, we sought to isolate $F_{W}$ from any transient effects
(e.g., viscous danov2000viscous  and wave burghelea2001onset ; diaz2022water
drag, inertia scholz2018inertial ) that could dampen the system’s evolution
and result in such a minuscule acceleration. We investigated $F_{W}$ alone by
restricting the boat’s motion to that of a simple pendulum without impeding
vibration (Fig. 3A), a method similarly employed to quantify water wave analog
Casimir forces denardo2009water . The boat was affixed along its central axis
$1.3$ cm above the waterline to a thin fishing line of length $L=1.4$ m via a
bowline knot. We calibrated the line such that when the pendulum angle
$\theta$ was zero, the tension force $F_{T}$ too was zero. For nonzero
$F_{W}$, the boat’s resultant displacement $\Delta x$ caused $F_{T}$ to
increase until reaching force balance (Fig. 3B). We measured $\Delta x$ for a
variety of $\omega$ ($0$-$42$ Hz) and $d^{\perp}_{0}$ ($1.9$-$3.8$ cm) and
observed typical values within $0$-$3$ mm. Since $L\gg\Delta x$, we assume the
boat undergoes negligible vertical displacement.

Figure 3: Hocking radiation field measurements. (A) Force diagram for pendulum
experiments used to directly measure $F_{W}$. (B) Archetypal boat displacement
plots for pendulum experiments at 33.5 Hz. Oscillations are attributed to the
interplay between $F_{W}$ and $F_{T}$. (C-D) $F_{W}$ versus $\omega$ and
$d^{\perp}_{\text{ss}}$. Red dotted lines denote the noise floor determined by
behavior far from boundaries. Simultaneous dependence on both parameters is
shown for 195 trials in (E).

By measuring $\Delta x$ in steady state, we can estimate the perpendicular
wave force $\overrightharp{F}_{W}=(m-\rho V)g\Delta\overrightharp{x}/L$, where
$\rho$ is the fluid density and $V$ is the liquid volume displaced by the
boat. We plotted $F_{W}$ as a function of $\omega$ and the steady-state hull-
wall distance $d^{\perp}_{\text{ss}}$ for 195 trials alongside a few fixed-
parameter slices (Figs. 3C-E). As expected, the qualitative behavior of
$F_{W}$ closely resembles the acceleration heatmap from Fig. 2C, with
increasing attraction below, light repulsion near, and near-zero effects above
$d^{\perp}_{T}(\omega)$. Despite the removal of transient effects, attractive
and repulsive forces remained small, respectively demonstrating
$\mathcal{O}(F_{W})=[10^{1},10^{2}]$ and $[10^{0},10^{1}]$ $\mu$N. We note
that most measurements with $d^{\perp}_{\text{ss}}<d^{\perp}_{T}$ fall outside
the experimental noise floor $2.9\pm 13.1$ $\mu$N.

Surface wave measurements – To better understand the role of the emanated
waves in generating a locomotive force, we employed the synthetic
SyntheticNote  Schlieren visualization technique Fast Checkerboard
Demodulation wildeman2018real  (FCD, see SI) to obtain quantitative
measurements of the wave field (Figs. 1C, 4A-B, Movie S3). For optimal
visualization quality, we minimized the water’s depth to $h_{\text{rest}}=5$
cm for all experiments. Imaging was performed with a high speed camera (AOS
X-PRI) at 500 FPS when the system had reached steady state and processed using
custom MATLAB code derived from Refs. wildeman2018real ; moisy2009synthetic .
We observed the wave train to possess $A\propto r^{-1/2}$ in accord with
established surface wave theory and follow the known dispersion relation for
GC waves:

$\omega^{2}=\left(gk+\frac{\gamma
k^{3}}{\rho}\right)\tanh{\left(h_{\text{rest}}k\right)},$ (1)

where $\gamma$ is the fluid’s surface tension and $g$ is the standard gravity
(Fig. S1C) lamb1924hydrodynamics .

FCD analysis of steady-state waves between the boat and wall reveals a net
field propagating outward from the boat (Fig. 4B). These waves share $\omega$
with those emitted on the boat’s far side and far from boundaries, but possess
reduced $A$ regardless of $\omega$ (Fig. 4C). We surmise that when the boat is
sufficiently close to the wall, reflected waves return with non-negligible
energy and modulate the free surface height at the hull. This modulation
impedes concurrent wave generation on the side nearest the wall while
minimally affecting the opposite side. Consequently, the steady-state
amplitude between the hull and wall is reduced. We liken these dynamics to the
reductions in height when jumping off a deformable medium
aguilar2016robophysical  or pumping a swing bae2005optimal  with poor timing.

Figure 4: Surface wave reconstructions reveal amplitude imbalance near
boundaries. (A-B) Archetypal reconstructions of 17.1 Hz waves far from and
near a boundary, respectively, with space-time heatmaps corresponding to
dotted yellow lines. $\eta(t,\overrightharp{r}$) describes the free surface
height with respect to $h_{\text{rest}}$. Dark gray regions were occupied by
solid objects (e.g., boat, wall). Light gray regions were deemed
unreconstructable (see SI). (C) FCD measurements reveal the net field near the
wall to have reduced $A(\omega)$.

A hydrodynamic model – Armed with an understanding of the wave fields both
near and far from boundaries, we motivate the boat’s locomotive behavior as it
relates to $d^{\perp}$ (Fig. 5A). The existence of a radiative force incident
on a wave emitter and proportional to square amplitude is a classical result
lamb1924hydrodynamics ; longuet1964radiation  observed in many systems with
asymmetric wave generation denardo2009water ; lee2017fluctuation ;
rhee2022surferbot . When the boat is far from boundaries, the generated waves
are spatially symmetric, leading to a net-zero $F_{W}$. Near boundaries,
reflected waves induce an amplitude imbalance, resulting in a finite $F_{W}$.
We postulate that at a perfect $d^{\perp}_{T}$, the reflected wave will have
energy insufficient to modulate the amplitude of emission but sufficient to
lightly force the boat away from the wall.

Though the amplitude dynamics successfully describe the boat’s attraction and
motionlessness for small and large $d^{\perp}$ respectively, they provide
insufficient reasoning for $F_{W}$’s observed frequency dependence. We
hypothesize that the unique properties of GC waves are relevant to these
complex hydrodynamics. Work by Hocking on the interactions of GC waves with
hard surfaces emphasizes the importance of wavenumber (alternatively
frequency) to radiation and reflection hocking1987reflection ;
hocking1987damping ; hocking1988capillary . Upon reflecting off a rigid
boundary, GC waves dissipate substantial energy through complex contact-line
dynamics. Within the accessible wavenumber range for our boat, the reflection
coefficient $R<0.22$ with $R\propto k^{3}$ and $k^{0.85}$ for $k\lesssim
7\leavevmode\nobreak\ \text{m}^{-1}$ and $k\gtrsim 20\leavevmode\nobreak\
\text{m}^{-1}$ respectively hocking1987reflection . Coupled with the
aforementioned amplitude modulation, this wavenumber dependence suggests that
higher frequency waves will have sufficient energy to induce attraction at
farther hull-wall distances.

Figure 5: Hocking radiation hypothesis and model prediction. (A) When
$d^{\perp}_{0}\gg d^{\perp}_{T}$, the reflected waves have insufficient energy
to affect the boat. When $d^{\perp}_{0}\ll d^{\perp}_{T}$, the reflected waves
perturb the free-surface height at the boat hull, yielding a reduced-amplitude
field. This amplitude asymmetry produces a net radiation force toward the
boundary. (B) Despite the non-monotonic nature of $A(\omega)$, the predicted
net radiation force near the wall compares favorably with experimental results
from Fig. 3C. (Inset) Radiation forces on boat sides due to emitted GC waves
as predicted by Eq. (3) and Fig. 4C.

Further, GC waves radiated by a vertically oscillating body have energy given
by the following hocking1988capillary :

$E_{R}=\frac{\pi}{2}\left(1+\frac{3\gamma k^{2}}{\rho g}\right)A^{2}.$ (2)

Considering the boat as two back-to-back, semi-circular wave emitters, this
expression implies the following radiation force incident upon one side:

$|F_{R}|=\frac{E_{R}k}{4\pi}=\left(\frac{k}{8}+\frac{3\gamma k^{3}}{8\rho
g}\right)A^{2}.$ (3)

The factor of $4\pi$ accounts for projecting the wave momentum normal to the
semi-circular boundary (see, e.g., Ref. lee2017fluctuation  for a more
detailed derivation). For our boat which has non-trivial $A(\omega)$, the
expected $F_{R}$ follows a power law with exponent between $3$ and $4$ (Fig.
5B Inset). We reiterate that the amplitude measurements were taken within the
attractive regime, and so we shift the origin of our power law to the observed
threshold frequency for attraction. The difference between $F_{R}$ on either
side of the boat yields a predicted $F_{W}$; this prediction matches well with
experimental results without any fitting parameters (Fig. 5B).

We summarize our postulated model of the boat’s boundary-driven locomotion in
four regimes. In all cases, when the boat first emits waves, the field is
symmetric, leading to a net-zero radiation force on the boat. When
$d^{\perp}_{0}\gg d^{\perp}_{T}$, the waves reflected off the boundary return
to the boat with negligible energy compared to emission. Consequently, the
boat experiences a force negligibly close to zero. When $d^{\perp}_{0}\ll
d^{\perp}_{T}$, the reflected waves hinder wave generation between the boat
and wall, leading to an observed amplitude reduction. Meanwhile, waves on the
far side remain unchanged; this broken symmetry yields a net force appearing
as a boat-wall attraction. When approaching $d^{\perp}_{T}$ from above,
reflected waves have insufficient amplitude to affect wave generation but
still carry non-negligible momentum. Symmetry is again broken and the boat
experiences a slight repulsive force. Since the energy of a reflected GC wave
increases with $k$, $d^{\perp}_{T}$ also increases with $k$ (and consequently
$\omega$). Should the original choice of $d^{\perp}_{0}$ be retained while
increasing $\omega$, the reflected waves will then have sufficient energy to
affect wave generation, causing the same result as when $d^{\perp}_{0}\ll
d^{\perp}_{T}$. Given the importance of the seldom-studied reflection and
generation properties of GC waves to the boat’s locomotion, we refer to these
asymmetric wave fields as “Hocking radiation fields.”

We note that our model can only explain the boat’s steady-state position using
wave amplitudes measured in that state. Additionally, Hocking’s theories on GC
waves require both the emitter and reflecting boundary to be stationary on
average. A much harder problem then is computing the system dynamics as the
wave field updates; how would one compute the position versus time of the
attracting boat in a dynamic environment? Indeed, we find the boat exhibits
complex attractive modes like “towing” in response to a moving boundary (see
SI, Movie S5). These dynamical experiments will help characterize transient
locomotive states owed to Hocking radiation in stationary and active
environments.

Conclusion – In this Letter, we revealed how a fluid surface swimmer can use
“Hocking radiation” to locomote without the need for a traditional propulsion
mechanism and made the first direct measurement of the corresponding force. In
doing so, we add to the growing list of transport phenomena that employ
surface wave fields both for propulsion and nonlocal interaction with fellow
substrate occupants lee2017fluctuation ; rhee2022surferbot ; ko2022small ;
harris2018interaction ; boersma1996maritime ; denardo2009water ;
saenz2021emergent ; ho2021capillary ; yuan2021wave ; wilcox1979sex ;
roh2019honeybees ; harris2015pilot ; benham2022gunwale . By symmetrically
generating waves near a boundary, our boat takes advantage of the reflection
dynamics unique to GC waves to self-propel with frequency- and distance-
dependent locomotive modes. Though the necessity of a nearby boundary
condition seemingly imposes restrictions on the relevance or usefulness of
Hocking radiation, manipulation of oscillation frequency and profile in
response to transient conditions may prove valuable in overcoming such
limitations.

###### Acknowledgements.

We thank Enes Aydin for designing and constructing the tank apparatus and for
helping during the boat design phase. We thank Ryan Hirsh for assisting with
pendulum experiments. We thank Paul Umbanhowar for helpful comments and
discussion. This work was funded by the Army Research Office GR00008673
(D.I.G.).

##  I Supplementary Information

Details of boat vibration – The robot boat vibrates in response to the
oscillation of an internally mounted eccentric motor (Vybronics Inc.
Cylindrical Vibration Motor VJQ24-35K270B). High-speed measurements of the
motor in motion reveal a consistent frequency $\omega_{0}$ in response to
power input $P_{0}$. We observe two response modes:
$P_{0}\propto\omega_{0}^{2}$ and $\omega_{0}^{3}$ for $\omega_{0}\lesssim$ and
$\gtrsim 30$ Hz respectively (Fig. S1A). The boat resonates at the crossover
between modes, which appears experimentally as both the maximum in
$A$-$\omega$ space (see main text) and the only significant deviation from
linearity in $\omega$-$\omega_{0}$ space (Fig. S1B).

Figure S1: Eccentric motor oscillation drives boat vibration and consequent
wave generation. (A) Motor frequency response falls into two distinct modes
with a resonance emerging at the crossover. (B) The boat’s vibration is damped
significantly due to coupling effects between the motor and boat and the fluid
surface and hull. (C) Wave amplitude measurements yielded a dispersion
relation comparable to established theory on GC waves lamb1924hydrodynamics .

We cast the boat’s vibratory response in 1D using the rotational analog of
Newton’s Second Law. By suspending the boat in midair on a string, we
eliminate the need to model the complex feedback mechanisms owed to surface
wave generation. Consequently, the relevant torques on the boat hull are
produced by gravity, air drag, and the eccentric motor, which we model as a
rotating unbalance palm2014system :

$\begin{split}I\ddot{\theta}=m_{0}\epsilon_{0}\omega_{0}^{2}R_{0}\sin{\left(\omega_{0}t-\theta\right)}-mgR_{B}\sin{\theta}\\\
-\frac{\pi}{10}\rho
c_{D}R_{B}^{5}\operatorname{sgn}\big{(}\dot{\theta}\big{)}\dot{\theta}^{2},\end{split}$
(4)

where $I$, $m$, $R_{B}$ and $c_{D}$ are the boat’s moment of inertia, mass,
radius, and drag coefficient respectively; $m_{0}$ and $\epsilon_{0}$ are the
rotating unbalance’s mass and radius respectively; $R_{0}$ is the distance
between the motor shaft and boat hull; and $\rho$ is the density of air. We
simulate Eq. (4) with an ordinary differential equation solver in MATLAB and
find an expected response frequency $\omega=0.99\omega_{0}$ (Fig. S1B).
However, physical measurement of the boat’s in-air vibration reveals a reduced
response driven at 69% the motor’s frequency. Placed in water, the boat’s
vibration drops further to 60%, with the surface wave frequency nearby at 57%.
We attribute these discrepancies to two sources of damping, namely the motor’s
non-idealized mounting to the boat and the coupling between the fluid surface
and the hull.

Figure S2: Boat vibrates primarily along the roll axis with small amplitude.
(A) Archetypal boat vibration at 38.1 Hz. (B) Amplitudes of oscillation along
roll, pitch, and vertical displacement axes. For nearly all accessible
frequencies, the primary boat oscillation occurs along the roll axis.

To better understand the boat’s vibratory response in 3D, we tracked the
oscillation with multiple high-speed cameras (OptiTrack) at 360 FPS (Fig. S2).
Unlike many established systems that employ periodic heaving (vertical)
motions to generate surface waves tatsuno1969transfiguration ; hulme1982wave ;
hocking1988capillary ; taneda1991visual ; benham2022gunwale , our boat
undergoes minimal vertical displacement. Instead, the eccentric motor induces
oscillations primarily along the fore-aft (roll) axis. Still, the boat’s
overall vibrational motion is miniscule, with a maximum roll amplitude
$\phi_{R}=0.20\degree\pm 0.02\degree$ corresponding to a vertical amplitude of
$0.21\pm 0.02$ mm.

Surface currents – As a further check on our hypothesis that the boat’s
attraction toward and repulsion from boundaries is the result of surface
waves, we investigated the possible existence of surface currents induced by
the boat’s wave generation. After fixing the boat’s lateral position such that
vibration was unimpeded, we seeded the fluid surface with lycopodium powder
(CAS number: 8023-70-9) for use with open-source Particle Image Velocimetry
(PIV) software in MATLAB thielicke2021particle . Results are shown in Movie
S4.

For $\omega<19$ Hz, no surface currents are produced, and seed particles trace
circular paths in the vertical plane as they bob over the waves (Fig. S3A).
For $19$ Hz $\leq\omega\leq 29$ Hz, vortices emerge as seed particles are
drawn in at the fore and aft, circulate along the boat perimeter, and eject in
jets at the left and right sides with maximum velocity $v\approx 8$ mm/s (Fig.
S3B). We rationalize the ingress and egress positions as the locations with
the weakest and strongest vibratory motion respectively, a result of the
orientation of the eccentric motor driving the oscillation. For $\omega>29$
Hz, a few smaller vortices with maximum velocity $v\approx 1$ mm/s emerge
inconsistently around the boat perimeter (Fig. S3C). All three regimes persist
when the boat is brought near a boundary.

The frequency dependence of these distinct surface current modes does not
correlate with that of Hocking radiation as described in the main text.
Furthermore, when orienting the boat such that the primary surface jets expel
toward a nearby boundary, the boat still experiences an attractive force where
mechanical intuition suggests a repulsion should emerge. For these reasons, we
rule out surface currents as a probable cause for Hocking radiation and
reaffirm our surface wave hypothesis.

Figure S3: Boat vibration generates minimal surface currents. (A-C) Surface
currents induced by wave generation at 6.3, 19.6, and 33.5 Hz respectively.
Currents shown are representative of behaviors within the three distinct
regimes: $\omega\in(0,19),[19,29],$ and $(29,42]$ Hz.

Synthetic Schlieren imaging SchlierenInCentipedesNote  – Before starting the
experiments, we captured a reference image of the background pattern
(checkerboard) as seen through a still free surface with the high speed
camera. During experiments, surface waves appeared as a distortion field
$\overrightharp{u}$ applied to the checkerboard. We compared the spatial
Fourier transform of the distorted checkerboard to that of the reference image
to find how the carrier peaks were modulated. When the free surface curvature
had focal length greater than the distance to the background pattern (i.e.,
the invertibility condition is met moisy2009synthetic ) we filtered the
modulated signal to extract $\overrightharp{u}(t,\overrightharp{r})$, which is
proportional to the gradient of the free surface height. Moisy and colleagues
quantify this invertibility condition as follows:

$h_{p}<h_{p,c}=\frac{\lambda^{2}}{4\pi^{2}\alpha\eta_{0}},$ (5)

where $h_{p}$ is the effective surface-pattern distance, $h_{p,c}$ is the free
surface focal length, $\lambda$ is the wavelength, $\alpha$ is the ratio of
indices of refraction given by $1-n_{\text{air}}/n_{\text{fluid}}$, and
$\eta_{0}$ is the wave amplitude moisy2009synthetic . Further, we adapted the
open-source code in Ref. wildeman2018real  for use with our apparatus,
incorporating a scale factor to account for additional interfaces between the
background pattern and the fluid free surface moisy2009synthetic .

Quantitatively identifying where the invertibility condition fails requires
knowledge of wave properties that are not known a priori and cannot be
reliably obtained from the reconstruction itself. However, we note that failed
reconstruction surface height data typically is highly discontinuous, both
with itself and with successfully reconstructed surface heights. We used this
characteristic to estimate regions where the reconstruction failed per video
frame with an autocorrelation method described by the following steps:

  1. 1.

Perform a 2D spatially-moving variance with square kernel given by the 8-way
nearest pixel neighbors.

  2. 2.

Compare the moving variance to a threshold value. We obtained our threshold
through trial-and-error but postulate that it is related to the effective
distance between the free surface and background pattern.

  3. 3.

Convert any pixels for which the variance exceeds the threshold to a mask.

  4. 4.

Perform minor cleanup on the mask using morphological operations. The result
is an estimate of all failed surface reconstructions in the frame.

Probing response to moving boundaries – Though the constrained pendulum system
enabled measurement of the Hocking radiation force, it restricted
phenomenological exploration to a firmly asymmetric-field regime. To probe the
existence of transition dynamics between the asymmetric and symmetric regimes
of the Hocking radiation phenomenon, we measured the response of a free-
floating boat with constant $\omega$ to a wall retreating with constant speed
$v_{\text{wall}}$ in 1D (Fig. S4, Movie S5). When initiated with a boat-wall
attraction, we posit the existence of a bifurcation in $v_{\text{wall}}$ at
which the wall would “tow” the boat with constant $d^{\perp}$.

Figure S4: Attractive Hocking radiation enables towing by moving boundary. (A)
Force diagram for towing experiments. (B) $d^{\perp}$ versus $t$ with $t_{0}$
corresponding to the wall’s initialization. For low $v_{\text{wall}}$, the
boat catches the wall within 20 s. For high $v_{\text{wall}}$, the wall
rapidly outpaces the boat. Slight variance of $d^{\perp}_{0}$ caused a
nonmonotonic trend with $v_{\text{wall}}$. (C) Lab- and (Inset) wall-frame
time series of towing experiment closest to $v_{\text{wall}}$ bifurcation.
Snapshots correspond to red points in (B).

We mounted the wall from previous experiments on a linear actuator (Firgelli®
FA-240-S-12-18) powered with constant current and pulse-width modulated
voltage. Velocimetry measurements thielicke2021particle  revealed minimal
surface currents ($v<5$ mm/s) near the wall’s center; consequently, we
performed experiments in this central region. We chose $d^{\perp}_{0}=1.18\pm
0.28$ cm and $\omega=41.9$ Hz such that the boat started firmly within the
attractive regime. Once $v_{\text{boat}}\approx
v^{\text{min}}_{\text{wall}}=2$ mm/s, we initialized the wall and recorded
$d^{\perp}(t)$ using an onboard Logitech C920 webcam (Fig. S4B).

When $v_{\text{wall}}\leq 4.9$ mm/s, acceleration induced by Hocking radiation
was sufficient for the boat to catch the retreating boundary. For large
$v_{\text{wall}}$, the wall swiftly outperformed the boat’s locomotion. We
most closely approached our expected bifurcation when $v_{\text{wall}}=7.7$
mm/s. During this trial, the wall towed the boat with $d^{\perp}<2$ cm for a
total distance of $43$ cm. Despite the stringent boundary requirement for
Hocking radiation to appear, proper choice of
$\left(\omega,v_{\text{wall}}\right)$ enabled the boat to travel over 10x
further than the corresponding $d^{\perp}_{T}(\omega)$.

## References

  * [1] Edward M Purcell.  Life at low reynolds number.  American Journal of Physics, 45(1):3–11, 1977.
  * [2] Ryan D Maladen, Yang Ding, Chen Li, and Daniel I Goldman.  Undulatory swimming in sand: subsurface locomotion of the sandfish lizard.  Science, 325(5938):314–318, 2009.
  * [3] Vasily Kantsler, Jörn Dunkel, Marco Polin, and Raymond E Goldstein.  Ciliary contact interactions dominate surface scattering of swimming eukaryotes.  Proceedings of the National Academy of Sciences, 110(4):1187–1192, 2013.
  * [4] Baxi Chong, Juntao He, Shengkai Li, Eva Erickson, Kelimar Diaz, Tianyu Wang, Daniel Soto, and Daniel I Goldman.  Self-propulsion via slipping: Frictional swimming in multilegged locomotors.  Proceedings of the National Academy of Sciences, 120(11):e2213698120, 2023.
  * [5] Kelimar Diaz, Baxi Chong, Steven Tarr, Eva Erickson, and Daniel I. Goldman.  Water surface swimming dynamics in lightweight centipedes.  arXiv preprint arXiv:2210.09570, 2022.
  * [6] Chris Roh and Morteza Gharib.  Honeybees use their wings for water surface locomotion.  Proceedings of the National Academy of Sciences, 116(49):24446–24451, 2019.
  * [7] Ki Yun Lee, Lu Wang, Jinhong Qu, and Kenn R Oldham.  Milli-scale biped vibratory water strider.  In 2019 International Conference on Manipulation, Automation and Robotics at Small Scales (MARSS), pages 1–6. IEEE, 2019.
  * [8] Eugene Rhee, Robert Hunt, Stuart J Thomson, and Daniel M Harris.  Surferbot: a wave-propelled aquatic vibrobot.  Bioinspiration & Biomimetics, 17(055001), 2022.
  * [9] Zhi-Ming Yuan, Minglu Chen, Laibing Jia, Chunyan Ji, and Atilla Incecik.  Wave-riding and wave-passing by ducklings in formation swimming.  Journal of Fluid Mechanics, 928, 2021.
  * [10] Andras Karsai, Deniz Kerimoglu, Daniel Soto, Sehoon Ha, Tingnan Zhang, and Daniel I Goldman.  Real-time remodeling of granular terrain for robot locomotion.  Advanced Intelligent Systems, page 2200119, 2022.
  * [11] Dave Gerr.  Propeller handbook.  International Marine Publishing, 1989.
  * [12] Teodor Burghelea and Victor Steinberg.  Onset of wave drag due to generation of capillary-gravity waves by a moving object as a critical phenomenon.  Physical Review Letters, 86(12):2557, 2001.
  * [13] R Stimson Wilcox.  Sex discrimination in gerris remigis: role of a surface wave signal.  Science, 206(4424):1325–1327, 1979.
  * [14] Graham P Benham, Olivier Devauchelle, Stephen W Morris, and Jerome A Neufeld.  Gunwale bobbing.  Physical Review Fluids, 7(7):074804, 2022.
  * [15] Anand U Oza, Rodolfo R Rosales, and John WM Bush.  A trajectory equation for walking droplets: hydrodynamic pilot-wave theory.  Journal of Fluid Mechanics, 737:552–570, 2013.
  * [16] Daniel Martin Harris.  The pilot-wave dynamics of walking droplets in confinement.  PhD thesis, Massachusetts Institute of Technology, 2015.
  * [17] Daniel M Harris, P-T Brun, Adam Damiano, Luiz M Faria, and John WM Bush.  The interaction of a walking droplet and a submerged pillar: From scattering to the logarithmic spiral.  Chaos: An Interdisciplinary Journal of Nonlinear Science, 28(9):096105, 2018.
  * [18] Pedro J Sáenz, Giuseppe Pucci, Sam E Turton, Alexis Goujon, Rodolfo R Rosales, Jörn Dunkel, and John WM Bush.  Emergent order in hydrodynamic spin lattices.  Nature, 596(7870):58–62, 2021.
  * [19] Sander Wildeman.  Real-time quantitative schlieren imaging by fast fourier demodulation of a checkered backdrop.  Experiments in Fluids, 59(6):1–13, 2018.
  * [20] Hendrik BG Casimir and Dirk Polder.  The influence of retardation on the london-van der waals forces.  Physical Review, 73(4):360, 1948.
  * [21] GL Klimchitskaya, U Mohideen, and VM Mostepanenko.  The casimir force between real materials: Experiment and theory.  Reviews of Modern Physics, 81(4):1827, 2009.
  * [22] LM Woods, Diego Alejandro Roberto Dalvit, Alexandre Tkatchenko, P Rodriguez-Lopez, Alejandro W Rodriguez, and R Podgornik.  Materials perspective on casimir and van der waals interactions.  Reviews of Modern Physics, 88(4):045003, 2016.
  * [23] Sipko L Boersma.  A maritime analogy of the casimir effect.  American Journal of Physics, 64(5):539–541, 1996.
  * [24] Alpha A Lee, Dominic Vella, and John S Wettlaufer.  Fluctuation spectra and force generation in nonequilibrium systems.  Proceedings of the National Academy of Sciences, 114(35):9255–9260, 2017.
  * [25] Bruce C Denardo, Joshua J Puda, and Andrés Larraza.  A water wave analog of the casimir effect.  American Journal of Physics, 77(12):1095–1101, 2009.
  * [26] C Parra-Rojas and R Soto.  Casimir effect in swimmer suspensions.  Physical Review E, 90(1):013024, 2014.
  * [27] D Ray, C Reichhardt, and CJ Olson Reichhardt.  Casimir effect in active matter systems.  Physical Review E, 90(1):013019, 2014.
  * [28] Ran Ni, Martien A Cohen Stuart, and Peter G Bolhuis.  Tunable long range forces mediated by self-propelled colloidal hard spheres.  Physical Review Letters, 114(1):018302, 2015.
  * [29] LM Hocking.  Reflection of capillary-gravity waves.  Wave Motion, 9(3):217–226, 1987.
  * [30] LM Hocking.  The damping of capillary–gravity waves at a rigid boundary.  Journal of Fluid Mechanics, 179:253–266, 1987.
  * [31] LM Hocking.  Capillary–gravity waves produced by a heaving body.  Journal of Fluid Mechanics, 186:337–349, 1988.
  * [32] Horace Lamb.  Hydrodynamics.  University Press, 1924.
  * [33] Michael S Longuet-Higgins and RW Stewart.  Radiation stresses in water waves; a physical discussion, with applications.  In Deep sea research and oceanographic abstracts, volume 11, pages 529–562. Elsevier, 1964.
  * [34] Krassimir D Danov, Rumiana Dimova, and Bernard Pouligny.  Viscous drag of a solid sphere straddling a spherical or flat surface.  Physics of Fluids, 12(11):2711–2722, 2000.
  * [35] Christian Scholz, Soudeh Jahanshahi, Anton Ldov, and Hartmut Löwen.  Inertial delay of self-propelled particles.  Nature Communications, 9(1):1–9, 2018.
  * [36] Classical Schlieren techniques rely on the precise alignment of two light-filtering masks. The replacement of a physical mask with a digitally-generated (synthetic) mask circumvents this non-trivial task (see Ref. [48]).
  * [37] Frédéric Moisy, Marc Rabaud, and Kévin Salsac.  A synthetic schlieren method for the measurement of the topography of a liquid interface.  Experiments in Fluids, 46(6):1021–1036, 2009.
  * [38] Jeffrey Aguilar and Daniel I Goldman.  Robophysical study of jumping dynamics on granular media.  Nature Physics, 12(3):278–283, 2016.
  * [39] Saebyok Bae and Yoon-Hwan Kang.  Optimal pumping in a model of a swing.  European Journal of Physics, 27(1):75, 2005.
  * [40] Hungtang Ko, Mathias Hadgu, Keyana Komilian, and David L Hu.  Small fire ant rafts are unstable.  Physical Review Fluids, 7(9):090501, 2022.
  * [41] Ian Ho, Giuseppe Pucci, Anand U Oza, and Daniel M Harris.  Capillary surfers: wave-driven particles at a fluid interface.  arXiv preprint arXiv:2102.11694, 2021.
  * [42] William J. Palm III.  System dynamics.  McGraw-Hill New York, third edition, 2014.
  * [43] M Tatsuno, S Inoue, and J Okabe.  Transfiguration of surface waves.  Rep. Res. Inst. Appl. Mech. Kyushu University, 17:195–215, 1969\.
  * [44] A Hulme.  The wave forces acting on a floating hemisphere undergoing forced periodic oscillations.  Journal of Fluid Mechanics, 121:443–463, 1982.
  * [45] Sadatoshi Taneda.  Visual observations of the flow around a half-submerged oscillating sphere.  Journal of Fluid Mechanics, 227:193–209, 1991.
  * [46] William Thielicke and René Sonntag.  Particle image velocimetry for matlab: Accuracy and enhanced algorithms in pivlab.  Journal of Open Research Software, 9(1), 2021.
  * [47] Wildeman’s open-source code [19] and author SWT’s modifications derived from Refs. [19, 37] as described in this section were also used in the analysis of Ref. [5]. The descriptive text is reproduced from the SI accompanying Ref. [5] with minor edits.
  * [48] Bruce R Sutherland, Stuart B Dalziel, Graham O Hughes, and PF Linden.  Visualization and measurement of internal waves by ‘synthetic schlieren’. part 1. vertically oscillating cylinder.  Journal of Fluid Mechanics, 390:93–126, 1999.


# A Monte Carlo Method for Evaluating Empirical Gyrochronology Models and its
Application to Wide Binary Benchmarks

Tomomi Otani Department of Physical Sciences, Embry-Riddle Aeronautical
University, 1 Aerospace Blvd, Daytona Beach, FL 32114, United States  Ted von
Hippel Department of Physical Sciences, Embry-Riddle Aeronautical University,
1 Aerospace Blvd, Daytona Beach, FL 32114, United States  Derek Buzasi
Department of Chemistry and Physics, Florida Gulf Coast University, 10501 FGCU
Blvd S, Fort Myers, FL 33965, United States  T. D. Oswalt Department of
Physical Sciences, Embry-Riddle Aeronautical University, 1 Aerospace Blvd,
Daytona Beach, FL 32114, United States  Alexander Stone-Martinez  Department
of Physical Sciences, Embry-Riddle Aeronautical University, 1 Aerospace Blvd,
Daytona Beach, FL 32114, United States  Department of Astronomy, New Mexico
State University, 1780 E University Ave, Las Cruces, NM 88003, United States
Patrice Majewski  Department of Physical Sciences, Embry-Riddle Aeronautical
University, 1 Aerospace Blvd, Daytona Beach, FL 32114, United States  Whiting
School & Engineering, Johns Hopkins University, 3400 North Charles Street,
Baltimore, MD 21218, United States

(Accepted Mar 21, 2022)

###### Abstract

Accurate stellar ages are essential for our understanding of the star
formation history of the Milky Way, Galactic chemical evolution, and to
constrain exoplanet formation models. Gyrochronology, a relationship between
stellar rotation and age, appears to offer a reliable age indicator for main
sequence (MS) stars over the mass range of approximately 0.6 to 1.3
$M_{\sun}$. Those stars lose their angular momentum due to magnetic braking
and as a result their rotation speeds decreases with age. Although current
gyrochronology relations are fairly well tested for young MS stars with masses
greater than 1 $M_{\sun}$, primarily in young open clusters, insufficient
tests exists for older and lower mass MS stars. Binary stars offer the
potential to expand and fill in the range of ages and metallicity over which
gyrochronology can be empirically tested. In this paper, we demonstrate a
Monte Carlo approach to evaluate gyrochronology models using binary stars. As
examples, we used five previously published wide binary pairs. We also
demonstrate a Monte Carlo approach to assess the precision and accuracy of
ages derived from each gyrochronology model. For the traditional Skumanich
models, the age uncertainties are $\sigma_{age}$/$age$ = 15-20% for stars with
$B-V$ = 0.65, and $\sigma_{age}$/$age$ = 5-10% for stars with $B-V$ = 1.5 and
rotation period P $\leq$ 20 days.

††journal: ApJ

##  1 Introduction

Ages are among the most difficult stellar properties to determine, yet they
provide key constraints to problems ranging from the formation and
habitability of exoplanets to the Galaxy’s chemical evolution and star
formation history. For decades it has been known that main sequence (MS) stars
of spectral type F through early M have outer convection zones with strong
magnetic fields and that they lose angular momentum due to magnetic braking
(Schatzman, 1962; Weber & Davis, 1967; Skumanich, 1972; Kawaler, 1988;
Pallavicini et al., 1981; Soderblom, 1985; Pinsonneault et al., 1989). This
magnetic braking among low mass MS stars was observed by Skumanich (1972), who
found that the rotation period was proportional to the square root of stellar
age. This is now referred to as the Skumanich law and this phenomenon of a
declining rotation rate with time is referred to as the slow rotator sequence.
Subsequently, both empirical/semi-empirical (Barnes, 2001, 2003, 2007, 2010;
Barnes & Kim, 2010; Brown, 2014; Angus et al., 2015; Gondoin, 2017; Angus et
al., 2019) and theoretical (Kawaler, 1988, 1989; Pinsonneault et al., 1989;
van Saders & Pinsonneault, 2013; Matt, et al., 2015; van Saders et al., 2016;
Spada & Lanzafame, 2020) frameworks for measuring stellar ages from rotation
periods, now called gyrochronology, have been developed. Numerous researchers
have made the case that gyrochronology could be developed into one of the most
precise stellar chronometric techniques for $\sim$0.6–1.0 M⊙ field stars
(Barnes, 2003; Soderblom, 2010; Epstein & Pinsonneault, 2014; van Saders et
al., 2016).

When a star arrives on the zero age main sequence, its rotational period is
typically between one and ten days (Irwin & Bouvier, 2009; Gallet & Bouvier,
2013). Rotation periods begin to follow the Skumanich law when stars are 500
to 700 Myr old, and by the age of 1 Gyr, all stars with masses between 0.6 and
1.3 $M_{\odot}$ will join this sequence. Although open clusters confirm that
the gyrochronology paradigm is valid for the Sun and lower mass stars, there
is ongoing concern about how age precision is affected by such factors as the
initial range of rotation rates (Barnes, 2010; Matt et al., 2012), the choice
of spin-down model (Aigrain et al., 2015), changing magnetic morphology
(Buzasi, 1997; Garraffo et al., 2018), and how mass loss rates may change for
stars older than the Sun (van Saders et al., 2016; Metcalfe & Egeland, 2019).
In addition, evidence is accumulating that the rate of rotational spin-down
may temporarily stall and does not follow the traditional Skumanich law
(Meibom et al., 2011b, 2015; McQuillan, Mazeh, & Aigrain, 2014; Rebull, 2017;
Curtis et al., 2019). Curtis et al. (2019) found that both early and late K
stellar rotation rates were the same in the open cluster NGC 6811 ($\sim$1
Gyr). Curtis et al. (2019) and Douglas et al. (2019) compared the open
clusters NGC 6811 and Praesepe and note that despite their age differences
both low mass star sequences below 0.8 $M_{\odot}$ seem to overlap, and
furthermore that low mass stars in NGC 6811 may not have spun down for the
past 300 Myr. Angus et al. (2020) and Spada & Lanzafame (2020) suggests that
the angular momentum transport from the core to the envelope of stars causes
this extended period over which rotation rates do not decrease. Additionally,
evidence is accumulating that field stars older than $\sim$1 Gyr rotate more
rapidly than predicted by the Skumanich law (Angus et al., 2015; van Saders et
al., 2016; van Saders, Pinsonneault, & Barbieri, 2018; Metcalfe & Egeland,
2019). Curtis et al. (2020) studied the dispersed star cluster Ruprecht 147
($\sim$2.7 Gyr) and found that its gyrochrone is flat compared to other
clusters with younger ages. Also, Agüeros et al. (2018) suggests that the slow
rotator sequence of the open cluster NGC 752 (1.3 Gyr) departs from the
classical Skumanich law. Both van Saders et al. (2016) and van Saders,
Pinsonneault, & Barbieri (2018) suggest that magnetic braking would stop when
stars become sufficiently old that their Rossby number ($R_{O}$), the ratio of
rotation period to convective turnover time, reaches approximately two.
Therefore, while the traditional Skumanich law may work well for younger F and
G stars, it may not apply to cooler and older stars. For these reasons, the
field of gyrochronology requires both a wider range of empirical tests and
further theoretical development.

Empirical tests based on ages derived from wide non-interacting binary stars
offer several advantages. First, components of each pair are coeval and should
exhibit the same rotation age. Many thousands of nearby pairs are potentially
useful for testing gyrochronology. Among the 1.3 million wide binaries found
in the $Gaia$ EDR3 by El-Badry, Rix, & Heintz (2021), about 22,000 are within
100 pc. The sheer number of stars in open clusters offers substantial
advantages, yet due to their distances, low-mass members of many clusters such
as NGC 6811, NGC 752, NGC 6819, Ruprecht 147, and M67, are too faint to
readily observe. Additionally, binary pairs span a well-sampled range of ages,
extending beyond the time most open clusters have lost most low mass members
and they span a wide range of metallicities. Finally, unlike asteroseismology,
which with current technology is limited to bright stars with a relatively
narrow range of masses, wide pairs contain stars that span the entire range of
masses over which gyrochronology is believed to apply.

Gyrochronology requires measured rotation periods. The $Kepler$ missions, with
their continuous monitoring of selected star fields, have obtained data useful
for deriving the rotation periods for thousands of stars. Unfortunately, the
vast majority of $Kepler$ and $K2$ stars do not have $B-V$ color measurements,
the most common proxy for mass upon which most current gyrochronology models
are based. For about 100 wide pairs in the original $Kepler$ field, Janes
(2017) constructed an empirical $B-V$ vs. $g-K$ relation. A more recent search
for wide pairs in the Kepler field was conducted by Godoy-Rivera (2018). This
study used $Gaia$ DR2 astrometry, parallaxes, radial velocities and
metallicities to identify candidate pairs.

In this paper we develop a Monte Carlo strategy for propagating uncertainties
in the observed stellar parameters and theoretical model coefficients in order
to quantify the resulting age uncertainties, then propose a model evaluation
method using wide binaries, under the assumption that binary components are
coeval. As an application of this strategy, we chose six gyrochronology models
that are based on rotation periods and colors (Barnes, 2007; Mamajek &
Hillenbrand, 2008; Meibom et al., 2011a; Angus et al., 2015, 2019; Spada &
Lanzafame, 2020). (Hereafter we refer to these models as BA07, MH08, M11, A15,
A19, and SL20, respectively).

Where age inconsistencies between binary components arise, they may be caused
by a variety of observational or theoretical reasons: [itemize] & one or both
members of the pair suffer from exceptionally large errors in color or
rotation period, one or both members of the pair lie outside the region of
parameter space in which gyrochronology is applicable, undetected tertiary
components modified the evolutionary history, mass, or rotation period of the
observed star, blending of the target star with another star along the line of
sight, the presence of multiple stellar spots or groups yielding an incorrect
rotation period, spots measured at different latitudes and/or differential
rotation, pulsation periods masquerading as rotation periods, and problems
with the gyrochronology models such as uncertainties, fitting parameters
inaccuracies, and missing physics in gyrochronology model. A precise
determination of both consistencies and inconsistencies among stars with
gyrochronological ages will help us investigate these issues in order to
improve both the observational techniques and the models.

The models that were used in this manuscript is shown in Section 2. The
uncertainty sensitivity of the five empirical models (B07, MH08, M11, A15, and
A19) are discussed in Section 3. The model evaluation method is discussed in
Section 4. Application of the method using the six gyrochronology models and
an evaluation of the age precision of those models are discussed in Section 5.
The conclusions drawn from these numerical experiments are described in
Section 6.

##  2 Gyrochronology Models

Barnes (2003) first proposed an empirical relationship between the rotation
period, color, and age of the form:

$P=\sqrt{t}\>f(B-V)$ (1)

where

$f(B-V)=\sqrt{B-V-0.5}-0.15(B-V-0.5)$ (2)

Here, $P$ is the rotation period in days, $t$ is the age of the star in Myr
and, $B-V$ is the color index. Later, Barnes (2007) updated the relationship
as:

$P=A^{n}\times a~{}(B-V-c)^{b}$ (3)

where $A$ is the age of the star in Myr, and $a$, $b$, $c$, and $n$ are
dimensionless free parameters. The four empirical models that are discussed in
this section (B07, MH08, M11, and A15) use this simple relationship. The
parameters used by each of these empirical models are shown in Table 1. To
estimate the age of a star, $P$ and $B-V$ must be obtained from observations.
Nearly continuous photometry from space telescopes ($Kepler$, $K2$, and
$TESS$) have provided rotation periods for thousands of stars. The majority of
$B-V$ color indices for these stars have come from ground-based photometry.

Table 1: Coefficients for four models. coefficients | B07 | MH08 | M11 | A15
---|---|---|---|---
a | 0.7725 | 0.407 | 0.700 | 0.40
$+$ $\sigma_{a}$ | 0.011 | 0.021 | 0.013 | 0.30
$-$ $\sigma_{a}$ | 0.011 | 0.021 | 0.013 | 0.05
b | 0.601 | 0.325 | 0.553 | 0.31
$+$ $\sigma_{b}$ | 0.024 | 0.024 | 0.052 | 0.05
$-$ $\sigma_{b}$ | 0.024 | 0.024 | 0.052 | 0.02
c | 0.4 | 0.495 | 0.472 | 0.45
$+$ $\sigma_{c}$ | N/A | 0.01 | 0.027 | N/A
$-$ $\sigma_{c}$ | N/A | 0.01 | 0.027 | N/A
n | 0.5189 | 0.566 | 0.52 | 0.55
$+$ $\sigma_{n}$ | 0.007 | 0.008 | N/A | 0.02
$-$ $\sigma_{n}$ | 0.007 | 0.008 | N/A | 0.09
Table 2: Model coefficients for A19. coefficients |  A19
---|---
$c_{A}$ | 0.65
$\sigma_{c_{A}}$ | 0.05
$c_{0}$ | -4.7
$\sigma_{c_{0}}$ | 0.5
$c_{1}$ | 0.72
$\sigma_{c_{1}}$ | 0.05
$c_{2}$ | -4.9
$\sigma_{c_{2}}$ | 0.2
$c_{3}$ | 29
$\sigma_{c_{3}}$ | 2
$c_{4}$ | -38
$\sigma_{c_{4}}$ | 4

Solving Equation 3 for stellar age, adding a subscript to age, $A_{g}$, to
acknowledge that this is not necessarily the true age of the star, but rather
the age provided via a gyrochronology relation, we have

$D(A_{g})=\left\\{\frac{G(P)}{D(a)\times\left[G(B-V)-D(c)\right]^{D(b)}}\right\\}^{\frac{1}{D(n)}},$
(4)

where $G(B-V)$ and $G(P)$ indicate that these observable parameters have
uncertainties characterized by Gaussian distributions using the observed
values and one-sigma uncertainties and $D(a,b,c,n)$ indicates that the model
parameter uncertainty distributions may not be Gaussian.

The Skumanich models do not always provide reliable results, particularly for
low mass main sequence stars (Meibom et al., 2011b, 2015; McQuillan, Mazeh, &
Aigrain, 2014; García et al., 2014; Douglas et al., 2017; Rebull, 2017; Curtis
et al., 2019). This may be due to angular momentum transport from stellar
interiors to their surfaces (Angus et al., 2020; Spada & Lanzafame, 2020).
Although stellar surfaces continuously lose angular momentum via magnetic wind
braking, this additional angular momentum from the interior compensates for
that loss, and it thus $appears$ that these star do not lose angular momentum
during this age range. In this paper, we use two models that include this
angular momentum transport (A19 and SL20).

The A19 model improves the age precision for F, G, and K dwarfs, including
stars that do not always follow the Skumanich-style slow rotator sequence. A19
also includes subgiant stars that were not included in previous gyrochronology
models. Because this model uses $Gaia$ $G_{BP}-G_{RP}$ colors, which are
available for more than a billion stars (Gaia Collaboration, 2018), it is
easier to find accurate and consistent color information than for those models
that rely on $B-V$ photometry. The A19 model was derived by combining both
isochrones and empirical gyrochronology models using data from Praesepe
($\sim$700 Myr) and the Sun. In order to incorporate A19 in a consistent
manner with the other models, we only used the gyrochronology component of the
full A19 model so only $B-V$ and P are required to estimate ages. A19 warns
that the gyrochronology models only works for late F, G, K and early M dwarfs
(0.56 $<$ $G_{BP}-G_{RP}$ $<$ 2.7) with Rossby numbers less than 2. The
relationship between rotation periods and age of the stars is:

$log_{10}(P)=c_{A}log_{10}(t)+\sum_{0}^{4}c_{n}[log_{10}(G_{BP}-G_{RP})]^{n}$
(5)

where P is rotation period in days, t is age of the star in years,
$G_{BP}$-$G_{RP}$ is the $Gaia$ color, and each coefficient is shown in Table
2. Ages from their model can be obtained from an open-source $Python$ package,
$stargate$, which we used (see details in A19).

SL20, updating their previous work (Lanzafame & Spada, 2015), construct
gyrochrones using the 0.4-1.3 $M_{\odot}$ star data from Praesepe ($\sim$700
Myr) and NGC6811 ($\sim$1 Gyr). Their model includes the reduced spin-down
observed in NGC 6811 due to the redistribution of angular momentum from the
stellar interior. Therefore it is built not only with the magnetic wind
braking law, but also with the semi-empirical mass dependence of the
rotational coupling time scale (Lanzafame & Spada, 2015; Somers &
Pinsonneault, 2016). This model agrees with the stalled surface spin-down
which was observed among low mass stars in Praesepe and NGC6811.

##  3 Uncertainty Sensitivity for Individual Stars

To propagate errors, we employ a Monte Carlo approach, allowing us to examine
the age precision achievable with these empirical models throughout the color-
period diagram. In this approach, we drew 500 values of color ($B-V$) and
rotation period ($P$) at grid points covering $B-V$ = 0.6-2.3 mag (or $G_{BP}$
\- $G_{RP}$ = 0.5-2.7) and $P$ = 0-60 days. For each color ($B-V$ or
$G_{BP}-G_{RP}$) and rotation period ($P$) combination, we obtained 10000
random normal distribution samples of color ($B-V$), rotation period ($P$),
and gyrochronology model coefficients using their mean values and the standard
deviations. For A15, two Gaussian distributions were created using both the
positive and negative uncertainties, to create this non-Gaussian distribution.
We did not include the effects of correlations between the parameters a and n
for A15. We also did not include the A19 and SL20 models in this uncertainty
analysis because these models did not include their uncertainties. We
calculated ages for each combination using the method discussed in Section 2,
then made distribution of the ages using the 10000 simulated data points
propagated through the appropriate model. From these distributions, we
obtained median ages ($age$) and standard deviations ($\sigma_{+}$ and
$-\sigma_{-}$). The mean of the absolute values of $\sigma_{+}$ and
$-\sigma_{-}$ are used to estimate the age uncertainties
$\frac{\sigma_{age}}{age}$.

Figure 1 presents the ranges of expected age precision in the form of this
normalized age uncertainty, $\frac{\sigma_{age}}{age}$, for the four
gyrochronology models of B07, MH08, M11, and A15. We emphasize that these are
age precision values, not age accuracy values, because they do not incorporate
any systematic errors in the models. The fractional rotation period
uncertainties, $\frac{\sigma_{P}}{P}$, for the primary and secondary were both
set at 2% of the observed period, which is the smallest uncertainties among
the wide binary and triple star examples that B07 published (36 Oph A, B, and
C). The $B-V$ uncertainties of 0.01, which B07 used in his analysis, was also
incorporated. For uncertainties of $G_{BP}$ \- $G_{RP}$, we set
$\sigma_{G_{BP}-G_{RP}}=0.033$ based on the $G$ mag of those targets (Evans et
al., 2018). For the coefficient uncertainties, the published values listed in
Table 1 were used. The uncertainties in model coefficients, observed rotation
period, and colors were propagated via this Monte Carlo approach to derive
uncertainties in the estimated ages.

Figure 1: Age uncertainty, $\frac{\sigma_{age}}{age}$, as a function of
stellar color, period, and their uncertainties. Top: The age uncertainty for
B07, MH08, and M11 models. Bottom: The age uncertainty for A15 and A19 models.
For A19, only the gyrochronology component of the model was used. Note that
the scale of A15 and A19 models (bottom) are different from the scale of the
other models (top). For the model uncertainties, we used the published values
in Table 1 and 2. For the observational parameter uncertainties we used
$\frac{\sigma_{P}}{P}=2\%$ and $\sigma_{B-V}=0.01$. We set
$\sigma_{G_{BP}-G_{RP}}=0.033$ based on the $G$ mag of those targets (Evans et
al., 2018).

Figure 1 shows how age uncertainties depend on the uncertainties in model
parameters, color, and rotation period. For B07, MH08, and M11 models, stars
with $\sigma_{P}$/$P$ = 0.10, the resulting age uncertainties are about 5-20%
for stars with $B-V\geq 1$ and above 20% for stars with $B-V\leq 1$. B07 used
an error propagation method to calculate the uncertainties of his model and
found that they are 20% for stars with $B-V$ = 0.5, 15% for stars with $B-V$ =
0.65, 13% for the stars with $B-V$ = 1.0 or 1.5, and always greater than 11%.
Our results are similar to his, however our results indicate slightly larger
age uncertainties (15-20%) for stars with $B-V$ = 0.65 and slightly smaller
uncertainties (5-10%) for stars with $B-V$ = 1.5 and $P$ $\leq$ 20 days.
Although these age uncertainties may be too large as the basis for some
research projects, e.g. refining the formation history of the Milky Way, they
may be sufficient for other projects such as interpreting the radiation
environment of an exoplanet. Using these models with high quality rotation
rates ($\sigma_{P}$/P = 0.05) is necessary within these models to keep age
uncertainties under 20%.

The coefficient uncertainties of the A15 and A19 models are significantly
larger than the coefficient uncertainties of the other two models, and these
models therefore yield larger age uncertainties. This does not mean that the
A15 and A19 models are of lower quality. Rather, they might simply be
incorporating more realistic uncertainties.111One of the reasons that A15
model obtained large coefficient uncertainties is that the model included old
asteroseismic stars that have stopped spinning down (van Saders et al., 2016).
They suggest that these large uncertainties only should apply to stars with
large Rossby numbers, $Ro$ = $P$/$\tau_{cz}$. We have chosen these five models
because their analytical form allows us to develop and demonstrate our Monte
Carlo approach to testing gyrochronology models.

##  4 Model Evaluation Method using Wide Binaries

As outlined in Section 1, we can use wide binaries to evaluate gyrochronology
models. We assume the two components are coeval (Kraus & Hillenbrand, 2009),
ratio the two gyrochronological ages, and expect this ratio to be close to
one. In this section, we propose a method to assess four Skumanich-based (B07,
MH08, M11, A15) and two non Skumanich-based (A19 and SL20) models.

For Skumanich-based models, assuming Equation 1 applies, the rotation periods
of the primary and secondary components of a binary system are

$P_{p}=(A_{g,p})^{n}\times a[(B-V)_{p}-c]^{b}$ (6)

and

$P_{s}=(A_{g,s})^{n}\times a[(B-V)_{s}-c]^{b}$ (7)

where a, b, c, and n are the parameters listed in Table 1, $P_{p}$ and $P_{s}$
are the primary and secondary rotation periods, and $[B-V]_{p}$ and
$[B-V]_{s}$ are the primary and secondary $B-V$ colors.

Although the model parameters (a, b, c, n) have uncertainties, for any given
value of those parameters the ages of the two stars must be the same (e.g.
$A_{p}=A_{s}$). Taking the ratio of equations 3 and 4 yields

$\frac{P_{p}}{P_{s}}=\left(\frac{A_{g,p}}{A_{g,s}}\right)^{n}\>\left[\frac{(B-V)_{p}-c}{(B-V)_{s}-c}\right]^{b}.$
(8)

Solving for the age ratio yields

$\frac{A_{g,s}}{A_{g,p}}=\left\\{\left[\frac{(B-V)_{p}-c}{(B-V)_{s}-c}\right]^{b}/\left(\frac{P_{p}}{P_{s}}\right)\right\\}^{1/n},$
(9)

which reduces the number of parameters, since the coefficient $a$ drops out of
the ratio. Therefore, we do not test the absolute value of the age parameter,
only the consistency of the ages derived for the two components. Nominally, we
expect $\frac{A_{g,s}}{A_{g,p}}$ = 1.

We incorporated an additional parameter, $\Delta$, as an age ratio
“tolerance”. This parameter represents the degree to which a gyrochronology
model is imperfect and yet still tolerably useful. Our goal is to use this age
tolerance parameter to explore the required precision of the models and
observational parameters in order to achieve useful gyrochronology ages. For
the remainder of this paper, we adopt $\Delta=0.1$, i.e. a 10% age difference
between the primary and secondary can be tolerated and still deemed
consistent. Other values of $\Delta$ could be chosen, based on the age
precision required for a particular project.

We again assumed Gaussian or other distributions in the observed quantities
and model parameters, as described above. Equation 9 is recast as

$D\left(\frac{A_{g,s}}{A_{g,p}}\right)=\left\\{\left[\frac{G(B-V)_{p}-D(c)}{G(B-V)_{s}-D(c)}\right]^{D(b)}/\left[\frac{G(P_{p})}{G(P_{s})}\right]\right\\}^{1/D(n)},$
(10)

where $G(B-V)$ and $G(P)$ are Gaussian distributions representing the observed
colors and periods, respectively, and their uncertainties; $D(b)$, $D(c)$, and
$D(n)$ are non-Gaussian distributions representing these model parameter
values and their uncertainties.

For non Skumanich-based models, we adopted a similar approach. We obtained the
following equation using $G(B-V)_{P}$, $G(B-V)_{S}$ (or $Gaia$-base colors
$G(G_{BP}-G_{RP})_{P}$, $G(G_{BP}-G_{RP})_{S}$ for A19), $G(P_{P})$, and
$G(P_{S})$;

$D\left(\frac{A_{g,s}}{A_{g,p}}\right)=\frac{f(G(B-V)_{P},G(P_{P}))}{f(G(B-V)_{S},G(P_{S}))}$
(11)

where $f(G(B-V),G(P))$ is the function to obtain ages for the A19 and SL20
models.

Figure 2 displays examples of Equation 6 & 11 via the age-ratio probability
distribution using each model for the color and period values for HD 155885
and HD 156026, two components of a triple star system listed in Barnes (2007).
In this figure the models parameters were drawn randomly from their quoted
error distributions and applied to Equation 6 to calculate
$\frac{A_{g,s}}{A_{g,p}}$. This process was repeated 10,000 times in order to
derive the age ratio distribution for this binary.

The rotation periods and colors of the stars displayed in Figure 2 are listed
in Table 3. For this binary, the age ratio probability distributions of all of
these models are inconsistent with $\frac{A_{g,s}}{A_{g,p}}$ = 1 $\pm$
$\Delta$ = 0.90 to 1.10. The estimated secondary star age is younger than the
primary age for all six models to a degree that exceeds uncertainties in the
observational data or model parameters.

Figure 2: An example of the age-comparison probability distribution from
Equation 6, using the models of B07, MH08, M11, A15, A19, SL20. HD 155885 and
HD156026 are used for this example with the published periods and $B-V$ colors
and their uncertainties (See in Table 3). The vertical bar shows the region
within 1 $\pm\Delta$, where $\Delta$ = 0.1.

To evaluate the degree to which the primary and secondary ages for any binary
are consistent within the tolerance $\Delta$, the fractions of the probability
distributions within 1 $\pm$ $\Delta$ were computed for each star, ratioed,
and converted to a percentage as follows:

${\rm
Coeval~{}probability}=100\times\left(\frac{\displaystyle\int_{1-\Delta=0.9}^{1+\Delta=1.1}\frac{A_{g,s}}{A_{g,p}}}{\displaystyle\int_{0}^{\infty}\frac{A_{g,s}}{A_{g,p}}}\right)$
(12)

This value corresponds to the likelihood that the binary pair is coeval within
the adopted tolerance (here 10%) and assumed model. Larger values from this
equation increase confidence in the consistency of the model and the veracity
of the observed binary properties. If this value is low, it is likely that the
stellar properties were not measured adequately or that there is a problem
with the gyrochronology model itself. Follow-up observations can determine
whether the former is the case. If not, adjustment to the empirical
gyrochronology model may be called for. We adopted 99.7% ($\geq$ 3 $\sigma$)
as the threshold for discriminating binary components with discrepant ages. In
the following section, we use the same formalism obtained above to examine a
related issue – the age precision possible with current and hypothetically
improved uncertainties in the empirical model constants and the observed
parameters.

##  5 Examples and Discussion

The wide binaries discussed in B07 were selected to illustrate the model
evaluation analysis described in Section 4. Among theses samples, B07 chose
pairs that had both rotation period and B-V information. He also used the
systems without rotational interaction between the components (except 36 Oph A
and B) The colors and rotation periods of each component are listed in Table
3. The dereddened $G_{BP}-G_{RP}$ was calculated using bolometric corrections
to $G_{BP}$ and $G_{RP}$ (Casagrande & VandenBerg, 2018). The rotation periods
and colors of these wide binary stars are shown in the color-period diagrams
in Figure 3 along with the gyrochrone grids of the B07, MH08, M11, A15, A19,
and SL20. As examples, we plot three cases (blue pairs in Figure 3) where the
ages of the two components appear to be consistent using all three
gyrochronology models and two cases (red pairs in Figure 3) in which the ages
of the two components appear to be inconsistent.

The lines connecting the components of age-consistent pairs are nearly
parallel to a gyrochrone while those of age-inconsistent pairs are clearly not
parallel to any gyrochrone. An eyeball assessment in this diagram provides a
helpful age consistency check, though to determine quantitative age
consistencies on must incorporate the observed and model uncertainties, the
latter of which are not apparent in this diagram.

Table 3: Wide binaries and a triple star system

used in the experiment presented in Figure set 4. The rotation period,
rotation period uncertainty, and de-reddened (B-V) are obtained from B07. 1,
2, and 3 are binary systems and 4, and 5 are in the same triple star system.
Components of 3 do not have Gaia $G_{BP}-G_{RP}$ because those apparent
magnitudes are brighter than the Gaia bright limit.  internal ID name period
$(B-V)_{0}$ $(G_{BP}-G_{RP})_{0}$ type (days) (mag) (mag) 1 $\xi$Boo A
(HD131156 A) 6.31 (0.05) 0.76 (0.01) 0.93 (0.06) G8V $\xi$Boo B (HD131156 B)
11.94 (0.22) 1.17 (0.01) 1.52 (0.06) K4V 2 61 Cyg A (HD 201091) 35.37 (1.3)
1.18 (0.01) 1.46 (0.06) K5V 61 Cyg B (HD 201092) 37.84 (1.1) 1.37 (0.01) 1.72
(0.06) K7V 3 $\alpha$Cen A (HD 128620) 28 (3.0) 0.67 (0.01) N/A G2V
$\alpha$Cen B (HD 128621) 36 (1.8) 0.87 (0.01) N/A K1V 4 36 Oph A (HD 155886)
20.69 (0.4) 0.85 (0.01) 1.06 (0.06) K1 36 Ohp B (HD 155885) 21.11 (0.4) 0.86
(0.01) 1.06 (0.06) K1 5 36 Oph B (HD 155885) 21.11 (0.4) 0.86 (0.01) 1.06
(0.06) K1 36 Oph C (HD 156026) 18.0 (0.4) 1.16 (0.01) 1.41 (0.06) K5

Figure 3: Color-period diagrams for example wide binaries (1, 2, 3) and a
triple system (4, 5) in the model grids of B07 (upper left), MH08 (upper
middle), M09 (upper right), A15 (lower left), A19 (lower middle), and SL 20
models (lower right). The age of each gyrochrone (in Myr) is listed in the
legend of the upper right panel. The crosses indicate the observed color and
rotation period of each target along with their uncertainties. Because the
components are coeval, the line connecting the two components of a binary
should lie on a gyrochrone. Blue symbols and lines indicate that both
components have consistent ages for these gyrochronology models. Red symbols
and lines indicate that both components have inconsistent ages. The A19 panel
does not include wide binary 3 because $Gaia$ photometry for these targets is
not available. The SL20 panel does not include binaries 2 or 3 because this
model is not available in the age range of those stars.

Figure set 4 displays the corresponding age ratio probability distributions of
the components of these binary and triple star systems. As listed in Table 3,
all targets are main sequence stars and all targets except $\xi$Boo A and
$\alpha$Cen A are K dwarfs. For each model, the fraction of the age ratio
distribution where the two components are consistent within the tolerance
level ($\Delta$ = $\pm$0.1), is listed in Table 4.

Table 4: Binary component consistency fractions for five pair-wise comparisons and six models. components | B07 | MH08 | M11 | A15 | A19 | SL20
---|---|---|---|---|---|---
$\xi$Boo A/B | 0.0000 | 0.0000 | 0.0579 | 0.0000 | 0.0000 | 0.0000
61 Cyg A/B | 0.4194 | 0.7389 | 0.4139 | 0.6991 | 0.3795 | N/A
$\alpha$Cen A/B | 0.2741 | 0.3392 | 0.2036 | 0.2789 | N/A | N/A
36 Oph A/B | 0.8388 | 0.8798 | 0.8309 | 0.8435 | 0.9064 | 0.9251
36 Oph B/C | 0.0000 | 0.0000 | 0.0000 | 0.0000 | 0.0000 | 0.0271

Binary no. 1 ($\xi$Boo A/B) is a wide binary consisting of components with
spectral types G8V and K4V. B07 obtained an age of 187 Myr and 265 Myr for
$\xi$Boo A and B, respectively. Our result eliminates the possibility that the
ages of the two component are consistent within the tolerance ($\Delta$ =
$\pm$0.1) for any models except M11. One of the possible reasons for these
inconsistent ages is that the stars are too young for the Skumanich law given
subsequent work and the measurement that $\xi$Boo A (P = 6.31 days) may not
yet be on the show rotator sequence.

Binary no. 2 (61 Cyg A/B) consists of a K5V+K7V wide binary pair. The
estimated ages by B07 are 2.12 and 1.87 Gyr, respectively. The probability
that the component ages are consistent within the tolerance is larger for A15
and MH08 models than other models.

Binary no. 3 ($\alpha$Cen A/B) consists of a G2V+K1V pair. The derived ages by
B07 are 4.6 and 4.1 Gyr, respectively. The component ages appear to be
consistent, though that may only be because the age ratio distributions are so
wide, which in turn is due to the large rotation period uncertainty of
$\alpha$Cen A.

Triple star system (no. 4 and 5) consists of two chromospherically active K1
dwarfs (36 Oph A/B) and a K5 dwarf at much greater separation (36 Oph C). B07
warns that A and B could potentially have interacted rotationally. The ages he
determined for A and B are both 1.43 Gyr, but the age of C is 590 $\pm$ 70
Myr. He mentioned the possibility that A and B may have interacted and thereby
decreased their period. Our result agrees with this age discrepancy between
A/B and C. A and B ages are consistent for all the models (0.83 to 0.93 of the
entire distribution are within the tolerance level for all models). However, B
and C ages are inconsistent for all models (0 to 0.03 of the entire
distribution are within the tolerance level for all models.

The color uncertainties are small for these example systems. The rotation
period uncertainties are also small for the majority of them. However, our age
comparison method does not work well (nor would any other method work well) if
either rotation period or color uncertainties are large. Figure set 4
demonstrates this in detail. Since the period uncertainties for $\alpha$Cen A
are large ($\sigma_{P}/P$ = 0.11), these distributions are broader than most
of the other distributions. Thus, even though it may appear that the ages of
the two components are consistent in the sense that they overlap 1.0, there is
actually a greater probability that the ages are inconsistent, i.e. outside
the 1 $\pm~{}\Delta$ range, than consistent and within 1 $\pm~{}\Delta$. Most
of the original $Kepler$ periods have reported uncertainties around
$\sigma_{P}/P$ = 0.05 (Aigrain et al., 2015). However, for $K2$, typical
uncertainties are larger, perhaps because the observation span is much shorter
than the original $Kepler$ data (see Fig. 4 of Reinhold & Hekker (2020)).
Therefore, K2 data may be insufficient for precise enough periods,
particularly for slower rotators. The best $B-V$ uncertainties that can
generally be obtained from ground-based telescopes are $\sim$ 0.02 mag (Tonry
et al., 2012), however most of the $TESS$ Input Catalog (Stassun et al., 2018)
entries have $\sigma_{B-V}\geq$ 0.05. On the other hand, colors obtained from
Pan-STARRS and $Gaia$ often have substantially smaller uncertainties. Pan-
STARRS $g-r$ uncertainties are around 0.04 and $Gaia$ $G_{BP}-G_{RP}$
uncertainties are 0.002 at $G$ $<$ 13, and 0.01 at $G$ = 17 (Tonry et al.,
2012; Evans et al., 2018). Therefore, we suggest that the next generation of
gyrochronology relations be based on photometry from $Gaia$ and Pan-STARRS, or
other large-field, high-quality surveys, for example as done by A19.

Figure 4: Comparison of age ratio probability distributions of components of
binary 1 ($\xi$Boo A/B) whose connecting line does not parallel a gyrochrone
which does not parallel a gyrochrone in Figure 3. Each gyrochronology model
(A15, B07, MH08, M11, A19, & SL20) is used for each subplot. The percentage of
each distribution that is within the tolerance area is shown in parentheses in
the lower right in each panel. A ten percent tolerance ($\Delta$ = $\pm$0.10)
in age agreement between components of such binary is illustrated by the
vertical black lines. Comparison of age ratio probability distributions of
components of each binary in shown in the Figure set (5 images), which is
available in the online journal.

Fig. Set4. Comparison of age ratio probability distributions of components of
binaries

##  6 Conclusion

Main sequence stars with masses of less than 1.3 $M_{\sun}$ lose their angular
momentum as they age and thus their rotation periods increase over time.
Although current gyrochronology concepts work relatively well for young MS
stars with masses between 1.0 and 1.3 $M_{\sun}$, we have learned in the past
ten years that improvements to the models are required for older and lower
mass MS stars. In this paper, we propose a numerical model evaluation method
using wide binaries under the assumption that binary components are coeval. We
test six models that rely only on color and period information (B07, MH08,
M11, A15, A19, and SL21) as examples. The use of binary stars supplements open
cluster data because the abundance of binaries allows us to cover and fill in
a wider range of stellar age and metallicity. We use our model evaluation
method to determine the probability that binary components have consistent
ages, finding examples where component stars are consistent as well as
examples where ages are inconsistent. Both consistent and inconsistent pairs
can help the community improve data quality and model fitting, and ultimately
the precision of gyrochronology ages. In addition, we used a Monte Carlo
approach to assess the precision and accuracy of ages derived from four
Skumanich-based models (B07, MH08, M11, A15)). The age uncertainties that we
obtained are similar to those obtained by B07 and are $\sigma_{age}$/$age$ =
5-20 % for $\sigma_{P}$/$P$ = 2% and $\sigma_{B-V}$ = 0.01. Our results show
that gyrochronological age uncertainty is slightly higher than derived by B07
for stars with $B-V$ = 0.65 ($\sigma_{age}$/$age$ = 15-20%) and slightly
smaller for the stars with $B-V$ = 1.5 and rotation period $P\leq$ 20 days
($\sigma_{age}$/$age$ = 5-10%). We also found that with the traditional
Skumanich-based gyrochronology models of B07, MH08, M11, and A15,
gyrochronology ages can be determined to approximately 5 to 20% if
$\sigma_{B-V}\leq$ 0.01 and $\frac{\sigma_{P}}{P}\leq$ 0.02, that are the
smallest uncertainties among binary samples of B07 (Please note that this is
the precision claimed by these specific models and this is not the actual
precision).

An anonymous referee has provided valuable insight through questions posed in
the report on the draft of this paper, and for which the authors are most
grateful.
This material is based upon work supported by the National Science Foundation
under Grant No. AST-1715718, AST-1910396, and AST-2108975 and by NASA under
Grant No. NNX16AB76G and 80NSSC21K0245.

## References

  * Agüeros et al. (2018) Agüeros, M. A., Bowsher, E. C., Bochanski, J. J., et al., 2018, ApJ, 862, 33
  * Aigrain et al. (2015) Aigrain, S. et al., 2015, MNRAS, 450, 3226
  * Angus et al. (2015) Angus, R., Aigrain, S., Foreman-Mackey, D., & McQuillan, A., 2015, MNRAS, 450, 1787
  * Angus et al. (2019) Angus, R., Morton, T. D., Foreman-Mackey, D., & van Saders, J. et al. 2019, AJ, 158, 173
  * Angus et al. (2020) Angus, R., Beane, A., Price-Whelan, A. M., & Newton, E. et al. 2019, AJ, 160, 90
  * Barnes (2001) Barnes, S. A., 2001, ApJ, 561, 1095
  * Barnes (2003) Barnes, S. A., 2003, ApJ, 586, 464
  * Barnes (2007) Barnes, S. A., 2007, ApJ, 669, 1167
  * Barnes (2010) Barnes, S. A., 2010, ApJ, 722, 222
  * Barnes & Kim (2010) Barnes, S. A., & Kim, Y., 2010, ApJ, 721, 675
  * Epstein & Pinsonneault (2014) Epstein, Courtney R., & Pinsonneault, Marc H., 2014, ApJ, 780, 2
  * Bravo et al. (2014) Bravo, J.-P., Roque, S., Estrella, R., Leão, I. C., & De Medeiros, J. R., 2014, A&A, 568, A34
  * Brown (2014) Brown, T. M., 2014, ApJ, 789, 101
  * Buzasi (1997) Buzasi, D. L., 1997, ApJ, 484, 2
  * Canto Martins et al. (2020) Canto Martins, B. L., Gomes, R. L., Messias, Y. S., De Lira, S. R., et al. 2020, ApJS, 250, 20
  * Curtis et al. (2019) Curtis, J. L., Agüreos, M. A., Douglas, S. T., & Meibom, S., 2019, ApJ, 879, 49
  * Casagrande & VandenBerg (2018) Casagrande, L., & VandenBerg, D., 2018, MNRAS, 479, 102
  * Douglas et al. (2017) Douglas, S. T., Agüeros, M. A., Covey, K. R., & Kraus, A. 2017, ApJ, 842, 83
  * Douglas et al. (2019) Douglas, S. T., Curtis, J. L., Agüeros, M. A., et al. 2019, ApJ, 879, 100
  * Curtis et al. (2020) Curtis, J. L., Agüreos, M. A., Matt, S. P., Covey, K. R. et al. ApJS, 250, 20
  * Evans et al. (2018) Evans, D. W., Riello, M., De Angeli, F., Carrasco, J. M. et al., 2018, “$Gaia$ Data Release 2: Photometric Content and Validation,” arXiv:1804.09368
  * El-Badry, Rix, & Heintz (2021) El-Badry, K., Rix, H. -W., Heintz, T. M., 2021, MNRAS, 506, 2269 El Badry et al. (2021)
  * Flewelling et al. (2016) Flewelling, H., Magnier, E., Chambers, K. et al. 2016, “The Pan-STARRS1 Database and Data Products,” arXiv:1612.05243
  * Gaia Collaboration (2018) Gaia Collaboration, Brown, A. G. A., Vallenari, A., & Prusti, T., et al. 2018, A&A, 616, A1
  * Gallet & Bouvier (2013) Gallet, F., & Bouvier, J. 2013, A&A, 556, A36
  * Garraffo et al. (2018) Garraffo, C., Drake, J. J., Dotter, A., Choi, J. M., et al., 2018, ApJ, 862, 1
  * García et al. (2014) García, R. A., Ceillier, T., Salabert, D., & Mathur, S., et al., 2014, A&A, 572, A34
  * Godoy-Rivera (2018) Godoy-Rivera, D., & Chaname, J., 2018, MNRAS, 479, 4440
  * Gondoin (2017) Gondoin, P., 2017, A&A, 599, A122
  * Guerrero et al. (2019) Guerrero, G., Zaire, B., Smolarkiewicz, P. K., de Gouveia Dal Pino, E. M., Kosovichev, A. G., & Mansour, N. N. 2019, ApJ, 880, 6
  * Irwin & Bouvier (2009) Irwin, J, & Bouvier, J. 2009, in IAU Symposium, Vol. 258, The Ages of Stars, ed. E. E. Mamajek, D. R. Soderblom, & R. F. G. Wyse, 363-374
  * Janes (2017) Janes, K., 2017, ApJ, 835, 75
  * Jenkins et al. (2016) Jenkins, J., Twichen, J., McCauliff, S., Campbell, J. et al., 2016, Proceedings of the SPIE, 99133E
  * Kawaler (1988) Kawaler, S, D., 1988, ApJ, 333, 236
  * Kawaler (1989) Kawaler, S, D., 1988, ApJ, 343, L65
  * Kraus & Hillenbrand (2009) Kraus, Adam, M., & Hillenbrand, Lynne, A., 2009, ApJ, 704, 531
  * Lanzafame & Spada (2015) Lanzafame, A, C., & Spada, F., 2015, A&A, 584, A30
  * Matt, et al. (2015) Matt, S. P., Brun, A. S., Baraffe, I., et al. 2015, ApJ, 695, 679
  * McQuillan et al. (2013) McQuillan, A., Aigrain, S. & Mazeh, T., 2013, MNRAS, 432, 1203
  * McQuillan, Mazeh, & Aigrain (2014) McQuillan, A., Mazeh, T., & Aigrain, S., 2014, ApJS, 211, 24
  * Mamajek & Hillenbrand (2008) Mamajek, E. E., & Hilenbrand, L. A., 2008, ApJ, 687, 1264
  * Mason et al. (2001) Mason, B., Wycoff, G., Hartkopf, W., Geoffrey, D., & Worley, C., 2001, AJ, 122, 3466
  * Matt et al. (2012) Matt, S., MacGregor, K., Pinsonneault, M., & Greene, T., 2012, ApJ, 754, 26
  * Metcalfe & Egeland (2019) Metcalfe, T., & Egeland, R., 2019, ApJ, 871, 39
  * Meibom & Mathieu (2009) Meibom, S., & Mathieu, R., 2009, ApJ, 695, 679
  * Meibom et al. (2011a) Meibom, S., Mathieu, R., Stassun, K., Liebesny, P., & Saar, S., 2011, ApJ, 733, 115
  * Meibom et al. (2011b) Meibom, S. et al. 2011, ApJ, 733, L9
  * Meibom et al. (2015) Meibom, S., Barnes, S., Platais, I., Gilliland, R., Latham, D., Mathieu, R. 2015, Nature, 517, 589
  * Pallavicini et al. (1981) Pallavicini, R., Golub, L., Rosner, R., Vaiana, G., Ayres, T., Linsky, J., 1981, ApJ, 248, 279
  * Pinsonneault et al. (1989) Pinsonneault, M. H., Kawaler, S. D., Sofia, S., & Demarque, P., 1989, ApJ, 338, 424
  * Rebull (2017) Rebull, L. M., Stauffer, J. R., Hillenbrand, L. A., Cody, A. M., Bouvier, J., Soderblom, D. R., Pinsonneault, M., & Hebb, L., 2017, ApJ, 839, 92
  * Reinhold & Hekker (2020) Reinhold, Timo, & Hekker, Saskia, 2020, A&A, 635, 43
  * Scargle (1998) Scargle, J. D., 1998, ApJ, 504, 405
  * Skumanich (1972) Skumanich, A., 1972, ApJ, 171, 565
  * Soderblom (1985) Soderblom, D.R., 1985, AJ, 90, 2103
  * Soderblom (2010) Soderblom, D.R., 2010, ARA&A, 48, 581
  * Somers & Pinsonneault (2016) Somers, G., & Pinsonneault, M.H., 2016, ApJ, 829, 32
  * Spada & Lanzafame (2020) Spada, F., & Lanzafame, A, C., 2020, A&A, 636, A76
  * Stassun et al. (2018) Stassun, K. G., Oelkers, R. J., Pepper, J. Paegert, M. et al., 2018, AJ, 156, 3
  * Schatzman (1962) Schatzman, E., 1962, Annales d’Astrophysique, 25, 18
  * Tonry et al. (2012) Tonry, J. L., Stubbs, C. W., Lykke, K. R., et al., 2012, ApJ, 750, 99
  * van Saders & Pinsonneault (2013) van Saders, J. L., & Pinsonneault, M. H., 2013, ApJ, 776, 67
  * van Saders et al. (2016) van Saders, J., Ceillier, T., Metcalfe, T., Aguirre, V., Pinsonneault, M., García, R., Mathur, S., & Davies, G., 2016, Nature, 529, 181
  * van Saders, Pinsonneault, & Barbieri (2018) van Saders, J., Pinsonneault, M., & Barbieri, M., 2018, ApJ, 872, 2
  * VanderPlas (2018) VanderPlas, J. T., 2018, ApJS, 236, 16
  * Weber & Davis (1967) Weber, E. J., & Davis, Jr., L., 1967, ApJ, 148, 217
  * Zacharias (2015) Zacharias, N. et al., 2015, AJ, 150, 101