Journal of Quantitative Spectroscopy & Radiative Transfer 195 (2017) 66–75

Contents lists available at ScienceDirect

Journal of Quantitative Spectroscopy & Radiative Transfer

journal homepage: www.elsevier.com/locate/jqsrt

Theory and practice of simulation of optical tweezers

Ann A.M. Bui, Alexander B. Stilgoe, Isaac C.D. Lenton, Lachlan J. Gibson,
Anatolii V. Kashchuk, Shu Zhang, Halina Rubinsztein-Dunlop, Timo A. Nieminen n
The University of Queensland, School of Mathematics and Physics, Brisbane, QLD 4072, Australia

art ic l e i nf o a b s t r a c t

Article history: Computational modelling has made many useful contributions to the field of optical tweezers. One as-
Received 30 September 2016 pect in which it can be applied is the simulation of the dynamics of particles in optical tweezers. This can
Received in revised form be useful for systems with many degrees of freedom, and for the simulation of experiments. While
15 December 2016
modelling of the optical force is a prerequisite for simulation of the motion of particles in optical traps,
Accepted 17 December 2016
Available online 20 December 2016
non-optical forces must also be included; the most important are usually Brownian motion and viscous
drag. We discuss some applications and examples of such simulations. We review the theory and
PACS: practical principles of simulation of optical tweezers, including the choice of method of calculation of
42.25.Fx optical force, numerical solution of the equations of motion of the particle, and finish with a discussion of
42.50.Wk
a range of open problems.
42.50.Tx
& 2016 Elsevier Ltd. All rights reserved.
87.80.Cc

Keywords:
Optical tweezers
Laser trapping
Optical force
Optical torque
Light scattering

Contents

1. Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 66
1.1. The need for simulations. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 67
1.2. Applications of simulations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 67
2. A recipe for simulation, part 1 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 68
3. Optical forces and torques. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 70
4. Viscous drag. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 72
4.1. Wall effects . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 73
5. A recipe for simulation, part 2: Brownian motion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 73
5.1. Nonspherical particles . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 73
6. Open questions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 74
Acknowledgements . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 74
References . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 74

1. Introduction computation of optical forces and torques is a light scattering

problem. While this is a challenging problem, and much work
The optical forces in optical tweezers result from the interac- remains to be done, there has been a great deal of progress, and for
tion of the trapping beam with the trapped particle. Thus, the many situations, it is straightforward to obtain the optical force
and torque. However, if we wish to simulate the behaviour of
n
Corresponding author. particles within optical tweezers, the optical force is only one of
E-mail address: [email protected] (T.A. Nieminen). the necessary ingredients. We will discuss some applications and

http://dx.doi.org/10.1016/j.jqsrt.2016.12.026
0022-4073/& 2016 Elsevier Ltd. All rights reserved.
 A.A.M. Bui et al. / Journal of Quantitative Spectroscopy & Radiative Transfer 195 (2017) 66–75 67

examples of such simulations, and review the theory and princi- freedom into feasibility. However, even with parallelisation, we
ples of simulation of optical tweezers. will still rapidly run into the limits of practicality due to the ex-
ponential growth of computational time with the number of de-
1.1. The need for simulations grees of freedom. Therefore, it can be necessary to resort to si-
mulation to obtain information we might prefer to find from a
Since it is usually straightforward to calculate the optical force complete force map. This will typically involve non-spherical
on a trapped particle, it is possible to characterise the trap by particles or multiple particles.
determining the optical force as a function of particle position (and On the other hand, even if it is feasible to calculate a complete
orientation if the particle is non-spherical). At first glance, this optical force map for the trap, we might still wish to perform si-
appears to provide complete information about the trap, and we mulations. In particular, a force map doesn't contain information
might ask what need there is to perform simulations. There are about the dynamics of a particle in the trap. While the optical force
two main answers to this question. First, it is not always feasible to —which the force map provides—is a key factor in the dynamics of
generate such a force map of the trap. Second, while a force map of the particle, the particle is also influenced by other forces: viscous
drag, thermal forces (driving Brownian motion), and possibly in-
this type does contain complete information about the trap in
teraction with other parts of the environment. If the dynamics are
some sense, it doesn't directly answer all questions we might have
of interest, we can use simulation to uncover it.
about the trap. In particular, the dynamics of a particle in the trap
To explore the dynamics of a particle in the trap, it can be
depend on its interaction with the surrounding environment as
possible, and advantageous, to use a pre-calculated force map. If it
well as the optical force. The dominant elements of that interac-
is feasible to calculate a complete force map with reasonable re-
tion are often Brownian motion and viscous drag, but other types
solution, the optical force at any position can be found by inter-
of interaction can also be important. Where the dynamics them-
polating between the points in the force map where the forces are
selves are the object of study (e.g., escape probabilities, synchro-
known. This interpolation can be performed very quickly (the
nised dynamics of trapped particles, etc.) or have a major impact
computational implementation should avoid copying the force
on the behaviour of interest (e.g., in the simulation of measure-
map to perform the interpolation). The required accuracy of the
ments to test calibration procedures), it is necessary to take these
interpolated force will determine the minimum resolution of the
non-optical forces into account.
force map. This resolution of the force map, along with the re-
The first of these cases results from situations with many degrees
quired spatial extent of the simulation, determines the number of
of freedom. To map the force as a function of position with useful
points required in the force map. If this exceed the number of time
(but not high) resolution typically requires about 30 steps along each steps required in the simulation, then direct calculation will be
degree of freedom (giving about 10 steps as forces change from zero more efficient. However, often the number of time steps will be
to a maximum value). If it takes 1 second to calculate the optical force much greater, and using a force map to find the optical force will
at a single position, this will give required computational times for be much more efficient. This will often be the case for optical traps
different degrees of freedom (DOF) of: with 2 or 3 degrees of freedom. An extreme case of this is where
the particle remains very close to its equilibrium position, and the
1 DOF Example: calculating axial and/or radial force–position trap can be represented in terms of a spring constant (which will
curves; finding equilibrium position along beam axis, generally be a diagonal tensor, with different spring constants in
and axial and radial spring constants. 30 to 60 points. different directions, or even a non-diagonal tensor).
Time: 0.5–1 minute.
2 DOF Example: mapping force for a spherical particle in a ro- 1.2. Applications of simulations
tationally symmetric trap (e.g., circularly polarised
Gaussian beam). 302 ≈ 1000 points. Time: E15 minutes. There are many possible applications of this. Most fall into
3 DOF Example: mapping force for a spherical particle in a trap three broad categories: simulations to understand experiments
lacking rotational symmetry (e.g., linearly polarised that have been performed, simulations to predict the results of
Gaussian beam). 303 ≈ 30, 000 points. Time: E8 hours. potential experiments, and simulations to explore optical traps
4 DOF Example: mapping force for a rotationally symmetric and the dynamics of trapped particles in ways that are not ac-
non-spherical particle in a rotationally symmetric trap. cessible experimentally.
304 ≈ 106 points. Time: E 10 days. The first of these, simulations of experiments that have been
5 DOF Example: mapping force for a rotationally symmetric performed, can be simply seeing if a simulated experiment mat-
non-spherical particle in a trap lacking rotational sym- ches measured results. This can be very useful if the experimental
metry; two spherical particles in a rotationally symmetry results are surprising. If agreement between simulated and mea-
trap. 305 ≈ 3 × 107 points. Time: E1 year. sured results is obtained, the physics and models used in the si-
6 DOF Example: mapping force for a non-spherical particle mulation adequately model reality. If agreement is not obtained,
lacking rotational symmetry in a trap lacking rotational then the model is either incomplete (e.g., physics not included
symmetry; two spherical particles in a trap lacking ro- significantly affect the measured results) or elements of the model
tational symmetry. 306 ≈ 109 points. Time: E30 years. are incorrect (invalid approximations, mathematical errors, in-
correct implementation in software, numerical errors).
Additional particles will add 2–3 translational degrees of free- For example, [79] observed the appearance of a third trapping
dom (depending on the symmetry of the trap) and 0, 2, or 3 ro- equilibrium position as two optical traps were moved close to-
tational degrees of freedom (depending on the symmetry of the gether. In this case, simulations were valuable for confirming that
particle). If the trapping beam varies in time, this adds another the third trap can be produced in this two-beam configuration,
degree of freedom, although if the time variation consists of even if the two trapping beams are not mutually coherent, i.e., the
switching between a small number of fixed positions, this will third trap doesn't depend on interference between the two trap-
only multiply the number of required calculations and the com- ping beams.
putational time by a small number. [33] used a combination of experiment and simulation to ex-
The above times do not take parallelisation of the calculations plore the transition from overdamped motion to underdamped
into account—this can readily bring one or two more degrees of motion as the size of trapped particles was reduced.
68 A.A.M. Bui et al. / Journal of Quantitative Spectroscopy & Radiative Transfer 195 (2017) 66–75

