Dariusz Kajtoch

Emilia Witkowska

Alice Sinatra

Spin-squeezed atomic crystal

Keywords: Dariusz Kajtoch, Emilia Witkowska, Alice PACS 03.75.Gg -Entanglement and decoherence in Bose-Einstein condensates PACS 03.75.Kk -Bose-Einstein condensation dynamic properties PACS 42.50.Dv -Squeezed states PACS 06.30.Ft -Time and frequency

HAL is

Introduction. -One of the most successful applications of cold atoms is atomic clocks that provide the best time standards. The intrinsic uncertainty in the measured clock transition frequency, so called standard quantum limit ∆ω ∝ (T √ N ) -1 where T is the interrogation time and N the number of atoms, has been already reached in atomic fountains using a microwave transition [1]. Increasing the atom number N and the interrogation time T in order to lower the value of the standard quantum limit is not straightforward in an atomic fountain. Indeed too large densities introduce an atom-number dependent frequency shift due to atomic interactions (collisional shift) and the interrogation time is limited by the fact that the cold cloud is in free expansion. Optical lattices that allow to confine and interrogate simultaneously a large number of atoms for longer times offer in this respect a crucial advantage [2][3][4][5][6][7] that, combined with the use of a larger reference frequency (optical transition) already outperforms the regular frequency standards [2,8,9]. Interestingly, many-body physics of cold atoms in an optical lattice opens new perspectives that can further improve atomic clocks. In particular, the configuration with one atom per site of an optical lattice is a mean to suppress the collisional shift, maximizing at the same time the density, both for bosonic [3][4][5] and also for fermionic atomic clocks [10]. On the top of these rapid developments, in this paper we envisage to introduce well designed quantum correlations among the atoms internal states, known as spin-squeezing [11][12][13][14][15][16], in order to beat the standard quantum limit and push even further the extraordinary precision of optical-lattice clocks in the long term. To this aim we generalize a spin-squeezing scheme successfully implemented in bimodal Bose-Einstein condensates [12,13] to the spatially multimode case of an optical lattice.

Our idea is to raise the optical lattice in an interacting bimodal Bose-Einstein condensate, adiabatically bringing the system from the superfluid to the Mott-insulator phase [17][18][19]. The correlations among the atomic spins build up in the superfluid phase where "each atom sees each atom" similarly to what happens in the absence of lattice where all the atoms share the same spatial mode. As the system approaches the Mott transition, the squeezing dynamics slows down to finally stop completely. While the condensate is destroyed as the atoms get localized in the lattice sites, the spin-squeezing survives and it is stored in the Mott-insulator phase. To the advantage of having a close packed ensemble atoms, where atomic interactions are highly suppressed meeting the basic requirements for an atomic clock, we then add the advantage of quantum correlations allowing for metrological gain with respect to independent atoms. We analyze in detail the scheme for a microwave transition relevant to microwave trappedatoms bosonic clocks [6,7,[20][21][22], for which a proof-ofprinciple experiment may be preformed with state of the art techniques.

Besides atomic clocks, the scheme can be used to "extract" the entanglement among the atoms that are initially in a common spatial mode [23], to obtain a state where (a) b) ( Fig. 1: (a) Initially, a Bose-Einstein condensate of atoms in an internal state a is prepared in a shallow 3D optical lattice. At time t = 0 an electromagnetic π/2-pulse puts each atom in a coherent superposition two internal states a and b, and the binary atomic interactions between cold atoms start the generation of squeezing in the system [27,28]. (b) Simultaneously, the lattice height is gradually increased, in such a way that the system enters the Mott-insulator phase at the "best squeezing time" t best for which squeezing is the largest, thus freezing the "best squeezing" ξ 2 best = ξ 2 (t best ) (see equation ( 5)) in the Mott-insulator phase.

each entangled atom can be individually addressed and manipulated with a quantum gas microscope [24,25]. Our "spin-squeezed crystal" is then a platform that allows a complete characterization of the entangled state [26] and could be used as a resource in the emerging field of quantum networks and multiparameter estimation.

The Model. -The protocol is sketched in Fig. 1. We consider a two component condensate with repulsive interactions, with symmetric coupling constants describing s-wave interactions between atoms in the two states g aa = g bb , and an adjustable interspecies coupling [13,29] in the phase-mixed regime g ab < g aa [30]; g ij = 4π 2 a ij /m where a ij is the s-wave scattering length for one atom in state i and one in state j, m is the atomic mass and is the Planck constant. The system Hamiltonian is

Ĥ = σ=a,b d 3 r Ψ † σ (r) - 2 2m ∇ 2 + V (r) Ψσ (r) + g σσ 2 d 3 r Ψ † σ (r) Ψ † σ (r) Ψσ (r) Ψσ (r) + g ab d 3 r Ψ † a (r) Ψ † b (r) Ψb (r) Ψa (r), (1)
where Ψa and Ψb are the bosonic field operators for atoms in the state a and b respectively satisfying usual commutation relations,

[ Ψσ (r), Ψ † σ ′ (r ′ )] = δ(r -r ′ )δ σσ ′ . ( 2
)
The system is confined in a 3D uniform optical lattice described by the periodic potential

V (r) = V 0 α=x,y,z sin 2 (kα) (3)
where k = 2π/λ is the lattice wavenumber and λ/2 is the lattice period, with the number of lattice sites M equal to the number of atoms N . Uniform systems can be realized in a flat bottom potential that is now possible to produce in the laboratory [31,32]. The fact that the potential walls are not infinitely steep releases the constraint of having exactly a filling factor of one to reach the Mott transition [33]. After each atom is prepared in a coherent superposition of internal states a and b by an electromagnetic π/2 pulse, the lattice depth is linearly increased in time

V 0 (t) = V init + (V c -V init )t/t best (4)
from an initial value V init to the critical value V c to reach the Mott transition, on a time scale given by the best squeezing time t best . We assume that the system dynamics is confined to the lowest Bloch band and it is described by the two-component Bose-Hubbard model [34] with timedependent hopping J(V 0 ) and interaction terms U σ (V 0 ), U ab (V 0 ). We quantify the spin-squeezing using the parameter [11,35]

ξ 2 = N ∆ Ŝ2 ⊥ min S 2 , ( 5
)
where S is the length of the mean collective spin and ∆ Ŝ2 ⊥ min is the minimal variance of the spin orthogonally to the mean spin direction. The collective spin operators, Ŝ± = Ŝx ± i Ŝy = d 3 r Ψ † a (r) Ψb (r), Ŝz = ( Na -Nb )/2 where Nσ is the atom number operator in the internal state σ, give access to spin-squeezing as a function of time.

Initial state and evolution. -Starting from a condensate in internal state a, a mixing π/2 pulse puts the system in a superposition of Fock states with N a atoms in the internal state a and

N b = N -N a atoms in b |Ψ 0 Na,N b (0 + ) : N a , N b
where we assume that each Fock state is in the spatial ground state for the corresponding atom number in the shallow lattice. We thus neglect the thermal excitations initially present in the system and the excitations created by the mixing pulse. A clean way to avoid any excitation in the experiment would be to perform the pulse for an ideal gas at zero temperature in a completely flat potential, and subsequently ramp up adiabatically the interactions [36] and the lattice. The very stringent condition of having no excitations at all is however not necessary in practice for our proposal. Based on our studies of spin-squeezing at finite temperature [37,38], we expect that the final squeezing will not be affected as long as the initial non-condensed fraction before the pulse is much smaller than the targeted squeezing, which can be achieved in the weakly interacting limit at sufficiently low temperature.

In the adiabatic approximation, during the evolution following the mixing pulse, each Fock state of the prepared superposition remains in an instantaneous ground state and picks up a time-dependent phase factor determined by the ground state energy E 0 (N a , N b , t) 6) for the system evolution is exact in the adiabatic approximation, as long as Ψ 0

