Thermometry of simulated Bose–Einstein condensates using machine learning

Since the first experimental realisation of Bose–Einstein condensation in dilute atomic gases [1–3]—one of the most precisely controllable quantum many-body systems—these systems have become important platforms for quantum simulation [4–6], precision metrology [7–9], and emerging quantum technologies [10, 11]. Central to these applications is the precise thermodynamical characterisation—particularly values of the temperature and chemical potential—which, in thermal equilibrium under the grand canonical ensemble, fully determine the system’s quantum statistical state and collective behaviour.

Conventionally, the temperature and chemical potential in ultracold atomic experiments are extracted using destructive time-of-flight imaging techniques [12–14]. In this method, the trapping potential is abruptly switched off, allowing the atomic cloud to expand freely. The temperature is then inferred from the spatial distribution of the cloud, assuming a Maxwell–Boltzmann velocity distribution. For a cloud with a Gaussian density profile, the $1/e^2$ width in the trap, σ₀, and the $1/e^2$ width after a time of flight t, σ, are related by $\sigma^2 = \sigma_0^2 + (k_\mathrm{B} T / m)\, t^2$ , where m is the atomic mass. However, this approach is inherently destructive and suffers from significant shot-to-shot variability, limiting experimental reproducibility and precision. Furthermore, the Maxwell–Boltzmann approximation is fundamentally unsuitable for describing quantum degenerate gases, which exhibit significant quantum statistical correlations and non-trivial interactions.

Minimally invasive thermometry approaches that use impurity-based polarons have been proposed and shown to be capable of nanokelvin and subnanokelvin precision by analysing impurity fluctuations in momentum and position [15–17].

In this paper, we demonstrate a proof-of-principle, minimally destructive thermometric technique using machine learning. Our method directly estimates temperature and chemical potential from in situ density profiles, facilitating repeated, rapid measurements on the same atomic sample without significantly perturbing it. Our approach is inspired by recent advances in machine learning to classify and characterise quantum fluid states [18–21]. Fundamentally, we learn the highly non-trivial mapping from a condensate’s atomic density to the thermodynamical parameters of interest—this is beyond what is possible in naïve fitting procedures currently used in cold atom experiments. We train our machine learning model using simulated density distributions generated from the stochastic Gross–Pitaevskii equation (SGPE), a theoretical framework for modelling finite-temperature dynamics and thermal fluctuations in Bose gases.

We validate our approach across a broad range of temperatures and chemical potentials, and evaluate its estimating capability using unseen density profiles with varied trapping geometries (despite only being trained on harmonically trapped condensates) and during thermalisation (despite only being trained on thermally equilibrated density profiles). Our results demonstrate the robustness and versatility of machine learning thermometry, suggesting its potential as a powerful, precise, and experimentally feasible tool for the investigation of ultracold atomic systems.

2.1. Dynamical treatment of the condensate

The theoretical treatment of a dilute Bose gas begins with the many-body Hamiltonian

Equation (1)

where

Equation (2)

incorporates the single-particle kinetic energy and the external trapping potential $V(\mathbf{x})$ , and $U(\mathbf{x-x^{\prime}})$ represents the two-body interatomic potential, where the factor 1/2 prevents double-counting of particle interactions.

Despite the strong short-range interactions characteristic of alkali atoms, ultracold gases of these species can be treated as weakly interacting systems [22]. At very low temperatures, the de Broglie wavelength λ greatly exceeds the range of the interatomic potential, permitting the approximation of interactions as elastic, point-like contacts characterised by the s-wave scattering length. This leads to the pseudopotential approximation $U(\mathbf{x}-\mathbf{x}^{\prime}) = g\delta(\mathbf{x}-\mathbf{x}^{\prime})$ , where the interaction strength $g = {4\pi\hbar^2a_s}/{m}$ is directly proportional to the scattering length a_s.

The Heisenberg equation of motion for the Bose field operator $\hat{\Psi}$ is given by

Equation (3)

We decompose the field operator into condensate and non-condensate contributions [23]:

Equation (4)

where $\hat{\Phi}(\mathbf{x},t)$ describes the condensate atoms and $\hat{\delta}(\mathbf{x},t)$ represents non-condensate atoms (thermal excitations, quantum fluctuations, or both); unlike the field operator $\hat{\Psi}(\mathbf{x})$ , condensate dynamics mean these may both have an explicit time-dependence. In the limit of large condensate occupation, the ensemble average reduces to a classical field: $\langle\hat{\Psi}(\mathbf{x})\rangle = \Phi(\mathbf{x},t)$ , such that $\hat{\Psi}(\mathbf{x}) = \Phi(\mathbf{x},t) + \hat{\delta}(\mathbf{x},t)$ .

Substituting this decomposition into (3) and taking the expectation value yields a form of nonlinear Schrödinger equation

Equation (5)

In the limit of negligible fluctuations, this reduces to the Gross–Pitaevskii equation, which provides an accurate description for weakly interacting gases with large atom numbers at temperatures typically below half the critical temperature for Bose–Einstein condensation.

2.2. Dynamics of Bose gases at finite temperatures

To investigate Bose gas dynamics across a broader temperature range, including near the phase transition, we require a treatment that incorporates thermal fluctuations beyond the mean-field approximation of (5).

We choose to model the growth of the condensate with the SGPE [24–27], which describes a static coupling of the condensate modes to a thermal bath. The growth is towards a state that corresponds to a thermodynamical equilibrium in the grand canonical ensemble defined by our choice of µ and T, which, by definition, must take real values. In essence, the SGPE is a phenomenologically damped Gross–Pitaevskii equation with additive noise. This approach implements a fluctuation–dissipation theorem that ensures the system relaxes to the correct equilibrium state, without spurious enhancement or depletion of either condensate or thermal components. The SGPE takes the form

Equation (6)

where

Equation (7)

is commonly referred to as the Gross–Pitaevskii ‘Hamiltonian,’ and the Gaussian noise correlations satisfy

Equation (8)

note $\eta(\mathbf{x},t)$ is a classical noise term, hence the angle brackets in (8) describe a purely statistical averaging, and not an expectation value over a quantum state.

We use a fourth-order Runge–Kutta scheme to solve equation (6) on a discretised spatio-temporal grid. There is an implicit projection out of high-momentum modes which are beyond the maximum momentum the spatial grid may represent. Under the numerical scheme as described in equations (6)–(8), we do not consider any interactions with the thermal cloud, i.e. we neglect an additive term $2g\tilde{n}(x)$ in equation (7) for the incoherent mean-field contribution (where $\tilde{n}$ is the non-condensate density and the factor of 2 arises from exchange symmetry under a Hartree–Fock theory). Such a term is important for quantitative agreement with experiments; for further discussion see e.g. [28, 29].

