## Abstract

Optically levitated nanoparticles in a vacuum offer a light–matter interface with broad and easy tunability of all key system parameters. However, the majority of previously reported experimental achievements in this area have only dealt with a single levitated object. Here, we demonstrate optical binding between multiple levitated objects confined in cross-polarized counter-propagating laser beams in a vacuum. We characterize the level of interparticle interaction, quantify its nonlinearity for various configurations of the system, and demonstrate its broad tunability. Our methodology for quantitative characterization of optically bound structures is supported by an extensive theoretical description and validated by numerical simulations. We believe the presented results represent a step toward the development of a framework of levitated optomechanics of complex coupled systems with a controlled level of coupling nonlinearity for experimental studies including, for example, mesoscopic entanglement.

© 2021 Optical Society of America under the terms of the OSA Open Access Publishing Agreement

## 1. INTRODUCTION

Optical levitation of particles in laser beams propagating in a vacuum was first demonstrated by Arthur Ashkin in his 1976 pioneering work [1]. The technique, referred to as levitated optomechanics, has recently attracted renewed interest owing to its applications in the emerging fields of macroscopic quantum mechanics [2] and stochastic thermodynamics [3]. These activities have been facilitated by the development of two strategies for cooling the center-of-mass motion of optically levitated, mesoscopic particles toward the quantum regime. The first one, independently proposed by three groups in 2010 [4–6], relies on coupling the trapped particle to an optical cavity, with the resulting efforts culminating in the very recent demonstration of ground-state cooling by Delić *et al.* [7]. The second approach is based on active feedback control [8–10], and is well suited to studying elementary statistical mechanics [11]. In particular, levitated optomechanics permits detailed studies of the Brownian motion of individual particles controllably coupled to the environment, whose positions can be tracked with high temporal and spatial resolution. This approach has yielded new insights into microscale thermodynamic processes taking place in the underdamped regime [3]; e.g., the first observation of Kramers turnover [11] or the emergence of nonconservative instabilities in thermally excited particle motion [12].

Extending levitated optomechanics to systems with multiple mechanical degrees of freedom offers new opportunities in research fields ranging from the quantum entanglement of “macroscopic” objects of increasing mass and size [13] to the statistical mechanics of coupled systems. Mechanisms for entangling levitated nanoparticles via nonlinear coupling have been considered theoretically [14,15], and experiments involving the selective cooling of individual phonon modes and redistribution of energy between nonlinearly coupled oscillators (in analogy with the famed Fermi–Pasta–Ulam–Tsingou problem [16]) can be envisaged [17,18].

When a collection of microobjects is immersed in an incident optical field, multiple light scattering events result in mutually coupled optical forces between the illuminated objects that depend on the configuration of the system as a whole. These optical binding forces [19] then drive self-organization of various micro- and nanoassemblies [19–23], dependent on the shape [24–26], material [27–29], and internal structure [30] of the objects comprising the ensemble, as well as on the polarization [31,32] and spatial distribution [33,34] of the incident optical field. Since the optical binding forces arise due to the deflection of momentum carried by the scattered light, the law of action and reaction between optically bound objects does not need to be satisfied [35]. As a consequence, optical binding forces have nonconservative components [36], which dramatically influence the motion of individual illuminated particles in the underdamped regime [12], while optically bound systems immersed in viscous liquids at low Reynolds numbers are known to display a rich variety of nonequilibrium dynamic behaviors [35,37,38]. Previous experimental studies of optical binding dynamics have most commonly been conducted at overdamped conditions; to date, only two studies have considered the behavior of optically bound matter in underdamped environments [39,40]. In particular, Bykov *et al.* reported long-range optical binding of several levitated microparticles mediated by intermodal scattering and interference inside the evacuated core of a hollow-core photonic crystal fiber [39], while Arita *et al.* demonstrated optical binding between two rotating chiral microparticles confined in a vacuum in independent circularly polarized optical traps [40].

In this paper, we investigate optical binding between multiple silica particles with identical submicron diameters confined in a vacuum using two noninterfering, counter-propagating Gaussian beams with mutually orthogonal linear polarizations. In this highly symmetrical particle–light field configuration, the motion of the illuminated particles is approximately conservative [41]. Our previous studies that were carried out with such an experimental system in the overdamped regime showed the formation of extended arrays of optically bound particles [42], which themselves exhibited exotic optical effects [43]. Here, we concentrate on the collective dynamics of optically bound matter in the underdamped regime at low ambient pressures. We observe the appearance of normal modes of vibration in the thermally driven particle motion, reconstruct the optical binding forces from experimentally recorded particle trajectories, quantify the level of nonlinear coupling between the particles, and investigate the influence of the beam intensity profiles on the particle dynamics. We compare our experimental findings with the predictions of a quantitative theoretical model of particle interaction via optical and hydrodynamic forces, and observe a very good agreement.

## 2. METHODS

#### A. Experimental Setup

