Iacopo Carusotto

Yvan Castin

Non-equilibrium and local detection of the normal fraction of a trapped two-dimensional Bose gas

Keywords: numbers: 67.85.-d, 47.37.+q, 37.10.Vz, 42.50.Gy

We propose a method to measure the normal fraction of a two-dimensional Bose gas, a quantity that generally differs from the non-condensed fraction. The idea is based on applying a spatially oscillating artificial gauge field to the atoms. The response of the atoms to the gauge field can be read out either mechanically from the deposited energy into the cloud, or optically from the macroscopic optical properties of the atomic gas. The local nature of the proposed scheme allows one to reconstruct the spatial profile of the superfluid component; furthermore, the proposed method does not require having established thermal equilibrium in the gas in the presence of the gauge field. The theoretical description of the system is based on a generalization of the Dum-Olshanii theory of artificial gauge fields to the interacting many-body context. The efficiency of the proposed measurement scheme is assessed by means of classical field numerical simulations. An explicit atomic level scheme minimizing disturbing effects such as spontaneous emission and light-shifts is proposed for 87 Rb atoms.

I. INTRODUCTION AND SUMMARY OF MAIN RESULTS

One of the most striking features of degenerate Bose gases in two dimensions is the possibility of having a superfluid behavior in the absence of a macroscopically populated Bose-Einstein condensate. The transition to the superfluid state is of the Berezinskii-Kosterlitz-Thouless (BKT) type, characterized by a sudden jump of the superfluid density from 0 to the universal value 4 (in units of the inverse square of the de Broglie thermal wavelength), independent from the details of the system [1,2]. At the transition point, the asymptotic behavior of the field correlation function changes from an exponential to a power-law decay at large distances. In contrast to the three-dimensional case, superfluidity is then not related to the appearance of a macroscopically occupied Bose-Einstein condensate in the thermodynamic limit.

Pioneering experiments have addressed the mechanical properties of two-dimensional layers of liquid Helium adsorbed on a substrate [3] and have characterized the universal jump of the superfluid fraction at the BKT critical point. On the other hand, liquid Helium experiments have limited access to the momentum distribution and the correlation functions of the fluid. The situation of ultracold atom experiments is almost the opposite: evidence of the BKT transition has been obtained from the coherence functions [4], the number of observed vortices [5], and the density profile after time-of-flight [6], while the macroscopic mechanical properties of the fluid have not been characterized yet.

Quite some effort has been recently devoted to the conceptual problem of how to experimentally detect genuine superfluidity in a quantum gas of ultracold atoms and not simply Bose-Einstein condensation [7]. A possibility explored in [8] is to look at the response of a gas in a toroidal trap to a static azimuthal artificial gauge field: a spectroscopic signature is proposed which should provide direct information on the total superfluid mass of the system. A different strategy proposed in [9] consists of looking at the evolution of the equilibrium density profile of a trapped gas when it is set into rotation.

In the present paper we propose two experimental protocols to measure the normal fraction of a gas in a local way, so to extract its spatial dependence in a trapped geometry. This feature is most relevant for atomic samples, as the superfluid core coexists with an external ring of normal gas [10]. In particular, the proposed diagnostic technique does not require to relate experimental observations after time of flight to in-trap quantities. Furthermore, in contrast to [9] our technique does not require thermodynamic equilibrium in the gas in presence of rotation [11] and may be applied to more general, nonequilibrium conditions.

The basic idea of our proposal is based on the definition of normal and superfluid fractions of a quantum fluid in terms of its current response to a transverse gauge field in the low-frequency and long-wavelength limit [START_REF] Pines | The Theory of Quantum Liquids[END_REF][START_REF] Dalfovo | [END_REF]. A spatially oscillating artificial gauge field [17][18][19][20] with a spatially localized envelope can be applied to the atomic gas using a suitable combination of laser beams. The response of the fluid to the gauge field can be detected either mechanically or optically. In the former case, one has to measure the amount of energy that is deposited in the atomic gas at the end of a suitable temporal sequence of gauge field. In the latter case, one can observe e.g. the phase shift that is experienced by the laser fields while crossing the atomic cloud.

The structure and the main results of the paper can be summarized as follows. In Sec.II we review the definition of the normal and superfluid fractions that we adopt throughout the whole paper. A strategy to generate the artificial gauge field with the suitable spatial geometry is presented in Sec.III using three laser beams, namely a coupling beam and two probe ones.

The first method to measure the normal fraction is by mechanical means. It is discussed in Sec.IV: a pulse of spatially modulated gauge field is suddenly applied to the gas and then slowly switched-off according to an exponential law in time. An analytical calculation within the linear response theory and local density approximation shows that the energy that is deposited in the gas at the end of the gauge field pulse is indeed proportional to the normal (total) density in the small fraction of the gas where the probe beam is focused, if the spatial modulation of the gauge field is orthogonal (parallel) to the coupling beam direction. Some of the complications that naturally occur in experiments are then numerically investigated: a Bogoliubov theory is used to assess the conditions to be imposed to the geometry and the temporal duration of the pulse; a classical field model is used to assess the conditions on the amplitude of the gauge field for the linear response theory to be valid. A main difficulty appears to be the relatively small amount of energy that can be deposited in the gas before nonlinear couplings become difficult to extrapolate out: the resulting figure is on the order of 1% of the total energy of the gas, which is however not far from the sensitivity of state-of-the-art thermodynamic measurements of the energy [21][22][23][24][25].

The second method to measure the normal fraction is by optical means. It is discussed in Sec.V A. The spatially modulated gauge field is imposed using the same laser beam configuration in a continuous wave regime to generate a stationary current pattern in the gas. This current pattern can be read out from the phase shift accumulated by the same probe laser beams after crossing the atomic cloud. Analytical calculations show that this phase shift is indeed proportional to the normal (total) density if the spatial modulation of the gauge field is orthogonal (parallel) to the coupling beam direction. For realistic configurations, the amount of the phase shift is anticipated to be of the order of a fraction of 10 -3 , i.e. small but still appreciable with present-day optical techniques.

Another optical set-up that is able to provide quantitative information on the normal fraction is discussed in Sec.V B: analytical calculations are used to relate the angular distribution of the scattered light off a single probe beam to the normal fraction of the gas. The main advantage of this last configuration is the rapidity of the measurement but the limiting factor is the relatively small number of photons (of the order of a few units) that are expected to be scattered in the useful directions.

In Appendix A we discuss how to cope with all those spurious effects that arise from a realistic configuration of atomic levels and laser fields, in particular spontaneous emission and the mechanical effect of undesired light-shifts; even though the underlying concepts are general, the discussion is mostly focused on the most promising case of 87 Rb atoms. The generalization of the Dum-Olshanii theory of artificial gauge fields to the many-body context is reported in Appendix B: this development is required to put the gauge field concepts on firm ground in presence of atomic interactions and to evaluate in a rigorous way the optical response of the atoms to the combined coupling and probe beams. Appendix C gives more details on the analytical derivation of the deposited energy and clarifies some issues related to the local density approximation. The framework for calculating the deposited energy using the Bogoliubov theory within the linear response regime is discussed in Appendix D. The last Appendix E discusses issues related to the statistical noise on the deposited energy both in the numerical calculation and in an actual experiment. In particular, it shows how useful information on the normal fraction could be extracted from the noise if a sufficiently precise determination of the initial energy was possible. We conclude in section VI.

II. DEFINITION OF SUPERFLUID AND NORMAL FRACTIONS

Our proposal to quantitatively assess the superfluidity of the two-dimensional atomic gas is based on the traditional definition of the normal fraction f n in terms of the response to a transverse gauge field coupling to the atomic current operator [START_REF] Pines | The Theory of Quantum Liquids[END_REF][START_REF] Dalfovo | [END_REF]. The Hamiltonian giving the coupling of the matter current to an arbitrary vector potential A(r) is

V = -d 2 r A(r) • j(r) (1)
with the current operator defined as usual as j(r) = 2im φ † (r)∇ φ(r) -h.c. (2) in terms of the bosonic field operator φ for the twodimensional gas. For a spatially homogeneous system, the linear response susceptibility relating the average current [27] to the applied gauge field can be easily written in momentum space and frequency domain as j (q, ω) = χ(q, ω) A(q, ω).

(3)

If the system is also invariant under reflection with respect to the direction of q, the susceptibility tensor χ(q, ω) turns out to be diagonal in the longitudinal/transverse basis with respect to q, with diagonal matrix elements χ L,T (q, ω), respectively. For a system of surface density ρ, the normal fraction f n of the system is then defined as the low-momentum, low-frequency limit of the susceptibility to transverse gauge fields:

f n = lim q→0 lim ω→0 m ρ χ T (q, ω). (4)
Note that the order of the limits is here important. A well-known sum rule based on gauge invariance imposes that the same limit for the longitudinal susceptibility χ L (q, ω) gives exactly unity,

1 = lim q→0 lim ω→0 m ρ χ L (q, ω). (5)
The definition (4) can then be extended to large but finite systems using the standard local density approximation.

III. HOW TO GENERATE THE GAUGE FIELD

We consider a three-dimensional gas of bosonic atoms in a strongly anisotropic, pancake-shaped trap. The axial confinement frequency ω z is much higher than the one ω along the xy plane; both the temperature T (times the Boltzmann constant k B ) and the chemical potential µ of the gas are assumed to be smaller than ω z . In this regime, the gas will be eventually described in terms of a two-dimensional Hamiltonian.

Building on an idea originally introduced in [17], an artificial gauge field coupling to the atomic current can be obtained by illuminating the atoms with several laser beams with suitably chosen frequencies, wavevectors, and waist profiles. Several schemes to generate artificial gauge fields for neutral atoms have been proposed in the last years [17][18][19]. The last proposal [19] was recently implemented on an atomic Bose-Einstein condensate: for sufficiently strong gauge fields, a disordered ensemble of vortices appeared in the gas [20]. . Right panel: side view. kc and k ± p are the wavevectors of the coupling and probe beams, respectively. In practice, k ± p = (k 2 p -q 2 /4) 1/2 ez ± (q/2) ex with q ≪ kp, where ei is the unit vector along axis i.

In the present paper, we shall focus our attention on the level configuration shown in Fig. 1. Three internal atomic levels in a Λ configuration are connected by three laser fields according to the sketch given in the left panel of Fig. 1: a coupling beam resonantly drives the |b → |e atomic transition, while a pair of probe beams resonantly drive the |a → |e transition with essentially the same detuning. The artificial gauge field originates from the spatial and temporal dependence of the resulting optically dark state [17]. All other atomic states are assumed to be far-off resonance. A discussion of their effect in the specific case of 87 Rb atoms is given in the Appendix A. A summary of the suggested experimental parameters for this specific atomic species is given in Table I.

The geometrical arrangement of the laser beams is sketched in the central and right panels of Fig. 1. The continuous-wave control beam propagates along the y direction with a wavevector k c and a carrier frequency ω c close to resonance with the |b → |e transition, with a detuning δ, and has a peak Rabi frequency Ω o c . Its waist profile is much wider than the size of the atomic cloud, so that it can be safely approximated by a plane wave.

The two probe beams share the same carrier frequency ω p close to resonance with the |a → |e transition. The carrier frequencies of the coupling and probe beams are chosen exactly on resonance with the Raman transition |a → |e → |b , that is ω p -ω c = ω b -ω a . The probe beams impinge on the atomic cloud with wavevectors k ± p ≃ k p e z ± q/2 close to the z direction and symmetrically located with respect to it. The difference q = k + pk - p lies along the xy plane and is in magnitude q ≪ k p = ω p /c. The probe beams spatial profile is taken to be a Gaussian with a waist w, centered at r 0 in the z = 0 plane.

In what follows, we will need the condition w ≫ q -1 , where q is small enough for the limit in Eq.( 4) to be almost reached. In the numerical examples to come, this requires qξ < 1, where ξ is the healing length of the gas. This automatically shows that the concept of normal fraction can be used only when the portion of the gas within the waist is "macroscopic" [28], here w ≫ ξ, which was physically expected. At the same time, the waist w is assumed to be much smaller than the atomic cloud radius, so to allow for a local measurement of the normal fraction. The peak Rabi frequencies of the two probe beams are Ω ± p (t), respectively (see Appendix B for the precise definition of the Rabi frequencies). The spatial dependence of the Rabi frequencies of both the coupling and the probe beams is then summarized by the following expressions [29],

Ω c (r, t) = Ω (0) c e -i∆ct e ikc•r (6)
Ω p (r, t) ≃ [Ω + p (t) e ik + p •r + Ω - p (t) e ik - p •r ] ×e -[(x-x0) 2 +(y-y0) 2 ]/w 2 . ( 7
)
where we have allowed for the coupling beam to have a small detuning ∆ c from Raman resonance on top of its carrier frequency at ω c . On the contrary, Ω p does not have a time dependent phase factor, but only contains the square root of a purely real non-negative switch-on and switch-off function f (t), common to the two probe beams:

Ω p (r, t) = Ω 0 p (r)[f (t)] 1/2 . (8)
In what follows, we shall restrict our attention to the non-saturating regime

|Ω c |, |Ω ± p | ≪ |δ + iΓ/2|
, where Γ is the decay rate of |e due to spontaneous emission, and δ is the common detuning of the probe and control beam carrier frequencies from the transitions |a → |e and |b → |e respectively. We shall also concentrate on the limit |Ω ± p | ≪ |Ω c | where the structure of the gauge field is the simplest.

As the transitions driven by the probe and control beams share the excited state |e , there exists for each spatio-temporal coordinates (r, t) an internal non-coupled state, for which the two excitation channels interfere destructively. In terms of the local Rabi frequencies Ω p (r, t) and Ω c (r, t), this non-coupled state reads

|N C(r, t) = |a -Ω p (r, t)/Ω c (r, t) |b (1 + |Ω p (r, t)| 2 /|Ω c (r, t)| 2 ) 1/2 . ( 9
)
Adiabatically eliminating the excited state |e , one sees that the bright orthogonal state, the so-called coupled state,

|C(r, t) = [Ω p (r, t)/Ω c (r, t)] * |a + |b (1 + |Ω p (r, t)| 2 /|Ω c (r, t)| 2 ) 1/2 (10) is separated from |N C(r, t) by a (complex) energy gap [δ ′ (r, t) -iΓ ′ (r, t)/2] ≡ [|Ω c (r, t)| 2 + |Ω p (r, t)| 2 ] 4(δ + iΓ/2) , (11)
where δ ′ and Γ ′ are the light-shift and the decay rate of the coupled state. If the energy gap is large enough as compared to both the motional coupling between |N C and |C due to the spatio-temporal dependence of Ω c and Ω ± p [30], and to the quantum of oscillation ω z along the tightly confined z direction, we can restrict the dynamics to the |N C internal state.

