# OptiTrap: Optimal Trap Trajectories for Acoustic Levitation Displays

Viktorija Paneva<EMAIL_ADDRESS>0000-0002-5152-3077 University of
BayreuthGermany ,  Arthur Fleig  University of BayreuthGermany
<EMAIL_ADDRESS>,  Diego Martínez Plasencia  University College
LondonUnited Kingdom<EMAIL_ADDRESS>,  Timm Faulwasser  TU Dortmund
UniversityGermany<EMAIL_ADDRESS>and  Jörg Müller  University
of BayreuthGermany<EMAIL_ADDRESS>
###### Abstract.

Acoustic levitation has recently demonstrated the ability to create volumetric
content by trapping and quickly moving particles along reference paths to
reveal shapes in mid-air. However, the problem of specifying physically
feasible trap trajectories to display desired shapes remains unsolved. Even if
only the final shape is of interest to the content creator, the trap
trajectories need to determine where and when the traps need to be, for the
particle to reveal the intended shape. We propose _OptiTrap_ , the first
structured numerical approach to compute trap trajectories for acoustic
levitation displays. Our approach generates trap trajectories that are
physically feasible and nearly time-optimal, and reveal generic mid-air
shapes, given only a reference path (i.e., a shape with no time information).
We provide a multi-dimensional model of the acoustic forces around a trap to
model the trap-particle system dynamics and compute optimal trap trajectories
by formulating and solving a non-linear path following problem. We formulate
our approach and evaluate it, demonstrating how _OptiTrap_ consistently
produces feasible and nearly optimal paths, with increases in size, frequency,
and accuracy of the shapes rendered, allowing us to demonstrate larger and
more complex shapes than ever shown to date.

ultrasonic levitation, minimum time problems, path following, volumetric
displays, phased arrays of transducers

††copyright: rightsretained††journal: TOG††journalyear: 2022††journalvolume:
1††journalnumber: 1††article: 1††publicationmonth: 1††price: 1††doi:
10.1145/3517746††ccs: Computing methodologies Computer graphics††ccs:
Computing methodologies Motion path planning††ccs: Computing methodologies
Computational control theory

##  1\. Introduction

Figure 1. While previously, shapes demonstrated on levitation displays were
limited to simple shapes with almost constant curvature (Fushimi et al., 2019;
Hirayama et al., 2019; Plasencia et al., 2020), our approach allows to render
generic complex paths. Shapes that have not been demonstrated before include
sharp edges as with the heart (B) and dolphin (D), as well as significant
changes of the curvature, as with the cat (A) and flower (C) (see
Supplementary Video). Photos of the physical particle are made with an
exposure time ranging from 0.2 to 1 s, such that each photo shows multiple
periods of the orbit.

Four photos labelled A, B, C and D show different graphics generated with the
levitation display. Image A shows an outline of a 2.4 centimetres wide cat,
rendered at 5 Hertz. Image B shows an outline of a 5.4 centimetres wide heart,
rendered at 10 Hertz. Image C shows an outline of a 5.6 centimetres wide
three-petal flower, rendered at 10 Hertz. Image D shows an outline of a 4.7
centimetres wide dolphin, rendered at 3.3 Hertz. A human index finger is
petting the dolphin. The volumetric dolphin graphics and the index finger are
of similar size.

Figure 2. _OptiTrap_ is an automated method to compute trap trajectories to
reveal generic mid-air shapes on levitation displays. The method accepts a
reference path (e.g., the shape of a heart) without any timing information as
an input (A). Our approach considers the capabilities of the device and its
trap-particle dynamics and combines these with a path following approach (B).
The approach produces feasible trap trajectories, describing when and where
the traps must be created (C). The resulting particle motion can be presented
on the actual device, yielding feasible paths and supporting complex objects
with sharp edges and/or significant changes of the curvature (D).

A flowchart with four elements sequentially connected with flow links. The
starting state A is a reference path q of theta, which captures the positions
of a reference shape of a heart. A flows to B, the OptiTrap algorithm,
composed of the trap-particle dynamics and optimal path following components.
The algorithm outputs a trap trajectory u of t, in state C, which contains
both the trap positions and timings. An image shows that the trap trajectories
do not overlap with the shape positions that were fed to the algorithm in
state A. Finally, the ending state D is the particle motion. A photo of the
rendered heart shape on the levitation display is shown above a stretched out
hand, indicating that the shape is three dimensional and rendered in mid-air.

Acoustic levitation has recently demonstrated the creation of volumetric
content in mid-air by exploiting the Persistence of Vision (PoV) effect. This
is achieved by using acoustic traps to rapidly move single (Hirayama et al.,
2019) or multiple (Plasencia et al., 2020) particles along a periodic
reference path, revealing a shape within $0.1s$ (i.e., the integration
interval of human eyes (Bowen et al., 1974)).

However, the way to define such PoV content remains unsolved. That is, there
are no automated approaches to compute the _trap trajectories_ , i.e., the
positioning and timing of the acoustic traps that will allow us to reveal the
desired shape.

Currently, content creators can only rely on trial and error to find
physically feasible trap trajectories resembling their intended target shapes,
resulting in a time consuming and challenging process. For instance, the
creator will need to define the timing of the path (i.e., not only _where_ the
particle must be, but also _when_). Such timing is not trivial and will affect
the overall rendering time and the accelerations applied to the particle,
which must be within the capabilities of the levitator. The way forces
distribute around the acoustic trap will also need to be considered to decide
where traps must be located.

Considering these challenges is crucial to design feasible trap trajectories,
as a single infeasible point along the trajectory typically results in the
particle being ejected from the levitator (e.g., approaching a sharp corner
too quickly). Thus, while it has been theorised that linear PoV paths of up to
4m and peak speeds of 17m/s are possible (Fushimi et al., 2020), actual
content demonstrated to date has been limited to simple shapes with almost
constant curvature along the path and much lower speeds (e.g., 0.72m/s in
(Ochiai et al., 2014); 1.2m/s in (Hirayama et al., 2019); 2.25m/s in
(Plasencia et al., 2020), combining six particles).

In this paper, we present _OptiTrap_ , the first structured numerical approach
to compute trap trajectories for acoustic levitation displays. _OptiTrap_
automates the definition of levitated PoV content, computing physically
feasible and nearly time-optimal trap trajectories given only a reference
path, i.e., the desired shape. As shown in Figure 1, this allows for larger
and more complex shapes than previously demonstrated, as well as shapes
featuring significant changes in curvature and/or sharp corners.

Our approach is summarised in Figure 2. _OptiTrap_ assumes only a generic
reference path $\textbf{q}(\theta)$ as an input, with no temporal information
(see Figure 2(A)). _OptiTrap_ formulates this as a path following problem (see
Figure 2(B)), computing the optimum timing in which a particle can traverse
such path. Our formulation considers the _Trap-Particle Dynamics_ of the
system, using a 3D model of acoustic forces around a trap. This results in a
non-invertible model, which cannot exploit differential flatness (Fliess et
al., 1995), on which most path following approaches rely. We instead provide a
coupling stage for the dynamics of our system and numerically invert the
system. This approach produces a trap trajectory $\textbf{u}(t)$ (Figure 2(C))
that results in the intended and nearly time-optimal particle motion. That is,
$\textbf{u}(t)$ defines the positions and timing of the traps that cause the
particle to reveal the target shape $\textbf{q}(\theta)$, according to the
capabilities of the actual device, cf. Figure 2(D).

In summary, we contribute the first structured numerical approach to compute
physically feasible and nearly time-optimal trap trajectories for levitation
displays, given only a reference path that the particle should follow. As a
core contribution, we provide a theoretical formulation of the problem in
terms of path following approaches, allowing optimum timing and device
properties (e.g., trap-particle dynamics) to be jointly considered. All the
details of the approach we propose are provided in Section 4.

We illustrate the potential of our approach by rendering shapes featuring
straight lines, sharp corners, and complex shapes, such as those in Figure 1.
We then provide an experimental validation of our approach, demonstrating
increases of up to 563% in the size of rendered objects, up to 150% in the
rendering frequency, and improvements in accuracy (e.g., shape revealed more
accurately). While the baseline shapes we compare against require trial and
error to determine optimal sizes or frequencies, our approach always yields
feasible paths (i.e., working in at least 9 out of 10 attempts) and makes
consistent use of accelerations very close to (but not exceeding) the maximum
achievable by the device, independently of the target shape.

These features allow _OptiTrap_ to render complex objects involving sharp
edges and significant changes in curvature that have never been demonstrated
before. Even more importantly, it provides a tool to systematically explore
the range of contents that levitation displays can create, as a key step to
exploit their potential.

##  2\. Background and Challenges

Figure 3.  Three different path timings for the cardioid shape with fixed
traversal time: (A) equidistant sampling of the arc length, (B) equidistant
sampling of the path parameter $\theta$, and (C) _OptiTrap_. The corresponding
acceleration magnitudes are depicted in (D). Strategy (A) is physically
infeasible due to infinite required acceleration at the corner ($\theta=\pi$).
Strategy (B) does not share this problem, but requires careful parametrization
of the path, while _OptiTrap_ (C) yields the timing automatically.

The figure consists of four subplots. The first three plots show a cardioid
shape in two dimentional space. The x-axis denotes the horizontal spatial
component from -4 to 4 centimetres and the y-axis the vertical from 8 to 16
centimetres, both at increments of two centimetres. Subplot D shows the
acceleration magnitude resulting from the timing strategies from A to C, on a
scale from 0 to 1000 metres per seconds squared as a function of the path
parameter theta, going from 0 to two pi. The most notable difference in the
timing strategies occurs at the sharp corner of the cardioid shape, where the
equidistant arc length strategy results in infinite acceleration.

Our goal is is to minimise complexity for content creators. Thus, given a
reference path (i.e., a shape defined by the content creator), our approach
must produce physically feasible trap trajectories, which accurately reveal
the desired path, while making optimal use of the capabilities of the device.

We formalise our content as a reference path $Q$, provided by the content
creator without any preassigned time information. It is an explicitly
parametrized curve

(1)
$Q:=\\{\xi\in\mathbb{R}^{3}\mid\theta\in[\theta_{0},\theta_{f}]\mapsto\textbf{q}(\theta)\\},$

where $\textbf{q}\colon\mathbb{R}\rightarrow\mathbb{R}^{3}$, as later
justified in Section 4.2, must be a twice continuously differentiable function
with respect to $\theta$, which is called the path parameter. Increasing
values of $\theta$ denote forward movement along the path $Q$. The starting
point on the path is $\textbf{q}(\theta_{0})$, while $\textbf{q}(\theta_{f})$
marks the end of the path. For periodic paths, such as a particle cyclically
revealing a PoV shape, $\textbf{q}(\theta_{0})=\textbf{q}(\theta_{f})$ and
$\dot{\textbf{q}}(\theta_{0})=\dot{\textbf{q}}(\theta_{f})$ hold. Please note
that while splines would provide a generic solution to this parametrization
(i.e., twice continuously differentiable fitting the points in $Q$), any other
parametric functions can be used. For example, taking $\theta\in[0,2\pi]$ and
$r>0$, a cardioid as in Figure 3, can be described by:

