Theory and simulation of the bistable behaviour of optically bound particles in the Mie size regime

Optical binding refers to the light-induced organization and ordering of microparticles. In this paper, we present a theoretical treatment of optical binding in the Mie size regime using the paraxial approximation for field propagation combined with the Lorentz force law. Experimental studies of the separation between two dielectric spheres in a counter-propagating (CP) geometry are compared to numerical predictions. We explore the bistable nature of the bound sphere separation and its dependency on the refractive index mismatch between the spheres and the host medium, with an emphasis on the fibre separation and adiabaticity of the system.


Introduction
The formation of large-scale self organized colloidal matter created and mediated solely by light-matter interaction is a very evocative and potentially powerful methodology. The area is of widespread interest to researchers in colloidal science, optical micromanipulation, photonic crystals and even bioengineering. In particular, optical binding of microspheres has been of recent interest in the field of optical micromanipulation and it is reasonable to assert that many unanswered and key questions remain in this field. The topic has a rich history going back to the work of Burns et al [1] and has recently gained major momentum with observations in counter-propagating (CP) geometries [2]- [4] and even links to waveguiding in photonic crystals [5]. In the domain of optically confined microparticles, two different forms of scattering can be distinguished as the underlying principle of such self-consistent ordering of matter: firstly, orthogonal scattering [1] from trapped objects leads to array formation of arrays of Rayleigh particles (where the particle size is smaller than the wavelength) and the interparticle spacing is typically a fraction of the incident wavelength. Conversely, CP beam geometry has led to the formations of arrays of particles from a few 100 nm upwards and particularly in the Mie size regime (where the particle size is greater than the wavelength) and is distinguished by large inter-particle separations: the spacing observed is actually several times that of the physical dimensions of each bound microparticle. It is this regime we concentrate upon in the present study.
Optical binding differs significantly form conventional optical tweezers where a predefined optical landscape can be used to align and order microparticles in a predefined structure [6]- [10]. In these studies, optical landscapes may be formed using interferometric or holographic trapping geometries and trapping of each object occurs with the gradient (or dipole) force. These forms of colloid arrangement do not exploit, or consider, the redistribution of light by the matter itself, which will undoubtedly be the key to large scale ordered aggregation of microparticles. The process of redistribution of the optical field via diffractive focusing and refocusing in a particle chain lies at the heart of optical binding [11]. Interestingly the system implies a coupling effect as the spheres position alters the light distribution and vice versa. Thus the system is nonlinear and can exhibit effects such as bistability and hysteresis [3]. Bistability between an optically bound state and a collapsed array (with touching particles) was previously noted by Singer et al [2].
In this paper, we present an in depth treatment of the underlying theory involved in the formation of optical binding of two dielectric spheres in a dual beam fibre trap in the Mie size regime and we support this with experimental observations. While theory has briefly been mentioned in previous studies [3,11], the purpose of the present paper is to show an in depth discussion of the theory for binding and discuss comparison with new experimental data for the dependency of the binding phenomena on pertinent parameters such as the refractive index differences of the media and sphere and the fibre separation, which has never been explored. Optical binding is a nonlinear process and may lead to bistability.
As we vary the fibre separation, we observe a bifurcation in the allowed equilibrium solutions and bistability may be observed. We show how the fibre separation can bias the choice of the spheres separation. Finally we explore how the sphere separation follows any branch of the bistable solution and define what we mean by adiabaticity in this system. The paper is organized as follows: in section 2 we layout the theoretical basis for our studies followed by a description of the experimental setup in section 3. We then describe experiments 3 Institute of Physics ⌽ DEUTSCHE PHYSIKALISCHE GESELLSCHAFT that elucidate the various forms of behaviour we have noted above for our optically bound system in section 4, and finish with our summary and conclusions in section 5.

