# Detection of Close Kuiper Belt Binaries with HST WFC3

Simon B. Porter Southwest Research Institute, 1050 Walnut St, Suite 300,
Boulder, CO 80302, USA  Susan D. Benecchi Planetary Science Institute, Tucson,
AZ, USA  Anne J. Verbiscer University of Virginia, Charlottesville, VA, USA
W. M. Grundy Lowell Observatory, Flagstaff, AZ, USA  Keith S. Noll NASA
Goddard Spaceflight Center, Greenbelt, MD, USA  Alex H. Parker SETI Institute,
Mountain View, CA, USA

(Received 5/31/2023; Revised 4/2/2024; Accepted 4/14/2024)

###### Abstract

Binaries in the Kuiper Belt are common. Here we present our analysis of the
Solar System Origins Legacy Survey (SSOLS) to show that using a PSF-fitting
method can roughly double the number of binaries identified in that dataset.
Out of 198 Kuiper Belt objects (KBOs) observed by SSOLS, we find 23 to be
visually separated binaries, while a further 19 are blended-PSF binaries
detectable with the method we present here. This is an overall binary fraction
of 21% for the SSOLS dataset of cold classical KBOs. In addition, we tested
our fitting methods on synthetic data, and while we were able to show it to be
very effective at detecting certain blended-PSF binary KBOs, fainter or closer
binary KBOs may easily be missed, suggesting that the close binary KBO
fraction could be even higher. These results strongly support the idea that
most (if not all) KBOs were formed through the Streaming Instability process,
and as a consequence, most KBOs were formed as near-equal mass binaries.

Classical Kuiper belt objects(250) — Trans-Neptunian objects(1705) — Asteroid
satellites(2207)

††journal: Planetary Science Journal††software:  Astropy (Astropy
Collaboration et al., 2022), Scipy (Virtanen et al., 2020), Photutils (Bradley
et al., 2022), Spiceypy (Annex et al., 2021), Tiny Tim (Krist, 1995),

##  1 Introduction

The Solar System is a highly evolved planetary system that betrays few clues
of its origin, with the key exception of the Kuiper Belt. The present day
Kuiper Belt consists almost exclusively of objects that formed outside of the
giant planets and were either pushed outward by giant planet migration (the
Neptune-resonant and scattered populations) or are still in the original
orbits where they formed around the Sun (the classical KBOs; McKinnon et al.,
2020). KBOs, and especially the cold classical KBOs (CCKBOs) which have the
least perturbed heliocentric orbits, have suffered very little bombardment
since their formation relative to asteroids interior of Jupiter (Mao et al.,
2021), and no real thermal processing like comets (Grundy et al., 2020). The
best constraints on the surfaces of KBOs come from the New Horizons mission,
which flew past both the Pluto system and the CCKBO Arrokoth. Pluto itself is
a highly-processed world with extensive surface-atmosphere interactions (Stern
et al., 2015). Pluto’s large satellite Charon and two of four small satellites
(Nix and Hydra) were imaged at sufficient resolution to show a small amount of
cratering (Robbins et al., 2017), but since all the satellites were likely
formed after a giant impact (Canup, 2011), that says little about Kuiper Belt
history as a whole. Arrokoth, a typical CCKBO (Porter et al., 2018), was much
more constraining. It is constructed of two distinct lenticular lobes,
connected on their equators, forming a contact binary (Spencer et al., 2020).
There is a single large crater on the smaller lobe, and all other craters are
significantly smaller (Spencer et al., 2020). This clean contact binary shape
is unlike any object previously encountered by a spacecraft, but provides
tantalizing clues to the formation processes of the solar system (McKinnon et
al., 2020).

A leading theory of the formation of the Outer Solar System is the Streaming
Instability (SI; Nesvorný et al., 2019). This postulates that solid bodies in
the circumsolar disk initially formed when small clumps of solid material
caused localized increases in gas drag, allowing the small clumps to very
rapidly grow to full sized KBOs within tens of years (Nesvorný et al., 2019,
and citations therein). An important consequence of this rapid growth (when
compared to more traditional slow accretion models) is that angular momentum
grows rapidly too. SI thus tends to produce binary systems (and occasionally
triple and higher systems), with most of the angular momentum in the mutual
orbits (Fraser et al., 2017). In addition, because the source angular momentum
is primarily from Keplerian shear, binaries produced by SI should be generally
prograde (Nesvorný et al., 2019). And because SI requires the gas disk to
happen, some fraction of SI-created binaries should have evolved inwards under
gas drag until they became contact binaries (as may have happened for
Arrokoth; McKinnon et al., 2020), while others evolved under mutual and solar
tides to become tight circular binaries (Porter & Grundy, 2012). SI thus
predicts a Kuiper Belt filled with binaries of varying mutual separations,
primarily prograde orbits, as well as a significant fraction of contact
binaries.

The Kuiper Belt actually observed does indeed contain many binaries, which are
primarily on prograde orbits (Grundy et al., 2019). Photometric lightcurve
studies appear to show that at least 15% of CCKBOs are contact binaries
(Thirouin & Sheppard, 2019), though this number could be an underestimate due
the difficulty of detecting contact binaries from lightcurves (Showalter et
al., 2021). Similar studies of the (closer and brighter) 3:2 Neptune resonant
“Plutinos” found a contact binary fraction of up to 50% (Thirouin & Sheppard,
2018). Arrokoth is particularly instructive, as its dual-lenticular contact
binary shape is best explained by formation as a binary by SI and slow
evolution to contact under gas drag (McKinnon et al., 2020). The case for SI
to be the leading formation mechanism for the Kuiper Belt, and the Outer Solar
System generally, is thus very strong. However, a key prediction of SI is
there should be a very large number of very tight ($<$3000 km separation)
binary KBOs. The overall fraction of KBOs represented by these tight binaries
is unknown, but must be significant to hide a large fraction of SI-produced
binaries below the angular resolution of HST and ground-based adaptive optics
systems. Measuring the separations of very tight binary KBOs (and especially
CCKBOs) thus provides critical constraints not only on the SI, but on the
processes that have sculpted the orbits of binary KBOs since then, including
solar perturbations (Porter & Grundy, 2012) and gas drag (McKinnon et al.,
2020).

The Wide Field Camera 3 (WFC3) on the Hubble Space Telescope (HST) offers the
best available tool to detect binary KBOs, as it is able to observe very faint
objects with very high angular precision. Here we show we can detect and
characterize binaries that are so close that they not clearly separated in
HST/WFC3 images. This is aided by the fact that the WFC3 point spread function
(PSF) is on average consistent and very well characterized, allowing for a
precise measurement of a binary KBO with a separation of at least 0.5 pixels,
as and bright as V=24. We apply this method to 198 CCKBOs observed with HST
program GO-15648 (the Solar System Origins Legacy Survey, SSOLS), and show
that it can roughly double the yield of binaries relative to those which can
be easily detected with visual inspection. We also use this model to create
1000 synthetic KBO cases implanted in real HST images, and show the limits of
this detection method.

Figure 1:  Examples of (top) a well-separated binary KBO, (middle) a blended
binary KBO, and (bottom) a KBO that appears to be single. The top separated
binary case accounted for 11.6% of the KBOs in the SSOLS dataset. For the
middle case, the single-PSF solution only subtracts the primary, leaving
behind a significant residual which is not there for the double-PSF solution,
strongly indicating that this KBO is a blended binary. The residuals for the
bottom case are almost identical for the single-PSF and double-PSF solutions,
indicating that it is not binary to the resolution of HST/WFC3. All of these
images are the F606W stacks for these objects; we performed the same analysis
for the F814W stacks and enforced consistency between them for objects that we
identified as blended binaries.

##  2 Target Selection and HST Imaging

The Cycle 26 Treasury Program ‘The Solar System Origins Legacy Survey’ (SSOLS)
observed 198 CCKBOs at optical wavelengths. This sample drew on three
successful ground-based discovery surveys, the Deep Ecliptic Survey (DES;
Elliot et al., 2005), the Canada France Ecliptic Plane Survey (CFEPS; Petit et
al., 2011), and the Outer Solar System Origins Survey (OSSOS; Bannister et
al., 2018). The DES was the first deep systematic survey to probe the
population of the Kuiper Belt, and discovered 320 KBOs and Centaurs, including
43 of the SSOLS targets (Elliot et al., 2005). CFEPS probed 321 degree2 of the
sky, discovering and characterizing 169 KBOs, discovering 4 of the SSOLS
targets (Petit et al., 2011). OSSOS was a followup of CFEPS, discovering over
800 KBOs and Centaurs, including 151 of the SSOLS targets (Bannister et al.,
2018). The SSOLS targets were chosen to be cold classical KBOs that did not
appear to be binary in the ground-based imagery, with the exception of 2005
EO304, the widest binary in the SSOLS target list. This sample set of CCKBOs
was designed to provide identifiable CCKBO binaries at a level that allows
discovery bias to be quantified and removed, to better understand the binarity
and luminosity function of the intrinsic CCKBO population. All of the targets
(listed in Table 1) had well-determined orbits with at least three oppositions
of data, to allow confidence in both their dynamical classification and that
would be observable by the Hubble Space Telescope (HST) with a small subarray.

The 198 KBOs were successfully observed with HST program GO-15648, ranging in
approximate apparent magnitude from V=21 to V=25. Observations of an
additional 14 objects were attempted, but failed due to HST not correctly
locking on to its guide star, leaving the observations unusable. The target
CCKBOs were between 35.98 and 53.43 astronomical units from Earth at the time
they were imaged, corresponding to a spatial resolution in the WFC3 images of
1044 to 1550 km/pixel. Each target KBO was observed with a single orbit, with
four images in the F606W filter (roughly V-band), and four images in the F814W
filter (roughly R-band), all using the UVIS2-C512C-SUB subframe. A four-point
box dither was used for each filter. These filters were chosen to provide the
most amount of color information while still maximizing the signal-to-noise
ratio (SNR) for binary detection. F606W and F814W have also been used for
previous HST observations of CCKBOs, particularly in support of the New
Horizons mission (Benecchi et al., 2019), allowing these observations to place
both the color distribution of the binaries in context of the broader cold
classical Kuiper Belt and the well-characterized New Horizons Distant KBO
subset. Here we only present the results of relative astrometry and binary
identification, while a following paper (in prep) will present the color
results. An additional paper (in prep) will use these binaries to present a
model for the size frequency distribution of binaries in the CCKBOs.

All of the data presented in this paper were obtained from the Mikulski
Archive for Space Telescopes (MAST) at the Space Telescope Science Institute.
The specific observations analyzed can be accessed via
https://doi.org/10.17909/6xbd-kd20 (catalog doi:10.17909/6xbd-kd20).

