\section{Introduction}
The interaction between ions and neutral atoms or molecules is of
central importance in atmospheric and astronomical processes and
environments. Prominent species in the interstellar environment
include H$_3^+$, CH$_2^+$ , HCO$^+$ and N$_2$H$^+$, among
others.\cite{tielens:2013} Additionally, ions are also considered to
play an important role in the formation of atmospheric
aerosols.\cite{hao:2015}\\

\noindent
Very recently,\cite{stutzki:2019} the HeH$^+$ ion, which was the first
molecule of the primordial universe,\cite{zyg98:151} has been detected
in interstellar space and means for the direct detection of H$_2^+$
have been discussed.\cite{black:2012} However, although H$_2^+$ is
most likely formed and present in space, e.g. through the HeH$^+$ + H
$\rightarrow$ H$_2^+$ + He reaction,\cite{stutzki:2019} (which is
believed to be the first atom-diatom reaction in the
universe\cite{lep02:R57}) collisions with H and H$_2$ are also
important loss channels of the ion. Nevertheless, with H$_2^+$ present
in the interstellar medium, it is also likely that the H$_2^+$--He
complex is formed. Hence, H$_2^+$--He plays an important role already
in the early stages of the Molecular Universe.\\

\noindent
The interaction between He and H$_2^+$ is also important for the
rotational cooling of H$_2^+$ through collisions with Helium as the
buffer gas.\cite{schiller:2017} This is an attractive way to generate
translationally and internally cold H$_2^+$ ions suitable for
precision measurements.\cite{schiller:2007,schiller:2014} With even
further increased precision and quantum state control of the ions,
fundamental natural constants such as the ratio of the electron to the
proton mass, $m_{\rm e} / m_{\rm p}$, can be determined with
unprecedented accuracy.\\

\noindent
Another process of interest which has been recently investigated is
the Penning ionization of $^3$S excited He colliding with
H$_2$.\cite{klein:2017} This process produces H$_2^+$--He with
sufficient energy to dissociate into ground state and rovibrationally
excited He and H$_2^+$ fragments. Such rovibrationally inelastic
half-collisions are particularly sensitive to the long-range part of
the intermolecular potential, which is dominated by polarization
interactions induced by the charge and quadrupole of
H$_2^+$. Furthermore, several long-range states for H$_2^+$--He have
been characterized from microwave spectroscopy and by using electric
field extraction.\cite{car96:395,gam02:6072} However, the
interpretation of these spectra has remained elusive, in part due to
the limited accuracy of the available potential energy surfaces.\\

\noindent
In the past, several PESs have been constructed at different levels of
theory to investigate the spectroscopy and dynamics of the H$_2^+$--He
complex.\cite{jos87:704,fal99:117,meu99:3418,pal00:1835,ram09:26,faz12:244306}
To characterize spectral transitions in the microwave region, an
accurate long-range potential is required.\cite{fal99:117,meu99:3418}
However, the level of theory used for the electronic structure
calculations in these earlier efforts was rather modest by today's
standards. Full configuration interaction (FCI) with the cc-pVQZ basis
has been used more recently but no explicit analytical representation
was included.\cite{ram09:26} Later, using the {\it ab initio} data of
Ref. \citen{ram09:26} a new PES was constructed by including an
explicit analytical formula only for the diatomic
potentials.\cite{faz12:244306}\\

\noindent
In the present work high-level electronic structure methods combined
with advanced representation techniques for global potential energy
surfaces and accurate representation of the long-range potential are
used. With these PESs quantum calculations of all bound states of
H$_2^+$--He with H$_2^+$ in its ground electronic and vibrational
state are then carried out. First, the computational methods are
presented, followed by the discussion of the bound states computed and
their interpretation in view of the near-dissociative states.\\



\section{Computational Methods}

\subsection{The Potential Energy Surfaces}
Two different levels of theory - a) multi reference configuration
interaction level including the Davidson correction
(MRCI+Q)\cite{wer88:5803,kno88:514} with the augmented Dunning-type
correlation consistent polarize hexaple zeta
(aug-cc-pV6Z)\cite{wil96:339} basis set and b) full configuration
interaction \cite{kno84:315,kno89:75} with the augmented Dunning
type correlation consistent polarized quintuple zeta
(aug-cc-pV5Z)\cite{dun89:1007,woo94:2975} basis set - are used in the
present work to calculate the {\it ab initio} energies. Initial
orbitals for the MRCI calculations were obtained using the complete
active space self-consistent field
(CASSCF)\cite{wen85:5053,kno85:259,wer80:2342} method with three $1s$
orbitals of H and He in the active space. The Molpro\cite{molpro}
software was used to perform all electronic structure calculations.\\

\noindent
The grids for the {\it ab initio} energy calculations are set up in
Jacobi coordinates ($R, r, \theta$).  Here, $r$ is the H$_2^+$ bond
length, $R$ is the distance between He and the center of mass of the
H$_2^+$ ion, and $\theta$ is the angle between $\vec{r}$ and
$\vec{R}$. The angular grid is defined by Gauss-Legendre quadrature
points chosen in the range between $0 \leq \theta \leq 90^\circ$ given
the spatial symmetry of the system. Details of the angular and radial
grids for the MRCI+Q and FCI calculations are given in Tables
S1 and S2 in the supporting
information.\\

