Imaging the collective excitations of an ultracold gas using statistical correlations

Advanced data analysis techniques have proved to be crucial for extracting information from noisy images. Here we show that principal component analysis can be successfully applied to ultracold gases to unveil their collective excitations. By analyzing the correlations in a series of images we are able to identify the collective modes which are excited, determine their population, image their eigenfunction, and measure their frequency. Our method allows to discriminate the relevant modes from other noise components and is robust with respect to the data sampling procedure. It can be extended to other dynamical systems including cavity polariton quantum gases or trapped ions.


I. INTRODUCTION
In the past few years, the degree of control of cold atom experiments has increased to an impressive level, from the control of the atomic interactions [1] and the trapping geometry [2], to the creation and observation of many-body correlated systems [3] or the control at the single atom level [4]. In order to extract quantitative measurements from such experiments one has to analyze a large number of images [5], which are fitted and compared to theoretical models [6]. For instance mean-field models describe remarkably well quantum gases at low temperature, including their dynamics [5,7,8]. However these simple models are far from exploiting all the information contained in the images.
This has motivated the development of alternative model-free approaches to analyze the experimental data. For example, with the minimal assumption that the image represents accurately the gas density profile, one can directly compute averaged observables to reveal the gas collective dynamics [9]. It is also quite efficient to represent the signal in the frequency domain, using Fourier transforms, to isolate the system response to a resonant excitation [10][11][12]. In some situations the noise itself contains a lot of information on the system [13] that can be recovered by studying the correlations within the images [14,15].
Here we show that a generic method of signal analysis, Principal Component Analysis (PCA) [16], provides a unique tool to extract all the relevant information from cold atom absorption images, without having to rely on a specific model. This tool has already been used to perform filtering [17][18][19], extract the phase in an interferometric signal [20,21] and identify the main noise sources in an experiment [22]. Recently it has been shown that PCA can be of interest to perform quantum state tomography [23]. As part of multivariate signal analysis methods PCA is widely used in numerous applications dealing with large amounts of data [16], to extract signals from a noisy background. * romain.dubessy@univ-paris13.fr The main result of this paper is that PCA can be extended to the study of the elementary excitations of an ultracold atomic gas and allow the direct observation of the system normal modes. Normal modes or Bogoliubov modes of ultracold atomic gases are the elementary low energy excitations of the system [6,24,25]. They provide a unique insight into the system properties. For example they can reveal the collective superfluid behaviour of Bose [26] and Fermi [27][28][29] gases or probe the system dimensionality [9,30]. Recently an analysis of a set of absorption images using time to frequency domain transformation [10][11][12] has been used to isolate a few low energy collective modes and study their damping. Having access to a method for data analysis which extracts the maximum information will be highly relevant for these studies.
This paper is organized as follows: the PCA method for noise filtering is discussed in section II. We then show in section III that the PCA enables a precise identification of the system low energy excitations. To support our analysis of the experimental data we compare our findings to the results of numerical simulations in section IV. Finally we discuss in section V the requirements for applying PCA to cold atom experiments and the possible improvements that may be achieved.

