# Spatial super-resolution in nanosensing with blinking emitters

Alexander Mikhalychev<EMAIL_ADDRESS>Atomicus
Sp. z o.o., Ul. Aleja Grunwaldzka 163, Gdansk, 80-266 Poland  Atomicus GmbH,
Amalienbadstr. 41C, Karlsruhe, 76227 Germany  Aleksandr Saushin  University of
Eastern Finland, Yliopistokatu 7, Joensuu, 80101 Finland  Alex Ulyanenkov
Atomicus GmbH, Amalienbadstr. 41C, Karlsruhe, 76227 Germany  Polina Kuzhir
University of Eastern Finland, Yliopistokatu 7, Joensuu, 80101 Finland

###### Abstract

We propose a method of spatial resolution enhancement in metrology
(thermometry, magnetometry, pH estimation and similar methods) with blinking
fluorescent nanosensors by combining sensing with super-resolution optical
fluctuation imaging (SOFI). To demonstrate efficiency of this approach, a
model experiment with laser diodes modeling fluctuating nanoemitters and
intentional blurring of the image is performed. The 2nd, 3rd, and 4th order
cumulant images provide improvement of the contrast and enable successful
reconstruction of smaller features of the modeled temperature (or any other
physical parameter) distribution relatively to the intensity-based approach.
We believe that blinking fluorescent sensing agents being complemented with
the developed image analysis technique could be utilized routinely in the life
science sector for recognizing the local changes in the spectral response of
blinking fluorophores, e.g. delivered targetly to the wanted cell or even
organelle. It is extremely useful for the local measurements of living cells’
physical parameters changes due to applying any external “forces”, including
disease effect, aging, healing or response to the treatment.

nanosensors, quantum dots, fluorescence, super-resolution optical fluctuation
imaging, spatial resolution, sub-Rayleigh imaging methods

##  I Introduction

Fluorescent nanosensors represent a powerful tool for measurement of
temperature, magnetic field, pH, and other parameters with high accuracy at
small scales [1, 2, 3, 4, 5, 6]. In particular, the dependence of their
fluorescent signal on the mentioned parameters can be utilized for
investigation of processes in living cells [2, 7, 8]. To ensure high spatial
resolution, one needs to localize the nanoemitters precisely and to
distinguish signals from closely located nanosensors. This requirement can be
fulfilled, for example, by placement of fluorescent nanodiomonds or quantum
dots at a manipulator needle [9, 10, 11, 12, 13]. This approach requires
precise mechanical control and does not allow to analyze environment
parameters at several points of the sample simultaneously. For a large number
of nanosensors, simultaneously attached to the sample, one can apply
deterministic super-resolution approaches based on selective (localized)
excitation of emitters (stimulated emission depletion (STED), ground state
depletion (GSD), reversible saturable optical fluorescence transitions between
spin levels (spin-RESOLFT), structured illumination imaging) [14, 15, 16, 17,
18], as well as stochastic methods (stochastic optical reconstruction
microscopy (STORM), photoactivated localization microscopy (PALM), points
accumulation for imaging in nanoscale topography (PAINT)) [19, 20, 21, 22],
ensuring independent localization of single emitters due to their blinking.
Those approaches are also quite sophisticated and require the nanosensors to
remain in the inactive state most of the time, thus increasing needed
exposition time.

Alternatively to the super-resolution approaches, based on scanning or
selective excitation (activation) of emitters, wide-field microscopy methods
exploiting specific temporal statistics of fluorescent objects, i.e.
stochastic blinking (in super-resolution optical fluctuation imaging — SOFI)
[23, 24, 25] or photon antibunching [26, 27, 28] have been successfully
implemented. The idea is based on postprocessing of time-dependent fluorescent
signal: its efficient exponentiation by considering correlations instead of
raw intensity, elimination of cross-terms from different emitters, and
subsequent enhancement of spatial resolution. Analysis of the fluorescence
blinking does not prevent usage of the “degrees of freedom” suitable for the
environment sensing, such as the dependence of the spectrum or fluorescence
lifetime on the temperature, acidity, pressure, etc. In the current
contribution we suggest to combine sensing capabilities of nanoparticles with
SOFI-based super-resolution for more accurate localization of the signal
source. Is is important to note that the information, used by SOFI to achieve
sub-Rayleigh spatial resolution, is already present in the collected
fluorescent signal and does not introduce any additional complications to the
measurement.

The idea of combined usage of blinking statistics of fluorophores for super-
resolution and certain other degree of freedom for getting additional
information is somewhat close to so called “color SOFI” [29, 30, 31, 32] —
super-resolution technique, which applies SOFI to several (typically two)
spectral channels. However, the tasks, solved by those two approaches, are
different. Color SOFI deals with discrimination between several (two or three)
discrete types of fluorescent emitters, serving as biological labels. In
contrast, the method proposed in this paper, deals with continuous-variable
sensing: the emitters are of the same type, but their fluorescent signal
differs due to variations of local environment.

To demonstrate principles of the proposed approach and show its feasibility,
we performed a model experiment mimicking blinking fluorescent nanoemitters.
We used classical photodiods with the controlled signal imitating random
switching between “bright” and “dark” states. The image was intentionally
blurred by inaccurate focusing on the detection plane so that the spots
produced by neighbor “emitters” were strongly overlapping. Model “temperature”
of the “emitters” was encoded in the “fluorescence lifetime” — the delay
between the “excitation pulse” and the emission response. The subsequent
application of the proposed SOFI-based sensing approach allowed to enhance
spatial resolution approximately twice.

The paper is organized as follows. In Section II, we discuss basic ideas of
the approach and list the procedures constituting the proposed technique of
spatially super-resolving sensing. Section III describes the model experiment,
data acquisition, and reconstruction results for the model “temperature”
distribution. Discussion of the results and possible directions of further
research are presented in Section IV.

##  II Basic ideas

###  II.1 Super-resolution optical fluctuations imaging

The commonly used model of wide-field fluorescence imaging includes a sample
composed of $N$ incoherent florescent emitters, located at positions
$\textbf{r}_{k}$ [23]. Fluctuations-based approaches assume that the emitters
are independently blinking, stochastically switching between “bright” and
“dark” states. The acquisition time is split into a finite number of frames.
The fluorescence signal in the frame $t$ ($t=1$, …, $N_{\text{f}}$, where
$N_{\text{f}}$ is the number of frames) can be described as [23]

$F(\mathbf{r},t)=\sum_{k}U(\mathbf{r}-\mathbf{r}_{k})\varepsilon_{k}s_{k}(t),$
(1)

where $U(\Delta\mathbf{r})$ is the point-spread function (PSF) of the optical
system, $\varepsilon_{k}$ is the constant molecular brightness of the $k$-th
emitter, and $s_{k}(t)$ is its time-dependent fluctuation ($s_{k}(t)=1$ if the
$k$-th emitters remains “bright” during the $t$-th frame; $s_{k}(t)=0$ for the
“dark” state).

The idea of resolution enhancement in the wide-field microscopy (including
quantum imaging, SOFI, antibunching imaging) is to construct signal with
power-law dependence on $U(\Delta\mathbf{r})$ and to take into account that,
for a localized PSF, $U^{n}$ is approximately $\sqrt{n}$ times narrower than
$U$ itself. For a single emitter ($N=1$), just raising the signal
$F(\mathbf{r},t)$ to certain power $n>1$ would be sufficient for achieving the
effect. For $N>1$ emitters, exponentiation of the signal leads to additional
cross-terms, hindering resolution of the emitters. For example, for $n=2$ the
squared signal equals

$F^{2}(\mathbf{r},t)=\sum_{k}U^{2}(\mathbf{r}-\mathbf{r}_{k})\varepsilon_{k}^{2}s_{k}^{2}(t)\\\
+\sum_{k,k^{\prime}:k^{\prime}\neq
k}U(\mathbf{r}-\mathbf{r}_{k})U(\mathbf{r}-\mathbf{r}_{k^{\prime}})\varepsilon_{k}\varepsilon_{k^{\prime}}s_{k}(t)s_{k^{\prime}}(t).$
(2)

To ensure resolution enhancement, one needs to construct a polynomial function
of the fluorescence signal with excluded cross-terms [23]. For example,
variance of the signal considered as a time series of frames both is quadratic
relatively to PSF and has zero cross-terms for independent blinking of the
emitters:

$\operatorname{Var}_{t}[F(\mathbf{r},t)]=\sum_{k}U^{2}(\mathbf{r}-\mathbf{r}_{k})\varepsilon_{k}^{2}\operatorname{Var}_{t}[s_{k}(t)].$
(3)

