\section*{Methods}

\subsection*{Model of rovibrational FFR in rare gas diatomic ions}

Rovibrational FFR in rare gas diatomic ions arise due to the spin-orbit interaction coupling different non-relativistic electronic states. 
The Hamiltonian describing the rovibrational motion of a diatomic molecule including relativistic and angular couplings\cite{BusseryJCP06} can be restricted here to the lowest $^2\Sigma^+$ and $^2\Pi$ states and reads
\begin{subequations}\label{eq:model} 
\begin{eqnarray}
\label{eq:H_Hund_a}
H = T_{rad} + V_{ad} + H_{SO} + H_{rot}
\end{eqnarray}
with $T_{rad}=-\frac{\hbar^2}{2\mu}\frac{d^2}{dR^2}\mathrm{diag}(1,1,1)$ the radial kinetic energy, 
\begin{eqnarray}
\label{eq:V_Hund_a}
V_{ad} + H_{SO} = 
\begin{pmatrix}
V^{^2\Sigma^+}(R) & -\frac{\sqrt{2}}{3}\Delta & 0 \\
-\frac{\sqrt{2}}{3}\Delta & V^{^2\Pi}(R) + \frac{1}{3}\Delta & 0 \\
0 & 0 & V^{^2\Pi}(R) - \frac{1}{3}\Delta 
\end{pmatrix}
\end{eqnarray}
with $V^\alpha(R)$ the adiabatic non-relativistic interaction potentials  and $\Delta$ the spin-orbit interaction, and 
\begin{eqnarray}
\label{eq:Hrot}
H_{rot} = \frac{1}{2\mu R^2}
\begin{pmatrix}
J(J+1)+\frac{9}{4}+\epsilon(J+\frac{1}{2})&
\sqrt{2}-\epsilon\sqrt{2}(J+\frac{1}{2}) & -\sqrt{2(J(J+1)-\frac{3}{4})}\\
\sqrt{2}-\epsilon\sqrt{2}(J+\frac{1}{2}) & J(J+1)+\frac{5}{4}&-\sqrt{J(J+1)-\frac{3}{4}} \\
-\sqrt{2(J(J+1)-\frac{3}{4})} & -\sqrt{J(J+1)-\frac{3}{4}}& J(J+1)-\frac{3}{4}
\end{pmatrix}
\end{eqnarray}
\end{subequations}
with $J$ denoting the total molecular angular momentum quantum number and 
$\epsilon$ the spectroscopic parity, $\epsilon=\pm 1$. In Eqs.~\eqref{eq:model}, we have neglected
the $R$-dependence of both the spin-orbit interaction and the angular couplings
and taken $\Delta$ to be equal to the value of
Ar$^+(^2P)$. This approximation is well justified since for HeAr$^+$ these couplings are estimated to differ from their asymptotic values by less than 1\%,\cite{GemeinCPL90} and a similar behavior is expected for NeAr$^+$.
We have also neglected the relativistic Cowan-Griffin term and the diagonal adiabatic correction, both of which are tiny. 

While all numerical results are based on Hamiltonian~\eqref{eq:model} without further approximations, we have uses the Hund's case (c) representation in our discussion which is appropriate to illuminate the rovibational structure of HeAr$^+$ and NeAr$^+$.
This implies that $J$ and $\Omega$, the body-fixed $z$-projection
of total electronic angular momentum, are good quantum numbers, with $\Omega=1/2$ for states $X$ and $A_2$ (corresponding to $^2\Sigma^+$ and $^2\Pi_{1/2}$), and $\Omega=3/2$ for state $A_1$ (corresponding to $^2\Pi_{3/2}$). 

\subsection*{Potential energy curves}

Empirical adiabatic potential curves obtained by fitting high-resolution spectroscopic data are available\cite{CarringtonJCP95} and have been used for 
the HeAr$^+$ ion, whereas
\textit{ab initio} calculations have been carried out for NeAr$^+$.
The lowest adiabatic $^2\Sigma^+$ and $^2\Pi$ electronic states of NeAr$^+$ have been obtained from 
supermolecular calculations employing a spin-unrestricted coupled cluster method
with single, double and noniterative triple excitations [UCCSD(T)], as implemented
in the Molpro package.~\cite{molpro} We have used the
correlation-consistent 
augmented aug-cc-pV6Z basis set for both atoms, and an extra set of diffuse
even-tempered functions has been added to describe the Ar atom. 

