Background-free 3D nanometric localisation and sub-nm asymmetry detection of single plasmonic nanoparticles by four-wave mixing interferometry with optical vortices

Single nanoparticle tracking using optical microscopy is a powerful technique with many applications in biology, chemistry and material sciences. Despite significant advances, localising objects with nanometric position accuracy in a scattering environment remains challenging. Applied methods to achieve contrast are dominantly fluorescence based, with fundamental limits in the emitted photon fluxes arising from the excited-state lifetime as well as photobleaching. Furthermore, every localisation method reported to date requires signal acquisition from multiple spatial points, with consequent speed limitations. Here, we show a new four-wave mixing interferometry technique, whereby the position of a single non-fluorescing gold nanoparticle is determined with better than 20 nm accuracy in plane and 1 nm axially from rapid single-point acquisition measurements by exploiting optical vortices. The technique is also uniquely sensitive to particle asymmetries of only 0.5% aspect ratio, corresponding to a single atomic layer of gold, as well as particle orientation, and the detection is background-free even inside biological cells. This method opens new ways of of unraveling single-particle trafficking within complex 3D architectures.

The ability to localize optically and eventually track the position of objects at the nanoscale requires ways to overcome the Abbe diffraction limit given by λ/(2NA), where λ is the wavelength of light and NA is the numerical aperture of the imaging lens. For visible light and high NA objectives this limit is roughly 200 nm. Current techniques for overcoming this limit can be divided into near-field and far-field optical methods. In near-field optics, super-resolution is achieved by localizing the light field (incident and/or detected) using optical fiber probes with sub-wavelength apertures [1], metal coated tips [2], or plasmonic nano-antennas [3]. These techniques can provide spatial resolutions down to about 10 nm, however are limited to interrogating structures accessible to optical tips and/or patterned substrates and require precise position control of the probe.
Far-field methods overcome this drawback and allow to localise single biomolecules inside cells. However, to achieve sufficient contrast and specificity against backgrounds, superresolution far-field field methods are dominantly fluorescence-based. They exploit either the principle of single emitter localisation [4,5] or optical point-spread function (PSF) engineering [6] and the fluorophore non-linear response in absorption or stimulated emission. Beside specific advantages and disadvantages of these two approaches, their implementation using fluorescence results in significant limitations. Fluorophores are single quantum emitters and thus only capable of emitting a certain maximum number of photons per unit time due to the finite duration of their excited-state lifetime, moreover they are prone to photobleaching and associated photo-toxicity.
Alternatively to using fluorescent emitters, far-field nanoscopy techniques have been shown with single metallic nanoparticles (NPs) optically detected owing to their strong scattering and absorbtion at the localised surface plasmon resonance (LSPR). These detection methods are photostable, and the achievable photon fluxes are governed by the incident photon fluxes and the NP optical extinction cross-section. Position localisation below diffraction similar to single emitter localisation can be achieved via wide-field techniques such as bright-field microscopy [7], dark-field microscopy [8], differential interference contrast [9], and interferometric scattering microscopy [10]. These techniques, however, are not background-free and either use large NPs (diameters ≥ 40 nm) to distinguish them against endogenous scattering and phase contrast in heterogeneous samples or work in optically clear environments. A more selective technique uses photothermal imaging where the image contrast originates from the refractive index change in the region surrounding the nanoparticle due to local heating following light absorption. This is a focused-beam scanning technique and has been used to track single 5 nm diameter gold NPs in two dimensions (however not 3D) [11]. By detecting the nanoparticle only indirectly via the photothermal index change generated in its surrounding, this method is also not free from backgrounds.
In fact, photothermal contrast has been shown in the absence of NPs due to endogeneous absorption in cells [12]. Moreover, as for every localisation method, it required signal acquisition from multiple spatial points to determine the position of the single NP, with consequent limitation in tracking speed.
Four-wave mixing (FWM), triply-resonant to the LSPR, was shown by us to be a very selective, high-contrast photostable method to detect single small gold NPs [13,14]. It is a third-order nonlinearity which originates from the change in the NP dielectric constant induced by the resonant absorption of a pump pulse and subsequent formation of a non-equilibrium hot electron gas in the metal [14]. It is therefore very specific to metallic NPs which are imaged background-free even in highly scattering and fluorescing environments [13]. In this work, we show theoretically and experimentally a new FWM detection modality which enables to determine the position of a single gold NP with < 20 nm accuracy in plane and < 1 nm axially from scanless single-point background-free acquisition in the 1 ms time scale, by exploiting optical vortices of tightly focussed light.

I. FOUR-WAVE MIXING INTERFEROMETRY TECHNIQUE
A sketch of the FWM technique is shown in Fig. 1a. The key developments in this new design compared to previous works are the epi-collection (reflection) geometry and the dualpolarisation heterodyne detection scheme. We use a train of femtosecond optical pulses with repetition rate ν L which is split into three beams, all having the same center wavelength, resulting in a triply degenerate FWM scheme. One beam acts as pump and excites the NP at the LSPR with an intensity which is temporally modulated with close to unity contrast by an acousto-optic modulator (AOM 1 ) driven at carrier frequency ν 1 with a square wave amplitude modulation of frequency ν m . The change in the NP optical properties induced by this excitation is resonantly probed by a second pulse at an adjustable delay time τ after the pump pulse. Pump and probe pulses of fields E 1 and E 2 , respectively, are recombined into the same spatial mode and focused onto the sample by a high NA microscope objective. The sample can be positioned and moved with respect to the focal volume of the objective by scanning an xyz sample stage with nanometric position accuracy. A FWM field E FWM (proportional to E 1 E * 1 E 2 ) is collected together with the probe in reflection (epi-direction) by the same objective, transmitted by the beam splitter (BS 1 ) used to couple the incident beams onto the microscope, and recombined in a second beam splitter (BS 2 ) with a reference pulse field (E R ) of adjustable delay. The resulting interference is detected by two pairs of balanced photodiodes. A heterodyne scheme discriminates the FWM field from pump and probe pulses and detects amplitude and phase of the field. In this scheme, the probe optical frequency is slightly upshifted via a second AOM (AOM 2 ), driven with a constant amplitude at a radiofrequency of ν 2 , and the interference of the FWM with the unshifted reference field is detected. As a result of the amplitude modulation of the pump at ν m and the frequency shift of the probe at ν 2 , this interference gives rise to a beat note at ν 2 with two sidebands at ν 2 ± ν m , and replica separated by the repetition rate ν L of the pulse train. A multichannel lock-in amplifier enables the simultaneous detection of the carrier at ν 2 − ν L and the sidebands at ν 2 ± ν m − ν L . Via the in-phase (Re) and in-quadrature (Im) components for each detected frequency, amplitude and phase of the probe field reflected by the sample (E 2r ) and of the epi-detected FWM field are measured (see sketch in Fig. 1a).
A key point of the technique is the use of a dual-polarisation balanced detection. Firstly, probe and pump beams, linearly polarised horizontally (H) and vertically (V) respectively in the laboratory system, are transformed into cross-circularly polarized beams at the sample by a combination of λ/4 and λ/2 waveplates. The reflected probe and FWM fields collected by the same microscope objective travel backwards through the same waveplates, such that the probe reflected by a planar surface returns V polarized in the laboratory system. The reference beam is polarised at 45 degree (using a polariser) prior recombining with the epi-detected signal via the non-polarizing beamsplitter BS 2 . A Wollaston prism vertically separates H and V polarizations for each arm of the interferometer after BS 2 . Two pairs of balanced photodiodes then provide polarization resolved detection, the bottom (top) pair detecting the current difference (for common-mode noise rejection) of the V (H) polarised interferometer arms. In turn, this corresponds to detecting the co-and cross-circularly polarised components of E 2r and E FWM relative to the incident circularly polarized probe, having amplitudes (phases) indicated as A ± 2r and A ± FWM (Φ ± 2r and Φ ± FWM ) in the sketch in Fig. 1a where + (−) refers to the co (cross) polarised component.