Since our focus lies in generating equilibrium thermal states rather than studying formation dynamics, the precise value of the dimensionless coupling parameter γ is not critical. This parameter controls how quickly the system equilibrates with its thermal environment—smaller values lead to slower equilibration, while larger values speed it up proportionally. However, γ must be chosen carefully: too small and the simulation becomes computationally impractical; too large and topological defects can form, actually slowing convergence. We use a spatially uniform rate γ = 0.01 throughout our simulations, which provides a good balance between computational efficiency and physical accuracy.

We evolve the system until we reach thermal equilibrium, as monitored by the relative change in condensate number:

Equation (9)

We assume equilibrium to have been achieved when $\Delta_N \lt 10^{-4}$ for five consecutive time steps.

For machine learning purposes, we partition the resulting atomic density profiles into three datasets: 80% for training, 10% for validation during each epoch, and 10% for final testing. The model is never exposed to the validation or test datasets during training. We use the training set to optimise the model parameters, the validation set to monitor generalisation and prevent overfitting during training, and the test set to provide an unbiased estimate of final performance on unseen data.

2.3. Atomic species and trap geometry

We consider ultracold, single-species atomic Bose gases of rubidium 87 with scattering length $a_\mathrm{s} = 5.29\,\mathrm{nm}$ [30] and atomic mass $1.44\times 10^{-25}\,\mathrm{kg}$ . The temperatures are sampled from $T\sim\mathrm{Uniform}(1200)\,\mathrm{nK}$ , encompassing both deeply degenerate and near-critical regimes. Chemical potentials are sampled from $\mu\sim\mathrm{Uniform}(20,80)\,\mathrm{nK}$ rather than from a regular grid to ensure robust model generalisation. Additional parameters are documented in our open-source code [31].

We use highly anisotropic trapping geometries. For harmonic traps, the potential is $V(x, y, z) = m(\omega_x^2 x^2+\omega_y^2 y^2+\omega_z^2 z^2)/2$ with transverse frequencies $\omega_x = \omega_y = 2\pi\times 25\,\mathrm{Hz}$ and axial frequency $\omega_z = 100\omega_x$ . We introduce a slight random anisotropy³ in the transverse dimensions to emulate experimental imperfections. We also consider toroidal traps with potential $V(x,y,z) = V_0 \{1-\exp( -\sigma^{-2}[\rho(x,y,z)-R]^2 ) \}$ , where V₀ is the trap depth, ρ is the radial distance from the torus centre, and σ and R are the minor and major radii, respectively.

In our model training we use quasi-two-dimensional (2D) condensates with an assumed Gaussian profile in the strongly confined z-dimension, except in section 4.2, where we test the model on column-integrated three-dimensional (but still highly anisotropic) condensates. The assumed tight axial confinement is equivalent to requiring $\hbar\omega_z \gg (\mu, k_{\mathrm{B}} T)$ and $\hbar\omega_x, \hbar\omega_y \ll \hbar\omega_z$ . The z-dimensional contribution is the harmonic oscillator ground state $\phi(z) = (\pi{\hbar}/{m\omega_{z}})^{-1/4}\exp(-m\omega_{z} z^2/2\hbar)$ , yielding the complete field $\Phi(x,y,z,t) = \Phi_\mathrm{2D}(x,y,t)\phi(z)$ . We can then evolve $\Phi_\mathrm{2D}(x,y,t)$ with a 2D equivalent to (6), where g and µ are replaced by the effective 2D parameters $g_\mathrm{2D} = g \sqrt{m\omega_{z}/2\pi\hbar}$ and $\mu_\mathrm{2D} = \mu - \hbar\omega_z/2$ .

3.1. Overview

We train our model on a set of M = 3000 atomic density profiles, $\rho(\mathbf{x}_\mathrm{2D})$ [where $\mathbf{x}_\mathrm{2D}\equiv(x,y)$ ], at thermal equilibrium, each of which is a square image of N rows of N pixels, where N = 256. The network architecture consists of an image processing and feature extraction pipeline followed by a prediction pipeline, as shown in figure 1. This architecture effectively processes the rich spatial information in atomic density profiles—extracting features at various orders of magnitude of the atomic density—to determine a pair of scalar values which we associate with the chemical potential and temperature.

Figure 1. Refer to the following caption and surrounding text. — **Figure 1.** Architecture of the convolutional neural network designed to take the density profile of an atomic Bose–Einstein condensate as an input, and from it determine the condensate’s chemical potential, µ, and temperature, T. The network processes 2D density profiles through three convolutional layers (extracting 12, 24, and 48 features, respectively), each followed by a maximum pooling layer with stride 2 (see appendix D for details). We apply global pooling to each feature before passing through two fully connected layers (FC1 and FC2), and use ReLU activation functions throughout.
Download figure:
Standard image High-resolution image

The image processing and feature extraction pipeline consists of three convolutional layers that increase the number of features (layer 1: 12, layer 2: 24, layer 3: 48) to be extracted in each layer, enabling the model to learn from different aspects of the input data (such as different orders of magnitude of the atomic density or any shape or curvature in the condensate). The reduction of spatial dimensions is not only for computational convenience⁴—by reducing the dimensionality of the input, we force the network to learn local, spatially coherent features and reduce the risk of overfitting (by using fewer parameters—referred to as weights and biases—in the network) [32].

The prediction pipeline transforms the spatial features extracted by the convolutional layers into scalar-valued estimations or ‘predictions’ of the chemical potential and temperature. By averaging each feature map in Layer 4, we reduce the matrix dimensions and capture global information from each feature. We then expand the model’s representational capacity from 48 to 512 neurons in the first fully connected layer (Layer 5 of the overall network)—this layer facilitates complex, highly non-linear combinations of the pooled features, which is essential for mapping the subtle details of the density profiles to the desired thermodynamic parameters. A large number of neurons in this layer improves the model’s predictive and generalisation capabilities—having fewer neurons in this layer was not as conducive to learning (and sufficiently few meant that the network could not determine the values of µ and T within any acceptable error) and more neurons did not appreciably improve our model’s predictive capabilities. There was no further processing in the output layer (e.g. an activation function); an in-principle consequence is the possibility of predicted negative values for the chemical potential (which can be meaningful, although not in the systems we consider) and the temperature (which is not meaningful). We never observed this to happen, and note that the validity of the SGPE model effectively assumes a finite temperature which is in some sense appreciable.

3.2. Scaled form of the stochastic SGPE

When generating the training data for the model, we consider a scaled system of units in order to avoid computing very small terms (e.g. of the order of $\hbar^2$ expressed in SI units). In particular, we introduce the following scaling factors: lengths are scaled by the harmonic oscillator lengths $\ell_{x,y,z} = \sqrt{\hbar}/{m\omega_{x,y,z}}$ , times by $1/\omega_x$ , energies by $\hbar\omega_{x}$ , and temperatures by $\hbar\omega_{x}/k_{\mathrm{B}}$ . The condensate field is rescaled as $\tilde{\Phi}(\mathbf{x}_\mathrm{2D}, t) = \ell_x \Phi(\mathbf{x}_\mathrm{2D}, t)$ , and the position and time variables become $\tilde{\mathbf{x}}_\mathrm{2D} = \mathbf{x}_\mathrm{2D} / \ell_x$ and $\tilde{t} = t / \tau$ respectively, where $\tau = 1/\omega_x$ . The Laplacian transforms as $\tilde{\nabla}_\mathrm{2D}^2 = \ell_x^2\nabla^2_\mathrm{2D}$ , while the effective 2D chemical potential and interaction strength are rescaled as $\tilde{\mu}_\mathrm{2D} = \mu_\mathrm{2D} / \hbar \omega_x$ and $\tilde{g}_\mathrm{2D} = g_\mathrm{2D} / (\hbar \omega_x \ell_x^2)$ . Neglecting the tildes for convenience, the resulting form of the SGPE is:

Equation (10)

with equivalently scaled Gaussian noise ensemble correlations

Equation (11)

It is these forms that we use to generate all atomic density samples in this work.

In our simulations, we consider a dimensionless spatial step size $\Delta x = 0.142\,62$ , which for the ⁸⁷Rb system and the $256\times 256$ grid we consider, corresponds to a total grid length of 80 µm. We evolve the system by a maximum of 100 000 time steps (although this is typically much smaller depending upon the dynamical time to thermalisation) using a dimensionless time step size $\Delta t = 10^{-3}$ . Other simulation parameters (in scaled and unscaled forms) can be found in our source code.

In the following sections, we will first describe the prediction pipeline, starting at layer 4 of the overall neural network, followed by a detailed explanation of the image processing and feature extraction pipeline. We outline the specific mathematical operations used in our convolutional network and contextualise what is often thought of as a ‘black box,’ in a presentation intended to be relevant and accessible for anyone with a physics background.

3.3. Prediction pipeline

The prediction pipeline is a fully connected, feedforward neural network, consisting of three layers, as shown in figure 2. This is preceded by the feature extraction pipeline, which reduces the dimensionality of the input images, through a sequence of three convolutional layers. The output of the feature extraction pipeline consists of 48 features $\Xi^3_i$ , which are $32\times 32$ square matrices, i.e. $\Xi^3_i \in \mathbb{R}^{N/8 \times N/8}$ , where $i \in \{1, \ldots, 48\}$ . Taking the average values over the rows and columns of the 48 features $\Xi^3_i$ produces 48 reduced feature representations $y_{i}^{4}$ . These are the input values for the first layer of the prediction pipeline, and the fourth layer overall.

Figure 2. Refer to the following caption and surrounding text. — **Figure 2.** The prediction pipeline architecture, showing in more detail the final part of the network architecture depicted in figure 1. Three fully connected, feedforward layers convert the spatial information from the feature extraction pipeline to a pair of scalar values associated with the chemical potential and temperature, highlighting the weights $w_{j,k}^{\ell}$ between layers and the biases $b_{j}^{\ell}$ associated with the neurons in a layer, for $\ell = 5$ .
Download figure:
Standard image High-resolution image

**Figure 2.** The prediction pipeline architecture, showing in more detail the final part of the network architecture depicted in figure 1. Three fully connected, feedforward layers convert the spatial information from the feature extraction pipeline to a pair of scalar values associated with the chemical potential and temperature, highlighting the weights $w_{j,k}^{\ell}$ between layers and the biases $b_{j}^{\ell}$ associated with the neurons in a layer, for $\ell = 5$ .
Download figure:
Standard image High-resolution image

In Layer 5, there are 512 neurons with ReLU activation (a necessarily nonlinear function that is commonly used in image recognition scenarios [33]). The pre-activation values are determined from the $y_{k}^{4}$ through

Equation (12)

where we define $w^\ell_{jk} \in \mathbb{R}$ as the weight connecting the kth neuron in layer $\ell-1$ to the jth neuron in layer $\ell$ , and $b^\ell_j \in \mathbb{R}$ as the bias for the jth neuron in layer $\ell$ (we use the same notation detailed in appendix A of [34]). We then apply the activation function to produce

Equation (13)

These values feed into Layer 6 (the final layer), through

Equation (14)

which gives estimations for the chemical potential $\mu = y^6_1$ and temperature $T = y^6_2$ , in nanokelvin. We initialise the weights and biases using the Kaiming scheme [35], as $w^5_{jk},\ b^5_j \sim U(-1/\sqrt{512}, 1/\sqrt{512})$ and $w^6_{jk},\ b^6_j \sim U(-1/\sqrt{2},1/\sqrt{2})$ (see appendix C in [34] for details). The final layer does not use an activation function, which is motivated by the physics. The chemical potential can be any real value, and the temperature can be any positive, real value, so no activation function is required (one might, in principle, use the $\mathrm{ReLU} = \max(0,x)$ activation for the temperature, but we did not observe negative temperatures being output from our model)⁵.

3.4. Image processing and feature extraction pipeline

The input atomic density can be thought of as a monochromatic image described by a matrix $\rho \in \mathbb{R}^{N \times N}$ . We show an example in figure 3(a)—while we display it using a false colour scale, in terms of information (each pixel has a single value assigned to it, describing the local density) it is functionally monochromatic. The image is processed through three convolutional layers, each followed by ReLU activation and maximum pooling (MaxPool). Figure 3 demonstrates intermediate steps to calculate a single feature from the first convolutional layer.

Figure 3. Refer to the following caption and surrounding text. — **Figure 3.** The first convolutional layer **Conv1**, elaborating in more detail the beginning of the image processing and feature extraction pipeline as depicted in figure 1. (a) Input of the atomic density, ρ. (b) Cross-correlation of the atomic density with the weights matrices (see equation (15)). (c) Application of the ReLU activation function and a bias (see equation (16)). (d) Carrying out maximum pooling (resulting in the halving of the image dimensions); this is one of j = 12 outputs of the first convolutional layer, as per equation (17). To the right we show zoomed-in sections of each step of the first convolutional layer as we construct the first layer feature maps.
Download figure:
Standard image High-resolution image

Layer 1: from a single input to 12 features

In the first layer ( $\ell = 1$ ), we extract C₁ = 12 features from the single input image ρ (such that the number of features in the ‘zeroth layer’ is C₀ = 1). The weights for this layer consist of C₁ = 12 filters connecting the single input channel (k = 1) to each output channel j, denoted as $\mathsf{W}^{\ell = 1}_{j,k = 1}(\sigma, \tau)$ for $j = 1, \dots, C_1$ . Each filter (in this layer and in all subsequent layers) is a square matrix of dimension $3 \times 3$ , and we initialise each weight element within these filters, independently, from a uniform distribution $\mathsf{W}^1_{j1}(\sigma, \tau) \sim U(-1/3, 1/3)$ (the range of the probability distribution function depends on the number of features in the previous layer—see appendix B).

We carry out cross-correlations (see appendix A for details) between the input image ρ and each filter $\mathsf{W}^1_{j1}$ to compute the intermediate feature maps $\zeta^1_j$ :

Equation (15)

