High precision wavefront control in point spread function engineering for single emitter localization

Point Spread Function (PSF) engineering is used in single emitter localization to measure the emitter position in 3D and possibly other parameters such as the emission color or dipole orientation as well. Advanced PSF models such as spline fits to experimental PSFs or the vectorial PSF model can be used in the corresponding localization algorithms in order to model the intricate spot shape and deformations correctly. The complexity of the optical architecture and fit model makes PSF engineering approaches particularly sensitive to optical aberrations. Here, we present a calibration and alignment protocol for fluorescence microscopes equipped with a spatial light modulator (SLM) with the goal of establishing a wavefront error well below the diffraction limit for optimum application of complex engineered PSFs. We achieve high-precision wavefront control, to a level below 20 m$\lambda$ wavefront aberration over a 30 minute time window after the calibration procedure, using a separate light path for calibrating the pixel-to-pixel variations of the SLM, and alignment of the SLM with respect to the optical axis and Fourier plane within 3 $\mu$m ($x/y$) and 100 $\mu$m ($z$) error. Aberrations are retrieved from a fit of the vectorial PSF model to a bead $z$-stack and compensated with a residual wavefront error comparable to the error of the SLM calibration step. This well-calibrated and corrected setup makes it possible to create complex `3D+$\lambda$' PSFs that fit very well to the vectorial PSF model. Proof-of-principle bead experiments show precisions below 10~nm in $x$, $y$, and $\lambda$, and below 20~nm in $z$ over an axial range of 1 $\mu$m with 2000 signal photons and 12 background photons.


Introduction
The diffraction limit to resolution is overcome in Single Molecule Localization Microscopy (SMLM) by estimating the location of individual molecules from sparsely distributed emission spots across the field of view of the camera [1,2,3,4]. The achievable resolution is limited by the localization precision and the density of fluorescent labels and can be on the order of several tens of nanometers in practice [5]. Several groups have extended this technique to 3D 1 localization. Estimation of the axial position of the molecules is made possible by adding a cylindrical lens to the optical path [6,7,8] or by using bi-plane or multi-focus imaging [9,10,11]. The elongation of the spot (astigmatism) or the difference in spot size (bi-plane imaging) encodes for the axial position. The addition of a Spatial Light Modulator (SLM) incorporated with a 4F relay system to the imaging light path enables the design of a broader class of Point Spread Functions (PSFs), opening up the possibility to optimize the localization performance of the microscope [12]. Notable proposals in the literature are single and double helix PSFs derived from Gauss-Laguerre modes [13,14,15,16], or made with annular zones with increasing helical charge [17,18], saddle point or tetrapod designs [19], which use higher orders of astigmatism in addition to primary astigmatism, phase-ramps [20] and self-bending beams [21]. Interferometric imaging can also be used to establish high-performance axial localization, but requires a complex 4π optical setup [22].
Other properties of the fluorescent molecules could be of interest next to the 3D-position, such as the emitter dipole orientation [23] and the wavelength of the emitted light. The latter is relevant for imaging multiple protein species in a specimen that are labeled with fluorophores with different (excitation and) emission spectra. Encoding the emission color into the PSF shape enables multicolor imaging with a single imaging light path and a single camera, operating at the full field of view. This has been proposed by our group in [24] for 2Dlocalization and generalized in [25] for 3D-localization. An alternative approach to '3D+λ' localization has been described in [26] by Shechtman et al.. The key idea of [24,25] is to have the SLM function as a Diffractive Optical Element (DOE), which splits the emission spot into two or three sub-spots. These sub-spots correspond to the diffraction orders generated by the DOE, and the distance between the spots and their relative intensity provide information on the emission wavelength of the fluorophore. The shape of the repetitive zones of the DOE can be designed for making the shape of the sub-spots change with the axial position of the fluorophore, thereby enabling 3D-localization.
A common characteristic of all engineered PSFs is their complexity compared to the simple 2D focused spots, which must be represented in the PSF model that is used in the parameter fitting algorithm for estimating the 3D-position (and possibly the emission color or molecular orientation). Simplified PSF models such as the Gaussian model [27], the scalar diffraction based Airy model, the Gibson-Lanni model [28], or effective models based on Hermite functions [29] cannot meet this requirement. A solution is the use of an experimental reference PSF, or a spline fit of such a PSF as model PSF [30,31,32,33], or the use of one or multiple Look Up Tables (LUTs) to estimate the z-position [34]. We have shown previously that a vectorial PSF model can also be used for complex 3D and 3D+λ engineered PSFs [25]. It is known that the vectorial PSF model is the physically correct model for image formation in high-NA fluorescence imaging systems. Another common characteristic of complex engineered PSFs is the sensitivity to aberrations that perturb the designed PSF shape and in this way negatively affect precision and accuracy. In order to achieve precisions down to the Cramér-Rao Lower Bound (CRLB), the best possible precision for an unbiased estimator, the aberration level of the optical system should be controlled to well within the diffraction limit (0.072λ root mean square wavefront aberration) [25], a condition which is often not met in practice. Correction of aberrations using a deformable mirror or with the SLM that is present anyway for producing the engineered PSF is therefore required. The control parameters of the adaptive optics component can be set using image based metrics [35,36,37] or via measurement of the to-be corrected aberrations. The latter may be done via phase retrieval algorithms based on the introduction of phase diversity, often in the form of a through-focus bead scan. This has been implemented in high numerical aperture microscope system [38], in localization microscopy [39] and used to improve the quality of a STED laser focus [40]. These algorithms rely on smoothing and apply constraints to suppress noise effects. The noise statistics however can be incorporated into a Maximum Likelihood Estimation (MLE) based phase retrieval algorithm, which makes it possible to assess the optimality of the retrieval process by means of a CRLB analysis [41].
Here, we present a calibration and alignment protocol for fluorescence microscopes equipped with an SLM with the goal of establishing a wavefront error well below the diffraction limit, which is needed for an optimum application of complex engineered PSFs in single emitter localization. This high-precision wavefront control is combined with the use of the vectorial PSF model, both in the MLE-based phase retrieval algorithm for estimating the aberrations, and in the localization algorithm for estimating the emitter position in 3D, signal photon count, and background photon level (and possibly emission color). We report on proof-of-principle experiments on fluorescent beads for different 3D+λ engineered PSFs for testing whether the proposed calibration protocol and fitting with a vectorial PSF model give rise to precisions matching the CRLB.
The outline of this paper is as follows. In section 2 the vectorial fitting routine and aberration retrieval is described, section 3 describes the experimental setup and the calibration and alignment techniques, section 4 shows the experimental results on the aberration retrieval and correction as well as on different 3D+λ engineered PSFs, and section 5 summarizes the main conclusions of the paper.

Vectorial PSF model for an emitter with non-zero size
The vectorial PSF model [42,25] is further refined by taking into account the non-zero size of the fluorescent beads we use in the experiments (175 nm and 200 nm), which induces a blurring on the scale of 1-2 pixels. We describe the PSF of such a bead as a PSF of a freely rotating single dipole emitter convoluted with a sphere ( r): with H bead ( r) the PSF, N the total number of signal photons detected at the camera, b the number of background photons per pixel and a the pixel size of the camera pixel, ( r) = 1 for | r| < r bead and 0 otherwise with r bead the radius of the bead and w lj ( r) the electric field component in the image plane, which is given by: with c n a normalization factor, A ( ρ) the amplitude, q lj ( ρ) the polarization vector components, defined elsewhere [43], and W ( ρ) the aberration function, which is the sum of a contribution from the microscope system and a contribution from the SLM. The wavevector k( ρ) is a function of the pupil coordinates: with n the refractive index of the medium and NA the numerical aperture of the objective lens. The expected photon count µ k for a given camera pixel is given by the integration of the PSF over the pixel area D k with size a × a: with r 0 the position of the emitter. The derivatives with respect to the fit parameters (θ = x, y, z, λ) which are needed for the fitting routine are similar to [25] but now involve a convolution: The derivatives of w lj with respect to the fit parameters remain the same: and: 2.2 Aberration retrieval from a through focus PSF scan Aberration retrieval is done by adapting the vectorial fitting routine to estimate the aberration coefficients from a through-focus image stack of a fluorescent 4 bead. The aberration function is expressed as a linear sum of root mean square (rms) normalized Zernike polynomialsZ m n ( ρ): where the appearing Zernike coefficients A m n are the fit parameters. The derivatives of the electric field components with respect to the Zernike coefficients needed for the MLE fitting routine are found to be: 3 Experimental setup and methods

Experimental setup
The setup, shown in Fig. 1, consists of a Nikon Ti-E microscope with a TIRF APO objective lens (NA = 1.49, M = 100), a 200 mm tube lens, a relay system with an SLM is built on one of the exit ports of the microscope. The relay system consists of two achromatic lenses (f 1 = 100 mm, Thorlabs AC254-100-A and f 2 = 200 mm, Thorlabs AC508-200-A), a nematic Liquid Crystal On Silicon (LCOS) SLM (Meadowlark, XY-series, 512x512 pixels, pixel size =15 µm, design wavelength = 532 nm) and a polarizing beam splitter to filter the x-polarized light which is not modulated by the SLM. The first achromatic lens relays the light on the SLM in a beam with a diameter of 3 mm and the second relay lens ensures Nyquist sampling of the fluorescent objects at the EMCCD (Andor iXon Ultra -X987, 512x512 pixels, pixel size = 16 µm, backprojected to object space 80 nm). The microscope is equipped with a set of lasers with wavelengths 405 nm (Coherent Cube), 488 nm (Coherent Sapphire), 561 nm (Coherent Sapphire), and 642 nm (MPB Communications). Either a dichroic filter set for the green (Ex: Semrock FF01-460/60-25, Di: Semrock Di02-R532-25X36, Em: Semrock FF01-545/55-25) or a quadband dichroic filter set (Chroma -TRF89902) is used. This standard setup is augmented with a novel, second light path for calibration of the SLM. This SLM calibration light path is designed for measuring the retardation difference between the x and y-polarized light incident on the SLM and consists of a laser (Thorlabs -CPS532, wavelength 532 nm, 0.9 mW) which illuminates the SLM via a beam expander, a polarizing beam splitter and a λ/2 wave-plate. The reflected and polarization-modulated light passes the λ/2 wave-plate and polarizing beam splitter again and is imaged onto a CMOS camera (Thorlabs -DCC1545M, 1280x1024 pixels, pixel size = 5.2 µm). A rotating diffuser is added to reduce speckle and two linear polarizing filters, which are aligned with the polarization axes of the polarizing beamsplitter, are added to reduce internal reflections. The intensity image captured by the CMOS camera is mapped to the intensity pertaining to specific SLM pixels I pxl . The polarization transfer through the λ/2 wave-plate and the SLM for the calibration light path is described by the Jones-matrix: with α the angle of the λ/2 wave-plate and φ pxl the retardance induced by an SLM pixel, and gives rise to an intensity: with I 0 the incident illumination intensity. The waveplate angle α is set to a small value, around 5 deg, in order to have a small maximum transmission. This is required for having a relatively long integration time, about an order of magnitude more than the rotation period of the rotating diffuser, for enabling good speckle reduction.

SLM calibration
In order to measure the modulation of a certain SLM pixel the mapping from the SLM onto the camera of the calibration path is needed. This mapping is obtained by applying a checkerboard pattern with increasing voltages to the SLM. The difference between the average captured image and the image when no voltage is applied is used as input for a corner detection algorithm (findcheckerboard from Matlab -Mathworks) to find the corner points. An affine transformation is fitted to these points and used to find the CMOS pixels corresponding to each SLM pixel. The calibration procedure for an SLM pixel is graphically explained in Fig. 2. First, the intensity response as function of applied voltage is measured in 256 steps, giving rise to a sequence of minima and maxima, which correspond to a retardation of π or 2π. All pixels inside the illuminated SLM plane appear to have three maxima, implying a total phase modulation of 4π or 1094 nm. The voltages for which these extrema occur are found by fitting parabolae to the three points near the extrema, which increases the precision and fully utilizes the 16 bit control of the SLM. The intensity is then divided into four segments which are scaled to [0 1] and converted to phase using the inverse of Eq. (11) over these segments. The phase response is used to construct an individual Look Up Table (LUT) for each SLM pixel, compensating the non-uniformity of the SLM. The LUT-parameters vary smoothly over the SLM and correspond roughly 6 with the Fabry-Perot fringes visible by eye, indicating that the differences in phase response are due to variations of the thickness of the liquid crystal layer. Additional pixel-to-pixel variations may arise from pixel-to-pixel variations in the underlying silicon switching circuitry. The complete calibration takes about 5 minutes (3 minutes scanning and 2 minutes computing time on a quadcore 3.3 GHz i7 processor), but can in principle be optimized to run faster.
The calibration is verified by applying a uniform retardation profile over the complete calibration range to measure the root mean square error of the reflected wavefront and compare the average intensity to the ideal intensity response. This verification is performed immediately after calibration and also after 45 minutes, and subsequently compared to the calibration provided by the manufacturer, which is only for 2π modulation, see Fig. 2(F). The rms error of the wavefront with the manufacturer's LUT is around 100 mλ and within specification, but too high for the stated goal of our research: to achieve wavefront control well below the diffraction limit in order to optimally apply complex engineered PSFs in single emitter localization. The individual pixel calibration method reduces the wavefront error by an order of magnitude to around 10 mλ immediately after the calibration and to 20 mλ after 45 minutes. The deterioration of performance over time is attributed to temperature fluctuations and the associated mechanical drift of the different optical components. We conclude that there is about a half hour window of opportunity to conduct experiments with a precision of the wavefront control within 20 mλ rms wavefront aberration. All subsequent experiments have been done within this time frame. We mention that dedicated temperature control of the setup can extend the total time with high-precision wavefront control.

Axial and lateral alignment of the SLM
The position of the Optical Axis (OA) of the emission path on the SLM is needed in order to project the phase profile at the correct position on the SLM. Furthermore the SLM should be aligned with the plane conjugate to the pupil plane of the objective lens, the Fourier plane, to ensure that every point inside the Field Of View (FOV) has the same phase modulation. A procedure is used to directly estimate this alignment by imaging a set of beads in the FOV and applying a defocus modulation, without the need for additional optical components. For a bead in the center of the FOV the chief ray of emitted light beam aligns with the OA and intersects the SLM plane at a position (ρ OA x ,ρ OA y ) with respect to the coordinate frame that has its origin at the center of the SLM area designated for use in the aberration correction and PSF engineering. Here ρ OA x and ρ OA y are dimensionless pupil coordinates (real space coordinates normalized by the pupil radius). If the SLM introduces defocus then this parabolic aberration profile will be decentered w.r.t. the OA: where A 0 2 is the Zernike fringe (not rms normalized) coefficient for defocus, and where W has length units. Apparently, not only defocus is introduced to the overall imaging system but also tip and tilt. As a consequence, the experimental PSF obtained from bead images will shift laterally, where the shift varies linearly with the applied defocus according to: (and similarly for the y-direction), as follows from comparing the tip/tilt in Eq. (12) to the k ( ρ) · r term in Eq. (2). Therefore, the lateral alignment of the SLM can be estimated from the shift of a bead in the center of the FOV and the axial alignment can be estimated from the shift of the beads at the rim of the FOV. The center of these beams will intersect the SLM at a different lateral position if the SLM is not aligned with the Fourier plane as illustrated in

Correction of oblique angle of incidence at SLM
The beam is reflected by the SLM under an oblique angle of θ = 20 degrees. The aberration function W (x, y) added by the SLM is related to the aberration function W slm (x , y ) under normal incidence, where (x, y) are the pupil coordinates normal to the beam axis and (x , y ) are the coordinates in the SLM plane, by: with θ the corresponding angle of incidence inside the SLM (sin θ = n sin θ , with n the refractive index inside the SLM, taken to be n = 1.5, and where the birefringence within the liquid crystal is neglected). Equation (14) describes two effects. The first is the scaling of the phase depth due to the oblique incidence, the second is the projection of the circular pupil onto the SLM plane, giving an anisotropic stretch and hence an elliptical cross-section (aspect ratio cos (20 deg) = 0.94). Implementing this transformation ensures that the phase profile contributed by the SLM corresponds to the aberration function for normal incidence used in the localization algorithm. The relative large pitch is sufficient to make a rough estimate of the emission spectra, but is not comparable to the spectral resolving power of a spectrometer. The spectral broadening is mainly due to the diffraction limited spot size on the camera, leading to an apparent non-zero emission in the bandstop regions of the dichroic. The spectral resolution is estimated to be on the order of 40 nm (estimate obtained from the ratio of the 0th order spot size to the distance between the 0th and 1st order times the peak emission wavelength). Blazed grating profiles are applied in the x and y-direction (the direction of incidence is tilted in the x-direction at the SLM) giving rise to substantially the same emission spectrum, thereby confirming that the anisotropic stretch of the SLM aberration function described by the obliquity correction of Eq. (14) is correct. The emission peaks are found at around 520 nm and 680 nm, and the weighted average emission wavelengths are 536 nm and 693 nm, for the green and red beads, respectively, matching the bead specifications. The path length modulation of the SLM pd is calibrated by applying a blazed grating with increasing phase depth. This results in diffraction orders m = 0, 1, . . . with amplitudes: with sinc (x) = sin (πx) / (πx), giving an intensity ratio between the zeroth and first diffraction order: By measuring this intensity ratio it is possible to estimate the phase depth modulated by the SLM and compare this to the expected applied phase depth modulation as shown in Fig. 4(C). This results in a phase depth modulation which is 3.8% higher for the green bead and can be regarded as a small correction of the 1/ cos(θ ) term in the oblique angle correction. The measured phase depth in the deep red is 74% of the phase depth expected from the SLM calibration experiments. This may be due to fringe field effects between the pixels [44,45,46], possibly in relation to the dielectric mirror and coatings that are optimized for green light, to inherent chromatic dispersion of the liquid crystal and/or to spurious reflections in the system contributing to the 0th order. Both differences in modulation are incorporated into the MLE-based localization/wavelength fitting routine by correcting the phase modulation with a factor 1.038 and 0.74 for the green and red bead, respectively. The observed variation of effective phase depth with wavelength is not incorporated into the calibration of the SLM for pixel-to-pixel variations with the additional calibration light path.
4 Experimental results

Aberration retrieval
The aberration retrieval and subsequent correction is tested experimentally for beads emitting in the green, using the dichroic filter set for green emission. In particular, the dominant Zernike aberration modes, the fit precision, and the goodness of fit is assessed. To this end through-focus PSF stacks of 21 slices in a 2 µm range are recorded, converted to photon counts with a gain calibration procedure, and fitted on a 31×31 pixel Region Of Interest (ROI), where the 3Dposition of the bead, the number of signal photons, the number of background photons per pixel, and a set of Zernike aberration coefficients are used as fit parameters in the MLE fitting routine. Measurements are repeated 5 times for determining the precision. The fitted signal photon count per focal slice is around 2.4×10 4 , the fitted background photon count around 19 photons/pixel. Figure 5(A) shows the retrieved Zernike aberration coefficients for fitted modes (n, m) satisfying n + |m| ≤ 2 (j + 1) with j = 1, 2, 3, 4, 5 the fit order (j = 1 includes primary astigmatism, coma, and spherical aberration, j = 2 includes secondary astigmatism, coma, and spherical aberration, and primary trefoil, etc.), Fig. 5(B) shows the corresponding retrieved wavefronts. The Zernike coefficients found at fit order j match well with the ones found at a lower fit order j − 1, indicating that there is little cross-talk between Zernike modes in the MLE fitting routine, and that there is no overfitting with large numbers of fitted aberration coefficients. Even for the case j = 5 for which 45 modes are retrieved, there seem to be no problems with convergence (typically only about 6 iterations are needed), and reproducibility of the fit. The dominant Zernike modes appear to be primary astigmatism (A 2 2 ), secondary coma (A 1 5 ), and two higher orders spherical aberration (A 0 6 , A 0 8 ). The correction collar of the objective lens is used to reduce the primary spherical aberration (A 0 4 ) as much as possible, but this does not seem to compensate for the higher orders of spherical aberration.
The experimental fit precision for all Zernike modes is below 1.5 mλ and typically around 0.6 mλ, somewhat higher than the CRLB, which is around 0.3 mλ [ Fig. 5(C)]. The impact of photon count on the fit precision of the aberration retrieval is further assessed with a simulation study. To this end through-focus PSF stacks of 21 slices in a 2 µm range are simulated where the total rms aberration level is kept constant at 10, 45 and 80 mλ. We use N cfg = 100 random instances per aberration level. The set of aberrations used consists of all Zernikes modes of radial order n ≤ 4 (except piston, tip, tilt, defocus). Next, shot noise is added corresponding to a range of signal photon counts N ph and a background of b = 10 photons per camera pixel, and these sets of noisy through-focus PSFs are used as input for the MLE aberration fitter. The fit precision is quantified by the average rms error of the wavefront: and the quality of the fit is evaluated by comparing this fit precision to the CRLB. Fig. 5(D) shows the result, and indicates that the precision scales as 1/ N ph , in agreement with expectations. The average residual wavefront error of the fit appears to be drop below 1 mλ for a through-focus stack with more than 10 4 signal photons, corresponding to the experimental conditions in the aberration retrieval tests. The goodness of fit is estimated with a chi-square test. The chi-square statistic is defined as: where n k is the measured photon count per pixel in each focal slice, and µ k is the photon count expected from the fit model. Here K is the total number of pixels in the fit region times the number of focal slices, i.e. the total number of statistically independent measurements. If the n k follow a Poissonian distribution with rates µ k then the mean and variance of the statistical distribution of χ 2 values follow as: var where the expectation values (n k − µ k ) 2 = µ k and (n k − µ k ) 4 = µ k +3µ 2 k for the Poisson-distribution are used. The statistical distribution of χ 2 values may be approximated by a normal distribution as K 1, even though the statistics of each measured pixel is Poissonian. The goodness of fit can then be quantified by the level of confidence found by comparing the experimental χ 2 -value to the mean and standard deviation of the expected normal distribution of χ 2 -values. Figure 5(E) shows the measured χ 2 -values in relation to the expected value K = 21×31 2 = 2.0×10 4 . It appears that the χ 2 -value converges to a level about 20% higher than the expected value, significantly more than the expected standard deviation 2×10 3 . This shows that there are still some model errors left, probably effects of photobleaching and illumination intensity variations (it is assumed in the fit that the same level of signal and background photon count applies to each focal slice). Other effects which possibly contribute to the discrepancy are the non-zero spectral bandwidth of the collected fluorescence emission, residual scattering at the SLM, and amplitude aberrations or apodization due to e.g. variations in transmission through the objective lens and the dichroics with pupil position [39].

Aberration correction
The overall rms wavefront error, as determined from the fitted Zernike coefficients converges to a value around 65 mλ, just below the diffraction limit, with increasing number of Zernike modes [ Fig. 5(F)], proving that aberration correction is needed for successful application of complex engineered PSFs. There is a trade-off between the number and type of Zernike modes that can be corrected and the peak-valley value of the wavefront that is needed for the desired engineered PSFs, in view of the limited phase dynamic range of the SLM. In the experiments on engineered PSFs we take all Zernike modes with radial order n ≤ 4 into account, similar to previous studies in the literature [18,41]. In order to evaluate the best possible performance of the aberration correction we also include second order coma (A 1 5 and A −1 5 ), second order spherical aberration (A 0 6 ), and third order spherical aberration (A 0 8 ). Figure 6 and Visualisation 1 and Visualisation 2 show the experimental and fitted through-focus PSF without and with aberration correction. The agreement between measured and fitted PSFs is quite well, especially in the 1 µm range around focus. It appears that the experimental through-focus PSF is more rotationally symmetric after correction, proving that asymmetry inducing aberrations are reduced. In addition, the asymmetry between spot shapes above and below focus is greatly reduced after correction, indicating that spherical aberration is largely eliminated. The wavefront error estimated with the aberration retrieval algorithm reduced from 59 ± 1 mλ to 13.4 ± 0.4 mλ [ Fig. 6(G)]. This residual error is close to the calibration precision of the SLM, indicating that the level of aberration correction is limited by the precision of the SLM control.
A final test of the aberration retrieval and correction procedure is performed by deliberately adding single Zernike modes (A m n = 60 mλ) on top of the corrected wavefront and subsequently feeding the aberrated through-focus stack to the aberration fitting routine, see Fig. 7 for the results. The fitting routine correctly retrieves each aberration with an estimated value of 67 ± 4 mλ, averaged over the 15 displayed aberrations, somewhat higher than the expected 60 mλ. All other retrieved aberrations remain at the level of 20 mλ or less, pointing to the specificity of the aberration retrieval and correction procedure. In particular, there is little crosstalk from added aberrations of radial order n to retrieved aberrations of order n − 1, which confirms the correct lateral alignment of the SLM phase profile with respect to the optical axis.

Analysis of engineered PSFs
Proof-of-principle experiments of engineered PSFs have been performed with green and red emitting beads, using the quadband dichroic filter set. Aberra-tions are corrected only up to radial order n ≤ 4 in order to save phase dynamical range on the SLM for the engineered PSFs. Primary spherical aberration is compensated by the correction collar of the objective lens. Three different designs have been tested. The first is a binary grating splitting the spot into two ±1st diffraction orders, where the grating zones are curved to induce astigmatism to the two orders (grating zone shape described by Zernike coefficients A −1 1,zone = 0.8λ 0 and A −2 2,zone = 0.15λ 0 , as in [25]). The nominal wavelength λ 0 is equal to 520 nm for the green emitting beads and 690 nm for the red emitting beads.The second engineered PSF is a blazed grating for splitting the spot into a 0th and +1st order (grating zones curved to induce astigmatism, the shape is described by Zernike coefficients A −1 1,zone = 1.4λ 0 and A −2 2,zone = 0.3λ) and an overall continuous astigmatic aberration profile with Zernike coefficient A −2 2,overall = −0.15λ 0 . The third engineered PSF is a double helix configuration (annular design [17,18] with four rings and exponent α = 1/2, and phase depth = λ 0 ). These design parameters are set by balancing the achievable precision with the axial range and with the footprint of the spot on the detector, but can in theory be improved by e.g. incorporating higher order astigmatism in case of the astigmatic profiles (mimicking the saddle-point or tetrapod PSF [12]). Through-focus image stacks of these engineered PSFs are recorded for different signal photon counts while keeping the background constant using the camera frame-time and the trans-illumination unit for conventional brightfield microscopy for providing extra background photons at shorter frame times as tuning parameters, and subsequently fitted using MLE and the vectorial PSF model. Figure 8 shows the phase profile for the three engineered PSF designs and examples of measured spots and the corresponding fits for two estimated signal photon counts (5 · 10 3 and 14 · 10 4 ). The overall experimental PSF, obtained by adding all recorded spots after first upsampling 3× and subsequently shifting with the fitted xy-position of the bead, is shown as well, along with the prediction of the vectorial PSF model. In this way the experimental PSF is built up from the cumulative signal of ∼ 10 8 photons. The agreement between experiment and the theoretical vectorial PSF is generally excellent, even the fringe structures at the largest defocus values match very well. The remaining discrepancies, mainly a slight broadening of the spots, is attributed to the non-zero spectral width of the light incident on the camera, due to the width of the emission spectrum and the width of the bandpass regions of the quadband dichroic. There is also a small asymmetry in the fringe structure, which is probably caused by the residual higher order spherical aberrations in the optical system.
The achieved precision in x/y/z/λ for the three engineered PSFs in focus as a function of signal photon count are shown in Figs. 9(A)-9(D), and for one signal level as a function of the axial position in Figs. 9(E)-9(H). The localization precision is obtained by fitting 25 acquisitions of a single bead at each z-position. The precisions follow the CRLB as a function of signal photon count and axial position, with the exception of the z-precision, which is somewhat worse than the CRLB. This is attributed to the higher orders of spherical aberration, which are not corrected for these datasets. The localization precision values appear to level off at values around x/y/z = 2/2/5 nm for high photon counts, probably due to effects of drift. The overall performance of the three engineered PSFs concerning precision appears to be quite similar. Figure 9(F) shows the average fitted zposition as a function of stage z-position. Both astigmatic engineered PSFs have a linear response with a fitted slope of 1.03±0.02 (blazed) and 1.01±0.02 (binary) in the green, but the double helix PSF and the blazed astigmatic PSF in the red underestimate the z-position slightly with a slope of 0.94±0.02 and 0.93±0.03, respectively. This bias may also be due to uncorrected higher order spherical aberration. The estimated wavelengths shown in Figs. 9(J)-9(K) are close to the measured weighted average emission wavelengths and resolve the green and red bead excellently. There seems to be a small overestimation of the wavelength for the green bead, which is constant over the entire axial range for the binary astigmatic PSF and the double helix PSF, but not for the blazed astigmatic PSF. The accuracy of the z-position is most likely affected by these small over and underestimations of the wavelength, as these parameters are closely coupled in the fitting routine.

Conclusion
In summary, we have shown how a dedicated calibration protocol for high-NA fluorescence microscopes equipped with an adaptive optical element such as an SLM can reduce the aberrations to around 20 mλ rms wavefront aberration. This high level of wavefront control enables single emitter localization with complex engineered PSFs, where the full vectorial PSF model is used in the fitting routine. A key ingredient in the calibration protocol is a separate SLM calibration light path for suppressing the effect of the non-uniformity of the SLM. In these well-controlled experimental circumstances the experimental PSFs conform very well to the predictions of the vectorial PSF model. The method is further tested with proof-of-principle experiments for fitting the 3D-position and the emission wavelength of single emitters. The precision in x/y/z/λ is similar to the CRLB in all four fit parameters, and a high accuracy in z and λ is obtained, without an experimental reference PSF or LUT.
An open issue is the residual model mismatch, as revealed by the χ 2 -test. A next step could be the expansion of the fit model with an intensity scaling that varies with the focal slice in order to take into account effects of photobleaching and illumination intensity fluctuations. Another inroad is to include apodization or amplitude aberrations, which model the dependence of the transmission of the different optical components (objective lens, dichroics) on pupil position. The optimum number of fit parameters can be assessed with statistical methods for model selection, in particular methods based on optimizing the Akaike information criterion, basically a weighted sum of the number of fit parameters and the maximum log-likelihood obtained by the fitting routine [47].
Another next step in the technical developments described in this paper would be to fit only the 3D-position keeping the emission wavelength fixed to the known weighted average fluorescence emission wavelength. This could possibly further improve the z-accuracy and precision. The non-zero spectral width could be taken into account by averaging the single-wavelength PSF over the spectrum weighted by the fluorescence emission and dichroic bandpass filter. Such an approach would enable identification of fluorescent labels with different emission spectra using a Generalized Likelihood Ratio Test (GRLT) [48] by comparing the likelihood of the fit for the known peak wavelengths of the different fluorescent labels.
An improvement of the current experimental setup for photon starved applications would be the replacement of the polarization sensitive LCOS-SLM by a polarization insensitive adaptive optical element such as a Deformable Mirror (DM). A similar calibration and alignment protocol as described here could be used, including the use of a calibration branch to suppress effects of thermal drift in the DM. Another direction of future research is to use the SLM or DM to mitigate the effects of sample induced aberrations. This can be done most effectively when the aberrations primarily originate from a sufficiently thin layer in the sample volume (on the order of the focal depth) by placing the adaptive optical element at a plane conjugate to the aberration inducing layer in the sample [49].
Data and software for aberration retrieval and '3D+λ' PSF fitting is available at [50].   Figure 4: A) The emission spectra of the red and green beads are measured by applying a blazed grating to the SLM. The distance between the zeroth and first order on the camera is then a measure for the wavelength, the spread of the first order spot is a measure for the emission spectrum. The spectra for gratings applied in the x and y direction (full and dashed lines) are identical, confirming the oblique angle correction Eq. (14). B) Illustration of the blazed grating and pitch p and path length step pd. C) The phase depth is calibrated by fitting the measured intensity ratio between the zeroth and first order with Eq. (16).

23
A B C D F n+|m|£12 n+|m|£10 n+|m|£8 n+|m|£6 n+|m|£4 Figure 5: A) Fitted Zernike aberrations coefficients (rms values) of the non aberration corrected microscope for different sets of modes taken into account. Fits have been done for modes n + |m| ≤ 2 (j + 1) with j = 1, 2, 3, 4, 5 the fit order. The coefficients found from higher order fits match reasonably well with the values found from the lower order fits. B) The retrieved aberrated wavefronts according to the fitted Zernike coefficients. C) The fit precision found in experiment (full lines) and according to the CRLB (dashed lines). D) Fit precision in simulation (data points) and CRLB of the fit (solid lines) as a function of photon count for different total rms aberration levels. E) The χ 2 value of the fit as a function of fit order j, flattening off at a value about 20% higher than the expected value (dashed line). F) The rms level of the fitted Zernike coefficients as a function of fit order j converging to a value just below the diffraction limit (dashed line).  Figure 7: A) Through-focus stacks of green emitting beads deliberately aberrated by single Zernike modes with coefficients A m n = 60 mλ and the corresponding theoretical through-focus stacks using the aberration coefficients found from the aberration fitting routine (through-focus range: [-1, 1] µm, 21 steps, estimated photon count was around 2.2×10 4 signal and 32 background photons). Scalebar indicates 1 µm and all exp/fit image pairs are contrast stretched with the same scale. B) The aberration retrieval appears to be mode specific, and estimates the Zernike coefficients as 67 ± 4 mλ, averaged over the 15 displayed Zernike modes.