Sebastiano Saccani

Sissa Moroni

E Vitali

M Boninsegni

Bose soft disks: a minimal model for supersolidity

Keywords: Supersolid, Cluster Solid, Bose Supersolidity, Cluster Solid, Bose Gas

The low temperature phase diagram of an assembly of Bose particles, interacting via a simple soft disk potential, is studied numerically by means of Monte Carlo simulations. At sufficiently high densities, the system features a cluster crystal phase, as observed in its classical counterpart. In the low temperature limit, and in a relatively wide region of parameter space, the quantum cluster crystal undergoes a transition to a supersolid phase, displaying both a finite superfluid fraction and crystalline order. We argue that soft-core Bosons ought be regarded as a "minimal model" of continuous-space supersolidity. With the aim of elucidating the interplay between superfluid response and solid order, we compute the excitation spectrum of the system. A comparison of that of the insulating cluster crystal with that of the supersolid, shows that the latter features an additional acoustic mode.

Introduction

The elusive supersolid phase of matter has intrigued condensed matter and manybody physicists for over half a century. The observation in 2004 by Kim and Chan [1] of non-classical rotational inertia in a crystal of 4 He rekindled the interest in such a phase of experimenters and theorists alike. Solid 4 He is universally believed to be the most likely candidate to display a supersolid phase, among naturally occurring condensed matter systems. At the time of this writing, however, the controversy over whether the anomaly observed by Kim and Chan indeed signals the transition of the system to a supersolid phase, has not yet been resolved [2]. Besides the feebleness of the observed superfluid response, as well its puzzling dependence on details of the measuring apparatus, it seems as if a significant source of ambiguity in the interpretation of existing data stems from the currently generally limited understanding of the supersolid phase.

In order to achieve greater insight, it would be beneficial to identify a theoretical many-body system whose supersolid properties could be established unambiguously, by means of rigorous calculations. That would allow one to characterize the phase not only through its most obvious qualities (namely superflow and crystalline order), but also some, like dynamical responses, that could be probed by different experiments than those carried out until now. Presumably, that would not only help resolve the controversy over solid helium, but also facilitate the detection of supersolid phases in other systems. Recently, Cinti et al. [3] investigated by Monte Carlo simulations a twodimensional model of bosons interacting via a purely repulsive pair-wise potential, featuring a dipolar 1/r 3 long-range tail but saturating to a constant finite value at short distance. This particular form of potential might be engineered in cold atomic systems, through a mechanism known as Rydberg blockade [4]. Such a system displays for sufficiently high density a crystal phase, i.e., one in which each lattice site is occupied by a number of particles K larger than one. In turn, such a crystal can turn superfluid at low enough temperature, phase coherence being established by quantum-mechanical tunnelling of particles across adjacent lattice sites.

Bounded pair potentials have been employed in classical systems as well, to describe the effective interactions between macromolecules. In that context too, it was found that a certain class of bounded repulsive pair potentials give rise, for sufficiently high density and low temperature, to multiple-occupancy crystal phases [5]. In particular, a criterion is known for the quantitative prediction of clustering, based on the form of the inter-particle potential [6]. We refer henceforth to the bunch of particles at each lattice site, as "cluster" or "droplet", and to the resulting phase as "cluster crystal" ("cluster supersolid" when appropriate).

The possible relevance of the long-range tail of the interaction, in stabilizing the supersolid phase, was recently studied [7] by considering a strictly short-range, soft disk type inter-particle potential,

v(r) = D if r ≤ a 0 if r > a (1)
where D is the potential barrier that two particles have to overcome in order to overlap (independent of the extent of the overlap), and a is the disk radius. The similarity between the phase diagram obtained with this simple interaction [7] and that found in the presence of a dipolar (or faster-decaying) tail [3], points to the fundamental irrelevance of the long-range tail of the potential. The soft-disk system can thus be considered a "minimal" underlying model for the presence of supersolidity in the continuum. As such, it is a worthwhile subject for further investigation for the reasons mentioned above, namely as a toy model embodying basic properties of the supersolid phase. With this aim, in this work we carry out a numerical study of several key properties of the soft disk model and its supersolid phase. Specifically, we compute the momentum distribution, the potential felt by a test particle (related to tunnelling of particles between different clusters), and the superfluid fraction. We also explore the evolution of the statistics of cycles of permutation of identical particles across the supersolid transition, and compare the excitation spectra of the supersolid and the insulating crystals.

The paper is organized as follows: in the next section we introduce the notation and briefly describe the computational methodology, devoting the rest of the manuscript to the presentation of the results.

Model and methodology

We consider a system of Bose particles of spin zero in two dimensions. The Hamiltonian of the system is

