\section{Supplemental material}
\subsection{Interaction time calculation}
We measured the atomic profile using time-of-flight absorption imaging. The atomic cloud is cigar-shaped with approximately Gaussian density, $n_{at}=n_0e^{-x^2/2\sigma_y^2-y^2/2\sigma_z^2-z^2/2\sigma_x^2}$ distribution in both radial, $\sigma_x=\sigma_y$=5 $\mu$m, and axial $\sigma_{z}$=70 $\mu$m axes. We measured the peak atomic density to be, $n_0=4\cdot10^{11}$ cm$^{-3}$ and the number of atoms in our trap is approximately 7,500.

The Langevin rate of collisions, $\Gamma=2 \pi n_{at} \sqrt{C_4/\mu}$, is independent of the collision energy. For our atomic experimental parameters of, $C_4=1.1\cdot10^{-56}$ J$\cdot$m$^4$ and $\mu\approx43.7$ amu, it yields a collision rate, $\Gamma$=0.9 kHz, for the atomic peak density.

In the first axial experiment, we prepare the atoms 50 $\mu$m vertically above the center of the ion trap. We then transport the atoms to overlap with the center of the trap using a piezo driven mirror which controls the position of the horizontal atoms' crossed dipole trap beam. The atoms' velocity during the transport is approximately constant, $v_{t}$=0.5 $\mu$m/ms. It takes the atoms 100 ms to reach the center of the trap. 

Due to the transport, the ion interacts first with the outskirts of the atomic profile which leads to much reduced and time dependent collision rate. To take into account the effect of the transport, we calculate the effective interaction time, $t_\textrm{eff}$, as if the atoms overlap the center of the trap from the beginning of the experiment,
\begin{equation}
	t_\textrm{eff}=\int\limits_{-\infty}^{t_t-x_t/v_t}e^{-\frac{(v_{t}t)^2}{2\sigma_y^2}}dt=t_0\left( \textrm{erf}\left( \frac{v_t t_t - x_t}{\sqrt{2}\sigma_y}\right) + 1 \right).
\end{equation}
Here, $0 \le t_t \le 100$ ms is the transport time, $x_t$=50 $\mu$m is the transport distance and $t_0=\sqrt{\frac{\pi}{2}}\frac{\sigma_y}{v_t}=9.2$ ms is the effective time for full transport, $t_t$=100 ms. The experimental interaction time, given in Fig. 2a-f, is then calculated by,
\begin{equation}
	t_\textrm{int.}=
	    \begin{cases}
			\begin{aligned}
				&t_\textrm{eff} & ,0 \le t_t \le x_t/v_t\\
				&t_0+(t_t-x_t/v_t) & ,t_t \ge x_t/v_t.
    		\end{aligned}
		\end{cases}
\end{equation}
For experiments with only single to few collisions (Fig. 2b-e in the main text) we stop the transport and release the atoms before they reach the center of the trap.

Our transport calculation takes into account the motion of the atoms with respect to the center of the ion trap. However, the ion is prepared in a classical coherent state and is preforming an oscillatory motion with an amplitude of $z_\textrm{ion}=71$ $\mu$m in the axial direction relative to the trap center. This motion amplitude is comparable to the atomic cloud size. We numerically calculate the effective reduced density due to this motion, $n_\textrm{eff}=0.79\cdot n_0$, which yields a slower Langevin collision rate, $\Gamma_\textrm{eff}$=0.71 kHz.

In the second, radial, experiment, we compare the experimental time to simulation time directly by implementing the motion of the atomic cloud into the simulation. 

\subsection{Calculation of the energy in the axial and radial modes of the trap}
The results of a numerical calculation for initial energy of 225K with different initial axial to radial energy distribution is given in Fig. \ref{fig:energy3d}. This is the same numerical calculation which was used throughout the main text. Here, we present the energy dynamics separated to the axial (z) and radial (x+y) modes of the trap. 