Table 1:  The SSOLS target CCKBOs, their F606W magnitudes, their combined relative probabilities of being binary (Equation 5), the best-fit separation of the double-PSF case, the best-fit flux ratio of the the double-PSF case, and our judgement if the object is a separated binary, blended binary, or single.  |  | F606W Mag | Bin. Prob. | Pix. Sep. | Flux Ratio | Binary?
---|---|---|---|---|---|---
(612087) | 1999 CU153 | 23.4 | 0.891 | 1.0 | 2.0:1 |
(439858) | 1999 ON4 | 23.6 | 0.425 | 1.1 | 2.5:1 |
(612157) | 2000 FG8 | 24.1 | 0.634 | 0.9 | 2.8:1 |
(455206) | 2001 FE193 | 23.0 | 0.082 | 1.5 | 5.7:1 |
(612213) | 2001 FK185 | 23.3 | 0.022 | 2.1 | 15.2:1 |
(385362) | 2002 PT170 | 23.3 | 0.014 | 2.7 | 9.0:1 |
(469509) | 2003 HC57 | 22.2 | 1.000 | 1.5 | 1.9:1 | Blended
(612549) | 2003 HG57 | 22.8 | 1.000 | 7.6 | 1.5:1 | Separated
(363401) | 2003 LB7 | 22.6 | 0.007 | 1.9 | 100.0:1 |
(612619) | 2003 SN317 | 22.9 | 0.998 | 0.8 | 2.0:1 | Blended
| 2003 UK293 | 23.4 | 0.186 | 1.2 | 10.2:1 |
| 2003 UV291 | 24.0 | 0.006 | 1.8 | 100.0:1 |
(612733) | 2003 YU179 | 23.0 | 1.000 | 1.1 | 1.7:1 | Blended
(444018) | 2004 EU95 | 23.0 | 0.001 | 6.2 | 100.0:1 |
| 2004 HD79 | 22.5 | 1.000 | 2.1 | 1.4:1 | Separated
| 2004 HE79 | 23.5 | 0.717 | 0.9 | 2.5:1 |
(469610) | 2004 HF79 | 22.9 | 1.000 | 3.2 | 1.4:1 | Separated
| 2004 HG79 | 23.1 | 0.629 | 1.0 | 2.6:1 |
(444025) | 2004 HJ79 | 23.3 | 0.958 | 1.3 | 2.6:1 | Blended
| 2004 HK79 | 23.4 | 1.000 | 3.1 | 1.1:1 | Separated
| 2004 KE19 | 22.7 | 1.000 | 3.3 | 2.2:1 | Separated
| 2004 KF19 | 22.7 | 0.043 | 4.1 | 54.2:1 |
| 2004 KG19 | 23.3 | 0.841 | 0.8 | 2.4:1 |
| 2004 MU8 | 23.0 | 1.000 | 4.3 | 1.3:1 | Separated
| 2004 PT117 | 23.1 | 0.984 | 1.0 | 2.3:1 | Blended
| 2004 PU117 | 23.2 | 0.089 | 1.1 | 6.6:1 |
| 2004 PV117 | 22.8 | 1.000 | 5.8 | 5.6:1 | Separated
| 2004 PW117 | 23.3 | 1.000 | 17.5 | 2.0:1 | Separated
| 2004 PX117 | 23.5 | 1.000 | 4.4 | 1.5:1 | Separated
(609222) | 2004 VB131 | 22.9 | 0.903 | 0.9 | 2.4:1 |
(609221) | 2004 VC131 | 22.0 | 0.294 | 1.3 | 2.2:1 |
| 2004 VD131 | 23.0 | 0.974 | 0.8 | 1.5:1 | Blended
| 2005 BW49 | 22.8 | 1.000 | 2.4 | 1.7:1 | Separated
| 2005 CE81 | 23.8 | 1.000 | 6.1 | 1.1:1 | Separated
(525461) | 2005 EN302 | 23.7 | 0.001 | 5.4 | 100.0:1 |
(525462) | 2005 EO304 | 22.6 | 1.000 | 43.1 | 4.2:1 | Separated
| 2006 CH69 | 23.2 | 1.000 | 23.0 | 1.7:1 | Separated
| 2006 JV58 | 23.3 | 1.000 | 3.6 | 1.7:1 | Separated
| 2006 QA181 | 23.5 | 0.010 | 3.0 | 9.9:1 |
(587670) | 2006 QE181 | 23.8 | 0.658 | 0.6 | 1.7:1 |
| 2006 QF181 | 23.4 | 0.999 | 0.8 | 1.6:1 | Blended
(523615) | 2006 UO321 | 23.7 | 0.088 | 1.4 | 9.4:1 |
| 2006 WF206 | 23.1 | 1.000 | 1.3 | 1.3:1 | Blended
| 2007 CQ79 | 23.4 | 0.884 | 1.1 | 3.6:1 |
| 2007 CS79 | 23.3 | 0.466 | 0.9 | 2.8:1 |
| 2007 DS101 | 22.8 | 1.000 | 2.1 | 1.3:1 | Separated
| 2013 EM149 | 23.0 | 1.000 | 1.3 | 2.0:1 | Blended
(500840) | 2013 GA138 | 23.9 | 0.483 | 0.7 | 1.4:1 |
| 2013 GB138 | 23.6 | 0.459 | 0.9 | 4.0:1 |
| 2013 GC138 | 23.6 | 0.219 | 1.9 | 12.4:1 |
| 2013 GD138 | 24.3 | 0.024 | 2.1 | 5.5:1 |
| 2013 GE138 | 24.3 | 0.053 | 1.5 | 5.9:1 |
| 2013 GF138 | 23.5 | 0.894 | 1.1 | 3.7:1 |
(500835) | 2013 GN137 | 22.9 | 0.997 | 1.2 | 2.5:1 | Blended
| 2013 GP137 | 23.8 | 0.005 | 1.7 | 100.0:1 |
(500836) | 2013 GQ137 | 23.5 | 0.076 | 1.5 | 5.9:1 |
| 2013 GR137 | 24.2 | 0.284 | 0.6 | 3.2:1 |
| 2013 GS137 | 23.9 | 0.317 | 0.9 | 3.1:1 |
(500837) | 2013 GT137 | 24.0 | 0.369 | 1.0 | 6.2:1 |
| 2013 GU137 | 23.7 | 0.668 | 0.7 | 1.5:1 |
(500838) | 2013 GV137 | 23.3 | 0.920 | 0.7 | 2.0:1 |
(500839) | 2013 GW137 | 23.8 | 0.070 | 1.6 | 4.5:1 |
| 2013 GX137 | 23.5 | 0.586 | 0.9 | 2.6:1 |
| 2013 GY137 | 23.7 | 0.130 | 0.8 | 2.5:1 |
(500856) | 2013 HT156 | 24.7 | 0.009 | 4.7 | 6.4:1 |
| 2013 SC101 | 24.1 | 0.036 | 0.9 | 100.0:1 |
| 2013 SD101 | 23.8 | 0.039 | 1.7 | 6.3:1 |
| 2013 SE101 | 24.1 | 0.340 | 1.0 | 2.8:1 |
| 2013 SF101 | 23.9 | 0.046 | 1.5 | 3.8:1 |
| 2013 SG101 | 24.5 | 0.042 | 2.4 | 11.4:1 |
| 2013 SJ100 | 24.5 | 0.628 | 0.9 | 6.8:1 |
| 2013 SL100 | 23.8 | 0.854 | 0.9 | 2.2:1 |
| 2013 SO100 | 24.1 | 0.073 | 2.1 | 12.6:1 |
(505446) | 2013 SP99 | 23.4 | 0.578 | 0.8 | 2.3:1 |
(505447) | 2013 SQ99 | 23.4 | 1.000 | 3.4 | 1.5:1 | Separated
| 2013 SV100 | 23.9 | 0.766 | 0.8 | 2.3:1 |
| 2013 SX100 | 24.2 | 0.028 | 2.3 | 4.2:1 |
| 2013 UB18 | 24.2 | 0.076 | 2.2 | 6.5:1 |
| 2013 UC18 | 24.5 | 0.040 | 1.8 | 3.3:1 |
| 2013 UD18 | 23.7 | 0.005 | 1.8 | 100.0:1 |
| 2013 UG18 | 24.3 | 0.090 | 0.5 | 5.3:1 |
(505476) | 2013 UL15 | 23.2 | 1.000 | 4.3 | 1.6:1 | Separated
| 2013 UL17 | 24.4 | 0.001 | 5.9 | 100.0:1 |
| 2013 UN15 | 23.8 | 0.712 | 0.6 | 2.2:1 |
| 2013 UN17 | 24.4 | 0.116 | 1.0 | 6.5:1 |
| 2013 UO15 | 23.3 | 0.796 | 0.7 | 1.8:1 |
| 2013 UP15 | 23.6 | 0.842 | 0.6 | 2.1:1 |
| 2013 UP17 | 24.3 | 0.137 | 0.8 | 2.2:1 |
| 2013 UR17 | 24.1 | 0.917 | 0.6 | 1.9:1 |
| 2013 UT17 | 24.0 | 0.337 | 0.8 | 5.0:1 |
| 2013 UW16 | 23.6 | 0.970 | 1.1 | 1.4:1 | Blended
| 2013 UW17 | 24.3 | 0.325 | 0.9 | 5.2:1 |
| 2013 UY16 | 24.1 | 0.664 | 0.8 | 2.8:1 |
| 2013 UY17 | 23.8 | 0.022 | 2.4 | 4.0:1 |
| 2014 UC228 | 24.2 | 0.394 | 1.3 | 1.6:1 |
| 2014 UC229 | 24.7 | 0.436 | 0.8 | 3.6:1 |
| 2014 UC230 | 24.5 | 0.129 | 1.7 | 8.9:1 |
(511551) | 2014 UD225 | 23.1 | 0.984 | 0.9 | 2.1:1 | Blended
(511552) | 2014 UE225 | 22.7 | 0.768 | 0.5 | 1.4:1 |
| 2014 UL228 | 24.0 | 0.113 | 1.0 | 4.1:1 |
| 2014 UL229 | 24.3 | 0.003 | 6.0 | 17.3:1 |
| 2014 UP228 | 24.0 | 0.009 | 3.3 | 10.7:1 |
| 2014 UY228 | 24.2 | 0.033 | 2.0 | 4.0:1 |
(523756) | 2014 WD509 | 22.5 | 1.000 | 6.4 | 2.0:1 | Separated
| 2015 GA57 | 24.2 | 0.121 | 3.5 | 30.7:1 |
| 2015 GB57 | 24.5 | 0.107 | 0.9 | 3.6:1 |
| 2015 GC57 | 24.0 | 0.480 | 1.2 | 2.7:1 |
| 2015 GC58 | 24.3 | 0.962 | 0.8 | 5.1:1 | Blended
| 2015 GD57 | 23.8 | 0.022 | 0.8 | 100.0:1 |
| 2015 GD59 | 24.2 | 0.048 | 0.6 | 19.4:1 |
| 2015 GE57 | 24.5 | 0.006 | 1.7 | 100.0:1 |
| 2015 GF56 | 24.3 | 0.018 | 1.9 | 12.3:1 |
| 2015 GF58 | 23.6 | 0.034 | 2.0 | 10.1:1 |
| 2015 GF59 | 24.1 | 0.355 | 1.6 | 4.7:1 |
| 2015 GG57 | 24.6 | 0.464 | 0.9 | 2.5:1 |
| 2015 GH58 | 23.8 | 0.992 | 1.5 | 9.7:1 | Blended
| 2015 GJ57 | 24.2 | 0.013 | 2.7 | 13.0:1 |
| 2015 GK58 | 24.6 | 0.383 | 1.1 | 4.9:1 |
| 2015 GL57 | 24.3 | 0.135 | 2.2 | 7.6:1 |
| 2015 GL58 | 24.3 | 0.036 | 2.1 | 4.7:1 |
| 2015 GM58 | 24.0 | 0.327 | 1.1 | 4.1:1 |
| 2015 GN58 | 24.3 | 0.154 | 2.3 | 9.4:1 |
| 2015 GO57 | 24.3 | 0.563 | 0.8 | 3.1:1 |
| 2015 GO58 | 24.6 | 0.452 | 1.2 | 2.5:1 |
| 2015 GP58 | 24.0 | 0.979 | 1.1 | 1.3:1 | Blended
| 2015 GR56 | 24.1 | 0.596 | 0.8 | 3.4:1 |
| 2015 GR57 | 24.0 | 0.008 | 3.3 | 12.3:1 |
| 2015 GS56 | 24.5 | 0.853 | 1.0 | 1.9:1 |
| 2015 GS57 | 24.0 | 0.150 | 1.1 | 3.2:1 |
| 2015 GT58 | 24.4 | 0.204 | 1.4 | 3.1:1 |
| 2015 GU56 | 23.5 | 0.952 | 0.8 | 1.7:1 | Blended
| 2015 GU57 | 23.9 | 0.006 | 1.6 | 100.0:1 |
| 2015 GU58 | 23.4 | 0.102 | 1.4 | 8.3:1 |
| 2015 GW56 | 24.4 | 0.037 | 1.4 | 6.5:1 |
| 2015 GW57 | 24.5 | 0.011 | 6.3 | 100.0:1 |
| 2015 GX56 | 24.3 | 0.152 | 2.0 | 7.7:1 |
| 2015 GY56 | 24.6 | 0.344 | 0.9 | 5.9:1 |
| 2015 GY57 | 24.6 | 0.102 | 1.8 | 6.8:1 |
| 2015 GZ56 | 23.5 | 0.007 | 6.6 | 100.0:1 |
| 2015 GZ57 | 24.6 | 0.018 | 0.5 | 100.0:1 |
| 2015 RA280 | 25.1 | 0.525 | 0.8 | 3.0:1 |
| 2015 RB280 | 24.3 | 1.000 | 5.6 | 1.7:1 | Separated
| 2015 RB281 | 23.5 | 0.029 | 1.8 | 12.0:1 |
| 2015 RC281 | 24.1 | 0.026 | 1.9 | 7.4:1 |
| 2015 RD280 | 24.4 | 0.015 | 0.9 | 100.0:1 |
| 2015 RE280 | 24.6 | 0.025 | 2.4 | 4.5:1 |
| 2015 RJ277 | 23.1 | 1.000 | 1.2 | 2.7:1 | Blended
| 2015 RO281 | 23.0 | 1.000 | 1.2 | 6.2:1 | Blended
| 2015 RP279 | 24.5 | 0.040 | 1.6 | 8.4:1 |
| 2015 RP280 | 24.0 | 1.000 | 3.9 | 1.8:1 | Separated
| 2015 RP281 | 25.2 | 0.439 | 0.5 | 1.4:1 |
| 2015 RT245 | 23.6 | 1.000 | 12.2 | 1.1:1 | Separated
| 2015 RZ279 | 24.2 | 0.017 | 3.2 | 5.0:1 |
| 2015 VA169 | 24.3 | 0.715 | 1.1 | 3.2:1 |
| 2015 VA172 | 23.6 | 0.045 | 2.0 | 10.9:1 |
| 2015 VB169 | 24.6 | 0.001 | 6.1 | 100.0:1 |
| 2015 VB170 | 23.7 | 0.874 | 1.0 | 2.3:1 |
| 2015 VB171 | 24.6 | 0.164 | 1.1 | 4.1:1 |
| 2015 VB173 | 24.3 | 0.005 | 1.6 | 100.0:1 |
| 2015 VC170 | 23.6 | 0.286 | 0.8 | 4.5:1 |
| 2015 VC172 | 23.8 | 0.003 | 6.2 | 100.0:1 |
| 2015 VD169 | 24.1 | 0.005 | 1.7 | 100.0:1 |
| 2015 VE169 | 24.0 | 0.019 | 2.3 | 6.1:1 |
| 2015 VF169 | 24.5 | 0.036 | 1.7 | 3.4:1 |
| 2015 VG169 | 24.7 | 0.040 | 1.4 | 7.0:1 |
| 2015 VH169 | 25.2 | 0.003 | 6.9 | 58.1:1 |
| 2015 VH171 | 23.7 | 0.143 | 1.1 | 4.4:1 |
| 2015 VH173 | 23.8 | 0.024 | 1.6 | 7.6:1 |
| 2015 VJ170 | 24.7 | 0.001 | 6.0 | 100.0:1 |
| 2015 VK169 | 24.9 | 0.032 | 0.4 | 100.0:1 |
| 2015 VK170 | 24.0 | 0.023 | 2.2 | 6.8:1 |
| 2015 VL171 | 24.5 | 0.014 | 0.6 | 100.0:1 |
| 2015 VM173 | 23.2 | 1.000 | 4.0 | 1.3:1 | Separated
| 2015 VN171 | 24.3 | 0.003 | 5.9 | 10.9:1 |
| 2015 VN172 | 24.3 | 0.265 | 0.9 | 2.8:1 |
| 2015 VO171 | 24.1 | 0.269 | 1.0 | 4.2:1 |
| 2015 VP168 | 24.4 | 0.547 | 0.8 | 1.8:1 |
| 2015 VP172 | 24.4 | 0.100 | 1.1 | 5.7:1 |
| 2015 VP173 | 24.3 | 0.012 | 3.7 | 9.2:1 |
| 2015 VQ168 | 24.4 | 0.002 | 4.6 | 100.0:1 |
| 2015 VQ169 | 23.9 | 0.712 | 0.8 | 2.6:1 |
| 2015 VQ172 | 23.6 | 0.121 | 1.8 | 9.1:1 |
| 2015 VQ173 | 23.7 | 0.824 | 1.0 | 2.8:1 |
| 2015 VR168 | 24.5 | 0.289 | 1.2 | 4.5:1 |
| 2015 VR172 | 23.6 | 0.315 | 1.0 | 8.6:1 |
| 2015 VS168 | 24.3 | 0.008 | 3.7 | 6.2:1 |
| 2015 VS172 | 24.3 | 0.232 | 0.8 | 2.0:1 |
| 2015 VT168 | 22.9 | 0.997 | 0.9 | 3.4:1 | Blended
| 2015 VU168 | 24.1 | 0.008 | 5.5 | 17.2:1 |
| 2015 VU169 | 23.8 | 0.064 | 1.6 | 7.6:1 |
| 2015 VU171 | 24.0 | 0.006 | 1.9 | 100.0:1 |
| 2015 VW168 | 24.7 | 1.000 | 4.8 | 1.8:1 | Separated
| 2015 VW170 | 24.4 | 0.005 | 9.7 | 5.0:1 |
| 2015 VW172 | 24.8 | 0.362 | 1.6 | 6.5:1 |
| 2015 VX169 | 24.2 | 0.088 | 1.5 | 5.0:1 |
| 2015 VY170 | 23.9 | 0.724 | 0.9 | 2.0:1 |
| 2015 VY172 | 24.1 | 0.069 | 1.6 | 8.1:1 |
| 2015 VZ169 | 24.6 | 0.051 | 1.6 | 4.4:1 |