Here, $s_1, t_1 \in \{1, 2, \dots, 256\}$ (assuming appropriate padding, see appendix C), and we have displayed the cross-correlation of an example atomic density with the final (learnt) weights matrices in figure 3 (b)) (see figure F1 in appendix F for the outputs of the cross-correlations of the same atomic density with all 12 weights matrices).

We then add a feature-map-specific bias term $b^1_j$ (these are always drawn from the same distribution, appropriate to the layer, as the weights) to each matrix element of the intermediate feature map, and subsequently also apply the ReLU activation function to each element:

Equation (16)

We show the activation function and bias applied to a single output of the cross-correlation in figure 3 c) (figure F2 in appendix F shows the results $\xi_{j}^{1}$ of adding the bias and applying ReLU activation to each of the intermediate feature maps $\zeta_{j}^{1}$ ).

Finally, we apply MaxPool with a $2 \times 2$ window and stride 2 (see appendix D for details) to obtain the first layer’s output feature maps $\Xi^1_j$ :

Equation (17)

This halves the matrix dimensions, so that $s_2, t_2 \in \{1, 2, \dots, 128\}$ . These C₁ = 12 feature maps $\Xi^1_j$ , each of size $128 \times 128$ , become the input to the next layer. We show a single, maximally pooled output in figure 3(d) (figure F3 in appendix F shows the full set of feature maps from the first layer for the particular sample density).

Layer 2: from 12 features to 24 features

In the second layer, we extract C₂ = 24 features from the C₁ = 12 input feature maps. The weights for this layer consist of $C_2 \times C_1 = 24 \times 12 = 288$ filters, denoted as $\mathsf{W}^2_{jk}(\sigma, \tau)$ for $j = 1, \dots, 24$ (output channel index) and $k = 1, \dots, 12$ (input channel index). We initialise each weight element within these filters independently, $\mathsf{W}^2_{jk}(\sigma, \tau) \sim U(-1/\sqrt{108}, 1/\sqrt{108})$ .

We compute the intermediate feature maps $\zeta^2_j$ by summing the cross-correlations over all input channels:

Equation (18)

See figure F4 in appendix F for displays of each output of these cross-correlations from the same sample input. In completely analogous fashion we again add bias terms $b^2_j$ , followed by application of the ReLU activation function:

Equation (19)

Figure F5 in appendix F equivalently shows the subsequent sample results for all values of j. Finally, we apply MaxPool with stride 2 to obtain the second layer’s output feature maps $\Xi^2_j$ :

Equation (20)

The matrix dimensions again halve, so $s_3, t_3 \in \{1, 2, \dots, 64\}$ . These C₂ = 24 feature maps $\Xi^2_j$ , each of size $64 \times 64$ , become the input to the next layer; figure F6 in appendix F displays each of these feature maps resulting from the same sample input density.

Layer 3: from 24 features to 48 features

In the third layer, we extract the final C₃ = 48 features from the C₂ = 24 input feature maps. These final features feed into the subsequent prediction pipeline. The weights for this layer consist of $C_3 \times C_2 = 48 \times 24 = 1152$ filters, denoted as $\mathsf{W}^3_{jk}(\sigma, \tau)$ for $j = 1, \dots, 48$ and $k = 1, \dots, 24$ . We initialise each weight element within these filters independently $\mathsf{W}^3_{jk}(\sigma, \tau) \sim U(-1/\sqrt{216}, 1/\sqrt{216})$ .

We carry out cross-correlations, add bias terms, and apply the ReLU activation function and MaxPool (with stride 2) in exactly the same way as in the second layer. Hence,

Equation (21)

Equation (22)

Equation (23)

Here, $s_4, t_4 \in \{1, 2, \dots, 32\}$ . These C₃ = 48 feature maps $\Xi^3_j$ , each of size $32 \times 32$ , are the final feature maps from the image processing and feature extraction pipeline that feed into the prediction pipeline, and can be seen in figure F9 in appendix F, with figures F7 and F8 showing the corresponding intermediate steps, for the same sample input density as for each of the previous figures in the appendix.

3.5. Computational considerations and learning

3.5.1. Cost function

Training a neural network involves optimising its parameters (the weights and biases) to minimise a cost function. For our model, which determines chemical potential and temperature, we consider the square error

Equation (24)

When determining $\mathcal{C}$ , we express both the temperature and chemical potential in nanokelvin; with regard to $\mu^\mathrm{predicted}$ and $\mu^\mathrm{true}$ , these values are strictly speaking chemical potentials divided by $k_{\mathrm{B}}$ , and then expressed in nanokelvin. The values making up $\mathcal{C}$ are thus comparable and of moderate magnitude.

3.5.2. Batching

For computational efficiency, we train the network using batches of atomic densities, each with independent associated temperature and chemical potential values, rather than processing them individually. We process input batches of size β (typically 16, 32, 64 or 128) as multidimensional arrays⁶, which we form by adding two additional channels, β and $C_\ell$ , to the matrix at a given layer $\ell$ . The cost function with batching is

Equation (25)

Weight matrices and bias terms in each convolutional layer are shared across all items in a batch, such that the same learned filters are applied identically to each input sample. This architectural choice does not reduce the number of parameters relative to processing individual inputs, but it avoids parameter duplication across the batch, thereby preserving model compactness and promoting generalisation. Moreover, while the data are different, because the operations applied to each sample are mathematically identical, the use of batches enables these computations to be vectorised and executed in parallel on GPU hardware. This leads to improved computational efficiency without altering the underlying functional form of the network. In practice, this parallelism is realised by arranging input data as multidimensional arrays with an added batch dimension, facilitating simultaneous convolution, activation, and pooling operations across multiple samples.

3.5.3. Training

To minimise the cost function in equation (25), we use the adaptive moment estimation (Adam) algorithm [37]—a widely used and generally effective adaptive learning rate method suitable for training large networks; see appendix D, section 4 of [34] for detailed derivations—using the default values for the learning rate, η = 0.001, and decay rates, $(\rho_1, \rho_2) = (0.9, 0.999)$ , and backpropagation for gradient calculation. We iteratively use Adam over 100 epochs—where each epoch processes the entire dataset in randomly selected batches (without replacement)—although we note a posteriori that using fewer epochs is possible since the cost saturates at around epoch 50, as shown in figure 5. After training, we save the learned weights and biases. Post-training, we can load these optimised parameters into a network with the same architecture to make rapid inferences from new data. We observe inference times for our model of the order of milliseconds in wall-clock time.

As described in section 3.4, in each layer the determination of the ζ matrices is associated with the introduction of the weights, and the determination of the ξ matrices with the introduction of the biases. Determination of the Ξ matrices is not part of the learning procedure, but processes the images for the next layer by halving the image dimensionality through MaxPool.

4.1. Interpretation of the feature maps

In figure 4 we show schematically the steps described mathematically in section 3.4 to produce from an initial density profile the feature maps from each of the three convolutional layers in the image processing and feature extraction pipeline. This displays subsets of the 12 feature maps produced in layer 1, the 24 feature maps produced in layer 2, and the 48 feature maps produced in layer 3, where plots, with colour axes, for the complete set of pre-activation images (the ζ matrices), post activation images (the ξ matrices), and feature maps (the Ξ matrices), for the same run, are displayed in appendix F.