Theory
Our model comprises two monochromatic laser fields of frequency ω CP along the z-axis which interact with a system of N transparent dielectric spheres of mass m, refractive index n s , and radius R, with centres at positions { r j (t)}, j = 1, 2, . . . , N, and which are immersed in a host medium of refractive index n h . The monochromatic electric field is expressed as a sum of positive and negative frequency components as wherex is the unit polarization vector of the field, ε ± ( r ) are the slowly varying electric field amplitudes of the right or forward propagating (+) and left or backward propagating (−) fields, and k = n h ω/c is the wavevector of the field in the host medium. The incident fields are assumed to be collimated Gaussians at coordinates z = 0 for the forward field and z = D f for the backward field where r 2 = x 2 + y 2 , w 0 is the initial Gaussian spot size, P 0 is the input power in each beam, and D f is the distance between the fibre ends. It is assumed that all the spheres are contained between the beam waists within the length D f R, and that the initial Gaussian spot size is larger than the optical wavelength, w o > λ/n h , λ being the free-space wavelength, so that the paraxial approximation is applicable to propagation in the host medium.
Consider first that the dielectric spheres are in a fixed configuration at time t specified by the centres { r j (t)}. Then the dielectric spheres provide a spatially inhomogeneous refractive index distribution which can be written in the form where θ(R − | r − rj(t)|) is the heaviside step function which is unity within the sphere of radius R centred on r = r j (t), and zero outside, and n s is the refractive-index of the spheres. Then, following standard approaches [12], the CP fields evolve according to the paraxial wave equations along with the boundary conditions in equation (2), where k = ω/c and ∇ 2 ⊥ = ∂ 2 /∂x 2 + ∂ 2 /∂y 2 is the transverse Laplacian describing beam diffraction. Thus, a given configuration of the dielectric spheres modifies the fields ε ± ( r ) in a way that can be calculated from the above field equations.