##  3 Image Analysis

The first step in our image analysis was to reproject all of the images, in
both filters, on to a common frame that cancelled out the dither pattern.
Pixels affected by saturation, crosstalk, or CTE tails (as identified in the
data quality part of the FLC files) were masked out here and in subsequent
analysis. For most orbits, a 300x300 pixel initial image was sufficient to
find the KBO, but a few required a 700x700 pixel window (i.e. the full frame
with ample margin for dither). We then manually selected the object from this
stacked image and recorded the offset to allow subsequent images to be
properly centered. We next used the offsets to reproject the images again now
centered on the objects with a smaller (181 pixel) window, upsampled from the
original image by a factor of two. These second stacks were used to check if
the object was a well-separated binary, as was the case for 23 of 198 KBOs. In
the case that it was a separated binary, we recorded the positions of the
assumed primary and secondary objects, and if no secondary was obvious, just
the primary object. We then produced the final reprojected images, centered on
the midpoint between the objects if a secondary was manually identified, or on
the primary object if not. These final images were separately stacked for each
of the two filters, F606W and F814W.

We then fit the separate stacks for each filter with a model Point Spread
Function (PSF) from the Tiny Tim package (Krist, 1995). We independently fit
both of the filters with both a single-PSF and double-PSF model, and compared
the resulting $\chi^{2}$. The image $\chi^{2}$ was estimated as the sum of the
square of the difference between the model and the stack of images. For the
single-PSF cases, we used the manual pick of the primary as the initial
location for the PSF, and then optimized the flux and x/y offset of the model
PSF to minimise the image $\chi^{2}$. If the object had been manually
identified as a binary in the prior steps, that was used as the initial
conditions for the location of the secondary for the double-PSF cases. If
there were not an identified secondary object, we used an initial guess of the
brightest pixel of the residuals from the single-PSF fit.

If the double-PSF model showed little improvement over the single-PSF model,
the two were held as equaly likely, while if the double-PSF model was a
significant improvement over the single-PSF model, we considered that to be
likely a binary. We quantified this as:

$\displaystyle\Delta_{prob}=exp(-0.5\times(\chi^{2}_{double}-\chi^{2}_{single}))$
(1) $\displaystyle rprob_{1}=\Delta_{prob}/(\Delta_{prob}+1.)$ (2)

Where $\Delta_{prob}$ is the relative likelihood of the double-PSF solution
versus the single-PSF solution, and $rprob_{1}$ is the likelihood that the KBO
is a binary based on the $\chi^{2}$ residuals. If $\chi^{2}_{double}$ is
almost the same as $\chi^{2}_{single}$, then $\Delta_{prob}$ is close to
unity, and $rprob_{1}$ is close to 0.5, as there is little evidence that
either solution is superior. However, if $\chi^{2}_{double}$ is much smaller
than $\chi^{2}_{single}$, then $\Delta_{prob}$ is much larger than unity, and
$rprob_{1}$ closer to 1.0, as there is strong evidence that the double-PSF
solution is superior.

While this comparison of the $\chi^{2}$ values is useful for filtering out the
cases that were not improved with the double-PSF model, it is susceptible to
low signal-to-noise ratio cases where the second PSF is fitting a coherent
noise spike. We addressed this by fitting both filters independently, and then
comparing the results for consistency. We deweighted the probability of being
a binary for any case where the location of the secondary was larger than 1/4
pixel between the filters, or where the difference in the primary and
secondary delta magnitudes was more than 1. The quarter-pixel test was based
on a few tests with both real and synthetic data. We quantified this as:

$\displaystyle rprob_{2}=0.25/|pixsep_{F606W}-pixsep_{F814W}|$ (3)
$\displaystyle rprob_{3}=1.0/|{\Delta}mag_{F606W}-{\Delta}mag_{F814W}|$ (4)
$\displaystyle rprob_{total}=rprob_{1}\times rprob_{2}\times rprob_{3}$ (5)

Where $pixsep_{F606W}$ and $pixsep_{F814W}$ are the primary-secondary
separation distance in pixels for each of the two filters,
${\Delta}mag_{F606W}$ and ${\Delta}mag_{F814W}$ are the primary-secondary
differences in magnitude for the two filters, and $rprob_{total}$ is the
combined estimated probability of the KBO being binary. The constraint for
${\Delta}mag$ consistency is weaker than for $pixsep$ to allow for the
secondaries to have different F606W-F814W colors than primaries, but still
useful to filter out inconsistent solutions.

##  4 Calibration of Detection Efficiency

Figure 2:  Cumulative distributions of the brightnesses and physical
separations of the SSOLS binaries in Table 1 (blue, solid) in comparison with
the known KBO binaries (orange, dashed) listed at
http://www2.lowell.edu/users/grundy/tnbs/status.html. The semimajor axes of
the SSOLS binaries is estimated as their apparent separation times $\sqrt{2}$.
The SSOLS binaries are both fainter and tighter than the known binary KBOs.
The only binary KBOs that are tighter than the SSOLS blended binaries are the
two detected by New Horizons (2011 JY31 and 2014 OS393; Weaver et al., 2022)
and the third by a stellar occultation (2014 WC510; Leiva et al., 2020).

We sought to calibrate the detection efficiency of very close binary KBOs (and
very close binary minor planets in general) with a set of simulations using
the same fitting techniques as described above on synthetic data.
Specifically, we selected four orbits from the SSOLS dataset (idy601, idy62h,
idy64j, and idy693) that had a corner that was clear of stars and any
diffraction spikes from the target KBO. For each of these four datasets, we
then generated 250 synthetic model binary KBOs and added them to the real data
in the F606W band, so that each synthetic binary KBO had real read noise and
cosmic rays. We varied brightnesses of the primaries of the binaries from V=21
to V=26, the separation of the binaries from 0 to 3 WFC3 pixels, and the flux
ratio between primary and secondary objects was tested at 1:1, 2:1, 5:1, and
10:1. The model assumed 0.1 pixels 1-$\sigma$ of jitter per exposure.

Our simulation results can be seen in Figure 3. Equal-brightness binaries (top
left of Figure 3) with V=21 combined brightness (roughly 400 km diameter for
both objects at 40 AU and 5% albedo, same for following estimates) can be
easily detected at separations as close as a half WFC3 pixel (roughly 600 km
at 40 AU). This drops off with brightness, with a full pixel of separation
being needed at V=22.5 ($\approx$200 km diameters, $\approx$1200 km
separation), and effectively no equal-brightness blended binaries were
detectable at a combined magnitude less than V=24.5 (roughly 80 km diameter
for both objects). HST/WFC3 is thus able to probe to the Roche limit for
equal-brightness binaries larger than roughly 400 km, but is ineffective at
detecting equal-brightness binaries smaller than 80 km diameter at separations
less than at least 90 radii apart. This is crucial to understanding the
missing tight binary population, as (523764) 2014 WC510 was observed by
stellar occultation to have a separation of 349 km, $<$4 radii apart (Leiva et
al., 2020), and 2011 JY31 was observed by New Horizons to have a separation of
198 km, a few radii apart (Weaver et al., 2022). In addition, the simulations
of Porter & Grundy (2012) showed there should be a large number of binary KBOs
at separations of a few radii, as that is where tidal forces circularized
orbits and stopped semimajor axis decay. HST/WFC3 is thus generally
insensitive to a potentially large population of binary KBOs that is both
known to exist from other detection methods, and predicted by theory.

Binary KBOs with a 2:1 brightness ratio fared only slightly worse, with a
similar binary recovery rate for V=21 (primary roughly 500 km diameter,
secondary roughly 300 km), down to half a pixel (roughly 600 km), and V=22.5
(primary roughly 240 km, secondary 170 km) also needing a full pixel of
separation (1200 km). However, the drop to a 2:1 flux ratio increased the
threshold of any blended binary detection at 3 pixels from V=24.5 to V=24.0
(120 and 85 km diameters). At 5:1 brightness ratio, the minimum combined
brightness for the detection of a satellite was V=23 (primary diameter about
210 km, secondary 95 km). At 10:1 brightness ratio, the minimum combined
brightness for the detection of a satellite was V=21.5 (primary diameter about
440 km, secondary 140 km).

