Kinetics of a single trapped ion in an ultracold buffer gas

The immersion of a single ion confined by a radiofrequency trap in an ultracold atomic gas extends the concept of buffer gas cooling to a new temperature regime. The steady state energy distribution of the ion is determined by its kinetics in the radiofrequency field rather than the temperature of the buffer gas. Moreover, the finite size of the ultracold gas facilitates the observation of back-action of the ion onto the buffer gas. We numerically investigate the system's properties depending on atom-ion mass ratio, trap geometry, differential cross-section, and non-uniform neutral atom density distribution. Experimental results are well reproduced by our model considering only elastic collisions. We identify excess micromotion to set the typical scale for the ion energy statistics and explore the applicability of the mobility collision cross-section to the ultracold regime.


I. INTRODUCTION
Trapped ion systems are among the most promising candidates for quantum information processing [1,2], precision measurements [3], and quantum chemistry [4]. For many of these applications it is required to cool the ions to low temperatures. To this end, various techniques such as laser cooling, resistive cooling, sympathetic cooling by other ions, or buffer gas cooling are routinely used. Ultracold atomic gases have recently become available in hybrid systems with trapped ions [5][6][7][8][9][10], extending the concept of buffer gas cooling to ultralow temperatures. Understanding this new regime and how it relates to conventional buffer gas cooling is essential for any future application to trapped ions.
Cooling of an ion in a buffer gas is caused by elastic collisions. They are dominated by the long range polarization interaction, and cross-sections are considerably larger than for collisions between two neutral atoms [11][12][13][14][15][16][17][18][19][20][21]. Therefore, collision rates are large and cooling is expected to be efficient [22]. Moreover, it has been proposed that internal degrees of freedom of ionic molecules could also be cooled by using an ultracold buffer gas [23]. Another specific feature of the polarization interaction is that collisions affecting the ion's mobility happen with rates independent of the collision energy [24]. This can lead to simplified system behaviour and has, for example, been applied in ion mobility spectrometry [25].
The motion of an ion in a radiofrequency (RF) trap can be decomposed into a fast driven motion, the micromotion, and a slow secular motion. In every collision the energy of the RF-field couples via the micromotion to the neutral atom's energy and the ion's secular energy. This can lead to ion energy removal or intake, sensitive to the RF-phase in the moment of the collision [26][27][28]. The average effect of many collisions results in cooling or heating and depends on parameters such as the ratio * cdz22@cam.ac.uk between the mass of the ion m i and the mass of the neutral atom m n . For very heavy neutrals runaway heating of the ion is expected [26], whereas very light neutrals enable efficient buffer gas cooling. In the mass ratio regime between the two extremes increased ion trap loss can be observed [29,30]. This effect has been explained by non-thermal energy distributions of the ion obtained from Monte Carlo simulations [31]. Such numerical simulations are a well established tool to model a trapped ion interacting with a buffer gas [32][33][34]. They can account for energy dependent scattering rates, complex electric field geometries, and other experimental parameters, which are difficult to treat analytically. In previous calculations, buffer gases at ambient temperatures have been assumed, and the cooling of the ion has been limited to the buffer gas temperature. However, in the recent experiments with ultracold neutral buffer gases a new energy scale related to excess micromotion has become dominant. The direct relation between the ion's mean energy and the excess micromotion has been observed in [8], for a system of Yb + − Rb.
Here, we investigate the kinetics of a single ion colliding elastically with an ultracold buffer gas, by applying Monte Carlo techniques. The effects on the ultracold neutral cloud are modelled using a semiclassical differential cross-section. The results on neutral atom loss and temperature increase, and the dependence of ion energy on excess micromotion are in good agreement with experiments.
The paper is organized as follows: In section II we describe the basic simulation procedure and the underlying physical model. The classical Langevin interaction model is applied in section III to derive the ion's energy statistics depending on the mass ratio, trap geometry, and scattering rate. In section IV we explain effects on the neutral atom cloud as a result of the energy dependent differential cross-section and the non-uniform neutral density distribution.

