Direct Monte Carlo simulation of the sympathetic cooling of trapped molecules by ultracold argon atoms

We demonstrate that cold molecules (H2 and benzene) can be created at temperatures below 1 mK by sympathetic cooling with laser-cooled rare gas atoms on timescales of seconds. The thermalization process is studied using the direct simulation Monte Carlo (DSMC) method, which allows a detailed analysis of the atomic and molecular spatial and energy distributions as a function of time. As part of this study, ultracold elastic cross sections for Ar–Ar and Ar–C6H6 are also calculated.


Introduction
There is currently considerable interest in the creation and application of cold molecular gases that range in temperature from a few kelvins down to the nK range. At temperatures below mK the thermal energy is sufficiently low so that the molecules can be trapped and manipulated by external electric [1], magnetic [2] or optical fields [3,4]. The ability to trap these gases allows measurements with much longer interaction times than was previously feasible, allowing highresolution spectroscopy [5] and measurements of the lifetimes of long-lived excited states [6]. Cold trapped molecules are seen as ideal candidates to measure parity violation at the molecular level [7] and to test extensions to the standard model by measurement of an electron dipole moment [8]. At still lower temperatures below 1 mK, the quantum mechanical nature of the centre-of-mass motion becomes important in the interactions between molecules in a trapped gas. In this regime, a variety of effects have been predicted including novel quantum phases in gases of polar molecules, and ultracold chemistry that is dominated by tunnelling and resonance phenomena.
Although simple cold molecules can be created by photoassociation and on Feshbach resonances, the means to cool more complex molecules to sub-mK temperatures is yet to be demonstrated since they cannot be laser cooled. Sympathetic cooling, where one species can be cooled by thermalizing collisions with a colder species, appears at first to be a promising route to achieve this because it is a well-established method for cooling atomic systems where the collision partners are below 1 mK. It has, however, not yet been demonstrated for mixtures of complex molecules and cold atoms where the molecules start in the mK regime. This is due to the many more degrees of freedom in molecular systems where collisions are more likely to be inelastic, leading to heating or loss of molecules and atoms from a trap. As a result, much of the cold molecule literature over the last few years has shown that this process is not likely to work over the wide range of collision energies even for many simple diatomic molecules and laser-cooled atoms [9]- [11].
In this paper, we show that even complex species such as benzene can be cooled to sub-mK temperatures on timescales of seconds using experimentally achievable parameters. To 3 do this, we use elastic thermalizing collisions between ultracold atomic rare gas atoms and molecules [12], which can be loaded into a deep optical trap. This scheme is attractive when compared to many other sympathetic cooling schemes because it reduces the possibility of inelastic collisions within the trap, which prevent thermalization. The choice of a rare gas atom as a collision partner minimizes the possibility of reaction between the collision partners while trapping in an optical field allows atoms and molecules to be trapped and cooled in their absolute ground state.
This paper reports on a detailed numerical study of the thermalization process of a mixture of atoms and molecules, based on the direct simulation Monte Carlo (DSMC) method [13]. Specifically, we have considered the thermalization of H 2 and C 6 H 6 , using Ar as a buffer gas. Although the experimental setup is quite general, we have decided to initially focus on these two molecules as they play crucial roles in several different scientific fields. For example, H 2 is the simplest molecule and has been extensively studied as a benchmark for comparing theory and experiments. Furthermore, it plays a major role in astrophysics. C 6 H 6 is representative of the large class of organic molecules. Essential to the simulation is the availability of reliable ultracold energy cross sections or, alternatively, accurate potential energy surfaces (PESs) from which the scattering observables can be deduced. In this respect the interaction between Ar and each of the two molecules has been the main focus of several studies. Recently, we calculated the elastic collision cross sections for rare gas collisions with H 2 and C 6 H 6 in the cold to ultracold regime. In this paper, this work is extended to include higher partial waves, and the elastic cross section for Ar-C 6 H 6 is also calculated for the first time in the cold regime (1 µK < E < 1 K).
The paper is organised as follows. Section 2 presents the DSMC method. In section 3 the microscopic parameters relevant to the simulation are analysed, with particular emphasis on the inter-species elastic cross sections. Section 4 discusses the macroscopic parameters for the simulations, such as the initial densities and energy distributions of the atoms and molecules, as well as the main approximations introduced in the simulation. Section 5 presents the main results of this work, and section 6 is devoted to the conclusions.

