Super-resolution optical fluctuation imaging—fundamental estimation theory perspective

The SOFI method is based on calculating temporal cumulants of measured intensity distribution in a number of time frames. It's often claimed, that the resolution can be increased by a factor $\sqrt{k}$ if the kth cumulant is computed, which is justified as follows. If the imaged sample consists of L independently fluctuating point emitters, then the light intensity observed in the image plane is:

$$\begin{align} I( \vec r, t) = \sum_{i = 1}^L P_i(t) U(\vec r - \vec{r_i}), \end{align} \tag{ 1 }$$

where the stochastic process P_i (t) represents the fluctuating brightness of ith emitter, $\vec{r_i}$ its position in the image plane, and $U(\vec r)$ is the PSF of the system. Since the emitters are independent, kth temporal cumulant of the signal (at a given $\vec r$) reads:

$$\begin{align} \kappa_k(\vec r) = \sum_{i = 1}^L \kappa_k [P_i(t)] U^k(\vec r - \vec {r_i}), \end{align} \tag{ 2 }$$

where $\kappa_k[P_i(t)]$ is the kth cumulant of the stochastic process P_i (t). The PSF is now replaced by its kth power. If the standard, Gaussian approximation of the PSF is used, U^k is narrowed by a factor $\sqrt k$ compared with U.

The described scheme can be improved in various ways. The most important modification is based on utilizing spatial correlations (cross-cumulants) in an image reconstruction algorithm [28]. Cross-cumulants of kth order are computed for all k-element subsets of all pixels, and the cumulant computed for a given set gives rise to the signal located in its centroid. Cross-cumulants corresponding to the same centroid can be summed directly, or with proper weights in order to maximize the signal-to-noise ratio (SNR) [29]. This approach not only allows to make use of information hidden in the spatial correlation, but also increases the number of pixels in the final image, which is significant from the practical point of view. The latter effect can be also achieved if recorded images are Fourier-interpolated before further processing [30].

Unfortunately, even after applying the described modifications, higher cumulants are more noisy, and it is not possible to achieve the unlimited resolution gain in practice—this effect reappears in all known experiments. Therefore, noise has to be taken into account in order to assess the maximal resolution gain achievable in SOFI. Some analysis of the impact of noise on the computed cumulants estimators have been made [31–33], but these studies have not employed estimation theory concepts such as the Fisher information (FI), and did not make an attempt to benchmark the performance of the methods against the fundamental limitations imposed by estimation theory. Our goal is to provide such a rigorous study.

Let's consider the simplest, yet representative case of imaging a binary object, which consists of two identical point emitters with fluctuating brightness. Those two emitters are assumed to lie on a known axis, perpendicular to the optical axis of the imaging system, so the problem becomes 1D, and only transverse resolution is studied. Moreover, we assume that the centroid of the object is also known, and only the distance between emitters (θ) needs to be estimated, see figure 1 (note that in case of emitters of different brightness all the reasoning will be basically unaltered provided one replaces the geometric distance between the sources by a quantity based on the second moment of intensity distribution, as this is the quantity that is subject to the Rayleigh curse [10, 24, 34]). Given a random vector $\boldsymbol{N}$ that represents the data, distributed according to a probability distribution which is a function of the estimated parameter $p_\theta(\boldsymbol{N})$, the variance $\textrm{Var} [ \tilde \theta]$ of any locally unbiased estimator of θ is lower bounded by $(\mathcal F_\textrm{meas})^{-1}$, where

$$\begin{align} \mathcal F_\textrm{meas}(\theta) = \int \frac{1}{p_\theta(\boldsymbol{N})} \left(\frac{ \partial p_\theta(\boldsymbol{N})}{\partial \theta } \right)^2 \, \textrm{d} \boldsymbol{N} \end{align} \tag{ 3 }$$

is the FI associated with the whole measurement [35]. For the purpose of comparing different strategies we will use the FI per photon $\mathcal{F}(\theta) = \mathcal{F}_\textrm{meas}(\theta)/\bar{N}$, where $\bar{N}$ is the mean number of photons involved in the experiment.

Zoom In Zoom Out Reset image size

From now on, we assume that the PSF is Gaussian with standard deviation σ. The FI per one photon for standard imaging (SI) of Poissonian sources with constant brightness [10, 23], $\mathcal F^\textrm{(SI)}$, is sketched as a function of θ in figure 2. A significant drop in estimation precision below the Rayleigh limit is visible. In super-resolution microscopy we are mostly interested in the sub-Rayleigh regime, i.e. we assume that $\theta \ll \sigma$. If no noise apart from shot noise is present, and the effect of finite spatial resolution of the detector is neglected, $\mathcal F^\textrm{(SI)}$ for small θ can be approximated as (see appendix A.1 for details)

$$\begin{align} \mathcal F^\textrm{(SI)}(\theta) = \mathcal F^\textrm{(M)}(\theta) = \theta^2 /8\sigma^4 + \mathcal O(\theta^4/\sigma^6), \end{align} \tag{ 4 }$$

where we have also indicated that SI is equivalent to the analysis based only on the mean of the total number of photons (M) collected over the whole duration of the experiment. Now, if the PSF is narrowed by a factor s, the FI in the limit $\theta \rightarrow 0$ increases by a factor s⁴—this observation allows us to connect PSF-size approach with an estimation theory approach. If a given super-resolution imaging scheme leads to an increase of FI from $\mathcal F^\textrm{(SI)}(\theta)$ to $\mathcal F(\theta)$, then for $\theta \ll \sigma$ this change is equivalent to narrowing of the PSF by a factor

$$\begin{align} \zeta = \lim_{\theta \rightarrow 0} \left(\mathcal{F}(\theta) /\mathcal{F}^\textrm{(SI)}(\theta) \right)^{1/4}. \end{align} \tag{ 5 }$$

The factor ζ will be called the resolution gain limit (RGL) and will serve us as a figure of merit to asses the performance of different super-resolution methods.

Zoom In Zoom Out Reset image size

Admittedly, it does not contain all the information on the performance of a given super-resolving technique, as it ignores the behaviour for larger θ (see figure 2), and the behaviour of the technique for more complicated objects. Nevertheless, this quantity captures in a simple way the essence of the super-resolving potential in a basic two point sources model, and allows to compare different methods in a well-defined way. Moreover, this quantity also provides us with a meaningful upper bound on the performance of a method in more complex imaging scenarios, as discrimination of binary objects is a prerequisite for resolving multiple sources. If a given super-resolving method yields a multiple-source image that can be equivalently regarded as the one obtained with a SI system with an appropriately narrowed PSF, then clearly when focusing on any two neighbouring points we can bound the performance of this method by our limit based on the simplified two-point imaging model. This limit may not be tight, but interestingly, even such an optimistic limit turns out to be lower than the resolution gain predicted by the naive PSF size analysis in some cases. It is worth pointing out, that in all cases examined in this work, the ratio $\left(\mathcal{F}(\theta) /\mathcal{F}^\textrm{(SI)}(\theta) \right)^{1/4}$ decreases with increasing θ. This implies that ζ remains a valid bound on the resolution gain enhancement irrespectively of the imaged points separation, and hence can be interpreted as a factor that reveals the maximal potential enhancement of a super-resolving method considered. Finally, maximization of ζ in a given protocol can be regarded as a rule of thumb prescription on the choice of parameters that is likely to lead to the optimal performance of the protocol in real-life scenarios.

In the definition of ζ (5), the same pixel size is used to calculate $\mathcal F$ and $\mathcal F^\textrm{(SI)}$. In order to examine the impact of the pixel size for different methods, we will use a slightly modified figure of merit defined as:

$$\begin{align} \zeta^\textrm{(pix)} = \lim_{\theta \rightarrow 0} \left(\mathcal{F}(\theta) /\left(\frac{\theta^2}{8~\sigma^4} \right) \right)^{1/4}. \end{align} \tag{ 18 }$$

