3D sub-diffraction imaging in a conventional confocal configuration by exploiting super-linear emitters

Sub-diffraction microscopy enables bio-imaging with unprecedented clarity. However, most super-resolution methods require complex, costly purpose-built systems, involve image post-processing and struggle with sub-diffraction imaging in 3D. Here, we realize a conceptually different super-resolution approach which circumvents these limitations and enables 3D sub-diffraction imaging on conventional confocal microscopes. We refer to it as super-linear excitation-emission (SEE) microscopy, as it relies on markers with super-linear dependence of the emission on the excitation power. Super-linear markers proposed here are upconversion nanoparticles of NaYF4, doped with 20% Yb and unconventionally high 8% Tm, which are conveniently excited in the near-infrared biological window. We develop a computational framework calculating the 3D resolution for any viable scanning beam shape and excitation-emission probe profile. Imaging of colominic acid-coated upconversion nanoparticles endocytosed by neuronal cells, at resolutions twice better than the diffraction limit both in lateral and axial directions, illustrates the applicability of SEE microscopy for sub-cellular biology.


Supplementary Note 1: Photo-stability
We studied the photo-stability of the UCNP emitters in two experimentally relevant scenarios. In both cases the particles were located inside the fixed neuronal cells which were used for the uSEE imaging in the paper.
In the first scenario, we scanned the sample through the excitation beam over a period of 5.5 hours at 0.4 mW excitation beam power at the back focal plane of the objective (2.2 mW µm −2 peak power density at the sample). This corresponds to imaging in uSEE microscopy mode with resolution two times better than the diffraction limit. We scanned under practically relevant conditions -dwell time of 3 ms, sampling rate of 100 × 100 pixels, field-of-view (FOV) of 1.67 × 1.67 µm. The emission from the particles in the FOV remained relatively stable, as only a minor downward trend in the emission intensity, on the order of 1% per each hour of imaging, was observed ( Supplementary  Figure 1 a).
In the second scenario, we scan at another practically relevant experimental setting, which produces diffraction limited images. This is achieved for excitation power of 2.2 mW at the back focal plane (11.8 mW µm −2 peak power density at the sample). We again scanned with the field-of-view of 1.67 × 1.67 µm, dwell time of 3 ms and a sampling rate of 100 × 100 pixels. Again, negligible deterioration of the emission intensity, on the order of 1% per each hour of imaging, was observed as seen in Supplementary Figure 1 b, which depicts the emission behaviour over a period of 7 hours.
It is also important to note that the imaging resolution remains the same throughout the whole duration of the imaging -215 nm for the uSEE imaging regime in Supplementary Figure 1 a and 420 nm for the diffraction-limited confocal imaging in Supplementary Figure 1 b. The alignment of the particle with respect to the beam has been verified at regular intervals. Supplementary Figure 1: The particle emission intensity decreases with a slow rate of 1% per hour, while the imaging resolution remains stable for practical imaging conditions. a Scanning an excitation beam of peak power density 2.2 mW µm −2 , in uSEE microscopy imaging conditions results in negligible drop of the emission intensity of the particle over a period of 5.5 hours. b Scanning an excitation beam of peak power density 11.8 mW µm −2 , in confocal mode corresponding to diffraction limited resolution, results in negligible decrease of the emission intensity of the particle over a period of 7 hours. In both (a) and (b) the resolution remains stable throughout the measurement.
Supplementary Note 2: Simulating the resolution of the optical system when imaging a non-linear emitter with an arbitrary excitation-emission curve Section A: General procedure for applying the theoretical framework In Supplementary Note 2, we detail on the theoretical framework, which we have developed, to numerically calculate the resolution of confocal microscopes employing non-linear emitters with a continually changing slope. The resolution is determined by the excitation beam shape and by the excitation-emission curve of the emitter. Due to the super-linearity of the emitter, small variations in these factors can lead to significant changes in the resolution. Thus, the excitation beam shape and the excitation-emission curve need to be precisely defined.
We adopt a fully vectorial model to calculate the excitation beam shape B ex (Supplementary Note 2, Section B). The model also accounts for key factors influencing the beam, namely the filling factor of the objective back-pupil and the polarization of the incident light. The resulting expression for the excitation beam shape (Eq. 15) is evaluated numerically in Matlab for the specific parameters of our setup (see Supplementary Note 2, Section F for details).
In addition to the excitation beam shape, a functional form of the excitation vs emission curve U is also needed to predict the resolution. We experimentally determine U m , as described in Supplementary Note 2, Section D and we find a functional form of U by fitting a smoothing spline in Matlab (see Supplementary Note 2, Section E for details).
Finally, equipped with the knowledge of U and the excitation beam shape B ex , we proceed to rule out the effect of the pinhole affecting the resolution in our system. The theory of the detection pathway filtering is presented in Supplementary Note 2, Section C, F along with the conclusion that for our specific experimental conditions, the excitation-emission curve U and the excitation beam shape fully determine the resolution of the system. Section B: Calculation of the excitation beam shape using fully vectorial approach We base our calculations on the established approach from Chapter 3 in Ref. [1] to calculate the electric field distribution in the focal region of an aplanatic, high numerical aperture (NA) microscope objective lens of focal length f (Supplementary Figure 2 a). In this section, we modify the approach in Ref. [1] to account for the circular polarization of the beam and for the specific filling factor of our microscope objective.
The light source in the form of a linearly polarized (along x direction), collimated Gaussian beam with the beam waist w 0 (see Supplementary Figure 2 a) can be expressed as where the negligible spreading of the beam along the propagation axis is omitted. Please note that we use (x , y ) to denote the coordinate system before the dichroic mirror (DM) and (x, y) for the coordinate system after the DM. The E 0 beam is passed through λ/4-waveplate with its optical axis at 45 • angle with respect to the direction of the incident polarization. Transforming the polarization of the E 0 beam from linear to circular is essential for the radial symmetry of the beam in the focus of high NA objectives. Using linear polarization leads to highly elliptic beams which results in asymmetric resolution in the lateral plane 1,2 . For this reason, we focus our subsequent analysis on the case of a circularly polarized beam and derive the equations describing the beam in the focus for this situation. The transformation of the beam E 0 by λ/4-waveplate can be mathematically expressed as Applying the transformation T λ/4 · E 0 leads to the following expression for E inc The circularly polarized beam E inc is then reflected on the dichroic mirror. If we assume that the dichroic mirror reflects the beam in a way that does not change the polarization state, the reflection is simply a coordinate mapping from (x , y ) → (x, y) and the beam incident on the objective lens becomes  Figure 2: Schematic diagram of a confocal microscope. In (a), we schematically show the geometry, parameters and coordinates used for the derivation of the excitation beam shape in the focus of an objective lens. The pinhole and the tube lens are displayed only for completeness as they have no role on the excitation beam path. In (b), we schematically show the formation of the dipole image in the detection plane for a dipole emitter positioned in the origin (0, 0, 0) and displaced by (ρ s , ϕ s , z s ) from the origin.
Here, we take advantage of the sine condition of geometrical optics (see Supplementary Figure 2 a and Ref. [1]) which relates the distance of the incident ray from the optical propagation axis (x 2 + y 2 ) and the angle of the corresponding refracted ray θ with respect to the same axis. As per Supplementary Figure 2 a, the incident and the refracted rays are related via In order to proceed further, we need to find an expression for E inc immediately to the outside of the spherical surface of radius f representing the objective lens. For the lens with a focal length f and after denoting c = (n x + in y + 0 n z ), with n x , n y and n z being the Cartesian coordinates unit vectors, we can express E inc in the spherical coordinates (f, θ, φ) with the origin in the focus of the aplanatic lens as Here, we have used the sine condition of geometrical optics x 2 + y 2 = f 2 sin 2 θ and projected E inc into vector components n ρ = cos φ n x + sin φ n y n φ = − sin φ n x + cos φ n y (7) as seen in Supplementary Figure 2 a. To find the expression of the field at the inside of the spherical surface we need to consider that the energy carried by the incident and the refracted ray must be the same. To satisfy this condition, the amplitude of the field at the inside of the lens has to be multiplied by a factor of n n √ cos θ, where n is the refractive index to the outside of the lens and n is the refractive index on the inside of the lens, respectively (see Ref. [1] for derivation of this condition). Here, we assumed that the permeability of the inside/outside media is equal to 1, which is a valid assumption for optical elements used in microscopes. We further assumed that the Fresnel transmission coefficients along n ρ and n φ are equal, which is a very good approximation for an air-glass interface of a lens. Finally, as seen in Supplementary  Figure 2 a, the spherical surface does the following mapping of the vectors (transition from outside to inside) where n θ = cos θ cos φ n x + cos θ sin φ n y − sin θ n z .
After taking the above into consideration, we can write the field immediately to the inside of the lens as: Note that we have not only expressed E lens on the spherical coordinate grid (f, θ, φ), but we have also transformed the vector components of E lens from cartesian to spherical coordinate system, centred in the focus of the aplanatic lens. Substituting from Eq. (7, 9) into Eq. (10) leads to Note that this final form of E lens is a vector field dependent on the spherical coordinate grid but with the vector components along the Cartesian directions (n x , n y , n z ). We now modify the expression to account for the filling factor f 0 of the back-aperture of the microscope objective lens, which influences both the lateral and the axial extend of the excitation beam shape. Large filing factors take advantage of the full NA of the microscope objective whereas small filing factors effectively decrease NA and therefore also decrease the resolution of the system. As the axial resolution improves roughly with the square of NA and lateral resolution improves linearly with NA, the filing factor can play a major role in the precise determination of the beam shape and the resolution of the microscope, respectively.
The aperture radius of the lens is given by f sin θ max , where θ max is the maximum acceptance angle of the cone of light entering the lens. The filling factor f 0 is defined as the ratio between the beam waist and the radius of the back-aperture of the lens f 0 = w 0 /f sin θ max . Inserting this into Eq. (11) leads to: where we defined the apodization function of the lens as The field E lens , just after the lens interface, is an angular spectrum representation of the field in the focus E exc , or in other words, E lens is an angular spectrum representation of the excitation beam E exc . The two fields are related via the Debye-Wolf diffraction integral 1 : where the origin of the cylindrical coordinates (ρ, ϕ, z) is in the focus of the aplanatic lens and k is the wave-number in the lens medium given by k = 2πn/λ 0 (λ 0 is the vacuum wavelength). After performing integration over φ, the focal field can be written as: where the functions I 0 , I 1 , I 2 are given by where J n is the n th -order Bessel function of the first kind. The integrals can be evaluated numerically to yield the excitation beam shape B ex (ρ, ϕ, z) = |E exc | 2 in the focus of a high NA microscope objective. We do this in the specific case of our setup in Supplementary Note 2, Section F.

