Daniel Comparat

A Study of molecular cooling via Sisyphus processes

Keywords: numbers: 37.10.Mn, 33.20.-t, 34.20.-b

We present a study of Sisyphus cooling of molecules: the scattering of a single-photon remove a substantial amount of the molecular kinetic energy and an optical pumping step allow to repeat the process. A review of the produced cold molecules so far indicates that the method can be implemented for most of them, making it a promising method able to produce a large sample of molecules at sub-mK temperature. Considerations of the required experimental parameters, for instance the laser power and linewidth or the trap anisotropy and dimensionality, are given. Rate equations, as well as scattering and dipolar forces, are solved using Kinetic Monte Carlo methods for several lasers and several levels. For NH molecules, such detailed simulation predicts a 1000-fold temperature reduction and an increase of the phase space density by a factor of 10 7 . Even in the case of molecules with both low Franck-Condon coefficients and a non-closed pumping scheme, 60% of trapped molecules can be cooled from 100 mK to sub-mK temperature in few seconds. Additionally, these methods can be applied to continuously decelerate and cool a molecular beam

In this article we study the Sisyphus method for cooling molecules. In such cooling, sketched in Fig. 1, external forces remove kinetic energy by transferring it into potential energy. A (absorption-)spontaneous emission step follows, creating non reversibility of the process which can be repeated by optical pumping the molecule back to its original state. The method is first compared to other ones, then illustrated with a simple 1D or 3D model, and finally detailed on specific cases, such as the NH molecule. We then conclude that the method is very versatile and is able to produce large samples of molecules at very low temperature.

FIG. 1: (Color online). Principle of Sisyphus cooling of molecules: 1) Removal of kinetic energy through motion in an external potential, 2) "Sisyphus transfer step": a dissipative process avoids the reverse motion. 3) This "one-way" (or "single photon") process can be repeated by bringing the molecule back in position using the trapping potential and (4) in its original internal state using light absorption in a "repumping step" (symbolized by the shaped spectrum for amplitude selection). r designs the radial coordinate. In 2D or 3D, due to the angular momentum the particles can miss the center, therefore rrepump could be non zero.

I. INTRODUCTION

A first simple idea to cool molecules is the use of evaporative cooling or collisions with dense and colder species, such as trapped laser cooled atoms or ions, in a so called sympathetic cooling scheme. Such a thermalisation technique has been demonstrated with molecular ions [1]. For a long time, reactive or inelastic collisions have strongly limited the efficiency of the process for neutral species [2,3]. Only very recently OH radicals have been cooled using evaporative cooling [4].

A second, straightforward idea is the laser cooling technique using transfer of photon momentum. While the first demonstration of light pressure on molecules occured in 1979 [5], it is only very recently that laser cooling of molecules has been observed [6]. To circumvent the general modification of the internal state occurring after the spontaneous emission step, this has indeed required choosing a very well suited molecule (SrF), which has a quasi-closed-level system with a very high Franck-Condon factor (see also Table I).

A third route, suitable for a larger number of molecular system, is to modify the kinetic energy into potential energy using external forces [7]. A so-called one-way [8][9][10][11], or irreversible, cooling has to be realized in order to avoid the reverse process which otherwise would heat the sample. This can be realized, for instance, by using an (absorption-)spontaneous emission step [12,13]. Finally, this can be repeated, in a so-called Sisyphus cooling process [14] (see Fig. 1). Compared to standard laser cooling based on photon momentum transfer, the reduction in temperature per spontaneous emission step is typically a few mK compared to a few µK. It can be a large fraction of the initial temperature. The process is therefore less sensitive to the modification of the internal state occurring after spontaneous emission. The Sisyphus cooling process, first proposed by Pritchard [15], has represented a milestone in the history of laser cooling by being responsible, through polarization induced light shifts, for breaking the Doppler limit in magneto-optical traps (MOT) [14]. The irreversible Sisyphus cooling cycle has then been realized for atoms [16] with the help of a gravity sag, both in magnetic [17] and in evanescent wave reflection fields [18], and with RF-induced transition both in an optical dipole trap [19] and in a magnetic trap [20]. A typical experiment, as reported in Ref. [19], concluded that "the final temperature achieved was (a somewhat disappointing) 17 µK". This explains why nowadays Sisyphus cooling is not heavily used by the cold atom community. But, such a low temperature would be a tremendous result for molecules. Indeed, a review, given by table II, of the produced cold molecular species indicates that most of them are produced at temperature of hundreds of mK. The Sisyphus effect has recently allowed the first efficient laser cooling of a polyatomic molecule from 400mK down to 30 mK with the phase-space density increased by a factor of 30 [21]. In such molecular Sisyphus cooling scheme some specific symmetric-top rotors, remain electrically trapped, and microwave transition between vibrational states provides the energy transfer [22]. Even if limited to a specific type of molecule, this extremely promising experiment already indicates the possibility to reduce the temperature below the mK range in tens or hundreds of seconds. One of the difficulties is to find molecules with short enough spontaneous emission times. A similar problem is found in Ref. [23], where it is suggested to put an optical cavity around a magnetically trapped OH sample in order to accelerate the spontaneous decay.

In this article, we suggest to use optical pumping and faster electronic transitions in order to generalize the method. We first discuss the required ingredients needed for efficient cooling such as the choice of the molecule, the trap and laser parameters. We then perform a detailed simulation of the whole process in realistic conditions in the case of NH molecule.

II. COOLING STRATEGY A. Choice of the molecular states

Each molecule has is own characteristics and the cooling strategy should be adapted to each of them by making choices for the trap, the electronic states and the lasers.

In Table I, we have listed some diatomic molecules [197] to illustrate their properties and explain why we think that the method can be implemented for most of the up to now produced cold molecular species. For simplicity, in this article, we do not consider any hyperfine structure, the splitting of which being typically in the MHz range and if needed can be resolved and manipulated by laser sidebands [6]. We simply note that the role of the hyperfine structure can be very important and, contrary to common thought, it can make the cooling simpler; for instance by creating the two trapping poten-tials U 1 , U 2 when without the hyperfine structure only a single one would exist.

For the choice of the lasers or trap, a key parameter is the AC or DC electric, magnetic or electromagnetic trapping capability of the molecule.

Electrostatic or magnetic trap can be used if the molecule present a strong enough energy shift ∆ due to Stark or Zeeman effect. In an external electric E or magnetic B field, the energy shift of the 2S+1 |Λ| Ω state depends on the rotational level J. For Hund's case a electronic state, it is respectively EdΩM/J(J + 1) and Bµ B (Λ -2Ω)ΩM/J(J + 1) where µ B is the Bohr magneton, d is the permanent electric dipole moment and M the projection of J/ along the quantization axis given by the local field at the particle position. Most of the molecules listed in Table I can thus be trapped in electric or magnetic trap. However some of them (with 1 Σ ground state for instance) cannot be electrically or magnetically trapped. In such case, one solution is to transfer the molecule to a state with sufficiently long lifetime where the method can then be applied [24,25]. Another solution, is to use a laser dipole trap [16,26]. But, because the initial temperature of the molecule is usually high (cf table II) a quite deep trap is required and so a close to resonance laser has to be used. For instance a 100 W laser focused on 1 mm diameter, detuned a fraction of a nm from a resonance leads to a typical trap depth of more than 10 mK and an oscillation frequency of ω ≈ 2π(500 Hz). Such laser can be realized through (fiber) laser amplification or using a Quasi-Continuous Wave laser [27]. Interestingly enough, the fact that the laser is close to resonance can be used to take benefit of the high diffusion rate (10 5 photons/s in the previous example) to perform the Sisyphus transfer between the two trapped state [16].

B. Removal elementary step strategy

Once the trap chosen, the cooling strategy is strongly linked to the energy removal elementary steps: step 1 and 2 in Figure 1.

In order to study the dynamics, we shall assume N particles initially at temperature T in potential U 1 (r). They are then transferred to the potential U 2 (r) and then, when reaching r repump = 0 (in the 1D case), are pumped back to U 1 . For simplicity, we shall illustrate first the strategies in the simple 1D picture and using two harmonic potentials U 1 (r) = mω 2 1 r 2 /2 and U 2 (r) = mω 2 2 r 2 /2 (where gravity is neglected). In such process there are no collisions to thermalize the sample, but for simplicity we shall use an effective temperature T (given by the average kinetic energy).

The Sisyphus transfer between state 1 and 2 can be obtained is several ways and using different strategies illustrated by 2 such as by Spontaneous or induced transfer that we want to discuss briefly:

• Spontaneous emission (figure 1)

Molecule Low state High state λ (nm) FC lifetime NH X 3 Σ - X 3 Σ - 3050 1 37 ms X 3 Σ - A 3 Π 336 1 400 ns X 3 Σ - b 1 Σ + 471 1 18 ms X 3 Σ - a 1 ∆ 794 1 12 s BH X 1 Σ + A 1 Π 433 1 160 ns YO X 2 Σ + A 2 Π 614 0.99 CH X 2 Π A 2 ∆ 431 0.99 530 ns SrF X 2 Σ + A 2 Π 651 0.98 23 ns CaH X 2 Σ + A 2 Π 693 0.98 CaF X 2 Σ + A 2 Π 606 0.98 20 ns YbF X 2 Σ + A 2 Π 552 0.95 BaF X 2 Σ + A 2 Π 860 0.95 50 ns MgH X 2 Σ A 2 Π 519 0.94 YbF X 2 Σ + A 2 Π 552 0.93 OH X 2 Π A 2 Σ + 308 0.91 700 ns ThO X 1 Σ + A 1 Σ + 943 0.86 VO X 4 Σ - B 4 Π 787 0.80 MnH X 7 Σ A 7 Π 556 0.78 CS X 1 Σ + A 1 Π 258 0.78 200 ns TiO X 3 ∆ A 3 Φ 709 0.72 CrH X 6 Σ + A 6 Σ + 866 0.68 CaS X 1 Σ + A 1 Σ + 658 0.59 CN X 2 Σ + A 2 Π 1097 0.50 10 µs CaO X 1 Σ + A 1 Σ + 866 0.43 X 1 Σ + A' 1 Π 1199 0 CO X 1 Σ + a 3 Π 206 0.31 10 ms X 1 Σ + A 1 Π 154 0.12 10 ns a 3 Π a' 3 Σ + 1453 0.04 3 µs a 3 Π d 3 ∆ 823 0.02 SrO X 1 Σ + A 1 Π 920 0.29 X 1 Σ + A' 1 Π 1074 0 SO X 3 Σ - A 3 Π 262 0.22 12 ns NO X 2 Π A 2 Σ + 227 0.17 200 ns X 2 Π a 4 Π 263 ? 100 ms SiO X 1 Σ + A 1 Π 235 0.15 9 ns N2 X 1 Σ + g a 1 Π 0.04 A 3 Σ + g B 3 Πg 8 µs PbO X 1 Σ + a 3 Σ + 629 0.02 H2 X 1 Σg B 1 Σ + g 111 0.01 1 ns c 3 Πg j 3 ∆g 567 0.99 LiH X 1 Σ + A 1 Σ + 385 0 X 1 Σ + B 1 Π 291 0.08 O2 X 3 Σ - g B 3 ∆ - g 286 0 predissociation
TABLE I: Some properties of some diatomic molecules. In italic are the molecules that have not been cooled up to now (cf. table II). The molecules are ordered by the rounded value of the v ′ = 0, v ′′ = 0 Franck-Condon factor. The main parameters needed to realize the Sisyphus cooling is a trapping capability (depending on the electronic state), a transition wavelength in the region cover by laser (typically near the visible region) and a good Franck-Condon ratio in order to simplify the repumping step. The spontaneous emission (lifetime) is mainly here to ensure the feasibility of the process.