Figure 4. Refer to the following caption and surrounding text. — **Figure 4.** Schematic showing the production of the final 48 feature maps from an initial density profile, via the *image processing and feature extraction pipeline*, which then feed into the *prediction pipeline* (see also figure 2) to determine estimated values of the chemical potential µ and temperature T. The initial density profile is the same as that of figure 3(a), the (1 of 12) pre-activation image is the same as in figure 3(b), the following post-activation image is the same as in figure 3(c), and the first of the following sample from 12 feature maps is the same as in figure 3(d). Relative to figure 1, this schematic depicts detailed progress, a single channel of each layer at a time, through the image processing and feature extraction pipeline (see section 3.5.2 for details). The end-of-layer outputs from all channels (individual ‘slices’ in figure 1) are inputs to the subsequent layer; the 48 feature maps output by layer 3 are individually globally maximally pooled, as per section 3.3. The complete set of pre-activation images, post-activation images, and feature maps, for each layer, can be seen in in appendix F.
Download figure:
Standard image High-resolution image

Increasingly abstract feature maps are learnt as we progress through the three layers. These features may include local variations in the density, patterns, or edges corresponding to the trapping potential, recognisable structures such as vortices, or other structures in the atomic cloud. The features may be local or global, with the convolutional layers trying to identify the fingerprints of the chemical potential and temperature values in the spatial structure of our input samples. The first set of feature maps, produced in the first convolutional layer (a subset of which are shown in figure 4), learn high-density features which we can reasonably infer to be associated with the chemical potential, since the chemical potential is fundamentally connected to the average atom number. For the complete set of pre-activation images, post activation images, and feature maps see appendix F, figures F1–F3. The third and final set of feature maps learn low-density features, which may be associated with thermal fluctuations at this scale. The role of temperature in the atomic density is only introduced in the thermal noise, given by the correlation function in equation (8). Accurate prediction of the temperature using our model is dependent upon being able to effectively capture thermal fluctuations, which the convolutional network is attempting to observe in this final convolutional layer. For the complete set of pre-activation images, post activation images, and feature maps see appendix F, figures F7–F9. Between the first and final layers, the second set of feature maps learn a mixture of large-scale features (associated with µ) and small-scale features (associated with T); for the complete set of pre-activation images, post activation images, and feature maps see appendix F, figures F4–F6. This hierarchical feature extraction is essential for accurately predicting both the chemical potential and temperature.

Recognising the distinct roles that the chemical potential and temperature have in our feature maps, we design a series of interrogations to ascertain how our model could be used in a range of experimentally relevant protocols. In particular, we ask: 1) does the predictive capability of the model depend on whether the density profile being observed is at or very close to thermal equilibrium (which strictly speaking is required for the chemical potential and temperature to be properly defined), e.g. could we predict the chemical potential and temperature during the thermalisation of a condensate? 2) If the model does not have a strong dependence on the equilibrium distribution, could we use the model on toroidally trapped condensates? 3) Can our model predict the chemical potential and temperature of an ensemble average of equilibrium states? We will address each of these questions in turn, immediately after presenting our choice of model, which we do in the next section.

4.2. Model accuracy

We considered a balance of computational cost (models extracting more features per layer have higher computational cost) and model accuracy (using a metric which we refer to as the accuracy within 5%, as defined below) to determine which model to use in our subsequent analysis. We considered models with batch sizes of $\beta = 16, 32, 64,$ or 128 (see section 3.5.2 for details), which extract either 6, 12 and 24 features, 12, 24 and 48 features or 16, 32 and 64 features in each convolutional layer.

Figure 5 shows the training and validation cost metrics for the models that extract 12, 24 and 48 features in the convolutional layers over a range of batch sizes. The absence of divergence between the training and validation cost metrics is often used as a heuristic to assess whether overfitting is occurring; in this case, it provides no such indication. Note that this figure alone is insufficient to ascertain the accuracy of the models. We observe no appreciable divergence in the training or validation cost metrics for all models (with fewer or more features)—all models approached final costs of the order of 10⁻¹, despite their accuracy varying significantly.

Figure 5. Refer to the following caption and surrounding text. — **Figure 5.** The training and validation cost metrics for batch sizes 16, 32, 64, and 128 from a model extracting 12, 24 and 48 features in the first, second, and third convolutional layers. A smaller batch size results in models with a lower overall training and validation cost.
Download figure:
Standard image High-resolution image

We define ‘accuracy within 5%’ as the proportion of samples with predictions within 5% of their true values—that is, the values which we set as input parameters to the stochastic evolution. Using this metric, figure 6 illustrates the predictive accuracy on the test dataset held back for this purpose, for models with 12, 24, and 48 convolutional features, across batch sizes $\beta = {16, 32, 64}$ and 128. We see that the standard deviation of absolute errors decreases with decreasing batch size, indicating more accurate predictions of both chemical potential and temperature. Smaller batch sizes generally improve generalisation [38] by introducing more stochasticity into the optimisation process, helping the model escape shallow local minima and find solutions that generalise better to unseen data, and this is what we appear to be observing here.

Figure 6. Refer to the following caption and surrounding text. — **Figure 6.** The histograms and corresponding normal distributions demonstrate the absolute error in the model’s prediction of the chemical potential (column 1) and temperature (column 2) of the held-out test data set to the true values input to our SGPE simulations. The plots share a common x-axis for easier comparison. These plots correspond to the models which extract 12, 24 and 48 features.
Download figure:
Standard image High-resolution image

As summarised in table 1, reducing the convolutional feature set to 6, 12, and 24 features has negligible effect on chemical-potential estimation but degrades temperature prediction, whereas the 12, 24, and 48 feature configuration gives the best balance of speed and accuracy. Models with 16, 32, and 64 features produced no measurable improvement over the model with 12, 24, and 48 features—accuracies unchanged to two decimal places. We have found that the model that extracts 12, 24, and 48 features, with a batch size of 16, offers a good balance between speed and accuracy. All analyses from this point on use this model.

Table 1. Test dataset accuracy for chemical potential µ and temperature T, defined as the proportion of predictions within 5% of the true values (see § 4.2). Results shown for three model groups (feature sets) across batch sizes β. The models with the best accuracy within 5% for the chemical potential and temperature are highlighted in bold.

Model (features)	Batch size β	Accuracy within 5% (µ)	Accuracy within 5% (T)
6, 12, 24	16	1.00	0.92
	32	0.99	0.92
	64	0.98	0.82
	128	0.98	0.73

12, 24, 48	16	1.00	0.96
	32	1.00	0.94
	64	0.99	0.90
	128	0.99	0.88

16, 32, 64	16	1.00	0.96
	32	1.00	0.94
	64	0.99	0.90
	128	0.99	0.88

