Dynamical symmetry and breathers in a two-dimensional Bose gas

TL;DR

This study explores the dynamical symmetry and breather phenomena in a two-dimensional Bose gas, demonstrating how scale invariance influences its evolution and revealing novel periodic breather states.

Contribution

The paper experimentally investigates the dynamical symmetry in a 2D Bose gas and introduces a new type of breather solution in the Gross--Pitaevskii equation.

Findings

Scale invariance relates different evolutions via a scaling transform.

Certain initial shapes lead to periodic breather states.

Experimental validation of Lorentz group symmetry in 2D Bose gases.

Abstract

A fluid is said to be \emph{scale-invariant} when its interaction and kinetic energies have the same scaling in a dilation operation. In association with the more general conformal invariance, scale invariance provides a dynamical symmetry which has profound consequences both on the equilibrium properties of the fluid and its time evolution. Here we investigate experimentally the far-from-equilibrium dynamics of a cold two-dimensional rubidium Bose gas. We operate in the regime where the gas is accurately described by a classical field obeying the Gross--Pitaevskii equation, and thus possesses a dynamical symmetry described by the Lorentz group SO(2,1). With the further simplification provided by superfluid hydrodynamics, we show how to relate the evolutions observed for different initial sizes, atom numbers, trap frequencies and interaction parameters by a scaling transform. Finally we show that some specific initial shapes - uniformly-filled triangles or disks - may lead to a periodic evolution, corresponding to a novel type of breather for the two-dimensional Gross--Pitaevskii equation.

Equations

$$\mathbf{r}\to \mathbf{r}/\lambda ,t\to t/{\lambda }^2.$$

$$E_{\mathrm{i}}=\frac{N^2{\hslash }^2}{2m}\tilde{g}\int |\psi (\mathbf{r}){|}^4{\mathrm{d}}^{2}r,$$

$$E_{\mathrm{p}}=\frac{N}{2}m{\omega }^2\int r^2|\psi (\mathbf{r}){|}^2{\mathrm{d}}^{2}r.$$

$$E_{\mathrm{k}}=\frac{N{\hslash }^2}{2m}\int |\nabla \psi {|}^2{\mathrm{d}}^{2}r.$$

$$\frac{\mathrm{d}{E}_{\mathrm{p}}}{\mathrm{d}t}$$

$$\frac{\mathrm{d}W}{\mathrm{d}t}$$

$$E_{\mathrm{p}}(t)=\frac{1}{2}{E}_{\mathrm{tot}}+\Delta E\cos (2\omega t),$$

$$\mathrm{i}\hslash \frac{\partial {\psi }_1}{\partial t}=-\frac{{\hslash }^2}{2m}{\nabla }^2{\psi }_1+\frac{{\hslash }^2\tilde{g}N}{m}|{\psi }_1{|}^2{\psi }_1+\frac{1}{2}m{\omega }_1^2{\mathbf{r}}^2{\psi }_1.$$

$${\psi }_2({\mathbf{r}}^{\prime },t^{\prime })=f(\mathbf{r},t){\psi }_1(\mathbf{r},t)$$

$$\lambda (t)={\left[\frac{1}{{\alpha }^2}{\cos }^2({\omega }_{1}t)+{\alpha }^2{\zeta }^2{\sin }^2({\omega }_{1}t)\right]}^{1/2},$$

$$\tan ({\omega }_2{t}^{\prime })=\zeta {\alpha }^2\tan ({\omega }_{1}t),$$

$$f(\mathbf{r},t)=\lambda (t)\exp \left(-\mathrm{i}\frac{m\dot{\lambda }{\mathbf{r}}^2}{2\hslash \lambda }\right),$$

$$\tan ({\omega }_2{t}^{\prime })={\omega }_{2}t,$$

$$\lambda (t)=(1+{\omega }_2^2{t}^2{)}^{1/2}.$$

$$(n_1|n_2)=\int n_1(\mathbf{r})n_2(\mathbf{r})d^{2}r$$

$$p[n_1,n_2,\lambda ]=\frac{(n_1^{(\lambda )}|n_2)}{||n_1^{(\lambda )}||||n_2||}$$

$$\mathcal{O}[n_1,n_2]=\underset{\lambda }{\max }p[n_1,n_2,\lambda ].$$

$${\partial }_{t}n+\nabla \cdot (n\mathbf{v})=0,$$

$$m{\partial }_t\mathbf{v}+\nabla \left(\frac{1}{2}m{\mathbf{v}}^2+\frac{{\hslash }^2}{m}\tilde{g}n+\frac{1}{2}m{\omega }^2{r}^2+P(n)\right)=0,$$

$$m{\partial }_t\mathbf{v}+\nabla \left(\frac{1}{2}m{\mathbf{v}}^2+\frac{{\hslash }^2}{m}\tilde{g}n+\frac{1}{2}m{\omega }^2{r}^2\right)=0.$$

$${\tilde{g}}_2{N}_2={\mu }^2{\tilde{g}}_1{N}_1,L_2=\alpha L_1,{\omega }_2=\zeta {\omega }_1,$$

$${\tilde{g}}_2{n}_2({\mathbf{r}}^{\prime },t^{\prime })$$

$${\mathbf{v}}_2({\mathbf{r}}^{\prime },t^{\prime })$$

$$\lambda (t)={\left[\frac{1}{{\alpha }^2}{\cos }^2({\omega }_{1}t)+{\left(\frac{\zeta \alpha }{\mu }\right)}^2{\sin }^2({\omega }_{1}t)\right]}^{1/2},$$

$$\tan ({\omega }_2{t}^{\prime })=\frac{\zeta {\alpha }^2}{\mu }\tan ({\omega }_{1}t).$$

$${\tilde{g}}_2{n}_2(\mu \mathbf{r},t_{\pi /2}^{\prime })={\tilde{g}}_1{n}_1\left(\mathbf{r},t_{\pi /2}\right).$$

$$\mathbf{r}\to \mathbf{r},t\to t+\beta ,$$

$$\mathbf{r}\to \mathbf{r}/\lambda ,t\to t/{\lambda }^2,$$

$$\mathbf{r}\to \frac{\mathbf{r}}{\gamma t+1},\mathbf{t}\to \frac{t}{\gamma t+1},$$

$$\mathbf{r}\to g(\mathbf{r},t)\equiv \frac{\mathbf{r}}{\gamma t+\delta }$$

$$\mathbf{t}\to h(t)\equiv \frac{\alpha t+\beta }{\gamma t+\delta },$$

$${\mathcal{P}}_0[{\psi }_1;\mathbf{r},t]=0$$

Full text

I Introduction

Symmetries play a central role in the investigation of a physical system. Most often they are at the origin of conserved quantities, which considerably simplify the study of the equilibrium states and the evolution of the system. For example, spatial symmetries associated with translation and rotation lead to the conservation of linear and angular momentum. More generally, it is interesting to determine the dynamical (or hidden) symmetries of the system under study, which can lead to more subtle conserved quantities. These symmetries are described by the group of all transformations of space and time that leave the action, therefore the equations of motion, invariant. A celebrated example is the $1/r$ potential in three dimensions, where there exists a dynamical symmetry described by the group O(4) for the bounded orbits Bander:1966 . When treated by classical mechanics, it leads to the conservation of the Laplace–Runge–Lenz vector, from which one deduces that the bounded orbits are actually closed trajectories.

Among the systems that display rich dynamical symmetries are the ones whose action is left invariant by a dilation transformation of space and time. Such scale-invariant systems were initially introduced in particle physics to explain scaling laws in high energy collisions Jackiw:1972 . We will consider here the non-relativistic version of scale-invariance, which applies to the dynamics of a fluid of particles. We consider the simultaneous change of length and time coordinates of each particle according to the scaling:

[TABLE]

In this dilation, the velocity of a particle is changed as $\boldsymbol{v}\to\lambda\boldsymbol{v}$. Therefore the kinetic energy of the fluid scales as $E_{\rm k}\to\lambda^{2}E_{\rm k}$, which ensures that the corresponding part of the action ($\propto\int E_{\rm k}\;{\rm d}t$) remains invariant in the transformation (1). If the interaction energy has the same scaling, $E_{\rm i}\to\lambda^{2}E_{\rm i}$, the total action of the fluid is invariant in the dilation. The simplest example of such a fluid is a collection of non-relativistic particles, either non-interacting ($E_{\rm i}=0$) or with pair-wise interactions described by a $1/r^{2}$ potential. A scale-invariant fluid possesses remarkably simple thermodynamic properties: for example its equation of state depends only on the ratio of chemical potential to temperature, instead of being an independent function of these two variables.

