2. Description of the detection scheme
The method relies on the ability of high-Q resonators to simultaneously optically trap [
21
K. Svoboda and S. M. Block, “Optical trapping of metallic Rayleigh particles,” Opt. Lett.
19(13), 930–932 (1994). [CrossRef]
[PubMed]
–
27
A. H. J. Yang, S. D. Moore, B. S. Schmidt, M. Klug, M. Lipson, and D. Erickson, “Optical manipulation of nanoparticles and biomolecules in sub-wavelength slot waveguides,” Nature
457(7225), 71–75 (2009). [CrossRef]
[PubMed]
] and detect the presence of NPs. Electrical field gradients caused by the evanescent fields of waveguides have been shown to allow trapping of small particles [
24
L. N. Ng, B. J. Luff, M. N. Zervas, and J. S. Wilkinson, “Propulsion of gold nanoparticles on optical waveguides,” Opt. Commun.
208(1-3), 117–124 (2002). [CrossRef]
,
26
T. J. Davis, “Brownian diffusion of nano-particles in optical traps,” Opt. Express
15(5), 2702–2712 (2007). [CrossRef]
[PubMed]
,
27
A. H. J. Yang, S. D. Moore, B. S. Schmidt, M. Klug, M. Lipson, and D. Erickson, “Optical manipulation of nanoparticles and biomolecules in sub-wavelength slot waveguides,” Nature
457(7225), 71–75 (2009). [CrossRef]
[PubMed]
]. In order to effectively trap a particle, the energy well associated to the trap has to be deep enough to overcome Brownian motion, i.e., the energy well has to be at least on the order of several
kBT, where
kB
is Boltzmann’s constant and
T is the temperature. However, even when field strengths are relatively modest, such that single NPs are only marginally trapped, NP
clusters are trapped for long periods of time once a critical NP cluster size is reached. This is due to the relative scaling of the trapping energy to the thermal energy associated with Brownian motion. Indeed, the thermal energy associated to Brownian motion is 3
kBT/2 irrespectively of the cluster size: The additional degrees of freedom for large clusters are internal; thermal energy associated to the movement of the center of mass remains unchanged. The trapping energy of the optical well on the other hand scales with the size of the cluster, since the trapping force is exerted on each of the aggregated NPs. As a result, the escape time of a NP cluster (i.e., the mean residence time inside the trap) scales as exp(|
NU
0|/
kBT) [
23M. Pelton, M. Liu, H. Y. Kim, S. Glenna, P. Guyot-Sionnest, and N. F. Scherer, “Optical trapping and alignment of single gold nanorods using plasmon resonances,” Proc. SPIE 6323, 63230E 1–9 (2006).
,
26
T. J. Davis, “Brownian diffusion of nano-particles in optical traps,” Opt. Express
15(5), 2702–2712 (2007). [CrossRef]
[PubMed]
], where
U
0 is the trapping energy for a single NP and
N is the number of NPs inside the cluster. It can be significantly larger for a cluster than for a single NP. For example, if the trapping energy of a single NP (
U
0) is 6
kBT, the associated escape time remains small, on the order of 50 ms in a typical ridge waveguide trap for 20 nm radius NPs (precisely, this is the “thermalization time” defined later in this section). The escape time of a two-particle cluster is larger by more than two orders of magnitude (~exp(6)) and is on the order of 10 seconds. This results in a large accumulation of NP dimers inside the optical trap, and a significant increase of NP concentration in the evanescent field of the resonator. This in turns modifies the refractive index of the solution around the resonator and induces a resonant frequency shift. In an assay where cluster aggregation is uniquely triggered by the presence of target molecules, this frequency shift is a unique signature of the target.
Figure 1
illustrates this scheme with the example of a slot waveguide ring resonator [
27
A. H. J. Yang, S. D. Moore, B. S. Schmidt, M. Klug, M. Lipson, and D. Erickson, “Optical manipulation of nanoparticles and biomolecules in sub-wavelength slot waveguides,” Nature
457(7225), 71–75 (2009). [CrossRef]
[PubMed]
,
28
T. Baehr-Jones, M. Hochberg, C. Walker, and A. Scherer, “High-Q optical resonators in silicon-on-insulator-based slot waveguides,” Appl. Phys. Lett.
86(8), 081101 1-3 (2005). [CrossRef]
] (slot waveguides are not essential, since trapping can also be obtained with other waveguide geometries such as ridge waveguides [
24
L. N. Ng, B. J. Luff, M. N. Zervas, and J. S. Wilkinson, “Propulsion of gold nanoparticles on optical waveguides,” Opt. Commun.
208(1-3), 117–124 (2002). [CrossRef]
,
26
T. J. Davis, “Brownian diffusion of nano-particles in optical traps,” Opt. Express
15(5), 2702–2712 (2007). [CrossRef]
[PubMed]
]).
Fig. 1 Waveguide coupled slot waveguide ring resonator capturing gold NP clusters after target molecule induced aggregation. Single NPs are only marginally trapped, NP accumulation in the evanescent field of the ring only occurs in the presence of clusters.
Ring resonators provide a practical technology platform for high-Q resonators for biodetection [
11
M. Iqbal, M. A. Gleeson, B. Spaugh, F. Tybor, W. G. Gunn, M. Hochberg, T. Baehr-Jones, R. C. Bailey, and L. C. Gunn, “Label-Free Biosensor Arrays Based on Silicon Ring Resonators and High-Speed Optical Scanning Instrumentation,” IEEE J. Sel. Top. Quantum Electron.
16(3), 654–661 (2010). [CrossRef]
]. As shown in
Fig. 1, a bus waveguide is coupled to a ring. The coupling strength is chosen such that light is fully dropped from the waveguide bus at resonance. This is the so-called critical coupling condition where ring coupling losses equal the roundtrip losses [
29
J. M. Choi, R. K. Lee, and A. Yariv, “Control of critical coupling in a ring resonator-fiber configuration: application to wavelength-selective switching, modulation, amplification, and oscillation,” Opt. Lett.
26(16), 1236–1238 (2001). [CrossRef]
]. Different types of waveguide geometries can be used for the ring itself, such as slot waveguides [
27
A. H. J. Yang, S. D. Moore, B. S. Schmidt, M. Klug, M. Lipson, and D. Erickson, “Optical manipulation of nanoparticles and biomolecules in sub-wavelength slot waveguides,” Nature
457(7225), 71–75 (2009). [CrossRef]
[PubMed]
] (
Fig. 2(a)
and
2(b)) or ridge waveguides [
24
L. N. Ng, B. J. Luff, M. N. Zervas, and J. S. Wilkinson, “Propulsion of gold nanoparticles on optical waveguides,” Opt. Commun.
208(1-3), 117–124 (2002). [CrossRef]
,
26
T. J. Davis, “Brownian diffusion of nano-particles in optical traps,” Opt. Express
15(5), 2702–2712 (2007). [CrossRef]
[PubMed]
] (
Fig. 2(c) and
2(d)). The waveguides shown in
Fig. 2 are assumed to be etched in 200 nm thick silicon-on-insulator (SOI) material. The slot is 200 nm wide with 250 nm wide silicon strips on either side. The ridge waveguide is 400 nm wide. The trapping energy distribution for single NPs (
U
0) is shown for waveguides carrying 2.5 mW of power and for 20 nm radius NPs, the default NP size assumed in this paper. Due to the high overlap of the investigated silicon waveguides with the solution, quality factors (
Q) and resonator finesse (
F) are limited by water absorption and are on the order of
Q = 10,000 (loaded
Q, i.e. including coupling losses to the waveguide) and
F = 30 (for a 50 μm radius) at 1580 nm. The quality factor is a measure of the resonance linewidth, and is an important metric to determine the minimum wavelength shift that can be detected. Typical trapping strengths on the order of
U
0 = 6
kBT are assumed for single NPs during numerical investigation of assay sensitivity (section 6). These trapping energy levels can be obtained with modest optical power: The finesse corresponds to the power ratio of the resonator waveguide to the bus waveguide, so that a trap strength of
U
0 =6
kBT can be reached for a bus waveguide power of 0.5 mW (0.5[mW]×
F/2.5[mW/
kBT]= 6
kBT).
Fig. 2 Waveguide cross-section and trapping energy (in units of kBT) for (a) a slot waveguide and (c) a ridge waveguide. (b) and (d) show the corresponding trapping energy through the center plane of the waveguides (x = 0). The mode in (a) is x-polarized, while the mode in (c) is y-polarized in order to maximize the E-field strengths in the solution.
While this form of assay is conceptually relatively straightforward, modeling its sensitivity limit is a complex problem. In particular, practical optical trapping strengths are limited by both thermal and temporal considerations. The optical field leads to heating, both due to direct absorption by the NPs, as well as absorption by water. Since excessive temperature leads to melting of the NP aggregates [
12
R. Elghanian, J. J. Storhoff, R. C. Mucic, R. L. Letsinger, and C. A. Mirkin, “Selective colorimetric detection of polynucleotides based on the distance-dependent optical properties of gold nanoparticles,” Science
277(5329), 1078–1081 (1997). [CrossRef]
[PubMed]
,
13
J. J. Storhoff, A. D. Lucas, G. Viswanadham, Y. P. Bao, and U. R. Müller, “Homogeneous detection of unamplified genomic DNA sequences based on colorimetric scatter of gold nanoparticle probes,” Nat. Biotechnol.
22(7), 883–887 (2004). [CrossRef]
[PubMed]
,
18
J.-S. Lee, M. S. Han, and C. A. Mirkin, “Colorimetric Detection of Mercuric Ion (Hg2+) in Aqueous Media using DNA-Functionalized Gold Nanoparticles,” Angew. Chem.
119(22), 4171–4174 (2007). [CrossRef]
] this is an important constraint. Also, the time required for the concentration of a cluster species (indexed by the NP number
N) to reach equilibrium in the trap (the “trap loading time”, related to the escape time) grows exponentially with the trapping energy
NU0
. In practice the time delay between turning on the optical field and reading out the result of the assay (the “assay stabilization time”) is finite, so that only concentrations of smaller cluster species are stabilized at readout, while the concentrations of larger clusters remains below their theoretical maximum. The mean residence time of NP clusters inside the trap (specifically, the “cluster thermalization time” as defined below) plays an important role for a second reason: The “cluster dissociation time” (i.e. the mean time before a cluster decomposes into smaller components due to target-NP bond dissociation) reduces the NP cluster concentration enhancement inside the optical trap if it is smaller than, or on the order of, the cluster thermalization time. This is due to the fact that equilibrium cluster concentrations assumed in this paper in the optical trap are calculated assuming that cluster capture and accumulation is only competing with cluster escape directly mediated by Brownian motion. If clusters can also dissociate (and subsequently escape due to reduced trapping strength), a second mechanism competes with cluster accumulation. However, if the thermalization time is smaller than the cluster dissociation time, clusters do not reside long enough inside the trap for dissociation to play a significant role. This is the limit in which reported numerical results are accurate, and is a limit verified by a variety of experimental conditions as discussed in section 3.
The outline of the paper is as follows: Section 3 contains a brief discussion of cluster dissociation times for several specific assay chemistries. A steady state model for the assay sensitivity is derived in section 4. This model is expanded to encompass time dynamics in section 5, a necessary step to obtain meaningful assay sensitivities, since they are limited by the assay stabilization time. In the interest of focusing on the actual assay sensitivity, the trap loading time is quickly estimated in section 5 based on a simplified model. However, a more accurate derivation is later provided in section 7. Section 6 is dedicated to a numerical investigation of the assay sensitivity. Section 7 derives the particle escape time based on the formalism introduced by Kramers [
30
H. A. Kramers, “Brownian motion in a field of force and the diffusion model of chemical reactions,” Physica
7(4), 284–304 (1940). [CrossRef]
]. In particular, trap loading times are derived in the 2D (silicon wire, or large radii rings) and 3D (smaller rings) limits. The experimental feasibility under the point of view of thermal constraints is verified in section 8.
In order to fully justify the time constants used in numerical evaluations, a complete theory of NP trapping dynamics is reported in section 7. Kramers derived in 1940 the lifetime of a particle captured in a one-dimensional (1D) potential well [
30
H. A. Kramers, “Brownian motion in a field of force and the diffusion model of chemical reactions,” Physica
7(4), 284–304 (1940). [CrossRef]
] (the Kramers escape time), albeit for a trap with a well-defined size, defined by a ridge or saddle point (
Fig. 3
). As previously pointed out [
26
T. J. Davis, “Brownian diffusion of nano-particles in optical traps,” Opt. Express
15(5), 2702–2712 (2007). [CrossRef]
[PubMed]
], the particle escape problem is different in the case of an evanescent optical field, as the potential energy asymptotically decays towards infinity: There is no clearly defined trap size (there is no saddle point). This raises important challenges in the theoretical determination of the particle escape time. If one defines particle escape as reaching a predetermined position (the trap boundary), the particle escape time is highly dependent on the position of the trap boundary, and does not converge when the distance from the center of the trap to the boundary is increased. In previous theoretical studies of particle escape times in a 1D evanescent optical trap [
26
T. J. Davis, “Brownian diffusion of nano-particles in optical traps,” Opt. Express
15(5), 2702–2712 (2007). [CrossRef]
[PubMed]
], the trap boundary was set at a point of predetermined potential energy level, on the order of
kBT. This boundary is justified by the fact that particle dynamics are dominated by Brownian motion (“free particles”) once the trapping energy is substantially less than the thermal kinetic energy, hence the particle is considered as having escaped past that point. This is an adequate escape criterion when comparing the cluster escape time to the cluster dissociation time (i.e., when calculating the “cluster thermalization time”): Indeed once a cluster reaches a trapping energy level small relative to
kBT, it reaches a point in the solution where,
at steady state, its ambient concentration corresponds to its bulk solution concentration (the “bulk concentration”) and is thus in thermodynamic equilibrium with smaller species. At this point cluster dissociation is already taken into account by the determination of the bulk concentration. However this escape criterion is not adequate to calculate the trap loading time: Indeed, during
transient trap loading the cluster concentration in the vicinity of the trap can get depleted below the bulk concentration value, slowing down the cluster capture rate. In other words, setting a trap boundary at a point of small trapping energy relative to
kBT is not sufficient to assume that the cluster concentration at this boundary is maintained at the bulk solution level during transient loading.
Fig. 3 Well shape used in the classic Kramers escape problem.
In the following, “thermalization time” is used specifically to refer to the mean time a particle remains inside the trap before returning to a point in the solution where its optical trapping energy is low (derived for an energy level of
kBT/4 in section 7) and where its concentration is close to the bulk concentration
in steady state conditions. The “trap loading time” refers to the time constant of
transient particle accumulation after the trap has been turned on. While the two are closely related, the trap loading time can exceed the thermalization time by one or several orders of magnitude and scales differently with the waveguide dimensions than the previously derived thermalization time [
26
T. J. Davis, “Brownian diffusion of nano-particles in optical traps,” Opt. Express
15(5), 2702–2712 (2007). [CrossRef]
[PubMed]
]. We will show that the particle escape time derived with Kramers’ formalism for a 3D evanescent trap scales simply as the sum of a well-defined trap loading time and a second term corresponding to the time needed for reaching the boundary by random-walk, while in lower dimensional systems the scaling law is more complicated, and transient loading effects have to be taken into account to derive a physical loading time. This is due to a fundamental property of random walks in 3D, in that the probability of a particle ever returning to the origin is smaller than 1, so that “particle escape” is a well defined condition even in the absence of a saddle point, while this is not the case in lower dimensional systems.
In order to model the assay response we introduce a cluster species dependent generalized overlap integral describing the overlap of the resonator field with the solution while taking into account the spatial variation of the cluster concentration throughout the optical trap. It is shown that the non-exponential (non-activated) energy dependence of the trap loading time is fundamentally linked to the trapping energy dependence of the generalized overlap, and that this relation is required to obtain physically meaningful results when evaluating the assay response.
4. Static model
In this section, a static expression of the effective index change seen by the resonator in the presence of NPs and targeted biomolecules is derived. Section 5 is dedicated to time varying aspects such as the finite trap loading time, and their impact on the assay sensitivity. The NP concentration near the surface of the resonator is a function of both the bulk solution concentrations of individual cluster species, determined by NP aggregation dynamics, as well as of the cluster concentration enhancements near the resonator as induced by cluster trapping. The induced shift in refractive index is then calculated based on a generalized overlap integral between the waveguide and the surrounding solution, that also takes into account the spatial variation of the concentration enhancements. The concentration enhancements are not instantaneous, but require a finite time to stabilize. This is taken into account in section 5 to obtain expressions for a realistic situation with finite stabilization times.
Optical trapping of small NPs relies on the dipole force induced in an electrical field gradient. The corresponding energy potential is given by
where
ε0
is the permittivity of free space,
εw
and
εg
are the relative dielectric constants of water (1.77) and gold (−90+10i at 1580 nm) [
33D. W. Lynch and W. R. Hunter, Handbook of Optical Constants of Solids, ed. Palik, E., Academic Press, Orlando, FL, 286–295 (1985).
],
N is the number of NPs in the cluster,
V the volume of a single NP,
the average squared electrical field, and
U0
the trapping energy of a single NP (defined as a positive number).
α is defined by the left portion of
Eq. (1) and is simply a compact notation for the expression inside the parenthesis. The optical trap leads to a local enhancement of the cluster concentration. Assuming that the cluster dissociation time is much larger than the cluster thermalization time, i.e., assuming that cluster dissociation time can be ignored, and assuming that the clusters are small enough to fit inside the trap, the equilibrium cluster concentration around the resonator is straightforward to derive from the diffusion equation
where
c(
N) is the concentration,
γ(
N) is the coefficient of friction and
D(
N)=
kBT/
γ(
N) is the diffusion constant, all for clusters of size
N. This results in a Boltzmann type distribution, with
at the bottom of the well, where
c
0(
N) is the concentration of clusters of size
N in the bulk of the solution (far from the trap). For
N > 1, it can be shown to scale as [
34
N. J. Lynch, P. K. Kilpatrick, and R. G. Carbonell, “Aggregation of ligand-modified liposomes by specific interactions with proteins. I: Biotinylated liposomes and avidin,” Biotechnol. Bioeng.
50(2), 151–168 (1996). [CrossRef]
[PubMed]
]
where
cNP
is the initial NP concentration prior to aggregation,
ϕ is a parameter describing the cluster formation statistics called the extent of aggregation [
34
N. J. Lynch, P. K. Kilpatrick, and R. G. Carbonell, “Aggregation of ligand-modified liposomes by specific interactions with proteins. I: Biotinylated liposomes and avidin,” Biotechnol. Bioeng.
50(2), 151–168 (1996). [CrossRef]
[PubMed]
],
υ(
Ν) is a coefficient in the Taylor expansion whose value depends on the type of cluster formation dynamics (e.g. diffusion or reaction limited [
34
N. J. Lynch, P. K. Kilpatrick, and R. G. Carbonell, “Aggregation of ligand-modified liposomes by specific interactions with proteins. I: Biotinylated liposomes and avidin,” Biotechnol. Bioeng.
50(2), 151–168 (1996). [CrossRef]
[PubMed]
,
35
S. Y. Park, J.-S. Lee, D. Georganopoulou, C. A. Mirkin, and G. C. Schatz, “Structures of DNA-linked nanoparticle aggregates,” J. Phys. Chem. B
110(25), 12673–12681 (2006). [CrossRef]
[PubMed]
]) and
O(
ϕN
) is a term of order
ϕN
. Physically,
ϕ is the number of bonds per particle if only 1 bond between any two particles is counted and cyclic connections are not allowed. In an extreme case where each target participates to a unique NP-NP bond
ϕ =
cT
/
cNP
, where
cT
is the initial target concentration. In general,
ϕ <
cT
/
cNP
since some targets contribute to redundant bonds (i.e., they bind a NP to a cluster to which it is already attached via another bond), or they might not participate to a bond at all (due to steric effects or due to dissociation). The specific relationship between
ϕ and
cT
/
cNP
is dependent on the assay chemistry, so that numerical results in section 6 are reported as a function of
cNP
and of
ϕ. Assuming
ϕ =
cT
/
cNP
gives a lower (best case) sensitivity for the assay and is a good approximation for dilute target concentration (small
ϕ, so that clusters are dominated by dimers and steric effects and cyclic or redundant bonds can be ignored) and long NP-NP bond lifetimes (longer than the particle collision time, so that all targets can be assumed to participate to a bond). In general,
υ(
Ν) depends on the cluster formation process, such as diffusion limited cluster formation (each NP collision leads to a bond formation) or reaction limited cluster formation (many collisions are required, probability of bond formation depends on cluster sizes). However
υ(
N) is close to 1 for
N = 1, 2 and 3 for both types of cluster formation dynamics. Precisely,
υ(
N) = 1 for all
N in the case of diffusion limited cluster formation, and
υ(
Ν) is respectively 1, 1 and 1.5 for
N = 1, 2 and 3 in the case of reaction limited cluster formation [
34
N. J. Lynch, P. K. Kilpatrick, and R. G. Carbonell, “Aggregation of ligand-modified liposomes by specific interactions with proteins. I: Biotinylated liposomes and avidin,” Biotechnol. Bioeng.
50(2), 151–168 (1996). [CrossRef]
[PubMed]
]. Since in the investigated regime the assay response is dominated by these three smallest cluster species (due to time domain limitations, as shown in section 5), the assay response is very close for both types of cluster formation dynamics at small and moderate
ϕ. The Taylor expansion,
Eq. (4), is only accurate for small
ϕ, the typical operating point since detection of low target concentrations is investigated. Complete analytical expressions derived from Smoluchowski kinetics [
34
N. J. Lynch, P. K. Kilpatrick, and R. G. Carbonell, “Aggregation of ligand-modified liposomes by specific interactions with proteins. I: Biotinylated liposomes and avidin,” Biotechnol. Bioeng.
50(2), 151–168 (1996). [CrossRef]
[PubMed]
] will be used to predict trends at large target concentrations in section 6. Here, decomposition in a Taylor series allows explicitly showing trends over
ϕ and
U
0.
By applying the Maxwell-Garnett formula [
36
C. R. Snyder and J. F. J. Douglas, “Determination of the Dielectric Constant of Nanoparticles. 1. Dielectric Measurements of Buckminsterfullerene Solutions,” Phys. Chem. B
104(47), 11058–11065 (2000). [CrossRef]
] for the dielectric constant of a solution with suspended particles, the refractive index in the vicinity of the resonator can be derived as
where
nA
is the Avogadro constant and
nw
the refractive index of water.
α has been defined in
Eq. (1). The first approximation is made by taking a first order approximation of the square root (always valid since index changes are small) while the second is obtained by replacing the cluster concentrations by its first order approximation, a step that is only valid for small
ϕ. It can be seen that an increasing number of terms contribute to the index change when
ϕ⋅exp(
U
0/
kBT) is approaching 1, thus potentially triggering a highly nonlinear response. As will be seen in the following, this nonlinearity is reduced by truncation of the series due to time dynamics and relatively linear responses can be obtained for the right choices of assay conditions (section 6).
Before converting
Eq. (5) into a wavelength shift, several corrective factors have to be taken into account. First, the cluster concentration enhancement does not occur to the same degree throughout the entire optical trap, since the depth of the energy well varies spatially. The refractive index change of the solution can be converted into a waveguide effective index change by means of an optical overlap integral [
37
J. Witzens, T. Baehr-Jones, and M. Hochberg, “Design of transmission line driven slot waveguide Mach-Zehnder interferometers and application to analog optical links,” Opt. Express
18(16), 16902–16928 (2010). [CrossRef]
[PubMed]
]. Here, the spatial variation of the concentration throughout the evanescent field of the waveguide is taken into account by a generalized overlap integral defined as
where Δ
neff
(
N) is the contribution of clusters of size
N to the effective index change of the waveguide, Γ is the generalized optical overlap (a function of
NU
0 /
kBT) and Δ
n
max(
N) is the refractive index change of the solution induced by clusters of size
N at the point of maximum cluster concentration, i.e., the bottom of the energy well. The generalized overlap is given by
where
H is the magnetic H-field and
Z
0 the impedance of free space. The exponential inside the integral is an addition to the usual overlap integral that takes into account the spatial dependence of the cluster concentration,
c(
x,
y)/max(
c). The final expression for the effective index change of the ring waveguide is
5. Time resolved model
For diffusion limited cluster formation,
υ(
N) = 1 for all
N.
Equation (5) is then diverging for
ϕ⋅exp(
U
0/
kBT) > 1 and obviously non-physical when time dynamics are ignored. In fact, the trapping energy and thus the trap loading time both increase for larger clusters, so that the concentrations of large cluster species do not have time to stabilize for finite assay stabilization times. The time-dependent effective index change can be expressed as
where
t is the assay stabilization time and
τload
(
N) is the trap loading time-constant for clusters of size
N. The latter can be estimated in the framework of a simple rate equation describing trapped particle dynamics (results will be rigorously derived using the framework of Kramers’ escape time derivation in section 7).
We consider a particle trap with a given cluster species dependent trap loading time
τload
(
N) and a cumulative cluster capture rate (2
πRw)
η such that the flux of captured clusters is given by (2
πRw)
ηc
0, where
R is the radius of the ring resonator,
w the width of the waveguide trap (i.e., the slot width or the ridge waveguide width),
η the cluster capture rate per trap area and
c
0 the bulk solution cluster concentration. Once the trap has been turned on the cluster capture dynamics are described by the rate equation
where
M is the number of clusters inside the trap. The cluster number at equilibrium
Meq
= (2
πRw)
τload
ηc0
and can be readily calculated by integrating the cluster concentration (
Eq. 3) over the optical trap
In
Eq. (11),
h is the physical height of the trap (depending on where the trap boundary is placed) and
htrap
(
NU
0/
kBT) is an effective trap height, such that the
excess number of particles captured after turning on the trap corresponds to the concentration at the bottom of the well times an effective volume 2
πRwhtrap
and is a well defined integral. For a natural trap boundary, i.e. relatively close to the waveguide, shortly after decay of the evanescent field, the number of clusters inside the trap is dominated by the excess concentration, so that
In order to quickly estimate η, we assume that it corresponds to a diffusion limited particle flow induced by placing an ideal particle sink in the solution. At t=0 the cluster concentration inside the trap (c
0) is several orders of magnitude smaller than the equilibrium concentration (~c
0exp(NU
0/kBT)) so that the trap is simply modeled as a particle sink with a zero concentration boundary condition. (2πRwη)c
0 is estimated as the steady state particle flow.
Figure 4(a)
shows the resulting concentration distribution as simulated in cylindrical coordinates. This is a well defined problem in a 3D environment, since given a fixed flux 2
πDr
0
c
0 (integrated over a solid angle 2
π) the concentration scales as
c
0(1-
r
0/
r) and converges to
c
0 far from the particle sink, where
r is the distance to the sink. Hence, the estimated particle flux converges to a non-zero value as the simulation domain size is increased (this is not the case in lower dimensional systems).
Fig. 4 Concentration distribution induced by a ring shaped particle sink with a 400 nm width and a 20 μm radius, as simulated in cylindrical coordinates. (b) Extracted value of the coefficient I as a function of the domain size d. (c) Coefficient I as a function of the estimated number of accumulated particles. The blue curve corresponds to the final number of accumulated particles for a 400 nm wide ridge waveguide (R = 20 μm) at NU
0 = 12.5 kBT. The red line corresponds to 10% of that number.
The total particle flux captured by the sink (
2πRw)
ηc
0 scales as
DR if
w/
R is kept constant (i.e. if the entire geometry of the particle sink is rescaled) since the solution to the diffusion equation can be simply rescaled. Hence it is natural to describe the particle capture rate as a function of
D/
w
where
I is a dimensionless constant that depends on
w/
R (and whose complete dependence on trap geometry is derived in section 7). This results in following expression for the trap loading time
In the following, the simulation domain is defined by constant concentration boundary conditions located at a distance
d above the sink, and at a distance
d to the left of the sink (i.e. the total domain size is given by (20 μm +
d) ×
d).
Figure 4(b) shows the estimated
I as a function of the simulation domain size for a 400 nm wide trap (the ridge waveguide) and a 20 μm ring radius. It can be seen that
I converges to
I
lim = 2.3 for an infinitely far constant concentration boundary. For a 200 nm wide trap (the slot) and a 20 μm ring radius
I converges to
I
lim = 2.5.
Three approximations entered the estimate of I, all of which correspond to an upper bound of the trap loading time (and thus to a pessimistic assessment of the assay sensitivity). The most important is the absence of fluid movement (i.e., a fully diffusion limited particle transport), a condition that might be hard to replicate in experimental conditions, since one would expect at least convection currents to be generated due to optically induced heating. The other two approximations are addressed in section 7 (and only result in a small overestimate here, on the order of 30% in the investigated regime).
One assumption is that the transient particle flux during trap loading equals the particle flux obtained with the steady state simulations under the effect of a permanent sink and infinitely far constant concentration boundary conditions.
In order for the particle flux to drop to the value predicted by a steady state simulation, a sufficient number of particles has first to be absorbed by the trap in order to deplete the surrounding solution down to the level seen in the steady state simulation. If this particle count is small relative to the total number of particles eventually accumulated by the trap, the trap loading time is accurately derived from the steady state particle flux, since the large majority of absorbed particles are limited by that flux. On the other hand, if this is not the case the trap loading will occur faster than predicted by the steady state particle flux, and is also non-exponential (since the incoming particle flux is time varying). A more accurate value of the particle flux can be estimated by setting the distance to the boundary conditions in the steady state simulation in such a way as to equate the cumulative particle depletion (the integral of the concentration depletion extracted from the steady state simulation) with the total number of particles accumulated by the trap (
Fig. 4(c)). It can be seen in
Fig. 4(c) that for the investigated geometry, 90% of the absorbed particles for a ridge waveguide trap with
NU
0 = 12.5
kBT are transported with a particle flux corresponding to
I > 1.9. In other words,
I
lim is fairly accurate for trapping energies larger than 12.5
kBT, while it corresponds to an upper bound of the loading time at smaller trapping energies. However, it should be noted that even though in general
I
lim is an upper bound at small energies, it is in practice an extremely robust estimate for order of magnitude calculations: As will be shown in section 7,
I reaches
half of its high energy limit
I
lim if
NU
0 is reduced down to 6
kBT, the latter corresponding to a
two orders of magnitude reduction of the exponentially energy activated term of the cluster concentration enhancement in
Eq. (14).
The second approximation consists in ignoring the effect of the trapping field on the particle flux, i.e., particles are considered to be “free” particles even close to the waveguide surface. The effect of the trapping field is taken into account in section 7 and results in a slight reduction of I
lim (corrected value ~1.9). This too is only a relatively small correction.
h
trap(
N) is dependent on the exact shape of the optical trap, and hence on the type of waveguide. For a ridge waveguide, the shape of the optical trap is given by
where
δ is the folding length of the evanescent field
intensity and is close to 100 nm for both the ridge and the slot waveguides. For the slot waveguide, the energy profile is locally parabolic around the bottom of the trap. The characteristic length scale for the center of the trap is defined as
where the second derivative is taken at the bottom of the well. For the slot waveguide shown in
Fig. 2(b),
δ’ = 320 nm.
htrap
can then be derived as
Ridge WG:
Slot WG:
In
Eq. (18), the integral was taken from
-∞ to ∞ since the cluster concentration is already very small by the time the trap gets truncated by the wafer substrate (e.g., the error is 5% for
NU
0 = 6
kBT and smaller at larger trapping energies). It can be seen that at equal well depth, a slot waveguide can capture more particles (
δ’ >
δ and the energy dependence scales as 1/sqrt(
NU
0) rather than 1/(
NU
0)). The latter is due to the fact that it has a flat bottom at the energy minimum. For these particular waveguide traps, the expression for the trap loading time can be derived as
where
δ’ has already been defined for the slot waveguide, and
δ’=
δ for the ridge waveguide. χ = 1 for the ridge waveguide and ½ for the slot waveguide. It should be noted that
Eq. (19) resulted from deriving
wh
trap from a 1D model. In the numerical examples given in section 6,
wh
trap is replaced by a 2D integral
S
trap calculated based on the 2D mode profile in order to increase the numerical accuracy.
Equation (19) shows that the non-activated (non-exponential) energy dependence of the cluster escape time is given by
S
trap/
wδ' ~
htrap
(
N)/
δ', an important trend that will be used in the following since the generalized optical overlap has the inverse energy dependence in the high trapping energy limit.
The same non-activated energy dependence is asymptotically verified by the generalized overlap integral at high trapping energies. This is due to the fact that at high trapping energies, the clusters all concentrate in the region of highest electrical field, so that
More accurate expressions can be derived for specific optical trap types as
Ridge WG:
Slot WG:
At this point, all the expressions needed to evaluate the assay response have been derived, and we can return to
Eq. (9), the high-level equation describing the assay. At the beginning of this section we observed that the series in
Eq. (9) has to be effectively truncated in order to avoid a non-physical, diverging result. This truncation automatically occurs at finite assay stabilization times, since the concentration of large cluster species does not have time to stabilize. For large cluster species, the trap loading time is much larger than the assay stabilization time, so that the exponential describing the time dynamics can be linearized. For large
N, the species dependent response then scales as
where
τ
0 and Γ
0 were respectively defined in
Eq. (19) and (21). This equation describes a regime where cluster accumulation in the optical trap is diffusion and transport limited (i.e., how fast clusters arrive at the trap) rather than retention limited (i.e., how long clusters remain captured). The series in
Eq. (9) converges since τ
0(
N) is an increasing function (due to an increasing coefficient of friction for larger clusters) and
ϕ < 1. As already pointed out, for cluster species whose concentration has time to stabilize, the generalized overlap integral is much more favorable for the slot waveguide than for the ridge waveguide. However it should be noted that for transport limited species described by
Eq. (24), the two types of waveguides are roughly equivalent. It is also apparent from
Eq. (24) that once the
transport limited cluster accumulation regime has been reached, two smaller clusters will contribute more to the assay response than one large fused cluster, due to the fact that the diffusion constant is decreasing with cluster size. This is an important aspect since it explains why the assay sensitivities shown in section 6 drop again past a critical extent of aggregation.
6. Estimated assay sensitivity
This section is dedicated to the numerical investigation of assay sensitivity. We start with a few considerations relating to an experimental realization: In order for this scheme to work, the optical trap needs to be continuously turned on while the resonator captures NP clusters. Since the resonance of the resonator shifts during cluster accumulation, a feedback system will have to be implemented in an experimental setup for the wavelength of a tunable laser to continuously track the resonance of the ring while it shifts. An alternate method could rely on a broadband light source that covers the whole wavelength range over which the resonance travels.
We start the numerical analysis with a few orders of magnitude: For 20 nm radius NPs,
τ
0 as defined by
Eq. (19) is on the order of 10 ms for both waveguide types (precisely 6.2 ms for the ridge and 9.9 ms for the slot waveguide). Based on the slot waveguide data, it would respectively take 1, 10, 60 and 600 seconds for the concentration of clusters corresponding to trapping energies
NU
0/
kBT = 5.5, 8, 9.9 and 12.3 to stabilize.
Figure 5
shows the response of a slot waveguide ring to 10 pM NPs in the absence of target (zero extent of aggregation) and to 10 pM NPs with an extent of aggregation (
ϕ) of 10%. For unambiguous experimental results, it is beneficial for the resonant frequency shift induced by the addition of target molecules to be much larger than the resonant frequency shift induced by the addition of NPs alone. This ensures that the expected wavelength shift is a clear proof of the presence of the target, rather than the result of variations of possibly ill-controlled initial NP concentration. It is apparent in
Fig. 5 that a 10% extent of aggregation leads to a clearly distinguished response to the target (for a proper choice of trapping energy
U
0 versus assay stabilization time).
Fig. 5 Example of a transmission spectrum for a slot waveguide resonator immersed in pure water (black), after addition of 10 pM of gold NPs (blue) and after addition of a sufficient amount of target to obtain 10% extent of aggregation (red). The trapping energy for single NPs is U
0/kBT = 6.2 and the assay stabilization time is 60s.
In order to maximize assay sensitivity, the well depth of the optical trap,
U
0, has to be adjusted as a function of the assay stabilization time. If the well depth is excessively large at a given assay stabilization time, such that the accumulation of NP dimers inside the optical trap is diffusion limited (in the regime described by
Eq. (24), a situation can ensue where the response to NP monomers is larger than the response to NP dimers. This is due to the fact that the diffusion constant of monomers is larger than for dimers, and due to the fact that that the momomer concentration (~
ϕ) is larger than the dimer concentration (~
ϕ
2) in typical conditions. Since the monomers are present independently on whether or not the target is present, this leads to a situation where the baseline response to unbound NP swamps the response to target induced NP clustering. The cross-over point depends on the assay stabilization time: The larger the stabilization time, the larger the trapping energies that can be applied before dimer accumulation becomes diffusion limited. In order to illustrate this,
Fig. 6
shows the effective NP concentration at the resonator corresponding to individual cluster species as a function of well depth for two different assay stabilization times (1 min. and 10 min.).
Fig. 6 Effective NP concentration at the assay readout time as a function of cluster size N and single NP trap energy U
0/kBT for (a) 1 minute assay stabilization time and (b) 10 minutes assay stabilization time. The NP concentration is 10 pM and the extent of aggregation is assumed to be 10%.
For diffusion limited cluster formation the species dependent bulk cluster concentrations take the form [
34
N. J. Lynch, P. K. Kilpatrick, and R. G. Carbonell, “Aggregation of ligand-modified liposomes by specific interactions with proteins. I: Biotinylated liposomes and avidin,” Biotechnol. Bioeng.
50(2), 151–168 (1996). [CrossRef]
[PubMed]
]
while for reaction limited cluster formation they take the form [
34
N. J. Lynch, P. K. Kilpatrick, and R. G. Carbonell, “Aggregation of ligand-modified liposomes by specific interactions with proteins. I: Biotinylated liposomes and avidin,” Biotechnol. Bioeng.
50(2), 151–168 (1996). [CrossRef]
[PubMed]
]
These expressions can be used to evaluate the assay response at large extents of aggregation, with results shown in
Fig. 7
. For low extents of aggregation, it can be seen that the assay responsivity first increases with trapping strength. Past a critical trapping strength (
U
0 ~6
kBT for 60s assay stabilization time) the response to monomers becomes dominant and swamps the response induced by cluster formation (to the point where the responsivity can become negative in
Fig. 7(a) due to clustering induced monomer depletion). It can be seen that for
U
0 = 6
kBT the responsivity is not only maximized, but also features better linearity relative to the extent of aggregation. Independently on the assumed cluster formation reaction dynamics (diffusion or reaction limited) the assay response drops past a certain extent of aggregation, and collapses for
ϕ =1. Initially, increasing the extent of aggregation leads to an increase of the concentration of small cluster species that contribute to the assay response. However, past a certain point large clusters form to the extent of depleting the smaller clusters. Since large clusters have a smaller diffusion constant, and their accumulation is diffusion rather than retention limited, they contribute less to NP accumulation inside the optical trap than several smaller clusters, with larger diffusion constants but with an identical cumulative NP count.
Fig. 7 Resonant wavelength shift as a function of target concentration and trapping strength for a 60s assay stabilization time. Curves correspond to U
0/kBT ranging from 2 to 14 in increments of 2. The initial NP concentration is 10 pM and the assay stabilization time 60s. (a) assumes diffusion limited cluster formation and (b) assumes reaction limited cluster formation.
Wavelength shifts of 1 pm have been shown to be readily tracked in a temperature compensated system with similar Q-factors [
11
M. Iqbal, M. A. Gleeson, B. Spaugh, F. Tybor, W. G. Gunn, M. Hochberg, T. Baehr-Jones, R. C. Bailey, and L. C. Gunn, “Label-Free Biosensor Arrays Based on Silicon Ring Resonators and High-Speed Optical Scanning Instrumentation,” IEEE J. Sel. Top. Quantum Electron.
16(3), 654–661 (2010). [CrossRef]
], so that a minimum wavelength shift of 1 pm is taken as baseline to derive sensitivity floors.
Figure 8
shows the assay sensitivity as afunction of stabilization time, assuming
ϕ =
cT
/
cNP
and assuming
ϕ = 0.1 (i.e., the initial NP concentration
c
NP is set to one order of magnitude above the target sensitivity). The well depth of the optical trap,
U
0/
kBT, is set to 0.765×log(
t)+3.05, where
t is the assay stabilization time in seconds. The sensitivity floors for specific assay chemistries with lower extents of aggregation can be derived by using
ϕ/[
cT
/
cNP
] as a corrective factor. The concentrations in
Fig. 8 can also be shown to correspond to trapping of just a few NP dimers (~5) by a typical ring resonator (
R = 50 μm) so that they also correspond to the ultimate sensitivity limit irrespectively of the detectable wavelength shift.
Fig. 8 Assay sensitivity as a function of stabilization time. The trapping energy is implicitly varied as described in the text, cNP
is assumed to be set at 10× the sensitivity floor, and ϕ is assumed to be equal to c
T
/c
NP
.
7. Model for cluster capture dynamics
In 1940 Kramers derived an analytical formula for the particle escape time out of a 1D energy well in the overdamped regime [
30
H. A. Kramers, “Brownian motion in a field of force and the diffusion model of chemical reactions,” Physica
7(4), 284–304 (1940). [CrossRef]
]. Kramers’ derivation assumes an energy well with a quadratic (harmonic) spatial energy dependence around the minimum (
U=-
U
0 +
κ
1
x
2/2 around the minimum) as well as a quadratic spatial energy dependence around the saddle point (
U = -
κ
2
x
2/2 around the saddle point). This applies to most 1D physical problems with well-defined, non-zero second-order derivatives at the well minimum and at the saddle point. However, the problem of a NP in an optical trap differs in several aspects. Most importantly, there is no saddle point, and the trap boundary is ill defined. In the case of a ridge waveguide, the bottom of the energy well is also non-harmonic and the second order derivative is ill defined. Nonetheless, the methodology used by Kramers can be applied to an optical trap, and has been previously applied to the case of a 1D evanescent field [
26
T. J. Davis, “Brownian diffusion of nano-particles in optical traps,” Opt. Express
15(5), 2702–2712 (2007). [CrossRef]
[PubMed]
]. However, there is an important difficulty associated to the absence of a saddle point in optical evanescent traps: In 1D and 2D, the escape time is a function of the location of the trap boundary, and diverges when increasing the distance from the center of the trap to the boundary. In order to alleviate this problem, it has been suggested in previous approaches [
26
T. J. Davis, “Brownian diffusion of nano-particles in optical traps,” Opt. Express
15(5), 2702–2712 (2007). [CrossRef]
[PubMed]
] to define the trap boundary as the location where the trapping energy becomes commensurate with the thermal kinetic energy. Here, we show that the ambiguity on particle escape time is an intrinsic property of systems with lower dimensionality and that the Kramers escape time can be unambiguously defined for an evanescent optical trap in 3D. By taking into account transient effects in lower dimensionality systems, physically meaningful trap loading times can be derived. We also derive a 2D transient model corresponding to an isolated silicon wire, or to a ring resonator in the limit of large radii and low trapping energies.
The section starts with a validation of the basic assumptions required in order for Kramers’ derivation of a particle escape time to be valid (overdamped limit). We then show that setting a boundary at an energy level commensurate to
kBT is a valid methodology to calculate the NP thermalization time, i.e., the time constant defined in section 2 as the characteristic time for a trapped NP cluster to reach, at steady state, a location in the solution where its local concentration is in thermodynamic equilibrium with other cluster species, and the relevant time constant for comparison with the cluster dissociation time. We proceed by deriving analytic expressions of the cluster thermalization times for the cases of both ridge and slot waveguides. This is followed by a discussion of trap loading times in lower dimensional systems and 3D systems, showing that the classic Kramers derivation yields a well defined trap loading time in 3D, but has to be refined with transient effects in lower dimensional optical traps to yield physically meaningful results. In 3D, the trap loading time is derived in two steps, as a locally 2D problem first analytically solved in the vicinity of the waveguide and then “stitched” to a 3D environment. In lower dimensional systems, a different criterion is introduced to place the trap boundary applied to the classic Kramers escape time derivation: The trap boundary is placed at a location such that the total particle concentration depletion inside the boundary during transient loading is commensurate with the total number of particles eventually captured by the trap. These methodologies result in two formulas for the trap loading time, one valid for 2D systems (
Eq. 44), e.g. straight waveguides or ring resonators in the limit of low trapping energy or large radius (in which case the problem reverts to a locally 2D problem), and another valid for small ring resonators at sufficiently high trapping energies (
Eq. (45). It should be further noted that these formulas predict trap loading times in the absence of fluid displacements, including convection currents, and are worst case estimates in the presence of fluid displacements.
Kramers’ derivation of a particle escape time is accurate in the overdamped regime. The NP dynamics are described by the Langevin equation corresponding to a harmonic oscillator with damping factor
where
m is the mass of the cluster and
ω
0 is the oscillation frequency. For the trap generated by the slot waveguide,
ω
0 is a well-defined quantity (
, where
m
0 is the mass of a single NP). In the case of the ridge waveguide, the highest oscillation frequency corresponds to maximum amplitude oscillations and is given by
. For
NU
0 = 12
kBT and
N = 2, this results respectively in damping factors of 60 and 85 for the ridge and the slot waveguides. In other words, cluster dynamics are well within the overdamped regime for both types of traps, and particle escape times can be derived with Kramers’ methodology.
We start the section by deriving the particle thermalization time in a 1D model [
26
T. J. Davis, “Brownian diffusion of nano-particles in optical traps,” Opt. Express
15(5), 2702–2712 (2007). [CrossRef]
[PubMed]
]. A 1D model is appropriate to calculate the cluster thermalization time when setting the boundary at an energy level on the order of
kBT. Indeed, this results in a boundary relatively close to the surface of the waveguide, so that escaping particles can be approximated as following a 1D flow. If the cluster escape time computed by setting the boundary at
εkBT (where
ε is a number smaller than 1) is much smaller than the cluster dissociation time, it can be concluded that the cluster concentration at the bottom of the well is exp([
NU
0-
εkBT]/
kBT)
cb
, where
cb
is the steady state cluster concentration at the boundary. For non-dissociating clusters,
cb
=
c
0exp(
ε). However, in any case
cb
>
c
0, resulting in a cluster concentration at the well minimum
c(
N) > exp(-
ε)exp(
NU
0/
kBT)
c
0. If 1-exp(-
ε) is small, the validity of
Eq. (3) is adequately verified. For example, if
ε = 1/4,
Eq. (3) is verified with an accuracy better than 22%. The boundary corresponding to the energy level
U =
kBT/4 is located at log(4
NU
0/
kBT)
δ = 390 nm from the waveguide surface for
NU
0 = 12
kBT. The distance between the boundary and waveguide remains on the order of the trap width, so that a 1D model remains relatively accurate (it actually underestimates the particle capture rate
η and provides an upper bound for the thermalization time).
The general expression for a particle escape time from a potential well is the ratio of the probability to find the particle inside the trap (the integral of the probability density function) to the particle flux escaping the trap. The probability density function is dominated by the bottom of the well, where it is close to the equilibrium concentration and is given by the Boltzmann statistics
Eq. (3). The escape flux on the other hand is a function of the particle concentration away from the well minimum. The latter is off-equilibrium due to the escaping flux. The particle flux can be written as the sum of the diffusion flux and drift flux:
In a 1D problem,
J can be approximated as a constant away from the bottom of the trap, since the latter is the primary particle source.
J can then be integrated as
where
ya
corresponds to the position of the bottom of the trap, and
yb
to the trap boundary. It should also be noted that
U(y) denotes an energy, i.e., with the sign convention of a physical energy, while
NU
0 is the well depth for a cluster of
N NPs, i.e.,
NU
0 = -
U(
ya
). As mentioned previously,
c(
ya
) is assumed to follow equilibrium statistics so that
c(
ya
) =
c
0exp(-
U(
ya
)/
kBT), but
c(
yb
) is set to 0, i.e. it is assumed that any particle reaching the trap boundary is absorbed by the boundary and never returns into the trap. The particle escape time is then expressed as
where the first integral corresponds to the particle flux and the second integral to the particle probability distribution (integrated starting from the slot or waveguide surface assumed to be at
y = 0). It has been suggested [
26
T. J. Davis, “Brownian diffusion of nano-particles in optical traps,” Opt. Express
15(5), 2702–2712 (2007). [CrossRef]
[PubMed]
] to use as a boundary the point where the trapping energy of the particle equals a given small fraction of the mean thermal energy, since the trapping force is then too small to play a significant role and the particle dynamics revert to those of a free, unbound particle. Here we set it to
kBT/4, half the mean 1D thermal energy. The first integral in
Eq. (30) corresponds to the 1D model for the inverse particle flux. It is dominated by the high-energy region at the outer edges of the trap, and is identical for both waveguides
Equation (31) is the equivalent of
Eq. (13), but the length scale
Iw has been replaced by
δ. Since
h
trap has already been derived in section 5 the escape times required for particles to reach the boundary at the energy level
kBT/4 can be directly derived as
This is the expression for the cluster thermalization time, and for the waveguide geometries explored here it is a factor
δ/
Iw = 4 (slot waveguide) to 8 (ridge waveguide) smaller than the trap loading time. For
NU
0 = 12
kBT this corresponds to a cluster thermalization time of respectively 10s and 170s for the ridge and for the slot waveguides. These numbers can be compared with the cluster dissociation times described in section 3. The particle escape time for ridge waveguides is identical to the escape time derived in [
26
T. J. Davis, “Brownian diffusion of nano-particles in optical traps,” Opt. Express
15(5), 2702–2712 (2007). [CrossRef]
[PubMed]
]. Here, the non-activated dependence on the well depth is also explicitly derived for both ridge and slot waveguides.
Next, we inspect the dependence of the calculated escape time on the boundary location and show that trends in lower dimensional systems fundamentally differ from 3D. The exact expression for the 1D escape time of the ridge waveguide is obtained by a change of variables and given by
where
z is the local well depth of the optical trap in units of
kBT. It is immediately apparent that both integrals diverge with vanishing
z
min =
NU
0exp(-
yb
/
δ)/
kBT (i.e. when the trap boundary is pushed out). The assumption that a particle never returns to the trap after crossing a given boundary is the crux of why the escape time is ill defined for a 1D (or 2D) evanescent trap: In the classic Kramers problem (
Fig. 3), the probability of a particle returning to the trap once it is past the saddle point and its potential energy has sufficiently decayed is vanishing, so that
c(
yb
)=0 is a natural assumption to make at the boundary. Similarly, in a 3D system the probability for a particle to ever return to the trap is smaller than 1 (for an unbounded solution). On the other hand, in a lower dimensional system without a saddle point, a particle returns to the core of the trap an increasing number of times before reaching a boundary that is moved away from the trap center. The significance of this divergence can be seen by deriving the asymptotic trend of the escape time, for example for the 1D evanescent field (ridge waveguide)
where
y
0 is the location where the trap energy reaches ~0.6
kBT. An additional refinement has been made in the particle concentration assumed in Eq. (34), made necessary because regions far from the optical trap are included. In
Eq. (30) and (33) the particle concentration was integrated assuming steady state statistics throughout the optical trap. This is an appropriate approximation when the integral is dominated by the particle concentration around the energy minimum as for a “natural” boundary close to the trap. Here, however, the asymptotic limit for remote trap boundaries is taken, so that the particle concentration far from the bottom of the trap contributes to a substantial amount. The latter is off-equilibrium. Since particle current in this region is dominated by diffusion, in a 1D system the particle concentration is described by
c
0(
yb
-
y)/(
yb
-
ya
) for large
y.
The first term in the right hand side of Eq. (34-b) can be immediately recognized as the time taken for a particle to freely diffuse over a distance
yb
in a 1D system. The second term is due to the fact that this is a 1D problem. In a random walk in a 1D lattice, a particle will on average return
X times to the origin before reaching a boundary located
X lattice steps away from the origin. The second term of Eq. (34-b) is the continuous space equivalent. Equivalently, it is also explained by the fact that in lower dimensional systems the flux generated by a particle sink does not converge for an increasingly remote boundary: In 1D, it scales as
Dc
0/(
yb
-
y
0) assuming a particle sink at
y
0.
Figure 9
compares the boundary dependent escape time to the asymptotic model.
Fig. 9 Particle escape time as a function of boundary position for NU
0/kBT = 12, as determined with the 1D model. The fit given by Eq. (35) is shown by the dashed line.
In a 3D model the concentration depletion is much less severe, as already discussed in section 5, due to the fact that under steady state conditions, given a fixed particle flux, the concentration in the solution far from the optical trap converges to a finite value. In a 3D environment, the probability for a particle to ever return to the optical trap via random walk is also less than unity, so that an escape time is a well-defined quantity even in the absence of a saddle point. If one defines r
0 such that the flux escaping from the 3D trap is 2πDr
0
c
0 (as in section 5), given the concentration at the center of the trap to be exp(NU
0/kBT)c
0, one finds that the particle concentration far from the trap scales as c
0(r
0/r-r
0/rb
) in order to maintain the constant flux and to satisfy the zero concentration boundary condition at rb
. This leads to (in spherical coordinates)
Even though the numerical value for r0 has not yet been determined, it can be immediately seen that the 3D escape time is the sum of the 3D diffusion term rb
2/6D and a finite offset corresponding to the trap loading time. Since h
trap has already been derived, only the escaping particle flux remains to be calculated in order to fully determine the trap loading time.
In section 5, the 3D trap loading time was derived assuming the trap to be equivalent to a flat, zero concentration boundary condition with the same dimensions than the waveguide surface (ridge waveguide) or than the slot footprint. A more accurate trap loading time taking into account the effect of the optical field on particle transport can be calculated by adapting
Eq. (28) and
(29) to a higher dimensional system. Instead of a fully analytical treatment of the complete 3D problem, it is broken down in a locally 2D problem (semi-analytically solved) with a numerically determined corrective factor taking into account the 3D environment. This is motivated by the fact that the local 2D problem around the waveguide is relatively straightforward to solve semi-analytically and that the global 3D particle transport farther from the waveguide, while more complex, is purely diffusion driven and independent of the waveguide specifics, and thus simply taken into account by numbers extracted from a 3D simulation of the diffusion equation.
Figure 10
shows the constant concentration contours around the particle sink (in cylindrical coordinates as in
Fig. 4(a)). Even though the diffusion equation is simulated with cylindrical coordinates, the concentration behaves as in a 2D system in the immediate vicinity of the trap (as seen by concentration contours centered on the waveguide), so that the escape time up to the dashed boundary in
Fig. 10(b) can be modeled as a 2D problem. On the other hand, past this boundary the evanescent optical field has sufficiently decayed to be ignored in the treatment of the surrounding 3D problem. The distance to this
intermediate boundary is denoted as
ri
and the concentration at the boundary as
ci
. The 3D environment is taken into account by numerically estimating
ci
as a function of the particle flux crossing the boundary. Since the estimated 2D trap loading time is slowed down by a factor
ci
/
c
0 there is a unique
ci
and particle flux pair that verifies self-consistency.
Fig. 10 Constant concentration contours generated by a 400 nm wide particle trap (trap boundary corresponding to c = 0) with R = 20 μm, simulated in cylindrical coordinates. The dashed line shows the ri
= 1 μm boundary of the 2D analytical calculation.
The relationship between
ci
and the particle flux crossing the intermediate boundary is simulated by placing a boundary condition at the dashed hemi-circle corresponding to
(so that (2
πR)
η
2D
wci
is the total flux) and by fitting
ci
as a function of
ηw. The subscript 2D in
η
2D serves as a reminder that this is the particle capture rate of the 2D problem referenced to the concentration at the intermediate interface
ci
rather than the bulk concentration
c
0. The interface between the two domains can be modeled this way since both the concentration and the particle flux are constant across the hemi-circular boundary (e.g. in
Fig. 10, the boundary coincides with a fixed particle concentration contour). For
R = 20 μm and a boundary located at
ri
= 1 μm from the center of the trap (
Fig. 10(b)) the concentration at the boundary is fitted as
Inside the intermediate boundary, the particle flux is analytically solved as a 2D problem. The conversion from
Eq. (28) to
Eq. (29) has to be adjusted in order to account for “flaring out” particle flow (
Fig. 12
). Indeed, the initial derivation relied on the particle current,
J, to be independent on
y in a 1D system. This is not true in a system with a higher dimensionality. In a 2D system, as locally around the waveguide, the integral of the particle flux over a contour encircling the waveguide is constant, but the local particle flux
J(
x,
y) itself is varying.
Equation (29) then needs to be recast into
where
c(
ya
) = exp(-
U(
ya
)/
kBT)
ci
(in equilibrium with the bulk concentration of the 2D problem),
c(
yi
) = 0 (absorbing boundary condition) and where
l(
y) is an effective length defined by
Fig. 12 (a) Constant particle flux contours and (b)
l(
y) as defined in
Eq. (38) as simulated from a flat 400 nm wide particle sink in cylindrical coordinates.
Close to the waveguide surface, , since the system can be approximated as a 1D system. Further from the waveguide, the problem acquires a cylindrical geometry. The particle flux fans out, and l(y)= πr where r is the distance to the waveguide.
Figure 11
shows the equi-energy levels for the trapping field. It can be seen that they also assume a cylindrical geometry away from the trap.
Figure 12 shows the particle flux as extracted from the simulation of the diffusion equation with a flat particle sink (section 5), i.e., ignoring the effect of the optical trapping field.
l(
y) has the expected trend, with
l(0) ~1.6
w >
w due to the edge effect on the sides of the particle sink (
Fig. 12(a)). In the following calculations,
Eq. (37) is evaluated assuming the real values for the optical trapping field, but assuming that
l(
y) can be approximated by the values corresponding to the flat particle sink.
l(
y) is fitted as
Fig. 11 Constant trapping energy contours for NU
0/kBT = 12. The contours correspond to kBT/4, kBT/2, kBT, 2kBT, 4kBT and 8kBT.
The first term in the right hand side of
Eq. (40) corresponds to the 2D problem with a boundary set at
ri
= 1 μm from the waveguide and 1/
l(
y) taking into account the 2D nature of the problem. The second term, 1.56, is a numerical correction taking into account the 3D environment for a 20 μm radius ring.
As can be seen in
Fig. 13
, the predicted values for
I
lim in the absence of an optical trapping field (
U
0 = 0) coincide with the values numerically estimated in section 5. For a non-zero
U
0,
I
lim is smaller since the optical trapping field effectively reduces the distance over which clusters have to freely diffuse before being captured by the trap. Another effect is that
I
lim is much less sensitive to
w once the optical trap has been turned on.
w impacts
l(
y) for small
y (
y <
w). However for
y smaller than
δ, exp(
U(
y)/
kBT) is small so that the impact on the integral in
Eq. (40) is reduced. I.e.,
I
lim is much less sensitive to
w in the region
w <
δ. Effectively, free particle diffusion only occurs outside the evanescent field of the waveguide where the particle flow already assumes a hemi-circular form when
w <
δ. For the typical optical trap dimensions explored in this paper, i.e.,
w on the order of 200 to 400 nm and R = 20 μm, the corrected value for
I
lim can be seen to be ~1.9.
Fig. 13 Coefficient
I
lim as calculated with
Eq. (40) in the absence of an optical trapping field (
U
0 = 0,
w=400 nm, dashed line) and with
U
0 = 12
kBT (continuous line,
δ=100 nm).
A general formula for I
lim will be derived later in this section as a function of the resonator dimensions, however an intuitive understanding of this formula requires a transient model for purely 2D systems (straight waveguides) to be derived first.
As already explained in section 5,
I
lim determines the trap loading time (it can be plugged in
Eq. (14) if
NU
0 is large enough such that the number of particles accumulated by the optical trap exceeds the number of particles that have to be captured in order for the particle flux to stabilize to its 3D limit. If the number of particles accumulated by the optical trap is smaller, the particle flux never reaches this stabilized limit during trap loading and remains larger, since the surrounding solution depletion is much less than predicted by the stabilized 3D model. For the geometry explored in section 5 (
R=20 μm) this occurs for
NU
0 > 12.5
kBT. For smaller trapping energies,
I
lim corresponds to an upper bound and one can estimate the trap loading time more precisely by making the right choice on the location of the boundary (
yb
): It should be placed such that the cumulative particle depletion inside the boundary corresponds to the number of particles accumulated by the trap. This boundary condition results in a particle flux corresponding to the lowest value the particle flux drops to during transient trap loading, and is a good approximation when comparing to fully numerical results.
For very small NU
0 the problem reverts to a 1D problem. The trap loading time then scales as τ~2δ2/D⋅exp(2NU
0/kBT)/(NU
0/kBT)2χ provided yb
~2δ ⋅exp(NU
0/kBT)/(NU
0/kBT)χ < w. Here, this only holds for very small trapping energies (NU
0/kBT < 2) due to the fact that w and δ are of the same order of magnitude.
In the intermediate energy range, the trap loading time corresponds to a transient 2D problem. As previously, an effective trap radius
r
0 is introduced to parameterize the particle flux, even though the definition of
r
0 is different here due to the fact that we are considering a 2D problem. If the domain boundary is set by a hemi-circle at a distance
rb
from the waveguide, at steady state the particle concentration scales as log(
r/
r
0)
c
0/log(
rb
/
r
0) given a steady state particle flux
πDc
0/log(
r
b/
r
0). From
Eq. (37),
r
0 can be shown to verify
In other words, r
0 is the radius of an ideal hemi-cylindrical particle sink such that the induced particle flux is equal to the particle flux generated by the actual optical trap. For example, for the ridge waveguide r
0 = 270 nm for NU
0 = 6 kBT and 330 nm for NU
0 = 12 kBT. The total particle depletion up to the boundary located at the distance rb
can be calculated as c
0
πrb
2/4/log(rb
/r
0). By equating it to the total number of particles eventually accumulated by the trap, c
0
δw⋅ exp(NU
0/kBT)/(NU
0/kBT)χ, one finds
For the ridge waveguide (
χ = 1) this implies log(
r
b/
r
0) ~
NU
0/2
kBT at large trapping energies (this is also true at small trapping energies if 2
δw/
πr
0
2~1). This results in a particle flux
wηc
0 ~2
πDc
0/(
NU
0/
kBT) and
for
NU
0 < 12
kBT, after which
I converges to its 3D limit,
I
lim. The expression for the transient 2D loading time is verified against time domain simulations of the diffusion equation (
Fig. 14(a)
).
Fig. 14 (a) 2D transient trap loading time simulated in the time domain for
δ = 100 nm and (blue crosses)
w = 125 nm,
r
0 = 150 nm (black dots)
w = 630 nm,
r
0 = 300 nm (red stars)
w = 3150 nm,
r
0 = 1100 nm. The continuous lines show the prediction given by
Eq. (43). It can be seen that it is accurate for large trapping energies in all three cases. For the case of the black dots, 2
δw/
πr
0
2 = 0.45 is close to 1 and it can be seen that
Eq. (43) is also accurate at low trapping energies. (b) shows the simulated 3D correction term (1.56 in
Eq. (40) as a function of
R when
ri
is held constant at 1 μm (black dots). The dashed line shows the prediction given by
Eq. (45) (with
r
0 set to 1 μm). The continuous curve shows the corresponding prediction for the loading time of the 400 nm wide ridge waveguide (with
r
0 set to 270 nm).
The transition between the 2D and the 3D regime is a function of the resonator radius R. In particular, the transition between the 2D and 3D regimes occurs when the depletion region estimated at the distance rb
from the trap center merges with the depletion region from the rest of the ring, i.e., when rb
becomes commensurate with the ring resonator radius R and the solution depletion region converts from a half-toroid to a half sphere. Thus it comes as no surprise that at NU
0 = 12 kBT, rb
~130 μm is on the order of 2R (40 μm).
It was numerically verified that the 3D corrective term (1.56 in
Eq. (40) for a 20 μm radius ring resonator) scales as log(6.7
R/ri
)/
π . The logarithmic dependence on
ri
can be a priori justified by the dependence of a 2D flux on the sink radius,
πDc
0/log(
rb
/
ri
).
Equation (40) can then be recast as
Equation (45) is verified against steady state simulations of the diffusion equation with varying ring resonator radius
R (
Fig. 14(b)).