Section C: Calculating the field in the detector plane
The precise knowledge of the excitation beam shape allows us to describe how strongly the particles get excited within each position of the beam ( Supplementary Figure 2 b). In the following we assume that the E exc (ρ, ϕ, z) induces an emitting dipole (emitter assumed infinitesimally small), with the dipole moment p aligned with the electric field of the excitation E exc : where α is the polarisability. Here, we normalized the magnitude of the dipole moment p with respect to the excitation intensity |E exc |, as we need to introduce the magnitude of the dipole emission at the later stage for the special case of a nonlinear emitter. The magnitude of the dipole emission depends on the excitation vs emission curve U , which encompasses the material interaction of the light with the particle. Therefore, we hypothesize that the amplitude of the radiating electric dipole moment for a case of a nonlinear emitter, having a continually changing slope, can be written as: and the nonlinear radiating dipole expression, with the emission strength taken into account, becomes The magnitude of the induced dipole moment and its orientation allows us to back-propagate the light emitted by the dipole though a microscope objective, and subsequently focus the light by a tube lens onto a detector plane.
The electric field in the detector space E 000 (ρ d , ϕ d , z d ) due to a radiating dipole in the focus of the aplanatic lens (0, 0, 0) (see Supplementary Figure 2 b) can be readily calculated using the dyadic Green's function ↔ G 0 (r, r 0 ) approach. The electric field E 000 (r) from a dipole p positioned at r 0 = (0, 0, 0) can be found via 1 : The procedure to find the dyadic Green's function for a case of radiating dipole embedded in medium of refractive index n, with the field collimated with a microscope objective lens of focal length f and focused by a tube lens of a focal length f onto a detector embedded in medium n is described in detail in Ref. [1]. For the purpose of brevity, here we only give the final result which is calculated for the assumption of f f (see Ref. [1] for more details). This assumption is generally valid in high resolution confocal systems with the tube lens having typically two orders of magnitude longer focal length than the microscope objective lens. With the above assumptions, the dyadic Green's function of the radiating dipole along the detection pathway can be written as: where functions K n (ρ d , ϕ d , z d ) are given by: and k is the wavevector of emission in the detector medium with refractive index n . Finally, the electric field in the detector space originating from the radiating dipole in the focus of the aplanatic lens can be written as where (ρ d , ϕ d , z d ) are coordinates in the detector space. The above equation along with Eq. (19) allows us to eliminate p and express the electric field in the detector plane as a function of the excitation field E exc , excitation vs emission curve U and the detection pathway Green's If we displace the dipole from the origin (0, 0, 0) (the sample with the emitter is experimentally scanned), to a new position given by r s = (ρ s , ϕ s , z s ) (see Supplementary Figure 2 b), the electric field in the image space will transform as: where is the transverse magnification of the detection optical pathway and (ρ p , ϕ p , z p ) are coordinates centred on the pinhole. The integrated detected signal through the pinhole due to a radiating dipole in the position (ρ s , ϕ s , z s ) can then be expressed as: where R is the radius of the pinhole. This equation describes the final resolution of the imaging system as registered by the detector for any radius R of the pinhole placed in front of the detector. When the pinhole diameter is larger than the image of the dipole, where the size of the dipole image is usually expressed in terms of an Airy unit, the above can be simplified to where B ex (ρ s , ϕ s , z s ) denotes the intensity beam profile of the excitation beam.
In reality, the integrated detected signal S is equal to the emission intensity profile of the particle I em and this equation can be represented as Eq. 1 in the main text.
This approximation must be considered with extreme care for the case of nonlinear emitters as the size of the dipole image in the pinhole plane depends on the excitation vs emission curve U . As a result, the size of the dipole image changes depending on the excitation power. In the saturation regime of U , and especially for the dipole out-of-focus, the size of the dipole's image can become significant and the above approximation might no longer be valid. However, the above approximation is valid for the range of parameters studied in this manuscript (see Supplementary Note 2, Section F for more details).
Section D: Measurement of the excitation vs emission curve U m The process of experimentally measuring U m is as follows. We first find a spatially isolated particle as described in Supplementary Methods (STED imaging) and ensure that the nearby particles do not contribute to the detected particle intensity for the whole excitation range. This is especially important at higher excitation powers where the saturation of U m leads to large particle profiles which can easily lead to overlap of emission intensities of neighbouring particles. We centre the particle in x,y and z by scanning the particle in both the xy and yz planes. Once the particle is at the centre of the field of view, we illuminate the particle with the excitation laser for a period of 10 s and we measure the accumulated counts on the APD over this period. This is repeated for each excitation power of the U m profile. Whenever the system is not actively acquiring emission signal, we place a beam stop in the excitation path to avoid optically induced damage of the particle, which could result in fading of its emission intensity (Supplementary Note 1; Supplementary Methods (Pre-illumination procedure for UCNPs)). Extreme care must be taken, especially at higher excitation powers, where the emission fading can affect the curve U m of the particle. After measuring U m , we sweep through the excitation range to double-check that the emission fading was negligible during the course of U m measurement.
We use a λ/2-waveplate in tandem with the polarising beam splitter cube and a set of neutral density filters (Supplementary Methods (Setup description)) to control the excitation power. The secondary pathway of the polarizing beam splitter cube is used for monitoring the power oscillations introduced to the system by the excitation laser (the power oscillation is around 10% due to the imperfect polarization extinction ratio of the polarization maintaining fibre). The active monitoring of such power oscillations is essential due to high nonlinearity slope of material which substantially magnifies any input power oscillation effect and can easily bias the slope measurement of U m . The attenuation factor of each neutral density filter used was measured experimentally for the relevant wavelength (NE20A-B, Thorlabs and NE10A-B, Thorlabs have an attenuation factor of 20.9 and 9.9, respectively). We ensure that the range of excitation powers for each neutral density filter configuration has sufficient overlap with neighbouring neutral density filter ranges to guarantee that the data is not biased due to interchange of optical components. We measure the background counts of the collection system without a particle present in the field of view and subtract this background from U m .
Finally, we independently verify -by taking images of particles at several excitation powers -that the emission values of U m measured on the static particle are identical to the case of the particle being scanned through the beam. The numerical procedure for the extraction of the peak emission intensity from the raw particle profile images is described in detail in Supplementary Methods (Measuring the resolution of the system for a non-linear emitter).
Section E: Spline interpolation U of a measured excitation vs emission curve U m As Eq. (28) depends on the local beam intensity rather than the excitation power, the raw excitation vs emission curve U m has to be expressed in terms of the peak beam intensity I peak vs emission, where I peak = max (|E exc (ρ, ϕ, z)| 2 ). The total excitation power P exc , as measured in the back-focal plane of the objective, is related to the peak intensity I peak in the focus of the beam as where T = 0.4 is the transmittance of our microscope objective at the 976 nm wavelength and c is the result of the numerical integration which is found after substituting from Eq. (15). The numerical integration is performed in Matlab over sufficiently large integration limits to properly account for the beam tail contributions to the total value of the integral. We numerically fit a spline U through the curve U m (expressed as I peak vs emission) in a log-log space using the Matlab function fit with optional parameters 'smoothingspline', and with the 'smoothingParam' parameter set to 0.97. The spline interpolation U expresses the relation I em = U (I peak ) in a functional form. The spline U is needed to calculate the resolution of the system from Eq. (28). This is done in Matlab using the function feval.
Section F: Numerical determination of the resolution First, we use Eq. (15), derived in Supplementary Note 2, Section B to numerically calculate the excitation beam profile for λ 0 = 976 nm excitation wavelength. We approximate our Olympus UPLSAPO 100XO, NA = 1.4, 100x, objective by an aplanatic lens with an effective focal length of f = 2 mm and with the back-aperture diameter of d = 8 mm. The refractive index of the objective lens, immersion oil, coverslip and the mounting medium where the particles are dispersed is set to n = 1.518. The beam waist of the collimated, circularly polarized Gaussian beam illuminating the back-focal plane of the objective is measured to be w 0 = 6 mm. This means that the filling factor of the beam is f 0 = 1.5. The numerical integration of Eq. (15) is done in Matlab using the integrate function.
The numerically calculated beam shape is subsequently used as an input in the functional form of U , as derived in Supplementary Note 2, Section E to evaluate the resolution of the entire system by using Eq. (28). We note that the spline fit of U for values of emission with SNR < 1 can lead to unphysical oscillations of U that can lead to the unphysical oscillatory character of the simulated emission profile. For this reason, we restrict the emission profile/resolution simulations to scenarios where the peak power of the excitation beam produces a signal above the noise floor, specifically at least 3σ from the noise floor mean value (σ denotes the standard deviation of the noise floor). Such a choice ensures that the unphysical oscillation appears in the tail of the simulated emission profile and therefore does not affect the full-width-at-half-maximum (FWHM) value, respectively resolution.
The transverse magnification M of the detection pathway, as defined by Eq. (26) is M = 76 (f = 2 mm, f = 100 mm, n = 1.518, n = 1). Considering M , the pinhole diameter of the collection fibre (105 µm) is larger than the effective Airy unit of the detected beam profile in the pinhole plane across the whole range of experimental studied cases. The pinhole used in our system has therefore negligible effect on the simulated resolution.
We verified the negligible effect of the pinhole numerically by scanning the particle through the excitation beam and using Eq. (27). We calculated the resolution of the system for two cases -for a pinhole with a diameter of 105 µm and for a pinhole diameter corresponding to the case of the pinhole fully open (integration range in Eq. (27) was set to 1 mm). Even at the highest excitation powers used experimentally, which lead to large effective Airy units due to the saturation of U , the effect of the pinhole was negligible and basically within the numerical accuracy of the simulation. Similarly, we verified experimentally that there is a sufficient leeway in the pinhole alignment in all three directions such that no change in the detection PSF is observed. The numerical and experimental verification of the negligible effect of the pinhole conclusively demonstrates that the observed resolution improvement in SEE microscopy is purely due to material effects and not a combination of pinhole filtering and material properties.