As shown by the figure 1, the simplest solution is simply to wait for spontaneous emission. This can be used when the first state (U 1 ) has a long enough lifetime so that the particle moves in the trap before its spontaneous decay [28][29][30].

Unfortunately, in many cases U 1 and U 2 will simply be two rotational or vibrational states of the electronical ground state having a very long spontaneous decay time. For such cases, another way is to control the Sisyphus step by applying an excitation laser, or any other light source (RF-microwave, ...), to the molecule initially in the potential U 1 to excited it to an excited state that will then decay towards U 2 .

• Continuous transfer (right and upper part of Figure 2): A simple way to mimic this case of a spontaneous emission with rate γ is to realize the Sisyphus transfer step by using a large and broadband laser covering the whole sample. This creates a quite uniform excitation rate γ toward an excited state that (quickly) decay toward U 2 . Usually γ is chosen to be smaller or on the same order of magnitude than the trap oscillation frequencies in order to avoid too fast transfer that would transfer the molecule before it has lost enough energy.

In such process (spontaneous emission or uniform excitation), in average, the particle loose kinetic energy simply by the fact that a slowing down molecule spends more time at the spacial points where the velocity is reduced. We can draw here some similarity with the optical pumping scheme in a blue optical molasses where atoms scatter less photons when they have less kinetic energy because there is less light present [31,32]. Hence, under the effect of its motion the molecule will be transferred toward U 2 when it has the smallest kinetic energy. More precisely, neglecting the photon recoil energies, and using survival probability described in the appendix B 4, we see that the average energy reduction is given by +∞ 0

(U 1 (r(t)) -U 2 (r(t)))γe -γt dt. Thus, in our simple harmonic trap example, using the typical trajectory r(t) for the average kinetic energy k B T , this energy reduction is [26]. This value can be easily improved by using a non uniform excitation rate γ(r) as graphically illustrated by the case 3 in Fig. 2. Moreover, we can always combine these ideas and, depending on the oscillation frequencies, laser power and cooling time available, find the optimum spatial distribution of the light power which optimizes the one way cooling.

k B ∆T = 0.4k B T (1 -ω 2 2 /ω 2 1 ) for γ = ω 2
• Induced spatial transfer (left and upper part of Fig-

ure 2):

When using a laser induced transfer, one idea is not to cover the whole sample by laser, but to catch the particle at the "top" of its trajectory that is when it has almost zero kinetic energy. This selection can be realized spatially by applying the laser only at this position. This simple realization is illustrated by the localized arrow in the left and upper part of Figure 2.

• Induced energy transfer (cf. middle and upper part of Figure 2): However, it can be sometimes simpler to apply the laser on the overall sample and to realize this selection spectroscopically. We can use the fact that the resonance condition is fulfilled only locally due to different potentials seen by the particles during the excitation (case 2 of Figure 2). For this last two cases, we can assume that the laser transition occurs only at location r = r res . Due to Boltzmann statistic the number of molecules reaching the transition point are N e -η , where η = U 1 (r res )/k B T . Assuming that the particles hardly move before being optically (re-)pumped, i.e. that all (re)pumping rate are faster than the trap frequencies, the total energy change for the sample is then

∆E = N k∆T = N e -η (U 2 (r res ) -U 1 (r res )) < 0.
The molecule has spent at most a quarter of the oscillation time in each potential before being transfer to the other one. Thus, a maximum cycling time is ∆t = π/2ω 1 + π/2ω 2 . This lead to the equation

∆T ∆t = - 2ω 1 ηe -η π 1 -ω 2 2 /ω 2 1 1 + ω 1 /ω 2 T (1)
As for evaporative cooling processes [33,34], we see that it is probably efficient to adjust the laser detuning (or the external trapping potential) as the particles are losing energy. This can be done, for instance, by keeping the ratio η of the transition energy to the thermal energy constant, with η = 1 being the best choice. The speed can be optimized by choosing ω 2 = ω 1 /2 with a temperature variation per cycle ∆T ≈ 0.28T and a (1/e) decay time of 17/ω 1 .

Independently of the chosen strategy we can drive the conclusion that a 1D Sisyphus cooling will provide a temperature variation per cycle on the order of ∆T ≈ 0.3T . Then the temperature should decays exponentially with a (1/e) decay time on the order of 20/ω 1 . Thus, cooling by a factor 10000, say from 100 mK to 10 µK, takes only 30 cycles. This has to be compared to the few thousands of cycles needed for standard laser cooling [6]. This simple study clearly highlights the potential of the Sisyphus method.

C. 1D, 2D or 3D cooling?

In real the particles moved in three dimension (3D) for instance being trapped in an harmonic trap with

U (r) = 1 2 mω 2 x x 2 + 1 2 mω 2 y y 2 + 1 2 mω 2 z z 2 .
The problem is that, even in such simple case, where the motion in one of the three axis x,y,z is decoupled from the others, the previous discussion, concerning 1D motion, cannot be generalized to the (2D or) 3D case. The reason is that, if Fig. 1 is still valid, where r designing the radial coordinate, now in 2D or 3D motion, due to the angular momentum, the particles can miss the center and r repump could be non zero. This is clearly indicated by the trajectories, calculated in the case of an harmonic 2D potential and shown in the lower part of Fig. 2. The dashed red line trajectory, corresponding to an isotropic trap case ω x = ω y , performs at constant r. Thus the process given in Fig. 1 does not work at all because r does not change, so the trajectory never crosses the desired transfer point (box in case 1 and solid circle in case 2 of Fig. 2) where the kinetic energy should be zero, simply because the kinetic energy is constant and never zero.

There is at least two ways to go around this important difficulty:

• Create asymmetries in the potential In order to avoid particles that orbit around the trap center with constant r for instance, it is possible to create asymmetries in the potential [17,23,28,35]. This is illustrated by the sec-ond trajectory (solid line in the lower part of Fig. 2), corresponding to an anisotropic trap case ω x = 0.8ω y . The trajectory is more complex now and present some points where the kinetic energy is almost zero and thus the trajectory crosses the line or box where the transfer should occur and the cooling will be efficient.

An important message is that a small modification of the trap can totally modify the efficiency of the cooling (see also Fig. 9).

• 3 × 1D cooling If we only look on one axis (let say the x axis), the 1D Sisyphus cooling picture given previously is perfectly valid and efficient with r = x. Indeed, the motion along x crosses the zero center and, at the edge of the trajectory, the velocity along x is zero as desired for an efficient Sisyphus transfer. The usual 1D process can thus be used to reduce first the velocity, or the temperature T x , along this axis. We then just have to repeat the process for the other axis y and z. Experimentally, the cooling along x can be done by using a spatial selection that affects only the x coordinate but not the other ones. For this we can choose the induced spatial transfer (case 1) strategy and send the laser beam only along a line (or a plane in a 3D vision) as shown on the left lower part of Fig. 2 by the hatched area.

D. Phase space density and optimized strategy

The temperature is not the only relevant parameter. Another important parameter, which is to be maximized, is the phase space density D = n 0 Λ 3 , where n 0 is the peak density and Λ = 2π 2 mkB T is the thermal de Broglie wavelength [34]. In the 3D isotropic harmonic case, for

U (r) = 1 2 mω 2 r 2 , we have D = N ω kB T 3 and n 0 = N (2πσ 2 r ) 3/2 where U (σ r ) = 1
2 k B T . A proper optimization of the cooling should be done by clearly defining the objective function to be optimized. Depending of the goal, we might thus prefer to optimize the time, the number of particles, the temperature, the phase space density or, as done experimentally for evaporative cooling strategies, the parameter Ḋ/D Ṅ /N [34,36,37]. The strategy has also to be chosen depending of constrains such as: vacuum limitation, available lasers, trapping possibilities, ... We have seen that the efficiency of the process strongly depends on its dimensionality, on the trap anisotropy and on the transfer strategy chosen. Comparison of all possible processes and dependance with trap or laser parameters is beyond the scope of this article. However, we have tried to initiate such discussion in the Appendix A and the results are summarized by the Fig. 10 which gives the time evolution of number of molecules N , Temperature T and phase space density D under different cooling strategies. These results should be taken only as rule of thumbs and not as accurate predictions. In fact, one important consideration is that the interaction time, the power broadening of the transition or the velocity dependance of the transition rate can lead to a much lower energy or spatial resolution than naively expected. We thus need to confirm these results by simulations. But, before doing so we first extract, from the appendix, few considerations for the efficiency for three different cooling strategies:

• one-way or single photon cooling.

In the first step (1 and 2) of Fig. 1, a single photon cooling decreases the temperature but the density stay the same because we simply remove the kinetic energy. However, it is possible to increase furthermore the phase space density by transferring the particles in a tighter trap without heating. This can be done by catching them in a tight trap as indicated by the dashed trap shown in the upper panel of case 1 of Fig. 2). In order to catch all particles it is even possible to have this small, tight trap U 2 (the black box in the lower panel of the figure) that dynamically follows the location of the transfer. In principle such single photon cooling processes looks extremely attractive [38]: starting with the more energetic molecules, and slowly sweeping the laser position (for case 1) or the transition frequency (for case 2) down, we can always transfer the molecules when they have almost zero kinetic energy. A simple optimization of the time evolution (trajectory) of this box is proposed in Ref. [39]. If the lower potential U 2 is almost flat we can also in principle remove all the energy with a single photon emitted per molecule. Several realizations of this one-way or single photon cooling already exist, but only for atoms [9,10,12,13,35,40].

These one-way or single photon cooling schemes has to be used if the repumping step 4) of Fig. 1 is difficult to realize. An obvious example is when the decay after step 2) populate a lot of levels, for instance due to bad Frank-Condon factors, because it will be quite difficult to optically pump all of them.