##  5 Results and Discussion

Figure 3:  Recovery efficiency of simulated binary KBOs with SSOLS-style
synthetic data, assuming jitter of 0.1 pixels and F606W-F814W color of 0.72.
Each panel corresponds to 1000 simulations with a different primary/secondary
flux ratio. The synthetic data was processed in the same way as the actual
SSOLS data. The separations of the synthetic binaries were varied uniformly
from 0 to 3 pixels, but binned and averaged here for clarity. Recovery
efficiency was as expected better for brighter satellites farther away, with
very few (but still some) 10:1 flux ratio binaries detectable. None of the
synthetic single-PSF cases were detected as a false-positive binary at 95%
confidence.

We have performed initial binary detection and for all the data in program
GO-15648. The results of our analysis of binary probability are listed in
Table 1, and compared to the known binary KBOs in Figure 2. Of the 198 objects
that were successfully observed, we identified 23 as being visually separated
binaries, a further 19 as blended binaries with a probability of being binary
of being $>$95% (per the formulation above), and 156 that appeared as singular
objects to our analysis. Our blended binary analysis was thus able to almost
double the yield of binary KBOs in the SSOLS dataset.

Our detection of a large number of blended-PSF binary KBOs implies that binary
KBOs at apparent separations at or just below the detection limit of HST/WFC3
may be very common. This argument is further strengthened by our synthetic KBO
analysis which shows that close binary KBOs dimmer than V=24 and with flux
ratios larger than 5:1 are very hard for WFC3 to detect, even with apparent
separations greater than one pixel. A full census of the binary fraction in
the Kuiper Belt is thus very difficult from direct imaging alone. Ground-based
adaptive optics systems are even more flux limited than HST (Grundy et al.,
2011) and cannot probe binaries that are fainter or closer than those
detectable with HST/WFC3. The James Webb Space Telescope (JWST) with the
NIRCam instrument only offers a marginal increase in angular resolution over
WFC3 ($\approx$30 milliarcseconds/pixel versus $\approx$40
milliarcseconds/pixel for WFC3), but may be able to detect fainter KBO
satellites due to its larger collecting area.

The strongest evidence that very tight binary KBOs may be common comes from
two unconventional binary detection methods, a stellar occultation and direct
spacecraft observations. The 3:2 Neptune resonant KBO (523764) 2014 WC510 was
observed by stellar occultation to have a separation of 349 km (Leiva et al.,
2020), which corresponds to an apparent separation of 13 milliarcseconds, or
less than a 1/3 of a WFC3 pixel. The New Horizons spacecraft observed five
KBOs with the LORRI camera at a distance of less than 0.1 astronomical units,
and two of them appeared to be very close binaries, 2011 JY31 and 2014 OS393
(Weaver et al., 2022). 2011 JY31 was resolved as a blended binary in both
LORRI epochs, and the combination of that with the assumption of the orbital
period being equal to the lightcurve period of 46 hours allowed the mutual
orbit separation to be fit as 199 km (Weaver et al., 2022). 2014 OS393 was
only observed to be binary in one of the two epochs, and Weaver et al. (2022)
estimated its apparent separation to be $\approx$150 km. These separations
would be $<$0.2 WFC3 pixels as seen by HST, and sure enough, neither 2011 JY31
nor 2014 OS393 appears obviously binary in HST WFC3 imaging (Weaver et al.,
2022). While these unconventional detections are hard to reproduce, they do
strongly argue that binary KBOs with separations below the detection threshold
of WFC3 should be common. These tight binaries may be a result of direct
formation from the protoplanetary disk by the streaming instability (e.g.
Fraser et al., 2017; Nesvorný et al., 2019), or the result of post-formation
orbital evolution (e.g. Porter & Grundy, 2012), or some combination of those
processes. The discovery of both of these very tight binary KBOs ($<$400 km)
with unconventional methods, and the large number of 1000-3000 km binaries
that we found in the SSOLS dataset, shows that there is likely to be a
substatial number of Kuiper Belt binaries with 300-1000 km separations.

The inclinations of binary KBOs to their orbits has been shown to be very
preferentially prograde (Grundy et al., 2019). This provides an important
constraint on the formation of binary KBOs, as it strongly supports the
formation of KBO binaries by SI (Nesvorný et al., 2019). However, this test
has only been performed for relatively wide KBO binaries (Grundy et al.,
2019), as very few tight binaries have fully determined orbits. Formation
models (e.g. Nesvorný et al., 2019) show that most binary KBOs formed wider
and more eccentric than the observed binary KBOs (Grundy et al., 2019), before
the mutual orbits evolved under perturbations to tighter, more circular orbits
(Porter & Grundy, 2012). Models imply that the prograde preference should be
preserved under that mutual orbit evolution (Porter & Grundy, 2012), but that
has not yet been proven, and can only be tested by measuring the orbits of
tight binary KBOs. Future observations of the blended binary KBOs discovered
in this program could fit their mutual orbits which would be very useful to
test if this prograde preference is indeed preserved to much tighter
separations.

##  6 Future Observations and Analysis

HST/WFC3 has historically been the high-resolution instrument available for
this work (e.g. Grundy et al., 2019, and citations therein), but JWST/NIRCam
does offer slightly higher resolution than WFC3 (Rieke et al., 2023), with a
much larger collecting area. However, JWST observations of KBOs is
challenging, as the design of the spacecraft prevents it from observing the
Kuiper Belt except at the quadratures, and after the micrometorite constraint
imposed in Cycle 2, only the post-opposition quadrature window is available.
Any KBO must therefore be observed by JWST within a narrow, roughly two week
window, away from the opposition surge that pushes faint KBOs to their maximum
brightness, and with slow moving background stars that may be confused for KBO
satellites. In addition, the narrow observing window of quadrature-only space
telescopes, like JWST and the Roman Space Telescope, makes recover of the
mutual orbits of binary KBOs effectively impossible, as they cannot cover the
required temporal range (Grundy et al., 2008). The Keck laser adaptive optics
system has been used extensively for recovery of binary KBOs at slightly lower
resolution than HST (Grundy et al., 2011). Future ground-based facilities like
the Extremely Large Telescope (ELT) and the Thirty Meter Telescope (TMT) may
be able to directly resolve even tighter binaries, with their larger
collecting areas being able to compensate for the throughput losses inherent
in adaptive optics systems. In addition, both of those future large telescopes
are planned to use multi-point laser adaptive optics systems (Diolaiti et al.,
2016; Boyer et al., 2014) which would enable them to observe KBOs without the
need of a stellar appulse for tip/tilt correction (as used by Keck),
increasing their operational flexibility.

The most promising pathway to surveying the very close binary population of
the Kuiper Belt may be through stellar occultations. As Leiva et al. (2020)
has shown, stellar occultation can detect binary KBOs down to arbitrary
separations, and the use of the Gaia DR3 star catalog (Gaia Collaboration et
al., 2021) allows these occultations to be targeted with great precision
(Porter et al., 2018). The advent of low cost cameras with GPS timing
precision enables stellar occultations to be observed with small, highly
mobile telescopes, greatly increasing the coverage possible for a given
occultation event. Given that 2/6 KBOs observed up close by New Horizons
(including Arrokoth) were apparently binaries with less than 200 km
separation, it would not take many stellar occultations to find additional
very close binary KBOs.

For future analysis of this and similar datasets, we plan to enhance binary
detection with the use of Machine Learning (ML). The current results are based
on shifting Tiny-Tim PSFs and assumptions about the noise properties of WFC3
images. To use ML, we would create a large amount of synthetic single and
binary KBOs (much larger than used for our calibration here) and then train a
neural network to distinguish between them. This could result in a much more
effective detection of barely-resolvable binary KBOs, and would enable the
application of this method to historical observations of KBOs by HST to find
even more blended-PSF binaries.

##  7 Conclusions

We performed a PSF-fitting process on 198 KBOs observed with HST/WFC3 and
found that while 11.6% appeared to be visually separated binaries, a further
9.6% were blended-PSF binaries. Our results show that binary KBOs at or below
the detection threshold of HST appear to be at least as common as binary KBOs
that appear to be well separated in WFC3/UVIS images. This was confirmed by
our application of the same fitting method to synthetic data, which found that
binary KBOs fainter than V=24 or with primary/secondary flux ratios larger
than 5:1 are very unlikely to be detectable with HST. These results appear to
support the idea that the Kuiper Belt was formed through the Streaming
Instability process, and that most if not all KBOs were born as binary
systems.

This work was supported by HST program GO-15648. HST data was obtained from
the Space Telescope Science Institute, which is operated by the Association of
Universities for Research in Astronomy, Inc., under NASA contract NAS 5–26555.

## References

  * Annex et al. (2021) Annex, A., Pearson, B., Seignovert, B., et al. 2021, AndrewAnnex/SpiceyPy: SpiceyPy 4.0.1, v4.0.1, Zenodo, Zenodo, doi: 10.5281/zenodo.4883901
  * Astropy Collaboration et al. (2022) Astropy Collaboration, Price-Whelan, A. M., Lim, P. L., et al. 2022, ApJ, 935, 167, doi: 10.3847/1538-4357/ac7c74
  * Bannister et al. (2018) Bannister, M. T., Gladman, B. J., Kavelaars, J. J., et al. 2018, ApJS, 236, 18, doi: 10.3847/1538-4365/aab77a
  * Benecchi et al. (2019) Benecchi, S. D., Borncamp, D., Parker, A. H., et al. 2019, Icarus, 334, 22, doi: 10.1016/j.icarus.2019.01.025
  * Boyer et al. (2014) Boyer, C., Adkins, S., Andersen, D. R., et al. 2014, in Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series, Vol. 9148, Adaptive Optics Systems IV, ed. E. Marchetti, L. M. Close, & J.-P. Vran, 91480X, doi: 10.1117/12.2056863
  * Bradley et al. (2022) Bradley, L., Sipőcz, B., Robitaille, T., et al. 2022, astropy/photutils: 1.5.0, 1.5.0, Zenodo, Zenodo, doi: 10.5281/zenodo.6825092
  * Canup (2011) Canup, R. M. 2011, AJ, 141, 35, doi: 10.1088/0004-6256/141/2/35
  * Diolaiti et al. (2016) Diolaiti, E., Ciliegi, P., Abicca, R., et al. 2016, in Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series, Vol. 9909, Adaptive Optics Systems V, ed. E. Marchetti, L. M. Close, & J.-P. Véran, 99092D, doi: 10.1117/12.2234585
  * Elliot et al. (2005) Elliot, J. L., Kern, S. D., Clancy, K. B., et al. 2005, AJ, 129, 1117, doi: 10.1086/427395
  * Fraser et al. (2017) Fraser, W. C., Bannister, M. T., Pike, R. E., et al. 2017, Nature Astronomy, 1, 0088, doi: 10.1038/s41550-017-0088
  * Gaia Collaboration et al. (2021) Gaia Collaboration, Brown, A. G. A., Vallenari, A., et al. 2021, A&A, 649, A1, doi: 10.1051/0004-6361/202039657
  * Grundy et al. (2008) Grundy, W. M., Noll, K. S., Virtanen, J., et al. 2008, Icarus, 197, 260, doi: 10.1016/j.icarus.2008.04.004
  * Grundy et al. (2011) Grundy, W. M., Noll, K. S., Nimmo, F., et al. 2011, Icarus, 213, 678, doi: 10.1016/j.icarus.2011.03.012
  * Grundy et al. (2019) Grundy, W. M., Noll, K. S., Roe, H. G., et al. 2019, Icarus, 334, 62, doi: 10.1016/j.icarus.2019.03.035
  * Grundy et al. (2020) Grundy, W. M., Bird, M. K., Britt, D. T., et al. 2020, Science, 367, aay3705, doi: 10.1126/science.aay3705
  * Krist (1995) Krist, J. 1995, in Astronomical Society of the Pacific Conference Series, Vol. 77, Astronomical Data Analysis Software and Systems IV, ed. R. A. Shaw, H. E. Payne, & J. J. E. Hayes, 349
  * Leiva et al. (2020) Leiva, R., Buie, M. W., Keller, J. M., et al. 2020, , 1, 48, doi: 10.3847/PSJ/abb23d
  * Mao et al. (2021) Mao, X., McKinnon, W. B., Singer, K. N., et al. 2021, Journal of Geophysical Research (Planets), 126, e06961, doi: 10.1029/2021JE006961
  * McKinnon et al. (2020) McKinnon, W. B., Richardson, D. C., Marohnic, J. C., et al. 2020, Science, 367, aay6620, doi: 10.1126/science.aay6620
  * Nesvorný et al. (2019) Nesvorný, D., Li, R., Youdin, A. N., Simon, J. B., & Grundy, W. M. 2019, Nature Astronomy, 3, 808, doi: 10.1038/s41550-019-0806-z
  * Petit et al. (2011) Petit, J. M., Kavelaars, J. J., Gladman, B. J., et al. 2011, AJ, 142, 131, doi: 10.1088/0004-6256/142/4/131
  * Porter & Grundy (2012) Porter, S. B., & Grundy, W. M. 2012, Icarus, 220, 947, doi: 10.1016/j.icarus.2012.06.034
  * Porter et al. (2018) Porter, S. B., Buie, M. W., Parker, A. H., et al. 2018, AJ, 156, 20, doi: 10.3847/1538-3881/aac2e1
  * Rieke et al. (2023) Rieke, M. J., Kelly, D. M., Misselt, K., et al. 2023, PASP, 135, 028001, doi: 10.1088/1538-3873/acac53
  * Robbins et al. (2017) Robbins, S. J., Singer, K. N., Bray, V. J., et al. 2017, Icarus, 287, 187, doi: 10.1016/j.icarus.2016.09.027
  * Showalter et al. (2021) Showalter, M. R., Benecchi, S. D., Buie, M. W., et al. 2021, Icarus, 356, 114098, doi: 10.1016/j.icarus.2020.114098
  * Spencer et al. (2020) Spencer, J. R., Stern, S. A., Moore, J. M., et al. 2020, Science, 367, aay3999, doi: 10.1126/science.aay3999
  * Stern et al. (2015) Stern, S. A., Bagenal, F., Ennico, K., et al. 2015, Science, 350, aad1815, doi: 10.1126/science.aad1815
  * Thirouin & Sheppard (2018) Thirouin, A., & Sheppard, S. S. 2018, AJ, 155, 248, doi: 10.3847/1538-3881/aac0ff
  * Thirouin & Sheppard (2019) —. 2019, AJ, 157, 228, doi: 10.3847/1538-3881/ab18a9
  * Virtanen et al. (2020) Virtanen, P., Gommers, R., Oliphant, T. E., et al. 2020, Nature Methods, 17, 261, doi: 10.1038/s41592-019-0686-2
  * Weaver et al. (2022) Weaver, H. A., Porter, S. B., Spencer, J. R., & The New Horizons Science Team. 2022, , 3, 46, doi: 10.3847/PSJ/ac4cb7


