Optimizing the localization precision in coherent scattering microscopy using structured light

2.1 Fisher information

To answer these questions, we will use the concept of the Fisher information (FI) [3], [7]. It accounts for the amount of information that the measurement data carry about an unknown parameter. Let us first consider the localization estimation using experimental iscat data. Our goal is to estimate the particle position R ₀ from the iscat images I ( x , y ; R 0 ) = I ( x , y ; θ ) that are noisy (Poisson noise) because of the limited number of photons collected [7]. We have followed the usual convention adopted in Fisher information theory and have denoted the estimation parameters with θ . In the following, we will only consider localization estimation, this is θ = R ₀, but our approach could be easily modified to deal with different estimation parameters as needed, e.g., in mass photometry or the discrimination of closely spaced scatterers [29].

Suppose that we have two different incoming fields ϵ _inc at hand. Through the Fisher information, we can determine which field provides greater sensitivity to the parameter of interest and should thus be used in noisy experiments. In the language of probability theory, p I ( x , y ; θ ) is a (Poisson) probability mass function whose expectation value gives I x , y ; θ . We next introduce the score

which accounts for the (normalized) change of p I when varying the estimation parameter θ _i , see also Figure 1. The Fisher information is defined as the variance of the score,

A large Fisher information indicates that the measurement data change significantly upon variation of θ _i , and thus a high localization precision can be obtained in experiment. The expectation value E in Eq. (5) has the form [7], [13]

where the sum extends over the pixels x, y of the image detector and Δx, Δy denote the pixel size. We have introduced a normalization factor N that is conveniently set to the intensity of the scattered fields in the backfocal plane [13] (proportional to the number of detected photons). In the study of structured light fields, we proceed somewhat differently and use instead the intensity of the incoming fields (proportional to the number of incoming photons). From the Fisher information, one can define the Cramér–Rao bound on the standard deviation of any unbiased estimator for the parameter θ _i per incoming photon as

The equality sign should be understood as an assignment of notation, where σ inc iSCAT sets a lower bound for the standard deviation achievable from any estimator saturating the Cramér–Rao bound. The corresponding bound per detected scattered photon is related through

where N_inc and N_det denote the numbers of incoming and detected photons, respectively.

2.2 Quantum Fisher information

We can extend the concept of Fisher information to the quantum regime, where the quantum Fisher information (QFI) quantifies the sensitivity of a quantum state to changes in the parameter of interest, independent of any specific measurement scheme. In other words, it doesn’t rely on iscat or any other variant of coherence microscopy but addresses the question what could be obtained under the best possible measurement conditions. As has been shown in Ref. [24], for coherent scattering measurements (as studied in this work), the quantum Fisher information can be obtained from the scattered fields in the backfocal plane through

where ρ_NA is the cutoff radius for a given numerical aperture (NA) of the lens. Below we will show that these optimized fields can be obtained in a computational approach through a matrix diagonalization, which makes the approach extremely simple and powerful. In the same way as for the Fisher information discussed above, we can introduce a quantum Cramér–Rao bound

which sets a lower bound on the standard deviation of the best possible estimator for the parameter θ _i and per incoming photon.

2.3 Normalization of Fisher information

When normalizing the Fisher information, we can do so in many ways. We can normalize per incoming photon, per scattered photon, per absorbed photon, or per detected photon. In order to choose the right normalization, we have to know what limits the accuracy in our measurements experimentally. Here are some examples:

1. If the number of photons per time available from the light source is limiting the experiments, then one should optimize per incoming photon. This is rarely a problem these days.

2. In fluorescence microscopy, it is known that the excitation, relaxation, and photobleaching of the fluorophores lead to reactive oxidative species that can harm live cells and organisms [30]. It is, therefore, often appropriate to maximize the Fisher information per excitation–relaxation cycle, i.e., per emitted fluorescence photon.

3. In coherent microscopy, phototoxicity is often less of a problem, and many experiments are limited by the low photon detection rate, i.e., the number of photons that can be recorded per pixel per frame. In this case, one should normalize per detected photon, which we did in a previous paper [13].

4. However, even if an imaging modality is based on elastic scattering, the illumination light can still cause phototoxicity, e.g., due to molecules that absorb light. A common example is the imaging of live-cells with blue or UV light, but even infrared light can inactivate bacteria at moderate intensities [31], which are commonly applied in coherent scattering microscopy. In this case, it is then again appropriate to normalize per incoming photon. With coherent scattering imaging being done at faster and faster framerates, requiring higher and higher excitation intensities [4], this regime will become more prominent in the future, which is why we chose this normalization for our manuscript.