H = - 1 2 N i=1 2 i + D i>j Θ(1 -r ij ), (2)
where r ij is the distance between particles i and j, all lengths are expressed in units of the soft-core diameter a, while all energies are expressed in units of • = 2 /ma 2 . The parameter D ≡ V / • can also be expressed as (a/ξ) 2 , where ξ is the quantummechanical penetration length of a potential barrier of height V . In the limit ξ → 0 the model (2) reduces to the hard-sphere gas. Path Integral Monte Carlo simulations of the system described by (2), enclosed in a cell with periodic boundary conditions of sides (L x , L y ) , have been performed by means of the Continuous-space Worm Algorithm in the grand canonical ensemble (i.e. at fixed temperature T , area A and chemical potential µ). Because this methodology is by now well-established, and is thoroughly described elsewhere, we omit here technical details, and refer the interested reader to the original references [8,9]. We denote by N the average number of particle, by N c the number of clusters contained in the simulation cell and by l the lattice parameter. We express the density ρ in terms of the dimensionless parameter r s = 1/ ρa 2 . The presence of a global superfluid response, extending to the whole system, is assessed numerically through the winding number estimator [10] of the superfluid fraction ρ s .

Results

In Figure 1 we summarize the main features of the phase diagram obtained in Ref. [7]. In the T = 0 limit, it is possible to identify two qualitatively different behaviors of the system in the strongly interacting (high D) and weakly interacting (low D) regimes. In the former case (D = 60 in Figure 1), for low values of the density the system is a superfluid gas, essentially behaving like an hard-disk fluid. This phase undergoes solidification into a triangular crystal on increasing µ. For sufficiently large values of D, the number of particles per cluster (unit cell) is K = 1. However, on further increasing the chemical potential, particles bunch into clusters which organize in a solid preserving the triangular structure. This phenomenon is classical in nature, being related to the energy cost associated to particles being at a distance less than the soft core diameter.

For lower values of D, the system can be kept superfluid up to a much higher density (r s < 1, meaning that particles easily overlap). Above a certain µ the system crystallizes directly into a cluster crystal with K > 1, the exact value depending on D. Moreover, even though density modulations appear in the system, the superfluid fraction in the system is not fully suppressed. This is the aforementioned cluster supersolid phase. This phenomenon is underlain by long cycles of exchanges of particles permutations between neighbouring clusters, eventually leading to the superfluidity of the system as a whole.

Upon further increase of the density, ρ s is progressively suppressed by the rising potential barrier between neighbouring lattice sites, which is associated with tunnelling probability. This eventually destroys phase coherence in the system, which goes into a normal (i.e. non superfluid) cluster phase. The same effect of tunnelling suppression can be obtained also by increasing D at constant r s , instead of µ, as can be seen in Figure 2, where the averaged potential felt by a "test" particle is plotted along the cristallographic direction [1, 0] passing trough the lattice sites. Note that the potential has minima in correspondence of the clusters. Higher potential barriers lead to increasingly compact clusters, i.e. particles pile up on a smaller spatial region (compare, e.g., the spread of the world lines in the two right panels of Figure 1). The pattern illustrated in Figure 2 suggests a correpondence with a soft-core Hubbard model, although in the present case the triangular lattice is self-assembled and not superimposed.

Finally, contraction of the world lines with ensuing demotion of tunnelling occurs when the temperature is increased. This is reflected in the distribution of the lengths of permutation cycles, shown in Figure 3, main panel. At low temperatures long permutation cycles can take place, involving particles pertaining to different clusters. At higher temperatures permutation cycles longer than the number of particles K inside a single droplet, i.e. longer than 6 or 7 in Figure 3, are progressively impaired. On the other hand shorter permutation cycles still occur, indicating that individual clusters remain superfluid even at temperatures where global phase coherence is lost. Long permutation cycles have a direct bearing on the superfluid fraction, shown in the inset of Figure 3 a function of the temperature for various system sizes. Data are consistent with a superfluid transition in the Kosterlitz-Thouless universality class, [11] as expected for a two-dimensional system.

The finite superfluid response is accompanied by a strong signature in the momentum distribution n(k), shown in Figure 4. Even though Bose condensation is suppressed in two dimensions at finite temperature, in the supersolid phase the calculated momentum distribution features a pronounced peak at k = 0, as well as a smaller peak at k = 2π/l. Neither is found in the insulating crystal. The former represents the integral of the occupation number of low-momentum states, corresponding to a slow power-low decay of the one-body density matrix (OBDM). The latter describes a substantial real-space modulation of the OBDM with the periodicity of the triangular lattice. While modulations in the OBDM are expected also in the insulating case [12], they are too small to be detected within the statistical precision of this calculation.

It is interesting to consider the limit on the superfluid fraction induced by density modulations, as pointed out by Leggett [13]. The bound reads as [14]:

ρ s ≤ min ϕ(r) 1 ρV v 2 0 ρ(r) | v s (r) | 2 dx v s (r) = m ϕ(r) ϕ(r + (n L x , n L y )) = ϕ(r) + 2πn (3)
where ϕ(r) is a phase function, ρ and v 0 are the average value over the simulation cell of ρ(r) and v s (r) respectively, and n, n and n are integers. This essentially states that if the system is uniform the superfluid fraction can be one, while if density modulations are present its value must be necessarily lower. Although the Leggett bound is very loose in strongly correlated systems such as solid 4 He or parahydrogen [14], it could provide an argument to explain why soft disks are supersolid under appropriate conditions while two-dimensional solid 4 He is never, should their density profiles look different. In order to verify this possibility we compared (Figure 5) the density profile of the supersolid soft disk system for D=5, T=0.1, µ=45, ρ s = 0.25 ± 0.01 with the one of solid 4 He in two dimensions at low temperature slightly above the melting density (ρ = 0.0765 Å-1 , T = 1K). Actually the two density profiles are remarkably similar, and the Leggett bound gives ρ s ≤ 0.39 for the soft disk and ρ s ≤ 0.37 for Helium. This shows that, within the constraint of Eq. 3 imposed by one-body properties, many-body effects have even qualitative effects on the superfluid response. In particular, the difference between hard-and soft-core systems is entirely encoded in two-body (and possibly In the former condition the system is supersolid, in the latter is a normal cluster solid.

higher) correlations. Indeed the pair distribution function of supersolid soft disks is peaked at contact [7], as opposed to the correlation hole induced by hard-core repulsion in Helium.

We finally discuss the dynamic structure factor S(k, ω) of the system, which is related to the elementary excitations such as phonons in the solid and rotons in the superfluid. It is not obvious in advance if the supersolid has a dominant mode of mixed character or distinct branches. The worm algoritm is a path integral Monte Carlo approach, and therefore one cannot have access to real time dynamics of the system. Instead the imaginary-time density-density correlation function can be computed

F (k, t) = T r(ρ -k e -tH ρ k e -(β-t)H ) N Z ρ k = j<N e ik•rj (4)
Z being the trace over the density matrix of the system. The F (k, t) is related to S(k, ω) by an inverse Laplace transform:

F (k, t) = dωe -tω S(k, ω) (5)
The numerical inversion of a Laplace transform is a classical ill-conditioned problem. To cope with this problem, various approaches have been developed, among which the Maximum Entropy (ME) method is probably the most used [15]. Here we resort instead to the Genetic Inversion via Falsification of Theories (GIFT) [16] approach, which offers a better resolution in the reconstructed spectra. For example, when applied to superfluid 4 He, GIFT is able to correctly separate the sharp quasi-particle excitations from the multiphonon branch, with appropriate spectral weights [16], whereas ME merges both structures [17,18]. A typical result of the GIFT inversion procedure is reported in Figure 6, for k = (0.25, 0) in units of reciprocal lattice spacing and for two different interaction stengths D at the same density r s = 0.421 and temperature T = 0.5. For D = 5 (solid line) the system is in the insulating cluster solid phase, and the spectral weight is concentrated in a single phononic peak. For D = 3 the system is in the supersolid phase, with a significant superfluid fraction (ρ s = 0.30 ± 0.01). The phononic peak is seen to shift to lower energies with respect to the insulating phase because the system is softer (see Figure 2). More interestingly, the spectrum displays a second peak with lower energy. Preliminary data on the dispersion of these spectral features suggests that both modes are acoustic. We thus conclude that a supersolid can sustain two distinct long-wavelength excitations, respectively related to the lattice dynamics and to the superfluid response.

Figure 1 .

 1 Figure 1. Qualitative low temperature phase diagram for high and low D as a function of µ. The panels show typical spatial configuration of the world lines resulting from simulations at low T , referring to the various phases. Results shown in the upper and the lower part of the figure correspond to simulations with D = 60 and D = 3, respectively.

  , r s =0.42 Cluster solid, D=5, r s =0.42

Figure 2 .

 2 Figure2. Potential felt by a test particle averaged over many configurations along the cristallographic direction[1, 0] passing trough the lattice sites, for two different heights of the interparticle potential. The minima of the potential correspond to the positions of the clusters, where the disk density is higher. The barrier between neighbouring lattice sites is increased for higher values of D, progressively supressing tunnelling probability.

Figure 3 .

 3 Figure 3. Average number of permutation cycles as a function of the permutation length for various temperatures (main panel), for D = 5, µ = 45, rs = 0.522. The average number of particles per droplet is between 6 and 7 in these simulations, therefore the number of permutation cycles longer than this value drops at all temperatures. However at low T permutation cycles involving a large number of droplets become possible: this allows for dissipation-less flow in the system. In the inset the superfluid fraction ρs as a function of the temperature is shown, for two different system sizes, containing Nc = 16 or Nc = 30 cristalline sites.

Figure 4 .

 4 Figure 4. Momentum distribution n(k) for D = 5, rs = 0.522, µ = 45 in the supersolid phase (T = 0.2, red disks) and in the insulating cluster solid (T = 4, blue circles). The peaks at k = 0 and k = 1 (in units of 2π/l) are due to (quasi)condensation and strong modulation of the OBDM with the lattice periodicity, respectively.

Figure 5 .

 5 Figure 5. Similarity of the density profiles ρ(r)/ρ for soft disks in the supersolid phase (D=5, T=0.1, µ=45,ρs = 0.25 ± 0.01), and for solid4 He near melting (ρ = 0.0765 Å-1 , T = 1K). These density profiles are obtained along the crystallographic direction[1, 0] passing trough a lattice site. r is in units of the lattice parameter.