# Harnessing re-programmable phase transitions to control the propagation of
sound waves

Audrey A. Watkins, Austin Eichelberg  Osama R. Bilal<EMAIL_ADDRESS>Department of Mechanical engineering, University of Connecticut, Storrs, CT
06269, USA

###### Abstract

Metamaterials can enable peculiar static and dynamic behavior (such as
negative effective mass density, dynamical stiffness, and Poisson’s ratio) due
to their geometry rather than their chemical composition. The geometry of
these metamaterials can be thought of as the phase of the material, which is
usually fixed once the material is fabricated. While there exist many
theoretical and numerical studies of metamaterials that can change phase, or
re-program, experimental realizations remain limited due to challenges in
manufacturability, the destructive nature of the re-programming and inherent
non-linearities. Through a combination of analytical, numerical and
experimental analyses, we utilize tunable, self-assembled, nonlinear magnetic
lattices to realize metamaterials with reversible phase transitions. Our
metamaterials are composed of free-floating disks, with embedded permanent
magnets, confined within magnetic boundaries. We exploit the non-destructive
nature of the adjustable magnetic boundaries to create a set of re-
programmable metamaterials to control the propagation of sound waves.
Furthermore, we demonstrate a robust, real-time tunable wave filter at ultra-
low frequencies. Our findings can expand the metamaterials horizon into
functional and tunable devices.

metamaterials, magnetic lattices, phase transitions, re-programmable,
reversible, nonlinear

Metamaterials, often composed of repeated, identical unit cells, are
structures with unique properties such as negative Poisson’s
ratios,lakes1987foam  vanishing shear modulus,milton1995elasticity ;
kadic2012practicability ; buckmann2014elasto  negative
stiffness,lakes2001extreme ; nicolaou2012mechanical  and polar
elasticity.bilal2017intrinsically  Such metamaterials could have potential
applications in wave control,bertoldi2008wave ; wang2013effects ;
shan2014harnessing ; shim2015harnessing ; raney2016stable ; bilal2017bistable
; bilal2017reprogrammable ; li2012switching ; airoldi2011design  and acoustic
cloaking.cummer2007one  The unique properties of these materials stem from
their geometric tessellations, rather than their chemical compositions. These
tessellations can be designed such that an applied external stimulus can
change the geometric arrangements of the unit cells, giving rise to new
domains, or phases, within the material.jin2020guided  Each of these phases
resembles a stable configuration of the material. The stability of the
metamaterial can be engineered such that multiple stable phases exist.
shan2015multistable ; chen2017harnessing ; jiang2019bifurcation ;
zareei2020harnessing ; schaeffer2015wave  The ability of metamaterials to
transform, or be programmed, into distinct stable phases correlates with the
presence of multiple local minima in the metamaterial energy
landscape.jin2020guided ; shan2015multistable ; benichou2013structures ;
restrepo2015phase ; frazier2017atomimetic ; zhang2019programmable ;
chen2019phase ; zhang2019energy ; yasuda2020transition ; yang2016phase  One
avenue to engineer phase transitioning metamaterials with multiple energy
minima is to capitalize on soft materials and their ability to withstand large
deformations without undergoing structural damage.li2012advances ;
li2012mechanics  This could facilitate the design of metamaterials with
reversible phase transitions.

Reversible transitions between phases, or energy minima, can be thought of as
re-programmbility of a material. Re-programmable metamaterials offer more
flexibility and adaptability in applications such as controlling
electromagnetic ma2014acoustic ; yang2014all ; bao2018design ; li2012flat ;
cui2016information ; arbabi2017planar ; lee2012reversibly  and acoustic waves.
janke2005applications  While there exist many theoretical and numerical
studies of programmable metamaterials, experimental realizations remain
limited wang2020tunable  due to challenges in
manufacturability,meloni2021engineering  the destructive nature of the re-
programming, and the inherent non-linearitylazarov2007low ;
nadkarni2014dynamics ; moleron2014solitary ; narisetti2010perturbation ;
fronk2017higher ; zaera2018propagation  within phase transitions. To overcome
the damaging nature of re-programmability and the common challenge of
manufacturing, we utilize easy-to-create, tunable, nonlinear metamaterials
capable of filtering sound waves through non-destructive, re-programmable, and
reversible phase transitions.

Figure 1:  Concept. Schematic of the metamaterials’ unit cell (opaque) used in
the experiments. The free floating disk is self-aligned in a) “plus sign” or
b) “square” orientation based on the relative position of the boundary walls.
c) Programmable transmission spectrum of the metamaterial when in phase 1
(orange) and phase 2 (blue).

In this paper, we utilize tunable, self-assembled lattices of repelling
magnets to create metamaterials that can control the propagation of sound
waves by means of phase transitions. We create our lattices by confining self-
aligned, free-floating disks with embedded permanent magnets to a fixed
boundary with a pre-engineered magnetic field. The magnetic repulsive forces
between the boundary magnets and the free-floating disks cause each disk to
self-align and rest within a magnetic potential well. Each magnetic potential
well and its respective floating disk act as a single building block, or a
unit cell, of the lattice. These identical unit cells repeat in space, forming
a one-dimensional metamaterial with unique wave propagation characteristics.
We design the geometry of the confining boundary such that the precise tuning
of the magnetic field controls the alignment of the disks. The relative
positions and orientation angles between the disks define the corresponding
metamaterial phase. By changing the boundary geometry and therefore the
resulting magnetic couplings within the unit cell, the transmittable
frequencies through the metamaterial can be programmed.

We create a unit cell that consists of a single disk bound by two sets of
adjustable boundary walls (Fig. 1). The outer set of the boundary walls
contains four fixed cylindrical magnets, each with a diameter of 3 $mm$ and
height of 3 $mm$, spaced by a distance $a$ horizontally and $b$ vertically.
The distances $a$ and $b$ define the size of the unit cell. The inner set of
the boundary walls contains four fixed cylindrical magnets, each with a
diameter of 2 $mm$ and height of 1 $mm$, spaced by the same distance $a$, both
horizontally and vertically (Fig. 1). The fixed magnets form a magnetic
potential well within the boundary, confining a floating disk that is embedded
with multiple magnets (e.g., four magnets) identical to those in the outer
boundary walls (i.e., 3 $mm$ height and diameter). Both inner and outer
boundaries can slide against each other in a parallel fashion. When the
magnets in both boundaries are vertically aligned with each other, the
resultant magnetic field causes the free-floating disks to self-align in their
unit cells with an orientation resembling a “plus sign.” This orientation is
phase 1 of our metamaterial (orange disk in figure 1 a). By sliding the set of
inner boundary walls by $\epsilon=a/2$, such that the inner boundary magnets
are halfway between the outer boundary magnets, the resultant magnetic field
causes each disk to reorient, resembling a “square shape.” This orientation is
phase 2 of our metamaterial (blue disk in figure 1 b). Each metamaterial phase
has its own wave transmission properties due to the different magnetic
couplings between the boundaries and the adjacent unit cells’ disks, allowing
for transmission by means of disk oscillation (Fig. 1 c). The boundary walls
can be adjusted an infinite number of times to reorient the floating disks
from phase 1 to phase 2 orientation and vice versa, allowing for the re-
programmability of the metamaterials.

We consider two types of disks (1) disks with two embedded magnets, and (2)
disks with four embedded magnets. The lattices composed of these disk types
are stable in two topologically distinct phases. In the case of two embedded
magnets per disk, the lowest energy states take place when the disks self-
align in an in-plane “vertical” orientation (phase 1) (Fig. 2 a(ii) inset) or
in a “slanted horizontal” orientation (phase 2) (Fig. 2 b(ii) inset). In the
case of four embedded magnets per disk, the lowest energy states exist when
the disks are in a “plus sign” orientation (phase 1) (Fig. 2 c(ii) inset) or
in a “square” orientation (phase 2) (Fig. 2 d(ii) inset). The transition
between the two low energy states (i.e., the two different phases) can be
achieved by simple sliding of the set of inner boundary walls, while the set
of outer boundary walls is fixed.

Figure 2:  Programmable transmission spectrum. Analytical dispersion curves
and experimental transmission ranges for two magnet per disk in a) phase 1 and
b) phase 2 orientation. Four magnet per disk behavior in the c) phase 1 and d)
phase 2 orientation. (a-d)(i) The analytical (red) and experimental (contour)
dispersion curves, (ii) the normalized transmission of a central disk and the
unit cell (inset), and (iii) a photo of three unit cells within the
metamaterial of the straight channel boundary in the phase 1 or phase 2
orientation (i.e., $\epsilon=0$ or $\epsilon=a/2$, respectively).

To analyze our metamaterials, we consider the unit cell as the basic building
block of each lattice. We use Bloch’s theorem bloch1929quantenmechanik  to
calculate the dispersion curves of our metamaterials by assuming that the unit
cell repeats infinitely in space. The dispersion curves represent a
correlation between frequency and wavenumber, revealing the expected
transmission regions within the frequency spectrum. The considered unit cell,
for each metamaterial phase, is composed of the free-floating disk, the fixed
boundary magnets, and the neighboring free-floating disks on each side (i.e.,
nearest neighbor interactions (Fig. A1)). The dispersion relation can be
written as an eigenvalue problem:

$[-\omega^{2}\textbf{M}+\textbf{K}(\bm{\kappa})]\bm{\phi}=0,$ (1)

where $\omega$ is the frequency, $\bm{M}=\begin{bmatrix}m_{i}&0\\\
0&m_{i}\end{bmatrix}$ is the mass matrix and $m_{i}$ is the mass of the disk
with $i$ embedded magnets, where $i=2$ or $4$ magnets. $\bm{K(\kappa)}$ is the
stiffness matrix (Eq. A1) as a function of the wavenumber, $\kappa$.
$\bm{\phi}=[u~{}v]^{T}$ is the Bloch displacement vector in the $x$ and $y$
directions, respectively.