VORTICES
To elucidate conceptually how the nanometric position accuracy arises from this dualpolarisation resolved FWM interferometry detection scheme, we simulated numerically the field distribution in the focal region of a 1.45 NA objective. The simulation parameters (wavelength, coverslip thickness, medium refractive index, back objective filling factor) were chosen to match the actual experimental conditions (see Methods). The amplitude and phase components of the field in the focal plane are shown in Fig. 1b for a left circularly (σ + ) polarised input field. Due to the high NA of the objective and the vectorial nature of the field, there is a significant cross-circularly polarised component that forms an optical vortex of topological charge l = 2, i.e. it has an amplitude (A − ) which is zero in the focus center and radially-symmetric non-zero away from the center, and a phase (Φ − ) changing with twice the in-plane polar angle. A point-like gold NP displaced from the focus center experiences this field distribution, and will in turn emit a field with an amplitude and a phase directly related to the NP radial and angular position. While this is the case for both reflected and FWM fields, only the FWM signal is background-free hence suited for tracking the NP in heterogeneous environments. In the following, we discuss the FWM field emitted from such NP and the corresponding coordinate retrieval.
We calculate the FWM field starting from the polarization of the NP at position r induced by the probe field, described by p(r) = 0 mα E 2 (r) where 0 is the vacuum permittivity, m is the dielectric constant of the medium surrounding the NP (glass and silicon oil in the experiment) andα is the particle polarizability tensor. We calculateα in the dipole approximation, valid for particle sizes much smaller than the wavelength of light (Rayleigh regime). To take into account non-sphericity of real particles, which will be relevant in the experiments as shown later, we adopted the model of a metallic ellipsoid with three orthogonal semi-axes of symmetry a, b and c. In the particle reference system,α is diagonal and its eigenvalues are given by [15] where is the dielectric permittivity of the particle, and L i with i = a, b, c are dimensionless quantities defined by the particle geometry (see Supplementary Information section S1.i). For an arbitrary particle orientation in the laboratory system, the polarizability tensor can be transformed usingα =M −1α M whereM is a rotation matrix.  ponents of the reflected probe field and FWM field A ± 2r and A ± FWM (Φ ± 2r and Φ ± FWM ) respectively, as a function of particle position in the sample focal plane (x, y), and in a section along the axial direction (x, z) through the focus, where + refers to the co-polarised component and -to the cross-polarised component relative to the left-circularly polarised incident probe. The calculation assumes a perfectly spherical gold NP in the dipole approximation. The inset shows a sketch of how the amplitude and phase of the FWM field ratio and the phase of the co-polarised FWM field can be used to locate in 3D the spatial position of the NP relative to the focus center. Linear grey scale from −π to π for all phases, and from m to M for field amplitudes, as indicated. The amplitude ratio of reflected probe and FWM are shown on a logarithmic greyscale over 5 orders of magnitudes.
The interference of the reflected probe field with the reference field is calculated using E ± 2r = (E ± R ) * ·α E + 2 whereα E + 2 is the particle induced polarisation for a σ + polarised incident probe field (we dropped the constant 0 m in p for brevity), and E ± R are reference fields equal to left (+) and right (-) circularly polarised input fields. The technique is configured such that the optical mode of probe and reference fields are matched, hence E R was calculated as the field distribution in the focal region in the same way as E 2 (see Methods), and backpropagated via time reversal. Similarly, the FWM interference is calculated as Here, δα is the pump-induced change of the particle polarizability which we have modeled as described in our previous work [14]. Briefly, δα arises from the transient change of the electron and lattice temperature following the absorption of the pump pulse by the NP. δα depends on the pump fluence at the NP, on the particle absorption cross-section, and on the delay time between pump and probe pulses (see also Supplementary Information Fig. S2 and S10). The simulations in Fig. 2 were performed to reproduce the experimentally measured FWM signal strength on a 30 nm radius gold NP at τ = 0.5 ps shown later (note that τ ∼ 0.5 ps is the delay for which the FWM amplitude reaches its maximum as a result of the ultrafast heating of the electron gas [14]).

