Vortex detection in atomic Bose–Einstein condensates using neural networks trained on synthetic images

We experimentally produce BECs with multiple vortices by rapidly cooling an atomic thermal gas, causing it to undergo a phase transition into the superfluid phase. During the quench process, the superfluid order parameter of the system develops independently in local areas, leading to the emergence of topological defects at the interfaces between these phase domains. This phenomenon is called the Kibble–Zurek mechanism [49, 50]. In the case of a superfluid system, these topological defects appear as quantum vortices.

The experiment begins with the preparation of a thermal gas of ⁸⁷Rb atoms in an optical dipole trap (ODT) with a highly oblate and elongated geometry, as described in [51]. The trapping potential of the ODT is expressed as

$$\begin{equation} V(x,y,z) = \frac{1}{2}m\omega^2 R_x^2 \bigg[ ( \frac{x}{R_x} )^2 +( \frac{y}{R_y})^{3.9} + ( \frac{z}{R_z} )^2 \bigg], \end{equation} \tag{ 1 }$$

where $\omega \approx 2\pi \times 7$ Hz is the trapping frequency for the x direction, m is the mass of a ⁸⁷Rb atom, and $\{R_x, R_y, R_z\} \approx \{65, 244, 2.8\}\,\mu\textrm{m}$ are the Thomas–Fermi radii of the final condensate. The ODT provides strong confinement along the z direction and this highly oblate geometry causes the vortex lines to be energetically aligned along the tightly confining z-axis.

We rapidly cool the thermal gas by decreasing the trap depth of the ODT from $1.15 U_c$ to $0.27 U_c$, where U_c is the critical trap depth at which the sample starts to undergo a phase transition. During the quench, the temperature of the sample is regulated by the depth of the trap through evaporation. To facilitate the formation of quantum vortices, we keep the sample in the ODT for an additional 1.25 s after the quench. We then release the trapping potential to allow the BEC to expand. During this expansion, the vortex cores expand faster, becoming visible in the subsequent imaging. After a 40 ms ToF expansion, we perform absorption imaging of the BEC along the z-axis to capture the density distribution of the BEC, which allows us to detect the presence and spatial arrangement of quantum vortices. The cooling speed is controlled by the decrease time of the trap depth, which ranges from 0.8 s to 11 s, resulting in a variable vortex number in the BEC. For our fastest quench, the vortex number exceeds 60 [30, 31]. The typical number of atoms in the condensate is approximately $N = 9.6\times 10^6$.

In figure 1(a), we display examples of experimental images in which the vortex cores are identified as localized areas of reduced density. The images are presented as 2D optical depth (OD) arrays and have been cropped to a size of 120 × 40 pixels, with each pixel being approximately 4.3 µm in size. This cropping ensures that the BEC sample is centered and enclosed within the image. To remove any outlier pixels that may have been caused by noise during the calculation of OD, we trim values that differ by more than four standard deviations from the surrounding values, replacing them with the average value of neighboring pixels. We then normalize the images so that the minimum pixel value is set to zero and the maximum is set to one, creating a consistent scale for analysis.

Zoom In Zoom Out Reset image size

The dynamics of weakly interacting atomic BECs can be effectively described by the GPE [4]. This equation takes the form

$$\begin{equation} \mathrm i\hbar \frac{\partial}{\partial t} \psi(x,y,z,t) = \bigg[ -\frac{\hbar^2}{2m}\nabla^2 + V(x,y,z) + \frac{4\pi \hbar^2 N a_\mathrm s}{m} \psi ^2 \bigg] \psi(x,y,z,t), \end{equation} \tag{ 2 }$$

where ψ is the condensate wave function and $a_\mathrm s$ is the scattering length of the atoms. The BEC dynamics in the z direction is highly suppressed due to tight confinement, so we simplify the equation by integrating out the z-dependent part of the wave function, resulting in an effective equation for a 2D complex wave function $\psi(x,y)$. The numerical solution of the 2D GPE is obtained on a grid space of 1024 × 512 with a cell size of 0.54 µm, which has a much higher resolution than our experimental images. The healing length of the condensate in the experiment was 0.26 µm and such high resolution is necessary for the GPE simulation. After determining the complex wave function of the system from the GPE, we calculate the density profile of the BEC by taking the squared norm of the wave function $|\psi(x,y)|^2$.