\noindent
The complete adiabatic surface for H$_2^+$--He can be expressed as a
many-body expansion\cite{var88:255}
\begin{equation}
\begin{split}
V_{\rm HeHH'}({r_{\rm HeH}},{r_{\rm HeH'}},{r_{\rm HH'}}) & =
V^{(1)}_{\rm He} + V^{(1)}_{\rm H} + V^{(1)}_{\rm H^+} \\ & +
V^{(2)}_{\rm HeH^+}(r_{\rm HeH}) + V^{(2)}_{\rm HeH'^+}(r_{\rm HeH'})
+ V^{(2)}_{\rm HH'^+}(r_{\rm HH'}) \\ & + V^{(3)}({r_{\rm
    HeH}},{r_{\rm HeH'}},{r_{\rm HH'}}),
\end{split}
\label{mb}
\end{equation}
where ${r_{\rm HeH}},{r_{\rm HeH'}}$ and ${r_{\rm HH'}}$ are the
distances between the respective atoms, and $V_{\rm HeHH'}({r_{\rm
    HeH}},{r_{\rm HeH'}},{r_{\rm HH'}})$ is the total energy of the
triatomic system at the corresponding geometry. The $V_i^{(1)}$ are
the atomic energies , whereas the $V_i^{(2)}(r_i)$ and
$V^{(3)}({r_{\rm HeH}},{r_{\rm HeH'}},{r_{\rm HH'}})$ are the two- and
three-body interaction energies, respectively, at corresponding
configurations.\\

\noindent
In general, two body interaction energies, i.e., the diatomic
potential, for a molecule AB can be expressed
as\cite{agu92:1265,faz12:244306}
\begin{equation}
 V^{(2)}_{\rm AB}(R_{{\rm AB}}) = \frac{c_0e^{-\alpha_{{\rm AB}}
     R_{{\rm AB}}}}{R_{{\rm AB}}} + \sum_{i=1}^{M}c_i\rho^i_{{\rm AB}}
 + V_{\rm long}(\tilde r),
 \label{ap2}
\end{equation}
with $c_0 > 0$ to ensure $V_{\rm AB}(R_{{\rm AB}}) \rightarrow \infty$
at $R_{{\rm AB}} \rightarrow 0$ and $\rho_{{\rm AB}} = R_{{\rm
    AB}}e^{-\beta_{{\rm AB}}^{(2)}R_{{\rm AB}}}$. The long range part,
$V_{\rm long}(\tilde r)$, can be written as\cite{fal99:117}
\begin{equation}
V_{\rm long}(\tilde r) = - \frac{\alpha_d q^2}{2\tilde r^4} -
\frac{\alpha_q q^2}{2\tilde r^6} - \frac{\alpha_o q^2}{2\tilde r^8} -
\frac{\beta_{ddq} q^3}{6\tilde r^7} - \frac{\gamma_d q^4}{24\tilde
  r^8},
 \label{dal}
\end{equation}
where $q$ is the charge, and $\alpha_d, \alpha_q$ and $\alpha_o$ are
the dipole, quadrupole and octopole polarizabilities for H and He,
respectively. $\beta_{ddq}$ and $\gamma_d$ are the first and second
hyperpolarizabilities, respectively. The values for the
polarizabilities of He and H are taken from
Refs. \citen{fal99:117,bis95:15} and $\tilde r$ is defined
as\cite{vel08:084307}
\begin{equation}
 \tilde r = r + r_l \exp^{(-(r-r_e))},
\end{equation}
to remove the divergence of the long range terms at short H-H and H-He
separations. Here, $r_l$ is a distance parameter and $r_e$ is the
equilibrium bond distance of the diatomic molecule. The parameters
used in this work to obtain the diatomic potentials are given in Table
\ref{tab:params}.\\





\begin{table}[ht]
\caption{Parameters used in the diatomic potentials. All values are in
  atomic units.}
\begin{tabular}{l|r|r}
    \hline
    \hline
         & H & He  \\
    \hline
    \hline
    Dipole polarizability $\alpha_d$&4.5&1.384\\
    Quadrupole polarizability $\alpha_q$ &15.0&2.275 \\
    Octopole polarizability $\alpha_o$&131.25&10.620 \\
    First hyperpolarizability $\beta_{ddq}$&159.75&20.41 \\
    Second hyperpolarizability $\gamma_d$&1333.125&37.56 \\
    \hline
    & H$_2^+$ & HeH$^+$  \\
    \hline
    \hline
    $r_l$ &10.0&8.0\\
    $r_{\rm eq}$&2.005815&1.4633\\
    \hline
    \hline
  \end{tabular}
  \label{tab:params}
\end{table}


\noindent
The linear parameters $c_i$ and the nonlinear parameters $\alpha_{{\rm
    AB}}$ and $\beta_{{\rm AB}}^{(2)}$ in Eq. \ref{dal} are determined