Supplementary Note 3: Axial resolution vs excitation and vs emission
The axial resolution is simulated from the excitation-emission curves in Figure 4  Supplementary Note 4: Comparison of 4% and 8% Tm uSEE images The 8% Tm doped particles enable better resolution, both in lateral and axial direction, than the 4% Tm doped particles -panels a and b, respectively in Supplementary Figure 4. The lateral uSEE microscopy resolution is 184 nm for the highly doped UCNPs, and 230 nm for the lower doped nanoparticles. The axial resolution is 390 nm for the 8% Tm sample and 607 nm for the 4% Tm sample. Additionally, the brightness of the 8% Tm is higher, so the images can be obtained faster (8 ms/pxl), compared to the 4% Tm doped sample (20 ms/pxl) and with a better signal-to-noise ratio (12 for the highly doped sample and 5 for the low doped sample).
To facilitate visual comparison between the images in the paper and the Supplementary Information, which are sometimes taken with a different scanning range, we have added a flat background, equal to the average image background, around the experimental images in Supplementary Figure 4. The background is added in the area outside of the dashed black boxes in Supplementary Figure 4. A similar visual aid is used in Figure 1, Figure 2 and Figure 3 in the main text.

Characterization of UCNPs
A TEM picture of the UCNPs is presented in Supplementary Figure 5 a.The particles have a shape of hexagonal prisms and they are lying on different facets in the presented Figure. This allows us to determine the average diagonal of the hexagon 45 ± 2 nm and the average height of the prism 47 ± 2 nm. Spectra of a dense UCNP sample of NaYF 4 , doped with 20% Yb and 8% Tm, taken at different excitation power are plotted in Supplementary Figure 5 b. Note that increasing the power does not simply re-scale the spectrum, but also changes the relative intensity of the peaks due to the complicated balance between different optical processes in the particles -non-radiative energy transfer between the Yb 3+ and Tm 3+ ions, cross-relaxation processes between neighboring Tm 3+ ions, non-radiative relaxation and phonon-assisted energy transfer 3 . Note that the 455 nm peak which is of interest in this paper increases dramatically its relative intensity compared to the other peaks. This strong non-linearity allows achieving super-resolution in this excitation power region via uSEE microscopy. The spectra are qualitatively similar to the literature reports on similar upconversion particles 4,5 .