We compute the analytical dispersion curves for a unit cell with dimensions
$a=27~{}mm$, $b=37~{}mm$, and $m_{2}=0.96$ grams (two embedded magnets per
disk) or $m_{4}=1.235$ grams (four embedded magnets per disk) (Fig. 2). The
generated dispersion curves account for both the longitudinal and shear modes,
which can also be described as the parallel polarization and perpendicular
polarization of the metamaterial, respectively. The calculated dispersion
curves for each phase are distinct from one another for both disk types.
Notably, in the phase 1 orientation for both types of disks (Fig. 2 a(i),
c(i)), neither of the transmission bands exist below a frequency of 1.8 Hz. In
the phase 2 orientation for both types of disks (Fig. 2 b(i), d(i)), the
emergent pass bands are almost grounded to a frequency of 0 Hz, indicating
lower frequency waves are transmittable through this phase. Therefore,
changing between phases for either disk type affects the lowest transmittable
frequencies. The calculated shear and longitudinal dispersion curves intersect
in all cases except for two magnets per disk in the phase 2 orientation, where
the shear and longitudinal dispersion curves are separated by a stop band
(Fig. 2 b(i)). In the case of two magnets per disk in phase 1, the calculated
transmission band spans 1.5 Hz with an upper cutoff frequency of 4 Hz (Fig. 2
a(i)). In contrast, in the case of four magnets per disk in phase 1, the
calculated transmission band spans 5.8 Hz ($\sim$ 4 times wider) with an upper
cutoff frequency of 7.6 Hz (Fig. 2 c(i)).

To experimentally validate the different phases of our metamaterials, we
realize two sets of adjustable boundary walls resembling those in figure 1.
The boundaries are cut from an acrylic glass sheet using a laser cutter (Full-
spectrum 24 pro-series). We fabricate the floating disks, with a diameter
$d=15~{}mm$, using a 3D printer (Connex 500). All embedded magnets in both the
boundary walls and the disks have the north pole facing upwards, replicating
the behavior of monopoles in plane. We place the boundary walls atop an air
bearing table (New way S1030002). We adhere a glass slide to the bottom side
of each disk to minimize friction, and place the disks within the sets of
boundary walls in random orientations. Then, we activate the air bearing table
and the resulting laminar flow beneath the disks creates a near frictionless
surface. The laminar layer of airflow beneath the disks has a fly height of 5
microns (about three orders of magnitude smaller than the longitudinal and
shear displacements of the disks), allowing us to consider the flexural motion
negligible. The disks move freely within their magnetic potential wells, self-
aligning into their lowest potential energy locations and orientations. The
observed low energy states, or phases, correspond to the analyzed lattices in
figure 2 (a-d)(iii). The re-programmability of the lattices, or the reversible
transition between these different phases, is achieved by sliding the
boundaries [see supplementary video 1].

To experimentally verify the analytically computed dispersion curves, we test
the dynamical characteristics of the different metamaterial phases once they
self-assemble (Fig. 2 (a-d)(iii)). We excite the right-most disk in each
channel longitudinally (i.e., along the direction of periodicity) with a chirp
signal between 0.2 Hz and 35 Hz using a mechanical shaker (Brüel and Kjær
4180) and a function generator (Keysight Technologies 33512B). Throughout the
duration of the chirp signal, we capture the motion of each disk using a
computer vision camera (Blackfly S USB3). We then post-process the images
using the digital image correlation software (DICe) and record the
displacement of each disk in time. The mechanical waves propagate through the
lattice in the form of an oscillatory motion of the free-floating disks. We
analyze the x-displacement (the direction of excitation) profiles of the
floating disks using 1D and 2D fast Fourier transform (2D-FFT). The 1D-FFT of
a central disk represents the transmission of the wave through the lattice at
different frequencies (Fig. 2 (a-d)(ii)). The resulting 2D-FFT is superimposed
as a contour map over the analytical dispersion curves for each metamaterial
(Fig. 2 (a-d)(i)).

The experimentally observed transmission ranges match well with the analytical
results. In addition, the transmission ranges obtained numerically through a
finite simulation using the Verlet method press2007numerical  is in agreement
with both the analytical and experimental data for phase 1 (Fig. A2). As noted
above, the increase in the number of embedded magnets in the disks causes both
the upper frequency cutoff values and the widths of transmittable frequencies
to increase. Similarly, the change in phase orientation from phase 1 to phase
2 lowers the range of transmittable frequencies. It is worth mentioning that
the shear mode is not excited in our experiments. However, we note the
presence of shear motion in the case of two magnets per disk in phase 2, due
to the coupling between the longitudinal and shear modes. The mode coupling is
a consequence of the “slanted horizontal” (as opposed to perfectly horizontal)
orientation of the disks in this phase, which we show numerically in figure
A3. We confirm the nature of the propagating modes in each metamaterial phase
by applying a pure shear excitation (i.e., perpendicular to the direction of
periodicity) (Fig. A4).

During the experiments, we observed a boundary effect on the left-most and
right-most floating disks. Upon shifting the boundary walls from phase 1 to
phase 2 and vice versa, the left-most and right-most floating disks did not
always autonomously change phase with the rest of the disks (i.e., sometimes
they remained in their initial phase orientations). These inconsistencies can
cause non-uniform re-programmability that has the potential to cause
discrepancies between the finite experimental and infinite analytical models
in more complex and large-scale systems.

Figure 3:  Circular metamaterial dynamics. Schematic of the circular channel
filled with six of the twenty-four floating disks embedded with four magnets
in each. Boundary alignments in the a) phase 1 and b) phase 2 orientation. The
displacements, scaled by a factor of 70, of each disk (green or red) and the
excitation source (blue) within the circular boundary in the phase 1
orientation when excited at c) 3.3 Hz and d) 1 Hz and e) in the phase 2
orientation when excited at 1 Hz.

To overcome the issues that could arise from the boundary effects, we design a
re-programmable metamaterial that has no left or right end: an infinite
circular channel. Similar to the inner and outer boundary walls in the
straight channel, we cut a boundary consisting of adjustable concentric
circles from acrylic glass. We embed the same 3 $mm$ x 3 $mm$ magnets in the
outer boundary walls and the 2 $mm$ x 1 $mm$ magnets in the inner boundary
walls, creating 24 magnetic potential wells (i.e., 24 unit cells). When the
floating disks are placed within the unit cells, they self-align in a manner
similar to the straight channel in either the phase 1 (Fig. 3 a) or phase 2
(Fig. 3 b) orientation, depending on the relative alignment between the
boundary walls. We note that in the phase 1 orientation the floating disks lay
against the outer circular boundary ring.

To analytically model the unit cells within the circular boundary, we utilize
the same Bloch analysis equations (Eqs. 1, A1) while accounting for the curved
boundary and the new equilibrium positions of the disks (Fig. A5). To
experimentally verify the dynamical characteristics of the different
metamaterial phases, we excite the metamaterial at the right vertical tangent
to the circular boundary. Similar to the straight channel, we post-process the
displacement of all the disks in time using the DICe software. The analytical
dispersion, experimental transmission, and experimental dispersion are plotted
in figure A6.

To demonstrate the programmability of our system in the infinitely periodic
circular channel, we excite the metamaterial in phase 1 at two different
frequencies corresponding to a pass band (at 3.3 Hz) and a stop band (at 1
Hz). Then, we induce a phase transition within our metamaterial by rotating
the outer boundary walls by $7.5^{\circ}$ (Fig. 3 a, b). As a result, all
disks change alignment from the “plus sign” orientation into the “square”
orientation (i.e., phase 2). Consequently, the dynamical characteristics of
the metamaterial change, and the frequency of 1 Hz is no longer within the
stop band, but rather is a pass band frequency. To visualize the behavior of
the metamaterials at different phases and excitation frequencies, we
superimpose the displacement profiles of all the disks on the circular channel
images (Fig. 3 c-e). The displacement of the excitation source is shown in
blue. The displacements of all other 23 freely floating disks are shown in
green (in case of pass band frequency) or red (in case of stop band
frequency). In case of phase 1 excitation at 3.3 Hz, the displacement profiles
show clear transmission of the wave through the metamaterial (Fig. 3 c). For
the same metamaterial phase excited at 1 Hz, the displacement profiles
indicate an exponential decay of the wave (Fig. 3 d). At the same excitation
frequency of 1 Hz, after the metamaterial transitions into phase 2, the wave
transmits through (Fig. 3 e).

Figure 4:  Re-programmability in time. a) Central unit cell of the straight
channel metamaterial with four magnets embedded in each floating disk
transitioning from the the (i) phase 2 to (ii) phase 1 back to (iii) phase 2
orientation. b) Displacement in time of the central disk as the metamaterial
is excited at a frequency $f=1.3$ Hz, the green highlight corresponding to the
phase 1 orientation and the red highlight corresponding to the phase 2
orientation. c) The 1D-FFT of the central disk in phase 1 (orange) and phase 2
(blue), and of the vibration source signal (black).

Inducing a phase transition in both the straight and circular channels is
reversible, re-programmable, and relatively simple. We harness these re-
programmable phase transitions as a real-time tunable wave filter. To validate
the concept of real-time tunability, we consider a metamaterial composed of
unit cells with four embedded magnets in each disk within the straight
channel. We induce cyclic phase transitions in the metamaterial (Fig. 4 a) as
we continuously excite it longitudinally at 1.3 Hz, corresponding to a pass
band in phase 2 and a stop band in phase 1 (Fig. 2 c, d). We track the
displacement of a central disk, while the lattice is being excited, through
the different phase transitions (Fig. 4 b). At time $t=0$ seconds, the
excitation is off and the metamaterial is in phase 2 (Fig. 4 a(i)). At $t=2$
seconds, the excitation source is on and the disk oscillates. At $t=11$
seconds, we start the phase transition into phase 1 (Fig. 4 a(ii)). Due to the
adjustment of the boundary walls during the phase transition, the disks
briefly register large displacements as they self-align into their phase 1
orientations. At $t=13$ seconds, the oscillatory motion of the disks begins to
decay exponentially and the disks reach their phase 1 positions and
orientations, approaching static equilibrium. At $t=33$ seconds, we start the
phase transition back into phase 2 (Fig. 4 a(iii)). The disks briefly register
large displacements during their shift back to the phase 2 orientation. At
$t=38$ seconds, the disks reach their steady-state behavior and show a
constant oscillatory motion, indicating the propagation of the sound waves
through the metamaterial. To identify the main frequency components within the
time signal, we apply 1D-FFT to the regions of attenuation (red) and
transmission (green) (Fig. 4 c). We also apply 1D-FFT to the signal of the
vibration source. The FFT shows clear transmission in the phase 2 orientation
(blue) and almost zero transmission in the phase 1 orientation (orange). The
time signal during the phase transitions demonstrates the behavior of a real-
time tunable wave filter.

We have utilized tunable, self-assembled, nonlinear magnetic lattices to
create metamaterials with reversible phase transitions. We harness these
controllable phase transitions to engineer the propagation of sound waves. Our
metamaterials are composed of free-floating disks with embedded permanent
magnets confined within magnetic boundaries. We exploit the non-destructive
nature of the adjustable magnetic boundaries to create a set of re-
programmable metamaterials within a straight channel. We further improve our
straight channel to mitigate the boundary effects by considering an infinitely
periodic circular arrangement of our metamaterials. Furthermore, we
demonstrate a robust, real-time tunable wave filter at ultra-low frequencies.
The applications for metamaterials with re-programmable phase transitions are
extensive, including wave guiding, cloaking, and vibration control. Our
findings can expand the capabilities of metamaterials used in functional and
tunable devices.