4
Institute of Physics ⌽ DEUTSCHE PHYSIKALISCHE GESELLSCHAFT We remark that even though the spheres move, and hence so does the refractive-index distribution, the fields will always adiabatically slave to the instantaneous sphere configuration.
The paraxial wave equation (4) is to be solved subject to the boundary conditions (2), and we have done this using the split-step beam propagation method. We remark that our equations are valid in the Mie size regime where the sphere diameter is larger than the wavelength 2R > λ [13], and for input fields that are well approximated as paraxial. Furthermore, in keeping with the experimental conditions considered here the refractive-index difference between the sphere and host medium is small, n = (n s − n h ) < 0.1, meaning that there is negligible backscatter from the spheres, and the scattering is dominantly in the forward direction, the so-called Mie effect [14]. More specifically, the focal length for a sphere in the small-angle approximation neglecting higher-order aberrations is f = R/(2 n) [15], and since the sphere only focuses rays that pass through it within the sphere radius R away from the axis, we can introduce an effective numerical aperture for the sphere NA = sin(θ max ) = R/f = 2 n, with θ max the maximum ray deflection angle due to the sphere. Since n < 0.1 for the experiments and simulations presented here, the spheres act as low NA < 0.2 focusing elements that can be well treated using scalar paraxial theory. This is true since the maximum deflection angle θ max = 1 due to the sphere is small, so that initial paraxial rays will remain paraxial and the incident polarization state of the field will be mainly unchanged. Furthermore, this analysis also demonstrates that the paraxial approximation employed here improves with decreasing refractive-index difference, but that it is expected to fail for larger values n > 0.2 where scalar theory fails.
To proceed we need equations of motion for how the sphere centres { r j (t)} move in reaction to the fields which requires calculation of the forces F j acting on the spheres. Assuming overdamped motion of the spheres with viscous damping coefficient γ and neglecting hydrodynamic coupling for the static situation to be considered here, the equations of motion for the sphere centres can be written Carrying through full three-dimensional simulations with modelling of the electromagnetic propagation in the presence of many spheres poses a formidable challenge, so here we take advantage of the symmetry of the system to reduce the calculation involved. First, for the cylindrically symmetric Gaussian input beams used here we assume that the applied laser fields supply strong enough transverse confinement that the sphere motion remains directed along the z-axis. This means that the positions of the sphere centres are located along the z-axis r j (t) = zz j (t), and that we may also concentrate on the component of the forces along the F j =ẑF j .
To calculate the forces acting on the spheres, we follow the approach of Zakharian et al [16] and express the cycle-averaged force acting on the jth dielectric using the Lorentz force law where T is taken as a few times the period of the optical wave 2π/ω, V j signifies that the volume integral is taken over the volume of the jth dielectric sphere centred at r j (t), ρ b is the charge density arising from bound charges at the surface of the sphere, and J b is the current density arising from bound charges in the volume of the sphere. The first term in the integrand ρ b E gives rise to a force density acting along the direction of polarization of the laser fields and perpendicular to the z-axis. The cylindrical symmetry of our problem dictates that the time averaged force arising from the volume integral of the force density ρ b E must be zero, and we therefore concentrate on the second term which we will see is directed along the z-axis. For the system of dielectric spheres, we have the optical polarization P( r, t) = ε 0 (n 2 ( r) − n 2 h ) E( r, t) which is only nonzero within the spheres, and the current density due to bound charges is given explicitly by J b = ∂ P/∂t. Bringing these results together, we obtain for the force directed along the z-axis.
By expanding the electric and magnetic fields in equation (7) into positive and negative frequency components as in equation (1), using Faraday's law iω B( r, ω) = ∇ × E( r, ω), dropping fast oscillating terms at 2ω by virtue of the time integral, we obtain the force directed along the z-axis expressed in terms of the slowly-varying electric field envelopes in equation (1) as In deriving these equations, we have assumed that the two CP beams are mutually incoherent thereby neglecting any interference between them, and we have also used equation (4) to obtain expressions for ∂|ε ± | 2 /∂z . The lower expression in equation (8) is used in the simulations to numerically determine the force on each sphere, and we note that this force is a combination of both gradient and scattering forces in general. In previous work [17], we employed an expression for the force on each sphere as the spatial gradient of the interaction energy between the spheres and field, and although this expression is formally identical to equation (8), we have found that equation (8) is much more accurate in simulations. A similar expression to the first line of equation (8) was previously derived by Tlusty et al [18], and Chaumet and Nieto-Vesperinas [19]. At this point we specialize to the case of two spheres N = 2 to illustrate our numerical approach. By virtue of the symmetry of the applied laser fields we assume the spheres are symmetrically placed around z = D f /2 with a separation D, and label the sphere at z = (D f − D)/2 as sphere 1, and that at z = (D f + D)/2 as sphere 2. Then for a given sphere separation D, we calculate the CP fields between z = [0, D f ] using the beam propagation method applied to equation (4) with the boundary conditions (2), and from the fields we numerically calculate the force F 1 = −F 2 for each sphere. By calculating the CP fields for a variety of sphere separations we can numerically find the sphere separations, where the force on each sphere is zero. By plotting the force F 2 (D) acting on the sphere at z = (D f + D)/2 versus sphere separation D, we can determine the stability of the solution, stable equilibria having a negative slope ∂F 2 /∂D < 0, indicating a restoring force. From the calculated forces, we may also numerically determine the effective potential for the sphere motion The potential will prove useful in intuiting the sphere motions. Optical binding arises from the fact that the force acting on a given sphere is composed of two components along the z-axis, one from the laser field whose beam waist is closest to the given sphere and a second oppositely directed force arising from the CP laser field that is partly refocused on to the given sphere by the other sphere [11]. Balancing of these two forces by the refocusing of the spheres provides an intuitive explanation of optical binding. On the other hand, bistability in optical binding is possible in the coupled light-sphere system due to feedback [3]: changing the sphere separation alters the electromagnetic field distribution via the focusing properties of the spheres, which in turn alters the forces on the spheres. Due to this feedback, the forces on the spheres viewed as a function of sphere separation can become highly nonlinear, and give rise to bistability.