To optically confine the studied particles inside a vacuum chamber and characterize their optical binding, we used a source of infrared laser light operating at the vacuum wavelength of 1064 nm with low intensity noise (Mephisto, Coherent Inc., Santa Clara, CA, USA). For the laser beam transformation, we used Thorlabs achromatic doublets with antireflection coating ACN254-XXX-C (L1–L6), dielectric mirrors PF10-03 (M1–M3), and aspheric lenses C240TME-C with antireflection coating (AS1). A collimated Gaussian beam from the laser was first expanded by a telescope formed by lenses L1 (${f_1} = 150\;{\rm mm} $) and L2 (${f_2} = 300\;{\rm mm} $) and projected on a spatial light modulator (SLM) (LCOS X10468-07, Hamamatsu Photonics K.K., Hamamatsu City, Japan). Phase mask encoded at the SLM diffracted the beam into the ${\pm}1$ diffraction orders that were used to generate the two counter-propagating trapping beams; the zeroth and higher orders were blocked by a stop placed in the focal plane of lens L3 (${f_3} = 400\;{\rm mm} $). The two transmitted ${1^{{\rm st}}}$-order beams were reflected from prisms P1 and collimated by lenses L4 (${f_4} = 200\;{\rm mm} $). These lenses formed telescopes with the lens L3, projecting the SLM plane on the mirrors M2. Telescopes consisting of lenses L5 (${f_5} = 100\;{\rm mm} $) and L6 (${f_6} = 150\;{\rm mm} $) then imaged the SLM plane onto the back focal planes of aspheric lenses AS1 ($f = 8\;{\rm mm} $, maximal NA = 0.5) that focused the two trapping beams into the vacuum chamber. The polarization of the beams was controlled with a pair of half-wave plates to form a cross-polarized dual-beam trap [12,33]. The focal planes of the two beams created in the trapping region by aspheric lenses AS1 were slightly displaced from each other along the beam propagation direction $z$ (by approximately $5\;\unicode{x00B5}{\rm m}$; see red lines in Fig. 1) to increase the axial trapping stability [20]. The beam waist radii ${w_0}$ in the trapping region could be adjusted in the range between $1 {-} 3\,\,\unicode{x00B5}{\rm m}$ by controlling the area of the diffraction grating imposed upon the SLM.

Silica (${{\rm SiO}_2}$) particles (mean diameter d = 600 nm, standard deviation of the particle size distribution 12%, Bangs Laboratories Inc., Fishers, IN, USA) were dispersed in isopropyl alcohol and after ${\sim}20$ min sonication of the suspension, droplets containing the particles were sprayed into the trapping region in the vacuum chamber employing an ultrasonic nebulizer (Beurer IH 50).

To enable position tracking of the optically trapped and bound particles, the sample was illuminated by an independent laser beam (Prometheus, vacuum wavelength 532 nm, Coherent Inc., Santa Clara, CA, USA) propagating along the $y$ direction perpendicular to the imaging $x - z$ plane. Large beam waist radius of ${w_0} = 40\,\,\unicode{x00B5}{\rm m}$ and low power (approximately 5 mW at the sample) of the green illuminating beam ensured its negligible contribution to the net optical force acting on the particles. To record the motion of the particles in the $x - z$ plane, a fast CMOS camera (Phantom V611, Vision Research, Inc., Wayne, NJ USA) was oriented along the $y$ direction; its exposure time was set to 1 µs and the frame rate was 300 kHz. Typically, we recorded at least 100,000 frames from the studied optically bound structures to obtain sufficiently long trajectories for the analysis of their motional dynamics.

The off-line tracking of the particle position from the high-speed video recordings was based on the determination of symmetries in the particle images [44]. Briefly, since a spherical particle produces an azimuthally invariant image, we used the shift property of the Fourier transform and looked for the best horizontal and vertical symmetries in the particle image, which provided us with the information about the in-plane $x$ and $z$ coordinates with the spatial resolution of ${\sim}10$ nm and the temporal resolution of ${\sim}3.3\,\,\unicode{x00B5}{\rm s}$.

#### B. Theoretical Model of Optically Bound Structures with Hydrodynamic Interactions

The detailed derivation of our theoretical model is presented in Section 1 of Supplement 1. Here, we summarize the most important elements of the model and its analytical results we directly used in the data processing and analysis.

The system of $N$ identical optically trapped and bound particles subject to hydrodynamic forces is described by the following set of Langevin equations in the frequency domain,

In the absence of hydrodynamic coupling, the thermal motion of a collection of optically bound spheres can be accurately described by conventional normal mode analysis. At this level of description, the normal modes of the system, given by the eigenvectors $\textbf{v}_i^{(0)}$ of the stiffness matrix, $\textbf{K}$, represent a set of $N$ independent harmonic oscillators. Each of these oscillators has a distinct resonant frequency determined by the respective eigenvalue $\lambda _i^K$ of $\textbf{K}$ and experiences identical viscous damping ${\xi _0} = 6\pi{\mu}a$, corresponding to the Stokes viscous drag coefficient of a spherical particle with radius $a$ moving in a fluid with viscosity $\mu$.

The presence of hydrodynamic interactions modifies the form of the normal modes and their damping, so that each quasi-normal mode has now a different effective damping coefficient associated with it. To approximate these effects, we introduce hydrodynamic coupling as a perturbation in the small parameter $\epsilon = a/r$, where $r$ is the mean separation between the spheres, and evaluate all relevant quantities to first order in $\epsilon$. This analysis is described in detail in Supplement 1, Section 1.3. To the first order, the quasi-normal modes are given by the eigenvectors ${\textbf{v}_i} \approx \textbf{v}_i^{(0)} + \epsilon \textbf{v}_i^{(1)}$ of the modified stiffness matrix of the hydrodynamically coupled system, which depends on the friction matrix, and the resonant frequency of mode $i$ is

As described in Section 1.4 of the Supplement 1, the stiffness matrix $\textbf{K}$ for a quasi-linear chain of optically bound spheres is approximated using a simple mass-and-spring type model, with each sphere interacting only with its nearest neighbors. In this approximation, each sphere feels the effects of an external optical trapping potential characterised by a stiffness $k$, and it is also subject to a harmonic binding force, characterized by an effective stiffness $b$. Hydrodynamic coupling between the neighboring spheres is then modeled using the Oseen tensor (see Section 1.5 of Supplement 1). In Tables 1–3, we provide the eigenvalues $\lambda _i^K$ of the stiffness matrix $\textbf{K}$, corresponding eigenvectors $\textbf{v}_i^{(0)}$ of $\textbf{K}$ representing the normal modes of a system without hydrodynamic coupling, first-order corrections $\textbf{v}_i^{(1)}$ to the eigenvectors of the quasi-normal modes of a hydrodynamically coupled system, as well as the correction factors ${c_i}$ that modify the effective modal friction, for 2, 3, and 4 identical interacting spheres. These quantities are expressed in terms of the stiffness coefficients $k$ and $b$, and the Stokes drag coefficient ${\xi _0}$. The case of two spheres can be solved exactly, without the use of perturbation, as shown in Table 1. The only remaining assumption is the use of the Oseen tensor. Specifically, oscillatory dynamics of two elastically and hydrodynamically coupled spheres can be fully characterized by two normal modes: the center-of-mass (COM) mode ($i = 1$), in which the two particles oscillate in phase within the external trapping potential while their separation distance remains constant, and the breathing (BR) mode ($i = 2$), in which the particles oscillate out-of-phase, so that the system’s center of mass does not move but the particle separation distance periodically changes [45]. Tables 2 and 3 show information for three and four spheres, respectively.

