Brownian dynamics simulations of sphere clusters in optical tweezers

Computationally modeling the behavior of wavelength-sized non-spherical particles in optical tweezers can give insight into the existence and stability of trapping equilibria as well as the optical manipulation of such particles more broadly. Here, we report Brownian dynamics simulations of non-spherical particles that account for detailed optical, hydrodynamic, and thermal interactions. We use a $T$-matrix formalism to calculate the optical forces and torques exerted by focused laser beams on clusters of wavelength-sized spheres, and we incorporate detailed diffusion tensors that capture the anisotropic Brownian motion of the clusters. For two-sphere clusters whose size is comparable to or larger than the wavelength, we observe photokinetic effects in elliptically-polarized beams. We also demonstrate that multiple trapping equilibria exist for a highly asymmetric chiral cluster of seven spheres. Our simulations may lead to practical suggestions for optical trapping and manipulation as well as a deeper understanding of the underlying physics.


Introduction
Modeling the behavior of wavelength-sized particles in optical tweezers can guide experiments in optical manipulation or assembly.Doing so requires calculating the optical forces and torques exerted on particles in a focused laser beam.These calculations are particularly challenging for wavelength-sized particles since neither the Rayleigh approximation nor ray optics can be validly used.Rather, detailed consideration of the electromagnetic fields incident on and scattered by the particle is necessary [1].T-matrices provide a well-established formalism for computing and describing electromagnetic scattering by non-spherical particles [2][3][4].To date, however, the use of T-matrices to calculate forces and torques in optical tweezers has mainly been limited to rotationally symmetric particles or particles that are smaller than the wavelength [5][6][7][8][9][10][11].
Even so, being able to compute optical forces and torques is not sufficient for understanding the behavior of non-spherical particles in optical tweezers.Key questions include the positions, orientations, and stability of any trapping equilibria that may exist.But as Bui et al. point out, finding equilibria solely from calculations of forces and torques at fixed positions and orientations is difficult for non-spherical particles because of the size of the parameter space [12].Consider for example a sphere in a circularly polarized beam.Rotational symmetry requires that the sphere center lie on the beam axis in equilibrium.Moreover, the sphere's orientation is irrelevant.So, to find the equilibrium trapping position, it is only necessary to calculate the axial component of the trapping force as a function of the sphere's position along the beam axis.No such simplifications are possible for particles lacking symmetry; mapping the trapping landscape for an asymmetric particle may require considering 3 positional and 3 orientational degrees of freedom [12].
An alternative strategy for finding equilibria is to perform dynamical simulations that use computations of optical forces and torques to model the actual motion of a particle in optical tweezers.If the goal is only to find trapping equilibria, it is not necessary to carefully consider other interactions that a trapped particle may experience but that vanish in static equilibrium, such as hydrodynamic resistance.In contrast, investigating the stability of any trapping equilibria or the dynamics of transitions between multiple equilibria requires more care.Optically-trapped particles suspended in a fluid experience Brownian motion involving both hydrodynamic interactions and random thermal interactions.These can both be highly anisotropic.Detailed consideration of Brownian motion is also needed to explore photokinetic effects arising from the angular momentum carried by an elliptically-or circularly-polarized beam since they can involve an interplay between optical and thermal interactions [13,14].
Here, we report the first Brownian dynamics simulations for wavelength-sized non-spherical particles in optical tweezers that account for optical interactions, hydrodynamic interactions, and thermal fluctuations in detail.We combine mstm, a well-established package for calculating T-matrices of clusters of spheres [4,15] with ott, a package that can calculate the optical forces and torques on a particle in a focused laser beam given a T-matrix [8,16].Furthermore, we use complete diffusion tensors for the clusters in order to realistically predict their anisotropic behavior in a fluid.Our work not only extends the size regime for which optical forces and torques on sphere clusters have previously been reported but also predicts their detailed motion in optical tweezers for the first time.While dynamical simulations of particles in optical tweezers have previously been reported [10,12,13,[17][18][19], our simulations incorporate both thermal fluctuations and anisotropic hydrodynamic resistance, and we do not linearize the optical forces and torques.We validate our simulations by considering single spheres and then explore the remarkably rich dynamical effects that occur for clusters of two spheres that are each comparable to or larger than the wavelength of the trapping laser.We then use our simulations to find multiple trapping equilibria for a highly asymmetric cluster of seven spheres.

Simulating Brownian motion with external interactions
Fig. 1.Coordinate systems used in simulations.In the laboratory coordinate system, a laser beam focused at (0, 0, 0) propagates in the +z direction.In the particle coordinate system, the spheres of a given cluster are located at fixed reference positions.The origin of the particle coordinate system is located at the cluster's center of mass.The orientation of the particle axes in the laboratory coordinate system describes the cluster's orientation.
We seek to simulate the trajectory of an arbitrary Brownian particle at temperature T in a fluid of viscosity η that is subject to external forces and torques.We assume that our particles are rigid and move in what we call the "laboratory frame."Consequently, we can describe the position and orientation of our particles by giving the laboratory coordinates of their center of mass and a rotation matrix.Specifically, for clusters of spheres, we define a reference configuration in which the position of each sphere is defined relative to the particle coordinate axes u 1 , u 2 , and u 3 , which we take to be unit vectors.We choose the origin of the particle coordinate axes to lie at the cluster center of mass (CM).The components of the particle coordinate axes in the laboratory frame are given by the columns of a 3 × 3 rotation matrix.Figure 1 illustrates the laboratory and particle coordinate systems.
In the overdamped limit in which the particle's inertia is negligible, we model the trajectory of a Brownian particle using a finite-difference approach by computing its generalized displacements during time steps of duration ∆t [18,20].Following Fernandes and García de la Torre [20], we consider the 6-component generalized displacement vector ∆q whose transpose is given by ∆q tr = (∆x, ∆y, ∆z, ∆φ 1 , ∆φ 2 , ∆φ 3 ) . (1) Here, ∆x, ∆y, and ∆z are displacements of the center of mass and ∆φ 1 , ∆φ 2 , and ∆φ 3 are infinitesimal rotations about the particle's axes u 1 , u 2 , and u 3 .The generalized displacement ∆q is computed in the particle coordinate system and subsequently needs to be transformed to laboratory coordinates.To find ∆q for a given time step, we separately consider two contributions: The first term is due to the particle's random Brownian motion, while the second term is due to the external interactions.The Brownian displacement ∆q B with components ∆q B i is simulated by generating normally distributed random numbers with covariance where the D i j are the elements of the particle's 6 × 6 diffusion tensor D [20].The displacement ∆q ext due to an external force F ext and torque n ext is given by where F ext and n ext have been combined into the 6-component vector F [20].Our main computational challenges are computing the external forces and torques F and finding the diffusion tensor D.

Calculating optical forces and torques with T-matrices
While the approach in Eqs.(3) and ( 4) is relevant to any external interaction, we now focus on the optical forces and torques exerted by optical tweezers.We consider a laser beam of vacuum wavelength λ 0 in a medium of refractive index n med that is focused by an objective lens with a given numerical aperture (NA).
Calculating the optical forces and torques is fundamentally a light scattering problem.As we briefly summarize here, we expand the position-dependent electric field of the incident beam E inc (kr) in vector spherical wave functions (VSWFs) [21]: Here, RgM nm (kr) and RgN nm (kr) are solutions of the vector Helmholtz equation that are finite at the origin and k is the beam's wavenumber.(See Nieminen et al. and references therein for details [21].)We use ott to truncate the expansion and calculate the expansion coefficients a nm and b nm for a focused Gaussian beam using an overdetermined least squares fit [8,22].The fitting procedure is necessary because a Gaussian beam is not a solution to the vector Helmholtz equation.
The M (1)  nm and N (1)  nm VSWFs asymptotically behave as outgoing spherical waves.The details of the scattering process are described by the T-matrix T, which relates the scattered coefficients to the incident coefficients: where we combine the incident and scattered field coefficients into single column vectors.ott can natively calculate T-matrices for homogenous spheres, whose elements are the Lorenz-Mie coefficients.However, ott cannot calculate T-matrices for sphere clusters.Instead, we use the multisphere superposition code mstm developed by Mackowski and Mishchenko [15].Our code calls mstm, reads the T-matrix it generates, and then uses ott to perform the matrix multiplication in Eq. ( 7) to find the scattered field expansion coefficients p mn and q mn .Subsequently, ott determines the optical forces and torques from the field expansion coefficients a nm , b nm , p nm , and q nm [8].In this approach, the computationally challenging light scattering problem only needs to be solved once in a general way when computing the T-matrix; subsequent force and torque calculations are fast.

Finding diffusion tensors
Brownian dynamics simulations require knowing the diffusion tensor D of a particle.Finding D requires solving the Stokes equations for creeping flow around the particle.For spheres of radius a, D is diagonal with two unique, nonzero elements: the translational diffusion constant D t = k B T/(6πηa) and the rotational diffusion constant D r = k B T/(8πηa 3 ).For clusters of two identical spheres, or dimers, D is also diagonal and has four unique elements.D t, describes translational diffusion along the dimer's long axis while D t,⊥ describes translational diffusion perpendicular to the long axis.Similarly, D r, describes rotational diffusion about the long axis and D r,⊥ describes rotational diffusion about the other two axes.We use the analytic solution of Nir & Acrivos to find D for dimers [23].However, for more complex sphere clusters, numerical methods are needed.We use the Fortran program BEST, which implements a boundary element solution to the Stokes equations over an arbitrary triangulated surface [24].As recommended, we perform repeated BEST calculations using increasing numbers of triangles and extrapolate to the limit of an infinite number of triangles.By comparing the results of BEST for single spheres and dimers to the analytic results, we estimate that using BEST results in uncertainties in the diffusion tensor elements of no more than 1%.

Implementation
We implement our Brownian dynamics algorithms in the package brownian_ot, which is freely available on Github [25].Our code integrates both a Matlab package (ott) and a Fortran 90 package (mstm) in a single interface.To make this interfacing straightforward, we implement brownian_ot in Python.The initial steps in simulating a given particle in a beam of a specified wavelength are determining the particle's diffusion tensor, determining the particle's T-matrix for scattering of light of that wavelength, and determining the beam's VSWF expansion coefficients.These calculations only need to be done once.Thereafter, we simulate the trajectory of a particle beginning from an arbitrary initial position and orientation as follows: 1. Use ott to compute the optical force and torque F on the particle in the laboratory coordinate system.
2. Correct the optical torque to be relative to the particle's center of diffusion [26] and transform all elements of F to the particle coordinate system.
3. Use Eq. ( 4) to calculate ∆q ext in the particle coordinate system.
4. Use Eq. ( 3) to calculate ∆q B in the particle coordinate system.
6. Transform the CM displacement (the spatial components of ∆q) to the laboratory coordinate system and update the CM position.
7. Calculate an unbiased rotation corresponding to the angular components of ∆q [27] and update the particle's orientation by rotation composition.

Repeat.
Internally, our orientation calculations use unit quaternions due to their compact storage and numerical stability.While our code runs on personal computers under Macintosh, Windows, or Linux operating systems, we typically run longer simulations using the Comet cluster at the San Diego Supercomputing Center, which is part of the Extreme Science and Engineering Discovery Environment (XSEDE) [28].The primary reason to use Comet is that we can run tens of independent simulations simultaneously, either to build up statistics or to explore parameter space.Run times depend primarily on the number of time steps to be simulated but also on the particle size since larger particles generally require a larger maximum order n max in the truncation of the VSWF expansions for the incident and scattered fields.On Comet, simulating 3 × 10 6 steps for a 1-µm-diameter sphere takes approximately 8 × 10 4 s, and simulating a trajectory of the same length for a cluster of seven 0.8-µm-diameter spheres takes approximately 1 × 10 5 s.

Simulation parameters
Throughout what follows, we consider a 1064-nm-wavelength beam focused by a 1.2 NA objective lens.Except when noted, the beam has a power of 5 mW and is left-circularly-polarized (LCP) with Jones vector (1, i).We also assume that the beam propagates in an isotropic medium with n med = 1.33, which corresponds to water.We assume that the medium has viscosity η = 1 × 10 −3 Pa s and that the particles are at T = 295 K.We consider non-absorbing particles with two different refractive indices: silica (Si; n = 1.45) and polystyrene (PS; n = 1.58).However, our code can be easily applied to particles with other refractive indices, including absorbing particles such as metallic nanospheres.
We choose a simulation time step ∆t = 1 × 10 −5 s throughout our simulations.Because optical forces and torques are strongly dependent on particle position and orientation, we choose the time step to be small enough that none of our particles translates a significant fraction of the wavelength or rotates through a significant angle during a single step.At the same time, ∆t is large enough that ballistic motion is negligible [18,29].Finally, except when noted, all simulations begin with the particle centers of mass at the origin and the particle reference axes aligned with the laboratory frame axes.

Spheres
To validate our simulations, we consider the behavior of a 1-µm-diameter PS sphere.The optical forces as a function of position are approximately linear near equilibrium (Fig. 2(a)):  9) and (10).
where κ x and κ z are the stiffnesses in the x and z directions, respectively.Because of radiation pressure, F z = 0 at z = z eq with z eq > 0 rather than at z = 0.The first 2 s of a simulated sphere trajectory shows that the sphere exhibits Brownian fluctuations about equilibrium (Fig. 2(b)).Mean-squared displacements (MSDs) calculated from the simulated trajectories demonstrate that the simulations perform as expected.The MSD for the coordinate x i is defined as MSD = [x i (t + τ) − x i (t)] 2 .We plot the x and z MSDs calculated from our simulations in Fig. 2(c).Since we know the stiffnesses κ x and κ z from the force calculations in Fig. 2(a) and also know the sphere radius, we can also predict the MSDs a priori.Solving the Langevin equation shows that a sphere of radius a moving in one dimension subject to a linear restoring force with stiffness κ has the MSD [ where the friction coefficient is given by γ = 6πηa.The dotted lines in the upper panel of Fig. 2(c) show that the predictions of Eq. ( 9) agree well with the values calculated from the trajectories.
In addition, because there are no significant optical torques, the sphere should undergo free rotational diffusion.We can similarly characterize the rotational dynamics of the sphere by calculating the mean-squared angular displacement (MSAD) of any of the reference axes in the laboratory frame.The MSAD for axis u i is defined as The lower panel of Fig. 2(c) shows the MSAD for u 1 since, for a sphere, all three axes are equivalent.A sphere with rotational diffusion constant D r has an MSAD given by [30,31] The agreement between the prediction of Eq. ( 10) (the dotted line in the lower panel of Fig. 2(c)) and the values determined from the simulated trajectories shows that the sphere indeed undergoes free rotational diffusion as expected.We next consider the simplest multi-sphere cluster: a dimer consisting of two identical spheres.Figure 1 shows the reference orientation for a dimer; we choose the long axis to be the u 3 axis.

Dimers
In order to verify that the dimers we consider can indeed be optically trapped, we compute the optical forces and torques on several characteristic dimers (Fig. 3).We perform all force calculations with the dimer long axis aligned with the direction of beam propagation.We observe three regimes.First, for a dimer whose constituent spheres are smaller than the wavelength of light, which we term "small," both F x and F z resemble the corresponding calculations for spheres (Fig. 3(a)).In particular, there is a linear regime near equilibrium and an equilibrium height z eq > 0. In addition, the lateral stiffness κ x is larger than the axial stiffness κ z .These results are similar to those reported by Borghese et al. for 220-nm-diameter PS spheres [6].Second, we consider an "intermediate" dimer whose spheres are comparable in size to the wavelength of light in the surrounding medium.For the dimer of 0.8-µm-diameter spheres shown in Fig. 3(b), we observe a qualitatively different F z graph with two maxima.Also, in contrast to the small dimer, the equilibrium trap stiffnesses κ x and κ z are comparable.Third, we consider a "large" dimer consisting of spheres that are larger than the wavelength, such as the 1.6-µm-diameter silica spheres in Fig. 3(c).Here, there are two heights z where F z = 0.In addition, the trap is stiffer axially than laterally: κ z > κ x .
The optical torques experienced by small, intermediate, and large dimers also differ.Since the incident beam is axisymmetric, we consider without loss of generality rotating the dimer by θ about the y axis and compute the y component of the optical torque n y .In all three cases, there is a restoring torque resulting in a stable angular equilibrium at θ = 0. Notably, for the small dimer in Fig. 3(a), the torque is approximately proportional to sin(2θ): where κ r is a rotational stiffness.However, Eq. ( 11) poorly describes the torques on intermediate and large dimers.Small, intermediate, and large dimers exhibit qualitatively different dynamical behavior in optical tweezers.Small dimers exhibit positional fluctuations similar to those for a trapped sphere and orientational fluctuations that would be expected for a particle experiencing a restoring torque.The beginning of a simulation trajectory for a small dimer is shown in Fig. 4, where we plot the CM coordinates and the laboratory-frame coordinates of all 3 particle-frame axes.Visualization 1 shows a rendering of the complete trajectory, and Fig. 5 is a still frame from the rendering.Once again, the CM position fluctuates about x = y = 0 and z = z eq .Unlike spheres, however, the restoring torque keeps the dimer axis u 3 pointing near the laboratory +z axis: u 3x and u 3y fluctuate near 0 and u 3z fluctuates near 1.The horizontal components of u 1 and u 2 range fully between -1 and +1, indicating that the dimer is free to rotate about its long axis.We analyze the behavior shown in the trajectory by computing MSDs and MSADs (Fig. 6).Because the diffusion tensor is anisotropic, we compute MSDs in the particle frame by resolving laboratory-frame displacements into components along the particle axes [31,32].We determine trap stiffnesses κ x and κ z from the linear behavior of F x and F z near equilibrium shown in Fig. 3(a).We also extract translational drag coefficients γ = k B T/D t, and γ ⊥ = k B T/D t,⊥ from the diffusion tensor.We can therefore compare the predictions of Eq. ( 9) to the MSDs calculated from the trajectories.
We also calculate the MSAD for the long axis u 3 since the orientation of this axis can be observed in experiments using techniques such as holographic microscopy [33].At short lag times τ, the optical torques do not affect the dimer: the MSAD is diffusive and is well-described by Eq. ( 10) with D r = D r,⊥ (Fig. 6(b)).At longer times, however, the effects of the optical tweezers become significant.Even approximating the torque to be proportional to sin(2θ) [Eq.( 11)], we cannot analytically solve the Einstein-Smoluchowski equation in order to determine the u 3 MSAD for arbitrary lag times.But assuming that the particle is in thermal equilibrium, we can compute the limiting value of the MSAD as τ → ∞. (See Supplemental Document for details.)The limiting value is given by lim τ→∞ where β = 1/(k B T) and F(x) is Dawson's integral.Since we can determine κ r by fitting the torque calculations shown in Fig. 3(a), we can plot the predictions of Eq. ( 12) in Fig. 6(b).
The results in Fig. 6 all differ from the analytical predictions by several percent.The discrepancies suggest that the translational and orientational degrees of freedom cannot be considered independently of each other.In particular, the optical forces on the dimer depend not only on position but also on orientation.Similarly, the optical torques depend not only on orientation but also on position.Moreover, the optical torque is only approximately described by Eq. (11).The trajectory of the intermediate dimer in Fig. 7 shows a marked difference from that of the small dimer: the dimer rotates deterministically about its long axis.This can be clearly seen in the behavior of u 1 and u 2 and in Visualization 2. The rotation is polarization-dependent: the dimer does not rotate when the incident beam is linearly polarized.Moreover, the direction of rotation reverses when the beam is right-circularly polarized (Visualization 3).We attribute this behavior to photokinetic spin-curl effects [14,34].While their analysis is only strictly valid in the Rayleigh limit, Ruffner and Grier predicted theoretically and demonstrated experimentally that the rotation frequency for particles experiencing spin-curl torques is proportional to S 3 /S 0 , where (S 0 , S 1 , S 2 , S 3 ) is the Stokes vector describing the incident polarization [14].We thus perform additional simulations with elliptically-polarized beams and determine the periods of rotation, averaged over 5 trajectories, by computing the u 1 MSAD and locating the first minimum (Fig. 8(a)).The rotation frequencies Ω we observe are indeed proportional to S 3 /S 0 (Fig. 8(b)).Here, a positive Ω corresponds to a counterclockwise rotation and a negative Ω corresponds to a counterclockwise rotation.
For the large dimer, we observe a similar deterministic rotation about the dimer axis (Fig. 9 and Visualization 4).However, there is an additional effect: the particle undergoes a wobble that is faster than the axial rotation.The effects of the wobble can be seen in the z component of u 1 and u 2 , in the lateral components of u 3 , and in the lateral position of the center of mass.Both the axial rotation and the wobble reverse direction when the beam is right-circularly polarized (Visualization 5).While experimentally detecting the axial rotation described in Fig. 8 for an optically homogeneous medium dimer would be challenging, it would be straightforward to  detect the wobble we predict here.We plan to explore these effects both computationally and experimentally in the future.Finally, our Brownian dynamics simulations reveal the existence of multiple trapping equilibria for a highly asymmetric sphere cluster.The 7-sphere cluster shown in Fig. 10 has no axes of rotational symmetry or planes of mirror symmetry, but it is the smallest rigid sphere packing that exhibits chirality [35].We simulate clusters of 0.8-µm-diameter silica spheres held in a beam that is linearly polarized in the x direction.As in the previous simulations, the cluster begins with its center of mass at the origin and in its reference orientation.
Rather than fluctuating or exhibiting deterministic motion about a single equilibrium, the simulated trajectory fluctuates about two different orientations (Fig. 11 and Visualization 6).This is most clearly visible in the behavior of the x and y components of the cluster-frame axes, although the two orientations also have different values of the CM position.
In order to determine whether the observed plateaus in the trajectory (e.g., between t = 14 and 19 s) in fact correspond to equilibria, we estimate candidate equilibrium positions and orientations from the plateaus and perform additional athermal simulations starting from the candidate equilibria.To estimate the equilibrium position, we average each CM coordinate in a plateau.To estimate the equilibrium orientation, we use the quaternion averaging algorithm of Markley and co-workers [36].We then perform additional simulations with the cluster starting from the candidate equilibria but without thermal fluctuations, and we confirm that the resulting trajectories approach asymptotic values.The refined equilibria obtained from the limiting values are shown in the horizontal lines of Fig. 11 and in Fig. 12(a)-(b).The two equilibria here correspond to a rotation of approximately 180 • about the z axis.
Remarkably, in a different simulation with the same initial conditions, we observe a different set of candidate equilibria (Fig. 13 and Visualization 7).We again perform additional simulations without thermal fluctuations and show the refined equilibria in Figs. 13 and 12(c)-(d).Once again, the two equilibria in this trajectory correspond to a 180 • rotation about the z axis.However, the equilibria of Fig. 12(c)-(d) involve a different sphere lying at the beam focus than in Fig. 12(a)-(b).
These results demonstrate the importance of considering thermal fluctuations in simulations of particles in optical tweezers.Here, thermal fluctuations allow the cluster to explore its trapping landscape and reveal multiple equilibria that we did not expect to find a priori.There may be additional trapping equilibria besides the ones we have reported.We intend to explore the trapping of this cluster in more detail as well as the transitions between equilibria.

Conclusions
Combining T-matrix-based calculations of optical interactions with a model of anisotropic Brownian motion has resulted in dynamical simulations that realistically capture the behavior of complex, wavelength-sized sphere clusters in optical tweezers.We have observed rich photokinetic effects even for the simplest case, clusters of two identical isotropic spheres, when the individual spheres are comparable to or larger than the wavelength.Our simulations have also demonstrated that multiple trapping equilibria exist for a wavelength-sized 7-sphere cluster with no symmetry.Finding this cluster's trapping equilibria solely from calculations of optical interactions at fixed positions and orientations would have been computationally challenging.Moreover, our incorporation of thermal fluctuations allowed us to more easily explore the trapping landscape.The fully general, on-demand computations of optical interactions needed for these simulations were only possible because of the speed of the T-matrix-based calculations.
Our work suggests two intriguing directions for future inquiry.First, systematically investigating effects such as the photokinetic axial rotation and wobble we have observed for intermediate and large dimers would be worthwhile.Our simulations allow continuous variation of parameters such as the particle size and refractive index and could usefully complement experimental investigations in which the particle properties cannot be as easily tuned.Moreover, since we can in principle determine the electromagnetic fields everywhere from the VSWF expansion of the incident beam and the T-matrix, more detailed investigation into the physical mechanism of these effects might be possible.This may be worthwhile since theoretical analyses based on the Rayleigh approximation cannot be strictly valid for particles of these sizes [34].
Second, it is also possible to extend this work to other systems.Using ott, we could investigate particle dynamics in beams that are not Gaussian, including Bessel beams [37] and Laguerre-Gaussian beams [38,39], both of which can carry orbital angular momentum.In addition, the generality of the T-matrix approach and the modularity of our code could allow us to consider other particles.ott natively supports T-matrix computations for small spheroids, and other packages reliably calculate the T-matrix for spheroids [40] or other axisymmetric particles [41].The discrete dipole approximation, which is implemented in ott, could also enable T-matrices to be computed for arbitrarily-shaped particles [42,43].The broad generality as well as detail of the Brownian dynamics simulations we have introduced here may make them useful both for guiding or interpreting experiments on optically-trapped particles as well as for gaining further insight into the underlying physics of optical trapping and manipulation.

Fig. 2 .
Fig. 2. (a) Optical forces on a 1-µm-diameter PS sphere in water in a 5 mW LCP beam as a function of particle position.F z is calculated as a function of z at x = y = 0. F x is calculated as a function of x at y = 0 and z = z eq (at which F z = 0).Dashed lines indicate linear approximation to optical forces near equilibrium giving stiffnesses κ x and κ z .(b) First 2 s of a 30-s-long simulated trajectory.(c) Laboratory-frame MSD and u 1 MSAD averaged from 5 independent trajectories with identical initial conditions.Shaded regions indicate standard deviations.Dotted lines are the predictions of Eqs.(9) and(10).

Fig. 3 .
Fig. 3. Optical forces and torques on small, intermediate, and large dimers in a 5 mW LCP beam.All force calculations are shown with the dimer long axis aligned with the z axis.F x and n y are calculated at z = z eq .In (c), the equilibrium position with z eq > 0 is chosen.Inset in (a) shows definition of θ.Dotted lines in force graphs indicate linear approximations near equilibrium.The dotted line in the torque graph for (a) shows a small-angle fit to Eq. (11).

Fig. 4 .
Fig. 4. First 2 s of the CM trajectory and rotation matrix elements for a small dimer composed of 0.4-µm-diameter silica spheres in a 5 mW LCP trap.The complete 30-s-long trajectory is shown in Visualization 1.

Fig. 5 .
Fig. 5. Rendering of an optically-trapped dimer composed of 0.4-µm-diameter silica spheres from Visualization 1, which shows the trajectory of Fig. 4. The stripes are drawn to enable visualizing rotations about all axes.Left: front view.The laser propagates upwards and is focused at the center of the box, whose sides are 6 µm long.Right: perspective view.The axes are 3.5 µm long.All subsequent visualizations have the same scale, and all visualizations are slowed 5× from real time.

Fig. 7 .
Fig. 7. First 2 s of the CM trajectory and rotation matrix elements for an intermediate dimer composed of 0.8-µm-diameter PS spheres in a 5 mW LCP trap.The complete 30-s-long trajectory is shown in Visualization 2. The dimer exhibits deterministic rotation about its long axis.

Fig. 8 .
Fig. 8. (a) u 1 MSADs for dimer of 0.8-µm-diameter PS spheres in 5 mW beams with left circular polarization, elliptical polarization with Jones vector (1/ √ 2, exp(iπ/6)/ √ 2), and linear polarization with Jones vector (1/ √ 2, 1/ √ 2).MSADs are averaged from 5 trajectories.The period of the dimer's deterministic rotation is determined from the first minimum of the MSADs (arrows).No rotation is observed for the linearly polarized beam.(b) Rotation frequency Ω as a function of the ratio of Stokes vector elements S 3 /S 0 .The blue and orange points correspond to the MSADs shown in (a).The sign of Ω corresponds to the direction of rotation.Dashed line: linear fit.

Fig. 9 .
Fig.9.First 2 s of the CM trajectory and rotation matrix elements for a large dimer composed of 1.6-µm-diameter silica spheres in a 5 mW LCP trap.The complete 30-s-long trajectory is shown in Visualization 4. The dimer exhibits both a slower deterministic rotation about its long axis (u 3 ) as well as a faster wobble of its long axis.

Fig. 10 .
Fig. 10.Reference orientation for chiral 7-sphere cluster.(a) and (b): Space-filling models viewed from two perspectives.(c): Ball-and-stick model viewed from same perspective as (b).The red sphere causes the cluster to be chiral.

Fig. S2. u 3
Fig. S2.u 3 MSADs calculated from simulations of a particle at T = 295 K with D r,⊥ = 5 s −1 and the labeled values of κ r .Each MSAD is averaged from 5 trajectories.Dashed lines indicate predictions of Eq. (S18), and the dotted line indicates the predictions of Eq. (S19).