• 3 × 1D cooling The 3 × 1D cooling is very efficient because it is simply three times the 1D cooling one. For instance the time needed to reach the same final temperature has to be multiplied by a factor 3 compared to the 1D case. It is thus a very fast cooling and it works with any kind of trap geometry. Another advantage of such 1D cooling, compared to the 3D one, is that the quantification axis is always the same (there is only one possible axis!) and transitions can thus be controlled using light polariza-tion, which can be extremely useful to excite only the states which are trapped (for instance by controlling the Zeeman sub-level in a magnetic trap). A last advantage is that the particles always cross the 0 coordinate along the considered axis and so the repumping light can then be focused only at the center of the trap in order to realize the ideal step 4 of Fig. 1. Finally we would like to mention the very useful rotation of phase space method used in Ref. [17] in order to cool always along the same direction. Quoting the authors: "First the vertical distribution is cooled. Then, by adiabatically changing the trap parameters, we produce a degeneracy ω x ∼ ω z . Anharmonicities cause the atom cloud to rotate in the x-z plane, exchanging the x and z distribution. After we adiabatically restore the trap to the initial field conditions. We then cool the vertical distribution again." Therefore only one laser direction is enough to cool the 3D sample!

The conclusion is that this three times 1D cooling may be simpler and more efficient that direct 3D cooling [17,22] but its practical implementation has to be done with care since the other axis should be spectator, not only in the cooling step but also in the repumping step.

• Standard 3D Sisyphus cooling scheme.

In the following we have chosen to study in more detail the Sisyphus cooling scheme because it keeps all molecules and does not require to move lasers. But , even without a detail simulation, we can directly see from Fig. 1 that, as the particle returned to the original potential (step 3), the reduction of temperature also induces an increase of the density and of the phase space density.

III. SISYPHUS COOLING ON NH: CASE WITH HIGH FRANCK CONDON VALUE

F C ∼ 1.
To catch more insight of the Sisyphus cooling process, we shall first study a case where the repumping step is simple. In order to have a fast cooling process we want to use an electronic transition. The simplest choice it to choose molecule with a Franck Condon factor close to unity creating a closed vibrational pumping scheme. With such close cycling scheme the Sisyphus cooling process is clearly advantageous compare to the one-way single-photon one because we can repeat the single photon process without losses.

As an example we shall investigate the NH (Imidogen) molecule. The A-X transition in NH is indeed almost perfectly diagonal and the v' = 0 -v" = 0 band has a Franck-Condon factor of better than 0.999 [40]. NH is a very well known molecule which has already been decelerated and cooled (see table II). It has been trapped both in electric [41] and magnetic traps, for more than 20 s [42], in the a 1 ∆ or its X 3 Σ -state respectively. The accumulation of Stark-decelerated NH molecules in a magnetic trap has been realized [40].

NH is probably one of the simplest systems to study cooling of molecule. A single photon cooling scheme has been envisioned [12]. For simplicity we shall concentrate here on NH in its X state trapped in magnetic field. But the case of the a 1 ∆ state and an electrostatic trap is also possible with the 325 nm c 1 Π ←a 1 ∆ Sisyphus transition and repumping on the 452 nm c 1 Π ←b 1 Σ + transistion.

A 471 nm b 1 Σ + ←X 3 Σ -transition is possible but a simple and efficient transition is the

A 3 Π; v ′ = 0 ←X 3 Σ -; v ′′ = 0 transition.
In order to perform the Sisyphus cooling we need to find two states with different magnetic moment and then couple them using laser light through an intermediate level. We choose not to use the absolute rotational ground state but the X 3 Σ -; N ′′ = 1, J ′′ = 1+. This allows to perform the Sisyphus cooling on the A 3 Π; N ′ = 1, J ′ = 0, -←X 3 Σ -; N ′′ = 1, J ′′ = 1+ transition that is, due to parity reasons, the only closed N ′ ← N ′′ rotational transition. A sketch of the relevant levels and transition strength, is shown in Figure 3 as calculated by a Program for Simulating Rotational Structure (PGO-PHER) [198] with the molecular constants taken from Ref. [43]. From Fig. 3, we can find two states with different magnetic moments to perform the Sisyphus cooling: U 2 being the potential of the M ′′ = 1 level (with a linear Zeeman shift) and U 1 being the M ′′ = 0 one (with a quadratic Zeeman shift). One drawback is the possible decay towards the untrapped M ′′ = -1 level, but, because of the fast electronic transition, a fast repumping from this state seems nevertheless possible.

Even if it is not the purpose of this article, we mention that using this transition standard laser cooling of NH should be quite easy. The fact that the transition is a J ′′ → J ′ = J ′′ -1 scheme and not the standard J → J + 1 one for magneto-optical trap (MOT), may lead to some difficulties because in 1D this produces a dark state that prevents the cooling [6,44]. However the dark state is not dark for the other lasers present in a 3D MOT. We have simulated such 3D cooling and found results comparable to the experimental demonstration of the so called 3D type-II magneto-optical trap on sodium atom [45]. However, as stressed in the Appendix, rate equations are not able to deal correctly with dark states so the simulation of a J ′′ → J ′ = J ′′ -1 transition using rate equation may be questionable.

A. Cooling in quadrupolar trap with energy degeneracy

For the Sisyphus cooling, we first study the simple case of a quadrupolar trap. We explain in detail our choice for the parameters in a way that can be generalized for other situations and then present the results.

The magnetic field is thus

B(R) = B ′ 0 (1 + α)Xe X + B ′ 0 (1 -α)Y e Y ) -2B ′ 0 Ze Z with B ′ 0 = 100 T/m
and α is a parameter that controls the anisotropicity of the trap. We choose the anisotropy parameter α = 0.17 (see figure 9) which optimize, after 20 ms, both the energy conversion for the pumping and repumping steps. This situation can easily be experimentally realized, for example using 2 orthogonal pairs of coils in anti-Helmholtz configuration. In a gradient of ∼ 100 T/m a typical molecule with potential energy of ∼ 150 mK (corresponds to 0.1 cm -1 ), sees a magnetic field slightly less than 0.2 T (see Fig. 3). This leads to a typical radius for the sample of σ r ≈ 2 mm. In order to affect most (say 90%) of the molecules, we choose an initial laser waist size of 5 mm and we start the Sisyphus process with a laser detuning corresponding to a transition with a potential energy of four times the energy: 500 mK (in temperature units).

The time evolution of the cooling is linked to two parameters: the typical "oscillation" time T 1 (2 ms in our case) and the number of "oscillations" required before the molecule is transferred. As indicated in the appendix, this number N d-1 0 (1 for a 1D, N 0 for a 2D, and N 2 0 for a 3D cooling) is linked to the selectivity ∆E in energy (or in spatial position) of the transition by ∆E ≈ E/N 0 in the 3D case. In our case, the study of the trajectories in the trap (see energies in figure 9) indicates that an accuracy on the order of ∆E = 20-70 mK is probably a good choice for a transfer after 20 ms. This suggests to take a value N 0 < 5 for our case of an initial average potential energy of ∼ 150 mK (temperature of 100 mK). Thus, we expect, a typical decay time for the temperature on the order of N 2 0 T 1 , i.e. below 50 ms. We thus choose in our simulation to reduce all lasers powers, waists and detunings (i.e. transition energy) exponentially with a time constant of τ = 50 ms.

The choice of the laser power is discussed in the Appendix B 3 e and the results is that the intensity I = 2P πw 2 should be on the order of kBT Γ v σr ΓL Γ (500nm) 3 λ 3 1 300 000 FC.AF . For transition at λ = 336 nm with a linewidth Γ = 1/(400ns) (see table I), a Franck Condon factor FC=1, an angular factor AF = 1/4 (ratio of the worst transition to the total one in Fig. 3), and a velocity v at the transition taken to be typicality (for reasonable N 0 value) tens of the thermal velocity; a standard Γ L = 2π(10M Hz) linewidth laser requires a power of nearly 100 mW. We choose such values for the simulation.

For this first study, but this will be refined later, we do not put any realistic repumping laser but we choose a constant and uniform repumping rate of 10 6 s -1 for all molecules having a potential energy bellow 15 mK as indicated by the green box in Fig. 3.

For the simulation, the absorption-emission processes are accurately calculated using rate equations taking into account the saturation, the laser detuning its linewidth and the Doppler effects [46]. The simulation used a Kinetic Monte Carlo algorithm, to solve exactly the rate equations, and a Verlet algorithm to drive the particles motion under the effects of gravity, magnetic field and recoil photons [47]. The molecular dipole moment is as- We choose 100 molecules and a 100 mW power laser propagating along Oz and circular left polarized (σ + ). The results of the cooling are presented in Figure 4. We first see that the repumping is efficient and that, despite the transient population of an untrapped potential curve, we do not loose any molecules during the process. The other important result is that the temperature is reduced, and therefore the phase-space density increases, but that the process stops, and even gets worse after 100 ms. If, by choosing better laser parameters, the final temperature can be made lower, we would like to stress that such limited cooling behaviour is quite general and that choosing better cooling strategies or parameters is not so obvious. Thus, this is an important result appearing in several situations that we have simulated when dealing with zero field trap at the center. By looking to the transitions radius and potential energy for the particles, we found that with degenerate energy levels near the trap centre (i.e. at low temperature), the spectral or spatial selectivity for the transitions is no longer ensured and that we excite undesirable levels during the process.

B. Cooling in trap without energy degeneracy

A simple way to avoid such nasty behaviour is to use a trap without energy degeneracy, that is with a non zero field at the center, where the spectral selectivity for the transitions is ensured. Another major advantage of such trap is that the field at the center defines a proper quantification axis than can be used to select even more efficiently the desired transition through proper laser polarizations.

Therefore, in order to lift degeneracies and to facilitate the spectral amplitude selection shaping for the optical pumping we choose here a Ioffe-Pritchard configuration. Keeping achievable values [48][49][50][51][52], the magnetic field is taken to be

B(R) = B ′ (Xe X -Y e Y ) + (B 0 + B ′′ Z 2 )e Z with B 0 = 0.1 T, B ′ = 1 T/cm and B ′′ = 0.1 T/cm 2 .
In order to avoid particles that orbit around the trap center, we need to create asymmetries in the potential. We choose α = 0.17 creating an asymmetric configuration with a 1.17 T/cm gradient along X and a 0.83 T/cm gradient along Y.

We still simulate particles initially at a temperature of T = 100 mK and randomly distributed using equipartition of the energy.

Because the trap is shallower, the sample is bigger. The Sisyphus laser has now a waist of 6 cm, a 5W power and a central frequency such as the resonant transition is for particles having a potential energy of 3k B T where T is the expected temperature. We expect an exponential decay of the temperature T , and consequently of the sample size. Thus the power, the detuning, as well the as the (initially 35 MHz) FWHM Lorentzian laser linewidth (this is to follow the interaction time and required resonance sharpness evolution) of the Sisyphus laser are decreased exponentially with the same time constant (170 ms). The bias field creates a very useful preferential quantization axis (the Z axis) and we thus choose to propagate the lasers along this axis with circular polarization. For this simulation we use a realistic repumping laser with 10W power, a 5 GHz laser FWHM linewidth in order to repump all needed untrapped particles at a waist of 3 cm. We did not find it useful and necessary to reduce the waist of the lasers and this is a clear simplification for experimentalists. In order to repump only useful particles, that is the particles having a low enough energy, the spectrum of repumping the laser is flat and cut at a frequency corresponding to particles in M ′′ = 1 with a potential energy of 0.4k B T . We just mention that spectral shaping for amplitude selection has experimentally been realized for optical pumping experiments [53][54][55][56], with similar laser power and shaping resolution using simple broadband diode lasers and diffraction gratings for the shaping [57][58][59].