[83] observed the trajectories of particles in a laser speckle performed on experimental data can be performed on the simu-
pattern, and compared the observed trajectories with simulated lated data. The simulated calibration or force measurement can be
trajectories in speckle patterns with the same average intensity. In compared with the actual optical force in the simulation. Examples
this case, the simulations do not aim to exactly replicate the ex- of simulations of this type include [86,84], and [28]. Such simu-
perimental situation, but to replicate it in a statistical sense. That lations can also readily include non-spherical particles ([10]). Si-
is, the speckle fields in the simulations have the same statistical milar comparison on known quantities in the simulation and si-
properties as the experimental speckle fields. Qualitative and mulated measurement of these quantities can be carried out for
statistical agreement between observed and simulated trajectories methods to measure properties of the surrounding medium, such
demonstrates that the observed behaviour is general, and does not as its viscoelasticity ([25]).
arise due to some abnormality in the experimental case. An application where the trapping beam varies in time is the
[89] determined the non-conservative force field from the simulation of control methods, where the position or power of the
motion of a trapped particles. Their experimental results were beam can be varied to achieve a desired effect ([5,1,46]). Variation
supported by simulations. In the experiment, the force field is of the beam power over time introduces an additional degree of
inferred from the motion of the particle, while in the simulations, freedom, and movement of the beam introduces 2 to 4 additional
the force field is known. This allows validation of the procedure degrees of freedom (time and 1 to 3 spatial degrees of freedom).
used to obtain the experimental force field. Similarly, improved trapping methods can be explored ([80]).
Similarly, the knowledge of the optical forces available in si- Simulations can be aimed at a more general exploration of the
mulations was used to validate escape force calibrations on chro- behaviour of optically-trapped particles. These can be specifically in-
mosomes by [41]. Measurements of the force required to pull vestigating the dynamics of trapped particles ([4,90,19,75,12,81]).
chromosomes free from an optical trap were performed in order to Another common goal is the study of the behaviour of non-spherical
estimate the forces exerted on chromosomes by a cell during cell particles, where the additional degrees of freedom motivate the use of
division (mitosis) from the power required to halt the motion simulations ([76,77,11]). These can include optically-driven micro-
([24]). Since the exact size and refractive index of the chromo- machines. One example is the use of simulations to determine the
somes were uncertain, a further series of experiments and simu- optimum illumination to drive a corrugated rotor with maximum
lations on the escape of spheres from optical traps due an applied torque efficiency, while retaining stable three-dimensional trapping
force were performed ([9]), revealing the dependence of the es- ([47]). Another example is an optical “wing”, consisting of a semi-cy-
cape trajectory, and the escape force, on the trapping power and lindrical rod ([2,78]), Such a structure can generate lift—an optical
rate of increase of the applied force. force acting normal to the direction of illumination—in addition to the
Simulations that deliberately differ from the experiments can
expected radiation pressure force. Simulations by [3] show that com-
show the effect of the difference on observations or measure-
plex rocking motion can occur. The simulations allow the effects of
ments. For example, [16] used Allan variance to quantify noise in
time-varying illumination, producing a parametrically driven non-
optical tweezers setups. Simulations were used to obtain (simu-
linear bistable oscillator, to be explored.
lated) data free of noise and long-term drift, providing a suitable
Finally, simulations can be useful for educational purposes
baseline for comparison with experimental results.
([85,66]).
As noted above, simulations are often necessary when the
For some of these simulations, it is not necessary to accurately
system has many degrees of freedom, such as when there are non-
model the dynamics of the particle. For example, to determine the
spherical or multiple particles. [8] used simulations to explore the
equilibrium position and orientation of a non-spherical particle
shape dependence of the trapping behaviour of non-spherical gold
within the trap, it is not necessary to correctly model the viscous
nanoparticles. [7] used simulations to support experiments rota-
drag. The translational and rotational drag tensors can be ap-
tional dynamics of multiple spheroidal particles in a dual beam
proximated by Stokes drag for a sphere, even if the particle is non-
trap. Following observations of optically-driven oscillations of el-
spherical. Brownian motion can be ignored, although it (or random
lipsoidal particles ([52]), [50] performed simulations to under-
jitter providing similar random motion) can be useful for pre-
stand the physical basis of the observed behaviour.
venting the particle from getting stuck in an unstable equilibrium.
In these examples above, simulations were performed to sup-
This can happen, for example, if the particle is a flat disc, which
port experiments. The opposite of this, where experiments are
would tend to align with its symmetry axis normal to the beam
performed to support simulations, is also common. The aim of the
axis ([6]), but will be in an unstable equilibrium if the simulation is
simulations can vary greatly, from demonstration of the feasibility
of a particular experiment before performing it, to using simula- begun with the disc on the beam axis, with its symmetry axis
tion as a tool to help design the experiment, through to a broad along the beam axis.
series of exploration via simulation with experiments being per- However, for many types of simulations, where actual or pro-
formed to validate the simulations. In this last case, the experi- spective experiments are being simulated, it is often important to
mental work might consist of only a small fraction of the range accurately model the dynamics, and to include all important de-
covered by the simulations. If the simulation method and im- tails of the interaction of the particle with its environment. At
plementation is already known to be reliable from previous vali- minimum, this can be expected to include viscous drag and
dation, then the reported work might consist purely of Brownian motion.
simulations.
Some of this more simulation-focussed work is similar to the
experiment-focussed work described above. For example, similar 2. A recipe for simulation, part 1
to the work of [89,67] also explored the non-conservative forces in
optical traps. The motion of a particle of mass m subject to a force F (t ) can be
As noted above, simulations allow the optical forces to be calculated from
known, and are therefore valuable for testing calibration methods.
d2 r
Simulated measurements, of the type that would be measured m = F (t ),
dt 2 (1)
experimentally in order to calibrate an optical trap, or to de-
termine the optical force field from the motion of a trapped par- given initial conditions for the position r (t = 0) and velocity
ticle, can be generated, and the same analysis that would be v (t = 0). The force F (t ) is the sum of contributions from various
 A.A.M. Bui et al. / Journal of Quantitative Spectroscopy & Radiative Transfer 195 (2017) 66–75 69