2.4 Computational approach

In our computational approach, we adopt the coherent scattering microscopy add-on [28] to the generic Maxwell solver toolbox nanobem [32], [33]. Contrary to our previous work, we are only interested in small and weak scatterers that can be described well in the dipole approximation [34] without solving the full Maxwell equations for the dielectric particle. All our simulations can be broken down into the following steps (see [28] for a more detailed discussion).

1. We start by specifying the incoming fields on a polar grid, ϵ _inc(ρ _m , φ _m ), which become focused toward the interface using the Richards–Wolf approach [35], [36] in order to obtain the fields E _inc( R ₀) at the scatterer position R _o .

2. The incoming fields induce an oscillating dipole moment p = α E of the scatterer, where α is the polarizability of a sphere [34], which produces outgoing waves e i kr r F sca ( r ̂ ) in the optical far-field [28], Eq. (6)] (see also far-field pattern in Figure 1). Note that in coherent scattering microscopy and for spherical scatterers, the dipole moment always points in the same direction as the excitation field.

3. The scattered fields are collected by a lens and converted to fields ϵ _sca(ρ _m , φ _m ) propagating in the backfocal plane along −z. Again, the fields are given on a polar grid.

4. Using the Richards–Wolf approach for imaging, we compute the scattered and reference fields in the image plane [28]. The intensity of the superposition of these fields then provides us with the iscat image I ( x , y ; R 0 ) , see Eq. (3).

We are now ready to compute both types of Fisher information. We first approximate the derivative of the scattered fields in the backfocal plane using finite differences.

with a similar expression for ∂ θ i I ( x , y ; θ ) . Here η is a small number and e ̂ θ i is the unit vector along the direction θ _i . This is depicted in Figure 1(b) for θ _i = X₀, where the scatterer is slightly displaced along x and the scattered far-fields are computed for the displaced scatterer. For the integral over the backfocal plane, see Eq. (9), we use a polar grid to get for a generic function F ( ρ , φ ) the expression

where w _m are the quadrature weights adopted in the numerical approach, such as Gauss–Legendre weights in this work. We next use that in linear response the scattered fields can be related to the incoming fields through a transition matrix T ( θ ),

where m labels the polar grid points and we have introduced the short-hand notation ϵ _m = ϵ (ρ _m , φ _m ). In computing the quantum Fisher information, we set the normalization factor to the intensity of the incoming fields, N = ∑ _m w _m | ϵ _inc,m|², and obtain from Eq. (9)

with the diagonal matrix w _mm′ = w _m δ_mm′. Equation (13) is a quadratic form, and the best possible fields are obtained from the solutions of the generalized eigenvalue problem (see [24] for a related approach)

The eigenvectors associated with the largest eigenvalue Λ = QFI(θ _i ) are the optimized fields ϵ inc ( o p t ) with the best possible localization precision. Instead of optimizing for a single localization parameter θ _i , in this work we replace the matrix F θ i with

where R ₀ = (X₀, Y₀, Z₀) is the particle position. Using this expression, in the diagonalization, we then optimize the localization precision for all three spatial directions rather than a single one, where γ is a scaling parameter that weights the relative importance for localization along X₀, Y₀ with respect to Z₀. Quite generally, one could also consider off-diagonal contributions to the Fisher matrix, but for the setups considered in this work, they turned out to be of minor importance and were thus neglected.

For iscat experiments, the optimization of the FI(θ _i ) in Eq. (5) cannot be expressed as a quadratic form in the incoming fields and cannot be converted to an eigenvalue problem or something related. We thus have to proceed differently. In our computational approach, we start from some initial guess and maximize the Fisher information in an iterative fashion using the fminunc function of matlab. More details about our approach will be given in Section 3.3.

3.1 A first glimpse on optimization

In order to get a better feeling of what we are aiming for, in Figure 2, we plot the field intensities and phase distributions for three incoming fields, namely a plane wave impinging from below on the glass–water interface (left, a–c), a focused Gauss beam (middle, a’–c’), and an optimized beam (right, a”–c”) that is obtained from the maximization of the quantum Fisher information in Eq. (14). For the plane wave polarized along x, we observe below the interface the interference between the primary incoming wave propagating in the +z direction and the secondary wave reflected at the interface propagating in the −z direction. In the middle panels, we show the results for a Gaussian incoming field with circular polarization