The results of the simulation, shown in Fig. 5, are impressive. With almost no losses, the temperature drops by a factor 1000 and the phase-space density increases by a factor 10 7 in only one second. Even better results are possible, but we may need to take into account collisional processes for such high phase space densities. We thus believe that further optimization of the parameters is useless.

IV. CASE WITH LOW FRANCK CONDON VALUE: OPTICAL PUMPING

With the previous NH case, we have treated a molecule with a very high Franck-Condon factor where vibration is frozen. However several molecules does not have such good properties (see Table I).

For such cases, with bad Frank-Condon factors, possible improvement of the Sisyphus cooling would consist of even further reducing the amount of (absorption)spontaneous emission cycles, for instance by using an external cavity [60] or coherent effects where the spontaneous emission step occurs only after several energy transfer processes [61][62][63][64][65]. But even, without using such coherent (often comlex) tricks, there are at least three possibilities to deal with such difficulties that we shall mention: :

• Use an optimized single photon cooling to avoid losses due to spontaneous emission.

We have already studied this solution and conclude that, if possible experimentally, three times 1D cooling may be more efficient that direct 3D cooling [17,22].

• Realize a Sisyphus cooling in the ground state using only ro-vibrational transitions as done in Ref. [21,22].

In light diatomic molecules the radiative lifetime of the excited ro-vibrational state is not too long and this solution can be efficiently implement. For instance the v ′ = 1 → v ′′ = 0 radiative lifetime of NH(X 3 Σ -) is determined to be 37 ms [66] (59 ms for OH [67]). Very powerful lasers, such as singlefrequency continuous-wave optical parametric oscillator (OPO), are nowadays commercially available to drive this transition (at 3.05 µm for NH).

Due to the long lifetime of the states it is preferable to work only with trappable states. This will often implies to work with M ′ = 2 (and thus J ′ = 2) upper state because it decays towards M ′′ = 1, 2, 3 lower states. Several lasers or microwave sources might be needed to repump these levels.

• Use an optical pumping scheme in order to repump the levels.

Obviously we can combined the three possibilities in order to find the best solution to the problem.

To help for further studies on Sisyphus cooling of such molecules, we deliberately choose to illustrate the third solution on a system, combining several "bad" points such as: no closed rotational neither vibrational transition scheme and quite poor Franck-Condon factors (≤ 0.5). Thus, we shall illustrate the method on a textbook example inspired by the CO case, but with scaled internuclear distance to get the Franck Condon factor of 0.5 towards the 3 Π 2 lower level. Quantities related to this level will be indicated by a double prime ′′ .

Then the 3 ∆ 3 upper level, indicated by a prime ′ , is preferred to the 3 Σ level because it forms with 3 Π 2 a closed system for spontaneous emission. Both potentials [199] are presented in Fig. 6 where the Franck-Condon factors, the rotational transitions and the trapping energy states are also displayed. They produce a Franck-Condon value near 0.5 for v ′ = 0 ← v ′′ = 0.

An obvious choice for the initial trapping U 1 is the lowest ro-vibrational level having the strongest trapped state that is the |v ′′ J ′′ M ′′ = |022 ′′ level. By choosing excitation only to |J ′ M ′ = |33 ′ , the rotational population can be locked in few trapped levels: |J ′′ M ′′ ∈ {|22, 32, 33, 42, 43, 44 ′′ }. Finally, in order to help the vibrational cooling, we do not use the best existing Franck-Condon factor (on v ′ = 1 → v ′′ = 0) but we deal with excitation towards v ′ = 0 to favor decay to low v ′′ values (see Fig. 6). Therefore for simplicity of both the explanation and the theoretical treatment, all lasers excitations end up to the sole |v ′ J ′ M ′ = |033 ′ .

We simulate particles initially at a temperature of 100 mK in an Ioffe trap with B 0 = 0.1 T,B ′ = 0.1 T/cm and B ′′ = 0.1 T/cm 2 .

FIG. 6: (Color online). Branching ratios, proportional to the broadening of the blue lines, between vibrational (left) and rotational (right) levels after transfer by the narrowband Sisyphus laser (dashed black) and the spectrally shaped amplitude selection broadband optical pumping (green hatched). J ′′ = 2 is not pumped for v ′′ = 0.

Several laser schemes or simple RF-microwave transfer can be used to make the Sisyphus or pumping transitions. We choose here a simple case with only two lasers: a single one to make the Sisyphus transition and a second one as the pumping laser near the trap center which shall also make optical pumping to bring back population from others ro-vibrational states. Thus, close to the trap center the (re-)pumping laser should bring the molecules back into the |022 ′′ state. A single laser would create a cylinder pumping zone. This would be a problem for particles moving mostly along its axis because they will never leave the cylinder and will be always under the effect of the repumping light. To be closer to the ideal situation of Fig. 1, with small r repump , we thus create a more spherical radial zone by dividing the pumping laser in three beams propagating along Ox, Oy and Oz. They are focused at the trap center on an initial 4 mm waist, all circularly left polarized (corresponding to σ + polarized for a quantization axis along the beam propagation) and their spectrum, 60 mW/cm -1 uniform power density, are shaped (amplitude selection) to drive only transitions towards |033 ′ . Finally, in order to accumulate population into the |022 ′′ level, all transitions from this state are removed by the spectral shaping.

Several possible experimental realizations exist for the repumping laser: spanning from femtosecond lasers, diode lasers, LED light, or even a supercontinuum laser, nowadays spanning 0.4-2.3 µm with 50 mW/nm uniform power density [68]. We suggest a simple experimental setup based on the use of several individual lasers driving individually the vibrational transitions. This would requires five diode lasers to cover the v ′′ = 0-4 → v ′ = 0 transitions respectively at 833, 971, 1163, 1449 and 1923 nm wavelength, that are all available commercially. The presence of the B 0 = 0.1 T bias field, lifts the energy levels through the Zeeman effect by 0.05 cm -1 . The required GHz range resolution is at the edge of what can be achieved by using a simple grating configuration, but higher resolution optical pulse shaping methods can be used [69]. We mention the optical arbitrary waveform line-by-line pulse shaping technique [59], based, for in-stance, on arrayed waveguide routers or liquid crystals on silicon [70][71][72], or the virtually imaged phased arrays that combine the high spectral resolution potential of etalons with the spectral disperser functionality of gratings [73]. Such devices have been applied to realize pulse shapers with sub-GHz spectral resolution [74][75][76]. But one simple method is the scanning diode method used in my group [56].

The Sisyphus transition, |022 ′′ → |033 ′ , is performed using a 1 W laser at 833 nm, 10 MHz linewidth, covering the whole sample with its 2 cm waist, left circularly polarized, propagating along Oz, with initial detuning chosen to initially transfer even very energetic particles with potential energy of 3.5k B 100 mK.

The waist size of repumping ans Sisyphus lasers, as well as their powers and the detunings are decreased exponentially with a time constant of 500 ms.

The results of the simulation is presented in Fig. 7 where the temperatures are calculated from the standard deviation of the histogrammed velocities. 25% of the |022 ′′ molecules are cooled from 100 mK to 1 mK within 1 s. Other molecular states are simply not repumped because they are in larger orbits, but no molecules are lost. Thus, repeating the same process every 1 s greatly improves the results with 60% of the molecules with temperature below 1 mK after 4 s. Clearly, there is plenty of room for optimization [200] and the lasers wavelengths, sizes, locations, polarizations and chosen transitions can be further optimized. For instance, simulation of a perfect Sisyphus scheme where particles are transferred and repumped when reaching their higher and lower radial coordinates leads to 100% of the molecules being cooled, from 100mK down to 1mK, within only 0.15 s.

V. GENERALIZATION OF THE METHOD

Before to conclude we would like to mention that the Sisyphus cooling method can be generalized with more complex potentials or absorption schemes for instance by during decelerators stages. We have illustrated in Fig. 8 such idea with molecules traveling repeatedly through uphill potentials. The idea is similar to the one proposed in Ref. [62] but now with a spontaneous emission step. The 3D simulation is done using a simple "toy model" spin 1/2 molecule (four levels A 2 Π 1/2 M ′ J = ±1/2 ← X 2 Σ, M ′′ J = ±1/2 transition) but is obviously more general. The global dynamics of the decelerator looks promising with negligible effect due to the transverse motion and only small Kinetic Energy dispersion due to photon recoil energy transfer. This scheme can be realized using permanent magnet and would allow continuous cooling.

FIG. 8: (Color online). Illustration of optical deceleration of a molecular beam using Sisyphus cooling when molecules travels through magnetic coils. Left Panel: illustration of the principle. Right Panel: Simulation using pairs of 3mm radius Helmholtz coils separated by 1.5cm and creating a 1T axial field. The circularly polarized lasers are 5 W focused with 1.5 mm waist, 1 GHz laser linewidth. The monokinetic molecular beams cover 1 mm radially.

VI. CONCLUSION

In conclusion, the fact that Sisyphus and optical pumping methods requires only a small number of spontaneous emission steps makes it especially useful for cold molecular systems with internal state modifications during the individual process. Similarly, an optical one-way pumping method can also help to realize transfer to a trap after a deceleration or a buffer gas cooling stage [9,10,12,13]. We stress here that a full and detailed simulation of the laser interaction is very important, a less detailed simulation can lead to incorrect results. We found that the laser power and its linewidth are parameters that are not so simple to tune. A relatively fast particle, a low power laser, an incorrect local polarization or a too small laser linewidth can result in transitions that are not effective. On the contrary a slow particle, a high power laser or a large linewidth can result in transitions occurring too early, out of resonance and not at the desired location nor energy location. We also found that a small modification of the trap can totally modify the efficiency of the cooling due to possible angular momentum quasiconservation. The spectral or spatial selectivity for the transitions have to be ensured and care has to be taken not to excite undesirable levels. A trap which lifts energy degeneracy allows simpler spectral selection of the transitions and in this case the cooling is found to be much easier by using preferential quantification axis allowing to control transition using light polarization.

We also show that three times 1D cooling may be simpler and more efficient that direct 3D cooling [17,22]. Indeed, in 1D the quantification axis is always the same and transitions can thus be controlled using light polarization, which can be extremely useful to excite only state which are trapped (for instance by controlling the Zeeman sub-level in a magnetic trap). Furthermore, rotation of phase space may be also advantageously used with only one laser direction to cool the 3D sample.

Obviously the Sisyphus method is very easy to implement for trapped quasi closed systems [77][78][79][80] and can then be a very impressive complement to standard laser cooling. For trapped species, the photon transition rate can be very slow and the scheme is particularly suitable for particles that require deep-UV lasers for electronic excitations. It also allows for the use of pure ro-vibrational transitions or of quasi-forbidden electronic transitions. Using state of the art (GHz) shaping resolution [56,59,69,75] the optical pumping step should not be a problem.