II. PRINCIPAL COMPONENT ANALYSIS
Let us briefly recall how PCA proceeds [16]. More detail (including formulas) is given in A. We start from a particular data set, which in our case is an ensemble of absorption images where the signal is proportional to the integrated atomic density. We first compute the average of the data set and subtract this mean image from all the images, thus obtaining an ensemble of centered images. We then compute the covariance matrix of this ensemble. The diagonal elements of the covariance matrix contain the variance of the pixels and off-diagonal elements quantify correlations between pixels. By diagonalizing this matrix we recover the eigenvectors, called Principal Components (PCs), which are thus uncorrelated. This statistical independence ensures that uncorrelated noise sources are associated to different principal components [16].
Our experiment is described in detail in reference [31]. Briefly, we produce a quantum degenerate gas of 87 Rb atoms confined in a radio-frequency (rf) dressed magnetic quadrupole trap. We can dynamically control the precise trap shape by varying the magnetic or rf fields, which results in selective excitations of the gas normal modes [30,31]. We measure the gas properties by performing in-situ absorption imaging along the strongly trapped vertical direction. The peak optical density is kept below 6 by repumping only a small fraction of the cloud from the F = 1 hyperfine ground state to the cycling transition. We carefully calibrate the imaging system following reference [32]. The gases we consider in this paper are in the quasi two-dimensional regime: the excitations along the imaging axis are frozen and the dynamics occurs only in the horizontal plane. In this plane the system is well described by a harmonic oscillator [31]. We apply the PCA to the study of the mode dynamics in an anisotropic quasi two-dimensional gas (µ/( ω z ) ∼ 1. 5) with ω x = 2π × 33 Hz and ω y = 2π × 44 Hz.  As an example of application, figure 1 displays the outcome of PCA applied to 27 images acquired in the same conditions. Due to variations of the stray light during image acquisition, fluctuations of atom number in the experiment or mechanical vibrations, the images are not exactly identical. The PCA decomposition identifies all these sources of noise, and we can identify them with a principal component: figure 1c probably accounts for atom number fluctuations while figure 1b and d indicates a small jitter of the camera position [17,22]. Higher order components, see figure 1e to h, reflect the presence of diffraction fringes on the probe beam intensity profile. For each of these components, the corresponding eigenvalue accounts for the fraction of the total variance due to the associated noise source.
Conversely, when the data set results from the variation of a parameter, PCA allows to probe the sensitivity of the system to this parameter. In particular, if the system is initially excited and evolves, measurements taken at different times allow to recover this variation as a principal component. In this paper, we exploit this possibility to directly measure the normal modes of a quantum degenerate gas confined in a harmonic trap.

III. EVIDENCING THE EXCITED MODES
We make use of our highly versatile trap potential to excite simultaneously several low energy eigenmodes of the gas. Namely we displace the trap minimum, we rotate the trap axes and change slightly the trap frequencies.
In the new trap the gas is strongly out of equilibrium and we record its evolution by taking images for different holding times in the trap, covering a time span of 100 ms. Figure 2 shows the result of PCA for this data set (133 images). Compared to figure 1 we see that the principal components have changed.
Let us now identify the first principal components. The first two PCs (see figure 2b and c) display a two-lobe pattern oriented respectively along the columns and the rows of the images: this is characteristic of a dipole oscillation of the cloud. This center of mass motion is due to the trap minimum displacement during the excitation process. The third PC (see figure 2d) indicates a global variation of the signal over the whole cloud, which can be interpreted as atom number fluctuations in the experiment. Part of these fluctuations are due to the fact that the lifetime in the trap is limited and atoms are lost as the holding time increases [41].The fourth PC (see figure 2e) possesses a striking spatial pattern with four lobes characteristic of a scissors excitation [33]. Note that the lobes are oriented at 45 • with respect to the trap axes (aligned with the first two PCs) as expected. The next two PCs (see figure 2f and g) look like compression modes of the gas, with a density depletion at the center of the cloud and a correlated augmentation of the density on the sides of the cloud.
The PCs are presented by decreasing eigenvalue, meaning that they account for less and less variations in the original data set. For this particular experiment the center of mass oscillation is the dominant excitation in the cloud, followed by the response to the rotation of the trap axis and marginally by the compression of the trap. In another experiment (not shown) where the trap rotation was not performed we have verified that no PC displayed the spatial pattern of a scissors excitation.
This analysis of the PCs is supported by the study of the time dependent oscillations of the associated weights, computed by projecting the centered original data set on  the PCs. The result of this computation is displayed on figure 3, for the dipoles and scissors components. Let us focus first on the first two weights: they exhibit sinusoidal oscillations at the expected trap frequencies (44 Hz and 33 Hz). This supports the fact that PCA has correctly identified as independent components the center of mass motion of the cloud along the trap axes. The scissors component displays a more complex oscillation pattern. We find that the best fit to the data is given by a sum of three sinusoides, at frequencies 12 Hz, 55 Hz and 77 Hz. This is related to the fact that the scissors component found by PCA is sensible both to the collective response of the superfluid part and to the collisionless oscillations of the normal part of the gas [33]. The simultaneous presence of these three frequencies has been evidenced in a three dimensional Bose-Einstein condensate [34] where simultaneous measurement of the superfluid and normal part of the cloud rotations were obtained by a bimodal distribution fit to the density profiles. Here we note that the same PC gives access to both the superfluid and the normal response to the rotation of the trap axes which might be used to measure their relative amplitudes. Let us stress that the PCA is able to separate the contributions of the different modes in a given experiment which could help to design better excitation pattern or focus on higher order modes [11]. In particular, being able to measure on the same data set the dipole mode frequencies gives access to the natural system clock [24]. Therefore PCA gives access to direct comparison between measured frequencies and predictions. Moreover, for the data set used in figure 2, we find that the simple hydro-dynamic models of references [7,8] fail to extract these frequencies present in the oscillation of the density, probably due to the fact that several collective modes are simultaneously excited. In this case it is really essential to use a model-free approach to analyze the data set.