sources: moving at terminal velocity, with viscous drag balancing the sum
of the other forces. When is this condition satisfied? If a particle is
F (t ) = Foptical + Fweight + Fbuoyancy + Fdrag + FBrownian + Fother , (2)
initially at rest in a fluid, and a force F is suddenly turned on, the
where we have explicitly listed the optical force, weight, buoyancy, approach to terminal velocity v0 = − D−1F will be characterised by
viscous drag, and thermal forces driving Brownian motion. We a time constant τ such that
have also included Fother to represent any other forces present. The v (t ) = v0 (1 − exp ( − t /τ )). (8)
weight and buoyancy are straightforward, with
Fweight = mg = ρparticle V g
The time constant τ is
(3)
τ = m|v0| /|F| , (9)
and
which, for a spherical particle, becomes
Fbuoyancy = − ρmedium V g, (4)
τ = 2ρparticle a2/( 9η), (10)
where V is the volume of the particle, ρparticle and ρmedium are the
densities of the particle and the surrounding medium, and g is the where a is the radius of the particle, and η is the (dynamic) visc-
local gravitational acceleration. These can almost always be trea- osity of the surrounding medium. Notably, this is independent of
ted as constant, and present no difficulty for numerical solution of the force. For a particle of radius a = 1 μm in water, this gives a
the differential Eq. (1). time constant of τ ≈ 2.4 × 10−7s. If we are calculating the motion of
The optical force Foptical is the force that typically has the most the particle using a time step Δt large compared to this (e.g.,
attention paid to it in discussion of optical tweezers and optical Δt > 10τ ), we can safely use Eq. (7). If we are using a time step
trapping. Depending on the particle in question (and the optical similar to or smaller than τ, we should, strictly speaking, use Eq.
trap), calculation of the optical force can vary from a formidable (1). In practice, if the force only changes by a small amount over
computational challenge to an already-solved problem with free- the time step Δt , the particle will already be moving at close to
ly-available implementations. We will discuss the calculation of terminal velocity, and Eq. (7) will yield acceptable results. [33]
optical forces in the following section. For the moment, we will consider a case where the difference between Eqs. (1) and (7)
consider the spatial scale over which the optical force varies, since matters over short times.
this directly affects the solutions of differential Eq. (1). In a typical Motion at the very low Reynolds numbers typical in optical
optical trap, the energy density varies from small to large values traps is outside our everyday experience. [71] gives an excellent
over a distance of half a wavelength or more. As a result, the op- and accessible overview.
tical forces will vary from small to large over a similar length scale, Brownian motion presents a serious difficulty: numerical so-
or over a distance comparable to the particle radius. For example, lution of differential equations such as (1) and (7) typically depend
when a large spherical particle is centred on the beam axis in on using a time step sufficiently short so that the time-varying
typical Gaussian beam optical tweezers, the radial force is zero, quantities in the equations (the forces, the position, and the ve-
and when the edge of the particle is on the beam axis, the radial locity) only vary by a small amount over the time step. However,
force is approximately maximum. For such a trap, we can assume no matter how short a time step is chosen, classical Brownian
that the length scale over which the optical force varies is the motion (i.e., Brownian motion in a continuous fluid [21]) always
larger of the two (i.e., the maximum of the half-wavelength and varies by a large amount. Therefore, it is simplest to remove
the particle radius). However, even for large particles, the force can Brownian motion from our differential equations, and treat it se-
vary over the half-wavelength scale, if interference effects are parately. We will return to this point after discussion of optical
important, such as when trapping in interference fringes produced forces, viscous forces, and Brownian motion.
by mutually-coherent counter-propagating beams. Knowledge of The final force in our differential equations, Fother , can represent
the length scale of variation in the optical force allows us to esti- many possible forces. For example, adhesion forces between par-
mate a suitable maximum distance to allow the particle to move ticles or particles and the microscope slide, electrostatic forces,
over a time step when numerically solving Eq. (1). magnetic forces, and more. Needless to say, some of these forces
The viscous drag force has major effects on the numerical so- can be challenging to model accurately. In some cases, a similar
lution of the differential equation describing the motion. We will approach to that suggested above for Brownian motion will be
discuss details of the calculation of the viscous drag later, and useful: remove the force from the differential equation, and treat it
restrict the current discussion to these effects on the solution. separately.
Since, typically, trapped particles are microscopic and are trapped Finally, for non-spherical particles and some spherical particles,
in a viscous environment, the interaction between the particle and it is necessary to consider rotational motion as well as transla-
the fluid is characterised by very low Reynolds numbers. In this tional motion. In this case, our differential equation will include
case, the viscous drag is linearly related to the velocity: optical torques, viscous drag torque (which will be linearly related
Fdrag = Dv, (5) to the angular velocity by a rotational drag tensor), and other
torques. These can be treated in a similar manner to their trans-
where D is a 3 × 3 drag tensor. Commonly, it is simply stated that lational counterparts. Rotational Brownian motion can be treated
since the Reynolds number is very low, we can neglect the inertial in a similar manner to translational Brownian motion. Weight and
term in Eq. (1) (recalling that the Reynolds number is the ratio of buoyancy can often be ignored, since for many particles, they
inertial effects to viscous effects in fluid motion), reducing Eq. (1) produce no torque about the centre of the particle. One compli-
to cation is that it is often desirable to perform the calculations of
0 = F (t ). (6) optical force and torque in a coordinate system fixed to the par-
ticle, necessitating transformations between the particle frame
Substituting (5), we obtain the first-order differential equation and the stationary frame. Due to the analogous nature of rotational
motion compared with translational motion, the rotational equa-
dr
−D = Foptical + Fweight + Fbuoyancy + FBrownian + Fother . tions of motion can be readily written following the translational
dt (7)
equations of motion, replacing masses with moments of inertia,
Note that differential Eq. (7) assumes that the particle is always forces with torques, and translational drag tensors with rotational
70 A.A.M. Bui et al. / Journal of Quantitative Spectroscopy & Radiative Transfer 195 (2017) 66–75