|ψ(t) = 2 -N/2 N Na=0 N N a e -i t 0 dt ′ E0(Na,N -Na,t ′ ) × |Ψ 0 Na,N -Na (t) : N a , N -N a . (6) Equation (
Na,N b and E 0 (N a , N b ) refer to the exact ground-state of the system with N a atoms in the internal state a and N b atoms in state b. The N -body ground state for two-component bosons on a lattice exhibits non-trivial correlations [39,40] that entangle the internal and external degrees of freedom in equation ( 6) which would render the calculation of the squeezing parameter a formidable task. We can however simplify the problem by introducing two approximations. (i) We assume that the squeezing is measured after raising the lattice to reach the deep Mott phase with one atom per site in the limit of vanishing tunnel coupling, which allows to disentangle the internal and external degrees of freedom in state (6). A consequence of this fact is that the action of the collective spin operators on the Fock states then obeys simple rules, e.g.

Ŝ+ |Ψ 0 Na,N b : N a , N b ≈ (N a + 1)N b |Ψ 0 Na+1,N b -1 : N a +1, N b -1
, as for the twomode Fock state where only internal degrees of freedom are considered. (ii) We evaluate the energies E 0 (N a , N b ) entering in the phase factors in (6) in the mean field approximation for the external degrees of freedom of the atoms, using the Gutzwiller method [34,[41][42][43]. Precisely, the total mean-field energy is minimized numerically in the subspace of fixed mean number of atoms in each component, i.e. Na = N a and Nb = N b , by a projection method similarly to Ref. [44]. This approximation is accurate enough in 3D and captures the transition to the Mott state marked by the suppression of on-site atom number fluctuations for integer filling and a threshold value of the U/J parameter of the Bose-Hubbard model. The results for the on-site correlation functions g

(2) aa = â † i â † i âi âi / â † i âi 2 and g (2) ab = â † i b † i âi bi / â † i âi b †
i bi and for the squeezing parameter ξ 2 as a function of time are shown in Fig. 2(a) and in Fig. 2(b) (green solid line) respectively. The time where the on-site g (2) functions go to zero, indicated with a vertical dotted line, marks the transition to the Mott phase. At this point the atoms cease to "feel each other" via the mean field interaction. As a consequence the squeezing dynamics stops and the spin correlations that built up in the superfluid phase are "frozen" in the Mott phase.

Effective χ Ŝ2

z model in the superfluid phase. -If the atom number is large, one can expand the ground state energy E 0 (N a , N b , t) in the phase factor of each Fock state in (6) around Na = Nb = N/2 up to the second order [45] 1

t 0 dt ′ E 0 (N a , N b , t ′ ) ≃ Φ 0 (N, t) + T (t) (N a -N b ) 2 4 . (7)
There is no linear term in (7) because of the a -b symmetric situation we consider, and the phase factor Φ 0 (N, t) depending on the total number of atoms does not play any role in the spin dynamics and can be neglected. Introducing the chemical potentials µ σ (t) = ∂ Nσ E 0 (N a , N b , t) and the parameter

χ(t) = 1 2 (∂ Na -∂ N b )[µ a (t) -µ b (t)
]/ , the function T (t) in equation ( 7) has the form

T (t) = t 0 dt ′ χ(t ′ ) . (8)
By changing the time variable from t to the dimensionless T (t), one then recovers the one-axis twisting (OAT) model [27] with an Hamiltonian proportional to Ŝ2 z . In particular, in the large N limit, the squeezing optimized over time

ξ 2 best = ξ 2 (t best ) is ξ 2 best ≃ 3 2/3 2 1
N 2/3 and the best squeezing time is

t best ≃ 3 1/6 N 2/3   1 V c -V init Vc Vinit dV 0 χ(V 0 )   -1 (9)
where we used T (t best ) ≃ 3 1/6 /N 2/3 and the linear ramp (4). We introduce the term "dynamic-OAT model" for the quadratic approximation (7) in the phase factors, as opposite to "static-OAT model" for which we directly take a time independent Hamiltonian Ĥ = χ(0) Ŝ2 z corresponding to the initial conditions after the pulse. Contrarily to the approach described in the previous section based on the numerical solution of equation ( 6), the dynamic-OAT model is valid only in the superfluid phase where the derivatives entering the definition of χ(t) are well defined. It allows however to get a physical insight and simple analytical results. In Fig. 2(b), giving the time evolution of the squeezing parameter, we show that the dynamic-OAT model (red circles) is in agreement with the numerical solution of equation ( 6) (green solid line) up to the phase transition already for N = 125. In Fig. 2(c) we show the corresponding time dependence of T (t). When approaching the transition, χ(t) tends to zero and T (t), that is the effective time of the effective χ Ŝ2 z model, tends to a constant. For optimal results, the Mott transition should occur at the best squeezing time, whose expression for large N is given in equation (9). As expected, the squeezing dynamics is slower in the dynamic-OAT model than in static one. We show the scaling of the best squeezing time with the total number of atoms N in Fig. 2(d). A fit gives t best ∝ N 0.35 , a slightly less favorable scaling than for the homogeneous static-OAT model, for which χ(0) homo = 8(a aa + a bb -2a ab )E R /πλ N and t best ∝ N 1/3 . The same figure confirms that the dependence of t best on the scattering lengths suggested by the homogeneous static-OAT model, approximately holds.

Adiabaticity and beyond mean field effects. -In our treatment based on equation ( 6), while we treat exactly the spin degrees of freedom that are in a quantum superposition in the initial state, we perform two main approximations concerning the external degrees of freedom of the atoms. First, we assume adiabaticity as the optical lattice lattice is raised, second we use the Gutwiller approximation to evaluate the ground state energy for a given spin Fock state. In this paragraph, we extend our analysis beyond mean field. As the squeezing in our scheme develops essentially in the superfluid phase, we first con-centrate on this phase and derive the validity conditions of the adiabatic approximation at heart of our treatment using the number-conserving Bogoliubov theory [46,47]. Then, in the end of the pararagraph, we discuss the effect of residual density fluctuations in the Mott state.

In the superfluid phase, where delocalization of atoms is energetically more favorable, a large majority of particles occupy the zero quasi-momentum Bloch state of the lowest band giving rise to condensation in momentum space. When the depletion of the condensate is small, the Bogoliubov method can capture the physical properties of the Bose-Hubbard model [48], although it cannot be pushed too far towards the phase transition boundary [48,49]. The starting point is the twocomponent Bose-Hubbard Hamiltonian [34] written in quasi-momentum representation. This Hamiltonian has the same form as for a two-component Bose-Einstein condensate in free space [50] with the field operators Ψσ (r) expanded in a plane wave basis, provided the kinetic energy ǫ q = 2 q 2 /2m and the coupling constants g σ , g ab are replaced by ǫ q (t) = -2J(V 0 ) γ cos(q • e γ ) and U σ (V 0 ), U ab (V 0 ) respectively, where e γ are primitive lattice vectors. One then introduces the number-conserving operators, Λq,a = â † 0 âq / √ N a and Λq,b = b † 0 bq / √ N b representing the non-condensed fields, and quadratizes the Hamiltonian in these fields. The quadratic Hamiltonian is diagonalized using the generalized Bogoliubov transformation Λq,σ (t) =u σ q,+ (t) βq,+ (t)

+ v σ q,+ (t) β † -q,+ (t)+ u σ q,-βq,-(t) + v σ q,-(t) β † -q,-(t), (10)
where βq,± (t) are Bogoliubov quasi-particle operators satisfying bosonic commutation relations and u σ q,± (t), v σ q,± (t) are known functions [30,50]. The Bogoliubov transformation in Eq. ( 10) is the same as in [50] with the same notations, provided that equations are expressed in terms of ∆E q (t) = ǫ q (t) -ǫ 0 (t).

In order to calculate the number of excitations created by the lattice rump starting at zero temperature with no excitation, we write the Heisenberg equation of motion for βq,± (t), and proceed similarly to appendix C of Ref. [36]. One can show that the Bogoliubov modes βq,+ (t) and βq,-(t) evolve independently, and are coupled to modes β-q,+ (t) and β-q,-(t) respectively

i d dt βq,± (t) β † -q,± (t) = ω q,± (t) -i Ω q,± (t) -i Ω q,± (t) -ω q,± (t) × βq,± (t) β † -q,± (t) , (11)
where the Bogoliubov spectrum is defined by ( ω q,± /∆E q ) 2 = 1 + ũa + ũb ± (ũ a -ũb ) 2 + 4ũ ab with ũσ = U σ n σ /∆E q and n σ = N σ /M for σ = a, b ,

ũab = U ab √ n a n b /∆E q and Ω q,± (t) = 1 2 d dt log ∆Eq(t)
ωq,±(t) .

p-4

Spin-squeezed atomic crystal

The total fraction of excitations 1 N q =0 n ex q (t) with n ex q (t) = ǫ=± Ψ bog (0)| β † q,ǫ (t) βq,ǫ (t)|Ψ bog (0) stays small, as long as the (sufficient) adiabaticity condition

|Ω q,± (t)| ≪ |2ω q,± (t)| (12)
is verified for each quasi-momentum q = 0. From the expression of Ω q,± (t) and for the linear ramp (4) one gets the expression of the adiabatic time for each mode

t q adiab,± = (V c -V init ) E R 4∆E q (V 0 ) d dV 0 ∆E q (V 0 ) ω q,± (V 0 ) (13)
which should be evaluated at q = (2π/N 1/3 l, 0, 0) with the lattice spacing l = λ/2, where the condition ( 12) is most stringent, and maximized over the ramp duration

V 0 ∈ [V init , V c ].
In the large atom number limit one gets

t adiab,± ≃ N →∞ N 1/3 2π (V c -V init )l 4J(V 0 ) d dV 0 J(V 0 ) c ± (V 0 ) , (14)
where c 2 ± = (l/ ) 2 J(V 0 )(( ω q,± /∆E q ) 2 -1)∆E q is the sound velocity of the phonon-like Bogoliubov excitation branches. We see from Eq. ( 14) that the adiabatic time shares almost the same scaling with N as the best squeezing time. Furthermore, the negative branch t adiab,-is always larger than t adiab,+ as long as one of the components is not completely depleted. This implies that t adiab,-alone sets the adiabatic time scale. We checked numerically that the condition t best ≫ t adiab,-holds for all parameters that we consider. For example, for N = 10 4 and parameters in Fig. 2, t adiab /t best ≈ 0.032.

Decoherence. -While two and three-body losses are suppressed in the Mott-insulator phase, they cannot be neglected in the superfluid phase, and play an important role as soon as the lost fraction of atoms at the best squeezing time becomes comparable to the squeezing ξ 2 best that one would have in the absence of decoherence [51]. For two particular configurations in Fig. 3 we show the lost fraction and the expected squeezing ξ 2 best in the absence of decoherence as a function of the initial atom number. As the lost fraction increases with N , while ξ 2 best decreases, the crossing of the two curves gives the maximum atom number N max for which the decoherence due to atom losses can be neglected. The lost fraction is obtained by solving rate equations for the mean atom number

d dt Nσ (t) = -γ (2) σ (t) g (2) σ (t) Nσ 2 (t) -γ (2)
ab (t) g (2)
ab (t) Na (t) Nb (t) (15) for two-body losses, and

d dt Nσ = -γ (3) σ (t) g (3) σ (t) Nσ 3 (16)
for three-body, where γ (m)

σ (t) = K (m) σ M m-1 d 3 r w 2m (r, t)
and γ

(2) a = a † i a † i a i a i / a † i a i 2 , calculated numerically in the adiabatic approximation using the Gutzwiller method. In the case of the transition |F = 1, m F = 1 ↔ |F = 2, m F -1 for 87 Rb atoms in Fig. 3(a), where the interspecies scattering length a ab can be tuned by Feshbach resonance [13], decoherence due to two-body losses is never negligible above N max ≈ 40 atoms, and limits the squeezing to 10dB already for a hundred atoms. On the other hand, for the transition |F = 1, m F -1 ↔ |F = 2, m F = -2 in Fig. 3(b), where the interspecies interaction can be tuned by slightly shifting the optical lattices for the two components [29] and two body losses are absent [56], the limit that is now imposed by three body losses is much less constraining and it does not influence the results up to large atom numbers N max ≈ 4 × 10 5 .

ab (t) = K (
In the Mott phase with one atom per site of a twocomponent system, two types of excitations exist: excitations that lead to double occupations, that are in general gapped, and "soft" excitations within the subspace of one atom per site that can be described by an effective spin model. In our scheme we need to maintain adiabaticity with respect to the first kind of excitations only, that is we need to remain in the manifold of single occupation. In this manifold, it is not important for us to be in the ground state, because the spin-dependent interaction energy is strongly reduced, scaling as J 2 /U [39,40,57,58] 1 , and consequently the squeezing dynamics is practically stopped. In particular, contrarily to what happens in the superfluid phase, the presence of a gap separating the low energy manifold with one atom per site from the states with double occupation, allows us to raise the lattice to the deep Mott regime J → 0 in a finite time, typically of order τ ≥ /U [59]. One can estimate the effect of the residual squeezing dynamics during this time by comparing χ resid τ , where χ resid scales as the spin-dependent interaction eneregy per particle divided by N , with the expression of the adimensional squeezing time T (t best ) ≃ N -2/3 below equation (9). As χ resid τ ≃ (J/U ) 2 N -1 ≪ N -2/3 we find that the residual squeezing dynamics during this adiabatic time is completely negligible. The only problem may come from the critical region where the gap closes as ∆/J ∝ (J/U ) c -(J/U ) in the thermodynamic limit, scaling as (∆/J) c ∝ N -1/3 at the transition in a finite size system [60]. At the critical point we then find the same scaling N 1/3 of the adiabatic time as in the superfluid case. However, even taking τ = N 1/3 /J, we would have χ resid τ ≃ (J/U )N -2/3 ≪ N -2/3 meaning that the effect of the residual squeezing dynamics would still be small, due to the small value of J/U near the critical point (J/U ) c ≃ 0.03.

Conclusions. -We study the formation of a spin-squeezed atomic crystal by bringing an interacting two-component Bose-Einstein condensate across the superfluid-to-Mott transition in a time-dependent optical lattice. The scheme could be directly used in lattice clocks using a microwave transition.

Fig. 2 :

 2 Fig. 2: (Color online) Numerical results: (a) On-site two-body correlation functions g (2) aa = â † i â † i âiâi / â † i âi 2 and g (2) ab = â † i b † i âi bi / â † i âi b † i bi calculated for the central Fock state |Na = N/2, N b = N/2 . (b) Spin-squeezing parameter ξ 2 as a function of t/tunit (bottom x-axis) and V0/ER (top x-axis) for the linear ramp; ER = 2π 2 2 /(mλ 2 ) is the recoil energy and λ/2 the lattice period. t -1 unit = (aaa + a bb -2a ab )ER/ λ. The green solid line shows numerical results with E0 in (6) calculated with the Gutzwiller method, while the red dots result from the dynamic-OAT model (7). Predictions of the static-OAT model (see text) are shown for comparison (blue dashed lines). (c) T (t) for the dynamic-OAT model (red solid line) and the static-OAT model (blue dashed line). Parameters are N = 125, aaa = a bb = 100.4a0 and a ab = 95a0, where a0 is the Bohr radius. The dashed horizontal lines in (b) and in (c) represent respectively ξ 2 best and T (t best ) of the OAT model for N = 125. By definition, the dynamic curves touch the dashed horizontal lines at the best squeezing time t best . (d) Scaling of the best squeezing time t best with N for aaa = a bb = 100.4a0 and different values of a ab . The dashed line is a fit t best ∝ N α with α = 0.353 ± 0.004.

Fig. 3 :

 3 Fig. 3: Best squeezing ξ 2 best = ξ 2 (t best ) in the absence of decoherence and lost fraction due to 2-body and 3-body losses at the time t best as function of the atom number in 87 Rb for the transitions (a) |1, 1 ↔ |2, -1 with Feshbach tuned a ab = 95a0 and dominant two body losses K (2) b = 8.1 × 10 -20 m 3 /s [52] and K (2) ab = 1.708 × 10 -19 m 3 /s [53] and (b) |1, -1 ↔ |2, -2 with three body losses only K (3) a = 5.4 × 10 -42 m 6 /s [54] and K (3) b = 1.8 × 10 -41 m 6 /s [55].

In our case the three interaction parameters U σσ ′ have the same order of magnitude, 0 < U ab Uaa, U bb , so that we will refer to a generic U in the following discussion.

We wish to thank Yvan Castin and Boris Svistunov for useful discussions. This work was supported by the Polish National Science Center Grants DEC-2015/18/E/ST2/00760, by the CNRS PICS-7403, and partially by the PL-Grid Infrastructure. D.K. acknowledges support from the French Government (BGF).


Dariusz Kajtoch

Emilia Witkowska

Alice Sinatra

Adiabaticity when raising a uniform 3D optical lattice in a bimodal Bose-Einstein condensate

Using the time-dependent Bogoliubov approach, we study adiabaticity for a two-component Bose-Einstein condensate in a 3D time-dependent optical lattice with unit filling, in the superfluid and weakly interacting regime. We show that raising the lattice potential height can couple the ground state of the Bogoliubov Hamiltonian to excited states with two quasiparticles of opposite quasimomenta. In the symmetric case for interactions and density in the two components these represent sound waves where the two components oscillate out of phase. We find an analytic expression of the adiabatic time, its dependence on the fraction of atoms in each component and its scaling with the system size.

I. INTRODUCTION

Spinor bosonic or fermionic atoms in optical lattices are playing an increasingly important role at the crossing of different fields such as statistical physics, condensed matter, atomic physics and quantum technologies. Actively investigated in present experiments, these systems display non trivial phase transitions and ground states stemming from the interplay between the spin and the external degrees of freedom of the atoms [1], and can be used to investigate novel superfluidity mechanisms [2]. Furthermore, they constitute a powerful platform for quantum computation [3,4], and offer fascinating perspectives for entangled state preparation and quantum metrology [5,6]. Among the different proposals using cold atoms in an optical lattice, several protocols require the possibility to adiabatically ramp up the optical lattice in the multi-component cold atoms system.

In close relation with the first experimental realizations, adiabaticity has been mainly studied, both experimentally [7] and theoretically [8,9] for a single species and in the presence of an external harmonic potential that brakes the translational symmetry of the lattice. In this paper we concentrate on the case of a uniform optical lattice, that is now possible to prepare in the laboratory thanks to the development of flat-bottom potentials [10,11]. We shall consider a single-component or a two-component condensate in the superfluid regime and study the adiabaticity condition when raising the lattice within the time-dependent Bogoliubov approach. For a quantum system with a discrete spectrum ĤΨ k = E k Ψ k , initially in its ground state Ψ 0 , the adiabaticity condition when a parameter of the Hamiltonian is varied in time starting from t = 0 takes the form [START_REF] Messiah | Mécanique quantique[END_REF]

Ψ k (t)| d dt Ψ 0 (t) |E k (t) -E 0 (t)| 1 ∀k = 0, ∀t > 0 (1)
Ψ k (t)| d dt Ψ 0 (t) is the coupling between the instantaneous ground state and other eigenstates, and E k (t) -E 0 (t) is the energy difference. The condition (1) can be equivalently rewritten as

Ψ k (t)|i d Ĥ dt |Ψ 0 (t) |E k (0) -E 0 (t)| 2 1 ∀k = 0, ∀t > 0 (2)
We will use this last formulation to interpret our result. Our paper is structured as follows: after recalling the Bose-Hubbard model and introducing all notations, the adiabaticity criterion for a single component Bose-Einstein condensate is derived in section II. The analysis is extended to the case of a two-components system in section III. While in this section we keep our formalism general, allowing to address the case where the atoms in the two components might have different masses, in the following we concentrate on the case of equal masses corresponding to different hyperfine states of the same atomic species. For this case, in section IV, we derive the adiabatic time for a linear ramp and we investigate the influence of a density imbalance between the two components. Conclusions are drown in section V.

II. ADIABATICITY CRITERION FOR A SINGLE COMPONENT

In the adiabatic evolution, when an external parameter of the Hamiltonian is changed in time, the quantum state remains an instantaneous eigenstate of the timedependent Hamiltonian at all times. On the contrary, if the change is too fast, the state will contain an admixture of excited states. For a Bose-Einstein condensate in an optical lattice in the weakly interacting regime, when the system is superfluid, the excitations are well described as Bogoliubov quasi-particles. We will regard the evolution as adiabatic as long as the total density of excited quasiparticles remains much less than one at all times.

A. Single component Bose-Hubbard model

We describe a system of N ultra-cold atoms, all in the same internal state and subject to an optical lattice potential, by the Bose-Hubbard Hamiltonian [START_REF] Jaksch | [END_REF]

Ĥ(t) = -J(t) i,j â † i (t)â j (t) + U (t) 2 i ni (t)[n i (t) -1],
(3) where â † i (t) creates a particle in the single-particle Wannier state w i (r, t) of the lowest energy band (l = 1) localized on the i-th site. The Bose-Hubbard model (3) considers only states in the lowest energy band, which is justified as long as the excitations energies to the higher bands are much larger than the energies involved in the system dynamics. The Wannier states w i (r, t) are conveniently constructed from the lowest band Bloch states ψ l=1,q (r, t) in the following way

w i (r, t) = 1 √ M q∈BZ e -iqR i ψ l=1,q (r, t), (4)
with

R i = d(i x u x + i y u y + i z u z )
, where d = λ/2 is the lattice spacing that we assume identical in the three spatial directions (λ is the optical lattice wavelength), u α are unit vectors and i x , i y , i z are integers. Summation in Eq. ( 4) extends over wave vectors belonging to the 1st Brillouin zone and M is the number of lattice sites. The Bloch states ψ l,q (r, t), labeled by the band index l and the quasi-momentum q, are eigenstates of the singleparticle Hamiltonian ĥ(t) = -

2 2m ∇ 2 + V 0 (t) α=x,y,z sin 2 (kα), (5)
where k = 2π/λ and V 0 (t) is the lattice potential height. If V 0 (t) is varied in time, the Wannier states and hence the creation and annihilation operators in the Bose-Hubbard Hamiltonian (3), as well as the hoping J(t) and interaction U (t) parameters, depend on time

U (t) = 4πa s 2 m d 3 r |w(r, t)| 4 , (6a)
J(t) = d 3 r w * i (r, t) - 2 2m ∇ 2 +V 0 (t) α=x,y,z sin 2 (kα) w j (r, t), (6b)
where a s is the s-wave scattering length characterizing binary short range interactions between cold atoms and m is the mass of an atom. In the limit in which V 0 E R where E R = 2 k 2 /(2m) is the recoil energy, the dependence of U and J on the lattice height V 0 can be approx-imated by [14,15]

U (V 0 ) = E R 8 π ka s V 0 E R 3/4 , (7a)
J(V 0 ) = E R 4 √ π V 0 E R 3/4 e -2 V 0 E R . (7b)
For brevity, we shall omit in the following to mark the band index l = 1 for the Bloch states of the lowest band.

B. Number-conserving Bogoliubov approach

We start from the Bose-Hubbard Hamiltonian (3) written in the quasi-momentum representation

Ĥ(t) = q q (t)â † q âq + U (t) 2M q1,q2,k â † q1-k â † q2+k âq1 âq2 , (8)
contrarily to the homogeneous case, the kinetic energy in the lattice takes the form

q (t) = -2J(t) α=x,y,z cos(q • du α ). (9)
In Eq. ( 8) and further we omit the explicit time dependence of the creation and annihilation operators to simplify the notation, but we keep in mind that, even in the Schrödinger picture, â † q and âq depend on time. In the number conserving Bogoliubov approach the amplitude of the field in the condensate mode is finally eliminated, and the interacting system is described as an ensemble of quasiparticles in the modes orthogonal to the condensate mode. The small parameter of the theory is the non-condensed fraction, and to the lowest non-zero order, that is the Bogoliubov order, the quasiparticles do not interact. The first step to find the Bogoliubov quasiparticles is to quadratize the Hamiltonian with respect to the non-condensed field. One then obtains

Ĥ(t) N 0 (t) + U (t) 2 n + q =0 â † q âq [ q (t) -µ(t)] + U (t) 2M q =0 â † 0 â † 0 âq â-q + â0 â0 â † q â † -q + 4N â † q âq ( 10
)
where n = N/M is the total atom density and the chemical potential is defined as

µ(t) = 0 (t) + U (t)n. (11)
To obtain Eq. ( 10), the relation â † 0 â0 = N -δ N , where δ N = q =0 â † q âq , was used. We now introduce the number conserving operators

Λq = 1 √ N â † 0 âq ( 12
)
in terms of which the Bogoliubov Hamiltonian (10) takes its final form

ĤBog (t) = H 0 (t) + 1 2 q =0 Λ † q , Λ-q σ z L q Λq Λ † -q , (13)
where the ground state energy is

H 0 (t) = N 0 (t) + U (t) 2 n - 1 2 q =0 [ q (t) -µ(t) + 2U (t)n] , (14)
the non-hermitian matrix L q has the form

L q = q (t) -µ(t) + 2U (t)n U (t)n -U (t)n -[ q (t) -µ(t) + 2U (t)n] , (15)
and σ z is the third Pauli matrix. In the derivation of (13) we used the following approximation

Λ † k Λq = â † k âq â0 â † 0 N ≈ â † k âq , (16)
consistent with the fact that we retain only terms that are at most quadratic in the non-condensed field. The Bogoliubov transformation

T q (t) bq b † -q = ūq (t) -v q (t) -v q (t) ūq (t) Tq(t) Λq Λ † -q , (17)
diagonalizes the quadratic Hamiltonian [START_REF] Jaksch | [END_REF] in terms of the Bogoliubov operators bq and b † q that satisfy bosonic commutation relations [ bq , b † q ] = δ q,q ĤBog (t) =H 0 + 1 2

q =0 ω q (t) + q =0 ω q (t) b † q bq , ( 18
)
where the Bogoliubov energy has the form

ω q (t) = ∆E q (t) [∆E q (t) + 2μ(t)] , (19)
∆E q (t) = q (t) -0 (t) , (20) μ
(t) = µ(t) -0 (t) = U (t)n , (21)
and the Bogoliubov modes are ūq (t) + vq (t) = ∆E q (t) ∆E q (t) + 2μ(t)

1/4 , (22)
ūq (t) -vq (t) = ∆E q (t) ∆E q (t) + 2μ(t) -1/4 . ( 23
)
In the limit of a small quasi-momentum dq 1 and for ∆E q U (t)n, the Bogoliubov energy ( 19) has a phonon-like dispersion

ω q (t) = c(t)q (24)
where c(t) = a 2J(t)U (t)n/ is the sound velocity. The Bogoliubov spectrum for a condensate in a uniform lattice was already found in [16] within the usual symmetrybreaking approach.

C. Time evolution in the Heisenberg picture

The column vector composed of Bogoliubov quasiparticle annihilation and creation operators

Bq (t) = bq (t) b † -q (t) , (25)
evolves according to the Heisenberg equation of motion

d dt Bq (t) = i [ ĤBog (t), Bq (t)] + ∂ Bq (t) ∂t H . ( 26
)
The first part on the right-hand side of ( 26) represents the free evolution of the quasi-particles

[ ĤBog (t), Bq (t)] = -ω q (t)σ z Bq (t) . ( 27
)
According to the Bogoliubov transformation ( 17), the second part on the right-hand side of ( 26) reads

∂ Bq (t) ∂t H = dT q (t) dt T -1 q (t) Bq (t)+T q (t)∂ t Λq (t) Λ † -q (t) (28
) By using the identity ū2 q (t)-v2 q (t) = 1 one can show that

dT q (t) dt T -1 q (t) = -Ω q (t)σ x , (29)
where

Ω q (t) = ūq (t) d dt vq (t) -vq (t) d dt ūq (t), (30)
and σ x is the first Pauli matrix. The coupling Ω q (t) can be expressed in terms of the quasi-particles energies, and takes the form

Ω q (t) = d dt log(ū q + vq ) = 1 2 d dt log ∆E q (t) ω q (t) . (31)
Let us now deal with the second term on the right-hand side of Eq. ( 28). The time derivative of the numberconserving operator ( 12) is

∂ t Λq (t) = 1 √ N ∂â † 0 (t) ∂t âq (t) + â † 0 (t) ∂â q (t) ∂t . (32)
By using the definition of the creation and annihilation operators in the first Bloch band

âq (t) = d 3 r ψ * q (r, t) Ψ(r), (33)
and expanding the field operator on the complete set of time-dependent Bloch states

Ψ(r) = k ψ k (r, t)â k (t) + l =1,k ψ l,k (r, t)â l,k (t) (34) one obtains ∂ t âq (t) = k C * k,q (t)â k (t) + l =1,k C * l,k,q (t)â l,k (t), (35)
where

C k,q (t) = d 3 r ψ * k (r, t) ∂ ∂t ψ q (r, t), (36)
C l,k,q (t) = d 3 r ψ * l,k (r, t) ∂ ∂t ψ q (r, t). ( 37
)
the first sum on the right-hand side of Eq. ( 35) runs within the lowest energy band, while the second one runs over all the other bands. For ramping times that are long with respect to the inverse recoil frequency /E R , one can neglect all inter-band couplings, i.e. C l,k,q = 0. Moreover, from the conservation of quasi-momentum we know that C k,q (t) = δ q,k C q,q (t) (see Appendix A). The time derivative (32) of the number conserving operators is then

∂ t Λq (t) Λ † -q (t) = C * q,q (t) Λq (t) Λ † -q (t) , ( 38
)
where we used the relations C -q,-q (t) = C * q,q (t) and C 0,0 (t) = 0 [START_REF]From the normalization of the Bloch states one can show that Re{Cq,q(t)} = 0. Moreover[END_REF] for the Bloch states in (36).

Gathering all the terms, the Heisenberg equation ( 26) for Bq (t) takes the form

i d dt Bq (t) = [ ω q (t)σ z -i Ω q (t)σ x ] Bq (t) + i C * q,q (t) Bq (t). ( 39
)
Notice that the second term in the above equation gives a global time-dependent phase factor, which can be removed by a gauge transformation. Finally, the solution of the Heisenberg equation ( 26) can be cast in the simple form

Bq (t) = A q (t) B q (t) bq (0) + B * q (t) A * q (t) b † -q (0), (40)
where

i d dt A q (t) B q (t) = ω q (t) -i Ω q (t) -i Ω q (t) -ω q (t) A q (t) B q (t) , (41)
with the initial conditions A q (0) = 1 and B q (0) = 0.

D. The adiabaticity parameter

During time evolution some Bogoliubov quasi-particles will be excited with respect to the quasi-particle vacuum state. The number quasi-particles created in mode q is determined by the coefficient B q (t):

n ex q (t) = Ψ Bog (0)| b † q (t) bq (t)|Ψ Bog (0) = |B q (t)| 2 , ( 42
)
where |Ψ Bog (0) is the Bogoliubov vacuum state at time t = 0. Excitations are avoided by suppressing the coupling between A q (t) and B q (t) in Eq. ( 41), which brings us to the adiabaticity condition

α q (t) = Ω q (t) 2 ω q (t) 1, ∀q, (43)
where Ω q is the coupling (31) and ω q is the Bogoliubov energy (19). Since both | Ω q (t)| and 1/ ω q (t) are monotonically decreasing functions of q, the left-hand side of Eq. ( 43) reaches its maximum value for minimal quasimomenta |q min | = q min . We thus introduce the adiabaticity parameter α as:

α = max 0≤t≤tramp α qmin (t), (44)
α qmin (t) ≡ 4∆E qmin d dt ∆E qmin ω qmin . ( 45
)
1. Limit of large N By expanding Eq. ( 9), ( 19) and (31) for a small wave vector q and non-zero interactions: one obtains

α qmin (t) qmin→0 = 4J(t) d dt J(t) c(t)q min . ( 46
)
The expression of the minimum wavenumber in the lattice

q min d = 2π N 1/3 , (47)
gives a scaling N 1/3 to the adiabaticity parameter. In the case of a linear ramp changing between V min and V max in a time t ramp

V 0 (t) = V max -V min t ramp t + V min , (48)
taking n = 1, and using the approximated formulas (7a)-(7b) for a deep enough lattice, we deduce the adiabatic time t adiab such that the evolution is adiabatic for t ramp t adiab ,

t adiab = E R V max -V min E R N 1/3 64( √ 2πka s ) 1/2 × max Vmin≤V0≤Vmax E R V 0 5/4 e V 0 E R . ( 49
)
Note that the fact the formula gives a diverging time in the limit V 0 /E R → ∞ (corresponding to J → 0) is not relevant here as our analysis is restricted to the superfluid regime were a condensate is present.

Ideal gas

When interactions tends to zero, i.e. U → 0 we cannot linearize the dispersion relation of quasi-particles. One rather has ω q (t) → ∆E q (t), (50) Ω q (t) → 0 .

(51)

In the ideal gas regime all the particles occupy the instantaneous q = 0 mode during time evolution and the evolution is always adiabatic. This effect is a consequence of quasi-momentum conservation, deriving from the fact the although the lattice height increases, the periodicity of the lattice is unchanged.

E. Interpretation of the adiabatic parameter

In this subsection we explicit the link between Eq. (43) and Eq. ( 2), gaining some physical insight into our result. Let us take the time derivative of the Bogoliubov Hamiltonian ( 18)

i d dt ĤBog (t) = i dH 0 dt +i q =0 d( ω q (t)) dt 1 2 + b † q bq + ω q (t) i d b † q dt bq + b † q d bq dt . (52)
The term in the first line of (52) cannot change the number of quasi-particles and will play no role. For the term in square brackets in the second line, using Eq (39) we find

i d b † q dt bq + b † q d bq dt = -i Ω q b-q bq + b † q b † -q , ( 53
)
showing that the time derivative of the Hamiltonian can couple the Bogoliubov ground state to states with two quasi-particles with opposite quasi-momenta

1 : q, 1 : -q| b † q b † -q |0 = 1, (54)
the energy difference being 2 ω q . Using Eq. ( 2), we then obtain the condition

α q (t) = |Ω q | 2 ω q 1, (55)
that coincides with the adiabaticity condition (43).

III. ADIABATICITY CRITERION FOR TWO-COMPONENTS

In this section we extend the reasoning presented in section II to derive an adiabaticity condition for the twocomponent system, including the case of different masses in the two components, corresponding to different isotopes [START_REF] Papp | [END_REF] or different atomic species [19,20].

A. Two-components Bose-Hubbard model

We consider a two-component Bose-Einstein condensate in an optical lattice potential.

As the atoms in the two components might have different masses, we introduce two tunneling parameters J a (t) and J b (t) of the Bose-Hubbard Hamiltonian, and two kinetic energies q,a , q,b and ∆E q,a , ∆E q,b defined as in equations (6b), (9) and (20). The interactions between atoms may also be different for different components. For convenience we introduce the interaction parameters

U σ (t) = a σ a s U (t), with σ = a, b, ab, (56)
with a σ=a,b,ab the s-wave scattering lengths respectively for two atoms in the state |a , two atoms in the state |b , or one atom in |a and the other in |b , where U (t) is still defined by Eq. (6a). We restrict here to the case of repulsive interactions U σ > 0 for σ = a, b, ab and to the miscible regime U ab < √ U a U b [21][22][23]. We note N a and N b the average number of atoms in each component, and N = N a + N b the total number of atoms.

B. Bogoliubov description for two-components

The Bose-Hubbard Hamiltonian for the twocomponent system in the quasi-momentum representation takes the form

Ĥ2C (t) = q q,a (t)ĉ † q,a ĉq,a + q,b (t)ĉ † q,b ĉq,b + q1,q2,k U a (t) 2M ĉ † q1-k,a ĉ † q2+k,a ĉq1,a ĉq2,a + U b (t) 2M ĉ † q1-k,b ĉ † q2+k,b ĉq1,b ĉq2,b + U ab (t) M ĉ † q1-k,a ĉ † q2+k,b ĉq1,a ĉq2,b , (57)
where ĉ † q,σ is the creation operator in the internal state |σ and quasi-momentum q. Similarly as in the singlecomponent case, we introduce a set of number conserving operators

Λq,a = 1 √ N a ĉ † 0,a ĉq,a , (58)
Λq,b = 1 √ N b ĉ † 0,b ĉq,b . (59)
In terms of the vector

Λ2C,q = Λq,a , Λ † -q,a , Λq,b , Λ † -q,b T , (60)
the two-component Bogoliubov Hamiltonian (57), once quadratized, reads Ĥ2C,Bog (t) = H 2C,0 (t) + 1 2

q =0 Λ † 2C,q Σ z L 2C,q Λ2C,q , (61)
where Σ z = σ z ⊕ σ z , and where the explicit form of the ground state energy H 2C,0 (t) and of the non-hermitian matrix L 2C,q , which are the two component equivalents of ( 14) and ( 15) respectively, is given in Appendix B. Equivalently to [START_REF]From the normalization of the Bloch states one can show that Re{Cq,q(t)} = 0. Moreover[END_REF], the Bogoliubov Hamiltonian (61) can be diagonalized using a Bogoliubov transformation: B2C,q (t) = T 2C,q (t) Λ2C,q ,

where we have introduced the vector B2C,q (t) = ( bq,+ (t), b † -q,+ (t), bq,-(t), b † -q,-(t)) T , (63) bq,± being annihilation operators of Bogoliubov quasiparticles, and the transformation matrix

T 2C,q (t) =     ūa q,+ (t) -v a q,+ (t) ūb q,+ (t) -v b q,+ (t) -v a q,+ (t) ūa q,+ (t) -v b q,+ (t) ūb q,+ (t) ūa q,-(t) -v a q,-(t) ūb q,-(t) -v b q,-(t) -v a q,-(t) ūa q,-(t) -v b q,-(t) ūb q,-(t)     .
(64) As in the homogeneous case without a lattice [START_REF] Castin | Les Houches LXXII:Coherent atomic matter waves[END_REF][START_REF] Tommasini | [END_REF][26][27], there are two excitation branches labeled ±:

ĤBog (t) = H 2C,0 (t) + 1 2 q =0 σ∈{+,-} ω q,σ (t) + q =0 σ∈{+,-} ω q,σ (t) b † q,σ (t) bq,σ (t), (65)
of energies

ω q,± = ω 2 q,a + ω 2 q,b 2
± ( ω 2 q,a -ω 2 q,b ) 2 4 + 4∆E q,a ∆E q,b U 2 ab n a n b   1/2 (66)
where ω q,a(b) are defined as in the single component case (19) ω q,a = ∆E q,a (∆E q,a + 2U a n a ), ( 67)

ω q,b = ∆E q,b (∆E q,b + 2U b n b ), (68)
with n σ = N σ /M . For single-particle unit filling M = N , n a and n b represent the fraction of atoms in component a and b respectively. The corresponding Bogoliubov amplitudes which appear in the transformation (62) are [26]

        ūa q,± va q,± ūb q,± vb q,±         =         2U ab √ n a n b ∆E q,b (∆E q,a + ω q,± ) 2U ab √ n a n b ∆E q,b (∆E q,a -ω q,± ) ( ω 2 q,± -ω 2 q,a )(∆E q,b + ω q,± ) ( ω 2 q,± -ω 2 q,a )(∆E q,b -ω q,± )         χ q,± , (69)
with the normalization coefficient χ q,± given by

χ q,± = [4∆E q,b ω q,± × 4∆E q,a ∆E q,b U 2 ab n a n b + ( ω 2 q,± -ω 2 q,a ) 2 -1/2 (70) is chosen to have σ=a,b |ū σ q,± | 2 -|v σ q,± | 2 = 1.

C. Heisenberg equation of motion

Following the procedure and arguments presented for the single component case in section II, the Heisenberg equation of motion for the operator B2C,q (t), defined in Eq. ( 63) takes the form

i d dt B2C,q (t) = [ ω q,+ (t)(σ z ⊕ 0 2×2 ) + ω q,-(t)(0 2×2 ⊕ σ z )] B2C,q (t) + i T 2C,q (t) dt T -1
2C,q (t) B2C,q (t) + i C * q,q B2C,q (t). ( 71)

The last term on the right-hand side of Eq. ( 71) gives a global time-dependent phase factor, which can be removed by a suitable gauge transformation. In general, the coupling term in the Eq. ( 71) is

T 2C,q (t) dt T -1 2C,q (t) =    0 -Ω q,+ (t) -Ω q,1 (t) -Ω q,2 (t) -Ω q,+ (t) 0 -Ω q,2 (t) -Ω q,1 (t) -Ω q,1 (t) -Ω q,2 (t) 0 -Ω q,-(t) -Ω q,2 (t) -Ω q,1 (t) -Ω q,-(t) 0    , ( 72
)
with the following couplings

Ω q,± (t) = σ=a,b ūσ q,± (t) d dt vσ q,± (t) -vσ q,± (t) d dt ūσ q,± (t) (73) Ω q,1 (t) = σ=a,b vσ q,- (t)
d dt vσ q,+ (t) -ūσ q,- (t)
d dt ūσ q,+ (t), ( 74
) Ω q,2 (t) = σ=a,b ūσ q,- (t)
d dt vσ q,+ (t) -vσ q,- (t)
d dt ūσ q,+ (t), (75)
that are explicitly calculated in appendix C, in equations (C7-C9). In the general case, the amplitudes bq,± (t) and b-q,± (t) are coupled in the time dependent optical lattice. In particular one has:

i d dt      bq,+ (t) bq,-(t) b † -q,+ (t) b † -q,-(t)      = M(t)      bq,+ (t) bq,-(t) b † -q,+ (t) b † -q,-(t)      (76)
with the evolution matrix M(t)

M(t) =    ω q,+ -i Ω q,1 -i Ω q,+ -i Ω q,2 -i Ω q,1 ω q,--i Ω q,2 -i Ω q,- -i Ω q,+ -i Ω q,2 -ω q,+ -i Ω q,1 -i Ω q,2 -i Ω q,--i Ω q,1 -ω q,-   
(77) By using the symplectic symmetry of the matrix M, that is of the form

M(t) = A B -B * -A * we can write      bq,+ (t) bq,-(t) b † -q,+ (t) b † -q,-(t)      =    A q,+ (t) C q,+ (t) B q,+ (t) D q,+ (t)    bq,+ (0) +    C q,-(t) A q,-(t) D q,-(t) B q,-(t)    bq,-(0) +     B * q,+ (t) D * q,+ (t) A * q,+ (t) C * q,+ (t)     b † -q,+ (0) +     D * q,-(t) B * q,-(t) C * q,-(t) A * q,-(t)     b † -q,-(0) (78)
where the evolution equation for each column vector in the right hand side of equation ( 76) is governed by the matrix M(t)

i d dt    A q,+ (t) C q,+ (t) B q,+ (t) D q,+ (t)    = M(t)    A q,+ (t) C q,+ (t) B q,+ (t) D q,+ (t)    (79)
etc., with and initial conditions A q,± (0) = 1 and B q,± (0) = C q,± (0) = D q,± (0) = 0.

D. Two components adiabaticity parameter

General case

According to the expansion (76), and assuming that no excitations were initially present, the number of Bogoliubov quasi-particle in mode q created by the ramp is

n ex q (t) = b † q,+ (t) bq,+ (t) + b † q,-(t) bq,-(t) (80) = σ=± |B q,σ (t)| 2 + |D q,σ (t)| 2 . ( 81
)
In order to minimize the total amount of excitations during the time evolution, one needs to reduce the offdiagonal couplings in the evolution equation of the amplitudes (79). This leads to the adiabaticity conditions :

Ω q,+ (t) 2 ω q,+ (t) 1 ; Ω q,-(t) 2 ω q,-(t) 1 (82) Ω q,1 (t) ω q,+ (t) -ω q,-(t) 1 (83) Ω q,2 (t) ω q,+ (t) + ω q,-(t)
1 .

(84)

Equal masses

Let us now concentrate on the case of equal masses, corresponding to two internal states of the same atomic species. In this case, with ∆E q,a = ∆E q,b ≡ ∆E q , J a (t) = J b (t) ≡ J(t), the couplings Ω q,1,2 (t) vanish as U a (t), U b (t) and U ab (t) share the same time dependence. It follows that the Bogoliubov modes of the "+" and "-" branch evolve independently and we are left with the adiabaticity conditions in equation (82), with Ω q,± (t) = 1 2

d dt log ∆E q (t) ω q,± (t) . (85)
The adiabatic conditions are most stringent at the minimal quasi-momenta |q min | = q min (47), and one introduces the two adiabaticity parameters α ±

α ± = max 0≤t≤tramp α qmin,± (t), (86)
α qmin,± (t) ≡ 4∆E qmin d dt ∆E qmin ω qmin,± , (87)
corresponding to the two excitation branches in the twocomponent system. In order for the evolution to be adiabatic the conditions α ± 1 must be satisfied.

In the limit of a large particle number N → ∞, q min → 0 and single-particle unit filling n a + n b = 1

α qmin,± qmin→0 = 4J(t) d dt J(t) c ± (t)q min , (88)
where c ± (t) are the sound velocities of the two excitation branches defined as follows

c 2 ± (t) = 1 2 c 2 a (t) + c 2 b (t) ± (c 2 b (t) -c 2 a (t)) 2 + 4 U 2 ab (t) U a (t)U b (t) c 2 a (t)c 2 b (t) , (89)
c a (t) = d 2J(t)U a (t)n a , (90)
c b (t) = d 2J(t)U b (t)n b . (91)
IV. ADIABATIC TIME

In this section, for the case of equal masses corresponding to two hyperfine states of the same atomic species, we derive the expression of the adiabatic time for a linear ramp and we investigate the influence of a density imbalance between two components.

For a linear ramp (48), by changing variables from t to V 0 in (86), the adiabaticity condition is t ramp t adiab,± , where

t adiab,± = max Vmin≤V0≤Vmax t ± (V 0 ), (92)
t ± (V 0 ) = (V max -V min ) 4∆E qmin (V 0 ) d dV 0 ∆E qmin (V 0 ) ω qmin,± (V 0 ) . (93)
Let us consider a situation in which V min = 4E R , V max = 13E R and the fraction n a of atoms in component a is varied from 0 to 1, with n b = 1 -n a .

In Fig. 1 we plot the times t ± (V 0 ) in Eq. ( 93) as a function of V 0 /E R , for different atom numbers and fractions n a . Within the selected range of V 0 , the maximum of t ± is reached for V 0 = V max . The maximum of the "minus" branch (bottom row) is always larger than that of the "plus" branch (top row), hence setting the minimal time scale for adiabatic evolution.

In Fig. 2 we plot the adiabatic time t adiab,-given by Eq. (92) as a function of n a for different atom numbers (colored solid lines). The black dashed line is the approximation for large N :

t adiab,± qmin→0 = max Vmin≤V0≤Vmax t lin ± (V 0 ), (94)
t lin ± (V 0 ) = (V max -V min ) 4J(V 0 ) d dV 0 J(V 0 ) c ± (V 0 )q min , (95)

A. Discussion of the results

To gain physical insight in the low energy excitations involved in the adiabaticty condition for the two component system, let us first consider the homogeneous case without the lattice [START_REF] Castin | Les Houches LXXII:Coherent atomic matter waves[END_REF]27]. A clear physical picture is obtained by linearizing the coupled Gross-Pitaevskii equations

i ∂ t ψ a = - 2 ∆ 2m + U a |ψ a | 2 + U ab |ψ b | 2 ψ a , (96) i ∂ t ψ b = - 2 ∆ 2m + U b |ψ b | 2 + U ab |ψ a | 2 ψ b , (97)
for the fields of the two components 96)-(97) gives

ψ a (r, t) = N a φ a (r, t) ψ b (r, t) = N b φ b (r, t), (98)
i ∂ t δψ a = - 2 ∆ 2m δψ a + n a U a (δψ a + δψ * a ) + √ n a n b U ab (δψ b + δψ * b ), ( 99
) i ∂ t δψ b = - 2 ∆ 2m δψ b + n b U b (δψ b + δψ * b ) + √ n a n b U ab (δψ a + δψ * a ), (100)
One then looks for eigenmodes in the form of plane waves

δψ σ δψ * σ = q =0 ūσ q vσ q e i(q•r-ωqt) . ( 101
)
where σ = a, b. Note that in the density-phase point of view where ψ a = √ ρ a e iθa and ψ b = √ ρ b e iθ b we obtain

δρ σ q = 1 2 (ū σ q + vσ q ), (102)
δθ σ q = 1 2i ūσ q -vσ q √ n σ . (103)
The linearized equations ( 99)-(100) decouple when taking the sum and difference. One then finds that

ūa q -va q = ω q ∆E q (ū a q + va q ), (104)
and similarly for b, where ∆E q = 2 q 2 /2m, is the kinetic energy, and the sum combination satisfies the eigenvalue problem ∆E q 1 -ω q ∆E q 2 I + 2M ūa q + va q ūb q + vb q = 0, (105) with I the identity matrix and In each case the maximum of t± is reached for the largest value of V0/ER. The plots where generated expanding the kinetic energy term ∆Eq min for a low quasi-momentum dqmin 1 and using the approximate formulas (7a)-(7b) for J and U a,b,ab with scattering lengths aaa = a bb = 100.4r0 and a ab = 95.0r0 where r0 is the Bohr radius, corresponding to the hyperfine sates F = ±1, mF = ∓1 of 87 Rb atoms with a Feshbach-tuned interspecies scattering length a ab [6,28]. The lattice wavelengths is λ = 800nm.

M = n a U a √ n a n b U ab √ n a n b U ab n b U b . ( 106
)
The eigenvalues of M are

λ ± = 1 2 U a n a + U b n b ± (U a n a -U b n b ) 2 + 4n a n b U 2 ab , (107)
and the two eigenenergies take the form

ω q,± = ∆E q [∆E q + 2λ ± ]. (108)
Looking for the eigenvectors of the 2x2 matrix M, one finds the ratio between the density perturbations in the two components for the two excitation branches

δρ b ± δρ a ± = ūb q,± + vb q,± ūa q,± + va q,± = λ ± -U a n a U ab √ n a n b = sgn(λ ± -U a n a ) λ ± -U a n a λ ± -U b n b 1/2 . ( 109
)
In the limit of a small wave vector q, they correspond to sound waves with the speed of sound given by

c homo ± = λ ± m . (110)
The relations in this subsection, derived for a homogeneous system, hold in the lattice, provided one reinterprets the kinetic energy term ∆E q according to Eq. (20). In particular Eq. (108) for the spectrum coincides with Eq. (66), Eq. ( 109) can be easily deduced by the modal functions (69), and Eq. (89) for the speed of sound reduces to (110) by performing the substitution

1 m → 2J(t)d 2 2
that is the low energy equivalent of the correspondence (20) for kinetic energy. Provided U b n b = 0, we can rewrite in the lattice

c 2 ± = d 2 J(t)U b (t)n b 2 × F, (111)
F =   1 + U a n a U b n b ± 1 - U a n a U b n b 2 + 4 n a n b U 2 ab U 2 b   . (112)
Note that when U a (t) and U b (t) have the same dependence on V 0 , the factor F in (111) gets out of the derivative in (95), and gives the dependence of the adiabatic time on the ratios between the atom numbers and scattering lengths in the two components. For the linear ramp one indeed has 92): N = 10 3 (black dash-dotted line), N = 10 4 (red dashed line), N = 10 6 (blue solid line), while the black dashed line is from the approximate expression (94) that holds in the large N limit. The plots where generated expanding the kinetic energy term ∆Eq min for a low quasi-momentum dqmin 1 and using the approximate formulas (7a)-(7b) for J and U a,b,ab with scattering lengths aaa = a bb = 100.4r0 and a ab = 95.0r0 where r0 is the Bohr radius. The lattice wavelengths is λ = 800nm.

t lin ± = N 1/3 (V max -V min ) 8π √ n b F 1 J d dV 0 J U b , (113)
which also shows that, for a linearized dispersion relation, t adiab,-is always larger than t adiab,+ .

Symmetric case

For U a = U b and n a = n b = 1/2, the factor √ n b F in the denominator of the t lin ± (113) is equal to (1 ± U ab /U a ) 1/2 showing the divergence of t lin -when U ab → U a . Equation (109) shows that the "minus" solution corresponds to a situation in which the two component oscillate out-ofphase while the "plus" solution corresponds to in-phase oscillations

δρ a -= -δρ b -, (114)
δρ a + = δρ b + . (115)
The out-of phase solution, whose speed of sound c -tends to zero as U ab approaches U a from below, announces the demixing instability for U 2 ab > U a U b .

Asymmetric case

In the asymmetric case, with U a = U b but n a n b , to the lowest order in n a , that is zero order, one has ω q,-= ∆E q , (116) ω q,+ = ∆E q (∆E q + 2U b n b ).

(117)

The "minus" branch ω q,-does not describe a sound wave, but is purely quadratic for small q, see (9) [START_REF]the Hartree-Fock (HF) limit[END_REF].

Correspondingly, for n a = 0 (and symmetrically for n a = 1) Ω q,-= 0 in equation ( 85) and the adiabatic time is zero as shown in the right panel of figure 2. We note however that the adiabatic time significantly increases where 0 < n a n b . This feature can be simply explained in the large N limit where, using equations ( 112) and (113), one gets

(t lin -) na n b (t lin -) na=n b na→0 ∼ 1 2 1 + U ab U b 1 √ n a . (118)
V. CONCLUSIONS

We found the adiabaticity condition when an optical lattice is raised in a uniform two-components Bose-Einstein condensate at single particle unit filling. We concentrated our analysis in the superfluid regime and we used the time dependent Bogoliubov approach. We find that the excitations that can brake adiabaticity are pairs of Bogoliubov excitations (sound waves) with opposite quasi-momenta, where, in the symmetric case the two components oscillate out of phase. The scaling of the adiabatic time with the system size is N 1/3 , given by the minimal quasi momentum that can be excited in the lattice. For large atom numbers we show that the adiabatic time is significantly larger in the strongly asymmetric case n a n b with respect to the symmetric case n a = n b , the ratio scaling as 1/ √ n a .

VI. ACKNOWLEDGMENTS

This work was supported by the Polish National Science Center Grants DEC-2015/18/E/ST2/00760 and by the CNRS PICS-7403. D.K. acknowledges the financial support of the French Government (BGF) and the French Embassy in Poland.

One has s a q,± = ∆E q,a ω ± 1 1 + x 2 q,± (C3) d a q,± = ω ± ∆E q,a 1 1 + x 2 q,± (C4) s b q,± = x q,± ∆E q,b ω ± 1 1 + x 2 q,± (C5)

d b q,± = x q,± ω ± ∆E q,b 1 1 + x 2 q,± . ( C6
)
Using equations (C3)-(C6), and equation ( 73)-( 75) we obtain

2Ω q,± = 1 1 + x 2 q,± d dt log ∆E q,a ω q,± + x 2 q,± 1 + x 2 q,± d dt log ∆E q,b ω q,± (C7)
-2(Ω q,1 + Ω q,2 ) = ω q,+ ω q,- 1 (1 + x 2 q,-)(1 + x 2 q,+ )

× d dt log ∆E q,b ∆E q,a |x q,-| |x q,+ | (C8)
-2(Ω q,1 -Ω q,2 ) = ω q,- ω q,+ 1 (1 + x 2 q,-)(1 + x 2 q,+ ) × d dt log ∆E q,a ∆E q,b |x q,-| |x q,+ | . (C9)

  around the uniform solutionsψ a = √ n a , ψ b = √ n b with chemical potentials µ a = n a U a + n b U ab , µ b = n b U b + n a U ab ,where φ a = φ b = 1/ √ V are the stationary condensate wave functions, n a = N a /V , n b = N b /V are the uniform atomic densities, and U σ = 4π 2 aσ m for σ = a, b, ab are the interaction constants. Linearization of Eq. (

N=10 6 FIG. 1 .

 61 FIG. 1. (Color online) Times t+(V0) (top row) and t-(V0) (bottom row) in equation (93) as a function of V0/ER, for three atom numbers: N = 10 3 , N = 10 4 and N = 10 6 from the left to the right. In each panel we show the result for different fractions of atoms in component a: na = 0.5 (black line), na = 0.1 (red dashed line) and na = 0.01 (blue dash-dotted line).In each case the maximum of t± is reached for the largest value of V0/ER. The plots where generated expanding the kinetic energy term ∆Eq min for a low quasi-momentum dqmin 1 and using the approximate formulas (7a)-(7b) for J and U a,b,ab with scattering lengths aaa = a bb = 100.4r0 and a ab = 95.0r0 where r0 is the Bohr radius, corresponding to the hyperfine sates F = ±1, mF = ∓1 of 87 Rb atoms with a Feshbach-tuned interspecies scattering length a ab[6,28]. The lattice wavelengths is λ = 800nm.

3 FIG. 2 .

 32 FIG. 2. (Color online) Adiabatic time t adiab,-as a function of the fraction na of atoms in component a, for different atom numbers, in linear scale (left panel) or in semi-log scale (right panel). The solid lines are from Eq. (92): N = 10 3 (black dash-dotted line), N = 10 4 (red dashed line), N = 10 6 (blue solid line), while the black dashed line is from the approximate expression (94) that holds in the large N limit. The plots where generated expanding the kinetic energy term ∆Eq min for a low quasi-momentum dqmin 1 and using the approximate formulas (7a)-(7b) for J and U a,b,ab with scattering lengths aaa = a bb = 100.4r0 and a ab = 95.0r0 where r0 is the Bohr radius. The lattice wavelengths is λ = 800nm.

Appendix A: Implications of discrete translational symmetry

The single particle Hamiltonian ( 5) is invariant under discrete translation R along any lattice vector. It follows that translation operator TR commutes with the Hamiltonian i.e. [ ĥ(t), TR ] = 0, at any time. Complete set of commuting observables theorem implies that both operators share the same set of eigenstates. The eigenstates of the operator TR are labeled with discrete quasimomentum q values [30]

with eigenvalue e iq•R . Equation (A1) suggests the following form of the eigenstate

where function u q (r, t) has the same periodicity as the lattice. Thus, single-particle Hamiltonian eigenstates are of the form ψ n,q (r, t) = e iq•r u n,q (r, t), with additional band index n. This is a statement of the Bloch theorem.

In a process where only the height of the potential varies in time and the periodicity of the lattice stays constant the time-evolution operator Û (t) defined as

where T denotes time-ordering operator, also commutes with TR . Consequently, the time-evolving state ψ(r, t) = Û (t)ψ n,q (r) will remain an eigenstate of the translation operator with the same eigenvalue as the Bloch state ψ n,q (r). In other words, quasi-momentum is conserved during time evolution. It follows from TR Û (t)ψ n,q (r) = Û (t) TR ψ n,q (r) = e iq•R Û (t)ψ n,q (r). (A4)

Although the quasi-momentum needs to be conserved, it does not exclude situation where other bands may be populated. In general one can write

Another consequence of translational invariance is the vanishing of the coupling term

when k = q. This can be shown once we consider the identity

and quasi-momentum conservation TR ∂ ĥ(t) ∂t ψ m,q (r, t) = ∂ ĥ(t) ∂t TR ψ m,q (r, t)

= e iq•R ∂ ĥ(t) ∂t ψ m,q (r, t). (A8)

Eq. (A8) shows that

This result combined with orthogonality of Bloch states gives a non-zero coupling term only if the quasimomentum of both states is the same.

Appendix B: Explicit formulas for the two-component system case

The non-hermitian matrix L is defined as follows

The energy introduced in Eq. (61) takes the form

Appendix C: Couplings among Bogoliubov modes

It is convenient to introduce the sum and differences of the Bogoliubov modal functions (69)

with σ = a, b, and express them using the parameters (C12) .

[30] The discrete nature of quasi-momentum is a consequence of Born-von Kármán boundary conditions in the finite size system.


(1) L'Institut National de l'Environnement Industriel et des Risques INERIS, Nancy, France. Email: Jannes-L.KINSCHER@ineris.fr, Tel. +33 3 55 66 28 81 (2) Institut de Physique du Globe de Paris IPGP, Sorbonne Paris Cité, UMR 7154 CNRS, Paris, France (3) University Paris Diderot, Paris and ANGE Team, INRIA, CETMEF, J.-Louis Lions, Paris, France (4) GFZ German Research Centre for Geosciences, Potsdam, Germany (5) Université de Lorraine, GeoRessources lab.,campus Artem, Nancy, F54042, France

In order to improve our understanding of hazardous underground cavities, the development and collapse of a ~200 m wide solution mining cavity was seismically monitored at a SOLVAY exploitation site at Cerville-Buissoncourt in the Lorraine basin, NE France (e.g. Kinscher et al. 2015). The recorded microseismic data set (~ 50.000 event files) represents two major microseismic crises, (i) one recorded in March to April 2008, related to a period of significant cavity roof gowth (without subsequent surface collapse), and (ii) one crisis in February 2009, associated with a 3 days lasting collapsing period, induced by controlled brine pumping operations. The numerous microseismic events (~ 50.000 event files) appeared as unusual swarming sequences with complex clusters of superimposed microseismic events (Fig. 1). Consequently, only a minor portion (~20%) of the data set could be processed through standard body wave phase based detection, location and source analysis routines. In this study, we present innovative methods able to generate a complete seismic catalogue of this unique microseismic data. First, a probabilistic, automatic detection location design has been developed that constrain the spatio-temporal characteristics of microseismic swarms in a first order (Kinscher et al. 2015). In this approach, mainly signal polarization attributes were used to identify P wave phases. Event locations were then obtained by using the information of the corresponding polarization angles, as well as signal amplitudes and a calibrated attenuation law describing the amplitude decay with distance. Moreover, seismic source mechanisms and standard source parameters were investigated by comparing observed and synthetic wave forms, amplitude spectra and peak-to-peak amplitude ratios, assuming different standard source models. Within this source study, it was possible to assess the common source mechanisms of many swarming events at once, using the observation from only one single three component station. The main results obtained by these developed methods, as well as by detailed statistical analysis of the spatio-temporal and energetic microseismic attributes are as follows: (i) the individual microseismic events are comparable to small natural tectonic earthquakes with moment magnitudes Mw ranging from around -3 to 1; (ii) source mechanisms for most microseismic events are remarkable stable and demonstrate a predominant thrust faulting (double-couple) regime with faults similarly oriented NW-SE, dipping 35° -55° , what might be related to the presence of systematically arranged pre-existing fractures; (iii) the origin of microseismic swarming is suggested in the incapacity to sustain larger strains and to release larger stresses, what seems to be related to the mechanical constitution of the rock strata overlying the cavity (i.e. low strength materials); (iv) caving and collapsing periods at the cavity roof are associated with systematic, self-reinforcing dynamics and have a distinct microseismic response, clearly observable from statistical analysis, which can be precisely described by empirical laws. These results underline that microseismic monitoring and the here developed methods are a powerful tool to constrain the mechanical and dynamical characteristics of an evolving and collapsing hazardous underground cavity. In addition, our results provide new important constraints for the design of adequate geomechanical models that will improve hazard assessment at underground openings.

Figure 1 :

 1 Figure 1: Example of a typical swarming sequence. Figure is taken from Kinscher et al. (2015).