This modification fixes the denominator to $\mathcal F^\textrm{(SI)}$ associated with the infinite spatial resolution of the detector. It allows us to observe, for example, how the SI resolution decreases when Δx is too large. The role of the pixel size becomes less trivial when auto-cumulants are used in the estimation (see figure 6). Pixels cannot be of course too large, but very small pixels are no longer the optimal choice because the information contained in the correlations between pixels is lost. In particular, higher auto-cumulants do not provide any extra information if one photon per frame per pixel is detected at most. The described problem disappears when cross-cumulants are used, and very small pixels again become advantageous. Notice the fact, that a similar trade-off is observed in the time domain, when one changes the frame time. Very short time frames would only be optimal if correlations between frames were used, but such schemes are not analyzed in this work.

Zoom In Zoom Out Reset image size

Even though the shot noise, resulting from the finite detection statistics, is the most fundamental one, other types of noise (e.g. camera noise) often play a significant role in SOFI. Let us study a simple model, in which background noise is Poissonian, uncorrelated, and its mean value µ_B is the same in each pixel and in each frame. In order to take such a noise into account, it is enough to repeat the whole reasoning described in section 3 with only one modification—the term associated with the background noise should be added to each random variable n_i,m (see appendix A.3 for more details).

Non-zero value of µ_B leads to ζ decrease in all examined cases (where SI without background noise is treated as a reference when ζ is computed). It turns out that cumulants based methods are more robust against this type of noise than the SI. As we show in figure 7 , the relative decrease of ζ due to the background noise is the highest, when only mean signal is used in the estimation. This result was to be expected, as the authors of SOFI method claim, that the background noise reduction is its important feature [6]. Nevertheless, the background noise affects the performance of each of the analyzed methods.

Zoom In Zoom Out Reset image size

So far, we have presented results for $\tau_\textrm{on} = \tau_\textrm{off}$ (realistic model) or p = 1/2 (simplified model). However, it is known that real emitters sometimes break this assumption, and favour one of the states. One can observe, that if off- state is more probable, the RGL becomes higher because more photons are emitted within frames in which two sources have different brightness. One should also take this asymmetry into account while dealing with optimizing the time frame for different estimation schemes—see figure 8.

Zoom In Zoom Out Reset image size

Let us examine the performance of an improved version of the image reconstruction scheme based on summing of the covariances of pixel pairs with the same centroid (M + XC2s). The authors of [29] propose to add weights to the summation procedure in order to maximize the SNR of the reconstructed image. Let $\kappa_1, \kappa_2,\ldots,\kappa_k$ be the covariances contributing to the same centroid. The signal located in this centroid is equal to $S = \kappa_1 + \kappa_2 + \cdots +\kappa_k$, when the basic scheme is considered. More generally, one can use the formula $S = w_1~\kappa_1 + w_2~\kappa_2 +\cdots+w_k \kappa_k$, where w_i are some arbitrary weights. As shown in [29], in order to maximize the $\textrm{SNR} = \frac{\left\langle S \right\rangle }{\sqrt{\textrm{Var}(S)}}$, one should choose w_i satisfying the following conditions:

$$\begin{align} \sum_{i = 1}^k w_i \left( \frac{A_{mi}}{\kappa_m} - \frac{A_{1i}}{\kappa_1}\right) = 0, \end{align} \tag{ 19 }$$

where m ∈ {2, 3, ..., k}, and $A_{ij} = \textrm{cov}(\kappa_i, \kappa_j)$. In reality, one obviously has no access to the exact values of κ_i and A_ij , so their estimators must be inserted into (19) in order to obtain the optimal weights. However, in the limit of large number of collected frames, the difference between the actual values and the estimators does not affect the weights w_i significantly. We will therefore assume, that the exact values of κ_i and A_ij are used in the weights computing procedure.

We can assume that w₁ = 1 without the loss of generality, as multiplying each weight by the same factor does not change the information content of the computed linear combination (in particular, the SNR remains the same). It's then straightforward to compute the rest of weights, by solving the set of equations (19). Then, the procedure for calculating $\boldsymbol{\mu}$, $\boldsymbol{\Sigma}$ , and afterwards $\mathcal F$ and ζ, is the same as the one described in appendix A.3—one only needs to replace the sum of the covariances with the proper linear combination.

As we show in figure 9, the improved procedure (M + XC2w) allows only a slight increase of ζ compared to the previously examined summation scheme (M + XC2s). The limit dictated by the approach in which covariances are treated as independent observables (M + XC2) is still far from being achieved. We therefore claim, that the problem of finding the optimal way to reconstruct the image using covariances, still remains open. It's partially due to the fact, that the SNR itself does not provide a full description of noise, when noise is correlated, as in the analyzed schemes.

Zoom In Zoom Out Reset image size

The photon number fluctuations utilized in SOFI are described by super-Poissonian distribution, and therefore can be explained with the help of semi-classical models. However, sub-Poissonian photon number distributions can be used to obtain super-resolution as well. As it is believed (and justified by the effective PSF analysis) [39], the resolution can be increased by a factor $\sqrt 2$ if two emitters are imaged, two-photon frames are observed, and due to anti-bunching phenomenon one can be sure, that at most one photon can be emitted from a single source within a single frame. Let's challenge this statement, and compute ζ for the described case and check how the fundamental shot noise affects the possibility of obtaining super-resolution. The PDF of measuring a frame consisting of photons at position $x_1, x_2$ is

$$\begin{align} & p_\theta(x_1, x_2) \nonumber\\ & \quad = \frac{ U (x_1+\theta/2) U(x_2 - \theta/2) + U(x_2+\theta/2) U(x_1 - \theta/2) }{2}. \end{align} \tag{ 20 }$$

After substituting the Gaussian form of U, expanding $p_\theta(x_1,x_2)$ into series around θ = 0, and using (3), we obtain the formula for the two-photon frame FI:

$$\begin{align} \mathcal F_{(2)} (\theta) = \frac{\theta ^2}{2~\sigma^4} + \mathcal O(\theta^4/\sigma^6). \end{align} \tag{ 21 }$$

The FI per one photon in this case is equal to $\mathcal F = \frac{1}{2} \mathcal F_{(2)}$, which after using (5) leads to $\zeta = \sqrt[4]{2}$. Despite optimistic assumptions (no noise apart from shot noise, only two-photon frames), our RGL again turns out to be lower than the resolution gain predicted by the PSF analysis.

The RGL introduced in this paper only applies to the transverse resolution. It is however known, that SOFI is a 3D method, and as such, allows to improve axial resolution as well. Nevertheless, the most important aspects of SOFI, which we want to study (i.e. how noise affects obtainable resolution) are all well illustrated with our model, in which both sources lie precisely in the image plane. To extend the whole reasoning to the 3D case, one should proceed analogously as in [40]. The information about the axial separation between the sources is encoded in the size of resulting PSFs, which are affected by the deviation from the image plane. When the deviation is very small, the FI associated with its estimation vanishes because the size of the PSF changes slowly with the axial displacement in this region. As with transverse resolution, the fluctuations will not result in a non-vanishing FI, but will rather result in a larger values of FI for small, yet non-zero deviations.

This case is particularly easy because subsequent photons are not correlated, and the FI per one photon can be calculated directly. Sources are equally bright, and the spatial resolution of the detector is infinite. Each photon position x is independently drawn from the PDF:

$$\begin{align} p_\theta (x) = \frac{U \left(x+\theta/2~\right) + U \left( x-\theta/2~\right)}{2}. \end{align} \tag{ A.2 }$$

Now the FI can be computed with the help of (3) in which vector $\boldsymbol{N}$ consists of just one element—a detected photon position x. The $\mathcal F$ can be therefore expressed as an integral

$$\begin{align} \mathcal F (\theta) = \int_{- \infty}^{\infty} \frac{1}{p_\theta(x)} \left( \frac{\partial p_\theta(x)}{\partial \theta} \right)^2 \, \textrm{d} x, \end{align} \tag{ A.3 }$$

which after substituting the Gaussian form of U(x) simplifies to