\subsection*{Calculation of lifetimes}

Resonance positions and widths have been obtained by adding a complex absorbing potential (CAP) to all three potential energy curves and subsequent diagonalisation of the resulting non-Hermitian Hamiltonian.~\cite{Riss1993}
The calculations were performed on a spatial grid extending from $3\,a_0$ to $200\,a_0$, using the Fourier method~\cite{RonnieGrid96} with $4095$ basis functions. Such a large number of grid points was necessary to properly converge the lifetimes of the narrow NeAr$^+$ resonances. 
The transmission-free CAP of Ref.~\onlinecite{Manolopoulos2002} 
with strength $E_{\mathrm{min}}=1295$~cm$^{-1}$ and absorption range $D=20\,a_0$ was used.
With this choice of parameters, the lifetimes are converged with a relative error of about $0.001$~\% for HeAr$^+$ and about $0.1$~\% for NeAr$^+$. 
$\Delta\phi$ used in Fig.~\ref{fig:lifetimes} refers to the short-range equivalent of the scattering phase shift that is defined asymptotically.

Additionally, the liftimes were calculated perturbatively using Fermi's golden rule and Hund's case (c), where 
the perturbation is given by the radial and rotational couplings. 
The Hund's case (c) representation was obtained by a unitary transformation of the full Hund's case (a) Hamiltonian. 
The FFR widths are determined as the transition rates of the $A_2$ rovibrational bound levels
to the continua of the $X$ and $A_1$ states. Since there are no radial couplings between the $A_1$ and $A_2$ states, the $A_1$ contribution to the FFR widths, respectively lifetimes, is only minor.
Bound and energy normalized continuum states were obtained by diagonalizing the unperturbed Hamiltonian using the same grid as described above. 
The transition rates were calculated for the continuum state with energy closest to the energy of the bound level. 
For NeAr$^+$, the perturbative calculations almost exactly reproduce the lifetimes obtained non-perturbatively using the CAP and full couplings. The largest deviations occur at the maximum lifetimes, where the resonance widths become very small and overlap integrals become quite sensitive to the energy and phase of the continuum state. Nevertheless, the positions of the peaks for phase protected states are in full agreement for all rovibrational levels.
For HeAr$^+$, perturbative and non-perturbative results also agree very well except for $v=0$ and low $J$ where quantitative deviations up to a factor of four occur. Except for $v=0$ and odd spectroscopic parity, the lifetimes decrease with $J$. This monotonous dependence is in disagreement with earlier calculations that predicted an oscillatory behavior with rotational excitation.\cite{GemeinCPL90} We attribute the disagreement mainly to a different treatment of the spin-orbit interaction, with Ref.~\onlinecite{GemeinCPL90} accounting for its effect only on the potential energy curves but not on the kinetic energy.

\subsection*{Estimating the lifetime scaling with $\mu$ from the Morse oscillator approximation}

For both HeAr$^+$ and  NeAr$^+$, the $A_2$ state interaction potential for $J=0$ is well approximated around the potential minimum by the Morse oscillator,
$V(R) = D \left( e^{-2 a (R-R_0)} - 2 e^{-a (R-R_0)} \right)$
with $D$ the well depth, $a$ a constant determining the width (with smaller $a$ resulting in a wider potential well) and $R_0$ the equilibrium distance. The decay rates are given by
\begin{align*}
  \gamma_v \propto \sin\left(\delta + \phi_v(k)\right)^2 \,
  \left| \tilde\psi_v(k) \right|^2
	\;,