Method
A computationally efficient method based on a Monte Carlo simulation of the microscopic collisions was pioneered by Bird [13], and is referred to as Bird's method. The method was originally developed to describe non-equilibrium gas flow in a number of problems in different scientific fields, such as aerospace, plasma modelling or engineering. An advantage of this method is that it can be applied to processes where the atoms/molecules have very large mean free paths that are outside the range of applicability of the Navier-Stokes equations. In molecular physics, the Boltzmann equation is the accepted mathematical model for studying gas flows. However, as it is a six-dimensional partial differential equation, it can become computationally demanding. It has been shown that Bird's method provides an accurate representation of the molecular flow regime, which is consistent with the Boltzmann equation. In recent years, this method has been successfully used to study the dynamics of evaporative cooling for Bose-Einstein condensation [14] and to model the collisional properties of ultracold atom mixtures [15]. Bird's method offers further advantages for studying these processes: it does not require the assumption of sufficient ergodicity, and it can be easily generalized to describe complex phenomena such as associative ionization, dissociation or trap loss, which are inaccessible to the Boltzmann formulation. In the following, we review the main features of the technique and its application to sympathetic cooling, in order to present our notation. The key idea is the decoupling of the in-trap molecular motion from the inter-molecular collisions over a time interval much smaller than the mean collision time. All atoms and molecules positions and velocities are mapped at the initial time t = 0 and then the particles are allowed to evolve freely according to the laws of classical dynamics. All particles experience an external potential U (x, y, z), which, for simplicity, in this paper is assumed to have a harmonic shape given by The parameter U 0 indicates the trap depth and 2 k the spring constant along the k-axis. In equation (1) the distances (x, y, z) are defined from the centre of the trap, which is located at (0, 0, 0).
The associated classical trajectories for this potential can be integrated exactly, and the motion along the different axes is completely decoupled. For example, along the x-axis, the motion is given by where x 0 and v 0 are the initial position and velocity at time t 0 , and ω x is the oscillation frequency, ω x = x √ U 0 /m, with m being the particle's mass. A time step t is introduced to discretize the time variable. At each time step the effect of inter-particle collisions is considered for all possible pairs. A random number is used to decide whether a specific pair collides, based on the following expression for the collision probability [16]: where v is the relative velocity in the centre-of-mass frame, µ is the reduced mass of the colliding pair, σ is the elastic cross section that depends on the collisional energy in the centreof-mass frame (E = µv 2 /2), and V is the enclosed volume. If the random number is less than the collision probability, P, the collision is accepted, and the centre-of-mass orientation of the post-collision momenta is assigned by using two further quantum numbers. It is customary to divide the trap into in macro-cells of volume V and to allow collisions only for particles occupying the same macro-cell. In fact in equation (3) there is no reference to the inter-particle distance, and a pair sitting at the opposite edges of the trap has the same probability of colliding as two atoms very close to each other, for the same relative speed. The free parameters V and t are chosen to satisfy P 1. Although all pairs within the same macro-cell are considered for collisions, the actual number of collisions is determined by an acceptance-rejection test. When a collision is accepted the particles maintain their positions, but experience a sudden variation of velocity as determined by two random angles in the centre-of-mass frame and by the conservation of momentum. When a particle's energy exceeds the trap height, the particle is removed from the simulation.

Microscopic parameters
Knowledge of the microscopic processes involved, such as the atom-molecule cross section, is an essential ingredient of the simulation. In the energy region of interest, i.e. below 100 mK,   elastic collisions are dominated by the s-wave contribution. Furthermore, in order to make the simulation numerically more efficient, the elastic cross section at low energy, restricted to the s-wave contribution, is described in [15] by an analytical expression that depends on two parameters, the scattering length a s and the effective range r eff : where k is the modulus of the relative wave vector (E =h 2 k 2 /2µ), E is the incident energy in the centre-of-mass frame and µ is the reduced mass of the complex. Figure 1 compares equation (4) to the full elastic cross section for Ar-Ar and Ar-H 2 , which were calculated including partial waves up to g waves (see below). It can be seen that as expected the s-wave contribution is dominant in the energy region of interest and that the cross section is well reproduced using equation (4). Figure 2 gives a comparison of the cross sections for different complexes, obtained with equation (4). The values of the parameters a s and r eff employed are summarized in table 1. The appendix describes in detail how the scattering length and effective range were obtained for each collision complex.