$$\begin{align} \sigma^2~\mathcal F (\theta) = \frac{1}{4} - \int^{\infty}_{-\infty} \frac{x^2~\exp\left(-\frac{1}{8~\sigma^2} (\theta -2 x)^2\right)}{2~\sigma^3~\sqrt{2~\pi } \left( \exp \left(\theta x/\sigma^2~\right)+1\right)}\, \textrm{d} x. \end{align} \tag{ A.4 }$$

To obtain the analytical form of the above integral for $\theta \ll \sigma$, one can expand the integrated function in the series around θ = 0, and perform the integration term by term to conclude that

$$\begin{align} \mathcal F (\theta) = \frac{1}{\sigma^2} \left( \frac{\theta^2}{8~\sigma^2}-\frac{\theta^4}{16~\sigma^4}+\frac{\theta^6}{24~\sigma^6}+\ldots \right), \end{align} \tag{ A.5 }$$

which is consistent with (4). In order to obtain the values of $\mathcal F(\theta)$ for larger θ, the introduced integral must be calculated numerically. To compute $\mathcal F$ in case of non-zero pixel size Δx, one needs to construct vector $\boldsymbol{N}$ which consists of mean values of the signal in different pixels only, and then proceed as in appendix A.3.

For the rest of this section, the assumption σ = 1 will be made. Let's consider the simplified model of fluctuations which is slightly more general than the one described in the main text. In each independent frame relative brightness of emitters placed at −θ/2 and θ/2, denoted by q₁ and q₂ respectively, is independently drawn from the same arbitrary probability distribution P(q_i ). The frame time and the mean emitters power are denoted by τ and $\bar P$ respectively—for the sake of simplicity we are going to use the quantity $\bar{n} = \bar P \tau$, which is the only relevant quantity as long as frames are independent, and is proportional to the mean number of photons detected per frame.

As mentioned in the main text, in order to compute the FI per one photon $\mathcal F$, we compute the FI for each fixed number of photons in a frame separately ($\mathcal F_{(n)}$), and then use the formula

$$\begin{align} \mathcal F = \frac{\left\langle \mathcal F_{(n)} \right\rangle _n}{\left\langle n\right\rangle _n}. \end{align} \tag{ A.6 }$$

$\mathcal F_{(n)}$ is the FI associated with the conditional PDF

$$\begin{align} p_\theta(x_1,\ldots,x_n|n) = \int \textrm{d} q_1~\int & \textrm{d} q_2\, p_\theta(x_1,\ldots,x_n|n,q_1,q_2) \nonumber\\ & \times P(q_1,q_2|n). \end{align} \tag{ A.7 }$$

Notice, that the above equation is a generalised form of (12). Equations (10), (11) and (13) are still valid in the general case, and can be used to compute conditional probabilities $p_\theta(x_1,\ldots,x_n|q_1,q_2,n)$ and $p_\theta(x_1,\ldots,x_n|n)$. We are now going to find an explicit expression for $p_\theta(x_1,\ldots,x_n|n)$ to compute $\mathcal F_{(n)}$ directly from the definition of FI. To do so, let us begin with inserting (11) and (A.1) into (10). Before performing the product in (10), we expand each factor into series around θ = 0. After keeping only the leading terms, we obtain

$$\begin{align} & p_{\theta}(x_1, \ldots, x_n | n, q_1, q_2)\nonumber\\ \quad =\ & (2~\pi)^{-n/2} \left( \prod_{i = 1}^n e^{-\frac{1}{8} x_i^2} \right) \nonumber\\ & \times \left( 1+ A_1~\theta+ A_2~\theta^2 + A_3~\theta^3+ A_4~\theta^4+\ldots \right), \end{align} \tag{ A.8 }$$

where

$$\begin{align} A_2 = \frac{1}{8} \sum_{i = 1}^n x_i^2 + \frac{Q_2}{4} \sum_{i\lt j} x_i x_j - \frac{1}{8}n , \end{align} \tag{ A.9 }$$

$$\begin{align} A_4 & = \frac{1}{384} \sum_{i = 1}^n x_i^4 - \frac{n}{64} \sum_{i = 1}^n x_i^2 + \frac{1}{128}n^2 + \frac{Q_2}{96} \sum_{i \neq j} x_i(x_j^3 - 3x_j) \nonumber\\ & \quad + \frac{1}{64} \sum_{i\lt j} x_i^2 x_j^2~ + \frac{Q_2}{32} \sum_{i\lt j,k \neq i, k \neq j} x_i x_j (x_k^2 - 1) \nonumber\\ & \quad + \frac{Q_4}{16} \sum_{i\lt j\lt k\lt m} x_i x_j x_k x_m , \end{align} \tag{ A.10 }$$

and the quantity Q_k is defined as

$$\begin{align} Q_k \equiv \left( \frac{q_1-q_2}{q_1+q_2} \right)^k. \end{align} \tag{ A.11 }$$

We do not specify the form of A₁ and A₃, which are not relevant as will be argued below. Now we can use (A.7) and (A.8) to obtain $p_\theta(x_1,\ldots,x_n|n)$. Let us denote the expected value of a function $X(q_1, q_2)$ with respect to $P(q_1, q_2 |n)$ by $\left\langle X \right\rangle _{q|n}$:

$$\begin{align} \left\langle X \right\rangle _{q|n} \equiv \int X(q_1, q_2) P(q_1, q_2 |n) dq_1 dq_2. \end{align} \tag{ A.12 }$$

Notice, that if we replace all Q_k terms in (A.9) and (A.10) by their mean values $\left\langle Q_k \right\rangle _{q|n}$, we obtain the mean values of coefficients—$\left\langle A_2~\right\rangle _{q|n}$ and $\left\langle A_4~\right\rangle _{q|n}$. Furthermore:

$$\begin{align} p_\theta(x_1,\ldots,x_n|n) & = (2~\pi)^{-n/2} \left( \prod_{i = 1}^n e^{-\frac{1}{8} x_i^2} \right) \nonumber\\ & \quad \times \left( 1+ \left\langle A_2~\right\rangle _{q|n} \theta^2 + \left\langle A_4~\right\rangle _{q|n} \theta^4 + \mathcal O (\theta^6) \right). \end{align} \tag{ A.13 }$$

We have just used the fact that odd coefficients $A_1,A_3,A_5,\ldots$ contain only terms proportional to Q_l , where l is odd. Moreover, from statistical identity of both sources it follows that $\left\langle Q_l \right\rangle _{q|n} = 0$ for odd l—that's the reason why odd coefficient could have been neglected from the beginning, and only even θ powers are present in (A.13). We are now going to calculate $\mathcal F_{(n)}$ using (3) which takes the form

$$\begin{align} \mathcal F_{(n)} = \int \frac{1}{p_\theta(x_1,\ldots,x_n|n)} \left( \frac{\partial p_\theta(x_1,\ldots,x_n|n)}{\partial \theta} \right)^2~\textrm{d} x_1~\ldots \textrm{d} x_n, \end{align} \tag{ A.14 }$$

which simplifies to

$$\begin{align} \mathcal F_{(n)} = (2~\pi)^{-n/2} \int& \prod_{i = 1}^n e^{-\frac{1}{8} x_i^2} \nonumber\\ & \times \left( 4~\left\langle A_2~\right\rangle _{q|n}^2~\theta^2 + \left(16~\left\langle A_2~\right\rangle _{q|n} \left\langle A_4~\right\rangle _{q|n} \right.\right. \nonumber\\ & \left.\left.- 4~\left\langle A_2~\right\rangle _{q|n}^3~\right) \theta^4 + \mathcal O (\theta^6) \right) \textrm{d}x_1\ldots \textrm{d}x_n. \end{align} \tag{ A.15 }$$

After inserting the formulas for $\left\langle A_2~\right\rangle _{q|n}$ and $\left\langle A_4~\right\rangle _{q|n}$, and performing the integration, we obtain