which becomes focused through a lens with a numerical aperture of NA = 1.4 (see Figure 1). E₀ is chosen such that the field strengths of the plane wave and the focused Gauss beam are the same at the scatterer position R ₀. When comparing the quantum Fisher information values of Eq. (9), reported at the bottom of Figure 2, we observe that the focused Gauss beam has a better localization precision. In the right panels, we show that the localization precision can be further boosted through an optimization of the quantum Fisher information. This shows that, with constant electric field at the scatterer position, a narrow illumination beam yields the highest quantum Fisher information regarding localization in the scattered light. This result has to be met with caution. First, with a tightly focused illumination beam, it will be difficult to tune the scattered and the reference fields separately (e.g., via defocusing for phase, or via the use of attenuators for amplitude tuning), potentially leading to the Fisher information being significantly smaller than the quantum Fisher information. Second, if a focused beam is scanned in order to image an extended field of view, the average quantum Fisher information will probably be lower for the focused beams.

Figure 2:

Intensity and phase maps for coherent scattering microscopy using a plane wave (left), a focused Gaussian beam (middle), and an optimized beam (right). The focusing lens has a numerical aperture of NA = 1.4. We report the intensity of the incoming beam in the (a) xy (z = 0) and (b) xz-plane. The glass–water interface is indicated by a dashed line, and we also show the scatterer with a diameter of 20 nm located on top of the interface. Panel (c) reports the phase map for the incoming field E_inc,x. Below the figures, we report the quantum Fisher information in arbitrary units. In the simulations, we use the same incoming intensities for the Gauss and optimized beam, the field strength of the plane wave is chosen such that it equals that of the Gauss beam at the position of the scatterer. The intensity of the optimized beam is scaled by a factor of 0.6.

3.2 Optimization of quantum Fisher information

In this section, we provide a detailed discussion of the optimized structured-light beams and the resulting quantum Fisher information values. We start from Eqs. (14a) and (15) with γ = 0.2 and comment on the importance of the scaling parameter γ in Appendix A. Figure 3 shows the optimized incoming modes in (a, a*) the backfocal plane and (b, b*) the plane where the scatterer is located. The left panels (a–c) report results for the scatterer on top of the interface, and the right ones (a*–c*) for a particle–interface distance of 1 μm. In the first rows of panels (a, a*) and (b, b*), we show results for a numerical aperture of NA = 1.2, and in the second rows for NA = 1.4.

Figure 3:

Modes with optimized quantum Fisher information in (a, a*) backfocal plane and (b, b*) plane where the scatterer is located, and for scatterer (a, b) on top of glass–water interface and (a*, b*) located 1 μm away from the interface. The colored quarter-disks in (a, a*) show the intensities of the incoming beams, and in the grayscale quarter-disks, we project the field-distribution onto radial, azimuthal, and circularly polarized (σ₊) unit vectors. The optimized modes exhibit (M1) radial, (M2) circular, and (M3) azimuthal polarization. For comparison, we also show a (G) Gaussian mode with circular polarization. The modes in the first and second rows of panels (a, b) show results for numerical apertures of 1.2 and 1.4, respectively. (c, c*) Eigenvalues of Eq. (14), which are proportional to the quantum Fisher information, as a function of NA for different modes. For the particle on top of the substrate and for large NA values, we scale the eigenvalues by a factor of 0.125 for better visibility.

Although in our optimization approach we do not enforce any symmetry constraints, the problem under study exhibits cylinder symmetry, mainly due to our consideration of small particles (note that things would strongly differ for larger particles [39]), and thus also the optimized solutions exhibit such symmetry. For this reason, we plot in Figure 3(a) the field intensities and polarization properties only in a quarter of a disk. The colored parts (upper-left quarter-disks) report the intensities of the optimized incoming fields ϵ _inc(ρ, ϕ) in the backfocal plane. In the remaining quarter-disks, we project the fields onto unit polarization vectors u _λ through

where the degree of polarization P is bound to values between zero and one. In the figure, we use (in clockwise direction) polarization vectors in the radial and azimuthal directions, as well u ± = ( x ̂ ± i y ̂ ) / 2 for the helicity basis σ_±. Black areas in the quarter-disks then indicate that the modes have an almost pure degree of the corresponding polarization state. Although not enforced by our optimization approach, all modes associated with the highest eigenvalues can be characterized in terms of (M1) radial polarization, (M2) circular polarization, and (M3) azimuthal polarization. All circular modes are twofold degenerate (only one of them is shown), and one could alternatively perform a superposition to obtain modes with horizontal and vertical linear polarization. In the figure, we additionally show results for a (G) Gaussian incoming beam with circular polarization and will comment on its properties later. Note that the radial and azimuthal unit vectors form a basis, but the vectors are not orthogonal to σ_±. For this reason, there is an approximately fifty percent σ₊ contribution to the radial (M1) and azimuthal (M3) modes. One observes that modes M1 and M3 have zero intensity at ρ = 0, whereas the circular modes M2 and G exhibit a rather weak intensity variation along ρ.