Figure 6 .

 6 Figure 6. S(k, ω) obtained with the described inversion procedure, for k = (0.25, 0) in units of reciprocal lattice spacing. Data are obtained for two values of D=3 and D=5 at rs = 0.421, for T = 0.5. In the former condition the system is supersolid, in the latter is a normal cluster solid.


M Chalony

J Barré

B Marcos
email: bruno.marcos@unice.fr

† A Olivetti

D Wilkowski
email: david.wilkowski@ntu.edu.sg

Long-range 1D gravitational-like interaction in a neutral atomic cold gas

Keywords: numbers: 37.10.De, 05.20.Jj, 04.80.Cc, 37.10.Gh, 04.40.-b

A quasi-resonant laser induces a long-range attractive force within a cloud of cold atoms. We take advantage of this force to build in the laboratory a system of particles with a one-dimensional gravitational-like interaction, at a fluid level of modeling. We give experimental evidences of such an interaction in a cold Strontium gas, studying the density profile of the cloud, its size as a function of the number of atoms, and its breathing oscillations.

I. INTRODUCTION

When interactions between the microscopic components of a system act on a length scale comparable to the size of the system, one may call them "long-range": For instance, the inverse-square law of the gravitational force between two point masses which is one of the most celebrated and oldest laws in physics. In the many particles world, it is responsible for dramatic collective effects such as the gravothermal catastrophe [1] or the gravitational clustering which is the main mechanism leading to the formation of the structure of galaxies in the present universe. Beyond gravitation, such long-range interactions are present in various physical fields, either as fundamental or as effective interactions: In plasma physics [START_REF] Elskens | Microscopic Dynamics of Plasmas and Chaos[END_REF], two dimensional (2D) fluid dynamics [START_REF] Chavanis | Dynamics and Thermodynamics of Systems with Long-Range Interactions[END_REF], degenerated quantum gases [START_REF] O'dell | [END_REF], ion trapping [5], to cite only these works. long-range interactions deeply influence the dynamical and thermodynamical properties of such systems. At the thermodynamic equilibrium, long-range interactions are at the origin of very peculiar properties, especially for attractive systems: The specific heat may be negative; canonical (fixed temperature) and microcanonical (fixed energy) ensembles are not equivalent. These special features have been known for a long time in the astrophysics community, in the context of self gravitating systems.

After the seminal works of Lynden-Bell and Wood [6] and Thirring [7], many contributions followed on this subject (see, for instance, [8] for a recent review), so that the equilibrium characteristics of attractive longrange interacting systems are theoretically well established. This situation is in striking contrast with the experimental side of the problem: There is currently no controllable experimental system exhibiting the predicted peculiarities. There have been some proposals to remedy to this situation: O'Dell et al. [START_REF] O'dell | [END_REF] have suggested creating an effective 1/r potential between atoms in a Bose-Einstein condensate using off-resonant laser beams; more recently, Dominguez et al. [9] have proposed taking advantage of the capillary interactions between colloids to mimic two-dimensional gravity, and Golestanian [10] has suggested experiments using thermally driven colloids. However, these proposals have not been implemented yet, and so far the dream of a tabletop galaxy remains elusive.

The key results of this paper are to show some experimental evidences of a gravitational-like interaction in a quasi-one-dimensional (hereafter 1D) test system consisting in a cold gas of Strontium atoms in interaction with two contra-propagating quasi-resonant lasers. To our knowledge, it is the first experimental realization of the 1D gravitational toy model, which can be compared with the theoretical predictions developed for more than 50 years by the astrophysical and statistical physics community. In the stationary regime, the cloud spatial distribution is in agreement with the well-known sech 2 law for the 1D self-gravitating gas at thermal equilibrium [11]. Moreover, the long-range attractive nature of the force is confirmed studying the cloud's size dependency as a function of the number of atoms. Out of equilibrium, the breathing oscillation frequency increases with the strength of the interaction as it should be for attractive interactions. Quantitatively our experimental results are in agreement with the expected 1/r α force with α = 0.

The paper is organized as follows. In Sec. II, we start from the radiation pressure exerted by the lasers and explain under which circumstances this force becomes analog to a 1D gravitational force. We then make some definite theoretical predictions on the size, density profile, and oscillation frequency of the interacting atomic cloud. The experimental setup is described in Sec. III. In the same section, the experimental results are compared with the theory.

II. MODEL AND THEORETICAL PREDICTIONS

The gravitational potential U (r) between two particles can be expressed through the Poisson equation ∇ 2 U (r) = A D Gmδ(r), where G is the coupling constant, m the mass of the particle and A D a numerical constant which depends on the dimension. The solution of the Poisson equation for the interpaticle potential U (r) in three dimension is the well-known

U (r) = Gm r , (1)
and in 1D,

U (r) = Gm|r| (2)
(for a review on 1D gravitational systems see, e.g., [12]). After using a mean-field approach (see below), we will show that such a potential should be at play in our experiment, under precise circumstances (see section II B).