We briefly note that column-integrated 3D atomic densities may also be used in our model with similar accuracy to the quasi-2D condensates. Given a three dimensional atomic density profile, $|\Phi(x,y,z)|^2$ , the column-integrated density (assuming the imaging is in the z-direction) is

Equation (26)

The column-integrated density may then be passed through the machine learning model as previously described with no further post-processing.

4.3. Prediction capability and robustness

4.3.1. Estimation capability during thermalisation

In our simulations the bath, with which there is formally both heat and particle exchange, is consistently parametrised throughout by a constant chemical potential, µ, and temperature, T. We do not consider dynamics following any form of quench, and the machine learning model is trained only on thermalised Bose gases. As described in equation (9), we use convergence of the atom number as a practical criterion for having achieved thermal equilibrium. However we note that estimated values for µ and T given by our model appear to trend to the final values more quickly, during the thermalisation process and after a relatively brief transient, as shown in figure 7. During the initial stages of evolution (approximately $t \lesssim 10\,\omega_x^{-1}$ , as shown in figure 7), the model significantly underestimates both µ and T. At these early times, the atomic density distribution has not yet developed the equilibrium spatial features that the convolutional layers have learned to recognise: specifically, the Thomas–Fermi-like bulk profile that encodes µ through the overall density scale and spatial extent, and the statistical properties of thermal fluctuations that encode T. However, once these characteristic structural features begin to emerge—which occurs whilst the condensate atom number is still evolving toward the equilibrium value defined by our convergence criterion (equation (9))—the model’s predictions converge rapidly to accurate values. This demonstrates that the model has learned to extract features indicative of the system’s eventual thermodynamic state before full equilibration is achieved, despite training exclusively on equilibrium density profiles.

Figure 7. Refer to the following caption and surrounding text. — **Figure 7.** The evolution of the estimates of the chemical potential and temperature over the thermalisation procedure on dual axes. The pink line is the absolute error in the chemical potential. The orange line is the absolute error in the temperature. The purple line is the atom number. Top: evolution of a harmonically trapped condensate. The temperature of the same is 30 nK. The chemical potential of the sample is also 30 nK. The absolute errors in the temperature and chemical potential plateau after 80 units of rescaled time.
Download figure:
Standard image High-resolution image

4.3.2. Toroidally trapped condensates

Toroidally trapped condensates permit metastable persistent currents (quantised flow of Bose–Einstein condensates in a multiply connected geometry) [39] and provide an appropriate geometry for experimental protocols such as weak links [40, 41] which are used in several atomtronic devices such as rotational sensors [10]. As an additional probe of the robustness of our model, which has been trained exclusively with harmonic trapping configurations, we apply it to toroidally trapped condensates. These have depth $V_0 = 60\mathrm{\,nK}$ , minor radius $\sigma = 20\,\mu$ m and major radius $R = 40\,\mu$ m; we choose V₀ such that $V_0/\mu \gt 1$ .

While, as shown in figure 8, the estimated µ and T trend towards stationary values relatively quickly, the temperature is generally closer to the input value than the chemical potential. This difference reflects that µ is encoded in the bulk density distribution, which differs qualitatively between harmonic and toroidal traps, whereas T is encoded in small-scale fluctuations that are more geometry-invariant.

Figure 8. Refer to the following caption and surrounding text. — **Figure 8.** The evolution of the estimates of the chemical potential and temperature over the thermalisation procedure on dual axes. The pink line is the absolute error in the chemical potential. The orange line is the absolute error in the temperature. The purple line is the atom number. Top: evolution of a harmonically trapped condensate. The temperature of the same is 22 nK. The chemical potential of the sample is also 22 nK. The absolute errors in the temperature and chemical potential plateau after 80 units of rescaled time. The depth of the trap is 60 nK. The minor and major radii are, respectively, 20 µm and 40 µm.
Download figure:
Standard image High-resolution image

Convolutional neural networks exploit weight sharing—each $3\times 3$ filter $\mathsf{W}_{jk}^\ell$ is applied identically across all spatial locations, so the same kernel contributes to every receptive field (i.e. each local $3\times 3$ patch of the input). During backpropagation, the gradient of a filter weight is accumulated over all receptive fields, so updating one weight modifies the response of that filter simultaneously at all positions. This enforces a translation-equivariant inductive bias [32], enabling the network to learn recurring local patterns independently of their absolute position. We interpret that convolutional filters scanning local receptive fields extract temperature-linked signatures that remain invariant under a change from harmonic to toroidal confinement.

Two additional factors accentuate the asymmetry: (i) the random transverse anisotropy injected during training (section 2.3) encourages robustness of local features to geometric perturbations, and (ii) the maximally global pooling before the fully connected layers discards absolute position, emphasising distributional summaries of fluctuation statistics rather than any geometry-specific structure.

4.3.3. Ensemble averages

Both average pooling and MaxPool are commonly used approaches in, e.g. image recognition. As discussed in appendix D, average pooling gives a general representation of smaller regions in the image whereas MaxPool emphasises the most prominent features in those regions. It is also possible to think of an equivalency between average pooling and taking an ensemble average over several noise trajectories, where individual trajectories are stochastic; change associated with averaging over these trajectories is smooth and deterministic. Akin to experiments, however, numerical runs obtained from a single noise realisation hold important physical information; one may interpret a single numerical run of the SGPE as an independent experimental realisation [24, 42–44].

As shown in figure 9, estimation of the chemical potential is quite comparable when using either maximum or average pooling, however when using average pooling information about the temperature is clearly lost. This is consistent with our observation, when considering toroidally trapped condensates, that estimation of the chemical potential is associated with the bulk distribution of the condensate (essentially similar regardless of the pooling type), whereas estimation of the temperature is associated with small-scale local density fluctuations, which will be averaged away by average pooling.

Figure 9. Refer to the following caption and surrounding text. — **Figure 9.** A comparison of the predictive capabilities of a network trained with maximum pooling and a network trained with average pooling. Both maximum pooling and average pooling are appropriate for learning and determining the chemical potential of Bose gases, but only maximum pooling is appropriate for learning and determining their temperature.
Download figure:
Standard image High-resolution image

Similarly, as we have tested, averaging over many trajectories to produce an averaged density washes out the temperature-dependent background fluctuations, and we observe that ensemble averages passed through our machine learning model can only be used to obtain the chemical potential and not the temperature of the sample.

5.1. Effective pixel size

To establish an effective pixel size in a typical experiment, we consider the imaging system of Wilson, et al [45], which—like our model—considers in situ imaging of condensates to obtain position distributions. Wilson, et al use a Point Grey Firefly MV CMOS camera with pixels of size $6.0\,\mu \times 6.0\,\mu$ m with a magnification of around 19.7, which leads to an effective pixel size of 6.0 µm/19.7 = 0.31 µm per pixel. We align our model’s theoretical maximum resolution to this effective pixel size based on this camera. For a quasi-2D ‘pancake’ condensate of diameter 80 µm in a harmonic trap, we need to use a grid size of $256\times 256$ in order to achieve an effective pixel size of about 0.31 µm.