#### C. Determination of Interparticle and External Forces and Trapping Stiffness from Particles’ Trajectories

To determine the spatial profiles of interparticle forces from the experimentally recorded particle positions, we first calculated the acceleration of breathing-mode oscillations along the actual mode trajectory, using numerical differentiation of the particle separation distance $\Delta z(t) = [z_2^{\rm a}(t) - z_1^{\rm a}(t)]$, where $z_1^{\mathrm{a}}, z_2^{\mathrm{a}}$ are the absolute coordinates of the particles [see Fig. 2(a)]. Applying Newton’s second law, we then obtained the averaged net interparticle force along the $z$ axis as $f_{z;21}^{{\rm net}}(\Delta z) = m{\langle \Delta \ddot z(t)\rangle _{[\Delta z,\Delta z + \delta z)}}$, where $\Delta \ddot z(t)$ is the breathing-mode acceleration and the operator ${\langle \rangle _{[\Delta z,\Delta z + \delta z)}}$ denotes averaging over all values of its argument observed for particle separation distances lying within the interval $[\Delta z,\Delta z + \delta z)$. This averaging procedure was applied within equally sized adjacent bins of width $\delta z \ll \langle {\Delta z} \rangle$ covering the relevant range of $\Delta z$ ($\langle \Delta z \rangle$ represents the mean separation between the particles), resulting in the profile of $f_{z;21}^{{\rm net}}(\Delta z)$ at a series of discrete values of $\Delta z$. Generally, $f_{z;21}^{{\rm net}}(\Delta z)$ contains deterministic contributions from the optical and viscous forces and a random contribution from the stochastic thermal force [40,45]. However, the averaged modal viscous drag $\xi {\langle \Delta \dot z(t)\rangle _{[\Delta z,\Delta z + \delta z)}}$, where $\xi$ is the coefficient of effective modal viscous friction (see Supplement 1 for details), is zero for all $\Delta z$, as during their motion, the particles are equally likely to be in the configuration with a given value of $\Delta z$ with a positive and negative modal velocity $\Delta \dot z(t)$ of the same magnitude. Similarly, the averaged modal Langevin thermal force ${\langle \Delta f_{z;21}^{\rm L}(\Delta z)\rangle _{[\Delta z,\Delta z + \delta z)}} = 0$ for all $\Delta z$. These assertions are rigorously true for a system in thermal equilibrium; hence, the fact they hold for the experimentally observed modal trajectories indicates we operated at near-equilibrium conditions. Since the viscous and thermal forces average out, $f_{z;21}^{\mathrm{net}}(\Delta z)$ directly represents the all-optical interparticle force $f_{z;21}(\Delta z)$. Using an analogical procedure, we determined the spatial profiles of the external optical trapping force as $f_{z;ext}^{{\rm net}}(\bar z) = m{\langle \ddot {\bar z}(t)\rangle _{[\bar z,\bar z + \delta z)}}$, where $\bar z(t) = [z_2^{\rm a}(t) + z_1^{\rm a}(t)]$ and the operator ${\langle \rangle _{[\bar z,\bar z + \delta z)}}$ denotes averaging over all values of its argument observed for particle separation distances lying within the interval $[\bar z,\bar z + \delta z)$.

The stiffnesses $k$ of the external trapping potential and $B,{B^{{\sin}}}$ of the interparticle potential, presented in Fig. 3(e) as functions of the mean particle separation distance $\langle {\Delta z} \rangle$, were determined from the experimental data by two independent methods:

- 1. The spatial profiles of the interparticle optical forces, $f_{z;21}^{{\rm net}}(\Delta z)$, were fitted with Eq. (8) to obtain the values of $B,{B^{{\sin}}}$ [using also Eq. (9) for ${B^{{\sin}}}$], whereas the values of $k$ were determined from the slopes of the spatial profiles of the external optical trapping forces, $f_{z;ext}^{{\rm net}}(\bar z)$ [symbols $\circ$ in Fig. 3(e)].
- 2. Power spectral densities of experimental modal trajectories were first fitted with Eq. (4) to extract the mode resonant frequencies, $\omega _{{\rm COM}}^{x,z}, \omega _{{\rm BR}}^{x,z}$, and subsequently, the values of $k$ and $B$ were determined using Eqs. (6) and (7) and the known mass of the used particles [symbols $\times$ in Fig. 3(e)].

## 3. RESULTS AND DISCUSSION

#### A. Modification of Optical Trapping Potential due to Light Scattering and Interference

In our optical binding experiments, we used the trapping geometry with two horizontal, counter-propagating laser beams with wave vectors $\textbf{k}$, laterally overlapping in a vacuum chamber so that the beam axes were aligned along the same line [33] [see Figs. 1 and 2(a) for illustration]. As the strength of transverse optical trapping in the $x - y$ plane is significantly higher than the axial trapping strength along the $z$ direction, the trapped particles tend to migrate to the common axis of the two beams and form a quasi one-dimensional (1D) structure oriented along the $z$ axis [20–22]. Figure 2(a) illustrates how the scattered photons modify the distribution of light intensity in the simplest case of two optically trapped particles. For each of the two counter-propagating incident beams, the light scattered in the beam propagation direction denoted by $\textbf{k}$ locally increases its optical intensity, whereas the backscattered wave interferes with the corresponding co-polarized incident beam, creating a weak standing wave along the optical axis [42].