\begin{figure}
	\centering
	\includegraphics[width=\linewidth,trim={0cm 0cm 1cm 0cm},clip]{Energy_AllModes_3D.eps}\\
	\caption{\textbf{Numerical calculation of the ion's energy dynamics during cooling}. The energy of the ion in the axial (z-axis) and radial (y-axis) is plotted as function of the time for reaching the Doppler cooling temperature (x-axis). Color represent different axial to radial energy ratios (legend).}
	\label{fig:energy3d}
\end{figure}

\subsection{Likelihood ratio test}
We calculate the likelihood, $\mathcal{L}_i^m$, of a single data instance, $x^m(t)$, to be described by each of the numerical trajectories, $\theta_i$:
\begin{equation}
    \mathcal{L}_i^m\left(x^m(t)|\theta_i\right)=\prod_t \frac{\theta_i(t)^{x^m(t)}e^{-\theta_i(t)}}{x^m(t)!}.
\end{equation}
Here, the trajectory, $\theta(t)$, is sampled at the experimental sampling times due to the 1 ms binning. The subscript, $i$, indicates which numerical calculation is used to calculate the likelihood. The superscript $m$ indicates what instance of data is evaluated. In the LRT we compare the likelihoods values from all the numerical calculations, $i$, for each experimental data, $m$, and pick the one with the maximal value: 
\begin{equation}
    \mathcal{L}_i^m\left(x^m|\theta_i\right)>\mathcal{L}_j^m\left(x^m|\theta_j\right) \forall j\neq i.
\end{equation}
We performed a stochastic numerical simulation to analyze the LRT method (see details in the following paragraph). The type I errors for wrongly detecting mode $i$ as mode $j$ are less than 1\% when the cooling trajectory is longer than 100 ms. We use the numerical calculation, $i$, chosen by the LRT to extract the energy of the ion for the experimental instance, $m$.

We used the results of the numerical calculation for different initial energy distributions between the axial and radial modes of the trap (Fig. 1b in the main text) to test our LRT analysis. We sampled the numerically calculated fluorescence trajectories according to the experimental photon-counting binning of 1ms and created a simulated stochastic data using a Poisson random distribution. We then performed the LRT analysis on the simulated data (100,000 trajectories for each numerical calculation) and determined the type I error due to Poisson noise.

When performing likelihood analysis to single-shot data we have the freedom to choose the time span on which we perform the likelihood. We usually choose a span from -100 ms to 10 ms which captures all the discrepancies of the different numerical calculations (see Fig. 1b of the main text). However, 100 ms cooling time corresponds to initial ion's energy of $\sim$100K. In order to perform LRT on colder events we need to reduce the likelihood span. In Fig. \ref{fig:likelihoodspan} we show the results of the stochastic simulation for different likelihood spans. We see that as the likelihood span reduces, the percentage of type I errors increases. For that, we choose a minimal span of 50 ms cooling time which limits our LRT analysis to $\sim$50 K. Below this energy we estimate the energy using the the numerical simulation with the initial energy ratio of Ez/Er=0.15.

\begin{figure}
	\centering
	\includegraphics[width=\linewidth,trim={0cm 0cm 1cm 0cm},clip]{TypeIErrors_span.eps}\\
	\caption{\textbf{Stochastic simulation of the LRT.} Type I error as function of the data span used for the LRT. Different colors indicates from which numerical simulation of different axial to radial energy ratio (legend) the stochastic data was created.}
	\label{fig:likelihoodspan}
\end{figure}

\subsection{Calculation of the scattering angle distribution}
We consider an ion with mass, m, and a velocity, v, which we choose to coincide with the z-axis of the trap, and an atom with mass, M, at rest in the lab-frame. For simplicity, we first consider the case in which m=M. In the center-of-mass frame, the ion's velocity is v/2 and the atom's velocity is -v/2 along the z-axis. After an elastic collision which preserves both momentum and energy, the ion's and atom's velocities remain v/2 with opposite directions, however, at angle $\theta$ from the z-axis. Since in head-on collision there is an equal probability to scatter into 4$\pi$, $\cos(\theta)$ is uniformly distributed between -1 to 1 where we implicitly choose the scattering angle $\phi$ to be zero.

In the lab frame, the velocity of the ion after a collision is given by:
\begin{equation}
    v/2(1+\cos(\theta))\hat{z}+v/2\sin(\theta)\hat{x}.
\end{equation}
Here, we rotate the radial axes such that the scattering products lie along the x-axis. The total energy of the ion before a collision, $E_b$, is:
\begin{equation}
    E_b=mv^2/2+mz^2\omega_z^2/2,
\end{equation}
where $z$ is the (unknown) position of the ion in the harmonic z-axis of the Paul trap and $\omega_z$ is the trap's harmonic frequency. The ion's energy after a collision, $E_a$, is given by:
\begin{equation}
    E_a=m(v/2)^2(1+\cos(\theta))^2/2+m(v/2)^2\sin^2(\theta)/2+mz^2\omega_z^2.
\end{equation}
Here, we assume that the collision is instantaneous such that the position of the collision, (0,0,z), is not changed during the collision. For that, there is no contribution from the radial trapping potential. To get rid of the dependence of the energy on the collision position, z, we look on the energy difference during a collision:
\begin{equation}\label{eq:Ediff}
    \Delta E=E_b-E_a=m(v/2)^2(1-\cos(\theta)).
\end{equation}
This expression depends on the ion's initial velocity, v, which is also unknown since the ion is in an harmonic potential (only the initial energy is known). To get rid of this dependence we divide the energy difference with the radial energy:
\begin{equation}
    \frac{\Delta E}{E_r}=2\frac{1-\cos(\theta)}{1-\cos^2(\theta)}.
\end{equation}
For the general case of an ion with mass m and an atom with mass M, this expression only slightly changes:
\begin{equation}\label{eq:1}
    \frac{\Delta E}{E_r}=2\frac{m}{M}\frac{1-\cos(\theta)}{1-\cos^2(\theta)}.
\end{equation}
We can solve this second order polynomial equation to get an expression for $\cos(\theta)$ which is given in Eq. \ref{eq:costheta} of the main text.

The red line in Fig. 3 of the main text is the Langevin scattering angle distribution given in terms of Eq. \ref{eq:1} (Fig. 3a) and $\cos(\theta)$ (Fig. 3b). We can reproduce this theoretical prediction exactly using our MD simulation as can be seen in Fig. \ref{fig:scattering_angle_sm} by the magenta circles. However, since each experimental measurement of the energy is destructive, we cannot obtain the exact value of the energy difference in a collision given in Eq. \ref{eq:Ediff}. Instead, we use the mean value of the initial energy, $\left<E_b\right>=183.7$ K, to get an approximated value of the energy difference,
\begin{equation}\label{eq:Ediffmean}
    \left<\Delta E\right>=\left<E_b\right>-E_a.
\end{equation}
The error for $\left<\Delta E\right>$ is given by the standard deviation of the initial energy, $\Delta \left<E_b\right>=2.3$ K. 

Another error source is due to the binning of the radial energy due to the LRT analysis. We bin the radial energy to the same binning used in the experimental data analysis. We calculate $E_z/E_a$ of each simulation run. For $E_z/E_a>0.9923$ we consider the event as no collision ($E_z/E_a=1)$. For $E_z/E_a$ values in the ranges given by the vector,  (0.9923, 0.9690, 0.9375, 0.7134, 0.3284, 0), we bin $Ez/E_a$=(0.9846, 0.9534, 0.9216, 0.5053, 0.1515). The radial energy is given by, $E_r=E_a(1-E_z/E_a)$. The results of this binning together with the mean energy effect is given by the blue line in Fig. \ref{fig:scattering_angle_sm} and Fig. 3 of the main text. The gray area indicates the values the simulation can vary when taking the edges instead of the mean value of each bin.

\begin{figure}
	\centering
	\includegraphics[width=\linewidth,trim={0cm 0cm 1cm 0cm},clip]{Figure4_sm.eps}\\
	\caption{\textbf{Scattering angle simulation.} Histograms of the energy difference divided by the radial energy after a collision (a) and the cosine of the scattering angle (b). Red lines are the theoretical prediction of Langevin scattering. Magenta circles are MD simulation results of exactly single collision of ion initialized in the axial direction with zero temperature atoms and no excess micromotion. The blue line is the same simulation, however, in the analysis we take into account the experimental limitations in the data analysis described in the text. Gray areas gives the uncertainty region for the data analysis.}
	\label{fig:scattering_angle_sm}
\end{figure}

\end{document}
/