The Fate of the Higgs Mode in a Trapped Dipolar Supersolid

TL;DR

This paper theoretically explores the excitation spectrum of a trapped dipolar quantum gas during the BEC-supersolid phase transition, revealing distinct Higgs and Nambu-Goldstone modes and their evolution.

Contribution

It provides a detailed analysis of the Higgs mode behavior in a finite-size dipolar supersolid, connecting symmetry breaking to excitation spectra.

Findings

Identification of Higgs and Nambu-Goldstone modes at the phase transition

Rapid increase in Higgs mode energy on the supersolid side

Strong coupling of Higgs mode to higher-lying modes

Abstract

We theoretically investigate the spectrum of elementary excitations of a trapped dipolar quantum gas across the BEC-supersolid phase transition. Our calculations reveal the existence of distinct Higgs and Nambu-Goldstone modes that emerge from the softening roton modes of the dipolar BEC at the phase transition point. On the supersolid side of the transition, the energy of the Higgs mode increases rapidly, leading to a strong coupling to higher-lying modes. Our study highlights how the symmetry-breaking nature of the supersolid state translates to finite-size systems.

Equations

$$i\hslash {\partial }_t\psi =H_{\mathrm{GP}}\psi ,$$

$$\left(\begin{array}{cc}{H}_{\mathrm{GP}}-\mu +\hat{X}& \hat{X}\\ -\hat{X}& -(H_{\mathrm{GP}}-\mu +\hat{X})\end{array}\right)\left(\begin{array}{c}u\\ v\end{array}\right)=\hslash \omega \left(\begin{array}{c}u\\ v\end{array}\right),$$

$$\rho (x,t)=|{\psi }_0(x){|}^2+2\lambda f(x){\psi }_0(x)\cos (\omega t),$$

$$\mathcal{I}[f]=\frac{{\left|⟨f⟩f^{\mathrm{Higgs}}\right|}^2}{⟨f⟩f⟨f^{\mathrm{Higgs}}⟩f^{\mathrm{Higgs}}},$$

$${\Phi }_{\mathrm{dip}}(\mathbf{r})=\int {\mathrm{d}}^3{r}^{\prime }{V}_{\mathrm{dd}}(\mathbf{r}-{\mathbf{r}}^{\prime })|\psi ({\mathbf{r}}^{\prime }){|}^2,$$

$$V_{\mathrm{dd}}(\mathbf{r})=\frac{3g_{\mathrm{dd}}}{4\pi }\frac{1-3{\cos }^2\vartheta }{|\mathbf{r}{|}^3},$$

$$\hat{X}\chi (\mathbf{r})={\psi }_0(\mathbf{r})\int {\mathrm{d}}^3{r}^{\prime }[V_{\mathrm{dd}}(\mathbf{r}-{\mathbf{r}}^{\prime })+g\delta (\mathbf{r}-{\mathbf{r}}^{\prime })]{\psi }_0^{\ast }({\mathbf{r}}^{\prime })\chi ({\mathbf{r}}^{\prime })+\frac{3}{2}{g}_{\mathrm{qf}}|{\psi }_0(\mathbf{r}){|}^3\chi (\mathbf{r})$$

$$\left(\begin{array}{cc}{H}_{\mathrm{GP}}-\mu +\hat{X}& -\hat{X}\\ \hat{X}& -(H_{\mathrm{GP}}-\mu +\hat{X})\end{array}\right)$$

$$(H_{\mathrm{GP}}-\mu )(H_{\mathrm{GP}}-\mu +2X)f={\hslash }^2{\omega }^{2}f$$

$$(H_{\mathrm{GP}}-\mu +2X)(H_{\mathrm{GP}}-\mu )g={\hslash }^2{\omega }^{2}g.$$