Trap
This section provides a brief outline of the experimental process and setup; a more detailed description has been given by Barker et al [23]. The main purpose of this section is to discuss the initial values for the DSMC simulation from an experimental perspective. A deep trap is required in order to hold the molecules in the mK range produced by Stark deceleration, while at the same time it must also be sufficiently deep to also trap the ground state rare gas atomic species. A quasi-electrostatic optical trap is an ideal universal trap in that  essentially all species can be trapped. However, in this type of trap, well depths in the 100 mK range are challenging and require high intensities. Such a trap can be created by utilizing a highfinesse build-up cavity to increase the intensity from tens of watts to the kW regime required. Trap depths of up to 150 mK are feasible for benzene and 25 mK for argon. The trapping potential inside a standing wave high-finesse cavity is an optical lattice, which is a periodic array of small cylindrical wells with a period of half the wavelength of the light used. In our case, this corresponds to a period of 532 nm with a radial size of 55 µm. The trapping potential is given below [23]: where w = 55 µm is the beam waist, and λ = 1.064 µm is the beam wavelength. While the potential shape is the same for all particles, its depth U 0 is proportional to the field intensity and the particle's polarizability, and it is therefore different for each species. Rare gas species in their ground state have been used, chosen because of their nonreactivity as compared to laser-cooled alkali metal atoms, as has previously been discussed [12,23]. All rare gas species must, however, be cooled in a long-lived metastable state where laser sources are available. Once cooled they can be quenched to their ground state. While helium is in many cases the best collision partner for cooling many molecular species, it cannot be easily quenched to its ground state. In addition, its recoil temperature following quenching would be too high for it to be useful as a collision partner. Collisions between ultracold ground state argon and molecular hydrogen have been shown to be favourable when compared to all other available rare gases except helium [12]. While all rare gas densities are limited by Penning ionization once quenched into their ground state, this process is not important and much higher trapped densities are feasible. The number of Ar atoms that can be loaded into each lattice site and their average final temperature are yet to be determined experimentally. However, in order to predict a characteristic thermalization time we estimate a range of atomic and molecular densities based on feasible loading schemes.
The ultracold argon atoms that are to be loaded into the optical trap are first created in their metastable state via a radiofrequency discharge. These atoms are then decelerated in a Zeeman slower and loaded on a magneto-optical trap. We can currently produce a metastable argon magneto-optical trap (MOT) that is approximately 1 mm in diameter with a density that is limited by Penning ionization to approximately 10 10 cm −3 . The density of ground state Ar available for sympathetic cooling is, however, not limited to this density, because when the metastable atoms are quenched to the ground state they are not subject to Penning ionization or light-assisted collisions.
One feasible method for loading atoms into the trap is to rapidly switch on the optical trap field and capture the atoms. This has the effect of heating the atoms and is thus not considered a viable loading scheme for sympathetic cooling. We instead propose to load the ground state argon atoms into a dipole trap that is overlapped by the MOT. This allows the metastable atoms to be loaded into the dipole trap while maintaining approximately the same temperature before they are quenched to the ground state. The large Stark shift in the deepest part of the trap will be used to selectively quench the coldest metastable atoms that are captured by the same trap. As the trap for ground state argon atoms is approximately 20 times less deep than the trap for the metastable atoms, we can also expect some cooling of the argon atoms due to conservation of phase-space density in the shallower trap. In addition, the atoms can be continually quenched to their ground state since the MOT can, in principle, be continually replenished. In this loading scheme the maximum densities are limited, therefore, by collisions with the ground and metastable atoms. Using similar loading schemes starting from the same initial MOT densities, trapped atom densities of 3-6 × 10 13 cm −3 have been achieved for rubidium atoms in a dipole trap overlapped by a MOT [24,25]. This number therefore seems a realistic prospect also for ground state argon accumulated in the trap. We therefore calculate the cooling dynamics in the trap based on densities that are approximately an order of magnitude lower and higher than this. We consider two well depths in use for sympathetic cooling. The first is the maximum value that can be produced by the maximum intensity available and corresponds to 24.7 mK for ground state Ar, 12 mK for H 2 and 150 mK for C 6 H 6 and a second lower set of values (6.4 mK-Ar, 2.4 mK-H 2 , 33 mK-C 6 H 6 ), which allows us to load in more atoms at the expense of the number of molecules that can be captured by the trap. Finally, once loaded with the cold atoms a rapid non-adiabatic loading of molecules into the same trap from the Stark decelerator would allow us to bring the cold atoms into contact with the hotter molecules without significantly heating the atoms. Table 2 reports, for each species, all relevant macroscopic parameters, including the trap height U 0 and the estimate number of particles initially trapped in each lattice site. For the Table 2. Parameters defining the initial conditions of the simulation. For each species the following is reported: the mass m, the trap height U 0 , the temperature T along each axis and the number N of particles trapped, assuming an initial density of 10 12 molecules cm −3 and 10 14 atoms cm −3 . Two different sets of trap depths are considered. atoms, the number of particles has been estimated by considering a gas of harmonic oscillators at equilibrium with temperature T = 100 µK and average density of 10 14 particles cm −3 . We consider the molecules as a free gas uniformly distributed in space, with a density of 10 12 particles cm −3 , and at equilibrium, but with different temperatures along the axes. This is a consequence of the slowing process, as the molecules are formed in a beam and subsequently slowed down by optical Stark deceleration. As a result, they will show a large spread in energy along the propagation axis and a much narrower one in the perpendicular directions. Their phase-space distribution can be described as The number of molecules trapped, N , is estimated by means of a six-dimensional integral of the phase-space distribution with the constraint that the total energy E has to be less than the trap height U 0 :