The approach can be extended to higher-order polynomials by using cumulants
instead of variance [23]:

$C_{t}^{(n)}[F(\mathbf{r},t)]=\sum_{k}U^{n}(\mathbf{r}-\mathbf{r}_{k})\varepsilon_{k}^{n}C_{t}^{(n)}[s_{k}(t)].$
(4)

The expression stems from linearity of cumulants relatively to a sum of
independent random variables:
$C^{(n)}[X_{1}+X_{2}]=C^{(n)}[X_{1}]+C^{(n)}[X_{2}]$ for independent $X_{1}$
and $X_{2}$.

Several lower-order cumulants for a random variable (process) $f(t)$ are
defined as

$C_{t}^{(2)}[f(t)]=\operatorname{E}_{t}\bigl{[}\bigl{(}f(t)-\operatorname{E}_{t}[f(t)]\bigr{)}^{2}\bigr{]},$
(5)
$C_{t}^{(3)}[f(t)]=\operatorname{E}_{t}\bigl{[}\bigl{(}f(t)-\operatorname{E}_{t}[f(t)]\bigr{)}^{3}\bigr{]},$
(6)
$C_{t}^{(4)}[f(t)]=\operatorname{E}_{t}\bigl{[}\bigl{(}f(t)-\operatorname{E}_{t}[f(t)]\bigr{)}^{4}\bigr{]}-3\bigl{(}C_{t}^{(2)}[f(t)]\bigr{)}^{2},$
(7)

and correspond to variance, asymmetry, and excess kurtosis respectively. Here,
$\operatorname{E}_{t}[f(t)]$ denotes time averaging of the random values
$f(t)$ (defined at each frame $t$) over an infinite number of frames
$N_{\text{f}}$. The random process of blinking is assumed to be stationary
(with negligible effect of bleaching) so that ergodicity can be applied to the
definition of cumulants. For finite acquisition time, expectation values are
estimated as finite-sample averages:

$\operatorname{E}_{t}[f(t)]\approx\frac{1}{N_{\text{f}}}\sum_{t=1}^{N_{\text{f}}}f(t).$
(8)

Strictly speaking, Eqs. (2) and (4) are valid for classical signals — in the
limit of a large number of detected photons, when the shot noise is small
relatively to the fluctuations caused by emitters’ blinking. Detailed
theoretical discussion of that assumption will be presented in Ref. [33].

To increase spatial resolution, SOFI deals with stochastic fluctuations of the
total fluorescence signal, integrated over each frame duration and certain
spectral range, rather than on spectral or temporal intensity distribution
inside a frame. The technique itself does not impose any specific constraints
on time or spectral windows, used during signal acquisition within frames.
That property is crucial for the proposed sensing approach and was efficiently
used in the color SOFI [29, 32].

###  II.2 Sensing with fluorescent nanoparticles

Sensing capabilities of fluorescent nanoparticles can be provided by
dependence of their fluorescence lifetime, spectrum, or total intensity on the
local parameters of the environment: temperature, pH, pressure, etc. [1, 2,
34, 35, 36, 37, 38, 39] (Fig. 1). Since blinking of the emitters introduces
significant noise to the total intensity of fluorescence, we focus on the
cases when the environment influences distribution of the fluorescent signal
over certain parameter (variable) $x$ within the frame. Let
$f(\mathbf{r},t,x)$ be the density of the signal over the variable $x$. Then,
the total fluorescent signal for the frame $t$ is defined as

$F(\mathbf{r},t)=\int f(\mathbf{r},t,x)dx.$ (9)

For spectral response of the fluorophores to local changes of environment, $x$
may correspond to the emitted radiation wavelength and $f(\mathbf{r},t,x)$ be
the spectral density of the signal. For sensing, based on fluorescence
lifetime like in the case of SiV centers in diamond needles [39], $x$ is the
delay between the excitation pulse and the detection of photons.

Figure 1: Sensing an environment parameter $\theta$ by its influence on
fluorescence signal: changes of the spectrum (a), fluorescence lifetime (b),
or total intensity (c). The latter case is incompatible with SOFI because of
blinking-induced noise of absolute intensity. Nonlinear dependence of the
estimator $\bar{\theta}=\Theta(\theta)$, defined by Eq. (15), can be taken
into account by appropriate calibration (d). Dotted arrows schematically
illustrate reconstruction of the environment parameter $\theta$ from the
estimator value $\bar{\theta}$ calculated from the fluorescence signal.

Eq. (1) can be rewritten in the following form for the density of the signal:

$f(\mathbf{r},t,x)=\sum_{k}U(\mathbf{r}-\mathbf{r}_{k})\varepsilon_{k}(x,\theta_{k})s_{k}(t),$
(10)

where the dependence on the $k$-th emitter intensity
$\varepsilon_{k}(x,\theta)$ on the variable $x$ (wavelength or delay of the
fluorescence signal) is affected by the local value of the environment
parameter $\theta=\theta_{k}$ in the vicinity of the $k$-th emitter. To infer
the local values $\theta(\mathbf{r})$ of the parameter in question, one needs
to specify certain model for the dependence of $\varepsilon_{k}(x,\theta)$ on
$\theta$ and, typically, to assume that all the emitters are described by the
same model
$\varepsilon_{k}(x,\theta_{k})=\varepsilon_{k}\varepsilon(x,\theta_{k})$ up to
a constant emitter-dependent factor $\varepsilon_{k}$. To eliminate additional
noise, caused by blinking of the emitters, it is worth normalizing the signal
by dividing it by its integral over $x$.

First, let us assume that the dependence of the emitter’s signal on the
parameter of environment is linear in the vicinity of the reference value
$\theta_{0}$:

$\varepsilon(x,\theta)=S_{0}(x)+S_{1}(x)(\theta-\theta_{0}).$ (11)

As it is shown in Appendix A, one can build an estimator

$\bar{\theta}=\theta_{0}+\frac{\int\varepsilon(x,\theta)Q_{1}(x)dx}{\int\varepsilon(x,\theta)Q_{0}(x)dx},$
(12)

where

$\int Q_{i}(x)S_{j}(x)dx=\delta_{ij},$ (13)

and ensure, by appropriate choice of the weights $Q_{0}(x)$ and $Q_{1}(x)$,
that the estimator retrieves all the information contained in the intensity
profile $\varepsilon(x,\theta)$ after its normalization (i.e., represents
sufficient statistics).

Eq. (12) indicates the fact that the profile $\varepsilon(x,\theta)$ can be
“compressed” to just two scalar characteristics (linear functionals)

$Y_{i}=\int Q_{i}(x)\varepsilon(x,\theta)dx,$ (14)

which are sufficient for retrieving all the available information about the
investigated parameter $\theta$. One can move further and consider the two
scalars $Y_{0}$ and $Y_{1}$ for some arbitrary weights $Q_{0}(x)$ and
$Q_{1}(x)$ without applying the constraints (13). Practically, those two
quantities may correspond to the signal accumulated over two spectral or
temporal windows. In that case (see Appendix B), one can still find two linear
combinations $\bar{Y}_{0}$ and $\bar{Y}_{1}$ of the scalars so that the
corresponding weights $\bar{Q}_{0}(x)$ and $\bar{Q}_{1}(x)$ satisfy Eq. (13).
However, such measurement-induced choice of the weights is, generally, non-
optimal and may lead to loss of information relatively to the full profile
analysis. The estimator for the environment parameter $\theta$ can be
constructed as

$\bar{\theta}=\theta_{0}+\frac{Y_{1}A_{00}-Y_{0}A_{10}}{Y_{0}A_{11}-Y_{1}A_{01}},$
(15)

where

$A_{ij}=\int Q_{i}(x)S_{j}(x)dx,\quad i,j=0,1,$ (16)

represent coefficients in decomposition of $Y_{i}$ over powers of the
deviation $\theta-\theta_{0}$:

$Y_{i}=A_{i0}+A_{i1}(\theta-\theta_{0}).$ (17)

In general case, the dependence of $\varepsilon(x,\theta)$ on $\theta$ is
nonlinear. Considering small deviations $|\theta-\theta_{0}|$, one can
linearize $\varepsilon(x,\theta)$ and define

$S_{0}(x)=\varepsilon(x,\theta_{0}),\quad
S_{1}(x)=\varepsilon^{\prime}_{\theta}(x,\theta_{0}).$ (18)

Then, following the procedures discussed above, one can construct the weights
$Q_{0}(x)$ and $Q_{1}(x)$ and calculate the quantities $Y_{i}$ according to
Eq. (14) for the given signal $\varepsilon(x,\theta)$. The dependence of the
linear estimator defined by Eq. (15) on the true value of the parameter
$\theta$ is, in general, also nonlinear and biased,
$\bar{\theta}=\Theta(\theta)\neq\theta$. Assuming that the dependence is
monotonically increasing in certain vicinity of the reference value
$\theta_{0}$, one can construct its inverse and build the estimator as