(2)
$\textbf{q}(\theta)=(0,r\sin(\theta)(1+\cos(\theta)),-r\cos(\theta)(1+\cos(\theta))+r)^{\top}.$

Whatever the approach used, revealing such reference paths typically involves
rapid particle movements, as PoV content must be revealed within about $0.1s$
(Bowen et al., 1974)) and must consider the feasibility of the path (i.e., the
trap-particle dynamics, the topology of the trap, and the capabilities of the
levitator). This results in the following two main challenges for _OptiTrap_ :
1) determining the optimal timing for the particle revealing the path; and 2)
computing the trap positions generating the required forces on the particle.

###  2.1. Challenge 1: Determining an Optimal Path Timing

The reference path $Q$ defines the geometry but not the timing (i.e., it
describes _where_ the particle needs to be, but not _when_ it needs to be
there). Our approach must compute such a timing, and the strategy followed
will have important implications on the velocity and accelerations applied to
the particle and, hence, on the physical feasibility of the content.

Different timing strategies and their effects on the feasibility can be
illustrated using the cardioid example from (2). Figure 3 depicts various
timing strategies applied to this shape, all of them traversing the same path
in the same overall time.

A straight-forward approach would be to sample equidistantly along the path,
as shown in Figure 3(A). This works well for lines and circles, where a
constant speed can be maintained, but fails at sudden changes in curvature due
to uncontrolled accelerations (see the acceleration at the corner of cardioid
in Figure 3(D), at $\theta=\pi$).

The second example in Figure 3(B) shows a simple equidistant sampling on the
path parameter $\theta$. As shown in Figure 3(D), this results in areas of
strong curvature naturally involving low accelerations, and can work
particularly well for low-frequency path parametrizations in terms of
sinusoidals (e.g., (Fushimi et al., 2020) showed such sinusoidal timing along
straight paths, to theorise optimum/optimistic content sizes). In any case,
the quality of the result obtained by this approach will vary depending on the
specific shape and parametrization used.

As an alternative, the third example shows the timing produced by our
approach, based on the shape and capabilities of the device (i.e., forces it
can produce in each direction). Please note how this timing reduces the
maximum accelerations required (i.e., retains them below the limits of the
device), by allowing the particle to travel faster in parts that were
unnecessarily slow (e.g., Figure 3(D), at $\theta=\pi$). For the same maximum
accelerations (i.e., the same device), this timing strategy could produce
larger shapes or render them in shorter times, for better refresh rates. This
illustrates the impact of a careful timing strategy when revealing any given
reference path. We address this challenge in Section 4.2.

###  2.2. Challenge 2: Computation of Trap Positions

In addition to the timing, the approach must also compute the location of the
traps. Previous approaches (Hirayama et al., 2019; Plasencia et al., 2020;
Fushimi et al., 2020) placed the traps along the shape to be presented, under
the implied assumption that the particle would remain in the centre of such
trap.

However, the location where the traps must be created almost never matches the
location of the particle. Acoustic traps feature (almost) null forces at the
centre of the trap and high restorative forces around them (Marzo et al.,
2015). As such, the only way a trap can accelerate a particle is by having it
placed at a distance from the centre of the trap. Such trap-particle distances
were measured by (Hirayama et al., 2019) to assess performance achieved during
their speed tests. Distortions related to rendering fast moving PoV shapes
were also shown in (Fushimi et al., 2019). However, none of them considered or
corrected for such displacements.

Please note the specific displacement between the trap and the particle will
depend on several factors, such as the particle acceleration required at each
point along the shape, as well as the trap topology (i.e., how forces
distribute around it). No assumptions can be made that the traps will remain
constrained to positions along or tangential to the target shape. We describe
how our approach solves this challenge in Section 4.3.

##  3\. Related Work

_OptiTrap_ addresses the challenges above by drawing from advances in the
fields of acoustic levitation and control theory, which we review in this
section.

###  3.1. Acoustic Levitation

Single frequency sound-waves were first observed to trap dust particles in the
lobes of a standing wave more than 150 years ago (Stevens, 1899). This has
been used to create mid-air displays with particles acting as 3D voxels
(Ochiai et al., 2014; Omirou et al., 2015, 2016; Sahoo et al., 2016), but such
standing waves do not allow control of individual particles. Other approaches
have included Bessel beams (Norasikin et al., 2019), self-bending beams
(Norasikin et al., 2018), boundary holograms (Inoue et al., 2019) or near-
field levitation (Andrade et al., 2016).

However, most display approaches have relied on the generic levitation
framework proposed by (Marzo et al., 2015), combining a focus pattern and a
levitation signature. Although several trap topologies (i.e., twin traps,
vortex traps or bottle beams) and layouts (i.e., one-sided, two-sided,
v-shape) are possible, most displays proposed have adopted a top-bottom
levitation setup and twin traps, as these result in highest vertical trapping
forces. This has allowed individually controllable particles/voxels (Marzo and
Drinkwater, 2019) or even particles attached to other props and projection
surfaces for richer types of content (Morales et al., 2019; Fender et al.,
2021). The use of single (Fushimi et al., 2019; Hirayama et al., 2019) or
multiple (Plasencia et al., 2020) fast moving particles have allowed for
dynamic and free-form volumetric content, but is still limited to small sizes
and simple vector graphics (Fushimi et al., 2020).

Several practical aspects have been explored around such displays, such as
selection (Freeman et al., 2018) and manipulation techniques (Bachynskyi et
al., 2018), content detection and initialisation (Fender et al., 2021) or
collision avoidance (Reynal et al., 2020).

However, no efforts have been made towards optimising content considering the
capabilities (i.e., the dynamics) of such displays, particularly for
challenging content such as the one created by PoV high-speed particles.
(Hirayama et al., 2019) showed sound-fields must be updated at very high rates
for the trap-particle system to engage in high accelerations, identifying
optimum control for rates above 10kHz. However, they only provided a few
guidelines (e.g., maximum speed in corners, maximum horizontal/vertical
accelerations) to guide the definition of the PoV content. (Fushimi et al.,
2020) provided a theoretical exploration of this topic, looking at maximum
achievable speeds and content sizes, according to the particle sizes and sound
frequency used. While the dynamics of the system were considered, these were
extremely simplified, using a model of acoustic trap forces and system
dynamics that are only applicable for oscillating recti-linear trajectories
along the vertical axis of the levitator. (Paneva et al., 2020) proposed a
generic system simulating the dynamics of a top-bottom levitation system. This
operates as a forward model simulating the behaviour of the particle given a
specific path for the traps, but unlike our approach it does not address the
inverse problem. Thus, _OptiTrap_ is the first algorithm allowing the
definition of PoV content of generic shapes, starting only from a geometric
definition (i.e., shape to present, no timing information) and optimising it
according to the capabilities of the device and the dynamics of the trap-
particle system.

###  3.2. Path Following and Optimal Control

The particle in the trap constitutes a dynamical system, which can be
controlled by setting the location of the trap. As such, making the particle
traverse the reference path can be seen as an optimal control problem (OCP).

In the context of computer graphics, optimal control approaches have been
studied particularly in the areas of physics-based character animation
(Geijtenbeek and Pronost, 2012) and aerial videography (Nägeli et al., 2017).
Such OCPs can be considered as function-space variants of nonlinear programs,
whereby nonlinear dynamics are considered as equality constraints. In
engineering applications, OCPs are frequently solved via direct discretization
(Von Stryk, 1993; Bock and Plitt, 1984), which leads to finite-dimensional
nonlinear programs. In physics-based character animation, such direct solution
methods are known as spacetime constraints (Witkin and Kass, 1988; Rose et
al., 1996).

If framed as _a levitated particle along a given geometric reference path
(i.e., the PoV shape)_ , the problem can be seen as an instance of a path
following problem (Faulwasser, 2012). Such problems also occur in aerial
videography, when a drone is to fly along a reference path (Nägeli et al.,
2017; Roberts and Hanrahan, 2016). In order to yield a physically feasible
trajectory, it is necessary to adjust the timing along that trajectory. This
can be done by computing feed-forward input signals (Faulwasser, 2012; Roberts
and Hanrahan, 2016) or, if the system dynamics and computational resources
allow, using closed-loop solutions, such as Model Predictive Control in
(Nägeli et al., 2017).

However, optimal control of acoustically levitated particles cannot be
approached using such techniques. All of the above approaches require
differential flatness (Fliess et al., 1995). That is, the underlying system
dynamics must be invertible, allowing for the problem to be solved by
projecting the dynamics of the moving object onto the path manifold (Nielsen
and Maggiore, 2008). This inversion approach also enables formulating general
path following problems in the language of optimal control, encoding desired
objectives (e.g., minimum time, see (Shin and McKay, 1985; Verscheure, 2009;
Faulwasser et al., 2014)).

As discussed later in Section 4.1, the distribution of forces around the
acoustic trap (i.e., our system dynamics) are not invertible, requiring
approaches that have not been widely developed in the literature. We do this
by coupling the non-invertible particle dynamics into a virtual system,
solvable with a conventional path following approach. We then use these
coupling parameters and our model of particle dynamics to solve for the
location of the traps.

##  4\. Optimal Control for Levitation Displays

This section provides a description of our _OptiTrap_ approach, which
automates the definition of levitated content, computing physically feasible
and nearly time-optimal trap trajectories using only a reference path as an
input. In Section 4.1, we first describe our specific hardware setup and the
general mathematical framework, then we consider existing models of the trap-
particle dynamics (Subsection 4.1.1), before introducing our proposed model
(Subsection 4.1.2). Next, we describe the two stages in our algorithm, which
match the challenges identified in Section 2. That is, Section 4.2 describes
the computation of optimum timing, approached as an optimal path following
problem, while Section 4.3 explains how to compute the trap locations. Please
note that Section 4 focuses on the general case of presenting levitated
content, that is, particles cyclically traversing a reference path (shape) as
to reveal it. Other cases, such as a particle accelerating from rest as to
reach the initial state (i.e., initial position and speed) required to render
the content can be easily derived from the general case presented here, and
are detailed in the Supplementary Material S1.

Figure 4. Overview of the components in our setup. We used two opposed arrays
of transducers at a distance of 23.9 cm and an OptiTrack system to track the
position of levitated particles in real time.

A photo of equipment. Two square array boards of transducers with dimension
16.8 centimetres are aligned to exactly face each other at a vertical distance
of 23.9 centimetres. Two optical motion cameras mounted on the side, are faced
towards the volume between the boards.