The amount of possible cooling schemes, using optical, magnetic, electric or gravitational forces combined with laser (or Radio-Frequency) transitions of any kind (Raman, coherent, stimulated, pulsed, shaped, multiphotonics, black-body ...) is so large and the method so versatile that it can be implemented to several species including ions [80,81], ranging from simple atoms to polyatomic molecules, that are difficult to laser cool, such as (anti-)hydrogen. It can be generalized to realize continuous cooling in new type of Stark or Zeeman decelerators. The low temperature achieved can then be a new starting point for further studies in precision measurements or in chemical control of highly correlated systems when collisions between molecules start to play a bigger role for lower temperatures.

Aknowlegments: I am indebted to N. Vanhaecke and H. Lignier for useful discussions. The research leading to these results has received funding from the European Research Council under the grant agreement n. 277762 COLDNANO.

Appendix A: Optimized strategy: one-way single photon cooling VS Sisyphus cooling

Trajectories in 3D traps

We describe (for simplicity) the trajectories in the case of harmonic trap but the discussion is general.

Trajectories in 3D harmonic traps are given, with (x, y, z) = (r 1 , r 2 , r 3 ) notation, by U 1 (r(t)) = 1 2 m i ω 2 i r i (t) 2 with r i (t) = r 0 i cos(φ i (t)) and φ i (t) = φ 0 i + ω i t. The total energy E = 1 2 m i ω 2 i r 0 i 2 is a conserved quantity. We want to reduce the kinetic energy E kin (t) = E -U 1 (r(t)). This constrains the trap geometry such that it should avoid angular momentum conservation where some particles can orbit around the trap center with a kinetic energy that never reduces: for instance if ω i = ω j and φ 0 i = φ 0 j + π/2 (dashed trajectory in Fig. 2). As studied in Ref. [82], we should look for a "chaotic" regime, with various different trajectories, increasing the likelihood of meeting the best condition for single photon cooling where kinetic energy is near zero: U 1 (r(t)) ≈ E that is when r 1 , r 2 , r 3 are all extremum. In the harmonic case, the coordinate r 1 becomes positive maximum when φ 1 (t) = 0[2π], at times t k1 = t 0 + k 1 T 1 , for k 1 = 0, 1, • • • , N 0 , spaced by the oscillation time T 1 = 2π/ω 1 . During these N 0 oscillations we want to have, at t k2 , r 2 (t k2 ) = r 0 2 cos(φ 0 2 + ω 2 t 0 + 2πk 2 ω 2 /ω 1 ) almost maximum. This optimal condition (r i maximum) should occur for all particles and thus for any trajectory independently of the initial condition φ 0 2 . A simple choice to realize such situation is ω 2 /ω 1 = 1/N 0 because during the N 0 periods, the phase φ 2 (t k1 ) spans the space [0, 2π] with interval of 2π/N 0 , and thus we found one φ 2 (t k2 )[2π] which approaches 0 with a maximum error of ∆φ = π/N 0 . This is the best possible result with N 0 points. Similarly, in N 2 0 oscillations, the choice ω 3 /ω 2 = 1/N 0 will allow φ 3 (t k2 ) to spans the space [0, 2π] with interval of 2π/N 0 . At a time t k3 where all positions are almost maximum, the maximum kinetic energy is

E -U 1 (r(t k3 )) = 1 2 m i ω 2 i (r 0 i ) 2 sin 2 (φ i (t k3 )) ≈ (∆φ) 2 E.
From this study we see that ratio of frequencies on the order of 1/N 0 are good. Coprime or incommensurable ratios may be also interesting but a simple (from the experimental point of view) choice is to modify slightly the trap geometry by a ratio 1 -1/N 0 . We have tested this hypothesis in the case of a simple quadrupole magnetic trap where the magnetic field is

B(R) = B ′ 0 (1+α)Xe X + B ′ 0 (1 -α)Y e Y ) -2B ′ 0 Ze Z with B ′ 0 = 100 T/m
and α is a parameter that controls the anisotropicity of the trap. We have looked for both for linear traps, i.e. with a state with a linear Zeeman effect, and for quadratic traps, i.e. with a state with a quadratic Zeeman effect: to be more precise the molecules used is NH in the exact same configuration than used for Fig. 3: using J ′′ = 1, M ′′ = 1 (linear)state or the J ′′ = 1, M ′′ = 0 (quadratic state). The typical "oscillation time" T 1 is 2 ms.

In order to see if the geometry allows trajectories to reach the best condition for cooling, i.e. with kinetic energy near zero for the Sisyphus pumping step and with potential energy near zero for the repumping step, we use 1000 molecules initially in thermal equilibrium. For each molecule we look for the minimal kinetic and potential energy reached during a given time: 5 ms or 20 ms, the first case corresponding to N 2 0 ∼ 5/2 and the second to N 2 0 ∼ 20/2. Then we took the worst molecules and plotted the result in Fig. 9 in function of the anisotropic parameter α. The simulation confirms that the anisotropicity of the trap is a crucial parameter and that α ∼ 1/N 0 and a minimal energy of E/N 2 0 are correct assumptions. In the case of an X, Y symetric trap (α = 0) it is clear than some molecules do not transfer their energy between kinetic and potential. Clearly the longer time we allow for the molecule to move the better the results is.

Energy or spatial selection in 3D traps

With this simulation and without entering into a more detail study [201], our general conclusion is that, in N 2 0 times the typical oscillation time T 1 along the fastest axis (say x), a good trap geometry, should allow the trajectories to reach the maximum energy

E within ∆E ≈ E/N 2 0 [202].
Finally, to be more precise we also have to take into account that the Sisyphus transition should occur when the particle stays in the E ± ∆E energy range, that is during a time ∆t ≈ T 1 /N 0 . This is N 3 0 times smaller than the overall time motion and, due to interaction time or power broadening of the transition, this can lead to lower energy resolution than the expected ∆E ≈ E/N 2 0 . The study of this effect is in fact the core of the simulation and will be done later numerically. However, a simple Lorentzian shape for the excitation rate (cf Eq. B3) leads to an energy resolution N 3/2 0 worse and thus ∆E ≈ E/N 1/2 0 . In conclusion, in order to have simple calculation, we assume an intermediate case of ∆E ≈ E/N 0 for all cooling strategies that we study. We have sketched them and the results in the Table 10.

The third column of Table 10 deals with the energy removal strategy (case 2 of Fig. 2). By looking at the typical case of E ∼ k B T , and with δE ≈ E/N 0 , we see that a single photon can lead to a reduction of the energy, or the temperature, by a factor N 0 . The depth (and shape) of the capture potential U 2 can thus be optimized, to be on the order of k B T /N 0 , i.e. in the harmonic case with frequencies ω ′ 2 i ∼ ω 2 i /N 0 . From an experimental point of view, such a flat potential can be realized using combined Zeeman and Stark effects, for instance [11]. In this case the density is unchanged and the phase space density D = N in a tight trap, cf black box in Fig. 2 Low Left Panel 1. This is the case study in the first column of table 10. The best spatial catch is when the particles reach the transfer zone of ∆E ≈ E/N 0 . However, the spatially located box is able to access only the turning points concentrated within its small region. For instance, from the figure Fig. 2, we clearly see that if the box was placed elsewhere with the same energy, the trajectory would never cross it. If we want to catch all particles, we have either to move the box to cover the all spatial transfer zone of ∆E ≈ E/N 0 (cf second and fifth column of table 10), or to move the trajectories to the box (using a small anharmonicity).

Here again in a 3D case we can cover this space in N 2 0 spots.

Another strategy is to use a second photon to pump the molecules towards a tighter trap when the molecules are close to the trap center, using for instance the Sisyphus cooling process (step 3 of Fig. 1, and first column of "Sisyphus step" of table 10). In this case, very similar to the dimple trap trick used in evaporative cooling [34,83], we need to transfer the molecule when the position is as close a possible to the center i.e. where the kinetic energy is maximal. The exact same discussion indicates that in N 2 0 times the typical oscillation time T ′ 1 ∼ √ N 0 T 1 in the new trap, the trajectories reach the maximal kinetic energy E ′ = E/N 0 within ∆E ′ ≈ E ′ /N 0 .

FIG. 10: (Color online). Typical evolution of number of molecules N , Temperature T (or energy E), phase space density D and time required t, for the sample under different d dimensional cooling strategies (pos. means a selection in position as in Fig. 2 1), Ener. a selection in energy as in Fig. 2 2) and Speed a selection with rate proportional to velocity as in Fig. 2 3). The values are given in the harmonic trap cases with angular frequencies initial ωin and final ω fin . The partial one-way ones takes only part of the population (with a box of given energy or spatial range). The Full one-way strategies move the previous box to catch all particles. The first Sisyphe step (step 3 of Fig. 1) put the population back in the original potential. The Sisyphus is thus simply the repetition of this step. The parameters, except the photons numbers, are given in power of N0 = kBT /∆E which is a parameter indicating the number of slices of the sample over which we catch particle by spacial (or energy) steps. The time evolution is given, in the "t" line, in function of the typical oscillation time T1 along the fastest axis motion: d -1 indicates for instance that the time evolve like t ∝ T1N d-1 0 .

The values given here are only indicative. But they can help to choose the parameters. For instance with T 1 ∼ 10 ms (ω 1 ∼ 2π(100 Hz)), a value of N 0 = 4 is possible even for a Full one way strategy spanning the all positions. Indeed, this gives in 3D (d = 3) a cooling time N 5 0 T 1 of 10 s which may be feasible. This should allow a huge gain in phase space density of 4 3d/2 ∼ 500 The first (1D) Sisyphus cooling scheme, leading to Eq. 1, was done with N 0 = 2, i.e. ω ′ = ω/2. When dealing with rotation, and angular momentum, several convention exist for the phase relation between coefficients, for the sens of rotations or for the name of the Euler angles. A good summary of the existing convention is given by [START_REF] Varshalovich | Quantum theory of angular momentum[END_REF]. We deal here with several different frames or axis: the lab fixed one (given by the vacuum chamber), the laser axis (given by the wavevector k); for a given molecule we have the local field B or E axis (depending on the position of the particle) and the molecular axis.

We use capital letter X,Y,Z for space fixed frame, with direction vectors (e X , e Y , e Z ). We use lower case letter for the local frame (e x , e y , e z ) with B (or E) field is along the Oz. We shall also use the (covariant) spherical and helicity basis vectors e 0 = e z , e ±1 = ∓ ex±iey √ 2 because they directly express light polarization. Another frame (e ′

x , e ′ y , e ′ z = k/ k ) is linked with the laser propagation. The original laser polarization is described in the helicity basis e

′ 0 = e ′ z , e ′ ±1 = ∓ e ′ x ±ie ′ y √ 2 .
We can write the laser polarization vector in the two frames as [START_REF] Varshalovich | Quantum theory of angular momentum[END_REF] Eq. 1-( 49):