Generalizing the single-particle theory of [17] to the many-body context in Heisenberg picture, one gets to an effective Hamiltonian for the component φ3D (r, t) of the three-dimensional atomic field operator in the (spatially and temporally-dependent) non-coupled state |N C(r, t) , φ3D (r, t) = N C(r, t)|a Ψa (r, t) + N C(r, t)|b Ψb (r, t) (12) in the simple form [31]:

H = d 3 r φ † 3D - 2 ∇ 2 2m + U 3D (r) + W 3D (r, t) φ3D -j 3D (r) • A 3D (r, t) + 1 2 g(r, t) φ † 3D φ † 3D φ3D φ3D ( 13
)
where U 3D is the three-dimensional trapping potential (supposed to be common to the internal states |a and |b ) and the vector gauge potential

A 3D (r, t) = i 2 N C(r, t)|[∇|N C(r, t) ] -c.c. ( 14
)
couples to the three-dimensional atomic current operator

j 3D (r) = 2im φ † 3D (r)∇ φ3D (r) -h.c. , (15)
and the scalar potential

W 3D (r, t) = - i 2 N C(r, t)| [∂ t |N C(r, t) ] -c.c. + + 2 2m i=x,y,z [∂ ri N C(r, t)|][∂ ri |N C(r, t) ] ( 16
)
couples to the three-dimensional density

n 3D (r) = φ † 3D (r) φ3D (r). ( 17
)
The derivation of the Hamiltonian ( 13) is based on the Quantum Stochastic Differential Equations formalism [START_REF] Gardiner | Quantum Noise[END_REF]; the details are given in the Appendix B.

The spatial and temporal dependence of the weights of the non-coupled state (9) in the internal states |a and |b reflects into a similar variation of the coupling constant describing the atomic interaction within the internal state |N C :

g 3D (r, t) = |Ω c | 4 g aa + 2|Ω p (r, t)| 2 |Ω c | 2 g ab + |Ω p | 4 g bb (|Ω c | 2 + |Ω p (r, t)| 2 ) 2 , (18)
where the coupling constants g aa , g ab and g bb originate from the a -a, a -b and b -b elastic s-wave interactions. In what follows, we shall be interested in isolating the response of the system to the gauge field A 3D . To this purpose, it will be useful to minimize the effect of all unwanted couplings to the density introduced by the scalar potential W 3D and by the spatio-temporal dependence of the interaction constant g 3D . This latter effect is minimized if one chooses states a, b with similar scattering properties g aa ≃ g ab ≃ g bb . In the limit |Ω p /Ω c | ≪ 1, one simply needs to have g aa ≃ g ab , as it is assumed from now on.

The atomic field along z is assumed to be frozen in the ground state of the harmonic confinement of wavefunction φ 0 (z). This allows to express the three-dimensional bosonic field φ3D in terms of the bosonic field for a twodimensional gas, setting φ3D (x, y, z) = φ 0 (z) φ(x, y).

(

) 19
Correspondingly, the two-dimensional coupling constant g has the expression

g = 2 m g = g 3D √ 2π a z ho ( 20
)
in terms of the three-dimensional coupling constant g 3D = 4π 2 a 3D /m, where a 3D is the s-wave scattering length in the state |a , and the size a z ho = /mω z of the ground state along z. Note that the dimensional reduction (19) does not require to be in the Lamb-Dicke limit k p,c a z ho ≪ 1, it is sufficient that the typical energy per particle corresponding to the various terms in H is ≪ ω z .

The effective two-dimensional gauge and scalar potentials A and W then result from an average of the Hamiltonian (13) over the motional ground state along z. Including in ( 14) and ( 16) the explicit form of the beam profiles and restricting ourselves to zeroth order in the small parameters q/k p,c , 1/wk p,c , and to second order in Ω ± p /Ω c , the resulting two-dimensional gauge potential turns out to be directed along the y axis and to have the form

A y (r) ≃ k c Ω + p e iq•r/2 + Ω - p e -iq•r/2 2 |Ω c | 2 e -2|r-r0| 2 /w 2 = k c |Ω + p | 2 + |Ω - p | 2 + [Ω + p Ω - * p e iq•r + c.c.] |Ω c | 2 × e -2|r-r0| 2 /w 2 . (21)
To the same level of approximation, the scalar potential has the form

W (r) = 2 (k 2 c + k 2 p ) 2m + ∆ c × Ω + p e iq•r/2 + Ω - p e -iq•r/2 2 |Ω c | 2 e -2|r-r0| 2 /w 2 , (22)
which can be made to vanish by choosing a detuning ∆ c that exactly compensates the recoil of the atoms after the Raman process a → e → b:

W ≡ 0 for ∆ c = - (k 2 c + k 2 p ) 2m . ( 23
)
It is worth pointing out that the temporal derivative df (t)/dt of the real switch-on and switch-off function f (t) of the probe beam Eq.( 8) has an exactly vanishing contribution to the ∂ t term in the right-hand side of Eq.( 16), so it does not contribute to the scalar potential W to all orders in Ω p /Ω c [START_REF]For the choice (23[END_REF]. The same conclusion holds for the gauge field, since there is no temporal derivative in (14). After expansion of the squared modulus as done in the second line of (21), two kinds of terms are immediately identified: (i) a slowly varying Gaussian term of size w and peak amplitude

|Ω + p | 2 + |Ω - p | 2 that
follows the laser envelopes, and (ii) an oscillating term at wavevector q with a Gaussian envelope of size w and peak amplitude |Ω + p Ω - p |. This spatially modulated term is indeed the one that we need to probe the normal fraction of the gas according to the definition (4): when q is taken along the x axis (y axis), it provides an almost purely transverse (longitudinal) contribution to the gauge field A. On the other hand, the slowly varying term always includes both longitudinal and transverse vector field components. Experimental procedures to subtract the effect of this unwanted term will be discussed in the next sections.

IV. DEPOSITED ENERGY MEASUREMENT

A. General idea

In this section we shall present a method to extract the value of the normal fraction from a measurement of the energy that is deposited in the system by a suitably designed gauge field sequence. The coupling beam is assumed to be always switched on. On the other hand, both probe beam intensities |Ω ± p | 2 are varied in time according to the (dimensionless) real envelope function f (t). This function is chosen to be 0 for t < 0 and to rapidly tend back to 0 at long times.

As already mentioned, we assume that the atomic interaction constants satisfy g aa ≃ g ab [START_REF]The spatio-temporal dependence of g3D(r, t) = gaa + ∆g(r, t) ≃ gaa + 2(g ab -gaa) |Ωp(r, t)/Ωc| 2 is responsible for spurious excitation of the gas via two different effects that can be captured by a Bogoliubov treatment[END_REF]. As soon as qw ≫ 1, the deposited energy is the sum of two independent contributions ∆E 1,2 corresponding to the decomposition (21) of the gauge field as the sum of a nonmodulated term and a modulated one at the wavevector q. Using standard linear response theory within the local density approximation as discussed in the Appendix C, the contribution ∆E 2 of the modulated term can be written in the simplified form

∆E 2 ≃ π 4 w 2 ǫ gauge 2 2 +∞ -∞ dω π ω |f (ω)| 2 Im[χ yy (q, ω)] (24)
where we have introduced the amplitude of the spatially modulated part of the gauge field,

ǫ gauge = 2 k c |Ω + p Ω - p |(t = 0 + ) |Ω c | 2 . (25)
A similar expression for the contribution ∆E 1 of the nonmodulated term is given in the Appendix C as (C8). The susceptibility χ yy appearing in (24) is the one (as defined in section II) of the fictitious spatially homogeneous twodimensional Bose gas that approximates the state of the trapped gas around the center r 0 of the probe beams' spot.

The expression (24) for the deposited energy involves the imaginary part of that susceptibility, while the normal fraction (4) involves the real part. To relate the two, one can make use of the well-known Kramers-Kronig relation of linear response theory,

lim ω→0 Re[χ(q, ω)] = +∞ -∞ dω ′ π Im[χ(q, ω ′ )] ω ′ . ( 26
)
For a suitably chosen envelope of the form

f (t) = e -γt Θ(t), (27)
where Θ(t) is the Heaviside step function, the Fourier transform is equal to f (ω) = i/(ω + iγ), and the integral in (24) indeed reduces to the real part of the susceptibility (26) in the γ → 0 limit. As a consequence, the deposited energy ∆E 2 for small q perpendicular (parallel) to k c can be related to the normal (total) density ρ n (ρ) of the trapped gas at position

r 0 by ∆E 2 ≃ π 4 w 2 m ǫ gauge 2 2 × ρ n (r 0 ) for q ⊥ k c ρ(r 0 ) for q k c . (28)
For intermediate angles α between q and k c , ∆E 2 is proportional to ρ cos 2 α + ρ n sin 2 α [36].

In an actual experiment, the undesired contribution ∆E 1 can be eliminated by noting its independence on the relative orientation of q and k c [see Eqs.(C8) and (C9)], as well as its different dependence on the probe amplitudes Ω ± p , proportional to [

|Ω + p | 2 + |Ω - p | 2 ] 2 rather than |Ω + p Ω - p | 2
. By measuring the deposited energy for at least two different values of the Ω + p /Ω - p ratio, one is able to isolate the relevant contribution (28).

B. How fast is the q, γ → 0 limit reached ?

An important point in view of experiments is to characterize how small the wavenumber q and the switchoff rate γ of the gauge field spatial modulation have actually to be taken to obtain a quantitatively accurate measurement of the normal fraction f n . To answer this question, we consider in this subsection the simplest case of a spatially homogeneous two-dimensional system in a square box of size L with periodic boundary conditions, excited by the gauge field in a plane wave form, i.e. in the limit w → ∞. We also limit ourselves to the case of a transverse gauge field with q = q e x perpendicular to k c = k c e y , A ideal (r, t) = e -γt Θ(t)e y ǫ gauge 2 e iqx + e -iqx . (29)

The deposited energy at the end of the gauge field sequence can be evaluated by means of the Bogoliubov theory of dilute Bose gases. The main steps of the calculation are sketched in the Appendix D. The final result reads

∆E 2 = ǫ gauge 2 2 N m f eff n ( 30
)
in terms of the wavevector-and γ-dependent effective normal fraction f eff n :

f eff n = 1 N k =0,-q 2 k 2 y m × Re n k -n k+q ǫ k+q -ǫ k -i γ (U k U k+q -V k V k+q ) 2 + 1 + n k + n k+q ǫ k+q + ǫ k -i γ (U k V k+q -V k U k+q ) 2 . ( 31
)
Here N is the total particle number, ρ = N/L 2 is the surface density, n . In (a), dependence of f eff n on the switch-off rate γ of the gauge field, for three different values of the gauge field spatial modulation wavevector q = (2π/L) ex (solid, qξ ≃ 0.1), q = 7 × (2π/L) ex (dotted, qξ ≃ 0.7) and q = 10 × (2π/L) ex (dashed, qξ ≃ 1). System parameters: square box (with periodic boundary conditions) of size L/ξ ≃ 63 containing N ≃ 40000 particles with interaction constant g = 0.1 2 /m at a temperature T /T d = 0.1. kBT d = 2π 2 ρ/m is the degeneracy temperature. The healing length ξ of the gas is defined by 2 /(mξ 2 ) = ρg and the Bogoliubov sound velocity is defined as cs = (ρg/m) 1/2 . The dashed curve corresponds to the same value of qξ ≃ 1 as in Fig. 4 and the vertical dotted lines indicate the values of γ/(csq) considered in that figure. The thin horizontal line is the prediction (34) of the quantum Bogoliubov theory in the thermodynamic limit. Interestingly, for the considered atom number, the finite-size thermodynamic formula Eq.( 39) evaluated with the quantum Bogoliubov theory gives already the same value at the scale of the figure [37]. The quadratic rather than linear dependence of f eff n on γ for small values of γ is a finite-size effect. The fact that in (a), the values of f eff n for the dotted curve are close to the arithmetic mean of the solid and dashed curves, for common values of γ/(csq), suggests that f eff n is roughly an affine function of (qξ) 2 for fixed γ/(csq). This guess is substantiated in panel (b) where the dependence of the effective normal fraction on the modulation wavevector q = q ex is shown for γ/(csq) = 0.2 [the crosses are the Bogoliubov predictions and the line is a linear fit of f eff n as a function of q 2 ].

ǫ k = 2 k 2 2m 2 k 2 2m + 2ρg 1/2 (32
is the usual Bogoliubov dispersion relation and the amplitudes of the Bogoliubov modes satisfy

U k + V k = 1 U k -V k = 2 k 2 /2m 2 k 2 /2m + 2ρg 1/4 . (33)
The n k are the thermal mean occupation numbers of Bogoliubov modes,

n k = 1/[exp(ǫ k /k B T ) -1].
The thermodynamic limit L → ∞ at fixed density N/L 2 can be worked out analytically by first taking the γ → 0 limit and then the q → 0 limit in the expression (31) for f eff n . In this way one recovers the usual Bogoliubov expression for the normal fraction, which in dimension two reads:

f n = 1 ρ d 2 k (2π) 2 2 k 2 y m (-∂ ǫ k n k ). ( 34
)
For a finite size system, the dependence of f eff n on γ and q is explored in Fig. 2. For the smallest non-zero wavevector value allowed by the chosen quantization box, the relative error on f eff n is already as small as 10% for γ/c s q = 0.15, where c s = (ρg/m) 1/2 is the Bogoliubov sound velocity. For γ/(c s q) fixed, in particular to the relevant value 0.2, and over the relevant range qξ ≤ 1, we have also found that f eff n is to a good approximation an affine function of (qξ) 2 , see Fig. 2b, which was not obvious.

Another interesting result of Bogoliubov theory applied to our system is a sufficient condition on the amplitude of the gauge field to be within the linear response regime. To this purpose, we can write the equations of motion for the Bogoliubov mode operators b k in the interaction picture in the presence of the timedependent gauge field, and impose that the amplitude change is small as compared to the initial value. The expression of this amplitude change, that we shall not give here, involves as usual energy denominators such as ǫ k -ǫ k±q + i γ.

In the limit γ → 0, the condition of linear response is thus most stringent for modes of wavevector k such that ǫ k = ǫ k±q , where the real part of the energy denominator can vanish. For the maximal value of k set by the thermal occupation, this leads to the sufficient condition

ǫ gauge (mk B T d ) 1/2 T d T 1/2 2 γ k B T d , (35)
in terms of the degeneracy temperature k B T d = 2π 2 ρ/m. This naive argument is however not able to determine to which extent this condition is actually necessary. This would require a higher order calculation which falls beyond the scope of the present work [47].

C. Numerical investigation

To further assess the validity and accuracy of our proposed scheme we have performed full scale numerical simulations of the response of a two-dimensional Bose gas at finite temperature to the complete gauge field (21), including the Gaussian envelope of the gauge field and a circular well trapping potential. A very useful tool to this purpose is the classical field model developed and applied in a number of recent works [14,38]. For this model both the thermal equilibrium state and the temporal dynamics can in fact be easily addressed with standard numerical tools and provide reliable results for the physics of the degenerate Bose gas.

We consider a classical (c-number) complex field defined on a square grid. The real space lattice constant b is chosen in terms of the thermal de Broglie wavelength λ = 2π 2 /mk B T as b/λ = π/(4ζ). The specific value ζ ≃ 0.80 [START_REF]G/π -2 ln 2, G being the Catalan constant G ≃ 0[END_REF] of the numerical coefficient is such that the classical field model correctly reproduces the total number of non-condensed particles for an ideal gas at zero chemical potential in the thermodynamic limit. This choice corresponds to setting the ultra-violet momentum cut-off

k max = π/b at 2 k 2 max /2m = ζk B T in the classical field theory.
In the canonical ensemble, the thermal probability distribution for the interacting classical field follows a Boltzmann law δ(

Ψ 2 -N ) exp(-E[ψ]/k B T ) with the norm- squared Ψ 2 = b 2
r |Ψ(r)| 2 fixed to the total atom number, and the discrete Gross-Pitaevskii energy functional [40] given by:

E[Ψ] = b 2 r Ψ * - 2 2m ∆ + U (r) + g 2 |Ψ| 2 Ψ, (36)
where U (r) is the trapping potential seen by the bidimensional gas. This probability distribution can be sampled by the long time limit of a Ito stochastic differential equation with drift terms and a noise term [START_REF] Mandonnet | [END_REF][START_REF] Sinatra | [END_REF], including projectors in order to keep the norm-squared constant, Ψ 2 = N :

dΨ = - 1 2 dτ Q Ψ - 2 2m ∆Ψ + U Ψ + g|Ψ| 2 Ψ + √ k B T b Q Ψ dξ - M -1 2N k B T dτ Ψ, ( 37
)
where M is the number of grid points, Q Ψ is the projector onto the subspace orthogonal to the classical field Ψ, dξ is a complex Gaussian, zero-mean, delta-correlated and temporally white noise such that dξdξ = 0 and

dξ * (r i ) dξ(r j ) = dτ δ ri,rj . (38)
Our numerical procedure simply consists in first generating a number n real of independent wavefunctions distributed according to the thermal Boltzmann law with energy functional (36) and then to let them evolve in real time according to the discrete time-dependent Gross-Pitaevskii equation including the gauge potential (21) [43].

As a first application of the classical field method, we have determined the normal fraction of a spatially homogeneous two-dimensional interacting gas at thermal equilibrium. This is done using the thermodynamic formula, FIG. 3: Classical field simulation of the normal fraction fn for a finite size two-dimensional, spatially homogeneous interacting Bose gas as a function of temperature. The calculation has been performed using the thermodynamic expression [START_REF]G/π -2 ln 2, G being the Catalan constant G ≃ 0[END_REF] for the normal fraction. The coupling constant is g = 0.1 2 /m. The different curves refer to simulations performed with different system sizes L, resulting in different numbers of grid points M = 16 2 (green, dotted), 32 2 (red, dashed), 64 2 (black, solid) [the grid spacing is b ≃ 0.99λ, where λ is the thermal de Broglie wavelength, see text]. The number of classical field realizations is n real = 1000. The dot-dashed line is the classical field prediction for the normal fraction fn of an ideal gas in the thermodynamic limit: as discussed in the text, the decrease at high temperatures is an artifact of the classical field model. The blue circle indicates the result of a numerical simulation of the deposited energy scheme as in Fig. 4 for a homogeneous system with an infinite beam waist w = ∞ (in which case we simply dropped the non-modulated term in the gauge field), a number of grid points M = 32 2 , a gauge field modulation wavevector q = (2π/L) ex ≃ (0.2/ξ) ex, and including a careful extrapolation of ǫgauge → 0 (numerics down to ǫgauge/(mkBT d ) 1/2 = 0.01) and of γ → 0 (numerics down to γ/(csq) = 0.05 where cs = (ρg/m) 1/2 is the Bogoliubov sound velocity).

applicable to finite size systems with periodic boundary conditions,

f thermo n = P 2 y N mk B T (39)
which involves the thermal variance of the total momentum P y of the gas [44,45]. The temperature dependence of the normal fraction is shown in Fig. 3 for increasing system sizes. The sudden variation around T /T d ≃ 0.13 becomes sharper and sharper as the system size is increased and should eventually correspond to a discontinuous jump in the superfluid fraction at the BKT transition [1,2,46]. The slow decrease for larger values of T /T d is instead an artifact of the ultraviolet cut-off that has to be imposed to the classical field model in any dimension d ≥ 2. Indeed, the same decrease is visible also in the case of an ideal gas, for which one can show that

f n = 1 -T /T d + O(e -T d /T T /T d ) in the thermodynamic limit.
The experimental estimation of the normal fraction obtained by the deposited energy method discussed in Sec.IV A is simulated in Fig. 4. The value of the deposited energy is extracted from the classical field simulation by taking the energy difference at the end of two evolutions using the same value of

|Ω + p | 2 +|Ω - p | 2 but different relative magnitudes of Ω ± p , Ω - p = 0 and |Ω + p | = |Ω - p | respectively.
This protocol aims at isolating the effect of the spatially modulated gauge potential: in the linear response limit, it is able to provide the exact value of ∆E 2 alone. The effective normal fraction is then extracted from the deposited energy via (28), and it is plotted in Fig. 4, as a function of the gauge field amplitude, for different values of the switch-off rate γ/(c s q) = 0.4, 0.2, 0.1 (in black, in red, in blue, from bottom to top) [c s = (ρg/m) 1/2 is the Bogoliubov sound velocity] and different geometries (empty squares and dashed lines for the spatially homogeneous case with periodic boundary conditions, filled circles and solid lines for the circular potential well).

The ǫ gauge -dependence in that figure allows to estimate the interval where the linear response approximation is reasonable, e.g. for ǫ gauge / √ mk B T d 0.08 the deviation due to non-linear effects is less than five percent. As expected, the estimate [START_REF] Stöferle | [END_REF] gives a more pessimistic bound around 0.04 for γ/(c s q) = 0.4. From an experimental point of view, we expect that values of the gauge field amplitude as high as ǫ gauge /(mk B T d ) 1/2 = 0.15 should be well achievable, see Table I. As shown in the figure, for ǫ gauge /(mk B T d ) 1/2 0.2, the effective normal fraction f n empirically has, within the error bars, a linear dependence with ǫ 2 gauge . This could be expected from the fact that the first order correction to the linear response theory scales as ǫ 4 gauge in the deposited energy ∆E, the term of order ǫ 3 gauge in ∆E being zero due to the odd parity of the matter-wave coupling to the gauge field. Assuming this linear dependence of the effective normal fraction, a fit was performed over the range (0.05) 2 ≤ ǫ 2 gauge /(mk B T d ) < (0.2) 2 , to extrapolate to ǫ gauge = 0, see the solid and dashed lines in Fig. 4. Note that this extrapolation procedure also eliminates the contribution to the effective normal fraction of the energy ∆E U p deposited by the undesired light-shifts, see the end of Appendix A. For larger values of the gauge field amplitude ǫ gauge , more serious non-linear effects set in in the effective normal fraction and it becomes less clear how to extrapolate them down to ǫ gauge = 0.

For the sake of completeness, it is important to note that for the weak gauge field amplitude that are required to be in the linear regime, or at least in the regime of the numerically suggested linear dependence of the effective normal fraction with ǫ 2 gauge , the deposited energy ∆E 2 is at most a few percent of the total energy of the system, which may be experimentally challenging to measure [49]. This value is however larger than the statistical uncertainty of the energy in the canonical ensemble FIG. 4: Classical field simulation of the deposited energy measurement scheme. Initial thermal equilibrium state at T /T d = 0.1 with g = 0.1 2 /m, and a system size L/ξ ≃ 63 (the number of grid points is M = 64 2 , and n real = 1000 independent realizations of the classical field are used). Gauge field modulation wavevector q ≃ (1/ξ) ex. The real time evolution is followed during τ = 3/γ. The complete form Eq.( 21) of the gauge field is considered, with r0 located in the center of the system, and a beam waist w = 30ξ. For each value of the gauge field intensity ǫgauge of Eq.( 25), two calculations are performed to extract the energy change ∆E2 due to the spatially modulated component, from which an approximate value of the normal fraction fn is obtained (see text): the first calculation with |Ω + p | 2 (0 + )/|Ωc| 2 = ǫgauge/( kc) and Ω - p ≡ 0, the second one with |Ω ± p | 2 (0 + )/|Ωc| 2 = ǫgauge/(2 kc). Empty squares with error bars: spatially homogeneous system of size L (with periodic boundary conditions), which corresponds to N ≃ 40000 atoms. Filled circles with error bars: system in a circular potential well (see text), with a total atom number N ≃ 30000 adjusted to have the same central density ρ, hence same degeneracy temperature T d and healing length ξ as in the homogeneous case. The suggested experimental values of Table I correspond to the circled point indicated by the oblique arrow and require the measurement of a relative energy change of ≃ 0.7%. Thin lines: quadratic fit (with no linear term, see text) of the effective fn as a function of ǫgauge, over the interval 0.05 ≤ ǫgauge/(mkBT d ) 1/2 < 0.2 (dashed lines for the spatially homogeneous system, solid lines for the system in the potential well). Black, red and blue lines and points (from bottom to top) correspond to an excitation sequence with γ/(csq) = 0.4, 0.2 and 0.1 respectively, with cs = (ρg/m) 1/2 the Bogoliubov sound velocity. The extrapolated zero-ǫgauge values of the effective normal fraction are indicated on the figure. The region between the horizontal dotted lines indicates the confidence interval of the thermodynamic prediction shown in Fig. 3. The residual 14% deviation of fn from the thermodynamic value (after linear extrapolation to γ = 0) is a finite qξ effect (see text).

with n real = 1000 realizations. As explained in the Appendix E, such a large number of realizations turned out to be necessary in the small ǫ gauge regime to compensate the fluctuations of the current j(t = 0) at the initial time. The resulting statistical error for the effective normal fraction is about 2% for the parameters of Table I.

In the numerical results of Fig. 4, this issue was circumvented by a simulation trick so that n real = 1000 is a sufficiently large number of realizations to keep error bars below 2% in the limit ǫ gauge → 0 [48]. It is also shown in the Appendix E that a measurement of the variance of the deposited energy (rather than its mean value) in principle also allows to access the normal fraction, if one knows experimentally the initial energy of the gas.

After linear extrapolation to γ = 0 from the lowest two values of γ in Fig. 4, which gives f n ≃ 0.176, the residual disagreement of f n with the thermodynamic result indicated by the dotted lines is about 14%, of the same order as the finite q correction predicted by Bogoliubov theory, see Fig. 2a. We have thus performed simulations for qξ = 1/2, for a spatially homogeneous system having the same parameters as in Fig. 4 and for a gauge field switch-off rate γ/(c s q) = 0.2 (not shown). It is found that the range of linearity of the effective normal fraction with ǫ 2 gauge is of the same order, up to ǫ gauge /(mk B T d ) 1/2 ≃ 0.27. Repeating the zero-ǫ gauge extrapolation procedure as in Fig. 4, we find for qξ = 1/2 an effective normal fraction of 0.164(3), rather than 0.146(2) for qξ = 1. Assuming a quadratic dependence of the effective normal fraction with q at fixed γ/(c s q), as predicted by Bogoliubov theory, see Fig. 2b, this leads to an additive correction to f n close to 0.024. We finally reach f n ≃ 0.20, which is within the statistical error bars of the thermodynamic value [50]. On a smaller system with M = 32 2 modes, we have reached a similar conclusion, taking even smaller values of qξ and γ/(c s q), see the blue circle in Fig. 3 and the figure caption. The deviation observed in Fig. 4 is therefore not a systematic error of the proposed method.

As a final check, we have tested the locality of the proposed measurement scheme by performing the simulation for two different geometries. Empty squares and dashed lines in Fig. 4 correspond to a spatially homogeneous system with periodic boundary conditions, while the filled disks and solid lines correspond to a system trapped in a circular well with steep walls of the form U (r) = ζk B T {tanh[(r -L/2)/(ξ/2)] + 1}, where the numerical parameter ζ determines the energy cut-off in the classical field model [START_REF]G/π -2 ln 2, G being the Catalan constant G ≃ 0[END_REF]. The probed region is at the center of the potential well, r 0 = 0. Even in the non-linear regime, where linear response theory fails, the effective normal fractions are almost the same in both geometries.

V. OPTICAL MEASUREMENT

The proposal that we have illustrated in the previous section was based on the measurement of atomic quantities, namely the deposited energy in the atomic cloud at the end of the gauge field sequence. The present section is devoted to the presentation and the characterization of an alternative, all-optical route to measure the normal fraction f n : information on the response of the atomic cloud to the gauge field can be retrieved from the trans-

(|Ω + p | 2 +|Ω - p | 2 ) t=0 + 2|Ωc | 2 0.09 gauge field switch-off rate γ = 0.2 csq with q = 1/ξ 1/γ = 4.8 ms q = 0.15kc
level schemes of Fig. 6: mitted probe beams once they have crossed the atomic cloud. Recent works have in fact pointed out that the strong frequency-dependence of the dielectric constant of an optically dressed medium in the electromagnetically induced transparency (EIT), already used experimentally to strongly reduce the light group velocity [START_REF] Fleischhauer | [END_REF][52][53], can be exploited for velocimetry experiments: Information on the current profile of an atomic cloud was predicted to be imprinted onto the phase of the transmitted probe beam [54,55].

first choice second choice level |a |F = 1, -1 |F = 2, -2 level |b |F = 2, -2 |F = 1, -1 level |e |F ′ = 2, -1 |F ′ = 2,
In the present case, the gas is illuminated by probe light of angular frequency ω p and coupling light of angular frequency ω c . As in Sec.IV the coupling light is a plane wave propagating along y axis. In the subsection V A, where it is proposed to measure the transmitted mean electric field amplitude, the probe light consists, as in Sec.IV, of two beams with wavevectors k ± p at a small angle with the z axis, as sketched in Fig. 1; differently from [54,55], the matter-wave current pattern is generated by the same laser beams that are then used for probing. The configuration is slightly different in the subsection V B, where it is proposed to measure the pattern of the scattered light intensity: in this case, the probe light consists of a single beam with a wavevector k p directed along the z axis, as sketched in Fig. 5. Even though no spatially modulated matter-wave current is generated by the laser beams, still one can extract information on f n from the zero-mean matter-wave currents induced by thermal fluctuations in the gas.

A. Extracting fn from the amplitude of transmitted light

The transmission and reflection of probe light from the two-dimensional atomic cloud can be described in terms of Maxwell equations. In particular, the dipole polarization of the atoms provides a source term for the probe electric field E p at angular frequency ω p : for the positive frequency parts, one has in the paraxial approximation with respect to the z axis:

(∆ + k 2 p ) E p = - k 2 p ǫ 0 P p , (40)
where k p = ω p /c and the Laplacian operator is threedimensional. In the considered Λ atomic configuration, see Fig. 1, one may think that the atoms occupy in each point of space the non-coupled state |N C(r) , leading to a vanishing mean atomic polarization P p . This is actually not exactly the case: Due to the atomic motion, the atomic internal state does not follow adiabatically and there exists a small coupling between the position dependent coupled state |C(r) and non-coupled state |N C(r) , the so-called motional coupling [30]. This leads to a small non adiabatic atomic polarization, that we now evaluate. Within a perturbative picture, we simply need to calculate the mean atomic polarization P p induced by the unperturbed laser fields. In terms of the three-dimensional atomic field operators, this reads

P p (r) = d ae Ψ † a (r) Ψe (r) . (41)
With the usual adiabatic elimination of state |e , as shown in the Appendix B, the atomic field in the ex-cited state can be written in terms of the atomic field operator χ3D in the coupled |C state, χ3D (r, t) = C(r, t)|a Ψa (r, t) + C(r, t)|b Ψb (r, t) [START_REF] Sinatra | [END_REF] as

Ψe ≃ Ω c 2(δ + iΓ/2) [1 + |Ω p /Ω c | 2 ] 1/2 χ3D + Γ 1/2 Be . (43)
From the explicit form of the noise term Be given in the Appendix B, it is immediate to see that it gives a zero contribution to the mean in Eq. [START_REF] Mandonnet | [END_REF]. Since the atoms are mostly in the uncoupled |N C state, we can approximate the atomic field Ψ a in the |a state by its |N C component. To lowest order in Ω p /Ω c we then have

P p (r) ≃ d ae Ω c 2(δ + iΓ/2) φ † 3D (r) χ3D (r) . ( 44
)
The next step is to perturbatively evaluate the field χ3D that is created in the coupled |C state by the motional coupling [30] of |C to |N C . The details of the procedure are given in the Appendix B. To first order in Ω p /Ω c and for the magic choice Eq.( 23) of ∆ c , one obtains after adiabatic elimination χ3D ≃ -

4 (δ + iΓ/2) m|Ω c | 2 ∇ φ3D • ∇ Ω p Ω c + Γ ′ 1/2 Bχ . ( 45
)
An explicit expression for the noise term Bχ is given in the Appendix B: again, the noise term has a zero expectation value and does not contribute to the optical polarization. The final form of the optical polarization in terms of the three-dimensional atomic density and current operators (17) and ( 15) reads, to zeroth order in q/k p :

P p (r) = - 4 |d ae | 2 |Ω c | 2 × × (k p -k c ) • j 3D (r) + 2im ∇ n 3D (r) E p (r), ( 46
)
where the Rabi frequency of the probe beam has been eliminated in favor of the electric field using the definition -d ea E p = Ω p /2 and the detuning δ has disappeared from the formula. The first term proportional to the atomic current operator has a simple semi-classical interpretation in terms of the reduced group velocity in the Electromagnetically Induced Transparency regime, as anticipated in [54,55]: in this regime, the refractive index strongly depends on the Raman detuning, which in turn depends on the atomic speed because of the Doppler effect. However, the expression (46) differs from the semiclassical one that was used in [55] in two ways. First, the current operator in Eq.(46) differs from the physical current of atoms by the gauge field

j phys = j 3D - 1 m n 3D A 3D . (47)
As the proposal in [55] addressed a pre-existing current profile and the weak probe beam induced a vanishingly small gauge field, the difference was irrelevant in that case. Here, on the contrary, the mean current is itself proportional to the gauge field so that the difference between the two operators really matters. Second, the expression (46) contains an extra term proportional to the average density gradient. In contrast to the first term, this one is purely imaginary. As a result, it only affects the intensity of the transmitted light via a combination of absorption and/or amplification effects. In particular, as it does not induce any phase shift on the light, it does not interfere with the proposal of [55].

In order to calculate the modification δE p induced by the atoms on the transmitted electric field of the probe, one has to insert the polarization (46) as a source term into the Maxwell equation (40). Within a standard approximation, we can neglect diffraction effects stemming from the in-plane part of the Laplace operator in (40) and integrate the z dependence across the atomic cloud. Taking into account the appropriate boundary conditions for δE p , this leads to the expression

δE p (x, y, z) ≃ e ikpz ik p 2ǫ 0 ∞ -∞ dz ′ e -ikpz ′ P p (x, y, z ′ ) (48)
for the transmitted field in the z > 0 region. In order for the approximation to be accurate, z has to be much larger than the thickness a z ho of the atomic pancake, but at the same time much smaller than the diffraction length k p /q 2 , where q is the characteristic wavevector of the inplane modulation of the atomic density and current.

Along z, the atomic field varies as the harmonic oscillator ground state wavefunction φ 0 (z), see Eq. (19). Performing the integral over z ′ , this gives the final expression for the variation of the transmitted field

δE p (r) = 2ik p |d ae | 2 ǫ 0 |Ω c | 2 k c • j(r) + 2im ∇ n(r) E p (r)
(49) in terms of the two-dimensional density n(r) and current j(r) operators. The first contribution proportional to the atomic current gives a phase shift, while the second contribution proportional to the atomic density gradient is responsible for absorption and amplification of the probe beam.

The atomic current profile created by the gauge field is evaluated using the linear response formulas (3) and (4) as discussed in detail in the previous sections. The gauge field is assumed to be switched on slowly enough as compared to the characteristic frequencies of all the excitation modes of the gas at wavevector q. Within the linear response regime, the contribution to the current due to the spatially modulated gauge field at q may be isolated by a suitable combination of measurements with different values of Ω ± p , which gives

j 2 (r) = ρk c mq|Ω c | 2 (q cos α + f n e z × q sin α) × Ω + p Ω - * p e iq•r + Ω + * p Ω - p e -iq•r e -2|r-r0| 2 /w 2 , ( 50
)
where r is now in the xy plane and α is the oriented angle that k c makes with q [36]. Inserting this expression into (49) and recalling the form (7) of the incident probe field, one can extract the phase shift experienced by the central part of the probe beams after crossing the atomic pancake,

∆φ ± 2 = 6π ρ m |Ω ∓ p | 2 |Ω c | 4 BΓ [cos 2 α + f n sin 2 α], (51)
where we have assumed k p ≃ k c . The total decay rate of the e state by spontaneous emission is indicated by Γ and B is the branching ratio for the decay to the a state, so that BΓ = |d ae | 2 k 3 c /(3π ǫ 0 ). As we have already mentioned above, the density gradient term in (49) only introduces an intensity modulation and is not responsible for any phase shift.

From a nonlinear optics point of view, the phase shift (51) can be interpreted as arising from a χ (3) optical nonlinearity of opto-mechanical origin similar to the one that was demonstrated in the experiment [56]: the nonlinear modulation of the optical response of the atoms is determined by the mechanical distortion of the cloud by the optical forces.

Inserting into (51) the values of Table I for the 87 Rb case, for the first choice of level scheme reported in the Appendix A, one has a branching ratio B = 1/4 and one finds a small, yet appreciable phase shift on the order of

∆φ choice 1 2 ≃ 6 • 10 -4 × [cos 2 α + f n sin 2 α]. (52)
For the second choice of level scheme in the Appendix A, the branching ratio is slightly larger, B = 1/3, but for the compromise choice Eq.(A25) |Ω c | 2 /Γ 2 is larger so that one finds a smaller phase shift

∆φ choice 2 2 ≃ 3 • 10 -4 × [cos 2 α + f n sin 2 α]. (53)
In addition to the phase shift of the transmitted beam that we have discussed so far, Bragg diffraction on the spatially modulated current profile produces a pair of additional beams of in-plane wavevector respectively ±3q/2 via a sort of four-wave mixing process. The relative intensity of these beams as compared to the incident probe beams is on the order of |∆φ 2 | 2 . For transverse gauge fields such that q • k c = 0 the contribution of the induced density gradient term of (49) vanishes by symmetry. In the case of longitudinal gauge fields, the relative correction is on the order of qξ.

B. Current fluctuations and the angular distribution of scattered light

All the calculations presented in the previous subsection aimed at evaluating the expectation value of the transmitted field amplitude. At this level of the theory, we were allowed to describe the probe beam as a coherent, classical field and we could neglect the fluctuations around the expectation value of both the light field amplitude and the atomic current and density operators. The formalism can be straightforwardly extended to quantum optical fields so to include the fluctuations of the atomic density and current. This is crucial when one aims at investigating the spontaneous scattering of light off the current fluctuations in the atomic gas. In this subsection, we shall in particular show how information on the normal fraction of the gas can be inferred from the angular distribution of scattered light. We shall make the approximation of replacing temporal derivatives of the electromagnetic field ∂ t E by -ick p E in Maxwell's equation. In particular, this misses retardation effects in the expression of the scattered fields in terms of the atomic dipoles, which is accurate since the system size is much smaller than c/Γ.

(k 2 p -Q 2 ) 1/2 ez + Q.
We consider the geometry sketched in Fig. 5: a single Gaussian probe beam is incident onto the atoms with a wavevector k p exactly orthogonal to the atomic plane, waist w centered at r 0 and a weak peak amplitude E 0 p . Inserting this form into (49) and taking the Fourier transform along the xy plane, one obtains the following operator equation for the scattered field component at in-plane wavevector

Q [57], E(Q, z) ≃ 2i k p |d ae | 2 ǫ 0 |Ω c | 2 π w 2 E 0 p d 2 q (2π) 2 e -q 2 w 2 /4 e -iq•r0 × k c • j Q-q + (Q -q) 2m n Q-q e ikz (Q) z . ( 54
)
Here, j Q and n Q are the spatial Fourier transforms of the two-dimensional current j(r) and density n(r) opera-tors; the Fourier transform of a product of two functions has been rewritten in terms of the convolution of their Fourier transforms. The z component of the propagation wavevector is determined by the photon dispersion as k z (Q) = (k 2 p -Q 2 ) 1/2 . The intensity of the scattered light at in plane wavevector Q is quantified by [57]

E † (Q)E(Q) = 2πk p |d ae | 2 w 2 |E 0 p | ǫ 0 |Ω c | 2 2 × d 2 q (2π) 2 d 2 q ′ (2π) 2 k c • j † Q-q + (Q -q) 2m n † Q-q × k c • j Q-q ′ + (Q -q ′ ) 2m n Q-q ′ × e i(q-q ′ )•r0 e -(q 2 +q ′ 2 ) w 2 /4 . ( 55
)
Since the system size is much larger than the waist w of the probe beam, we can for simplicity assume an effective translational symmetry along the xy plane. As a consequence, the correlation function that appears in ( 55) has a delta-distribution shape around equal wavevectors Qq = Qq ′ (see e.g. the next Eq.( 56)).

In contrast to the schemes proposed in the previous sections, where the duration 1/γ of the experiment had to be at least on the order of 1/(c s q), where c s = (ρg/m) 1/2 is the Bogoliubov sound velocity, the light scattering experiment discussed here can be performed on a much faster time scale, only limited by the characteristic rate Γ ′ of the internal atomic evolution time, Eq.(B7). As a result, the experiment can be performed in the small wavevector region Qξ ≪ 1 where the contribution to (55) of the terms involving the density fluctuations n Q is negligible [58]. Of course, efficient isolation of the scattered light from the incident beam requires that the scattering angle θ sc ≃ Q/k p be much larger than the diffraction cone of the probe beam, i.e. Q ≫ 1/w.

The instantaneous correlation function of the current in the y direction parallel to k c can be evaluated applying the fluctuation-dissipation relation (E5) to the current operator j q in an infinite space geometry. This gives

j † Q,y j Q ′ ,y = (2π) 2 δ 2 (Q -Q ′ )k B T Re[χ yy (Q ′ , ω = 0)] ≃ Q→0 (2π) 2 δ 2 (Q -Q ′ ) k B T ρ(r 0 ) m × [cos 2 φ sc + f n sin 2 φ sc ], (56)
where φ sc is now the azimuthal angle between k c and Q. Inserting this expression into (55) and taking the thermodynamic limit, one gets to the final expression for the scattered intensity in the momentum Q-space [57],

E † Q E Q = k p k c |d ae | 2 ǫ 0 |Ω c | 2 2 2πw 2 |E 0 p | 2 × k B T ρ(r 0 ) m [cos 2 φ sc + f n sin 2 φ sc ]. ( 57
)
To estimate the relative intensity of scattered light, it is useful to rewrite the expression (57) for the momentum space intensity E † Q E Q in terms of physically more transparent quantities such as the angular distribution I(θ sc , φ sc ) of scattered intensity. For small scattering angles |θ sc | ≪ 1, the infinitesimal solid angle and momentum space volume elements are related by dΩ = sin θ sc dθ sc dφ sc ≃ θ sc dθ sc dφ sc ≃ d 2 Q/k 2 p , so that

I sc (θ sc , φ sc ) ≃ k 2 p (2π) 2 E † Q E Q . ( 58
)
This immediately leads to the final expression for the angular distribution of the scattering intensity [57]

I sc (θ sc , φ sc ) I inc = k 2 p k c |d ae | 2 πǫ 0 |Ω c | 2 2 k B T ρ(r 0 ) m × [cos 2 φ sc + f n sin 2 φ sc ] (59)
in units of the incident intensity,

I inc = d 2 r |E p (r)| 2 = πw 2 2 |E 0 p | 2 . ( 60
)
From this expression, it is immediate to see that information on the normal fraction of the gas can be retrieved from the azimuthal dependence of the scattered intensity. In terms of the total spontaneous emission decay rate Γ of the e state and the |e → |a branching ratio B (we recall that BΓ = |d ae | 2 k 3 c /(3π ǫ 0 )), Eq.( 59) can be rewritten in the more transparent form [57],

I sc (θ sc , φ sc ) I inc = 9k 2 B T T d 2π 2 B 2 Γ 2 |Ω c | 4 [cos 2 φ sc + f n sin 2 φ sc ]. (61)
To estimate the feasibility of the proposed light scattering experiment, we now derive an upper bound on the number of useful scattered photons in a single shot of duration τ . Calculating the Poynting's vector of the probe beam, and using (60), we find an incident flux of probe photons Φ inc = (k p w) 2 |Ω 0 p | 2 /(12BΓ). Integrating in (59) the term proportional to f n over solid angles in the cone θ sc ≤ 1/(k p ξ), we obtain the flux Φ use sc of useful scattered photons. As a maximal duration, we take τ = 1/Γ non-ad fluo where the fluorescence rate of the atoms due to motional coupling between the non-coupled and the coupled states is given by Eq.(A2) (with 2|Ω + p | 2 replaced here with |Ω 0 p | 2 ). The number of single shot useful scattered photons is thus bounded by

N use ph ≤ 3πB 16 k B T ω z ρw 2 f n 1 (k c ξ) 2 . ( 62
)
Remarkably the Rabi frequencies Ω c and Ω 0 p have canceled out in the ratio of the scattered flux to the fluorescence rate. One recognizes in the right-hand side of (62) the effective mean number of atoms N eff at = πρw 2 f n /4 in the normal component illuminated by the probe beam, as in Eq.( 28). There is however a severe geometrical reduction factor, 1/(k c ξ) 2 , due to the small aperture of the useful scattering cone. For the parameters of Table I, with B = 1/3, and taking a waist w = 30ξ and f n = 0.2 as in Fig. 4, we find N eff at ≃ 1400, 1/(k c ξ) ≃ 0.15, which leads to N use ph ≤ 5. This remains accessible to current quantum optics experiments. For k B T / ω z fixed, the upper bound in Eq.( 62) scales as ρ 2 , since 1/ξ 2 scales as ρ, so that larger values of photon numbers for a given waist may be obtained by increasing the density ρ of the two-dimensional Bose gas.

VI. CONCLUSIONS AND PERSPECTIVES

In this paper we have proposed and validated two methods to measure the superfluid fraction of a quantum fluid of ultracold atoms. The idea is to apply an artificial gauge field to the atoms with spatial oscillations within a localized envelope and to detect, within the linear response regime, the matter-wave current pattern that is generated in the fluid. This can be done in a mechanical way by measuring the energy that is deposited in the fluid at the end of a gauge field sequence. This can also be done in an all-optical way by observing the phase shift experienced by the same laser beams that are used to generate the artificial gauge field or the angular pattern of scattered light. We have shown that, by a careful choice of the parameters and of the atomic level scheme for 87 Rb, two experimental obstacles, the spontaneous emission and (in the case of the mechanical method) undesired light-shifts, can be put to an acceptable level. The experimental challenge remains in the required high sensitivity of the measurements, that is a detection of small energy changes in the mechanical method, and a detection of small optical phase shifts or small photon numbers in the all-optical methods.

The interest of the proposed methods is twofold: they do not require that the gas reaches thermal equilibrium in presence of the gauge field, and furthermore they give the possibility of reconstructing in a local way the spatial profile of the superfluid fraction of a trapped gas, independently from the presence or the absence of a Bose-Einstein condensate. This last feature is attractive in the study of the Berezinskii-Kosterlitz-Thouless transition to a superfluid state in two-dimensional Bose gases and of the superfluidity properties of Bose gases in disordered environments. It would also be interesting to extend the method to the study of superfluidity in multi-component atomic fermionic gases, which may require identification of suitable level schemes.

Appendix A: Experimental issues

In this Appendix we review some issues that may hinder an experimental implementation of our proposal. Our attention will be concentrated onto the most relevant case of 87 Rb atoms considered in the experiment of [59,60] and, in a two-dimensional context, in [4]. The level structure of this atomic species is sketched in Fig. 6. A possibly important advantage of this atom in view of an experimental implementation of the present proposal is that the singlet and triplet scattering lengths of ground state atoms are equal within a few percent [61]. We thus expect that all scattering lengths between arbitrary F = 1 or F = 2 sublevels of the ground state have all almost the same values [62], which leads to g aa ≃ g ab and therefore to a suppression of the spatio-temporal variation of the effective interaction constant g 3D (r, t) defined in (18). In the deposited energy method proposed in Sec.IV, this is important to reduce the emission of phonons in the atomic gas by the temporal modulation of the interaction constant. In the optical detection scheme of Sec.V, this is also important to suppress the contribution of the interaction term to the field χ3D and then to the optical polarization. Other atomic species such as Yb [18] or metastable He [63, 64] and/or different laser beam configurations [19,20,26] are expected to be useful for other purposes, e.g. to suppress spontaneous emission and/or generate artificial gauge fields with different geometries [65].

Two possible choices for the three states |a , |b , |e forming the Λ system on the D1 line of 87 Rb are considered, as sketched in the two panels of Fig. 6. For each choice, we determine the undesired effects (spontaneous emission, light-shifts, Raman leaks) stemming from deviations from the perfect adiabatic following of the noncoupled state by the moving atoms and from optical transitions to other levels not included in the Λ system. An eye will also be kept on trying to maximize the |e to |a branching ratio so as to reinforce the optical signal of Sec.V. To minimize spontaneous emission within the low saturation regime, we shall allow for a small detuning δ of both the probe and coupling beam carrier frequencies from the |a → |e and |b → |e transitions, respectively. The Raman detuning of the two beams is taken in a way to always fulfill the magic Raman condition (23).

a. First choice

The coupling beam propagates along the y axis and is taken with a σ + polarization with respect to the y quantization axis. The probe beam is taken as linearly polarized along y. For the three atomic levels forming the Λ system, we take

F' = 2 F' = 1 F = 2 F = 1 ∆ hf,g ∆ hf,
|a ≡ |F = 1, m F = -1 , |b ≡ |F = 2, m F = -2 and |e ≡ |F ′ = 2, m F ′ = -1 [66].
This scheme of levels and lasers is illustrated in the left panel (a) of Fig. 6. To estimate the importance of the non-adiabatic coupling between the |N C and the |C states due to the atomic motion, we can evaluate the ratio ρ C /ρ N C of the two-dimensional densities in the two states close to the center of the laser spot. This is done using the explicit formula (45) for the field in the |C state [69]. Using the fact that for quasi-2D samples the gradient is mostly along the harmonically trapped z direction and assuming Ω + p = Ω - p , we obtain for k p ≃ k c :

ρ C ρ N C ≃ 4(δ + iΓ/2)k c Ω + p Ω 3 c 2 ω z m . (A1)
Accuracy of the adiabatic approximation requires that this ratio is much smaller than unity.

The finite population that is present in the |C state as a consequence of non-perfect adiabaticity is responsible for the spontaneous emission of photons at a single atom rate:

Γ non-ad fluo = Γ ′ ρ N C ρ C = 4Γ |Ω + p | 2 |Ω c | 4 k 2 c m ω z , ( A2
)
where the fluorescence rate Γ ′ of the coupled state |C is defined by Eq.( 11) and a spatial average has been performed. Remarkably, Γ non-ad fluo does not depend on the detuning δ.

Other contributions to the fluorescence rate come from non-resonant excitation processes. Dominating among these are the excitation of the |a state to the states |F ′ = 1 or 2, m F ′ = 0 by the coupling beam at a rate

Γ c fluo = Γ |Ω c | 2 4∆ 2 hf,g , (
Γ Λ ′ fluo = 2 3 Γ |Ω + p | 2 ∆ 2 hf,e . ( A4
)
In these expressions, we have taken into account the tabulated hyperfine dipole matrix elements of the various optical transitions, and we have introduced the hyperfine splittings ∆ hf,g and ∆ hf,e given for 87 Rb in the caption of Fig. 6. Limiting ourselves to the most relevant regime where |Ω + p /Ω c | 2 > 1/100, we see that in the present case of 87 Rb atoms, Γ Λ ′ fluo ≫ Γ c fluo . The total fluorescence rate can then be approximated as the sum of Γ Λ ′ fluo and Γ non-ad fluo . For a given value of the gauge field (proportional to |Ω + p /Ω c | 2 ), the total fluorescence rate is minimized to

Γ min fluo ≃ 4 √ 6 3 Ω + p Ω c 2 Γ ∆ hf,e k 2 c ω z m 1/2 (A5)
by a careful choice of the coupling beam Rabi frequency

Ω opt c 2 = ∆ hf,e 6 k 2 c ω z m = √ 3∆ hf,e 2π 1/2 k c ma 3D g. (A6)
Here, we have expressed ω z in terms of the threedimensional scattering length a 3D and the reduced twodimensional coupling constant as given by Table I and Eq. (20). Inserting the actual parameters of the 87 Rb atom, and taking g = 0.1, we obtain |Ω opt c | 2 /Γ 2 ≃ 0.21. It remains to adjust the detuning δ to be in the weak saturation regime,

s ≡ |Ω c | 2 /2 |δ + iΓ/2| 2 1 10 . ( A7
)
Let us take the same values of the gauge field sequence as in Fig. 4: γ/(c s q) = 0.2 and qξ = 1 and ǫ gauge ≃ 0.15(mk B T d ) 1/2 , c s = (ρg/m) 1/2 being the Bogoliubov sound velocity. This choice of ǫ gauge leads to

|Ω + p (0 + )/Ω c | 2 ≃ 0.03(ρλ 2 c ) 1/2
: as this quantity has to be much smaller than 1 in order for the gauge field description of Sec.III to be valid, it is safe to impose ρλ 2 c < 10. Integrating over the exponential switch-off ramp of Ω p and eliminating g in terms of ω z and a 3D , this gives for the total fluorescence probability per atom,

P fluo = √ 6 2 Γ ∆ hf,e k B T d ω z ρ 2 g 2 1/2 = √ 6 4 Γ ∆ hf,e 1 (ρa 2 3D ) 1/2 .
(A8) For 87 Rb with the choice ρλ 2 c = 9, we obtain the not very impressive result,

P fluo ≃ 0.22 . (A9)
For the sake of completeness, it is important to note that for this choice, k B T = 0.1k B T d remains smaller than ω z , so that the Bose gas retains a two-dimensional character. The existence of other atomic levels in addition to the ones strictly needed to create the gauge field is responsible not only for dissipative effects such as fluorescence, but also creates reactive effects such as a spatially and temporally-dependent light-shift of the non-coupled state |N C . Among the dominating processes, the parasitic Λ ′ scheme creates a modulated light-shift potential

U Λ ′ (r) = |Ω p (r)| 2 3∆ hf,e . (A10)
A shift of the same order of magnitude arises from the coupling of the |a state to the |F ′ = 2 or 1, m F ′ = 0 by the coupling beam. An estimate of the energy deposited in the system by the U Λ ′ term as compared to the one ∆E 2 due to the gauge field can be obtained with Bogoliubov theory: Using Eq.(D15) with U 0 = (Ω + p Ω - * p )(0 + )/(3∆ hf,e ), η = γ and Q = q, one gets for γ ≪ ǫ q ≪ ρg:

∆E U Λ ′ ∆E 2 = |Ω c | 4 9f n (k c c s ) 2 ∆ 2 hf,e . (A11)
For actual parameters, the energy change due to ∆U Λ ′ turns out to be non-negligible. For the optimal value of the Rabi frequency Ω opt c from Eq.(A6), the ratio is

∆E U Λ ′ ∆E 2 = 2 ω z 3f n ρg . (A12)
For 87 Rb, one finds the discouraging result ∆E U Λ ′ /∆E 2 ≃ 600g/(ρλ 2 c f n ), which remains much larger than unity even for ρλ 2 c = 9. Even if the deposited energy by the spurious potential U Λ ′ is much larger than the desired one of the gauge field, a suitable extrapolation procedure may take advantage of the different dependence on the laser intensities to isolate the effect of the gauge potential. An alternative possibility is to exploit the fact that the ∆E U Λ ′ contribution does not depend on the direction of q: Within the regime of linear response, this contribution can therefore be eliminated by taking the difference of the energy changes for respectively longitudinal and transverse gauge fields.

Another possible nuisance is the existence of stimulated Raman processes that may out couple the noncoupled C state to atomic ground state sublevels |c other than |a and |b , via excited state sublevels other than |e [70]. One may however check that for the proposed scheme these leaky Raman couplings are detuned from resonance by a frequency amount at least ∆ hf,g in absolute value and therefore harmless.

b. Second choice

Another possible choice for 87 Rb atoms is to take a σ -polarization for the coupling beam propagating along the y axis. The probe beam is again linearly polarized along y. The atomic levels forming the Λ system are now

|a ≡ |F = 2, m F = -2 , |b ≡ |F = 1, m F = -1 and |e ≡ |F ′ = 2, m F ′ = -2 .
The strong two-body losses that are generally experienced by the upper hyperfine manifold of the ground state are here suppressed by the choice of a maximal m F state for |a : as collisions between ultra-cold atoms mostly occur in the s-wave scattering channel, conservation of the sum of the m F 's then prevents transition to the lower hyperfine manifold.

The fluorescence rate per atom due to the motional coupling between |N C and |C is still given by Γ non-ad fluo as defined in Eq.(A2). As there is no longer any parasitic Λ ′ system, the fluorescence due to laser excitation of |a or |b to excited state sublevels other than |e is now dominated by the transitions |b → |F ′ = 1 or 2, m F ′ = -1 due to the probe beam. Thanks to the reduced occupation probability ≃ |Ω p | 2 /|Ω c | 2 of sublevel |b in the atomic state |N C and to the larger hyperfine splitting ∆ hf,g of the ground state, the fluorescence rate is strongly suppressed. After spatial averaging it amounts to

Γ p fluo = 3 2 Γ |Ω + p | 4 ∆ 2 hf,g |Ω c | 2 . (A13)
Other fluorescence processes on the D2 line (e.g. the transition |a → |J ′ = 3/2, F ′ = 3, m F ′ = -3 excited by the coupling beam) are several orders of magnitude weaker than Γ p fluo thanks to the huge fine structure splitting of 2π × 7 THz.

Since the two terms in the sum Γ non-ad fluo + Γ p fluo experience different switch-off functions e -γt and e -2γt , we integrate over time to calculate the total fluorescence probability:

P fluo = 4Γ γ |Ω + p (0 + )| 2 |Ω c | 4 k 2 c m ω z + 3Γ 4γ |Ω + p (0 + )| 4 |Ω c | 2 ∆ 2 hf,g . (A14)
After optimization over the coupling beam intensity, the minimal fluorescence probability

P min fluo = 2 √ 3 γ Γ ∆ hf,g Ω + p (0 + ) Ω c 3 k 2 c m ω z 1/2 (A15)
is obtained for a coupling beam Rabi frequency such that

|Ω opt c | 4 = 16 3 |Ω c | 2 |Ω + p | 2 ∆ 2 hf,g k 2 c m ω z . (A16)
Introducing the reduced quantities γ = γ/(ρg) and ǫgauge = ǫ gauge /(mk B T d ) 1/2 , and eliminating ω z in terms of g and a 3D , we finally obtain

P min fluo = 9 128π 1/4 ǫ3/2 gauge γ Γ ∆ hf,g λ c a 3D (ρλ 2 c ) -1/4 . (A17)
For the parameters of Table I, in particular ρλ 2 c = 9, one finds |Ω opt c | 2 /Γ 2 ≃ 5.5 so that a detuning |δ| > 5Γ is required to remain in the weak saturation regime. This resulting probability of spontaneous emission per atom in the deposited energy measurement is very small, P fluo ≃ 0.008.

(A18)

As compared to the first one, this second choice then provides a strong reduction of the spontaneous emission rate by a factor almost 30. Another advantage of this second choice is that lightshift effects are potentially smaller thanks to the absence of the parasitic Λ ′ scheme. The probe beam on the |b → |F ′ = 1 or 2, m F ′ = -1 transitions produces a lightshift which, after average in the |N C state, leads to the spurious potential

U p (r) = - |Ω p | 4 4|Ω c | 2 ∆ hf,g . (A19)
The amount of energy that is deposited in the gas by this spurious potential can be estimated using twice Eq.(D15), first with

U 0 = -|Ω + p (t = 0 + )| 4 /(|Ω c | 2 ∆ hf,g ), η = 2γ, Q = q,
and second with U 0 four times smaller, η = 2γ, Q = 2q. Neglecting 2 γ with respect to ǫ q , and taking qξ = 1, we obtain the following estimate for the spurious deposited energy

∆E U p = 133 160 N ρg 2 |Ω + p (0 + )| 8 |Ω c | 4 ∆ 2 hf,g . (A20)
For the coupling beam Rabi frequency (A16) minimizing spontaneous emission, the ratio of the energies deposited by the gauge field and the spurious potential amounts to:

∆E U p ∆E 2 = 133 √ 2 960π 3/2 g ǫgauge f n λ c ρ 1/2 a 2 3D .
(A21)

For g = 0.1 and ǫgauge = 0.15 and using the 87 Rb parameters summarized in Table I, one finds the still quite unfortunate result

∆E U p ∆E 2 ≃ 4 f n ≫ 1. (A22)
A possibility to overcome this difficulty and separate ∆E from ∆U p is to use the same strategy proposed to separate ∆E 2 from ∆E 1 by exploiting the different dependence of the two quantities on the ratio |Ω + p /Ω - p |(0 + ). This can be combined with the choice of a compromise value of |Ω c | 2 /Γ 2 that allows to strongly suppress the light-shift potential without introducing a too large spontaneous emission rate.

To this purpose, we fix ǫgauge = 0.15, γ = 0.2ρg/ , qξ = 1, g = 0.1 and we take as free parameters X = ρλ 2 c and Y = |Ω c | 2 /Γ 2 . Inserting the relevant parameters for 87 Rb as in Table I, we obtain

P fluo = 6.86 • 10 -2 X 1/2 Y + 7.46 • 10 -4 Y, (A23)
∆E U p ∆E 2 = 0.132Y 2 f n . (A24)
A reasonable compromise between the two competing effects is to choose X = 9 and Y = 0.5, which corresponds to

ρλ 2 c = 9 and |Ω c | 2 Γ 2 = 0.5 . ( A25
)
As a result, for the same parameters T /T d = 0.1 and f n = 0.2 used in the classical field simulations of section IV, we obtain the quite encouraging values

P fluo ≃ 0.045 (A26) ∆E U p ∆E 2 ≃ 0.16 . (A27)
We have checked that at the resulting temperature k B T / ω z = 0.62 the Bose gas retains a two-dimensional character and that the validity of the gauge field model of Sec.III is guaranteed by the resulting probe beam Rabi frequency

|Ω + p (0 + )| 2 /|Ω c | 2 = 0.03X 1/2 = 0.09 ≪ 1.
To conclude, we remark that the energy deposited by the spurious light-shift is suppressed in the extrapolation procedure proposed in subsection IV C to eliminate effects beyond the linear response regime to the gauge field, that is an extrapolation of the effective normal fraction that is linear in ǫ 2 gauge . Since ∆E U p scales as |Ω + p (0 + )| 8 with the probe beam Rabi frequency, it varies as ǫ 4 gauge with the amplitude of the gauge field, so it indeed contributes to the effective normal fraction as ǫ 2 gauge , which extrapolates to zero. This works of course for ǫ gauge small enough for the mechanical effect of the spurious lightshift to be treatable in the linear response regime.

Appendix B: Dum-Olshanii theory for many-body systems

In a seminal work [17], Dum and Olshanii have shown that an effective gauge field appears in the theoretical description of a three-level atom interacting with a laser field on a Λ transition. Here we use the formalism of the Quantum Stochastic Differential Equations (see e.g. §8.3.2 of [START_REF] Gardiner | Quantum Noise[END_REF]) to extend this idea to an interacting Bose gas in second quantized form.

We start with the master equation for the density operator σ of the many-body system, assuming for simplicity that spontaneous emission corresponds to a net loss of atoms,

d dt σ = 1 i [H, σ] + Γ d 3 r Ψe (r)σ Ψ † e (r) - 1 2 { Ψ † e (r) Ψe (r), σ} , (B1)
where [, ] ({, }) stands for the commutator (anticommutator). The Hamiltonian H is the sum of the single-particle kinetic and trapping terms, of the interaction terms (that we formally model as local Dirac-delta interactions), of the internal energy of the atomic excited state, and of the coupling terms of the atoms to the laser fields. Since the loss rate Γ greatly exceeds the kinetic, trapping and interaction energies, we can neglect the external dynamics of the excited state and write

H ≃ d 3 r α=a,b - 2 2m Ψ † α ∆ Ψα + U 3D Ψ † α Ψα + d 3 r g aa 2 Ψ †2 a Ψ2 a + g bb 2 Ψ †2 b Ψ2 b + g ab Ψ † a Ψ † b Ψb Ψa + d 3 r (-δ) Ψ † e Ψe + d 3 r Ω p 2 Ψ † e Ψa + Ω c 2
∂ t Ψe = 1 i [ Ψe , H] - 1 2 Γ Ψe + Γ 1/2 Fe (r, t) . (B3)
Here, the quantum noise term Fe is δ-correlated in position and time, e.g. [ Fe (r, t), F † e (r ′ , t ′ )] = δ(rr ′ )δ(t -t ′ ), and we recall that the expectation value of normally-ordered products of noise operators vanish, e.g.

F † e Fe = 0, since the bath does not provide an incoming flux of e atoms.

The only non-zero contributions to the commutator in Eq.(B3) originate from the excited state internal energy and from the atom-laser coupling term. This latter term can be expressed solely in terms of the atomic field operator χ3D in the coupled internal state |C , as defined in Eq. [START_REF] Sinatra | [END_REF]. Along the lines of [71], we formally integrate [∂ t + (-iδ + Γ/2)] Ψe = Ŝ neglecting a transient of duration 1/Γ as Ψe (r, t) = +∞ 0 dτ e -(-iδ+Γ/2)τ Ŝ(r, t -τ ).

(B4)

The Rabi frequencies Ω c,p and the atomic field χ3D have a negligible variation during 1/Γ and may be replaced by their values at time t in the integrand. This leads to Eq.( 43) of the main text, where the noise term is defined as Be (r, t) = +∞ 0 dτ e -(-iδ+Γ/2)τ Fe (r, t -τ ). As explained in section III, we are in a regime where the atoms are mostly in the non-coupled state and the field χ3D in the coupled state is small and a perturbation expansion in powers of χ3D can be performed. The gauge field formalism discussed in Sec.III for the evolution of the atomic field φ3D in the non-coupled state [defined in Eq.( 12)] is already recovered at zeroth order in χ3D . From this zeroth order approximation of φ3D , it is then easy to obtain the first order contribution to the field χ3D that is required in Sec.V to evaluate the optical polarization of the moving atoms.

From Eq.( 12) the equation of motion for φ3D is

∂ t φ3D = 1 i [ φ3D , H] + α=a,b Ψα ∂ t N C|α . (B5)
By the very definition of non-coupled state, the excited state internal energy and the atom-laser coupling terms give an exactly vanishing contribution to the commutator. In the other terms of the Hamiltonian as well as in the last sum in Eq.(B5), we can perform the approximation Ψα ≃ α|N C φ3D , which is accurate at zeroth order in χ3D . After an integration by part and noting that N C|∂ t |N C and N C|∇|N C are purely imaginary quantities, we find that up to this order φ3D follows a purely Hamiltonian evolution governed by Eq.( 13).

The equation of motion of χ3D has the form

∂ t χ3D = 1 i [ χ3D , H] + α=a,b Ψα ∂ t C|α . (B6)
The commutator with the internal excited state energy term introduces a Ψe term, that we replace with Eq.( 43): in this way, both a noise term and a complex, position dependent energy term -(iδ ′ + Γ ′ /2) χ3D appear in the equation. The real quantities δ ′ and Γ ′ are given by Eq.( 11) and correspond to light-shift and damping effects, respectively. Since Γ ′ is much larger than the kinetic, trapping, interaction and recoil energies of the atoms, we can neglect these latter terms in the evolution equation of the coupled state, and only keep the coupling to φ3D . This amounts to keeping in Eq.(B6) only the contributions to the kinetic, trapping and interaction terms of the Hamiltonian H that contain one single factor χ † 3D and an arbitrary number of φ3D and φ † 3D factors. In this way, we obtain

∂ t χ3D ≃ -(iδ ′ + Γ ′ 2 ) χ3D + φ3D C|[-∂ t + i 2m ∆]|N C + i m ∇ φ3D • C|∇|N C + 1 i G φ † 3D φ2 3D + Γ ′ 1/2
Fχ (B7)

where we have introduced a complex position and time dependent coupling constant

G = C|a a|N C [(g aa -g ab ) N C|a a|N C -a ↔ b] (
F † χ (r ′ , t ′ )] ≃ |δ+iΓ/2| 2 Γ e i(t-t ′ )δ-Γ|t-t ′ |/2 δ(r -r ′
). Since we are working in a low saturation regime in which Γ ≫ Γ ′ , the time dependent factor in front of δ(rr ′ ) may be replaced with a Dirac of t -t ′ , so that Fχ is in practice a spatio-temporally delta-correlated noise.

The last step is to expand Eq.(B7) to first order in Ω p /Ω c . We also limit ourselves to zeroth order in q/k p and in 1/(k p w), and we neglect the temporal derivative of the switch-off function f (t). Then (i) for g aa ≃ g ab , the fourth contribution in the right-hand side of Eq.(B7) vanishes [72], and (ii) for the magic choice Eq.( 23), the second contribution in the right-hand side of Eq.(B7) vanishes. With the same adiabatic elimination technique adopted for Ψe and taking into account the fact that δ ′ and Γ ′ vary very slowly on the scale of 1/Γ ′ , we are finally led to the final equation Eq.(45) with a noise term defined by Bχ (r, t) = +∞ 0 dτ e -[iδ ′ +Γ ′ /2](r,t)τ Fχ (r, t -τ ). We complete the discussion by giving the back-action of the field χ3D on the field φ3D , a back-action that was already considered for a specific single atom geometry in [30]. The linear coupling of φ3D to χ3D originates from terms in the Hamiltonian that are linear in χ3D , leading to

∂ t φ3D back = χ3D N C|[-∂ t + i 2m ∆]|C + i m ∇ χ3D • N C|∇|C + 1 i [G χ † 3D φ2 3D +2G * φ † 3D φ3D χ3D ]. (B9)
Expression of the back-action solely in terms of φ3D and noise operators is obtained by replacing χ3D in the resulting equations of motion with its adiabatic approximation derived from Eq.(B7). This leads in general to a lengthy formula. For simplicity, we give the result for g aa = g ab , to leading order in Ω p /Ω c , we also neglect the contribution to C|∂ t |N C of the time-derivative of the switch-off function f (t), and we use the specific form Ω p /Ω c considered in this paper, see Eqs. (6,7), restricting to zeroth order in q/k p and 1/(k p w) so that

∂ t φ3D back ≃ - i 2 m 2 4(δ + iΓ/2) |Ω c | 2 Ω p Ω c 2 × × [(k p -k c ) • ∇] 2 φ3D + noise terms. (B10)
After reduction to the xy plane, the deterministic term gives rise to two corrections to the evolution of φ3D : (i) a complex position dependent energy shift,

(δ ′′ -iΓ ′′ /2) = -2(δ + iΓ/2)ω z 2 k 2 c m |Ω p | 2 |Ω c | 4 , (B11)
and (ii) a complex correction to the mass along y, δm

y = 8 k 2 c |Ω p | 2 (δ + iΓ/2)/|Ω c | 4 .
The quantity δ ′′ is the lightshift potential experienced by the non-coupled bidimensional field. The spatial average of the fluorescence rate Γ ′′ of the non-coupled field coincides with the Γ non-ad fluo fluorescence rate previously discussed in (A2), as it should be. For the parameters of Table I, the reactive corrections δ ′′ and δm y are small provided that the detuning is not too large, |δ/Γ| < 5. For instance, an estimate for the undesired energy deposited by the δ ′′ potential can be obtained from the Bogoliubov theory for a homogeneous system, see Eq.(D15), leading to

∆E δ ′′ ∆E 2 ≈ 5 • 10 -5 f n (δ/Γ) 2 (B12)
where ∆E 2 is the desired deposited energy giving access to the normal fraction f n .

the exact relations:

∆E ≡ +∞ -∞ dt d H(t) dt (C2) = +∞ -∞ dt d 2 rA(r, t) • d dt j(r) (t), (C3)
where the second equality comes from a time-dependent Hellmann-Feynman theorem and a temporal integration by parts. Calculating j(r) (t) by linear response theory gives

∆E ≃ R dω 2π d 2 r d 2 r ′ ω Im   α,β χ ex αβ (r, r ′ ; ω)A α (r, ω) * A β (r ′ , ω)   , (C4)
where χ ex is the exact current susceptibility in real space taking into account the spatial inhomogeneity of the trapped cloud. Note that, contrarily to Eq.(C2), Eqs.(C3),(C4) still hold when f (t) = 0 for t < 0 and has a discontinuous jump in t = 0. We now use the particular form (C1) for A and consider the relevant limiting case qσ ≫ 1, q min(ξ, λ) ≪ 1, where ξ is the healing length of the gas and λ is the thermal de Broglie wavelength. We also assume that σ is much smaller than the radius of the trapped cloud, so that the density variation within a region of radius σ around r 0 may be neglected.

Within a local density approximation, we then replace χ ex with the susceptibility χ of a spatially homogeneous system with a density equal to the one of the trapped gas at position r 0 and with the same temperature, χ ex (r, r ′ ; ω) ≃ χ(rr ′ ; ω).

This local density approximation leads to

∆E ≃ R dω 2π ω |f (ω)| 2 d 2 R d 2 u |c + | 2 + |c -| 2 + c * + c -e -iq•(R+u/2) + c + c * -e iq•(R+u/2) |c + | 2 + |c -| 2 + c * + c -e -iq•(R-u/2) + c + c * -e iq•(R-u/2) Im[χ yy (u; ω)] e -|R-r0| 2 /σ 2 e -u 2 /(4σ 2 ) (C5)
where we have performed the change of variables r = R + u/2, r ′ = Ru/2. As we work in the qσ ≫ 1 regime, we have for example

d 2 R e -2iq•R e -|R-r0| 2 /σ 2 = e -2iq•R0 πσ 2 e -q 2 σ 2 ≪ πσ 2 (C6)
so that all the oscillating terms in R may be neglected.

Introducing the Fourier transform of χ yy (k; ω), which is an even function of k due to parity or rotational invariance, we obtain

∆E = ∆E 1 + ∆E 2 (C7) ∆E 1 ≃ +∞ -∞ dω 2π ω|f (ω)| 2 (2πσ 2 ) 2 d 2 k (2π) 2 Im χ yy (k; ω) |c + | 2 + |c -| 2 2 e -k 2 σ 2 (C8) ∆E 2 ≃ +∞ -∞ dω 2π ω|f (ω)| 2 (2πσ 2 ) 2 d 2 k (2π) 2 Im χ yy (q + k; ω)2|c + | 2 |c -| 2 e -k 2 σ 2 . (C9)
The second contribution ∆E 2 comes from the spatially modulated gauge field at q, while the first contribution ∆E 1 is due to the non-modulated term which follows the broad Gaussian envelope. The expression (24) in the main text is obtained from (C9) by noting that the integration over k is effectively limited by the Gaussian factor to a small region of radius 1/σ in which one is allowed to neglect the k-dependence of the susceptibility. Naively, one could guess that a necessary condition for the accuracy of our local density approximation is that the switch-off time of the gauge field γ -1 is short as compared to the characteristic time 2R/v for the induced mechanical perturbation to cross the whole cloud, to be reflected by its boundaries and to turn back to the excitation zone where it can interfere with the excitation process, v being the fastest between the sound and thermal speeds in the cloud of radius R.

This condition is actually sufficient, but not necessary within linear response theory. We now show for f (t) = Θ(t)e -γt , as in Eq.( 27), that the γ → 0 limit for the deposited energy scheme exists and coincides with the perturbation induced by the gauge field in the thermodynamic equilibrium state. As one can show by inserting the explicit form of the temporal Fourier transform of the gauge field into (C4) and performing the integral over ω, the deposited energy can be written in the form

∆E ≃ 1 2 d 2 r d 2 r ′ α,β A α (r, t = 0 + ) A β (r ′ , t = 0 + ) Re χ ex αβ (r, r ′ ; ω = iγ) (C10)
where we have introduced the Kubo formula for the exact current-current susceptibility

χ ex αβ (r, r ′ ; ω) = λ,λ ′ (π λ -π λ ′ ) λ|j α (r)|λ ′ λ ′ |j β (r ′ )|λ E λ ′ -E λ -ω -i0 + (
C11) in terms of the thermal equilibrium population π λ of quantum state λ. This quantity can be simply related to the susceptibility at thermodynamic equilibrium,

χ th αβ (r, r ′ ) = lim γ→0 Re[χ ex αβ (r, r ′ ; ω = iγ)] + + 1 k B T λ,λ ′ ;E λ =E λ ′ π λ λ|j α (r)|λ ′ λ ′ |j β (r ′ )|λ . (C12)
We recall that the thermodynamic susceptibility relates the mean current in a thermal equilibrium state at temperature T to the applied (weak) static gauge field via j = χ th * A, where * is the spatial convolution product. In the present case of an interacting gas, one can safely expect that the second line in (C12) gives a negligible contribution as there is no systematic degeneracy and the current operator j has no diagonal matrix elements since the eigenstate wavefunctions may be taken real. Within linear response theory, we thus obtain

∆E → γ→0 1 2 d 2 rd 2 r ′ α,β χ th αβ (r, r ′ )A α (r, 0 + )A β (r ′ , 0 + ).
(C13) Since local density approximation in the thermodynamic equilibrium state is a standard procedure, we expect that it can be applied to evaluate the right-hand side of Eq.(C13) [73]. As a consequence, within linear response theory, our deposited energy method (that measures the left-hand side of Eq.(C13)) should be accurately described by the local density approximation down to the γ → 0 limit.

Appendix D: Some results of linear response theory and the Bogoliubov expression of the deposited energy

A system of time independent Hamiltonian H 0 experiences, at times t > 0, a time dependent weak perturbation of Hamiltonian -ǫf (t)V, where ǫ → 0, the dimensionless time dependent factor f (t) is zero for t < 0 and tends rapidly to zero for t → +∞, and the operator V is time independent. At time t = +∞, the system is free again, with a mean energy modified by the perturbation. The question is to calculate the mean energy change to second order in ǫ.

Suppose first that, at t = 0 -, the system is prepared in the eigenstate |λ of H 0 of eigenenergy E λ . The energy change δE between time 0 and time +∞ is

δE = lim t→+∞ ψ(t)|(H 0 -E λ )|ψ(t) (D1)
where |ψ(t) is the system state vector at time t. The usual time dependent perturbation theory gives an expansion in powers of ǫ:

|ψ(t) = |ψ 0 (t) + ǫ|ψ 1 (t) + ǫ 2 |ψ 2 (t) + . . . ( D2
)
where |ψ 0 (t) = exp(-iE λ t/ )|λ ,

|ψ 1 (t) = - t 0 dτ i e -iH0(t-τ )/ f (τ )Ve -iE λ τ / |λ , (D3)
and the expression of higher order contributions is not needed. Using (H 0 -E λ )|ψ 0 (t) = 0 and ψ 0 (t)|(H 0 -E λ ) = 0, we find to second order in ǫ:

δE ≃ lim t→+∞ ǫ 2 ψ 1 (t)|(H 0 -E λ )|ψ 1 (t) . (D4)
In this paper, f (t) = Θ(t) exp(-γt), with γ > 0 and Θ(t) is the Heaviside step function, see Eq.( 27). Also, the system is prepared initially in a statistical mixture of eigenstates of H 0 with a probability distribution π λ . After explicit integration of (D3) over τ and then average over |λ , the expression for the signal to be detected experimentally is

Signal (V) ≡ lim ǫ→0 δE ǫ 2 = 1 2 Re λ,λ ′ (π λ -π λ ′ ) E λ ′ -E λ -i γ | λ ′ |V|λ | 2 . (D5)
The sum may be restricted to E λ = E λ ′ since the contributions with E λ = E λ ′ are zero. This also shows that the signal has a finite limit for γ → 0 + . Note that in a thermal equilibrium state π λ = Z -1 exp(-E λ /k B T ), the signal is necessarily positive.

The calculation of the noise on the experimental signal can be performed along the same lines. One defines δE 2 ≡ lim t→+∞ ψ(t)|(H 0 -E λ ) 2 |ψ(t) with the initial state vector |ψ(0) = |λ , and one finds after average over the initial state:

[Noise (V)] 2 ≡ lim ǫ→0 δE 2 ǫ 2 = γ→0 λ,λ ′ ,E λ =E ′ λ π λ | λ|V|λ ′ | 2 ≃ λ π λ [ λ|V 2 |λ -λ|V|λ 2 ], (D6)
where the approximate equality is based on the assumption that there are no systematic degeneracies in the many-body spectrum. The ǫ 2 scaling of the variance in Eq.(D6) shows that the typical value of the energy change at the end of the excitation sequence is of order ǫ. This scaling is to be contrasted with the ǫ 2 one of the expectation value that is suggested by Eq.(D5). We now apply the general formula Eq.(D5) to the Bogoliubov analysis of subsection IV B. In this case, ǫ = ǫ gauge /2 and V = V q + V -q with

V q = [0,L] 2 d 2 r e iqx j y (r) = k k y m a † k+q a k (D7)
where a k is the annihilation operator of a particle of the gas of wavevector k, and we have set q = q e x . In a translationally invariant system, the eigenstates |λ can be taken of well defined total momentum; as the action of V ±q changes this total momentum by ± q, the two operators V q and V -q cannot interfere in the signal and thus Signal (V) = 2 Signal (V q ). In terms of the annihilation operators b k of Bogoliubov quasiparticles, we can split

V q = V (0) q + V (2) q + V (-2) q (D8)
in terms of

V (0) q = k =0,-q k y m (U k U k+q -V k V k+q ) b † k+q b k (D9) V (2) q = ′ k =0,-q k y m (V k U k+q -U k V k+q ) b † k+q b † -k (D10) V (-2) q = ′ k =0,-q k y m (U k V k+q -V k U k+q ) b k b -(k+q) .(D11)
The primed sum ′ indicates restriction of the sum over wavevectors such that k y > 0. In Bogoliubov theory, the eigenstates |λ may be taken in the form of Fock states of quasiparticles. Since V (n) q changes the total number of quasiparticles by the amount n, the terms in the right hand side of Eq.(D8) cannot interfere in the signal and Signal (V q ) = Signal (V (0) q )+Signal (V (2) q )+Signal (V (-2) q

). (D12) Thanks to the clever writing of V (n) q with the constraint k y > 0 [74], there are no interferences in Eqs.(D9,D10,D11) between the different terms of the sums over k. As a result, the whole signal is the sum over the contribution of the different k's. The last trick is to express the ratios π λ ′ /π λ in terms of the mean occupation numbers n k of the Bogoliubov modes of energy ǫ k and to make use of the identity

e βǫ k = n k + 1 n k (D13)
satisfied by the Bose law. A little bit of rewriting taking advantage of the relation Signal(V) = (N/m)f eff n and of the remarks in [74], finally leads to Eq.( 31).

An alternative procedure is to calculate the currentcurrent susceptibility (C11) for a spatially homogeneous system within the Bogoliubov theory, which for q ⊥ e y gives in dimension d:

χ yy (q; ω) = 1 L d k =0,-q 2 k 2 y m 2 n k -n k+q ǫ k+q -ǫ k -ω -i0 + U 2 k+q U 2 k -U k+q V k+q U k V k - 1 + n k + n k+q -ǫ k+q -ǫ k -ω -i0 + V 2 k+q U 2 k -U k+q V k+q U k V k + 1 + n k + n k+q ǫ k+q + ǫ k -ω -i0 + U 2 k+q V 2 k -U k+q V k+q U k V k + n k+q -n k ǫ k -ǫ k+q -ω -i0 + V 2 k+q V 2 k -U k+q V k+q U k V k . (D14)
From Eq.(C4) one then recovers expression Eq.(31) for the effective normal fraction f eff n . Another application of Eq.(D5) is to calculate the energy deposited by the external potential U(r, t) = (U 0 e iQ•r +c.c)Θ(t)e -ηt in the spatially homogeneous case, for η > 0. This is useful for Appendix A and Appendix B to estimate the effect of undesired light-shifts. In second quantized form, and to leading order in Bogoliubov theory, one then has ǫV = N 1/2 (U Q +

V Q )[U 0 (b † Q + b -Q ) + h.c.
]. These terms do not interfere in Eq.(D5). For nonzero temperature, using (D13), we then obtain a temperature independent result

∆E U ≃ 2N |U 0 | 2 (U Q + V Q ) 2 Re 1 ǫ Q -i η . (D15)
Remarkably, this also allows to calculate the energy change ∆E g due to the switch-on-and-off of a spatially modulated coupling constant, δg(r, t) = (δg 0 e iq•r + c.c.)Θ(t)e -γt . For a spatially homogeneous system, to leading order of Bogoliubov theory, the relevant operator is ǫV = N 1/2 (U q + V q )[ρδg 0 (b † q + b -q ) + h.c.], so that one can formally apply Eq.(D15) with U 0 = ρδg 0 . This can be applied to the variation of the coupling constant due to g aa = g ab in Eq. (18). In this case δg 0 = 2g[(g ab -g aa )/g aa ](Ω In the section IV C, while presenting the numerical results on the deposited energy measurement, we mentioned that the statistical noise on the deposited energy was larger for smaller values of the gauge field amplitude ǫ gauge .

To understand this feature, it can be useful to rewrite the deposited energy for a single realization of the classical field simulation in the form δE = -d 2 r j(r, 0) • A(r, 0 + ) -+∞ 0 + dt d 2 r j(r, t) ∂ t A(r, t). (E1)

The first term comes from the abrupt switch-on of the gauge field. For each realization, it is of order ǫ gauge but averages to zero in the limit of an infinite number of realizations of the experiment as j(r, 0) = 0. In any actual calculation, an average over a finite number n real of realizations is taken, which gives a non-zero random value for δE scaling as ǫ gauge / √ n real . The relevant signal δE is given by the second term in (E1), obtained from the classical Hamiltonian identity dH/dt = ∂ t H. For small values of the gauge field switch-off rate γ, this term is of order O(ǫ 2 gauge ) as in this limit j adiabatically follows the thermal equilibrium value for the instantaneous value of the gauge field. As a result, the number of realizations that are needed to extract the signal out of the statistical noise due to the first term grows as |ǫ gauge | -2 , which perfectly explains the numerical observation.

Furthermore, it has been demonstrated by a number of recent cold atom experiments that noise is not always just an hindrance but can be also a source of useful physical information [75,76]. As a simple example, we consider here the amplitude of the noise on the energy that is deposited in the system at each realization of the experiment. This quantity is quantified by the average δE 2 of the square of the deposited energy in the γ → 0 limit. Looking at (E1), it is immediate to see that in the small ǫ gauge limit the dominant contribution comes from the square of the first term, which suggests that the noise on the deposited energy is related to the variance of the instantaneous fluctuations of the current operator. From the fluctuation-dissipation theorem, this quantity can then be related to the normal fraction of the gas.

This idea can be put on solid grounds by developing a full quantum calculation. The linear response theory calculation performed along the lines of Appendix D leads to the expression in Heisenberg picture

lim γ→0 [H 0 (+∞) -H 0 (0)] 2 ≃ [ d 2 r j(r) • A(r, 0 + )] 2 (E2)
where H 0 is the unperturbed Hamiltonian, that is the Hamiltonian without coupling to the gauge field. This relation connects the variance of the quantum equivalent of the deposited energy to the instantaneous fluctuations of the current operator and confirms our expectation based on the classical field model. Moreover, combined with the fluctuation-dissipation theorem, it can be the starting point for another proposal to measure f n .

For a generic hermitian operator V with vanishing diagonal matrix elements in the eigenbasis of H 0 , the fluctuation-dissipation theorem of linear response theory relates the imaginary part of the susceptibility χ to the Fourier transform of the correlation function S V V (t) = V (t) V (0) ,

Im[χ V V (ω)] = 1 2 S V V (ω) 1 -e -ω/kBT . ( E3
)
The Fourier transform S V V (ω) of the correlation function in the thermodynamic equilibrium state is defined as usual as

S V V (ω) = ∞ -∞ dt e iωt V (t) V (0) . (E4)
Under the assumption that most of the spectral weight of the V operator lies in the low-energy region ω ≪ k B T , we can approximate 1 -e -ω/kBT ≃ ω/k B T . This is a quite standard approximation of many-body theory and is generally accurate in the small q limit [START_REF] Pitaevskii | Bose-Einstein Condensation[END_REF]. After a few manipulations, it leads to the general expression

S V V (t = 0) = dω 2π S V V (ω) ≃ dω 2π 2k B T ω Im[χ V V (ω)] = k B T Re[χ V V (ω = 0)], (E5)
where the equivalent of Eq.( 26) was used to obtain the last identity. An application of the fluctuationdissipation relation (E5) to the susceptibility and the fluctuations of the mass current in liquid He can be found in [START_REF] Dalfovo | [END_REF].

The link between the variance of the deposited energy and the normal fraction is immediately obtained by applying (E5) to the specific operator

V = -d 2 r j(r) • A(r, 0 + ) (E6)
and isolating the contribution of the spatially modulated gauge field proportional to Ω + p Ω - p . In this way, using Eq.(C10), one is led to the final expression It is however crucial to note that a measurement of f n based on the relation (E7) requires taking expectation values of the Hamiltonian operator at different times. This may be experimentally challenging as it requires either a non-destructive measurement of the initial energy of the system at t = 0 before switching on the gauge field, or a very precise a priori knowledge of its value in a sort of microcanonical ensemble [77].

FIG. 1 :

 1 FIG.1: Scheme of the set-up under consideration. Left panel: generic sketch of the Λ configuration of atomic levels and laser beams involved in the optical processes. Center panel: view from above of the two-dimensional atomic pancake (lying within the dashed circle) and of the laser beams (the hatched disk is the spot of the probe beams in the xy plane). Right panel: side view. kc and k ± p are the wavevectors of the coupling and probe beams, respectively. In practice, k ± p = (k 2 p -q 2 /4) 1/2 ez ± (q/2) ex with q ≪ kp, where ei is the unit vector along axis i.

2D density ρλ 2 c = 9 ρ

 9 = 14µm -2 degeneracy temperature kBT d ≡ 2π 2 ρ/m T d = 500 nK temperature T = 0.1 T d T = 50 nK 2D interaction constant g = mg/ 2 g = 0.1 transverse confinement ωz = 0.23 2 k 2 c /m = 0.16 kBT d ωz 2π = 1.65 kHz healing length ξ ≡ (ρg) -1/2 ξ = 0.84µm reduced gauge field amplitude ǫgauge ≡ ǫgauge (mk B T d ) 1/2 ǫgauge = 0.15 probe beam Rabi frequencies

  ∆EU /∆E2 (for fn = 0.2) 33 0.16 TABLE I: Suggested values of the physical parameters for an experimental measurement of the normal fraction of a twodimensional Bose gas of 87 Rb using an artificial gauge field produced by laser (coupling and probe beam) excitation on the D1 line with an optical wavelength λc = 795 nm. The first block characterizes the thermal equilibrium of the gas. The second block determines the gauge field. The third block deals with the issues of spontaneous emission and spurious lightshift, for two atomic level schemes: the "first choice" is Fig.6a, the "second choice" is Fig.6b. For this second choice, the indicated value |Ωc| 2 /Γ 2 = 0.5 corresponds to a compromise between minimization of the fluorescence and of the spurious light-shift, the truly minimal P fluo being < 0.01, see (A18). Note that P fluo and ∆EU /∆E2 are basically unchanged if one takes qξ = 1/2 and γ/(csq) = 0.4. The three-dimensional scattering length a3D ≃ 100 Bohr radii is related to the twodimensional coupling constant by Eq.(20). A useful relation is mΓ/( k 2 c ) ≃ 792. cs = (ρg/m) 1/2 is the Bogoliubov sound velocity.

FIG. 5 :

 5 FIG.5: Scheme of the scattering geometry under examination. Probe light is incident at wavevector kp = kpez and the scattered light is collected at a wavevector ksc = (k 2 p -Q 2 ) 1/2 ez + Q.

2 FIG. 6 :

 26 FIG.6: Scheme of the internal levels of 87 Rb atom involved in the D1 transition (J = 1/2 → J ′ = 1/2). The hyperfine splitting of the ground (excited) state is ∆ hf,g = 2π × 6.834 GHz (∆ hf,e = 2π × 814 MHz). The natural linewidth of the excited state is Γ = 2π × 5.75 MHz. Thin (thick) lines indicate transitions that are induced by the probe (coupling) beams. Solid (dashed) lines indicate the desired (main undesired) transitions. The relevant dipole matrix elements are shown in units of the reduced dipole element J = 1/2||er||J ′ = 1/2 of the D1 line, in boldface for the desired transitions and in italic for the undesired ones. Our first choice for the Λ system |a , |e , |b states is shown in the left (a) panel, with|a = |F = 1, mF = -1 , |e = |F ′ = 2, m F ′ = -1 , |b = |F = 2, mF = -2 . The second choice is shown in the right (b) panel, with |a = |F = 2, mF = -2 , |e = |F ′ = 2, m F ′ = -2 , |b = |F = 1, mF = -1. Note that we have taken the y axis as the quantization axis of angular momenta.

  A3) and the excitation of the non-coupled state |N C to |e ′ = |F ′ = 1, m F ′ = -1 by the total probe plus coupling field with an effective Rabi frequency -2Ω p / √ 3, which results in the fluorescence rate on the parasitic Λ ′ configuration |a → |e ′ → |b :

  Ψ † e Ψb + h.c. . (B2) As previously defined, δ is the common value of the detuning of the probe and coupling beams from the |a → |e and |b → |e transitions. In a Heisenberg picture for the open atomic system, the ground state atomic field operators Ψα=a,b satisfy the usual evolution equations i ∂ t Ψα = [ Ψα , H]. On the other hand, conservation of the canonical commutation relations of the fields and of the Hermitian conjugation relation between Ψα and Ψ † α requires including a quantum Langevin term Fe in the evolution equation for the excited state field Ψe ,

  B8) with the convention g ba = g ab . The noise term is defined by Fχ = -i|δ + iΓ/2|(Ω * c /|Ω c |) Be . Its correlation properties are determined by the commutation relation [ Fχ (r, t),

  Appendix E: The noise on the deposited energy

  (+∞) -H 0 (0)] 2 , (E7) which demonstrates an alternative way of extracting the value of the normal fraction f n from a measurement of the statistical variance of the deposited energy in a series of experiments.

  + p Ω - * p )(0 + )/|Ω c | 2 so that for qξ = 1 and γ → 0, For the values of Table I and |g ab -g aa | 0.1|g aa | as expected for 87 Rb, this gives ∆E g /∆E 2 ≈ 7 • 10 -4 /f n which is negligible.

	∆E g ∆E 2	≃	16 5f n	mρg ( k c ) 2	g ab -g aa g aa	2	.	(D16)

Acknowledgments

We acknowledge useful discussions with G. La Rocca, S. Stringari, F. Gerbier, C. Salomon and J. Dalibard. The work of Y.C. was done as part of the ERC Project FERLODIM N.228177. Y.C. is a member of IFRAF and acknowledges financial support from IFRAF. I.C. acknowledges financial support from ERC through the QGBE grant.

Appendix C: Derivation of the expression for the deposited energy

We start from a two-dimensional system at thermal equilibrium with no average current and we apply a gauge field of the form A(r, t) = f (t)e y |c + e iq•r/2 + c -e -iq•r/2 | 2 e -(r-r0) 2 /(2σ 2 ) , (C1) where e α is the unit vector along direction α. The derivable dimensionless envelope function f (t) is assumed to be zero for t < 0 and to rapidly tend to zero for t → +∞. The time-independent coefficients c ± have the dimension of the square root of a momentum.

We are interested in evaluating the energy change of the system from t = 0 -to t = +∞ at the lowest order in c ± . We work in Schrödinger picture and we first use L 0 dy h(x, y, t) = 0, ∀x, ∀t. This requires a bit more care as the evolution operator is diagonal neither in the realspace nor in the k-space, nor in any mixed representation (real-space along some direction, Fourier-space along the other). In a first stage, we performed a local gauge transform on Ψ of the form U Ψ = e -iH(x,y,t)/ Ψ with ∂yH(x, y, t) = h(x, y, t). In this way, the operator that one has to exponentiate to obtain the evolution operator is simply proportional to p 2 y . Since h(x, y, t) has a zero spatial mean along y axis, H(x, 0, t) = H(x, L, t), ∀x, ∀t, so that the gauge transform is compatible with the periodic boundary conditions. In a second stage, we realized however that this gauge transform trick is not exact for our lattice model, since py is not simply ( /i)∂y (contrarily to continuous space). As a consequence, relation (E1) was not well satisfied numerically. Since the non-zero-spatial-mean version h(x, y, t) is factorized in a function f (x, t) of x and t times exp[-y 2 /(2σ 2 )], with w = 2σ, we simply diagonalized numerically once the operator W = py exp[-y 2 /(2σ 2 )]+h.c. (in practice, a 64 × 64 hermitian matrix for Fig. 4) and evaluated the evolution operator exp[idtf (x, t)W ] in the corresponding eigenbasis (where it is diagonal).

[44] A. J. Leggett, Rev. Mod. Phys. 71, S318 (1999).

[45] I. Carusotto [47] Above the critical temperature, there is no superfluid phononic branch, and one may expect that the condition qξ 1 for f eff n to reach the zero-q limit is replaced by the condition that qλ 1, where λ is the thermal de Broglie wavelength, which is less stringent in the high temperature limit. To test this intuitive idea, Bogoliubov theory is not applicable, but it is tempting to use the ideal gas model as a first guess. Then for qλ 1 and T ≪ T d , one obtains for the ideal quantum gas in the thermodynamic limit

and the chemical potential µ ≃ -kBT exp(-T d /T ). The condition that this is much smaller than fn, which is here equal to unity, thus does not only involve a condition on qλ, it also involves a comparison of 2 q 2 /2m to the chemical potential, which is much more stringent for T ≪ T d . This could be expected from the fact that, for the occupation numbers of the low momenta modes, the Bose-Einstein formula reduces to the Lorentzian shape

), which has a narrow momentum width ∝ |µ| 1/2 . As an example, to have f eff n = 0.9 for the ideal gas in the zero γ limit, one has to take qλ ≃ 0.48 for T /T d = 0.18, and qλ ≃ 0.08 for T /T d = 0.1 (these values of q correspond to qξ = 0.64 and qξ = 0.08 for the interacting gas with g = 0.1 2 /m, respectively).