###  4.1. Modelling the Trap-Particle Dynamics

We start by describing the model of the trap-particle dynamics used for our
specific setup, shown in Figure 4. This setup uses two opposed arrays of 16×16
transducers controlled by an FPGA and an OptiTrack tracking system (Prime 13
motion capture system at a frequency of 240Hz). The design of the arrays is a
reproduction of the setup in (Morales et al., 2021), modified to operate at
20Vpp and higher update rates of 10kHz. The device generates a single twin-
trap using the method described in (Hirayama et al., 2019), allowing for
vertical and horizontal forces of $4.2\cdot 10^{-5}$N and $2.1\cdot 10^{-5}$N,
respectively, experimentally computed using the linear speed tests in
(Hirayama et al., 2019) (i.e., 10cm paths, binary search with 9 out 10 success
ratios, with particle mass $m\approx 0.7\cdot 10^{-7}$ kg). Note the
substantial difference between the maximum forces in the vertical and
horizontal directions. Our approach will need to remain aware of direction as
this determines maximum accelerations.

In general, the trap-particle dynamics of such system can be described by
simple Newtonian mechanics, i.e.,

(3)
$m\ddot{\textbf{p}}(t)=F(\textbf{p}(t),\dot{\textbf{p}}(t),\textbf{u}(t)).$

Here, $\textbf{p}(t)=(p_{x},p_{y},p_{z})^{\top}\in\mathbb{R}^{3}$ represents
the particle position in Cartesian $(x,y,z)$ coordinates at time
$t\in\mathbb{R}_{0}^{+}$, and $\dot{\textbf{p}}(t)$ and $\ddot{\textbf{p}}(t)$
are the velocity and acceleration of the particle, respectively. The force
acting on the particle is mostly driven by the acoustic radiation forces and,
as such, drag and gravitational forces can be neglected (Hirayama et al.,
2019). Therefore, the net force acting on the particle will depend only on
$\textbf{p}(t)$ and on the position of the acoustic trap at time $t$ denoted
by $\textbf{u}(t)=(u_{x}(t),u_{y}(t),u_{z}(t))^{\top}\in\mathbb{R}^{3}$.

As a result, our approach requires an accurate and ideally invertible model of
the acoustic forces delivered by an acoustic trap. That is, we seek an
accurate model predicting forces at any point around a trap only in terms of
$\textbf{u}(t)$ and $\textbf{p}(t)$. Moreover, an invertible model would allow
us to analytically determine where to place the trap as to produce a specific
force on the particle given the particle location.

For spherical particles considerably smaller than the acoustic wavelength and
operating in the far-field regime, such as those used by our device, the
acoustic forces exerted can be modelled by the gradient of the Gor’kov
potential (Bruus, 2012). This is a generic model, suitable to model acoustic
forces resulting from any combination of transducer locations and transducer
activations, but it also depends on all these parameters, making it inadequate
for our approach.

Our case, using a top-bottom setup and vertical twin traps is much more
specific. This allows for simplified analytical models with forces depending
only on the relative position of the trap and the particle, which we compare
to forces as predicted from the Gor’kov potential in Figure 5.

####  4.1.1. Existing Models of the Trap-Particle Dynamics

Figure 5. Analytical models of horizontal (A) and vertical (B) acoustic forces
around a trap. The plots on the left show forces along the main axes X and Z,
and analytical models provide a good fit. The plots on the right show
horizontal (A) and vertical (B) forces across 2D slices through the trap
centre along the X and Z axes. Spring and Sinusoidal models fail to predict
forces outside the main axes, while our proposed model provides accurate
reconstruction within the region of interest (highlight).

Two subplots A and B showing the acoustic force on the y-axis going from F
maximum to F minimum, as a function of the displacement between the acoustic
trap and the particle, ranging from -0.01 to 0.01 metres.

Simple _spring_ models have been used extensively (Paneva et al., 2020;
Fushimi et al., 2018, 2020), modelling trapping forces according to a
_stiffness_ parameter $\mathcal{K}_{i}$, that is, with forces being
proportional to the distance of the particle to the centre of the trap:

$F_{i}(\textbf{p},\textbf{u}):=\mathcal{K}_{i}\cdot|u_{i}-p_{i}|,\quad
i\in\\{x,y,z\\}.$

Such models are usually refined by providing specific stiffness values for
each dimension, but they are only suitable for particles remaining in close
proximity to the centre of the trap, i.e., the linear region near the centre
of the trap. As a result, _spring_ models are only accurate for systems moving
particles slowly, requiring low acceleration and forces (so that particles
remain within the linear region).

_Sinusoidal_ models have been proposed as an alternative (Fushimi et al.,
2018, 2020), providing accurate fitting from the centre of the trap to the
peaks of the force distribution in Figure 5, according to the peak trapping
force $\mathcal{A}_{i}$ and characteristic frequency $\mathcal{V}_{i}$:

$F_{i}(\textbf{p},\textbf{u}):=\mathcal{A}_{i}\cdot\sin(\mathcal{V}_{i}\cdot(u_{i}-p_{i})),\quad
i\in\\{x,y,z\\},$

However, both of these models (i.e., _spring_ and _sinusoidal_) are only
suitable for particles placed along one of the main axes of the acoustic trap,
not for particles arbitrarily placed at any point around it. This is
illustrated on the right of Figure 5, which provides an overview of how forces
distribute on a horizontal and vertical 2D plane around a trap, according to
each model (i.e., Gor’kov, _spring_ , _sinusoidal_ , and _Ours_ , detailed in
the next subsection).

Please note how the three previous models show good matching along the
horizontal and vertical axes (represented as blue and red lines). However, the
_spring_ and _sinusoidal_ models are of one-dimensional nature (e.g., force
$F_{x}$ only depends on X distance $(u_{x}-p_{x})$) and become inaccurate at
points deviating from the main axes.

####  4.1.2. Our Model

The complexity of the force distribution modelled by the Gor’kov potential
increases as the distance to the centre of the trap increases. However, we
only need to derive a model allowing us to predict where to place the trap to
produce a specific force on the particle. This allows us to limit our
considerations to the region corresponding to the peaks designated by
$\mathcal{A}_{r}$ and $\mathcal{A}_{z}$ in Figure 5.

For points within this region, forces around twin traps distribute in a mostly
_axis-symmetric_ fashion, which can be approximated as:

(4a) $\displaystyle F_{r}(\textbf{p},\textbf{u}):=\ $
$\displaystyle\mathcal{A}_{r}\cdot\cos\left(\mathcal{V}_{z}\cdot{(u_{z}-p_{z})}\right)\cdot$
$\displaystyle\sin\left(\mathcal{V}_{xr}\cdot{\sqrt{(u_{x}-p_{x})^{2}+(u_{y}-p_{y})^{2}}}\right),$
(4b) $\displaystyle F_{z}(\textbf{p},\textbf{u}):=\ $
$\displaystyle\mathcal{A}_{z}\cdot\sin\left(\mathcal{V}_{z}\cdot(u_{z}-p_{z})\right)\cdot$
$\displaystyle\cos\left(\mathcal{V}_{zr}\cdot\sqrt{(u_{x}-p_{x})^{2}+(u_{y}-p_{y})^{2}}\right),$

where $\mathcal{A}_{r},\mathcal{A}_{z}$ represent peak trapping forces along
the radial and vertical directions of the trap, respectively.
$\mathcal{V}_{z}$, $\mathcal{V}_{xr}$, $\mathcal{V}_{zr}$ represent
characteristic frequencies of the sinusoidals describing how the forces evolve
around the trap. The resulting forces can be converted into acoustic forces in
3D space from these cylindrical coordinates, with azimuth
$\phi=\arctan((u_{y}-p_{y})/(u_{x}-p_{x}))$:

(5) $F(\textbf{p},\textbf{u})=\begin{pmatrix}F_{x}(\textbf{p},\textbf{u})\\\
F_{y}(\textbf{p},\textbf{u})\\\
F_{z}(\textbf{p},\textbf{u})\end{pmatrix}:=\begin{pmatrix}F_{r}(\textbf{p},\textbf{u})\cos\phi\\\
F_{r}(\textbf{p},\textbf{u})\sin\phi\\\
F_{z}(\textbf{p},\textbf{u})\end{pmatrix}.$

We validated our model by comparing its accuracy against the forces predicted
by the gradient of the Gor’kov potential. More specifically, we simulated 729
single traps homogeneously distributed across the working volume of our
levitator (i.e., 8x8x8cm, in line with (Hirayama et al., 2019; Fushimi et al.,
2019)), testing 400 points around each trap for a total of 400x729 force
estimations.

Table 1. Fit parameters for the Spring, Sinusoidal and Axis-symmetric models and the respective average relative errors when compared to Gor’kov. Model | Spring | Sinusoidal | Axis-symmetric
---|---|---|---
Fit | $\mathcal{K}_{x}$ | -0.0071 | $\mathcal{A}_{x}$ | 0.00009 | $\mathcal{A}_{r}$ | 0.0004636
$\mathcal{K}_{y}$ | -0.0071 | $\mathcal{A}_{y}$ | 0.00009 | $\mathcal{A}_{z}$ | 0.0002758
$\mathcal{K}_{z}$ | -0.94 | $\mathcal{A}_{z}$ | -0.0019 | $\mathcal{V}_{z}$ | 1307.83
|  | $\mathcal{V}_{x}$ | -68.92 | $\mathcal{V}_{xr}$ | -476.49
|  | $\mathcal{V}_{y}$ | -68.92 | $\mathcal{V}_{zr}$ | 287.87
|  | $\mathcal{V}_{z}$ | 1307.83 |  |
Error | 63.3% | 35.9% | 4.3%

Table 1 summarises the error distribution achieved by these three models when
compared to Gor’kov, showing an average relative error as low as 4% for our
model and much poorer fitting for the other two models. Full details and raw
data used in this validation can be found in the Supplementary Material S2.

As a final summary, this results in a model for our trap-particle dynamics
that only depends on p and u and which fits the definition of a second-order
ordinary differential equation:

(6) $m\ddot{\textbf{p}}(t)=F(\textbf{p}(t),\textbf{u}(t)).$

While the resulting model is accurate (4% relative error), we note that it is
not invertible, which will complicate the formulation of our solution
approach. For instance, for a particle at $\textbf{p}=(0,0,0)$, force
$F(\textbf{p},\textbf{u})=(0,0,\mathcal{A}_{z}\cdot\cos\left(\mathcal{V}_{zr}\cdot
x\right))$ can be obtained with either
$\textbf{u}=(x,0,\pi/(2\mathcal{V}_{z}))$ or
$\textbf{u}=(-x,0,\pi/(2\mathcal{V}_{z}))$.

###  4.2. Open-Loop Optimal Path Following