Setup description
We realize uSEE microscopy on a home-built conventional confocal setup (Supplementary Figure 6). Several additions have been introduced, which are not critical for SEE microscopy, but help to facilitate focusing/sample navigation (optical paths in Supplementary Figure 6 d, e) and imaging of conventional biologically targeted dyes, such as WGA-Alexa Fluor R 647 (optical path in Supplementary Figure 6 b). Additionally, optical path in Supplementary  Figure 6 b allows STED super-resolution mode for an auxiliary purpose of finding isolated particles. Two laser lines of 976 nm and 808 nm can be used to illuminate the sample. The system allows illumination with individual laser lines, but also simultaneously with the two co-aligned lasers.
976nm laser delivery arm (Supplementary Figure 6 a) The excitation source is a 976 nm laser diode (Thorlabs, BL976-PAG900), coupled into a polarization maintaining fibre (PM). The output from the PM fibre is collimated by an achromatic doublet lens L 1 (Thorlabs, AC254-035-B, f = 35 mm) to a waist of 6 mm. The polarization state of the beam is cleared up by passing the light through a linear polarizer LP (Thorlabs, LPVIS100-MP).
The linearly polarized light passes through a half-wave plate λ 1 /2 (Thorlabs, WPMH10M-980) and a polarising beam splitter (PBS) cube (Thorlabs, CCM1-PBS252/M). The PBS and the half-wave plate tandem enable control of the laser power by rotating the half-wave plate. The part of the beam reflected at the PBS is focused by an achromatic doublet lens L 2 (Thorlabs, AC254-060-B-ML, f = 60 mm) onto the power meter (Thorlabs S121C head with PM100D acquisition unit), which is used for active monitoring of the laser power on the excitation beam path.
The part of the beam transmitted through the PBS passes through a set of interchangeable neutral density filters (ND set ) (this enables access to a wider range of excitation powers) and the beam is subsequently redirected by two silver mirrors (Thorlabs, PF10-03-P01) towards the rest of the setup.
808nm laser delivery arm (Supplementary Figure 6 b) The 808 nm beam is delivered by a laser diode (Lumics, LU0808M250, 250 mW). The beam is collimated by an air-spaced achromatic doublet lens L 3 (Thorlabs, ACA254-060-B) to the same beam waist as the excitation beam (6 mm).
The beam then passes through a zero-order half-wave plate λ 2 /2 (Thorlabs, WPH10M-808) and a polarising beam splitter (PBS) cube (Thorlabs, CCM1-PBS252/M) to enable control of the optical power passing through the system.
If the STED mode is desired, the beam transmitted through the PBS then passes through a vortex plate VP (Holo/Or, VL-209-5-Y-A). The vortex plate, consisting of 16 discrete phase-steps, transforms the Gaussian beam into a Laguerre-Gaussian (LG) beam with a topological charge l = 1.
Afterwards, the laser beam is redirected by a set of two silver mirrors M (Thorlabs, PF10-03-P01) towards the rest of the setup. Figure 6 c) A dichroic mirror DM 1 (Chroma, ZT860lpxr) transmits the 976 nm laser and reflects the 808 nm laser to combine and re-direct them towards the rest of the setup. This allows using only one beam at a time, or both co-aligned beams for simultaneous sample illumination.