[48] To avoid the rapid increase of the error bars on fn for ǫgauge → 0, due to the effect discussed in Appendix E, we performed, for each stochastic field realisation of the thermal ensemble, two independent temporal evolutions, one in presence of the gauge field A and the other in presence of the opposite gauge field -A. We then took the average of the two corresponding deposited energies. This in particular exactly eliminates all the noisy terms in the deposited energy that are linear in the gauge field.

[49] Measuring the gas energy is in general reduced to measuring the gas spatial density profile, before time of flight (to measure the trapping energy) and after time of flight (to measure the initial kinetic plus interaction energy, which is the expansion energy). These density measurements have reached a precision giving the absolute equation of state of the gas, the pressure for example, within five percent [21,22] or even at the percent level [23],

but not yet at the sub-percent level. Our scheme however requires a measurement of an energy change, not an absolute measurement of the energy, giving hope for a better precision.

Recent experiments have actually measured the change of the momentum [24] and of the energy [25] of a gas to determine its structure factor with good precision, not for a gauge field excitation as proposed here, but for a Bragg laser excitation.

[50] We have checked by performing simulations (for qξ = 1, γ/(csq) = 0.4 and ǫgauge/(mkBT d ) 1/2 = 0.02) with a duration τ = 6/γ that the truncation of the time evolution at τ = 3/γ as done in Fig. 4 introduces an error on the normal fraction which is within the statistical error bars.