This section computes the timing for the particle, so that it moves along the
given reference path $Q$ from (1) in minimum time and according to the
dynamics of the system. We approach this as an open-loop (or feed-forward)
path following problem, as typically done in robotics (Faulwasser, 2012). That
is, we design an optimal control problem that computes the timing
$t\mapsto\theta(t)$ as to keep the levitated particle on the prescribed path,
to traverse the path in optimum (minimum) time while considering the trap-
particle dynamics (i.e., minimum _feasible_ time).

Our non-invertible model of particle dynamics calls for a more complex
treatment of the problem, which we split in two parts. In the first part we
derive a virtual system for the timing law, which we describe as a system of
first-order differential equations, solvable with traditional methods. In the
second part we couple the timing law with our non-invertible dynamics,
introducing auxiliary variables that will later enable pseudo-inversion (i.e.,
to compute trap placement for a given force) as described in Section 4.3.

####  4.2.1. Error Dynamics and the Timing Law.

The requirement that the particle follows the path $Q$ exactly and for all
times means that the deviation from the path equals zero for all
$t\in\mathbb{R}^{+}_{0}$, i.e.,

$\textbf{e}(t):=\textbf{p}(t)-\textbf{q}(\theta(t))\equiv 0.$

If the path deviation $\textbf{e}(t)$ is $0$ during the whole interval
$[t_{0},t_{1}]$, this implies that the time derivatives of $\textbf{e}(t)$
also have to vanish on $(t_{0},t_{1})$. Considering this (i.e.,
$\dot{e}(t)\equiv 0$ and $\ddot{e}(t)\equiv 0$) results in:

(7a) $\displaystyle\textbf{p}(t)$ $\displaystyle=\textbf{q}(\theta(t)),$ (7b)
$\displaystyle\dot{\textbf{p}}(t)$
$\displaystyle=\dot{\textbf{q}}(\theta(t))=\frac{\partial\textbf{q}}{\partial\theta}\dot{\theta}(t),$
(7c) $\displaystyle\ddot{\textbf{p}}(t)$
$\displaystyle=\ddot{\textbf{q}}(\theta(t))=\frac{\partial^{2}\textbf{q}}{\partial\theta^{2}}\dot{\theta}(t)^{2}+\frac{\partial\textbf{q}}{\partial\theta}\ddot{\theta}(t).$

Thus, for particles on the path $Q$, system dynamics of the form (6), and
provided we are able to express u as a function of p and $\ddot{\textbf{p}}$
(i.e., the system inversion described in Section 3.2), the position, the
velocity, and the acceleration of the particle can be expressed via
$\theta,\dot{\theta},\ddot{\theta}$, respectively. For details, a formal
derivation, and tutorial introductions we refer to (Faulwasser, 2012;
Faulwasser et al., 2017).

Observe that in (7) the partial derivatives
$\frac{\partial^{2}\textbf{q}}{\partial\theta^{2}}$ and
$\frac{\partial\textbf{q}}{\partial\theta}$, as well as $\dot{\theta}$ and
$\ddot{\theta}$ appear. This leads to two additional observations: First, the
parametrisation q from (1) should be at least twice continuously
differentiable with respect to $\theta$, so that the partial derivatives are
well-defined. This justifies our constraint in Section 2. Please note this
does not rule out corners in the path $Q$, as shown by the parametrisation (2)
of the cardioid. Second, the time evolution of $\theta$ should be continuously
differentiable, as otherwise large jumps can occur in the acceleration.

To avoid said jumps, we generate the timing $t\mapsto\theta(t)$ via the double
integrator

(8) $\ddot{\theta}(t)=v(t),$

where $v(t)\in\mathbb{R}$ is a computational degree of freedom, used to
control the progress of the particle along $Q$. The function $v(t)$ enables us
to later cast the computation of time-optimal motions along reference paths as
an OCP, see (16).

The in-homogeneous second-order ordinary differential equation (8) takes care
of jumps, but needs to be augmented by conditions on $\theta$ and
$\dot{\theta}$ at initial time $t=0$ and final time $t=T$. This is done to
represent the periodic nature of our content (i.e., the particle reveals the
same path many times per second):

(9)
$\theta(0)=\theta_{0},\quad\theta(T)=\theta_{f},\quad\dot{\theta}(0)=\dot{\theta}(T),$

where $\theta_{0}$ and $\theta_{f}$ are taken from (1) and the total time $T$
will be optimally determined by the OCP. The first two equations in (9) ensure
that the path $Q$ is fully traversed, while the third one ensures the speeds
at the beginning and end of the path match.

With these additional considerations, we define the (state) vector
$\textbf{z}(t):=(\theta(t),\dot{\theta}(t))^{\top}$, which finally allows us
rewrite our second-order differential equations in (8) as a system of first-
order differential equations:

(10) $\dot{\textbf{z}}(t)=\begin{pmatrix}0&1\\\
0&0\end{pmatrix}\textbf{z}(t)+\begin{pmatrix}0\\\
1\end{pmatrix}v(t),\quad\textbf{z}(0)=\textbf{z}_{0},\leavevmode\nobreak\
\textbf{z}(T)=\textbf{z}_{T}.$

Please note that (10) is an equivalent _virtual_ system to (8) and (9),
solvable using standard Runge-Kutta methods (Butcher, 2016). The system (10)
is suitable to generate the timing law, but the particle dynamics (6) still
need to be included, as described next.

####  4.2.2. Coupling of non-invertible particle dynamics:

To include the particle dynamics (6) in the computation of the timing, we need
to couple them with the system (10). The typical approach is to rewrite (6)
as:

(11)
$M(\ddot{\textbf{p}}(t),{\textbf{p}}(t),{\textbf{u}}(t)):=m\ddot{\textbf{p}}(t)-F(\textbf{p}(t),\textbf{u}(t))=0$

and make sure that this equation locally admits an inverse function:

(12) ${\textbf{u}}(t)=M^{-1}(\ddot{\textbf{p}}(t),{\textbf{p}}(t)),$

This would allow us to easily compute the location of our traps. However, such
inversion is not straightforward for our model, and to the best of our
knowledge, standard solution methods do not exist for such cases.

We deal with this challenge by introducing constraints related to our particle
dynamics. These will allow us to compute feasible timings while delaying the
computation of the exact trap location to a later stage in the process.

More specifically, we introduce auxiliary variables $\zeta_{1},...,\zeta_{6}$
for each trigonometric term in the force $F$, and we replace
$p_{i}=q_{i}(\theta)$ for $i\in\\{x,y,z\\}$ (i.e., particle position exactly
matches our reference path):

(13a) $\displaystyle\zeta_{1}$
$\displaystyle=\sin\left(\mathcal{V}_{xr}{\sqrt{(u_{x}-q_{x}(\theta))^{2}+(u_{y}-q_{y}(\theta))^{2}}}\right),$
(13b) $\displaystyle\zeta_{2}$
$\displaystyle=\cos\left(\mathcal{V}_{z}\cdot{(u_{z}-q_{z}(\theta)}\right),$
(13c) $\displaystyle\zeta_{3}$
$\displaystyle=\sin\left(\mathcal{V}_{z}\cdot{(u_{z}-q_{z}(\theta)}\right),$
(13d) $\displaystyle\zeta_{4}$
$\displaystyle=\cos\left(\mathcal{V}_{zr}{\sqrt{(u_{x}-q_{x}(\theta))^{2}+(u_{y}-q_{y}(\theta))^{2}}}\right),$
(13e) $\displaystyle\zeta_{5}$ $\displaystyle=\sin\phi,$ (13f)
$\displaystyle\zeta_{6}$ $\displaystyle=\cos\phi.$

With this, we are now able to formally express the force (5) in terms of
$\zeta:=(\zeta_{1},...,\zeta_{6})$, i.e.,

(14)
$\tilde{F}(\zeta):=\begin{pmatrix}\mathcal{A}_{r}\zeta_{1}\zeta_{2}\zeta_{6}\\\
\mathcal{A}_{r}\zeta_{1}\zeta_{2}\zeta_{5}\\\
\mathcal{A}_{z}\zeta_{4}\zeta_{3}\end{pmatrix}.$

Using these auxiliary variables, we can now couple the trap-particle dynamics
(6) with the virtual system (10). To this end, similar to (11), we define the
following constraint along path $Q$:

(15)
$\widetilde{M}(\theta(t),\dot{\theta}(t),v(t),\zeta(t)):=m\ddot{\textbf{q}}(\theta(t))-\tilde{F}({\zeta}(t))=0.$

Observe that i) we need $\dot{\theta}(t)$ due to (7c); and ii)
$\ddot{\theta}(t)$ can be replaced by $v(t)$ due to (8).

Finally, combining this with the prior first-order differential system allows
us to conceptually formulate the problem of computing a minimum-time motion
along the path $Q$ as:

(16) $\displaystyle\min_{v,T,\zeta}$ $\displaystyle\leavevmode\nobreak\
T+\gamma\int_{0}^{T}v(t)^{2}\mathrm{d}t$ subject to
$\displaystyle\dot{\textbf{z}}(t)$ $\displaystyle=\begin{pmatrix}0&1\\\
0&0\end{pmatrix}\textbf{z}(t)+\begin{pmatrix}0\\\
1\end{pmatrix}v(t),\quad\textbf{z}(0)=\textbf{z}_{0},\leavevmode\nobreak\
\textbf{z}(T)=\textbf{z}_{T},$ $\displaystyle 0$
$\displaystyle=\widetilde{M}(\textbf{z}(t),v(t),\zeta(t)),$
$\displaystyle\zeta(t)$ $\displaystyle\in[-1,1]^{6}.$

The constraint on $\zeta(t)$ is added since the trigonometric structure of
(13) is not directly encoded in the OCP, while $\gamma\geq 0$ is a
regularisation parameter.

Figure 6. Effect of the regularisation parameter $\gamma$ on the position of
the optimised traps.

Three subplots A, B and C show a cardioid shape in two dimensional space. The
x-axis shows the horizontal spatial component from -5 to 5 centimetres, at
steps of 5, while the y-axis shows the vertical from 8 to 16 centimetres, at
steps of 4. Subplot A shows the computed trap positions when the
regularisation parameter gamma is low. Sudden jumps at the highest and lowest
vertical position, as well as at the sharp edge of the shape are highlighted.

The case $\gamma=0$ corresponds to strictly minimising time, which typically
leads to an aggressive use of the forces around the trap. The example in
Figure 6(A) shows the trap accelerating the particle before arriving at the
corner, then applying aggressive deceleration before it reaches the corner. At
the top and bottom of the shape, again the particle is accelerated strongly,
then it is suddenly decelerated by the traps placed behind it, in order to
move around the curve. This would be the optimum solution in an ideal case
(i.e., a device working _exactly_ as per our model), but inaccuracies in the
real device can make them unstable. Moreover, we observed the overall
reductions in rendering time to usually be quite small.