In contrast to previous work on optical levitation of nanoparticles in a stationary, multistable optical trapping potential with an oscillatory standing wave component [46], optically bound nanoparticles in our experiments did not feel any preexisting, static periodic potential landscape. Instead, modulation of the instantaneous trapping potential sensed by a given particle dynamically changed with the changing positions of the other simultaneously illuminated particle (or particles). This dependence of the trapping potential on the positions of individual illuminated particles then resulted in coupling of the particles’ motion and in their collective dynamical response. Since the modulation depth of the optical interference field produced by backscattering was rather low, the thermally driven particles could frequently hop between different trapping positions and, thus, sampled the net anharmonic trapping potential containing multiple potential wells [see inset in Fig. 2(b) for illustration].

#### B. Linearly Coupled Harmonic Oscillators

As discussed in detail in Section 1 of Supplement 1 and summarized in the Methods section of this paper, a system of $N$ optically and hydrodynamically interacting particles moving in an external optical trapping potential can be described as a set of $N$-coupled damped harmonic oscillators. With an appropriate coordinate transformation, this set can be represented by an equivalent system of $N$-uncoupled linear oscillators, whose resonant frequencies define the normal modes of the system. For two coupled particles, oscillatory dynamics can be fully characterized in terms of the COM and BR normal modes. Along the $z$ axis, the transformations from the absolute positions of the particles $z_1^{\rm a},z_2^{\rm a}$ to the normal-mode coordinates ${z_{{\rm COM}}},{z_{{\rm BR}}}$, with the corresponding resonant frequencies $\omega _{{\rm COM}}^z,\omega _{{\rm BR}}^z$, read as (see the Methods section and Sections 1.4 and 1.6 in Supplement 1 for details):

Figures 2(c) and 2(d) show representative sequences of video frames recorded from two trapped, coupled particles in independent experimental runs. Depending on the actual initial conditions (particle position and velocity) at the time of entering the trap, the particles were observed to oscillate either in phase in the COM mode [see Fig. 2(c)] or out of phase in the BR mode [see Fig. 2(d)]. In general, the thermally driven motion of two coupled particles is a linear combination of the COM and BR normal modes with different relative weights. This notion can be illustrated by the power spectral densities (PSDs) of particle positions calculated from the particle trajectories [$x_1^{\rm a}(t),z_1^{\rm a}(t)$] and [$x_2^{\rm a}(t),z_2^{\rm a}(t)$] recorded in the $x - z$ plane. As shown in Fig. 2(e), the position PSDs of both particles 1 and 2 along the $x$ and $z$ directions feature distinct peaks corresponding to $\omega _{{\rm COM}}^x/2\pi$ ($\omega _{{\rm COM}}^z/2\pi$) and $\omega _{{\rm BR}}^x/2\pi$ ($\omega _{{\rm BR}}^z/2\pi$). The differences between the transverse ($x$) and axial ($z$) values of the mode resonant frequencies stem from the different values of $k$ and $b$ along these two directions [22]. Upon carrying out the coordinate transformations indicated in Eqs. (6) and (7), the two normal modes can be separated, as verified by the PSDs of $({x_{{\rm COM}}}, {z_{{\rm COM}}})$ and $({x_{{\rm BR}}}, {z_{{\rm BR}}})$ plotted in Fig. 2(f).

The resonant frequency $\omega _{{\rm BR}}^z$ of the axial breathing mode reflects the overall optical interparticle force characterized by the spring constants $k$ of the external optical potential and $b$ of the optical binding interaction [see Eq. (7)]. To determine the profile of this force from the experimental data, we first calculated the acceleration of breathing-mode oscillations along the actual mode trajectory, using numerical differentiation of the particle separation distance $\Delta z(t) = [z_2^{\rm a}(t) - z_1^{\rm a}(t)]$, and subsequently obtained the averaged net interparticle force along the $z$ axis as $f_{z;21}^{{\rm net}}(\Delta z) = m{\langle \Delta \ddot z(t)\rangle _{[\Delta z,\Delta z + \delta z)}}$ (see the Methods section for details). Since the hydrodynamic and thermal forces average out, $f_{z;21}^{{\rm net}}(\Delta z)$ directly represents the all-optical axial interparticle force $f_{z;21}(\Delta z)$. Figure 2(b) compares the experimentally determined interparticle force, $f_{z;21}(\Delta z)$ (blue line), with its theoretical counterpart $f_{{\rm z;21}}^{{\rm th}}(\Delta z) = ({f_{z;2}} - {f_{z;1}})$ (black line), calculated as the difference of the axial optical forces ${f_{z;2}}, {f_{z;1}}$ acting on the two illuminated, optically bound particles. For the actual experimental values of system parameters, the model based on the multiple-scattering Mie theory [42,47] agrees very well with the experiments. A similar conclusion applies to the experimental (orange line) and theoretical (black line) interparticle potentials ${U_{{\rm z;21}}}(\Delta z)$, calculated by the integration of the corresponding forces as ${U_{{\rm z;21}}}(\Delta z) = - \int {f_{{\rm z;21}}}(\Delta z){\rm d}\Delta z$.

#### C. Non-linearly Coupled Particles

Within the experimentally relevant range of $\Delta z$, the axial optical interparticle potential resembles a harmonic potential well with a superimposed sinusoidal modulation [see the inset in Fig. 2(b)]. Consequently, the associated axial optical interparticle force ${f_{{\rm z;21}}}(\Delta z)$ is a superposition of a linear term and a nonlinear oscillatory component, so

The spatial arrangement and motional dynamics of optically trapped and bound particles strongly depend on the distribution of light intensity of the illuminating beams. As illustrated in Fig. 3(a), the interparticle separation distance can be sensitively and reversibly tuned by changing the beam waist radius ${w_0}$ [34]. This behavior can be intuitively explained in terms of the interplay between the optical binding and trapping forces. Specifically, the particles coupled by optical springs with stiffness $b$ are simultaneously subject to the trapping forces of the external optical potential with stiffness $k$, which tend to push the particles together [20,45]. Hence, the interparticle springs are, in effect, always partially compressed (see Supplement 1, Section 1.4, for details). For larger beam waists, the stiffness of the external trapping potential––and, thus, the magnitude of the trapping forces––is reduced and the interparticle springs expand toward their relaxed length.