Combined delivery path 980nm and 808nm lasers (Supplementary
After the dicroic mirror, the laser beam passes through a position with a removable achromatic doublet lens (Thorlabs, AC254-0400-B-ML, f = 400 mm). The RL lens is placed in the beam path only when navigating within the sample. The removable lens (RL) is placed in such a way to generate a partially focused spot at the back-pupil of the microscope objective. This essentially generates a highly defocused beam in the sample plane illuminating a large field-of-view (FOV) improving the navigation within the sample. A slight adjustment to the RL lens positioning displaces the partially focused beam to the outer rim of the back-focal plane of the objective, which allows to image a relatively large FOV in a simplistic dark-field mode. While imaging in the SEE microsopy or STED microscopy mode, the RL lens is not present in the beam path.
After the RL position, the laser beam is reflected by a dichroic mirror DM 2 (Chroma, T750spxrxt) and passes through a removable pellicle RP (Thorlabs, BP145B2). The RP is in place only during sample navigation mode to re-direct the reflected laser light or luminescence towards the camera arms and when acquiring the point spread function (PSF) of the laser(s). Otherwise the RP is removed to maximize the amount of luminescent signal reaching the detector.
Next, the laser beam passes through a quarter-wave plate λ/4 (Thorlabs, AQWP10M-980), which is oriented in such a way to produce a circularly polarized beam. Afterwards, the beam is reflected by a silver mirror M (Thorlabs, PF10-03-P01) into the rigidly fixed microscope objective MO (Olympus, UPLSAPO 100XO, NA = 1.4, 100x) with a back pupil radius of 4 mm (this leads to a filling factor of f 0 = 1.5). Finally, the MO focuses the beam on the sample mounted on a piezo-electric 3-axis stage XYZ (Thorlabs, MAX311D/M with BPC203 piezo controller).
The light reflected, scattered or emitted from the sample is collected through the same objective MO. Then, it can be send to either one of the camera arms for focusing purposes and sample navigation or towards the APD for PSF measurement. Figure 6 d) We use the laser beams reflected or scattered from the sample for focusing and coarse alignment of the excitation and STED beams with respect to the sample. The faint reflections of the beams are visible whenever the beam encounters an interface with even a very small refractive index mismatch. This helps to focus on a proper position within the sample. Additionally, with the removable lens RL placed in position, we use this beam path to find cells or to find particles in the sample.

