Excitation spectrum of a trapped dipolar supersolid and its experimental evidence

TL;DR

This paper investigates the excitation spectrum of a dipolar Bose gas transitioning from superfluid to supersolid, revealing two distinct excitation branches linked to crystal and superfluid orders, supported by experimental evidence.

Contribution

It provides a theoretical and experimental analysis of the elementary excitations in a dipolar supersolid, demonstrating the coexistence of crystal and superfluid excitation modes.

Findings

Two excitation branches appear in the supersolid phase.

Experimental two-frequency response confirms broken symmetries.

Finite-size effects are minimal in the excitation spectrum.

Abstract

We study the spectrum of elementary excitations of a dipolar Bose gas in a three-dimensional anisotropic trap across the superfluid-supersolid phase transition. Theoretically, we show that, when entering the supersolid phase, two distinct excitation branches appear, respectively associated with dominantly crystal and superfluid excitations. These results confirm infinite-system predictions, showing that finite-size effects play only a small qualitative role, and connect the two branches to the simultaneous occurrence of crystal and superfluid orders. Experimentally, we probe compressional excitations in an Er quantum gas across the phase diagram. While in the Bose-Einstein condensate regime the system exhibits an ordinary quadrupole oscillation, in the supersolid regime we observe a striking two-frequency response of the system, related to the two spontaneously broken symmetries.

Equations

$$\mathrm{i}\hslash \frac{\partial \psi (\mathbf{r},t)}{\partial t}$$

$$\displaystyle\hskip 56.9055pt+\Delta\mu[n]\Big{)}\psi(\mathbf{r},t),$$

$$\psi (\mathbf{r},t)=\left({\psi }_0(\mathbf{r})+\eta \delta \psi (\mathbf{r},t)\right)e^{-i\mu t},$$

$$\delta \psi (\mathbf{r},t)=u_l{e}^{-i{ϵ}_{l}t/\hslash }+v_l^{\ast }{e}^{i{ϵ}_{l}t/\hslash }.$$

$$S(\mathbf{k},{\omega }^{\prime })$$

$$\psi (\mathbf{r},t)e^{i\mu t}\approx \sqrt{|{\psi }_0(\mathbf{r}){|}^2+2\eta \delta {\rho }_l(\mathbf{r})\cos {\omega }_{l}t}{e}^{-i\eta \delta {\phi }_l(\mathbf{r})\sin {\omega }_{l}t},$$

$$V_{\text{in(out)}}$$

$$\displaystyle\langle f(y)\rangle_{R_{j}}=\int_{y\in R_{j}}f(y)\,dy\bigg{/}\int_{y\in R_{j}}\,dy.$$

$$n_i(k_x,k_y)=\bar{n}(k_x,k_y)+\sum _{p=1}^{N_s}{w}_{p,i}{\mathcal{C}}_p(k_x,k_y),$$

$$n(\mathbf{r},t)\approx n_0(\mathbf{r},t)+2\sum _l\eta \delta {\rho }_l(\mathbf{r})\cos \left({\omega }_{l}t+{\varphi }_l\right),$$

$$n^{(p)}(k_x,k_y,t)=\bar{n}(k_x,k_y)+W_p(t){\mathcal{C}}_p(k_x,k_y).$$

Full text

Supersolidity – a paradoxical quantum phase of matter that combines crystal rigidity and superfluid flow – was suggested more than half a century ago as a paradigmatic manifestation of a state in which two continuous symmetries are simultaneously broken Boninsegni and Prokof’ev (2012). In a supersolid, the spontaneously broken symmetries are the gauge symmetry, associated with the phase coherence in a superfluid, and the translational invariance, signalizing crystalline order. The striking aspect is that, in a supersolid of indistinguishable bosons, the same particles are participating in developing such two apparently antithetical, yet coexisting, orders. Originally predicted in quantum solids with mobile bosonic vacancies Andreev and Lifshitz (1969); Chester (1970); Leggett (1970), the search for supersolidity has fueled research across different areas of quantum matter from condensed matter to atomic physics, including quantum gases with nonlocal interparticle interactions  Henkel et al. (2010); Cinti et al. (2010); Boninsegni (2012); Lu et al. (2015); Macia et al. (2016); Cinti and Boninsegni (2017); Wenzel et al. (2017); Baillie and Blakie (2018); Chomaz et al. (2018); Roccuzzo and Ancilotto (2019); Böttcher et al. (2019a); Tanzi et al. (2019a); Chomaz et al. (2019); Kora and Boninsegni (2019); Zhang et al. (2019).

