Superstatistical energy distributions of an ion in an ultracold buffer gas

TL;DR

This paper develops a superstatistical framework to analytically describe the non-thermal energy distributions of ions in RF traps interacting with ultracold buffer gases, revealing how Tsallis statistics naturally emerge from physical parameters.

Contribution

It introduces a first-principles superstatistical model for ion energy distributions in RF traps with ultracold gases, linking Tsallis statistics to controllable experimental parameters.

Findings

Analytic ion energy distributions derived from first principles.

Tsallis statistics emerge from physical mechanisms involving heating and fluctuations.

Distributions depend on experimentally controllable parameters.

Abstract

An ion in a radiofrequency ion trap interacting with a buffer gas of ultracold neutral atoms is a driven dynamical system which has been found to develop a non-thermal energy distribution with a power law tail. The exact analytical form of this distribution is unknown, but has often been represented empirically by q-exponential (Tsallis) functions. Based on the concepts of superstatistics, we introduce a framework for the statistical mechanics of an ion trapped in an RF field subject to collisions with a buffer gas. We derive analytic ion secular energy distributions from first principles both neglecting and including the effects of the thermal energy of the buffer gas. For a buffer gas with a finite temperature, we prove that Tsallis statistics emerges from the combination of a constant heating term and multiplicative energy fluctuations. We show that the resulting distributions essentially depend on experimentally controllable parameters paving the way for an accurate control of the statistical properties of ion-atom hybrid systems.

Equations

$$e_q(x)=(1+(1-q_T)x{)}^{\frac{1}{1-q_T}}$$

$${\ddot{r}}_j(\tau )+[a_j-2q_j\cos (2\tau )]r_j=0,$$

$${\mathbf{v}}^{\prime }=\frac{1}{1+\tilde{m}}\mathbf{v}+\frac{\tilde{m}}{1+\tilde{m}}{\mathbf{v}}_{\mathbf{n}}+\frac{\tilde{m}}{1+\tilde{m}}\mathrm{R}\cdot (\mathbf{v}-{\mathbf{v}}_{\mathbf{n}})$$

$$E^{\prime }=\eta E+c_1\sqrt{Eϵ}+c_2ϵ,$$