Figure 3(b) and (b*) shows the intensities of the focused fields E _inc in the plane of the scatterer that is located in the center of the panels. The radial mode M1 has a maximum intensity at the particle position R ₀, where the electric field points in the z-direction, contrary to the azimuthal mode that has zero electric field intensity there. The intensity maxima of M2 and G are located at R ₀, and the electric field lies in the xy-plane and has a circular polarization. In the panels, we also report the relative importance of the eigenvalues D = 100 Λ/Λ_max, with a factor of D = 100 corresponding to the largest eigenvalue.

The dependence of the eigenvalues (that equals the quantum Fisher information) on the numerical aperture is shown in Figure 3(c) and (c*). For NA values below say 1.3, the circular mode M2 performs best; however, the Gaussian mode with our ad-hoc choice of a unit variance performs almost equally well. This is remarkable. We have developed a framework for optimizing the quantum Fisher information, enabling the computation of fields that maximize localization precision, and have not imposed any constraints on the shape of the incoming fields ϵ _inc. Also their polarization properties have been left arbitrary. Nevertheless, the diagonalization approach comes up with four “boring” modes to be classified in terms of radial, circular, and azimuthal polarization. And for NA ≤ 1.3, the performance is even comparable to that of a simple Gaussian beam. Alas, we had hoped for better.

Things change for larger NA values, in particular for the scatterer on top of the substrate. Here, the radial mode M1 has the largest eigenvalue, which further increases for NA values that are large enough to support incoming waves with propagation angles above the cutoff value for total internal reflection (TIR). These wave components generate evanescent fields above the glass–water interface, which induce a strong dipole moment of the scatterer in z-direction and correspondingly an efficient scattering into the substrate. Indeed, when comparing panels (c) and (c*) for the scatterer on top and above the interface, we observe that for the large gap distance of 1 μm evanescent waves play no role and the eigenvalues remain bound to those NA values where TIR waves can be excited. It is thus the appearance of evanescent waves that boosts the localization precision. We note that evanescent waves not only play a role for radial polarization but also for all other modes, see Figure 3(c), although the localization enhancement is by far the strongest for the M1 mode.

In Appendix A, we show that the optimized fields also perform well for particles located in the vicinity of the spot for which they have been optimized. We also demonstrate that the precise value of the scaling parameter γ is not overly important.

Figure 4 provides a deeper insight into the contributions to the quantum Fisher information. The transition matrix in Eq. (12) can be conveniently decomposed into two contributions

where T _inc( θ ) accounts for the transition from the incoming fields to the scatterer, and T _sca( θ ) for the transition from the scatterer to the backfocal plane. The terms entering Eq. (15) can then be decomposed into the contributions

Figure 4:

Incoming, scattered, and remaining contributions of quantum Fisher information for different modes and NA values, see Eq. (19), and for (M1) radial, (M2) circular, and (G) Gaussian mode. We show results for the scatterer located (a) on top and (b) 1 μm away from the interface. For details, see text.

The expression in the first line accounts for the change of the incoming fields upon variation of θ _i while keeping the scattering properties fixed (incoming contribution). This contribution is expected to be large for excitation fields with high spatial gradients, where a slightly displaced particle is excited differently. The expression in the second line accounts for the modified scattering properties for fixed incoming fields (scattering contribution). Likewise, the expressions in the third and fourth lines account for mixed contributions (rest). The relative importance of the respective contributions is reported in Figure 4. Importantly, for the best modes (M2 for small NA values, M1 for larger NA values), the largest contribution is the scattering one. This shows that the high localization precision of these modes is not due to the fact that the scatterer is excited by a highly structured light field where the excitation properties change significantly when varying the particle position. Rather the scattering channel plays the important role, where a change of the particle position R ₀ leads to far-field modifications that carry the information about R ₀. For large NA values and the particle on top of the substrate, also the evanescent fields at R ₀ are transformed into waves propagating toward the lens, which then carry additional information about R ₀.