Recent experiments have revealed that axially elongated dipolar quantum gases can undergo a phase transition from a regular Bose-Einstein condensate (BEC), possessing a homogeneous density in the local-density-approximation sense, to a state with supersolid properties, where density modulation and global phase coherence coexist Böttcher et al. (2019a); Tanzi et al. (2019a); Chomaz et al. (2019). Such experiments, complementing the ones with BECs coupled to light Léonard et al. (2017a); Li et al. (2017); Léonard et al. (2017b), have opened a whole set of fundamental questions, covering the very real meaning of superfluidity in a supersolid state, its shear transport, and phase rigidity.

Of particular relevance is the study of the spectrum of elementary excitations, which governs the system’s response to perturbations Landau (1947); Bogoliubov (1947); Pitaevskii and Stringari (2016). Typically, phase transitions occur in concomitance with drastic modifications of the excitation spectra, – as in the case of the emergence of roton excitations in He II or the phononic dispersion for BECs – and similar dramatic changes are expected when crossing the superfluid-supersolid transition. Theoretical studies of uniform (infinite) gases with periodic boundary conditions and soft-core Saccani et al. (2012); Macrì et al. (2013); Rossotti et al. (2017) or dipolar interactions  Macia et al. (2012); Bombin et al. (2017); Roccuzzo and Ancilotto (2019), have shown two distinct branches appearing in the excitation spectrum of a supersolid state – one for each broken symmetry. Their coexistence has been identified as an unambiguous proof of supersolidity, being the direct consequence of the simultaneous presence of superfluid and crystalline orders Andreev and Lifshitz (1969); Pomeau and Rica (1994); Saccani et al. (2012); Macrì et al. (2013).

An important issue is to understand if these trademarks survive – and can be measured – in the experimentally relevant regimes of a finite-size quantum gas, confined in all three spatial dimensions. In this Letter, we address these points by performing full spectrum calculations and by experimentally exciting collective modes in an erbium quantum gas. Both theory and experiment show the existence of two distinct classes of excitations, one connected to crystal modes and the other to phase modes, providing the finite-size equivalent of the two-branches spectrum for infinite systems.

In our study, we consider a three-dimensional dipolar quantum gas confined in an axially elongated ($y$) harmonic trap with transverse orientation ($z$) of the atomic dipoles. These systems are well described by an extended Gross-Pitaevskii equation (eGPE), including nonlinear terms, accounting for contact interactions depending on the scattering length $a_{\rm s}$, the anisotropic long-range dipole-dipole interaction (DDI), and quantum fluctuations in the form of a Lee-Huang-Yang type of correction  Wächtler and Santos (2016a); Bisset et al. (2016); Wächtler and Santos (2016b); Schmitt et al. (2016); Chomaz et al. (2016); Baillie and Blakie (2018); Roccuzzo and Ancilotto (2019); Tanzi et al. (2019a); Chomaz et al. (2019); Böttcher et al. (2019a); Zhang et al. (2019); see also Sup . We calculate ground-state wavefunctions, $\psi_{0}(\mathbf{r})$, by minimizing the energy functional resulting from the eGPE using the conjugate-gradients technique Ronen et al. (2006). As shown in Fig. 1 (insets), the ground state evolves with decreasing $a_{\rm s}$ from a regular BEC (a,b) to a supersolid state with axial density-wave modulation (c-f) and finally to an insulating array of independent droplets (ID) (g-h) Boninsegni (2012); Macrì et al. (2013); Roccuzzo and Ancilotto (2019); Böttcher et al. (2019a); Chomaz et al. (2019).