E ±
2r and E ± FWM have amplitudes and phases as a function of particle position in the focal region as shown in Fig. 2 for the case of a perfectly spherical NP. Similar to the spatial distributions of the focussed field shown in Fig. 1b, E − 2r and E − FWM form optical vortices of l = 2 topological charge. The ratio E − 2r /E + 2r is also shown with its amplitude A − 2r /A + 2r and phase Φ − 2r − Φ + 2r . The FWM field distribution has a narrower PSF than the reflection, as expected from the third-order nonlinearity. The phase of E FWM is shifted compared to E 2r due to the phase difference between δα andα. Slices along the axial direction (z) are shown for E FWM (a similar behavior is observed for E 2r , see Supplementary Information Fig. S1).
Of specific importance for the localization of the NP in 3D are the amplitude and phase of the FWM field ratio E − FWM /E + FWM and the phase of the co-polarized E + FWM as highlighted by the red panels and sketch in Fig. 2.
The co-polarized FWM amplitude is much stronger (∼ 90 fold) than the cross-polarized FWM, and has a phase which is independent of the lateral position over the PSF width.
Conversely, along the z-axis this phase is linear in the displacement between particle and center of the focus and can be used to determine the particle z coordinate. This can be easily understood as due to the optical path length difference between the particle and the observation point. For a plane wave of wavevector k = 2πn/λ with n refractive index in the medium the phase would be 2kz, the factor of 2 accounting for the back and forth path in reflection geometry reflection. We find a linear relationship between Φ + FWM and z (see Supplementary Information Fig. S3) with a slope ∂z/∂Φ = 38.8 nm/rad, slightly larger than λ/(4πn) = 28.8 nm/rad. This is due to the propagation of a focussed beam with high NA where a Gouy phase shift occurs, reducing the wavevector in axial direction due to the wavevector spread in lateral direction.
For the in-plane radial coordinate R of the NP position relative to the focus position (see sketch in Fig. 2), we find that the FWM amplitude ratio A − FWM /A + FWM scales quadratically with R up to R ∼ 60 nm, such that this coordinate can be calculated as Notably, by using the FWM ratio the retrieved R is independent of the pump, probe and reference powers, and of the NP size (in the dipole approximation). Conversely, R 0 is specific to the NA of the microscope objective used and the probe beam fill factor (see Supplementary Information Fig. S4) and decreases with increasing NA, showing that high NAs are required to localize the NP in R down to small distances. Finally, the angular position coordinate ϕ can be taken from the phase of the FWM ratio Supplementary Information Fig. S3). In essence, using these relationships we can locate the NP in 3D (z, R and ϕ) by scanless polarisation-resolved and phase-resolved FWM acquisition at a single spatial point.

III. BACKGROUND-FREE FOUR-WAVE MIXING DETECTION
Prior to quantifying experimentally the nanometric localisation accuracy, it is important to emphasize that our FWM detection is background-free even in scattering and/or autofluorescing environments, making it applicable to imaging single small NPs inside cells, surpassing other methods reported in the literature. Localizing single plasmonic NPs with nanometric accuracy even at sub-millisecond exposure times can be achieved in an optically clear environment with simpler techniques such as dark-field microscopy [8] and interferometric scattering microscopy [10]. However, since these techniques use the linear response of the NP they are substantially affected by endogenous scattering and fluorescence, which severely limit their practical applicability in heterogeneous biological environments.
To exemplify this point, Fig. 3(top) shows fixed HeLa cells that have internalized gold NPs plane image of a 25 nm radius gold NP in 5% agarose gel, via the co-circularly polarised reflection amplitude A + 2r (e), cross-circularly polarised reflection amplitude A − 2r (f) and phase Φ − 2r (g), cocircularly polarised FWM amplitude A + FWM (h) and phase Φ + FWM (i). FWM was acquired with 0.5 ms pixel dwell time, pixel size in plane of 38 nm. Grey scales are linear from −π to π for all phases, and from m to M for field amplitudes, as indicated.
of 20 nm radius, imaged with FWM using a 1.45 NA oil-immersion objective. High-resolution DIC microscopy was available in the same instrument (for details see Supplementary Information section S2.iii and S2.vii). Fig. 3a shows the DIC image of a group of HeLa cells on which reflection and FWM imaging was performed in the region highlighted by the dashed frame. The co-circularly polarized reflection image A + 2r shown in Fig. 3b correlates with the cell contour seen in DIC, and shows a spatially varying contrast due to thickness and refractive index inhomogeneities in the sample. Even with a particle diameter as large as 40 nm, gold NPs are not distinguished from the cellular contrast neither in DIC nor in the A + 2r reflection image. Detecting the cross-polarized reflection A − 2r (Fig. 3c) which has been suggested as a way to improve contrast [16] is still severely affected by the cellular scattering background. On the contrary, the co-circularly polarised FWM amplitude A + FWM shown in  With the reflection amplitude scaling as the volume of the particle, a 15 nm radius NP would be indistinguishable from the background in the present case. Furthermore, the interference with the background limits the accuracy of particle localization. Notably, the phase Φ − 2r is also severely affected by the scattering from the gel. Conversely, the FWM amplitude A + FWM and phase Φ + FWM shown in Fig. 3h,i are background-free (as can be seen by the random phase outside the particle) and clearly resolve the NP despite the heterogeneous surrounding.