In the particular experiment illustrated in Fig. 3(a), the particle separation was varied between 4 and $13\,\,\unicode{x00B5}{\rm m}$ by increasing ${w_0}$ from 1.6 to $2.6\,\,\unicode{x00B5}{\rm m}$. The experimentally determined optical interparticle forces, ${f_{{\rm z;21}}}(\Delta z)$ (colored lines), together with the corresponding fits to Eq. (8) (black lines), are shown in Fig. 3(b) for three different values of ${w_0}$. The corresponding fitted values of $\alpha$ significantly increase with reduced interparticle separation distance, indicating the increasing nonlinearity of the optical binding [42]. Larger modulation amplitudes $\alpha$ then lead to the appearance of multiple stable particle configurations indicated by dotted vertical lines in Figs. 3(b)–3(d). Depending on the value of ${w_0}$, various anharmonic interparticles potentials, ${U_{{\rm z;21}}}(\Delta z)$, were observed, including double wells at small interparticle separations [see yellow line in Fig. 3(d)] and multistable potentials for larger separations [see blue line in Fig. 3(d)]. Due to weak barriers between the local minima in the potential landscape, the particles frequently jumped between multiple stable configurations, with a transition rate comparable to $\omega _{{\rm BR}}^z/2\pi = \sqrt {B/(4{\pi ^2}m)}$ [see Fig. 3(c)].

Figure 3(e) compares the stiffness $k$ of the external trapping potential with the harmonic ($B$) and sinusoidally modulated (${B^{{\rm \sin}}}$) interparticle stiffness for different mean particle separation distances, $\langle {\Delta z} \rangle$. Individual stiffness values were found using two independent procedures. In the first one, the stiffness was calculated from the spatial profiles of the interparticle force ($B,{B^{{\sin}}}$) and external trapping force ($k$) acting on the system of optically bound particles; these force profiles were determined from the analysis of experimental particle trajectories (see the discussion above and in the Methods section). Alternatively, the stiffness was estimated from the position of the characteristic spectral peaks in the PSDs of normal-mode coordinates and known particle mass (see Methods section). As illustrated in Fig. 3(e), the values of $k$ and $B$ obtained from the two procedures agree very well for both the $x$ and $z$ directions. Due to the small oscillation amplitude of particles confined in the secondary optical traps in multistable potentials and broadening of spectral peaks caused by anharmonicity of these potentials and relatively large viscous damping, we were not able to observe clear evidence of oscillatory peaks corresponding to $\omega _{{\rm BR}}^{z,\sin} = \sqrt {{B^{{\rm \sin}}}/m}$ in the modal PSDs (see Supplement 1, Section 2, for details).

#### D. Pressure Dependence and Collision of Particles

In the previous studies of the motion of individual optically levitated particles [8,9], it was observed the traps became unstable at pressures around 1 mbar, so an active feedback control scheme was required for trapping at such lower pressures. To study the dynamics of optically bound levitated particles under these ambient conditions, we gradually reduced the pressure from atmospheric level down to the critical value and continuously monitored the particles’ response. In doing so, we took advantage of our position detection method based on an ultrafast CMOS camera, which gave us quantitative, properly calibrated information about the amplitude of particles’ motion without any prior knowledge of the particles’ properties [12].

In general, ambient pressure affects the profile of internal temperature within the levitated particles, which subsequently modifies their density and refractive index and influences the optical forces acting on the particles [49]. This effect was qualitatively studied using a theoretical model presented in Section 3 of the Supplement 1. However, by analysing particle trajectories recorded at different pressures, we found that the optical interparticle forces and potentials were almost independent of the pressure [see Figs. 4(a) and 4(b)]. Consequently, any changes in the optical properties of the particles were rather small. We conclude that the observed increase in the amplitude of particles’ motion was primarily caused by increasing the kinetic energy of the ambient gas heated due to the absorption of the trapping laser light in the levitated particles [49–51].

Moreover, silica particles are porous and contain traces of condensed water, which, when heated, can evaporate, strongly affecting the particles’ motion [52]. Using the width of the probability density functions for velocity components $\nu$, ${\rm PDF}(\nu) \propto \exp [- m{{\nu}^2}/(2{k_{\rm B}}T)]$, where $m$ is the mass of the particle and ${k_{\rm B}}$ is the Boltzmann constant [see examples in Fig. 4(c)], we estimated the effective temperatures $T$ of the axial and transverse BR and COM normal modes of optically bound particles [Fig. 4(d)]. We observed different effective temperatures for the particle motion along and perpendicular to the optical axis. We attribute this difference, at least in part, to anisotropic absorption of the incident laser light within the sphere volume [50]. While reducing the ambient pressure, the amplitudes of particle oscillations were observed to increase and, at a pressure around 1 mbar, the particles were observed to collide, as shown in Fig. 4(e).

To investigate the possible effects of nonconservative optical forces on the stability of optically levitated [53] and bound [25,35] structures, we carried out full 3D stochastic simulations of the system, using the multiple Mie scattering theory to correctly describe the optical forces acting on the system (see Section 4 of the Supplement 1). These simulations did not show any evidence of instability in the optically bound structures for pressures down to ${10^{- 7}}\;{\rm mbar}$ and incident powers corresponding to the typical experimental values. Thus, nonconservative optical forces can be ruled out as the main source of the observed instability.

#### E. Optically Bound Structures with More Than Two Particles