Most physical systems exhibiting scale invariance also possess a more general conformal invariance, where time and space are modified by conformal transformations instead of the simple dilations given in Eq.(1) Nakayama:2013 . In the non-relativistic domain, this conformal invariance exists for the Schrödinger equation describing the motion of the two systems mentioned above, free particles Hagen:1972 ; Niederer:1972 and particles interacting with a $1/r^{2}$ potential Alfaro:1976 . In both cases, the dynamical symmetry group associated to this scale/conformal invariance is the Lorentz group SO(2,1). This is also the case for the three-dimensional pseudo-spin $1/2$ Fermi gas in the unitary regime (for a review, see e.g. Zwerger:2011 ). There, the scattering length between the two components diverges, ensuring the required disappearance of a length scale related to interactions. In addition to the existence of a universal equation of state, this dynamical symmetry leads to a vanishing bulk viscosity Son:2007 ; Elliott:2014 , and also to general relations between the moments of the total energy and those of the trapping energy in a harmonic potential Werner:2006 .

In this article we consider another example of a scale/conformal invariant fluid with the SO(2,1) dynamical symmetry, the "weakly interacting" two-dimensional (2D) Bose gas. The concept of "weak interaction" means in this context that the state of the gas is well described by a classical field $\psi(\boldsymbol{r},t)$. This field is normalized to unity ($\int|\psi|^{2}\;{\rm d}^{2}r=1$), so that the density of the gas reads $n(\boldsymbol{r},t)=N|\psi(\boldsymbol{r},t)|^{2}$ where $N$ is the number of particles. In the scaling of positions, the 2D matter-wave field changes as $\psi(\boldsymbol{r})\to\lambda\,\psi(\lambda\boldsymbol{r})$, which guarantees that the norm is preserved and that the dynamical part of the action $\propto{\rm i}\hbar\int{\rm d}t\int{\rm d}^{2}r\;\psi^{*}\,\partial_{t}\psi$ is invariant. The interaction energy of the gas then reads for contact interaction

[TABLE]

where $m$ is the mass of a particle and $\tilde{g}$ the dimensionless parameter characterizing the strength of the interaction. One can immediately check that $E_{\rm i}$ obeys the $\lambda^{2}$ scaling required for scale invariance, which can be viewed as a consequence of the dimensionless character of $\tilde{g}$. The classical field description used here is valid if one restricts to the case of a small coupling strength $\tilde{g}\ll 1$ Svistunov:2015 . This restriction is necessary because of the singularity of the contact interaction $\frac{\hbar^{2}}{m}\tilde{g}\delta(\boldsymbol{r})$ in 2D when it is treated at the level of quantum field theory. Note that the condition $\tilde{g}\ll 1$ does not constrain the relative values of the interaction and kinetic energies. Actually in the following we will often consider situations where $E_{\rm i}\gg E_{\rm k}$ (Thomas–Fermi regime).

So far, the scale/conformal invariance of the weakly interacting 2D Bose gas has been mainly exploited to measure its equation of state Hung:2011 ; Yefsah:2011 . Also, one of its dynamical consequences in an isotropic 2D harmonic potential of frequency $\omega$ has been explored: the frequency of the breathing mode was predicted to be exactly equal to $2\omega$ for any $\tilde{g}$ Kagan:1996b ; Pitaevskii:1997a ; gritsev:2010 , as tested in Refs.  Chevy:2002 ; Vogt:2012 . Note that in the presence of a harmonic potential, the whole system is not scale-invariant anymore, but it still possesses a dynamical symmetry described by the group SO(2,1), as shown in Pitaevskii:1997a . Recently, deviations from this prediction for $\tilde{g}\gtrsim 1$, an example of a quantum anomaly Olshanii:2010 , have been observed Peppler:2018 ; Holten:2018 .

The purpose of our work is to go beyond static properties of the weakly interacting 2D Bose gas and its single-mode oscillation in a harmonic potential, and to reveal more general features associated with its dynamical symmetry. To do so, we study the evolution of the gas in a 2D harmonic potential of frequency $\omega$, starting from a uniformly-filled simple area (disk, triangle or square). Here we use $\tilde{g}\leq 0.16$ so that the classical field description is legitimate. We first check (Sec. II) the prediction from Pitaevskii:1997a that $E_{\rm k}+E_{\rm i}$ should have a periodic evolution in the trap with the frequency $2\omega$. We then investigate the transformations linking different solutions of the equations of motion. These transformations are at the heart of the dynamical symmetry group SO(2,1). In practice we first link the evolution of clouds with the same atom number and homothetic initial wave functions in harmonic potentials with different frequencies (Sec. III). Then, restricting to the case where superfluid hydrodynamics is valid, we derive and test a larger family of transformations that allows one to connect the evolutions of two initial clouds of similar shapes with different sizes, atoms numbers, trap frequencies and interaction strengths (Sec. IV). Finally in Sec. V we explore a property that goes beyond the symmetry group of the system and that is specific to triangular and disk-shaped distributions in the hydrodynamic limit: we find numerically that these distributions evolve in a periodic manner in the harmonic trap, and we confirm this prediction over the accessible range for our experiment (typically two full periods of the trap $4\pi/\omega$). These particular shapes can therefore be viewed as two-dimensional breathers for the Gross–Pitaevskii (non-linear Schrödinger) equation in the hydrodynamic limit dauxois2006physics . They also constitute a novel example of universal dynamics in a quantum system prepared far from equilibrium Eigen:2018aa ; Ern:2018 ; Prufer:2018aa .

II Evolution of potential energy

Our experiment starts with a 3D Bose-Einstein condensate of 87Rb that we load around a single node of a vertical ($z$) standing wave created with a laser of wavelength 532 nm. The confining potential along $z$ is approximately harmonic with a frequency $\omega_{z}/(2\pi)$ up to 4.9 kHz. The interaction parameter is $\tilde{g}=\sqrt{8\pi}\,a_{s}/\ell_{z}$, where $a_{s}$ is the 3D s-wave scattering length and $\ell_{z}=(\hbar/m\omega_{z})^{1/2}$. The interaction energy per particle and the residual temperature are both smaller than $\hbar\omega_{z}$ so that the vertical degree of freedom is effectively frozen Ville:2017 . The initial confinement in the horizontal $xy$ plane is ensured by "hard walls" made with a light beam also at 532 nm. This beam is shaped using a digital micromirror device (DMD), and a high-resolution optical system images the DMD pattern onto the atomic plane Aidelsburger:2017b , creating a box potential on the atoms. The cloud fills uniformly this box potential, and it is evaporatively cooled by adjusting the height of the walls of the box. For all data presented here, we kept the temperature low enough to operate deep in the superfluid regime, $T/T_{c}<0.3$, where $T_{c}$ is the critical temperature for the Berezinskii-Kosterlitz-Thouless transition. At this stage the atoms are prepared in the $F=1,m_{F}=0$ hyperfine (ground) state, which is insensitive to magnetic field.

Once the gas has reached equilibrium in the 2D box, we suddenly switch off the confinement in the $xy$ plane and simultaneously transfer the atoms to the field-sensitive state $F=1,m_{F}=-1$ using two consecutive microwave transitions, via the intermediate state $F=2,m_{F}=0$. Most of the experiments are performed in the presence of a magnetic field that provides the internal state $F=1,m_{F}=-1$ with an isotropic harmonic confinement in the $xy$ plane, with $\omega/2\pi$ around $19.5$ Hz. We estimate the anisotropy of the potential to be $\lesssim 2\%$. We let the cloud evolve in the harmonic potential for an adjustable time before making an in-situ measurement of the spatial density $n(\boldsymbol{r})=N|\psi(\boldsymbol{r})|^{2}$ by absorption imaging.

The measurement of $n(\boldsymbol{r})$ gives access to both the interaction energy (2) and the potential energy in the harmonic trap

[TABLE]

Since the gas is an isolated system, we expect the total energy $E_{\rm tot}=E_{\rm k}+E_{\rm i}+E_{\rm p}$ to be conserved during the evolution, where the kinetic energy $E_{\rm k}$ reads:

[TABLE]

The SO(2,1) symmetry for a 2D harmonically trapped gas brings a remarkable result: $E_{\rm k}+E_{\rm i}$ and $E_{\rm p}$ should oscillate sinusoidally at frequency $2\omega$ Pitaevskii:1997a . More precisely, using the 2D Gross–Pitaevskii equation one obtains the relations

[TABLE]

where we have defined $W=\omega m\int\boldsymbol{r}\cdot\boldsymbol{v}\,n\ {\rm d}^{2}\boldsymbol{r}$, and the velocity field $\boldsymbol{v}(\boldsymbol{r})=\frac{\hbar}{m}\mbox{Im}\left[\psi^{*}(\boldsymbol{r})\boldsymbol{\nabla}\psi(\boldsymbol{r})\right]/|\psi(\boldsymbol{r})|^{2}$. Initially the gas is prepared in a steady state in the box potential, so that $\boldsymbol{v}=0$ hence $W(0)$ is also null. Therefore the potential energy evolves as

[TABLE]

where $\Delta E=\frac{1}{2}[E_{\rm p}(0)-E_{\rm k}(0)-E_{\rm i}(0)]$ can be positive or negative. A similar result holds for the sum $E_{\rm k}+E_{\rm i}$ (with $\Delta E$ replaced by $-\Delta E$), but not for the individual energies $E_{\rm k}$ or $E_{\rm i}$.

We show in Fig. 1(a) the evolution of the potential energy per particle for an initially uniformly-filled square. Although the density distribution is not periodic (see Fig. 1(b)), the potential energy $E_{\rm p}$ evolves periodically and is well fitted by a cosine function with a period that matches the $2\omega$ prediction and the expected zero initial phase. For a better adjustment of the data, we added a (small) negative linear function to the fitting cosine. Its role is likely to account for the residual evaporation rate of atoms from the trap ($\sim 0.1\$\mathrm{s}$^{-1}$).

This simple dynamics can be viewed as a generalisation of the existence of the undamped breathing mode at frequency $2\omega$ that we mentioned in the introduction Kagan:1996b ; Pitaevskii:1997a . We emphasize that this result is a consequence of the SO(2,1) symmetry and would not hold for the Gross–Pitaevskii equation in 1D or 3D.

III General scaling laws

An important consequence of the dynamical symmetry of the 2D Gross–Pitaevskii equation is the ability to link two solutions $\psi_{1,2}$ of this equation corresponding to homothetic initial conditions: one can relate $\psi_{1}(\boldsymbol{r},t)$ and $\psi_{2}(\boldsymbol{r}^{\prime},t^{\prime})$, provided they evolve with the same parameter $\tilde{g}N$ and the same trap frequency $\omega_{1}=\omega_{2}$. By using a simple scaling on space and time, this link can be further extended to the case $\omega_{1}\neq\omega_{2}$.

The general procedure is presented in Appendix A and we start this section by summarizing the main results. Consider a solution of the Gross–Pitaevskii equation $\psi_{1}(\boldsymbol{r},t)$ for the harmonic potential of frequency $\omega_{1}$:

[TABLE]

Using scale/conformal invariance, we can construct a solution $\psi_{2}(\boldsymbol{r}^{\prime},t^{\prime})$ of the Gross–Pitaevskii equation with the frequency $\omega_{2}=\zeta\omega_{1}$ using:

[TABLE]

where space is rescaled by $\boldsymbol{r}^{\prime}={\boldsymbol{r}}/{\lambda(t)}$ with

[TABLE]

where the dimensionless parameter $\alpha$ is the homothetic ratio between the initial states. The relation between the times $t$ and $t^{\prime}$ in frames 1 and 2 is

[TABLE]

and the multiplicative function $f$ is

[TABLE]

where $\dot{\lambda}\equiv\frac{{\rm d}\lambda}{{\rm d}t}$. The two solutions $\psi_{1,2}(t)$ correspond to the evolution of two clouds with the same parameter $\tilde{g}_{1}N_{1}=\tilde{g}_{2}N_{2}$. At $t=0$, these two wave functions correspond to the ground states of the Gross-Pitaevskii equation in the box potentials with characteristic lengths $L_{1,2}$, with $L_{2}=\alpha L_{1}$. Both initial wave functions $\psi_{1,2}(0)$ can be chosen real, and the scale invariance of the (time-independent) 2D Gross-Pitaevskii equation ensures that they are homothetic: $\alpha\psi_{2}(\alpha\boldsymbol{r},0)=\psi_{1}(\boldsymbol{r},0)$. For example in the limit $E_{\rm i}\gg E_{\rm k}$, $\psi(0)$ corresponds to a uniform density in the bulk, and goes to zero at the edges on a scale given by the healing length $\xi\equiv[N\hbar^{2}/(2mE_{\rm i})]^{1/2}$. For two box potentials of homothetic shapes filled with the same number of particles, the ratio $\xi_{2}/\xi_{1}$ is equal to the ratio $L_{2}/L_{1}$.

We explore experimentally this mapping between two evolutions in the particular case $L_{1}=L_{2}$ and $\omega_{1}\to 0$, i.e., $\alpha=1$ and $\zeta\to+\infty$. This corresponds to comparing the evolution of clouds with the same shape and the same size either in a harmonic potential or in free (2D) space. The choice of the initial shape is arbitrary; here we start from a uniform triangle of side length $40.2(3)\$\mathrm{\SIUnitSymbolMicro m} with $3.9(3)\times 10^{4}$ atoms and let it evolve either in a harmonic potential of frequency $\omega_{2}/(2\pi)=19.7(2)\,$\mathrm{Hz}, or without any potential ($\omega_{1}=0$). In both cases, we record images of the evolution, examples of which are given in Fig. 2(a) and 2(b). These two evolutions should be linked via Eq. (9). The relation (11) between $t$ and $t^{\prime}$ reads

[TABLE]

and the relation (10) becomes:

[TABLE]

The relation (13) indicates that the scaling transformation maps the first quarter of the oscillation period in the harmonic trap $\omega_{2}t^{\prime}\leq\pi/2$ onto the ballistic expansion from $t=0$ to $t=\infty$. In the absence of interactions, this result has a simple physical interpretation: after the ballistic expansion between $t=0$ to $t=\infty$, the asymptotic position distribution reveals the initial velocity distribution of the gas, whereas the evolution in the harmonic trap during a quarter of oscillation period exchanges initial positions and initial velocities. We emphasize that the mapping (13) also holds for an interacting system, as a consequence of the SO(2,1) symmetry underlying the Gross–Pitaevskii equation 111 One may question the validity of the Gross-Pitaevskii equation (GPE), hence of scale invariance, after a long expansion time when the gas occupies a large area $R^{2}$. If we were interested in the ground state of a box of size $R\to\infty$ and a given $\tilde{g}N$, we would indeed expect deviations with respect to GPE because the relevant momenta $k\sim R^{-1}$ would tend to 0, and logarithmic corrections in $k$ to the coupling constant would become significant Petrov:2000a . Here, this issue is absent because the initial interaction energy is converted into kinetic energy at the beginning of the expansion. The relevant atomic momenta thus remain $\sim(2mE_{\rm int}/N\hbar^{2})^{1/2}$ at all times, which validates the use of the GPE. .

In order to reconstruct the scaling laws (13) and (14) from the measured evolutions, we compare each image $n_{1}(\boldsymbol{r},t)$ for the free evolution with the set of images $n_{2}(\boldsymbol{r}^{\prime},t^{\prime})$ obtained for the in-trap evolution. More precisely, we start by defining the overlap ${\cal O}[n_{1},n_{2}]$ between two images in the following way:

• •

We introduce the scalar product $(n_{1}|n_{2})$ between two images

[TABLE]

and the norm of an image $||n_{1}||=\sqrt{(n_{1}|n_{1})}$.

• •

In order to relate two images that differ by a spatial scaling factor $\lambda$, we introduce the quantity

[TABLE]

