arXiv CSopen access

Can machine learning for quantum-gas experiments be explainable?

2026 · arxiv_cs

arXiv CS · Papers · License: Open Access · 2026

machine learning, deep learning, neural networks

Can machine learning for quantum-gas experiments be explainable? I. B. Spielman 1, 2, 3, ∗ and Justyna P. Zwolak 1, 2, 4, † 1

National Institute of Standards and Technology, Gaithersburg, MD 20899, USA Department of Physics, University of Maryland, College Park, MD 20742, USA 3 Joint Quantum Institute, University of Maryland, College Park, MD 20742, USA 4 Joint Center for Quantum Information and Computer Science, University of Maryland, College Park, MD 20742, USA (Dated: May 19, 2026)

arXiv:2605.18689v1 [cond-mat.quant-gas] 18 May 2026

Virtually all aspects of many-body atomic physics are challenging: experiments are technically demanding, datasets have become enormous, and the memory and CPU requirements for classical simulation of generic quantum systems often scale exponentially with system size. Machine learning (ML) methods are already assisting in each of these areas and are poised to become transformative. Here, we focus on two specific applications of ML to cold-atom-based quantum simulators. These devices generally generate data in the form of images; we first showcase denoising of raw images and then identify solitonic waves in Bose-Einstein condensates. In both of these examples, we comment on the interplay between performance, model complexity, and interpretability.

INTRODUCTION

Today’s atomic quantum simulators and computers are experimentally elaborate [1, 2], produce increasingly large datasets, and are computationally bottlenecked: classical simulation of generic quantum systems scales exponentially with system size [3]. Machine learning (ML) methods already help navigate this combination of experimental complexity, data volume, and computational intractability, and are poised to have a transformational impact [4]. The preceding chapter by H. Schlömer and A. Bohrdt provides an up-to-date overview of the application of ML techniques to cold-atom physics, with both experimental and theoretical use cases. In this chapter, we take a complementary tack by explicitly showcasing the use of ML for two specific experimental tasks: the construction of “reference images” used in the first step of experimental data processing [5], and the identification of solitary waves [6, 7] in images of Bose–Einstein condensates (BECs). Together, these highlight the importance of high-quality training data [8] and the trade-offs in model architecture between performance and scientific transparency [9]. Almost all neutral atom quantum devices operate in a cyclic manner [2, 5]. Atoms are first laser cooled from tens to thousands of Kelvin and captured in a magnetooptical trap (MOT). A carefully scripted control sequence then performs further cooling and state preparation to yield the atomic ensemble under study. After a series of classical control steps and one or more measurements, the ensemble is discarded, and the process repeats. Thus, unlike solid-state quantum devices (and, to a lesser extent, ion-based systems), neutral-atom quantum systems are short-lived. Every stage of this process can be improved using ML, including hardware design, control se-

∗ [email protected] † [email protected]

quence discovery [10], measurement [11], and data processing [12, 13]. Quantum devices offer new capabilities for sensing weak signals, simulating quantum phenomena, and performing a range of computational tasks [14–16]; however, their performance is often limited by measurement. For example, in neutral-atom and ion arrays, the overall speed is limited by measurement times [2]; atomic clocks and magnetometers require sub-Poissonian atom-number detection for quantum-enhanced sensing [17]; and in coldatom quantum simulators, the identification of quantum correlations is often limited by technical noise [18]. Most neutral-atom sensors, simulators, and computers acquire data as images; in particular, quantum gas experiments often rely on bright-field techniques with coherent laser light, such as absorption or phase-contrast imaging [5], schematically illustrated in Fig. 1a. As seen in Fig. 1b, coherent-light imaging such as this presents unique challenges because each unwanted source of optical scattering, for example, dust on optical elements or etaloning, adds its own interference pattern to the detected image. The first section of this chapter contrasts traditional linear-algebra and ML techniques for removing these artifacts (Fig. 1c). Processing such as this is essential, as each stage of the overall analysis pipeline impacts the next: for example, improvements in initial data processing yield higher-quality datasets, which benefit all downstream analysis and enable further refinements of the initial stage. Even with high-quality denoised data, extracting a desired signal can be challenging. Our second section demonstrates this stage of data analysis through the feature identification task of locating excitations, called solitons, in BECs [6]. Solitons, robust shape-preserving solitary waves, are present in many physical systems and arise from an interplay between nonlinearity and dispersion. In experiments with ultracold gases, solitonic excitations provide a vivid window into nonlinear, nonequilibrium physics, and (excluding systems where the spin degree of freedom is relevant) are either “bright” or

2 (a) Schematic

Aperature Imaged dust

Quantum gas

Camera

Dust Probe laser Image plane

Object plane Fourier plane

Probe image

Dark image (×20)

200

100

200 x, pixels

CCD Counts

(b) Raw data Atoms image

y, pixels

FIG. 1. Bright-field imaging for cold atoms. (a) Schematic of the imaging system. A dust mote scatters light near the object plane, which is refocused near the image plane. (b) Raw image data showing a complex combination of interference patterns. The left panel shows the probe with atoms present, the center panel shows the probe alone, and the right panel shows the background image. The red curves denote the reconstruction window. (c) Optical depth. The left panel uses no reference recovery; the middle panel uses standard linear-algebra techniques; and the right panel uses a ResNet ML model.

ML Recovery 2.0

150

1.5 1.0

125

0.5 100

y, pixels

No recovery 175

0.0

75 150

200 x, pixels

“dark”. Bright solitons are stable for attractively interacting systems, where the tendency for the BEC to collapse to a single point is counterbalanced with the kinetic energy resulting from bending the macroscopic wave function. The reverse holds for dark solitons, which are stable in repulsively interacting systems. Owing to the comparative ease of preparing large repulsively interacting BECs, experiments with dark solitons are more common and are the focus of this chapter. Recent experimental developments have enabled the preparation of single solitons, or ensembles of solitons, with key properties—such as position and velocity—tunable by experimental design [19, 20]. In elongated BECs, dark solitons often appear as narrow density depletions in absorption or phase-contrast images, with physical information such as velocity and local background density reflected in the soliton depth and width. Large collections of images acquired across experimental parameters can therefore reveal soliton statistics, formation mechanisms, and interaction dynamics. In practice, however, much of this physics is locked inside high-volume image datasets. A first obstacle is surprisingly basic: before we can study soliton physics quantitatively, we must reliably determine which images contain solitonic excitations, the number of excitations,

150

200 x, pixels

150

200 x, pixels

−0.5

and their positions. This “localization-and-classification” step quickly becomes the dominant bottleneck when experiments produce thousands to tens of thousands of images [21, 22]. The core difficulty in localizing and classifying solitonic excitations is the lack of a fixed, rigid visual template for “a soliton in a BEC image.” The observed depth and width of a dark soliton depend on its velocity and the imaging resolution, while shot-to-shot atom-number fluctuations, thermal effects, and imaging artifacts can create soliton-like dips or obscure genuine excitations. The situation becomes even more challenging when multiple excitations coexist: solitons can be close together, merge, decay, or move between frames, complicating the implementation of simple rule-based detectors. These realities make naive automated heuristics brittle. One effective response is to use convolutional neural networks (CNNs) trained directly on images. CNNs can learn robust visual features for classification (e.g., “no soliton” versus “soliton present”) and can be extended to localization or object detection when the task demands soliton positions. CNN-based approaches have been shown to dramatically speed up dark-soliton detection and enable high-throughput analysis that would be impractical by eye [21].

3 CNNs introduce two practical challenges for scientific use. First, as supervised models, they require substantial labeled data for training and validation; in the absence of reliable simulations, labels often come from manual annotation of experimental image collections, which is slow, labor-intensive, and can be inconsistent across annotators or across experimental conditions [23]. Second, even when accuracy is high, CNN decision rules are typically opaque: the model may memorize subtle correlates (e.g., imaging artifacts or background structure) that are not part of the underlying soliton physics, complicating verification and failure-mode analysis. Common explanation tools for CNNs, such as saliency maps [24] or Gradientweighted Class Activation Mapping (Grad-CAM) [25], provide heuristic visualizations rather than a transparent model of how each physical feature influences the decision. This is a critical issue when the goal is not just a label, but a trustworthy scientific inference workflow that includes an interpretable link between the label and the underlying physically meaningful cues [22]. In this spirit, we consider alternatives—with CNNlevel performance—that yield human-readable explanations of the evidence used for each decision. Specifically, we compare a high-accuracy image-native CNN classifier with an explainable boosting machine (EBM) model that operates on physically motivated, low-dimensional image features and produces transparent, additive explanations [26]. EBMs are a family of so-called glass-box models: their predictions decompose exactly into a sum of learned feature-shape functions (and, optionally, a small number of pairwise interaction terms), enabling global visualizations and per-image attributions without posthoc explanation methods [27]. In our context, the fully intelligible EBMs can not only match but actually surpass CNN classifiers’ performance, addressing the “what did the model decide?” as well as the “what evidence did it use, and does that evidence align with the physics?” questions. With this backdrop, we can now articulate the most essential message of this chapter: ML works best starting with large quantities of high-quality data, a well-defined figure of merit, and the simplest viable model. Data is the bedrock of ML: a good dataset (labeled or unlabeled) contains a dense and complete set of “examples” so that the right trained model can learn to interpolate between them [8]. This underscores the importance of understanding the dataset’s content: what are typical signals? What is the scale of noise or uncertainty? What features of the data are of interest? How often do these features appear? The answers to this sort of question will guide the assembly and curation of the dataset. ML is not a free lunch: good performance without spending considerable effort assembling the dataset is rare. The “complexity” of an ML model depends on a range of interrelated considerations, including the size of the data vector traversing the model, the sophistication of the model architecture, the hyper-parameters required to describe the model (such as the number of layers in

a neural network), as well as the overall number of parameters in the model. However, in practice, the most common failure mode in applying ML to experimental physics is not an exotic choice of architecture, but an illposed learning objective: before collecting data or training a model, one must specify what is to be learned and how success will be measured. This means translating a scientific goal into a concrete prediction task (classification, regression, denoising, segmentation), defining the target quantity and its domain (e.g., a “reference image” vs. an artifact-free density or excitation count vs. excitation locations vs. excitation type), and choosing labels or surrogate targets that are consistent with the physics and the measurement process. Equally important is deciding what variability the model should be invariant to (shot-to-shot atom number, global translation, fringe patterns) and what variability it should remain sensitive to (true excitations, physically meaningful contrast changes), since these choices determine whether the model learns the intended signal or a convenient correlate. Finally, a well-defined figure of merit turns “it seems to work” into a quantitative statement and exposes failure modes early. As we will see, these choices influence training efficiency, performance, and interpretability. II.

EXPERIMENTAL IMAGING OF COLD-ATOM SYSTEMS

Our examples employ bright-field techniques such as absorption or phase-contrast imaging. As depicted in Fig. 1a, an illuminating probe laser beam is first attenuated or phase-shifted by the atomic sample [5, 28]; then, the resulting optical field traverses the imaging system and the final 2D intensity I(x, y) is imaged as a function of position (x, y). Relevant atomic properties, such as the spatial density, are encoded in the fractional transmission ≃ I(x, y)/Ip (x, y), and therefore require knowledge of the probe intensity Ip (x, y) in the absence of atoms. In practice, reference images are acquired shortly before or after the data image; during this interval, coherent-imaging artifacts often change slightly, leading to imperfect division. For example, in Fig. 1a, the probe laser is shown to travel into a vacuum system (square) where it encounters a dust mote, leading to an outgoing scattered optical field. This pattern of light is faithfully imaged to the “imaged dust” just prior to the camera, where the interference pattern between the scattered light and the probe, as attenuated by the atomic ensemble, is measured. In realistic imaging systems, several such sources can be simultaneously present, leading to a highly structured probe laser, as in Fig. 1b (left). Mitigating these artifacts is the first focus of this chapter, in Sec. III. In what follows, we exclusively employ resonant absorption imaging to obtain the 2D column density of atoms Z n(x, y) = dz ρ(x, y, z), (1)