With a growing number of constituent particles, oscillatory dynamics of optically bound structures becomes significantly more intricate, containing contributions from multiple normal modes of increasing complexity. Figure 5 summarizes the analysis of motion of three identical optically bound particles. We calculated the pairwise optical interparticle forces from the time-averaged accelerations of the particles as ${f_{{\rm z;21}}}(\Delta z) = {m}\langle \ddot z_2^{\rm a} - \ddot z_1^{\rm a}\rangle$ and ${f_{{\rm z;32}}}(\Delta z) = {m}\langle \ddot z_3^{\rm a} - \ddot z_2^{\rm a}\rangle$, where $\Delta z = (z_2^{\rm a} - z_1^{\rm a}),\Delta z = (z_3^{\rm a} - z_2^{\rm a})$ are the respective interparticle distances. According to the plots presented in Fig. 5(a), the profiles of ${f_{{\rm z;21}}}(\Delta z)$ and ${f_{{\rm z;32}}}(\Delta z)$ are almost identical, indicating a high degree of symmetry of the studied experimental configuration.

As predicted by the theoretical model of coupled, damped harmonic oscillators, which is developed and discussed in detail in Supplement 1 and summarized in the Methods section, the general motion of the particles along both axial $z$ direction and transverse $x$ direction can be decomposed into three normal modes. In particular, the particles display the COM mode, in which they all oscillate in phase with identical amplitudes and the frequency ${\omega _{{\rm COM}}} \approx \sqrt { k/m}$, symmetric breathing (s-BR) mode, in which the middle particle remains stationary and the peripheral particles oscillate out of phase with identical amplitudes and the frequency ${\omega _{{\rm s} - {\rm BR}}} \approx \sqrt {(k + b)/m}$, and antisymmetric breathing (a-BR) mode, in which the peripheral particles oscillate in phase with identical amplitudes, while the middle particle oscillates out of phase with them with the frequency ${\omega _{{\rm a} - {\rm BR}}} \approx \sqrt {(k + 3b)/m}$. After carrying out the coordinate transformations described in Supplement 1 and in the Methods section of this paper, the PSDs of the normal-mode coordinates can be calculated. Plots of these PSDs shown in Figs. 5(b) and 5(c) for the axial and transverse normal modes, respectively, clearly indicate distinct resonant frequencies for the three modes. The observed differences between the transverse and axial values of ${\omega _{{\rm COM}}}$, ${\omega _{{\rm s} - {\rm BR}}}$, and ${\omega _{{\rm a} - {\rm BR}}}$ result from different values of the trapping and binding stiffness along the two directions [22]. To compare the experimentally observed resonant frequencies with the model predictions, we determined the stiffness $k$ of the external optical trapping potential and the binding stiffness $b$ from ${\omega _{{\rm COM}}}$ and ${\omega _{{\rm s} - {\rm BR}}}$ measured along both $x$ and $z$ directions, and, subsequently, used these values of $k$ and $b$ to calculate the expected ${\omega _{{\rm a} - {\rm BR}}}$. Comparison of thus calculated ${\omega _{{\rm a} - {\rm BR}}}$ with the directly measured value obtained from the corresponding modal PSD revealed an acceptable 10 % difference.

When working in overdamped conditions in aqueous environments, it is possible to experimentally study the dynamics of formation and motion of large assemblies of optically bound particles over extended periods of time [42]. On the other hand, at low ambient pressures associated with low viscous damping, we often observed collisions between individual trapped particles and irreversible formation of small particle clusters, which typically occurred during the particle loading process or during the evacuation of the sample chamber (see Supplement 1 for examples). Despite these experimental challenges, we were still able to form extended chains containing as many as 10 optically bound particles at atmospheric pressure. Usually, instabilities of the bound structures arose when the ambient pressure was reduced, resulting in the loss of the particles from the trap. However, we succeeded in resolving the normal modes of a chain of four optically bound particles confined at the pressure of 20 mbar (see Supplement 1, Section 5, for illustration).

## 4. CONCLUSION

In summary, we demonstrated stable optical binding of up to four microscopic silica particles confined in the underdamped regime, at pressures lower than 20 mbar, using cross-polarized counter-propagating Gaussian beams. We then presented a methodology that determines the parameters of such coupled oscillators––in particular, the trapping and binding stiffness––from the camera record of particles’ trajectories. We were able to clearly resolve collective vibrational modes of optically bound structures characteristic of linearly coupled harmonic oscillators moving in the presence of weak viscous damping and hydrodynamic coupling. In addition, we observed and characterized sinusoidal modulations of the optical binding potentials describing a weak nonlinear interparticle coupling, which is caused by interference of the incident and backscattered light [42] and results in multistability of such optically bound matter [48]. A similar type of externally controllable coupling nonlinearity could facilitate optically induced entanglement between the trapped mesoscopic particles [14,15], and is of significant interest for the systematic studies of energy redistribution between coupled nonlinear oscillators [17,18]. Currently, the lowest achievable gas pressure at which stable optically bound structures can be formed and maintained is limited to ${\sim}1\,\,{\rm mbar}$ by heating effects, which substantially increase the center-of-mass motion of the particles, causing unwanted collisions and instabilities. However, by implementing a parametric feedback cooling scheme [8,9] or cavity-assisted cooling [7], both recently developed for individual optically levitated nanoparticles, it is feasible to selectively cool the collective vibrational modes of the structures, possibly all the way down to the quantum ground state. With efficient cooling protocols in place, the number of coupled degrees of freedom of the system could potentially be further expanded, which represents an important step toward the development of a mesoscopic experimental model of solid-state matter.

## Funding

Grantová Agentura České Republiky (GA18-27546S); European Commission (CZ.1.05/2.1.00/01.0017); Akademie Věd České Republiky (RVO:68081731); Ministerstvo Školství, Mládeže a Tělovýchovy (LO1212); European Regional Development Fund (CZ.02.1.01/0.0/0.0/16_026/0008460).

## Disclosures

The authors declare no conflicts of interest.

## Supplemental document

See Supplement 1 for supporting content.

## REFERENCES

**1. **A. Ashkin and J. Dziedzic, “Optical levitation in high vacuum,” Appl. Phys. Lett. **28**, 333–335 (1976). [CrossRef]

**2. **Z.-Q. Yin, A. A. Geraci, and T. Li, “Optomechanics of levitated dielectric particles,” Int. J.Mod. Phys. B **27**, 1330018 (2013). [CrossRef]