The regularisation ($\gamma>0$) enables _nearly_ time-optimal solutions. More
specifically, this is a strictly convex regularisation that penalises high
magnitudes of the virtual input $v=\ddot{\theta}$. This is most closely
related to the accelerations applied to the particle, hence avoiding
aggressive acceleration/deceleration as shown in Figures 6(B) and (C). Several
heuristics could be proposed to automate selection of a suitable value for
$\gamma$ (e.g., use smallest $\gamma$ ensuring that the dot product of
$\dot{\textbf{p}}(t)$ and $\dot{\textbf{u}}(t)$ remains always positive), but
this step is considered beyond the scope of the current paper.

Note that the OCP in (16) yields the (nearly) optimal timing along the path
$Q$, respecting the particle dynamics and hence solving _Challenge 1_.
Moreover, it yields the required forces through $\zeta(t)$. However, it does
not give the trap trajectory $\textbf{u}(t)$ that generates these forces. To
obtain the trap trajectory and thus solve _Challenge 2_ , we refine (16) in
the following section.

###  4.3. Computing the Trap Trajectory

The second stage in our approach deals with _Challenge 2_ , computing the trap
trajectory $\textbf{u}(t)$ based on the solution of (16). In theory, we would
use the values for $\zeta_{i}(t)$, $i=1,...,6$ and $q_{j}(\theta(t))$,
$j\in\\{x,y,z\\}$ to determine the particle position and force required at
each moment in time, obtaining the required trap position $\textbf{u}(t)$ by
solving (13).

In practice, solving $\textbf{u}(t)$ from (13) numerically is not trivial,
particularly for values of $\zeta_{i}$ close to $\pm 1$, where numerical
instabilities could occur, particularly for $\zeta_{1}$ and $\zeta_{4}$.

We attenuate these difficulties by providing further structure to the OCP
(16), specifically regarding $\zeta$:

(17) $\zeta_{2}^{2}+\zeta_{3}^{2}=1,\quad\zeta_{5}^{2}+\zeta_{6}^{2}=1.$

We also constrain the solvability of (13) by using a constant back-off
$\varepsilon\in\interval[open]{0}{1}$ in order to avoid numerical
instabilities:

(18) $-1+\varepsilon\leq\zeta_{i}\leq 1-\varepsilon,\quad i=1,...,6.$

For the sake of compact notation, we summarise the above constraints in the
following set notation

(19) $\mathcal{Z}:=\left\\{\zeta\in\mathbb{R}^{6}\leavevmode\nobreak\
|\leavevmode\nobreak\ \eqref{eq:zeta_trig_pythagoras}\text{ and
}\eqref{eq:consistConstraints}\text{ are satisfied}\right\\}.$

The final OCP is then given by

(20) $\displaystyle\min_{v,T,\zeta}$ $\displaystyle\leavevmode\nobreak\
T+\gamma\int_{0}^{T}v(t)^{2}\mathrm{d}t$ subject to
$\displaystyle\dot{\textbf{z}}(t)$ $\displaystyle=\begin{pmatrix}0&1\\\
0&0\end{pmatrix}\textbf{z}(t)+\begin{pmatrix}0\\\
1\end{pmatrix}v(t),\quad\textbf{z}(0)=\textbf{z}_{0},\leavevmode\nobreak\
\textbf{z}(T)=\textbf{z}_{T},$ $\displaystyle 0$
$\displaystyle=\widetilde{M}(\textbf{z}(t),v(t),\zeta(t)),$
$\displaystyle\zeta(t)$ $\displaystyle\in\mathcal{Z}.$

Finally, given $\zeta(t)\in\mathcal{Z}$ and $\theta(t)$ from the solution of
(20), we numerically solve (13) for $u_{i}$, $i\in\\{x,y,z\\}$, using Powell’s
dog leg method (Mangasarian, 1994). Please note that higher values of
$\varepsilon$ will limit the magnitudes of the forces exploited by our
approach. A simple solution is to perform a linear search over $\varepsilon$,
only increasing its value and recomputing the solution if Powell’s method
fails to converge to a feasible trap location. For details on discretising the
OCP (20), we refer to the Supplementary Material S3.

##  5\. Evaluation

This section evaluates the capability of our algorithm to support content
creation by generating physically feasible trap trajectories given only a
reference path (i.e., a shape) as an input. We select a range of shapes and
analyse the effects each step in our approach has on the final achievable
size, rendering frequency, and reconstruction error, for each shape. We also
inspect how these steps influence the presence of visual distortions in the
end result. Finally, we analyse the resulting acceleration profiles.

Figure 7. Test shapes used for evaluation, as rendered by our _OptiTrap_
approach. The circle (A) provides a trivial timing case. The cardioid (B) and
squircle (C) represent simple shapes featuring sharp corners and straight
lines. The fish (D) is used as an example of composite shape featuring
corners, straight lines, and curves.

The image consists of four photos A, B, C and D. The photos show illuminated
outlines of the different shapes rendered by the levitation interface in mid-
air.

###  5.1. Test Shapes

We used four test shapes for the evaluation: a circle, a cardioid, a squircle,
and a fish, all shown in Figure 7.

The circle is selected as a trivial case in terms of timing. It can be
optimally sampled by using a constant angular speed (i.e., constant
acceleration), and is there to test whether our solution converges towards
such optimum solutions. The cardioid and squircle represent simple shapes
featuring sharp corners (cardioid) and straight lines (squircle), both of them
challenging features. Finally, the fish is selected as an example composite
shape, featuring all elements (i.e., corners, straight lines, and curves).

All of these shapes can be parameterised by low-frequency sinusoids (see
Supplementary Material S4), ensuring the path parameter $\theta$ progresses
slowly in areas of high curvature. This provides a good starting point for our
baseline comparison that we explain in Section 5.2. While our approach works
for content in 3D, we perform the evaluation on shapes in a 2D plane of the 3D
space to facilitate the analysis of accelerations in Section 5.6, which we
will perform in terms of horizontal and vertical accelerations (i.e., the
independent factors given our axis-symmetric model of forces), allowing us to
validate _OptiTrap_ ’s awareness of particle directions and how close it can
get to the maximum accelerations allowed by the dynamics of acoustic traps, as
discussed in Section 4.1.

Notice in Figure 7, that due to the timing differences, different shape
elements exhibit different levels of brightness. To obtain homogeneous
brightness along the shape, please refer to the solution proposed by (Hirayama
et al., 2019), where the particle illumination is adjusted to the particle
speed.

###  5.2. Conditions Compared

We compare three approaches to rendering levitated shapes, where each approach
subsequently addresses one of the challenges introduced in Section 2. This
will help us assess the impact that each challenge has on the final results
obtained.

The first condition is a straight-forward _Baseline_ , with homogeneous
sampling of the path parameter, which matches the example strategy shown in
Figure 3(B), and placing traps where the particle should be. The _Baseline_
still does not address any of the challenges identified (i.e., optimum timing
or trap placement), but it matches approaches used in previous works (Hirayama
et al., 2019; Plasencia et al., 2020; Fushimi et al., 2020) and will
illustrate their dependence on the specific shape and initial parametrisation
used.

The second condition, _OCP_Timing_ , makes use of our _OptiTrap_ approach to
compute optimum and feasible timing (see Subsection 4.2.1), but still ignores
_Challenge 2_ , assuming that particles match trap location (i.e., it skips
Subsection 4.2.2). The third condition is the full _OptiTrap_ approach, which
considers feasibility and trap-particle dynamics and deals with both the
timing and trap placement challenges.

These conditions are used in a range of comparisons involving size, frequency,
reconstruction error as well as comparative analysis of the effects of each
condition on the resulting particle motion, detailed in the following
sections.

###  5.3. Maximum Achievable Sizes

This section focuses on the maximum sizes that can be achieved for each test
shape and condition, while retaining an overall rendering time of 100ms111The
circle rendered at 100ms would exceed the size of our device’s working volume.
Thus, this shape was rendered at 67ms. (i.e., PoV threshold). For each of
these cases, we report the maximum achievable size and reliability, which were
determined as follows.

Maximum size is reported in terms of shape width and in terms of meters of
_content per second_ rendered (Plasencia et al., 2020), and their
determination was connected to the feasibility of the trap trajectories,
particularly for the _Baseline_ condition.

Table 2. Maximum shape width and meters of content per second achieved with the _Baseline_ , _OCP_Timing_ , and _OptiTrap_ approach at constant rendering frequency. Successful trials out of 10. | _Baseline_ | _OCP_Timing_ | _OptiTrap_ | Increase
---|---|---|---|---
Shape | Freq. | Width | Content per | Success | Width | Content per | Success | Width | Content per | Success | in Width w.r.t.
| (Hz) | (cm) | Second (m) | Rate | (cm) | Second (m) | Rate | (cm) | Second (m) | Rate | the _Baseline_
Circle | 15 | 7.00 | 3.30 | 10/10 | 7.00 | 3.30 | 10/10 | 7.00 | 3.30 | 10/10 | 0%
Cardioid | 10 | 8.05 | 2.48 | 10/10 | 9.09 | 2.80 | 10/10 | 9.09 | 2.80 | 10/10 | 12.9%
Squircle | 10 | 0.80 | 0.29 | 9/10 | 5.30 | 1.90 | 10/10 | 5.30 | 1.90 | 10/10 | 562.5%
Fish | 10 | 7.77 | 2.42 | 10/10 | 8.76 | 2.75 | 10/10 | 8.76 | 2.75 | 10/10 | 12.7%

More specifically, we determined maximum feasible sizes for the _Baseline_
condition by conducting an iterative search. That is, we increased the size at
each step and tested the reliability of the resulting shape. We considered the
shape feasible if it could be successfully rendered at least 9 out of 10 times
in the actual levitator. On a success, we would increase size by increasing
the shape width by half a centimetre. On a failure, we would perform a binary
search between the smallest size failure and the largest size success,
stopping after two consecutive failures and reporting the largest successful
size. This illustrates the kind of trial and error that a content designer
would need to go through without our approach and it was the most time
consuming part of this evaluation.

For the other two conditions (i.e., _OCP_Timing_ and _OptiTrap_), maximum
sizes could be determined without needing to validate their feasibility in the
final device. Again, we iteratively increased the shape size using the same
search criteria as before (i.e., 5mm increases, binary search). We assumed any
result provided would be feasible (which is a reasonable assumption; see
below) and checked the resulting total rendering time, increasing the size if
the time was still less than 100ms. At each step, a linear search was used,
iteratively increasing $\varepsilon$, until feasible trap locations could be
found (i.e., equation (13) could be solved without numerical instabilities).
Adjusting the value of the regularisation parameter $\gamma$ was done by
visually inspecting the resulting trap trajectories, until no discontinuities
could be observed (see Figure 6). Please note that this is a simple and quick
task, which, unlike the _Baseline_ condition, does not involve actual testing
on the levitation device and could even be automated. All solutions provided
can be assumed feasible, and the designer only needs to choose the one that
better fits their needs.