4 the 3D density ρ(x, y, z) integrated along the imaging line of sight. In the simplest case of low intensity, each atom contributes equally to the attenuation of the probe beam. This leads to Beer’s law exponential behavior characterized by the optical depth I(x, y) − Ibg (x, y) , (2) OD(x, y) = − log Ip (x, y) − Ibg (x, y) giving n(x, y) =

OD(x, y) , σ0

(3)

where σ0 = 3λ2 /(2π) is the resonant absorption cross section for a probe laser of wavelength λ, and Ibg (x, y) is an additional image taken with the probe laser off to account for any common-mode background signal. For brevity, this discussion omits several considerations that are irrelevant for reference reconstruction, but are important for accurate determination of n(x, y); see Refs. [5] and [28] for details. It is common to release the atomic ensemble from its confining potential for a brief timeof-flight (TOF) prior to imaging. In many cases, this period of nominally ballistic evolution enables direct access to the momentum distribution; however, for the data presented in Sec. IV, TOF expansion is dominated by interaction effects. In most experiments, dark solitons are too small to be resolved in-situ; they are made visible by TOF expansion. Figure 1b shows an example with atoms (left), probeonly (center), and dark (right) images; these data were taken in situ (i.e., with no TOF). These data were selected to illustrate the consequences of a highly contaminated probe laser: numerous overlapping interference structures are visible. When the scattering sources leading to this pattern spatially drift by even a fraction of the optical wavelength between these exposures (for example, due to acoustic vibrations of optical elements), coherentimaging artifacts change noticeably, leading to imperfect division. The red squircle (a function that interpolates between a square and a circle) encloses the expected location of the absorption signal, which is barely visible to the eye. Simply computing the OD using Eq. (2) yields the image in Fig. 1c (left) with artifacts comparable in magnitude to the underlying signal. The artifacts are greatly reduced using both conventional linear algebra (middle) and ML methods (right), finally revealing the pattern of four separate groups of atoms. These data were taken from Ref. [29], which followed the dynamics of these groups as they moved in a larger background BEC (not visible). III.

REFERENCE RECONSTRUCTION

With this basic understanding of bright-field imaging’s how and why as it applies to cold-atom experiments, we return to the reference recovery task in more detail. At

a high level, the task is to take data as in Fig. 1b and predict the intensity Ip (x, y) that actually illuminated the atomic ensemble; in practice, this amounts to using the intensity exterior to the red squircle (where no atoms are present) to predict the interior intensity. This can be accomplished by collecting a large “training set” of reference images (no atoms present) that reflect the full variability present in the probe beam. For relatively clean beams, this might require 10 to 100 images; the depicted probe beam is a worst-case scenario, requiring a training set of ≈ 104 images. Using these data, both linear-algebra (LA) and ML techniques can be employed to construct an optimal reference image for each image containing atoms [30–32]. LA methods are established tools that greatly reduce coherent-imaging artifacts, but because the observed intensity is a quadratic function of the underlying complex field, intensity-domain summation can require many linear combinations to eliminate residual artifacts. When framed in the ML context, these are very simple, highparameter-count models that lead to a highly compressed intermediate vector representation, i.e., a latent space. In this context, such models are highly interpretable, perform well, and are trivial to “train.” Here, we contrast LA with more traditional ML techniques and explore the trade-offs between interpretability, performance, and development effort. LA methods operate in a highly intuitive way. They construct a collection of independent templates for the interior and exterior regions; then, for an image with atoms, they use the exterior-region templates to determine the appropriate linear combination of interiorregion templates. In this way, LA templates can have the physical interpretation patterns that tend to change together, ideally from individual scattering sources. This physical interpretation breaks down for non-linear contributions resulting from interference between multiple sources. A key secondary benefit of this approach is noise filtering: white noise from all sources (photon/photoelectron shot noise, read noise, ...) will generally have little in common with the templates and will be removed from the reconstructed reference. By contrast, published ML-based tools for reference recovery use U-Net architectures where down-sampling CNN layers are followed by up-sampling with ConvTranspose or PixelShuffle layers, along with bypass connections [33–35]. These models only slightly outperform LA approaches, are opaque black-box tools, and offer no in-principle reduction of white detection noise outside the windowed region. This section begins by introducing conventional LAbased reference-recovery methods [30–32], which serve as a benchmark for our ML implementation. We then demonstrate that a complex, ResNet-based autoencoder yields only slightly improved performance at the expense of interpretability, efficiency, and complexity.

5 A.

Linear algebra formalism

Since the following analysis does not rely on the images’ two-dimensional (2D) spatial structure, we work with flattened P -dimensional vectors, where P equals the number of pixels. Vectors and linear operators are denoted in boldface, such as v and O, with components vp and Opq . The conjugate transpose is denoted ∗ by [O† ]pq ≡ Oqp (index reversal with conjugation); and the identity operator is I. The basic approach uses a “training” set T = {tα } of independently measured reference vectors, collected as the columns of T, a P × |T | matrix with elements Tpα . In practice, T should be sufficiently large that its variability is representative of the underlying distribution. We adopt the convention that Latin indices (p, q, r, . . . ) denote image-space (pixel) coordinates, while Greek indices (α, β, γ, . . . ) label set elements (training samples, basis vectors, . . . ). We then use these vectors to optimally recover each data vector d in regions unaffected by the atomic system and extrapolate to the remainder, optimized with respect to a windowed least-squares criterion. This process requires that the signal of interest be concentrated within an “interior” region, isolated by a window W, a square matrix with eigenvalues ∈ [0, 1] (typically diagonal in the pixel basis, but this need not be the case); the remaining “exterior” region is correspondingly isolated by a mask M with the same spectral bounds (often M = I − W, but this is not required). For any mask or window O, we define a metric GO ≡ O† O (positivesemidefinite Hermitian), which induces the inner product ⟨u, v⟩G = u† Gv and norm ||v||2G = ⟨v, v⟩G . In particular, we use the interior, exterior, and Euclidean metrics GW , GM , and GE = I, respectively; furthermore, unannotated inner products and norms (e.g., ⟨u, v⟩) are taken to be Euclidean. This framing is particularly useful when uncertainties are known in the form of an inverse covariance matrix Σ−1 , which can be expressed as Σ−1 = R† R in terms of the Cholesky factor R. For precision-scaled residuals, the following arguments follow through with each metric changed to GO → GOR = R† O† OR. In practice, when uncertainties are known, we reconstruct precision-scaled data Rd, and use precision-scaled training data. This removes the need for further discussion of uncertainties and sidesteps the nuance that, while W and M are fixed, the per-observation uncertainties in general differ. For every data vector d, our task is to optimally reconstruct a reference r using a basis of linearly independent vectors B = {bα }, collected as the columns of B, a P×|B| matrix with elements Bpα . Since the basis vectors will be generated from the training set T , we require |B| ≤ |T | to avoid underdetermination. Conventional linear techniques employ a matrix product r = B a(d)

(4)

with amplitudes given by a linear map a(d) = A† d, de-

fined by the P × |B| matrix A. The elements of B can further be taken to be orthonormal with respect to a chosen metric, i.e., B† GB = I, since any right-factor used to orthonormalize B can be absorbed into A.

Basis sets

Both the basis B and amplitude map A must be determined from the training data T to best address the question: which set of vectors best describes the interior region (where signal is present) using the amplitudes A† d determined from the exterior? By contrast, the two most common approaches simply reconstruct complete reference images. The first of these approaches directly uses all |T | training vectors as the basis, i.e., setting B = T. By construction, these fully span the training dataset and are operationally linearly independent because |T | ≲ 103 is, in practice, much smaller than the dimension P ≈ 106 of the vectors derived from typical 1024×1024 images.

PCA basis set

The second common approach uses principal component analysis (PCA) to generate an orthonormal basis BPCA with a number of elements |BPCA | ≤ |T | sufficient to capture the structure present in the training set while suppressing contributions from stochastic noise. In this simple example, we focus on the case of no mask or window GM,W = I, giving the Euclidean orthonormality condition B† B = I. As shown below, this approach often yields substantial dimensionality reduction, with |BPCA | ≃ 16 even when |T | ≃ 1000. In this context, we seek an orthonormal basis of dimension |BPCA | and transfer matrix A that together most effectively describe the training set, with residual vectors δtα ≡ tα − BA† tα . Collecting these residuals yields the P ×|T | residual matrix δT = T − BA† T. Given these definitions, the optimal basis minimizes a loss function given by the total variance of the residual vectors L=

X α

||δtα ||2 = ||δT||2 = Tr δT† δT ,

(5)

where ||O|| denotes the Frobenius norm of the matrix O. Expanding and using the cyclic property of the trace yields the quadratic loss function L = Tr TT† I − BA† − AB† + AA† ,

(6)

where the first term is the total variance of the training set, and the remainder quantifies the variance explained by the PCA basis. The required orthonormality constraint adds a Lagrange term Tr Λ B† B − I to L, with the Lagrange multiplier matrix Λ.

6 Taking matrix derivatives with respect to A† and B† leads to the coupled equations 0 = TT† A − TT† B

and

0 = BΛ − TT† A

(7)

where the orthonormality constraint implies Hermitian Lagrange multipliers ΛW , giving the required real valued eigenvalues. In direct analog with the PCA case, we solve this problem by left-multiplying both expressions by T† , obtaining the coupled equations

the first is simply solved by setting A = B. In the second, we identify Λ as a diagonal matrix whose entries are the largest |BPCA | eigenvalues of TT† (a positive semidefinite Hermitian matrix with non-negative eigenvalues), and taking the corresponding eigenvectors as the columns of B. When every eigenvalue is included, L → 0, indicating that every training vector can be perfectly recovered. In practice, this indicates overfitting because each probe image has a nonnegligible contribution from stochastic noise sources such as photon shot noise. This formulation obscures the fact that, although TT† is a P ×P matrix, it has at most |T | nonzero eigenvalues. We address this by left-multiplying by T† to obtain a new eigenvalue problem

in terms of the |T | × |T | Hermitian matrices CO = T† GO T, the |T | × |B| basis vectors −1/2 B̄O = CO T† GO BO , orthonormality constraint B̄†W B̄W = I, and with O ∈ {W, M}. Again, following the PCA case, the first equation is −1/2 1/2 quickly solved with ĀM = CM CW B̄W leaving behind the eigenvalue equation

0 = B̄Λ − CB̄,

0 = B̄W Λ − CW B̄W .

(8) †

with the |T | × |T | Hermitian matrix C = T T, a |T |×|BPCA | matrix of basis vectors B̄ = C−1/2 T† B (with C−1/2 defined on the support of C), and the orthonormality condition B̄† B̄ = I. This is conveniently solved by the singular value decomposition (SVD) of T = UΣV† , giving the diagonal singular-value matrix Σ → Λ1/2 , the right singular vectors V → B̄, and the left singular vectors U. Thus, the final basis vectors are obtained via the inverse transformation B = TC−1/2 B̄ = U. D.

Optimal basis sets

We now go beyond the standard PCA case by asking: which exterior basis provides the most information about the interior region, and which basis best describes the interior? Relative to the PCA calculation in the preceding section, addressing these questions requires two changes. (i) Including the mask-metric GM to the transfer matrix, i.e., giving coefficients a = A† GM d in Eq. (4). (ii) Using windowed norms ||δtα ||GW in Eq. (5). The PCA solution employed the orthonormality condition B† B = I, which rendered the loss function quadratic in B, enabling a simple linear-algebra solution. Here the use of the windowed norm suggests the use of an interior basis matrix BW with orthonormality constraints B†W GW BW = I. In the masked-extrapolation setting, we minimize the interior reconstruction loss X L= ||δtα ||2GW = Tr δT† GW δT , (9) α