We start considering a quasi-1D {cold atomic gas + 1D quasi-resonant laser beams} system; an atomic gas, with a linear density n(z), is in interaction with two contrapropagating laser beams. The two beam intensities I + (z) and I -(z), where I + (-∞) = I -(+∞) ≡ I 0 , respectively propagating in the positive and negative direction, are much smaller than the atomic line saturation intensity I s . Thus the atomic dipolar response is linear. The radiation pressure force of the lasers on an atom, having a longitudinal velocity v z , is given by [13]:

F ± (z, v z ) = ± k Γ 2 Γ 2 4(δ ∓ kv z ) 2 + Γ 2 I ± (z) I s , (3)
where is the reduced Planck constant, Γ the bare linewidth of the atomic transition, k the wave number, and δ the frequency detuning between an atom at rest and the lasers. For a cloud of N atoms, the attenuation of the laser intensity is given by:

dI ± = ∓ σ ± 2πL 2 ⊥ N I ± n(z)dz, (4)
where n(z) is the normalized linear density profile and

σ ± = 6π k 2 Γ 2 g(v z ) 4(δ ∓ kv z ) 2 + Γ 2 dv z (5)
is the average absorption cross-section for a single atom. g(v z ) is the normalized longitudinal velocity distribution and 2πL 2 ⊥ is the transverse section of the cloud. At equilibrium g(v z ) is an even function so σ -= σ + ≡ σ. The optical depth is defined as:

b = σ 2πL 2 ⊥ N +∞ -∞ n(z)dz = σN 2πL 2 ⊥ . (6)
Atoms also experience a velocity diffusion due to the random photon absorptions and spontaneous emissions:

This is modeled by a velocity diffusion coefficient D introduced in Eq. (7). In experiments, δ < 0 such that the force, given in Eq. ( 3), is a cooling force counteracting the velocity diffusion. We now describe the N atoms by their phase space density in 1D, f (z, v z , t). As in [14], we write a Vlasov Fokker-Planck equation

∂f ∂t + v z ∂f ∂z -ω 2 z z ∂f ∂v z + 1 m ∂ ∂v z [(F + (z, v z ) + F -(z, v z ))f ] = D ∂ 2 f ∂v 2 z . (7)
which is, for most of the cases, a reasonable modeling of long-range force systems in the mean-field approximation (see, e.g., [START_REF] Binney | Galactic Dynamics[END_REF]). The second term in Eq. ( 7) is an inertial one, whereas the third one describes a harmonic trapping force being a good approximation of the dipolar trap used in the experiment [START_REF] Chu | [END_REF]. Indeed the dipolar potential, in the longitudinal axe of interest, can be written as

U dip (z) = -U 0 1 + z zR 2 (8)
with z R = 1.2 mm, U 0 = 1 2 k B T trap , and T trap = 20µK. The observed rms longitudinal size being L z 400µm, it is reasonable to perform a Taylor expansion around z = 0 to get the harmonic approximation:

U dip (z) ≈ -U 0 1 - z z R 2 . (9)
having a characteristic frequency

ω z = k B T trap mz 2 R 1/2 . ( 10
)
The fourth term of Eq. ( 7) contains the mean-field force F ± divided by the atomic mass m. The right hand side describes a velocity diffusion. The use of a one dimensional model is justified by the fact that the ratio between the rms transverse L ⊥ and longitudinal L z size of the cloud measured in the experiment is

L ⊥ /L z ≈ 2 × 10 -2 .
Eq. ( 7) neglects atomic losses and dependencies in position and velocity of the velocity diffusion coefficient. One notes that the attractive force coming from the beams absorption (Eq. ( 3) and ( 4)) is known since the early days of laser cooling and trapping [17]. However, in an usual 3D setting this attractive force is dominated by the repulsive force due to photons reabsorption [18], which, in the small optical depth limit, may be seen as an effective repulsive Coulomb force. By contrast, in a 1D configuration with an elongated cloud along the cooling laser beams, the probability of photons reabsorption is reduced by a factor of the order of L ⊥ /L z , in comparison with the isotropic cloud having the same longitudinal optical depth. In our experiment, the reduction factor is about 2 × 10 -2 , so that the repulsive force can be safely ignored. Similar but weaker reduction of the probability of photons reabsorption is also expected for the 2D geometry, which opens the possibility of experimental systems analogous to 2D self-gravitating systems.

A. Fluid approximation

In order to solve Eq. ( 7) we assume that the system can be described using a fluid approach: The velocity distribution at time t does not depend on the position, except for a macroscopic velocity u(z, t). We write then the one point distribution function f as

f (z, v z , t) = mN n(z, t) 1 ∆(t) g v z -u(z, t) ∆(t) . ( 11
)
The velocity distribution g(v z ) is even, centered around u; the velocity dispersion is characterized by a time modulation ∆(t). Integrating Eq. ( 7) over dv z and over v z dv z , we obtain the fluid equations:

∂n ∂t + ∂ ∂z (nu) = 0 (12) ∂(nu) ∂t + ∂ ∂z u 2 + ∆(t) 2 v 2 z g(v z )dv z n + ω 2 z zn - 1 m n (F + + F -)g v z -u(z, t) ∆(t) ∆(t)dv z = 0. (13)

B. Stationary solution

We first look for a stationary solution; this imposes u = 0 and ∆ = 1. Eq. ( 12) is then automatically satisfied; Eq. ( 13) for the stationary density n(z) reads:

v2 z ∂n ∂z + ω 2 z zn - 1 m n [F + + F -]g(v z )dv z = 0 ,( 14
)
where we have used the notation v 2 z g(v z )dv z = v2 z . Eq. ( 4) is easily integrated, yielding:

I + (z) = I 0 e -b z -∞ n(s)ds (15) I -(z) = I 0 e -b +∞ z n(s)ds (16)
The exponentials are expanded up to first order, according to the small optical depth hypothesis b ≪ 1:

I + (z) ≃ I 0 1 -b z -∞ n(s)ds (17)
I -(z) ≃ I 0 1 -b +∞ z n(s)ds . ( 18
)
Introducing these expressions for I ± into Eq. ( 14), we obtain finally

v2 z ∂n ∂z + ω 2 z zn -N Cn +∞ -∞ sgn(s -z)n(s)ds = 0 , (19)
where

C = 3 Γ 2mkL 2 ⊥ I 0 I s σ k 2 6π 2 . (20)
Eq. ( 19) is equivalent to an equation describing the stationary density of an assembly of N trapped particles of mass m, with gravitational coupling constant G, in an external harmonic trap of frequency ω z , in a heat bath at temperature T , with the correspondence:

v2 z ↔ k B T m (21a) C ↔ Gm , (21b)
where k B is the Boltzmann constant. Two characteristic lengths are identified;

L ni = k B T mω 2 z ( 22
)
is the characteristic size of the non-interacting gas in its external harmonic holding potential. Using Eq. ( 10) we get

L ni = T T trap z R . (23)
The other characteristic length L i is associated with the interaction strength:

L i = k B T N Cm . (24)
Using these notations we write Eq. ( 19) as

∂n ∂z + zn L 2 ni - n L i +∞ -∞ sgn(s -z)n(s)ds = 0 . ( 25
)
The first term of (25) favors the density spreading. In contrast with the 2D and 3D cases, it always prevents the collapse of the cloud [19]. The second term describes an external harmonic confinement coming from the dipole trap in the experiment. The third term is the attractive interaction due to laser beams absorption. It corresponds to an 1D gravitational potential expresses in Eq. ( 2). If the inequality L i ≪ L ni is fulfilled, Eq. ( 25) is the one expected for a 1D self-gravitating gas at thermal equilibrium [11]. It yields the profile:

n(z) = 1 4L i sech 2 z 2L i . (26)
A generalization of Eq. ( 25) is written as:

∂n ∂z + 1 k B T ∂U dip ∂z n-An +∞ -∞ |s-z| -α sgn(s-z)n(s)ds = 0 , (27)
including the exact form of the dipole trap (8), and the variation of the interaction exponent α of a 1/r α attractive force. This expression is used to compare theory with experiments in section III. A is a free parameter controlling the interaction strength, and thus the width of the equilibrium profile.

C. Breathing oscillations

To probe the dynamics of the system, we now go back to Eqs. (12,13), linearizing these equations with respect to u and ∆ -1: For small amplitude oscillations. One notes that this approximation is much less restrictive than linearizing with respect to the velocity v z . We then compute [F + + F -]f dv z :

[F + + F -]f dv z ≃ c 1 (I + -I -)n + c 2 (I + + I -)nu + c 3 (∆ -1)(I + -I -)n , (28)
where the constants c i involve integrations with respect to v z . The first term is the gravitational-like force, as in (19) with n(z) replaced by the time-dependent density n(z, t). The second one is a friction, which a priori depends weakly on z through I + + I -. Since I + -I - is of order b ≪ 1, the third term, of order b(∆ -1), is neglected. We assume that the dynamics is captured by a single parameter λ(t), using the ansatz [20]:

f (z, v z , t) = mN n(z/λ)g(λv z -λz). ( 29
)
When compared with (11), this amounts to assume: u = λ λ z, ∆ = 1/λ. We introduce the notations . and . 0 for the spatial average of a quantity over the density at time t and the stationary density respectively. Then

z 2 = λ 2 z 2 0 , zu = λ λ z 2 0 , u 2 = λ2 z 2 0 . (30)
We note that Eq. ( 12) is automatically satisfied by the ansatz (29). To obtain an equation for λ, we integrate Eq. ( 13) over z dz. We obtain, for λ close to 1 (small amplitude oscillations):

λ + κ λ + ω 2 (λ -1) = 0 (31)
with κ an effective friction and a breathing oscillation frequency:

ω br = ω z (3(p -1) + 4) 1 2 . ( 32
)
p measures the compression of the cloud:

p = L 2 ni L 2 z . (33)
In experiments where the effective friction is rather small, Eq. ( 32) is expected to be a fair approximation for the breathing oscillation frequency. More generally, assuming a power law two-body interaction force in the gas 1/r α , the simple relation for ω br in the weak damping limit becomes [20]:

ω br = ω z ((3 -α)(p -1) + 4) 1 2 . ( 34
)
This formula relates ω br to α and p, and will be used in section III D. Eq. ( 34) was derived in [20] assuming a velocity independent interaction term, which would be obtained by linearizing the radiation pressure force [START_REF] Chavanis | Dynamics and Thermodynamics of Systems with Long-Range Interactions[END_REF] in velocity. This is not a reasonable approximation in our experiments [START_REF]For our experimental data points shown in Fig[END_REF], but we have shown here that (34) is still expected to provide a reasonable approximation for the breathing frequency in the limit of small optical depth.

III. EXPERIMENTS A. Experimental setup

The sample preparation is done in the same way as depicted in [START_REF] Chalony | [END_REF]. More details about laser cooling of Strontium in a magneto-optical trap (MOT) can be found in [23]. After laser cooling, around 10 5 atoms at T ≃ 3 µK are loaded into a far detuned dipole trap made of a 120 mW single focused laser beam at 780 nm. Analysis are performed using in-situ images taken with a CCD at different instant of the experimental sequence. The longitudinal profile is obtained averaging over the irrelevant remaining transverse dimension. We directly measure the longitudinal trap frequency ω z = 6.7(0.5) Hz from relaxation oscillations of the cold cloud (see example of temporal evolutions in Fig. 1). The radial trap frequency ω ⊥ = 470(80)Hz is deduced from cloud size measurements. The beam waist is estimated at 23( 2) µm leading to a potential depth of T trap ≃ 20 µK.

50 ms after loading the dipole trap (corresponding to t = 0 in Fig. 1), a contra-propagating pair of laser beams, red-detuned with respect to the 1 S 0 → 3 P 1 intercombination line at 689 nm (radiative lifetime: 21 µs), is turned on for 400 ms. These beams, aligned with respect to the longitudinal axis of the cloud, generate the effective 1D attractive interaction. When the 1D lasers are on, we apply a B = 0.3 G magnetic bias field, for two important reasons: First, the Zeeman degeneracy of the excited state is lifted such that the lasers interact only with a twolevel system made out of the m = 0 → m = 0 transition which is insensitive to the residual magnetic field fluctuation. Second, the orientation of magnetic field bias, with respect to the linear polarization of the dipole trap beam, is tuned to cancel the clock (or transition) shift induced by the dipole trap on the transition of interest [START_REF] Chalony | [END_REF].

The temperature along the 1D laser beams, in our experimental runs, is found to be in the range of 1 -3 µK. Even at such low temperatures, and in sharp contrast with standard broad transitions, the frequency Doppler broadening kv z remains larger than Γ. As a direct consequence, the optical depth b depends on the exact longitudinal velocity distribution g(v z ) (see Eqs. ( 5) and ( 6)) which are not necessarily gaussian [START_REF] Chalony | [END_REF]. Since we measure only the second moment of the distribution g(v z ), namely: vz or T , one has enough control to assert the b ≪ 1 limit, thus the occurrence of the self-gravity regime. However we can perform only qualitative tests of our theory described in section II B.

At t = -50ms, the MOT cooling laser beams are turned off, leaving the trapped atomic cloud in an outof-equilibrium macroscopic state. Without the 1D lasers, we observed a weakly damped oscillation of the breathing mode and of the center of mass position (blue circles in Fig. 1). One notes that damping is caused by anharmonicity of the dipole trap and not by thermalization of the gas which is negligible on the experimental timescale. In presence of the 1D laser beams, overdamped or underdamped oscillations of the cloud are observed.

B. Stationary state's density profile

Let us first consider the stationary state in the overdamped situation (red circles in Fig. 1). After the transient phase (t < 30 ms), the rms longitudinal size of the atomic gas reaches a plateau at a minimal value of L z ≃ 120 µm with T ≃ 2 µK. The slow increase of the cloud's size after the plateau (t > 150 ms) goes with an increase of the temperature up to 4 µK at the end of the time sequence. The origins of the long time scale evolution are not clearly identified, but it is most likely due to coupling of the longitudinal axis with the uncooled transverse dimensions because of imperfect alignment of the 1D laser beams with the longitudinal axis of the trap and nonlinearities of the trapping forces. At the plateau where temperature is around 2 µK the noninteracting gas is expected to have a rms longitudinal size of L z = L ni ≃ 370 µm. Hence, a clear compression of the gas by a factor of three is observed. It is due to the attractive interaction induced by the absorption of the 1D laser beams. Moreover, the estimated optical depth is b < 0.6. We then approach the two previously mentioned conditions -b ≪ 1 and L z ≪ L ni -for being in the 1D self-gravitating regime as discussed in section II B. In Fig. 2, where b ≃ 0.4, we test the effective interaction in the gas by assuming a power law two-body interaction force in the gas 1/r α and fitting the experimental linear density distribution for different values of α in the presence of a dipolar trap; α = 0 corresponds to 1D gravity. We see that the best fit seems to be for α ∈ [0, 1/2].