Once the maximum achievable sizes were determined, these were tested in the
actual device. We determined the feasibility of each shape and condition by
conducting ten tests and reporting the number of cases where the particle
succeeded to reveal the shape. This included the particle accelerating from
rest, traversing/revealing the shape for 6 seconds and returning to rest.

Table 2 summarises the results achieved for each test shape and condition.
First of all, it is worth noting that all trials were successful for
_OptiTrap_ and _OCP_Timing_ approaches, confirming our assumption that the
resulting paths are indeed feasible and underpinning _OptiTrap_ ’s ability to
avoid trial and error on the actual device during content creation.

The results and relative improvements in terms of size vary according to the
particular condition and shape considered. For instance, it is interesting to
see that all conditions yield similar final sizes for the circle, showing that
_OptiTrap_ (and _OCP_Timing_) indeed converge towards optimum solutions.

More complex shapes where the _Baseline_ parametrisation is not optimal (i.e.,
cardioid, squircle, and fish), show increases in size when using _OCP_Timing_
and _OptiTrap_. More interestingly, the increases in size vary greatly between
shapes, showing increases of around $12\%$ for the fish and cardioid, and up
to $562\%$ for the squircle. This is the result of the explicit
parametrisation used, with reduced speeds at corners in the fish and cardioid
cases, but not in the case of the squircle.

It is also interesting to see that _OCP_Timing_ and _OptiTrap_ maintain high
values of _content per second_ rendered, independently of the shape. That is,
while the performance of the _Baseline_ approach is heavily determined by the
specific shape (and parametrisation) used, the OCP-based approaches
(_OCP_Timing_ and _OptiTrap_) yield results with consistent _content per
second_ , determined by the capabilities of the device (but not so much by the
specific shape). Finally, please note that no changes in terms of maximum size
can be observed between _OCP_Timing_ and _OptiTrap_.

###  5.4. Maximum Rendering Frequencies

Table 3. Maximum rendering frequencies achieved with the _Baseline_ , _OCP_Timing_ , and _OptiTrap_ approach, for shapes of equal size. Successful trials out of 10. | _Baseline_ | _OCP_Timing_ | _OptiTrap_
---|---|---|---
Shape | Freq. | Success | Freq. | Success | Freq. | Success
| (Hz) | Rate | (Hz) | Rate | (Hz) | Rate
Circle | 15 | 10/10 | 15 | 10/10 | 15 | 10/10
Cardioid | 8 | 10/10 | 10 | 10/10 | 10 | 10/10
Squircle | 4 | 10/10 | 10 | 10/10 | 10 | 10/10
Fish | 9 | 9/10 | 10 | 10/10 | 10 | 10/10

In this second evaluation, we assessed the effect of the timing strategy on
the maximum achievable rendering frequencies. To do this, we selected the
maximum achievable sizes obtained in the prior evaluation for each shape. We
then reproduced such sizes with the _Baseline_ approach, searching for the
maximum frequency at which this approach could reliably render the shape
(using the same searching and acceptance criteria as above).

That is, while we knew _OCP_Timing_ and _OptiTrap_ could provide 10Hz for
these shapes and sizes, we wanted to determine the maximum frequency at which
the _Baseline_ would render them, as to characterise the benefits provided by
optimising the timing.

The results are summarised in Table 3. As expected, no changes are produced
for the circle. However, there is a 25% increase for the cardioid, 150% for
the squircle, and 11% increase for the fish using the _OptiTrap_ method.
Again, please note that no differences can be observed in terms of maximum
frequency between _OCP_Timing_ and _OptiTrap_.

###  5.5. Reconstruction Accuracy

To evaluate the reconstruction accuracy, we recorded each of the trials in our
maximum size evaluations (see Section 5.3) using an OptiTrack camera system.
We then computed the Root Mean Squared Error (RMSE) between the recorded data
and the intended target shape, taking 2s of shape rendering into account
(exclusively during the cyclic part, ignoring ramp-up/ramp-down). For a fair
comparison across trials, we normalise the RMSE with respect to the traversed
paths, by dividing the RMSE by the total path length of each individual test
shape. The results are summarised in Table 4.

Table 4. RMSE and Path-normalised (PN) RMSE with respect to the total path length of each test shape, for the _Baseline_ , _OCP_Timing_ , and _OptiTrap_ approach, at constant rendering frequency. | _Baseline_ | _OCP_Timing_ | _OptiTrap_
---|---|---|---
Shape | RMSE | PN | RMSE | PN | RMSE | PN
| (cm) | RMSE | (cm) | RMSE | (cm) | RMSE
Circle | 0.196 | 0.893 | 0.146 | 0.666 | 0.114 | 0.521
Cardioid | 0.142 | 0.575 | 0.153 | 0.545 | 0.120 | 0.429
Squircle | 0.056 | 1.956 | 0.0824 | 0.434 | 0.080 | 0.421
Fish | 0.111 | 0.458 | 0.147 | 0.536 | 0.0581 | 0.212

On average, both _OCP_Timing_ and _OptiTrap_ provide better results in terms
of accuracy, when compared to the _Baseline_. It is particularly worth noting
that the accuracy results, in terms of the raw RMSE for _OptiTrap_ and
_OCP_Timing_ , are better for the cardioid when compared to the _Baseline_ ,
even though a larger shape is being rendered. While the raw RMSE of
_OCP_Timing_ and _OptiTrap_ is slightly larger for the squircle, we need to
take into account that the size rendered by _OCP_Timing_ and _OptiTrap_ is
almost 6 times larger. Comparing accuracy in terms of normalised RMSE shows
consistent increases of accuracy for _OptiTrap_ compared to the _Baseline_ ,
with overall decreases in the RMSE of $41.7\%$ (circle), $25.4\%$ (cardioid),
$78.5\%$ (squircle), and $53.8\%$ (fish). _OCP_Timing_ shows a decrease of
$25.4\%$, $5.14\%$, and $77.8\%$ of the normalised RMSE for the circle,
cardioid, and squircle, with the exception of the fish, where the normalised
RMSE increased by $17.1\%$, when compared to the _Baseline_.

Comparing the reconstruction accuracy of _OCP_Timing_ and _OptiTrap_
highlights the relevance of considering trap dynamics to determine trap
locations (i.e., Section 4.3). _OptiTrap_ consistently provides smaller RMSE
than _OCP_Timing_ , as a result of considering and accounting for the trap-to-
particle displacements required to apply specific accelerations. We obtain a
$21.8\%$, $21.3\%$, $2.9\%$ and $60.5\%$ decrease in the normalised RMSE, for
the circle, cardioid, squircle, and fish, respectively. Such differences are
visually illustrated in Figure 8, showing the effects on the cardioid and
fish, for the _OCP_Timing_ and _OptiTrap_ approaches. It is worth noting how
placing the traps along the reference path results in the cardioid being
horizontally stretched, as the particle needs to retain larger distances to
the trap to keep the required acceleration (horizontal forces are weaker than
vertical ones). This also results in overshooting of the particle at corner
locations, which can be easily observed at the corner of the cardioid and in
the fins of the fish. For completeness, the visual comparison for the
remaining two test shapes is provided in the Supplementary Material S5.

Figure 8. Visual comparison between shapes rendered with traps located
according to particle dynamics using _OptiTrap_ (left) and traps placed along
the reference path using _OCP_Timing_ (right) for the cardioid (top) and the
fish (bottom) test shapes. Please note undesired increases in size and error
in sharp features, such as corners.

The image consists of a two-by-two photo matrix. There are two photos of a
rendered cardioid in the first row, and fish in the second row. The cardioid
and the fish in the first column are rendered using the OptiTrap algorithm,
and in the second column the two shapes are rendered using OCP_Timing. The
levitated graphics in the second column are a bit larger in size and less
precise than those in the first column.

###  5.6. Analysis of Acceleration profiles

Figure 9. Particle acceleration and speed for the fish rendered at 10 Hz using
_OptiTrap_ (solid black) and the _Baseline_ (dash-dotted purple). While the
speed profiles (D) of both trap trajectories are similar, the trap trajectory
generated by our approach is feasible, whereas the one generated by the
_Baseline_ is not. The main reason is that the _Baseline_ exceeds the feasible
horizontal acceleration, while our approach caps it to feasible values (B).
Our approach compensates by using higher (available) vertical acceleration
(C). Note that the times where the _Baseline_ applies a higher total
acceleration (A) than our approach are those where our approach respects the
constraints on the feasible horizontal acceleration (B).

The figure consists of four subplots. Subplot A shows the total acceleration
in metres per seconds squared from 0 to 400, in steps of 200. Subplot B shows
the horizontal acceleration in metres per seconds squared from -200 to 400, in
steps of 200. Subplot C shows the vertical acceleration in metres per seconds
squared from -500 to 500, in steps of 500. Subplot D shows the total speed in
metres per seconds from 0 to 4, in steps of 2. All quantities are plotted as
functions of time, going from 0 to 100 milliseconds, at steps of 20.

Finally, we examined the acceleration profiles produced by _OptiTrap_ and how
these differ from the pre-determined acceleration profiles of the _Baseline_.
Please note that we only compare and discuss _OptiTrap_ and _Baseline_ in
Figure 9, and only for the fish shape. _OCP_Timing_ is not included, as it
results in similar acceleration profiles as _OptiTrap_. Only the fish is
discussed, as it already allows us to describe the key observations that can
be derived from our analysis. For completeness, the acceleration profiles for
all remaining test shapes are included in the Supplementary Material S6.

As introduced above, Figure 9 shows the particle accelerations and speeds of a
particle revealing a fish shape of maximum size (i.e., a width of 8.76cm),
rendered over 100ms using the _OptiTrap_ and _Baseline_ approaches. It is
worth noting that while _OptiTrap_ succeeded in rendering this shape,
_Baseline_ did not (the maximum width for _Baseline_ was 7.77cm).

A first interesting observation is that although _OptiTrap_ provides lower
total accelerations than the _Baseline_ during some parts of the path (see
Figure 9(A)), it still manages to reveal the shape in the same time. The key
observation here is that the regions where the total acceleration is lower for
_OptiTrap_ match with the parts of the path where the horizontal acceleration
is very close to its maximum; for example, note the flat regions in Figure
9(B) of around $\pm 300m/{s}^{2}$). This is an example of _OptiTrap_ ’s
awareness of the dynamics and capabilities of the actual device, limiting the
acceleration applied as to retain the feasibility of the path.