by fitting the expression with the {\it ab initio} energies using the
Levenberg-Marquardt nonlinear multidimensional fitting
method.\cite{pre92} The optimized linear and nonlinear parameters for
the diatomic potentials calculated via fitting are given in Tables
S3 and S4.\\

\noindent
The three-body interaction energies, $V^{(3)}({r_{\rm HeH}},{r_{\rm
    HeH'}},{r_{\rm HH'}}) \equiv V^{(3)}(r, R, \theta)$ are calculated
from Eq. \ref{mb}. For a particular configuration of He-H$_2^+$,
$V^{(3)}(r, R, \theta)$ can be calculated using the reproducing kernel
Hilbert space\cite{ho96:2584} (RKHS) approach.\\

\noindent
The procedure for computing the analytical energy of a given
configuration from a set of known {\it ab initio} energies is briefly
described here.  According to the RKHS theorem, the value of a
function $f({\bf x})$ can be evaluated from a set of known values
$f({\bf x}_i)$ at positions ${\bf x}_i$ as a linear combinations of
kernel products
\begin{equation}
\tilde{f}({\bf{x}}) = \sum_{i=1}^{N} c_i K ({\bf{x}}, {\bf{x}}_i),
\end{equation}
where $c_i$ are the coefficients and $K({\bf{x}}, {\bf{x}}_i)$ are the
reproducing kernels. The coefficients are calculated from the known
values by solving a set of linear equations
\begin{equation}
 f({\bf{x}}_j)= \sum_{i=1}^{N} c_i K ({\bf{x}}_i, {\bf{x}}_j).