\end{align*}
where  $\phi_v(k) = \mathrm{arg}[\tilde\psi_v(k)]$ and $\tilde\psi_v(k)$  the Fourier transform of the Morse oscillator eigenstate. The absolute value squared $\left| \tilde\psi_v(k) \right|^2$ thus determines the minimum lifetime that can be expected. On top of a smooth variation  with $k$ (or energy) due to
$\left|\tilde\psi_v(k)\right|^2$, the decay rates $\gamma_v$ oscillate, with the roots corresponding to phase protection, cf. Eq.~\eqref{eq:condition}.
Exploiting the fact that $\tilde\psi_v(k)$ is known analytically, in terms of complex   $\Gamma$-functions,~\cite{BancewiczJPA98} we obtain  the scaling of the lifetimes, resp. decay rates, with the parameters of the potential and the reduced mass. Specifically, we find that the number of roots in $\gamma_v$ and thus the chance for phase protection increase with $D$, $R_0$, $a$, and $\mu$.
This agrees well with our observation of phase protection for NeAr$^+$
and its absence for HeAr$^+$, since the reduced mass and the well depth are larger for NeAr$^+$ than for HeAr$^+$, while the position of the minimum and the potential width are very similar for the two molecules.
More generally, our estimate predicts an isotope effect on predissociation lifetimes which agrees well with experimental observations for N$_2^+\,\,$ \cite{GoversJPB73} and Ne$_2^+$.\cite{GluchEPJD07}


\subsection*{Velocity Map Imaging after Merged Beam Penning Ionization }
Cold Penning ionization reactions of He$^*$-Ar and Ne$^*$-Ar are realized by the merged beam technique described by Henson et. al.\cite{HensonSci12}. The pulsed supersonic  beams of He/Ne and Ar are created by two Even-Lavie valves\cite{Even2000} positioned at a 10 degree angle. The He/Ne beam is excited to a metastable paramagnetic state by dielectric barrier discharge (DBD)\cite{Luria2009} and are deflected by a curved magnetic guide such that it is co-propagating with the Ar beam. The short opening duration of the Even-Lavie valves allow collision energies as low as 10$\,$mK $\times k_B$ to be reached.\cite{Shagam2013}

The two beams meet in the reaction region inside the velocity map imaging (VMI) detector.
The beams are skimmed by two parallel razors before the VMI detector to limit their distribution to a 1$\,$mm plane. The velocity map imaging setup consists of 8 separate plates electrodes positioned perpendicular to the beam path, followed by a 9$^{th}$ grounded electrode. These are followed by a flight tube approximately 1$\,$m long. The electrode immediately above the beam (second plate) has an aperture of 1$\,$cm in diameter for the ions to pass through, limiting the overall region from which the ions can be imaged to a 1$\,$mm by 10$\,$mm diameter cylinder. The remaining 6 electrode rings have an aperture of 40$\,$mm diameter.

Argon is ionized upon collision with He($^3S_1$)/Ne($^3P_2$), leading to Ar$^+$ as well as HeAr$^+$/NeAr$^+$ for associative ionization. The ions are then accelerated towards the micro-channel plate (MCP) detector. The MCP is followed by a phosphor screen that facilitates the imaging of the product ions by a charge couple device (CCD) camera. The center of the individual product ions is found as these are accumulated.

Since the time of ionization is not precisely known for our reaction, we must adjust the VMI method to allow for mass selectivity. To this end, we have used two  VMI operation modes: (i)  constant voltages mode  where particles with different mass-to-charge ratio  are spatially separated by an external magnetic field and (ii) a pulsed mode \cite{Mikosch2007,Mikosch2008}  that  selects particles with a certain mass-to-charge ratio by time of flight (TOF). In our experiments, the products from Penning (Ar$^{+}$) and associative ionization (HeAr$^{+}$ and NeAr$^{+}$) for the two reactions were separated using these methods.

(i) Constant voltages with bias magnetic field: In this mode, all 8 plates of VMI have constant voltages throughout the measurement. By applying a uniform magnetic field parallel to the neutral beam collision axis, the product ions get shifted on the detector according to their respective mass-to-charge ratios. The VMI resolution is not affected since the magnetic field is uniform. Figure~\ref{fig:vmi}(c) is acquired with a magnetic field of approx. 10$\,$Gauss along the entire VMI.

 (ii) Pulsed mode: This approach combines the idea of velocity mapping with TOF mass spectrometry allowing the imaging of a single mass/charge on the detector. In this method, all the plates are simultaneously pulsed to their respective imaging voltages during interaction time. The imaging pulse to the VMI electrodes follows a cleanup pulse with a period defining the `product accumulation time'. To keep the VMI resolution optimal, this timing is limited by the expansion of the products outside the imaging volume. When the imaging pulse is applied, all the ions formed during `accumulation time' fly up to the detector and get separated during their flight times, allowing TOF mass spectrometry to be performed. The MCP at the end of the flight tube is pulsed on for 50 ns selecting a single mass product for imaging. This method is used for collecting the data for Fig.~\ref{fig:vmi}(b), which shows Ar$^+$ obtained from He*-Ar collisions. Note that the actual size of the image formed on the detector is 5$\,$mm in diameter.