To create synthetic images of BECs with quantum vortices, we first calculate the ground state of the system using the imaginary-time evolution of the GPE [52, 53]. Then, we imprint phase singularities onto the condensate wave function. The phase singularity mask is given by $\theta(x,y) = \sum_{i = 1}^{N_v} s_i \arctan(x-x_i,y-y_i)$, where $(x_i,y_i)$ is the position of the ith phase singularity and $s_i = \pm 1$ is the sign of phase winding randomly chosen. The number of singularities, N_v , imprinted on the condensate ranges from 0 to 100, and their locations are randomly distributed throughout the sample area. After phase imprinting, we evolve the condensate with the new phase profile for a short hold time using the real-time evolution of the GPE, during which phase singularities develop into vortices in the system. Then, to account for the effects of the ToF imaging process, where the density profile around the vortex cores are slightly altered, we allow the sample to expand by suddenly releasing the trapping potential. In the free expansion, we gradually reduce the interaction strength for the same ToF time to simulate the effect due to the fast expansion along the z direction [5, 6]. Finally, we create a density image through the squared norm of the wavefunction. Figure 1(b) shows typical simulated BEC images.

After obtaining a synthetic BEC image, we label it with the positions of the vortices in the simulated BEC sample. We know the initial positions of the phase singularities on the condensate, but the vortices may have shifted and dissipated during the evolution after the phase-imprinting sequence. To identify the vortex positions, we use a brute-force algorithm based on the phase information of the wave function [45]. The algorithm searches for all local density minima as potential vortices and calculates the integral of the phase gradient $\nabla (\textrm{arg}(\psi) )$ along the circumference of each candidate point. If the integral is non-zero, the point is identified as a vortex. To ensure that only the vortices within the BEC are accurately labeled, we define a sample area criterion using a cutoff density, denoted by ρ_cut, and exclude any singularity in the outer region of the sample with density $\rho \unicode{x2A7D} \rho_\textrm{cut}$ from the labeled vortex list. Figure 2(a) shows the change of the sample area enclosed by the contour ($\rho = \rho_\textrm{cut}$) as the value of ρ_cut is adjusted. Here, the sample density is normalized so that the maximum density, $\textrm{max}( |\psi(x,y)|^2)$, is set to one.

Zoom In Zoom Out Reset image size

The last step in creating a synthetic training set is to pre-process the simulated images to make them look like the ones obtained in the experiment. Initially, the simulated images are cropped to 960 × 320 pixels and then downsized to 120 × 40 pixels, so that they have the same resolution and size as the experimental images. After that, Gaussian noise is added to each pixel of the image with a uniform noise strength across the image. The Gaussian noise has a zero mean and a standard deviation denoted as σ_n. By adjusting the tuning parameter σ_n, the noise level of the images can be changed, as shown in figure 2(b). Finally, the image is normalized to set the minimum and maximum pixel values to zero and one, respectively, completing the pre-processing sequence (figure 1(b)).

Our vortex detection algorithm is based on the CNN model proposed by Metz et al [45]. This CNN architecture consists of seven convolutional layers and three maxpool layers. The stride of the layers is 1, except for the first two maxpool layers, which have a stride of 2 to reduce the output size. The CNN takes greyscale density images as a single input channel, and it has three output channels, each reduced both in height and width by a factor of four compared to the input channel. Each 4 × 4 grid cell of the input, serving as the unit window for vortex detection, is coarsened into a set of three values for the output channels, which represent the probability of a vortex core being in the grid cell and the rescaled x and y positions of the vortex core within the grid cell, respectively. When the size of the input image is $\mathrm{H} \times \mathrm{W}$, the resulting output tensor has dimensions of $\mathrm{H}/4 \times \mathrm{W}/4 \times 3$, which are denoted as $\widetilde{Y}_{ijk}$, where i and j are the position indices of a 4 × 4 grid cell, and k is the index of the output channels. For ground-truth data of the vortex positions, the array Y_ijk is created, where $Y_{ij1} = 1$ when a grid cell contains a vortex and $Y_{ij1} = 0$ when a grid cell contains no vortex. The loss function used for training the CNN is defined as follows:

$$\begin{align} L = \sum_{\mathrm{batch}}\sum_{i,j} \Bigg[ -w_{1}Y_{ij1}\log(\widetilde{Y}_{ij1}) -(1-Y_{ij1})\log(1-\widetilde{Y}_{ij1}) + w_{2}Y_{ij1} ( {(Y_{ij2}-\widetilde{Y}_{ij2})}^2+{(Y_{ij3}-\widetilde{Y}_{ij3})}^2 ) \Bigg], \end{align} \tag{ 3 }$$

where w₁ and w₂ are hyperparameters that control the weighting of each term. In this study, w₁ and w₂ are set to 10, as in [45].