Camera sample navigation -laser reflection/scattering (Supplementary
The reflected beam propagates back through the microscope objective and is partially reflected on the removable pellicle RP. The removable mirror RM is not in place in this mode and the beam propagates to a fixed pellicle P1 (Thorlabs, BP145B2). Here, it is partially reflected towards the PSF collection arm and partially transmitted towards a second identical pellicle P2 (Thorlabs, BP145B2). At the second pellicle P2 part of the light is transmitted and used for practical alignment purposes (not shown here). The light reflected from P2 passes through a neutral density filter ND 1 (Thorlabs, NE20A-A). Finally, the reflection from the sample plane is imaged by a tube lens TL 2 (Thorlabs, AC254-200-B-ML, f = 200 mm) onto a camera CAM2 (CCD camera, Olympus, DP72). Figure 6 e) The coarse navigation within a luminescent sample is facilitated by imaging the luminescence onto a second camera -a low-noise camera CAM1 (Thorlabs, DCC3260M). The luminescence is collected through the same objective MO. The removable pellicle RP reflects part of the light reflected, scattered or emitted by the sample. To send the light towards CAM1, a removable silver mirror RM (Thorlabs, PF10-03-P01) is positioned after the pellicle. The luminescent signal is filtered from the excitation light by a shortpass filter F 1 (Thorlabs, FESH0750) and imaged onto the camera by a tube lens TL 1 (Thorlabs, AC254-200-A-ML, f = 200 mm). Figure 6 f ) To measure the PSF of the lasers, we use a standard procedure of imaging the scattered/reflected light by a sample with dispersed gold nanoparticles (size 80 nm) -Supplementary Methods (Measuring the excitation beam shape using a gold nanoparticle).

