Dynamical symmetry and breathers in a two-dimensional Bose gas

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.


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 [1].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 scaleinvariant systems were initially introduced in particle physics to explain scaling laws in high energy collisions [2].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: In this dilation, the velocity of a particle is changed as v → λv.Therefore the kinetic energy of the fluid scales as E k → λ 2 E k , which ensures that the corresponding part of the action (∝ E k dt) remains invariant in the transformation (1).If the interaction energy has the same scaling, E i → λ 2 E 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 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) [3].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 [4,5] and particles interacting with a 1/r 2 potential [6].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.[7]).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 [8,9], and also to general relations between the moments of the total energy and those of the trapping energy in a harmonic potential [10].
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 ψ(r, t).This field is normalized to unity ( |ψ| 2 d 2 r = 1), so that the density of the gas arXiv:1903.04528v2[cond-mat.quant-gas]26 Apr 2019 reads n(r, t) = N |ψ(r, t)| 2 where N is the number of particles.In the scaling of positions, the 2D matter-wave field changes as ψ(r) → λ ψ(λr), which guarantees that the norm is preserved and that the dynamical part of the action ∝ i dt d 2 r ψ * ∂ t ψ is invariant.The interaction energy of the gas then reads for contact interaction where m is the mass of a particle and g the dimensionless parameter characterizing the strength of the interaction.
One can immediately check that E i obeys the λ 2 scaling required for scale invariance, which can be viewed as a consequence of the dimensionless character of g.The classical field description used here is valid if one restricts to the case of a small coupling strength g 1 [11].This restriction is necessary because of the singularity of the contact interaction 2 m gδ(r) in 2D when it is treated at the level of quantum field theory.Note that the condition g 1 does not constrain the relative values of the interaction and kinetic energies.Actually in the following we will often consider situations where So far, the scale/conformal invariance of the weakly interacting 2D Bose gas has been mainly exploited to measure its equation of state [12,13].Also, one of its dynamical consequences in an isotropic 2D harmonic potential of frequency ω has been explored: the frequency of the breathing mode was predicted to be exactly equal to 2ω for any g [14][15][16], as tested in Refs.[17,18].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 [15].Recently, deviations from this prediction for g 1, an example of a quantum anomaly [19], have been observed [20,21].
The purpose of our work is to go beyond static properties of the weakly interacting 2D Bose gas and its singlemode 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 ω, starting from a uniformly-filled simple area (disk, triangle or square).Here we use g ≤ 0.16 so that the classical field description is legitimate.We first check (Sec.II) the prediction from [15] that E k + E i should have a periodic evolution in the trap with the frequency 2ω.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π/ω).These particular shapes can therefore be viewed as two-dimensional breathers for the Gross-Pitaevskii (non-linear Schrödinger) equation in the hydrodynamic limit [22].They also constitute a novel example of universal dynamics in a quantum system prepared far from equilibrium [23][24][25].