drag tensors. It should be noted that for particles with a chiral p = Ta, (13)
shape, rotational and translational motion can be coupled through
viscous drag [56]; in this case, an addition coupling tensor will be where T is the T-matrix, or transition matrix, or system transfer
included in both the translational and rotational equations of matrix. This assumes that the electromagnetic properties of the
motion. For descriptions of simulations involving rotational mo- particle are linear and constant (i.e., the particle doesn't change
tion, including equations of motion, see [11,10]. over time). With this description of scattering, the scattering
properties of the particle and the details of the incident field are
separated, leading to high efficiency for repeated calculation. Once
the T-matrix T has been calculated, it can be used repeatedly for
3. Optical forces and torques
calculation under different illumination conditions (as long as the
wavelength remains the same). Thus, as the particle moves within
The computational modeling and simulation of optical twee-
the optical trap, the T-matrix T remains constant, and the incident
zers is essentially a light scattering problem ([61,59]). The trapping
field, described by a changes. The changes in a can be found by
beam interacts with the trapped particle—this is the scattering
using the transformation properties of the basis functions ψn(inc)
aspect of the problem—and a force results. As there are a large
under translation and rotation. It is important to note that the
number of computational approaches to light scattering ([39]),
requirement that the particle not change over time means that it is
there are a large number of computational approaches to calcu-
necessary to use a coordinate system in which the particle is fixed.
lating optical forces ([59,38]). A complete review of all of the
Thus, calculations of the optical force and the torque are per-
methods would be a monumental (and book-length, if not multi- formed in the particle rest frame.
volume) task, and we will not attempt it here. Instead, we will Further efficiency results from analytical integration of the
discuss the elements of calculation of optical forces that are most momentum flux, using known results for products of integrals of
important for deciding which method will be used for such cal- the basis functions over a sphere. This reduces the formulae for
culations, and refer readers to appropriate technical literature for optical forces and torques from integrals to sums of products of
particular methods. the mode amplitudes ([59]). This avoids the need to calculate the
We will begin with an overview of the T-matrix method and fields over a grid of points in order to perform such integral nu-
why it is often the method of first choice for simulations. Note that merically, and also avoiding numerical error due to the resolution
for a spherical particle, the T-matrix method is essentially of the computational grid.
equivalent to generalised Lorenz–Mie theory (GLMT) ([29]). Then, However, it is not always feasible to calculate the T-matrix of
after noting cases where the T-matrix method might not be the the particle. The most common methods for calculating the
best choice, or even feasible, we review some basic principles of T-matrix of a particle are generalised Lorenz–Mie theory (GLMT),
calculating optical forces that can affect the choice of alternative when the particle is a uniform isotropic sphere (but note that
methods. GLMTs exist for non-spherical particles as well ([31])) and the
In general, it is safe to conclude that where calculation of the extended boundary condition method (EBCM), also known as the
T-matrix is feasible, the T-matrix method appears to be the ideal null-field method, developed by [87,88]. However, other methods
method. The T-matrix method is not a method of calculating light are possible ([51,40,62,29,55,49,48]). A comparison of EBCM,
scattering by a particle, but a formalism in which the already cal- point-matching ([62]), and the discrete dipole approximation
culated scattering properties of the particle can be expressed in the (DDA) ([48]) is given by [72].
form of the T-matrix. The extended boundary condition method Where it is impractical or impossible to calculate the T-matrix,
(EBCM) is widely used to calculate the T-matrix, being the original other methods can be sought. In many ways, the calculation of
method used by [87,88]. Thus, “T-matrix method” is often used optical forces and torques is a simple scattering problem. Often,
synonymously with EBCM, but the distinction between them there is a single particle, comparable in size to the wavelength, and
should be recognised ([59,29,31]). sufficiently far from other particles and surfaces so that multiple
In the T-matrix method, we represent the incident and scat- scattering can be ignored (it should be noted that “multiple scat-
tered fields in terms of discrete sets of vector-valued basis func- tering” is a rather artificial concept ([54]), and it is possible to treat
tions ψn(inc) and ψn(scat) , where n is a mode index labelling the scattering by a single object in a multiple scattering formalism
functions, each ψn being a solution of the vector Helmholtz (e.g., using DDA) or scattering by a group of objects in a single
equation. Using these bases, we can write the incident field am- scattering formalism ([30])). The incident field is monochromatic
plitude as and coherent. The complication is that the incident field is not a
∞ plane wave, but a focussed beam. This, and the desired outputs
E (0inc) = ∑ an ψn(inc), being the force and torque rather than the fields, or scattering
n (11) cross-sections, or scattering patterns, means that existing com-
putational implementations of particular methods might be
where an are the mode amplitudes (or expansion coefficients) of the
unsuitable.
incident wave, and the scattered wave amplitude as
There are some general theoretical points that merit discussion,
∞
since they can affect the choice of computational method or details
E 0(scat) = ∑ pk ψk(scat),
(12) of how a method is implemented. In order to calculate the forces,
k
there are two different approaches that we can take. First, we can
where pk are the mode amplitudes of the scattered wave. For use conservation of momentum, and find the difference between
computational practicality, these sums must be truncated at some the incoming momentum flux and the outgoing momentum flux
finite nmax . For a basis set of vector spherical wavefunctions, as of the light. This difference is the rate at which momentum is
usually used in the T-matrix method, the truncation criterion gi- transferred to the particle—that is, the force exerted on the par-
ven by Brock is suitable, giving a relative error due to truncation of ticle. Second, we can directly calculate the force using the Lorentz
about 10 6 ([60]). force law (or the Helmholtz force law, or other force law). These
With truncation, the mode amplitudes of the incident and two approaches are summarized in Fig. 1.
scattered waves can be written as column vectors a and p , and At this point, one might be surprised to discover that there are
their relationship can be expressed in matrix form as multiple expressions given in the literature for the momentum
 A.A.M. Bui et al. / Journal of Quantitative Spectroscopy & Radiative Transfer 195 (2017) 66–75 71

Fig. 1. The optical force can be calculated from the momentum flux or by an
electromagnetic force law.

flux of light, and also multiple electromagnetic force laws. This is

the Abraham-Minkowski controversy, where we encounter com-
peting expressions for the momentum density of an electro-
magnetic field ([68]). With more than one possible expression for
the momentum density, how can we choose the correct one, or at
least the best one to use?
The key is to note that while there are different expressions for
the electromagnetic momentum density in material media, the Fig. 3. Choices of surface over which to integrate momentum flux. We can choose a
different versions all correspond to identical expressions for the surface conforming to the surface of the particle, or a surface of simple geometry
total momentum density in the material medium associated with enclosing the particle. Alternatively, we can choose a spherical surface in the far
field, which can simplify the calculation by allowing us to make far-field
the electromagnetic wave. Where the electromagnetic momenta
approximations.
differ, the difference is matched by opposing differences in what is
labelled material momentum or interaction momentum. Since we
case, we obtain a volume charge density from the variation of the
must calculate the total force, whether or not it is described as
dielectric polarisation in the particle, as shown in Fig. 2. The
purely electromagnetic or the sum of an electromagnetic and a
Helmholtz force law gives the force acting on the dipole moment
material force, the difference in how the total momentum is di-
per unit volume, and the Lorentz force law the force acting on the
vided into electromagnetic and material (and possibly other)
equivalent charges. In both cases, the total force is identical.
components is not fundamental. However, it is convenient to be
If we choose to find the force from conservation of momentum,
able to calculate a single quantity rather than multiple quantities we need to choose a closed surface over which to integrate the
that, when added, equal that single quantity. Noting that for cases momentum flux. There are three principle choices: a surface
where the electromagnetic properties of the medium can be de- conforming to the surface of the trapped particle, a surface of
scribed completely with a constant permittivity and permeability, simple geometry close to the particle enclosing it, and a spherical
the Minkowski momentum is the total momentum ([68]), this is surface in the far field, as shown in Fig. 3. The latter two of these
the simplest choice of momentum density. are often the best choices. If a surface in the far field is chosen, it
With each possible choice of momentum density, there is an can be possible to make far-field approximations to simplify cal-
associated electromagnetic force law. If we begin by choosing a culation of the momentum flux. If a nearby surface of simple
force law, we can derive an expression for the momentum density geometry is chosen, the same surface can be used for particles of
and momentum flux of the electromagnetic field. Doing this in different shapes, simplifying implementation.
reverse, we can begin with an expression for the momentum, and If there are multiple particles within the trap, we will usually
obtain an electromagnetic force law. The most commonly en- need to calculate the optical force acting on each particle. This is
countered force laws are the Lorentz force law, giving the force important, for example, when considering optical binding
acting on charges and currents, and the Helmholtz force law which ([14,15]). In this case, we cannot find the individual forces by in-
includes forces acting on induced dipole moments. These con- tegration of the total field in the far field. We can instead use
nection between these force laws can be seen if we consider the surface surrounding each individual particle, integration over each
induced dielectric polarisation in a particle. We can represent this of which will yield the force acting on the enclosed particle. In a
either by the dipole moment per unit volume, or by equivalent multiple-particle T-matrix formulation, it is still possible to use the
charges. For a uniformly polarized sphere, the dipole moment per usual single-particle summation formulae ([59]) to find the force,
unit volume is uniform throughout the sphere, but we can replace if we use incident and scattered field mode amplitudes for each
this by an equivalent surface charge density. In the more general particle individually.
While there are many possible methods, they largely fall into
three groups: finite element methods (FEM), finite difference
methods, of which the most notable variants are the finite-dif-
ference time-domain method (FDTD) and the finite-difference
frequency-domain method (FDFD), and approximate methods,
notably Rayleigh scattering and geometric optics or ray optics.
In the finite element method (FEM), the computational space is
subdivided into finite elements (Volakis et al., 1998). The values
relevant to the PDE inside or at the surface of each element are
approximated by some known function, perhaps linear or a higher
order function ([82]). The interaction between each element or an
element and its surrounding elements is described by a matrix
Fig. 2. The optical force can be calculated from the momentum flux or by an that depends on the particular definition and the chosen division
electromagnetic force law. of the computational space. For scattering problems, FEM might be
72 A.A.M. Bui et al. / Journal of Quantitative Spectroscopy & Radiative Transfer 195 (2017) 66–75