IV. EXPERIMENTAL LOCALISATION: THE ROLE OF NANOPARTICLE EL-LIPTICITY
To experimentally quantify the nanometric localisation and its accuracy due to photon shot-noise, we started by examining a single nominally spherical 30 nm radius gold NP drop cast onto a glass coverslip and immersed in silicon oil, using an index-matched 1.45 NA oil-immersion objective.
The experimental data shown in Fig. 4 are scans of the (x, y, z) sample position. They reveal all main features seen in the calculations, namely a ring-like spatial distribution inplane of the amplitude ratio A − FWM /A + FWM and a phase of the cross-polarized component, and in turn of the FWM ratio, twice rotating from 0 to 2π along the in-plane polar angle.
A main difference to the calculations is the observation of two displaced nodes, rather than a single central node, in the amplitude of the cross-polarized component, and in turn a nearly constant phase of the cross-polarized component in the central area. The FWM ratio resolved across different axial planes is also shown in Fig. 4. We observe that, besides these minima, the axially resolved distributions agree with the calculated behavior of a linear relationship between Φ + FWM and z and a z-dependent Θ 0 manifesting as a rotated phase Notably, by measuring on > 50 particles in the sample we could not find any pattern forming a l = 2 optical vortex as simulated in Fig. 2. Different particles showed minima with different relative displacement, or single displaced minima, or no minima within the spatial range of sufficient signal to noise (see Supplementary Information Fig. S12), suggesting that these minima are related to physical differences between particles. Indeed, we can explain these experimental findings assuming a small particle asphericity. This is shown in the calculations in Fig. 4 where we used an ellipsoid nanoparticle with semi-axes a = 30.135 nm and b = c = 30 nm along the x, y and z axis in the particle reference system, which is rotated in plane relative to the laboratory system (see also see Supplementary Information Fig. S5 and S6). It is remarkable that an asymmetry of only 0.5% in aspect ratio, or about one atomic layer of gold, manifests as a significant perturbation of the cross-polarized field patterns compared to the spherical case. With such sensitivity to asymmetry, the lack of experimentally observed pattern corresponding to a perfectly spherical particle is not surprising, considering the real particle morphology as observed in transmission electron microscopy (TEM) for this sample (see top inset in Fig. 4 and Supplementary Information   Fig. S8).
The shot-noise in the experiment as well as the deviation from perfect sphericity affect the localisation accuracy. This is analyzed in Fig. 5. Firstly, we performed simulations as in (in both the in-phase Re and in-quadrature Im components) in a spatial region away from the particle, where no FWM is detectable. The standard deviation σ of this distribution was deduced, and was found to be identical in both components, and for the co-polarised and cross-polarised components, as expected for an experimental noise dominated by the shot-noise in the reference beam (see Supplementary Information Fig. S9 for the dependence of σ on the power in the reference beam). A relative noise figure was defined as σ/A 0 with A 0 being the maximum measured value of the co-polarized FWM field amplitude.
The simulations were performed using a statistical distribution of the FWM field values (Re and Im components) at each spatial pixel having the same relative noise σ/A 0 as the experimental data.
To quantify the uncertainty in the localisation of the particle coordinates, we then cal- observe that adding the experimental noise results in a localization uncertainty of less than 20 nm. Increasing the relative noise by one order of magnitude to 0.7% results in a localization uncertainty of ∼ 50 nm, as shown in Fig. 5c. We emphasize that ∆ is dominated by the In fact, considering that the axial direction z is determined directly by the slope of Φ + FWM and does not involve the cross-polarized FWM, we find z localization accuracies as small as 0.3 nm for 0.7% noise in Φ + FWM . Note that the relative shot-noise in the experiments scales as 1/ V √ tI 1 √ I 2 where V is the particle volume (in the Rayleigh regime), I 1 (I 2 ) is the intensity of the pump (probe) beam at the sample, and t is the acquisition time. It can therefore be adjusted according to the requirements. For example, a 15 nm radius NP imaged under the experimental conditions as in Fig. 4 would give rise to an 8-fold increase in relative noise, hence still maintaining the localization uncertainty to below 50 nm. These results elucidate that a localization accuracy of below 20 nm in-plane and below 1 nm in axial direction is achievable with the proposed method with a realistic shot-noise level as in the experiment.
However, the lack of particle symmetry is a limitation for the in-plane localization. Introducing the particle asymmetry, without shot-noise, as calculated in Fig. 4, results in a significant deviation ∆ ∼ 100 nm in the central area, due to the lack of a central node in the cross-polarized FWM amplitude (see Fig. 5d). On the other hand, it is remarkable how sensitive the described FWM technique is to particle asymmetry, which can be used as a new tool to detect particle ellipticity down to a/b − 1 = 10 −4 (corresponding to atomic accuracy comparable to TEM) as well as particle orientation. We find that the FWM amplitude ratio in the focus center scales linearly with the particle ellipticity and that the phase of the FWM ratio in the focus center scales with the in-plane particle orientation angle (which for the data in Fig.4 was found to be γ = 150 • , see Supplementary Information Fig. S5 and S6). We also note that when examining gold NPs of 5 nm radius (see Supplementary Information   Fig. S8 and S13), we found a large proportion (∼ 70%) having a FWM amplitude ratio 02 in the center of the focal plane, in this case limited by the signalto-noise ratio rather than the particle asymmetry, suggesting that smaller NPs might be intrinsically more mono-crystalline hence symmetric.
Importantly, when NPs are not immobilized onto a surface but freely rotating (a more relevant scenario for particle tracking applications) we can expect that the rotational averaging over the acquisition time (the rotational diffusion constant in water at room temperature of a 30 nm radius nanoparticle is ∼ 10 4 rad 2 /s) would result in an effective symmetry, and would lead to a localization accuracy only limited by shot-noise. This case is discussed in the next section.

V. ROTATIONAL AVERAGING
To mimic a relevant biological environment such as the cytosolic network, as well as having single NPs freely rotating but not diffusing out of focus, gold NPs of 25 nm radius were embedded in a dense (5% w/v) agarose gel in water (see Supplementary Information section S2.ii) and measured using an index-matched 1.27 NA water-immersion objective. For these experiments the focused size of the pump beam was increased by a factor of two (by under-filling the back focal aperture), to enlarge the region where a single NP could diffuse while still being excited by the pump field hence giving rise to FWM (this also increases the maximum cross-polarized FWM amplitude). Conversely, the probe beam was tightly focused to exploit the full NA of the objective an in turn exhibit an l = 2 optical vortex in the cross-circularly polarised component as calculated in Fig. 2. Supplementary Information Section S1.vii). Conversely, it is not possible to reproduce the experimental findings by assuming a non-rotating randomly-oriented asymmetric particle (see Supplementary Information Section S2.xii), hence free rotation is key to the pattern observed in Fig. 6a. Since the only particle asymmetry that does not average upon rotation is 3D chirality, we suggest that this response is a manifestation of chirality (albeit a detailed theoretical understanding of the nonlinear optical response of a chiral particle is beyond the scope of this work). We remind that these quasi-spherical NPs have irregularities of few atoms clusters which can lead to symmetry breaking, as already seen for their ellipticity, and thus also 3D chirality (see TEM in Supplementary Information Fig. S8). Notably, an l = 1 optical vortex offers a simpler dependence for the retrieval of the NP position coordinates, with a FWM amplitude ratio A − FWM /A + FWM scaling linearly with R and a phase Θ = Φ − FWM − Φ + FWM directly changing with the polar angle ϕ (see Supplementary Information Fig. S14). Using these dependencies (alongside the axial position z given by Φ + FWM as previously shown), we have retrieved the NP position coordinates in 3D for the scan shown in Fig. 6a and compared them with the position coordinates recorded from the scanning piezoelectric sample stage. This is shown in Fig. 6b. We observe deviations below 50 nm in x, y between the retrieved NP positions and the coordinates from the stage and a remarkable accuracy in z better than 10 nm, limited by systematic drifts rather than shot-noise. From the shot noise σ/A 0 = 0.007 in these measurements and the coordinate retrieval parameters, we deduce the uncertainty to be below 20 nm in x, y and 1 nm in z, using only 1 ms acquisition time.