$\bar{\theta}_{\text{nonlin}}=\Theta^{-1}(\bar{\theta}).$ (19)

Practically, the procedure means plotting the calibration curve
$\bar{\theta}=\Theta(\theta)$ and using it for finding the parameter value
$\theta=\Theta^{-1}(\bar{\theta})$ for the value $\bar{\theta}$ calculated
from the input data $\varepsilon(x,\theta)$ (Fig. 1(d)).

###  II.3 Combining sensing and SOFI

The proposed idea of combining spatial resolution with sensing consists in
applying SOFI separately to several channels relatively to the variable $x$,
which defines signal distribution (spectral or temporal) within each frame.
Naively, one can just replace the integrated signal $F(\mathbf{r},t)$ by its
density $f(\mathbf{r},t,x)$ in Eq. (4) and get spectrally or temporarily
resolved cumulant image:

$C_{t}^{(n)}[f(\mathbf{r},t,x)]=\sum_{k}U^{n}(\mathbf{r}-\mathbf{r}_{k})\varepsilon_{k}^{n}\varepsilon^{n}(x,\theta_{k})C_{t}^{(n)}[s_{k}(t)].$
(20)

Practically, such approach would mean that one takes the binned signal density
$f(\mathbf{r},t,x_{i})$, specified for a discrete set of values $x_{i}$, and
performs analysis of frame-wise temporal statistics to each bin $x_{i}$ by
considering the time series $\\{f(\mathbf{r},t,x_{i}\\}_{t}$. However, such
time series is unlikely to satisfy the requirement of having enough photons in
each frame, which is important for reliable application of SOFI.

Such issue can be resolved by the above-discussed “compression” of a full
profile over $x$ into several scalar functionals. Following Eq. (14), one can
define the signals, integrated over two spectral or temporal windows (indexed
as $j=0,1$) within each frame, as follows:

$F_{j}(\mathbf{r},t)=\int f(\mathbf{r},t,x)Q_{j}(x)dx\\\
=\sum_{k}U(\mathbf{r}-\mathbf{r}_{k})\varepsilon_{k}Y_{kj}s_{k}(t),$ (21)

where $Y_{kj}$ is calculated according to Eq. (14) for $\theta=\theta_{k}$.

In contrast to the binned density $f(\textbf{r},t,x_{i})$, the integrated
signals $F_{j}(\mathbf{r},t)$ are intense enough to be exploited for SOFI. The
cumulants can be constructed for each window separately:

$C_{j}^{(n)}(\mathbf{r})=C_{t}^{(n)}[F_{j}(\mathbf{r},t)]\\\
=\sum_{k}U^{n}(\mathbf{r}-\mathbf{r}_{k})\varepsilon_{k}^{n}Y_{kj}^{n}C_{t}^{(n)}[s_{k}(t)].$
(22)

The quantities $C_{0}^{(n)}(\mathbf{r})$ and $C_{1}^{(n)}(\mathbf{r})$ have
the meaning close to $Y_{0}$ and $Y_{1}$ and can be used for estimation of the
local parameter value $\theta(\mathbf{r})$ similarly to Eq. (15). If the
deviations $|\theta_{k}-\theta_{0}|$ of the parameter values $\theta_{k}$ at
the emitters’ positions from the reference value $\theta_{0}$ are small, the
factor $Y_{kj}^{n}$ in Eq. (22) can be decomposed in the way similar to Eq.
(17):

$Y_{kj}^{n}\approx A_{j0}^{n}+nA_{j0}^{n-1}A_{j1}(\theta_{k}-\theta_{0}).$
(23)

Comparison of Eqs. (17) and (23) suggests that one can construct the estimator
for the local value $\theta(\mathbf{r})$ of the environment parameter as

$\bar{\theta}^{(n)}(\mathbf{r})=\theta_{0}\\\
+\frac{C_{1}^{(n)}(\mathbf{r})A_{00}^{n}-C_{0}^{(n)}(\mathbf{r})A_{10}^{n}}{n[C_{0}^{(n)}(\mathbf{r})A_{10}^{n-1}A_{11}-C_{1}^{(n)}(\mathbf{r})A_{00}^{n-1}A_{01}]}.$
(24)

Substituting Eqs. (22) and (23) into Eq. (24), one can see that the
constructed estimator indeed represents the parameter distribution, blurred by
the diffraction:

$\bar{\theta}^{(n)}(\mathbf{r})=\theta_{0}\\\
+\frac{\sum_{k}U^{n}(\mathbf{r}-\mathbf{r}_{k})\varepsilon_{k}^{n}C_{t}^{(n)}[s_{k}(t)](\theta_{k}-\theta_{0})}{\sum_{k}U^{n}(\mathbf{r}-\mathbf{r}_{k})\varepsilon_{k}^{n}C_{t}^{(n)}[s_{k}(t)]}.$
(25)

The effect of blurring reduction for higher-order cumulants, characteristic
for SOFI, is still present in Eq. (25): the PSF $U$ is raised to the power
$n$.

For larger deviations $\theta-\theta_{0}$, the nonlinearity of the response
should be taken into account similarly to Eq. (19). For a uniform distribution
of the parameter $\theta$ ($\theta(\mathbf{r})=const=\theta$), one has
$C_{j}^{(n)}(\mathbf{r})\propto Y_{j}^{n}$ and can find the (nonlinear)
dependence $\bar{\theta}^{(n)}=\Theta_{n}(\theta)$ of the estimator
$\bar{\theta}^{(n)}$ on the parameter value $\theta$, calculated according to
Eq. (24) with the substitution $C_{j}^{(n)}(\mathbf{r})\mapsto Y_{j}^{n}$.
Then, in the general case of non-uniform parameter distribution, the nonlinear
estimator will be

$\bar{\theta}^{(n)}_{\text{nonlin}}(\mathbf{r})=\Theta_{n}^{-1}(\bar{\theta}^{(n)}(\mathbf{r})),$
(26)

where $\bar{\theta}^{(n)}(\mathbf{r})$ is calculated according to Eq. (24).
General properties of the estimator are discussed in Appendix C.

###  II.4 Reconstruction of parameter distribution for noisy signal

While for ideal noiseless data the constructed estimator (26) would provide
reasonable local values for the investigated parameter distribution
$\theta(\mathbf{r})$, its application to the noisy data (especially for
sparsely located emitters) deserves additional attention. The structure of
expression (22) implies that the cumulant values $C_{j}^{(n)}(\mathbf{r})$
estimated from finite-sample data are more reliable near emitters
($\mathbf{r}\approx\mathbf{r}_{k}$ for certain $k$), but may contain
significant relative shot noise contribution in the region between the
emitters where the signals are weak. The effect becomes more pronounced for
higher-order cumulants, since the effective PSF
$U^{n}(\mathbf{r}-\mathbf{r}_{k})$ is narrower for larger $n$, while the
contribution of shot noise is generally growing for larger cumulant orders.

The inaccuracy $\Delta\bar{\theta}(\mathbf{r})$ of the reconstructed parameter
value
$\bar{\theta}(\mathbf{r})=\bar{\theta}_{\text{nonlin}}^{(n)}(\mathbf{r})$ is
determined by the relative contribution of the shot noise to the quantities
$C_{0}^{(n)}(\mathbf{r})$ and $C_{1}^{(n)}(\mathbf{r})$, which is, roughly,
inversely proportional to the square root of the detected signal:

$\Delta\bar{\theta}(\mathbf{r})\propto\frac{1}{\sqrt{\left|C_{0}^{(n)}(\mathbf{r})+C_{1}^{(n)}(\mathbf{r})\right|}}.$
(27)

The likelihood of getting the estimated distribution
$\bar{\theta}(\mathbf{r})$ for the true distribution $\theta(\mathbf{r})$ can
be represented as

$L\left[\bar{\theta}(\mathbf{r})|\theta(\mathbf{r})\right]\propto\exp\left(-\int
d^{2}\mathbf{r}\frac{[\bar{\theta}(\mathbf{r})-\theta(\mathbf{r})]^{2}}{2\Delta\bar{\theta}^{2}(\mathbf{r})}\right).$
(28)

In the regions with low signal, the inaccuracy
$\Delta\bar{\theta}(\mathbf{r})$ can be large, making the estimate
$\bar{\theta}(\mathbf{r})$ unreliable. The estimate can be refined in the
spirit of Tikhonov’s regularization [40]. Very sharp features of the
reconstructed distribution of the environment parameter $\theta$ from a
blurred image with finite density of fluorescent emitters are likely to be
artifacts rather than to provide useful information about the sample. The
statement can be captured as a reduced prior probability of the distributions
with large values of second derivatives:

$p_{0}[\theta(\mathbf{r})]\propto\exp\left(\frac{1}{2D^{2}}\int
d^{2}\mathbf{r}\left|\nabla^{2}\theta(\mathbf{r})\right|^{2}\right),$ (29)

where the parameter $D$ represents the scale for the derivatives and
$\nabla^{2}=\partial^{2}/\partial x^{2}+\partial^{2}/\partial y^{2}$ is
2-dimensional Laplacian.

According to Bayes’ formula [41, 42], the posterior probability of the true
distribution $\theta(\mathbf{r})$ conditioned by the estimated distribution
$\bar{\theta}(\mathbf{r})$ equals

$p\left[\theta(\mathbf{r})|\bar{\theta}(\mathbf{r})\right]\propto
L\left[\bar{\theta}(\mathbf{r})|\theta(\mathbf{r})\right]p_{0}[\theta(\mathbf{r})].$
(30)

The distribution $\theta(\mathbf{r})$ can be reconstructed according to
maximum a posteriori method:

$\theta(\mathbf{r})=\operatorname*{argmax}\limits_{\theta(\mathbf{r})}p\left[\theta(\mathbf{r})|\bar{\theta}(\mathbf{r})\right].$
(31)

Practically, the parameter values are estimated at a discrete grid of points
$\mathbf{r}_{ij}$. In that case, the derivatives in Eq. (30) can be replaced
by finite differences, and the values $\theta(\mathbf{r}_{ij})$ can be found
by unconstrained quadratic minimization (Appendix D).

###  II.5 Sensing method summary

To realize the proposed sensing approach with spatial super-resolution
capabilities, the following procedures should be performed.

1\. Model construction. The models for parameter-dependent spectral or
temporal profiles $\varepsilon(x,\theta)$ of the emitters’ fluorescence and
for the emitters’ blinking should be constructed (either theoretically or from
empirical calibration procedure).

2\. Sensing calibration. The optimal weights $Q_{0}(x)$ and $Q_{1}(x)$ of the
two detection channels are defined in Appendix A. If the experimental
capabilities do not allow to perform such optimal measurement, the actual
weights $Q_{0}(x)$ and $Q_{1}(x)$ are defined by the available detection setup
(the profiles of the channels’ sensitivity). For sensing based on the $n$-th
order cumulant, the linear estimator is defined according to Eq. (24). The
nonlinear relation $\bar{\theta}^{(n)}=\Theta_{n}(\theta)$ between the
estimator and the true value of the parameter is found by theoretical
consideration or direct measurements of the sample response for a uniform
distribution of the parameter $\theta(\mathbf{r})=const$. The final nonlinear
estimator is constructed according to Eq. (26).

3\. Data acquisition. The fluorescent response of the sample is collected in
terms of 2-channel framed images $F_{k}(\mathbf{r}_{ij},t)$, where $k=0$, 1 is
the channel index, $t$ is the frame index, and $\\{\mathbf{r}_{ij}\\}$ is the
set of detection positions (pixel positions for an array detector or the
scanning positions).

4\. Data processing. For each channel $k$ and each detection position
$\mathbf{r}_{ij}$ separately, the cumulant value
$C_{k}^{(n)}(\mathbf{r}_{ij})=C_{t}^{(n)}[F_{k}(\mathbf{r}_{ij},t)]$ is
calculated. For $n=2$, 3, 4 the expressions are provided by Eqs. (5)–(7) with
the expectation values replaced by finite-sample averages. The local estimates
of the parameter $\bar{\theta}_{\text{nonlin}}^{(n)}(\mathbf{r})$ are
calculated according to Eq. (26).

5\. Data interpolation. To account for the difference of the estimate
reliability at different points, the corrected parameter distribution is
constructed according to Eq. (31) with the method described in Appendix D.

##  III Model experiment

###  III.1 Setup description

To verify the proposed technique, we performed the following model experiment.
We used a paper screen as a model of fluorescent sample. The screen was
irradiated with 13 KY-008 lasers controlled by the microcontroller board
Arduino MEGA [43]. The image of the screen was recorded with Thorlabs
DCC1545M-GL camera with 50 mm lens (Fig. 2). To simulate the diffraction
limit, the lens was placed 2 mm closer to the camera than it was necessary to
get a clear picture. The images of clear and unfocused laser spots are shown
in Fig. 3. Each laser played the role of a fluorescent particle.
Correspondingly, each particle had prescribed “relaxation time” meaning the
time between “excitation pulse” and “fluorescent emission”. To denote
“excitation pulse” in the recording data, Arduino made all lasers blink two
times. Since the minimal distinguishable time gap is 1 frame (about 43 ms), we
used frames as a measure of time. Therefore, the “excitation pulse” was a
series of 2 one-frame flashes of lasers separated by 1 frame of darkness. The
“fluorescent emission” was 1-frame flashing of a laser. To simulate the
fluorescent particles with different temperature, the lasers were divided into
“hot” and “cold” ones. For the hot lasers, the “relaxation time” was 5 frames
with standard deviation of 2 frames, while the cold lasers had “relaxation
time” of 10 frames with the same standard deviation.

Figure 2: The scheme of model experiment. The lasers controlled with Arduino
microcontroller board irradiate the paper screen. The light scattered by the
screen is collected with lens to the Thorlabs camera. Figure 3: The images of
laser spots on the paper screen: (a) the lens is placed to get a clear image;
(b) the lens is 2 mm closer to the camera than it is necessary for a clear
image.

###  III.2 Data acquisition

The data were recorded in series of 10000 frames with the frame rate of 24
fps. For each “excitation cycle”, each laser was randomly assigned its state
(active or dark) and the time delay between the “excitation pulse” and its
relaxation (“fluorescent emission”). The delay was sampled according to normal
distribution with the mean determined by the model temperature of the laser
(“hot” or “cold”). The random sampling was performed by the Arduino
microcontroller. Then, Arduino produced an “excitation pulse” and made each
active laser to flash after the time delay assigned to it. After modeling
“fluorescent emission” of all the active lasers withing the current cycle
(when all the lasers are switched off), Arduino waited for 1 frame, and
repeated the cycle. In the end, we got a video containing intensity
distribution in time and space.

###  III.3 Data processing and results

First, the recorded video was divided into frames with `ffmpeg` program [44]
and processed with Python script. The script calculated total intensity for
each frame and identified “excitation pulses” according to specific
alternation of high and low intensity of frames. The frames designating the
“excitation pulse” were excluded from further calculations. According to the
notations, introduced in Section II, the “excitation cycles” are indexed by
$t=1$, …, $N_{\text{f}}$.

The signal for each “excitation cycle” was split into two channels and
accumulated pixelwise. The “hot” signal $I_{0}(i,j;t)$ is the pixelwise sum of
intensities for the 1st - 7th frames after the $t$-th “excitation pulse”,
where indices $i$ and $j$ enumerate pixel positions in each frame. The “cold”
signal $I_{1}(i,j;t)$ is defined in the same way starting from the 8th frame
after the “excitation pulse” (Fig. 4).

Figure 4: Example of ”excitation cycle” timeline, where two pulses appear in
the “hot” time window and two pulses - in the “cold” time window.

The intensity (mean signal value) and the 2nd, 3rd, and 4th order cumulants
for the two channels where defined according to Eqs. (5)–(8). The obtained
intensity and cumulant images are shown in Fig. 5. For visualization purposes,
the “cold” and “hot” channels were assigned to the blue and red channels of
RGB palette respectively. Since the cumulant of the 3rd order can be negative,
the absolute value of the cumulant was used to create the image. One can see
that the intensity distribution (Fig. 5(a)) is blurred like original unfocused
image (Fig. 3(b)). However, it is possible to distinguish the light of “hot”
and “cold” lasers because of different relaxation times. The light spots are
more visible in the second-order cumulant image (Fig. 5(b)), but the
background caused by blurring is still strong. The third-order cumulant image
(Fig. 5(c)) is clear enough and the background is hardly visible. The best
quality image is given by the fourth-order cumulant (Fig. 5(d)). It allows one
to distinguish clearly spots of “hot” and “cold” lasers. To summarize, the
experiment is in a good agreement with the theory, and the calculation of
cumulants of intensity fluctuations allows one to visualize the particles with
different relaxation times.

Figure 5: The two-channel signal for intensity (a) and cumulants of the 2nd
(b), 3rd (c), and 4th (d) order.

###  III.4 Reconstruction of local “temperature”

To fit the above-described model experiment into the theoretical framework,
discussed in Section II, we assign the normalized dimensionless “temperature”
equal to $\theta=-1$ to “cold” lasers and $\theta=+1$ to “hot” lasers. The
emission profile variable $x$ corresponds to the camera frame number within
the “hyperframe” (excitation cycle) of SOFI in the considered case. For such
definitions, the model emission profile is described as

$\varepsilon(x,\theta)=\frac{1}{\sqrt{2\pi}\sigma}\exp\left(-\frac{(x-7.5+2.5\theta)^{2}}{2\sigma^{2}}\right)+b(x),$
(32)

where $\sigma=2$ and the additional parameter-independent background signal is
modeled as
$b(x)=\eta\exp\\{-(x-7.5)^{2}/(2\sigma_{\text{b}}^{2})\\}/(\sqrt{2\pi}\sigma_{\text{b}})$
with $\sigma_{\text{b}}=5$ and $\eta=0.1$.

A simple threshold filter is used to split the two detection channels:

$Y_{0}=\sum_{x\leq 7}\varepsilon(x,\theta),$ (33) $Y_{1}=\sum_{x\geq
8}\varepsilon(x,\theta).$ (34)

The decomposition coefficients $A_{ij}$ (Eq. (17)) can be approximated as

$A_{j0}=\frac{1+\eta}{2},\quad A_{j1}=\pm\frac{1}{\sqrt{2\pi}},$ (35)

where $j=0$ and 1 correspond to the “$+$” and “$-$” signs respectively. The
linear estimator $\bar{\theta}^{(n)}(\mathbf{r})$ is calculated according to
Eq. (24) with $\theta_{0}=0$, while its more accurate nonlinear counterpart is
constructed according to Eq. (26) using the procedure outlined before it (Fig.
6). Further, the reconstructed distribution can be refined using the method
outlined in the Section II.4 and Appendix D.

Figure 6: The nonlinear estimator
$\bar{\theta}_{\text{nonlin}}^{(n)}=\Theta^{-1}(\bar{\theta}^{(n)})$ for the
model experiment and analysis of intensity ($n=1$) and cumulants ($n=2$, 3,
4). The dashed diagonal line is shown for visual guidance.

The results are presented in Fig. 7. As predicted, usage of higher-order
cumulants enhances spatial resolutions and enables analysis of finer features
of the sample. According to the standard performance of SOFI, using fourth
order cumulant instead of intensity leads to approximately twice narrower
efficient PSF.

Figure 7: Reconstructed distribution of the model normalized temperature
$\theta$ for analysis of intensity (a), the second (b), the third (c) and the
fourth (d) order cumulants. Only the central part of the image is shown.

##  IV Discussion

To get quantitative estimates of the proposed approach performance, we applied
several sharpness measures to the reconstructed distribution of the local
“temperature”. The results shown in Fig. 8 indicate that the most blurred
distribution is obtained for the traditional intensity-based approach.
Derivative-based measures (Tenengrad and high-pass filtering [45]) give
preference to the second-order cumulant (after application of noise-reduction
algorithm from subsection II.4), while eigenvalue-based [46] and histogram
entropy based [45] methods assign the best score to the fourth-order cumulant.

Figure 8: Sharpness of the reconstructed distribution of the model normalized
temperature $\theta$ for different orders of cumulants used for the analysis.
Solid (dashed) lines correspond to the distributions after (before) the noise
correction described in subsection II.4.

We believe that the proposed approach and its proof-of-concept experiment can
pave the way for:

  * •

detailed analysis and optimization of the combined sensing and imaging
approach in terms of classical and quantum Fisher information;

  * •

reconstruction of the sample characteristic distribution by combined usage of
several orders of cumulants simultaneously;

  * •

extension of the approach to multimodal sensing, when several characterisics
of the sample are reconstructed simultaneously;

  * •

implementation of Fourier spectrum re-weighting for additional de-blurring and
cumulant orders combination;

  * •

application of the “sliding window” approach [47] for more efficient
processing of many-pixel images.

##  V Conclusions

We proposed a simple and elegant technique capable of combining sensing with
super-resolution imaging, based on non-trivial temporal statistics typical for
blinking fluorescent nanoemitters. The approach is universal since the
dependence of the fluorophores response on the characteristic to be
reconstructed can take different form: changes of the spectrum, fluorescence
lifetime, etc. It is shown that splitting the fluorescence signal into just
two detection channels at each spatial detection point is sufficient for
reliable reconstruction of the local characteristic distribution.

The conducted model experiment proofs soundness and efficiency of the approach
and can be used as a guidance for the real-world experiments with, for
example, living cells and structured fluorescent nanodiamonds [48] or blinking
graphene quantum dots ensembles [49]. In agreement with the theoretical
predictions, usage of cumulants instead of emission intensity enhances
sharpness of the reconstructed parameter distribution and enables analysis of
smaller features of the sample.

Proposed methodology will be useful for local measurements of living cells’
physical parameters, such as temperature, pH, electric field gradient,
presence of reactive oxygen spices, i.e. the characteristics which could be
recognized as spectral changes of appropriate fluorophore delivered targetly
to the wanted cell or even organelle. In its turn, this advanced local sensing
tool could give an important insights into fundamental life science, which
still lacks instruments capable of monitoring the cellular balance with
nanometre spatial resolution to guarantee measurements consistency for that
dynamic system.

###### Acknowledgements.

This study was accomplished through the financial support of Research Council
of Finland (Flagship Programme PREIN, decision 320166; research project,
decision 357033); QuantERA project EXTRASENS, in Finland funded through
Research Council of Finland, decision 361115), Horizon Europe MSCA FLORIN
Project 101086142.