[56] F. Brennecke, S. Ritter, T. Donner, T. Esslinger, Science 322, 235 (2008).

[57] Eqs.(54,55,57,59,61) have been obtained from the operatorial form of (49), where the atomic field operators have already been reduced to 2D. This is an approximation. The exact mean intensity of the light field involves the expectation value of products of φ3D and φ † 3D that are not normally ordered, and one has to put the atomic field operators in normal order before the 2D reduction. In the operatorial form of (49), this amounts to replacing kc • [j + 2im ∇n] with (kckp) • J, where J(x, y) = +∞ -∞ dz im φ † 3D (x, y, z)∇ φ3D(x, y, z). Then one finds additive corrections to the above mentioned equations. For example, for kp ≃ kc, in Eq.( 59) one has to add to the expression k B T ρ(r 0 ) m [cos 2 φsc + fn sin 2 φsc] the correction ωzρ(r0)/(2m). This correction, corresponding to spontaneous emission due to the motional coupling between |N C and |C along z, is a flat function of φsc, that can thus be distinguished experimentally from the relevant fn sin 2 φsc component.

[58] As a consequence of the continuity equation, the fluctuations of the density and of the longitudinal current are related by ωnq = q • jq. For small wavevectors, fluctuations are mostly sound-like with ω ≃ cs q. Inserting this result into the continuity equation, one finds that qnq/2m ≪ jL,q if qξ/2 ≪ 1. [72] At the lowest order in Ωp/Ωc, the effect of a nonzero gaa -g ab for the measurement scheme of section V A is to introduce a smooth polarization profile Pp = [4|dae| 2 /( 2 |Ωc| 2 )] (gaa -g ab ) φ †2 3D φ 2 3D Ep. As the corresponding phase shift does not depend on Ωp/Ωc, it can be isolated in the experiment. The contribution of higher order terms O[|Ωp/Ωc| 3 ] have a non-trivial spatial structure and may interfere with the signal to be measured. Their relative value as compared to ∆φ2 is however a factor roughly [(gaa -g ab )/gaa] g (ρλ 2 c )/(2π 2 fn) weaker.

[73] We have checked that the right-hand side of Eq.(C13) is simply -∆F (within linear response theory), where ∆F is the free-energy of the system at thermal equilibrium in presence of the static gauge field A(r, 0 + ) minus the free energy in the absence of gauge field, for fixed temperature and atom number. This establishes a formal link between our proposal and the one of [9].

[74] We used the fact that, for a generic function F satisfying F [-(k + q)] = F (k) for all k, one has