with residuals δT = T − BW A† GM T. Taking matrix derivatives of L with respect to A† and B†W yields the coupled equations 0 = GM TT† GW BW − GM TT† GM A 0 = GW BW ΛW − GW TT† GM A

1/2

0 = CW B̄W − CM ĀM 1/2 1/2 0 = B̄W Λ − CW CM ĀM

(10) (11)

(12)

This is solved by introducing the interior and exterior † SVD expansions of OT = UO ΣO VO , with left singular vectors UO , right singular vectors VO and singular values ΣO . All together this gives solutions: Λ = Σ2W , † B̄W = VW , and ĀM = VM Σ−1 M VM VW ΣW . Or in the initial representation WBW = UW † MA = UM Σ−1 M VM VW ΣW .

(13) (14)

If dimension reduction is desired it is implemented by truncating ΣW and the columns of UW and VW . E.

Linear algebra implementation

We now move from formal derivation to the implementation, performance and interpretability of the two linear methods—“standard PCA” and “optimal linear”— before comparing them with more complex deep-learning models. At the level of numerical linear algebra, both of these algorithms consist of just a few calls to the BLAS or LAPACK libraries. While it is true that these libraries do the numerical heavy lifting, considerable effort is required to properly utilize the dataset comprising ≳ 4 × 103 resonant absorption images acquired over about a week of experimental operation. We therefore begin by exploring the limitations of a simple-minded deployment, then introduce the decomposition of datasets into training, validation, and test sets, and finally introduce k-fold crossvalidation. The first results of both reference recovery implementations are shown in Fig. 2a; both standard PCA (blue) and optimal linear (red) each yield a spectrum that, aside from the j = 0 component, approximately follows a power law with exponent −2 (dotted black) before reaching a noise floor[60]. Typically, components in the powerlaw regime contain useful information, while those in the

Eigenvalue

(a) 100

(b) j=0

j=1

j=2

j=3

j=4

j=5

j = 73

j = 74

10−2

10−4

10−6

100

(c)

101 Component index

2.97

Loss

102

2.95

2.93

101

FIG. 2. Linear algebra-based reference reconstruction. (a) Normalized standard PCA eigenvalues (blue), the optimal linear spectrum Σ2W (red), and a power-law with exponent −2. (b) The first 6 and last 2 standard PCA components, with windowed and masked regions delineated by the red squircle. (c) Reconstruction loss from 8-fold cross-validation for standard PCA (blue) and the optimal linear (red) methods. Training and validation losses are denoted by solid and dashed curves, respectively. The inset expands the nearsaturation regime and adds 1-sigma error bands. Both validation curves exhibit local minima marked by circles.

100

30 50 # retained

101 Components retained

plateau regime are predominantly noise. This overall impression is confirmed by the first and last few principal components for the standard PCA method shown in Fig. 2b, where the j = 0 component is essentially the average probe; j = 1 shows a frequently appearing diffraction ring along with additional features; and j = 2 includes a phase-shifted diffraction ring with even more artifacts. By contrast, the final components are dominated by white noise. This interpretability is possible because the principal components are themselves simple images, i.e., templates, whose importance is given by the spectrum in Fig. 2a. Each eigenvalue can be interpreted as the reduction in the training loss [Eqs. (5) or (9)] resulting from that component (this is usually called the explained variance). Finally, Fig. 2c plots the training loss resulting from the first j components for both the standard PCA (blue) and optimal linear (red) methods[61] In both cases, the training loss decreases monotonically as more components are included, suggesting that the best strategy is to retain every component. This is not the case. In ML language, these examples can be viewed as purely linear models; for example, retaining 75 components yields a “model” with ≈ 12 × 106 parameters[62]. Models of this size are susceptible to overfitting and memorization; in this case, each element of the training dataset contains a non-generalizable component consisting of photon shot noise, which will be learned when a sufficiently large number of components are retained. Data science resolves these types of issues by first randomly splitting the overall dataset into k ≈ 10 subsets, designating one subset for validation, and leaving the remainder for training. When repeated for all k possible

validation sets, this is called k-fold cross-validation. The dashed curves in Fig. 2c (only visible in the inset) show 8-fold cross-validation losses for both the standard (blue) and optimal (red) methods along with error bands reflecting the variation across the folds. Unlike the training loss, which monotonically decreases, the validation loss in the inset first saturates, and then increases: this is typical over-training behavior. This analysis identifies the optimal number of components to retain, 37 for the standard PCA method and 24 for the optimal linear method. In this case, the added sophistication of the optimal-linear approach pays off, not in improved ultimate performance, but rather a reduction in overall parameter count. These linear models serve as baselines for the non-linear models introduced next. In a more complete ML workflow, one would first set aside a randomly selected subset, ≈ 10 % of the full dataset, as a test set that is not used during model design, training, or validation. For demonstration purposes in this chapter, we did not set aside a test set.

ML model modalities

From an ML viewpoint, both of our linear methods are minimal encoder-decoder models with two dense layers (one input and one output), no activation functions, and with equal input and output dimensions; in other words, they are purely linear autoencoders, schematically shown in Fig. 3a. In abstract terms, these models linearly map a P -dimensional input vector d (flattened image) to a lower-dimensional latent vector a of dimension |B|, and then employ a second linear transform to produce a re-

8 FIG. 3. ML model modalities for reference recovery. (a) A simple encoderdecoder. When no nonlinear activation functions are present, this reduces to the linear models introduced in Secs. III C and III D. (b) A U-Net with shortcut connections. In both panels, rectangles denote data vectors and trapezoids denote layers, oriented to indicate the typical direction of dimension reduction or expansion.

(a) Encoder-decoder

(b) U-Net autoencoder

Input

Decoder

Encoder

Decoder

Encoder

Decoder

N repeats

Encoder

Decoder

Latent Decoder Output

Latent

constructed P -dimensional output vector [e.g., Eq. (4)]. The overall encoder-decoder architecture is agnostic to the specific contents of the “encoder” and “decoder blocks.” For example, using the notation introduced in Sec. III A, one such block could be a standard dense layer ! X ap = ϕ([Ad + b]p ) = f Apq dq + bp (15) q

which combines a linear transformation A, an additive bias b, and a pointwise nonlinear activation function ϕ applied to each output (see Sect. IV E for common choices). A single such layer is called a perceptron, and stacking two or more perceptrons gives the multilayer perceptron (MLP) architecture [39], one of the canonical deep-neural-network architectures. The first stage of these models flattens 2D images into simple vectors, discarding implicit spatial structure and forcing the model to relearn any relevant spatial relations. This is inefficient, often requiring a large number of parameters and long training times when direct methods such as SVD are unavailable. CNNs overcome both limitations by employing small convolutional kernels that operate directly in image space. This replaces the “encoder” stage of an encoder-decoder with one or more convolutional layers. Likewise, the “decoder” stage can be replaced by convolutional-like layers such as ConvTranspose or PixelShuffle. ConvTranspose is an approximate deconvolution layer, and PixelShuffle operates by first using a CNN to grow its 1D depth-wise channels and then reshaping them into a 2D image. Published ML-based solutions to the reference recovery task typically employ the more sophisticated U-Net architecture shown in Fig. 3b. Originally developed for image segmentation, the U-Net extends the autoencoder pattern with shortcut connections and often gives highquality results even with modest-sized training sets [40]. While the U-Net pattern can, in principle, use generic encoder and decoder blocks, it typically uses CNNs for encoding and ConvTranspose or PixelShuffle layers for decoding.

ML implementation

With this backdrop, we finally turn to an ML implementation of the reference recovery problem. For brevity, this demonstration augments the encoder-decoder pattern [63] with a latent space ResNet [42] to deliver a modest performance improvement. This comes at the expense of drastically increased implementation overhead, including model design, regularization, data augmentation, hyperparameter selection, and training. Model design Several model architectures were explored before converging on the ResNet encoder-decoder shown in Fig. 4. The ResNet encoder-decoder model organizes individual encoders and decoders (trapezoids) into high-level functional blocks with two distinct low-dimensional data pathways. The “latent” path uses purely linear encoders and decoders (in the encoder-embedder and deembedderdecoder stages); in isolation, it recovers the performance of the linear models. The “residual” path is transformed to and from image space by single-layer perceptrons (in the encoder-embedder and deembedder-decoder stages). ResNet blocks then exchange information between these two data paths, where the latent path carries the already highly predictive linear reference-recovery information, while the residual path corrects for nonlinear artifacts left behind by the linear model. At the architecture’s level, this model introduces just three hyperparameters: the number of ResNet blocks NUM RESNET, and the dimensions of the data paths DIM LATENT and DIM RESIDUAL.

Regularization and data-augmentation Larger or more complicated models are often prone to over-training and generally require large datasets. In this context, regularization methods combat overtraining and accelerate training, while data augmentation expands the effective size of the dataset (which can also reduce overtraining).

9 Input Encoder Latent Annotator

Encoder Embedder

Embedder Residual

Residual Encoder

ResNet block

Residual Decoder Extractor Latent Decoder Latent Encoder

Decoder

Deembedder Decoder

Deembedder

Output

FIG. 4. ResNet encoder-decoder. Data-flow is indicated by lines, while dimension-changing layers are denoted by trapezoids. Functional groups are denoted by larger rectangles.

Our example ResNet encoder-decoder model-level implements regularization with both a Dropout layer (prior to the encoder-embedder block) and BatchNorm layers (following the residual encoder and residual decoder layers). A Dropout layer randomly zeros elements of the data vector between layers (and rescales them to maintain the overall signal level). BatchNorm layers rescale data vectors so that each component has a unit standard deviation across the training dataset. The optimal location of these layers in the model is determined empirically. Data augmentation can be quite sophisticated, for example, generating data with auxiliary models. For simplicity, we augment by including small spatial displacements of the images prior to flattening. This introduces two more hyperparameters: the dropout fraction DROPOUT FRAC, and the maximum augmentation shift MAX SHIFT.

Training In contrast to the linear models, the addition of nonlinearities and hidden layers makes the loss landscape complex, with many local minima. This is a generic problem

in deep learning, and has spawned many gradient-based optimizers, with well-known examples including Adam, AdamW, and Lion. These optimizers themselves contain hyperparameters, such as the learning rate, and, in the case of AdamW, additional regularization parameters. AdamW was found to outperform other options for this task; the weight decay present in AdamW controls the growth of large eigenvalues in the linear Encoder/Decoder stages of our ResNet model. In addition, training follows a schedule in which one or more parameters are adjusted throughout training. For example, here we begin with a OneCycle stage for 1000 epochs that quickly increases and decreases ETA, followed by a linear decrease for 9000 epochs, and concluding with a power-law with exponent −1 for the remainder of training. All together, this introduces several more hyperparameters: the overall learning rate ETA and the weight decay LAMBDA, in addition to the many parameters describing the training schedule. All hyperparameters were found empirically.

Hyperparameter selection All hyperparameters must be individually tuned for optimal training and final model performance. This requires a second stage of optimization wrapped around the minimization problem solved by training. In practice, one finds coarse values for these using individual training runs, optimizing on the validation loss, and then switches to longer k-fold cross-validation runs to identify the final values. For example, DIM LATENT and DIM RESIDUAL were initially set to the ideal number of components retained in the linear models, and then optimal performance was found by searching away from this point. Searching across the full hyperparameter space can be timeconsuming; in this case, it requires testing 100+ hyperparameter combinations to achieve acceptable performance and reasonable training speed. The 8-fold cross-validation training curves for the final ResNet encoder-decoder model are shown in Fig. 5, with training and validation data indicated by the solid and dashed curves, respectively. In both cases, the loss drops by orders of magnitude in the first ≈ 100 epochs before slowing down; by epoch 1000, training has drastically slowed, entering into a fine-tuning phase. As seen in the inset, the validation loss continues to decrease after 105 epochs and is slightly lower than that of the linear models (markers); the final training loss is well below the validation loss, indicating some degree of memorization. The ResNet uses the same reconstruction loss [Eq. (9)], preprocessing, image resolution, and cross-validation folds as the linear benchmarks, so differences in Figs. 2 and 5 are attributable to the model rather than the data pipeline. Careful inspection of the images in Fig. 1c confirms