##  Appendix A Optimal estimator for linear-response sensing

Here, we show that the information about the environment parameter $\theta$,
contained in the full signal density profile $\varepsilon(x,\theta)$,
described by Eq. (11), can be efficiently compressed into two scalar
characteristics $Y_{i}$ (Eq. (14)), which are used for building the estimator
$\bar{\theta}$ according to Eq. (12). The analysis is performed in the
vicinity of the reference value $\theta_{0}$ of the parameter, i.e. the
deviation $|\theta-\theta_{0}|$ is assumed to be small.

We assume that the signal profile is measured at a discrete set of bins
$\\{x\\}$ and contains shot noise satisfying Poisson statistics:

$\varepsilon(x,\theta)=S_{0}(x)+S_{1}(x)(\theta-\theta_{0})+\delta(x),$ (36)

where $\operatorname{E}[\delta(x)]=0$,
$\operatorname{E}[\delta(x)\delta(x^{\prime})]\approx
S_{0}(x)\delta_{xx^{\prime}}$. To reduce the influence of the noise caused by
the emitters’ blinking on the final estimates of the environment parameter, it
is useful to consider the normalized signal
$\varepsilon(x,\theta)/\sum_{x^{\prime}}\varepsilon(x^{\prime},\theta)$. The
Fisher information of that normalized density profile relatively to the
parameter $\theta$ at the point $\theta=\theta_{0}$ can be calculated as