used to refer to a number of methods that involve solving either

sparse or dense systems of equations. The most important FEM as
far as optical tweezers is concerned is the the discrete dipole ap-
proximation (DDA). The physical interpretation of DDA involves
representing a large scattering particle by multiple smaller inter-
acting dipoles whose polarisability is known ([54]). The interaction
between each dipole is described by a dense matrix; the resulting
linear system approximates the scattering by the combined object.
[92] provide a good overview of DDA including recent develop-
ments and comparisons to other methods. Unlike other FEM, DDA
doesn't require the space surrounding the scatter to be discretised,
unless the surrounding space is inhomogeneous or contains other
objects. DDA is more suitable for smaller isolated particles and
particles with smaller (relative) refractive indexes due to the re-
quirement to solve a dense linear system.
The finite-difference time-domain (FDTD) method refers to a
method described by [91] for solving systems of coupled partial
differential equations. Although it was originally formulated for
solving the Maxwell equations, FDTD can also be applied to other
systems of differential equations. The original formulation of FDTD
for the Maxwell equations involved calculating the electric and
magnetic fields at locations on a structured grid spanning the
computational space. Spatial derivatives in the Maxwell equations
are calculated using second order finite difference approximations
involving the adjacent locations on the structured grid. The fields
are advanced through time using a leapfrog scheme with second
order accuracy, where the electric and magnetic fields are updated
at alternate half integer time steps. Since Yee's original method,
there have been numerous improvements and specialisations such
as unconditionally stable methods or single step methods
([37,73]). One disadvantage of the FDTD method is difficulty in
representing objects with smooth surfaces using a structured
Cartesian grid; when the surface doesn't conform to the Cartesian
grid, this introduces staircasing error ([37]). The simplest approach
is to increase the grid resolution. However, this results in a major
increase in the computational requirements. Other alternatives
include using non-Cartesian grids such as spherical or circular
grids, sub-gridding certain regions or incorporating a local dis-
tortion near curved object boundaries ([35]). Use of non-Cartesian
Fig. 4. Comparison of different computational methods for simulating optical
grids requires calculation of the Jacobian and correcting the FDTD tweezers. Further from the centre is better; e.g., the Rayleigh approximation, GLMT
update equations appropriately, special attention should also be and EBCM are the best of these methods for the calculation of forces on very small
given to discontinuities in the mesh. particles. The particle complexity includes both the particle geometry and com-
The finite-difference frequency-domain (FDFD) method is very position (inhomogeneity, anisotropy, non-linearity). The accuracy and computa-
tional requirements depend on the type of problem being solved, so this compar-
similar to FDTD except it assumes time harmonic solutions for the
ison should be treated as a qualitative guide, rather than an exact quantitative
incident and scattered fields ([49]). FDFD is very similar to FEM comparison.
where only interactions between adjacent elements are con-
sidered. This results in a linear system describing the scatter. Un- methods depend on the type of problem being solved, so this
like FDTD, FDFD performs the scattering calculation for only a comparison should be treated as a qualitative guide, rather than an
particular frequency, while FDTD is able to calculate the scattering
exact quantitative comparison. This comparison is based on our
of multiple frequencies simultaneously.
own experience with these methods.
The choice to use a particular computational method to model
optical tweezers greatly depends on the regime the problem falls
into. For very large and very small particles the ray optics and
Rayleigh approximations are able to model particles with rea- 4. Viscous drag
sonable accuracy. DDA approaches the Rayleigh limit for small
particles but is able to simulate larger particles with fairly high Viscous drag is a key factor in the dynamics of a particle in an
accuracy but scales relatively poorly with memory and time. The optical trap. It determines the speed at which the optical forces
FDTD and FDFD methods are able to simulate an extended range of and torques will move or rotate the particle, and it also affects
particles but rely on being able to discretise the computational Brownian motion. For simulation of optical tweezers, we wish to
space to describe fine details of the scatterer or rapidly changing obtain the translational viscous drag tensor for the particle (and
fields. FDFD assumes a particular form for the wave solutions, but the rotational drag tensor, if we need to include rotation in the
is only able to simulate a single frequency; in comparison, FDTD simulation). For the case of a spherical particle, this is straight-
can easily deal with illumination such as short pulses, and can forward, since there is a simple analytical solution: Stokes drag on
include non-linear effects including frequency doubling, frequency a sphere. This gives
mixing, etc. A summary of these comparisons is presented gra-
phically in Fig. 4. The performance and capabilities of the different D = 6πηaI, (14)
 A.A.M. Bui et al. / Journal of Quantitative Spectroscopy & Radiative Transfer 195 (2017) 66–75 73

where I is the identity tensor. If the rotational drag tensor is re-

quired, this is equal to
Drot = 8πηa3I. (15)

For cylinders, convenient formulae are given by [18]. For more

general cases, it can be necessary to resort to solving the fluid flow
and calculating the drag tensors. At low Reynolds numbers, the
fluid flow is described by the Stokes equation:
η∇2v = ∇p, (16)

where v is the velocity field of the fluid flow, and p is the pressure
field. The two most promising approaches appear to be using the
general solution in spherical coordinates ([43,65]) and direct fi-
nite-difference solution.

4.1. Wall effects Fig. 5. Convergence of final position of particle undergoing Brownian motion in an
optical trap as a function of time step.

It should also be noted that nearby surfaces affect the viscous

water at 300 K, and the distance scales with the square root of Δt ,
drag on a particle. In general, this is a difficult problem ([34]).
we would need a time step of approximately 10 4 s if we want the
However, the simple case of a sphere near a plane wall has a known
distance moved to be less than 1% of the wavelength.
solution. The approximate solution by [23] is often used. However, it
We can investigate the effect of our choice of time step quan-
is only accurate when the particle is a large distance away from the
titatively. We can generate a series of discrete Brownian steps for a
wall, and fails when the distance between the wall and the closest
a time step Δt0 , and calculate the motion of the particle. Then, we
part of the particle is less than 1 particle radius. That is, it cannot be
can double the time step, and sum successive pairs of Brownian
used when it is most needed. The exact solution presents serious
steps, to obtain half the number of steps, each twice as long in
difficulties in calculation. Fortunately, a simple and very accurate
time. This can be repeated, allowing us to investigate the con-
approximation formula is available ([13]).
vergence of the calculation with decreasing time step. The dis-
placement over the time step due to Brownian motion can be used
to find an average velocity due to Brownian motion over the time
5. A recipe for simulation, part 2: Brownian motion
step; this can be expressed as an average force over the time step
and included in a predictor–corrector method such as Runge–
Brownian motion in a viscous fluid, in the absence of other
Kutta. An example is shown in Fig. 5 for a particle of radius 1 μm ,
forces, can be easily modelled. The probability distributions for
comparing the convergence of trajectories as the time step is re-
displacements of a spherical particle in each of the x, y, and z di-
duced. The comparison includes both Euler's method and a fixed-
rections over a time interval Δt are normal (Gaussian), with var-
step 4th-order Runga–Kutta method. The results indicate that a
iance equal to 2DΔt , where
time step of 10 4 s (for a distance less than 1% of the wavelength)
kB T gives a reasonably small error. If higher accuracy is desired, a
D=
6πηa (17) shorter time step can be chosen. Due to the square root depen-
dence of the distance, the scaling is relatively poor.
is the diffusion coefficient, and kB is Boltzmann's constant, T is the We recommend that a similar analysis of convergence be per-
absolute temperature, η is the (dynamic) viscosity, and a is the formed, especially if simulations are sufficiently lengthy so that a
radius of particle. Thus, it is straightforward to simulate Brownian just-small-enough for acceptable error time-step should be chosen.
motion using a Monte Carlo method. To calculate a displacement For very short time steps, the Stokes drag formula can be in-
over Δt , we can generate 3 normally distributed random numbers, appropriate ([27,42]), and for very short time steps, the transition
Rx, Ry, and Rz with variance equal to 1, giving us to ballistic motion becomes apparent ([36]). As long as the time
Δx = (2DΔt )1/2Rx (18) step is not too short, this will not present any difficulty. If simu-
lations of a particle trapped in gas are being performed, then the
for the displacement in the x-direction, and similar results in the y and transition between ballistic and continuous regimes is important
z directions. This can then be repeated for subsequent time steps, with ([45]) at larger time steps.
new random numbers Rx, Ry, and Rz generated for each time step.
Notably, classical Brownian motion of this type is self-similar 5.1. Nonspherical particles
across all time scales, i.e., fractal ([21,58]), and consequently, the
accuracy of a Monte Carlo simulation like this is independent of If the particle is nonspherical, the translation and rotational drag
the choice of step size Δt . Therefore, if we aim to simulate a series tensors will be different in different directions, and the variance of
of measurements of particle position, it is sufficient to calculate the Brownian motion will be different in different directions.
the particle position at only the times at which the position is It is simplest to calculate the random Brownian steps in the rest
measured, and Δt is the time interval between the measurements. frame of the particle, and then transform the motion to the sta-
There is no need to calculate the position for intermediate times. tionary frame.
However, in the presence of other forces, this changes. It be- The question of suitable time steps was left open earlier. From
comes necessary for the distance the particle moves in a single the discussion of Brownian motion above, we see that Brownian
time step to be small enough so that the other forces do not motion can be the main factor limiting our choice of time step.
change too much. Since the optical force can change greatly over However, it is wise to check that the optical force will not move
half a wavelength, the distance must be a small fraction of the the particle too far during the time step. A maximum time step
wavelength. Noting that a particle of radius 1 μm will move, on based on the optical force can be found, and another based on the
average, a distance of 1 μm in 2.1 s due to Brownian motion in Brownian motion. The smaller of the two can then be used as the
74 A.A.M. Bui et al. / Journal of Quantitative Spectroscopy & Radiative Transfer 195 (2017) 66–75