Laser PSF collection pathway (Supplementary
The reflected beam from the particle follows the same trajectory as above until it hits the fixed pellicle P1. At pellicle P1, the beam is reflected towards the PSF collection arm, where it is focused through an achromatic doublet lens L 4 (Thorlabs, AC254-100-B-ML, f = 100 mm). Afterwards, the beam is reflected at a dichroic mirror DM 3 (Chroma, T750spxrxt) and coupled into multimode fibre MMF (Thorlabs, M43L01, core diameter 105 µm).
During the laser PSF measurements, the interchangeable filter IF is not present. Finally, the light exiting the MMF is detected by a single photon avalanche photodiode APD (Excelitas, SPCM-AQRH-14-FC) which is synchronized with the piezo-electric 3-axis stage XYZ (Thorlabs, MAX311D/M with BPC203 piezo controller) to generate the beam PSF profiles.
If this mode is not used, this beam path is closed by a beam block immediately after pellicle P1.
Luminescence PSF collection pathway (Supplementary Figure 6 g) The luminescence is collected through the same objective MO. The removable pellicle RP is removed to maximize the collected luminescent signal. The signal is separated from the laser line by a dichroic mirror DM 2 (Chroma, T750spxrxt) and directed by a set of two silver mirrors M to an achromatic doublet lens L 5 (Thorlabs, AC254-100-A-ML, f = 100 mm). The signal then passes through a dichroic mirror DM 3 (Chroma, T750spxrxt) and a narrow fluorescence filter IF (Semrock FF01-448/20-25) to isolate the 455 nm emission line. The light is then coupled into a multimode fibre MMF (Thorlabs, M43L01, core diameter 105 µm) and detected by a single photon avalanche photodiode APD (Excelitas, SPCM-AQRH-14-FC) which is synchronized with the piezo-electric 3-axis stage XYZ (Thorlabs, MAX311D/M with BPC203 piezo controller) to generate the luminescent image. For the bio-measurements presented in the paper, the cells are stained with the cellular membrane marker WGAbiotin coupled to streptavidin tagged with Alexa Fluor R 647 and the nuclear intercalating dye DAPI. We visualized the Alexa Fluor R 647 with an identical procedure, where the excitation is done with the 808 nm laser and the detection is through an IF isolating the red emission from the dye (Chroma, ET660/60M-2P).
This pathway is also used to collect spectra of a dense sample. During the spectral acquisition, the fibre collecting the luminescent signal (PSF collection block) is connected to a spectrometer (Ocean Optics USB 2000) rather than an APD. During spectral acquisition, the beam is not scanned but remains static. When not used, this path is blocked after the dicroic mirror DM 2 by a beam block. Measuring the excitation beam shape using a gold nanoparticle

RP
We experimentally measure the excitation beam shape by scanning a gold particle with a diameter of 80 nm through the excitation beam and by detecting the scattered light intensity for each particle position. The scanning of the sample is realized using a commercially available software ImSpec which drives the stage and acquires the particle images. The field-of-view of scanning is set sufficiently large to account for the tails of the excitation beam shape. For all images presented in the paper, we sample above the Nyquist rate to avoid sampling artefacts. The detector pinhole/fibre diameter is sufficiently large to not influence the excitation beam shape by the detection pathway filtering (Supplementary Note 2, Section F). In our case, the pinhole size/collection fibre diameter is 105 µm and the transverse magnification M = 76. This leads to an expected Airy unit in the detection plane equal to 60 µm for λ = 976 nm wavelength. As a result, the detection pinhole can be considered fully open in the following (even for an out-of-focus gold particle) and the measured beam shape can be considered a valid excitation beam of our optical system.
We verify the above assumption by calculating the excitation beam shape using Eq. (15) and compare the simulation with the experimentally measured excitation beam shape profiles (Supplementary Figure 7 a,b). In Supplementary Figure 7 a, the simulated xz cross-section and the corresponding measurement is presented. Apart from a small tilt of the experimental beam, the two profiles are practically identical. In Supplementary Figure 7 b, x and z line cross-sections going through the centre of the beam in Supplementary Figure 7 a are presented, with calculated and measured FWHM essentially identical for the two cases.
The FWHM along the lateral (x) and the axial (z) direction is determined from experimental images by using the Matlab function improfile. For each direction x and z, we used improfile three times to extract pixel values along three adjacent lines. Next, the three improfile outputs were averaged to reduce the noise level. Finally, we used Matlab's function interp1 to increase the sampling rate of the averaged profile ten times to facilitate fitting. We fitted the resulting profile using Matlab's fit function with 'smoothingspline' parameter and with 'SmoothingParam' set to 0.02. Such fit of the experimental data was subsequently used for the precise determination of the FWHM value. Pre-illumination procedure for UCNPs As shown in Supplementary Note 1, the UCNP particles we are using show high photo-stability for hours, under relevant imaging conditions of uSEE microscopy and diffraction-limited confocal microscopy. However, to verify the validity of the developed theoretical framework over a broad region of excitation powers, covering different slopes of the excitation-emission curve, we sometimes needed to subject the particle to prolonged periods of intense illumination, outside of typical imaging conditions. Such illumination can lead to a drop of the emission intensity of the particle, which can significantly alter the excitation-emission curve and compromise the presented results.
To avoid such effects, we have employed a pre-illumination procedure. It is illustrated for the same pair of particles as in Supplementary Figure 1 and was applied to the particles after the measurements in Supplementary  Figure 1. The UCNP in the middle of the image was centered to the beam and was left there without scanning for a period of 2.5 hours. The excitation power was 2.2 mW at the back focal plane of the objective (11.8 mW µm −2 peak power density at the sample). The emission of the particle declined to about 40% of the initial emission value before stabilizing - Supplementary Figure 8. Evidently from the inset images, this procedure has a local effect and only influences the particle in the middle of the beam. We have used this procedure whenever we wanted to ensure a long-term emission stability under extreme excitation conditions over prolonged periods of time, for example we applied the pre-illumination step to the particle studied in Figure 1 and Figure 3 in the main text. Under typical experimentally relevant conditions, used for obtaining the rest of the data in the paper, this procedure was not needed. The pre-illumination procedure reduces the emission intensity of the particle, but allows reaching a stable emission level after illumination at 11.8 mW µm −2 peak power density for 2.5 hours.
Measuring the resolution of the system for a non-linear emitter Experimentally, the resolution of the system for the case of a luminescent particle is obtained by scanning the particle through the excitation beam, acquiring the corresponding image of the particle and measuring the FWHM of the emission profile. Similar to Supplementary Methods (Measuring the excitation beam shape using a gold nanoparticle), the resulting image is used to determine the FWHM along the lateral (x) and the axial (z) direction by using the Matlab function improfile averaged over three lines going through the centre of the imaged particle. We again increased the sampling rate of the averaged line profile ten times by using linear interpolation between the original profile points by using interp1 function in Matlab. We fit the resulting profile averaged over three lines using Matlab's fit function with 'smoothingspline' parameter and with 'SmoothingParam' set to 0.02. This leads to a sufficiently smooth fit even for highly noisy images at low excitation powers, with very good estimators of the profile's peak and half intensity points -the crucial parameters for determining the FWHM of the non-linear particle. The background level is subtracted from the images to remove bias to the FWHM evaluation.

STED imaging
The stimulated emission depletion (STED) super-resolution modality of the setup is not necessary to achieve uSEE imaging. In our studies we use the STED capability for navigation and verification purposes. STED allows imaging of the UCNPs used in this paper (NaYF 4 , doped with 20% Yb and 8% Tm) with resolution close to the size of the particles, which helps identifying isolated UCNPs - Supplementary Figure 9 a. Additionally, the 808 nm laser line which is used for depletion in the STED mode can be used for excitation of the WGA-Alexa Fluor R 647 dye and help confirm the position of the UCNPs with respect to the cells.