The spectrum of elementary excitations is calculated by numerically solving the Bogoliubov de Gennes equations, which are obtained from an expansion of the macroscopic wavefunction as $\psi(\mathbf{r},t)=\left[\psi_{0}(\mathbf{r})+\eta(u_{l}e^{-i\epsilon_{l}t/\hbar}+v_{l}^{*}e^{i\epsilon_{l}t/\hbar})\right]e^{-i\mu t}$ with $\eta\ll 1$ and linearizing the eGPE around $\psi_{0}$ Pitaevskii and Stringari (2016); Ronen et al. (2006); van Bijnen et al. (2010); Chomaz et al. (2018). Here, $\mu$ is the ground-state’s chemical potential. By solving the resulting eigenvalue problem, we find a set of discrete modes, numbered by $l$, of energy $\epsilon_{l}=\hbar\omega_{l}$ and amplitudes $u_{l}$ and $v_{l}$. We calculate the dynamic structure factor (DSF), $S(k,\omega)$, which informs on the system’s response when its density is perturbed at a given modulation momentum $k$ and with an energy $\hbar\omega$ Blakie et al. (2002); Brunello et al. (2001); Pitaevskii and Stringari (2016). Whereas in the absence of an external trap the spectrum is continuous and the DSF is a $\delta$-peak resonance at the Bogoliubov mode $(\omega_{l},k_{l})$, the confining potential yields instead a discretization of the excitation spectrum and a $k$ broadening in $S(k,\omega)$. For a given energy (i. e. a single mode), finite-size effects may even yield several peaks in $k$, see e. g. three-peak structures at large energy in Fig. 1 (a,b). For the considered parameters, these finite-size effects are more pronounced in Er than Dy, since the latter exhibits a larger number of maxima in the density-modulated phases, rendering its excitation spectrum more reminiscent of the infinite-system case; see Fig 1.

Figure 1 shows the calculated excitation spectrum for ground-states in the regular BEC, the supersolid, and the ID phases for a Dy (upper row) and Er (lower row) quantum gas. In the BEC regime close to the supersolid transition (a, b), the spectrum of excitations shows a single excitation branch with the characteristic phonon-maxon-roton dispersion of a BEC O’Dell et al. (2003); Santos et al. (2003); Blakie et al. (2012); Bisset et al. (2013); Jona-Lasinio et al. (2013), as recently measured Petter et al. (2019). When the roton fully softens (at $a_{\rm s}=a_{\rm s}^{*}$), the ground-state becomes density modulated with a wave number close to the roton one, $k_{\rm rot}$. Here, the excitation spectrum develops additional structures, marked by the appearance of nearly degenerate modes; (c, d). When lowering $a_{\rm s}$, we find that these modes start to separate in energy, where some harden and the others soften, and two excitation branches become visible; (e, f). This result resembles that of infinite systems, where the broken translational and gauge symmetry are each associated with the appearance of one excitation branch Saccani et al. (2012); Macrì et al. (2013); Roccuzzo and Ancilotto (2019). Additionally, we observe that the spectrum acquires a periodic structure, reminiscent of Brillouin zones in a crystal, with reciprocal lattice constant $k\simeq k_{\rm rot}$. Modes with an energy higher than the maxon (energy maximum at $k<k_{\rm rot}$) seem to have a single-droplet-excitation character, and they will be the subject of future investigations. When further decreasing $a_{s}<a_{\rm s}^{*}$, the lower-lying branch decreases both in energy and in DSF values, whereas the opposite occurs for the higher branch. Eventually, when reaching the ID regime, the lower branch progressively vanishes, underlying the disappearance of global superfluidity; (g, h).

We focus on the properties of the excitation spectrum in the supersolid regime. The interesting question is how the two branches relate to the two orders in the systems, crystal and superfluid. To gain insight, we study the system’s dynamics when a single mode $l$ is excited with amplitude $\eta\ll 1$ by writing $\psi(t)e^{i\mu t}\approx\sqrt{|\psi_{0}|^{2}+2\eta\delta\rho_{l}\cos\omega_{l}t}e^{-i\eta\delta\varphi_{l}\sin\omega_{l}t}$, in terms of density perturbations $\delta\rho_{l}=(u_{l}+v_{l}^{*})\left|\psi_{0}\right|$ and phase perturbations $\delta\varphi_{l}=(u_{l}-v_{l}^{*})/\left|\psi_{0}\right|$. The subsequent time evolution of the axial density profile is shown in Fig. 2(a-c) for three relevant cases. For simplicity, only the two extremes of the mode oscillation are shown. The mode character can be understood by noting that phase gradients correspond to mass currents. Large gradients inside a density peak imply motion of the density peak (e.g. Fig. 2(a)) and relate to crystal modes. Large phase gradients between density peaks, signify a superfluid current of particles tunneling from one density peak to another (e. g. Fig. 2(b)), and are associated with phase modes. However, in our system, the phase/crystal mode classification is not strict and we find that these two characters mix; see Fig. 2 (a-c). Particularly, we observe both behaviors simultaneously in Fig. 2 (c). Such a mixing is expected from the long-range nature of the DDI, coupling density, and position of the peaks Saccani et al. (2012); Macrì et al. (2013). Note that the character of the mode can change with $a_{\rm s}$. For instance, the mode in Fig.2 (c) develops an almost pure crystal character for decreasing $a_{\rm s}$. To quantify a mode’s character, we plot in Fig. 2 (d) the DSF spectrum at a fixed $a_{\rm s}$, colored according to the ratio $C$ of phase variances inside, and between the density peaks Sup . This allows us to differentiate the dominant character of the two branches, being phase type for the lower branch, and crystal type for the upper one.