actual time step for the computation. If Brownian motion can be convective flows, and the convective flows can alter the tempera-
neglected, then it can be convenient to use an adaptive step size ture distribution, and also the position of the particle within the
solver for initial value problems. If we wish to calculate simulated trap (thus altering the absorption of light and the temperature
measured at specific times, we can choose time steps so that these distribution). This coupling makes the solution of the problem
specific times match times at which we calculate the position of difficult; an iterative method might be required. The time scales
the particle in the trap. Alternatively, it may be possible to inter- involved can be investigated—if conduction dominates energy
polate between the calculations. transport, then it may be possible to ignore the effect of convec-
tion on the temperature distribution, and the problem, while still
challenging, is greatly simplified.
6. Open questions

A number of open questions and unsolved problems in simu- Acknowledgements

lation of optical tweezers remain. We present a selection of them
here, in the spirit of presenting useful and interesting challenges This research was supported under Australian Research Council's
to those who wish to tackle them. Interesting work has been Discovery Projects funding scheme (project number DP140100753).
performed on some of these topics, giving a hint of many inter-
esting results yet to be uncovered.
References
 Optical force on complex particles. While optical forces and
torques exerted on a wide range of particles can be readily cal- [1] Aguilar-Ibañez C, Suarez-Castanon MS, Rosas-Soriano LI. A simple control
culated, those on large and complex particles remain challenging. scheme for the manipulation of a particle by means of optical tweezers. Int J
 Nonlinear particles. Particles with non-linear electromagnetic Robust Nonlinear Control 2011;21(3):328–337.
[2] Artusio-Glimpse AB, Peterson TJ, Swartzlander GA. Refractive optical wing os-
properties have the potential for many interesting behaviours in cillators with one reflective surface. Opt Lett 2013;38(6):935–937.
optical traps ([69,70,20]). [3] Artusio-Glimpse AB, Schuster DG, Gomes MW, Swartzlander GA. Rocking mo-
 Deformable particles. These present a double challenge. First, it tion of an optical wing theory. Appl Opt 2014;53(31):I1–I9.
[4] Banerjee A, Balijepalli A, Gupta SK, LeBrun TW. Generating simplified trapping
is necessary to calculate not just the optical and viscous forces probability models from simulation of optical tweezers system. J Comput Inf Sci
acting on the particle, but also the stresses and consequent Eng 2009;9(2):021003.
deformation of the particle. Second, the deformation results in [5] Banerjee AG, Chowdhury S, Losert W, Gupta SK. Real-time path planning for
coordinated transport of multiple particles using optical tweezers. IEEE Trans
change in the optical force. The deformable particles of most Autom Sci Eng 2012;9:669–678.
interest are red blood cells ([44,17,74]) and other cells ([32]), [6] Bayoudh S, Nieminen TA, Heckenberg NR, Rubinsztein-Dunlop H. Orientation of
but simpler objects such as vesicles, which are sometimes used biological cells using plane polarised Gaussian beam optical tweezers. J Mod
Opt 2003;50(10):1581–1590.
as simple analogs of cells, are also of interest ([63]).
[7] Brzobohatý O, Arzola AV, Šiler M, Chvátal L, Jákl P, Simpson S, et al. Complex
 Wall effects on viscous drag on non-spherical particles. The rotational dynamics of multiple spheroidal particles in a circularly polarized,
movement of non-spherical particles near surfaces is important dual beam trap. Opt Express 2015;23(6):7273–7287.
[8] Brzobohatý O, Šiler M, Trojek J, Chvátal L, Karásek V, Zemánek P. Non-spherical
in many biological systems. One interesting example is the
gold nanoparticles trapped in optical tweezers shape matters. Opt Express
motion of sperm, which dramatically change in their swimming 2015;23(7):8179–8189.
behaviour near surfaces. ([22,64]). Optical tweezers offers an [9] Bui AAM, Stilgoe AB, Khatibzadeh N, Nieminen TA, Berns MW, Rubinsztein-
opportunity to explore this behaviour, either by trapping sperm Dunlop H. Escape forces and trajectories in optical tweezers and their effect on
calibration. Opt Express 2015;23(19):24317–24330.
and measuring swimming forces (as done for free-swimming [10] Bui AAM, Stilgoe AB, Nieminen TA, Rubinsztein-Dunlop H. Calibration of
sperm by [57]) or by trapping and moving analogs near nonspherical particles in optical tweezers using only position measurement.
surfaces. Simulations would be very helpful for identifying Opt Lett 2013;38(8):1244–1246.
[11] Cao Y, Stilgoe AB, Chen L, Nieminen TA, Rubinsztein-Dunlop H. Equilibrium
changes in motion that result from changed behaviour of sperm orientations and positions of non-spherical particles in optical traps. Opt Ex-
near surfaces; such simulations would need to account for wall press 2012;20:12987–12996.
effects on the motion. [12] Cao Y, Zhu T, Lv H, Ding W. Spin-controlled orbital motion in tightly focused
 Interaction between trapped particles and complex biological high-order laguerre-gaussian beams. Opt Express 2016;24(4):3377–3384.