## References

  * (1) R. Lakes,  _Science_ 1987, _235_ 1038\.
  * (2) G. W. Milton, A. V. Cherkaev,  _Journal of engineering materials and technology_ 1995, _117_ , 4 483.
  * (3) M. Kadic, T. Bückmann, N. Stenger, M. Thiel, M. Wegener,  _Applied Physics Letters_ 2012, _100_ , 19 191901.
  * (4) T. Bückmann, M. Thiel, M. Kadic, R. Schittny, M. Wegener,  _Nature communications_ 2014, _5_ , 1 1.
  * (5) R. S. Lakes, T. Lee, A. Bersie, Y.-C. Wang,  _Nature_ 2001, _410_ , 6828 565.
  * (6) Z. G. Nicolaou, A. E. Motter,  _Nature materials_ 2012, _11_ , 7 608.
  * (7) O. R. Bilal, R. Süsstrunk, C. Daraio, S. D. Huber,  _Advanced Materials_ 2017, _29_ , 26.
  * (8) K. Bertoldi, M. C. Boyce,  _Physical Review B_ 2008, _78_ , 18 184107.
  * (9) P. Wang, J. Shim, K. Bertoldi,  _Physical Review B_ 2013, _88_ , 1 014304.
  * (10) S. Shan, S. H. Kang, P. Wang, C. Qu, S. Shian, E. R. Chen, K. Bertoldi,  _Advanced Functional Materials_ 2014, _24_ , 31 4935\.
  * (11) J. Shim, P. Wang, K. Bertoldi,  _International Journal of Solids and Structures_ 2015, _58_ 52\.
  * (12) J. R. Raney, N. Nadkarni, C. Daraio, D. M. Kochmann, J. A. Lewis, K. Bertoldi,  _Proceedings of the National Academy of Sciences_ 2016, _113_ , 35 9722.
  * (13) O. R. Bilal, A. Foehr, C. Daraio,  _Proceedings of the National Academy of Sciences_ 2017, _114_ , 18 4603.
  * (14) O. R. Bilal, A. Foehr, C. Daraio,  _Advanced materials_ 2017, _29_ , 39 1700628.
  * (15) J. Li, J. Shim, J. Deng, J. T. Overvelde, X. Zhu, K. Bertoldi, S. Yang,  _Soft Matter_ 2012, _8_ , 40 10322.
  * (16) L. Airoldi, M. Ruzzene,  _New Journal of Physics_ 2011, _13_ , 11 113010.
  * (17) S. A. Cummer, D. Schurig,  _New Journal of Physics_ 2007, _9_ , 3 45.
  * (18) L. Jin, R. Khajehtourian, J. Mueller, A. Rafsanjani, V. Tournat, K. Bertoldi, D. M. Kochmann,  _Proceedings of the National Academy of Sciences_ 2020, _117_ , 5 2319.
  * (19) S. Shan, S. H. Kang, J. R. Raney, P. Wang, L. Fang, F. Candido, J. A. Lewis, K. Bertoldi,  _Advanced Materials_ 2015, _27_ , 29 4296.
  * (20) T. Chen, O. R. Bilal, K. Shea, C. Daraio,  _arXiv preprint arXiv:1710.04723_ 2017.
  * (21) Y. Jiang, L. M. Korpas, J. R. Raney,  _Nature communications_ 2019, _10_ , 1 1.
  * (22) A. Zareei, B. Deng, K. Bertoldi,  _Proceedings of the National Academy of Sciences_ 2020, _117_ , 8 4015.
  * (23) M. Schaeffer, M. Ruzzene,  _Journal of Applied Physics_ 2015, _117_ , 19 194903\.
  * (24) I. Benichou, S. Givli,  _Journal of the Mechanics and Physics of Solids_ 2013, _61_ , 1 94.
  * (25) D. Restrepo, N. D. Mankame, P. D. Zavattieri,  _Extreme Mechanics Letters_ 2015, _4_ 52\.
  * (26) M. J. Frazier, D. M. Kochmann,  _Advanced Materials_ 2017, _29_ , 19 1605800.
  * (27) Y. Zhang, B. Li, Q. Zheng, G. M. Genin, C. Chen,  _Nature communications_ 2019, _10_ , 1 1.
  * (28) T. Chen, X. Zhang, X. Yan, B. Zhang, J. Jiang, S. Zhang, C. Guo, D. Huang, M. Qi,  In _Behavior and Mechanics of Multifunctional Materials XIII_ , volume 10968. International Society for Optics and Photonics, 2019 109680W.
  * (29) Y. Zhang, D. Restrepo, M. Velay-Lizancos, N. D. Mankame, P. D. Zavattieri,  _Scientific reports_ 2019, _9_ , 1 1.
  * (30) H. Yasuda, L. Korpas, J. Raney,  _Physical Review Applied_ 2020, _13_ , 5 054067.
  * (31) D. Yang, L. Jin, R. V. Martinez, K. Bertoldi, G. M. Whitesides, Z. Suo,  _Extreme Mechanics Letters_ 2016, _6_ 1\.
  * (32) S. Li, B. Sun,  _Advances in soft matter mechanics_ ,  Springer, 2012.
  * (33) B. Li, Y.-P. Cao, X.-Q. Feng, H. Gao,  _Soft Matter_ 2012, _8_ , 21 5728.
  * (34) G. Ma, M. Yang, S. Xiao, Z. Yang, P. Sheng,  _Nature materials_ 2014, _13_ , 9 873.
  * (35) Y. Yang, I. I. Kravchenko, D. P. Briggs, J. Valentine,  _Nature communications_ 2014, _5_ , 1 1.
  * (36) L. Bao, Q. Ma, G. D. Bai, H. B. Jing, R. Y. Wu, X. Fu, C. Yang, J. Wu, T. J. Cui,  _Applied Physics Letters_ 2018, _113_ , 6 063502.
  * (37) X. Li, S. Xiao, B. Cai, Q. He, T. J. Cui, L. Zhou,  _Optics letters_ 2012, _37_ , 23 4940.
  * (38) T.-J. Cui, S. Liu, L.-L. Li,  _Light: Science & Applications_ 2016, _5_ , 11 e16172.
  * (39) A. Arbabi, E. Arbabi, Y. Horie, S. M. Kamali, A. Faraon,  _Nature Photonics_ 2017, _11_ , 7 415.
  * (40) S. Lee, S. Kim, T.-T. Kim, Y. Kim, M. Choi, S. H. Lee, J.-Y. Kim, B. Min,  _Advanced Materials_ 2012, _24_ , 26 3491.
  * (41) L. Janke, C. Czaderski, M. Motavalli, J. Ruth,  _Materials and Structures_ 2005, _38_ , 5 578.
  * (42) Y.-F. Wang, Y.-Z. Wang, B. Wu, W. Chen, Y.-S. Wang,  _Applied Mechanics Reviews_ 2020, _72_ , 4.
  * (43) M. Meloni, J. Cai, Q. Zhang, D. Sang-Hoon Lee, M. Li, R. Ma, T. E. Parashkevov, J. Feng,  _Advanced Science_ 2021, 2000636.
  * (44) B. S. Lazarov, J. S. Jensen,  _International Journal of Non-Linear Mechanics_ 2007, _42_ , 10 1186.
  * (45) N. Nadkarni, C. Daraio, D. M. Kochmann,  _Physical Review E_ 2014, _90_ , 2 023204.
  * (46) M. Molerón, A. Leonard, C. Daraio,  _Journal of Applied Physics_ 2014, _115_ , 18 184901\.
  * (47) R. K. Narisetti, M. J. Leamy, M. Ruzzene,  _Journal of Vibration and Acoustics_ 2010, _132_ , 3\.
  * (48) M. D. Fronk, M. J. Leamy,  _Journal of Vibration and Acoustics_ 2017, _139_ , 5\.
  * (49) R. Zaera, J. Vila, J. Fernandez-Saez, M. Ruzzene,  _International Journal of Non-Linear Mechanics_ 2018, _106_ 188\.
  * (50) F. Bloch,  _Journal of Physics_ 1929, _52_ , 7-8 555.
  * (51) W. H. Press, S. A. Teukolsky, W. T. Vetterling, B. P. Flannery,  _Numerical recipes 3rd edition: The art of scientific computing_ ,  Cambridge university press, 2007.
  * (52) A. A. Watkins, A. Eichelberg, O. R. Bilal,  _Under review_ 2021.
  * (53) W. Jiao, S. Gonella,  _arXiv preprint arXiv:1911.10410_ 2019.
  * (54) A. A. Watkins, O. R. Bilal,  _Frontiers in Materials_ 2020, _7_ 410\.

##  Appendix A Straight Boundary Analysis

In our infinite model analysis, i.e., the dispersion curve calculations, we
consider the nearest neighbor interactions between the disks. We consider a
periodically repeating unit cell with disk index $j$. Each disk can have two
or four embedded magnets. The stiffness matrix in the utilized model (Eq. A1)
accounts for the repulsion forces between each embedded magnet ($n_{i}$),
boundary magnet ($\beta_{ii}$) in the same unit cell, and the disk magnet
($\alpha_{iii}$) in the adjacent unit cells $j+1$ and $j-1$. The stiffness
matrix, $\bm{K(\kappa)}$, consists of the static repulsive forces between the
disk magnets and the boundary magnets in the lowest potential energy locations
(Fig. A1) and can be written as watkinseichelbergbilal2021frustration :

$\mathbf{\textbf{K}}=\sum_{n=1}^{N}\Biggl{(}\sum_{\alpha=1}^{2N}\biggl{[}f_{,d}(d_{n,\alpha})_{l}\textbf{e}_{n,\alpha}\otimes\textbf{e}_{n,\alpha}\left[cos(\kappa\cdot\textbf{R}_{n,\alpha})-1\right]\\\
+\frac{f(d_{n,\alpha})_{l}}{(d_{n,\alpha})}(\textbf{I}-\textbf{e}_{n,\alpha}\otimes\textbf{e}_{n,\alpha})\left[cos(\kappa\cdot\textbf{R}_{n,\alpha})-1\right]\biggr{]}\\\
-\sum_{\beta=1}^{4}\biggl{[}f_{,d}(d_{n,\beta})_{s}\textbf{e}_{n,\beta}\otimes\textbf{e}_{n,\beta}+\frac{f(d_{n,\beta})_{s}}{(d_{n,\beta})}(\textbf{I}-\textbf{e}_{n,\beta}\otimes\textbf{e}_{n,\beta})\biggr{]}\\\
-\sum_{\beta=5}^{\Gamma}\biggl{[}f_{,d}(d_{n,\beta})_{l}\textbf{e}_{n,\beta}\otimes\textbf{e}_{n,\beta}+\frac{f(d_{n,\beta})_{l}}{(d_{n,\beta})}(\textbf{I}-\textbf{e}_{n,\beta}\otimes\textbf{e}_{n,\beta})\biggr{]}\Biggr{)}$
(A1)

where $N$ is the number of magnets embedded in a single floating disk, $n$ is
the index of the magnets in unit cell $j$, $\alpha$ is the index of the
magnets in the nearest neighboring unit cells (e.g. $j+1$, $j-1$), $\otimes$
is the dyadic product, $\bm{R_{n,\alpha}}$ is the unit vector between magnet
$n$ and $\alpha$, $\beta$ is the index of the boundary magnets, and $\Gamma$
is the number of boundary magnets. In the phase 1 orientation (i.e., when
$\epsilon=0$), $\Gamma=8$ which accounts for the four inner boundary magnets
and the four outer boundary magnets of the unit cell. In the phase 2
orientation (i.e., when $\epsilon=a/2$), $\Gamma=6$ which accounts for all
four inner boundary magnets and two outer boundary magnets (positioned halfway
between the inner boundary magnets) of the unit cell. We model the repulsive
forces measured between the embedded disk magnets and the boundary magnets as
an inverse power law $f(d)_{(.)}=A_{(.)}d^{\gamma_{(.)}}$.
$A_{(.)}=3\mu\beta^{2}/4\pi$ represents the interaction constant between the
disk magnets and the smaller inner boundary magnets ($A_{s}=2.210\times
10^{-10}$) or the larger outer boundary magnets ($A_{l}=9.617\times 10^{-11}$)
where $\mu$ is the permeability of air and $\beta$ is the magnetic moment,
$f_{,d}()$ is the first derivative of $f()$, $\gamma_{s}=-3.036$, and
$\gamma_{l}=-4$. jiao2019distinctive ; watkins2020demultiplexing

Figure A1 depicts a schematic of the unit cell vectors that appear in the
stiffness matrix in equation A1 for the straight channel boundary. Unit cells
in the phase 1 orientations (Fig. A1 a, c) have eight fixed boundary magnets
(i.e., $\Gamma=8$) while unit cells in the phase 2 orientations (Fig. A1 b, d)
have six fixed boundary magnets (i.e., $\Gamma=6$). The number of boundary
magnets within the unit cell corresponds to the $\Gamma$ value used in the
stiffness equation (Eq. A1). The index of each disk magnet is labeled. An
example unit vector between disk magnet $n_{i}$, where $i=2$, and boundary
magnet $\beta_{ii}$, where $ii=2$ and the same disk magnet $n_{2}$ and
adjacent disk magnet $\alpha_{iii}$, where $iii=4$ is shown in figure A1 a.
The coupling between disk magnets $n$ in unit cell $j$ and disk magnets
$\alpha$ in unit cells $j+1$ and $j-1$ resembles the nearest neighbor
interactions accounted for in the dispersion analysis.

Figure A1: Straight channel unit force vectors. Labeling and representation of
the unit force vectors within a unit cell, $j$ (opaque). a) Example unit
vector between disk $j$ magnets ($n$), boundary magnets ($\beta$), and the
adjacent disks’ embedded magnets $\alpha$ (transparent) used in the stiffness
matrix. Labeling and indexing of two magnets per disk b) phase 2 and four
magnets per disk c) phase 1 and d) phase 2.