IV. COMPARISON TO NUMERICAL SIMULATIONS
We pursue our investigation numerically in order to compare the principal components to normal modes. We use a zero temperature two-dimensional mean field model of our cloud and perform a numerical time-dependent simulation which mimics the experimental sequence. We then extract the simulated density profiles using a regular time sampling, thus obtaining a data set of 152 computed images. We finally compare the PCA of this data set to the actual normal modes of the trap, computed using the Bogoliubov-de Gennes equations. Details on the simulations are given in B.    files quantitatively, we compute the scalar product between the principal components and the eigenmodes images. We find a high degree of overlap for the largest five principal components (dipoles: 99.7% and 99.4%, scissors: 98.5%, monopole-like 98.8% and quadrupolelike 89.2%) when projected onto the corresponding eigenmode. This supports the idea that the largest principal components can indeed be identified with a well defined normal mode (see also C).
To confirm this result, we compare the oscillation frequency of the principal components (obtained by fitting a sinusoidal function to the time dependent weight of the simulated density profiles) ω pca to the frequency of the mode given by the Bogoliubov-de Gennes equations ω diag and to an analytic hydrodynamical model ω th [24]. The results obtained with the data of figure 4 are reported in table I. We find that for the largest principal components the simple sinusoidal behavior correctly fits the data and gives a value compatible with the Bogoliubov-de Gennes theory, within the numerical uncertainty [42].
For the collective modes, we expect the correct value for the mode frequency to be given by the diagonalization procedure, as the hydrodynamical model is only approximate. There is an excellent agreement between ω pca and ω diag for the values reported in table I, thus validating our experimental findings. However this is not true for the principal components with a small variance, which exhibit complex temporal behaviors. We observe that these components do not have a significant overlap with one of the modes: PCA is not able to identify them.
We conclude from this example that PCA provides a robust way of evidencing the dominantly excited modes in an out of equilibrium ultracold gas. Once the relevant components are isolated, PCA allows to extract the mode time dependence, without having to rely on a model for the atomic response.

V. DISCUSSION
We now discuss the requirements for PCA to be efficient and compare it with Fourier analysis. PCA is a statistical method: the data set has to span a sufficiently large number of configurations for the correlations between two different normal modes to average to zero. In particular, to resolve two different modes with close frequencies, the total acquisition time has to be larger than the beat note period. However it is not necessary to use an even sampling during this time period. In addition if the populations in the two excitations are very different, resulting in very different contributions to the variance, PCA separates them efficiently, even for shorter observation times, see the discussion in C.
Fourier transformation methods can also be quite efficient for identifying collective modes [10][11][12]. However they come with stronger constraints: the time sampling must be regular and the total time should be a multiple of the signal time period. It supposes an a priori knowledge of the signal frequency which may have to be determined iteratively. Moreover Fourier transform gives only access to frequencies which are multiple of a fundamental frequency, which complicates the analysis for systems with multiple excitations (see supplementary data). Finally, a white noise contributes to each Fourier component, whereas it is naturally filtered out in PCA. PCA is not subject to such constraints: if we reduce the size of the data set used in section IV, for example by keeping only one image out of ten, PCA is still able to find out PCs close to the excited eigenmodes (dipoles, scissors and monopole-like with 95% fidelity, but the quadrupolelike component is absent, see supplementary data), even if the Nyquist-Shannon sampling theorem is not verified any longer [43]. In that sense PCA is more efficient than Fourier methods.
In conclusion we have shown that, beyond noise filtering [17][18][19], PCA provides a powerful statistical tool to analyze experimental as well as numerical data sets. When applied to time-dependent systems, it allows for a model-free discrimination of the normal modes and to the measurement of their populations and frequencies.
We expect PCA to be particularly relevant for the study of samples where fluctuations play a major role in the physics. Examples include the random creation of defects in the Kibble Zurek mechanism [35] or the correlations between vortices and anti-vortices in a two dimensional superfluid [36]. We note that PCA is a very general method and would be suitable for other systems with time dependent signals. In cavity polariton quantum gases, where images can be taken in real time, PCA will allow to extract the relevant information inside a large data set [37]. Finally, cold trapped ions systems behave as crystals supporting many collective modes, which could be studied using PCA [38].
We envision that PCA is suited to perform Bogoliubov spectroscopy, in the spirit of the method used in references [10,11]. A mode largely excited by a resonant excitation will be easily identified by PCA and its frequency precisely determined by measuring its eigenvalue with respect to the modulation frequency. In contrast to Fourier methods, PCA can identify the dominantly excited mode using samples covering only one oscillation period of the mode either by recording the time evolution or by varying the excitation phase. This property should prove useful in particular to study over-damped modes.