**3. **J. Gieseler and J. Millen, “Levitated nanoparticles for microscopic thermodynamics—a review,” Entropy **20**, 326 (2018). [CrossRef]

**4. **D. Chang, C. Regal, S. Papp, D. Wilson, J. Ye, O. Painter, H. Kimble, and P. Zoller, “Cavity opto-mechanics using an optically levitated nanosphere,” Proc. Natl. Acad. Sci. USA **107**, 1005–1010 (2010). [CrossRef]

**5. **P. F. Barker and M. N. Shneider, “Cavity cooling of an optically trapped nanoparticle,” Phys. Rev. A **81**, 023826 (2010). [CrossRef]

**6. **O. Romero-Isart, M. Juan, R. Quidant, and J. I. Cirac, “Toward quantum superposition of living organisms,” New J. Phys. **12**, 033015 (2010). [CrossRef]

**7. **U. Delić, M. Reisenbauer, K. Dare, D. Grass, V. Vuletić, N. Kiesel, and M. Aspelmeyer, “Cooling of a levitated nanoparticle to the motional quantum ground state,” Science **367**, 892–895 (2020). [CrossRef]

**8. **T. Li, S. Kheifets, and M. G. Raizen, “Millikelvin cooling of an optically trapped microsphere in vacuum,” Nat. Phys. **7**, 527–530 (2011). [CrossRef]

**9. **J. Gieseler, B. Deutsch, R. Quidant, and L. Novotny, “Subkelvin parametric feedback cooling of a laser-trapped nanoparticle,” Phys. Rev. Lett. **109**, 103603 (2012). [CrossRef]

**10. **F. Tebbenjohanns, M. Frimmer, V. Jain, D. Windey, and L. Novotny, “Motional sideband asymmetry of a nanoparticle optically levitated in free space,” Phys. Rev. Lett. **124**, 013603 (2020). [CrossRef]

**11. **L. Rondin, J. Gieseler, F. Ricci, R. Quidant, C. Dellago, and L. Novotny, “Direct measurement of Kramers turnover with a levitated nanoparticle,” Nat. Nanotechnol. **12**, 1130–1133 (2017). [CrossRef]

**12. **V. Svak, O. Brzobohatý, M. Šiler, P. Jákl, J. Kaňka, P. Zemánek, and S. H. Simpson, “Transverse spin forces and non-equilibrium particle dynamics in a circularly polarized vacuum optical trap,” Nat. Commun. **9**, 5453 (2018). [CrossRef]

**13. **M. Arndt and K. Hornberger, “Testing the limits of quantum mechanical superpositions,” Nat. Phys. **10**, 271–277 (2014). [CrossRef]

**14. **H. Nguyena and F. Bernards, “Entanglement dynamics of two mesoscopic objects with gravitational interaction,” Eur. Phys. J. D **74**, 69 (2020). [CrossRef]

**15. **S. Qvarfort, S. Bose, and A. Serafini, “Mesoscopic entanglement through central-potential interactions,” J. Phys. B **53**, 235501 (2020). [CrossRef]

**16. **E. Fermi, P. Pasta, S. Ulam, and M. Tsingou, “Studies of the nonlinear problems,” Technical Report LA-1940 (Los Alamos Scientific Lab, 1955).

**17. **J. Ford, “The Fermi-Pasta-Ulam problem: paradox turns discovery,” Phys. Rep. **213**, 271–310 (1992). [CrossRef]

**18. **R. Livi, M. Pettini, S. Ruffo, M. Sparpaglione, and A. Vulpiani, “Equipartition threshold in nonlinear large Hamiltonian systems: the Fermi-Pasta-Ulam model,” Phys. Rev. A **31**, 1039–1045 (1985). [CrossRef]

**19. **M. M. Burns, J.-M. Fournier, and J. A. Golovchenko, “Optical matter: crystallization and binding in intense optical fields,” Science **249**, 749–754 (1990). [CrossRef]

**20. **S. A. Tatarkova, A. E. Carruthers, and K. Dholakia, “One-dimensional optically bound arrays of microscopic particles,” Phys. Rev. Lett. **89**, 283901 (2002). [CrossRef]

**21. **W. Singer, M. Frick, S. Bernet, and M. Ritsch-Marte, “Self-organized array of regularly spaced microbeads in a fiber-optical trap,” J. Opt. Soc. Am. B **20**, 1568–1574 (2003). [CrossRef]

**22. **K. Dholakia and P. Zemánek, “Gripped by light: optical binding,” Rev. Mod. Phys. **82**, 1767–1791 (2010). [CrossRef]

**23. **T. Čižmár, L. C. D. Romero, K. Dholakia, and D. L. Andrews, “Multiple optical trapping and binding: new routes to self-assembly,” J. Phys. B **43**, 102001 (2010). [CrossRef]

**24. **O. Brzobohatý, A. V. Arzola, M. Šiler, L. Chvátal, P. Jákl, S. Simpson, and P. Zemánek, “Complex rotational dynamics of multiple spheroidal particles in a circularly polarized, dual beam trap,” Opt. Express **23**, 7273–7287 (2015). [CrossRef]

**25. **S. H. Simpson, P. Zemánek, O. M. Maragò, P. H. Jones, and S. Hanna, “Optical binding of nanowires,” Nano Lett. **17**, 3485–3492 (2017). [CrossRef]

**26. **M. G. Donato, O. Brzobohatý, S. H. Simpson, A. Irrera, A. A. Leonardi, M. J. Lo Faro, V. Svak, O. M. Maragò, and P. Zemánek, “Optical trapping, optical binding, and rotational dynamics of silicon nanowires in counter-propagating beams,” Nano Lett. **19**, 342–352 (2019). [CrossRef]

**27. **V. Demergis and E.-L. Florin, “Ultrastrong optical binding of metallic nanoparticles,” Nano Lett. **12**, 5756–5760 (2012). [CrossRef]

**28. **Z. Yan, S. K. Gray, and N. F. Scherer, “Potential energy surfaces and reaction pathways for light-mediated self-organization of metal nanoparticle clusters,” Nat. Commun. **5**, 3751 (2014). [CrossRef]