We apply upsampling to the experimental and synthetic images to feed them to the CNN. Nearest-neighbor upsampling is used, where each pixel in the original image is repeated $n_\textrm{up}$ times in both the row and column directions without any interpolation. The original size of the images is 120 × 40 pixels, and after upsampling, they become $120n_\textrm{up}\times40n_\textrm{up}$ in size. The value of $n_\textrm{up}$ can be adjusted to change the resolution of the images, but it should be noted that the images with different $n_\textrm{up}$ values look the same because no interpolation is used during the upsampling process.

The CNN is trained on a dataset of 2000 synthetic images, with a batch size of 100. The training process is conducted for 300 epochs, with an initial learning rate of $\eta = 10^{-3}$. After 300 epochs, the learning rate is decreased to 10⁻⁴. To ensure a fair comparison of the performance of CNNs trained with various parameter combinations, the early stopping method is used with a patience value of 50 epochs, without fixing the training epoch. This means that, while monitoring the loss function, if it does not show further reduction over the course of 50 epochs, the algorithm regards it as having reached its minimum and concludes the training process.

The performance of the CNN is evaluated by comparing the number of vortices identified by the machine to the number counted by a human observer. To quantify the discrepancy between the machine and human counts, we calculate the root mean squared error (RMSE) value that is given by

$$\begin{equation} \mathrm{RMSE} = \bigg( \frac{\sum_{i}^{n}{[ N_{\mathrm{m}}(i)-N_{\mathrm{h}}(i) ]^2}}{n} \bigg)^{1/2}, \end{equation} \tag{ 4 }$$

where $N_\textrm{m}(i)$ is the machine count for the ith image, $N_\textrm{h}(i)$ is the corresponding human count, and n = 350 is the number of experimental images used for the test. Out of the 350 experimental images, the distribution of $N_\textrm{h}$ is as follows: 56.2% have 0 to 20 vortices, 18.9% have 20 to 40 vortices, 20.6% have 40 to 60 vortices, and 4.3% have 60 to 80 vortices. The RMSE value provides an estimate of the average difference between the machine and human counts per image. The confidence threshold is set to 0.5 for the CNN's vortex detection.

The performance of the machine for vortex detection is affected by the tuning parameters such as the upsampling size $n_\textrm{up}$ for imaging resolution, the cutoff density ρ_cut, and the noise strength σ_n in the formation of synthetic training sets. In this section, we explain how the machine performance is changed by the three parameters. We train various CNNs with different parameter settings and compare their RMSE values to determine the optimal parameter combination for accurate vortex detection.

We first investigate the effect of image resolution on machine performance by adjusting the upsampling size $n_\textrm{up}$ from one to five, while keeping $\sigma_\textrm{n} = 0.07$ and $\rho_\textrm{cut} = 0.01$ fixed. The results of the evaluation are shown in figure 3. As the size of the upsampling increases from $n_\textrm{up} = 1$ to 2, the machine performance improves significantly, reducing the RMSE value by a factor of about 3. Further increases in $n_\textrm{up}$ result in gradual improvements in detection performance, and for $n_\textrm{up}\gt3$, the RMSE value converges to approximately 3. The scatter plots of figure 3(a) compare the vortex number counts $N_\textrm{m}$ and $N_\textrm{h}$ for various $n_\textrm{up}$ to show the details of the machine performance. When $n_\textrm{up} = 1$, i.e. without upsampling, the machine exhibits poor performance, particularly for cases with a high vortex number, indicating that the image resolution is not sufficient to resolve vortices when clustered with high density (figure 3(b)). It is interesting to note that since the 4 × 4 grid cell is used as the unit window for single vortex detection in the CNN model, $n_\textrm{up} = 4$ corresponds to having a resolution of 1 pixel in the bare experimental image. This seems to explain the RMSE convergence for $n_\textrm{up} \unicode{x2A7E} 3$. The results for various $n_\textrm{up}$ suggest that we can utilize the CNN detection algorithm effectively for low-resolution experimental images by simple upsampling without interpolation. The required upsampling size may be smaller in experimental settings with better resolution.

Zoom In Zoom Out Reset image size

Next, we explore the effect of the parameters σ_n and ρ_cut on machine performance. We measure the RMSE value for a range of noise strength σ_n from 0.01 to 0.15, for different ρ_cut values and a fixed $n_\textrm{up} = 4$. The results shown in figure 4(a), demonstrate a U-shaped dependence on the noise strength σ_n, with an optimal noise strength of around 0.07, regardless of the density cutoff value used. At $\sigma_\textrm{n} = 0.07$ and $\rho_\textrm{cut} = 0.01$, the RMSE has a minimum value of approximately 3. Interestingly, a CNN trained with low-noise images tends to underestimate the number of vortices, although it often overestimates for images with a few vortices (figure 4(c)). This behavior appears counterintuitive because it is expected that, trained with low-noise images, the machine might interpret strong noises as vortices, leading to overcounting. When the machine is trained with high-noise images, it shows an underestimate of the vortex number for the entire experimental image test set, reflecting its low sensitivity to vortices, ignoring them as noises (figures 4(b) and (c)).