Experiment
A dual beam fibre optical trap [20] was used for our experimental studies (shown in figure 2) and discussed in detail elsewhere [2,3,11]. It was operated by light from a continuous wave Ytterbium Fibre Laser (IPG Photonics) at λ = 1070 nm which was coupled into two single mode fibres (F1 and F2) with equal field distribution to ensure symmetric array formation at half the fibre separation. By choosing a path length mismatch of the two fibres larger than the laser coherence length (∼1 mm), standing wave and interferometric effects were avoided [2,21]. To obtain a variable waist separation of the fields, one fibre (F2) could be displaced in z via a motorized micromanipulator, so a fibre separation (D f ) of between 40 to 140 µm could be obtained; additionally the separation could be cycled between these points with a preset speed (v).
A combination of the microscope and single beam tweezers allowed us to initiate and view simultaneously the array via a 60× microscope objective and a CCD camera. The camera was connected to a computer with frame grabber card to capture the images of the optically bound array with a sample rate of 25 frames s −1 .
Data analysis was performed by utilizing a Lab VIEW script with a pattern recognition algorithm to determine the position of the beads (error: better than ±0.5 µm) and the fibre separation (error: better than ±4%) within each frame of a captured movie. The overall mean values of the sphere separation were taken from up to 12 datasets of about 100 to 400 measurements each with typically three different sphere pairs. The results are shown as blue crosses in the following figures and the associated spread in the measurements are indicated by vertical bars. The associated error in our measurements is dependent on material parameters, namely, the refractive index n h of the host medium (with an accuracy of ±0.001), and the refractive index n s of the sphere for which measurements by Bangs laboratories indicated a spread of the refractive index of ±1%. Furthermore, absorption of water in the sphere can alter the refractive index [22] and therefore was estimated to be n s = 1.41 with errors of ±0.5%. Finally the sphere radius R has a size distribution (Bangs laboratories) such that the STD of the mean diameter for 3 µm is <10%.
The term n = (n s − n h ) denotes the refractive index difference between the sphere and the host medium. Silica microspheres (Duke Scientific n = 1.40 to 1.46@589 nm) were diluted in a de-ionized water and sucrose solution [2] to give a variable refractive index mismatch. The refractive index of the solution was measured with a Brix refractometer recalibrated for 1070 nm [23,24]. This enabled us to perform measurements for a variable n of between 0.03 to 0.09.
It is important to note that for low n, the optical forces acting are weak and the array creation (after loading by the helper tweezers) was observed after a timescale of several seconds before we achieved equilibrium separations. The low refractive index difference between the spheres and the host medium made an increase of the waiting time between the initiations of the array and measurement necessary due to a longer response time of the array. During measurements a fluctuation in the sphere separation may be caused due to the open sample cell used in our experiment, where a slight flow can perturb the array formation in contrast to a closed cell. These and other noise sources such as laser fluctuations (maximal 3% of output power) and external vibrations were compensated by long integration times of up to 1 min for each data set acquired. To account for evaporation of the host medium, which changes the concentration (by <4% over the duration of data acquisition of 5 min) and hence the refractive index, the host medium was washed off after each measurement and new medium and spheres added.
The fibre separation affects the distribution of the optical field and its variation can have significant influence on the bistable behaviour of an optical bound array, which is shown in the following simulations and experiments.