Second, it is worth noting that maximum horizontal accelerations are
significantly smaller than vertical accelerations. As such, it does make sense
for horizontal displacements to become the limiting factor. In any case,
neither acceleration exceeds its respective maximum value. The _OptiTrap_
approach recovers any missing time by better exploiting areas where
acceleration is unnecessarily small. It is also worth noting that the value of
the horizontal and vertical accelerations are not simply being capped to a
maximum acceleration value per direction (e.g., simultaneously maxing out at
$\pm 300m/{s}^{2}$ and $\pm 600m/{s}^{2}$ in the horizontal and vertical
directions), as such cases are not feasible according to the dynamics of our
acoustic traps.

Third, it is interesting to note that the final acceleration profile is
complex, retaining little resemblance with the initial parametrisation used.
This is not exclusive for this shape and can be observed even in relatively
simple shapes, such as the cardioid in Figure 3(D) or the remaining shapes,
provided in the Supplementary Material S6.

The complexity of these profiles is even more striking if we look at the final
speeds resulting from both approaches (see Figure 9(D)). Even if the
acceleration profiles of _Baseline_ and _OptiTrap_ are very different, their
velocity profiles show only relatively subtle differences. But even if these
differences are subtle, they mark the difference between a feasible path
(i.e., _OptiTrap_) and a failed one (i.e., _Baseline_). These complex yet
subtle differences also illustrate how it is simply not sensible to expect
content designers to deal with such complexity, and how _OptiTrap_ is a
necessary tool to enable effective exploitation of PoV levitated content.

##  6\. Discussion

This paper presented _OptiTrap_ , an automated approach for optimising timings
and trap placements, as to achieve feasible target shapes. We believe this is
a particularly relevant step for the adoption of levitation PoV displays, as
it allows the content creator to focus on the shapes to present, with feasible
solutions being computed automatically, while making effective usage of the
capabilities of the device. As such, we hope _OptiTrap_ to become an
instrumental tool in helping explore the actual potential of these displays.

However, _OptiTrap_ is far from a complete content creation tool. Such a tool
should consider the artist’s workflows and practices. Similarly, testing and
identifying most useful heuristics to tune our approach (e.g., the
regularisation) or visualisation tools identifying tricky parts of the shapes
(i.e., requiring high accelerations) should be included as a part of this
process. Our goal is simply to provide the base approach enabling this kind of
tools.

Even this base approach can be extended in a variety of ways. Our levitator
prototype is built from off-the-shelf hardware, and is still subject to
inaccuracies that result in distortions in the sound-fields generated (Fushimi
et al., 2019). As such, more accurate levitation hardware or, alternatively, a
model reflecting the dynamics of the system in a more accurate manner would be
the most obvious pathway to improve our approach. It is worth noting that both
factors should be advanced jointly. A more accurate model could also be less
numerically stable, potentially leading to worse results if the hardware is
not accurate enough.

The high update rates of 10kHz required by the levitator and the millisecond
delays introduced by optical tracking systems indicate that closed-loop
approaches can be both promising and challenging avenues to explore. Assuming
a tracking device synchronised with the levitation device (as to map current
trap locations with real particle positions), _OptiTrap_ could be combined
with learning-based approaches. These learning-based approaches will require
example paths (i.e., with initial timing and trap placement), and their
achievable complexity and convergence will be limited by the examples
provided. In such cases, _OptiTrap_ can be used to always generate feasible
initial trajectories (i.e., to avoid system restarts on failure). Thus
learning-based approaches could be used to further refine _OptiTrap_ beyond
our current model, as to account for device inaccuracies such as those
discussed by (Fushimi et al., 2019).

However, higher gains can be obtained from more radical changes in the
approach. For instance, our approach optimises the timing and trap placement,
but it does not modify the target shapes. As illustrated in Figure 9, slight
changes in speed have significant effects on the acceleration profiles (and
feasibility) of the shapes. This is even more prominent for position, where
small changes can heavily influence the acceleration and feasibility of target
shapes. As such, approaches exploiting subtle modifications to the shape could
lead to significant gains in rendering performance.

Another interesting possibility would be extending our approach to use several
particles. As shown in (Plasencia et al., 2020), while the use of several
particles does not increase the overall power that can be leveraged, it does
allow for increased flexibility. That is, the intensity/forces of each trap
can be individually and dynamically adjusted, as to match the needs of the
region of the path that each particle is revealing. Also, particles can each
be rendering specific independent features, so the content is not limited to a
single connected path, and the particles do not waste time/accelerations
traversing parts of the path that will not be illuminated (i.e., visible).
This approach, however, entails significant challenges. The first obvious
challenge is the reliability of the intensity control of the traps. (Plasencia
et al., 2020) demonstrate accurate control of the stiffness at the centre of
the traps, but the effects of multiple (interfering) traps in each trap’s
topology (i.e., how forces distribute around the trap) is yet to be studied. A
second challenge is that each particle is not forced to traverse the path at
the same speed/rates, with such independent timing progression becoming an
additional degree of freedom to account for.

Finally, further extensions to our work can come from its application to
domains other than PoV displays. An obvious next step would be to adapt
_OptiTrap_ to photophoretic displays, which trap particles using optical traps
instead of acoustic traps (Smalley et al., 2018; Kumagai et al., 2021). This
would involve including a model of the dynamics of such optical traps, but it
can also involve further challenges, such as modelling the response times of
galvanometers and LC panels involved in creating the trap.

Our approach can also be adapted to applications requiring objects to be
transported quickly and accurately. For example, contactless transportation of
matter has a wealth of applications in areas such as the study of physical
phenomena, biochemical processes, materials processing, or pharmaceutics
(Foresti et al., 2013). Our method can help solving such problems by directly
computing a rest-to-rest solution of the matter to be transported, given an
estimation of the mass of that matter, the acoustic force acting on the
matter, and a path from start to target position.

##  7\. Conclusion

In this paper we proposed the first structured numerical approach to compute
trap trajectories for acoustic levitation displays. _OptiTrap_ automatically
computes physically feasible and nearly time-optimal trap trajectories to
reveal generic mid-air shapes, given only a reference path. Building on a
novel multi-dimensional approximation of the acoustic forces around the trap,
we formulate and show how to solve a non-linear path following problem without
requiring or exploiting differential flatness of the system dynamics. We
demonstrate increases of up to 563% in size and up to 150% in frequency for
several shapes. Additionally, we obtain better reconstruction accuracy with up
to a 79% decrease in the path-normalised RMSE. While previously, feasible trap
trajectories needed to be tuned manually for each shape and levitator, our
approach requires calibration of each individual levitator just once. We are
confident that the ideas in this paper could form the basis for future content
authoring tools for acoustic levitation displays and bring them a key step
closer to real-world applications.

###### Acknowledgements.

This research has received funding from the European Union’s Horizon 2020
research and innovation programme under grant agreement #737087 (Levitate) and
from the AHRC UK-China Research-Industry Creative Partnerships (AH/T01136X/2).