VI. SINGLE PARTICLE TRACKING
A demonstration of the practical applicability of the method for scanless single particle tracking in 3D is shown in Fig. 7 on gold NPs of 25 nm radius embedded in a dense (5% w/v) agarose gel in water, which provides a heterogeneous environment as shown in Fig. 3(bottom). For practical purposes, we implemented a simple in situ calibration of the in-plane NP position coordinates without the need for prior knowledge and/or characterization of the particle optical response. As shown in Fig. 7a, we apply a small oscillation of known amplitude (16 nm) and frequency (25 Hz) to one axis of the sample stage. When the FWM field ratio from the NP encodes the NP in-plane position as discussed in the previous section, this oscillation is detected in the measured E − FWM /E + FWM with an amplitude and a phase which we can use to directly calibrate the FWM field ratio in terms of in-plane position coordinates of known size and direction (see also Supplementary Information section S2 x.iii). Fig. 7a shows the power spectrum of the Fourier transform of E − FWM /E + FWM which exhibits a peak at this modulation frequency. The correspondingly calibrated particle position coordinates over time are shown Fig. 7b. Conversely, if for a NP we do not observe this oscillation in E − FWM /E + FWM , the FWM field ratio is a measure of the NP in-plane asymmetry and orientation. This case is shown in Fig. 7c-d (for the power spectrum see Supplementary   Information section S2 x.iii). Importantly, we can still accurately measure the axial position coordinate, which is directly given by Φ + FWM and does not involve the cross-polarized FWM, while now E − FWM /E + FWM encodes information on the NP asymmetry and orientation, rather than its in-plane position. Fig. 7c shows examples of traces of the NP z position acquired over several tens of seconds (with 2 ms point acquisition) from which we can see a 'jumping behavior' around two preferred axial locations, suggesting the presence of two pockets in the agar gel separated by about 100 nm. Fig. 7d shows the z position together with E − FWM /E + FWM (as real and imaginary parts) as a zoom over a time window (from the top trace in (c)) during which the NP got stuck in a corner, below the center of the lower pocket.
From the strong and slowly varying FWM field ratio, the corresponding time evolution of the NP in-plane orientation angle and aspect ratio is given in Fig. 7f. Finally, Fig. 7e  Shape asymmetry of the NP of as little as 0.5% ellipticity, corresponding to about one atomic layer of gold, also induces a cross-circularly polarized field. This exceptional sensitivity to asymmetry eventually limits accurate position retrieval in-plane. This can be overcome by observing a NP freely rotating, such that rotational averaging of the asymmetry occurs over the acquisition time.
Experimentally, we show that the FWM detection is completely background-free in scattering environments such as biological cells, and outperforms existing methods such as reflectometry, scattering and differential interference contrast. A localization accuracy of below 50 nm in-plane, and below 10 nm axial, limited by systematic drifts, is shown with a single NP of 25 nm radius freely rotating in an agarose gel, at an acquisition time of only 1 ms. The shot-noise limited accuracy is found to be below 20 nm in plane and 1 nm axially.
Smaller NPs can be measured by correspondingly increasing the intensity of the incident beams and/or the acquisition time. Notably, we find that, during free rotation, the crosscircularly polarized FWM field distribution in the focal plane is an l = 1 optical vortex, as opposed to the predicted l = 2 for a perfectly spherical particle, and we suggest that this is a manifestation of NP chirality not averaged during rotation.
To demonstrate the practical applicability of the method for single particle tracking, we provide two examples of NPs diffusing in a dense agarose gel network. In one case, we show particle position coordinates retrieved using an in situ calibration procedure and the corresponding tracking in 3D. In the second case, we show that when the cross-circularly polarized component is dominated by the effect of shape asymmetry, it can be used for tracking the in-plane asymmetry and orientation, while the axial position coordinate provides information on the particle diffusion.
Ultimately, this method paves the way towards a new form of single particle tracking, where not only the NP position, but also its asymmetry, orientation, and chirality are detected with sub-millisecond time resolution, revealing much more information about the NP and its complex dynamics (e.g. hindered rotation) while moving and interacting within a disordered 3D environment. For the calculations in Fig. 6 the simulation parameters were NA=1.27, n = 1.333, β = 2.15 for the pump and β = 0.97 for the probe beam. The calculation was performed for a linear polarization of the incident field along the x axis. This results in a vectorial field at the focus called E x (r). For the orthogonal linear incident polarization, the results were rotated counterclockwise to obtain E y . To simulate circular incident polarization, the calculated field maps E x and E y were combined with complex coefficients, namely:    Background-free 3D nanometric localisation and sub-nm asymmetry detection of single plasmonic nanoparticles by four-wave mixing interferometry with optical vortices -Supplementary Information S1. CALCULATIONS

i. Polarisability of a non-spherical gold nanoparticle
We describe a non-spherical nanoparticle (NP) as a metallic ellipsoid with three orthogonal semi-axes of symmetry a, b and c. In the particle reference system the polarisability tensorα is diagonal, and its eigenvalues are given by [15] α where is the dielectric constant of the NP, m is the dielectric constant of the surrounding medium, and L i with i = a, b, c are dimensionless quantities defined by the NP geometry as Only two out of three geometrical factors are independent, since for any ellipsoid L a + L b + L c = 1. For a spherical particle a = b = c, and For prolate spheroids (i.e. cigar-shaped) with a > b = c or oblate (i.e. pancake-shaped) spheroids with a = b > c the geometrical factors can be expressed analytically [15]. In the most general case, we calculated the integrals in Eq. S2 using the numerical integration function integral of Matlab.
ii. Reflected probe field Fig. S1 shows the reflected probe field as function of the NP position, calculated as in Fig. 2 of the main manuscript, including its axial dependence through the focus. iii. Pump-induced change of particle polarisability: Spatial distribution To calculate the pump-induced change of the particle polarizability δα we used the model developed in our previous work [14].