To test our predictions, we experimentally study the collective excitations in an erbium quantum gas across the BEC-supersolid-ID phases. We prepare a BEC at $a_{\rm s}=64\,a_{0}$. The atoms are confined in an axially elongated optical-dipole trap of harmonic frequencies $2\pi\times(\nu_{x},\nu_{y},\nu_{z})=2\pi\times(259(2),30(1),170(1))$ Hz and polarized along $z$ by an external magnetic field; see Chomaz et al. (2018, 2019). To probe our system, we perform standard absorption imaging after 30 ms of time-of-flight expansion, yielding measurements of the momentum space density $n(k_{x},k_{y})$ Sup . Using the tunability of the contact interaction via magnetic Feshbach resonances Chin et al. (2010), we can prepare the system at desired locations in the phase diagram – in the BEC, supersolid, or ID phase – by linearly ramping down $a_{\rm s}$ in 20 ms to the target value. We then allow the system to stabilize for 10 ms. At this point, we record an atom number of typically $5\times 10^{4}$ for the supersolid regime. We confirmed the relevant $a_{\rm s}$ ranges by repeating the matter-wave interferometric analysis of Ref. Chomaz et al. (2019). While in the BEC region the momentum distribution shows a regular, nearly Gaussian single peak, in the supersolid regime the in-trap density modulation gives rise to coherent interference patterns along $k_{y}$, consisting of a central peak with two lower-amplitude side peaks; see Fig. 3(a).

After preparing the system in the desired phase, we excite collective modes in the gas by suddenly reducing the axial harmonic confinement to 10$\%$ of its initial value (i. e. $\nu_{y}\approx 3\,$Hz) for 1 ms, before restoring it again. The atomic cloud is subsequently held for a variable time $t_{\rm h}$, before releasing it from the trap and recording the time evolution of $n(k_{x},k_{y})$. As the lifetime of the supersolid state is limited to around 40 ms Chomaz et al. (2019), we focus on $t_{\rm h}\leq 30$ ms. As expected, in the BEC phase, we predominantly observe an oscillation of the axial width, connected to the lowest-lying quadrupole mode Pitaevskii and Stringari (2016). In the supersolid regime, the situation is more complex; see Fig. 1(c-f). Here, multiple modes, of both crystal and phase character, can be simultaneously populated, resulting in a convoluted dynamics of the interference pattern.

We therefore employ a model-free statistical approach, known as principal component analysis (PCA) Jolliffe (2002), to study the time evolution of the measured interference patterns at a fixed $a_{\rm s}$. This method has been successfully used to study e. g., matter-wave interference Segal et al. (2010) and collective excitations Dubessy et al. (2014) in ultracold-gas experiments. The PCA analyzes the correlations between pixels in a set of images, decomposes them into statistically independent components, and orders these principal components (PCs) according to their contributions to the overall fluctuations in the dataset.

In a dataset probing the system dynamics after an excitation, the PCA can identify the elementary modes with the PC’s weights in the individual images exhibiting oscillations at the mode frequencies Sup ; Dubessy et al. (2014). We apply the PCA to the time-evolution of the interference patterns after the trap excitation. Figure 3 (b) shows the PCA results in the supersolid regime at $a_{\rm s}=49.8\,a_{0}$. We identify two leading PCs, which we label as PC1 and PC2. Their weights oscillate with different amplitudes and at distinct frequencies, namely 41(1)Hz for PC1 and 52(5)Hz for PC2. The comparison between the measured frequencies and the theoretically calculated mode energies indicates that, following our trap excitation, the second and third lowest lying even-mode are simultaneously populated. As shown in Fig. 2 (b) and (c), these modes possess a phase and a mixed character, respectively. Note that we apply an overall shift of $-4.3\,a_{0}$ to the $a_{\rm s}$ value for the experimental data; for more details see the discussion in Refs. Petter et al. (2019); Böttcher et al. (2019b).

