# Fast, accurate, and predictive method for atom detection in site-resolved images of microtrap arrays

Marc Cheneau, Romaric Journet, Matthieu Boffety, Francois Goudail, Caroline Kulcsár, Pauline Trouvé-Peloux

► To cite this version:

Marc Cheneau, Romaric Journet, Matthieu Boffety, Francois Goudail, Caroline Kulcsár, et al. Fast, accurate, and predictive method for atom detection in site-resolved images of microtrap arrays. 2025. [hal-05239297v1](https://hal HAL Id: [hal-05239297](https://hal.science/hal-05239297v1)
# Fast, accurate, and predictive method for atom detection in site-resolved images of microtrap arrays

Marc Cheneau,* Romaric Journet, Matthieu Boffety, François Goudail, and Caroline Kulcsár

Université Paris-Saclay, Institut d'Optique Graduate School,
CNRS, Laboratoire Charles Fabry, 91127, Palaiseau, France

Pauline Trouvé-Peloux

DTIS, ONERA, Université Paris-Saclay, 91123, Palaiseau, France

(Dated: September 4, 2025)

We introduce a new method, rooted in estimation theory, to detect individual atoms in site-resolved images of microtrap arrays, such as optical lattices or optical tweezers arrays. Using labelled test images, we demonstrate drastic improvement of the detection accuracy compared to the popular method based on Wiener deconvolution when the inter-site distance is comparable to the radius of the point spread function. The runtime of our method scales approximately linearly with the number of sites, and remains well below 100 ms for an array of 100 × 100 sites on a desktop computer. It is therefore fully compatible with a real-time usage. Finally, we propose a rigorous definition for the signal-to-noise ratio of the problem, and show that it can be used as a predictor for the detection error rate. Our work opens the prospect for future experiments with increased array sizes, or reduced inter-site distances.

## I. INTRODUCTION

Over the past 15 years, site-resolved imaging of atoms or molecules in arrays of optical microtraps, whether optical lattices or optical tweezers arrays, has become a central technique for probing spatial correlations in quantum systems [1–4]. In a typical experimental setting, a one- or two-dimensional array loaded with atoms, with an inter-site distance ranging from a few hundred nanometers to a few micrometers, is imaged onto a digital camera with a high-resolution objective. After calibration of the coordinates of the sites in the image, and of the profile of the point spread function (PSF) of the optical system, the image is processed to detect the presence or absence of an atom in each site. The reconstruction of the site occupancies is a classical estimation-detection problem [5, 6], where the signal to be estimated before detection consists in the ‘brightness’ of each site.

When the inter-site distance (a) is much larger than the PSF radius ($r_{PSF}$), a simple binning approach is often sufficient to accurately estimate the site brightnesses [7–10]. Otherwise, one has to resort to more elaborate estimation methods, which we broadly categorize as local non-linear least squares [11–14], local iterative [15, 16], deconvolution [17–24], local maximum likelihood [25, 26], and supervised [27, 28] or unsupervised [29] neural networks. Recently, the authors of ref. [30] benchmarked several methods with a Cramér-Rao bound using labelled test images simulating an optical tweezer array experiment. They showed that a non-linear least-squares approach achieves a better accuracy than, for instance, the popular Wiener or Richardson-Lucy deconvolution, but at the cost of a prohibitive computation time.

In this article, we employ a generalized Wiener filter approach—which provides the optimal linear estimate of the signal, given the statistical properties of the signal and the noise [5]—to the reconstruction of the site occupancies. Using labelled test images, we benchmark our method with the standard Wiener deconvolution, and show that our method drastically improves the detection accuracy when the sites are not optically resolved, that is $a \lesssim r_{PSF}$. This means that some settings which were deemed improper for single-atom detection—because of a short inter-site distance or a low site brightness, for instance—can actually be processed efficiently.

Importantly, our method is applicable in real-time, as its runtime scales approximately linearly with the number of sites, and remains well below 100 ms for an array of 100 × 100 sites on a desktop computer. To help the reader come to grips with the numerical implementation, we provide them with a tutorial example in the form of a Jupyter notebook [31].

Last but not least, we propose a rigorous definition of the signal-to-noise ratio for the reconstruction problem, i.e. a synthetic indicator of the intrinsic ‘difficulty’ of the estimation-detection problem under given experimental conditions. Furthermore, we show that this indicator is a good predictor for the detection accuracy, and can therefore be used to rationalize the design of the imaging system for future experiments. For instance, one could determine the minimum numerical aperture (NA) necessary to achieve a target detection error rate; since a reduced NA usually comes with an increased field of view, larger arrays could be processed in neutral-atom based platforms for quantum computation or quantum simulation [32–34]. Another possibility is to minimize the inter-site spacing for a given NA, offering a direct benefit to experiments on collective light scattering in tweezer arrays [35], or dipolar gases in optical lattices [36], which require short inter-site distances to make the dipolar in-

<footnote>* marc.cheneau@institutoptique.fr</footnote>
teractions sizable.

## II. MODEL AND METHODOLOGY

To quantitatively evaluate the performance of different reconstruction methods, we generate labelled test images using an idealized model [37]. We define a square array with $N_s$ sites in an image with $N_p$ pixels, and denote a the inter-site distance (measured in pixels). The top left site coincides with the center of a pixel unless a global offset is applied to the site positions.

Each site has the same occupancy probability $p$. The brightness of the empty sites is zero, while that of the occupied site follows a normal distribution with mean $\mu$ and variance $\sigma^2$. The finite variance accounts for site-to-site variations of the scattering rate due to non-uniform trap depths or illumination, for instance. The atoms are treated as independent scatterers, which is a good approximation as long as $a$ is much larger than the imaging wavelength divided by $2\pi$.

Once the brightness of each site has been defined, we build a noiseless image by summing the PSF of each site, multiplied by the corresponding brightness, and adding a uniform background $k$. We have chosen a Gaussian PSF for simplicity [38]; it is characterized by its half width at half maximum (HWHM), identical for all sites, and normalized such that the area integral of the PSF is equal to unity. We truncate the PSF when the distance to the center is larger than 3 times the HWHM.

Finally, we model the shot noise by treating each pixel of the noiseless image as the mean value of a Poisson distribution. We also add a normally distributed noise with zero mean and uniform variance $r^2$ to account for the read noise of the camera. This noise model is well suited for CMOS cameras, assuming that the noise induced by the dark current and the digitization are negligible.

The relationship between the site brightnesses and the pixel values in the test image can be summarized in the following matrix equation:

$$y \sim \mathcal{P}(Mx + k) + \mathcal{N}(0, r^2). \quad (1)$$

Here, and in the following, small bold letters denote column vectors, capital bold letters matrices, and (small or capital) normal letters scalars. The vector $y$ represents the pixel values, $x$ the brightnesses, $k = k1$ the background, and $r^2 = r^21$ the read noise. $M$ is the measurement matrix, with $(M)_{ij}$ the integral over the $i$th pixel of the PSF of the $j$th site. The normalization of the PSF implies that $\sum_i (M)_{ij} = 1$ for all sites $j$. $\mathcal{P}(u)$ and $\mathcal{N}(u, v)$ denote, respectively, the multivariate Poisson distribution with mean $u$, and the multivariate normal distribution with mean $u$ and diagonal covariance matrix diag($v$). $1$ and $0$ denote the column vectors of ones and zeros.

An important parameter is the ratio of the inter-site distance ($a$) to the PSF HWHM ($r_{PSF}$), which determines

<figure><img src="image_1.png" /><figcaption>FIG. 1. **Example of test image.** The model parameters used to generate the image are given in Table I. The inset shows a zoom into the white square, with the white dots indicating the site positions. The signal-to-noise ratio of this image is about 15 dB, see Section VIII. Our reconstruction method detects the atoms with an error rate of 0.2 ± 0.1%, compared to 1.2 ± 0.2% with a standard deconvolution estimator, see Section VII A. The uncertainties represent the standard deviation over 1000 images.</figcaption></figure>

<table><caption>TABLE I. **Model parameters for the image in Fig. 1.** The unit of length is the pixel size. These parameters make a signal-to-noise ratio of about 15 dB, see Section VIII.</caption><tr><td>number of sites</td><td>N<sub>s</sub></td><td>50 × 50</td></tr><tr><td>occupancy probability</td><td>p</td><td>0.6</td></tr><tr><td>inter-site distance</td><td>a</td><td>3 pixel</td></tr><tr><td>PSF HWHM</td><td>r<sub>PSF</sub></td><td>2 pixel</td></tr><tr><td>brightness mean</td><td>μ</td><td>200 count / site</td></tr><tr><td>brightness std.</td><td>σ</td><td>20 count / site</td></tr><tr><td>background</td><td>k</td><td>0 count / pixel</td></tr><tr><td>read noise std.</td><td>r</td><td>1 count / pixel</td></tr></table>

if neighboring sites are optically resolved ($a \gg r_{PSF}$) or not ($a \lesssim r_{PSF}$). The first case is typically encountered in experiments with optical tweezer arrays, while the second characterizes most optical lattice experiments. We treat both cases in this article, although our primary focus and motivation is for optical lattices.

Fig. 1 shows one typical test image, which we will use throughout the article to illustrate our method. The parameters used to generate this image are given in Table I. They were chosen such that the reconstruction is barely possible using conventional techniques. The corresponding signal-to-noise ratio, as defined in Section VIII, is about 15 dB.

We quantify the detection accuracy of any reconstruction method using the detection error rate:

$$\text{DER} \triangleq (\text{FP} + \text{FN})/N_s, \quad (2)$$
where FP and FN stand for the number of false positives and false negatives, and the terms 'positives' and 'negatives' refer to sites labelled as occupied or empty, respectively. The sites are labelled by comparing their estimated brightness to a uniform threshold. In this work, we have determined the threshold by minimizing the DER, knowing the ground truth for each test image [39]. For an image like in Fig. 1, our reconstruction method detects the atoms with an error rate of 0.2 ± 0.1 %, compared to 1.2 ± 0.2 % with a standard deconvolution estimator, see Section VII A. The uncertainties represent the standard deviation over 1000 images.

### III. OPTIMAL LINEAR ESTIMATOR

We now turn to the description of our estimation method. It starts with the choice of a linear model to describe the relationship between the data ($\mathbf{y}$) and the variables ($\mathbf{x}$):

$$ \mathbf{y} = \mathbf{M}\mathbf{x} + \mathbf{k} + \mathbf{n}. \quad (3) $$

Here, $\mathbf{y}$, $\mathbf{x}$, $\mathbf{k}$ and $\mathbf{M}$ have the same meaning as in Eq. (1). The vector $\mathbf{n}$ accounts for the shot noise (originating from both the atomic and background signals) and for the camera read noise; it has a zero mean value. Eq. (3) is relevant in most experimental contexts. It is even exact for our test images, because the read noise is additive, and it is always possible to split a random vector with a Poisson statistics into its mean value and uncorrelated fluctuations with zero mean value. In the following, we will always subtract the background from the image, and replace $\mathbf{y} - \mathbf{k}$ with $\mathbf{y}$ for simplicity.

Ideally, the problem of estimating $\mathbf{x}$ given $\mathbf{y}$ should be treated in the maximum likelihood sense with an accurate statistical model for the signal and the noise [40]. Given the bimodal nature of the brightness probability distribution (empty sites have zero brightness, and occupied sites have a mean brightness $\mu > 0$), the resulting estimation algorithm would be non-linear, therefore time-consuming, and prone to fall into local minima.

To minimize the computation time and ensure the global convergence of the method, we want to limit ourselves to linear estimators. The optimal linear estimator (OLE), also called generalized Wiener filter, is the matrix $\mathbf{H}$ which minimizes the mean squared error between the estimated and true values of the variables:

$$ \mathbf{H}_{\text{opt}} \triangleq \arg\min_{\mathbf{H}} [\text{MSE}(\mathbf{H})], \quad (4) $$

with

$$ \text{MSE}(\mathbf{H}) \triangleq \langle \|\hat{\mathbf{x}}(\mathbf{H}) - \mathbf{x}\|^2 \rangle, \quad (5) $$

$$ \hat{\mathbf{x}}(\mathbf{H}) \triangleq \mathbf{H}(\mathbf{y} - \mathbf{M}\langle\mathbf{x}\rangle) + \langle\mathbf{x}\rangle. \quad (6) $$

Here, $\|\cdot\|$ is the Euclidean norm, and $\langle\cdot\rangle$ denotes the expectation value over the probability distributions of $\mathbf{x}$ and $\mathbf{n}$, which are treated as random variables.

We show in Appendix A that the solution to Eq. (4) can be put in the form

$$ \mathbf{H}_{\text{opt}} = (\mathbf{M}^{\top}\mathbf{\Sigma}_{n}^{-1}\mathbf{M} + \mathbf{\Sigma}_{x}^{-1})^{-1}\mathbf{M}^{\top}\mathbf{\Sigma}_{n}^{-1}, \quad (7) $$

where $\mathbf{\Sigma}_{x} \triangleq \langle(\mathbf{x} - \langle\mathbf{x}\rangle)(\mathbf{x} - \langle\mathbf{x}\rangle)^{\top}\rangle$ is the covariance matrix of the variables, and $\mathbf{\Sigma}_{n} \triangleq \langle\mathbf{n}\mathbf{n}^{\top}\rangle$ is the covariance matrix of the noise. To compute $\hat{\mathbf{x}}(\mathbf{H}_{\text{opt}})$ in practice, we solve the linear system

$$ \mathbf{A}(\hat{\mathbf{x}} - \langle\mathbf{x}\rangle) = \mathbf{b}, \quad (8a) $$

with

$$ \mathbf{A} = \mathbf{M}^{\top}\mathbf{\Sigma}_{n}^{-1}\mathbf{M} + \mathbf{\Sigma}_{x}^{-1}, \quad (8b) $$

$$ \mathbf{b} = \mathbf{M}^{\top}\mathbf{\Sigma}_{n}^{-1}(\mathbf{y} - \mathbf{M}\langle\mathbf{x}\rangle), \quad (8c) $$

using the conjugate gradient method with an incomplete LU preconditioner. See Appendix B for a brief description of the numerical implementation. The full code is available online in the form a tutorial example [31].

### IV. BAYESIAN INTERPRETATION

The OLE derived in the previous section depends on the mean and covariance matrix of $\mathbf{x}$ and $\mathbf{n}$, and therefore on the probability distributions of these quantities. The more knowledge we have about these distributions, the more accurate will the estimate be.

Let us start by using the model described in Section II, and suppose that the parameters $p$, $\mu$ and $\sigma^2$ are known. (A method to estimate these parameters using only the recorded images is presented in Section VI A.) The mean and covariance matrices can then be computed using the relations

$$ \langle(\mathbf{x})_i\rangle = p\mu, \quad (9a) $$

$$ (\mathbf{\Sigma}_x)_{ij} = [p(1-p)\mu^2 + p\sigma^2]\delta_{ij}, \quad (9b) $$

$$ (\mathbf{\Sigma}_n)_{ij} = [r^2 + k + p\mu(\mathbf{M}\mathbf{1})_i]\delta_{ij}, \quad (9c) $$

where $\delta_{ij}$ is the Kronecker delta. Here, each site is assumed to have the same occupancy probability $p$, which is as much as we can tell before an image is recorded. We call $p$ the prior probability.

Once an image is recorded, we gain some information about the probability for each site to be occupied in this particular sample: a site which appears bright is more likely to be occupied than a site which appears dark. Let us assume that we can attribute a posterior occupancy probability $(\mathbf{p})_i$ to each site $i$, depending on the apparent site brightness. (We give an explicit protocol to define the posterior occupancy probability in Section VI B). The mean and covariance matrices can then be computed using the relations

$$ \langle(\mathbf{x})_i\rangle = (\mathbf{p})_i\mu, \quad (10a) $$

$$ (\mathbf{\Sigma}_x)_{ij} = \{(\mathbf{p})_i[1 - (\mathbf{p})_i]\mu^2 + (\mathbf{p})_i\sigma^2\}\delta_{ij}, \quad (10b) $$

$$ (\mathbf{\Sigma}_n)_{ij} = [r^2 + k + \mu(\mathbf{M}\mathbf{p})_i]\delta_{ij}. \quad (10c) $$
In the following, we will refer to the estimator obtained from Eq. (9) as being a priori optimal, and to the estimator obtained from Eq. (10) as being a posteriori optimal. We will see that the a posteriori OLE is indeed more accurate than the a priori OLE.

## V. INTERPRETATION IN TERMS OF REGULARIZED LEAST-SQUARES PROBLEM

The OLE can also be interpreted as the solution to a weighted and regularized least-squares problem:

$$ \hat{x} = \arg \min_{x} \left[ (y - Mx)^{\top} \Sigma_{n}^{-1} (y - Mx) + (x - \langle x \rangle)^{\top} \Sigma_{x}^{-1} (x - \langle x \rangle) \right]. \quad (11) $$

The first term of the objective function quantifies the fidelity to the data, weighted by the variance of the noise on each pixel. The second term regularizes the problem by penalizing the solutions which spread too far away from the expectation $\langle x \rangle$. It becomes important when the PSFs of neighboring sites strongly overlap, and the standard least-squares problem is ill-posed. In the regularized least-squares method, one usually treats the regularization operator as a scalar hyperparameter, and optimizes its value based on some criterion (for instance cross-validation). Here, the Wiener filter approach directly provides us with the optimal operator.

Adopting the point of view of a regularized least-squares problem sheds new light on the difference between the a priori and a posteriori OLEs. Imagine that the posterior occupancy probability $(p)_i$ of some site $i$ is close to zero or one. For this site, the corresponding variance $(\Sigma_x)_{ii}$ becomes small, and the regularization forces the estimate $(x)_i$ to be close to the expectation $(p)_i \mu$. The result is an effective decoupling of the site $i$ from its neighbors, which in turn can be estimated more accurately.

## VI. LEARNING ABOUT THE PROBABILITY DISTRIBUTIONS FROM THE DATA

We now explain how one can determine the probability distributions of $x$ and $n$ directly from the recorded images.

### A. Prior distributions

The first thing to do is estimate $\langle x \rangle \triangleq p\mu$ from the sum of all pixel values:

$$ \langle x \rangle = \sum_{i} (y)_{i} / N_{s}. \quad (12) $$

We can then compute $\langle y \rangle = \langle x \rangle M1$, and remove it from the recorded image.

Next, we replace $\Sigma_n$ and $\Sigma_x$ by scalar matrices $\Sigma_n I$ and $\Sigma_x I$ ($I$ is the identity matrix). Within this approximation, the a priori OLE takes the simpler form $H_{opt} = (M^T M + (\Sigma_n/\Sigma_x)I)^{-1}M^T$. We treat $\Sigma_n/\Sigma_x$ as a hyperparameter, denoted $\gamma$, and tune its value to maximize the separation between empty and occupied sites in the distribution of estimated brightnesses, $\hat{x}(\gamma)$. The rationale behind this strategy comes from the equivalence with a regularized least-square problem put forward in Section V: if $\gamma$ is too small, the estimator will tend to amplify the noise, resulting in a smeared distribution; if $\gamma$ is too large, the estimator will squeeze the solution around zero. The optimal value should instead provide the most accurate estimate, and preserve the bi-modal character of the brightnesses, leading to a more accurate detection.

To quantify the separation between empty and filled sites, we use the fourth standardized moment (also called kurtosis) of $\hat{x}(\gamma)$. The kurtosis is a convenient alternative to Fisher's linear discriminant for quantifying the separation between the empty and filled sites in the distribution of estimated brightnesses, because it does not require fitting the distribution with a particular model. We have verified numerically that these two quantities are univocally related once the value of the occupancy probability $p$ is fixed.

Fig. 2 illustrates the optimization of $\gamma$, and confirms that the value which maximizes the separation between the empty and filled sites, denoted $\gamma_{opt}$, simultaneously minimizes the DER. This optimal value matches with the expected value $\Sigma_n/\Sigma_x$, with $\Sigma_x = (\Sigma_x)_{ii}$, $\Sigma_n = \sum_i (\Sigma_n)_{ii}/N_p$, and $\Sigma_x, \Sigma_n$ given by Eqs. (9b) and (9c).

Finally, we fit a Gaussian mixture model to $\hat{x}(\gamma_{opt})$:

$$ \pi(x) = (1 - \phi) \mathcal{N}(x; \mu_0, \sigma_0) + \phi \mathcal{N}(x; \mu_1, \sigma_1), \quad (13) $$

and determine $p$, $\mu$ and $\sigma^2$ (see the dashed lines in the top panel of Fig. 3). The first two parameters fall off immediately: $p$ coincides with the mixture parameter $\phi$, and $\mu$ follows from Eq. (12) [41]. For $\sigma^2$, we can assume that the variance of the Gaussian corresponding to occupied sites $(\sigma_1^2)$ is equal to the true variance $(\sigma^2)$ plus the estimation errors, which we take equal to the variance of the Gaussian corresponding to empty sites $(\sigma_0^2)$. This reasoning is valid as long as the estimation errors of different sites are not strongly correlated.

### B. Posterior distributions

To determine the posterior distributions, we need to attribute a posterior occupancy probability to each site. To this end, we exploit the Gaussian mixture model fitted to the a priori OLE, and define

$$ (p)_i \triangleq \frac{\phi \mathcal{N}[(\hat{x})_i; \mu_1, \sigma_1]}{(1 - \phi) \mathcal{N}[(\hat{x})_i; \mu_0, \sigma_0] + \phi \mathcal{N}[(\hat{x})_i; \mu_1, \sigma_1]} \cdot (14) $$
<figure><img src="image_2.png" /><figcaption>FIG. 2. Optimizing the hyperparameter of the a priori OLE. The optimal value $\gamma_{opt}$ is found by minimizing the kurtosis of $\hat{x}(\gamma)$, which simultaneously minimizes the DER. Each point is the average over 500 test images. The error bars represent the standard deviation. The same set of images is used for all values of $\gamma$. The top panel compares the distributions of $\hat{x}(\gamma)$ for a single test image, and $\gamma = \gamma_{opt}/20$, $\gamma_{opt}$, and $\gamma_{opt} \times 20$. The three plots are drawn using the same scale for each axis.</figcaption></figure>

Then, using the values of $\mu$ and $\sigma^2$ determined earlier, we can compute the first two moments of the posterior distributions using Eq. (10).

In Fig. 3, we compare the estimated brightness distributions of the *a priori* and *a posteriori* OLEs for the test image displayed in Fig. 1. Qualitatively, the distribution obtained with *a posteriori* OLE seems to be closer to the true distribution (shown in gray in the background) than with the *a priori* OLE. Quantitatively, the *a posteriori* OLE reaches $0.2 \pm 0.2\%$, compared to $0.7 \pm 0.2\%$ for the *a priori* OLE, where the uncertainty is the standard deviation over 1000 test images generated with the same parameter set. For comparison, a linear estimator based on the popular Wiener deconvolution filter (see Section VII A) gives a DER of $1.2 \pm 0.2\%$.

## VII. BENCHMARKING

We now benchmark the *a priori* and *a posteriori* OLEs against a linear estimator based on the popular Wiener deconvolution filter, which we call the ‘deconvolution estimator’.

<figure><img src="image_3.png" /><figcaption>FIG. 3. OLEs for the test image in Fig. 1. The top and bottom panels show the distribution of estimated brightnesses for the *a priori* and *a posteriori* OLEs, respectively. In both panels, the subplots on the right show the underlying occupancy probability and brightness variance, and the gray histogram in the background represents the true brightness distribution. In the top panel, the dashed lines are the two modes of the most likely Gaussian mixture distribution, with weights, means and variances $(1-\phi)\mu_0, \sigma_0^2$ (empty sites), and $\phi, \mu_1, \sigma_1^2$ (filled sites), see Eq. (13).</figcaption></figure>

### A. Deconvolution estimator

The deconvolution estimator proceeds as follows: The image is first centered around zero by removing the mean pixel value, and filtered using the Wiener deconvolution filter with a uniform (in reciprocal space) regularization operator $\lambda \in \mathbb{R}^+$.

To extract the site brightnesses from the filtered image, we perform the convolution by a disk-shaped kernel of radius $d \in \mathbb{R}^+$, and compute the linear interpolation of the result at the site coordinates. Since the Wiener deconvolution filter, the convolution, and the linear interpolation are all linear operations, the deconvolution estimator as a whole is a linear estimator.

For each set of parameters used in the article, we have performed a joint optimization of the hyperparameters $\lambda$ and $d$ by maximizing the kurtosis of $\hat{x}(\lambda, d)$, as discussed in Section VI A.
<figure><img src="image_4.png" /><figcaption>FIG. 4. Detection error rate and signal-to-noise ratio. Top panel: The colored lines are the DER level curves at 0.1% as a function of $\mu$ and $a$, with the other parameters as in Table I. Each point is averaged over 100 test images. The gray lines are the SNR level curves at 10, 15, 20, 25 and 30 dB, see Section VIII. Lower panel: Cuts through the upper panel along the lines $\mu = 1000$ (A), and $a = 1.5 r_{\text{PSF}}$ (B). The SNR is shown in gray.</figcaption></figure>

## B. Detection error rate

Fig. 4 shows the level curves at 0.1% of each estimator in the space spanned by $\mu$ and $a$, with the other parameters fixed to the values of Table I. One can clearly distinguish two regimes. When the sites are optically resolved ($a \gg r_{\text{PSF}}$), all estimators have exactly the same accuracy. In this regime, the DER is limited by the fluctuations of the measured signal—due to the variability of the number of scattered photons, or the noise—, rather than by the ability of the estimator to correctly estimates the brightnesses. In contrast, when the sites are not optically resolved ($a \lesssim r_{\text{PSF}}$), the detection accuracies arrange as follows: *a posteriori* OLE > *a priori* OLE > deconvolution, and our method drastically improves the accuracy compared to the standard method.

## C. Runtime

Fig. 5 shows the runtime (wall-clock time) of the different estimators for the parameters of Table I and an increasing number of sites. While the deconvolution

<figure><img src="image_5.png" /><figcaption>FIG. 5. Runtime. The runtime is the averaged over 5 test images generated with the parameters of Table I, except for the varying number of sites.</figcaption></figure>

method is always faster, the runtime of the OLEs approximately linearly with the number of sites, and remains well below 100 ms for an array of 100 × 100 sites, which makes them both perfectly suitable for real-time usage in all existing experiments.

The runtimes were measured on a desktop computer with an Intel Core i9-10900K CPU. Note that the optimization of the hyperparameters and the computation of the Gram matrix are *not* included, because these objects can be cached and reused for all images recorded under similar experimental conditions. Note also that the runtime of the *a posteriori* OLE includes that of the *a priori* OLE, since the *a priori* OLE is required to determine the posterior occupancy probability.

## D. Robustness against calibration errors

All estimators used or proposed so far for the reconstruction of the site occupancy assume that the site coordinates and the PSF have been correctly calibrated beforehand. In this last section, we study how robust the OLEs and the deconvolution estimator are against two possible calibration errors: an error in the site positions (global offset along one axis of the site array) and an error in the PSF HWHM (global scaling factor). By calibration error, we mean that the value of the parameter used to compute the estimate is not equal to the value of the parameter which was used to generate the test image.

We have computed the DER of each estimator for both types of calibration errors and varying error amplitudes, see Fig. 6. One sees that the OLEs display a comparable or even better robustness against the calibration errors compared to the deconvolution estimator. They should therefore perform well under realistic experimental conditions.
<figure><img src="image_6.png" /><figcaption>FIG. 6. Robustness against calibration errors. In the left panel, the coordinates of the sites are affected by a global offset along one axis. The offset is given relative to the inter-site distance a. In the right panel, the PSF HWHM is affected by a uniform scaling factor. Each point is an average over 500 test images generated with the parameters of Table I. The error bars represent the standard deviation over the test images.</figcaption></figure>

VIII. SIGNAL-TO-NOISE RATIO

A valuable by-product of our approach is a rigorous definition of a signal-to-noise ratio for the estimation-detection problem, namely

$$ \text{SNR} \triangleq 10 \log_{10} \left( \frac{N_s \mu^2}{\text{MSE}(\mathbf{H}_{\text{opt}})} \right), \quad (15) $$

with (see Appendix A)

$$ \text{MSE}(\mathbf{H}_{\text{opt}}) = \text{trace}[(\mathbf{I} - \mathbf{H}_{\text{opt}}\mathbf{M})\mathbf{\Sigma}_x]. \quad (16) $$

In the limit $a \gg r_{\text{PSF}}$, the weighted Gram matrix $\mathbf{M}^\top\mathbf{\Sigma}_n^{-1}\mathbf{M}$ is diagonal, and no longer depends on $a/r_{\text{PSF}}$. The optimal MSE then reduces to

$$ \text{MSE}(\mathbf{H}_{\text{opt}}) = \sum_{i=1}^{N_s} \left[ 1 - \frac{(\mathbf{M}^\top\mathbf{\Sigma}_n^{-1}\mathbf{M})_{ii}(\mathbf{\Sigma}_x^{-1})_{ii}}{(\mathbf{M}^\top\mathbf{\Sigma}_n^{-1}\mathbf{M})_{ii} + (\mathbf{\Sigma}_x^{-1})_{ii}} \right]. \quad (17) $$

Fig. 4 shows level curves of the SNR at 10, 15, 20, 25 and 30 dB (gray lines) alongside the DER level curves. This graph confirms that the SNR is a good predictor for the DER of both OLEs as their DER level curves are superimposed with SNR level curves in the regimes $a \gg r_{\text{PSF}}$ and $a \lesssim r_{\text{PSF}}$. In-between these two regimes, a transition occurs.

Being able to predict the DER based on the SNR will be extremely useful to design new experiments as it gives the guarantee of reaching a target detection accuracy without the need to over-design the imaging system.

IX. CONCLUSION

We have used the generalized Wiener filter approach to solve the problem of individual atom detection in arrays

of optical microtraps. By design, this method is expected to be the best possible linear estimator for the site brightness. Here, we have shown that it indeed surpasses the standard method based on Wiener deconvolution, drastically reducing the detection error rate when the sites are not optically resolved.

In contrast with maximum likelihood estimators, the computation time remains fully compatible with real-time applications. It is also robust against calibration errors. Last but not least, it offers the possibility (not illustrated in this article) to easily account for spatial inhomogeneities in the distribution of atoms, site brightness, read noise, or PSF. For instance, one can easily disregard ‘hot’ camera pixels by removing the corresponding lines of the measurement matrix $\mathbf{M}$, or account for the variations of the PSF across the field of view when the latter is so wide that optical aberrations cannot be neglected.

Benchmarking against supervised or unsupervised machine learning remains to be done. Artificial neural networks have the ability to precisely model the bimodal brightness probability distribution, but their versatility comes at the cost of a delicate training phase. On the contrary, our method depends on a few hyperparameters only; these hyperparameters have a clear physical interpretation, and they can be optimized using only the images to be analyzed.

Finally, the Wiener filter theory provides us with a rigorous definition of a signal-to-noise ratio for the estimation-detection problem. This SNR aggregates all experimental parameters in a single number, and completely determines the detection error rate. This opens the prospect for the joint optimization of optical elements and image processing algorithms to improve the global performance of the imaging system in future experiments [42, 43]. Such approach should prove especially useful for experiments aiming at increasing the number of sites or reducing the inter-site distance beyond what is currently feasible.

ACKNOWLEDGMENTS

The authors thank Félix Faisant, Thomas Chalopin, David Clément, and Isabelle Bouchoule for their feedback on the manuscript, and Anaïs Molineri for early work on the subject.

Appendix A: Derivation of the OLE

In this section we explain how we have derived the expression (7) of the OLE. We take $\langle \mathbf{x} \rangle = \mathbf{0}$ for simplicity. We start by writing the squared norm as the trace of the
outer product and expand all terms:

$$ \text{MSE}(\mathbf{H}) = \langle \text{trace}[(\mathbf{H}\mathbf{y} - \mathbf{x})(\mathbf{H}\mathbf{y} - \mathbf{x})^{\top}] \rangle \quad (A1) $$

$$ = \text{trace}[\mathbf{H}\langle\mathbf{y}\mathbf{y}^{\top}\rangle\mathbf{H}^{\top} - \mathbf{H}\langle\mathbf{y}\mathbf{x}^{\top}\rangle - \langle\mathbf{x}\mathbf{y}^{\top}\rangle\mathbf{H}^{\top} + \langle\mathbf{x}\mathbf{x}^{\top}\rangle]. \quad (A2) $$

The OLE is such that the derivative of the MSE with respect to $\mathbf{H}$ is identically zero, which gives [44, Eqs. (100), (104) and (111)]:

$$ \left. \frac{\mathrm{d}\mathrm{MSE}}{\mathrm{d}\mathbf{H}} \right|_{\mathbf{H}_{\mathrm{opt}}} = 2[\mathbf{H}_{\mathrm{opt}}\langle\mathbf{y}\mathbf{y}^{\top}\rangle - \langle\mathbf{x}\mathbf{y}^{\top}\rangle] = \mathbf{0}, \quad (A3) $$

hence

$$ \mathbf{H}_{\mathrm{opt}} = \langle \mathbf{x}\mathbf{y}^{\top} \rangle (\langle \mathbf{y}\mathbf{y}^{\top} \rangle)^{-1} . \quad (A4) $$

We proceed further by replacing the image vector $\mathbf{y}$ with $\mathbf{M}\mathbf{x} + \mathbf{n}$:

$$ \langle \mathbf{y}\mathbf{x}^{\top} \rangle = \langle (\mathbf{M}\mathbf{x} + \mathbf{n})\mathbf{x}^{\top} \rangle \quad (A5) $$

$$ = \mathbf{M}\langle\mathbf{x}\mathbf{x}^{\top}\rangle, \quad (A6) $$

and

$$ \langle \mathbf{y}\mathbf{y}^{\top} \rangle = \langle (\mathbf{M}\mathbf{x} + \mathbf{n})(\mathbf{M}\mathbf{x} + \mathbf{n})^{\top} \rangle \quad (A7) $$

$$ = \mathbf{M}\langle\mathbf{x}\mathbf{x}^{\top}\rangle\mathbf{M}^{\top} + \langle\mathbf{n}\mathbf{n}^{\top}\rangle, \quad (A8) $$

In deriving these results, we have used the noise properties $\langle\mathbf{x}\mathbf{n}^{\top}\rangle = \langle\mathbf{x}\rangle\langle\mathbf{n}^{\top}\rangle$ ($\mathbf{x}$ and $\mathbf{n}$ are uncorrelated), and $\langle\mathbf{n}\rangle = \mathbf{0}$.

The matrices $\langle\mathbf{x}\mathbf{x}^{\top}\rangle$ and $\langle\mathbf{n}\mathbf{n}^{\top}\rangle$ are the covariance matrices of $\mathbf{x}$ and $\mathbf{n}$, denoted $\Sigma_x$ and $\Sigma_n$ in the main text. Combining Eqs. (A4), (A6) and (A8) gives

$$ \mathbf{H}_{\mathrm{opt}} = \Sigma_x \mathbf{M}^{\top} (\mathbf{M} \Sigma_x \mathbf{M}^{\top} + \Sigma_n)^{-1} . \quad (A9) $$

This expression turns out to be inconvenient for efficiently computing the solution to the problem. We obtain the more useful expression in Eq. (7) by applying the so-called Woodbury identity [44, Eq. (158)] (which holds because $\Sigma_x$ and $\Sigma_n$ are both positive definite by definition).

Finally, we can also compute the MSE corresponding to the OLE by combining Eqs. (A2), (A6), (A8) and (A9):

$$ \text{MSE}(\mathbf{H}_{\mathrm{opt}}) = \text{trace}[\mathbf{H}_{\mathrm{opt}}(\mathbf{M}\Sigma_x\mathbf{M}^{\top} + \Sigma_n)\mathbf{H}_{\mathrm{opt}}^{\top} - \mathbf{H}_{\mathrm{opt}}\mathbf{M}\Sigma_x - \Sigma_x\mathbf{M}^{\top}\mathbf{H}_{\mathrm{opt}}^{\top} + \Sigma_x] \quad (A10) $$

$$ = \text{trace}[(\mathbf{I} - \mathbf{H}_{\mathrm{opt}}\mathbf{M})\Sigma_x]. \quad (A11) $$

## Appendix B: Numerical implementation

We solve the linear system defined by Eq. (8) using the preconditioned conjugate gradient (CG) method, which is an iterative algorithm very well suited to large, sparse linear least-square problems. Here, the sparsity of the matrix $\mathbf{A}$ follows from the fact that the sites are only coupled to their closest neighbors, and the covariance matrices $\Sigma_x$ and $\Sigma_n$ are assumed to be diagonal.

The use of a preconditioner is crucial to ensure fast and reliable convergence of the CG algorithm because the condition number of the matrix $\mathbf{A}$ quickly increases with the ratio $a/r_{\mathrm{PSF}}$. We have tested different preconditioners (diagonal, incomplete Cholesky and incomplete LU decompositions), and found the Crout version of the incomplete LU (ILU) decomposition to be the most efficient and reliable. We adapt the number of nonzero entries to keep in the LU matrices to the ratio $a/r_{\mathrm{PSF}}$ by checking the convergence of the CG algorithm. The CG algorithm is stopped when $\|\mathbf{A}\mathbf{x} - \mathbf{b}\|/\|\mathbf{b}\| < 10^{-2}$ for all values of $a/r_{\mathrm{PSF}}$.

As the covariance matrix of the noise depends on the occupancy probability (see Eqs. (9c) and (10c)), one should in theory compute the weighted Gram matrix $\mathbf{M}^{\top}\Sigma_n^{-1}\mathbf{M}$ using either the prior or the posterior occupancy probability. In practice, we found that both the a priori and a posteriori OLEs obtained with the scalar approximation $\Sigma_n = \Sigma_n\mathbf{I}$, $\Sigma_n = \sum_i(\Sigma_n)_{ii}/N_p$, are almost indistinguishable from those obtained without this approximation. We therefore compute only the non-weighted Gram matrix $\mathbf{M}^{\top}\mathbf{M}$, and use it for both the a priori and a posteriori OLEs, which saves computation time.

The numerical implementation is in Python, and relies mainly on the standard libraries NumPy [45] and SciPy [46]. The Wiener deconvolution filter is from the scikit-image library [47]. The computation of the incomplete LU decomposition uses the `ilupp` packages [48], which itself is based on the `ILU++` software [49]. If this package is not available, one may instead rely on the routine provided by Scipy, which itself uses the SuperLU library [50]. We use the Intel Math Kernel Library to speed up sparse matrix-matrix products, with Python bindings provided by the `sparse dot mkl` package [51].

The full code is available online in the form of a tutorial notebook [31].

[1] A. Browaeys and T. Lahaye, Many-body physics with individually controlled rydberg atoms, Nature Physics 16, 132 (2020).

[2] L. Henriet, L. Beguin, A. Signoles, T. Lahaye, A. Browaeys, G.-O. Reymond, and C. Jurczak, Quantum computing with neutral atoms, Quantum 4, 327 (2020).

[3] C. Gross and W. S. Bakr, Quantum gas microscopy for single atom and spin detection, Nature Physics 17, 1316 (2017).
(2021).

[4] A. M. Kaufman and K.-K. Ni, Quantum science with optical tweezer arrays of ultracold atoms and molecules, Nature Physics **17**, 1324 (2021).

[5] S. M. Kay, *Fundamentals of statistical signal processing, Vol. I: Estimation Theory* (Prentice-Hall, 1993).

[6] S. M. Kay, *Fundamentals of statistical signal processing, Vol. II: Detection Theory* (Prentice-Hall, 1998).

[7] M. Kwon, M. F. Ebert, T. G. Walker, and M. Saffman, Parallel low-loss measurement of multiple atomic qubits, Physical Review Letters **119**, 180504.

[8] M. A. Norcia, A. W. Young, and A. M. Kaufman, Microscopic Control and Detection of Ultracold Strontium in Optical-Tweezer Arrays, Physical Review X **8**, 041054 (2018).

[9] A. Cooper, J. P. Covey, I. S. Madjarov, S. G. Porsev, M. S. Safronova, and M. Endres, Alkaline-Earth Atoms in Optical Tweezers, Physical Review X **8**, 041055 (2018).

[10] I. S. Madjarov, Entangling, controlling, and detecting individual strontium atoms in optical tweezer arrays (2021).

[11] W. S. Bakr, J. I. Gillen, A. Peng, S. Folling, and M. Greiner, A quantum gas microscope for detecting single atoms in a Hubbard-regime optical lattice, Nature **462**, 74 (2009).

[12] M. F. Parsons, F. Huber, A. Mazurenko, C. S. Chiu, W. Setiawan, K. Wooley-Brown, S. Blatt, and M. Greiner, Site-Resolved Imaging of Fermionic $^6$Li in an Optical Lattice, Physical Review Letters **114**, 213002 (2015).

[13] A. Omran, M. Boll, T. A. Hilker, K. Kleinlein, G. Salomon, I. Bloch, and C. Gross, Microscopic Observation of Pauli Blocking in Degenerate Fermionic Lattice Gases, Physical Review Letters **115**, 263001 (2015).

[14] A. Alberti, C. Robens, W. Alt, S. Brakhane, M. Karski, R. Reimann, Artur Widera, and D. Meschede, Super-resolution microscopy of single atoms in optical lattices, New Journal of Physics **18**, 053010 (2016).

[15] J. F. Sherson, C. Weitenberg, M. Endres, M. Cheneau, I. Bloch, and S. Kuhr, Single-atom resolved fluorescence imaging of an atomic Mott insulator, Nature **467**, 68 (2010).

[16] P. Schauss, *High-resolution imaging of ordering in Rydberg many-body systems*, Ph.D. thesis, Ludwig-Maximilians-Universität München (2015).

[17] M. Miranda, R. Inoue, Y. Okuyama, A. Nakamoto, and M. Kozuma, Site-resolved imaging of ytterbium atoms in a two-dimensional optical lattice, Physical Review A **91**, 063414 (2015).

[18] M. F. Parsons, *Probing the Hubbard Model With Single-Site Resolution*, Ph.D. thesis, Harvard University (2016).

[19] L. W. Cheuk, *Quantum Gas Microscopy of Strongly Correlated Fermions*, Ph.D. thesis, Massachusetts Institute of Technology (2017).

[20] R. Yamamoto, H. Ozawa, D. C. Nak, I. Nakamura, and T. Fukuhara, Single-site-resolved imaging of ultracold atoms in a triangular optical lattice, New Journal of Physics **22**, 123028 (2020).

[21] K. Kwon, K. Kim, J. Hur, S. Huh, and J.-y. Choi, Site-resolved imaging of a bosonic Mott insulator of $^7$Li atoms, Physical Review A **105**, 033323 (2022).

[22] A. La Rooij, C. Ulm, E. Haller, and S. Kuhr, A comparative study of deconvolution techniques for quantum-gas microscope images, New Journal of Physics **25**, 083036 (2023).

[23] J. Mongkolkiattichai, L. Liu, D. Garwood, J. Yang, and P. Schauss, Quantum gas microscopy of fermionic triangular-lattice mott insulators, Physical Review A **108**, L061301 (2023).

[24] S. Buob, J. Höschele, V. Makhalov, A. Rubio-Abadal, and L. Tarruell, A strontium quantum-gas microscope, PRX Quantum **5**, 020316 (2024).

[25] M. Martinez-Dorantes, W. Alt, J. Gallego, S. Ghosh, L. Ratschbacher, Y. Völzke, and D. Meschede, Fast non-destructive parallel readout of neutral atom registers in optical potentials, Physical Review Letters **119**, 180503 (2017).

[26] A. H. Burrell, D. J. Szwer, S. C. Webster, and D. M. Lucas, Scalable simultaneous multiqubit readout with 99.99% single-shot fidelity, Physical Review A **81**, 040302 (2010).

[27] L. R. B. Picard, M. J. Mark, F. Ferlaino, and R. v. Bijnen, Deep learning-assisted classification of site-resolved quantum gas microscope images, Measurement Science and Technology **31**, 025201 (2019).

[28] J. Verstraten, K. Dai, M. Dixmerias, B. Peaudecerf, T. d. Jongh, and T. Yefsah, In-situ imaging of a single-atom wave packet in continuous space, arXiv:2404.05699 [quant-ph] (2024).

[29] A. Impertro, J. F. Wienand, S. Häfele, H. von Raven, S. Hubele, T. Klostermann, C. R. Cabrera, I. Bloch, and M. Aidelsburger, An unsupervised deep learning algorithm for single-site reconstruction in quantum gas microscopes, Communications Physics **6**, 1 (2023).

[30] J. Winklmann, A. Alberti, and M. Schulz, Comparison of atom detection algorithms for neutral atom quantum computing, in *2024 IEEE International Conference on Quantum Computing and Engineering (QCE)*, Vol. 01 (2024) pp. 1048–1057.

[31] Available at https://gitlab.in2p3.fr/gaz-quantiques-lcf/strontium/atom-detection.

[32] F. Gyger, M. Ammenwerth, R. Tao, H. Timme, S. Snigirev, I. Bloch, and J. Zeiher, Continuous operation of large-scale atom arrays in optical lattices, Physical Review Research **6**, 033104 (2024).

[33] G. Pichard, D. Lim, E. Bloch, J. Vanecloo, L. Bourachot, G.-J. Both, G. Mériaux, S. Dutartre, R. Hostein, J. Paris, B. Ximenez, A. Signoles, A. Browaeys, T. Lahaaye, and D. Dreon, Rearrangement of individual atoms in a 2000-site optical-tweezer array at cryogenic temperatures, Physical Review Applied **22**, 024073 (2024).

[34] N.-C. Chiu, E. C. Trapp, J. Guo, M. H. Abobeih, L. M. Stewart, S. Hollerith, P. Stroganov, M. Kalinowski, A. A. Geim, S. J. Evered, S. H. Li, L. M. Peters, D. Bluvstein, T. T. Wang, M. Greiner, V. Vuletić, and M. D. Lukin, Continuous operation of a coherent 3,000-qubit system, arXiv:2506.20660 [quant-ph] (2025).

[35] B. Hofer, D. Bloch, G. Biagioni, N. Bonvalet, A. Browaeys, and I. Ferrier-Barbut, Single-atom resolved collective spectroscopy of a one-dimensional atomic array, arXiv:2412.02541 [quant-ph] (2024).

[36] L. Su, A. Douglas, M. Szurek, R. Groth, S. F. Ozturk, A. Krahn, A. H. Hébert, G. A. Phelps, S. Ebadi, S. Dickerson, F. Ferlaino, O. Marković, and M. Greiner, Dipolar quantum solids emerging in a hubbard quantum simulator, Nature **622**, 724 (2023).

[37] More elaborate models exist to simulate specific experimental contexts [52], but their number of parameters
make them less practical to use.

[38] We have also tested our method with an Airy disk, and obtained very similar results. One can convert the Gaussian HWHM to the radius of the Airy disk by multiplying the former with a factor 2.4, assuming that both functions have the same HWHM.

[39] In an experimental context, one would rather determine the threshold by first fitting the distribution of estimated brightnesses with a mixture model, and then find the brightness corresponding to an equal likelihood of being empty or filled.

[40] G. Demoment and Y. Goussard, Inversion within the probabilistic framework, in *Bayesian Approach to Inverse Problems* (John Wiley & Sons, Ltd, 2008) Chap. 3, pp. 59–78.

[41] Eq. (12) provides a more accurate estimate of $\mu$ than the Gaussian mixture parameter $\mu_1$ because of the 'squeezing' effect of the regularization.

[42] F. Diaz, F. Goudail, B. Loiseaux, and J.-P. Huignard, Increase in depth of field taking into account deconvolution by optimization of pupil mask, Opt. Lett. **34**, 2970 (2009).

[43] M.-A. Burcklen, H. Sauer, F. Diaz, and F. Goudail, Joint digital-optical design of complex lenses using a surrogate image quality criterion adapted to commercial optical design software, Applied Optics **57**, 9005 (2018).

[44] K. B. Petersen and M. S. Pedersen, The matrix cookbook (2012), version 20121115.

[45] C. R. Harris, K. J. Millman, S. J. van der Walt, R. Gommers, P. Virtanen, D. Cournapeau, E. Wieser, J. Taylor, S. Berg, N. J. Smith, R. Kern, M. Picus, S. Hoyer, M. H. van Kerkwijk, M. Brett, A. Haldane, J. F. del Río, M. Wiebe, P. Peterson, P. Gérard-Marchant, K. Shep-

pard, T. Reddy, W. Weckesser, H. Abbasi, C. Gohlke, and T. E. Oliphant, Array programming with NumPy, Nature **585**, 357 (2020).

[46] P. Virtanen, R. Gommers, T. E. Oliphant, M. Haberland, T. Reddy, D. Cournapeau, E. Burovski, P. Peterson, W. Weckesser, J. Bright, S. J. van der Walt, M. Brett, J. Wilson, K. J. Millman, N. Mayorov, A. R. J. Nelson, E. Jones, R. Kern, E. Larson, C. J. Carey, I. Poblat, Y. Feng, E. W. Moore, J. VanderPlas, D. Laxalde, J. Perktold, R. Cimrman, I. Henriksen, E. A. Quintero, C. R. Harris, A. M. Archibald, A. H. Ribeiro, F. Pedregosa, P. van Mulbregt, and SciPy 1.0 Contributors, SciPy 1.0: Fundamental Algorithms for Scientific Computing in Python, Nature Methods **17**, 261 (2020).

[47] S. van der Walt, J. L. Schönberger, J. Nunez-Iglesias, F. Boulogne, J. D. Warner, N. Yager, E. Gouillart, T. Yu, and the scikit-image contributors, scikit-image: image processing in Python, PeerJ **2**, e453 (2014-06).

[48] C. Hofreither, ilupp – ILU algorithms for C++ and Python, https://github.com/c-f-h/ilupp.

[49] J. Mayer, ILU++: A new software package for solving sparse linear systems with iterative methods, PAMM **7**, 2020123 (2007).

[50] X. S. Li, An overview of SuperLU: Algorithms, implementation, and user interface, ACM Trans. Math. Softw. **31**, 302 (2005).

[51] C. A. Jackson, Python wrapper for Intel Math Kernel Library (MKL) matrix multiplication, https://github.com/flatironinstitute/sparse Dot.

[52] J. Winklmann, D. Tsevas, and M. Schulz, Realistic neural atom image simulation, in *2023 IEEE International Conference on Quantum Computing and Engineering (QCE)*, Vol. 01 (2023) pp. 1349–1359.