iv. Parameters for NP position localisation
From the calculations in Fig. 2 of the main manuscript we find a linear relationship between Φ + FWM and z as shown in Fig. S3 with a slope ∂z/∂Φ = 38.8 nm/rad, slightly larger than λ/(4πn) = 28.8 nm/rad. This is due to the propagation of a focussed beam with high NA where a Gouy phase shift occurs, effectively increasing the wavelength in axial direction. Note that the measured Φ + FWM has an offset, due to the phase of the reference field. In the shown representation, the value of Φ + FWM at the focal plane was taken as offset and subtracted, such that Fig. S3 directly gives ∂z/∂Φ. v. Dependence of R 0 on objective NA As described in the previous section and in the manuscript, the four-wave mixing amplitude ratio A − FWM /A + FWM scales quadratically with the radial coordinates R (for small values of R), such that in first approximation this coordinate can be calculated as R = Fig. S4 shows the dependence of the parameter R 0 in the focal plane (z = 0) on the objective numerical aperture using a filling factor β = 0.83 as in the experiment, and in a medium of n = 1.5185. We find that R 0 scales quadratically with 1/NA. The width w 0 for the in-plane distribution of A + FWM fitted as a Gaussian function using the curvature near the center, i.e. A + FWM = A 0 e −(R/w 0 ) 2 , is also shown, representing an effective point-spread function which scales linearly with 1/NA. vi. Amplitude and phase of FWM ratio versus particle ellipticity Although being a limitation to achieve nanometric position accuracy, it is remarkable how sensitive the described FWM technique is to particle asymmetry, which could be used as a new tool to detect particle ellipticity down to 10 −4 (corresponding to atomic accuracy comparable to TEM) as well as particle orientation. This is shown in Fig. S5. The FWM amplitude ratio in the focus center scales linearly with the particle ellipticity a/b − 1 (where we assumed b = c and a, b and c aligned along the x, y and z axis respectively). Furthermore the distance between the two displaced amplitude minima scales with the square root of the ellipticity. Fig. S5a shows these dependencies (red lines) and the experimental values (dashed lines) observed in Fig. 4 of the main manuscript, which indicate an ellipticity of about 0.5%. The occurrence of higher order modes in the measured data could also be used with appropriate modelling (beyond the the Rayleigh regime) to determine more details about the specific shape of the measured particle. Fig. S5b shows that the finite value of the phase of the FWM ratio in the focus center is given by -2γ plus an offset, with the particle orientation angle γ defined as the angle between the longer particle axis and the x-axis (see sketch in Fig. S5). Conversely, this phase is almost independent on the ellipticity (see Fig. S6). Furthermore, the orientation angle of the internodal axis relative to the x-axis Firstly, we note that for a circularly polarized incident field, the longitudinally polarized field component (E z ) in the focal plane exhibits an l = 1 optical vortex, while the crosscircularly polarized component is an l = 2 vortex. This is shown in Fig. S7 for the case of a 1.45 NA objective, as calculated in the main manuscript in Fig. 1b.
In order to reproduce the experimental findings in Fig. 6a of the main manuscript, we assumed a particle polarisability tensor representing a rotationally-averaged elliptical particle (hence an effectively spherical particle with equal semiaxis) plus a contribution coupling the longitudinal field component into the (x, y) plane. Therefore, the particle polarisability tensor was described as:α with α 0 given by Eq. S1 for a = b = c = 25 nm in a surrounding medium consisting of water with index n = √ m = 1.333, and C being a complex coefficient, the amplitude of which was adjusted to reproduce the strength of the FWM amplitude ratio in the experiment (for the calculations in Fig. 6c we used |C|=0.37). ii. Gold nanoparticles in agarose gel Before use, glass slides and coverslips were cleaned with acetone. 500 mg of agar powder (high molecular biology grade agarose, Bioline Cat. 41025) was placed in a conical glass beaker filled with 10 mL of pure water, and the mixture was boiled in a microwave for about 60 s until a homogeneous liquid was formed. Nominally spherical gold NPs of 25 nm radius (BBI Solutions) in pure water were added to the mixture to achieve a typical concentration of 10 9 NPs/mL. A chamber was formed by attaching an adhesive imaging spacer of 0.12 mm thickness with a 13 mm diameter hole (GraceBioLabs) onto a glass coverslip. 20 µL volume of the gold NP agar mixture was pipetted into the chamber and sealed with a glass slide.

iii. Gold nanoparticles in fixed cells
HeLa cells were seeded onto gridded coverslips and then loaded with 20 nm radius gold NPs via clathrin mediated endocytosis, using the transferrin (Tf) ligand attached onto the NP via a biotin-streptavidin conjugate. Tf-biotin was purchased from Sigma Aldrich and 20 nm radius gold NPs covalently bound to streptavidin were purchased from Innova Biosciences. Cells were placed on ice for 10 min to inhibit endocytosis, incubated with 50 µg/mL Tf-biotin in ice-cold serum-free medium for 8 min and then washed 3 times in ice-cold phosphate-buffered saline (PBS) at pH 7.4. Cells were then incubated with a dilution of gold NP-streptavidin conjugate in ice-cold serum-free medium, followed by washing 3 times in PBS pH 7.4. Cells were subsequently incubated in pre-warmed imaging medium for 6 hours to enable endocytosis, and then fixed in 3% paraformaldehyde for 10 mins at room temperature. They were washed 3 times in PBS at room temperature and mounted onto a glass slide using a DAKO (Dako UK Ltd) mounting medium at 80%.