## References

  * (1)
  * Andrade et al. (2016) Marco A. B. Andrade, Anne L. Bernassau, and Julio C. Adamowski. 2016.  Acoustic levitation of a large solid sphere.  _Applied Physics Letters_ 109, 4 (2016), 044101\.  https://doi.org/10.1063/1.4959862 arXiv:https://doi.org/10.1063/1.4959862
  * Bachynskyi et al. (2018) Myroslav Bachynskyi, Viktorija Paneva, and Jörg Müller. 2018.  LeviCursor: Dexterous Interaction with a Levitating Object. In _Proceedings of the 2018 ACM International Conference on Interactive Surfaces and Spaces_ (Tokyo, Japan) _(ISS ’18)_. ACM, New York, NY, USA, 253–262.  https://doi.org/10.1145/3279778.3279802
  * Bock and Plitt (1984) Hans Georg Bock and Karl-Josef Plitt. 1984.  A multiple shooting algorithm for direct solution of optimal control problems.  _IFAC Proceedings Volumes_ 17, 2 (1984), 1603–1608.
  * Bowen et al. (1974) Richard W. Bowen, Jordan Pola, and Leonard Matin. 1974\.  Visual persistence: Effects of flash luminance, duration and energy.  _Vision Research_ 14, 4 (1974), 295 – 303.  https://doi.org/10.1016/0042-6989(74)90079-0
  * Bruus (2012) Henrik Bruus. 2012\.  Acoustofluidics 7: The acoustic radiation force on small particles.  _Lab Chip_ 12 (2012), 1014–1021.  Issue 6. https://doi.org/10.1039/C2LC21068A
  * Butcher (2016) John Charles Butcher. 2016\.  _Numerical methods for ordinary differential equations_.  John Wiley & Sons.
  * Faulwasser (2012) Timm Faulwasser. 2012\.  _Optimization-based solutions to constrained trajectory-tracking and path-following problems_.  Ph. D. Dissertation. Otto-von-Guericke-Universität Magdeburg, Universitaetspl. 2.
  * Faulwasser et al. (2014) Timm Faulwasser, Veit Hagenmeyer, and Rolf Findeisen. 2014\.  Constrained Reachability and Trajectory Generation for Flat Systems.  _Automatica_ 50, 4 (2014), 1151–1159.  https://doi.org/10.1016/j.automatica.2014.02.011
  * Faulwasser et al. (2017) Timm Faulwasser, Tobias Weber, Pablo Zometa, and Rolf Findeisen. 2017.  Implementation of Nonlinear Model Predictive Path-Following Control for an Industrial Robot.  _IEEE Trans. Contr. Syst. Techn._ 25, 4 (2017), 1505–1511.  https://doi.org/10.1109/TCST.2016.2601624
  * Fender et al. (2021) Andreas Rene Fender, Diego Martinez Plasencia, and Sriram Subramanian. 2021\.  _ArticuLev: An Integrated Self-Assembly Pipeline for Articulated Multi-Bead Levitation Primitives_.  Association for Computing Machinery, New York, NY, USA.  https://doi.org/10.1145/3411764.3445342
  * Fliess et al. (1995) Michel Fliess, Jean Lévine, Philippe Martin, and Pierre Rouchon. 1995.  Flatness and defect of non-linear systems: introductory theory and examples.  _Int. J. Contr._ 61, 6 (1995), 1327–1361.
  * Foresti et al. (2013) Daniele Foresti, Majid Nabavi, Mirko Klingauf, Aldo Ferrari, and Dimos Poulikakos. 2013.  Acoustophoretic contactless transport and handling of matter in air.  _Proceedings of the National Academy of Sciences_ 110, 31 (2013), 12549–12554.  https://doi.org/10.1073/pnas.1301860110 arXiv:https://www.pnas.org/content/110/31/12549.full.pdf
  * Freeman et al. (2018) Euan Freeman, Julie Williamson, Sriram Subramanian, and Stephen Brewster. 2018.  Point-and-Shake: Selecting from Levitating Object Displays. In _Proceedings of the 2018 CHI Conference on Human Factors in Computing Systems_ (Montreal QC, Canada) _(CHI ’18)_. ACM, New York, NY, USA, Article 18, 10 pages.  https://doi.org/10.1145/3173574.3173592
  * Fushimi et al. (2020) Tatsuki Fushimi, Bruce W. Drinkwater, and Thomas L. Hill. 2020\.  What is the ultimate capability of acoustophoretic volumetric displays?  _Applied Physics Letters_ 116, 24 (2020), 244101\.  https://doi.org/10.1063/5.0008351 arXiv:https://doi.org/10.1063/5.0008351
  * Fushimi et al. (2018) Tatsuki Fushimi, Bruce W. Drinkwater, Asier Marzo, and Thomas L. Hill. 2018.  Nonlinear trapping stiffness of mid-air single-axis acoustic levitators.  _Applied Physics Letters_ 113, 3 (2018), 034102\.  https://doi.org/10.1063/1.5034116 arXiv:https://doi.org/10.1063/1.5034116
  * Fushimi et al. (2019) Tatsuki Fushimi, Asier Marzo, Bruce W. Drinkwater, and Thomas L. Hill. 2019.  Acoustophoretic volumetric displays using a fast-moving levitated particle.  _Applied Physics Letters_ 115, 6 (2019), 064101\.  https://doi.org/10.1063/1.5113467 arXiv:https://doi.org/10.1063/1.5113467
  * Geijtenbeek and Pronost (2012) Thomas Geijtenbeek and Nicolas Pronost. 2012.  Interactive Character Animation Using Simulated Physics: A State-of-the-Art Review.  _Comput. Graph. Forum_ 31, 8 (Dec. 2012), 2492–2515.  https://doi.org/10.1111/j.1467-8659.2012.03189.x
  * Hirayama et al. (2019) Ryuji Hirayama, Diego Martinez Plasencia, Nobuyuki Masuda, and Sriram Subramanian. 2019\.  A volumetric display for visual, tactile and audio presentation using acoustic trapping.  _Nature_ 575, 7782 (2019), 320–323.
  * Inoue et al. (2019) Seki Inoue, Shinichi Mogami, Tomohiro Ichiyama, Akihito Noda, Yasutoshi Makino, and Hiroyuki Shinoda. 2019.  Acoustical boundary hologram for macroscopic rigid-body levitation.  _The Journal of the Acoustical Society of America_ 145, 1 (2019), 328–337.
  * Kumagai et al. (2021) Kota Kumagai, Shun Miura, and Yoshio Hayasaki. 2021.  Colour volumetric display based on holographic-laser-excited graphics using drawing space separation.  _Scientific reports_ 11, 1 (2021), 1–9.
  * Mangasarian (1994) Olvi L Mangasarian. 1994\.  _Nonlinear programming_.  SIAM.
  * Marzo and Drinkwater (2019) Asier Marzo and Bruce W. Drinkwater. 2019.  Holographic acoustic tweezers.  _Proceedings of the National Academy of Sciences_ 116, 1, 84–89.  https://doi.org/10.1073/pnas.1813047115 arXiv:https://www.pnas.org/content/116/1/84.full.pdf
  * Marzo et al. (2015) Asier Marzo, Sue Ann Seah, Bruce W. Drinkwater, Deepak Ranjan Sahoo, Benjamin Long, and Sriram Subramanian. 2015.  Holographic Acoustic Elements for Manipulation of Levitated Objects.  _Nature communications_ 6 (Oct. 2015), 8661.  http://dx.doi.org/10.1038/ncomms9661
  * Morales et al. (2021) Rafael Morales, Iñigo Ezcurdia, Josu Irisarri, Marco A. B. Andrade, and Asier Marzo. 2021\.  Generating Airborne Ultrasonic Amplitude Patterns Using an Open Hardware Phased Array.  _Applied Sciences_ 11, 7 (2021).  https://doi.org/10.3390/app11072981
  * Morales et al. (2019) Rafael Morales, Asier Marzo, Sriram Subramanian, and Diego Martínez. 2019.  LeviProps: Animating Levitated Optimized Fabric Structures Using Holographic Acoustic Tweezers. In _Proceedings of the 32Nd Annual ACM Symposium on User Interface Software and Technology_ (New Orleans, LA, USA) _(UIST ’19)_. ACM, New York, NY, USA, 651–661.  https://doi.org/10.1145/3332165.3347882
  * Nägeli et al. (2017) Tobias Nägeli, Lukas Meier, Alexander Domahidi, Javier Alonso-Mora, and Otmar Hilliges. 2017\.  Real-time planning for automated multi-view drone cinematography.  _ACM Transactions on Graphics (TOG)_ 36, 4 (2017), 1–10.
  * Nielsen and Maggiore (2008) Christopher Nielsen and Manfredi Maggiore. 2008.  On local transverse feedback linearization.  _SIAM Journal on Control and Optimization_ 47 (2008), 2227–2250.
  * Norasikin et al. (2019) Mohd Adili Norasikin, Diego Martinez-Plasencia, Gianluca Memoli, and Sriram Subramanian. 2019\.  SonicSpray: A Technique to Reconfigure Permeable Mid-Air Displays. In _Proceedings of the 2019 ACM International Conference on Interactive Surfaces and Spaces_. 113–122.
  * Norasikin et al. (2018) Mohd Adili Norasikin, Diego Martinez Plasencia, Spyros Polychronopoulos, Gianluca Memoli, Yutaka Tokuda, and Sriram Subramanian. 2018\.  SoundBender: Dynamic Acoustic Control Behind Obstacles. In _Proceedings of the 31st Annual ACM Symposium on User Interface Software and Technology_ (Berlin, Germany) _(UIST ’18)_. Association for Computing Machinery, New York, NY, USA, 247–259.  https://doi.org/10.1145/3242587.3242590
  * Ochiai et al. (2014) Yoichi Ochiai, Takayuki Hoshi, and Jun Rekimoto. 2014\.  Pixie Dust: Graphics Generated by Levitated and Animated Objects in Computational Acoustic-Potential Field.  _ACM Trans. Graph._ 33, 4, Article 85 (July 2014), 13 pages.  https://doi.org/10.1145/2601097.2601118
  * Omirou et al. (2015) Themis Omirou, Asier Marzo, Sue Ann Seah, and Sriram Subramanian. 2015.  LeviPath: Modular Acoustic Levitation for 3D Path Visualisations. In _Proceedings of the 33rd Annual ACM Conference on Human Factors in Computing Systems_ (Seoul, Republic of Korea) _(CHI ’15)_. Association for Computing Machinery, New York, NY, USA, 309–312.  https://doi.org/10.1145/2702123.2702333
  * Omirou et al. (2016) Themis Omirou, Marzo Asier Perez, Sriram Subramanian, and Anne Roudaut. 2016.  Floating charts: Data plotting using free-floating acoustically levitated representations. In _2016 IEEE Symposium on 3D User Interfaces (3DUI)_. 187–190.  https://doi.org/10.1109/3DUI.2016.7460051
  * Paneva et al. (2020) Viktorija Paneva, Myroslav Bachynskyi, and Jörg Müller. 2020\.  Levitation Simulator: Prototyping for Ultrasonic Levitation Interfaces in Virtual Reality. In _Proceedings of the 2020 CHI Conference on Human Factors in Computing Systems_ _(CHI ’20)_. ACM, New York, NY, USA.
  * Plasencia et al. (2020) Diego Martinez Plasencia, Ryuji Hirayama, Roberto Montano-Murillo, and Sriram Subramanian. 2020.  GS-PAT: High-Speed Multi-Point Sound-Fields for Phased Arrays of Transducers.  _ACM Trans. Graph._ 39, 4, Article 138 (July 2020), 12 pages.  https://doi.org/10.1145/3386569.3392492
  * Reynal et al. (2020) Maxime Reynal, Euan Freeman, and Stephen Brewster. 2020\.  Avoiding Collisions When Interacting with Levitating Particle Displays. In _Extended Abstracts of the 2020 CHI Conference on Human Factors in Computing Systems_ (Honolulu, HI, USA) _(CHI EA ’20)_. Association for Computing Machinery, New York, NY, USA, 1–7.  https://doi.org/10.1145/3334480.3382965
  * Roberts and Hanrahan (2016) Mike Roberts and Pat Hanrahan. 2016.  Generating Dynamically Feasible Trajectories for Quadrotor Cameras.  _ACM Trans. Graph._ 35, 4, Article 61 (July 2016), 11 pages.  https://doi.org/10.1145/2897824.2925980
  * Rose et al. (1996) Charles Rose, Brian Guenter, Bobby Bodenheimer, and Michael F. Cohen. 1996.  Efficient Generation of Motion Transitions Using Spacetime Constraints. In _Proceedings of the 23rd Annual Conference on Computer Graphics and Interactive Techniques_ _(SIGGRAPH ’96)_. Association for Computing Machinery, New York, NY, USA, 147–154.  https://doi.org/10.1145/237170.237229
  * Sahoo et al. (2016) Deepak Ranjan Sahoo, Takuto Nakamura, Asier Marzo, Themis Omirou, Michihiro Asakawa, and Sriram Subramanian. 2016.  JOLED: A Mid-Air Display Based on Electrostatic Rotation of Levitated Janus Objects _(UIST ’16)_. Association for Computing Machinery, New York, NY, USA, 437–448.  https://doi.org/10.1145/2984511.2984549
  * Shin and McKay (1985) Kang Shin and Neil McKay. 1985.  Minimum-time control of robotic manipulators with geometric path constraints.  _IEEE Trans. Automat. Contr._ 30, 6 (1985), 531 – 541.
  * Smalley et al. (2018) DE Smalley, E Nygaard, K Squire, J Van Wagoner, J Rasmussen, S Gneiting, K Qaderi, J Goodsell, W Rogers, M Lindsey, et al. 2018\.  A photophoretic-trap volumetric display.  _Nature_ 553, 7689 (2018), 486.
  * Stevens (1899) Walter Le Conte Stevens. 1899\.  A Text-Book of Physics–Sound.
  * Verscheure (2009) Diederik Verscheure. 2009\.  _Contribution to contact modeling and identification and optimal robot motion planning_.  Ph. D. Dissertation. Katholieke Universiteit Leuven, Faculteit Ingenieurswetenschappen, Department Werktuigkunde.
  * Von Stryk (1993) Oskar Von Stryk. 1993\.  Numerical solution of optimal control problems by direct collocation.  In _Optimal control_. Springer, 129–143.
  * Witkin and Kass (1988) Andrew Witkin and Michael Kass. 1988.  Spacetime Constraints. In _Proceedings of the 15th Annual Conference on Computer Graphics and Interactive Techniques_ _(SIGGRAPH ’88)_. Association for Computing Machinery, New York, NY, USA, 159–168.  https://doi.org/10.1145/54852.378507