The velocity detected by the VMI is calibrated by creating a cold beam of Rydberg argon atoms at a range of velocities by DBD on a beam of argon. The beam is then field ionized when the VMI is pulsed and its position is detected for every velocity. The velocity of the argon beam is measured by another MCP positioned in front of the argon valve. Both the UV light from the discharge and the metastable argon atoms are detected by the MCP giving the time of flight.

\subsection*{Measurement of decay rates}
NeAr$^{+}$ ions were formed in crossed molecular beams of argon and neon, produced by Even-Lavie valves. The neon beam was excited to the metastable $^{3}P_2$ level with a dielectric barrier discharge (DBD),  with all ions generated during the excitation deflected away. The beams crossed at 12 degrees, their angles respective to the axis of the skimmer. Assuming translational beam temperatures of 1$\,$K, we obtain a collision energy of around 16$\,$meV, which corresponds to a reaction temperature of 190$\,$K. The collision products were skimmed before being accelerated to 4.2 keV, directed through several beam manipulation elements and injected into the EIBT trap. The setup is depicted in Fig.~\ref{fig:lifet_exp}(a), for further details we refer to Ref. \onlinecite{rahinov2012}. Mass-sensitive detection of the trapped particles by a pick-up electrode was used to optimize the setup parameters such as manipulation voltages and trap timing. Neutral products stemming from dissociation of the ions or neutralizing collisions with the background gas are lost from the trap and were detected on a MCP located on-axis behind the trap. NeAr$^{+}$ could be trapped with minimal contaminating ions, as confirmed by FFT analysis of the pick-up signal and the distinct oscillation observed in the neutral count, as depicted in Fig.~\ref{fig:lifet_exp}(b). The experimental data is binned in 200$\,$ns-intervals to bring out the oscillations on the scale of the NeAr${^+}$ frequency in the trap (Fig.~\ref{fig:lifet_exp}(b)), and in 1$\,$ms-intervals for showing the overall decay in Fig.~\ref{fig:lifet_exp}(c).


The measured data is fitted to a multi-exponential decay,
$f(t)=\sum_i a_i\exp(-t/\tau_i)$. The best agreement,
corresponding to a correlation coefficient of $R^2=0.995$, 
is found when including six terms in the fit, cf. Table~\ref{tab:fit},
where the longest decay time corresponds to collisions with the background gas.

\begin{table}[tbp]
\centering
\begin{tabular}[c]{|c|c|c|c|c|c|c|}
\hline
$\tau_i$ & 50.0$\,\mu$s & 1.50$\,$ms & 2.03$\,$ms & 3.44$\,$ms & 15.8$\,$ms & 1.72$\,$s \\ \hline
$\quad a_i\quad$ & $1.01\times 10^{4}$ & $7.48\times 10^{3}$ & $6.06\times 10^{3}$ & $8.31\times 10^{3}$ & 406 & 40.9 
\\\hline
\end{tabular}
\caption{Decay times $\tau_i$ with weights $a_i$ obtained when fitting the experimental data to $\sum_i a_i \exp(-t/\tau_i)$.}
\label{tab:fit}
\end{table}


\begin{acknowledgments}
CPK is grateful for a Rosi and Max Varon Visiting Professorship. 
Financial support from the German-Israeli Foundation, grant
no. 1254, is gratefully acknowledged. 
\end{acknowledgments}

\subsection*{Author contributions}

AB, WS and DMR carried out the theoretical calculations. JWR designed and YS built the VMI apparatus, and YS, PP and NB performed the VMI measurements. NJH carried out the lifetime measurements with the help of OH. EN and CPK planned and supervised the project. All authors contributed to the discussion of the results and the writing of the manuscript.