Approximations
In the description of the evaporative and sympathetic cooling of the mixture of atoms and molecules, a number of approximations are introduced. First, we treat the confining potential as a harmonic trap, with a cutoff at U 0 : This approximation allows for a significant reduction of the computational burden of the simulation. It fails when particles have an energy close to the trap height and as a consequence have a much longer oscillation time. This approach could in principle lead to an underestimation of the thermalization time. Table 3 reports the time of oscillation along the xand z-coordinates for Ar, H 2 and C 6 H 6 , for several energies expressed as a percentage of the trap height. Table 3. Oscillation time τ inside the trap, along the radial and axial directions, for the potential of equation (5) for different energies expressed as a percentage of the trap height. The last row reports the harmonic oscillation times as obtained from equation (8). Note the different timescales for r and z.  . The table shows that the harmonic approximation is reasonable for energies as high as half the trap height, i.e. 6, 12 and 75 mK for H 2 , Ar and C 6 H 6 , respectively. On the top of the trap, the deviation is within a factor of 3. However, as will be shown later, although the molecules have an initial energy comparable to the trap height, they cool very rapidly below the half line, and the total cooling time is not affected by overestimating the first part of the process. The number of particles trapped in each lattice site is, however, estimated using the exact trapping potential of equation (5). A second key assumption is that we neglect any macroscopic quantum effect of the trap. That is, while the inter-particle dynamics is governed by the Schroedinger equation, we consider the interaction between the particles and the external field U (r, z) to be purely classical. At the same time, we consider an MB distribution for the gas. Table 4 shows the quantum of energy for each species in the trap. It can be seen that the classical approximation is particularly poor at describing the motion along the z-axis. On the other hand, a quantum mechanical description of the particle's motion is extremely difficult to achieve computationally, and most previous works have followed the same classical assumption. Bose-Einstein statistics was considered in [14,16]. However, applications of the DSMC method to describe Rb-Cs mixtures have shown good agreement with experiment to within a few per cent for the thermalization rates [15].
A final comment must be made on the statistical properties of the DSMC simulations. Such simulations were originally intended for very large systems (N ≈ 10 10 ). The statistical assumptions underlying the method rely on a large number of collisions. In previous works similar to our work [16,26], samples of more than 5000 particles were considered and even higher densities were simulated introducing the concept of the macro-particle, that is, the description of several atoms or molecules at once as a bunch. From table 2, however, it can be noted that our simulations deal with at most 10 3 particles. In particular, the optimal number of molecules is constrained by imposing a large atom/molecule ratio in each lattice site. Therefore, simulations with as few as ten molecules per site are presented and discussed. Although the statistical character is partly recovered by averaging over a large number of such simulations, it must be stressed that such low-density simulations may not be fully reliable, and according to our experience a minimum number of about 50 particles per species should be included in each simulation. This is because if the number of particles is too low, parameters such as the macro-cell and the time step have a tremendous effect on the simulation, which is contrary to the assumption that the simulation results should be a weak function of those parameters.