To numerically validate the behavior of our metamaterials, we simulate the
frequency response function (FRF) of a finite metamaterial array composed of
20 disks in phase 1 using the Verlet method. In the case of two magnets per
disk (Fig. A2 a(ii) inset), the simulated FRF resembles the experimental
1D-FFT relatively well, accounting for a slight difference in the width of the
range of transmission between the two plots (Fig. A2 a(ii)). The upper
frequency value of the numerically simulated FRF is 4.1 Hz, whereas in the
experimentally obtained 1D-FFT, it is 4.6 Hz. The lower frequency value of the
numerically simulated FRF is 2.8 Hz, where as in the experimentally obtained
1D-FFT, it is 2.6 Hz. In the case of four magnets per disk (Fig. A2 b(ii)
inset), the numerical FRF and 1D-FFT match relatively well (Fig. A2 b). The
upper frequency value of transmission of the numerically simulated FRF is 7.4
Hz, whereas in the 1D-FFT it is 8 Hz. The lower frequency value of
transmission of the numerically simulated FRF is 2.9 Hz, whereas in the 1D-FFT
it is 3.1 Hz. Overall, the numerically simulated FRFs are a good depiction of
the behavior of the metamaterials in the phase 1 orientations when excited
longitudinally. The slight discrepancy in the numerical FRFs and experimental
ones is due to the variation in the size of the embedded magnets in the two
different boundary sets, as these magnets do not have the same thickness and
therefore are not in the same plane.

Figure A2:  Longitudinal Excitation. Numerical FRF of phase 1 metamaterial
within the straight channel boundary with a) two magnets per disk and b) four
magnets per disk. (a, b)(i) The analytical (red) and experimental 2D-FFT
(contour map) dispersion curves. (a, b)(ii) The numerical FRFs (green)
superimposed on the experimental transmission (black) through the
metamaterial. (iii) photo of the self-assembled unit cells within the phase 1
metamaterials

Within the main text of the paper, we consider the longitudinal excitation
(i.e., along the direction of periodicity) of the straight channel (Fig. 2).
This longitudinal excitation causes the displacement profiles of the disks to
correspond with the longitudinal dispersion curves only, as long as the
dispersion modes are not coupled (Fig. 2 a, c, d). When the modes are coupled
and the metamaterial is longitudinally excited, shear motion presents itself
as well (Fig. 2 b). For instance, in the case of two magnets per disk in phase
2, the “slanted horizontal” orientation of the floating disks causes the shear
and longitudinal modes to couple. To validate our observation, we numerically
simulate three different phases of the metamaterials with two embedded magnets
per disk using the Verlet method (Fig. A3). We compare the longitudinal versus
shear motion of (1) vertically aligned disks (Fig. A3 a(i)), (2) “slanted
horizontal” disks (Fig. A3 b(i)), and (3) perfectly aligned horizontal disks
(Fig. A3 c(i)). We excite each metamaterial longitudinally and plot the
simulated longitudinal and shear displacements for each configuration (Fig. A3
(a-c)(ii)). The longitudinal excitation causes shear motion to occur only when
the floating disks are in the “slanted horizontal” orientation. We note that
the perfectly horizontal alignment of the disks is not physical and such phase
is not stable for these unit cell dimensions, however, we simulate its
response to longitudinal excitation to show that the presence of the shear
mode only exists when the disks are aligned at an angle (i.e., “slanted”).

Figure A3:  Shear presence in longitudinal excitation. Schematic of the two
magnet per disk metamaterial in the a)(i) phase 1 “vertical,” b)(i) phase 2
“slanted horizontal,” and c)(i) phase 2 “perfectly horizontal” orientation.
(a-c)(ii)The numerically obtained FRF including both the longitudinal (black)
and shear (pink) modes of transmission corresponding to the disk orientation.
Figure A4:  Shear Excitation. Two magnet per disk behavior when the shear mode
is excited in the a) phase 1 and b) phase 2 orientation. Four magnet per disk
behavior when the shear mode is excited in the c) phase 1 and d) phase 2
orientation. (a-d)(i) The analytical (red) and shear experimental (contour)
dispersion curves, (ii) the normalized transmission of a central disk and the
unit cell (inset), and (iii) a photo of three unit cells within the
metamaterial of the straight channel boundary in the phase 1 or phase 2
orientation (i.e., $\epsilon=0$ or $\epsilon=a/2$, respectively).

To experimentally confirm the nature of the propagating modes through our
straight channel metamaterials (i.e., longitudinal versus shear), we excite
the self-aligned lattices (Fig. A4 (a-d)(iii)) with chirp signal perpendicular
to the direction of periodicity. The resulting motions of the floating disks
are tracked using DICe, and 1D- and 2D-FFTs are used to post-process the
y-displacements of the disks (Fig. A4). In the case of two magnets per disk in
the phase 1 orientation (Fig. A4 a), the analytical shear dispersion curve
(Fig. A4 a(i)) spans from 2.5 Hz to 2.9 Hz and the transmission (Fig. A4
a(ii)) spans from 2.2 Hz to 2.8 Hz. In the case of two magnets per disk in the
phase 2 orientation (Fig. A4 b), the analytical shear dispersion curve (Fig.
A4 b(i)) spans from 3.8 Hz to 5.1 Hz and the transmission (Fig. A4 b(ii))
spans from 2.9 Hz to 4.3 Hz. The transmission pass band that is present in
phase 2 of two magnets per disk is shifted slightly lower than the shear
dispersion curve, but the width of frequency values transmitted is very
similar, 1.3 Hz versus 1.4 Hz for the analytical shear (Fig. A4 b(i) (red))
and experimental data (Fig. A4 b(ii)), respectively. In the case of four
magnets per disk in the phase 1 orientation (Fig. A4 c), the analytical shear
dispersion curve (Fig. A4 c(i)) spans from 1.8 Hz to 3.9 Hz and the
transmission (Fig. A4 c(ii)) spans from 1.3 Hz to 3.7 Hz. In the case of four
magnets per disk in the phase 2 orientation (Fig. A4 d), the analytical shear
dispersion curve (Fig. A4 d(i)) spans from 5.5 Hz to 6.4 Hz and the
transmission (Fig. A4 d(ii)) spans from 4.5 Hz to 5.7 Hz. In both phases with
four magnets per disk, the analytical shear dispersion curves and transmission
match relatively well, noting a slight downward shift of the transmission pass
band.

Supplementary Video 1: This video shows the re-programmability and
reversibility of the straight metamaterial with four magnets embedded in each
disk. First, we align the boundary walls such that the metamaterial is in
phase 1 (each disk resembles a “plus sign”). Then, we initiate the phase
transition such that the disks reorient themselves and the metamaterial is in
phase 2 (each disk resembles a “square”). Lastly, we initiate another phase
transition and the metamaterial returns to phase 1.

Supplementary Video 2: This video shows the re-programmability and
reversibility of the metamaterial within the circular boundary with four
magnets embedded in each disk. First, we align the boundary walls such that
the metamaterial is in phase 1 (each disk resembles a “plus sign”). Then, we
initiate the phase transition such that the disks reorient themselves and the
metamaterial is in phase 2 (each disk resembles a “square”). Lastly, we
initiate another phase transition and the metamaterial returns to phase 1.

##  Appendix B Circular boundary analysis

To experimentally validate the self-alignment of the disks within the circular
boundary into two topologically distinct phases, we cut concentric circular
boundary walls out of acrylic glass sheets resembling those in figure 3 a, b.
The outer diameter of the outermost ring is 300 $mm$, while the inner diameter
is 285 $mm$. The thickness of the second largest ring is 7.5 $mm$. The inner
diameter of the innermost ring is 190 $mm$, while the outer diameter is 205
$mm$. The thickness of the second smallest ring is 7.5 $mm$. The embedded
magnets in the inner and outer sets of boundary walls resemble those in the
straight channel boundary (i.e., smaller magnets embedded in the inner set,
larger magnets embedded in the outer set) all with the north pole facing
upwards. Each unit cell is fabricated in $15^{\circ}$ increments, totalling 24
unit cells within the circle. We use the same method to fabricate and prepare
the disks that were used in the straight channel boundary. We place the
boundary on the air bearing table and the disks within the boundary walls in
random orientations. We activate the air bearing, creating a layer of laminar
flow beneath the disks. The disks self-align into their lowest potential
energy locations and orientations. Based on the relative angle between the
boundary walls, the disks align in the phase 1 or phase 2 orientations (Fig.
A6 (a-d)(iii)).

Similar to the analysis of the unit cell within the straight boundary, the
dispersion analysis for the unit cell within the circular boundary considers
the nearest neighbor interactions (Fig. A5). The labeling and indexing is the
same for the circular boundary as that for the straight boundary (Fig. A1).
The unit cells in the phase 1 orientation (Fig. A5 a, c) have eight boundary
magnets ($\Gamma=8$), while the unit cells in the phase 2 orientation (Fig. A5
b, d) have six boundary magnets ($\Gamma=6$). An example unit vector between
disk magnet $n_{i}$, where $i=2$, and boundary magnet $\beta_{ii}$, where
$ii=2$ and the same disk magnet $n_{2}$ and adjacent disk magnet
$\alpha_{iii}$, where $iii=4$ is shown in figure A5 a. While the dispersion
model is similar to that of the straight boundary, in the circular case, the
separation distances between the top and the bottom boundaries are different,
taking into account the curvature of the unit cell.

Figure A5: Circular channel force unit vectors. Labeling and representation of
the unit force vectors within a unit cell, $j$ (opaque). a) Example unit
vector between disk $j$ magnets ($n$), boundary magnets ($\beta$), and the
adjacent disks’ embedded magnets $\alpha$ (transparent) used in the stiffness
matrix. Labeling and indexing of two magnets per disk b) phase 2 and four
magnets per disk c) phase 1 and d) phase 2.

The analytically computed dispersion curves for the circular boundary are
plotted in figure A6 (a-d)(i). Similar to the straight boundary (Fig. 2), the
shift from phase 1 to phase 2 has a measurable effect on the lower frequencies
of the dispersion curves. In the phase 1 orientation, for both types of disks,
neither of the transmission bands exist below 1 Hz (Fig. A6 a(i), c(i)). In
the phase 2 orientation, for both types of disks, both transmission bands
exist at or below 0.15 Hz (Fig. A6 b(i), d(i)). This phenomenon is consistent
with that seen in the straight boundary case (i.e., the adjustment of the
boundary walls and the consequential change in locations and orientations of
the floating disks from phase 1 to phase 2 causes lower expected transmission
frequencies). In addition, increasing the number of embedded magnets in each
disk has a measurable effect on the transmission frequencies. The calculated
transmission band (i.e., the dispersion curve frequency range) in the case of
two magnets per disk in phase 1 (Fig. A6 a(i)) spans 1.1 Hz with an upper
frequency cutoff of 2.4 Hz. In contrast, the calculated transmission band the
case of four magnets per disk in phase 1 (Fig. A6 c(i)) spans 3.3 Hz (3 times
wider) with an upper frequency cutoff of 4.3 Hz.

Figure A6:  Transmission through circular channel. Analytical dispersion
curves and experimental transmission ranges for two magnet per disk in a)
phase 1 and b) phase 2 orientation. Four magnet per disk behavior in the c)
phase 1 and d) phase 2 orientation. (a-d)(i) The analytical (red) and
experimental (contour) dispersion curves, (ii) the normalized transmission of
a the tenth disk from the source and the unit cell (inset), and (iii) a photo
of three unit cells within the metamaterial of the circular channel boundary
in the phase 1 or phase 2 orientation (i.e., shifted $0^{\circ}$ or
$7.5^{\circ}$, respectively).

To experimentally verify the analytically computed dispersion curves, we
excite the metamaterials with a tangential chirp signal. Throughout the
duration of the chirp signal, we capture the motion of each floating disk
using a computer vision camera. We post-process the captured images using the
DICe software to extract the displacement of each disk over time. We analyze
the displacements of the floating disks in each of the four metamaterial
phases (characterized separately) using a 1D- and 2D-FFT. The 1D-FFT is the
measured transmission through the lattice at the tenth disk from the wave
source (Fig. A6 (a-d)(ii)). The 2D-FFT is the experimental dispersion curve,
plotted as a contour map, which we superimpose on the analytical dispersion
curve for each metamaterial phase (Fig. A6 (a-d)(i)). The individual disk
transmissions (1D-FFT) and the post-processed dispersion curves (2D-FFT) match
well. Similar to the straight boundary case, for two magnets per disk in the
phase 2 orientation, we note the emergence of shear motion, even with only
longitudinal excitation present, due to the coupling of the shear and
longitudinal modes.