$$\begin{align} \mathcal F_{(n)} & = \theta^2~\left[ \frac{n}{8} + \frac{\left\langle Q_2~\right\rangle _{q|n}^2}{8} n(n-1) \right] \nonumber\\ & \quad + \theta^4~\left[ -\frac{n}{16} - \frac{\left\langle Q_2~\right\rangle _{q|n}^2}{16} n (n-1) \right. \nonumber\\ & \quad \left. \qquad\qquad -\frac{\left\langle Q_2~\right\rangle _{q|n}^3}{16} n(n-1)(n-2) \right] + \mathcal O (\theta^6). \end{align} \tag{ A.16 }$$

Notice, that now $\theta \ll 1$ condition is not sufficient to ensure that the term with θ⁴ is negligible compared to the θ² term, because the powers of n are different in both terms. The approximation

$$\begin{align} \mathcal F_{(n)} \simeq \theta^2~\left[ \frac{n}{8} + \frac{\left\langle Q_2~\right\rangle _{q|n}^2}{8} n(n-1) \right] \end{align} \tag{ A.17 }$$

is nevertheless justified provided $\theta^4 n^3~\ll \theta^2 n^2$, which is equivalent to the condition $\theta \ll n^{-1/2}$. The above inequality holds for all n that give relevant contribution to the final result if $\theta \ll \bar{n}^{-1/2}$. Equation (A.6) allows us to compute the one-photon FI:

$$\begin{align} \mathcal F = \frac{\theta^2}{8} \left[ 1 + \frac{\left\langle \left\langle Q_2~\right\rangle _{q|n}^2 n(n-1) \right\rangle _n}{ \left\langle n \right\rangle _n} \right] + \mathcal O (\theta ^4~\bar{n} ^2). \end{align} \tag{ A.18 }$$

Note that if emitters do not fluctuate, q₁ is always equal to q₂, so $\left\langle Q_2~\right\rangle _{q|n} = 0$, and we recover (A.5) up to the 2nd order of θ—using (A.16) one can additionally check that the coefficient at θ⁴ for non-fluctuating case is also correctly retrieved. The quantity $\left\langle Q_2~\right\rangle _{q|n}$ can be regarded as a measure of fluctuations intensity—it becomes larger, if the normalized difference between q₁ and q₂ takes large values with high probability. As intuitively expected, the FI per one photon increases with $\left\langle Q_2~\right\rangle _{q|n}$, as well as RGL defined in (5), which in our case has a form

$$\begin{align} \zeta_\textrm{max} = \left(1 + \frac{\left\langle \left\langle Q_2~\right\rangle _{q|n}^2 n(n-1) \right\rangle _n}{ \left\langle n \right\rangle _n}\right)^{1/4}. \end{align} \tag{ A.19 }$$

In order to obtain a more specific expression for $\zeta_\textrm{max}$, let us consider the two-level model of emitters mentioned in the main text, i.e.

$$\begin{align} P \left(q_i = q_\textrm{off} \right) &= p \end{align} \tag{ A.20 }$$

$$\begin{align} P\left(q_i = q_\textrm{on} \right)& = 1-p \end{align} \tag{ A.21 }$$

for i ∈ {1, 2}. Recall that $q_\textrm{off}+q_\textrm{on} = 1$, and the fluctuation strength is defined as $\alpha = 1-q_\textrm{off}/q_\textrm{on}$. It's easy to show that the mean number of photons detected per frame is

$$\begin{align} \left\langle n \right\rangle _n = 2~\bar{n} \left( p^2 q_\textrm{off} + p(1-p) + (1-p)^2 q_\textrm{on} \right). \end{align} \tag{ A.22 }$$

The variable Q₂ takes a non-zero value only in two equally probable cases when $q_1~\ne q_2$, so

$$\begin{align} \left\langle Q_2~\right\rangle _{q|n} = 2~\left( q_\textrm{on} - q_\textrm{off} \right)^2 P(q_1 = q_\textrm{on}, q_2 = q_\textrm{off}|n). \end{align} \tag{ A.23 }$$

The probability $P(q_1 = q_\textrm{on},q_2 = q_\textrm{off}|n)$ is calculated using (13), and the fact that

$$\begin{align} P(n|q_1, q_2) = \frac{(\bar{n} (q_1+q_2))^n \exp(-\bar{n}(q_1+q_2))}{n!}. \end{align} \tag{ A.24 }$$

The expression present in (A.19) can be written as an infinite sum

$$\begin{align} \left\langle \left\langle Q_2~\right\rangle _{q|n}^2 n(n-1) \right\rangle _n = \sum_{n = 0}^{\infty} \left\langle Q_2~\right\rangle _{q|n}^2 n(n-1) P(n), \end{align} \tag{ A.25 }$$

which after inserting (A.23) and making some simplifications (e.g. changing variables $q_\textrm{on}, q_\textrm{off}$ to α) takes the form

$$\begin{align} \left\langle \left\langle Q_2~\right\rangle _{q|n}^2 n(n-1) \right\rangle _n = 4 p^2 (1-p)^2~\left(\frac{\alpha}{2-\alpha} \right)^4~\bar{n}^2~S, \end{align} \tag{ A.26 }$$

where

$$\begin{align} S = \sum_{n = 0}^{\infty} \frac{e^{-\bar{n}} \bar{n} ^n}{n!} \left( B e^{A \bar{n}}(1-A)^n+C+D e^{-A \bar{n}} (1+A)^n\right)^{-1}, \end{align} \tag{ A.27 }$$

and the following definitions are used: $A = \frac{\alpha}{2-\alpha}$, $B = p^2 (1-A)^2$, C = 2p(1 − p), $D = (1-p)^2(1+A)^2$. Joining together (A.19), (A.22), and (A.26), we obtain the following expression:

$$\begin{align} \zeta_\textrm{max} = \left(1 + G(p,\alpha, \bar{n}) \bar{n} \right)^{1/4}, \end{align} \tag{ A.28 }$$

where

$$\begin{align} G(p,\alpha, \bar{n}) = \frac{2 p^2 (1-p)^2~\alpha^4}{(2-\alpha)^3 (1-p \alpha)} S. \end{align} \tag{ A.29 }$$

To obtain the value of $\zeta_\textrm{max}$ for arbitrary parameters, one needs to approximate the infinite sum S numerically. However, it is possible to prove that (see appendix A.4)