$F_{\text{full}}(\theta)=\sum_{x}\frac{S_{1}^{2}(x)}{S_{0}(x)}-\frac{1}{\sum_{x}S_{0}(x)}\Bigl{(}\sum_{x}S_{1}(x)\Bigr{)}^{2},$
(37)

where the second term describes “information loss” caused by the normalization
[50].

For a low level of noise, the estimator $\bar{\theta}$ (constructed from the
noisy signal according to Eq. (12)) is a random variable obeying approximately
normal distribution with $\operatorname{E}[\bar{\theta}]=\theta$ and
$\operatorname{Var}[\bar{\theta}]\approx\sum_{x}Q_{1}^{2}(x)S_{0}(x)$. The
Fisher information after “compression” of the profile $\varepsilon(x,\theta)$
into the scalar $\bar{\theta}$ equals

$F_{\bar{\theta}}(\theta)=\Bigl{[}\sum_{x}Q_{1}^{2}(x)S_{0}(x)\Bigr{]}^{-1}.$
(38)

According to general properties of Fisher information [51, 52], the following
inequality holds:

$F_{\bar{\theta}}(\theta)\leq F_{\text{full}}(\theta).$ (39)

On the other hand, one can prove by direct calculation that the choice

$Q_{i}(x)=\frac{q_{i}(x)}{\sum_{x^{\prime}}q_{i}(x^{\prime})S_{i}(x^{\prime})},\quad
i=0,1,$ (40)

where