iv. Lock-In detection of the modulation
A digital lock-in amplifier Zurich Instruments HF2LI was used to detect the FWM signal at both sidebands ν 2 ± ν m − ν L = (2 ± 0.4)MHz together with the reflected probe at the carrier ν c = ν 2 − ν L = 2MHz (ν 2 =82 MHz is the frequency up-shift of the probe field via AOM 2 , ν m = 0.4 MHz is the intensity modulation of the pump via AOM 1 , and ν L = 80 MHz is the pulsed laser repetition rate, see also Fig. 1 in main paper). To produce a modulation of the pump intensity, we supplied a modulation at ν m to the digital input of the AOM 1 driver (Intraaction DFE-834C4-6), generated by the internal oscillator of the lock-in. To provide an electronic reference to the lock-in at the carrier frequency of 2 MHz, we assembled an external electronic circuit using coaxial BNC components (Minicircuits), mimicking the optical mixing. The laser repetition frequency was taken from the 80 MHz photodiode signal output of the Ti:Sa laser (Newport/Spectra Physics MaiTai). After filtering for the first harmonics using a 50 MHz high pass (BHP-50+) and a 100 MHz low-pass (BLP-100+), the signal was mixed with the frequency ν 2 from the +10dBm reference output of the AOM 2 driver (Intraaction DFE-774C4-6) using an electronic mixer (ZP-1LH+) and a 10 MHz low pass filter (BLP-10.7+) to produce an electric sine wave at the difference frequency ν 2 −ν L = 2 MHz. The latter was then converted from sine to TTL using a comparator (Pulse Research Lab -PRL-350TTL-220), to be suited as external reference for the HF2LI using the DIO0 input (the two available high-frequency analog inputs were used for the two balanced-diode signals).
To discuss the analysis of the dual sideband detection by the lock-in, we start by describing the detected voltage of the balanced-diode signal as the real part of U = a c e iωct (1 + a m cos(ω m t + ϕ)) . (S5) with the complex carrier amplitude a c , which is measured directly by the ZI-HF2LI via the in-phase (Re) and in-quadrature (Im) components at the carrier frequency ν c = ω c /(2π), and the complex relative modulation of the carrier a m with the frequency 0 < 2πν m = ω m < ω c and the phase ϕ. The modulation is due to the intensity modulation of the pump field, which is a real quantity and can thus be described by a cosine. This can be rewritten as U = a c e iωct + a c a m 2 e iϕ e i(ωc+ωm)t + a c a m 2 e −iϕ e i(ωc−ωm)t .
Accordingly, the complex amplitudes of the lower and upper sidebands are a ± = a c a m e ±iϕ /2, which are detected together with a c by the ZI-HF2LI using the dual sideband modulation detection. The modulation phase ϕ can be determined modulo π by and is due to the delay τ between the modulation output of the lock-in and the detection of the amplitude modulation by the lock-in, specifically ϕ = −ω m τ , resulting in the modulation term cos(ω m (t − τ )). The modulation of the detected probe field, which is the FWM signal, is given by The determination of ϕ modulo π leads to an uncertainty of the sign of the modulation.
The phase shift ϕ can be determined including its sign by directly detecting the reflected pump beam which is frequency up-shifted via AOM 1 by the amount ν 1 = 83 MHz, using ν 1 − ν L instead of ν 2 − ν L in the lock-in reference, and knowing that a m > 0. Since the phase shift is given by electronic and acousto-optic delays it is stable to a few degrees. Once ϕ is known, the four-wave mixing field a m a c is determined using Eq.(S8).
v. Shot-noise limited detection: Dependence on reference power The experimental noise was evaluated by taking the statistical distribution of the measured FWM field (in both Re and Im components) in a spatial region away from the particle where no detectable signal is present. The standard deviation σ of this distribution was deduced for each component. For the measurement conditions as in Fig. 4 in the main paper, σ was found to be the identical in the Re and Im component, and for the co-polarised and cross-polarised components, as expected for an experimental noise given by the shot-noise of the reference beam and the electronic noise. Fig. S9 shows σ as a function of the total power P R of the reference beam entering the balanced detector, which is split onto the 4 We used P R ∼ 3 mW in the experiments shown.
vi. Pump-induced change of particle polarisability: Pump power and delay dependence An example of the co-polarised FWM field amplitude A + FWM for a 30 nm radius NP at the focus center as a function of the pump-probe delay time and power is shown in Fig.S10.
As discussed in detail in our previous work [14], the pump-probe delay dependence reflects Fig. S10 follows the expected third-order nonlinear behavior for small powers, and the onset of saturation at higher power. Pump power (P 1 ) and probe power (P 2 ) were increased while keeping a constant ratio 2P 2 = P 1 . At small powers, in the limit of a pump-induced change of the particle polarisability linear with the pump intensity, the FWM field amplitude scales as P 1 √ P 2 ∝ P 3/2 1 .

vii. Background-free FWM detection in cells
To exemplify that our FWM detection is very specific to gold NPs and free from background even in highly scattering and fluorescing environments, Fig. 3

5.5).
Notably, FWM acquisition can be performed simultaneously with confocal fluorescence microscopy for correlative co-localisation analysis. An example of this simultaneous acquisition is shown in Fig. S11. Confocal epi-fluorescence detection was implemented in the same microscope set-up, by conjugating the sample plane onto a confocal pin-hole of adjustable opening in front of a photomultiplier detector (Hamamatsu H10770A-40). Excitation occurred via the same laser beam used for FWM; fluorescence was collected via the same objective (in this case 1.45 NA oil-immersion) and spectrally separated using a dichroic beam splitter for pick-up and a bandpass filter transmitting in the 600-700 nm wavelength range in front of the photomultiplier. Fig. S11 shows 3T3 cells loaded with 10 nm radius gold NPs using a protocol similar to the one described above for HeLa cells. Fixed cells were image with DIC, reflection, FWM and simultaneous confocal fluorescence. The co-circularly polarized reflection image A + 2r in Fig. S11b correlates with the cell contour seen in DIC in the region highlighted by the dashed frame, and shows a spatially varying contrast due to thickness and refractive index inhomogeneities in the sample. Epi-fluorescence in Fig. S11c correlates with the cell morphology seen in DIC and reflection image, and is dominated by autofluorescence near the cell nucleus. The co-circularly polarised FWM amplitude A + FWM shown in Fig. S11d is a maximum intensity projection over a 3 µm z-stack and clearly indicates the location of single gold NPs in the cell. FWM imaging is background-free (throughout the z-stack) despite the significant scattering and autofluorescing cellular contrast seen in  NPs. The heterogeneity of the observed patterns is consistent with the heterogeneous NP shape and asymmetry observed in TEM.