Results
For an array of two spheres, and fixed fibre separation D f and sphere refractive-index n s , the stable equilibrium sphere separations depend on the refractive-index difference n between the spheres and host medium. In the first set of experiments, the influence of n on the array spacing D for 5 µm and 1.28 µm diameter spheres with a fixed fibre spacing of D f = 90 µm was explored. The experimental data and comparison with our numerical modelling is shown in figure 3.
For the 5 µm spheres (left panel in figure 3), we see that the equilibrium separation decreases with n, whereas for the 1.28 µm spheres (right panel in figure 3), the equilibrium separation is seen to increase with n, both in agreement with the numerical modelling. Furthermore, for the smaller sphere size, the numerics predicts a bistable region at around n = 0.06. Experimentally the bistability could not be resolved for this example, but still a trend for the experimental sphere separation D to increase with n is evident.
An intuitive picture of the dependence of the equilibrium spacing D on the refractive index difference n is as follows: the small angle approximation to the focal length of the spheres obtained in section 2 is f = R/(2 n), which shows that the focal length of the sphere decreases with increasing refractive index difference. Since optical binding arises physically due to refocusing of the light fields by the two spheres on to each other, we expect the equilibrium sphere separations to follow the same trend as the sphere focal length, namely that it would decrease with increasing n. This trend is clearly seen in figure 3 for the 5 µm diameter spheres (left panel in figure 3), but not for the 1.28 µm diameter spheres (right panel in figure 3). We shall elucidate the reason why this difference can occur below.
The second set of experiments was designed to more fully explore the bistability in optical binding that we reported previously [3]. In particular, the stable equilibrium separation D of two 3 µm diameter silica spheres was measured as a function of the refractive index difference for fixed fibre separations of D f = 70, 90, 100 µm ±4%. The experimental results are shown in blue in the left-hand side panels of figures 4(a)-(c) respectively along with the corresponding numerical model results indicated by the red dots. In all of these experimental plots, the blue crosses represent the overall mean value for a series of realizations, and the blue vertical bars delineate the spread in measured values. For example, in figure 4(b) the overall mean values of the sphere separation was taken over on average 12 datasets of about 300 measurements each with typically three different sphere pairs. We see that there is good overall agreement between the experimental results and the stable solutions from the numerical model, the regions of negative slope in the D versus n plot being found to be unstable. In particular, the experiment shows the bistability predicted by the theoretical modelling.
For the fibre spacing D f = 70 µm in figure 4(a), we see that the sphere separation hovers around D = 11 µm, and a bifurcation point appears around n = 0.09 beyond which a new stable upper branch appears. This is also seen in the plot of the numerically generated effective potential U from equation (9) which is plotted on the right-hand side of figure 4(a) as a function of D and n. In particular, for n = 0.06 we see a global potential minimum at around D = 10 µm and marked (4), but for n > 0.09 two potential minima, marked (1) and (3), are evident (along with an unstable potential maximum marked (2)). Furthermore, we note that the lower branch, which exists below the bifurcation point, exhibits the expected trend that the sphere separation decreases with increasing n.
In contrast, for the fibre spacing D f = 90 µm in figure 4(b), we see that the bifurcation point is reduced to n ≈ 0.077, and both the upper and lower stable branches are equally evident in the sphere separation and potential plots. Note that whereas the sphere separation for the lower bistable branch has a tendency to decrease with n, the sphere separation for the upper bistable branch has a tendency to increase with n. Thus, in the vicinity of bistability of the optical binding the simple argument based on focusing that the sphere separation must decrease with increasing refractive index difference is negated, and this underlies the differences seen for the two cases in figure 3 with and without bistability.
As the fibre spacing is further increased to D f = 100 µm in figure 4(c), we see that the bifurcation point has increased to n ≈ 0.087, but in this case it is the lower branch that appears only beyond the bifurcation point (in contrast to the case in figure 4 for D f = 70 µm where it is the upper branch that only appears beyond the bifurcation point). The plot of sphere separation versus refractive index difference in this case is mainly dominated by the upper branch with the trend for D to now increase with n.
We also note that, as expected physically, the experimental fluctuations indicated by the blue vertical bars in figures 4(a)-(c), are largest closest to the bifurcation points where new bistable branches appear. Furthermore, it is also seen that the deviation between theory and experiment in figures 4(a)-(c) is largest for smaller values of n. This is understood by realizing that as the index mismatch decreases, the net optical forces acting on the spheres also get smaller, so the equilibria are created by cancellation of ever smaller forces due to the CP fields. In this situation, the numerical equilibria become more and more sensitive to the precise material parameters, whereas for larger index mismatches the equilibria are more robust against the precise parameters. In particular, for low values of the index difference the optical forces acting on the bead are smaller, and the corresponding effective potentials are shallower, meaning that the bead positions are more susceptible to systematic and environmental fluctuations in the experiment. This means that for lower values of the index difference, there is a less strong correlation between the theory and experiment in spite of the fact the paraxial approximation underlying our simulations is better for smaller index differences. Nevertheless, there is very good overall agreement between the numerical and experiment results in figures 4(a)-(c) and excellent evidence for bistability.
As indicated by the numerical simulations and experimental findings, the potential landscape of the bistability is strongly dependent on the fibre separation. In particular, the observed bistability in optical binding suggests that hysteresis in optical binding could be observed if the fibre separation was slowly cycled in the bistable regime so that the sphere separation would follow its local stable equilibrium value. For this experiment, the separation between the fibre ends was slowly cycled from 140 to 40 µm and back to 140 µm. For a relatively slow cycle velocity of the fibre ends, 1 µm s −1 , the spheres were experimentally found to adiabatically follow the changing fibre separation in the system and hysteresis was observed. This is shown as the blue data in figure 5 (left panel) for which the experimentally measured sphere separation D is plotted parametrically as a function of the fibre separation D f . The red data in figure 5 (left panel) shows the theoretical sphere separation as a function of fibre separation, and there is excellent overall agreement between the theoretical calculations and experiment.
In contrast, if the fibre spacing is changed too quickly, the hysteresis loop becomes washed out and eventually vanishes. This is illustrated in figure 6 (left panel) where we plot the  experimentally measured sphere separation parametrically as a function of fibre separation for velocities v black = 30 µm s −1 , v blue = 10 µm s −1 , and v red = 4 µm s −1 , as D f is cycled between 60-100 µm. If the speed is reduced to v red = 2 µm s −1 , as in the right panel in figure 6, the hysteresis was observed to extend to 130 µm, close to the theoretical value in figure 5. We remark, however, that even with a fixed pair of spheres and fibre velocity the 'switch up' point in the bistable loop could vary from shot to shot, but the 'switch down points' were much more robust. This is shown in the blue and black curves in the right-hand side panel in figure 6, where a fibre velocity of 1 µm s −1 was used. Such shot-to-shot variation is associated with 'critical slowing down' in the vicinity of switching points where the effective potential U is flat and the force is close to zero, and the system is therefore more prone to the effects of noise. Figure 5 (left panel) shows U as a function of sphere spacing D and fibre spacing D f , and we see that around the upper switch point D f ≈ 110 µm the potential is relatively flat as a function of D, so we expect shot-to-shot variations. In contrast, near the switch down point D f ≈ 60 µm the potential has much more structure as a function of D, and we expect that the switch down dynamics should be much more deterministic as seen in the experiments.

Summary and conclusions
We have presented an in-depth experimental and numerical investigation of the dependency of the optical binding of two spheres as a function of the refractive index mismatch and the fibre separation with an emphasis on the bistable nature of the system. Physically, bistability can occur in optical binding as the light is modified (refocused) by the spheres and the light in turn tells the spheres how to move. It is this feedback that allows the force acting on the spheres viewed as a function of sphere separation to become nonlinear and to display several zero crossings (see above figures), and thus to show bistability. Our model is applicable directly to sphere sizes in the Mie and Lorentz-Mie regime and shows good agreement with experiment. Optical binding is in many respects in its infancy and is likely to play a pivotal role in self assembly of crystalline structures using optical forces.