where $n_{1}^{(\lambda)}(\boldsymbol{r})=\lambda^{2}n_{1}(\lambda\boldsymbol{r})$ is the image rescaled from $n_{1}(\boldsymbol{r})$ by the factor $\lambda$, with the same atom number: $N_{1}=\int n_{1}(\boldsymbol{r})\;{\rm d}^{2}r=\int n_{1}^{(\lambda)}(\boldsymbol{r})\;{\rm d}^{2}r$. Note that the definition of the norm given above entails $||n_{1}^{(\lambda)}||=\lambda||n_{1}||$. By construction the quantity $p[n_{1},n_{2},\lambda]$ is always smaller or equal to 1, and it is equal to 1 only when the image $n_{1}^{(\lambda)}$ is identical to $n_{2}$ up to a multiplicative factor.

• •

Finally, for a couple of images $(n_{1},n_{2})$ we vary $\lambda$ and define their overlap as

[TABLE]

In practice, for each image $n_{1}(t)$ acquired at a given time $t$, we determine the time $t^{\prime}_{\rm opt}$ where the overlap between $n_{1}(t)$ and $n_{2}(t^{\prime})$ is optimal. We denote $\Lambda(t)$ the value of the scaling parameter $\lambda$ for which the value ${\cal O}[n_{1}(t),n_{2}(t^{\prime}_{\rm opt})]$ is reached (see Supplemental material for more details). Since the center of the images may drift during the evolution, we also allow for a translation of $n_{2}$ with respect to $n_{1}$ when looking for the optimum in (16-17).

The result of this mapping between the two evolutions is shown in Fig. 2(c) and 2(d). In Fig. 2(c), we plot $t^{\prime}_{\rm opt}$ as a function of $t$. The prediction (13) is shown as a continuous line and is in good agreement with the data. In Fig. 2(d), we show the variation of the corresponding optimal scaling parameter $\Lambda(t)$. Here again the prediction (14) drawn as a continuous line is in good agreement with the data. The overlap between the density distributions at the corresponding times is shown in the inset of Fig. 2(d) and is always around $0.95$, confirming that these density distributions have very similar shapes. Indeed, the overlap between two images averaged over a few experimental realizations, taken in the same conditions ranges from $0.98$ to $0.99$ due to experimental imperfections.

Finally we note that here we connected solutions of the Gross-Pitaevskii equation (8) with the same atom number $N_{1}=N_{2}$. Actually the results derived above also apply to pairs of solutions with $\tilde{g}_{1}N_{1}=\tilde{g}_{2}N_{2}$, since only the product $\tilde{g}N$ enters in the Gross-Pitaevski equation (8).

IV Scaling laws in the hydrodynamic regime

In the previous section, we have linked the evolution of two clouds with the same atom number $N$ (or the same $\tilde{g}N$). We show now that it is also possible to link evolutions with different $N$’s and $\tilde{g}$’s, provided we restrict to the so-called hydrodynamic (or Thomas–Fermi) regime, where the healing length $\xi$ is very small compared to the size of the gas.

IV.1 General formulation

The Gross–Pitaevskii equation (8) can be equivalently written in terms of the density and the velocity fields as

[TABLE]

where $P(n)=-\hbar^{2}/2m\ (\boldsymbol{\nabla}^{2}\sqrt{n})/\sqrt{n}$ is the so-called quantum pressure. When the characteristic length scales over which the density and velocity vary are much larger than the healing length $\xi$, one can neglect the contribution of the quantum pressure in (19):

[TABLE]

This approximation, corresponding to the Thomas–Fermi limit, leads to the regime of quantum hydrodynamics for the evolution of the density $n$ and the irrotational velocity field $\boldsymbol{v}$ Pitaevskii:2016 . It enriches the dynamical symmetries of the problem, as we see in the following. For our experimental parameters, this approximation is legitimate since the healing length is a fraction of micrometer only, much smaller than the characteristic size of our clouds (tens of micrometers).

We consider two homothetic shapes, e.g. two box-like potentials with a square shape, with sizes $L_{1,2}$ and filled with $N_{1,2}$ atoms. We assume that we start in both cases with the ground state of the cloud in the corresponding shape, so that the initial velocity fields are zero. Note that contrarily to the case of Sec. III, the ratio between the healing lengths $\xi_{2}/\xi_{1}$ is not anymore equal to $L_{2}/L_{1}$ so that the initial wave functions are not exactly homothetic, but this mismatch occurs only close to the edges over the scale $\sim\xi_{1,2}\ll L_{1,2}$. As before, at time $t=0$ we switch off the potential creating the shape under study, and switch on a harmonic potential with frequency $\omega_{1,2}$. Our goal is to relate the two evolutions with parameters $(\tilde{g}_{1}N_{1},L_{1},\omega_{1})$ and $(\tilde{g}_{2}N_{2},L_{2},\omega_{2})$.

The general transformation involves three dimensionless constant parameters $\mu,\alpha,\zeta$:

[TABLE]

and reads:

[TABLE]

with $\dot{\lambda}=\frac{{\rm d}\lambda}{{\rm d}t}$. The spatial variables are rescaled as $\boldsymbol{r}^{\prime}=\boldsymbol{r}/\lambda(t)$ with the function $\lambda$ now given by

[TABLE]

and the relation between the times $t$ and $t^{\prime}$ in frames 1 and 2 is:

[TABLE]

With a calculation similar to that detailed in Appendix A, one can readily show that if $(n_{1},\boldsymbol{v}_{1})$ is a solution of the hydrodynamic equations (18,20) for the frequency $\omega_{1}$, then $(n_{2},\boldsymbol{v}_{2})$ is a solution for the frequency $\omega_{2}$. If $\mu=1$, these equations also apply beyond the Thomas–Fermi limit, as shown in Sec. III. More strikingly, they show that, in the quantum hydrodynamic regime, the evolution of any cloud is captured by a universal dynamics that depends only on its initial geometry.

IV.2 Connecting evolutions with a fixed trap frequency, a fixed size and different $\tilde{g}N$

We present here the experimental investigation of the scaling described above, focusing on the case $L_{1}=L_{2}$ and $\omega_{1}=\omega_{2}$, i.e., $\alpha=\zeta=1$. In other words, we compare the evolution of two clouds with the same initial shape and density distribution, different atom numbers and different interaction strengths in a given harmonic trap. For simplicity, we consider the result of the evolution at times $t$ and $t^{\prime}$ such that $\omega_{1}t=\omega_{2}t^{\prime}=\pi/2$, which satisfies the constraint (25). In this case $\lambda(t)=1/\mu$, so that the general scaling (22) reads

[TABLE]

We start with a cloud in a uniform box potential with the shape of an equilateral triangle of side length $L=38.2(3)\$\mathrm{\SIUnitSymbolMicro m}. At $t=0$ we transfer the atoms in the harmonic trap of frequency $\omega/2\pi=19.6\ $\mathrm{Hz} and remove the box potential. At $t=\pi/(2\omega)$ we image the cloud. We perform this experiment for different values of $\tilde{g}$ – and slightly different atom numbers – corresponding to the product $\tilde{g}N$ between 200 and $4000$. This leads to a ratio $\xi/L$ always smaller than $0.03$, ensuring that we stay in the quantum hydrodynamic regime. The variation of $\tilde{g}$ is achieved by changing the intensity $I$ of the laser beams creating the vertical confinement, with $\tilde{g}\propto I^{1/4}$. The values of $\tilde{g}$ are obtained from the measurement of the vertical frequency $\omega_{z}$ (see Supplemental Material).

We analyze the series of images using the same general method as in Sec. III. We select arbitrarily one image as a reference point (here, the one corresponding to $\tilde{g}N\approx 2000$, shown as a red square on Fig. 3). Then, we calculate the best overlap between this reference point and all other images obtained for different $\tilde{g}N$’s, and extract an optimal scaling parameter $\Lambda$. The results of this analysis are displayed on Fig. 3. The inset shows that the overlap is close to 1 for all values of $\tilde{g}N$, indicating that the clouds all have the same shape, as expected from (26). On the main graph of Fig. 3, we show the variations of $\Lambda^{-2}$ with $\tilde{g}N$. The scaling law (24) predicts that $\Lambda^{-2}=\mu^{2}\propto\tilde{g}N$, which is indicated by the solid line passing by the origin and the reference point. Here again, this prediction is in excellent agreement with the data. Note however that the result for the largest $\tilde{g}N$ is slightly lower than the theoretical prediction, which we attribute to the fact that, for larger powers in the vertical confining laser beam, the local defects of the potential that it creates start to play a significant role.