10 flow. In this language, a solitonic excitation is not merely a dip or a bump in density; it is a coupled structure in both density and phase, stabilized by a balance between dispersion (which tends to broaden a wave packet) and nonlinear interactions (which can self-steepen or selffocus the wave).

Loss

103

Loss

2.97

2.95

102 2.93

103

101 101

102

103 Epoch

104 Epoch

104

105

FIG. 5. ResNet training with 8-fold cross-validation. Training and validation data are denoted by solid and dashed curves, respectively, and include 1-sigma error bands. The inset zooms in on the fine-tuning stage of training, and the markers show the results for the linear models in Fig. 2 for conventional PCA (blue) and optimal linear (red).

that the background is slightly less noisy for the ML case. Still, this improvement came at the expense of vastly increased effort (days rather than weeks), both at the human and computational levels. Furthermore, the limited explainability of the linear models has been sacrificed. The main conclusion is that ML methods can offer improvements upon standard numerical methods, but the gains can be modest, and the final value proposition should be considered. For example, it is much more straightforward to integrate LA-based reference recovery into daily laboratory operations, where real-time data analysis informs decision-making. Further, experiments focusing on the overall density distribution, such as the Gaussian thermal distributions in Fig. 1c, would seldom benefit from the often negligible improvements offered by ML models. However, experiments studying the pixel-bypixel fluctuations of the atom signal away from its mean (i.e., atom shot noise) are an ideal use case for ML tools, as they often operate in regimes where the imaging noise scale exceeds the signal of interest.

IV.

LEARNING WITH UNDERSTANDING: CASE OF SOLITONIC EXCITATIONS

Solitons are robust, localized wave structures that, despite waves’ typical tendency to spread, preserve their shape as they move through their host medium. They appear in many nonlinear media, such as water [43, 44], optical fibers [45], and plasmas [46]. Bose-Einstein condensates offer an unusually simple, clean, and controllable environment in which to study solitonic excitations. This is because over a wide range of conditions, the many-body system can be described by a single macroscopic wavefunction with a well-defined amplitude and phase [47]. The amplitude encodes the atomic density, while spatial variations of the phase describe superfluid

Repulsively interacting 1D (or quasi-1D) BECs support dark (or gray) solitons, of which the static black soliton is a special case. They can be thought of as moving “phase defects” that produce a localized depletion in the atomic density [48]. A black soliton is a static solution that fully depletes the condensate density and is accompanied by a full π phase jump. Gray solitons are moving solutions of increased width that only partly deplete the density, and have a phase jump that drops to zero (or increases to 2π) as their velocity increases to the condensate speed of sound [49]. Seminal experiments demonstrated that such solitons can be created in a controlled way—for example, by imprinting a spatial phase pattern onto a BEC (“phase imprinting” or “phase engineering”) and then watching the resulting excitation propagate [6]. As noted in the introduction, attractively interacting BECs have bright, rather than dark, soliton solutions. In elongated (cigar-shaped) condensates—where the axial dynamics are often much slower than the tightly confined radial dynamics—these excitations often appear as narrow, stripe-like density depletions in TOF imaging (see Sec. II). As we mentioned in the Introduction, several practical challenges complicate the analysis of such images: the observed contrast and apparent width depend on the soliton’s velocity and on the experimental imaging resolution; thermally excited atoms can obscure the signal, while interference fringes add distortions; and there is no single reliable visual template for identifying solitonic excitations. Moreover, the acquired 2D images are typically projections of a 3D cloud, meaning that a localized 3D structure can appear as a deceptively simple 2D stripe. As a result, the same physical excitation can look different across experimental conditions, and distinct excitations can look superficially similar in a single-shot image. For high-throughput analysis, it is often useful to adopt labels that reflect what is reliably distinguishable in the imaging modality at hand. In TOF images of elongated BECs, one can define coarse classes based on the presence or absence of well-defined features as “no excitation,” “single excitation,” and “other excitations” (e.g., multiple stripes, low-contrast events, or otherwise ambiguous structures) [21, 23]. Within the single-excitation category, using additional data-processing approaches, such as modeling and physics-informed fitting, allows for introducing finer distinctions between longitudinal solitons[64] and (transverse) solitonic vortices [22]. In the following section, we briefly overview the most common types of observed excitations in elongated BECs.

11 x

1 ms

2 ms

5 ms

7 ms

10 ms

FIG. 6. Images showing evolution of the BEC (measured after TOF) from 1 ms (panel A) to 10 ms (panel E) after imprinting a phase step of about 1.5π on the top half of the condensate. A positive density disturbance is moving rapidly in the +x direction, while a dark solitonic excitation (most likely a vortex ring) is moving significantly slower in the −x direction (opposite to the direction of the applied force). Reproduced with permission from Ref. [6].

A short guide to soliton types in elongated BEC experiments

Since many experiments with BECs are conducted in highly elongated traps, it is natural to first model the excitation as effectively one-dimensional motion along the long axis. However, real condensates are only quasi 1D rather than strictly 1D, and the residual threedimensional structure can matter qualitatively. Planar dark solitons in three dimensions can be dynamically unstable and can evolve into other long-lived solitarywave structures [51]. A useful way to organize the landscape is therefore to distinguish solitons that are well approximated as quasi-1D kinks in density from 3D solitary waves that retain soliton-like density signatures but carry vorticity and more complex geometry, such as solitonic vortices, vortex rings, and more complex structures [51, 52]. Dark and gray kink solitons—These are canonical “soliton” excitations: a localized depletion accompanied by a change in the phase that becomes a step for a static black soliton In images of elongated BECs, they often appear as a stripe-like depletion roughly transverse to the long axis, see Fig. 6. Slow solitons tend to appear deeper and more localized, while faster solitons are typically shallower and harder to distinguish from background structure. In our experimental pipeline, the term kink soliton is used to emphasize this phase-change character, even when the imaging signature is simply a density notch [48, 53]. More broadly, once three-dimensional effects become relevant, it is best to think in terms of a family of related solitary-wave excitations rather than a single “soliton” archetype [54]. Solitonic vortices—A solitonic vortex is a vortex line embedded in an elongated BEC whose density depletion can resemble a soliton sheet, but whose phase twists around a core due to the phase gradient. Experimentally, solitonic vortices can leave distinctive signatures after expansion, such as a twisted planar density depletion, and interferometric measurements can reveal dislocations associated with circulation [52]. For image-based classification, the important lesson is that a stripe-like depletion does not uniquely identify a kink soliton: vortical excitations can masquerade as solitons in projection, see Fig. 7.

Vortex rings and other 3D solitary waves—Vortex rings are closed-loop vortex lines that can also appear as localized depletions, sometimes evolving from unstable planar solitons [51]. Depending on geometry and imaging axis, a vortex ring may appear as a pair of depletions, a partially filled notch, or a more complex structure. To summarize, solitonic excitations are rich but visually heterogeneous: their appearance varies with velocity, dimensionality, and imaging artifacts, and different excitations can project to similar stripe-like patterns. This is exactly the setting where one might want to maintain the balance between high-throughput classification and localization—CNNs excel here—and model decisions that can be audited against physics—which is where interpretable models become attractive. In the following sections, we discuss these two settings in more depth. We also introduce the simplified dark solitons dataset we use for all demonstrations. B.

The dark solitons in BECs dataset 2.0

The dataset used for experiments discussed later in this chapter is a subset of the Dark solitons in BECs dataset (v2.0) [55]. This dataset was originally created to support the development and benchmarking of ML methods for cold-atom experiments and to provide a concrete testbed at the interface of ML and quantum physics [23]. It consists of more than 1.6 × 104 pre-processed experimental images of BECs after TOF, acquired under a range of conditions, some with and some without solitonic excitations present. Since the images are produced by a real experimental apparatus, they naturally include the variability and imperfections that make automated analysis challenging and scientifically relevant: fluctuations in atom number and cloud shape, imaging noise, and genuine shot-to-shot differences in the excitations themselves. As introduced in Sec. II, these data were obtained via absorption imaging after TOF[65], yielding the raw data as in Figs. 8(a)–(c). The condensate, visible in the red rectangle in Fig.8a, occupies only a small fraction of the field of view. Figures 8(d)–(f) show the optical depth computed using Eq. (2) after rotating by ϕ ≈ 45◦ to align with the long axis of the BEC, with rotated coordinates denoted by [iϕ , jϕ ]. To produce standardized inputs

Opt. density

PRL 113, 065302 (2014) (a) (a)

(b) (b)

week ending 8 AUGUST 2014

PHYSICAL REVIEW LETTERS

FIG. 7. (a)–(b) Experimental images of the integrated density distribution of two BECs after TOF of 120 ms, revealing the presence of a solitonic vortex with opposite circulation. (c) Theoretical 3D calculation for the experimental configuration, with clockwise circulation. In all panels, arrows indicate the atomic flow. Reproduced with permission from Ref. [52].

FIG. 1 (color online). Integrated triaxial density distribution of a BEC after a time of flight of 120 ms in presence of a solitonic vortex aligned along y. For each condensate we report absorption images along the two radial directions, horizontal y (left) and vertical z (top), and the residuals of axial (x) imaging after subtracting the Thomas-Fermi profile fit. (a)–(b) Experimental snapshots of two condensates with opposite circulation. (c) Theoretical 3D calculation with clockwise circulation and a μtheo ≃ μexp =3 (white noise has been added in order to better compare the theoretical calculation to the experimental results). Arrows indicate the atomic flow. Other examples of solitonic vortices with different orientation and shape can be seen in the Supplemental Material [16].

suitable for both manual labeling and ML training, each

optical density image is then fit to a column-integrated Thomas-Fermi profile