Results
In order to study the thermalization of the molecular gas, we have considered a range of possible initial spatial distributions for both atoms and molecules. To achieve a sensible reduction of the temperature of the molecules, the number of atoms in the trap must be at least one order of magnitude more than the number of molecules, as the final temperature can be estimated as the weighted mean of the two initial temperatures. In fact, if the number of molecules n m is taken as a fraction of the number of atoms n a , n m = xn a , and the aim is to reach a final temperature T , it can be easily estimated that x T /T m , where T m is the initial average energy of the molecules. However, it must be considered that sympathetic cooling is not the only process in the trap but also elastic losses play a role in reducing the gas temperature. This phenomenon, as we will see, is particularly enhanced for benzene.
When presenting our results we prefer to refer to the number of particles per lattice site rather than the macroscopic atomic or molecular density. In fact, when the trap is loaded, the gas density goes from being uniform to becoming strongly inhomogeneous, as the particles are concentrated in many small clouds that alternate with empty spaces. The macroscopic density thus depends on which volume is chosen to enclose the particles. The number of particles per lattice site, in contrast, is an indisputable quantity. A second subjective quantity is the thermalization time. In our analysis, it has been defined as the first time when the atom's and molecule's temperatures become closer than 10 µK. However, after this time, both temperatures oscillate around the thermalization value, due to the statistical fluctuations associated with the smallness of the samples. The error bars on the thermalization temperature reflect the amplitude of the oscillations. Furthermore, the concept of thermalization refers more properly to an equilibrium state of the gas, where the number of particles in a small energy window is constant over time and corresponds to the MB value. During the cooling, both atoms and molecules, and particularly the latter, are almost never energetically at equilibrium. We define the gas temperature T as that of a Boltzmann gas with the same mean energyÊ; therefore T =Ê/3k B , where k B is the Boltzmann constant.
Finally, the computer code used to produce all the results presented in this section will be the subject of a separate future publication, where the numerical aspects of the simulations will be explained in detail [27]. For the simulations presented in this paper, we have divided the trap into about 20 cells. The cell distribution of the molecules varies during the simulation. As they cool they tend to concentrate in the inner cells, while initially they are distributed more uniformly. The atoms are nearly always in the central cell. The time step is chosen of the order of 10-50 µs, and the simulations are run for about 1000 000 time steps. For completeness, all the inputs necessary to reproduce the results of this paper using COOL [27] are given in the supplementary material, available at stacks.iop.org/NJP/12/113002/mmedia.

Ar-H 2
Three initial spatial densities have been considered for both the atoms and the molecules for each pair of trap depths, corresponding to the trapping of 100, 250 and 500 atoms and 10, 25 and 50 molecules in each lattice site for U 0 = 24.7 mK (12 mK) for atoms (molecules) and 50, 500 and 3000 atoms and 10, 75 and 300 molecules for U 0 = 6.372 mk (2.397 mK). For each pair of initial densities, 20 different simulations were run, using for each a different seed for the initialization of the random sequence. In all simulations, the atoms were given random energies around 300 µK, whereas the molecules were initially distributed with a truncated MB with temperatures of 400 µK along the xand z-axes and 1 K along the y-axis, truncated at U 0 . This corresponds to an initial molecular average energy of 2.6 mK for the deeper trap and 0.6 mK for the shallower trap. Table 5 reports estimates for the final temperature T (in µK) the thermalization time τ (in s), and the percentage of atoms and molecules left in the trap, for the two sets of trap depths considered. As the trap depth for Ar is higher than that for H 2 , all atoms remain in the trap, whereas a small percentage of molecules evaporate due to inter-molecule collisions. The effectiveness of the cooling process depends on the atom/molecule ratio; the larger this ratio the lower the equilibrium temperature of the system. When a small number (10) of molecules are considered, thermalization requires a much longer time than the one considered in the simulations. Apart from this inconvenient situation, in most of the cases considered, the final temperature reaches the sub-mK regime. The thermalization time is normally less than 10 s. Figure 3 shows the atomic and molecular temperatures as a function of time, for the case of 250 atoms and 50 molecules. It also shows the percentage of molecules left in the trap as a fraction of the initial sample. Figure 4 shows the atomic and molecular spatial (upper row) and energy (lower row) distributions at three different times for the same case as figure 3. The spatial distribution shows the particle's position along the Rand z-axis at the selected time averaged over one oscillation period, and not the instantaneous position. This is because the instantaneous position can be misleading, while the average position gives a much more faithful idea of the amplitude of the oscillation for each particle and does not depend on the specific time when the 'observation' is made. For harmonic motion the average position is directly connected to the particle's energy, so the plot gives an indication of the energy distributions of the gases, where the greater the distance from the origin, the larger the particle's energy. For a similar reason there are very few particles near the origin, as the density of states becomes very small, and decays as E 2 as E → 0 for a gas of harmonic oscillators. Furthermore, this effect is enhanced by the reduction of the spatial volume near the origin, which goes as R 2 . The energy distribution plots (bottom row) show histograms of the fraction of the total number of particles as a function of their energies. The plots also show the equivalent MB energy distributions, that is, the ones with the same mean energy of the particles. The two panels in Table 5.