ix. FWM measurements on 5 nm radius nanoparticles
The FWM amplitude ratio A − FWM /A + FWM in the focus center was measured on nominally spherical 5 nm radius gold NPs of the type shown in Fig. S8b, and compared with the ratio obtained on nominally spherical 30 nm radius gold NPs of the type shown in Fig. S8a.
The results are summarized in Fig. S13. The shot-noise limit σ/A 0 in these experiment is indicated by the dashed lines. Notably, we observe that for the 5 nm radius NPs, a large proportion (∼ 70%) has a FWM amplitude ratio A − FWM /A + FWM ≤ 0.02 in the center of the focal plane, in this case limited by the signal-to-noise ratio rather than the particle asymmetry. Conversely, almost all NPs of nominal 30 nm radius have A − FWM /A + FWM > 0.02 which is well above the shot-noise limit and is due to the particle asymmetry.
x. Parameters for position localisation of a freely rotating single NP From the measurements on a freely rotating NP in an agarose gel shown in the main manuscript in Fig. 6 we have deduced the relationship between the FWM amplitude ratio A − FWM /A + FWM and phase Φ − FWM − Φ + FWM and the position coordinates in the focal plane, as measured from the sensor positions in the scanned piezoelectric sample stage. This is shown in Fig. S14. The (x, y) experimental scans in Fig. 6 were first converted in 2D images along the polar coordinates R, ϕ. The data for A − FWM /A + FWM were averaged over the full range of the ϕ coordinate to obtain the one-dimensional dependence versus R shown in Fig. S14. Similarly, the data for Φ − FWM − Φ + FWM were averaged in R over a region of about 60 nm < R < 100 nm to obtain the dependence versus ϕ. The FWM field phase Φ + FWM versus z also shown in Fig. S14 was taken from rapid 3D scans. We find that Φ + FWM varies linearly with z, as discussed in Section S1.iv, with the slope ∂z/∂Φ = 37 nm/rad (note that the experiments in the main manuscript in Fig. 6 were performed using a 1.27 NA waterimmersion objective). We also find that A − FWM /A + FWM has a simple linear dependence in R given by R/R 0 = A − FWM /A + FWM with R 0 = 1.33 µm in the focal plane. Moreover, the phase of the cross to co-polarised FWM ratio Θ = Φ − FWM − Φ + FWM versus polar angle follows the relationship −ϕ = Θ − Θ 0 with Θ 0 = 0.22 for the data shown in Fig. S14. in comparison with those observed during a scan where the highest value of A − FWM in the focus is measured. This is shown in Fig. S15. The behaviour is quite similar on both traces, dominated by systematic effects in the z coordinate due to a cross-coupling to the xy scanning, and a slow drift in z -note that the z range shown is only 50 nm. Superimposed, we find random variations in z which are larger than the shot noise (estimated to be ≤ 3 nm under the experimental conditions in Fig. S15). Based on our calculation of the effect of NP asymmetry on the cross-circularly polarized FWM (see also section S1.vi), the trace corresponding to the large A − FWM is attributed to the NP not rotating but stuck in a position where the NP ellipticity determines A − FWM in the focus center (from which we estimate the NP ellipticity to be a/b − 1 = 0.08). For the trace corresponding to the observation of the l = 1 mode in A − FWM , we do observe larger jumps in the z-coordinate and the range spanned is slightly larger, consistent with the idea that the NP is freely rotating in 3D in this case, and hence also bouncing about within the agarose gel pocket in which it is trapped.
To estimate the effect of restricted diffusion of a NP confined in a cage in agar, as sketched in Fig. S16, we assume a NP radius of R P and a cage radius of R C . The free diffusion of the NP in water is given by for each dimension of space, where k B is Boltzmann's constant, T is the temperature, and η is the dynamic viscosity of water. Here we are interested only in the z coordinate. The boundary of the cage will limit this diffusion, confining the particle motion to a range 2(R C − R P ), which the particle explores in the characteristic time If the measurement integration time τ is smaller than τ C , the measured position fluctuation will explore the full cage size. If instead τ > τ C , the measured position will be an averaged position, with less fluctuations, an effect also known as motional narrowing. The standard deviation of the position fluctuations scales as the inverse root of the number of averages, hence is given by σ = 2(R C − R P ) 12(1 + τ /τ C ) (S11) where the 1/ √ 12 factor originates from the standard deviation of a uniform distribution.
The resulting σ as function of the cage radius is given in Fig. S17. We find that due to the motional narrowing, σ is suppressed to below 5 nm for R C 45 nm, for which the NP has 40 nm space to move and thus also to rotate. The observed fluctuations in Fig. S15, which have a σ in in the order of 5 nm, are therefore consistent with the assumption of free rotational diffusion.
xii. Randomly oriented non-rotating particle We cannot reproduce the experimental finding of an l = 1 optical vortex in the crosscircularly polarized FWM by assuming a non-rotating randomly-oriented asymmetric particle. This is shown in Fig. S18 where we have calculated the statistical distribution of the FWM ratio A − FWM /A + FWM in the focus center assuming an asymmetric particle ellipsoid with axis a = b(1+e 1 ) and c = b(1+de 1 ) using the ellipticity e 1 = 0.0825 and considering random orientations in 3D. Note that by varying the parameter d from 1 to 0 we change the particle shape from oblate to prolate. The histograms shown in Fig. S18 are significantly different  For practical purposes, we implemented a simple way to achieve in situ calibration of the in-plane NP position coordinates without the need for prior knowledge and/or characterization of the particle optical response. This is achieved by applying a small oscillation of known amplitude (here 16 nm) at 25 Hz frequency to the x axis of the sample stage, which is accurately measured by the position sensor in the stage. When the FWM field ratio from the NP encodes the NP in-plane position as discussed in the main manuscript, this oscillation is detected in the measured E − FWM /E + FWM . To find the amplitude A of this oscillation (as a complex number i.e. including its phase) we then Fourier transform E − FWM /E + FWM and multiply it with a Gaussian filter centered at 25 Hz (with 0.5 Hz FWHM). We also Fourier transform the detected x position sensor from the sample stage to deduce the corresponding complex amplitude x 0 . We can then calibrate E − FWM /E + FWM in terms of NP in-plane positions as X = Re(( x 0 / A)E − FWM /E + FWM ) and Y = Im(( x 0 / A)E − FWM /E + FWM ). Conversely, if a NP exhibits a FWM field ratio dominated by the shape asymmetry, we do no longer observe this oscillation in E − FWM /E + FWM . This case is shown in the paper in Fig.7c-f. The power spectrum for this case is shown in Fig. S19.