$q_{0}(x)=S_{0}(x)-S_{1}(x)\frac{\sum_{x^{\prime}}S_{0}(x^{\prime})S_{1}(x^{\prime})}{\sum_{x^{\prime}}S_{1}^{2}(x^{\prime})}$
(41)

and

$q_{1}(x)=\frac{S_{1}(x)}{S_{0}(x)}-\frac{\sum_{x^{\prime}}S_{1}(x^{\prime})}{\sum_{x^{\prime}}S_{0}(x^{\prime})}$
(42)

ensures that the maximal value of Fisher information is reached:
$F_{\bar{\theta}}(\theta)=F_{\text{full}}(\theta)$.

##  Appendix B Practical estimator for sensing with two spectral or temporal
windows

In this appendix, we show how scalar characteristics $Y_{i}$, defined by Eq.
(14) without any particular constraints imposed on the weights $Q_{j}(x)$, can
be modified so that Eq. (13) are satisfied.

In general case, the values $Y_{i}$ are decomposed according to Eq. (17). One
can construct the following linear combinations of those quantities:

$\bar{Y}_{0}=\frac{Y_{0}A_{11}-Y_{1}A_{01}}{A_{00}A_{11}-A_{10}A_{01}}\mapsto
1,$ (43)
$\bar{Y}_{1}=\frac{Y_{1}A_{00}-Y_{0}A_{10}}{A_{00}A_{11}-A_{10}A_{01}}\mapsto\theta-\theta_{0}.$
(44)

Linearity of Eq. (17) implies that the weights $\bar{Q}_{i}(x)$, corresponding
to the quantities $\bar{Y}_{j}$ are related to the initial weights $Q_{i}(x)$
by the same relations:

$\bar{Q}_{0}(x)=\frac{Q_{0}(x)A_{11}-Q_{1}(x)A_{01}}{A_{00}A_{11}-A_{10}A_{01}},$
(45)
$\bar{Q}_{1}(x)=\frac{Q_{1}(x)A_{00}-Q_{0}(x)A_{10}}{A_{00}A_{11}-A_{10}A_{01}}.$
(46)

One can check by direct substitution that the conditions (13) are satisfied
for the constructed weights $\bar{Q}_{i}(x)$.

Fig. 9 shows an example of weights calculated for the model signal densities
depicted in Fig. 1(a,b).

Figure 9: Weights for sensing by fluorescence spectrum (a,c,e) and lifetime
(b,d,f) for the model signal densities depicted in Fig. 1(a,b): optimal
weights (a,b) defined by Eq. (40); practical weights (c,d) corresponding to
Gaussian spectral filters and threshold-based treatment of fluorescence
lifetime; optimal linear combinations (e,f) of practical weights accoring to
Eqs. (45) and (46).

##  Appendix C General properties of cumulant-based estimator for noiseless
data

To be physically reasonable, the estimator
$\bar{\theta}_{\text{nonlin}}(\mathbf{r})$, defined by Eq. (26) should satisfy
several general requirements when applied to noiseless data.

1\. For spatially uniform distribution $\theta(\mathbf{r})=const=\theta$ the
estimator should reconstruct the correct value:
$\bar{\theta}_{\text{nonlin}}(\mathbf{r})=\theta$ for all $\mathbf{r}$. In
that case, according to Eq. (22), one has
$C_{1}^{(n)}(\mathbf{r})/C_{0}^{(n)}(\mathbf{r})=Y_{1}^{n}/Y_{0}^{n}$ and,
therefore, Eq. (24) retains its value after the replacement
$C_{j}^{(n)}(\mathbf{r})\mapsto Y_{j}^{n}$. Then, according to the definition
(26), the equality $\bar{\theta}_{\text{nonlin}}(\mathbf{r})=\theta$ holds for
all $\mathbf{r}$.

2\. If the parameter values $\theta_{k}$ in the vicinity of the emitters are
bounded as $\theta_{\text{min}}\leq\theta_{k}\leq\theta_{\text{max}}$, the
estimator value is also bounded:
$\theta_{\text{min}}\leq\bar{\theta}_{\text{nonlin}}(\mathbf{r})\leq\theta_{\text{max}}$
for all $\mathbf{r}$. According to Eq. (15), the estimator
$\bar{\theta}=\Theta(\theta)$ is a function of the ratio $Y_{1}/Y_{0}$, which
itself is a function of $\theta$. The monotonicity of the dependence
$\Theta(\theta)$ implies that both of the above-mentioned functions are also
monotonic. Therefore, for each emitter, the ratio $Y_{k1}/Y_{k0}$ is bounded
by its values at $\theta_{k}=\theta_{\text{min}}$ and
$\theta_{k}=\theta_{\text{max}}$. The same conclusion is valid for the ratio
$Y_{k1}^{n}/Y_{k0}^{n}$, since $n$ is positive.

Assuming that the blinking statistics is the same for all emitters (i.e.,
$C_{t}^{(n)}[s_{k}(t)]$ depends only on $n$, but not on $k$) and taking into
account the definition (22), one can conclude that the ratio
$R(\mathbf{r})=C_{1}^{(n)}(\mathbf{r})/C_{0}^{(n)}(\mathbf{r})$ is bounded by
its values for uniform distributions with $\theta_{k}=\theta_{\text{min}}$ and
$\theta_{k}=\theta_{\text{max}}$ for all $k$, too. The derivative of
$\bar{\theta}^{(n)}$, defined by Eq. (24), over the ratio $R$ has constant
sign, independent of $R$:

$\frac{\partial\bar{\theta}^{(n)}(\mathbf{r})}{\partial
R(\mathbf{r})}=\frac{A_{00}}{nA_{01}}\frac{(A_{10}/A_{00})^{n-1}A_{11}/A_{01}-(A_{10}/A_{00})^{n}}{[R(\mathbf{r})-(A_{10}/A_{00})^{n-1}A_{11}/A_{01}]^{2}}.$
(47)

Therefore, the values $\bar{\theta}^{(n)}(\mathbf{r})$ are also bounded:

$\Theta(\theta_{\text{min}})\leq\bar{\theta}^{(n)}(\mathbf{r})\leq\Theta(\theta_{\text{max}}).$
(48)

Finally, taking into account the definition (26) and monotonicity of
$\Theta^{-1}$, we have

$\theta_{\text{min}}\leq\bar{\theta}^{(n)}_{\text{nonlin}}(\mathbf{r})\leq\theta_{\text{max}}.$
(49)

3\. If the emitters are evenly distributed (i.e., for each point of interest
$\mathbf{r}$ there is at least one emitter with
$U(\mathbf{r}-\mathbf{r}_{k})\varepsilon_{k}\geq C>0$ for some constant $C$),
the parameter estimation by $\bar{\theta}^{(n)}_{\text{nonlin}}(\mathbf{r})$
is local:

$\frac{\partial\bar{\theta}^{(n)}_{\text{nonlin}}(\mathbf{r})}{\partial\theta_{k}}\rightarrow
0\text{ for }|\mathbf{r}-\mathbf{r}_{k}|\rightarrow\infty.$ (50)

The statement is valid if the PSF of the imaging system is local:

$U(|\mathbf{r}-\mathbf{r}_{k}|)\rightarrow 0\text{ for
}|\mathbf{r}-\mathbf{r}_{k}|\rightarrow\infty.$ (51)

Then, the relative contribution of a distant emitter to the cumulant value
becomes negligible for $|\mathbf{r}-\mathbf{r}_{k}|\rightarrow\infty$:

$\left|\frac{1}{C_{j}^{(n)}(\mathbf{r})}\frac{\partial
C_{j}^{(n)}(\mathbf{r})}{\partial\theta_{k}}\right|\propto\frac{U^{n}(|\mathbf{r}-\mathbf{r}_{k}|)\varepsilon_{k}^{n}}{C^{n}}\rightarrow
0.$ (52)

Therefore, it does not influence the value of the estimator
$\bar{\theta}^{(n)}_{\text{nonlin}}(\mathbf{r})$.

##  Appendix D Expressions for regularized parameter estimation

For parameter values, defined at a regular grid $\mathbf{r}_{ij}$, the
derivatives in Eq. (29) can be approximated as finite differences. Then, the
expression for log-probability takes the form

$-\log
p\left[\theta(\mathbf{r}_{ij})|\bar{\theta}(\mathbf{r}_{ij})\right]\propto\sum_{i,j}\frac{[\bar{\theta}(\mathbf{r}_{ij})-\theta(\mathbf{r}_{ij})]^{2}}{2\Delta\bar{\theta}^{2}(\mathbf{r}_{ij})}\\\
+\sum_{i,j}\frac{[2\theta(\mathbf{r}_{ij})-\theta(\mathbf{r}_{i-1,j})-\theta(\mathbf{r}_{i+1,j})]^{2}}{2D^{2}}\\\
+\sum_{i,j}\frac{[2\theta(\mathbf{r}_{ij})-\theta(\mathbf{r}_{i,j-1})-\theta(\mathbf{r}_{i,j+1})]^{2}}{2D^{2}}.$
(53)