$$f_{\eta }(\eta )=\frac{1}{\delta }\frac{a_1{a}_2}{a_1+a_2}\{\begin{array}{ll}{\left(\frac{\delta }{\eta }\right)}^{a_1+1}& \eta \ge \delta \\ {\left(\frac{\eta }{\delta }\right)}^{a_2-1}& 0<\eta <\delta \end{array}$$

$$f_{E_0}(E_0)=\frac{E_0^k{\beta }_0^{k+1}}{\Gamma (k+1)}{e}^{-E_0{\beta }_0},$$

$$f_{E^{\prime }}(E^{\prime })$$

$$={\int }_{\eta =0}^{\eta =\infty }\frac{1}{\eta }\frac{(E^{\prime }/\eta {)}^k{\beta }_0^{k+1}}{\Gamma (k+1)}{e}^{-(E^{\prime }/\eta ){\beta }_0}{f}_{\eta }(\eta )\mathrm{d}\eta .$$

$$f_{\eta }(\eta )=\delta (\eta -{\eta }_c),$$

$$f_{E^{\prime }}(E^{\prime })=\frac{E^{\prime k}{\beta }_0^{k+1}}{{\eta }_c^{k+1}\Gamma (k+1)}{e}^{-\frac{E^{\prime }{\beta }_0}{{\eta }_c}}.$$

$$f_{E^{\prime }}(E^{\prime })={\int }_{{\beta }^{\prime }=0}^{{\beta }^{\prime }=\infty }\frac{E^{\prime k}{\beta }^{\prime k+1}}{\Gamma (k+1)}{e}^{-E^{\prime }{\beta }^{\prime }}\frac{{\beta }_0}{{\beta }^{\prime 2}}{f}_{\eta }\left(\frac{{\beta }_0}{{\beta }^{\prime }}\right)\mathrm{d}{\beta }^{\prime }.$$

$$f_{{\beta }^{\prime }}({\beta }^{\prime })=\frac{{\beta }_0}{{\beta }^{\prime 2}}{f}_{\eta }\left(\frac{{\beta }_0}{{\beta }^{\prime }}\right),$$

$$f_{E^{\prime }}(E^{\prime })={\int }_{{\beta }^{\prime }=0}^{{\beta }^{\prime }=\infty }\frac{E^{\prime k}{\beta }^{\prime k+1}}{\Gamma (k+1)}{e}^{-E^{\prime }{\beta }^{\prime }}{f}_{{\beta }^{\prime }}({\beta }^{\prime })\mathrm{d}{\beta }^{\prime }.$$

$$f_{{\beta }^{\prime }}({\beta }^{\prime })={\int }_{\eta =0}^{\eta =\infty }\eta f_{\beta }(\eta {\beta }^{\prime })f_{\eta }(\eta )d\eta .$$

$${\beta }_i={\beta }_{i-1}/{\eta }_i.$$

$${\beta }_n={\beta }_0\prod _{i=1}^{n}1/{\eta }_i.$$

$$f_{{\beta }_n}({\beta }_n)=\frac{1}{\sqrt{2\pi n}\sigma {\beta }_n}\exp [-\frac{(\ln {\beta }_n-\ln {\beta }_0+n\mu {)}^2}{2n{\sigma }^2}],$$

$$f_{E_n}(E_n)$$

$$\times \frac{1}{\sqrt{2\pi n}\sigma {\beta }_n}\exp [-\frac{(\ln {\beta }_n-\ln {\beta }_0+n\mu {)}^2}{2n{\sigma }^2}]\mathrm{d}{\beta }_n.$$

$$f_{E_n}(E_n)=\frac{{\hat{\beta }}^3{E}_n^2}{4\sqrt{\hat{\beta }{E}_{n}n{\sigma }^2+1}}\exp \left(-\hat{\beta }{E}_n\right)$$

$$\times \left(\text{erf}\left(\sqrt{\frac{\hat{\beta }{E}_{n}n{\sigma }^2+1}{2n{\sigma }^2}}\right)+1\right)\exp \left(-\frac{n{\sigma }^2}{2}{\left(\hat{\beta }{E}_n-2\right)}^2\right),$$

$$\mathcal{L}[{\beta }^{k+1}{f}_{\beta }(\beta )]=(-1{)}^{k+1}\frac{d^{k+1}}{d{E}^{k+1}}\mathcal{L}[f_{\beta }(\beta )],$$

$$T_i={\eta }_i{T}_{i-1}+\kappa T_a,$$

$$f_{\beta }(\beta )=\frac{1}{\beta \Gamma (n_T)}{e}^{-\frac{\beta n_T}{⟨\beta ⟩}}{\left(\frac{\beta n_T}{⟨\beta ⟩}\right)}^{n_T}.$$

$$f_{E,T}(E)={\left(\frac{n_T}{⟨\beta ⟩}\right)}^{-k-1}\frac{\Gamma (k+n_T+1)}{\Gamma (k+1)\Gamma (n_T)}\frac{E^k}{{\left(\frac{⟨\beta ⟩E}{n_T}+1\right)}^{k+n_T+1}}.$$

$${\int }_{\eta =0}^{\eta =\infty }{f}_{\eta }(\eta ){\eta }^{n_T}d\eta =1.$$

$$n_T=a_1-a_2=\frac{⟨\eta ⟩-4⟨{\eta }^2⟩+3⟨\eta ⟩⟨{\eta }^2⟩}{⟨\eta ⟩-2⟨{\eta }^2⟩+⟨\eta ⟩⟨{\eta }^2⟩}.$$

$$⟨T⟩=⟨\eta ⟩⟨T⟩+\kappa T_a=\frac{\kappa T_a}{1-⟨\eta ⟩},$$

$$⟨T⟩=\frac{1}{k_B⟨\beta ⟩}\frac{n_T}{n_T-1}.$$

$$⟨\beta ⟩=\frac{1}{k_B\kappa T_a}\frac{n_T}{n_T-1}(1-⟨\eta ⟩).$$

$${\ddot{r}}_j(\tau )+[a_j-2q_j\cos (2\tau )]r_j=0,$$

$$r_j(\tau )=A_j\cos ({\varphi }_{j0})\mathrm{ce}(a_j,q_j,\tau )-A_j\sin ({\varphi }_{j0})\mathrm{se}(a_j,q_j,\tau ),$$

Full text

Ion-atom hybrid systems are realized by superimposing cold ions in a radiofrequency (RF) trap with trapped ultracold atoms Härter and Hecker Denschlag (2014); Sias and Köhl (2014); Willitsch (2015). RF traps use rapidly oscillating electric fields to dynamically confine the ions. In an adiabatic regime Gerlich (1992), the resulting motion of an ion can be represented as a thermal component (“secular motion”) superimposed by small-amplitude oscillations at the RF frequency (“micromotion”) Gerlich (1992). In a hybrid trap, the ion undergoes frequent collisions with neutral atoms which disrupt its motion and lead to energy exchange between the secular motion and the RF field Major and Dehmelt (1968); Moriwaki et al. (1992); DeVoe (2009); Zipkes et al. (2011); Chen et al. (2014); Höltkemeier et al. (2016).

These processes lead to a distortion of the ion’s secular-energy distribution from thermal (Boltzmann) to one better described by a power law at high energy Levy and Solomon (1996); DeVoe (2009); Zipkes et al. (2011); Chen et al. (2014). The precise knowledge of the ion energetics is crucial for understanding the properties and dynamics of hybrid systems and their derived applications. Consequently, this problem has been the subject of intense recent research DeVoe (2009); Zipkes et al. (2011); Chen et al. (2014); Höltkemeier et al. (2016, 2016). In the high-energy limit, expressions for the mean energy and the power-law exponent have been derived Chen et al. (2014). The complete ion energy distribution has often been modeled DeVoe (2009); Biesheuvel et al. (2016); Meir et al. (2016a) by Tsallis (q-exponential) functions Tsallis (1988, 2009),

[TABLE]

for $q_{T}>1$ where $C_{q}$ is a normalization factor. $q_{T}$ is a parameter which characterizes the deviation from a standard exponential function which is recovered in the limit $q_{T}\rightarrow 1$. However, the application of q-exponentials has remained empirical Biesheuvel et al. (2016); Höltkemeier et al. (2016); Meir et al. (2016b) since their first introduction for fitting numerical energy distributions DeVoe (2009).

Using the formalism of superstatistics Beck (2001); Beck and Cohen (2003), we introduce a framework for the statistical mechanics of ion-atom hybrid systems. We derive analytic ion secular-energy distributions both neglecting and including the thermal energy of the ultracold buffer gas and confirm their validity by comparison with numerical simulations. For a buffer gas at zero Kelvin, we obtain an energy distribution with no steady-state and an exponential decay at high energies. For a buffer gas at finite temperature, we prove from first principles the emergence of Tsallis statistics thus vindicating its application in the present context. The energy distributions derived here depend on experimentally adjustable parameters which opens the door for a rational experimental control of the statistical properties of ion-atom hybrid systems.

The motion in each direction $r_{j},j\in(x,y,z)$, of an ion in a quadrupole RF trap is given by the Mathieu differential equations,

[TABLE]

where $\tau=\Omega t/2$ and $q_{j},a_{j}$ are the Mathieu stability parameters Major et al. (2005). In a hybrid system, the ion interacts with neutral atoms through a polarization potential. We treat the dynamics as a series of elastic collisions in the Langevin approximation with an energy-independent rate Gioumousis and Stevenson (1958); Chen et al. (2014). The velocity $\mathbf{v}$ of the ion after a collision is Zipkes et al. (2011); Chen et al. (2014),

[TABLE]

where $\tilde{m}=m_{n}/m_{i}$ is the ratio of the atom’s to the ion’s mass, $\mathrm{R}$ is a rotation matrix Arvo (1992) and $\mathbf{v_{n}}$ is the velocity of the neutral atom. Primes refer to post-collision quantities. As we have no control over the instantaneous velocities of the particles at the time of collisions, $\mathbf{v}$ and $\mathbf{v_{n}}$ are random variables.

From Eq. (3), the ion’s secular energy after a collision can be derived to be (see supplemental material),

[TABLE]

where $\epsilon=\frac{m_{n}}{2}|\mathbf{v_{n}}|^{2}$ is the kinetic energy of the neutral atom and $\eta,c_{1},c_{2}$ are coefficients (see supplemental material). Assuming that the buffer gas density is uniform, these coefficients are independent of the values of $E$ and $\epsilon$. For an ion much hotter than the buffer gas $(E\gg\epsilon)$, we approximate $E^{\prime}\approx\eta E$. The stability of the ion motion in the buffer gas with respect to runaway heating is determined by the distribution of $\eta$. As a rule, the motion is stable for a mass ratio $\tilde{m}\lesssim 1.4$ for Mathieu parameters $q\ll 1$ Moriwaki et al. (1992); DeVoe (2009); Zipkes et al. (2011); Chen et al. (2014); Höltkemeier et al. (2016).

Fig. 1 shows numerical distributions $f_{\eta}(\eta)$ for the energy-transfer parameter $\eta$ for $q=0.1,\tilde{m}=0.75$ and $q=0.5,\tilde{m}=1.25$ which correspond to stable and unstable ion motions, respectively. The numerical simulations were performed following DeVoe’s approach DeVoe (2009).

The solid lines in Fig. 1 correspond to log-Laplace distributions of the form Kozubowski and Podgórski (2003),

[TABLE]

with $a_{1},a_{2}>0$ which have previously been used as to model processes involving multiplicative fluctuations Kozubowski and Podgórski (2003). The parameter $\delta$ representing the maximum of the distribution was found to be $\approx 1$, reflecting the fact that most collisions result in little changes to the ion’s energy. The values of $a_{1}$ and $a_{2}$ may be estimated by calculating the first and second moment $\langle\eta\rangle$ and $\langle\eta^{2}\rangle$, respectively, of the distribution using Eq. (5) and matching them to the expressions found numerically from Eq. (4).

Let us now assume that the ion is initially prepared in a thermal state at temperature $T_{0}$, as may be the situation after Doppler laser cooling Metcalf and van der Straten (1999); Honerkamp (2002). The resulting distribution for the ion’s initial energy $E_{0}$ is,

[TABLE]

where $\beta_{0}=1/(k_{B}T_{0})$, $\Gamma$ is the Gamma function and the pre-exponential factor represents the density of states ($k=2$ for a three-dimensional harmonic oscillator Pethick and Smith (2001)).

We now consider the effects of collisions with the neutral atoms. Initially, we neglect their thermal energy and set $\epsilon=0$ in Eq. (4) such that $E^{\prime}=\eta E$. The resulting energy distribution can be written as Riley et al. (2010),

[TABLE]

We first consider the case in which every collision multiplies the energy by a fixed amount, $\eta_{c}$. The distribution for $\eta$ is then given by a Dirac $\delta$ function,

[TABLE]

so that

[TABLE]

This is still a thermal distribution, except that it can now be written in terms of $\beta^{\prime}=\beta_{0}/\eta_{c}$.

We now generalise this approach to an arbitrary $f_{\eta}(\eta)$ by making the change of variables $\beta^{\prime}=\beta_{0}/\eta$ in Eq. (Superstatistical energy distributions of an ion in an ultracold buffer gas),

[TABLE]

The energy distribution after a collision can thus be represented by a superposition of thermal states. This problem can be treated within the formalism of superstatistics, i.e., the superpositions of several statistics as in our case the ones of $\eta$ and $E$ in Eq. (Superstatistical energy distributions of an ion in an ultracold buffer gas) Beck (2001); Beck and Cohen (2003); Rouse and Willitsch (2015).

We can now define a distribution for $\beta^{\prime}$,

[TABLE]

which is used to recast Eq. (10) into the form

[TABLE]

Eq. (12) has the form of a Laplace transform $\mathcal{L}$. For general distributions $f_{\beta}(\beta),f_{\eta}(\eta)$ one gets

[TABLE]

Repeated application of Eq. (13) and substitution into Eq. (12) can then be performed to obtain the energy distribution of an ion after $n$ collisions.

Thus, we formulate a recurrence relation for $\beta$ after collision number $i$,

[TABLE]

Since the ion is initially in a thermal state, we take $\beta_{0}$ to be constant. After $n$ collisions starting from $\beta_{0}$, we get

[TABLE]

Each value of $\eta$ is assumed to be independently and identically distributed, and so by applying the central limit theorem the product $\prod_{i=1}^{n}1/\eta_{i}$ is log-normally distributed for large $n$ Riley et al. (2010). Hence, from Eq. (11) we write,

[TABLE]

where $\mu=\langle\ln\eta\rangle$ and $\sigma^{2}=\langle(\ln\eta)^{2}\rangle-\langle\ln\eta\rangle^{2}$.

We now return to the energy distribution. By inserting Eq. (16) into Eq. (12), we obtain,

[TABLE]

We use the Laplace integration method Asmussen et al. (2016) to find an approximate analytical solution for $k=2$. We obtain

[TABLE]

where $\hat{\beta}$ is the point at which the integrand of Eq. (B.3) is maximal. In the high-energy limit for $k=0$, Eq. (B.3) has been shown to exhibit an exponential decay Touchette and Beck (2005); Rabassa and Beck (2014). From the general property of the Laplace transform,

[TABLE]

follows that if the high-energy behavior for $k=0$ is an exponential decay, then this holds true for any integer value of $k$. Thus, we conclude that a purely multiplicative model of the heating process does not lead to Tsallis statistics which is characterized by a power-law tail for the distribution at high energies.

In order to test the validity of Eq. (B.3), a series of simulations were performed at a buffer gas temperature $T=0$ K and varying the mass ratio or number of collisions. The results are plotted in Fig. 2 along with the distributions computed from Eq. (B.3). The $\mu$ and $\sigma$ parameters were computed from numerical distributions $f_{\eta}(\eta)$ such as the ones shown in Fig. 1. At low collision numbers, the agreement is generally poor, which is expected due to the assumption in the derivation of Eq. (B.3) that the central limit theorem can be applied. Moreover, for all collision numbers, the agreement is less good at low energies due to the Laplace integration method being valid only in the limit $E\rightarrow\infty$. However, for higher energies and numbers of collisions, Eq. (B.3) becomes an increasingly better representation of the simulated data.

For comparison, the numerical data for 25 collisions at a mass ratio of 1.0 is presented in Fig. 3 together with the distribution predicted using Eq. (B.3). The red dashed line represents a Tsallis distribution obtained from a maximum-likelihood estimation (MLE) to the numerical data. It can be clearly seen that Tsallis statistics is a poor match for a buffer gas at zero Kelvin, while Eq. (B.3) provides much better agreement.

Neither the energy nor the $\beta$ distributions, Eq. (B.3) and Eq. (16), respectively, converge to a steady state with increasing $n$. This is a known property of an unbounded multiplicative random walk, and in the present case results from the neglect of the temperature of the buffer gas allowing the ion to reach an arbitrarily low temperature Levy and Solomon (1996); Sornette and Cont (1997).

For a buffer gas at a finite temperature, we have to adopt a different procedure as the change of the ion energy following a collision is no longer a purely multiplicative process, see Eq. (4). Because the buffer gas velocity distribution is isotropic, the $c_{1}$ coefficient averages to zero such that it can be neglected. We are thus left with $E^{\prime}=\eta E+c_{2}\epsilon$. Assuming again that the ion’s energy distribution can be represented as a superposition of thermal states as in Eq. (10), it follows that $\langle E\rangle=(1+k)k_{B}\langle T\rangle$. This suggests that we can rephrase the problem of finding an energy distribution to one of finding the underlying temperature distribution. We approximate that the contributions from $\epsilon$ in Eq. (4) can be treated as a constant source of heating proportional to the temperature of the buffer gas $T_{a}$ which ensures that the ion’s steady-state temperature is non-zero. This is a good approximation if the thermal fluctuations of the buffer gas are much smaller than the ones of the ion. The ion temperature after a collision is then,

[TABLE]

where $\kappa$ is a heating coefficient (see supplemental material).

To find the required temperature distribution, we solve the recurrence relation Eq. (20). The mathematical solution of this problem has been outlined in Refs. Sornette and Cont (1997); Biró and Jakovác (2005) and leads to a gamma distribution for $\beta$:

[TABLE]

Multiplying by the density of states and applying the Laplace transform we obtain the ion energy distribution,

[TABLE]

The parameter $n_{T}$ can be obtained from the condition Sornette and Cont (1997)

[TABLE]

This integral may be solved numerically, or alternatively we make use of the empirical distribution Eq. (5). From substituting Eq. (5) into Eq. (23), we obtain,

[TABLE]

assuming $\delta=1$ in Eq. (5). To fully characterize Eq. (22), we also require the value for $\langle\beta\rangle$. From Eq. (20), it follows that

[TABLE]

Averaging $T=1/(k_{B}\beta)$ over Eq. (21), we get,

[TABLE]

Equating Eqs. (25) and (26) we find,

[TABLE]

This derivation is only valid for $n_{T}>1$ and $\langle\eta\rangle<1$. If either of these conditions is not met, the mean temperature diverges because the ion motion becomes unstable.

The distribution Eq. (22) has the form of a q-exponential Eq. (1) multiplied by a $E^{k}$ term. For $k=0$ (one-dimensional), it reduces to the standard q-exponential, and for $k=2$ (three-dimensional) it is equivalent to the form used in Ref. Meir et al. (2016b), if we set their exponent $n=n_{T}+3$. We have therefore shown that Tsallis statistics are physically meaningful for the present problem under the condition that the variance of the thermal fluctuations are sufficiently small so that the additive noise due to the thermal energy of the atoms can be approximated as a constant.

Fig. 4 shows a comparison of the MLE values of the parameters $n_{T}$ and $1/\langle\beta\rangle$ extracted from numerical simulations with their predictions from Eqs. (24) and (27), respectively. Below the critical mass ratio given by the intersection of the curves with the grey horizontal line in Fig. 4 (a), the ion motion is stable. Up to near this point, the predictions for both parameters are very close to the values extracted from numerical data, vindicating the assumptions leading to the derivation of Eq. (21). Above the critical mass ratio, the predicted mean $\langle\beta\rangle$ becomes increasingly inaccurate as a result of energy correlations between different coordinate axes not accounted for in the present model (see Chen et al. (2014) and supplemental material).

From Eq. (22) and (23), it becomes clear that the energy distribution and therefore the statistical properties of the ion depend on the buffer gas temperature and the distribution $f_{\eta}(\eta)$. The latter depends on system parameters such as the atom-ion mass ratio and the Mathieu parameters of the trap which are defined in advance by the experimenter. By varying these parameters, $f_{\eta}(\eta)$ and therefore the Tsallis distribution Eq. (22) can be tuned in a deterministic manner allowing for a control of the statistical properties of the system.

Beyond the current application, the formalism developed here represents a general framework for describing the statistical mechanics of an ion in a buffer gas which can be used to, e.g., compute thermodynamic functions Beck (2011). The present treatment can also be extended to localized buffer gases. These developments will be reported elsewhere.

Appendix A Supplemental Material

A.1 Ion motion

The treatment of the motion of an ion in a quadrupole RF trap is given in standard texts, e.g., Ref. Major et al. (2005). For ready reference, we repeat here the salient points. The equation of motion in each direction $r_{j},j\in(x,y,z)$, of an ion in a quadrupole RF trap is given by the Mathieu differential equation ,

[TABLE]

where $\tau=\Omega t/2$ and $q_{j},a_{j}$ are the Mathieu parameters for axis $j$. For stable ion trajectories, the solution can be written as a sum of even (ce) and odd (se) Mathieu functions,

[TABLE]

with amplitude $A_{j}$ and initial phase $\phi_{j0}$. The velocity of the ion is obtained to be

[TABLE]

with,

[TABLE]

and the equivalent result for $\dot{\mathrm{se}}$. Here, the $c_{2m,j}$ are the coefficients of a Fourier expansion of ce and $\beta_{j}\approx\sqrt{(a_{j}+q_{j}^{2}/2)}$ for small $q_{j}$Major et al. (2005). By considering only the lowest order term $m=0$, we obtain the secular position of the ion as a function of $t$,

[TABLE]

which represents a harmonic oscillation of amplitude $\tilde{A}_{j}=A_{j}c_{0}$, frequency $\omega_{j}=\beta_{j}\Omega/2$ and initial phase $\phi_{j0}$. The secular velocity can be found through differentiation of the secular position, and from these we can define the secular energy,

[TABLE]

We define the total secular energy of the ion to be,

[TABLE]

which is the energy of a 3D harmonic oscillator with frequencies $\omega_{x,y,z}$.

Appendix B Derivation of post-collision ion energy

Here, we outline a proof that the ion’s energy after a micromotion-interrupting collision is given by Eq. (4) in the main text, i.e., $E^{\prime}=\eta E+c_{1}\sqrt{E\epsilon}+c_{2}\epsilon$, and a method to derive the coefficients $\eta,c_{1},c_{2}$ in this expression. We shall also discuss how the random variables involved may be generated so that these coefficients may be sampled for comparisons with numerical simulations.

The velocity after a collision is given by Eq. (3) of the main text, i.e.,

[TABLE]

In terms of the amplitudes $A_{j}$ and the initial secular phases $\phi_{j0}$, the components of the velocity are given by,

[TABLE]

with the equivalent for the post-collision velocity, where we use $A_{j}^{\prime},\phi_{j0}^{\prime}$ to indicate the post-collision values. To proceed, we eliminate $\phi_{j0}^{\prime}$. We approximate that the ion does not move during the collision and so by equating the positions before and after we can find an expression for $\phi_{j0}^{\prime}$. Thus, we take,

[TABLE]

then expand the Mathieu functions in a Fourier series, and apply standard trigonometric addition formulae to produce,

[TABLE]

Next, we apply the harmonic addition theorem to rewrite this as a single trigonometric function Oo and Gan (2012),

[TABLE]

where,

[TABLE]

and

[TABLE]

By equating the position before and after the collision we find,

[TABLE]

Dividing by $C_{j}$ and rearranging for the post-collision secular phase gives,

[TABLE]

We may now substitute this into Eq. (36) to find the post-collision velocities as a function of $A_{j}^{\prime},A_{j},\phi_{j0}$,

[TABLE]

where

[TABLE]

This result is substituted into Eq. (35) for the components of $\mathbf{v^{\prime}}$, and likewise Eq. (36) is substituted in for the components of $\mathbf{v}$. Thus we obtain a vector equation linking the post-collision amplitudes of motion to the pre-collision amplitude and secular phases, the neutral velocity components, and the random rotation matrix $\mathrm{R}$ in Eq. (35). The equation can then be simplified such that each component of the left-hand side is given by $s_{j}$ and so we write,

[TABLE]

where $\mathbf{s}^{T}=(s_{x},s_{y},s_{z})$ and $\mathbf{v_{R}}$ contains the remaining terms, none of which depend on $A_{j}^{\prime}$. Next, we must convert this to a set of equations for the secular energy. These energies are related to the amplitudes by Eq. (33); by inspection we will need to square the equations in order to extract the secular energy, since they are presently linear in the amplitude. We take the outer product of $\mathbf{s}$ with itself resulting in,

[TABLE]

The diagonal terms of Eq. (47) are given by,

[TABLE]

where the proportionality factor can be found from Eq. (33). By equating these terms to the corresponding term in the matrix given by $\mathbf{v_{R}}\otimes\mathbf{v_{R}}$ we obtain an expression for $\tilde{E_{j}^{\prime}}$ in terms of the amplitudes and secular phases for each axis before the collision, the elements of the rotation matrix, $\tau$ and the velocity of the neutral atom. We may use these to find the steady-state energy by averaging over the collision parameters (e.g. $x_{1},\phi_{j},\tau$) and requiring that in the steady state $\langle\tilde{E_{j}^{\prime}}\rangle=\langle\tilde{E_{j}}\rangle$ to generate three simultaneous equations of the general form,

[TABLE]

where $\alpha E_{n}$ is the fraction of the neutral energy transferred to the ion, and the $\eta_{ij}$ are coefficients describing the combination of the transfer of energy between the axes and the random fluctuation of the energy due to micromotion interruption. The set of three equations then may then be solved to find the three mean steady-state energy components $\langle\tilde{E_{j}}\rangle$. This approach leads to essentially the same set of equations as is found in Ref. Chen et al. (2014), but in terms of the secular rather than the time-averaged energies. Finally, we note that the off-diagonal elements may be averaged over in the same manner to determine quantities such as $\langle A_{x}A_{z}\rangle$ which contain information about the correlations between the motion along each axis and may be of use in future investigations.

Instead of directly averaging over these quantities to obtain the mean energies, we may also use them to prove that the form of Eq. (4) in the main text is correct, and to extract the multiplicative coefficient $\eta$. Returning to the set of three equations defined by the diagonal elements of the matrices, we convert the system to a form of spherical coordinates defined by,

[TABLE]

with the the two angles defined in the interval $[0,\pi/2)$, and an equivalent transformation applied to the primed quantities. The advantage of this coordinate system is that it simplifies the factoring of the total energy given by $E=\frac{m_{i}}{2}\rho^{2}$ from the expressions, since each secular energy component $\tilde{E_{j}}$ is proportional to $\rho^{2}$ and a function of $\phi_{\rho},\theta_{\rho}$. For $v_{n}$, standard spherical coordinates may be used,

[TABLE]

with $\epsilon=\frac{m_{n}}{2}|v_{n}|^{2}$. Performing the conversion to spherical coordinates is then followed by rearranging each equation such that the terms proportional to $\rho^{\prime}$ are on the left hand side, and summing the three equations together. This produces an equation of the form,

[TABLE]

In terms of the energies this is,

[TABLE]

where the coefficients have been redefined to include the factors of $m_{i}$ and $m_{n}$ to match the result given in Eq. (4) of the main text.

The factor $\eta$ is a function of nine random variables – the three initial secular phases $\phi_{0,j}$, the time $\tau$, the three variables $x_{1},x_{2},x_{3}$ used in the random rotation matrix given in Ref. Arvo (1992), and $\theta_{\rho},\phi_{\rho}$ which describe the relative distribution of the secular energy between the three axes. Of these, the $\phi_{0,j}$ are uniformly distributed on $[0,2\pi)$ for a homogenous buffer gas, and the three rotation matrix variables $x_{1},x_{2},x_{3}$ are uniformly distributed on $[0,1)$. The angles $\theta_{\rho}$ and $\phi_{\rho}$ are given by,

[TABLE]

and,

[TABLE]

In the ideal case, the temperature is a constant and is equal for each axis. Under these conditions, the probability distribution may be found analytically by starting from the thermal distribution for $A_{j}$ (see Eq. (63) below) and applying the standard methods for finding functions of random variables Riley et al. (2010). The final results are,

[TABLE]

and

[TABLE]

Random sampling of these distributions may be achieved through the inverse transform method by taking a random variable $u_{i}$ uniformly distributed in $[0,1)$ :

[TABLE]

and

[TABLE]

In practice, especially at higher mass ratios and values of the Mathieu $q$ parameter, this assumption breaks down, and the temperature for each axis is defined by a separate distribution. For linear RF traps with radial symmetry Major et al. (2005), the distributions for the radial $x$ and $y$ axes are identical and $f_{\phi_{\rho}}(\phi_{\rho})$ is approximately unchanged, but $f_{\theta_{\rho}}(\theta_{\rho})$ must be corrected to take into account the breakdown of equipartition between the $(x,y)$ and $z$ axes Chen et al. (2014). We will take each axis to have a different inverse temperature $\tilde{\beta_{j}}$, leading to the expected different mean energy for each axis. We approximate that, for a linear trap, $\tilde{\beta_{x}}=\tilde{\beta_{y}}$ and define $\xi=\tilde{\beta_{x}}/\tilde{\beta_{z}}=T_{z}/T_{x}$. Under these conditions, $f_{\theta_{\rho}}(\theta_{\rho})$ can be re-derived as ,

[TABLE]

For low mass ratios, the breakdown from equipartition is small and it suffices to set $\xi=\langle\xi\rangle\approx\tilde{E_{z}}/\tilde{E_{x}}$. At higher mass ratios, the increased correlation between the energy along each axis is such that $\langle\xi\rangle$ decreases less rapidly than predicted, and the higher-order moments of $\xi$ must be taken into account.

We may make use of these distributions and the expression for $\eta$ to calculate $\langle\eta\rangle$ and $\langle\eta^{2}\rangle$ by averaging it in turn over each of these distributions. The integrations over $\phi_{x,y,z},x_{1},x_{2},x_{3},\phi_{\rho},\theta_{\rho}$ may be performed analytically, leaving only the integration over $\tau$ to be performed numerically. We find that the remaining function of $\tau$ is periodic and so integrating over a single period is sufficient to calculate $\langle\eta\rangle$ and $\langle\eta^{2}\rangle$ in terms of the Mathieu parameters and the mass ratio. As shown in Fig. 5, the mean value calculated using this procedure and the stated distributions for the random variables involved is in excellent agreement with the values found from simulations.

B.1 Derivation of the heating coefficient $\kappa$

This approach may also be used to calculate the rate of transfer between the thermal energy of the atoms and the ion, Eq. (20) in the main text. By assuming an ion initially at the centre of the trap at rest, we may calculate the energy gained in a collision and apply equipartition to translate this to an increase in the temperature of the ion. Briefly, we take the Eq. (53) and set $E^{\prime}=3k_{b}T^{\prime},E=0,\epsilon=\frac{3}{2}k_{b}T_{a}$ to obtain $T^{\prime}=\kappa T_{a}$, where $\kappa$ is a function of the random variables $x_{1},x_{2},x_{3},\phi_{n},\theta_{n},\tau$ as defined earlier, the mass ratio $\tilde{m}$ and the Mathieu parameters $a_{j},q_{j}$. As outlined in the main text, we assume that $\kappa$ can be set to its mean value and so we average over the three uniform variables $x_{j}$ and the angular distributions of $\phi_{n},\theta_{n}$. The result is,

[TABLE]

where $w_{j}^{2}$ is the Wronskian defined by $(\dot{\mathrm{ce}}(a_{j},q_{j},\tau)\mathrm{se}(a_{j},q_{j},\tau)-\mathrm{ce}(a_{j},q_{j},\tau)\dot{\mathrm{se}}(a_{j},q_{j},\tau))^{2}$, which is time-independent Chen et al. (2014). Furthermore, when averaged over a complete period, $\langle\mathrm{ce}(a_{j},q_{j},\tau)^{2}+\mathrm{se}(a_{j},q_{j},\tau)^{2}\rangle=1$, and so we obtain,

[TABLE]

where the approximation applies in the limit $q\rightarrow 0$, i.e. a time-independent trap. In Fig. 6, we plot the calculated value of $\kappa$ for a time-dependent trap and the value extracted from numerical simulations for $q=0.5$ over a range of mass ratios, finding an excellent agreement between the two (black solid line). The approximate function for a time-independent trap (red dashed line) does not adequately describe $\kappa$ at this high value of $q$, but captures the dependence on the mass ratio.

B.2 Numerical methods

In order to generate numerical energy distributions for comparison with our statistical-mechanical models, we employ the Monte-Carlo simulation approach described by DeVoe DeVoe (2009). The ion’s initial state is calculated by first generating a value for the secular amplitude of motion, $\tilde{A}_{j}$, for each axis from the thermal distribution,

[TABLE]

where $\omega_{j}=\beta_{j}\Omega/2$, and $\beta_{0}$ is the initial inverse temperature. The inverse transform method may be applied to generate random values from this distribution by taking a uniform random number $u\in[0,1)$ and solving Abramowitz and Stegun (1970),

[TABLE]

The secular phase $\phi_{j}$ is taken from a uniform distribution $[0,2\pi)$ and used with the Mathieu amplitude $A_{j}=\tilde{A}_{j}/c_{0,j}$ to calculate the initial position and velocity of the ion from Eqs. (29) and (30). In this way, we ensure that the initial states are generated according to a fixed initial temperature and that the initial secular velocity is correctly correlated to the initial position. Since collisions occur at an energy-independent rate, they are described by Poisson statistics, and so the time between collisions may be sampled from an exponential distribution. The simulation is advanced directly from one collision to the next through use of the matrix propagator method to avoid errors introduced by numerical integration DeVoe (2009). Collisions are simulated by updating the ion’s velocity components according to Eq. (35). Velocities of atoms are drawn from a normal distribution at a fixed temperature, or set to zero to investigate only the multiplicative effects. Three uniformly distributed random numbers $\in[0,1)$ are required for the random rotation matrix, which is calculated as described in Ref. Arvo (1992). The simulation is then advanced to the next collision and the process repeats until a certain number of collisions have been simulated at which point the final secular energy is calculated.

The first two moments of the unknown distribution $f_{\ln\eta}(\ln\eta)$, i.e. $\mu$ and $\sigma^{2}$, are required as parameters in Eqs. (16) and (18) of the main text. We may estimate these parameters from the empirical distribution for $\eta$, Eq. (5) in the main text, or we may extract them from numerical simulations. The estimated value of $\sigma^{2}$ is typically only accurate to within $\pm 20\%$ and so these are instead extracted from simulations of single collisions with buffer-gas atoms at 0 K, for which $\eta=E^{\prime}/E$. Performing a large number of such simulations allows the estimation of $\mu$ and $\sigma^{2}$ from the resulting distribution of $\eta$. These estimates do not take into account the deviation from equipartition, since they are calculated from a thermal state, but are sufficiently accurate to allow comparison of the derived distribution to the numerical results.

For comparisons to values estimated with Eq. (24) and (27) in the main text, $\langle\beta\rangle$ and $n_{T}$ may also be found by maximum likelihood estimation (MLE) of numerical data. This has been shown to be more accurate than curve-fitting to binned data for Tsallis functions, avoids the issues of choosing an appropriate bin size, and requires only numerical root finding for one variable rather than a 2D optimization process Shalizi (2007). The procedure for standard Tsallis functions with $k=0$ is described in Ref. Shalizi (2007) and may be applied straightforwardly for other values of $k$.

Unless stated otherwise, numerical simulations are performed assuming an initial ion temperature $T_{0}$ of 1 mK, an arbitrary collision rate of 1000 s*-1*, an RF frequency of $8\times 2\pi$ MHz in a linear Paul trap with $q_{x}=0.2$, $q_{y}=-0.2$, $q_{z}=0.0$. $a_{z}$ is set such that the axial frequency is 100 kHz, and $a_{x,y}$ are equal to $-\frac{1}{2}a_{z}$. All computations were performed in Mathematica 10.2 using the internal implementation of the Mathieu functions and calculation of $\beta_{j}$. When necessary, coefficients for the Fourier series representation of the Mathieu functions were calculated using Miller’s algorithm Früchting (1969). Figures presented are log-binned and represent the results of at least 10000 trials.

B.3 Approximate calculation of the ion energy distribution for a zero-temperature buffer gas

To evaluate the ion energy distribution for a zero-temperature buffer gas (Eq. (17) in the main text),

[TABLE]

we use the Laplace method to find an approximate solution for $k=2$ following Ref. Asmussen et al. (2016). Briefly, the integrand of Eq. (B.3) has a maximum at the point $\beta_{n}=\hat{\beta}$,

[TABLE]

where $\mathcal{W}$ is the Lambert-$\mathcal{W}$ function Asmussen et al. (2016). We define $g(E_{n},\beta_{n})$ to be the logarithm of the integrand of Eq. (B.3) such that,

[TABLE]

and then replace $g(E_{n},\beta_{n})$ with its Taylor series to second order around the point $\beta=\hat{\beta}$. This leads to a Gaussian integral which can be analytically evaluated,

[TABLE]

which is asymptotically correct for $E_{n}\rightarrow\infty$, since as $E_{n}$ increases, the integral becomes more sharply peaked around $\hat{\beta}$ and the approximation becomes more precise Asmussen et al. (2016). The same method can be applied for an arbitrary value of $k$.

B.4 Derivation of the $\beta$ distribution for an ion in a buffer gas at finite temperature

Briefly, following Ref. Sornette and Cont (1997), we take Eq. (20),

[TABLE]

and rewrite this as,

[TABLE]

Using $(T_{i}-T_{i-1})/T_{i-1}\approx d\ln T/dt$, this expression can be converted into an overdamped Langevin equation for $x=\ln T$,

[TABLE]

where $\mu=\langle\eta\rangle-1$ and $\eta$ has been decomposed into its mean $\langle\eta\rangle$ and a fluctuating part $\hat{\eta}(t)$. This Langevin equation can be approximated by a Fokker-Planck equation which in steady state is given by,

[TABLE]

The boundary conditions are fixed by $f_{T}(0)\rightarrow 0$ and $f_{T}(\infty)\rightarrow 0$, corresponding to $f_{x}(x)\rightarrow 0$ for $x\rightarrow\pm\infty$. The corresponding solution of Eq. (72) is then given by,

[TABLE]

where $A$ is a normalization constant. Proceeding directly to $\beta=e^{-x}$ and normalizing for $\mu<0$ we find,

[TABLE]

where $\nu=-\frac{2\mu}{\sigma^{2}}$ and $b=\frac{-\mu}{k_{B}T_{a}\kappa}$. This is a gamma distribution, in agreement with the result obtained in Ref. Biró and Jakovác (2005) for multiplicative fluctuations with an additive constant. Due to the approximations used in the derivation, the forms given for the two parameters $\nu$ and $b$ correspond to those expected for a log-normal distribution of $\eta$, which is not the case for the present system. We therefore replace $\nu$ and $b$ with $n_{T}$ and $\langle\beta\rangle$ and calculate values for these appropriate for the observed form of $f_{\eta}(\eta)$ as described in the main text.