Acknowledgments
LPL is UMR 7538 of CNRS and Paris 13 University. We acknowledge helpful discussions with Aurélien Perrin.

Appendix A: Principal Component Analysis
We provide here a short recipe to apply PCA to the analysis of density profiles. Other examples of applications are given in references [16,17]. We stress that the mathematical formalism is quite simple and that most data analysis softwares provide standard implementation of PCA. We are interested in density profile images and assume that the pixels of each image are stored (row wise) in a single vector. The first step is to center the data set by computing the average image and subtract it to each image. The whole data set can then be stored in an N ×P matrix, denoted B, where N is the number of images and P is the number of pixels. Thus B i,j contains the j-th pixel of the centered image i.
Next we want to compute the eigenvalues of the covariance matrix S = B T B/ (N −1) where B T is the transpose of matrix B. This P × P matrix is in general quite large so it is hard to diagonalize it directly. However it is quite simple to show that its rank is at most N . Indeed, assuming that X is an eigenvector of S with eigenvalue λ (meaning SX = λX), it is straightforward to verify that Y = BX is an eigenvector of the square N × N symmetric matrix Σ = BB T /(N −1) with the same eigenvalue λ. Therefore S and Σ have the same spectrum, of at most N real eigenvalues. Knowing an eigenvector Y of Σ, the corresponding eigenvector of S is simply X = B T Y . Finally let us stress that these vectors are orthogonal since the S matrix is real and symmetric. We define the associated PCs by normalizing the eigenvectors to unity. In the case of both a large number of pixels and a large number of images the diagonalization of S and Σ is hard to compute. However since we are a priori only interested in the PCs with the largest variance they can be efficiently computed by iterative methods [39].
The PCs provide an orthonormal basis spanning the subspace of the data set and therefore each original image can be represented as a sum of the mean image and the weighted contributions from each principal component. These weights are obtained by projecting the centered image onto the corresponding principal component. By selecting only relevant principal components the noise can be partially filtered out of the reconstructed images [16].
where t is expressed in units of ω −1 x , x and y in units of a x = /(M ω x ), and ψ ≡ ψ(x, y, t) in units of a −1 x . M is the atomic mass, N the number of atoms and g 2D = √ 8πa/a z is the reduced coupling constant, where a is the contact interaction scattering length and a z = /(M ω z ) is the size of the vertical harmonic oscillator ground-state. The potential reads: where ǫ = ω 2 y /ω 2 x quantifies the trap in plane anisotropy and the arbitrary angle θ allows to rotate the trap axes. The auxiliary parameters x 0 , y 0 and α can be used to induce a trap displacement and compression. Table II details the value of the parameters appearing in (B2) before and after the excitation.
The numerical wave function is represented on a square 128 × 128 grid with an equivalent full width of 15a x in both x and y directions. For the computations we used g 2D N = 1000, matched to the experimental conditions. We use this model to compare the outcome of two numerical computations. On the one hand we mimic the experiment described in section III by a) computing the system ground state for the initial potential using an imaginary time evolution algorithm; b) using this result as the input of a real time evolution in the final potential; and c) performing PCA on regularly sampled density profiles during the evolution. The evolution algorithm relies on a time splitting spectral method, from t = 0 to t = 37.7 (in dimensionless units, corresponding to 6 periods of the weakest trap axis) using a time step of 10 −3 . The total time is chosen to be close to a multiple of both dipole modes oscillation period (6 periods and 8 periods respectively): this ensures that the average density profile computed in PCA is not skewed. The sampling is performed every T s = 0.126. The result of this procedure is shown in the left panel of figure 4.
On the other hand we directly compute the small excitation spectrum using Bogoliubov-de Gennes equations obtained from the linearisation of (B1) around the system ground state in the final trap. This implies the ability to diagonalize a square 2 15 × 2 15 matrix which is quite challenging. Fortunately this matrix is sparse and we are interested only in the lowest energy part of the spectrum, which means we do not have to compute all the eigenstates. We designed a fast custom C program that uses a combination of an iterative method [39] with an efficient sparse matrix library [40] to compute the relevant eigenvectors. The result of this procedure is displayed in the right panel of figure 4.
Appendix C: Identification of the principal components with the normal modes We have shown that PCA is very efficient to identify the normal modes of an excited ultracold gas. This may be surprising but can be understood in the framework of small excitations. Using a hydrodynamic model [24], the gas out of equilibrium density profile may be expanded as: ρ(r, t) = ρ 0 (r) + k c k cos [ω k t + φ k ] f k (r), where k labels the normal mode of frequency ω k , f k (r) describes the mode normalized density profile and c k is related to the mode population.
In the experiment we observe the gas only at discrete times {t n } n∈ [1,N ] and positions {r i } i∈[1,P ] and therefore we can write the i-th pixel of the n-th image as ρ(r i , t n ) = ρ 0 (r i ) + k c k cos [ω k t n + φ k ] f k (r i ) + ε(r i , t n ), where we added a pixel and time dependent noise contribution ε(r i , t n ). PCA starts with the evaluation of the centered data set, by averaging over the sampling times {t n }: B n,i = k c k cos [ω k t n + φ k ] f k (r i ) + ε(r i , t n ) + δ(r i ), where the δ(r i ) term is close to zero for a total sampling time T ≫ (min k ω k ) −1 .
Then the S matrix elements can be written as: where the term ∆(r i , r j ) is the effective noise covariance between pixels i and j, due to the initial noise distribution and finite sampling size induced errors. Providing that the ∆(r i , r j ) term is small enough it is straightforward to verify that the principal components of matrix S are the vectors {f k (r i )} i , with eigenvalue ∼ c 2 k /2. In particular this is true for T ≫ (min k =k ′ |ω k − ω k ′ |) −1 . This constraint on T is more stringent than the previous one, especially when two normal modes are close to degeneracy. However if these two modes have small populations they have a small contribution to the ∆(r i , r j ) term.
The conclusion of this analysis is twofold. On the one hand PCA correctly identifies the most excited[44] eigenmodes of the system. On the other hand, the total time sampling should be large enough to resolve the beat note between these modes. For practical purposes we empirically found that taking T equal to one beat note period is sufficient, see for example figure 3.