$$\mathcal{P}=\{\begin{array}{ll}1,& \underset{x}{\max }q(x)>1-ϵ\\ 0,& \underset{x}{\max }q(x)<ϵ\end{array},$$

$$S(\mathbf{k},\omega )=\sum _j{\left|\int {\mathrm{d}}^{3}r(u_j^{\ast }(\mathbf{r})+v_j^{\ast }(\mathbf{r}))e^{i\mathbf{k}\cdot \mathbf{r}}{\psi }_0(\mathbf{r})\right|}^2\delta (\omega -{\omega }_j).$$

$$\tilde{S}(\mathbf{k},\omega )=\sum _{a_{\mathrm{s}}=94.00a_0}^{a_{\mathrm{s}}=100a_0}\mathcal{I}[f_{a_{\mathrm{s}}}](\omega )S_{a_{\mathrm{s}}}(\mathbf{k},\omega ),$$

$${\mathbf{X}}_j=\mathbf{M}+\sum _{k=1}^{N_t}{Y}_{kj}{\mathbf{A}}_k,j=1,...,N_t.$$

$${\mathbf{X}}_j\approx \mathbf{M}+\sum _{k=1}^{\tilde{N}}{Y}_{kj}{\mathbf{A}}_k,j=1,...,N_t.$$

Full text

Superfluids that spontaneously break $U(1)$ symmetries occur in many fields of physics ranging from condensed matter to high energy physics to cosmology. Under certain conditions, such superfluids can undergo an additional quantum phase transition, breaking a continuous translational symmetry leading to the counter-intuitive state of matter known as supersolid Leggett (1970); Prokof’ev (2007); Boninsegni and Prokof’Ev (2012). While the experimental proof of this state in helium has remained elusive Chan et al. (2013), it has been induced as a density modulation via the invocation with light fields in spin-orbit coupled Bose-Einstein condensates (BECs) Li et al. (2017) and BECs coupled to crossed optical cavities Léonard et al. (2017a, b). In contrast to these experiments, recent work on arrays of dipolar quantum droplets Roccuzzo and Ancilotto (2019); Tanzi et al. (2019a); Böttcher et al. (2019); Chomaz et al. (2019); Tanzi et al. (2019b); Guo et al. (2019); Natale et al. (2019) has shown that the supersolid state of matter can arise purely from internal interactions of a quantum many-body system. In general, the additional breaking of the continuous translational symmetry gives rise to two collective modes: a phonon like Nambu-Goldstone mode as well as an amplitude Higgs mode Nambu (2015); Pekker and Varma (2015). In this letter, we study the behavior of the Higgs mode in a trapped system at the transition from the superfluid to a supersolid phase.

Starting with the pioneering work on the behavior of the amplitude Higgs mode in condensed matter systems Sooryakumar and Klein (1981); Littlewood and Varma (1982), the study of this collective mode has been a subject of intense theoretical Huber et al. (2007, 2008); Huber and Blatter (2009); Podolsky et al. (2011); Bruun (2014) and experimental Bissbort et al. (2011); Endres et al. (2012); Sherman et al. (2015); Leonhard (2016); Léonard et al. (2017b); Behrle et al. (2018) investigation. In contrast to particle physics, where Lorentz invariance and particle-hole symmetry provide a stable Higgs particle Higgs (1964), this collective mode is in general strongly damped due to its coupling to Goldstone modes. So far its observation required systems, where the effective low energy action exhibits an at least approximate particle-hole symmetry, yielding Lorentz-invariant equations of motion and a decoupling of the amplitude Higgs mode from the Nambu-Goldstone modes Huber et al. (2008); Pekker and Varma (2015); Leonhard (2016).

We start our analysis by describing the expected behavior of the collective low energy modes across the superfluid to supersolid quantum phase transition in the thermodynamic limit. In the superfluid regime, we expect a single Goldstone mode characterized by the superfluid sound velocity. Then, the quantum phase transition is triggered by the softening of the roton, which eventually vanishes at the transition; at this critical point, the system exhibits three gapless modes. In the supersolid phase, two gapless modes remain and correspond to the coupled Goldstone mode between the superfluid phase mode associated with the $U(1)$ symmetry breaking and the phonon mode associated with the translational symmetry breaking. In turn, the amplitude Higgs mode of the translational symmetry breaking becomes gapped and strongly damped due to the absence of particle-hole symmetry.

The focus of this manuscript is on the modification of this scenario in trapped systems, where modes and their energies are discretized due to the finite size. We map out the collective excitation spectrum across the superfluid-supersolid phase transition in a trapped system as experimentally realized with dipolar atomic gases Tanzi et al. (2019a); Böttcher et al. (2019); Chomaz et al. (2019); Tanzi et al. (2019b); Guo et al. (2019); Natale et al. (2019) and identify how both the Goldstone and Higgs modes emerge. The analysis is based on the numerical investigation of the Gross-Pitaevskii equation (GPE) and the Bogoliubov-de Gennes (BdG) equations, which is a well established method to describe the quantum phase transition as well as the low-lying collective modes. We identify a characteristic behavior of the Nambu-Goldstone mode and the amplitude Higgs mode close to the phase transition. Remarkably, the amplitude Higgs mode is stable close to the quantum phase transition even in absence of particle-hole symmetry due to the energetic splitting of the modes in finite-size systems of trapped geometries. Away from the critical point, we find that the increasing gap of the amplitude Higgs mode, then leads to a strong hybridization and avoided level crossings, which is the corresponding phenomenon to the strong damping of the amplitude Higgs mode in the thermodynamic limit. We emphasize that this characteristic fingerprint of the collective modes can be used as a tool to identify the precise position of the quantum phase transition in trapped systems and will provide new incentives to experimentally discover the Higgs mode in a system, where it arises purely from intrinsic interactions.

In view of the recent discovery of supersolidity in ultracold dipolar atomic gases Roccuzzo and Ancilotto (2019); Tanzi et al. (2019a); Böttcher et al. (2019); Chomaz et al. (2019); Tanzi et al. (2019b); Guo et al. (2019); Natale et al. (2019), we consider a system of $N=30\,000$ ${}^{162}\mathrm{Dy}$ atoms, in an harmonic potential with trapping frequencies $\omega_{\mathrm{trap}}=2\pi\,(30,\,90,\,110)\,\mathrm{Hz}$, where the magnetic field is aligned along the $y$-direction. It is the same system we considered in our previous work Guo et al. (2019), where we investigated in detail the low-energy Goldstone mode and proved the system’s supersolidity, including phase rigidity.

Such a system is described by the extended Gross-Pitaevskii equation (eGPE) Ronen et al. ; Wenzel et al. (2017); Roccuzzo and Ancilotto (2019); sup

[TABLE]

where we define $H_{\mathrm{GP}}\coloneqq H_{0}+g|\psi|^{2}+\Phi_{\mathrm{dip}}[\psi]+g_{\mathrm{qf}}|\psi|^{3}$ and $\psi$ is normalized to the atom number $N=\int\mathrm{d}^{3}r\,|\psi(\bm{r})|^{2}$. $H_{0}=\hbar^{2}\nabla^{2}/2m+V_{\mathrm{ext}}$ contains the kinetic energy and trap confinement $V_{\mathrm{ext}}(\bm{r})=m(\omega_{x}^{2}x^{2}+\omega_{y}^{2}y^{2}+\omega_{z}^{2}z^{2})/2$. The contact interaction strength $g=4\pi\hbar^{2}a_{\mathrm{s}}/m$ is given by the scattering length $a_{\mathrm{s}}$ and the dipolar mean field potential $\Phi_{\mathrm{dip}}$ by the dipolar length $a_{\mathrm{dd}}=\mu_{0}\mu_{\mathrm{m}}^{2}m/(12\pi\hbar^{2})$ ($a_{\mathrm{dd}}\approx 130\,a_{0}\ \text{for }^{162}\mathrm{Dy})$, where $\mu_{\mathrm{m}}$ is the magnetic moment. Within the local density approximation, quantum fluctuations are taken into account by the term $g_{\mathrm{qf}}|\psi|^{3}$, which is the Lee-Huang-Yang (LHY) correction Lee et al. (1957) for dipolar systems Lima and Pelster (2011, 2012), necessary to prevent the dipolar collapse Petrov (2015); Kadau et al. (2015). First, we perform imaginary time evolution of Eq. (1) to find the ground state $\psi_{0}$. We obtain, as a function of the scattering length a regular BEC, a supersolid or an array of isolated droplets, as schematically shown at the top of Fig. 1.

In order to study the elementary excitations, we follow Refs. Morgan et al. (1998); Huepe et al. (2003); Ronen et al. (2007); Baillie et al. (2017); Chomaz et al. (2018); Roccuzzo and Ancilotto (2019) and shortly review the BdG equations here. We linearize the wavefunction $\psi(\bm{r},t)=[\psi_{0}(\bm{r})+\lambda(u(\bm{r})e^{-i\omega t}+v^{*}(\bm{r})e^{i\omega t})]e^{-i\mu t/\hbar}$ around its ground state $\psi_{0}$ with the Bogoliubov amplitudes $u$ and $v$, oscillating at a frequency $\omega/2\pi$ as a small pertubation with amplitude $\lambda$. $\mu$ is the chemical potential and $u,\ v$ are normalized by $\int\mathrm{d}^{3}r\,(u(\bm{r})^{2}-v(\bm{r})^{2})=1$. Substituting the expression for $\psi$ in Eq. 1 and retaining terms up to linear order in $\lambda$ yields a system of linear equations that can be neatly expressed in matrix form

[TABLE]

where the operator $\hat{X}$ Baillie et al. (2017); Chomaz et al. (2018) is defined by its action on either function $\chi=u,\,v$ in the supplementary material sup . $\hat{X}$ contains the contact and dipolar interactions, as well as quantum fluctuations. It is essential to note here that our system intrinsically requires all three contributions, as a BEC with only contact and dipolar interactions in our parameter region of scattering lengths would collapse. Only through the stabilization with beyond mean-field effects given by the LHY correction, the droplet phase can exist Petrov (2015); Kadau et al. (2015); Ferrier-Barbut et al. (2016).

By solving Eq. (2) for different scattering lengths across the BEC-supersolid phase transition, we gain insight into the energetic and spatial character of the elementary excitations in the self-organized supersolid regime. The excitation energies as a function of the scattering length are shown in Fig. 1. Here we focus on the low-lying modes, in particular the eight lowest, which are decoupled from the $y$- and $z$-direction and mainly feature an excitation along the droplet array. It is instructive to label the energy branches by the corresponding mode’s parity as even (odd), if the corresponding function that describes the density fluctuation is a symmetric (anti-symmetric) function along the direction of the droplet array sup . The mode parity is essential in understanding the spectrum. As a function of one interaction parameter (in our case the scattering length $a_{\mathrm{s}}$), energy levels with modes of opposite parity are allowed to cross each other, while energy levels with modes of equal parity hybridize, leading to avoided crossings and level repulsion von Neumann and Wigner (1929); Landau and Lifshitz (1981).

To distinguish the superfluid and supersolid phases, we look at the droplet density overlap, as defined in our previous works Böttcher et al. (2019); Guo et al. (2019) through the ratio of density in between central and side droplets and peak density in the central droplet, which is shown as the blue background in Fig. 1 (noted as min/max ratio). One can also think of this quantity as an inverse measure of the density modulation strength, which is connected to the superfluid fraction. In Fig. 1, we indicate the phase transition by a dashed line for clarity. Let us now begin by considering the BEC phase at $a_{\mathrm{s}}=100\,a_{0}$ on the right hand side of Fig. 1 and move in our discussion systematically towards smaller scattering lengths. The lowest branch is the anti-symmetric dipole mode, located at $30\,\mathrm{Hz}$, exactly the trap frequency. It is not affected by a change of scattering length, which is an instance of Kohn’s theorem Khon (1961) that demands the dipole mode to be always completely decoupled from all interactions, and thus from all other modes. This can serve as a test of the numerical accuracy of our calculations Ronen et al. , and our simulation recovers the dipole mode exactly for all considered scattering lengths. The next higher mode is the lowest quadrupole mode, which is followed by other modes, that have been studied for non-dipolar Stringari (1996) and dipolar O’Dell et al. (2004); van Bijnen et al. (2010) systems in various trap geometries. We find that the parity of the modes is alternating between odd and even as a function of the excitation energy on the BEC side, which energetically separates the modes well, as the kinetic energy increases in alternating parity.

Moving towards smaller scattering lengths, one can recognize clearly that a bundle of branches quickly decreases in energy (see $a_{\mathrm{s}}\in(97.0,98.5)\,a_{0}$ in Fig. 1). It consists of the modes induced by the roton softening Giovanazzi and O’Dell (2004); Corson et al. (2013); Chomaz et al. (2018) and is always composed of one symmetric and one anti-symmetric mode, which are degenerate until the scattering length approaches the phase transition (around $a_{\mathrm{s}}=96.8\,a_{0}$). We note that the roton instability has previously been predominantly studied in 2D- systems San (2003); Ronen et al. (2007); Wilson et al. (2008); Jona-Lasinio et al. (2013); Natu et al. (2014); Shen and Quader (2018) and recently in 1D-systems Petter et al. (2019); Chomaz et al. (2018). One can clearly see the anti-crossings between modes of equal parity as the roton mode softens.

Further towards smaller scattering lengths, around $a_{\mathrm{s}}=96.8\,a_{0}$, the system approaches the superfluid-supersolid phase transition, where these two modes continuously split in energy, losing their degeneracy. This can be intuitively understood, as, in contrast to the normal modes in the BEC that are defined by the confinement, the roton is defined by the intrinsic interactions and is localized only in the center of the BEC, where the density is almost flat. Even in our finite system, we remarkably observe that this leads to a degeneracy of the even and odd roton modes. Crossing the phase transition, a density modulation emerges in the ground state and one is no longer able to freely choose the position at which to place the modulation. This lifts the degeneracy of even and odd modes that match their wavelength to the spacing of the droplets. Thus, we emphasize that the significant change at the phase transition is the energetic splitting of the roton mode into the Goldstone and Higgs modes of the supersolid. In fact the splitting point of the roton bundle into Higgs and Goldstone branches may serve as a rigorous definition for the transition point from the superfluid to a supersolid state. It is interesting to note that this splitting point coincides exactly with the simple measure of overlap between the droplets that has been used in recent works Böttcher et al. (2019); Guo et al. (2019), and is shown as the blue background in Fig. 1. It underlines the validity of this simple measure to characterize the quantum critical point, as it was done in Böttcher et al. (2019); Guo et al. (2019); see also Zhang et al. (2019).

This splitting of the roton bundle on the BEC side into the Higgs and the Goldstone mode of the supersolid can be seen as the main feature in Fig. 1. The even mode that increases its energy with decreasing scattering length features an oscillation between the droplet crystal and the superfluid background and therefore can be identified as the Higgs mode (as illustrated further below, see Fig. 2(b)). On the other hand, the decreasing mode is the low-energy Goldstone mode that has been recently described in Guo et al. (2019). Despite the drastic change in energy, we observe that the mode pattern of these modes does not strongly change throughout the phase transition sup .

The Higgs mode in an infinite system is theoretically expected to be gapped, with a vanishing gap only exactly at the quantum critical point Huber et al. (2007, 2008); Podolsky et al. (2011); Pekker and Varma (2015). Due to the finite size of the here considered trapped system, the Higgs gap at the phase transition is lifted to a finite frequency of $\Delta_{\mathrm{min}}=16.97\,\mathrm{Hz}$ at $a_{\mathrm{s}}=96.75\,a_{0}$. In systems closer to the thermodynamic limit, such as with higher atom number and lower longitudinal confinement we expect the minimum Higgs gap at the quantum critical point to decrease, consistent with the predictions of a vanishing gap in the thermodynamic limit Huber et al. (2007, 2008); Huber and Blatter (2009); Podolsky et al. (2011); Pollet and Prokof’Ev (2012).

Contrary to the gapped Higgs mode, the Goldstone mode in an infinite system is always gapless Huber et al. (2007, 2008); Huber and Blatter (2009); Podolsky et al. (2011); Bruun (2014). This again translates to a finite frequency of the Goldstone mode in our finite-size system, seen as the lowest anti-symmetric mode in Fig 1. As the superfluid fraction decreases towards smaller scattering lengths, the frequency approaches zero, consistent with the decrease of the associated sound velocity of this mode. This is when the first symmetric mode also reaches energies below the dipole mode energy. We can identify this branch as the higher-$k$ Goldstone mode, with the lowest asymmetric branch being the $k=0$-Goldstone mode (labeled ”Goldstone” in Fig. 1).

The phase transition from BEC to supersolid breaks the translational symmetry, in addition to the already broken $U(1)$ symmetry, which is essential for the existence of the Higgs (amplitude) mode Huber et al. (2008). The broken translational symmetry means that the system is no longer Galilean invariant, thus we expect the Higgs mode to be overdamped Bissbort et al. (2011); Pollet and Prokof’Ev (2012); Podolsky and Sachdev (2012); Pekker and Varma (2015). We illustrate this further in Fig. 2, where we show the coupling strength of the Higgs mode to other modes and provide a view into the spatial modulation pattern. However, exactly at the phase transition, in a narrow range of scattering lengths $a_{\mathrm{s}}\in(96.45,96.65)\,a_{0}$ the Higgs mode exists in an isolated state and is not coupled to any other mode, which is remarkable. This is a finite-size feature, since in the present system size, the energetic spacing between the lowest and other symmetric modes, especially at the quantum critical point, is significant.

To visualize how the Higgs mode couples to the other modes, we display the time evolution of these modes near the coupling points, as shown in Fig. 2(b-g). The density time evolution Ronen et al. ; Martin and Blakie (2012) is given by

[TABLE]

where $f(x)=u(x)+v(x)$ and $\psi_{0}(x)$ are 1D-cuts in the $y=0$ plane at $z=0$ through the density fluctuation $f$ and the ground state $\psi_{0}$, respectively. Defining a measure on the strength of the Higgs character within each mode will allow us to follow its coupling and decay into the other symmetric modes more quantitatively. Here, we use the density mode pattern $f^{\mathrm{Higgs}}(x)$ at the scattering length $a_{\mathrm{s}}=96.65\,a_{0}$ (Fig. 2(b)) as a reference, since at this scattering length, the droplets have clearly developed, yet the energy of the Higgs branch is furthest from other symmetric modes, thus minimizing its coupling. We then calculate its fidelity (the modulus squared overlap)

[TABLE]

for all mode patterns $f$ and map it to a color code on the points shown in Fig. 1, which results in the signature of Fig. 2(a).

The time evolution of the Higgs density pattern (Fig. 2(b)) shows that this mode is characterized by an oscillating exchange of amplitude between crystal and fluid, which is in line with the illustration of an amplitude oscillation with respect to the order parameter in the valley of the famous Mexican-hat potential Podolsky et al. (2011); Endres et al. (2012); Pekker and Varma (2015); Guo et al. (2019); Léonard et al. (2017b, a). For slightly lower scattering lengths, we observe that the branch on which we identify the most typical Higgs character now decreases in energy and develops an additional breathing character while maintaining part of the BEC background oscillation character, see Fig. 2(c). As this occurs, the two higher lying symmetric modes increase in energy and develop the characteristic amplitude fluctuation of the Higgs mode, see Fig. 2(d, f). In combination with an increase of the Higgs overlap (4) on these branches, we can identify that the actual Higgs branch steadily increases in energy away from the quantum critical point after its minimum has been surpassed, schematically indicated by the red line in Fig. 2(a). Effectively, this means that the Higgs branch hybridizes with the other symmetric modes, as has been noted in a previous work considering strongly interacting lattice Bosons Bissbort et al. (2011), and as a consequence decays into the other symmetric modes due to its overdamped nature Bissbort et al. (2011); Pollet and Prokof’Ev (2012); Podolsky and Sachdev (2012); Pekker and Varma (2015). Fig. 2 shows remarkably well the coupling and decay of the Higgs mode to the higher symmetric modes, indicated by the green color, but also the loss of Higgs character within each branch. Following along the lowest symmetric branch in Fig. 2 towards smaller scattering lengths, the strength of the Higgs character within this branch decays rapidly, as it develops a more dominant breathing character, see Fig. 2(c). With a further decrease of scattering length, the mode develops into the higher-$k$ Goldstone mode, that manifests as an oscillation of atoms between the central and the two side droplets sup . We additionally show points of strongest Higgs character within the two higher symmetric branches in Fig. 2(d, f) and emergence of breathing character and decay along the branches in Fig. 2(e, g). Towards the isolated droplet regime at even smaller scattering lengths, the Higgs mode has disappeared, as its character decays with the exponentially decreasing overlap between the droplets. We can relate this to the results obtained in two-dimensional superfluids Pollet and Prokof’Ev (2012); Endres et al. (2012); Bruun (2014), where the Higgs mode can only be seen close to the quantum critical point and its character disappears towards the Mott-insulating state of matter.

In view of the existence of an isolated Higgs branch at the quantum critical point and its traceability into other symmetric modes within a range in scattering length of about $1\,a_{0}$, we propose the Higgs mode to be experimentally observable with current experimental setups Chomaz et al. (2018); Guo et al. (2019); Böttcher et al. (2019). To this end, we have performed real time evolution simulations of the eGPE (1) using the standard experimental protocol of starting in the BEC regime with a subsequent linear ramp of the scattering length down into the supersolid state of matter. We systematically verified sup that in these simulations, the strongest density response is found to resemble the Higgs density pattern from our BdG calculations as together with a weaker density response that resembles the breathing patterns (as in the mixed state Fig. 2(c)), proving the possibility of the experimental observation of the Higgs mode in the future.

In conclusion, we have mapped out the spectrum of collective excitations across the BEC-supersolid phase transition and identified the Goldstone and Higgs branches that emerge from the two roton modes. We thereby improve the understanding of the Goldstone mode, as measured in Guo et al. (2019), and predict the existence of an isolated and stable Higgs mode in a narrow range of interaction strengths. The Higgs mode is gapped, its energy increases further away from the transition and it strongly couples to higher lying modes. Let us emphasize again: here we considered a system, exclusively subject to its own intrinsic interactions and a trapping potential. By varying only one of the intrinsic interaction parameters, the system undergoes a supefluid-supersolid phase transition. In the supersolid, the confinement enhances energetic splitting between Higgs and other modes, which leads to an isolated and decoupled Higgs mode. Due to the clear framework of the here considered system, it provides an excellent opportunity to study a correspondence between models in high energy physics Arkani-Hamed et al. (1998); Randall and Sundrum (1999) and condensed matter systems.

As a next step, we plan to experimentally confirm the here presented findings and map out the coupling of the Higgs mode to higher lying modes.

I Supplementary material

I.1 Simulation details

In this section, we quickly summarize simulation details not further explained in the main text.

Our system is described by the eGPE Eq. (1). The dipolar mean field potential is given by

[TABLE]

where

[TABLE]

is the dipolar interaction with the dipolar interaction parameter $g_{\mathrm{dd}}=4\pi\hbar^{2}a_{d}d/m$ given by the dipolar length $a_{\mathrm{dd}}=\mu_{0}\mu_{\mathrm{m}}^{2}m/(12\pi\hbar^{2})$ defined by the magnetic moment $\mu_{\mathrm{m}}$ and $\vartheta$ is the angle between $\bm{r}$ and the magnetic field axis.

Furthermore, $g_{\mathrm{qf}}=32/(3\sqrt{\pi})ga^{3/2}Q_{5}(\epsilon_{\mathrm{dd}})$ is taking into account quantum fluctuations within the local density approximation for dipolar systems Lima and Pelster (2011, 2012), where $\epsilon_{\mathrm{dd}}=g_{\mathrm{dd}}/g=a_{\mathrm{dd}}/a_{\mathrm{s}}$ is the relative dipolar strength. In our simulations, we use the approximation $Q_{5}(\epsilon_{\mathrm{dd}})=1+\frac{3}{2}\epsilon_{\mathrm{dd}}^{2}$ Lima and Pelster (2012); Ferrier-Barbut et al. (2016); Wenzel et al. (2017); Bisset et al. (2016); Baillie et al. (2017). The mean field dipolar potential Eq. (5) is effectively calculated using a Fourier transform, where we employ a spherical cutoff for the dipolar potential Eq. (6), and choose its radius such that there is no spurious interaction between periodic images Góral and Santos (2002); Ronen et al. ; Lu et al. (2010). To find the ground state in all simulations, we use the method of imaginary time propagation.

In order to study the elementary excitations, we employ the BdG formalism introduced in the main text. The coupling operator $\hat{X}$ in the BdG matrix representation Eq. (2) of the main text is given by Baillie et al. (2017); Chomaz et al. (2018)

[TABLE]

Note that we follow the sign convention ($+u$, $+v$) Morgan et al. (1998); Ronen et al. in the expansion of $\psi$, instead of ($+u$, $-v$) Baillie et al. (2017); Chomaz et al. (2018), which inverts the sign on the off-diagonal elements of the BdG matrix Eq. (2). The particle (hole) components are usually, in our convention, identified by the functions $u$ ($v$) Ronen et al. . If we set $v\rightarrow-v$, we can convert the matrix representation used in our convention ($+u$, $+v$) to the matrix

[TABLE]

previously used in the ($+u$, $-v$) convention Baillie et al. (2017); Chomaz et al. (2018).

In practice, we employ the same method described in Ronen et al. , which is to rewrite $f=u+v$ and $g=v-u$, corresponding to the density and phase fluctuations, respectively, which enables the reduction of Eq. (2) after a short calculation into the equations

[TABLE]

and

[TABLE]

We solve the equation for $f$ and calculate the respective $g$ as in Ronen et al. . With the solutions for $f$ and $g$, $u=(f-g)/2$ and $v=(f+g)/2$ can be recovered. The solution for the equation of $f$ is obtained by using the Arnoldi method to calculate the lowest eigenvalues and eigenfunctions Huepe et al. (2003); Ronen et al. ; Martin and Blakie (2012). Note that the BdG equations (2) as presented in the main text always have the solution $g=\psi_{0}$ with $\omega=0$, as one can easily verify from Eq. (9). This solution is sometimes referred to as neutral mode Morgan et al. (1998); Huepe et al. (2003); Ronen et al. and has its origin in fixing the phase of the wavefunction by introducing $\mu$ in the linearization. However, we emphasize that this is not a physical degree of freedom and refrain from calling it a ”mode” since $(H_{\mathrm{GP}}-\mu)g=0$ for $g=\psi_{0}$ is merely stating that the ground state fulfills the eGPE, which is redundant.

We point out that the conventional BdG description, as we present it in this paper, is linearizing the the wavefunction around its ground state with quasiparticles that are not orthogonal to the condensate Morgan et al. (1998). In a previous study on nonlinear mixing of these quasiparticles Morgan et al. (1998), it was shown that transfer of the quasiparticle population from one mode into another is inversely dependent on the the inital finite amount of excitation. The mixing processes are enhanced by a finite temperature and dominated by energy conserving conversion processes Morgan et al. (1998). This will provide an additional theoretical challenge in the study of the coupling of the Higgs mode and its dependence with considerations of finite-temperature effects. In particular, a detailed study of the Higgs mode in the coupled regime and at higher energies may require an orthogonal quasiparticle basis.

We classify the excitation modes in terms of the parity of the density response. For this, we choose $\epsilon=0.1$ as the criterion that determines the parity,

[TABLE]

where $q(x)=|f(x)-f(-x)|/(2\,\underset{x}{\max}\ |f(x)|)$, $\mathcal{P}=1,\,0$ for odd or even modes, respectively.

I.2 Elementary excitation spectra

Since we have a finite sized system, we obtain a discrete spectrum of elementary excitations for every scattering length, rather than a continuous one for an infinite system. In order to visualize the momentum dependence of the obtained modes (which we will label by an index $j$ in the following), we calculate the zero-temperature dynamic strucure factor Brunello et al. (2001); Blakie et al. (2002), given by

[TABLE]

It is also known as the density correlator Huber et al. (2007), typically measured as a response function for Bragg spectroscopy, since it correlates and Fourier transforms the density fluctuation $f=u+v$ with the ground state $\psi_{0}$ as $\mathcal{F}[f^{*}\cdot\psi_{0}](\bm{k})$, where $\mathcal{F}$ is the Fourier transform. For illustration purposes in Fig. S1, we convolve (11) with a Gaussian along the $\omega$-axis.

In Fig. S1, we show 4 different scattering lengths, $a_{\mathrm{s}}=\{97.50,\,96.65,\,96.45,\,95.50\}\,a_{0}$ corresponding to a BEC state close to the supersolid phase transition, a state that displays the isolated Higgs mode, a state that displays the mixed Higgs mode, and a state further towards the isolated droplet regime, respectively. The left hand panels show the dynamic structure factor Eq. (11), and the right hand side shows Higgs character $\mathcal{I}$ (Eq. 4 in the main text), of the modes at the correponding energies. In Fig. S1(a), we see a discrete version of the typical dispersion relation of dipolar BECs at an interaction strength that features the roton minimum. As the roton further softens towards smaller scattering lengths, the system approaches the phase transition and finally, as the quantum critical point has been crossed, splits into the Goldstone and Higgs modes, see Fig. S1(b). The lowest mode is the low-energy Goldstone mode that we briefly discussed in the main text and was the main subject of our previous work Guo et al. (2019). The symmetric mode above the Goldstone mode is the isolated Higgs mode, energetically well separated from other symmetric modes which all display a very small Higgs character $\mathcal{I}$, as one can see in the side panel of the spectrum. In contrast, at $a_{\mathrm{s}}=96.45\,a_{0}$ (Fig. S1(c)), as the Higgs mode rises in energy, it hybridizes with the higher-lying symmetric modes and mixes its character with them, resulting in a decreased Higgs strength on the lowest symmetric mode and an increased Higgs strength on the higher symmetric modes. Even further towards the isolated droplet regime (Fig. S1(d)), the Higgs mode has disappeared (there is only residual overlap $<30\,\%$ with other symmetric modes).

We emphasize that the $k_{x}$ in Fig. S1(b-c) is to be understood by folding the roton momentum $k_{\mathrm{rot}}$ back to zero, since after the density modulation has emerged, the dispersion relation of the system is best described within first Brillouin zone $k_{x}\in[0,k_{\mathrm{rot}})$. One can see in Fig. S1(b) that the splitting of the roton mode into Goldstone and Higgs modes takes place at $k_{\mathrm{rot}}$ as well, so that the prominent isolated Higgs mode we mainly discuss in the main text is the $k=0$-Higgs mode. We observed the spectra such as the ones shown in Fig. S1 for all scattering lengths and checked for the $k$-dependence of the Higgs mode. In order to do so, we weighted every mode’s structure factor with its Higgs overlap, and superimposed the weighted structure factors of all scattering lengths,

[TABLE]

in order to filter out modes with smaller and keep modes with larger Higgs character. We find that the maximum strength of $\tilde{S}$ shifts slightly away from $k_{\mathrm{rot}}$, but exactly for the scattering lengths at which the Higgs mode hybridizes with the breathing modes. In this sense one cannot distinguish anymore whether $\tilde{S}$ probes partly the varying $k$ of the breathing modes or not.

I.3 Experimental observability of the Higgs mode

In order to verify that the isolated Higgs mode should be experimentally observable, we performed real time evolution simulations of the eGPE and followed a standard experimental protocol. Starting from a BEC ground state at $a_{\mathrm{s}}=112\,a_{0}$, we linearly ramp the scattering length within $30\,\mathrm{ms}$ down into the supersolid regime at scattering lengths $a_{\mathrm{s}}=96.65,\ 96.45\,a_{0}$ and observe the subsequent time evolution up to $300\,\mathrm{ms}$.

For later comparison, we show the BdG mode patterns corresponding to the pure Higgs and breathing density variation functions $f$ in Fig. S2. These are the 1st and 3rd lowest symmetric modes of Fig. 1 in the main text, at $a_{\mathrm{s}}=96.65\,a_{0}$. The Higgs pattern corresponds to the time evolution shown in Fig. 2(a). The breathing pattern is similar to Fig. 2(c), but without a mixed Higgs character.

By performing the experimental scheme for different final scattering lengths, we obtain the full density evolution shown in Figs. S3, S4 for the scattering length at which we expect to find an isolated Higgs mode, or a mode that displays mixed Higgs and breathing character, respectively. After the $30\,\mathrm{ms}$ ramp and another approximate $15\,\mathrm{ms}$, the droplets form and begin an oscillatory motion.

In order to decompose this oscillatory motion into the elementary excitations, as we have studied them in the BdG formalism, we employ a model-free statistical analysis, the principal component analysis (PCA) Jolliffe (2002); Segal et al. (2010); Bertram et al. (2014); Dubessy et al. (2014); Ferrier-Barbut et al. (2018); Cao et al. (2019). In short, PCA provides a way to map a timeseries of density distributions to a set of density variation distributions, whose weight varies in time. PCA has a wide range of general applications Jolliffe (2002), from image analysis, to dimensional reduction of large datasets, but it turns out as well that there is a direct correspondence Dubessy et al. (2014); Barr (2015) between the variation distributions, or principal components (PCs), and the density variation function $f$ obtained by the BdG formalism, so that we can identify a specific PC with a specific BdG mode and we may see the weight of this mode as a function of time. One of PCA’s properties is that the signal (in our case the timeseries shown in Fig. S3(a), S4(a)) can be reconstructed exactly using a superposition of all PCs, oscillating at their respective weights. However, a very small subset out of all PCs accounts for most of the density variation (or information) contained within the timeseries. This is essentially how PCA is used for dimensional reduction. More explicitly stated, it means the following Segal et al. (2010); Dubessy et al. (2014). Consider a set of $N_{t}$ images of $p$ pixels each. Typically for high resolution images or simulation grids, $N_{t}<p$, which we will assume in the following. We represent the $j$’th image by a vector $\bm{X}_{j}$, the mean image by $\bm{M}=\sum_{j=1}^{N_{t}}\bm{X}_{j}/N_{t}$, write the dataset in a $p\times N_{t}$ matrix $B=[\bm{X}_{1}-\bm{M},\ ...\ ,\bm{X}_{N_{t}}-\bm{M}]$ and apply PCA to this dataset. PCA diagonalizes the covariance matrix of $B$ and thus essentially finds a lower dimensional linear transform $Y$ (an $N_{t}\times N_{t}$ matrix) that represents the data $\{\bm{X}_{j}\}_{j=1}^{N_{t}}$ by a new set of $N_{t}$ orthonormal basis vectors $\{\bm{A}_{j}\}_{j=1}^{N_{t}}$ (the PCs), which are pairwise uncorrelated, sorted in decreasing variance and yield the original data exactly as

[TABLE]

PCA is a total-variance-preserving operation, meaning that the total variance of the signal is captured exactly by the sum of the variances of the individual PCs. But since we ordered the PCs in decreasing variance, it means that we can account for most of the variance, considering only a small subset $\{\bm{A}_{j}\}_{j=1}^{\tilde{N}}$ with $\tilde{N}\ll N$, of all PCs. We normalize the variance of each PC by the total variance as a relative measure of total variance accounted for. Reconstruction of the original dataset with only partial information is then written as

[TABLE]

In the simulation, we take around $N_{t}=1500$ snapshots of the wavefunction and we find that merely considering $\tilde{N}=2$ PCs, as we show further below, accounts for more than $90\,\%$ of the total variation.

From Eqs. (13, 14), we see that the weight of each PC for the $j$’th image is contained in the $j$’th column of $Y$. Hence, in the following we will label the 1D-cuts along the $x$-axis of the $k$’th PC as $\mathcal{A}_{k}(x)$, corresponding to $\bm{A}_{k}$ and its weight as a function of time as $\mathcal{Y}_{k}(t)$, corresponding to $\{Y_{kj}\}_{j=1}^{N_{t}}$.

In the comparison of the PCs shown in Figs. S3, S4, the similarity between the BdG Higgs and breathing patterns shown in Fig. S2 is very clear, so that we can identify that the Higgs mode has been clearly, and since it is the first PC, most strongly, excited for both values of the final scattering length $a_{\mathrm{s}}=96.65,\ 96.45\,a_{0}$. In both cases we find the breathing pattern to be second most strongly excited and again we identify the strong similarity between the BdG mode pattern and the PCs.

However, if we observe the weights of the PCs as a function of time, we find for the scattering length $a_{\mathrm{s}}=96.65\,a_{0}$, where we expect to find the isolated Higgs mode, that both the weights of Higgs and breathing patterns (Fig. S3(c, e)) are oscillating at a single frequency of $17.09\pm 0.01\,\mathrm{Hz}$ and $40.46\pm 0.01\,\mathrm{Hz}$, respectively. The frequencies of the PCs match the BdG results of $21.82\,\mathrm{Hz}$ and $40.50\,\mathrm{Hz}$ at this scattering length quite well, which, in addition to the very similar shape of the modes density variation distributions, underlines that the modes have been correctly identified and that they are excited while being decoupled from each other. The slight difference of the frequencies for the Higgs mode might be due to different numerical grids that were used to carry out the BdG calculation and the full real time evolution. Additionally, we observe a clear envelope of the breathing mode’s weight in Fig. S3(e), which may indicate that with the ramping procedure, even higher lying modes have been excited. Over the timespan of $300\,\mathrm{ms}$, these higher lying modes decay into the lower lying breathing mode, thus increasing its weight in the timeseries. In contrast, the amplitude of the Higgs weight in Fig. S3(c) is approximately constant, indicating that neither the population of this mode increases (no quasiparticle decays into the Higgs mode) nor that it decreases (no quasiparticle is excited from the Higgs to higher lying modes). This can be understood from Fig. 1(a), since at $a_{\mathrm{s}}=96.65\,a_{0}$, the Higgs branch is energetically well separated from the more densely distributed, higher lying symmetric modes. Thus, it is less likely to exchange the quasiparticle occupation with any of the higher symmetric modes.

In contrast, for $a_{\mathrm{s}}=96.45\,a_{0}$, where we expect to find the coupled Higgs mode, there are multiple frequency components in the signal, which we can see again from the weights (Fig. S4(c, e)). This means that, though we again have clearly excited both the Higgs and breathing mode, they are now coupled to other modes. As a result, the identification of the Higgs mode in this regime is complicated by the fact that despite the main density variation in a timeseries such as Fig. S4(a) being described by the two PCs in Fig. S4(b, d), their weight’s frequencies can not be assigned to a single elementary exciation anymore.

Another way to identify the excitation of a specific mode in the experimental protocol, is to reconstruct the signal only using specific subsets of PCs, which must not necessarily be the ones with highest variance, and compare it to the time evolution one obtains from the BdG equations, as in Fig. 2(b-g). E.g. if one is interested in whether 2 arbitrary PCs $m_{1},\,m_{2}$ correspond to a specific (or a superposition of) BdG mode(s), one may study the expression $\bm{M}+Y_{{m_{1}}j}\bm{A}_{m_{1}}+Y_{{m_{2}}j}\bm{A}_{m_{2}}$ for $j=1,\ ...,\ N_{t}$ instead of Eq. 14, where $m_{1},\,m_{2}\neq 1,\,2$. Calculating a measure on how close to the BdG time evolution (Eq. 3 in the main text) is to the reconstructed signal using the $m_{1}$ and $m_{2}$ PCs can yield better correspondence between BdG and PCA for more complex PCs.

As a final remark, ramping the scattering length further into the isolated droplet regime, $a_{\mathrm{s}}\lesssim 96\,a_{0}$, the first PCs do not indicate the spatial shape of the Higgs or breathing pattern as shown in Fig. S3, S4(b,d) anymore, but rather a collection of other highly complex breathing modes. This means that, as expected from the discussion of the decaying Higgs character towards the isolated droplets in the main text, the Higgs mode is not excited by this protocol, or decays into other symmetric modes on timescales that are not accessible in our simulation and in current experimental setups.