To visualize the role of each PC on the interference-pattern dynamics, we apply a partial recomposition of the images, accounting only for the PC of interest; see Sup . The effect of PC1 on the axial dynamics is shown in Fig. 3(c), mainly being an axial breathing of the central peak, accompanied by weaker in-phase breathing of the side peaks. Instead, PC2 exhibits a dominant variation of the side-peak amplitude; see Fig. 3 (d). These results show a good agreement with the calculated time evolutions of the interference patterns for the second and third even modes, shown in Fig. 3 (e-f).

Finally, we study the evolution of the modes across the BEC to supersolid and ID phases. We repeat the collective excitation measurements for various $a_{\rm s}$, and, using the PCA, we extract the oscillation frequencies of all the leading PCs. Figure 4 shows our experimental results together with the mode tracking from the spectrum calculations. For a give elementary mode $l$, we plot $\omega_{l}$ as well as the response amplitude $R_{l}=m\omega_{y}^{2}\langle\ell|\hat{y}^{2}|0\rangle/2\hbar\omega_{l}$, which indicates the probability to be excited by our trap-excitation scheme. For completeness, the figure shows both even and odd modes, although only even modes are coupled to our trap-excitation scheme. Here, $|0\rangle$ and $|\ell\rangle$ denote, respectively, the ground and excited states of interest, and $\hat{y}$ is the axial position operator.

In the BEC regime, besides the roton mode that progressively softens with decreasing $a_{\rm s}$, the other modes show a regular spacing in energy and are nearly constant with $a_{\rm s}$. In both the theory and experiment, we observe that just one mode couples to the trap-excitation scheme. This mode has a compressional, axial breathing character. Experimentally, we observe that all the leading PCs oscillate at the same frequency, suggesting that they account for the same mode Sup . In this regime, both the PCs frequencies, $\omega_{l}$, and $R_{l}$ remain rather constant. At the supersolid phase transition, reached around $a_{\rm s}=50.6\,a_{0}$, the numerical calculations reveal that different modes undergo an abrupt change and can mix with each other. Their energy and phase/crystal character exhibits a strong dependence on $a_{\rm s}$. Here, several modes respond to the trap-excitation scheme, as shown by the value of $R_{l}$. In the PCA we observe that the leading PCs now oscillate at distinct frequencies and have different characters (see also Fig. 3). One set of PCs reduces their frequency when lowering $a_{\rm s}$, indicating (at least) one phase mode that softens strongly in the supersolid regime, even below the trap frequency $\nu_{y}$. Another set of PCs shows a frequency that remains hard when decreasing $a_{\rm s}$. Calculations of $C$ show that this mode changes character along the phase diagram and eventually becomes crystal type.

In conclusion, the overall agreement between the experiment and theory confirms the calculations in the supersolid regime, revealing two distinct branches with respective crystal and superfluid characters. The trademarks of supersolidity expected in infinite systems thus carry over to the finite-size ones currently available in laboratories. The knowledge of the excitation spectrum will provide the base for future investigations related to the superfluid properties and phase rigidity in a supersolid state.

• Correspondence and requests for materials should be addressed to Francesca.Ferlaino@uibk.ac.at.

.1 Calculation of the Bogoliubov Spectrum

Our theory is based on an extended version of the Gross-Pitaevskii equation (eGPE)

[TABLE]

where $\psi(\mathbf{r},t)$ is the dipolar quantum-gas’ wave function $\psi(\mathbf{r},t)$. The eGPE includes the kinetic energy, external trap potential and the mean-field effect of the interactions Pitaevskii and Stringari (2016); Baranov (2008). The first three terms of Eq. (.1) account for the kinetic energy, the external harmonic trapping potential, and the mean-field interactions, respectively. The latter includes the contact and the dipolar interactions. In order to study the supersolid phase, it is fundamental to also include a beyond-mean-field corrections in order to stabilize the supersolid state against the roton instability. This is done by adding a term in the form of the Lee-Huang-Yang correction, $\Delta\mu[n]$ Wächtler and Santos (2016a); Bisset et al. (2016); Wächtler and Santos (2016b); Schmitt et al. (2016); Chomaz et al. (2016); Baillie and Blakie (2018); Roccuzzo and Ancilotto (2019); Tanzi et al. (2019a); Chomaz et al. (2019); Böttcher et al. (2019a); Zhang et al. (2019); see also Gammal et al. (2000); Bulgac (2002); Lu et al. (2015); Petrov (2015). This is typically included as a correction to the chemical potential obtained under the assumption of local density approximation Lima and Pelster (2011); Lee and Yang (1957). However, recent experimental results have raised the questions about the range of validity of such a treatment since quantitative disagreements at a level of few $\%$ have been observed when comparing the theory results with the experimental findings Schmitt et al. (2016); Chomaz et al. (2018); Ferrier-Barbut et al. (2018); Cabrera et al. (2018); Petter et al. (2019); Böttcher et al. (2019b). To the best of our knowledge, this is still an open question, which will need future additional theoretical investigations. To compensate for this effect, throughout this letter, we shift $a_{\rm s}$ by $-4.3a_{0}$. To calculate the ground-state (GS) wave-function, $\psi_{0}(\mathbf{r})$, we then minimize the energy functional resulting from the eGPE using the conjugate-gradients technique Ronen et al. (2006).