\end{equation}
Here it is worth mentioning that the RKHS approach exactly reproduces
the input data at the reference points. The derivatives of
$\tilde{f}({\bf{x}})$ can be calculated analytically from the kernel
functions $K({\bf{x}}, {\bf{x}}')$. For a multidimensional function
the $D$-dimensional kernel can be constructed as the product of $D$
1-dimensional kernels $k(x, x')$
\begin{equation}
 K({\bf{x}}, {\bf{x}}') = \prod_{d=1}^{D}k^{(d)}(x^{(d)}, x'^{(d)}),
\end{equation}
where $k^{(d)}(x^{(d)}, x'^{(d)})$ are the 1-dimensional kernels for
$d$-th dimensions.\\

\noindent
For the radial dimensions ($r$ and $R$) a reciprocal power decay
kernel
\begin{equation}
\label{k24}
 k^{[2,4]}(x,x') = \frac{2}{15}\frac{1}{x^5_{>}} -
 \frac{2}{21}\frac{x_<}{x^6_>},
\end{equation}
is used in the present work where, $x_>$ and $x_<$ are the larger and
smaller values of $x$ and $x'$. The value of this kernel smoothly
decays to zero according to $x^{-4}$ as the leading term in the
asymptotic region, which gives the correct long-range behavior for
atom-diatom type interactions. For the angular dimension, a Taylor
spline kernel
\begin{equation}
 k^{[2]}(z,z') = 1 + z_<z_> + 2z^2_<z_> - \frac{2}{3}x^3_<,
\end{equation}
is used, where $z_>$ and $z_<$ are analogous to $x_>$ and $x_<$. Here,
the variable $z$ is defined as
\begin{equation}
z = \frac{1 - {\rm cos} \theta}{2},
\end{equation}
so that the values of $z$ are always in the interval [0,1].\\

\noindent
Finally, the 3-dimensional kernel is
\begin{equation}
\label{3dk}
 K({\bf{x}}, {\bf{x}}') = k^{[2,4]}(R,R')k^{[2,4]}(r,r')k^{[2]}(z,z'),
\end{equation}
where, ${\bf{x}}, {\bf{x}}'$ are $(R, r, z)$ and $(R', r', z')$,
respectively. A computationally efficient toolkit is used in this work
to calculate the coefficients and in evaluating the
function\cite{unk17:1923}. Adding a small regularization parameter
(here $\lambda = 10^{-20}$ for the MRCI+Q data) to the diagonal
elements provides additional numerical stability. In practice,
$\lambda$ is increased until a regular solution is obtained for the
inversion. For FCI no regularization is required.\\

\noindent
To represent the long range part of the H$_2^+$--He interaction
the analytical form from Ref. \citen{fal99:117}
 \begin{equation}
 \begin{split}
  V_{\rm long}(R, r, \theta) = & - \frac{\alpha_d q^2}{2R^4} -
  \frac{\alpha_q q^2}{2R^6} - \frac{\alpha_0 q^2}{2R^8} -
  \frac{\beta_{ddq} q^3}{6R^7} - \frac{\gamma_d q^4}{24R^8}
  \\ &-\frac{3\alpha_dq\Theta(r)P_2(\rm{cos}\theta)}{R^6}
  -\frac{5\alpha_dq\Phi(r)P_4(\rm{cos}\theta)}{R^8}
  -\frac{6\alpha_qq\Theta(r)P_2(\rm{cos}\theta)}{R^8}
  \\ &-\frac{C_6^0(r) + C_6^2(r)P_2(\rm{cos}\theta)}{R^6} -
  \frac{C_8^0(r) + C_8^2(r)P_2(\rm{cos}\theta)}{R^8}
  \end{split}
  \label{tl}
 \end{equation}
is used. Here, the first five terms represent the charge+induced
multipole interactions, the sixth term represents the
dipole+quadrupole induction interaction and the seventh and eighth
terms represent the higher order induced-dipole+hexadecapole and
induced-quadrupole+quadrupole interactions, respectively.  Here,
$\Theta(r)$ and $\Phi(r)$ are the quadrupole and hexadecapole moments
of H$_2^+$, respectively. The last two terms in Eq. \ref{tl} are the
contributions from dispersion interactions. The $r-$dependence of the
moments and dispersion coefficients is included by representing them
as a second degree Taylor series. All parameters used in Eq. \ref{tl}
are those from Ref. \citen{fal99:117}. The parameters for He are
given in Table \ref{tab:params}.\\

\noindent
For configurations with $R \gtrsim 9.0$ a$_0$ the MRCI+Q/aug-cc-pV6Z
calculations are discontinuous along the $R-$coordinate, see Figure
\ref{fig1}, where MRCI+Q/aug-cc-pV6Z energies are given fixed values
of $\theta$ and $r$.  For large values of $R$, MRCI+Q/aug-cc-pV6Z
energies are discontinuous which originates from the Davidson
correction of the MRCI energies because the order of the states along
a potential energy scan can swap. This then leads to discontinuities
in the Davidson-corrected energies. Hence, for the long range part of
the PES the explicit analytical long-range expression (see
Eq. \ref{tl}) was used to construct the full 3D PES which is referred
to as MRCI+Q+LR in the following. In order to smoothly connect the
short- and long-range parts of the MRCI+Q PES a Fermi (switching)
function is used, see Figure \ref{fig1a}:
\begin{equation}
  f_s(R) = \frac{1}{{\rm exp}\frac{R-R_0}{\delta R}+1}
\label{eq:fermi}
\end{equation}
where $R_0 = 8.5$ a$_0$ and $\delta R = 0.2$ a$_0$.  The function has
a value of 0.5 at $R = 8.5$ a$_0$. The total potential $(V_{\rm tot})$
is then calculated as
\begin{equation}
 V_{\rm tot} = f_s V_{\rm RKHS} + (1-f_s)V_{\rm long},
\label{sfunc}
\end{equation}
where $V_{\rm RKHS}$ is the short range part of the interaction
potential obtained from RKHS interpolation using the many body
expansion.\\

\begin{figure}
\includegraphics[scale=1.4]{fig1.pdf}
\caption{MRCI+Q/aug-cc-pV6Z and FCI/aug-cc-pV5Z energies as a function
  of $R$ for fixed $\theta = 172.935^\circ$ and $r = 2.0$ a$_0$. The
  black line is the RKHS interpolation of the FCI energies}
\label{fig1}
\end{figure}

\noindent
Full CI calculations are smooth out to $R \sim 50$ a$_0$, contrary to
MRCI+Q, see Figure \ref{fig1}. Hence, the full 3-dimensional PES was
also calculated using FCI using a somewhat smaller basis set, i.e.,
aug-cc-pV5Z. This PES, called FCI in the following, was again
represented as a RKHS. Although the FCI energies are smooth in the
$R-$long range, a third PES (FCI+LR) was constructed by using the same
long range expression used for the MRCI+Q+LR PES. For the FCI+LR PES
the parameter values in the switching function were $R_0 = 13.5$ a$_0$
and $\delta R = 0.25$ a$_0$ in Eq. \ref{eq:fermi}.\\

\begin{figure}
\includegraphics[scale=1.6]{fig2.pdf}
\caption{Upper panel: Energies obtained from MRCI+Q/aug-cc-pV6Z, RKHS,
  analytical long-range and RKHS+analytical are plotted as a function
  of $R$ for fixed $\theta = 172.935^\circ$ and $r = 2.0$ a$_0$. Lower
  panel: the weights for the short range interaction energies and long
  range interaction energies as a function of $R$.}
\label{fig1a}
\end{figure}


\subsection{Bound state calculations}
Ro-vibrational bound state calculations for different $J$ states with
$e$ and $f$ symmetries are carried out in scattering coordinates using
the 3D discrete variable representation (DVR) method with the DVR3D
program suite.\cite{ten04:85} The radial Gauss-Laguerre quadrature
grids consist of 86 and 32 points along $R$ and $r$ coordinates,
respectively. For the Jacobi angle $\theta$, a grid of 36
Gauss-Legendre points was used and for the radial grids ($r$, $R$) the
wavefunctions were constructed using Morse oscillator functions. For
the diatom (H$_2^+$), $r_e = 2.5$ a$_0$, $D_e = 0.1026$ $E_{\rm h}$
and $\omega_e = 0.018$ $E_{\rm h}$ are used and with these parameters
the $r-$grid covered points between 0.92 to 3.8 a$_0$. As the wave
functions for the near-dissociation states need to cover large values
along $R$, the corresponding values were $R_e = 11.5$ a$_0$, $D_e =
0.08$ $E_{\rm h}$, and $\omega_e = 0.00065$ $E_{\rm h}$ which defined
the $R$ grid between 1.82 and 20.87 a$_0$. The $r_2$
embedding\cite{ten04:85} is used to calculate the rotationally excited
states, where the $z-$axis is parallel to $R$ in body-fixed Jacobi
coordinates. For the $J > 0$ calculations, the Coriolis couplings are
included. In the $r_2$ embedding, calculations with $ipar = 1$ and 0
correspond to the $ortho$ and $para$ H$_2^+$, respectively. The $e$
and $f$ symmetries are assigned by the parity operator $p$.\\

\noindent
Another method by which we calculated the bound states is the
coupled-channels variational method (CCVM). It is similar to a
coupled-channels (CC) scattering calculation, but instead of
propagating the radial coordinate $R$ to solve the CC differential
equations it uses a basis also in $R$ and obtains the desired number
of eigenstates of the Hamiltonian matrix with the iterative Davidson
algorithm \cite{davidson:75}. For the angular motion of H$_2^+$ in the
H$_2^+$--He complex we used a free rotor basis with $j_{H_2^+}$
ranging from 0 to 14 (or 16, in tests). The basis in the H$_2^+$
vibrational coordinate $r$ contains the $v = 0 - 7$ eigenfunctions of
the free H$_2^+$ Hamiltonian for $j_{H_2^+} = 0$ on a grid of 110
equidistant points with $r = 0.25 - 5.5~a_0$. The basis in $R$ was
obtained by solving a one-dimensional (1D) eigenvalue problem with the
radial kinetic energy and a potential $V_{\rm eff}(R)$. This potential
is a cut through the full 3D potential of He-H$_2^+$ with $\theta$ and
$r$ fixed at the equilibrium values, to which we added a term linear
in $R$ with a slope that was variationally optimized by using the $R$
basis in full 3D calculations of the lower He-H$_2^+$ levels. The 1D
radial eigenvalue problem was solved with sinc-DVR
\cite{groenenboom:93b} on a 357-point grid with $R = 2 - 50~a_0$. In
order to converge also near-dissociative states we finally included
120 radial basis functions in the 3D full direct product basis.\\

\section{Results and Discussion}

\subsection{Quality of the PESs}
First, the quality of the {\it ab initio} calculations and their RKHS
representation is considered. In Figure \ref{fig2} the analytical
energies are compared with the {\it ab initio} energies for a few
selected Jacobi angles at $r = 2.0$ a$_0$ for the MRCI+Q+LR PES.  A
similar comparison is also shown for the FCI PES in Figure
S1. Excellent agreement between the two sets of data is
found, see Figure \ref{fig2}. Figure S2 presents the
contour plot of the analytical energies for the He-H$_2^+$ system for
$r = 2.0$ a$_0$.\\

\begin{figure}
\includegraphics[scale=0.74]{fig3.pdf}
\caption{The analytical energies obtained from MRCI+Q+LR PES (solid
  lines) and the MRCI+Q/aug-cc-pV6Z {\it ab initio} energies as a
  function of $R$ for several Jacobi angles and at fixed $r = 2.0$
  a$_0$.}
\label{fig2}
\end{figure}

\noindent
The quality of the RKHS representation of the MRCI+Q+LR and FCI PES is
reported in Figure S3. For the grid points used to
generate the RKHS representation, the agreement between reference
points and the reproducing kernel is excellent with $R^2$ values of
$(1-2\times10^{-9})$ and $(1-3\times10^{-12})$ for MRCI+Q+LR and FCI
PESs, respectively. In addition, {\it ab initio} energies were also
calculated at the MRCI+Q/aug-cc-pV6Z and FCI/aug-cc-pV5Z level of
theory for off-grid geometries. They are also reported in Figure
S3 together with the RKHS energies evaluated at these
geometries. Again, the agreement between the electronic structure
calculations and the RKHS representation is good with $R^2$ value of
$(1-2\times10^{-7})$ for the MRCI+Q+LR PES.\\

\noindent
The equilibrium geometry of the FCI/aug-cc-pV5Z surface is a linear
He--H-H configuration ($r = 2.0749$ a$_{0}$ and $R=2.9712$ a$_{0}$),
with an energy of --2735.1 cm$^{-1}$ below the H$_{2}^{+}$
asymptote. This compares with the MRCI+Q+LR calculations for which
$r=2.0745$ a$_{0}$, $R=2.9713$ a$_{0}$ and depth --2736.2 cm$^{-1}$
and the earlier QCISD(T)/aug-cc-pvQz PES\cite{meu99:3418} ($r=2.0750$
a$_{0}$, $R=2.9720$ a$_{0}$ and depth --2717.0 cm$^{-1}$)
values. Hence, the structures of all PESs differ by less than 0.01
a$_0$ but the energetics varies over a range of $\sim 20$ cm$^{-1}$
whereby the dissociation energies for the two PESs from the present
work only differ by 1.1 cm$^{-1}$.\\

\subsection{Bound States}
The ground state energy of H$_2^+$($v=0,j=0)$--He computed from the
MRCI+Q+LR surface for $ortho-$H$_2^+$--He using DVR3D and CCVM are
--1795.1567 and --1795.3328 cm$^{-1}$, respectively. The same
energies, are obtained from the FCI PES as --1793.7632 and --1793.9067
cm$^{-1}$, using DVR3D and CCVM, respectively. These values are $\sim
40$ cm$^{-1}$ lower compared to those reported
previously\cite{meu99:3418} (--1754.269 cm$^{-1}$) on the QCISD(T)
PES. In Ref.\cite{meu99:3418} only 3 quadrature points along $r$ were
used for the 3D bound state calculations, which may not be sufficient
to fully converge the energies. The ground state energy is also
calculated in the present work following a time dependent wave packet
approach\cite{kon16:034303} on a 2D potential fixing $r$ at 2.0
a$_0$. These ground state energies are in fair agreement with previous
results (--1603 vs. --1593 cm$^{-1}$).\cite{meu99:3418}\\

\noindent
For {\it para-} H$_2^+$--He the ground state energies obtained from
the MRCI+Q+LR PES using DVR3D and CCVM are --1795.1575 and --1795.3352
cm$^{-1}$, respectively.  For the FCI surface the ground state
energies of {\it para-} H$_2^+$--He are calculated as --1793.7639 and
--1793.9091 cm$^{-1}$) using DVR3D and CCVM, respectively.  The
difference between DVR3D and CCVM is less than 0.17 cm$^{-1}$ for
both, $ortho$ and $para-$H$_2^+$--He. For the FCI+LR PES all the bound
states obtained from different methods for both, $ortho$ and
$para-$H$_2^+$--He are within 0.02 cm$^{-1}$ or less of the FCI PES
results. Hence, only the results obtained from the MRCI+Q+LR and FCI
PESs are reported.\\


\begin{figure}
\centering
\includegraphics[scale=0.5]{fig4.pdf}
\caption{Bound state energies for $J=0 \dots 6$ for {\it
    ortho-}H$_2^+$-He with $e$ symmetry). Results on the MRCI+Q+LR
  surface (black) and FCI (red) are shown. The ground state energy
  (--1793.7632 cm$^{−1}$) for the FCI/aug-cc-pv5Z surface is marked
  with a dashed line.}
\label{fig:fig4}
\end{figure}


\noindent
A direct comparison for all and the near-dissociative (within 20
cm$^{-1}$ of dissociation) bound $e$ states for {\it
  ortho-}H$_2^+$--He and $J=0$ to $J=6$ from DVR3D calculations and
using the MRCI+Q+LR (red) and FCI (black) PESs is given in Figures
\ref{fig:fig4} and \ref{fig:fig5}. All states up to the dissociation
limit of $v=0, j=1$ state of H$_2^+$ are reported. The level pattern
for the two PESs is nearly identical. The distribution of the energy
difference $\Delta E$ between the MRCI+LR and the FCI PESs for
calculations with DVR3D or CCVM is given in Figure S5.\\

\begin{figure}
    \centering
    \includegraphics[scale=0.50]{fig5.pdf}
    \caption{Bound state energies for ($J = 0 \dots 6$)
      (\textit{ortho}) with $e$ symmetry. With MRCI+Q+LR surface
      (black) and FCI (red). States within $\sim 20$ cm$^{-1}$ of the
      dissociation are reported.}
\label{fig:fig5}
\end{figure}

\noindent
The transitions that were probed by the microwave experiments lie
close to dissociation. Hence, a particular focus here is on accurately
computing these stationary states and to determine whether any
candidate transitions can be identified from using the MRCI+Q+LR and
the FCI PESs. A tentative assignment in particular for the 15.2 GHz
and 21.8 GHz transitions has been given previously based on
experiments using electric field dissociation.\cite{gam02:6072} They
were analyzed using an effective Hamiltonian. The 15.2 GHz transition
was assigned to a low-$N$ transition (in the terminology of
Ref.\cite{gam02:6072}, $N$ is the spin-free angular momentum which is
$J$ in the present work) with $\Delta N = 0$ with $N=3$ or $N=4$ in
{\it ortho-}H$_2^+$--He. In the following, $N$ is used when referring
to the analysis of the experiments\cite{gam02:6072} whereas $J$ is
used when discussing the present calculations. The fine and hyperfine
splittings due to coupling of electron and total nuclear spin, and
coupling of the resultant to the rotational angular momentum of the
nuclei are both less than 100 MHz, so are several orders of magnitude
smaller than the separations between rotational levels of the
complex. Thus, identification of $N$ with $J$ is a meaningful
approximation.\\

\noindent
For the 21.7 GHz transition on the other hand the analysis led to an
assignment involving $\Delta N = 1$ with $N = 11$ and $N' = 10$ in
{\it para-}H$_2^+$--He. While the analysis leading to a $\Delta N = 1$
transition involving {\it para-} H$_2^+$ is based on physical grounds,
that to a high-$N$ state involves fitting of the Zeeman pattern which
is more approximate. The selection rules for these transitions are $e
\leftrightarrow f$ for $\Delta J = 0$ and $e \leftrightarrow e$ or $f
\leftrightarrow f$ for $\Delta J = \pm 1$, respectively.\\

\noindent
First, the near-dissociative states for {\it ortho-}H$_2^+$--He are
discussed. All near-dissociative states from the MRCI+Q+LR and FCI
PESs using the DVR3D and CCVM methods are reported in Figures
\ref{fig:fig6}a to d. All energies for $J=0$ to 6 are also reported in
Tables \ref{tab:mrci} and \ref{tab:fci}. There is one $e/f$ parity
doublet with $\Delta J = 0$ with a transition frequency between 10 and
18 GHz, involving the $J=2$ state for {\it ortho-}H$_2^+$--He. Using
DVR3D the transition frequency is 14.4(5) GHz whereas with the CCVM
code the transition is at 9.3 GHz. The parity doublet in both cases is
within 2 cm$^{-1}$ of dissociation which makes it a pair of
near-dissociative states. This is also confirmed by considering the
expectation value for the $R-$coordinate for the two states involved
which are $\langle R \rangle = 13.1$ a$_0$ for the $e$ state and
$\langle R \rangle = 12.7$ a$_0$ for the $f$ state which confirms
their long range character as suggested from the experiments.\\

\begin{figure}
\centering
\includegraphics[scale=0.75]{fig6.pdf}
\caption{Near dissociation {\it ortho-}H$_2^+$--He states (in
  cm$^{-1}$) and predicted transition frequencies in (GHz) using the
  MRCI+Q+LR (left, panels a and c) and FCI/aug-cc-pv5 (right, panels b
  and d) PESs. States with {\it e} (red) and {\it f} symmetry (black)
  are reported separately. Results from DVR3D and CCVM are in the top
  and bottom row, respectively. The 14.4(5) and 9.3 GHz parity doublet
  is a candidate for the 15.2 GHz line observed experimentally which
  had been assigned to a parity doublet with $\Delta N =
  0$.\cite{gam02:6072}}
\label{fig:fig6}
\end{figure}

\noindent
Next, the near-dissociative states for {\it para-}H$_2^+$--He are
discussed for the two PESs and the two methods to compute bound
states, see Figures \ref{fig:fig7}a to d. The only near-dissociative
states involving either an $e/e$ or an $f/f$ transition with a
transition frequency around 20.1(2) and 16.7(6) GHz from DVR3D and
CCVM involves a $J=0$ and a $J=1$ state. The $\langle R \rangle =
15.0$ a$_0$ for the $J = 0, e$ state and $\langle R \rangle = 16.0$
a$_0$ for the $J = 1, f$ state, show the long range nature of the wave
functions. The potential candidate for the 21.8 GHz transition is not
found for the high $J$ states sufficiently close to dissociation to be
part of a suitable candidate transition.\\


\begin{figure}
\centering
\includegraphics[scale=0.75]{fig7.pdf}
\caption{Near dissociation {\it para-}H$_2^+$--He states (in
  cm$^{-1}$) and predicted transition frequencies in (GHz) using the
  MRCI+Q+LR (left, panels a and c) and FCI/aug-cc-pv5 (right, panels b
  and d) PESs. States with {\it e} (red) and {\it f} symmetry (black)
  are reported separately. Results from DVR3D and CCVM are in the top
  and bottom row, respectively.}
\label{fig:fig7}
\end{figure}


\begin{table}[ht]
\caption{Near dissociation states calculated using MRCI surface with
  DVR3D and CCVM for He-H$_2^+$ in cm$^{-1}$. Zero is set to the
  energy of H$_2^+$($v=0,j=0$) state. The dissociation limit for
  $ortho$ is at 58.2336 cm$^{-1}$ and 58.3124 cm$^{-1}$ for the DVR3D
  and CCVM respectively, which corresponds to H$_2^+$($v=0,j=1$),
  whereas the dissociation limit for {\it para-} is at 0 cm$^{-1}$
  which corresponds to H$_2^+$($v=0,j=0$).}
\begin{tabular}{c|cc|cc|rr|rr}
\hline
\hline
    \multicolumn{5}{c}{\it ortho}& \multicolumn{4}{c}{\it para}\\
    \hline
    \multicolumn{3}{c}{$e$} & \multicolumn{2}{c}{$f$}& \multicolumn{2}{c}{$e$} & \multicolumn{2}{c}{$f$} \\
    \hline
    $J$ & CCVM   &  DVR3D  & CCVM  &  DVR3D  & CCVM    &    DVR3D &   CCVM  & DVR3D \\
   \hline
    0 & 39.896 & 40.321 &        &        & --15.641 & --15.719 \\
      & 55.689 & 55.876 &        &        &  --1.509 &  --1.137 \\
    1 & 43.199 & 43.668 & 54.629 & 54.464 & --13.965 & --14.035 \\
      & 54.514 & 54.582 &        &        &  --0.954 &  --0.430 \\
      & 56.959 & 57.273 &        &        &          &           \\
    2 & 49.345 & 49.755 & 56.329 & 56.439 & --10.712 & --10.766 \\
      & 56.639 & 56.923 &        &        &  --0.061 &    0.895 \\
      & 58.272 & 59.186 &        &        &          &           \\
    3 & 44.938 & 44.941 & 41.798 & 41.783 &  --6.082 & --6.109   \\
      & 55.443 & 55.588 &        &        &          &           \\
    4 & 46.022 & 45.971 & 50.718 & 50.712 &  --0.587 & --0.540   \\
      & 54.375 & 54.375 &        &        &          &           \\
    5 & 40.655 & 40.711 & 40.682 & 40.680 &  --7.877 & --7.799   \\
      & 55.518 & 55.507 &        &        &          &           \\
    6 & 55.985 & 55.959 &        &        &          &            & --13.897 & --13.746 \\
\hline
\hline
\end{tabular}
\label{tab:mrci}
\end{table}

\begin{table}[ht]
\caption{Near dissociation states calculated using FCI surface with
  DVR3D and CCVM for He-H$_2^+$ in cm$^{-1}$. Zero is set to the
  energy of H$_2^+$($v=0,j=0$) state. The dissociation limit for
  $ortho$ is at 58.2320 cm$^{-1}$ and 58.3124 cm$^{-1}$ for the DVR3D
  and CCVM respectively, which corresponds to H$_2^+$($v=0,j=1$),
  whereas the dissociation limit for $para$ is at 0 cm$^{-1}$ which
  corresponds to H$_2^+$($v=0,j=0$).}
  \begin{tabular}{c|cc|cc|rr|rr}
    \hline
    \hline
     \multicolumn{5}{c}{$ortho$}& \multicolumn{4}{c}{$para$}\\
    \hline
     \multicolumn{3}{c}{$e$} & \multicolumn{2}{c}{$f$}& \multicolumn{2}{c}{$e$} & \multicolumn{2}{c}{$f$} \\
    \hline
    \hline
      $J$ & CCVM   &  DVR3D & CCVM    & DVR3D  & CCVM      &  DVR3D \ & CCVM  & DVR3D\\
   \hline
    0 & 39.575 & 40.036 &         &        & --16.005 & --16.067\\
      & 55.662 & 55.860 &         &        &  --1.527 &  --1.150\\
    1 & 42.934 & 43.438 & 54.601  & 54.648 & --14.315 & --14.368\\
      & 54.494 & 54.575 &         &        &  --0.969 &  --0.448\\
      & 56.916 & 57.232 &         &        &          &         \\
    2 & 49.213 & 49.647 & 56.314  &  56.437 & --11.021 & --11.058\\
      & 56.623 & 56.916 &         &        &  --0.061 &    0.881\\
      & 58.271 & 59.162 &         &        &          &         \\
    3 & 44.497 & 44.513 & 41.473  & 41.473 & --6.311  &  --6.320\\
      & 55.422 & 55.577 &         &        &          &         \\
    4 & 45.652 & 45.611 & 50.538  & 50.547 & --0.680  &  --0.619\\
      & 54.067 & 54.084 &         &        &          &         \\
    5 & 40.565 & 40.626 & 40.854  & 40.860 & --8.092  &  --8.005\\
      & 55.219 & 55.220 &         &        &          & \\
    6 & 55.626 & 56.610 &         &        &          &  & --11.890  & --11.742\\
    \hline
    \hline
  \end{tabular}
  \label{tab:fci}
\end{table}


\section{Conclusions}
Two new PESs at the MRCI+Q+LR and FCI level of theory with large basis
sets and represented as a reproducing kernel have been used to
determine all bound and near-dissociative states for {\it ortho-} and
{\it para-}H$_2^+$--He. It is found that at both levels of theory the
bound states compare to within fractions of a wavenumber when
stationary states are determined from the same nuclear quantum
code. Moreover, the stationary states on one and the same PES
determined from two different quantum bound state codes (DVR3D and
CCVM) also agree closely, typically within less than fractions of one
cm$^{-1}$. This provides stringent benchmarks on potential transitions
that have been observed experimentally. One such assignment is for the
15.2 GHz transition which corresponds to {\it ortho-}H$_2^+$--He. The
transition found in the present work involves an $e/f$ parity doublet
with $J=2$. This compares with a tentative assignment to an $e/f$
parity doublet involving either a $J=3$ or $J=4$ state. For the 21.8
GHz transition which had been tentatively assigned to an $e/e$ or
$f/f$ transition in {\it para-}H$_2^{+}$--He the candidate,
near-dissociation states are $J=0$ and $J=1$, both of which are within
less than 2 cm$^{-1}$ of dissociation and the transition frequencies
range from 16 to 22 GHz. However, no high $J-$ candidate states
suitable for transition were found.\\

\noindent
The present work presents two high-accuracy, fully dimensional PESs
for H$_2^+$--He together with quantum bound state calculations that
provide potential assignments of experimentally characterized,
near-dissociation states. The results from both, MRCI+LR and full CI
PESs, using two different approaches for calculating the quantum bound
states are largely consistent. It will be interesting to use the
present PESs in future inelastic scattering calculations.\\

\section*{Conflicts of interest}
There are no conflicts of interest to declare.

\section*{Acknowledgments}
This work was supported by the Swiss National Science Foundation
through grants 200021-117810, the NCCR MUST and the AFOSR (to MM). The
authors thank Prof. J. Tennyson for exchange on the DVR3D program.