II. SIMULATION MODEL
The time evolution of a trapped single ion colliding with ultracold atoms is modelled using a simulation consisting of an analytical and a numerical part. In the time between collisions, trajectories are analytically described using the pseudo-potential approximation, while elastic collisions are taken into account using Monte Carlo techniques.
This formula, describing a three-dimensional harmonic oscillator, will be used throughout the following calculations to approximate the ion's position. The total secular energy E x + E y + E z will be referred to as the ion energy.
The motion of the ion is affected by collisions with the neutral atoms. We assume them to be instantaneous, meaning, the timescale of the collision is shorter than Ω −1 T , which is the shortest timescale of the motion of the ion. This assumption implies that every collision is sensitive to the momentary relative velocity, including the micromotion [37]. Therefore, we consider a number of effects causing contributions to the micromotion. Firstly, the intrinsic micromotion described in Eqn. (3) is proportional to the distance of the ion from the centre of the RF-quadrupole field. Secondly, static offset electric fields displace the minimum of the ion trapping potential by a distance (∆x, ∆y) from the symmetry axis of Eqn. (1). Thirdly, RF pickup on end-cap electrodes can lead to micromotion along the trap symmetry axis. These contributions are summed in the expression where c z parameterizes the micromotion along the trap symmetry axis [38]. The ion velocity considered for collisions is with v sec = d dt r sec . This is similar to the time derivative of Eqn. (3) but includes all the excess micromotion terms.

B. Collision dynamics
The simulation uses classical trajectories for the motion of the ion. Therefore, its validity is restricted to ion energies well above the energy quanta of secular motion ω. The temperature of the ultracold buffer gas is assumed to be well below this energy scale, and in the collision the neutral atom's initial energy is neglected. Due to conservation of energy and momentum, the elastic scattering process is defined by the scattering angles (θ, φ). The ion's velocity changes according to with v ion,i ( v ion,f ) being the initial (final) velocity given by Eqn. (6) at the time of the collision, β = mn mi+mn and R is the rotation matrix determined by θ and φ, with respect to the direction of v ion,i . From v sec,f and r sec a new set of ϕ j and E j can be determined, which describes the ion's trajectory after the collision.
To illustrate its impact on the motion of the ion, Eqn. (7) can be rewritten in terms of the secular velocity, yielding This expression shows how v mm couples to the secular velocity in every collision. Note that v mm is, by definition of Eqn. (5), the same before and after the collision since it only depends on the position r sec and time t of the instantaneous collision.

C. Scattering rate
The probability dP c for the ion to collide with a neutral atom within a short time interval dt defines the scattering rate It is proportional to the neutral atom density n( x) at the ion's position. The cross-section σ(E c ) is usually a function of the collision energy E c , which, neglecting the energy of the neutral atom, is given by E c = β mi 2 v 2 ion . The ion's position changes on a timescale of ω −1 j while the velocity of the ion v ion changes on a timescale of Ω −1 T . In general, Γ(t) will therefore be time dependent in a non-trivial way. The sampling method used to efficiently choose the time of collision is explained in Appendix A.
A technical description of the main simulation loop is given in Appendix B.