II. EVOLUTION OF POTENTIAL ENERGY
Our experiment starts with a 3D Bose-Einstein condensate of 87 Rb 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 ω z /(2π) up to 4.9 kHz.The interaction parameter is g = √ 8π a s / z , where a s is the 3D s-wave scattering length and z = ( /mω z ) 1/2 .The interaction energy per particle and the residual temperature are both smaller than ω z so that the vertical degree of freedom is effectively frozen [26].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 [27], 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 ω/2π around 19.5 Hz.We estimate the anisotropy of the potential to be 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(r) = N |ψ(r)| 2 by absorption imaging.
The measurement of n(r) gives access to both the interaction energy (2) and the potential energy in the har-monic trap Since the gas is an isolated system, we expect the total energy E tot = E k + E i + E p to be conserved during the evolution, where the kinetic energy E k reads: The SO(2,1) symmetry for a 2D harmonically trapped gas brings a remarkable result: E k + E i and E p should oscillate sinusoidally at frequency 2ω [15].More precisely, using the 2D Gross-Pitaevskii equation one obtains the relations where we have defined W = ωm r • v n d 2 r, and the velocity field v(r) = m Im [ψ * (r)∇ψ(r)] /|ψ(r)| 2 .Initially the gas is prepared in a steady state in the box potential, so that v = 0 hence W (0) is also null.Therefore the potential energy evolves as where ] can be positive or negative.A similar result holds for the sum E k + E i (with ∆E replaced by −∆E), but not for the individual energies E k or E 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 p evolves periodically and is well fitted by a cosine function with a period that matches the 2ω 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 (∼ 0.1 s −1 ).
This simple dynamics can be viewed as a generalisation of the existence of the undamped breathing mode at frequency 2ω that we mentioned in the introduction [14,15].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 ψ 1,2 of this equation corresponding to homothetic initial conditions: one can relate ψ 1 (r, t) and ψ 2 (r , t ), provided they evolve with the same parameter gN and the same trap frequency ω 1 = ω 2 .By using a simple scaling on space and time, this link can be further extended to the case 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 ψ 1 (r, t) for the harmonic potential of frequency ω 1 : Using scale/conformal invariance, we can construct a solution ψ 2 (r , t ) of the Gross-Pitaevskii equation with the frequency ω 2 = ζω 1 using: where space is rescaled by r = r/λ(t) with where the dimensionless parameter α is the homothetic ratio between the initial states.The relation between the times t and t in frames 1 and 2 is and the multiplicative function f is where λ ≡ dλ dt .The two solutions ψ 1,2 (t) correspond to the evolution of two clouds with the same parameter g1 N 1 = g2 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 = αL 1 .Both initial wave functions ψ 1,2 (0) can be chosen real, and the scale invariance of the (timeindependent) 2D Gross-Pitaevskii equation ensures that they are homothetic: αψ 2 (αr, 0) = ψ 1 (r, 0).For example in the limit E i E k , ψ(0) corresponds to a uniform density in the bulk, and goes to zero at the edges on a scale given by the healing length ξ ≡ [N 2 /(2mE i )] 1/2 .For two box potentials of homothetic shapes filled with the same number of particles, the ratio ξ 2 /ξ 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 ω 1 → 0, i.e., α = 1 and ζ → +∞.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) µm with 3.9(3) × 10 4 atoms and let it evolve either in a harmonic potential of frequency ω 2 /(2π) = 19.7(2)Hz, or without any potential (ω 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 reads and the relation (10) becomes: The relation (13) indicates that the scaling transformation maps the first quarter of the oscillation period in the harmonic trap ω 2 t ≤ π/2 onto the ballistic expansion from t = 0 to t = ∞.In the absence of interactions, this result has a simple physical interpretation: after the ballistic expansion between t = 0 to t = ∞, 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 [28].
In order to reconstruct the scaling laws ( 13) and ( 14) from the measured evolutions, we compare each image n 1 (r, t) for the free evolution with the set of images n 2 (r , t ) obtained for the in-trap evolution.More precisely, we start by defining the overlap O[n 1 , n 2 ] between two images in the following way: • We introduce the scalar product (n 1 |n 2 ) between two images and the norm of an image ||n 1 || = (n 1 |n 1 ).
• In order to relate two images that differ by a spatial scaling factor λ, we introduce the quantity where n (λ) 1 (r) = λ 2 n 1 (λr) is the image rescaled from n 1 (r) by the factor λ, with the same atom number: Note that the definition of the norm given above entails ||n is identical to n 2 up to a multiplicative factor.
• Finally, for a couple of images (n 1 , n 2 ) we vary λ and define their overlap as In practice, for each image n 1 (t) acquired at a given time t, we determine the time t opt where the overlap between n 1 (t) and n 2 (t ) is optimal.We denote Λ(t) the value of the scaling parameter λ for which the value O[n 1 (t), n 2 (t 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 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 Λ(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.Optimal rescaling factor between the corresponding images n1(t) and n2(t opt ).On the two graphs (c) and (d), the solid lines are the theoretical predictions given by ( 13) and ( 14).The inset of figure (d) presents the overlap between the corresponding images of the two series.On (c) and (d), the error bars indicate the confidence intervals within two standard deviations of the fits used to reconstruct the scaling laws.
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 g1 N 1 = g2 N 2 , since only the product gN 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 gN ).We show now that it is also possible to link evolutions with different N 's and g's, provided we restrict to the so-called hydrodynamic (or Thomas-Fermi) regime, where the healing length ξ is very small compared to the size of the gas.

A. General formulation
The Gross-Pitaevskii equation ( 8) can be equivalently written in terms of the density and the velocity fields as where √ 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 ξ, one can neglect the contribution of the quantum pressure in (19): 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 v [29].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 ξ 2 /ξ 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 ∼ ξ 1,2 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 ω 1,2 .Our goal is to relate the two evolutions with parameters The general transformation involves three dimensionless constant parameters µ, α, ζ: and reads: with λ = dλ dt .The spatial variables are rescaled as r = r/λ(t) with the function λ now given by and the relation between the times t and t in frames 1 and 2 is: With a calculation similar to that detailed in Appendix A, one can readily show that if (n 1 , v 1 ) is a solution of the hydrodynamic equations (18,20) for the frequency ω 1 , then (n 2 , v 2 ) is a solution for the frequency ω 2 .If µ = 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.
B. Connecting evolutions with a fixed trap frequency, a fixed size and different gN We present here the experimental investigation of the scaling described above, focusing on the case L 1 = L 2 and ω 1 = ω 2 , i.e., α = ζ = 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 such that ω 1 t = ω 2 t = π/2, which satisfies the constraint (25).In this case λ(t) = 1/µ, so that the general scaling (22) reads We start with a cloud in a uniform box potential with the shape of an equilateral triangle of side length L = 38.2(3)µm.At t = 0 we transfer the atoms in the harmonic trap of frequency ω/2π = 19.6Hz and remove the box potential.At t = π/(2ω) we image the cloud.We perform this experiment for different values of gand slightly different atom numbers -corresponding to the product gN between 200 and 4000.This leads to a ratio ξ/L always smaller than 0.03, ensuring that we stay in the quantum hydrodynamic regime.The variation of g is achieved by changing the intensity I of the laser beams creating the vertical confinement, with g ∝ I 1/4 .The values of g are obtained from the measurement of the vertical frequency ω 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 gN ≈ 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 gN 's, and extract an optimal scaling parameter Λ.The results of this analysis are displayed on Fig. 3.The inset shows that the overlap is close to 1 for all values of gN , 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 Λ −2 with gN .The scaling law (24) predicts that Λ −2 = µ 2 ∝ gN , 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 gN 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 = π/(2ω), i.e., t = ∞  25) and (24) where the values of the parameters α and µ are measured independently.The uncertainty of these values are represented as a shaded area.On (d), this area is too narrow to be discernable.On (c)-(e), the error bars indicate the confidence intervals within two standard deviations of the fit that we use to reconstruct the scaling laws.They are too small to be seen on (d).
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 [30].Note that we also observe such a diamond-like shape at t = π/(2ω) starting from a square box, albeit with a non-uniform density (see Supplemental Material).

C. Connecting evolutions with a fixed trap frequency, different sizes and different gN
Finally we compare the evolution of two clouds with homothetic shapes and α, µ = 1, ζ = 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)µm, contains N 1 = 3.7(3) × 10 4 atoms, and its initial density distribution is shown on Fig. 4(a).The second one has a side length L 2 = 36.8(5)µm and contains N 2 = 5.4(3) × 10 4 atoms (Fig. 4(b)).The ratio ξ/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 g, and take pictures after different evolution times.We expect that the two evolutions n 1 (r, t) of the first cloud and n 2 (r , t ) of the second cloud are linked via Eqs.(22,24,25), with parameters α = L 2 /L 1 = 1.36(4) and µ = 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 of the second frame and the corresponding time t of the first frame is shown on Fig. 4(d), and the best scaling factor Λ(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 α, µ and ζ, 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 p is periodic with period T /2 ≡ π/ω for an arbitrary initial state ψ(r, 0) [see Eq. ( 7)].Of course, the existence of this periodicity does not put a strong constraint on the evolution of ψ(r, t) itself.Due to the non-linear character of the Gross-Pitaevskii equation, the evolution of ψ 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 uniformlyfilled squares, pentagons or hexagons, we indeed observe that even though E p (jT /2) = E p (0) for integer values of j, the shapes n(r) = N |ψ(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 ψ(r, t) that undergoes a periodic evolution in an isotropic harmonic trap of frequency ω (for a generalisation to different settings, see e.g.[22,31]).According to this definition, the simplest example of a breather is a steady-state ψ ss (r) of the Gross-Pitaevskii equation, e.g. the ground state.Other breathers are obtained by superposing ψ ss with one eigenmode of the Bogoliubov-de Gennes equations, resulting from the linearization of the Gross-Pitaeveskii equation around ψ ss .In principle (with the exception of the breathing mode [15]), 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 [32].Note that for a negative interaction coefficient g in 1D, a bright soliton forms a stable steady-state of the Gross-Pitaevskii equation (even for ω → 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.[33] 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 sim-ulating the evolution of a wave function ψ(r, t) with the Gross-Pitaevskii equation [34].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 | ψ(0)|ψ(t) | between the wave functions at times 0 and t.The calculation was performed on a square grid of size N s × 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 gN = 25600.Note that by contrast to the "scalar product between images" introduced above, the quantity | ψ(0)|ψ(t) | 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 max ≈ T /2, for increasing grid sizes N s = 64, . . ., 1024.The product gN is adjusted such that the healing length ξ = [N 2 /(2mE i )] 1/2 = a , where is the grid spacing and a 2 = 0.5, 1, 2, 4, 8.The condition a N s ensures that ξ is much smaller than the size of the triangle (Thomas-Fermi regime), while having a 1 provides an accurate sampling of the edges of the cloud.The overlap between |ψ(0) and |ψ(t max ) increases with the grid size, and reaches 0.995 for the largest grid.
In the simulation, the trapping frequency ω is adjusted such that |∆E| E 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 ψ(r, t = 0) in traps with different frequencies.In particular, if the evolution starting from ψ(r, 0) in a trap of frequency ω 1 is periodic with period π/ω 1 , the evolution in another trap with frequency ω 2 will be periodic with period π/ω 2 (see Eq. ( 11)).Two simulations with the same ratio a/N s ∝ ξ/L, where L = 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), in-creasing 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 (ξ/L 1) is necessary for obtaining a periodic evolution of the shape with period T /2.Indeed, in the ideal gas case (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 gN .Indeed this periodicity holds in both limiting cases g = 0 (ideal gas) and gN 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 |∆E| E 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 ≈ 2T /7 for the evolution of the overlap between n(r, 0) and n(r, t).To illustrate this, Fig. 5(c) displays four density distribution at times between 0 and 2T .Let us assume that this periodicity 2T /7 is exact when ∆E = 0.For a diskshaped initial distribution with any value of ∆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 | ψ(0)|ψ(t) |, starting from the ground state in a disk-shaped box potential centered on a 512 × 512 grid.The disk diameter is chosen equal to half the grid size and the simulation is run for gN = 12800.This simulation shows that the overlap between ψ(r, 0) and ψ(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, . . ., 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 0 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 ∝ ξ/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ω [14,15], 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 r 2 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.In-deed 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 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 ω.The period of these breathers are π/ω and 4π/ω 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 boxlike 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 [35,36].
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 g 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" [19,37,38].For example, the frequency of the breathing mode of a gas in a harmonic potential then differs from its classical value 2ω.It remains to be understood if a similar "quantum anomaly" shows up for the breathers described in this work.and σ is the scattering cross-section of the atoms.
On Fig. S10, the image at ωt = π/2 corresponds to t = 12 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 ω.The frequency ω of the harmonic potential is chosen such that E p (0) ≈ E tot /2.
Here the calculations are performed on a grid 512×512, with an initial state in the Thomas-Fermi regime.The healing length is ξ ≈ 1.4 for the triangular case and ξ ≈ 2 for the other cases, where 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, . . ., 39.
We show in figure S16 the evolution of the overlap | ψ(0)|ψ(t) | for the set of data of Figs.S11-S15 (grid 512×512), together with the evolution of the same quantity calculated on a smaller grid (256 × 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.

FIG. 1 .
FIG. 1.Time evolution of the potential energy per particle of a 2D gas of 87 Rb atoms in an isotropic harmonic potential of frequency ω for a square of side length 27.6(5) µm with 4.1(2) × 10 4 atoms.(a): Evolution of the potential energy per particle.Each point is an average of 7 to 10 realizations, and the error bars indicate the standard deviation of these different realizations.The frequency of the trap is measured with the oscillation of the center of mass: ω/2π = 19.3(1)Hz.The oscillations of Ep are fitted with a cosine function and an additional linear slope (continuous line).This slope is −0.25(4)Hz/ms and accounts for the loss of particles from the trap.The fitted frequency is 38.5(1)Hz, which is compatible with ω/π, as predicted by the SO(2,1) symmetry of the gas.(b): Density distribution of an initially uniform gas after the evolution in a harmonic potential at times ωt = 0, π, 2π, 3π, 4π, corresponding to the first periods of the potential energy indicated by the labels from 1 to 5. The horizontal black lines represent 10 µm.

FIG. 2 .
FIG. 2. Evolution of a gas with triangular shape (side length 40.2(3) µm, 3.9(3) × 10 4 atoms) for two different values of the harmonic trapping frequency.(a) and (b): Averaged images of the density distribution after a variable evolution time in the harmonic potential of frequency ω1 = 0 and ω2/2π = 19.7(2)Hz respectively.The images result of an average over 5 to 10 realizations, and the horizontal black lines represent 10 µm.Pairs of images with approximately corresponding evolution times have been chosen.(c): Optimal time t opt (t) for which the overlap between images of the first and the second evolutions is maximum.(d):Optimal rescaling factor between the corresponding images n1(t) and n2(t opt ).On the two graphs (c) and (d), the solid lines are the theoretical predictions given by (13) and (14).The inset of figure (d) presents the overlap between the corresponding images of the two series.On (c) and (d), the error bars indicate the confidence intervals within two standard deviations of the fits used to reconstruct the scaling laws.

FIG. 3 .
FIG. 3. Scaling factor at ωt = π/2 for different values of gN .(a): Initial density distribution of the cloud.(b)-(d): Density distributions of the cloud after an evolution during t = π/(2ω) in the harmonic trap for different values of gN .For (a)-(d), the horizontal black lines represent 10 µm.Main graph: Best scaling factor Λ −2 as a function of gN .The red square corresponds to the reference image and its ordinate is fixed to 1.The solid line represents the prediction (26).The shaded area represents its uncertainty, due to the one in the atom number of the reference point.The vertical error bar represent the precision at two standard deviations of the fit that determines Λ −2 .(e): Value of the overlap between the density distributions and the reference point.The error bars due to the fit are smaller than the black points.

FIG. 4 .
FIG. 4. Mapping between two clouds with the same shape, different sizes and different atom numbers.(a), (b): Initial density distribution of the two clouds.The horizontal black lines represent 10 µm.(c) Best overlap between each image of the first series of images and the images of the second one.(d): Optimal time t opt of the second evolution as a function of the time t of the first evolution.(e): Optimal scaling factor Λ(t) between the first and second evolutions.On (d) and (e), the solid lines are the predictions (25) and(24) where the values of the parameters α and µ are measured independently.The uncertainty of these values are represented as a shaded area.On (d), this area is too narrow to be discernable.On (c)-(e), the error bars indicate the confidence intervals within two standard deviations of the fit that we use to reconstruct the scaling laws.They are too small to be seen on (d).

FIG. 5 .
FIG. 5. (a): Density distributions of an initially triangular-shaped cloud at t/T = 0, t/T = 0.08, t/T ≈ 1/4 and t/T ≈ 1/2.The first and last distributions are close one to another.(b): Scalar product between the initial density distribution of a triangular-shaped cloud (red square) and the density distributions during its evolution in the harmonic trap.The first point is fixed at 1.The dashed lines indicate where t/T is a multiple of 1/2.The shape seems to be periodic of period T /2.(c): Density distributions of an initially disk-shaped cloud at t/T = 0, t/T ≈ 2/7, t/T ≈ 1 and t/T ≈ 2. The first two and the last distributions are close one to the others.(d): Scalar product between the initial density distribution of a disk-shaped cloud (red square) and the density distributions during its evolution in the harmonic trap.The first point is fixed at 1.The dashed lines indicate where t/T is a multiple of 2/7.The shape seems to be periodic of period 2/7.On (a) and (c), the horizontal black lines represent 10 µm.On (b) and (d), the black arrows indicate the point corresponding to density distributions shown on (a) and (c) respectively.The error bars represent the statistical error of the measurement.

FIG. 6 .FIG. 7 .
FIG. 6. (a): Calculated density distributions at times t/T = 0, 1/8, 1/4, 1/2 and (b): calculated time evolution of | ψ(0)|ψ(t) |, starting from the ground state in a triangular box.The numerical integration of the Gross-Pitaevskii equation is performed on a 512 × 512 grid.The triangle is centered on the grid, with a side length equal to half the grid size.We chose gN = 25600 corresponding to an initial healing length ξ ≈ , where is the grid step.(c): Calculated density distributions at times t/T = 0, 2/7, 1, 2 and (d): calculated time evolution of | ψ(0)|ψ(t) |, starting from the ground state in a disk-shaped box.The numerical integration of the Gross-Pitaevskii equation is performed on a 512 × 512 grid.The disk is centered on the grid, with a diameter equal to half the size of the grid.We chose gN = 12800 leading to an initial healing length ξ ≈ 2 , where is the grid step.On (b) and (d), the black arrows indicate the times corresponding to the snapshots presented on (a) and (c).
FIG. S15.Time evolution of the density |ψ(r, t)| 2 calculated from the solution of the Gross-Pitaevskii equation for an initial disk box potential (see text for details).Each image correspond to time t = jT /8, with j = 0, . . ., 39.The evolution is approximately periodic with period 2T .