In a next step, we study the Bogoliubov de Gennes (BdG) excitation spectrum of a dipolar Bose-Einstein condensate trapped in a harmonic cigar shaped potential Pitaevskii and Stringari (2016); Ronen et al. (2006). Our calculations are obtained by expanding the wavefunction $\psi(\mathbf{r},t)$ around $\psi_{0}(\mathbf{r})$. Here, we write:

[TABLE]

where $\eta\ll 1$, $\mu$ is the chemical potential of the ground state and

[TABLE]

The spatial modes $u_{l}$ and $v_{l}$ are oscillating in time with the corresponding frequency $\omega_{l}=\epsilon_{l}/\hbar$. We then linearize the eGPE around $\psi_{0}$ at first order in $\eta$. By solving the set of coupled linear equations, we obtained the discrete modes, numbered by $l$, of energy $\epsilon_{l}$ and amplitudes $u_{l}$ and $v_{l}$. We define the (odd) even parity of the mode from their amplitude $u_{l}$ and $v_{l}$ being (anti-)symmetric in $y$.

In order to illustrate the spectrum, we compute the dynamic structure factor (DSF), since it directly gives information about the density response of the system when perturbed at specific energies and momenta. At $T=0$ the DSF is defined as Blakie et al. (2002); Chomaz et al. (2018):

[TABLE]

where the sum is over the different spatial modes and $\mathbf{k}$ is the wave vector. In Fig. 1 and Fig. 2 we plot the DSF of Eq. (2). For better visualization, we use an energy broadening of $0.09\,h\nu_{y}$ and $0.12\,h\nu_{y}$ for Fig. 1 and Fig. 2, respectively, similar to what was done in Ref. Blakie et al. (2002).

.2 Defining the mode character

Within the Bogoliubov theory and in the linear regime, the effect of the population of the mode $l$ on the global state dynamics can be studied using the following expression Pitaevskii and Stringari (2016)

[TABLE]

where the density fluctuations $\delta\rho_{l}=(u_{l}+v_{l}^{*})\left|\psi_{0}\right|$ and phase fluctuations $\delta\varphi_{l}=(u_{l}-v_{l}^{*})/\left|\psi_{0}\right|$ have been separated.

In order to evaluate the dominant character of each mode $l$, we introduce the quantity $C$. As discussed in the main text, the crystal and phase mode differentiate from each other by the spatial region where $\delta\varphi_{l}$ varies the most. For crystal modes, this is inside the density peaks, resulting e. g.  in a center-of-mass motion of one individual peak, which leads to a change of the crystal structure. Differently, for phase modes, $\delta\varphi_{l}$ changes the most between neighboring peaks, signalizing a particle exchange between peaks and thus a modification of the atom numbers in the peaks. We quantify these two types of character by computing the spatial variance of $\delta\varphi_{l}(\mathbf{r})$ inside the density peaks, $V_{\textrm{in}}$, and in between them, $V_{\textrm{out}}$. The quantities $V_{\textrm{in}}$ and $V_{\textrm{out}}$ are defined as follow.

For a given axial density cut of the GS wave function $|\psi_{0}(0,y,0)|^{2}$, we first define the region inside (between) the density peaks by identifying the different density maxima (minima) and number them by $j\in[1,N_{\textrm{in(out)}}]$. In a next step, we compute the mean distance $d$ between all density minima to their neighbouring maxima. Finally, we isolate the region $R_{j}=[-d/3,+d/3]$ of space centered around each maxima (minima) and calculate:

[TABLE]

The mean $\langle\cdot\cdot\cdot\rangle_{R_{j}}$ is defined for a generic function $f$ as

[TABLE]