Ar-C 6 H 6
The analysis of the thermalization of benzene follows closely the discussion of the previous subsection. Two sets of trap depths are considered, and for each nine pairs of initial atomic and molecular densities, corresponding to 10, 25 and 50 molecules in the deeper trap, and 10, 75 and 300 in the shallower trap. The thermalization process was simulated using 20 runs with a different seed for the sequence of random numbers for all 18 pairs of initial atomic and molecular densities. The atoms were initially given random energies around 300 µK.
The molecules were initialized with truncated MB distributions with temperatures indicated in Table 2. The trap depth for benzene is much higher (150 mK) than that for H 2 , which is reflected in the rather large initial energy of the benzene molecules. In particular, the average benzene energy is 35 mK for the deeper trap and 7.5 mK for the shallower trap. This means that most molecules have an initial energy greater than the Ar trapping potential, resulting in fast evaporation of Ar when the hot molecular gas is co-trapped with the cold atomic gas. We have therefore investigated two different trapping schemes, one in which the atoms and molecules are simultaneously co-trapped at time t = 0 and a second in which first the molecules alone are trapped and then the atoms are added to the trap with a time delay of a few seconds. Table 6 shows the effect of loading pure benzene in the trap, for a time of 5 s. The table shows the initial and final temperatures and numbers of molecules for each of the six different densities considered. It can be seen that pre-thermalization has a very small effect when the molecular density is low, due to the small probability of molecule-molecule collisions. On the other hand, for the 300 molecule case, the molecular temperature is greatly reduced after 5 s. The price to pay is a significant reduction of the number of molecules. The initial and final energy distribution for this case is shown in figure 5. Not only the mean energy is sensibly reduced, but also the energy distribution becomes smoother and more symmetric around its mean value.
The results are reported in table 7. As expected, for low atom/molecule density ratios the final temperature is quite high. Also, the thermalization time becomes quite large. In some cases thermalization is not completed within the time range of the simulation (10-40 s). When the two Table 6. Pre-thermalization of pure benzene in the trap. The thermalization time is 5 s for all cases. The first column indicates the number of molecules in each site, and all the results refer to a set of 20 runs. The other columns report the initial and final temperatures, and the initial and final numbers of molecules. species are co-trapped at time t = 0, a large portion of the Ar atoms are expelled from the trap, reducing the effectiveness of the sympathetic cooling, as the Ar density becomes too small. In more than one case all Ar atoms are expelled from the trap. All the atomic densities considered for set 1 are probably too low to effectively cool benzene to sub-mK temperatures. In fact, due to the very high initial energy of benzene an atom/molecule ratio of 10 is probably too small. As the molecular densities for set 1 are also small, the effect of pre-thermalization is also negligible. A more efficient process is obtained for set 2, where larger densities are considered for both atoms and molecules. In particular, figures 6 and 7 describe the case of 300 molecules/500 atoms with pre-thermalization. The right panel of figure 6 shows the atomic and molecular The same as table 5 but for Ar-benzene. The last two rows report the results for pre-thermalized benzene.
The case highlighted in bold is illustrated in detail in figures 6 and 7. In the last two rows, the thermalization time is considered from the moment the two gases are co-trapped.   Figure 7. Spatial (upper row) and velocity (lower row) distributions for the Ar atoms (black) and benzene molecules (red), at time t 0 = 0, t 1 = 1 s and t 2 = 2 s, from left to right. In the velocity distribution plots, the histograms show the fraction of the total number of particles with energy within a chosen range, while the curves represent the Boltzmann distributions with the corresponding mean energy.
temperatures and numbers of particles as a function of time. The left panel shows the same quantities when atoms and molecules are loaded simultaneously in the trap. In the latter case, all atoms are expelled from the trap and the molecular temperature stabilizes at 2000 µK.
In the former case instead the atomic evaporation is limited and the molecular temperature reaches the sub-mK regime. Figure 7 shows the spatial (upper row) and energy (lower row) distributions at three different times, namely t 0 = 0, t 1 = 5 s and t 2 = 15 s, for the scheme with pre-thermalization. While initially the molecular and gas clouds are very well separated in space, after a short time they occupy the same volume. Also the energy distributions show a typical Maxwell-Boltzmann shape at time t = 15 s.