Interestingly, the shape for $t^{\prime}=\pi/(2\omega)$, i.e., $t=\infty$ for an evolution without any trap, is close to a uniformly filled triangle but inverted compared to the initial one (see insets of Fig. 3). The emergence of such a simple form after time-of-flight is reminiscent of the simple diamond-like shape obtained for the 3D expansion of a uniform gas initially confined in a cylindrical box Gotlibovych:2014 . Note that we also observe such a diamond-like shape at $t=\pi/(2\omega)$ starting from a square box, albeit with a non-uniform density (see Supplemental Material).

IV.3 Connecting evolutions with a fixed trap frequency, different sizes and different $\tilde{g}N$

Finally we compare the evolution of two clouds with homothetic shapes and $\alpha,\mu\neq 1,\zeta=1$, which means clouds with different initial sizes, different atom numbers and evolving in the same harmonic trap. We perform an experiment where the initial shape is a square with a uniform density. The first cloud has a side length $L_{1}=27.0(5)\$\mathrm{\SIUnitSymbolMicro m}, contains $N_{1}=3.7(3)\times 10^{4}$ atoms, and its initial density distribution is shown on Fig. [4](#S4.F4)(a). The second one has a side length $L_{2}=36.8(5)\ $\mathrm{\SIUnitSymbolMicro m} and contains $N_{2}=5.4(3)\times 10^{4}$ atoms (Fig. 4(b)). The ratio $\xi/L$ is around $0.01$ for these two clouds. We let them evolve in the same harmonic potential described above and with the same interaction parameter $\tilde{g}$, and take pictures after different evolution times. We expect that the two evolutions $n_{1}(\boldsymbol{r},t)$ of the first cloud and $n_{2}(\boldsymbol{r}^{\prime},t^{\prime})$ of the second cloud are linked via Eqs. (22, 24, 25), with parameters $\alpha=L_{2}/L_{1}=1.36(4)$ and $\mu=\sqrt{N_{2}/N_{1}}=1.21(8)$. We analyze the two series of images with the same procedure as in Sec. III and determine the scaling laws that link the two evolutions one to the other. The best overlaps between the images of the first and second series are shown in Fig. 4(c). They are all above 0.97, indicating that the two evolutions are indeed similar. The relation between the time $t^{\prime}$ of the second frame and the corresponding time $t$ of the first frame is shown on Fig. 4(d), and the best scaling factor $\Lambda(t)$ is shown on Fig. 4(e). The solid lines show the theoretical predictions (25) and (24), which are in very good agreement with the experimental data.

With the three experiments described in Sec. III and Sec. IV, we have tested the scaling laws (22)-(25) independently for the three parameters $\alpha$, $\mu$ and $\zeta$, demonstrating that, in the quantum hydrodynamic regime, the evolution of a cloud initially at rest depends only on its initial shape, up to scaling laws on space, time and atom density.

V Two-dimensional breathers

In Sec. II, we have seen that due to the SO(2,1) symmetry, the evolution of the potential energy $E_{\rm p}$ is periodic with period $T/2\equiv\pi/\omega$ for an arbitrary initial state $\psi(\boldsymbol{r},0)$ [see Eq. (7)]. Of course, the existence of this periodicity does not put a strong constraint on the evolution of $\psi(\boldsymbol{r},t)$ itself. Due to the non-linear character of the Gross-Pitaevskii equation, the evolution of $\psi$ is not expected to be periodic, as illustrated on Fig. 1(b) for a square initial shape. When looking experimentally or numerically at various initial shapes like uniformly-filled squares, pentagons or hexagons, we indeed observe that even though $E_{\rm p}(jT/2)=E_{\rm p}(0)$ for integer values of $j$, the shapes $n(\boldsymbol{r})=N|\psi(\boldsymbol{r})|^{2}$ at those times are notably different from the initial ones. We found two exceptions to this statement, which are the cases of an initial equilateral triangle and a disk. This section is devoted to the study of these very particular states that we call "breathers".

In the present context of a fluid described by the Gross–Pitaevskii equation, we define a breather as a wave function $\psi(\boldsymbol{r},t)$ that undergoes a periodic evolution in an isotropic harmonic trap of frequency $\omega$ (for a generalisation to different settings, see e.g. Bishop:1980 ; dauxois2006physics ). According to this definition, the simplest example of a breather is a steady-state $\psi_{\rm ss}(\boldsymbol{r})$ of the Gross–Pitaevskii equation, e.g. the ground state. Other breathers are obtained by superposing $\psi_{\rm ss}$ with one eigenmode of the Bogoliubov–de Gennes equations, resulting from the linearization of the Gross-Pitaeveskii equation around $\psi_{\rm ss}$. In principle (with the exception of the breathing mode Pitaevskii:1997a ), the population of this mode should be vanishingly small to avoid damping via non-linear mixing. Extending this scheme to the superposition of several modes in order to generate more complex types of breathers seems difficult. Indeed the eigenmode frequencies are in general non-commensurable with each other, therefore the periodicity of the motion cannot occur as soon as several modes are simultaneously excited 222For the ground state of a harmonically confined 2D gas in the Thomas–Fermi limit, the mode frequencies are $\omega\left(2n^{2}+2nm+2n+m\right)^{1/2}$ with $n,m$ positive or null integers Stringari:1998 ; Ho:1999 .. Note that for a negative interaction coefficient $\tilde{g}$ in 1D, a bright soliton forms a stable steady-state of the Gross–Pitaevskii equation (even for $\omega\to 0$) and thus also matches our definition. In that particular 1D case, a richer configuration exhibiting explicitly the required time-periodicity is the Kuznetsov–Ma breather, which is obtained by superposing a bright soliton and a constant background (see e.g. Zhao:2018 and refs. in).

Here, we are interested in 2D breathers that go well beyond a single mode excitation and we start our study with the uniform triangular shape. In this case, for experiments performed with a gas in the Thomas–Fermi regime, we find that the evolution of the shape is periodic with period $T/2$ within the precision of the measurement. As an illustration, we show in Fig. 5(a) four images taken between $t=0$ and $T/2$. The scalar product $(n(0)|n(t))$ between the initial distribution and the one measured at times $T/2$, $T$, $3T/2$ and $2T$, shown in Fig. 5(b), is indeed very close to 1. We could reproduce the same result for various initial atom numbers.

We did not find an analytical proof of this remarkable result, but we could confirm it numerically by simulating the evolution of a wave function $\psi(\boldsymbol{r},t)$ with the Gross–Pitaevskii equation 333In spite of the fact that our parameters are well in the Thomas–Fermi regime, we perform the numerical analysis using the Gross–Pitaevskii equation (8) and not the quantum hydrodynamic equations (18,20). The reason is that the discontinuity of the density that appears in the latter case on the edge of the sample may lead to numerical singularities in the subsequent dynamics.. We show in Fig. 6(a) a few snapshots of the calculated density distribution, and in Fig. 6(b) the evolution of the modulus of the (usual) scalar product $|\langle\psi(0)|\psi(t)\rangle|$ between the wave functions at times [math] and $t$. The calculation was performed on a square grid of size $N_{s}\times N_{s}$ with $N_{s}=512$. The initial wave function is the ground state of a triangular box with the side length $N_{s}/2$ centered on the grid, obtained by imaginary time evolution for $\tilde{g}N=25600$. Note that by contrast to the "scalar product between images" introduced above, the quantity $|\langle\psi(0)|\psi(t)\rangle|$ is also sensitive to phase gradients of the wave functions. Its evolution shows clear revivals approaching unity for $t$ close to multiples of $T/2$.

We show in Fig. 7(a) the finite-size scaling analysis of the value of the first maximum of this scalar product occurring at $t_{\rm max}\approx T/2$, for increasing grid sizes $N_{s}=64,\ldots,1024$. The product $\tilde{g}N$ is adjusted such that the healing length $\xi=[N\hbar^{2}/(2mE_{\rm i})]^{1/2}=a\ell$, where $\ell$ is the grid spacing and $a^{2}=0.5,1,2,4,8$. The condition $a\ll N_{s}$ ensures that $\xi$ is much smaller than the size of the triangle (Thomas–Fermi regime), while having $a\gtrsim 1$ provides an accurate sampling of the edges of the cloud. The overlap between $|\psi(0)\rangle$ and $|\psi(t_{\rm max})\rangle$ increases with the grid size, and reaches $0.995$ for the largest grid.

In the simulation, the trapping frequency $\omega$ is adjusted such that $|\Delta E|\ll E_{\rm tot}$ in Eq. (7); the cloud then keeps an approximately constant area over time, which is favourable for the numerics. Note that this choice does not restrict the generality of the result, since the scaling laws seen in Sec. III allow one to connect the evolution of a given $\psi(\boldsymbol{r},t=0)$ in traps with different frequencies. In particular, if the evolution starting from $\psi(\boldsymbol{r},0)$ in a trap of frequency $\omega_{1}$ is periodic with period $\pi/\omega_{1}$, the evolution in another trap with frequency $\omega_{2}$ will be periodic with period $\pi/\omega_{2}$ (see Eq. (11)).

Two simulations with the same ratio $a/N_{s}\propto\xi/L$, where $L=\ell N_{s}/2$ is the size of the initial cloud, describe the same physical system with a better accuracy as $a$ and $N_{s}$ are increased. For the results in Fig. 7(a), increasing the number of pixels $N_{s}$ for a fixed $a/N_{s}$ makes the scalar product closer to $1$. If this result could be extended as such to arbitrary large values of $N_{s}$, this would demonstrate that the ground state of a triangular box evolves periodically in a harmonic potential. However, a closer look at the results of this finite-size scaling analysis seems to indicate that $a$ should either be kept constant or increased at a slower rate than $N_{s}$ to have the scalar product approaching 1 in an optimal way. Of course this conjecture deduced from our numerical analysis needs to be further explored with analytical tools, which is out of the scope of the present paper.

The requirement for the Thomas-Fermi regime ($\xi/L\ll 1$) is necessary for obtaining a periodic evolution of the shape with period $T/2$. Indeed, in the ideal gas case ($\tilde{g}=0$), the evolution over $T/2$ corresponds to an inversion of the initial shape with respect to the origin, i.e., a triangle pointing upwards for the case of interest here (Fig. 5(a)). One may then wonder about the existence of a periodicity $T$ for the triangular shape, irrespective of the product $\tilde{g}N$. Indeed this periodicity holds in both limiting cases $\tilde{g}=0$ (ideal gas) and $\tilde{g}N$ large (Thomas–Fermi regime). However numerical simulations show unambiguously that the evolution is not periodic in the intermediate case.

We have also run the same simulations for other simple regular polygons (square, pentagons, hexagon). We did not observe a similar revival of the initial wave function over the time period $(0,5T)$ (see Supplemental Material for details).

Finally we turn to the case of a disk-shaped initial cloud (Fig. 5(c)). The experiment was performed with a cloud prepared such that $|\Delta E|\ll E_{\rm tot}$ in Eq. (7), so that the potential energy is approximately constant over time. In this particular case, the experimental result shown in Fig. 5(d) seems to indicate a periodicity $\approx 2T/7$ for the evolution of the overlap between $n(\boldsymbol{r},0)$ and $n(\boldsymbol{r},t)$. To illustrate this, Fig. 5(c) displays four density distribution at times between [math] and $2T$. Let us assume that this periodicity $2T/7$ is exact when $\Delta E=0$. For a disk-shaped initial distribution with any value of $\Delta E$, the evolution cannot be $2T/7$-periodic. Indeed the potential energy of the cloud is only $T/2$-periodic, which is not a submultiple of $2T/7$. However, all the disk-shaped clouds should have a $2T$ periodicity, which is the least common multiple of $T/2$ and $2T/7$. As we show now, this $2T$ periodicity is well supported by a numerical analysis.

We show in figure 6(c) snapshots of the the calculated density distribution, and in figure 6(d) the time evolution of the overlap $|\langle\psi(0)|\psi(t)\rangle|$, starting from the ground state in a disk-shaped box potential centered on a $512\times 512$ grid. The disk diameter is chosen equal to half the grid size and the simulation is run for $\tilde{g}N=12800$. This simulation shows that the overlap between $\psi(\boldsymbol{r},0)$ and $\psi(\boldsymbol{r},t)$ indeed recovers values close to 1 at times close to multiples of $2T/7$, as observed experimentally.

A closer inspection of figure 6(d) indicates that the time evolution of the overlap is in good approximation periodic with period $2T$, with a symmetry around $t=T$ as well as around $t=2T$. If the evolution is effectively periodic with period $2T$, the symmetry around these points is expected. Indeed the wave function is chosen real for $t=0$, and will thus be real also at $2T$ (up to a global phase). Therefore the evolution must be symmetric around those points thanks to time-reversal symmetry. On the other hand, this symmetry does not show up around the other local maxima $j2T/7$ ($j=1,\ldots,6$), indicating that one does not expect a full overlap with the initial state for those points.

In order to investigate further the revival around $2T$, we have run a finite-size scaling analysis for the same grid sizes as for the triangles, and for $a^{2}=1,2,4,8,16$ (Fig. 7(b)). We find that the numerical results are compatible with a full recovery of the initial wave function at time $2T$, with a scalar product between the wave functions at times [math] and $2T$ attaining a maximum of $0.9986$ for the largest grid size $N_{s}=1024$ and $a^{2}=8$. In this case, the optimal value of $a$ for a given $N_{s}$ (marked with a dot in Fig. 7) increases with $N_{s}$; note that the optimal ratio $a/N_{s}\propto\xi/L$ decreases when $N_{s}$ increases, which guarantees that the cloud remains in the Thomas-Fermi regime.

To conclude this section, we emphasize that the phenomenon described here is notably different from the existence of a breathing mode at frequency $2\omega$ Kagan:1996b ; Pitaevskii:1997a , mentioned in the introduction and explored in Sec. II. Here, we observe a periodic motion of the whole cloud, not just of the second moment $\langle r^{2}\rangle$ of the position. We also note that the observed phenomenon is a genuine non-linear effect, which cannot be captured by a linearization of the motion of the cloud around an equilibrium position. Indeed the state of the gas at an intermediate time may dramatically differ from the state at initial time or after a full period, both in terms of size and shape. A proper analysis of these breathers may require a multimode approach, with the observed phenomenon resulting from a mode synchronization effect via non-linear couplings.

VI Summary and outlook

In this paper, we investigated experimentally some important consequences of the dynamical symmetries of the two-dimensional Gross-Pitaevkii equation, describing the evolution of a weakly interacting Bose gas in a harmonic potential. Firstly, we showed that the SO(2,1) symmetry leads to a periodic evolution of the potential energy and to scaling laws between the evolution of clouds with the same atom number and the same interaction parameter. Secondly, we showed that in the quantum hydrodynamic regime, more symmetries allow one to describe the evolution of the gas by a single universal function irrespective of its size, atom number, trap frequency, and interaction parameter $\tilde{g}$. This universal evolution depends only on the initial shape and velocity field of the cloud. Thirdly, we identified two geometrical box-like potentials, equilateral triangle and disk, which lead to a periodic motion of the wave function when one starts with a gas uniformly filling these shapes and releases it in a harmonic potential of frequency $\omega$. The period of these breathers are $\pi/\omega$ and $4\pi/\omega$ for the triangles and the disks, respectively. This result was confirmed by a numerical simulation for a cloud initially in the Thomas-Fermi regime of the box-like potential, giving an overlap respectively larger than 0.995 and 0.998 between the initial state and the state after one period.

The existence of these breathers raises several interesting questions. First it is not immediate that their existence is a direct consequence of the dynamical symmetries of the system. If this is the case, such breathers could appear also for other systems exhibiting the SO(2,1) symmetry, like a three-dimensional unitary Fermi gas or a cloud of particles with a $1/r^{2}$ interaction potential. Remarkably the later case can be approached using classical (Newton) equations of motion; we have performed a preliminary numerical analysis with up to $10^{5}$ particles, which indicates that an initial triangular (resp. disk) shape with uniform filling also leads to an approximate periodic evolution in a harmonic potential with same period $T/2$ (resp. $2T$) as the solution of the Gross–Pitaevskii equation. We also note that in the 1D case, the spectrum of the Hamiltonian of a gas of particles interacting with a repulsive $1/r^{2}$ potential is composed of evenly spaced energy levels, ensuring a periodic evolution of the system for any initial state Calogero:1971 ; Sutherland:1972 .

The allowed shapes for such breathers is also an intriguing question. In our exploration (both experimental and numerical), we only found this behavior for triangles and disks but one cannot exclude that complex geometrical figures can show a similar phenomenon. Another issue is related with thermal effects. For all studies reported here, we operated with a gas deeply in the degenerate regime, which is well approximated by the zero-temperature Gross–Pitaevskii formalism. A natural extension of our work is therefore to study to which extent the present findings will subsist in the presence of a significant non-superfluid component. For our experimental setup, this will require a significant increase in the vertical trapping frequency, so that the vertical degree of freedom remains frozen for the thermal component of the gas.

Finally we recall that the SO(2,1) symmetry is only an approximation for the description of a two-dimensional Bose gas. It is valid when the gas can be modelled by a classical field analysis, hence for a small interaction parameter $\tilde{g}\ll 1$. For stronger interactions, one has to turn to a quantum treatment of the fluid. This breaks the scale invariance and the SO(2,1) symmetry that exist at the classical field level, providing an example of a "quantum anomaly" Holstein:1993 ; Cabo:1998 ; Olshanii:2010 . For example, the frequency of the breathing mode of a gas in a harmonic potential then differs from its classical value $2\omega$. It remains to be understood if a similar "quantum anomaly" shows up for the breathers described in this work.

Acknowledgements.

This work is supported by DIM NanoK, ERC (Synergy UQUAM), QuantERA ERA-NET (NAQUAS project) and the ANR-18-CE30-0010 grant. We thank Yvan Castin, Cheng Chin, Ignacio Cirac, Lei Feng, Jörg Schmiedmayer and Steven Simon for stimulating discussions.

Appendix A Symmetry groups of the Schrödinger and 2D Gross–Pitaevskii equations

For completeness we summarize in this appendix the main properties of the transformations that leave invariant the Schrödinger equation (i) for a free particle and (ii) for a particle confined in a harmonic potential. The ensemble of these transformations forms a group called the maximal kinematical invariance group, which is parametrized in the 2D case by 8 real numbers. In what follows, we are interested only in the subgroup that is relevant for scale/conformal invariance. For example in the case of a free particle, 5 parameters are related to space translations, changes of Galilean frame and rotations, which do not play a role in our study. We are then left with 3 parameters, corresponding to time translations, dilations and special conformal transformations. These transformations also leave the 2D Gross-Pitaevskii equation invariant. In the following we identify their generators and show that they obey the SO(2,1) commutation algebra. We will follow closely the approach of Niederer:1972 ; Niederer:1973 , which was developed for the Schrödinger equation describing the motion of a single particle, but also applies with little modifications to the case of the non-linear Gross-Pitaevskii equation. In this appendix we set $\hbar=1$ to simplify the notations.

A.1 Free particles

Although we are interested ultimately in the case where the particles evolve in a harmonic potential, we start by a brief summary of the free particle case, for which the algebra is slightly simpler, while involving transformations of a similar type. In Niederer:1972 , it was shown that in addition to space translations, rotations, and Galilean transformations, the three following transformations leave invariant the free-particle Schrödinger equation:

• •

The translations in time

[TABLE]

since the Hamiltonian has no explicit time dependence.

• •

The dilations

[TABLE]

already introduced in Eq. (1) of the main text.

• •

The so-called "expansions":

[TABLE]

which correspond to a special conformal transformation for the time.

The combination of these transformations forms a 3-parameter group with the most general transformation written as:

[TABLE]

with the constraint $\alpha\delta-\beta\gamma=1$. The dilation (28) is obtained by setting $\beta=\gamma=0$, $\delta=\lambda$ and $\alpha\delta=1$.

Let us consider a function $\psi_{1}(\boldsymbol{r},t)$ which is a solution of the Gross-Pitaevskii equation in free space:

[TABLE]

with

[TABLE]

Starting from $\psi_{1}(\boldsymbol{r},t)$, we define the function $\psi_{2}(\boldsymbol{r}^{\prime},t^{\prime})$ as

[TABLE]

with $\boldsymbol{r}^{\prime},t^{\prime}$ set as

[TABLE]

and

[TABLE]

With a tedious, but straightforward calculation, one can check that $\psi_{2}(\boldsymbol{r}^{\prime},t^{\prime})$ is also a solution of the Gross–Pitaevskii equation:

[TABLE]

for any value of the parameters $\alpha,\beta,\gamma,\delta$ with the constraint $\alpha\delta-\beta\gamma=1$. The group of transformations (30, 31) thus allows one to generate an infinite number of solutions of the Gross-Pitaevskii equation. We could pursue this analysis by determining the generators associated with the action of these transformations on the wave functions $\psi(\boldsymbol{r},t)$, but we postpone it to the case of a harmonically confined system which is more relevant for our physical system. The two studies are anyway very similar and the symmetry groups of the two systems have the same Lie algebra Niederer:1972 ; Niederer:1973 .

A.2 Particles in a harmonic trap

In the presence of an isotropic harmonic potential of frequency $\omega$, the general transformations on position and time leaving invariant the Schrödinger equation are also defined by a set of four numbers $(\alpha,\beta,\gamma,\delta)$ with the constraint $\alpha\delta-\beta\gamma=1$ Niederer:1973 . Setting

[TABLE]

the change in position is

[TABLE]

with

[TABLE]

while the transformation on time $t\to t^{\prime}=h(t)$ reads

[TABLE]

Note that the time translations belong to this set of transformations, as expected for a time-independent problem. They are obtained by taking $\alpha=\delta=\cos(\omega t_{0})$ and $\beta=-\gamma=\sin(\omega t_{0})$.

We start with a solution $\psi_{1}$ of the Gross-Pitaevskii equation in the trap:

[TABLE]

with

[TABLE]

Using this group of transformations, we can generate another function $\psi_{2}(\boldsymbol{r}^{\prime},t^{\prime})$ satisfying

[TABLE]

following the definitions (34-35) with now

[TABLE]

The fact that $\psi_{2}$ is a solution of the Gross-Pitaevskii equation was proven for the non-interacting case in Niederer:1973 , and one can check that the contribution of the interaction term proportional to $|\psi|^{2}\psi$ cancels in the 2D case thanks to the scaling $f\propto\lambda$.

In the main text, we use a specific version of the transformation $(\boldsymbol{r},t)\to(\boldsymbol{r}^{\prime},t^{\prime})$ that (i) maps the time $t=0$ onto the time $t^{\prime}=0$, and (ii) is such that $\dot{\lambda}(0)=0$ since we want to relate a real solution $\psi_{1}$ onto another real solution $\psi_{2}$ ($\psi_{1}$ and $\psi_{2}$ are both ground state wave functions in a box-like potential). These two conditions, in association with $\alpha\delta-\beta\gamma=1$, impose $\beta=\gamma=0$ and $\delta=1/\alpha$, hence:

[TABLE]

and

[TABLE]

Finally we note that the simple dilation transformation $\boldsymbol{r}^{\prime}=\boldsymbol{r}/\sqrt{\zeta}$, $t^{\prime}=t/\zeta$ allows one to relate a solution of the Gross-Pitaevskii equation $\psi_{1}(\boldsymbol{r},t)$ in a trap with frequency $\omega_{1}$ to a solution

[TABLE]

in a trap with frequency $\omega_{2}=\zeta\omega_{1}$:

[TABLE]

We can thus combine this dilation with the transformation (46-47) in order to obtain the transformation that links two (initially real) solutions $\psi_{1}(\boldsymbol{r},t)$ and $\psi_{2}(\boldsymbol{r}^{\prime},t^{\prime})$ of the Gross–Pitaevskii equation for a given $\tilde{g}N$, obtained in harmonic traps with frequencies $\omega_{1,2}$ and starting with homothetic initial conditions with characteristic lengths $L_{1,2}$. This transformation reads:

[TABLE]

with

[TABLE]

and $\alpha=L_{2}/L_{1}$, $\zeta=\omega_{2}/\omega_{1}$. This corresponds to the scaling (10) used in the main text.

A.3 Generators and SO(2,1) symmetry

We now look for the infinitesimal generators of the transformation $\psi_{1}\to\psi_{2}$ in the presence of a harmonic potential (Sec. A.2), and show that they fulfill the commutation algebra characteristic of the SO(2,1) group. We focus here on the transformation (39-41) which relates solutions of the Gross–Pitaevskii equation for the same non-linear coefficient $\tilde{g}N$ and the same trap frequency $\omega$.

We first note that the set of four numbers $(\alpha,\beta,\gamma,\delta)$ with the constraint $\alpha\delta-\beta\gamma=1$ actually forms a set of three independent parameters, as for the free-particle case (Sec. A.1). To this set of numbers we can associate a matrix

[TABLE]

of the group SL(2,R). In order to simplify our discussion, we consider the following three subgroups of SL(2,R), each parametrized by a single parameter $s_{j}$, $j=1,2,3$:

[TABLE]

and

[TABLE]

We will obtain three independent generators by considering a small displacement from the unit matrix for each subgroup ($|s_{j}|\ll 1$). In all three cases, we will write the passage from $\psi_{1}$ to $\psi_{2}$ as:

[TABLE]

where we introduce the time-dependent generator $\hat{L}_{j}(t)$. The goal is to determine explicitly these operators and their commutation relation, in order to check that they satisfy the SO(2,1) algebra.

Generator associated to $s_{1}$.

We have in this case

[TABLE]

where the $\hat{\sigma}_{j}$, $j=x,y,z$, are the Pauli matrices. We first get $\lambda(t)=1-\frac{s_{1}}{2}\cos(2\omega t)$ so that

[TABLE]

and the infinitesimal changes in $\boldsymbol{r},t$ are

[TABLE]

This allows one to determine the passage from $\psi_{1}$ to $\psi_{2}$ as in (55) with

[TABLE]

Generator associated to $s_{2}$.

We find

[TABLE]

In this case $\lambda(t)=1+\frac{s_{2}}{2}\sin(2\omega t)$, and

[TABLE]

It also provides the transformation of space and time coordinates:

[TABLE]

This corresponds to a transformation similar to the one considered above in Eq. (58), with the time translation $t\to t+\pi/(4\omega)$. The associated operator for the passage from $\psi_{1}$ to $\psi_{2}$ is thus

[TABLE]

Generator associated to $s_{3}$.

Finally we have for the third case

[TABLE]

We simply have $\lambda(t)=1$, $f(\boldsymbol{r},t)=1$, and this case corresponds to the time translations mentioned above, for which we have

[TABLE]

The operator $\hat{L}_{3}(t)$ is thus

[TABLE]

From the expressions of the three generators $\hat{L}_{j}$ we easily find the commutations relations valid at any time

[TABLE]

which are characteristic of the Lorentz group SO(2,1). As explained in Pitaevskii:1997a , this set of commutation relations allows one to construct in particular families of solutions with an undamped breathing motion.

S1: Measurement of the vertical frequency

The value of the interaction parameter $\tilde{g}=\sqrt{8\pi}a_{s}/\ell_{z}$ is determined by the measurement of the vertical frequency $\omega_{z}$, since $\ell_{z}=\sqrt{\hbar/(m\omega_{z})}$. We measure this frequency for different values of the intensity of the laser producing the vertical confinement with the following procedure. We prepare a cloud with a high intensity of the laser, which means a high vertical confinement ($\omega_{z}/(2\pi)\approx 7\,$\mathrm{kHz}), and in $1\,$\mathrm{ms} we reduce this intensity to a final value $I$, giving a vertical kick to the cloud due to gravity and exciting a dipolar oscillation along the vertical axis. We let it evolve during a variable amount of time, then remove all confinements to let the cloud expand freely during $12\,$\mathrm{ms}. We measure its vertical position that varies sinusoidally with time. A fit of this vertical position provide us an estimate of $\omega_{z}$. On Fig. [S8](#A1.F8) are shown the results of these measurements as a function of $\sqrt{I}$. The points lie on a line as expected, and a linear fit gives $\omega_{z}/(2\pi)=4.95\,$\mathrm{kHz}$\times\sqrt{I}-0.03\,$\mathrm{kHz}, where $\sqrt{I}$ is expressed in the arbitrary unit used on Fig. S8. We are then able to extrapolate the value of the vertical frequency, and therefore of $\tilde{g}$ for any intensity $I$ in the range explored on the figure.

S2: Details on the determination of the scaling laws

We develop here with an example the procedure to reconstruct the scaling laws that we have used in sections III and IV C. We consider the density distribution $n_{1}(t_{1})$ in the series presented on Fig. 2(a) at time $t_{1}=5.9\,$\mathrm{ms}, and we search for the time $t^{\prime}_{\mathrm{opt}}(t_{1})$ and the rescaling parameter $\Lambda(t_{1})$ that relate this density distribution to the series presented on Fig. 2(b). We show on Fig. [S9](#A1.F9)(a) the overlap between the image $n_{1}(t_{1})$ and all the images $n_{2}(t^{\prime})$ of the second series. The optimal time is found to be $t_{\mathrm{opt}}(t_{1})=4.3(7)\,$\mathrm{ms}, indicated by the maximum of the parabolic fit. The time $t^{\prime}_{1}$ that is the closest to $t^{\prime}_{\mathrm{opt}}(t_{1})$ for which we have an image is $t^{\prime}=4\,$\mathrm{ms}$$. On Fig. S9(b) we show the variation with $\lambda$ of the scalar product $p\left(n_{1}(t_{1}),n_{2}(t^{\prime}_{1}),\lambda\right)$. The optimal parameter is $\Lambda(t_{1})=0.81(5)$, indicated by the parabolic fit.

S3: Time evolution of density distributions in a harmonic trap (experiment)

The time evolution of the cloud whose potential energy is shown on Fig. 1 is displayed on Fig. S10. The imaging is performed by transfering a fraction of the atoms in the state $F=2,m_{F}=0$ and an absorption imaging is performed with a laser beam resonant with the $F=2\to F^{\prime}=3$ transition. The fraction of the transfered atoms is chosen so that the optical density of the imaged atoms doesn’t exceed 1. The value of the atomic density is then obtained with $n=\mathrm{OD}/(\varepsilon\sigma)$, where $\mathrm{OD}$ is the measured optical density, $\varepsilon$ is the fraction of the transfered atoms and $\sigma$ is the scattering cross-section of the atoms.

On Fig. S10, the image at $\omega t=\pi/2$ corresponds to $t=12\,$\mathrm{ms}$$, and we recover the observation that an initial square shape evolves as a diamond-shaped cloud, as discussed in the main text (Sec. IV B).

S4: Time evolution of density distributions in a harmonic trap (numerical)

We show in figures S11-S15 the time evolution of the density distribution, calculated using the numerical solution of the time-dependent Gross–Pitaevskii equation. The initial distribution corresponds to the ground state in a box-like potential. The box has the shape of a regular polygon (triangle, square, pentagon, hexagon) or a disk. It is centered on the calculation grid, and its summits are located on a circle of diameter $L/2$, where $L$ is the total size of the grid. As in the experiment, the box confinement is removed at time $t=0$ and the wave function then evolves in an isotropic harmonic confinement with frequency $\omega$. The frequency $\omega$ of the harmonic potential is chosen such that $E_{\rm p}(0)\approx E_{\rm tot}/2$.

Here the calculations are performed on a grid $512\times 512$, with an initial state in the Thomas–Fermi regime. The healing length is $\xi\approx 1.4\,\ell$ for the triangular case and $\xi\approx 2\,\ell$ for the other cases, where $\ell$ is the grid spacing. Each panel contains 40 images, which display (from left to right and from top to bottom) the density distribution at time $t=jT/8$, with $j=0,\ldots,39$.

We show in figure S16 the evolution of the overlap $|\langle\psi(0)|\psi(t)\rangle|$ for the set of data of Figs. S11–S15 (grid $512\times 512$), together with the evolution of the same quantity calculated on a smaller grid ($256\times 256$, dotted line). The periodicity $T/2$ (resp. $2T$) appears clearly for an initially triangular (resp. disk-shaped) cloud. For these two cases, the local maxima of the overlap get closer to unity as the grid size increases, and the evolution of the overlap is approximately symmetric around the times $t=jT/4$ for the triangles ($t=jT$ for the disks), as expected from time reversal symmetry arguments for a periodic evolution (see main text). On the opposite, no such hint for a periodic behavior shows up on this time interval for a cloud whose initial shape is a uniform square, pentagon or hexagon.