5.2. Further considerations

Whilst our model’s resolution is appropriate for experimental data analysis, real BEC images often include additional complexities such as imaging artefacts, focus variations, and dark counts. In the work of Wilson et al, the optical resolution of the microscope (with a numerical aperture of 0.25) is actually diffraction-limited to ${\approx}1.9\,\mu$ m when illuminated with light of wavelength 780 nm. More recent advances in in situ imaging of BECs by [46] and [47] report diffraction-limited resolutions of 0.7 µm and 0.9 µm, respectively. In practice this means that experimental images are effectively our simulated atomic density profiles convolved with the point-spread function of width given by the diffraction-limited resolution of the imaging system. Taking this known complexity into account—alongside augmentation with the other aforementioned experimental artefacts—would enhance the model’s robustness for real data. We also note that, while our model has been trained on in situ imaged position distributions, training on expanding time-of-flight velocity distributions is an in-principle straightforward extension, although significantly more demanding to generate if the full dynamics are included.

We also draw attention to very recent independent work, developed in parallel by De Sousa et al [48], for machine learning thermometry of a dilute, laser cooled (non-quantum-degenerate) vapour of potassium 39 atoms, from experimental data.

We have introduced a proof-of-principle machine learning model which can, with a single image, estimate the values of the chemical potential and temperature of a Bose-condensed cloud of atoms. We use convolutional neural networks—the foundation of most image recognition models—to take an atomic density profile and estimate important thermodynamic parameters. We have demonstrated that our model for a harmonically trapped condensate can estimate to some degree the chemical potential and temperature for systems that it has not previously been trained on, including during thermalisation, and on toroidally trapped condensates. This work joins recent applications of machine learning in the quantum fluids literature, such as the identification of topological defects such as vortices and solitons, and the reconstruction of vortex filaments [19–21].

We thank K E Wilson for insightful discussions about the experimental considerations of this work. We thank the EPSRC, Grant No. EP/T015241/1 and EP/T518001/1, for funding this project.

The data that support the findings of this study are openly available at the following URL/DOI: https://github.com/0jg/bec-ml-thermometry [31] and https://github.com/0jg/sgpe-rs [49]

In section 3.4, we introduced the cross-correlation function. We define the cross-correlation function (between y and w) as

Equation (A.1)

Relative to the more generally known convolution operation, the cross-correlation operation is effectively a clockwise rotation of the kernel by 180^∘. We extend equation (A.1) to matrices by introducing another index, σ,

Equation (A.2)

The cross-correlation operation may be interpreted as a Frobenius inner product of a sub-matrix of $\mathsf{Y}$ and the kernel, $\mathsf{W}$ [50].

For historical and conventional reasons, it is the cross-correlation that is calculated in machine learning libraries such as PyTorch [51] and Tensorflow [52]—‘convolutional’ neural networks are a technical misnomer, but during training, convolutions and cross-correlations should converge on the same result.

Since most of the data we are interested in are 2D matrices, it is equation (A.2) that we use throughout this paper.

In a convolutional layer $\ell$ (where $\ell \in \{1, 2, 3\}$ ), we aim to extract $C_\ell$ features (output channels) from the $C_{\ell-1}$ features (input channels) provided by the previous layer (layer $\ell-1$ ). For the first layer ( $\ell = 1$ ), the input is the single map ρ, so C₀ = 1.

To compute the jth output feature map (where $j \in \{1, \dots, C_\ell\}$ ), we use a set of weights (filters). Specifically, for each input feature map k (where $k \in \{1, \dots, C_{\ell-1}\}$ ), there is a 2D filter $\mathsf{W}^\ell_{jk}$ of spatial size $k_W \times k_W$ . For $\ell = 1$ , k can only be 1, so the filters are $\mathsf{W}^1_{j1}$ . The collection of $C_{\ell-1}$ filters $\{\mathsf{W}^\ell_{j1}, \dots, \mathsf{W}^\ell_{j,C_{\ell-1}}\}$ is used together (summing their individual cross-correlation results with corresponding input maps) to produce the jth intermediate output map $\zeta^\ell_j$ .

We initialise each individual weight element $\mathsf{W}^\ell_{jk}(\sigma, \tau)$ (where $\sigma, \tau$ are spatial indices within the $k_W \times k_W$ filter) independently from the uniform distribution [51]:

Equation (B.1)

Here, $C_{\ell-1}$ is the number of input channels (‘fan-in’ channels) to the layer, and $k_W^2$ is the spatial area of the filter. The term $C_{\ell-1} k_W^2$ represents the total number of inputs that contribute to a single element in the output map before activation (the fan-in). We sample the biases $b_j^\ell$ from the same distribution U.

The total number of distinct weight matrices (filters), N_filters, in a convolutional neural network with L layers is

Equation (B.2)

The machine learning problem is to find the appropriate weights $\mathsf{W}^\ell_{jk}$ and biases $b^\ell_j$ which optimise the mapping from the input data to the desired output.

To preserve the spatial dimensions of the feature maps during the cross-correlation operation (before pooling), we apply zero padding around the input feature maps of layer $\ell$ . For a filter (kernel) of size $k_W \times k_W$ , where k_W is odd (e.g. $k_W = 3$ as used in the main text), we add κ zeros to each side (top, bottom, left, right) of the input maps, where

Equation (C.1)

For $k_W = 3$ , κ = 1. This ‘same’ padding ensures that the output of the cross-correlation, $\zeta^\ell_j$ , has the same matrix dimensions as the input feature maps, $\Xi^{\ell-1}_k$ , before pooling is applied.

Pooling operations reduce the spatial dimensions (width and height) of feature maps. This provides a degree of translation invariance and reduces the computational cost in subsequent layers. Pooling operates independently on each feature map. A pooling kernel of size $k_\mathrm{P} \times k_\mathrm{P}$ slides across the input map with a stride Σ. Common types are average pooling (AvgPool), which computes the mean of the values within the kernel window, and MaxPool, which takes the maximum value.

Figure D1. Refer to the following caption and surrounding text. — **Figure D1.** Illustration of pooling using a kernel of size $k_\mathrm{P}\times k_\mathrm{P} = 2\times 2$ and stride length $\Sigma = 2$ on a $4\times 4$ grid. The output is a square matrix with reduced width and height $\lfloor (N_\mathrm{in}-k_\mathrm{P})/\Sigma \rfloor + 1 = \lfloor (4-2)/2 \rfloor + 1 = 2$ . In the main text, we use MaxPool with $k_\mathrm{P} = 2$ and $\Sigma = 2$ . The convolutional filters have size $k_W \times k_W = 3 \times 3$ .
Download figure:
Standard image High-resolution image