**29. **Y. Li, H. Xin, Y. Zhang, H. Lei, T. Zhang, H. Ye, J. J. Saenz, C.-W. Qiu, and B. Li, “Living nanospear for near-field optical probing,” ACS Nano **12**, 10703–10711 (2018). [CrossRef]

**30. **O. Brzobohatý, R. J. Hernàndez, S. Simpson, A. Mazzulla, G. Cipparrone, and P. Zemánek, “Chiral particles in the dual-beam optical trap,” Opt. Express **24**, 26382–26391 (2016). [CrossRef]

**31. **C. D. Mellor, T. A. Fennerty, and C. D. Bain, “Polarization effects in optically bound particle arrays,” Opt. Express **14**, 10079–10088 (2006). [CrossRef]

**32. **F. Nan and Z. Yan, “Probing spatiotemporal stability of optical matter by polarization modulation,” Nano Lett. **18**, 1396–1401 (2018). [CrossRef]

**33. **T. Čižmár, O. Brzobohatý, K. Dholakia, and P. Zemánek, “The holographic optical micro-manipulation system based on counter-propagating beams,” Laser Phys. Lett. **8**, 50–56 (2011). [CrossRef]

**34. **O. Brzobohatý, V. Karásek, T. Čižmár, and P. Zemánek, “Dynamic size tuning of multidimensional optically bound matter,” Appl. Phys. Lett. **99**, 101105 (2011). [CrossRef]

**35. **S. Sukhov, A. Shalin, D. Haefner, and A. Dogariu, “Actio et reactio in optical binding,” Opt. Express **23**, 247–252 (2015). [CrossRef]

**36. **S. Sukhov and A. Dogariu, “Non-conservative optical forces,” Rep. Prog. Phys. **80**, 112001 (2017). [CrossRef]

**37. **V. Karásek, M. Šiler, O. Brzobohatý, and P. Zemánek, “Dynamics of an optically bound structure made of particles of unequal sizes,” Opt. Lett. **42**, 1436–1439 (2017). [CrossRef]

**38. **L. Chvátal, O. Brzobohatý, and P. Zemánek, “Binding of a pair of Au nanoparticles in a wide Gaussian standing wave,” Opt. Rev. **22**, 157–161 (2015). [CrossRef]

**39. **D. S. Bykov, S. Xie, R. Zeltner, A. Machnev, G. K. L. Wong, T. G. Euser, and P. St.J. Russell, “Long-range optical trapping and binding of microparticles in hollow-core photonic crystal fibre,” Light Sci. Appl. **7**, 22 (2018). [CrossRef]

**40. **Y. Arita, E. M. Wright, and K. Dholakia, “Optical binding of two cooled micro-gyroscopes levitated in vacuum,” Optica **5**, 910–917 (2018). [CrossRef]

**41. **S. Divitt, L. Rondin, and L. Novotny, “Cancellation of non-conservative scattering forces in optical traps by counter-propagating beams,” Opt. Lett. **40**, 1900–1903 (2015). [CrossRef]

**42. **O. Brzobohatý, L. Chvátal, and P. Zemánek, “Optomechanical properties of optically self-arranged colloidal waveguides,” Opt. Lett. **44**, 707–710 (2019). [CrossRef]

**43. **O. Brzobohatý, L. Chvátal, A. Jonáš, M. Šiler, J. Kaňka, J. Ježek, and P. Zemánek, “Tunable soft-matter optofluidic waveguides assembled by light,” ACS Photon. **6**, 403–410 (2019). [CrossRef]

**44. **I. T. Leite, S. Turtaev, X. Jiang, M. Šiler, A. Cuschieri, P. St.J. Russell, and T. Čižmár, “Three–dimensional holographic optical manipulation through a high–numerical–aperture soft–glass multimode fibre,” Nat. Photonics **12**, 33–39 (2018). [CrossRef]

**45. **B. R. Slezak and B. D’Urso, “A microsphere molecule: the interaction of two charged microspheres in a magneto-gravitational trap,” Appl. Phys. Lett. **114**, 244102 (2019). [CrossRef]

**46. **R. Diehl, E. Hebestreit, R. Reimann, F. Tebbenjohanns, M. Frimmer, and L. Novotny, “Optical levitation and feedback cooling of a nanoparticle at subwavelength distances from a membrane,” Phys. Rev. A **98**, 013851 (2018). [CrossRef]

**47. **D. Mackowski, *The Extension of Mie Theory to Multiple Spheres* (Springer, 2012), pp. 223–256.

**48. **V. Karásek, T. Čižmár, O. Brzobohatý, P. Zemánek, V. Garcés-Chávez, and K. Dholakia, “Long-range one-dimensional longitudinal optical binding,” Phys. Rev. Lett. **101**, 143601 (2008). [CrossRef]

**49. **E. Hebestreit, R. Reimann, M. Frimmer, and L. Novotny, “Measuring the internal temperature of a levitated nanoparticle in high vacuum,” Phys. Rev. A **97**, 043803 (2018). [CrossRef]

**50. **J. Millen, T. Deesuwan, P. Barker, and J. Anders, “Nanoscale temperature measurements using non-equilibrium Brownian dynamics of a levitated nanosphere,” Nat. Nanotechnol. **9**, 425–429 (2014). [CrossRef]

**51. **J. Millen, T. S. Monteiro, R. Pettit, and A. N. Vamivakas, “Optomechanics with levitated particles,” Rep. Prog. Phys. **83**, 026401 (2020). [CrossRef]

**52. **F. Ricci, “Levitodynamics toward force nano-sensors in vacuum,” Ph.D. thesis (Institut de Ciencies Fotoniques, 2019).

**53. **Y. Roichman, B. Sun, A. Stolarski, and D. G. Grier, “Influence of nonconservative optical forces on the dynamics of optically trapped colloidal spheres: the fountain of probability,” Phys. Rev. Lett. **101**, 128301 (2008). [CrossRef]