D. Inelastic collisions
Inelastic processes like charge exchange, spin exchange or molecule formation have been predicted to occur in the hybrid system [11,22,[39][40][41]. In experiments [6][7][8][9], charge exchange, which is typically signaled by the loss of the ion, has been observed at rates many orders of magnitude lower than the elastic collision rate, in the non-resonant case. Spin exchange collisions can occur with rates comparable to elastic scattering [22], and energy from internal states can be transferred to the external degrees of freedom. In a spin-stretched configuration, however, spin exchange is suppressed.
In our simulation, we can include inelastic effects by introducing additional, competing rates, defined as in Eqn. (9), but with inelastic cross-sections σ i (E c ) instead. The effect of any inelastic event on the hybrid system depends primarily on the question whether the original ion still exists after the process. If this is not the case, the simulation can be stopped at the first occurrence. If the ion continues to exist, the amount of energy released or absorbed by the internal states of the colliding particles needs to be considered in a modified version of Eqn. (7). In either case, our simulation is able to predict the rate at which inelastic events occur, given the inelastic crosssection σ i (E c ). On the other hand, the simulation can be used, if a rate is measured experimentally, to determine σ i (E c ).
In the following sections we consider elastic processes only, assuming inelastic processes either to happen very rarely, in line with the experimens [7][8][9], or to involve only small amounts of internal energy, which do not significantly affect the system.

III. LANGEVIN SCATTERING
For the motion of an ion in a neutral gas, mainly large angle scattering is considered relevant, as small deflections do not significantly change the ion's trajectory. This assumption leads to the Langevin scattering model, which successfully describes the ion's mobility in previous experiments with ions in a neutral buffer gas. Here, we apply the Langevin scattering model to the trapped ion system including the effects of micromotion. We investigate the properties of the ion's energy and will later compare these results with those from a more complete semiclassical scattering model, in section IV. We will indeed find good agreement between the two models in describing the energy distribution of the ion, confirming the above assumption, that large angle scattering events determine the evolution of the ion's energy.
The ion-neutral interaction is dominated by the attractive polarization interaction, which is of the form with C 4 = α 0 Q 2 /(4πǫ 0 ) 2 being proportional to the neutral particle polarizability α 0 . R is the internuclear separation. Classically one can define a critical impact pa- [42]. Collisions with impact parameter b > b c lead to small deflections and are neglected. Impact parameters b < b c result in inwardspiralling trajectories, which lead to almost uniformly distributed scattering angles into all directions as in hardsphere scattering. The resulting cross-section for large angle scattering, σ L = π b 2 c is proportional to the inverse collision velocity, and leads to a scattering rate independent of the collision energy [24].

A. Energy scale
Our aim is to determine the energy scale of the ion on its classical trajectory after many collisions such that the initial conditions for E j can be neglected. Neither the Langevin scattering at its energy-independent rate, nor the ultracold neutral buffer gas, which we assume at T = 0, introduce such a scale. As a consequence the only energy scale in the system is defined by the excess micromotion, see Eqn. (5). This is in contrast to the case of a buffer gas with non-negligible temperature, where the energy scale is rather set by the temperature of the neutral gas [31].
To associate the excess-micromotion with an energy scale we define with v mm,e being the full velocity amplitude of the micromotion for an ion in the centre of the ion trapping potential. v mm,e depends on the displacement (∆x, ∆y) caused by the uncompensated offset electric field. Note that when the offset electric field is compensated using a photon correlation measurement [43,44], photon shot noise usually limits the lowest achievable E mm,e . Since E mm,e is the only energy scale in the colliding system it is sufficient to express all energies in units of E mm,e . This gives general results for any amplitude of uncompensated micromotion. It also implies that any statistical measure of ion energy has to scale with E mm,e and therefore with the square of the excess micromotion amplitude. Such quadratic dependence of the mean ion energy has been experimentally observed in [8].

B. Energy spectrum
In order to treat scattering in the presence of micromotion in its most general form we choose all trap frequencies (ω x,y,z , Ω T ) to have irrational ratios. We assume a very low homogeneous neutral atom density such that the scattering rate Γ is much smaller than the trap frequencies. This ensures that two consecutive collisions happen at uncorrelated positions.
We obtain the energy spectrum from the simulation by binning the secular energy after each collision on a logarithmic scale. This measures a logarithmic energy probability distribution dP (E)/d log(E). Fig. 1 shows such energy spectra for three different mass ratios of the atom and the ion. For heavier neutral atoms the mean energy of the ion increases and a tail in the spectrum towards higher energies becomes dominant. Even for equal Energy spectra of the ion for three different mass ratios. 10 8 energies are sampled into logarithmically spaced bins with an energy resolution of 100 bins per decade. We choose ∆x = ∆y and cz = 0 and trap parameters q 2 a = 50. In the case where the ion mass is twice the neutral mass (black) the ion's mean energy is 0.8 Emm,e, in the case for equal masses (red) 5 Emm,e. The red dashed line corresponds to a thermal energy distribution with 5 Emm,e. The spectrum for a lighter ion, mn m i = 1.7 (blue), contains a significant contribution of very high energies, typically leading to quick ion loss due to finite trap depth. masses (β = 0.5) the tail towards high energies is evident when compared to the thermal distribution with the same mean energy. For any mass ratio, the obtained spectrum distinctly differs from a thermal distribution, also because very low energies (E ≪ E mm,e ) are extremely rare. This supports the validity of the classical treatment of trajectories and instantaneous collisions (E > Ω T ).
A power law, dP (E)/d log(E) ∝ E α nicely fits the tail in the spectrum towards high energies for m n > m i or β > 0.5 [31]. As β is increased further, towards even heavier neutrals, the negative exponent α approaches 0, at which point the spectrum does not converge with time anymore and runaway heating starts to dominate the evolution of the ion's energy. The critical mass ratio parameter β crit can be found by extrapolating the exponent α(β) towards α(β crit ) = 0. The quantity β crit is not a universal number but is a function of the trap geometry. It depends on the ratio ωp ωz or, expressed in the ion trap parameters a and q, on q 2 a = ( 2ωp ωz ) 2 . The dependency is explained by the fact that axial confinement leads to radial deconfinement and thereby to an increase of the ratio between average micromotion and secular energy. For three different cases the extrapolation to β crit is shown in Fig. 2. Our data are compatible with the critical mass ratio previously found for a specific trap geometry [31]. For very elongated traps (ω z ≪ ω p ) we find β crit = 0.685, corresponding to a mass ratio mn mi = 2.17 . The exponents α obtained from fitting power laws to the high energy tails of the ion energy spectra are plotted against the mass ratio coefficient β. This is done for three different trap geometries, for a spherical trap ( q 2 a = 6, black rectangles), for an elongated trap ( q 2 a = 50, red circles) and for the extreme case with negligible axial confinement ( q 2 a = 10 10 , blue triangles). The data is fitted with second order polynomials to extrapolate to βcrit for which the exponent α becomes 0.

C. Average energy and lifetime of the ion
To evaluate the effectiveness of buffer gas cooling for different mass ratios we calculate the average of the energy spectrum of the ion. As a physical measure, the arithmetic mean comes close to the definition of a temperature, albeit the clearly non-thermal distribution. We show the arithmetic mean in Fig. 3 for the case q 2 a = 50. Although the arithmetic mean diverges for α ≥ −1, the energy spectrum can still be normalized for α < 0. In this range the median continues to be a well defined statistical measure for the energy of the ion.
In experiments, a large probability density at high energies leads to a rapid ion loss due to a finite trap depth E td [30,31]. We have numerically evaluated the required trap depth E td to limit the ion loss probability per collision to P loss (Fig. 3). For large β we can approximate the required trap depth by E td = E mm,e (P loss ) 1/α . The results show that mass ratios with light neutrals are preferred for efficient buffer gas cooling. However, even for the heavy neutral scenario stable trapping and buffer gas cooling are possible for carefully chosen trap geometry, trap depth and micromotion compensation. The arithmetic mean is expected to diverge for β > 0.554. An example for the required trap depth for P loss < 10 −5 is given (blue triangles) and compared to E td = Emm,e (P loss ) 1/α (dotted line) using α obtained from power law fits to the tails of the energy spectra. All data are for q 2 a = 50.

D. Higher collision rates
For all the results discussed so far the collision rate has been assumed very low, and as long as the condition Γ ≪ ω is fulfilled the results remain unchanged. As the collision rate approaches or even exceeds the trap frequency ω x,y , the probability to have consecutive collisions at correlated positions increases. Under this condition we observe a reduced median energy and an increased power law tail.

IV. SEMICLASSICAL SCATTERING
Up to this point the classical Langevin model has been used to explore the ion's energy spectrum and its dependence on mass ratio, trap geometry and collision rates.
Here we make use of a semiclassical description of the interaction process. The solution to the quantum mechanical scattering problem can be found by expanding the wavefunction into partial waves. In the regime where many partial waves contribute, the total elastic cross section scales like E −1/3 c [39]. The resulting energy dependent scattering rate and angular dependence lead to additional effects as compared to the Langevin model. These are necessary to explain the back action on the neutral cloud.
We model the interaction potential by the long range polarization interaction of Eqn. (10) plus repulsion at short distances. The full differential cross-section is calculated using [45] dσ The angular momentum of a partial wave is l and k = √ 2µE c is the collision momentum. The scattering phase η l can be obtained by solving the radial Schrödinger equation which involves the centrifugal potential 2 l(l+1) 2µR 2 . The resulting centrifugal barrier increases in height with angular momentum (∼ ( l) 4 ). Partial waves with l < l 0 = 1/ 2µ √ 2C 4 E c have a collision energy larger than the height of the centrifugal barrier, probe the deep potential well and are reflected from the hard core. The exact form of the potential, relevant to determine η l , is typically not known. Therefore phase shifts η l for l < l 0 are assumed to be uniformly distributed within [0, 2π) [11]. In this approximation each partial wave contributes with σ l = 2πl k 2 to the total crosssection. Summing σ l up to l 0 reproduces the Langevin cross-section.
The full quantum mechanical cross-section includes additional contributions from partial waves with l > l 0 . As the centrifugal barrier is higher than the collision energy, these partial waves are scattered from the centrifugal barrier, if tunnelling effects are neglected. The phase shifts can be semiclassically approximated by [11,39] with R 0 = l+1/2 k .
A. Modeling the differential cross-section The probability distribution for the deflection angle θ is numerically calculated using Eqn. (12), Eqn. (13) and randomly distributed η l for l < l 0 . We sum Eqn. (12) for l up to 20000 and average over 100 different random sets of η l . The differential cross-section calculated in this way depends on the reduced mass µ, the collision energy E c = 2 k 2 2 µ and C 4 only. For the case of a 174 Yb + ion colliding with a 87 Rb atom, four probability distributions of the form Eqn. (14) are shown in Fig. 4. The main feature of I(θ, E c ) is a forward scattering peak, which gets more pronounced as the energy increases [41]. The integral of the differential cross-section reproduces the expected E −1/3 c energy dependence, and its magnitude is in agreement with [11]. The four curves are for different collision energies Ec. The forward scattering peak at small θ is more pronounced for higher energies. Its shape is emphasized in (b) where the scattering angle θ is displayed logarithmically. The normalized I(θ, Ec) (smooth lines) for the four different energies are compared with logarithmically binned output from the random θ-generating function (step like lines). The approximation is optimized to reproduce the height and position of the forward scattering peak.
To implement the differential cross-section in the Monte Carlo simulation a parameterization of the normalized I(θ, E c ) is used to create a function that returns a random θ for a given collison energy. The distribution I(θ, E c ) is modelled in four intervals using two power laws (∝ θ p ), a flat top and a flat background. The parameters peak height, background offset and interval limits are energy dependent and are well approximated by power laws (∝ E p ′ c ). These power laws are obtained from fits to differential cross-sections for more than 30 different energies E c in the range between k B ×1 µK and k B ×100 K. The θgenerating function uses these parameter functions and inverse transform sampling. Fig. 4b compares sampled output of the θ-generating function with the normalized differential cross-sections.
B. Effects of the energy dependent scattering rate on the ion energy spectrum We have simulated the kinetics of the ion in an ultracold buffer gas using the parametrized differential crosssection to investigate how this affects the ion energy spectrum. Different from the Langevin case, the collision rate depends on the instantaneous ion energy. Therefore efficient simulation relies on the collision time sampling described in Appendix A. The ion energies are binned and weighted by the time the ion remains at the specific energy.
We have performed simulations for an ion trap with frequencies ω x,y,z = 2π × {151, 153, 42} kHz, Ω T = 2π × 42.5 MHz, excess micromotion parameters ∆x = ∆y = 2 µm and neutral density n = 10 18 m −3 for the system 174 Yb + -87 Rb, reproducing conditions comparable to [8]. We find good agreement with the previous results from the Langevin model and conclude that the Langevin model is sufficient to describe the ion's energy statistics. Formally, however, the system does not necessarily scale only with the excess-micromotion energy E mm,e anymore, as the differential cross-section introduces its own energy scale. In the next section we will demonstrate that the full differential cross-section is necessary to predict effects on the cold neutral atoms.

C. Neutral cloud evolution
Ultracold neutral atomic clouds have atom numbers ranging up to 10 9 . Compared to a room temperature buffer gas, this limited number of atoms and the good isolation from the environment allow the observation of collision effects on the neutral gas. The main observables are the number of neutral atoms N a and their temperature T a . In experiments, these values can be obtained from time-of-flight imaging and are suitable to verify the simulation model.
The back-action of the ion onto the neutral gas is a result of the energy transfer per collision. For general two-body elastic collisions it is given by depending on the scattering angle θ. For very small deflections θ ≪ 1, resulting from the forward scattering peak in the differential cross-section, only very little energy is transferred to the neutral atom. The distribution of transferred energies E t , shown in Fig. 5, can be understood as a convolution of the collision energy distribution and the energy dependent differential cross-section.  . The semiclassical distribution of E c is also slightly shifted towards higher energies with respect to the Langevin E c , since collisions are more likely to happen at higher energies. The transferred energies differ only little between to two models for E t k B ×0.03 K. This reflects the statement that the ion's mobility is well described by Langevin type collisions. However, the significant peak at low E t of the semiclassical distribution causes most of the effects on the cold neutral gas.
So far all the results were obtained assuming a uniform density distribution of neutral atoms. This will now be replaced by a spatial distribution for a thermal gas in a harmonic trap with temperature T a and atom number N a . Considering a finite trap depth E td,a for the neutral atoms, every collision with E t > E td,a will lead to an atom loss, whereas every E t < E td,a will increase the temperature T a . The simplified model used in the simulation assumes immediate thermal equilibration of the neutral gas. Then a collision with E t < E td,a simply increases T a by Et 3 kB Na . The loss of an atom for E t > E td,a will decrease N a by 1 but also affect the temperature de-pending on the atoms energy. The new temperature is calculated using the total energy of the neutral cloud, and the potential and average kinetic energy of the lost atom known from the position r of the collision. This can lead to evaporative cooling or heating effects depending on the position of the ion in the neutral gas. For ion trajectories larger than the size of the buffer gas the ion can only collide in the centre of the trap where its motional energy is mostly related to the secular motion rather than the micromotion. This suppresses the power law tail of the ion energy distribution and reduces the ion's average energy. Hence, using tight traps to confine the neutral atoms might help to overcome the constraints on ion trap depth, trap geometry and mass ratio. Here we compare the simulation predictions for the effects on neutral atoms with experimental data from Yb + − Rb [8]. The measured quantities are the loss of neutral atoms and temperature increase of a cold thermal cloud, after 8 s of interaction and for different excessmicromotion energies E mm,e . Fig. 6 displays the data together with the simulation results.
Initial conditions for the neutral 87 Rb cloud in the simulation are given by T a,0 = 250 nK and N a,0 = 2.25×10 6 . Neutral trap frequencies of 2π × {28, 28, 8} Hz result in an initial central density of n(0) = 1.9 × 10 18 m −3 . The neutral trap depth is E td,a = k B × 8 µK.
A single 172 Yb + ion is trapped with parameters given in section IV B. The excess micromotion parameters ∆x and ∆y are varied between 0 and 15 µm, and additionally there is micromotion along the trap symmetry axis with c z = 2 µm. The excess-micromotion energy scale, E mm,e /k B thus ranges from 90 mK to 7 K.
The semiclassical simulation predicts both the shape and the magnitude of the experimental results well. In contrast, the Langevin scattering is not suited to model these measurements because the ultracold atoms are highly sensitive to small energy transfer E t . They correspond to collisions with small deflection angles, which are neglected by the Langevin model, cf. Fig. 5.

V. CONCLUSION
We have numerically investigated the kinetics of a single trapped ion interacting with an ultracold neutral gas. Our results explain the effects of the mass ratio, trap geometry and excess micromotion on the ion's energy spectrum. We have applied two different collision models of the atom-ion interaction, the Langevin and the semiclassical scattering model. Both yield similar ion energy statistics. The greater simplicity of the Langevin model introduces a characteristic energy scale for the ion energy statistics. However, experimentally observed effects on the neutral cloud can only be explained by the semiclassical model. Forward scattering events with small energy transfer are affecting both the neutral cloud temperature and the atom loss rate. Cold atom-ion collisions could then be used to locally remove atoms resulting in efficient cooling of quantum gases.

Appendix A: Collision time sampling
A scattering rate Γ sets the probability for a collision to take place within a time interval dt. In our case with n being the neutral atom density, σ the cross-section and the velocity v ion as defined in Eqn. (6). Usually the density is a function on the ion's position, the crosssection depends on the collision energy and v ion oscillates rapidly in time, leading to a Γ(t) with non-trivial time dependence. In the following we explain the method used to randomly generate collision times t with the exact distribution defined by Γ(t).
In general, the process of an event (collision) taking place with a rate Γ(t) can be modelled using a differential equation for the probability Q(t) that the event has not yet happened after the time t, The probability distribution for an event to take place after the time t is defined by P (t) = −dQ/dt. In the simple case with constant Γ(t) = Γ 0 the solution is and a random time can be obtained using inverse transform sampling, with r being a uniformly distributed random number in the interval (0, 1]. For time dependent Γ(t) the analytic solution for the probability distribution function is Non-trivial time dependence of Γ(t) usually requires a numerical approach to sample t from Eqn. (A5). One straight-forward method would be to discretize time into small steps, calculate Γ(t) and its contribution to P (t) for every step, thereby numerically integrating the function P (t) up to a randomly chosen trigger value, at which point the event takes place. However, another method is employed here, which proves to be much faster and does not suffer from discretization errors. It works for Γ(t) that have an upper bound Γ m = sup(Γ(t)), or where such a condition can be enforced by introducing a cutoff. In the specific case of the trapped ion, all factors in Eqn. (A1) are easily limited by considering the energies E j defining the trajectory and the excess micromotion parameters and using the peak neutral density. This upper bound can be adjusted after each collision event having affected E j and n( x).
The algorithm works by advancing the system by a time t according to Eqn. (A4) with Γ 0 = Γ m . Then the rescaled rate is calculated for the resulting state of the system after the time t and an event takes place if γ(t) > r, with r being another uniformly distributed random number in the interval [0, 1). If the event does not take place (γ(t) ≤ r) the algorithm iteratively loops back to advance the system by an additional t, again according to Eqn. (A4). The method is exact in that it reproduces the probability distribution function Eqn. (A5). A proof follows below. The efficiency of the method is the ratio of the average of Γ(t) and Γ m , ǫ = Γ(t) /Γ m = γ(t) . This means that in order to generate N events it can be expected that Γ(t) needs to be evaluated N/ǫ times. Proof: We start from writing an expression for the probability distribution function P s (t) obtained with the suggested method. Since the final t can be the result of any number of iterations, P s (t) is a sum of all these possibilities, P s (t) = P s,1 + P s,2 + P s,3 + . . ., where P s,i is the probability that t results as the time of event after i iterations. The first few terms are given below.
Note that the products of P Γm (as defined in Eqn. (A3)) in the integrals always combine to Γ (i−1) m P Γm (t). Therefore P Γm (t) is taken out of the sum as a common prefactor, Now the upper boundaries of all the partial integrals can be set equal to t, since they only induce ordering to the time series {t 1 , t 2 , . . . , t n }. This rescales the terms by the number of possible orderings (n!). Then the partial integrals reduce to a single integral to the power of n, The following is a short description of how the main simulation loop has been implemented. 1) setup system configuration: ω j , Ω T and excess micromotion, β, neutral trap frequencies, T a , N a , n( x), dσ dΩ , ...