ǫ = p ′ =-1,0,+1 ǫ ′ * p ′ e ′ p ′ = p=-1,0,+1 ǫ * p e p , with ǫ p = p ′ =0,±1 (-1) p ′ ǫ ′ p ′ D * (1)
pp ′ (0, θ, ϕ). We find:

ǫ +1 = 1 + cos θ 2 e -iϕ ǫ ′ +1 + sin θ √ 2 ǫ ′ 0 + 1 -cos θ 2 e iϕ ǫ ′ -1 ǫ 0 = - sin θ √ 2 e -iϕ ǫ ′ +1 + cos θǫ ′ 0 + sin θ √ 2 e iϕ ǫ ′ -1 ǫ -1 = 1 -cos θ 2 e -iϕ ǫ ′ +1 - sin θ √ 2 ǫ ′ 0 + 1 + cos θ 2 e iϕ ǫ ′ -1
where θ = ( Oz, Oz ′ ) = ( e z , k) is the angle between the local field and the laser propagation axis. When dealing with pure (circular) polarization and because the transition probability depend only on the modulus square of the transition amplitude the value of φ becomes irrelevant and we choose φ = 0.

Transition Matrix Elements and Line Strengths

The electromagnetic field, due to the lasers L, can be written

E(r, t) = 1 2 L E L e i(kL.r-ΦL(t)) + E * L e -i(kL.r-ΦL(t))
For each laser L the polarization vector, defined by E L (r, t) = E L (r, t)ǫ L , is (written directly in the local axis) ǫ L = p=-1,0,+1 ǫ * p e p and the irradiance, called improperly intensity, is

I L = ε 0 |E L | 2 c/2 (B1)
The light-molecule interaction hamiltonian is -μ.E, where μ is the electric dipole moment operator.

The Doppler effect, and the laser linewidth, are taken into account by writing

Φ L (t) = (ω L -k L .v)t + Φ L (t).
Where Φ L (t) is a fluctuating phase. As shown below, or through the Wiener-Khinchin-(Einstein-Kolmogorov) theorem, its statistical average is linked to the spectral irradiance distribution I ω (ω). In the simple Lorenzian case

I ω (ω) = IL π ΓL/2 (ω-ωL) 2 +(ΓL/2) 2 with FWHM Γ L .
We shall only describe the notation for Hund's case a, but for more complex cases I use the known molecular parameter of the molecule and then the transition strength and Zeeman or Stark effect shifts (∆E =

E 0 + Sign(C)[-∆/2 + ∆ 2 /4 + C 2 F 2 ]
for the of each levels where F is the field) given by PGOPHER (Program for Simulating Rotational Structure) see also [START_REF] Brown | Rotational Spectroscopy of Diatomic Molecules[END_REF].

We use the simplified notation

′ | μ.ǫ * | ′′ = Ω ′ v ′ J ′ M ′ | μ.ǫ * |Ω ′′ v ′′ J ′′ M ′′
where the electronical part is simply noted |Ω , the vibrational |v and the rotational |JM Ω . The integration on the electronical coordinates r i gives first the electronic dipole moment µ el (R) = Ω ′ q e i=1 r i .e * q Ω ′′ (R) depending on the internuclear distance R. With the short notation that q refers to vector quantized on the molecular frame whereas p refer to vector on the local (field) frame. The integration on the vibrational coordinate R gives the dipole moment of the transition

µ v ′ v ′′ q = v ′ µ el (R) v ′′ . If µ el (R)
is not constant we can assume a linear variation in R which leads to

µ v ′ v ′′ q = v ′ |v ′′ µ el (R v ′ v ′′ )
where the R v ′ v ′′ = v ′ |R|v ′′ centroid is the "distance" where the transition takes (classically speaking) place. When no other information is available we take the dipole moment of the transition from Einstein's formula for the lifetime τ ′ of the upper state τ ′-1 ≈

ω 3 Ω ′ Ω ′′ 3ǫ0c 3 π µ el 2 .
The integration over the rotational coordinates leads to the final results for the Rabi Frequency Ω of the transition:

Ω = ′ | μ.E L | ′′ = E L µǫ * p (B2) = E L ǫ * p v ′ |v ′′ µ el (R v ′ v ′′ ) × J ′ Ω ′ , 1q|J ′′ Ω ′′ J ′ M ′ , 1 -p|J ′′ M ′′
where p = M ′′ -M ′ and q = Ω ′′ -Ω ′ are the only values leading to non zero amplitude.

The intensity of a transition is thus proportional to the

Franck-Condon factor | v ′ |v ′′ | 2 .
FIG. 11: (Color online). Schematics of a levels i and j interacting with a broadband laser and with other levels.

Laser interaction

The optical Bloch equations for the density matrix element

ρ ij = i|ρ|j = ρ * ji are [86]: ρii = j Γ ji ρ jj - j Γ ij ρ ii + i j μij .Eρ ji -ρ ij μji .E ρi =j = -iω ij ρ ij - Γ coll + k (Γ ik + Γ jk ) 2 ρ ij + i k μik .Eρ kj -ρ ik μkj .E
Where for completeness we have added some possible collisional dephazing through Γ coll . The rate equation we want to derive can be justified only using several approximations [START_REF] Cohen-Tannoudji | Processus d'Interaction entre Photons et Atomes[END_REF][START_REF] Höppner | [END_REF] ; especially when coherence (non diagonal) terms are small and thus a dominant part is played by the population (diagonal) terms. We thus keep only the terms ρ kj or ρ ik when k = i or k = j, even if this prevent to simulate some dark states coherent schemes [203] Under this approximation we have:

ρii = j Γ ji ρ jj - j Γ ij ρ ii - j Im ρji L Ω L ij e i(Φij-ΦL) ρij = -i (ω ij -dΦ ij /dt) ρij - Γ coll + k (Γ ik + Γ jk ) 2 ρij + i 2 (ρ jj -ρ ii ) L ′ Ω L ′ ij e i(Φij-Φ L ′ ) + Ω L ′ ij * e i(Φij+Φ L ′ )
The Rabi frequency Ω L ij characterizes the strength of the transition between the states |i and |j and is defined by

Ω L ij = i| μ.E L |j e ikL.r
For more generality we have written these equations for ρij = ρ ij e iΦij (t) , where Φ ij (t) would have to be later chosen cleverly for instance to remove the time dependent exponential terms (Φ ii = 0). The next step [88] consists to formally integrate the second equation ρij -

b(t)ρ ij = L ′ c L ′ (t) in ρij (t) = L ′ t 0 c L ′ (t ′ )e t t ′ b dt ′ .
The standard choice of Φ ij = Φ L is not possible when several lasers L are present, thus we report the formula in the first equation which then does not depend on the choice of Φ ij . We finally make a statistical average (over several experiments) noted ρ ij =< ρ ij >.

The situation may be complex with, for instance, several lasers coherently driving the same transition to realize a travelling optical lattice plus incoherent lasers on the same or other transitions for optical pumping. Following [88], we shall model the phase Φ L (t) as a Poissonian random process with T L = 2/Γ L as the average time between successive phase jumps. The statistical average

< e i[Φ L (t)-Φ L ′ (t ′ )]
> is non zero only if there is no phase jump between t and t ′ . For the laser L this occurs with a probability P 0 = e -(t-t ′ )/TL . Thus, if L and L ′ are coherent, that is if Φ L -Φ L ′ have no more fluctuating phases (Φ L = Φ L ′ ), we have < e i[Φ L (t)-Φ L ′ (t ′ )] >= e -(t-t ′ )ΓL/2 . The net result of this statistical average is to remove the Φ L fluctuating phases and simply to add the laser linewidth Γ L , as done for the collisional terms Γ coll , to the coherence relaxation rate which becomes half of

Γ L ij = Γ coll + k (Γ ik + Γ jk ) + Γ L a. Rate equation
Even if coherences, or Bloch equations, are tractable with a Monte Carlo method as described in [89], it is much simpler to use rate equations for each particle. We thus use the method of steepest descent or stationary phase, by integrating by part to write

t 0 dt ′ Ω L ij (t ′ )e -[i(ωij±ωL)+Γ L ij /2](t-t ′ ) ∼ Ω L ij (t) i(ωij±ωL)+Γ L ij /2 .
This is valid only for low saturation that is when the saturation parameter

s = 2|Ω L ij | 2 4(ω ij -ω L ) 2 + Γ L ij 2 ≪ 1
The resonant terms, that is when |δ L,0 ij = ω L -ω ij | ≪ ω ij , are the dominant one and terms containing e i[ΦL+Φ L ′ ] , will be dropped (so called rotating-wave approximation) because they are much smaller than resonant ones and because the statistical average leads to exponential decreasing terms. A slow time variation due to small frequency difference between coherent lasers L, L ′ driving the ij transition could still occur but, in order to avoid remaining oscillating terms, we shall suppose it smaller than 1/Γ L ij . We define the detuning

δ L ij = ω L -ω ij -k L .v = δ L,0 ij -k L .v + ∆ j -∆ i
where ∆ i and ∆ j are the energy shifts created by the external trapping potentials. We combined lasers in class "L c " of mutually coherent lasers by writing

L = Lc Lc
L . An important parameter is the total Rabi Frequency for the class L c :

Ω Lc ij = Lc L Ω L ij e -iδ L ij t = Lc L i| μ.E L |j e ikL.r(t) e -i(ωL-ωij)t
Where r(t) = r + vt, which in the simulation is taken to be the current location of the particle.

We finally obtain the rate equations used in our simulation and illustrated in Fig. 11:

ρii = j [Γ ji ρ jj -Γ ij ρ ii + γ ij (ρ jj -ρ ii )] (B3) γ ij = Lc |Ω Lc ij | 2 Γ Lc ij 2 + 4δ Lc ij 2 Γ Lc ij (B4)
The physical interpretation is obvious, γ ij is the rate of excitation, but also of stimulated emission, of the transition and is given simply by the incoherent sum over the individual rates.

With i| μ.E L |j = | i| μ.e p |j E L ǫ p |e iφ L p we see that the rate is proportional to the total laser taken into account interferences due to lasers polarizations:

|Ω Lc ij | 2 = Lc L |Ω L ij | 2 + (B5) 2 Lc L ′ >L |Ω L ij Ω L ′ ij | cos[(ω L -ω L ′ )t -(k L -k L ′ ).r(t) -(φ L p -φ L ′ p )]
where

| Ω L ij | = |µǫ p E L |.

b. Bloch equation with laser linewidth

For completeness, we mention that in the case of a single class of coherent lasers L c , we can choose one laser L 0 of this class and Φ ij = Φ L0 to derive generalization of the Bloch equation. Time derivative of the statistical average equations, before the integration by part, leads to Bloch equations with laser linewidth:

ρii = j Γ ji ρ jj - j Γ ij ρ ii - j Im (ρ ji Ω ij ) ρij = iδ Lc ij ρij - Γ Lc ij 2 ρij + i 2 (ρ jj -ρ ii )Ω ij
where -ΦL) . We obviously recover the rate equation (B3) using ρij = 0 and we find ρij = A full quantum theory of light shift and heating in an optical trap can be found in Ref [90], but we shall here follow a simple deviation based on Ehrenfest's theorem. Atomic motion is affected by the force F = μ .∇E where μ = T r( ρ μ) = ij ρji μij . We first write it in the single coherent lasers class case, because general result is simply the sum over the classes. The statistical average leads to F = i>j Re ρ * ij ∇Ω ij where we have separated the sum using levels with ω i > ω j . This can be written as:

Ω ij = L Ω L ij e i(Φij
F = -i>j (ρ jj -ρ ii ) γ ij 2δij Γ Lc ij Re ∇Ωij Ωij -Im ∇Ωij Ωij
where a clear separation between the dipolar and the scattering force appear F = F dip + F scat .

In the general case we found:

F dip = -∇U dip ; U dip = Lc i>j γ ij (r) δ ij Γ Lc ij (ρ jj -ρ ii )
In our Monte Carlo simulation, a molecule is always in a given state |j . Thus, in addition to the previous shifts created by the magnetic and electric fields, the dipolar potential indicates that we have to shift the j energy level by

Lc i>j δ Lc ij γ Lc ij (r) Γ Lc ij - Lc k<j δ Lc kj γ Lc kj (r) Γ Lc kj (B6)
We recover the standard result for a two level system that for red detuning the upper state is up-shifted and the lower one is down-shifted by the laser. This picture is not the same that the usual dressed state one [14] where particles are in dressed state and in a superposition of |i and |j levels with a light shift given by δ ln(1 + s). However we recover the same average shift in our Monte Carlo simulation because we oscillate in time between |j and |i . Finally, we mention that we use this shift only to calculate the forces but we do not add it as a real shift in the levels. The energy levels are thus not shifted by lasers and no new laser detuning are calculated. The main reasons is to avoid unphysical accumulation of dipole shifts. Indeed, if we calculate a new detuning with this formula, at the next update of the energy levels a novel detuning will be added to this one and we would have accumulation of the detuning that would not correspond any-more to reality.

The scattering force is given by

F scat = Lc i>j γ Lc ij (ρ jj -ρ ii )Im ∇Ω Lc ij Ω Lc ij |Ω Lc ij | 2 Im ∇Ω Lc ij Ω Lc ij = Lc L k L |Ω L ij | 2 + 2 Lc L ′ >L k L + k L ′ 2 × |Ω L ij Ω L ′ ij | cos[(ω L -ω L ′ )t -(k L -k L ′ ).r(t) -(φ L p -φ L ′ p )]
The photon absorption rate is given by γ Lc ij and the photon momentum can be chosen depending on the "participation" of the laser beams to the total irradiance by comparison with formula B5. In the simple example of a retro-reflected laser creating a standing wave lattice (Ω L ij = Ω L ′ ij and k L = -k L ′ ) no force is present but the absorption rate is non zero.

More general cases, especially in the saturated regime, can be treated with ad-hoc formulas (see [46]) but we have not implemented them.

d. General laser spectrum

We have, up to now, followed a standard procedure to derive the absorption rate (B3): that is statistical average (and deconvolution approximation) followed by an integration by part [START_REF] Höppner | [END_REF]88]. However, the result seems non-physical in the far-off resonance case, when

|δ Lc ij | ≫ Γ Lc ij .
We found an absorption rate that is proportional to the laser linewidth, where obviously, because the laser is very far from resonance, only the laser intensity, not its linewidth, should play a role [91]. The problem arises in a Cauchy-lorentzian laser irradiance spectrum because the spectrum does not decrease fast enough and has still a large amplitude even far from its center. But, the problem is more general and even for other laser spectrum it arises from the fastest oscillating term in e -[i(δ L ij +Φ L )+Γ L ij /2](t) that is now due to the detuning and not any more to the laser phase fluctuation. Thus, in such far detuned case the integration by part has to be performed before the statistical average. The net result, for this far off resonant case, is that the phase fluctuation plays no more role and thus Γ L has to be removed from Γ L ij which becomes Γ L ij = Γ coll + k (Γ ik + Γ jk ) = Γ ij in the absorption rate formula (B4).

Fortunately, the near and far off resonant results can be combined together in a single formula. Indeed Eq. (B4) can be written, in the single laser case to simplify, as

γ ij = L ω ⊗ I ω (ω ij + k L .v)π|µǫ p | 2 /(ǫ 0 c) that is pro- portional to the convolution of the Cauchy-Lorentzian natural linewidth L ω (ω) = 1 π Γ ij /2 ω 2 +(Γ ij /2)
2 by a Cauchy-Lorentzian laser irradiance spectrum I ω (ω). However for more realistic laser spectral shape such as a Gaussian one the formula

γ ij = [L ω ⊗ I ω ](ω ij + k L .v)π|µǫ p | 2 /(ǫ 0 c) (B7)
agrees with our expectation even in the far-off resonant case where the formula agrees with the well known re-

sult γ ij = |Ω L ij | 2 4δ L ij 2 Γ ij .
As another example, for very broadband lasers with a spectrum including the resonance, the Lorentzian natural linewidth can be considered as a Dirac peak and γ ij is simply proportional to the irradiance at the transition wavelength I ω (ω ij + k L .v). The physical interpretation of Eq. (B7) is obvious: in an incoherent (rate equation) model, a laser lineshape can be seen as a sum of several narrowband (delta or narrow Gaussian) lasers lines of intensity I ω (ω), so, in a low saturation and perturbative approach, the total rate of excitation is thus simply the sum of the incoherent rates.

In our simulation, we keep the same spirit, and for each coherent lasers class, without making a precise integration of the formulas over the laser spectrum, we calculate the rates using Eq. (B7) and then forces simply by replacing Γ L ij by Γ ij in the formulas.

e. Parameters needed

We want here to discuss the irradiance needed to efficiently perform the transition. This is usually defined by the saturation irradiance (intensity) I sat defined as γ = I Isat Γ 2 (where we have dropped the ij subscript). With a laser linewidth at least as large as the transition linewidth we have

Γ L ij ≈ Γ L and γ ≈ Ω 2 π 2 Iω(ωij) I . At resonance γ ∼ µ 2 E 2 L ΓL 2 Γ 2 L
. Including the Franck-Condon (FC) and angular factors (AF), Einstein's formula for the natural lifetime 1/Γ leads to γ/FC.AF ∼ 6πc 2 ω 3 I Γ ΓL and so γ/FC.AF ∼ 300 000 s -1 λ 500 nm

3 I W/m 2 Γ Γ L
Depending on the laser transition, the needed transition rate is not the same. Clearly for the Sisyphus transition the transition probability should be near unity when the particle arrive at the transition point; whereas for the repumping transition we may have much more time to do it (except if it is needed to repump the particle during the Sisyphus transfer), typically a fraction of the trapping period.

For the narrow linewidth Sisyphus laser, in order to have an efficient transition, the time γ -1 to realize the transition should match the time where the particle spend at resonance (that is when |δ| < Γ L ). Due to the (thermal) velocity v = k B T /m ∼ 1 -10 m/s in the trap (depth k B T and size σ r ∼ 1 mm), the time spend at the resonance condition is only ΓL kBT σr v . We thus find γ ∼ kBT Γ v σr . For an harmonic trap σ r ∝ T /ω 1 and the irradiance needed is thus proportional to the temperature. This explain why we have to reduce the power during the cooling process. As example, a (150 ns lifetime, 500nm wavelength) 1 MHz transition natural linewidth require γ ∼ 10 7 s -1 and an irradiance of I ∼ 100W/m 2 .

For the case of a laser used for optical (re-)pumping, we can simply choose its broadband linewidth Γ L such that the particle is always at resonance, i.e. Γ L ∼ k B T . In practice we use a much broader, and thus a correspondingly more intense, laser in order to cover several rotational or vibrational levels. For T ∼ 100 mK this require a linewidth of Γ L ∼ 2π(2 GHz) and for a transition rate of γ -1 ∼ 1 ms an irradiance of 7 mW/cm 2 .

Kinetic Monte Carlo exact solution of the master or rate equation

We recall here the method, which have been published in Ref. [47], used to perform the simulation and especially the Kinetic Monte Carlo method for solving rate equations and the N-body integrator.

The method is able to solve any system driven by a master equation

dP k dt = N l=1 Γ kl P l - N l=1 Γ lk P k (B8)
This equation describes the time evolution of the probability P k of a system to occupy each one of a discrete set of states numbered by k. Each process occurs at a certain average rate Γ lk (t). With, P k = ρ kk we clearly see that our system is described with such a master equation. The states ensemble contains all molecule internal states and the rate are the the spontaneous emission plus the lasers absorption and stimulated emission rates described in Eq (B3).

One of the simplest algorithm to solve this equation, sometimes called fixed time step algorithm, is based on the first-order formula P k (t + dt) ≈ P k (t) + N l=1 Γ kl (t)P l (t)dt -N l=1 Γ lk (t)P k (t)dt. The main disadvantage of this fixed time step method [92] is that dt has to be small enough to maintain accuracy and such that at most one reaction occurs during each time step: meaning Γ lk (t)dt ≪ 1. The Kinetic Monte Carlo (KMC) algorithm solve this problem by choosing optimal time step evolution of the system. Furthermore, the KMC method makes exact numerical calculations and cannot be distinguished from an exact molecular dynamics simulation, but is orders of magnitude faster. KMC method is indistinguishable from the behavior of the real system, reproducing for instance all possible data in an experiment including its statistical noise. Surprisingly enough, up to now it has been more or less limited to the study of chemical reactions, surface or cluster physics (diffusion, mobility, vacancy motion, transport process, epitaxial growth, dislocation, coarsening, ...).

The Kinetic Monte Carlo algorithm uses the fact that the system has a Markovian behavior. For a system initially at time t in state k, the probability that the system has not yet escaped from state k at time t ′ is given by exp t ′ t l Γ lk (τ )dτ . At time t ′ a reaction takes place, so just before t ′ the system is still in state k, we therefore have to generate a new configuration l by picking it out of all possible new configurations with a probability proportional to Γ lk (t ′ ).

The KMC algorithm is then the iteration of the following steps.

• Initializing the system to its given state called k at the actual time t.

• Creating the new rate list Γ lk for the system, l = 1, . . . , N .

• Choosing a unit-interval uniform random number generator [204] r: 0 < r ≤ 1 and calculating the first reaction rate time t ′ by solving t ′ t N l=1 Γ lk (τ )dτ = -ln r.

• Choosing a unit-interval uniform random number generator r ′ : 0 < r ′ ≤ 1 and searching for the integer l for which R l-1 < r ′ R N ≤ R l where R j = i=1,j Γ ik (t ′ ) and R 0 = 0. This can be done efficiently using a binary search algorithm.

• Setting the system to state l and modifying the time to t ′ . Then go back to the initial step.

Simple method to solve the N-body equations of motion

Discussion of possible solutions to solve the N-body equations of motion is given in Ref. [47]. We simply recall here the main conclusions.