The mode character is then evaluated by considering the ratio $C=\nicefrac{{V_{\textrm{in}}}}{{V_{\textrm{{out}}}}}$. $C$ is large for modes with prevalent crystal and small for the ones with dominant phase character. In Fig. 2 (d) we encode the information on $C$ as a color scale on the DSF spectrum. The same color map is used to illustrate the modes of the panels (a-c) in Fig. 2, confirming their correct assignment. For completeness, we also illustrate in Fig. S1 the modes’ character on the spectrum of a 164Dysupersolid, using the parameters of Fig. 1 (e) of the main text.

.3 Applying the principal component analysis to our data

.3.1 Dataset for applying the PCA

To identify the excited modes from our experimental data, we apply a general statistical method called principal component analysis (PCA) Jolliffe (2002); Segal et al. (2010); Dubessy et al. (2014) to a set of measured density distributions after a time-of-flight expansion. For our trap-excitation measurement, a dataset for the PCA is composed as follow. For each target value of $a_{\rm s}$, we record the time evolution of the density distribution for holding time, $t_{\rm h}$, between 0 and 30 ms. For each $t_{\rm h}$, we record between 15 and 30 repeated images, all together yielding a dataset of $N_{m}\gtrsim 200$ images. Each experimental run $i$ yields a two-dimensional density distribution $n_{i}(k_{x},k_{y})$. By performing a simple two dimensional Gaussian fit, we extract $71\times 71$ pixels region-of-interest (ROI) centered on the atomic cloud (the pixel’s width in $k_{x,y}$ is $0.32\mu{\rm m}^{-1}$). In addition, we post-select the shots in which the atom number, the axial cloud size and the transverse cloud size vary by less than 20%, 30% and 15% than their mean values, respectively.

.3.2 PCA’s working principle

To apply the PCA, we represent each ROI of a dataset as a vector $\rho_{i}(s)$ where $s$ represent the index of the pixel ($s\in[1,N_{p}]$, $N_{p}$ is the number of pixels in one image). We compute the mean vector image $\bar{\rho}(s)=\sum_{i=1}^{N_{m}}\rho_{i}(s)/N_{m}$ and consider the variations of the pixel values in each vector image compared to $\bar{\rho}$, $\delta\rho_{i}(s)=\rho_{i}(s)-\bar{\rho}(s)$. Finally, we consider the covariance matrix of these variations $Cov(p,s)=\sum_{i=1}^{N_{m}}\rho_{i}(s)\rho_{i}(p)/(N_{m}-1)$, which is real symmetric. By simply diagonalizing the covariance matrix, the PCA constructs a new basis of $N_{p}$ vector-images, called principal components (PCs) and written $\mathcal{C}_{p}(s)$ ($p\in[1,N_{p}]$) in the original pixel basis, that are uncorrelated one from an other. The PCs satisfy $Cov\mathcal{C}_{p}=\lambda_{p}\mathcal{C}_{p}$ where $\lambda_{p}$ is the eigenvalue of the covariance matrix associated to the PC $p$. The original vector images can be all rewritten in this new basis as $\rho_{i}(s)=\bar{\rho}(s)+\sum_{p=1}^{N_{s}}w_{p,i}\mathcal{C}_{p}(s)$, where $w_{p,i}=\sum_{s=1}^{N_{p}}\mathcal{C}_{p}(s)\rho_{i}(s)$ is the weight of the component $p$. We note that, by converting back the pixel representation to the original two-dimensional momentum space, the above decomposition means

[TABLE]

where $\mathcal{C}_{p}(k_{x},k_{y})$ encompasses now the density-distribution change induced by the PC $p$. The fact that the covariance matrix is diagonal in the PC basis indicates that the PCs correspond to uncorrelated sources of variations in the dataset. More explicitly, the coefficients $w_{p,i}$ show no correlations in between different $p$. This feature makes the PCA a powerful tool, e. g. to identify and discriminate between elementary modes of different frequencies when applied to time evolution data, as used in Ref. Dubessy et al. (2014). An example of the obtained two leading PCs in the supersolid region is given in Fig. S2.

.3.3 Identifying the elementary modes of a quantum gas via the PCA

We quickly remind the working principle, of the identification of modes via the PCA. In the linear regime, the contribution of each mode to density oscillations is expected to decouple and separate temporal and spatial variations as:

[TABLE]