**Figure D1.** Illustration of pooling using a kernel of size $k_\mathrm{P}\times k_\mathrm{P} = 2\times 2$ and stride length $\Sigma = 2$ on a $4\times 4$ grid. The output is a square matrix with reduced width and height $\lfloor (N_\mathrm{in}-k_\mathrm{P})/\Sigma \rfloor + 1 = \lfloor (4-2)/2 \rfloor + 1 = 2$ . In the main text, we use MaxPool with $k_\mathrm{P} = 2$ and $\Sigma = 2$ . The convolutional filters have size $k_W \times k_W = 3 \times 3$ .
Download figure:
Standard image High-resolution image

In this paper, we apply MaxPool after the activation function in each convolutional layer ( $\ell \in \{1, 2, 3\}$ ). The inputs to the pooling operation are the activated feature maps $\xi^\ell_j$ . We use a pooling kernel size $k_\mathrm{P} = 2$ and a stride $\Sigma = 2$ . If an input map has matrix dimensions $N_\mathrm{in} \times N_\mathrm{in}$ , the output map dimensions $N_\mathrm{out} \times N_\mathrm{out}$ are given by:

Equation (D.1)

For $k_\mathrm{P} = 2, \Sigma = 2$ , this simplifies to $N_\mathrm{out} = N_\mathrm{in}/2$ .

The output feature maps after pooling in layer $\ell$ , denoted $\Xi^\ell_j$ , are computed from the activated maps $\xi^\ell_j$ . Using $s^{{\prime}}, t^{{\prime}}$ for the output map, $s, t$ for the input map $\xi^\ell_j$ ), the MaxPool operation is:

Equation (D.2)

With $k_\mathrm{P} = 2$ and $\Sigma = 2$ , this matches the formulae used in the main text, e.g. equation (17).

The average pooling equivalent (not used in the main text) would be:

Equation (D.3)

For conceptual comparison, we can draw an analogy between the operations in a convolutional layer (without pooling) and a fully connected neural network layer. Consider a simplified 1D case, where the output $y_j^\ell$ of neuron j in layer $\ell$ of a fully connected neural network is typically computed as:

Equation (E.1a)

Equation (E.1b)

Here, $\mathcal{A}$ is the activation function, $z_j^\ell$ is the pre-activation value, $y_k^{\ell-1}$ are the activations from the previous layer ( $\ell-1$ ), $w_{jk}^\ell$ is the weight connecting input neuron k to output neuron j, $b_j^\ell$ is the bias for neuron j, and $n_{\ell-1}$ is the number of neurons in layer $\ell-1$ .

In the convolution neural network context (equations (18) and (19)), ignoring pooling and 2D structure for simplicity), the computation for the jth feature map involves:

Equation (E.2a)

Equation (E.2b)

The key difference lies in equation (E.2b) using a local cross-correlation ( $\star$ ) operation with shared weights within the filter $\mathsf{W}^\ell_{jk}$ , rather than a simple weighted sum (dot product) over the entire previous layer output as in equation (E.1b). However, both involve a form of weighted summation over inputs from the previous layer (summing over index k up to the number of input features/neurons, $C_{\ell-1}$ or $n_{\ell-1}$ ), followed by a bias and activation function, highlighting the conceptual link.

The figures in this appendix show the complete feature extraction process for the same atomic density profile used in figure 3. Each set of three figures corresponds to one convolutional layer: pre-activation maps (cross-correlation outputs), post-activation maps (after ReLU and bias), and final feature maps (after MaxPool).

Figures F1–F3 show the first layer outputs (12 features), which primarily capture high-density regions and large-scale structural information. Figures F4–F6 show the second layer outputs (24 features), which learn intermediate-scale patterns. Figures F7–F9 show the third layer outputs (48 features), which are sensitive to fine-scale fluctuations that encode temperature information.

Figure F1. Refer to the following caption and surrounding text. — **Figure F1.** The cross-correlation, ζ¹, of a single image of an atomic density, $\rho(x,y)$ , with 12 weights kernels, $\mathsf{W}^1_i$ , $i\in\{1,\cdots,12\}$ , as determined by equation (15). The cross-correlation may result in positive or negative values, as indicated by the colour bars.
Download figure:
Standard image High-resolution image

**Figure F1.** The cross-correlation, ζ¹, of a single image of an atomic density, $\rho(x,y)$ , with 12 weights kernels, $\mathsf{W}^1_i$ , $i\in\{1,\cdots,12\}$ , as determined by equation (15). The cross-correlation may result in positive or negative values, as indicated by the colour bars.
Download figure:
Standard image High-resolution image

Figure F2. Refer to the following caption and surrounding text. — **Figure F2.** ξ¹, the result of applying an activation function as determined by equation (16).
Download figure:
Standard image High-resolution image

Figure F3. Refer to the following caption and surrounding text. — **Figure F3.** The maximally pooled features, and therefore the output of the first convolutional layer, Ξ¹, as determined by equation (17).
Download figure:
Standard image High-resolution image

Figure F4. Refer to the following caption and surrounding text. — **Figure F4.** The cross-correlation, ζ², of the features extracted from the previous layer, Ξ¹, with 24 weights kernels, $\mathsf{W}^2_i$ , $i\in\{1,\cdots,24\}$ , as determined by equation (18). The cross-correlation may result in positive or negative values, as indicated by the colour bars.
Download figure:
Standard image High-resolution image

**Figure F4.** The cross-correlation, ζ², of the features extracted from the previous layer, Ξ¹, with 24 weights kernels, $\mathsf{W}^2_i$ , $i\in\{1,\cdots,24\}$ , as determined by equation (18). The cross-correlation may result in positive or negative values, as indicated by the colour bars.
Download figure:
Standard image High-resolution image

Figure F5. Refer to the following caption and surrounding text. — **Figure F5.** The post-activations, ξ², as determined by equation (19).
Download figure:
Standard image High-resolution image

Figure F6. Refer to the following caption and surrounding text. — **Figure F6.** The maximally pooled features, and therefore the second layer outputs, Ξ², as determined by equation (20).
Download figure:
Standard image High-resolution image

Figure F7. Refer to the following caption and surrounding text. — **Figure F7.** The cross-correlation, ζ³, of the features extracted from the previous layer, Ξ², with 48 weights kernels, $\mathsf{W}^3_i$ , $i\in\{1,\cdots,48\}$ , as determined by equation (21). The cross-correlation may result in positive or negative values, as indicated by the colour bars.
Download figure:
Standard image High-resolution image

**Figure F7.** The cross-correlation, ζ³, of the features extracted from the previous layer, Ξ², with 48 weights kernels, $\mathsf{W}^3_i$ , $i\in\{1,\cdots,48\}$ , as determined by equation (21). The cross-correlation may result in positive or negative values, as indicated by the colour bars.
Download figure:
Standard image High-resolution image

Figure F8. Refer to the following caption and surrounding text. — **Figure F8.** The post-activations, ξ³, as determined by equation (22).
Download figure:
Standard image High-resolution image

Figure F9. Refer to the following caption and surrounding text. — **Figure F9.** The maximally pooled features, and therefore the final output of the image processing pipeline, Ξ³, as determined by equation (23).
Download figure:
Standard image High-resolution image