Zoom In Zoom Out Reset image size

The sample-area criterion, represented by the cutoff density ρ_cut in the synthetic images, also has a significant impact on the CNN performance. The density cutoff ρ_cut determines how close to the sample boundary is a vortex labeled in the synthetic image data. For noise strengths of 0.03 or higher, the CNN's performance improves as the cutoff density cutoff decreases, even to 0.01 (figure 4(a)). As shown in figure 4(d), the CNN trained for higher ρ_cut tends not to count density dips near the sample boundary as vortices. It is worth noting that Metz et al [45] found that a CNN works efficiently with $\rho_\textrm{cut} = 0.15$ to analyze simulated BEC images. The optimal value of the cutoff density may differ under various task conditions.

The peculiar dependence of the machine performance on noise strength σ_n observed in the previous section motivates further investigation to find a more effective implementation of signal noise in the synthetic training dataset. Our current approach to creating a synthetic dataset involves adding Gaussian noise with a uniform strength across the entire image. However, it should be noted that noise in actual experimental images is the result of a complex combination of multiple factors, and its strength can vary spatially over the sample and fluctuate for different realizations. The primary source is atom-shot noise, which is caused by the quantized nature of atoms and leads to fluctuations in atom density. In the simplest scenario of a non-interacting ideal gas, the counting of atoms, N, within a specific volume follows a Poisson distribution, resulting in a variance of atom number $(\Delta N)^2 \sim \bar{N}$ with $\bar{N}$ being the mean value of N. Even when interactions among particles are taken into account, atom-shot noise remains a significant source of noise for atom image data [54, 55]. Additionally, the inherent Poisson statistics associated with photon counting in the imaging process introduce another unavoidable source of noise as well as other technical noises.

To gain insight into the noise characteristics of our experimental images, we measure the noise strength and its spatial variations over the image. For each pixel, we calculate the difference δ, between its value and the average value of its four nearest and four next-nearest neighboring pixels. We then estimate the noise strength σ_n in a region of interest by taking the standard deviation of δ in the region. Figures 5(a)–(c) show the probability distributions of δ obtained from three different regions in the image; a high-atom-density center region, a low-atom-density sample boundary region, and a background region without atoms, respectively. The noise distributions are found to be well modeled by Gaussian profiles, and the measured σ_n values increase with increasing atom density, which is in agreement with the predictions from atom-shot and photon-shot noises. In particular, the noise strength in the high-density center region is $\sigma_\textrm{n} = 0.073$, close to the optimal noise strength of 0.07 for synthetic images, which explains the optimization. Additionally, figure 5(d) displays the histogram of the noise strength in the high-density region across all experimental images, showing that the noise strength fluctuates from one image to another, centered around 0.073.

Zoom In Zoom Out Reset image size

The proposition that the best training dataset is one collected from real task situations suggests that it is desirable to create a synthetic image dataset that includes all the noise characteristics of experimental images. However, it is difficult to incorporate every noise factor into synthetic images through a comprehensive theoretical analysis. As an alternative, we construct a training set with synthetic images that have different noise strengths while maintaining uniform Gaussian noise. We train the CNN using a set of 2000 training images with varying noise levels, ranging from 0.00 to 0.07, where $n_\textrm{up} = 4$ and $\rho_\textrm{cut} = 0.01$. The resulting RMSE value is 2.49, which is approximately 15% lower than the minimum RMSE obtained when using the optimally selected uniform noise strength $\sigma_\textrm{n} = 0.07$ (figure 6(a)). This improvement indicates that the non-uniform noise strength of the experimental images can be better addressed by the noise strength mixing method.

Zoom In Zoom Out Reset image size

Turbulent BECs contain a large number of vortices distributed in a disorganized manner (figure 6(b)) and they lead to experimental images with intricate patterns, such as linear chains of density dips created by multiple vortices and faint cracks at the sample boundaries that resemble vortices (figure 6(c)). This complexity makes it hard to accurately count the vortices, even for a human observer. Nevertheless, our vortex detection algorithm performs well, with an RMSE of ≈ 3, which is more than satisfactory for the purposes of this study.