[13] Chaoui M, Feuillebois F. Creeping flow around a sphere in a shear flow close to
environments. The work on deformation of red blood cells noted a wall. Q J Mech Appl Math 2003;56(3):381–410.
above can be considered a special case of this. More generally, a [14] Chaumet PC, Nieto-Vesperinas M. Optical binding of particles with or without
the prescence of a flat dielectric surface. Phys Rev B 2001;64:035422.
trapped particle can interact with membranes, macromolecules,
[15] Chvatal L, Brzobohaty O, Zemanek P. Binding of a pair of Au nanoparticles in a
cells, complex fluids, etc. Modelling its interaction with such an wide Gaussian standing wave. Opt Rev 2015;22:157–161.
environment can be challenging. Where the behaviour of living [16] Czerwinski F, Richardson AC, Oddershede LB. Quantifying noise in optical
cells needs to be included (e.g., swimming behaviour of bacteria or tweezers by Allan variance. Opt Express 2009;17(15):13255–13269.
[17] Dao M, Lim C, Suresh S. Mechanics of the human red blood cell deformed by
sperm, ingestion of the trapped particle by a macrophage, etc.), optical tweezers. J Mech Phys Solids 2003;51:2259–2280.
realistic models of the behaviour are required. This is a very [18] de la Torre JG, Bloomfield VA. Hydrodynamic properties of complex, rigid,
complex problem that remains largely untouched. biological macromolecules: theory and applications. Q Rev Biophys 1981;14
 Heating and thermal effects, including convective flow. While
(1):81–139.
[19] Deng Y, Bechhoefer J, Forde NR. Brownian motion in a modulated optical trap.
heating is often ignored in optical trapping simulations—wave- J Opt A 2007;9(8):S256.
lengths and beam powers are often chosen to be such that ab- [20] Devi A, De AK. Theoretical investigation on nonlinear optical effects in laser
trapping of dielectric nanoparticles with ultrafast pulsed excitation. Opt Ex-
sorption and consequent heating in minimal, to avoid damage to press 2016;24(19):21485.
live biological specimens—there can be significant heating when [21] Einstein A. Investigations on the theory of the Brownian movement.New York:
absorbing particles are trapped. Heating introduces a wide range Dover; 1956.
[22] Elgeti J, Kaupp UB, Gompper G. Hydrodynamics of sperm cells near surfaces.
of effects, from changes in the viscosity of the surrounding fluid
Biophys J 2010;99(4):1018–1026.
due to increased temperature, convection currents, and effects [23] Faxén H. Der Widerstand gegen die Bewegung einer starren Kugel in einer
such as thermophoresis ([26]). If there are liquid–liquid or liquid– zähen Flüssigkeit, die zwischen zwei parallelen ebenen Wänden einges-
gas interfaces present, the dependence of surface tension on chlossen ist. Ann der Phys 1922;373(10):89–119.
[24] Ferraro-Gideon J, Sheykhani R, Zhu Q, Duquette ML, Berns MW, Forer A.
temperature can produce strong flows due to Marangoni convec- Measurements of forces produced by the mitotic spindle using optical twee-
tion ([53]). In general, the temperature distribution drives the zers. Mol Biol Cell 2014;24:1375–1386.
 A.A.M. Bui et al. / Journal of Quantitative Spectroscopy & Radiative Transfer 195 (2017) 66–75 75