Numerical methods such as the ordinary Runge-Kutta methods are not ideal for integrating Hamiltonian systems because they do not conserve energy. On the contrary symplectic integrators does conserve energy.

Because of the computational round-off error and stability domain issues, algorithms are generally not good to go beyond third or fourth order. One very popular N-body integrator is the fourth (local) order "Hermite" predictor-corrector scheme by Makino and Aarseth [START_REF] Heggie | The gravitational million-body problem[END_REF][START_REF] Aarseth | Gravitational N-body Simulations: Tools and Algorithms[END_REF][START_REF] Zwart | [END_REF]. We use here a simpler but still efficient algorithm, the so called velocity leapfrog-Verlet-Störmer-Delambre algorithm: It is of O((∆t) 3 ) accuracy for both position r and velocity v for a ∆t timestep. It has also the big advantage that accuracy can be improved by using higher order symplectic integrators such as the one by Yoshida [96].

In our case the velocity can be modified directly by photon recoils in absorption or in emission and the acceleration is simply calculated using the gradient of the potential. A typical time-scale is given by the motion in laser (or trapping) fields. Thus the time-step should be small fraction of the ratio of the waist, or the trap size, over the thermal velocity of the particle.

Finally, we combine the KMC and the N-body integrator in the following way: we first calculate an expected (KMC) reaction time t ′ -t. Then, if t ′ -t < ∆t, that is a reaction should occurs before the motion, we move, through Eq. (B9) the particles, using the timestep t ′ -t and the reaction takes place at time t ′ . On the contrary, if ∆t < t ′ -t, a dynamical (N-body) evolution of the system is made, no reaction takes place but due to its markovian probabilistic behaviour, the system is still govern by Eq. (B8). Then, after each change, of the position or of the internal state, the laser fields and potentials are recalculated and new transition rates are calculated. It is convenient to choose a ∆t timestep such that the calculated laser excitation rates are almost constant over it, in order to allow us to calculate the reaction time t ′ -t = -ln r N l=1 Γ lk (t) . We use these rules of thumb to start the simulation but we finally reduce the timestep ∆t until we obtain convergence of the results which usually occurs for γ∆t ∼ 1. [125], KRb [126], RbCs [127], NaCs [128], LiCs [129] 100 200 000

Method

Three body collision Rb 2 [118] Li 2 [106,130] 0.2 2 000 000 laser cooling SrF [6], YO [131], CaF [132] 0.3

Cryogeny (Buffer Gas) [133] CaH [134], VO [135] CaF [136], PbO [137,138], O 2 [139] NH [140], ND, CrH, MnH [141], ND 3 , H 2 CO [142], YbF [143], NH, NH 3 , O 2 , ThO, Naphthalene [START_REF] Doyle | Cold, Trapped Molecules via Cryogenic Buffer Gas Methods[END_REF], BaF, SrO, YbF, YO [133] 400 000 10 12

Field Slowing: Stark [START_REF] Van De Meerakker | [END_REF] CO [24], NH 3 , ND 3 [146,147], OH [148,149], OD [150], H [199] We calculate them using CO reduced mass, harmonic potentials and rigid rotor energy levels given by ω(v + 1/2) + BJ(J + 1). Equilibrium point, vibrational frequency and rotational constant for respectively lower and higher potential are: 1.2a0, 1700 cm -1 , 1.7 cm -1 and, 1.37a0, 1200 cm -1 , 1.3 cm -1 . The v ′′ = 0, v ′ = 0 transition being at 12000 cm -1 .

[200] In our simple scheme a lot of particles are transferred into |v ′′ = 0, J ′′ = 2, M ′′ = 2 levels, having same trapping potential than the initial one and not a less stiff one as desired for good Sisyphus cooling. This is a waste of time because a pumping step is required before another cycle.

[201] A proper study of a general, non harmonic, trap case should be done using revolving orbits, orbital precession, recurrence or Poincaré analysis for the investigation of this dynamical systems. This should illustrate in more detail which trap geometry is an efficient one.

FIG. 2 :

 2 FIG. 2: (Color online). Illustration of possible strategies to remove the kinetic energy from the sample. The upper part illustrate the process in a 1D picture whereas the lower part indicates it in a 2D picture. The transfer position is indicated by an hatched area. Two kinds of trajectories, but with the same kinetic energy, are shown: the dashed red line correspond to an isotropic trap case ωx = ωy and the solid line with ωx = 0.8ωy. 1) Left panel: the selection is done spatially (in the hatched region shown by the line or box) to catch the molecule where it has almost zero kinetic energy. 2) Central panel: the selection is done spectrally depending on the potential energy of the particle when it has almost zero kinetic energy. 3) right panel: the selection is done depending on the velocity through the lifetime of the (absorption or spontaneous emission) transition.

FIG. 3 :

 3 FIG. 3: (Color online). Principle of NH Sisyphus cooling of molecules. The width of the arrows are proportional to the branching ratios.

FIG. 4 :

 4 FIG. 4: (Color online). Results of NH Sisyphus cooling of molecules in a linear (quadrupole) trap with energy degeneracy at center. Molecules number N (left axis) in the M ′′ = 1 (square) and M ′′ = 0 (circle) states as well as the temperature T (triangle, right axis) in function of time.

FIG. 5 :

 5 FIG. 5: (Color online). Results of NH Sisyphus cooling of molecules in a quadratic (Ioffe-Pritchard) anisotropic trap. Molecules number N in the M ′′ = 1 and M ′′ = 0 states and temperature T are display on the left axes. The relative phase space density D is display in right axes in function of time.

FIG. 7 :

 7 FIG.7:(Color online). Time evolution of the (stacked) number N of particles in |v ′′ J ′′ state and the (three axis) averaged temperature T of the particles in |022 ′′ . The sequence is repeated every 1 s.

  in the d dimensional case. It is possible to increase furthermore the phase space density by transferring the particles in a tighter trap without heating. This can be done by catching them

FIG. 9 :

 9 FIG. 9: (Color online). Effect of the trap anisotropy, on the minimum kinetic E kin (lower panel for a quadratic trap) and potential energy Epot (higher panel for a linear trap) reached after 5 ms (black square) or 20 ms (red circle) evolution time of N = 1000 molecules initially at 100 mK. The worst energy case among the N = 1000 molecules is plotted, in temperature units, in function of α. 1 -α can be seen as the aspect ratio of the trap (see text for detail).

  Appendix B: Rate equations, Kinetic Monte Carlo and N-body integrator used 1. Polarization and Rotation matrices

  ρ jj ) which, as expected, has a smaller value, by the saturation parameter s, than the populations ρ iiρ jj . c. Dipolar potential and light shift

  [202] By removing the positivity condition on the coordinates, in the previous discussion, we have twice more maximum and ∆φ becomes π/2N0 ∼ 1/N0.[203] It is sometimes possible to use the rate equations, even with coherent dark states, simply by changing from the initial basis states to a new basis which includes the bright and dark states such as (|i ± |j )/ √ 2.[204] In our case we use the free implementations by GSL (GNU Scientific Library) of the Mersenne twister unitinterval uniform random number generator of Matsumoto and Nishimura.

  85,87 Rb 2[101][102][103], Cs 2 [104], 40 K 2 [105], Li 2 [106-108], Na 2 [109], 40 KRb [110], 41 K 87 Rb [111], Cr 2[112], Li 3[113], NaK[114] Cs 2[116], H 2 [117], Rb 2 [118], Li 2 [119], Na 2 [120], K 2 [121-123], He * 2 [124], Ca 2

	Molecule	T(µK)	N
	Feshbach, RF [99, 100]	0.1	100 000
	Photoassociation [115]

  2 CO[151], NH[152], SO 2[153], C 7 H 5 N[154], CO[151], ND 3[173], D 2 O [174], CH 3 F [175], CF 3 H [176], CH 3 CN [177], H 2 O , D 2 O, HDO [178], NH 3 , CH 3 I, C 6 H 5 CN, C 6 H 5 Cl [179]BeH + , YbH + [183], AF350 + =C 16 H 14 N 2 O 9 S + [184], MgH + [1], O + 2 , MgO + , CaO + [185], H + 2 , H + 3 [186], BaO + [187], NeH + , N + 2 , OH + , H 2 O + , HO + 2 , ArH + , CO + 2 , KrH + , C 4 F + 8 , R6G + [182], Cyt 12+ , Cyt 17+ [188] GAH + =C 30 H 46 O 4 [189] HCN, Ne -, Kr -, ArHF, tetracene-Ar, HCN-H 2 , HD, D 2 , Ag8-Ne N , Ar N , Kr N , Xe N (N = 1 -135), NaCs, LiCs, HF-(H 2 ), OCS-(H 2 ) N (N = 1 -17), Cs2 , Rb 2 , CO, HCl, amino acid, 3-hydroxyflavone, xanthine, [Na(H 2 O) N ] + (N = 6 -43) [196], HF-N 2 O, Mg-HF, Mg-(HF) 2 , CH 3 -H 2 O, Cs + N (H 2 O) M , CF 3 I, CH 3 I, ...

			10 000 1 000 000
		YbF [155],
		LiH [156], CaF [157]
	Rydberg	H 2 [158-161]
	Optical	C 6 H 6 [162], NO [163]
	Zeemann	O 2 [164]
	Beam collision	NO [165], KBr (13K) [166],	400 000
		ND 3 [167]
	Beam dissociation	NO [168]	1 600 000
	Rotating Nozzle	O 2 , CH 3 F, SF 6 [169, 170], CHF 3 [171],	1 000 000
		perfluorinated C 60 [172]
	Velocity filtering	H 2 1 000 000	10 9
	Sympathetic cooling [180-182]		20 000	1000
	Nanodroplet [190-195]	Mg 1-3 1 000 000	10

TABLE II :

 II Slow and cold molecule list with T < 1 K with an estimation of their temperature T and number N . For each species only the reference corresponding to the pioneering work is given. Program for Simulating Rotational Structure, C. M. Western, University of Bristol, http://pgopher.chm.bris.ac.uk.

	e
	and r ′′ e , leads for the Franck-Condon (in SI units) 2 √ ω ′′ e ω ′ e ω ′ e +ω ′′ e -µ(r ′ e -r ′′ e ) 2 ω ′′ e ω ′ e / (ω ′ e +ω ′′ e ) . e [198] A

Appendix C: List of formed cold molecules

The following table presents a list of the cold molecular species with T < 1 K. For each species only the reference corresponding to the pioneering work is given. The combination of the buffer gas pre-cooling with the Stark velocity filtering is listed either under the "Cryogeny" or under the "Velocity filtering" category. The occurrence of some "exotic" molecules has been detected only as a loss process in a degenerate gas but not really produced and isolated. This is the case of Cs 3 and Cs 4 [97] and thus they do not appear in the table. Similarly, several photoassociation experiments have produced electronically excited cold molecules (as single example Yb [98]) but we do not include them in the list as they have a lifetime limited to some tens of nanoseconds. Furthermore, we do not always detail the isotope difference, for instance NH 3 can be the same as 15 NH 2 D.