11) loop back to 3)
The setup of the system configuration in 1) contains mainly parameters which do not change during the collisions, such as trap frequencies or the atom-ion mass ratio. Exceptions are the neutral atom number N a and temperature T a , which act as initial conditions. The initial state for the ion energy in 2) is mostly unimportant, as the simulation will iterate over many collisions and the information of the initial state is lost after a few collisions. When looking at steady state statistics of the ion energy, the values after the first few collisions can simply be ignored, thus effectively letting the system evolve for a short time before measuring its properties. The points 3) to 7) implement the collision time sampling algorithm described in Appendix A. The maximally possible collision rate Γ m is calculated from the peak density n max (T a , N a ) of the neutral atoms and on the maximum value of σ(E c ) v ion (t). The latter is limited by the highest possible ion velocity, which depends on the secular energy E j and the excess micromotion. The calculation of γ(t) in 6) uses the v ion and r ion obtained in 5) to determine Γ(t) with Eqn. (9). If a collision takes place at the chosen time t, 7) gives the scattering angles according to the differential cross-section. 8) simulates the back-action on the neutral atoms. 9) modifies the velocity of the ion. In 10) the new trajectory parameters E j and ϕ j are determined, to represent the motion of the ion up to the next collision. Information on any system parameter can be retrieved at user-defined points within the simulation loop. The ion energy is typically registered after 10), the collision and transferred energies, E c and E t after 8). Random sampling is used in 4) and 6), to determine the time of collision, and in 7), for the scattering angles.