The above discussion explains why certain modes are better suited for particle localization than others: the optimized modes essentially maximize the field strength (and correspondingly the induced dipole moments) at the scatterer positions. Depending on the NA value of the focusing lens, the induced dipole moments are either located in the interface plane (M2) or perpendicular to it (M1). At the same time, the distribution of the incoming fields is surprisingly simple and could be approximated by even simpler Gaussian mode profiles with circular or radial polarization without significantly deteriorating the localization precision.

3.3 Optimization of Fisher information

So far we have computed the incoming fields that optimize the quantum Fisher information of Eq. (14). We should stress that our approach is quite general, and correspondingly the localization precision bounds obtained are generally valid for any type of coherent localization microscopy. However, it is not clear whether these bounds can indeed be reached in actual experiments and, if this is the case, which technique and localization protocol should be used. In this section, we make a first step in answering these questions and evaluate the localization performance of iscat. We start by comparing the quantum and iscat Cramér–Rao bounds, Eqs. (7) and (10), using the fields obtained from the optimization of the quantum Fisher information, and then continue to optimize directly the iscat Fisher information of Eq. (6). In our following discussion, we only present results for a scatterer located on top of the interface. As has been shown above, for smaller NA values, say below 1.3, the results for scatterers on top and away from the interface are very similar. In setups with larger NAs additionally evanescent waves can be excited, which are of importance only for scatterers on top of the interface.

Figure 5 shows results for the iscat Cramér–Rao bounds, which are obtained by starting from the optimized solutions of the quantum Fisher information and maximizing the Fisher information in Eq. (5) using an iterative procedure. Similar results were also obtained when starting from a random guess for the incoming fields. Let us first concentrate on the dashed lines in panels (c), which report the NA dependence of σ inc iSCAT for the localization precision along x (left) and z (right), using as incoming modes the solutions of the eigenvalue problem of Eq. (14) shown in Figure 3(a). With the exception of the dip around NA = 1.16, which is marked with an arrow and can be associated with waves close to the Brewster angle (a detailed discussion can be found in Ref. [28]), for NA values below 1.3 both the M2 mode with circular polarization and the Gaussian mode perform best. Things change for larger NA values, where the mode M1 with radial polarization achieves the best localization precision along x, however, at the expense of a lower precision along z. As has been discussed previously, this can be attributed to the excitation of evanescent waves at the glass–water interface.

Figure 5:

Optimization of the Fisher information for particle localization in iscat, see Eq. (6), and for a scatterer located on top of a glass–water interface. Electric field intensity in (a) backfocal plane and (b) plane of scatterer, and for the different NA values reported in the panels. (c) Localization precision along x (left panel) and z (right panel) as a function of the NA. The solid lines report the quantum Cramér–Rao bounds and the dashed lines the iscat Cramér–Rao bounds, as obtained from Eq. (7) for the optimized fields computed from the solutions of Eq. (14). The symbols show the bounds as obtained from the optimization of the Fisher information (γ = 0.2). For details, see text.

When comparing the quantum (solid lines) and iscat (dashed lines) Cramér–Rao bounds, one observes for small NA values a ratio of about three for the localization precision along x (left) and a significantly smaller ratio for z (right). This finding is in agreement with those for plane-wave excitations [13], [40]. For NA values above 1.3, the ratio of the localization precision along x decreases, with values of 2.2 for NA = 1.4 and 1.7 for NA = 1.45. At the same time, the localization precision along z deteriorates.

We additionally performed optimizations for the iscat Fisher information given in Eq. (6). In the optimizations, we expanded the incoming fields in a basis obtained from the solutions of Eq. (14), keeping the 25 modes with highest eigenvalues, and optimized the expansion coefficients in an iterative procedure such that the Fisher information becomes maximized. We started the optimization either with the mode associated with the highest eigenvalue or with random expansion coefficients, obtaining similar results in both cases. Quite generally this doesn’t mean that our optimization approach provides us with the global maximum; however, it shows that it is not overly sensitive to the initial guess for the incoming fields.

The symbols in Figure 5 show the Cramér–Rao bounds obtained from these optimizations. One observes that σ inc iSCAT ( x ) can be significantly reduced through optimization while keeping the precision along z almost unaltered, at least for NA values below 1.3. The optimized incoming field intensities are depicted in (a) the backfocal plane and (b) the plane of the scatterer. Closer inspection shows that for all NA values only a few basis modes with identical polarizations become mixed. As we don’t enforce any symmetry constraints in our optimizations, the modes originating from the iterative maximization procedure are slightly distorted and rotated. We have currently no clear explanation why these optimized modes perform better and cannot exclude that modes with an even better performance exist and will address these points in future work. A brief discussion of the scaling parameter γ in the context of the optimization results is given in Appendix A.