$$\begin{align} \lim_{\bar{n} \rightarrow \infty} S = \left\{\begin{array}{ll} C^{-1}& 0\lt \alpha \leqslant 1 \\ (B+C+D)^{-1} & \alpha = 0 \end{array} \right. , \end{align} \tag{ A.30 }$$

which allows us to provide an analytical expression for $\zeta_\textrm{max}$ scaling in the infinitely bright sources regime:

$$\begin{align} \lim_{\bar{n} \rightarrow \infty} \zeta_\textrm{max} \bar{n} ^{-1/4} = \left(\frac{p(1-p) \alpha^4}{(2-\alpha)^3(1-p \alpha)} \right)^{1/4}. \end{align} \tag{ A.31 }$$

Numerical analysis show, that the replacement of S by C⁻¹ in (A.28), which leads to equation

$$\begin{align} \zeta_\textrm{max} \simeq \left( 1 + \frac{p(1-p) \alpha^4}{(2-\alpha)^3(1-p \alpha)} \bar{n} \right)^{1/4} , \end{align} \tag{ A.32 }$$

becomes a good approximation for large enough $\bar{n}$. The comparison between $\zeta_\textrm{max}$ computed numerically for finite $\bar{n}$ using (A.28) and its analytical approximation valid for large $\bar{n}$ (A.32) is shown in figure A1.

From now on the spatial resolution of the camera is not assumed to be infinite, and the whole detection area is divided into $M_\textrm{pix}$ pixels of size Δx. The positions of the centroids of subsequent pixels are denoted by $x_1,\ldots,x_{M_\textrm{pix}}$. The detection time is divided into $M_\textrm{fr}$ frames, mth frame covers the time interval $\left[ (m-1) \tau, m \tau \right]$, where $m \in \{1,2,\ldots,M_\textrm{fr} \}$. We are going to consider only two-level blinking model, and stay with its simplified version with independent frames for a while. Our goal is to compute the FI associated with different choices of vector $\boldsymbol{N}$. Let's consider the case studied in the main text, in which $\boldsymbol{N}$ is described by (6). Subsequent frames are independent, so the central limit theorem can be directly used to prove that $\boldsymbol{N}$ is normally distributed.

One needs to use slightly more subtle arguments to extend the above reasoning to the estimation based on 2nd auto-cumulant only (AC2). $\boldsymbol{N}$ consists of 2nd auto-cumulant (variance) estimators for each pixel

$$\begin{align} \boldsymbol{v}_m = \left[n_{1,m}^2 - \left\langle n_1~\right\rangle ^2, \ldots, n_{M_\textrm{pix},m}^2 - \left\langle n_{M_\textrm{pix}} \right\rangle ^2~\right]^T, \end{align} \tag{ A.33 }$$

where $\left\langle n_i \right\rangle = \frac{1}{M_\textrm{fr}}\sum_{m = 1}^{M_\textrm{fr}} n_{i,m}$ denotes the mean value estimator. Unfortunately, this estimator depends on detected photon numbers from different frames, so vectors $\boldsymbol{v}_m$ are not mutually independent anymore. However, in the limit $M_\textrm{fr} \rightarrow \infty$ the mean value estimator becomes very accurate compared to the variability of the number of photons in a given pixel in a single frame because the variance of n_i,m does not depend on $M_\textrm{fr}$, and the variance of $\left\langle n_i \right\rangle$ scales as $1/M_\textrm{fr}$. That means, that the replacement of the mean value estimator with its exact value in $\boldsymbol{v}_m$ does not affect the distribution of $\boldsymbol{N}$ in the limit $M_\textrm{fr} \rightarrow \infty$. Vectors $\boldsymbol{v}_m$ become independent after making the described replacement, which allows us to conclude, that $\boldsymbol{N}$ is normally distributed.

In order to compute the FI associated with the normally distributed vector $\boldsymbol{N}$ in the most general case in which both the mean vector $\boldsymbol{\mu}$, and the covariance matrix $\boldsymbol{\Sigma}$ depend on the estimated parameter θ, one can use the formula [35]

$$\begin{align} \mathcal F_\textrm{(meas)} = \frac{\partial \boldsymbol{\mu}^\top}{\partial \theta} \boldsymbol{\Sigma}^{-1} \frac{\partial \boldsymbol{\mu}}{\partial \theta} + \frac{1}{2} \textrm{tr} \left( \boldsymbol{\Sigma}^{-1} \frac{\partial \boldsymbol{\Sigma}}{\partial \theta} \boldsymbol{\Sigma}^{-1} \frac{\partial \boldsymbol{\Sigma}}{\partial \theta} \right). \end{align} \tag{ A.34 }$$

If vectors $\boldsymbol{v}_m$ are independent and identically distributed (i.i.d.), the mean vector and the covariance matrix for each $\boldsymbol{v}_m$ are denoted by $\boldsymbol{\mu}_{(1)}$ and $\boldsymbol{\Sigma}_{(1)}$, then $\boldsymbol{\mu} = \boldsymbol{\mu}_{(1)}$, and $\boldsymbol{\Sigma} = \frac{1}{M_\textrm{fr}} \boldsymbol{\Sigma}_{(1)}$. We therefore see, that the 2nd term in (A.34) is negligible compared to the 1st term in the limit $M_\textrm{fr} \rightarrow \infty$ (provided the first term is non-zero), and hence in this limit we may write

$$\begin{align} \mathcal F_\textrm{(meas)} = \frac{\partial \boldsymbol{\mu}^\top}{\partial \theta} \boldsymbol{\Sigma}^{-1} \frac{\partial \boldsymbol{\mu}}{\partial \theta} = M_\textrm{fr} \frac{\partial \boldsymbol{\mu}_{(1)}^\top}{\partial \theta} \boldsymbol{\Sigma}_{(1)}^{-1} \frac{\partial \boldsymbol{\mu}_{(1)}}{\partial \theta}. \end{align} \tag{ A.35 }$$

In order to compute the FI per one photon $\mathcal F$ one needs to compute the elements of $\boldsymbol{\mu}$ and $\boldsymbol{\Sigma}$, use (A.35), and then divide $\mathcal F _\textrm{(meas)}$ by the average total photon number. Let $v_{1,m}, v_{2,m},\ldots,v_{n,m}$ be the elements of the vector $\boldsymbol{v}_m$. In some cases we are going to use a short-hand notation $v_{i,1} \equiv v_i, n_{i,1} \equiv n_i$, because the 2nd index can be omitted in many situations when single frame statistics are considered. Then:

$$\begin{align} \boldsymbol{\mu} & = \left[ \left\langle v_1~\right\rangle , \left\langle v_2~\right\rangle , \ldots,\left\langle v_n \right\rangle \right]^T, \end{align} \tag{ A.36 }$$

$$\begin{align} \boldsymbol{\Sigma} & = \frac{1}{M_\textrm{fr}} \left( \left[ \begin{array}{cccc} \left\langle v_1~v_1~\right\rangle & \left\langle v_1~v_2~\right\rangle & \ldots & \left\langle v_1~v_n \right\rangle \\ \left\langle v_2~v_1~\right\rangle & \left\langle v_2~v_2~\right\rangle & \ldots & \left\langle v_2~v_n \right\rangle \\ \vdots & \vdots & \ddots & \vdots \\ \left\langle v_n v_1~\right\rangle & \left\langle v_n v_2~\right\rangle & \ldots & \left\langle v_n v_n \right\rangle \end{array} \right] - \boldsymbol{\mu} \boldsymbol{\mu}^\top \right). \end{align} \tag{ A.37 }$$

In every considered case each v_i (i ∈ {1, 2, ..., n}) can be written as a linear combination of elements of the form $n_j^{k_1} n_l^{k_2}$ where $j,l \in \{1,2,\ldots,M_\textrm{pix} \}$, and $k_1,k_2$ are natural exponents (possibly zero). Therefore, it is enough to be able to compute expected values of products $\left\langle n_{j_1}^{k_1} n_{j_2}^{k_2} \ldots n_{j_r}^{k_r} \right\rangle$ for $r \leqslant 4$ to reconstruct all terms of $\boldsymbol{\mu}$ and $\boldsymbol{\Sigma}$. The procedure used to compute these expected values is as follows. PSFs of both sources are numerically integrated over different pixels—we construct variables

$$\begin{align} U_{j,1} & = \int_{x_j-\Delta x/2}^{x_j+\Delta x/2} U(x + \theta/2) dx ,~ \nonumber\\ U_{j,2} & = \int_{x_j-\Delta x/2}^{x_j+\Delta x/2} U(x - \theta/2) dx, \end{align} \tag{ A.38 }$$

where $j\in \{1,2,\ldots,M_\textrm{pix}\}$ denotes the pixel label. Now we use the fact, that the number of photons detected in each pixel, when the sources brightness are fixed, is described by a Poisson distribution with a mean value

$$\begin{align} \left\langle n_j | q_1, q_2~\right\rangle = \left( q_1~U_{j,1} + q_2~U_{j,2} \right) \bar{n}, \end{align} \tag{ A.39 }$$

where $q_1, q_2~\in \{q_\textrm{off},q_\textrm{on} \}$ denote relative brightness of the 1st and the 2nd emitter respectively. When external, Poissonian noise (studied in section 5.2) is present, the above equation should be modified to

$$\begin{align} \left\langle n_j | q_1, q_2~\right\rangle = \left( q_1~U_{j,1} + q_2~U_{j,2} \right) \bar{n} + \mu_B, \end{align} \tag{ A.40 }$$

where µ_B is the mean value of this noise in each pixel in each frame. Further steps remain valid, as conditional random variables $n_j|q_1,q_2$ are Poissonian and mutually independent with and without external noise. The expected value of a product of their powers is

$$\begin{align} \left\langle n_{j_1}^{k_1} n_{j_2}^{k_2} \ldots n_{j_r}^{k_r}|q_1,q_2~\right\rangle = \prod_{i = 1}^r M_{k_i}\left(\left\langle n_{j_i} | q_1, q_2~\right\rangle \right), \end{align} \tag{ A.41 }$$

where M_k (v) denotes a kth raw moment of Poisson distribution with a mean value v, e.g. M₀(ν) = 1, M₁(ν) = ν, $M_2(\nu) = \nu^2 +\nu$, $M_3(\nu) = \nu^3+3~\nu^2+\nu$, $M_4(v) = \nu^4+6~\nu^3 + 7~\nu^2 + \nu$. Equation (A.41) is only valid if indices $j_1,\ldots,j_r$ are mutually different—if any index repeats, one should replace an expression of the form $n_j^{k_1} n_j^{k_2}$ with an expression $n_j^{k_1+k_2}$, repeat such a procedure as long as there are any repetitions left, and only then use (A.41) directly. Already described steps allow us to compute conditional expected values. In order to compute the desired expected values $\left\langle n_{j_1}^{k_1} n_{j_2}^{k_2} \ldots n_{j_r}^{k_r} \right\rangle$ one only needs to average the conditional ones over four different configurations of the emitters using the formula

$$\begin{align} \left\langle X\right\rangle & = p^2~\left\langle X|q_\textrm{off}, q_\textrm{off} \right\rangle +p(1-p) \left( \left\langle X|q_\textrm{off}, q_\textrm{on} \right\rangle \right. \nonumber\\ & \quad \left.+\left\langle X|q_\textrm{on}, q_\textrm{off} \right\rangle \right) +(1-p)^2~\left\langle X|q_\textrm{off}, q_\textrm{off} \right\rangle. \end{align} \tag{ A.42 }$$

By performing the described steps numerically, we obtain $\mathcal F(\theta)$ (and consequently ζ) associated with different image reconstruction algorithms in the simplified blinking model.

We will now show, how to extend this scheme to the case of Markov process based realistic model. Let's first remind, that the relative brightness of each emitter is described by a Markov process with two possible states with different relative brightness: $q_\textrm{on}$ and $q_\textrm{off}$. The brightness of each emitter is a function of time P_i (t) (i ∈ {1, 2}), which takes two possible values: $q_\textrm{off} \bar P$, $q_\textrm{on} \bar P$. During a short time interval $\left[t, t+\delta t\right]$ ith emitter emits P_i (t)δt photons on average. Lifetimes of on- and off- states are equal to $\tau_\textrm{off}$ and $\tau_\textrm{on}$ respectively. The probability that a given emitter remains in a fixed state with a lifetime τ_i for a time period t is proportional to exp(−t/τ_i ). Let's now introduce a more formal description of the Markov process, which leads to such an exponential behaviour. At any time t, the state of an emitter is described by a vector $\begin{bmatrix} p_\textrm{off} \\ p_\textrm{on} \end{bmatrix}$, where $p_\textrm{off}$ and $p_\textrm{on}$ denote the probabilities of finding the emitter in off- and on- state respectively. The time evolution of the emitter state is given by

$$\begin{align} \begin{bmatrix} p_\textrm{off} \\ p_\textrm{on} \end{bmatrix}(t + \Delta t) = \boldsymbol{T} (\Delta t) \begin{bmatrix} p_\textrm{off} \\ p_\textrm{on} \end{bmatrix}(t), \end{align} \tag{ A.43 }$$

where $\boldsymbol{T} (\Delta t)$ is a transition matrix defined as

$$\begin{align} \boldsymbol{T}(\Delta t) = \begin{bmatrix}t_{00} & t_{10} \\ t_{01} & t_{11} \end{bmatrix}(\Delta t) = \exp \left( \Delta t \begin{bmatrix} -\tau_\textrm{off}^{-1} & \tau_\textrm{on}^{-1} \\ \tau_\textrm{off}^{-1} & -\tau_\textrm{on}^{-1} \end{bmatrix} \right). \end{align} \tag{ A.44 }$$

It's easy to check, that after a long evolution the state always converges to

$$\begin{align} \begin{bmatrix} p_\textrm{off} \\ p_\textrm{on} \end{bmatrix}(\Delta t \rightarrow \infty) = \begin{bmatrix} \tilde p_\textrm{off} \\ \tilde p_\textrm{on} \end{bmatrix} = \frac{1}{\tau_\textrm{off} + \tau_\textrm{on}}\begin{bmatrix} \tau_\textrm{off} \\ \tau_\textrm{on} \end{bmatrix}. \end{align} \tag{ A.45 }$$

$\begin{bmatrix} \tilde p_\textrm{off} \\ \tilde p_\textrm{on} \end{bmatrix}$ is a stationary state of the process, and will be used as an initial state in our considerations—if no information about the previous run of the process is available, the probability of finding an emitter in a given state is proportional to its lifetime. Using the transition matrix $\boldsymbol{T}(\Delta t)$ elements, let us define the $\boldsymbol{S}(\Delta t)$ matrix:

$$\begin{align} \boldsymbol{S}(\Delta t) = \begin{bmatrix}t_{00}(\Delta t) q_\textrm{off} \bar P & t_{01}(\Delta t) q_\textrm{off} \bar P \\ t_{10}(\Delta t) q_\textrm{on} \bar P & t_{11}(\Delta t) q_\textrm{on} \bar P \end{bmatrix}, \end{align} \tag{ A.46 }$$

which allows us to write down formulas for temporal brightness correlations in a compact way:

$$\begin{align} & \left\langle P_i(t_1) P_i(t_2) \ldots P_i(t_r) \right\rangle \nonumber \\ & \quad = \begin{bmatrix} \tilde p_\textrm{off} & \tilde p_\textrm{on} \end{bmatrix} \boldsymbol{S}(t_2-t_1) \boldsymbol{S}(t_3-t_2)\ldots \boldsymbol{S}(t_r-t_{r-1}) \begin{bmatrix} q_\textrm{off} \bar P \\ q_\textrm{on} \bar P \end{bmatrix}. \end{align} \tag{ A.47 }$$

In the above formula $t_1 \leqslant t_2 \leqslant \cdots \leqslant t_r$, and $\left\langle \bullet \right\rangle$ denotes averaging over Markov processes P₁(t), P₂(t). Two emitters are independent, so in order to compute a product in which P₁ and P₂ terms are mixed, one can use the formula

$$\begin{align} \left\langle \mathcal G_1~\left [P_1(t) \right] \mathcal G_2~\left [P_2(t) \right] \right\rangle = \left\langle \mathcal G_1~\left [P_1(t) \right] \right\rangle \left\langle \mathcal G_2~\left [P_2(t) \right] \right\rangle , \end{align} \tag{ A.48 }$$

which is true for all functionals $\mathcal G_1$, $\mathcal G_2$. Further on, the following integrals of correlations over detection time frames will be useful:

$$\begin{align} \chi_1 & = \int_0^\tau \left\langle P_i(t)\right\rangle \textrm{d}t = \left\langle P_i\right\rangle \tau \end{align} \tag{ A.49 }$$

$$\begin{align} \chi_{2,m} & = \int_0^\tau \textrm{d}t_1~\int_{(m-1) \tau}^{m \tau} \textrm{d}t_2~\left\langle P_i(t_1) P_i(t_2) \right\rangle \end{align} \tag{ A.50 }$$

$$\begin{align} \chi_{3,m} & = \int_0^\tau \textrm{d}t_1~\int_0^\tau \textrm{d}t_2~\int_{(m-1) \tau}^{m \tau} \textrm{d}t_3~\left\langle P_i(t_1) P_i(t_2) P_i(t_3) \right\rangle \end{align} \tag{ A.51 }$$

$$\begin{align} \chi_{4,m} & = \int_0^\tau \textrm{d}t_1~\int_0^\tau \textrm{d}t_2~\int_{(m-1) \tau}^{m \tau} \textrm{d}t_3~\int_{(m-1) \tau}^{m \tau} \textrm{d}t_4~ \nonumber\\ & \quad \times \left\langle P_i(t_1) P_i(t_2) P_i(t_3) P_i(t_4) \right\rangle. \end{align} \tag{ A.52 }$$

At this point, we are prepared to attack the problem of computing $\mathcal F$. Although the frames are now correlated, we still consider vectors $\boldsymbol{N}$ that can be written as in (6), and we can use the central limit theorem in its extended version [38] because correlations between frames decay exponentially with time. Therefore, it is again enough to calculate $\boldsymbol{\mu}$ and $\boldsymbol{\Sigma}$ associated with $\boldsymbol{N}$, and then use (A.34). The scaling of both terms in this equation remains the same, so the 2nd term again disappears in the limit $M_\textrm{fr} \rightarrow \infty$. Equation (A.36) is still valid, and can be used to calculate $\boldsymbol{\mu}$, but in order to compute $\boldsymbol{\Sigma}$ elements one needs to take into account correlations between frames:

$$\begin{align} \boldsymbol{\Sigma}_{ij} = \frac{1}{M_\textrm{fr}^2} \sum_{m,m^{^{\prime}} = 1}^{M_\textrm{fr}} \textrm{cov}(v_{i,m}, v_{j,m^{^{\prime}}}). \end{align} \tag{ A.53 }$$

In the limit $M_\textrm{fr} \rightarrow \infty$, using the homogeneity of the Markov processes, we can simplify our formula:

$$\begin{align} \boldsymbol{\Sigma}_{ij} = \frac{1}{M_\textrm{fr}} \left( \textrm{cov}(v_{i,1},v_{j,1})+2~\sum_{m = 2}^{\infty} \textrm{cov}(v_{i,1}, v_{j,m}) \right). \end{align} \tag{ A.54 }$$

Analogously to the previous case, photon numbers in different pixels and time frames are uncorrelated and described by a Poisson distribution if functions $P_1(t), P_2(t)$ are fixed, and we have:

$$\begin{align} \left\langle n_{j,m} |P_1(t), P_2(t) \right\rangle \equiv \mu_{j,m} = \int_{(m-1)\tau}^{m \tau} \left( P_1(t) U_{1,j} + P_2(t) U_{2,j}\right) dt. \end{align} \tag{ A.55 }$$

Conditional products of variables n_j,m are computed with the help of the formula

$$\begin{align} \left\langle n_{j_1, m_1}^{k_1} n_{j_2, m_2}^{k_2} \ldots n_{j_r, m_r}^{k_r}|P_1(t), P_2(t) \right\rangle = \prod_{i = 1}^r M_{k_i}\left(\mu_{j_i,m_i} \right) \end{align} \tag{ A.56 }$$

valid if pairs $(j_i, m_i)$ mutually differ on at least one position. All terms of $\boldsymbol{\mu}$ and $\boldsymbol{\Sigma}$ are linear combinations of expectation values (where averaging over Markov processes is made) $\left\langle n_{j_1, m_1}^{k_1} n_{j_2, m_2}^{k_2} \ldots n_{j_r, m_r}^{k_r} \right\rangle$. We want restrict ourselves to $\boldsymbol{N}$ which consist of 2nd order correlations at most, so from now on we assume that $k_1+k_2+ \cdots +k_r \leqslant 4$. Then, after expanding the RHS of (A.56), and using formulas for Poisson distribution moments, we see that all conditional expected values are linear combinations of products of µ_j,m with at most four terms. Expectation values required to reconstruct $\boldsymbol{\mu}$ and Σ are linear combinations of similar products averaged over Markov processes. Such products can be written with the help of variables χ, defined in (A.49)–(A.52),

$$\begin{align} \left\langle \mu_{j,m} \right\rangle = (U_{1,j} + U_{2,j}) \chi_1, \end{align} \tag{ A.57 }$$

$$\begin{align} \left\langle \mu_{j,0} \mu_{j^{^{\prime}},m} \right\rangle & = (U_{1,j}U_{1,j^{^{\prime}}} + U_{2,j}U_{2,j^{^{\prime}}}) \chi_{2,m} + (U_{1,j} U_{2,j^{^{\prime}}} \nonumber\\ & \quad + U_{2,j} U_{1,j^{^{\prime}}}) \chi_1^2, \end{align} \tag{ A.58 }$$

$$\begin{align} \left\langle \mu_{j,0} \mu_{j^{^{\prime}},0} \mu_{j^{^{\prime\prime}},m} \right\rangle & = U_{1,j}U_{1,j^{^{\prime}}}U_{1,j^{^{\prime\prime}}} \chi_{3,m} +U_{1,j}U_{1,j^{^{\prime}}}U_{2,j^{^{\prime\prime}}} \chi_{2,0}\chi_1~\nonumber\\ & \quad +(U_{1,j}U_{2,j^{^{\prime}}}U_{1,j^{^{\prime\prime}}}+U_{1,j}U_{2,j^{^{\prime}}}U_{2,j^{^{\prime\prime}}}) \chi_{2,m} \chi_1 \nonumber\\ & \quad + (1~\leftrightarrow 2), \end{align} \tag{ A.59 }$$

$$\begin{align} & \left\langle \mu_{j,0} \mu_{j^{^{\prime}},0} \mu_{j^{^{\prime\prime}},m} \mu_{j^{^{\prime\prime\prime}},m} \right\rangle \nonumber\\ & = U_{1,j}U_{1,j^{^{\prime}}}U_{1,j^{^{\prime\prime}}}U_{1,j^{^{\prime\prime\prime}}} \chi_{4,m} +(U_{1,j}U_{1,j^{^{\prime}}}U_{1,j^{^{\prime\prime}}}U_{2,j^{^{\prime\prime\prime}}}\nonumber\\ & \quad+U_{1,j}U_{1,j^{^{\prime}}}U_{2,j^{^{\prime\prime}}}U_{1,j^{^{\prime\prime\prime}}}+U_{1,j}U_{2,j^{^{\prime}}}U_{1,j^{^{\prime\prime}}}U_{1,j^{^{\prime\prime\prime}}}\nonumber\\ & \quad+U_{1,j}U_{2,j^{^{\prime}}}U_{2,j^{^{\prime\prime}}}U_{2,j^{^{\prime\prime\prime}}})\chi_{3,m} \chi_1~\nonumber \\ & \quad+ U_{1,j}U_{1,j^{^{\prime}}}U_{2,j^{^{\prime\prime}}}U_{2,j^{^{\prime\prime\prime}}} \chi_{2,0}^2 + (U_{1,j}U_{2,j^{^{\prime}}}U_{1,j^{^{\prime\prime}}}U_{2,j^{^{\prime\prime\prime}}}\nonumber\\ & \quad+U_{1,j}U_{2,j^{^{\prime}}}U_{2,j^{^{\prime\prime}}}U_{1,j^{^{\prime\prime\prime}}}) \chi_{2,m}^2 +(1~\leftrightarrow 2). \end{align} \tag{ A.60 }$$

Notation $+(1~\leftrightarrow 2)$ means that terms with swapped indices 1 and 2, that correspond to the 1st and 2nd emitter, should be added.

Let's summarize the procedure used to compute $\boldsymbol{\mu}$ and $\boldsymbol{\Sigma}$ for the Markov process based model. First, (A.36) and (A.54) are applied, and all terms are expressed as linear combinations of expected values $\left\langle n_{j_1, m_1}^{k_1} n_{j_2, m_2}^{k_2} \ldots n_{j_r, m_r}^{k_r} \right\rangle$. Then, conditional expected values $\left\langle n_{j_1, m_1}^{k_1} n_{j_2, m_2}^{k_2} \ldots n_{j_r, m_r}^{k_r}|P_1(t), P_2(t) \right\rangle$ are computed with the help of (A.56). Averaging over Markov processes P₁(t), P₂(t) is made after expanding the RHS of (A.56), formulas ((A.57)–(A.60)) are then utilized to express $\left\langle n_{j_1, m_1}^{k_1} n_{j_2, m_2}^{k_2} \ldots n_{j_r, m_r}^{k_r} \right\rangle$ using variables U_1,j , U_2,j , and χ (defined in (A.49)–(A.52)). Notice, that the sum present in (A.57) is infinite, but in our case all sums converge. Moreover, it is possible to express all the terms of $\boldsymbol{\mu}$ and $\boldsymbol{\Sigma}$ in the considered cases with the help of variables $U_{1,j}, U_{2,j}$ (computed numerically), χ₁, χ_2,1, χ_3,1, χ_4,1, and the following infinite sums (all of them converge):

$$\begin{align} S_1 & = \sum_{m = 2}^\infty\left( \chi_{2,m} - \chi_1^2~\right), \end{align} \tag{ A.61 }$$

$$\begin{align} S_2 & = \sum_{m = 2}^\infty\left( \chi_{3,m} - \chi_{2,1} \chi_1~\right), \end{align} \tag{ A.62 }$$

$$\begin{align} S_3 & = \sum_{m = 2}^\infty\left( \chi_{4,m} - \chi_{2,1}^2~\right), \end{align} \tag{ A.63 }$$

$$\begin{align} S_4 & = \sum_{m = 2}^\infty\left( \chi_{2,m}^2 - \chi_1^4~\right) \end{align} \tag{ A.64 }$$

where χ variables are expressed as functions of the setup parameters analytically with the help of (A.47), (A.49)–(A.52), and then sums $S_1, S_2, S_3, S_4$ are also computed analytically (they can be expressed as sums of geometric series).

In this section we provide a rigorous proof of (A.30). For α = 0, we have A = 0, the value of the sum does not depend on $\bar{n}$ and can be easily computed, and the proof becomes trivial. For $0\lt \alpha \leqslant 1$ we need to prove that

$$\begin{align} \lim_{\bar{n} \rightarrow \infty}\!\sum_{n = 0}^{\infty} \frac{e^{-\bar{n}} \bar{n} ^n}{n!} \!\left( B e^{A \bar{n}}(1\!-A)^n{+}C+D e^{-A \bar{n}} (1{+}A)^n\right)^{-1}\!\!=\!C^{-1}, \end{align} \tag{ A.65 }$$

where $0\lt A\leqslant 1$, $B,D \geqslant 0$, $C\gt 0$. Let us define

$$\begin{align} X_n(\bar{n}) : = \frac{e^{-\bar{n}} \bar{n} ^n}{n!} \left( B e^{A \bar{n}}(1-A)^n+C+D e^{-A \bar{n}} (1+A)^n\right)^{-1}. \end{align} \tag{ A.66 }$$

We are going to use the following property of Poisson distribution:

$$\begin{align} \forall _{\delta \gt 0, \epsilon \gt 0} ~ \exists_{N_0} ~ \forall_{\bar{n} \gt N_0}\!\!: ~ 0 \lt \sum_{n = 0}^\infty \frac{e^{- \bar{n}} \bar{n} ^n}{n!} -\sum_{n = (1-\epsilon) \bar{n}}^{(1+\epsilon) \bar{n}} \frac{e^{- \bar{n}} \bar{n} ^n}{n!} \lt \delta. \end{align} \tag{ A.67 }$$

Intuitively, for large enough $\bar{n}$, all probable values of Poisson distribution are located in the range $\bar{n} (1~\pm \epsilon)$, because the standard deviation of Poisson distribution with mean $\bar{n}$ is $\sqrt{\bar{n}}$. Let's notice that $X_n(\bar{n} ) \leqslant C^{-1} \frac{\bar{n}^n e^{- \bar{n}}}{n!}$, so the following inequality is true:

$$\begin{align} & 0\lt \sum_{n = 0}^\infty X_n(\bar{n}) -\sum_{n = (1-\epsilon) \bar{n}}^{(1+\epsilon) \bar{n}} X_n( \bar{n}) \leqslant C^{-1} \nonumber\\ & \quad \times\left( \sum_{n = 0}^\infty \frac{e^{- \bar{n}} \bar{n} ^n}{n!} -\sum_{n = (1-\epsilon) \bar{n}}^{(1+\epsilon) \bar{n}} \frac{e^{- \bar{n}} \bar{n} ^n}{n!} \right), \end{align} \tag{ A.68 }$$

which after using (A.67) allows us to conclude that:

$$\begin{align} \forall _{ \epsilon \gt 0} ~ \lim_{\bar{n} \rightarrow \infty} \sum_{n = 0}^\infty X_n(\bar{n}) = \lim_{\bar{n} \rightarrow \infty} \sum_{n = (1-\epsilon) \bar{n}}^{(1+\epsilon) \bar{n}} X_n(\bar{n}). \end{align} \tag{ A.69 }$$

Let's now fix $\epsilon , \delta_1 \gt 0$ satisfying

$$\begin{align} A+(1-\epsilon) \log (1-A) \lt - \delta_1, \end{align} \tag{ A.70 }$$

$$\begin{align} -A+(1+\epsilon) \log (1+A) \lt -\delta_1. \end{align} \tag{ A.71 }$$

For A = 1 the first inequality may be omitted. It's easy to show, that the described choice of positive constants and δ₁ is always possible. Such a choice allows us to conclude, that for $n \in \left[ \bar{n} (1-\epsilon), \bar{n} (1+\epsilon) \right]$:

$$\begin{align} e^{A \bar{n}} (1-A)^n \leqslant e^{\bar{n} (A+(1-\epsilon) \log (1-A))} \lt e^{- \delta_1~\bar{n}}, \end{align} \tag{ A.72 }$$

$$\begin{align} e^{-A \bar{n}} (1+A)^n \leqslant e^{\bar{n} (-A+(1+\epsilon) \log (1+A))}\lt e^{- \delta_1~\bar{n}}. \end{align} \tag{ A.73 }$$

Therefore,

$$\begin{align} 1 \leqslant \left( \frac{e^{- \bar{n}} \bar{n} ^n}{C n!} \right) / X_n(\bar{n}) \leqslant 1 + \frac{B+D}{C} e^{- \delta_1~\bar{n}}, \end{align} \tag{ A.74 }$$

and consequently

$$\begin{align} & \left( 1+\frac{B+D}{C} e^{-\delta_1~\bar{n}} \right)^{-1} \sum_{n = (1-\epsilon) \bar{n}}^{(1+\epsilon) \bar{n}} \frac{e^{- \bar{n}} \bar{n} ^n}{C n!} \leqslant \sum_{n = (1-\epsilon) \bar{n}}^{(1+\epsilon) \bar{n}} X_n(\bar{n}) \nonumber \\ & \quad \leqslant \sum_{n = (1-\epsilon) \bar{n}}^{(1+\epsilon) \bar{n}} \frac{e^{- \bar{n}} \bar{n} ^n}{C n!}. \end{align} \tag{ A.75 }$$

Since $e^{-\delta_1~\bar{n}} \rightarrow 0$ for $\bar{n} \rightarrow \infty$, and sum $\sum_{n = (1-\epsilon) \bar{n}}^{(1+\epsilon) \bar{n}} \frac{e^{- \bar{n}} \bar{n} ^n}{C n!}$ converges to C⁻¹ (as a consequence of (A.67)), we can use the so called sandwich theorem to finally conclude that

$$\begin{align} \lim_{\bar{n} \rightarrow \infty} \sum_{n = 0}^\infty X_n(\bar{n}) = \lim_{\bar{n} \rightarrow \infty} \sum_{n = (1-\epsilon) \bar{n}}^{(1+\epsilon) \bar{n}} X_n(\bar{n}) = C^{-1}. \end{align} \tag{ A.76 }$$

Zoom In Zoom Out Reset image size