Conclusions
This paper presents a detailed analysis of the thermalization process for a molecular gas (either H 2 or C 6 H 6 ) in a buffer of ultracold Ar atoms. The macroscopic parameters utilized in these simulations are based upon achievable experimental conditions [23]. The simulations allow us to study the thermalization process as a function of time for different initial densities of atoms and molecules. We obtain density distributions and energy profiles of the two gases as a function of time, which can then be compared directly with future experimental observations. The two molecular species considered have allowed us to investigate two physical scenarios quite different from each other: one where the molecule polarizability is smaller than the atom polarizability (i.e. H 2 ), and the opposite situation (benzene), where the molecule is more polarizable than the atom. In the first case, we found that the molecules reach sub-mK temperatures for a wide range of initial densities. This result can now be used to guide the design of the experiment and to give an indication of the cooling times expected for a variety of different atomic and molecular densities in the trap. The other situation, where the molecular polarizability is larger, makes the cooling more difficult to achieve, as the initially large molecular energy results in significant atomic evaporation, which strongly reduces the effectiveness of the cooling process. The combination of initial large densities and a two-step procedure allows one to improve the cooling process. The final temperatures are estimated to be 330 ± 30 µK for H 2 and 600 ± 100 µK for C 6 H 6 , in a time of less than 10 s.
The DSMC method employed here is rather general, and only requires an a priori knowledge of the inter-particle cross sections. In this application, we have considered only elastic scattering, but the method can be easily generalized to include inelastic channels. Under other experimental conditions, where both elastic and inelastic channels are available, it is difficult to estimate the effectiveness and feasibility of sympathetic cooling. In those cases, DSMC calculations can provide a quantitative analysis of the cooling process. At the same time, however, in order to simulate such complex experiments, where several hundreds of particles interact among themselves and the confining field, a number of approximations were introduced. These include the approximation of a harmonic trapping potential and the use of an analytic expression for the cross section with collision energy. These could be relaxed if necessary, but would result in time-consuming calculations. The key approximations employed are the use of classical trajectories for the particles inside the trap, and that we neglect the role of the strong optical fields on the interactions between the particles. The study of a gas with a quantized centre-of-mass motion induced by the tight binding of a steep trap is extremely computationally demanding, but was pioneered by Zoller and Gardiner [28]. Similar studies for Bose-Einstein condensates have shown very good agreement between theory and experiment [15]. The effect of strong optical fields on scattering observables can be analysed and will be the subject of our future work.
The study of collision cross sections between all possible colliding pairs is an essential step for the simulation of evaporative cooling of molecules and sympathetic cooling of molecules by ultracold atoms. In this paper, we reviewed the Ar-Ar interactions, and we presented new results for the Ar-H 2 and Ar-C 6 H 6 cross sections, whereas the H 2 -H 2 cross section was taken from the work of Quéméner et al [22]. The C 6 H 6 -C 6 H 6 cross section was estimated. The detailed study of Ar-H 2 and Ar-C 6 H 6 collisions will be the subject of our future works. Due to the resonant character of the Ar 2 dimer, comparable PESs yield significantly different values for the scattering length and therefore the ultracold elastic cross section. However, as they are orders of magnitude larger than the Ar-molecule cross section, this uncertainty should not significantly affect the reliability of simulations because Ar-Ar collisions are much more probable than Ar-molecule collisions. In fact, a first experimental application of cooling apparatus could be used to resolve the uncertainty in the Ar-Ar PES by measuring the thermalization of a pure gas of Ar atoms. Molecule-molecule collisions are highly unlikely owing to the low trapped density of the molecules and therefore an accurate knowledge of the elastic cross section associated with this process was not essential for our simulations. A completely different case, however, is when pre-thermalization is considered: in this scheme, the magnitude of the molecule-molecule cross section is fundamental for the process to be effective. However, we have determined that a scattering length of 10 Å is sufficient for pre-thermalizing benzene over a timescale of 1 s.

Appendix. Inter-particle elastic cross section
The next subsections describe in detail how the elastic cross section was obtained for each complex.