in the whole nodal plane, while in a solitonic vortex it has density profile and the velocity field would be isotropic opposite sign in the two half planes separated by the vortex [Fig. 2(a)]. Conversely, if the vortex is oriented along the line. A consequence is that, if a condensate containing a radial direction of an elongated trap, as in our case, the (" is released soliton equiphase surfaces, which originate from the vortex and overall density 2 # )3/2 from the2trap, the structure is preserved, must be orthogonal to the condensate outer surface, are iϕ −while i0 in the casejϕof−a jsolitonic 0 TF vortex the forced to bend and to gather in the shortnradial = direction − the two density , 0 n0 max 1 −asymmetric flow induces + δn, (16) i,j depleted sides to twist [Fig. 2(b)]. The phase gradient is then concentrated in a Riin opposite directions. RjFrom the standard definition of superfluid velocity, v ¼ ðℏ=mÞ∇ϕ, it planar region whose width is much smaller than the axial follows that the atomic flux is oriented towards increasing length of the condensate. The high phase gradient correphases; hence, the depleted regions move in the opposite sponds to a high atomic velocity field and the system tends direction. If at least one extreme of the vortex line is to reduce its energy by depleting the density in the plane. oriented along y or z a remarkable structure appears (see The resulting structure is a solitonic vortex [7,8]. Far from Fig. plane twisted physics-informed around a hollow vortex the density depleted region the phase pattern is similar to excitation (PIE) classifier, which uses providing estimates of the cloud center core. [i0 ,1): j0a],solitonic the peak Such a feature was already clearly visible in [1]. The that of a gray soliton [Fig. 2(c)], the two sides of the parameters of a Ricker wavelet (i.e., the Mexican hat 2Dcondensate density n0different , thephases Thomas–Fermi radiisign[Rofia, R along vortex, was previously observed exciting having (Δϕ ¼ π in the case of j ], which collective modes [17], following the vortex core precession a stationary solitonic offset vortex and soliton). A for key the wavelet) fitted to the residual density (condensate with an overall δna dark to correct small changes [18,19], or using interferometric techniques [20], is here difference is that in a soliton the phase gradient is the same density profile with Thomas-Fermi fit subtracted) to in probe intensity between images. The areorientation extractedimages from the twist in single absorption images. centered at assign one of six excitation categories within class-1. cropped to a fixed size of 164 × 132 pixels We compare our observations with the solution π ℏ2 2 [i0 , j0 ]. Finally, an elliptical mask (with ofsemi-major and the Gross-Pitaevskii (GP) equation The iℏ∂ t ψ ¼ − physics-informed feature representation used 2m ∇ ψþ 1 2 2 2 2 2 mω ðx þ α r Þψ þ gjψj ψ, where α ¼ ω =ω ⊥ x and g ¼ x 2 minor radii determined from the Thomas-Fermi fit) is throughout this section is illustrated in Fig. 9. Starting (a) 4πℏ2 as =m with as being the scattering length. The trapped applied to suppress background noise outside theiscloud. stationary state obtained by means of an the imaginary from ODtime image of a BEC containing a lone solitonic π of the GP equation, where the phase pattern of a Together, these steps enforce consistent evolution geometry across excitation (Fig. solitonic vortex is initially imprinted [7,8,21]. This state is 9a), we extract a 1D density profile along (b) used as the initialthe condition for simulation of axis the shots, reduce irrelevant pixel variation,thenand make thethe analysis and subtract the Thomas-Fermi backfree expansion. In order to speed up the simulation time and remaining intensity structure more directly comparable to isolate the localized depletion (Fig. 9b). This (c) make the calculation feasible also ground for long expansions, we dynamically rescale the GP equation by means of the is then fit with a simplified Mexican-hat across experimental conditions. residual profile Thomas-Fermi scaling ansatz in the transverse direction FIG. 2 (color online). Theoretical solutions of in situ density 9b, red). The model’s parameters provide [22–28]. A typical Dark result is shownmodel in Fig. 1(c)(Fig. for a threeand phase profiles for (a) a presence vortex aligned alongof the excitations—The axis of an Coarse labels: dimensional (3D) simulation of a condensate expanding for physically meaningful description of the axially symmetric potential, a solitonic vortex (b) and a soliton a compact and solitons BECs todataset provides several 10ωx and μtheo ¼ 10ℏω⊥ . The expand120 ms with ω⊥ ¼types (c) orientedin perpendicularly the axis of an(v2.0) elongated axially excitation, see ing for solitonic vortex clearly develops a planar density Fig. 9c. In particular, the amplitude A, indicate the atomic flow. ML models of symmetric labelstrap. toArrows support training automated

analysis. The first category of labels is a deliberately sim065302-2 ple set of human assigned coarse labels that capture the most robust visual distinctions available [21]: class-0 (Fig. 8d), indicating no excitations; class-1 (Fig. 8e), indicating a single, well-defined excitation; and class-2 (Fig. 8f), indicating shots that are neither unambiguously empty nor clean single-excitation cases (e.g., multiple stripes, additional density modulations, or ambiguous structures). This three-way split mirrors the practical decision an experimentalist typically makes when scanning large image sets and provides a stable foundation for both benchmarking and downstream, more detailed labeling. During data curation, an additional class was introduced, class-8, to indicate data flagged during the process as potentially mislabeled [23]. Roughly 60 % of the v2.0 is unlabeled and assigned to class-9. Fine labels: physics-informed excitation categories—The class-1 data, i.e., with exactly one identified excitation, are annotated with a finer, physically motivated taxonomy stored as the excitation PIE label [22]. This additional label is provided by a

width σ, position xc , and symmetric a and b asymmetric shoulder parameters derived independently from the top and bottom half of the image capture the key geometric and structural properties that distinguish different excitation types. The ten Mexican hat fit parameters are used by the PIE classifier. Figure 10 shows examples of BECs from all six PIE classes.

Problem setup

In the remainder of this chapter, we focus on the class-1 subset of the dataset. The fine-grained excitation PIE labels are strongly imbalanced, with the dominant longitudinal soliton category containing 2,229 images and the smallest clockwise vortices category containing only 28 images. Thus, we apply several curation steps common in data science to reduce the extreme imbalance while preserving the physics. First, we remove images labeled as canted, as a combination of orientation variability and small class size

j (px)

(b)

1.0 200

0.5

jφ (px)

250

500

(d)

100

250

500

(e)

250

500 i (px) 1.5

(f)

1.0

0.5

0.0 0

100

100 (g)

100

(h)

100

iφ (px)

100

iφ (px)

(i)

0 0

makes the category difficult to learn reliably. Second, we merge physically symmetric categories into consolidated classes, i.e., we combine top and bottom partial soliton classes into a single partial category, and clockwise and counterclockwise solitonic vortices into a single vortex category. To ensure consistency within the merged classes, the bottom partial and counterclockwise vortex images are transformed by symmetry operations so that their visual appearance is always top partial and clockwise vortex, respectively (i.e., the dataset contains no bottom partial or counterclockwise vortex examples). We refer to the resulting three-class subset as the reduced soliton dataset. The reduced set provides a realistic, experimentally grounded benchmark for comparing interpretable EBM pipelines trained on low-dimensional, physics-informed image representations with image-native CNN models. Our goal is to train a classifier that maps the observed data to a probability distribution over labels while remaining robust to common experimental imperfections, such as shot-to-shot variations in atom number, small translations, and residual imaging artifacts.

1.5

(c)

Optical Density

(Data scaled by 100)

(a)

400

OD (unitless)

FIG. 8. (a)–(c) Raw data showing (a) probe with atoms I, (b) probe only Ip , and (c) background Ibg (scaled by 100×). As in Fig. 1, red boxes outline the region where the BEC is located. This data was used to derive the optical depth (e). (d)–(f) Representative preprocessed images from the Dark solitons in BECs dataset 2.0 showing a BEC (d) without excitations (class-0), (e) with a lone excitation (class-1), and (f) with multiple excitations (class-2). The green arrows indicate all solitonic excitation positions. (g)–(i) Overall column-integrated optical depth (gray) and corresponding fluctuations away from the nominal ThomasFermi profile (black) for a BEC (g) without excitations, (h) with a lone excitation, and (i) with multiple excitations. The green dashed lines mark the location of the depletions in the density fluctuations corresponding to the solitonic excitation indicated in (d)–(f).

CCD Counts (×103 )

A brief introduction to Explainable Boosting Machines

PCA enforces ℓ2 -orthogonality among the principal components, making it a powerful tool for compactly representing high-dimensional images and capturing dominant modes of variance. However, orthogonality in pixel space does not necessarily translate to semantic complementarity of the underlying physical structures. Even the leading components might mix multiple sources of variation, such as condensate shape fluctuations and diffrac-

100

tion or fringe patterns, as seen in the first six principalcomponent images in Fig. 2b. This mixing complicates interpretation. Thus, while PCA provides an efficient coordinate system, it does not by itself yield a transparent explanation of which physically meaningful cues drive a downstream classification decision. To introduce interpretability, we turn to EBMs. EBMs are accurate, inherently interpretable predictive models that combine the structure of generalized additive models (GAMs) with the training power of modern boosting. GAMs are an extension of linear and generalized linear models[66], in which linear terms, such as βi xi , are replaced with learned 1D shape functions fi (xi ), yielding an additive score of the form FGAM (x) = β0 +

P X

fi (xi ),

(17)

i=1

where P is the number of features [58]. Since each feature’s contribution can be visualized directly by plotting fi , this structure preserves interpretability while still allowing nonlinear relationships. A limitation of purely additive models is that they cannot capture interactions, i.e., situations where the effect of one feature depends on the value of another. A standard extension is to augment the GAM score with a small number of pairwise interaction terms, FGA2 M (x) = β0 +

P X i=1

fi (xi ) +

fij (xi , xj ),

(18)

(i,j)∈I

where I is a small set of interacting feature pairs and each fij is a learned two-dimensional interaction function. This form, sometimes called a generalized additive model with pairwise (i.e., second-order) interactions

14 Optical Density 0.0

0.5

1.0

1.5

jφ (px)

Width

100

Asymmetric Shoulder Height

(a)

50 0

OD (unitless)

Symmetric Shoulder Height

100

Amplitude

(b)

f (x) = A[1 −

0 0

(x−xc )2 a (x − xc )2 + b(x − xc )]e− 2σ2 2 σ

100

(GA2 M, where “2” indicates the highest allowed interaction order), retains interpretability since the model can still be inspected through one-dimensional shape plots fi and a small number of two-dimensional interaction surfaces fij [26]. Training GAM and GA2 M models amounts to learning from data functions {fi }, and optionally {fij }, that minimize an appropriate loss function. In principle, one could fit these functions using splines or other smooth function classes; however, in practice, the functions must be flexible enough to capture sharp thresholds and non-smooth structure common in experimental data. Moreover, fitting many functions simultaneously can lead to different terms competing to explain the same variation, reducing stability and complicating interpretation. EBMs address these issues by learning the GAM/GA2 M form using modern boosting, constructing the functions {fi } and (optionally) {fij } through an iterative error-correction process. Given a dataset of N labeled examples D = {(X(i) , y (i) )}N i=1 , where X(i) ∈ RP is a P -dimensional feature vector (in our case, PIE-derived fit parameters and descriptors) and y (i) is the corresponding class label, we want to learn the score function F(x) of the GAM/GA2 M. The additive, linear sum of feature contributions is connected with the labels via the so-called link function g(·): y = g (F(X))

(19)

FIG. 9. Physics-informed feature extraction for solitonic excitations. (a) Image of a BEC with a solitonic excitation identified by the red arrows. (b) 1D density profiles, both with (black) and without (grey) Thomas-Fermi background subtraction. The red curve plots a fitted Mexican-hat model capturing the localized depletion structure. (c) Schematic of the fit-model illustrating the key five parameters: amplitude A, center xc , width σ, and symmetric and asymmetric shoulder heights a and b. These parameters provide a compact, physically interpretable description of the excitation and form the basis of the feature representation used for classification.

Link functions enable EBM to handle various response types, including classification probabilities, counts, and continuous regressions. For binary classification, one typically uses the logistic sigmoid, p̂(y = 1 | X) = σ(F (X)) =

2 Lsq (y, F(X)) = 12 y − F (X) .

(23)

Starting from a simple baseline model, the algorithm

(20)

For continuous regression, common choices include the identity link µ(X) = F(X) and positive-valued links such as µ(X) = exp(F(X)) (log link) or µ(X) = log 1 + eF (X) (softplus link). For count data (e.g., Poisson regression), a standard choice is the log link µ(X) = exp(F(X)), which enforces nonnegative predicted mean counts. The model parameters, i.e., the functions fi and fij , are fit by minimizing the average loss (also called an empirical risk ), N 1 X (i) L y , F(X(i) ) , Remp (F) = min F ∈G N i=1

(21)

where G denotes the GAM/GA2 M model class and L is a loss function. For binary classification, a standard choice is the logistic (cross-entropy) loss,

Llog (y, F(X)) = −y log σ(F (X)) − (1 − y) log(1 − σ(F (X))) ,

while for regression, a common choice is the squared loss,

1 . 1 + e−F (X)

(22)

repeatedly identifies what the current model still gets wrong and adds a small correction that improves the fit. In practice, residuals are often used in regression problems, and loss gradients are used in classification tasks.

j (px) (a)

100

(c)

1.5

(e)

1.0

0.5

0.0 0

100

(b)

100

(d)

100 1.5

(f)

1.0

0.5

0.0 0

100

In boosting-based training, the “error-correction” step can be made precise: at iteration t one computes pseudoresiduals (negative loss gradients) ∂ (i) (i) rt = − ℓ y ,z , ∂z z=Ft−1 (X(i) ) (i)

(24)

and then fits a weak learner to predict rt from the inputs. Crucially, each correction is built using a shallow decision tree and the weak learner is constrained to depend only on a single feature (when updating a main effect fi ) or a single feature pair (when updating an interaction term fij ); the correction is then absorbed into the corresponding fi or fij . Thus, each boosting step modifies only one interpretable component of the additive model. By accumulating many such small, constrained corrections, EBMs can learn highly nonlinear shapes and sharp thresholds while preserving the GAM/GA2 M additive structure, ensuring that the final prediction remains fully decomposable into interpretable per-feature and pairwise contributions. Importantly, EBMs enable not only global explanations in the form of plots of fi and fij , but also local, per-example breakdowns of contributions without relying on post-hoc interpretability methods. The same structure that makes EBMs interpretable also limits their utility: EBMs operate on explicit features, so their explanations are only as meaningful as the tabular representation provided to the model. While many ML methods can learn hierarchical features directly from pixels, EBMs do not “look” at raw images in a native way; they instead require semantically meaningful, low-dimensional descriptors whose individual coordinates can be named, plotted, and interpreted. Fortunately, in the solitonic-excitation setting, such a compact, yet human-intelligible, representation is given by the Mexican hat wavelet fitted to the condensate density profile [22]. This motivates the complementary use of CNNs, which learn feature representations directly from images without requiring a predefined tabular encoding.

100

Optical Density

FIG. 10. Example images of data from the six PIE classes: (a) a longitudinal soliton, (b) canted excitation, (c) top partial excitation, (d) bottom partial excitation, (e) clockwise vortex, and (f) counterclockwise vortex. The vorticity in (e) and (f) is assigned based on how the density fronts curve around the excitation axis. The green arrows indicate all excitation positions.

Optical Density

iφ (px)

A brief introduction to convolutional neural networks

CNNs are an image-native alternative to EBMs. Rather than requiring an explicit tabular representation, CNNs learn a hierarchy of spatial features directly. As an image passes through the layers, a CNN repeatedly transforms it into new arrays called activations or feature maps. These feature maps—the model’s internal representation—are learned intermediate encodings of the input that are useful for the task. Early layers tend to capture simple local patterns such as edges, blobs, or contrast changes; intermediate layers combine these into higher-level motifs (textures, repeated stripes, localized notches); and deeper layers form task-specific representations that are most informative for separating classes. This flexibility comes at a price. In a trained CNN, no single neuron or filter typically corresponds to a single human-meaningful concept; instead, evidence is distributed across many channels and layers, making it difficult to attribute decisions to a small set of interpretable features. CNNs construct intermediate feature maps by convolving an input image X with a bank of learned filters. Let F 0 ≡ X, and denote the output of layer ℓ by F ℓ ∈ RHℓ ×Wℓ ×Cℓ , where Hℓ × Wℓ are the spatial dimensions of the feature map and Cℓ is the number of channels. A typical convolutional block takes the form F ℓ+1 = φ(Wℓ ∗ F ℓ + bℓ ) ,

(25)

where ∗ denotes convolution, Wℓ and bℓ are learned weights and biases, and, similar to the dense layer in Eq. (15), φ(·) is a pointwise activation function [39]. Common interior activation functions include: Rectified linear unit (ReLU) ReLU(x) = max(0, x);

(26)

A common default in CNNs, ReLU is simple, fast, encourages sparse activations, and helps mitigate vanishing

16 gradients compared to saturating nonlinearities. However, for negative inputs, units sometimes get stuck at 0 (so-called “dying ReLU”); moreover, because ReLU is unbounded above, large activations can occur.

vector. This is often done using global pooling, which averages or maximizes each channel across its spatial dimensions, H

h(X)c = Leaky ReLU LReLU(x) = max(αx, x), where 0 < α ≪ 1;

(27)

Leaky ReLU introduces a small negative slope to reduce dying-ReLU behavior while retaining ReLU-like simplicity. The trade-off is a new hyperparameter (the negative slope); the output remains unbounded above.

ℓ X ℓ 1 X Fℓ (m, n, c), Hℓ Wℓ u=1 v=1

(31)

producing a feature vector h(X) ∈ RCℓ . This vector summarizes the presence and strength of the learned features across the entire image. For classification tasks, a final linear layer maps this feature vector to a set of output scores, or logits, z(X) = Wout h(X) + bout .

(32)

These logits are then converted into predictions using an appropriate output activation function: Gaussian error linear unit (GELU) GELU(x) = x Φ(x),

sigmoid (28)

where Z x Φ(x) = −∞

−1

P r(y = 1|X) = σ(z(X)) = [1 + exp (−z(X))]

; (33)

The sigmoid is used for binary classification and maps a scalar logit to a value between 0 and 1. √1 exp(−t2 /2) dt. 2π

(29)

GELU is a smooth nonlinearity widely used in Transformer-based models. It preserves a nonzero gradient for negative inputs and can improve optimization in some settings, at the cost of slightly higher computational complexity (the exact form involves the Gaussian cumulative distribution function and is often approximated).

In addition to convolution, CNN architectures commonly include pooling operations that reduce the spatial dimensions of feature maps while retaining the most important information. Pooling acts independently on each channel and replaces small spatial neighborhoods with a single summary value. For example, in max pooling with a window size k × k, Fℓ+1 (m, n, c) = max Fℓ (km + i, kn + j, c), i=1,...,k j=1,...,k

softmax P r(y = k|X) = softmax(z(X))k exp (zk (X)) = PK ; j=1 exp (zj (X))

The softmax is used for multi-class classification; it converts the logit vector into a probability distribution over K mutually exclusive classes.

linear (identity) I(X) = X.

(30)

where m and n denote the vertical and horizontal spatial indices, respectively, and c indexes the channel, reduces the spatial resolution by a factor of k in each direction. By progressively decreasing spatial resolution, pooling helps CNNs build representations that are more robust to small translations and reduces the computational cost of deeper layers. In modern architectures, pooling is sometimes replaced or supplemented by strided convolutions that achieve a similar downsampling effect. After several convolutional and pooling layers, the feature maps F ℓ encode increasingly abstract, spatially compressed representations of the input. To convert these representations into a prediction, the network first aggregates the spatial information into a compact feature

(34)

(35)

The trivial identity activation function is often used in regression tasks to predict continuous numerical values directly. Thus, CNNs map the input image X to a prediction through a sequence of learned feature transformations followed by a task-appropriate output link. Given a labeled dataset D = {(X(i) , y (i) )}N i=1 , the network parameters are learned by minimizing an empirical loss function. For classification tasks, a common choice is the cross-entropy loss, Lce (θ) = −

N 1 X log P rθ y (i) | X(i) , N i=1

(36)

17 where θ denotes all learned weights and biases. Optimization algorithms such as stochastic gradient descent adjust the parameters to minimize this loss, enabling the network to learn feature representations and classification rules directly from the image data. However, unlike EBMs, where the learned functions correspond directly to named input features, the features learned by CNNs are implicit and must be interpreted indirectly. Post-hoc tools such as gradient-based saliency maps [24] or GradCAM [25] can highlight pixels that influence the output, but these methods provide heuristic visualizations rather than the exact additive decomposition available in EBMs. We use CNNs as a high-performance baseline for comparing EBM performance and interpretability.

Results: classification benchmarking

Each image in the reduced soliton dataset is represented as a 2D array X∈RH×W , with H=164 and W =132, and the corresponding Mexican hat fit parameters introduced in Sect. IV B stored as a 1D array Y∈RP , with P =14. In the context of supervised ML, where each data point is accompanied by a corresponding label, we represent the dataset as D = {(X(i) , Y(i) , y (i) )}N i=1 , with labels y (i) ∈{0, 1, 2} indicating excitation type (longitudinal soliton, partial, and vortex, respectively). The workflow for EBM-based classification progresses in two stages. The first uses the reduced soliton dataset as defined in Sect. IV C and serves as a baseline for assessing how well an interpretable, physics-informed model can separate longitudinal solitons, partial excitations, and vortices. The second stage uses the same reduced dataset and augments the vortex class with physically meaningful, symmetry-preserving transformations to test whether improved representation of the rarest class leads to better classification performance. These two stages allow us to separate the baseline capability of the EBM from the gains achievable through targeted augmentation of the most underrepresented excitation type. The EBM performs strongly in both settings, and targeted augmentation leads to a clear improvement in aggregate performance in terms of both accuracy and the F1 score. Precision measures how often predictions assigned to a given class are correct; recall measures how often true examples of that class are successfully identified; and the F1-score summarizes the balance between precision and recall by taking their harmonic mean, which penalizes large discrepancies between the two; as a result, the score is high only when both precision and recall are simultaneously high. The balanced accuracy increases from 92.2 % on the baseline reduced dataset to 94.3 % after vortex augmentation, while the average F1-score rises from 94.0 % to 96.2 %. Taken together, these results indicate that the main limitation of the baseline model is not a lack of expressive power, but the scarcity of representative vortex examples in the training data. Detailed quantitative results are summarized in Ta-

TABLE I. Classification report for the reduced soliton dataset, comparing baseline performance of the EBM trained on the reduced dataset (Base.) with performance obtained after applying physically meaningful augmentations to the vortex class (Augm.). Metrics are grouped by soliton class. Metric

Longitudinal Partial Vortex Base. Augm. Base. Augm. Base. Augm. Precision 0.98 0.93 0.92 0.97 0.95 0.96 Recall 0.97 0.95 0.85 0.98 0.92 0.93 F1-score 0.97 0.94 0.88 0.98 0.94 0.94 Support 446 160 13 446 159 27 ROC-AUC 0.997 0.986 0.936 0.997 0.981 0.976 PR-AUC 0.993 0.994 0.998 0.994 0.993 0.999

ble I. For each case, we report the standard per-class classification metrics: precision, recall, F1-score, and support (the number of test examples in the class). We also report one-vs-rest[67] ROC-AUC (receiver operating characteristic–area under the curve) and PR-AUC (precision-recall–area under the curve), which quantify how well the model separates the target class from all others across classification thresholds. On the baseline reduced dataset, the EBM performs strongly on the two majority classes, achieving F1-scores of 0.97 for longitudinal solitons and 0.94 for partial excitations. Performance on the vortex class is also encouraging, but noticeably weaker, with precision 0.92, recall 0.85, and F1-score 0.88. This behavior is consistent with the dataset’s strong class imbalance: while the model learns reliable decision boundaries for the dominant classes, it has fewer examples from which to infer a stable representation of vortices. Applying physically meaningful augmentation to the vortex class leads to clear improvements across all vortex metrics, with recall increasing from 0.85 to 0.93, and the F1-score from 0.88 to 0.95. At the same time, performance on the longitudinal and partial classes remains essentially unchanged, indicating that augmentation enhances minority-class recognition without degrading performance on the majority classes. These results suggest that the EBM is not fundamentally limited by its expressive power on this task; rather, a significant portion of the remaining error can be attributed to the scarcity of representative vortex examples in the training data. This overall improvement is driven almost entirely by better identification of the vortex class, rather than by changes in the already strong performance on the two majority classes. A key advantage of EBMs is that their decision process can be inspected directly in terms of the features used for training. Figure 11 shows both global and local aspects of this interpretability: the top row ranks features by mean absolute score, indicating which fitted parameters matter most for each class decision, while the bottom row shows representative shape functions that quantify how changes in a given feature shift the model score, where the ver-

Ordered features

(a) longitudinal

(b) partial

A(t)

A(b)

b(t)

(t)

A(t)

(t)

b(b)

(t)

a(t)

(b)

σ (b)

(b)

0.0 1.0 2.0 Mean absolute score

Score

(d) 5

0.0 1.0 2.0 Mean absolute score

(e) 5

(f)

-5 -1.5 A(t)

2.0

1 0

-1

-5 -5.0

0.0 1.0 2.0 Mean absolute score

-5.0

-1.5 A(b)

2.0

-2.5

tical axis corresponds to the additive contribution of the feature to the class-specific logit (score). For this analysis separate Mexican-hat fits were performed to density profiles from the top and bottom portions of the image (extending the single-profile fits in Fig. 9); the resulting fit parameters are labeled by (t) and (b). For all three classes, the most prominent features are the excitation amplitudes A(t/b) together with the shoulder parameters a(t/b) and b(t/b) , indicating that the model relies primarily on the depth and asymmetry of the density depletion when distinguishing excitation types. The first notable difference is in the global featureimportance rankings in Fig. 11a–c. For the longitudinal and partial classes, the most important parameters are A(t) , a(t) , and b(t) , indicating that classification is driven primarily by the amplitude and detailed shape of the excitation in the top part of the condensate, together with the A(b) parameter capturing the vertical asymmetry in the condensate. By contrast, the vortex class is dominated by A(b) , b(t) , A(t) , and b(b) , showing that vortex identification relies most strongly on the combination of excitation amplitude and asymmetric shoulder structure. In the language of Fig. 9c, this suggests that vortices are distinguished by the interplay between depletion depth and asymmetric distortion, consistent with their more complex density signatures. The shape-function in Fig. 11d–e provides a particularly clear interpretation of the model’s decision logic. The unequal contributions of A(t) and A(b) reflect the underlying bias in the dataset—with only top partial solitons—capturing the directional asymmetry in the excitation. Longitudinal solitons exhibit approximately symmetric density depletions, whereas partial excitations and vortices are intrinsically asymmetric. Consequently, the relative imbalance between A(t) and A(b) provides a key discriminative signal, and the model correctly assigns different contributions to these features to capture this structure. As shown in Fig. 11d–e, increasing

0.0 xc(b)

2.5

FIG. 11. EBM results. Top row: global feature importance, quantified by the mean absolute contribution of each feature to the model score for (a) longitudinal solitons, (b) partial excitations, and (c) vortices. The ordering highlights which parameters are most influential in distinguishing each excitation type. Bottom row: representative feature curves for (b) (d) A(t) , (e) A(b) , and (f) xc , illustrating how variations in individual parameters affect the model score. The vertical axis (“Score”) represents the contribution of the corresponding feature to the model output: positive values increase the likelihood of the class, while negative values decrease it. The three colors in panels (d)– (f) correspond to the longitudinal solitons (green), partial excitations (orange), and vortices (green) labels, respectively.

A(t) strongly favors the partial-excitation class while suppressing the longitudinal class, whereas increasing A(b) favors the longitudinal and, to some extent, vortex class, while penalizing partial excitations. Importantly, the model does not rely on these amplitudes independently, but rather on their combined and nonlinear contributions, allowing it to distinguish symmetric from asymmetric excitation profiles. In contrast to the amplitude-based features, the po(b) sitional feature xc exhibits a nearly flat contribution across the range of interest, indicating that it plays a negligible role in classification. This is also consistent with the underlying physics: while solitonic excitations can appear at different locations within the condensate, their classification is determined primarily by their shape and symmetry rather than their absolute position. Thus, (b) the EBM correctly identifies xc as a non-discriminative feature in this setting, further supporting the alignment between the learned model and physical intuition. We tested two CNN-based classification architectures. The first is a simple CNN with three convolutional layers followed by a fully connected layer established a baseline for image-native learning on the reduced soliton dataset. The second is a more advanced architecture specifically designed for solitonic-excitation classification, adapted from Ref. [21]. These two configurations allow us to assess how architectural choices impact performance. Since the models used in this chapter are not hyperparameteroptimized, all experiments are performed with a fixed training set and a held-out test set, without a separate validation set. For the simple CNN, we report average performance across 10 independent training runs to account for variability from random initialization and stochastic optimization. The model achieves a mean balanced accuracy of 78.3(1.3) % and an F1-score of 0.82(1). Unsurprisingly, performance is strongest for the longitudinal class, with recall 0.96(1) and F1-score 0.94(1), and moderate

19 TABLE II. Classification report for CNN-based models on the reduced soliton dataset, comparing the mean performance over 10 runs for a simple CNN architecture (Simp.) with the best observed performance for the CNN architecture adapted from Ref. [21] (Mod.). Metrics are grouped by soliton class. Metric

Longitudinal Partial Vortex Simp. Mod. Simp. Mod. Simp. Mod. Precision 0.92(1) 0.94 0.83(4) 0.86 0.90(8) 0.79 Recall 0.96(1) 0.96 0.77(3) 0.80 0.62(4) 0.85 F1-score 0.94(1) 0.95 0.80(2) 0.83 0.73(2) 0.82 Support 446 446 159 159 27 27 ROC-AUC 0.955(5) 0.953 0.942(6) 0.943 0.984(5) 0.983 PR-AUC 0.972(2) 0.969 0.888(13) 0.894 0.824(24) 0.819

icant additional labeling effort, especially for the rarest excitation types. Moreover, CNN’s internal decision rules are not directly accessible. Saliency maps [24] and GradCAM [25] can highlight pixels influential for a given prediction, but they do not provide an explicit, lowdimensional model ofwhich physical cues drive the decision or how different cues trade off. The EBM, by contrast, incorporates domain knowledge through its feature representation, reducing the reliance on large labeled datasets and enabling strong performance in a dataconstrained setting.

for partial excitations, with recall 0.77(3) and F1-score 0.80. The vortex class remains the most challenging, with a recall 0.62(4) and an F1-score 0.73(2). This behavior reflects the dataset’s class imbalance: the network learns robust representations for the dominant classes but struggles to generalize to the relatively scarce vortex examples. The architecture adapted from Ref. [21] delivers substantially improved performance. The best observed model achieves a balanced accuracy of 87.0 % and an F1-score of 0.87. Both minority classes benefit from the improved architecture, with particularly notable gains for the vortex class, whose recall increases from 0.62(4) to 0.85 and F1-score from 0.73(2) to 0.82. These results indicate that, when provided with sufficient architectural capacity and appropriate inductive bias, CNNs can learn discriminative representations directly from image data, even for underrepresented classes. Despite these improvements, the performance of the CNN models is consistently below that of the EBM-based approach for every metric considered, including ROCAUC and PR-AUC, with the largest advantage in PRAUC for the minority classes. The uniformly high values of the ROC-AUC and PR-AUC for all classes, as shown in Table II, indicate that the model assigns well-separated scores to each excitation type even before a specific decision threshold is chosen. The improvement in vortex ROC-AUC from 0.936 to 0.976 after augmentation shows that the augmented EBM ranks vortex examples more reliably relative to non-vortex examples, consistent with the observed increase in vortex recall and F1-score. However, the lower PR-AUC for vortices, about 0.82 regardless of the CNN architecture, shows that this separability does not translate into equally strong precision–recall behavior for the rare vortex class. This gap highlights that, for this dataset, the physicsinformed feature representation provides a more effective basis for classification than learning directly from pixels alone. In some sense, identifying the feature representation is part of the overall architecture of an EBM. While CNN performance would likely improve with substantially larger labeled datasets, achieving such improvements would require targeted data collection and signif-

DISCUSSION AND OUTLOOK

A central message of this chapter is that the effectiveness of machine learning methods in scientific settings depends critically on the quality, structure, and size of the underlying dataset. In both examples considered here, even relatively large datasets benefited from carefully designed augmentation strategies that improved coverage of relevant physical features. Equally important is the use of validation—and ideally independent test-sets—to ensure that reported performance reflects genuine generalization rather than overfitting. Our first example of reference reconstruction illustrates that relatively simple linear-algebra-based methods can perform competitively while remaining interpretable. In these “models”, the learned components can be directly visualized as images, providing insight into correlated variations in the data. More complex deep learning approaches yielded only modest performance gains, at the expense of interpretability. Our second example, soliton classification, provides an even stronger contrast: the physics-informed, interpretable model not only matched but exceeded the performance of the more complex CNNbased approach. So, “can machine learning be explainable?” is a resounding yes. Results in this section demonstrate that machine learning models can be both accurate and interpretable, and that the commonly assumed “tradeoff” between performance and explainability need not hold in all settings. In particular, when domain knowledge can be incorporated into the representation—as in the case of physics-informed features—interpretable models may achieve strong performance even in data-limited regimes. At the same time, it is important to recognize that different modeling approaches scale differently with data. Deep learning methods are likely to benefit more strongly from large, well-labeled datasets, whereas interpretable, feature-based models can achieve robust performance with fewer labeled examples by leveraging prior knowledge but may require tabular representations. In practice, the choice of model therefore depends not only on achievable accuracy, but also on the availability of labeled data, the cost of data curation, and the need for interpretability and validation. In such cases, you, the scientist, must make a judgment call weighing the rela-

20 tive value of raw performance against understanding the machine learning decision process, accepting that, on average, high-performing black-box solutions often lead to unexpected outcomes in seemingly random situations.

ACKNOWLEDGMENTS

Advancing machine learning for scientific applications will require continued attention to dataset design, feature representation, and evaluation protocols (including choosing appropriate loss functions and determining relevant success measures) that emphasize not only performance but also reliability, interpretability, and consistency with underlying physical principles. Such considerations will be essential for establishing trust in machine learning systems deployed in scientific discovery and measurement-driven domains.

IBS thanks E. Gvozdiovas, A. M. Piñeiro and M. Zhao for assistance in assembling datasets for the reference recovery section. This work was partially supported by the National Institute of Standards and Technology. The views and conclusions contained in this paper are those of the authors and should not be interpreted as representing the official policies, either expressed or implied, of the U.S. Government. The U.S. Government is authorized to reproduce and distribute reprints for Government purposes notwithstanding any copyright noted herein. Any mention of commercial products is for information only; it does not imply recommendation or endorsement by NIST.

[1] C. D. Bruzewicz, J. Chiaverini, R. McConnell, and J. M. Sage. Trapped-ion quantum computing: Progress and challenges. Appl. Phys. Rev., 6(2):021314, 2019. [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] A. W. Harrow and A. Montanaro. Quantum computational supremacy. Nature, 549(7671):203–209, 2017. [4] G. Carleo, I. Cirac, K. Cranmer, L. Daudet, M. Schuld, N. Tishby, L. Vogt-Maranto, and L. Zdeborová. Machine learning and the physical sciences. Rev. Mod. Phys., 91:045002, 2019. [5] W. Ketterle, D. S. Durfee, and D. M. Stamper-Kurn. Making, probing and understanding Bose-Einstein condensates, pages 67–176. Proceedings of the International School of Physics ”Enrico Fermi”. IOS Press, 1999. [6] J. Denschlag, J. E. Simsarian, D. L. Feder, C. W. Clark, L. A. Collins, J. Cubizolles, L. Deng, E. W. Hagley, K. Helmerson, W. P. Reinhardt, et al. Generating Solitons by Phase Engineering of a Bose-Einstein Condensate. Science, 287(5450):97–101, 2000. [7] R. Carretero-González, D. J. Frantzeskakis, and P. G. Kevrekidis. Nonlinear waves in Bose–Einstein condensates: physical relevance and mathematical techniques. Nonlinearity, 21(7):R139, 2008. [8] D. Zha, Z. P. Bhat, K.-H. Lai, F. Yang, Z. Jiang, S. Zhong, and X. Hu. Data-centric artificial intelligence: A survey. ACM Comput. Surv., 57(5), 2025. [9] C. Rudin. Stop explaining black box machine learning models for high stakes decisions and use interpretable models instead. Nat. Mach. Intell., 1(5):206–215, 2019. [10] Z. Vendeiro, J. Ramette, A. Rudelis, M. Chong, J. Sinclair, L. Stewart, A. Urvoy, and V. Vuletić. Machinelearning-accelerated Bose-Einstein condensation. Phys. Rev. Res., 4:043216, 2022. [11] A. Impertro, J. F. Wienand, S. Hä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. Commun. Phys., 6(1):166, 2023. [12] A. Bohrdt, C. S. Chiu, G. Ji, M. Xu, D. Greif, M. Greiner, E. Demler, F. Grusdt, and M. Knap. Classifying snap-

shots of the doped Hubbard model with machine learning. Nat. Phys., 15(9):921–924, 2019. [13] N. Käming, A. Dawid, K. Kottmann, M. Lewenstein, K. Sengstock, A. Dauphin, and C. Weitenberg. Unsupervised machine learning of topological phase transitions from experimental data. Mach. Learn.: Sci. Technol., 2(3):035037, 2021. [14] C. L. Degen, F. Reinhard, and P. Cappellaro. Quantum sensing. Rev. Mod. Phys., 89:035002, 2017. [15] E. Altman, K. R. Brown, G. Carleo, L. D. Carr, E. Demler, C. Chin, B. DeMarco, S. E. Economou, M. A. Eriksson, K.-M. C. Fu, et al. Quantum simulators: Architectures and opportunities. PRX Quantum, 2:017003, 2021. [16] Y. Alexeev, D. Bacon, K. R. Brown, R. Calderbank, L. D. Carr, F. T. Chong, B. DeMarco, D. Englund, E. Farhi, B. Fefferman,et al. Quantum computer systems for scientific discovery. PRX Quantum, 2:017001, 2021. [17] L. Pezzè, A. Smerzi, M. K. Oberthaler, R. Schmied, and P. Treutlein. Quantum metrology with nonclassical states of atomic ensembles. Rev. Mod. Phys., 90:035005, Sep 2018. [18] M. F. Parsons, A. Mazurenko, C. S. Chiu, G. Ji, D. Greif, and M. Greiner. Site-resolved measurement of the spincorrelation function in the Fermi-Hubbard model. Science, 353(6305):1253–1256, 2016. [19] G. Theocharis, A. Weller, J. P. Ronzheimer, C. Gross, M. K. Oberthaler, P. G. Kevrekidis, and D. J. Frantzeskakis. Multiple atomic dark solitons in cigarshaped Bose-Einstein condensates. Phys. Rev. A, 81(6):063604, 2010. [20] A. R. Fritsch, Mingwu Lu, G. H. Reid, A. M. Piñeiro, and I. B. Spielman. Creating solitons with controllable and near-zero velocity in Bose-Einstein condensates. Phys. Rev. A, 101(5):053629, 2020. [21] S. Guo, A. R. Fritsch, C. Greenberg, I. B. Spielman, and J. P. Zwolak. Machine-learning enhanced dark soliton detection in Bose–Einstein condensates. Mach. Learn.: Sci. Technol., 2(3):035020, 2021. [22] S. Guo, S. M. Koh, A. R. Fritsch, I. B. Spielman, and J. P. Zwolak. Combining machine learning with physics: A framework for tracking and sorting multiple dark solitons. Phys. Rev. Research, 4(2):023163, 2022. [23] A. R. Fritsch, S. Guo, S. M. Koh, I. B. Spielman, and J. P.

21 Zwolak. Dark solitons in Bose–Einstein condensates: a dataset for many-body physics research. Mach. Learn.: Sci. Technol., 3(4):047001, 2022. [24] K. Simonyan, A. Vedaldi, and A. Zisserman. Deep inside convolutional networks: Visualising image classification models and saliency maps. arXiv:1312.6034, 2013. [25] R. R. Selvaraju, M. Cogswell, A. Das, R. Vedantam, D. Parikh, and D. Batra. Grad-CAM: Visual Explanations from Deep Networks via Gradient-Based Localization. Int. J. Comput. Vis., 128(2):336–359, 2020. [26] Y. Lou, R. Caruana, J. Gehrke, and G. Hooker. Accurate intelligible models with pairwise interactions. In Proceedings of the 19th ACM SIGKDD International Conference on Knowledge Discovery and Data Mining, pages 623– 631, 2013. [27] H. Nori, S. Jenkins, P. Koch, and R. Caruana. InterpretML: A Unified Framework for Machine Learning Interpretability. arXiv:1909.09223, 2019. [28] E. Altuntas and I. B. Spielman. Self-Bayesian aberration removal via constraints for ultracold atom microscopy. Phys. Rev. Research, 3:043087, 2021. [29] M. Zhao, J. Tao, and I. B. Spielman. Kolmogorov scaling in turbulent 2d Bose-Einstein condensates. Phys. Rev. Lett., 134:083402, 2025. [30] X. Li, M. Ke, B. Yan, and Y. Wang. Reduction of interference fringes in absorption imaging of cold atom cloud using eigenface method. Chin. Opt. Lett., 5(3):128–130, 2007. [31] S. R. Segal, Q. Diot, E. A. Cornell, A. A. Zozulya, and D. Z. Anderson. Revealing buried information: Statistical processing techniques for ultracold-gas image analysis. Phys. Rev. A, 81:053601, 2010. [32] L. Niu, X. Guo, Y. Zhan, X. Chen, W. M. Liu, and X. Zhou. Optimized fringe removal algorithm for absorption images. Appl. Phys. Lett., 113(14):144103, 2018. [33] G. Ness, A. Vainbaum, C. Shkedrov, Y. Florshaim, and Y. Sagi. Single-exposure absorption imaging of ultracold atoms using deep learning. Phys. Rev. Applied, 14:014011, 2020. [34] Da-Wei Ying, Si-Hui Zhang, Shu-Jin Deng, and Hai-Bin Wu. Single shot imaging for cold atoms based on machine learning. Acta Phys. Sin., 72(14):144201, 2023. [35] K. Lee and Y. Shin. Dual-species atomic absorption image reconstruction using deep neural networks. arXiv:2508.12120, 2025. [36] In our experience, this power-law scaling is ubiquitous, but we are not aware of the underlying mechanism. [37] For conventional PCA, this curve can be obtained by cumulatively summing the eigenvalues. This does not generalize to the optimal linear method, which is not built from a simple sum over orthogonal components. [38] To use standard SVD tooling with our large dataset, the initial 644 × 484 images were down-sampled to 322 × 242. [39] I. Goodfellow, Y. Bengio, and A. Courville. Deep learning. MIT press Cambridge, 2016. [40] O. Ronneberger, P.Fischer, and T. Brox. U-net: Convolutional networks for biomedical image segmentation. In Medical Image Computing and Computer-Assisted Intervention (MICCAI), volume 9351 of LNCS, pages 234– 241. Springer, 2015. [41] As opposed to the U-Net with shortcut connections, this structure is better able to filter noise processes such as shot noise and readout noise. [42] K. He, X. Zhang, S. Ren, and J. Sun. Deep residual

learning for image recognition. In 2016 IEEE Conference on Computer Vision and Pattern Recognition (CVPR), pages 770–778, 2016. [43] J. S. Russel. Report of the Committee on Waves, pages 417–468. British reports VI, 1837. [44] M. Lakshmanan. Tsunamis and Oceanographical Applications of Solitons, pages 8506–8521. Springer New York, New York, NY, 2009. [45] A. Hasegawa. Soliton-based optical communications: an overview. IEEE Journal of Selected Topics in Quantum Electronics, 6(6):1161–1172, 2000. [46] M. Kono and M. M. Škorić. Nonlinear Physics of Plasmas, volume 62 of Springer Series on Atomic, Optical, and Plasma Physics. Springer Berlin, Heidelberg, 2010. [47] F. Dalfovo, S. Giorgini, L. P. Pitaevskii, and S. Stringari. Theory of Bose-Einstein condensation in trapped gases. Rev. Mod. Phys., 71(3):463–512, 1999. [48] D. J. Frantzeskakis. Dark solitons in atomic BoseEinstein condensates: from theory to experiments. J. of Phys. A, 43(21):213001, 2010. [49] T. Busch and J. R. Anglin. Motion of Dark Solitons in Trapped Bose-Einstein Condensates. Phys. Rev. Lett., 84(11):2298–2301, 2000. [50] This class includes any feature whose projection is a clean density depletion, potentially including: kink solitons, longitudinally aligned solitonic vortices, and large-radius vortex-rings to name a few. [51] B. P. Anderson, P. C. Haljan, C. A. Regal, D. L. Feder, L. A. Collins, C. W. Clark, and E. A. Cornell. Watching Dark Solitons Decay into Vortex Rings in a Bose-Einstein Condensate. Phys. Rev. Lett., 86(14):2926–2929, 2001. [52] S. Donadello, S. Serafini, M. Tylutki, L. P. Pitaevskii, F. Dalfovo, G. Lamporesi, and G. Ferrari. Observation of solitonic vortices in Bose-Einstein condensates. Phys. Rev. Lett., 113(6):065302, 2014. [53] L. M. Aycock, H. M. Hurst, D. K. Efimkin, D. Genkina, H.-I. Lu, V. M. Galitski, and I. B. Spielman. Brownian motion of solitons in a Bose-Einstein condensate. Proc. Natl Acad. Sci., 114(10):2503–2508, 2017. [54] A. M. Mateo and J. Brand. Stability and dispersion relations of three-dimensional solitary waves in trapped BoseEinstein condensates. New J. Phys., 17(12):125013, 2015. [55] J. P. Zwolak, S. Guo, A. R. Fritsch, and Ian B. Spielman. Dark solitons in BECs dataset 2.0. National Institute of Standards and Technology, 2021. [56] After accounting for the magnification of our imaging system, the 648 × 488 pixel raw images (Point Grey FL3, 5.6 µm pixel pitch) have an effective pixel size of 0.93 µm, appreciably smaller than our 2.8 µm optical resolution. [57] Generalized linear models are a broad class of models that extend ordinary linear regression to outcomes that are not well-modeled by a Gaussian with constant variance. [58] T. Hastie and R. Tibshirani. Generalized Additive Models. Stat. Sci., 1(3):297–310, 1986. [59] In the one-vs-rest setting, each class is treated as the positive class in turn, with all remaining classes grouped as negatives. [60] In our experience, this power-law scaling is ubiquitous, but we are not aware of the underlying mechanism. [61] For conventional PCA, this curve can be obtained by cumulatively summing the eigenvalues. This does not generalize to the optimal linear method, which is not built from a simple sum over orthogonal components.

22 [62] To use standard SVD tooling with our large dataset, the initial 644 × 484 images were down-sampled to 322 × 242. [63] As opposed to the U-Net with shortcut connections, this structure is better able to filter noise processes such as shot noise and readout noise. [64] This class includes any feature whose projection is a clean density depletion, potentially including: kink solitons, longitudinally aligned solitonic vortices, and large-radius vortex-rings to name a few. [65] After accounting for the magnification of our imaging

system, the 648 × 488 pixel raw images (Point Grey FL3, 5.6 µm pixel pitch) have an effective pixel size of 0.93 µm, appreciably smaller than our 2.8 µm optical resolution. [66] Generalized linear models are a broad class of models that extend ordinary linear regression to outcomes that are not well-modeled by a Gaussian with constant variance. [67] In the one-vs-rest setting, each class is treated as the positive class in turn, with all remaining classes grouped as negatives.

Record · ID 200483 · SHA-256 6006b40e76f78f82
Retrieved via Conceptio — every document is proof-bundled with source, license, and retrieval metadata.