The sums representing the second and the third terms do not include boundary
values of $i$ and $j$ respectively. The parameter $D$ scales as
$D=const\times\langle\Delta\bar{\theta}(\mathbf{r}_{ij})\rangle w$, where $w$
is the PSF width expressed in the grid steps and representing the expected
minimal size of the reconstructed distribution features, which are not
artefacts yet. The proportionality coefficient can be calculated from detailed
theoretical investigation of the noise models, but, in practice, it can be
just chosen empirically.

Unconstrained minimization of the quantity, defined by the right-hand side of
Eq. (53), is trivial and includes just linear inversion according to the
least-squares method. One of the approaches to obtaining physically reasonable
constrained estimates
$\theta_{\text{min}}\leq\theta(\mathbf{r}_{ij})\leq\theta_{\text{max}}$ is to
perform constrained quadratic optimization, which, however, becomes much more
computationally hard for a large number of grid points $\mathbf{r}_{ij}$. A
more efficient way to get reasonable results is to constrain the initial
estimate $\bar{\theta}(\textbf{r}_{ij})$ to the range
$[\theta_{\text{min}},\theta_{\text{max}}]$ before applying the minimization
of the expression (53). If all the values satisfy the constraint
$\bar{\theta}(\textbf{r}_{ij})\in[\theta_{\text{min}},\theta_{\text{max}}]$,
the minimizing set $\\{\theta(\textbf{r}_{ij}\\}$ will also belong to the
physical range with certainty.

## References

  * Norman et al. [2021] Victoria A Norman, Sridhar Majety, Zhipan Wang, William H Casey, Nicholas Curro, and Marina Radulaski.  Novel color center platforms enabling fundamental scientific discovery.  _InfoMat_ , 3(8):869–890, 2021.
  * Wu and Weil [2022] Yingke Wu and Tanja Weil.  Recent developments of nanodiamond quantum sensors for biological applications.  _Advanced Science_ , 9(19):2200059, 2022.
  * Radtke et al. [2019] Mariusz Radtke, Ettore Bernardi, Abdallah Slablab, Richard Nelz, and Elke Neu.  Nanoscale sensing based on nitrogen vacancy centers in single crystal diamond and nanodiamonds: Achievements and challenges.  _Nano Futures_ , 3(4):042004, 2019.
  * Xu et al. [2023] Ying Xu, Weiye Zhang, and Chuanshan Tian.  Recent advances on applications of nv- magnetometry in condensed matter physics.  _Photonics Research_ , 11(3):393–412, 2023.
  * Baranov et al. [2011] Pavel G Baranov, Anna P Bundakova, Alexandra A Soltamova, Sergei B Orlinskii, Igor V Borovykh, Rob Zondervan, Rogier Verberk, and Jan Schmidt.  Silicon vacancy in sic as a promising quantum system for single-defect and single-photon spectroscopy.  _Physical Review B_ , 83(12):125203, 2011.
  * Anisimov et al. [2016] AN Anisimov, Dmitrij Simin, VA Soltamov, SP Lebedev, PG Baranov, GV Astakhov, and V Dyakonov.  Optical thermometry based on level anticrossing in silicon carbide.  _Scientific Reports_ , 6(1):33301, 2016.
  * Fu et al. [2007] Chi-Cheng Fu, Hsu-Yang Lee, Kowa Chen, Tsong-Shin Lim, Hsiao-Yun Wu, Po-Keng Lin, Pei-Kuen Wei, Pei-Hsi Tsao, Huan-Cheng Chang, and Wunshain Fann.  Characterization and application of single fluorescent nanodiamonds as cellular biomarkers.  _Proceedings of the National Academy of Sciences_ , 104(3):727–732, 2007.
  * Hui et al. [2010] Yuen Yung Hui, Bailin Zhang, Yuan-Chang Chang, Cheng-Chun Chang, Huan-Cheng Chang, Jui-Hung Hsu, Karen Chang, and Fu-Hsiung Chang.  Two-photon fluorescence correlation spectroscopy of lipid-encapsulated fluorescent nanodiamonds in living cells.  _Optics express_ , 18(6):5896–5905, 2010.
  * Blakley et al. [2015] SM Blakley, IV Fedotov, S Ya Kilin, and AM Zheltikov.  Room-temperature magnetic gradiometry with fiber-coupled nitrogen-vacancy centers in diamond.  _Optics letters_ , 40(16):3727–3730, 2015.
  * Blakley et al. [2016] SM Blakley, IV Fedotov, LV Amitonova, EE Serebryannikov, H Perez, S Ya Kilin, and AM Zheltikov.  Fiber-optic vectorial magnetic-field gradiometry by a spatiotemporal differential optical detection of magnetic resonance in nitrogen–vacancy centers in diamond.  _Optics Letters_ , 41(9):2057–2060, 2016.
  * Chang et al. [2017] K Chang, Alexander Eichler, Jan Rhensius, Luca Lorenzelli, and Christian L Degen.  Nanoscale imaging of current density with a single-spin magnetometer.  _Nano letters_ , 17(4):2367–2373, 2017.
  * Pelliccione et al. [2016] Matthew Pelliccione, Alec Jenkins, Preeti Ovartchaiyapong, Christopher Reetz, Eve Emmanouilidou, Ni Ni, and Ania C Bleszynski Jayich.  Scanned probe imaging of nanoscale magnetism at cryogenic temperatures with a single-spin quantum sensor.  _Nature nanotechnology_ , 11(8):700–705, 2016.
  * Malykhin et al. [2022] Sergei Malykhin, Yuliya Mindarava, Rinat Ismagilov, Fedor Jelezko, and Alexander Obraztsov.  Control of nv, siv and gev centers formation in single crystal diamond needles.  _Diamond and Related Materials_ , 125:109007, 2022.
  * Hell and Wichmann [1994] Stefan W Hell and Jan Wichmann.  Breaking the diffraction resolution limit by stimulated emission: stimulated-emission-depletion fluorescence microscopy.  _Optics letters_ , 19(11):780–782, 1994.
  * Rittweger et al. [2009a] Eva Rittweger, Dominik Wildanger, and SW Hell.  Far-field fluorescence nanoscopy of diamond color centers by ground state depletion.  _Europhysics Letters_ , 86(1):14001, 2009a.
  * Rittweger et al. [2009b] Eva Rittweger, Kyu Young Han, Scott E Irvine, Christian Eggeling, and Stefan W Hell.  STED microscopy reveals crystal colour centres with nanometric resolution.  _Nature Photonics_ , 3(3):144–147, 2009b.
  * Maurer et al. [2010] PC Maurer, Jeronimo R Maze, PL Stanwix, Liang Jiang, Alexey Vyacheslavovich Gorshkov, Alexander A Zibrov, Benjamin Harke, JS Hodges, Alexander S Zibrov, Amir Yacoby, et al.  Far-field optical imaging and manipulation of individual spins with nanoscale resolution.  _Nature Physics_ , 6(11):912–918, 2010.
  * Classen et al. [2017] Anton Classen, Joachim von Zanthier, Marlan O Scully, and Girish S Agarwal.  Superresolution via structured illumination quantum correlation microscopy.  _Optica_ , 4(6):580–587, 2017.
  * Hess et al. [2006] Samuel T Hess, Thanu PK Girirajan, and Michael D Mason.  Ultra-high resolution imaging by fluorescence photoactivation localization microscopy.  _Biophysical journal_ , 91(11):4258–4272, 2006.
  * Rust et al. [2006] Michael J Rust, Mark Bates, and Xiaowei Zhuang.  Sub-diffraction-limit imaging by stochastic optical reconstruction microscopy (storm).  _Nature methods_ , 3(10):793–796, 2006.
  * Schnitzbauer et al. [2017] Joerg Schnitzbauer, Maximilian T Strauss, Thomas Schlichthaerle, Florian Schueder, and Ralf Jungmann.  Super-resolution microscopy with dna-paint.  _Nature protocols_ , 12(6):1198–1228, 2017.
  * Gu et al. [2013] Min Gu, Yaoyu Cao, Stefania Castelletto, Betty Kouskousis, and Xiangping Li.  Super-resolving single nitrogen vacancy centers within single nanodiamonds using a localization microscope.  _Optics express_ , 21(15):17639–17646, 2013.
  * Dertinger et al. [2009] Thomas Dertinger, Ryan Colyer, Gopal Iyer, Shimon Weiss, and Jörg Enderlein.  Fast, background-free, 3D super-resolution optical fluctuation imaging (SOFI).  _Proceedings of the National Academy of Sciences_ , 106(52):22287–22292, 2009.
  * Dertinger et al. [2010] Thomas Dertinger, Ryan Colyer, Robert Vogel, Jörg Enderlein, and Shimon Weiss.  Achieving increased resolution and more pixels with superresolution optical fluctuation imaging (SOFI).  _Optics express_ , 18(18):18875–18885, 2010.
  * Dertinger et al. [2013] Thomas Dertinger, Alessia Pallaoro, Gary Braun, Sonny Ly, Ted A Laurence, and Shimon Weiss.  Advances in superresolution optical fluctuation imaging (SOFI).  _Quarterly reviews of biophysics_ , 46(2):210–221, 2013.
  * Schwartz and Oron [2012] Osip Schwartz and Dan Oron.  Improved resolution in fluorescence microscopy using quantum correlations.  _Physical Review A_ , 85(3):033812, 2012.
  * Schwartz et al. [2013] Osip Schwartz, Jonathan M Levitt, Ron Tenne, Stella Itzhakov, Zvicka Deutsch, and Dan Oron.  Superresolution microscopy with quantum emitters.  _Nano Letters_ , 13(12):5832–5836, 2013.
  * Monticone et al. [2014] D Gatto Monticone, K Katamadze, Paolo Traina, E Moreva, Jacopo Forneris, Ivano Ruo-Berchera, Paolo Olivero, IP Degiovanni, Giorgio Brida, and M Genovese.  Beating the abbe diffraction limit in confocal microscopy via nonclassical photon statistics.  _Physical Review Letters_ , 113(14):143602, 2014.
  * Gallina et al. [2013] Maria Elena Gallina, Jianmin Xu, Thomas Dertinger, Adva Aizer, Yaron Shav-Tal, and Shimon Weiss.  Resolving the spatial relationship between intracellular components by dual color super resolution optical fluctuations imaging (sofi).  _Optical Nanoscopy_ , 2(1):1–9, 2013.
  * Grußmayer et al. [2020] Kristin S Grußmayer, Stefan Geissbuehler, Adrien Descloux, Tomas Lukes, Marcel Leutenegger, Aleksandra Radenovic, and Theo Lasser.  Spectral cross-cumulants for multicolor super-resolved sofi imaging.  _Nature Communications_ , 11(1):3023, 2020.
  * Liu et al. [2020] Zhihe Liu, Jie Liu, Zhe Zhang, Zezhou Sun, Xuyang Shao, Jintong Guo, Lei Xi, Zhen Yuan, Xuanjun Zhang, Daniel T Chiu, et al.  Narrow-band polymer dots with pronounced fluorescence fluctuations for dual-color super-resolution imaging.  _Nanoscale_ , 12(14):7522–7526, 2020.
  * Glogger et al. [2021] Marius Glogger, Christoph Spahn, Jörg Enderlein, and Mike Heilemann.  Multi-color, bleaching-resistant super-resolution optical fluctuation imaging with oligonucleotide-based exchangeable fluorophores.  _Angewandte Chemie_ , 133(12):6380–6383, 2021.
  * Horoshko et al. [2024] D. Horoshko, M. Kolobov, A. Mikhalychev, and P. Kuzhir.  Fisher information analysis for localisation of blinking sources of light.  To be published, 2024.
  * Bushuev et al. [2016] EV Bushuev, V Yu Yurov, AP Bolshakov, VG Ralchenko, EE Ashkinazi, AV Ryabova, IA Antonova, PV Volkov, AV Goryunov, and A Yu Luk’yanov.  Synthesis of single crystal diamond by microwave plasma assisted chemical vapor deposition with in situ low-coherence interferometric control of growth rate.  _Diamond and Related Materials_ , 66:83–89, 2016.
  * Maletinsky et al. [2012] Patrick Maletinsky, Sungkun Hong, Michael Sean Grinolds, Birgit Hausmann, Mikhail D Lukin, Ronald L Walsworth, Marko Loncar, and Amir Yacoby.  A robust scanning diamond sensor for nanoscale imaging with single nitrogen-vacancy centres.  _Nature nanotechnology_ , 7(5):320–324, 2012.
  * Arcizet et al. [2011] Olivier Arcizet, Vincent Jacques, Alessandro Siria, Philippe Poncharal, Pascal Vincent, and Signe Seidelin.  A single nitrogen-vacancy defect coupled to a nanomechanical oscillator.  _Nature Physics_ , 7(11):879–883, 2011.
  * Babinec et al. [2010] Thomas M Babinec, Birgit JM Hausmann, Mughees Khan, Yinan Zhang, Jeronimo R Maze, Philip R Hemmer, and Marko Lončar.  A diamond nanowire single-photon source.  _Nature nanotechnology_ , 5(3):195–199, 2010.
  * Lesik et al. [2013] Margarita Lesik, Piernicola Spinicelli, Sébastien Pezzagna, Patrick Happel, Vincent Jacques, Olivier Salord, Bernard Rasser, Anne Delobbe, Pierre Sudraud, Alexandre Tallaire, et al.  Maskless and targeted creation of arrays of colour centres in diamond using focused ion beam technology.  _physica status solidi (a)_ , 210(10):2055–2059, 2013.
  * Golubewa et al. [2022] Lena Golubewa, Yaraslau Padrez, Sergei Malykhin, Tatsiana Kulahava, Ekaterina Shamova, Igor Timoshchenko, Marius Franckevicius, Algirdas Selskis, Renata Karpicz, Alexander Obraztsov, et al.  All-optical thermometry with NV and SiV color centers in biocompatible diamond microneedles.  _Advanced Optical Materials_ , 10(15):2200631, 2022.
  * Tikhonov and Arsenin [1977] Andrey N Tikhonov and Vasiliy Y Arsenin.  _Solutions of ill-posed problems_.  Washington, DC: John Wiley & Sons, New York, 1977.
  * Gelman et al. [2013] A. Gelman, J. B. Carlin, H. S. Stern, and D. B. Dunson.  _Bayesian data analysis_.  Chapman and Hall/CRC, 2013.
  * Koch [2007] K.-R. Koch.  _Introduction to Bayesian statistics_.  Springer Science & Business Media, 2007.
  * [43] Arduino hardware and software.  https://www.arduino.cc.
  * [44] A complete, cross-platform solution to record, convert and stream audio and video.  https://ffmpeg.org.
  * Krotkov [1988] Eric Krotkov.  Focusing.  _International Journal of Computer Vision_ , 1(3):223–237, 1988.
  * Wee and Paramesran [2007] Chong-Yaw Wee and Raveendran Paramesran.  Measure of image sharpness using eigenvalues.  _Information Sciences_ , 177(12):2533–2552, 2007.
  * Mikhalychev et al. [2019] AB Mikhalychev, B Bessire, IL Karuseichyk, AA Sakovich, M Unternährer, DA Lyakhov, Dominik L Michels, André Stefanov, and D Mogilevtsev.  Efficiently reconstructing compound objects by quantum imaging with higher-order correlation functions.  _Communications Physics_ , 2(1):134, 2019.
  * Quarshie et al. [2024] Mariam M Quarshie, Sergei Malykhin, Alexander Obraztsov, and Polina Kuzhir.  Nano-and micro-crystalline diamond film structuring with electron beam lithography mask.  _Nanotechnology_ , 35(15):155301, 2024.
  * Belko et al. [2022] Nikita Belko, Lena Golubewa, Vyacheslav Chizhevsky, Sopfy Karuseichyk, Dmitry Filimonenko, Marija Jankunec, Hamza Rehman, Tatsiana Kulahava, Polina Kuzhir, and Dmitri Mogilevtsev.  Hysteresis and stochastic fluorescence by aggregated ensembles of graphene quantum dots.  _The Journal of Physical Chemistry C_ , 126(25):10469–10477, 2022.
  * Mikhalychev et al. [2021] Alexander Mikhalychev, Konstantin Zhevno, Svetlana Vlasenko, Andrei Benediktovitch, Tatjana Ulyanenkova, and Alex Ulyanenkov.  Fisher information for optimal planning of x-ray diffraction experiments.  _Journal of Applied Crystallography_ , 54(6):1676–1697, 2021.
  * Lehmann and Casella [2006] Erich L Lehmann and George Casella.  _Theory of point estimation_.  Springer, 2006.
  * Cohen [1968] M Cohen.  The Fisher information and convexity.  _IEEE Trans. Inf. Theory_ , 14(4):591–592, 1968.