A.1. Ar 2
The Ar 2 molecule is considered a benchmark system, and has been extensively studied, both theoretically and experimentally, with an emphasis on characterizing the Ar-Ar PES. The HFDID1 empirical potential of Aziz [17] is considered the most accurate available. Recently, two highly accurate ab initio potentials have also been developed, by Patkowski et al [29] and Slavicek and et al [30]. Furthermore, Tang and Toennies produced an extensive collection of semi-empirical rare gas dimer potentials, including Ar-Ar [31]. Table A.1 gives our calculated values for the scattering length and for the effective range for these PESs. The table also shows the number N of vibrational bands supported by each potential. The scattering length varies over a wide range, due to the resonant character of the Ar 2 molecule. With the exception of the Tang and Toennies PES, the other three surfaces predict a state very close to dissociation and therefore produce a very large scattering length. Under these circumstances, a small difference between two PESs can result in very different values for a s . In particular, the PESs from Slavicek et al [30] and Patkowski et al [29] support a ninth vibrational state very close to dissociation, resulting in a large positive scattering length, whereas the HFDID1 PES [17] only supports eight states, with a ninth state just above dissociation, resulting in a large negative scattering length. Figure A.1 shows the zero-energy wavefunction for Ar 2 calculated with different PESs. It can be noted that the four curves are in very good agreement in the short-range region, and differ Comparison of four different Ar-Ar interactions. The quantities considered are the zero r 0 of the potential, its minimum r min and the depth V min of the potential well, the number N of vibrational states supported, the energy E 9 of the ninth vibrational state, the scattering length a s and the effective range r eff .
Reference  Zero-energy wavefunction u 0 for the Ar 2 complex calculated with the different Ar-Ar interactions discussed in table A.1. From top to bottom; using the Ar-Ar interaction of (a) Tang and Toennies [31], (b) Aziz [17], (c) Slavicek et al [30] and (d) of Patkowski et al [29]. At large R, where R is the Ar-Ar distance, the difference between the Tang and Toennies PES and the other three PESs is particularly noticeable. All the interactions are very similar in the shortrange region, yielding nearly indistinguishable wavefunctions. markedly for R > 10 Å. This behaviour can be easily explained as all potentials support a state very close to dissociation, and as a consequence the zero-energy wavefunction is sensitive to the smallest change of the potential. Table A.1 also reports three quantities characterizing the PES, namely the zero of the curve r 0 , the minimum r min , and the well depth V min . These quantities also appear to be very similar for the four PESs. A detailed analysis of where the PESs deviate most from each other or which part affects the scattering length most is beyond the scope of this work.
Despite the large spread in the scattering lengths, we can conclude that three PESs are in good agreement with each other for this work, while the PES of Tang and Toennies yields a particularly small a s compared to the other three. In particular, the cross sections associated with each pair of values (a s , r eff ) are all much larger than the cross sections of the other processes, as is shown in figure 2. Therefore, the probability of Ar atoms colliding with each other is higher than the probability of their colliding with the molecules. This fact is further supported by the larger number of atoms in the trap compared to the molecules, as will be discussed in section 5. It is thus safe to conclude that the particular value of the Ar-Ar scattering length is not critical to the thermalization process, provided that it yields a cross section much larger than the Ar-molecule one. We have therefore assumed the values corresponding to Aziz's PES, as this semi-empirical PES has been calibrated over a large number of experimental data.

A.2. Ar-H 2
The low-energy cross section for the Ar-H 2 complex was studied in detail by Barletta [19]. The main findings of that work are that, unlike the Ar-Ar case, five different PESs available in the literature yield comparable results, with a spread of only 20% in the scattering length. Furthermore, the three-body Ar-H 2 complex can be effectively modelled as a two-body system interacting with an effective potential. The calculations of [19] were restricted to the s-wave contribution. For this work, we have calculated higher partial waves, up to l = 4, using a onedimensional model as indicated in [19]. Those higher waves, as shown in the right panel of figure 1, contribute very little in the energy region where most of the collisions take place, namely 100 µK < E < 10 mK. In this work, the values of a s = 9.930 Å and r eff = 7.00 Å were used, which correspond to those given by the PES of Bissonnette et al [18].

A.3. Ar-C 6 H 6
The elastic Ar-C 6 H 6 cross section was determined using MOLSCAT package [32] with the PES from Pirani et al [20]. These calculations closely parallel our previous calculations of the Ne-C 6 H 6 cross section at low temperatures [33].

A.4. H 2 -H 2
The elastic cross section for H 2 -H 2 collisions has been analysed in detail recently [22,34,35], using different available PESs. Particularly useful for the purpose of our work is figure 7(a) of Quéméner et al [22], which shows the cross section for several orders of magnitude of the colliding energy. There is a strong shape resonance at E ≈ 3 K, well above the trap height. In fact, in the cold regime, the cross section can be considered constant. The value of the H 2 -H 2 scattering length is not reported explicitly in the cited works, but it can be estimated from the cited figure to be a s ≈ 12.5 Å, corresponding to σ (E = 0) = 1960 Å 2 . The effective range can be estimated to be very small as the curve is flat, and therefore we set it to zero.
A.5. C 6 H 6 -C 6 H 6 The determination of the ultracold C 6 H 6 -C 6 H 6 cross section is a very difficult task, both experimentally and theoretically. Due to the low density used for benzene in our simulations, the precise value of the cross section is not important as the dominant process will be Ar-benzene collisions. In this work, we assumed the value a s = 10 Å and set the effective range to zero.