with $\phi_{l}$ an arbitrary phase for the mode $l$. This relation should also hold for the density distribution after the gas’s free-expansion. If one considers that the image index $i$ encloses a time dependence ($t_{i}$), the equations (3) and (4) have a very similar structure, associating $\mathcal{C}_{p}(k_{x},k_{y})$ and $w_{p,i}$ to $\rho_{l}(\mathbf{r})$ and $\cos\left(\omega_{l}t_{i}\right)$, respectively. Thus the PCA-based identification of uncorrelated components in the time-evolution of the density profiles should enable to identify the elementary modes of the system. The corresponding PCs’ weights are then expected to oscillate in time at the frequency $\omega_{l}$ of the modes. In particular, the PCA should separate the modes oscillating at different frequencies and differentiate them from other sources of fluctuations or of dynamics (e.g. dissipation). Following Ref. Dubessy et al. (2014), we note that modes can be properly distinguished if the period associated to their beating is smaller than the total time for which the time-evolution is recorded, or, even for shorter probe time, if they have different enough amplitudes of oscillations (i.e. excitation probability).

From our dataset with repeated realizations of each hold time $t_{\rm h}$, we thus consider, for each PC $p$, the mean weights at time $t_{\rm h}$, $W_{p}(t_{\rm h})=\sum_{i/t_{i}=t_{\rm h}}w_{p,i}/\sum_{i/t_{i}=t_{\rm h}}1$. We then fit $W_{p}(t_{\rm h})$ to a sine function $A_{0}+A_{s}\cos\left(\omega t_{\rm h}+\phi\right)$ and extract the PC’s frequency ($\omega$) and amplitude $A_{s}$ of oscillation. We then consider as relevant the PCs that show oscillation of amplitude $A_{s}>8\times 10^{-4}$, frequency $\nu>20\,$Hz, and where the oscillation frequency can be extracted with a precision $<10\%$. Examples of the time evolution of $W_{p}$ and of their fits are shown in Fig. 3 (b).

We note that the PCA does not always assemble in a unique PC all the correlations in the pixel values that follow the same time dependence, and a single mode can be artificially split into several components in the analysis process. To better understand this behavior, we performed tests on theoretical calculations and compare them to the experiments. Theoretically, we specifically populate a single Bogoliubov mode on top of the ground-state, we then compute the interference patterns as a function of the hold time $t_{\rm h}$, similar to what is done in Fig. 3 (e-f), and finally we apply the PCA. For each mode considered, both for regular (BEC) or density-modulated ground-states, several leading PCs are found to oscillate. Their frequencies match the mode frequency while their oscillation amplitude decreases with the PC’s index. Typically the ratio in the oscillation amplitudes between the first and the second PC is about 10, and the amplitude of the larger-index PCs are negligible. Therefore, in the cases where the modes are the most strongly excited, i.e. mainly in the BEC regime (see $R_{l}$ scaling in Fig. 4 for the excitation amplitude), one can indeed expect that several PCs are sensitive to a single mode in experiment, matching our observation.

From those theory tests, we can also better understand the origin of this artificial splitting of one mode in several PCs. Indeed, the oscillations of the different PCs are found to have the same frequencies but different phases, typically shifted by about $\pi/2$. As it treats the pixels independently, the PCA gets confused by such $\pi/2$ phase shifts in the oscillations occurring in different regions of space, i.e. pixels’ values that distinctly oscillate, starting from their extremal or medial values. The PCA then artificially splits the oscillations occurring in these different regions into several components while they correspond to the same mode. Finally, this effect can be further favored in typically imperfect experimental settings, by the addition of experimental noise as well as other technical (e.g. imaging artifacts) or physical (e.g. dissipation) effects, which yield differences in the pixel values. Our theory tests, however, show that those additions are not the main reasons for the observed splitting.

Based on the conclusions of those tests, in the experiment (see discussion of Fig. 4 of the main text), we interpret as probing distinct modes only the PCs showing different frequencies, while PCs whose frequencies match within their error bars are interpreted as probing a single elementary excitation of the system.

.3.4 Partial recomposition

To isolate the effect of each PC on the complex time-evolution of the interference patterns, we use partial recomposition of the images inspired from Eq. (3). In particular we define

[TABLE]

This is equivalent to consider that a single PC is ”excited”, similarly to what can be done in theory for the individual excited modes of the BdG spectrum (see Fig. 2) and its description in the main text and Supp. Mat.). In Fig. 3 (c-d), we show examples of the axial cuts of $n^{(p)}(k_{x},k_{y},t)$ for two of the leading PCs. We note that here, as well as for all experimental data shown in this manuscript, the axial cuts correspond to the average of the density distributions for $|k_{x}|<1.6\,\mu{\rm m}^{-1}$.