In absence of the 1D laser beams, we have checked that the experimental linear density distribution have the expected profile of a non interacting gas in our dipole trap having a z R = 1.2(1) mm Rayleigh length.

C. Cloud's longitudinal size

In the self-gravitating regime a 1/N dependency of L z is expected at fixed temperature (see (26) and the definition of L i ). Fig. 3 shows that the cloud's size L z is in agreement with this prediction for two temperature ranges: 1.5(2) µK (blue circle) and 2.1(2) µK (red star). Fits correspond to the blue dashed line for 1.5(2) µK and the red dashed line for 2.1(2) µK. The fitting expression is

N = a(1/L z -L z /L 2 ni ) , (35)
where a and L ni are free parameters depending on the temperature of the gas. If L ni ≫ L z , the self-gravitating regime is recovered in the fitting expression. However Eq. (35) takes into account the presence of a harmonic trap. Eq. ( 35) can be simply derived using the generalized virial theorem (see Eq. ( 11) in Ref. [24]) and it is in perfect agreement with numerical integrations of Eq. ( 25). The fits give L ni ≃ 0.5 mm, slightly larger than the expected value of L ni at these temperatures. The 1/N dependency of L z in the self-gravitational regime is consistent with a long-range interaction with α = 0. Unfortunately as discussed above, the residual Doppler effect prevents a quantitative comparison with the pre-diction of our model.

D. Breathing oscillations

Let us now consider the evolution of the trapped cold cloud in the underdamped situation (as an example see green circles in Fig. 1). Without the 1D lasers, the ratio of the eigenfrequencies of the breathing mode ω br and the center of mass ω z is found to be close to two, as expected for a non-interacting in a harmonic trap. As an example blue curve, shown in Fig. 1, gives ω br /ω z = 1.9 (1). If now the attractive long-range interaction is turned on, ω br is expected to follow Eq. (34) whereas ω z should remain unchanged.

Fig. 4 summarizes the comparisons between the measured ratio (ω br /ω z ) 2 and the predictions deduced from the relation (34). p is computed from the experimental data in the stationary state. We expect α = 0, however to judge the nature of the long-range attractive interaction, three plots respectively for α = 0, 1 and 2 are shown. If the α = 2 case can be excluded, the experimental uncertainties does not allow to clearly discriminate between α = 0 and α = 1. In conjunction with Fig. 2, we conclude that the system is reasonably well described by a gravitational-like interaction, α = 0.

IV. CONCLUSION AND PERSPECTIVES

In this paper, we give strong indications of an 1D gravitational-like interaction in a Strontium cold gas induced by quasi resonant contra-propagating laser beams. First, we show that in the self-gravitating limit, the density distribution follows the theoretically expected profile. Moreover, the scaling of the cloud size with the number of atoms follows the predicted 1/N law. Finally, the modification of breathing frequency of the cloud, due to the long-range interaction, is correctly described by a self-gravitating model.

Other phenomena can also be investigated e.g. in relation with plasma physic; Landau damping should be observed studying the return to equilibrium of the system after various perturbations. Moreover, the actual experimental system could be easily extended to 2D geometry suggesting interesting consequences: By contrast with the 1D case, a 2D self-gravitating fluid undergoes a collapse at low enough temperature, or strong enough interaction. Hence, it is conceivable that an experiment similar to the one presented in this paper, in a pancake geometry, would show such a collapse [25].

FIG. 1 :

 1 FIG. 1: Upper part: (Color online) Typical temporal evolutions of Lz the rms longitudinal size of the atomic cloud for three different 1D beam intensities. The laser detuning is δ = -5Γ for all curves. Lower part: The center of mass (Cdm) position of the atomic gas without the 1D lasers (I = 0). The y axis origin is arbitrary.

1 FIG. 2 :

 12 FIG.2:(Color online) Density linear distribution for N = 10 5 . The black circles are the experimental data with I = 0.02Is, δ = -6Γ and b ≃ 0.4. The profiles were symmetrized to improve the signal to noise ratio. The curves are least square fits of the data using Eq. (27) containing the exact form of the dipole trap and a two-body interaction force 1/r α . The fits are performed for each α by fixing the normalization and varying the interaction strength.

FIG. 3 :

 3 FIG. 3: (Color online) Dependency of the longitudinal size of the cloud with the number of atoms for δ = 5.7(5)Γ and I = 0.3Is. The blue circle (red star) data points correspond to temperature 1.5(2) µK (2.1(2) µK). The optical depth is in the range of 0.6-0.2 according to atoms number variations. The blue and the red dashed lines are fits using Eq. (35).

FIG. 4 :

 4 FIG. 4: (Color online) Comparison for α = 0, 1 and 2 of the experimental ratio (ω br /ωz) 2 and the predictions deduced from the relation (34). The values of p are measured on the experiment.

Acknowledgements: This work is partially supported by the ANR-09-JCJC-009401 INTERLOP project. DW wishes to thank Frédéric Chevy for fruitful discussions.