[25] Fischer M, Berg-Sørensen K. Calibration of trapping force and response func- [58] Nelson E. Dynamical theories of Brownian motion. Princeton University Press,
tion of optical tweezers in viscoelastic media. J Opt A: Pure Appl Opt 2007;9: Princeton, NJ; 1967.
S239–S250. [59] Nieminen TA, du Preez-Wilkinson N, Stilgoe AB, Loke VLY, Bui AAM, Ru-
[26] Flores-Flores E, Torres-Hurtado SA, Páez R, Ruiz U, Beltrán-Pérez G, Neale SL, et al. binsztein-Dunlop H. Optical tweezers theory and modelling. J Quant Spectrosc
Trapping and manipulation of microparticles using laser-induced convection Radiat Transf 2014;146:59–80.
currents and photophoresis. Biomed Opt Express 2015;6(10):4079–4087. [60] Nieminen TA, Loke VLY, Stilgoe AB, Heckenberg NR, Rubinsztein-Dunlop H. T-
[27] Franosch T, Grimm M, Belushkin M, Mor FM, Foffi G, Forró L, et al. Resonances matrix method for modelling optical tweezers. J Mod Opt 2011:528–544.
arising from hydrodynamic memory in Brownian motion. Nature 2011;478:85–88. [61] Nieminen TA, Rubinsztein-Dunlop H, Heckenberg NR. Calculation and optical
[28] Gong Z, Chen H, Xu S, Li Y, Lou L. Monte-Carlo simulation of optical trap measurement of laser trapping forces on non-spherical particles. J Quant
stiffness measurement. Opt Commun 2006;263(2):229–234. Spectrosc Radiat Transf 2001;70:627–637.
[29] Gouesbet G. T-matrix formulation and generalized Lorenz-Mie theories in [62] Nieminen TA, Rubinsztein-Dunlop H, Heckenberg NR. Calculation of the
spherical coordinates. Opt Commun 2010;283:517–521. T-matrix general considerations and application of the point-matching
[30] Gouesbet G, Grehan G. Generalized Lorenz-Mie theory for assemblies of method. J Quant Spectrosc Radiat Transf 2003;79–80:1019–1029.
spheres and aggregates. J Opt A: Pure Appl Opt 1999;1:706–712. [63] Noguchi H, Takasu M. Structural changes of pulled vesicles a brownian dy-
[31] Gouesbet G, Lock JA. On the electromagnetic scattering of arbitrary shaped namics simulation. Phys Rev E 2002;65:051907.
beams by arbitrary shaped particles: a review. J Quant Spectrosc Radiat Transf [64] Nosrati R, Driouchi A, Yip CM, Sinton D. Two-dimensional slither swimming of
2015;162:31–49. sperm within a micrometre of a surface. Nat Commun 2015;6:8703.
[32] Guck J, Schinkinger S, Lincoln B, Wottawah F, Ebert S, Romeyke M, et al. [65] Pak OS, Lauga E. Generalized squirming motion of a sphere. J Eng Math
Optical deformability as an inherent cell marker for testing malignant trans- 2014;88(1):1.
formation and metastatic competence. Biophys J 2005;88(5):3689–3698. [66] Perkins TT, Malley CV, Dubson M, Perkins KK. 2010. An interactive optical
[33] Haghshenas-Jaryani M, Black B, Ghaffari S, Drake J, Bowling A, Mohanty S. tweezers simulation for science education. Proceedings SPIE 7762, pp. 776215.
Dynamics of microscopic objects in optical tweezers experimental determi- [67] Pesce G, Volpe G, Chiara De Luca A, Rusciano G, Volpe G. Quantitative assessment
nation of underdamped regime and numerical simulation using multiscale of non-conservative radiation forces in an optical trap. EPL 2009;86:38002.
analysis. Nonlinear Dyn 2014;76:1013–1030. [68] Pfeifer RNC, Nieminen TA, Heckenberg NR, Rubinsztein-Dunlop H. Momentum
[34] Happel J, Brenner H. Low reynolds number hydrodynamics.Dordrecht: of an electromagnetic wave in dielectric media. Rev Mod Phys 2007;79
Kluwer; 1991. (4):1197–1216.
[35] Hastings FD, Schneider JB, Broschat SL, Thorsos EI. An FDTD method for ana- [69] Pobre R, Saloma C. Single Gaussian beam interaction with a Kerr microsphere
lysis of scattering from rough fluid-fluid interfaces. IEEE J Ocean Eng 2001;26 characteristics of the radiation force. Appl Opt 1997;36:3515–3520.
(1):94–101. [70] Pobre R, Saloma C. Radiation force exerted on nanometer size non-resonant Kerr
[36] Huang R, Chavez I, Taute KM, Lukić B, Jeney S, Raizen MG, et al. Direct ob- particle by a tightly focused Gaussian beam. Opt Commun 2006;267:295–304.
servation of the full transition from ballistic to diffusive Brownian motion in a [71] Purcell EM. Life at low Reynolds number. Am J Phys 1977;45(January (1)):3–11.
liquid. Nat Phys 2011;7:576–580. [72] Qi X, Nieminen TA, Stilgoe AB, Loke VLY, Rubinsztein-Dunlop H. Comparison of
[37] Inan US, Marshall RA. Numerical Electromagnetics: the FDTD Method.Cam- T-matrix calculation methods for scattering by cylinders in optical tweezers.
bridge: Cambridge University Press; 2011. Opt Lett 2014;39:4827–4830.
[38] Jones PH, Maragò OM, Volpe G. Optical tweezers: principles and applications. [73] Raedt HD, Michielsen K, Kole JS, Figge MT. Solving the Maxwell equations by
Cambridge: Cambridge University Press; 2015. the Chebyshev method a one-step finite-difference time-domain algorithm.
[39] Kahnert FM. Numerical methods in electromagnetic scattering theory. J Quant IEEE Trans Antennas Propag 2003;51(11):3155–3160.
Spectrosc Radiat Transf 2003;79–80:775–824. [74] Rancourt-Grenier S, Wei M-T, Bai J-J, Chiou A, Bareil PP, Duval P-L, Sheng Y.
[40] Kahnert FM, Stamnes JJ, Stamnes K. Surface-integral formulation for electro- Dynamic deformation of red blood cell in dual-trap optical tweezers. Opt
magnetic scattering in spheroidal coordinates. J Quant Spectrosc Radiat Transf Express 2010;18(10):10462–10472.
2003;77:61–78. [75] Ren Y, Wu J, Zhong M, Li Y. Monte-Carlo simulation of effective stiffness of
[41] Khatibzadeh N, Stilgoe AB, Bui AAM, Rocha Y, Cruz GM, Loke V, et al. De- time-sharing optical tweezers. Chin Opt Lett 2010;8(2):170–172.
termination of motility forces on isolated chromosomes with laser tweezers. [76] Simpson SH, Hanna S. Holographic optical trapping of microrods and nano-
Sci Rep 2014;4:6866. wires. J Opt Soc Am A 2010;27(6):1255–1264.
[42] Kheifets S, Simha A, Melin K, Li T, Raizen MG. Observation of Brownian motion [77] Simpson SH, Hanna S. Orbital motion of optically trapped particles in La-
in liquids at short times instantaneous velocity and memory loss. Science guerre-Gaussian beams. J Opt Soc Am A 2010;27(9):2061–2071.
2014;343:1493–1496. [78] Simpson SH, Hanna S, Peterson TJ, Swartzlander GA. Optical lift from dielectric
[43] Lamb H. Hydrodynamics. 5th Edition. Cambridge: Cambridge University Press; semicylinders. Opt Lett 2012;37(19):4038–4040.
1924. [79] Stilgoe AB, Heckenberg NR, Nieminen TA, Rubinsztein-Dunlop H. Phase-
[44] Li J, Dao M, Lim C, Suresh S. Spectrin-level modeling of the cytoskeleton and transition-like properties of double-beam optical tweezers. Phys Rev Lett
optical tweezers stretching of the erythrocyte. Biophys J 2005;88:3707–3719. 2011;107:248101.
[45] Li T, Kheifets S, Medellin D, Raizen MG. Measurement of the instantaneous [80] Taylor MA, Waleed M, Stilgoe AB, Rubinsztein-Dunlop H, Bowen WP. Enhanced
velocity of a Brownian particle. Science 2010;328:1673–1675. optical trapping via structured scattering. Nat Photonics 2015;9:669–673.
[46] Li X, Cheah CC, Hu S, Sun D. Dynamic trapping and manipulation of biological [81] Trojek J, Chvátal L, Zemánek P. Optical alignment and confinement of an el-
cells with optical tweezers. Automatica 2013;49(6):1614–1625. lipsoidal nanorod in optical tweezers a theoretical study. J Opt Soc Am A
[47] Loke VLY, Asavei T, Stilgoe AB, Nieminen TA, Rubinsztein-Dunlop H. Driving 2012;29(7):1224–1236.
corrugated donut rotors with laguerre-gauss beams. Opt Express 2014;22 [82] Volakis JL, Chatterjee A, Kempel LC. 1998. Finite element method for electro-
(16):19692–19706. magnetics: antennas, microwave circuits, and scattering applications. IEEE
[48] Loke VLY, Nieminen TA, Heckenberg NR, Rubinsztein-Dunlop H. T-matrix Press, New York.
calculation via discrete dipole approximation, point matching and exploiting [83] Volpe G, Kurz L, Callegari A, Volpe G, Gigan S. Speckle optical tweezers mi-
symmetry. J Quant Spectrosc Radiat Transf 2009;110:1460–1471. cromanipulation with random light fields. Opt Express 2014;22(15):18159–
[49] Loke VLY, Nieminen TA, Parkin SJ, Heckenberg NR, Rubinsztein-Dunlop H. FDFD/ 18167.
T-matrix hybrid method. J Quant Spectrosc Radiat Transf 2007;106:274–284. [84] Volpe G, Petrov D. Torque detection using Brownian fluctuations. Phys Rev Lett
[50] Loudet J-C, Mihiretie BM, Pouligny B. Optically driven oscillations of ellipsoidal 2006;97:210603.
particles. Part II ray-optics calculations. Eur Phys J E 2014;37:125. [85] Volpe G, Volpe G. Simulation of a Brownian particle in an optical trap. Am J
[51] Mackowski DW. Discrete dipole moment method for calculation of the T Phys 2013;81:224–230.
matrix for nonspherical particles. J Opt Soc Am A 2002;19(5):881–893. [86] Volpe G, Volpe G, Petrov D. Brownian motion in a nonhomogeneous force field
[52] Mihiretie BM, Snabre P, Loudet JC, Pouligny B. Optically driven oscillations of and photonic force microscope. Phys Rev E 2007;76:061118.
ellipsoidal particles. Part I: experimental observations. Eur Phys J E 2014;37:124. [87] Waterman PC. Matrix formulation of electromagnetic scattering. Proc IEEE
[53] Miniewicz A, Bartkiewicz S, Orlikowska H, Dradrach K. Marangoni effect visua- 1965;53(8):805–812.
lized in two-dimensions Optical tweezers for gas bubbles. Sci Rep 2016;6:34787. [88] Waterman PC. Symmetry, unitarity, and geometry in electromagnetic scat-
[54] Mishchenko MI. Electromagnetic scattering by particles and particle groups. tering. Phys Rev D 1971;3:825–839.
Cambridge: Cambridge University Press; 2014. [89] Wu P, Huang R, Tischer C, Jonas A, Florin E-L. Direct measurement of the
[55] Mishchenko MI, Travis LD, Mackowski DW. T-matrix method and its appli- nonconservative force field generated by optical tweezers. Phys Rev Lett
cations to electromagnetic scattering by particles a current perspective. J 2009;103:108101.
Quant Spectrosc Radiat Transf 2010;111:1700–1703. [90] Xu S-H, Li Y-M, Lou L-R, Sun Z-W. Computer simulation of the collision fre-
[56] Moffatt HK. 1977. Six lectures on general dluid dynamics and two on hydro- quency of two particles in optical tweezers. Chin Phys 2005;14(2):382–385.
magnetic dynamo theory. in: Balian R, Peube, J-L. (Eds.), Fluid dynamics. [91] Yee KS. Numerical solution of initial boundary value problems involving
Gordon and Breach, London, pp. 149–233. Maxwell's equations in isotropic media. IEEE Trans Antennas Propag 1966;14
[57] Nascimento JM, Shi LZ, Chandsawangbhuwana C, Tam J, Durrant B, Botvinick (3):302–307.
EL, et al. Use of laser tweezers to analyze sperm motility and mitochondrial [92] Yurkin MA, Hoekstra AG. The discrete dipole approximation an overview and
membrane potential. J Biomed Opt 2008;13:014002. recent developments. J Quant Spectrosc Radiat Transf 2007;106:558–589.