Supplementary Material
We provide here supplementary material to support the results and conclusions of our paper. Figure 5 shows the result of a Thomas-Fermi fit [7,8] to the experimental data for the data set of figure 2. While the center of mass oscillations are correctly extracted, the angle of the cloud shows no clear oscillating time variation, as would be expected for a scissors excitation [33]. For this particular experiment, PCA is able to extract more information on the system than a fitting procedure, see figure 3.  Figure 6 displays the time dependent weights of the simulated density profiles on the first five principal components. We empirically find that components with a variance above 1% of the total variance follow a sinusoidal oscillation, while components with a smaller variance are not well identified and do not oscillate at a single frequency. Figure 7 shows a Fourier analysis of the simulated density profiles. Because of the finite time sampling only discrete frequencies are present in the analysis. The scissorslike and monopole-like density profiles are spread over several frequency components: this FFT method is thus less efficient to decorrelate the modes than PCA, for this particular time sampling. Obviously by choosing a better time sampling the FFT signal could be optimized. However PCA provides a more straightforward method to identify multiple excitations, with less constraints on the data sampling. Figure 8 displays the result of PCA applied to the simulated density profiles where the result of the simulation has been undersampled by keeping only one image out of ten. PCA is still efficient and identifies four modes with 96% fidelity out of the five dominantly excited modes.  In particular for the monopole-like PC (figure e), there is on average less than one sampling point per period.