\section{Introduction}\label{section:1-Introduction}
Condensation can occur when one or more condensable gas species are exposed to temperature differences around their dewpoint temperature.
In particular, environments connected to the atmosphere at standard ambient conditions are subject to the condensation of water vapor within the volume or at comparatively cool surfaces.
Studying the interaction between condensation and flow is therefore highly relevant to a wide range of applications, from meteorological flows to residential and vehicular ventilation.

Phase transitions in turbulent flows connect a variety of effects concerning both fluid dynamics and thermodynamics in the interplay between both phases.
While correlations and models can describe the heat transfer on a global level for specific fluid mixtures and parameter ranges\cite{Herranz1998, Bonner2013}, inclusion of phase transition in fully resolved simulations is possible as well.
The interaction can be realized via a two-way coupling, where a carrier phase determines the trajectories of the dispersed particles of the condensed phase, and the mass loading due to the particles in turn influences the flow of the carrier fluid. Another option is a four-way coupling, where particle-particle interactions are included in addition to the two-way coupling\cite{Kuerten2015}.
The particle motions can be calculated using Lagrangian methods \cite{Mashayek1998} or Eulerian models \cite{Fevrier2005,Masi2010}.
In these systems, the high computational costs associated with the simulation of the dispersed phase are unavoidable in order to capture the full dynamics of the flow.

In flows at ambient conditions in systems open to the atmosphere, the vapor concentration is typically small compared to the concentrations found in technical applications.
Additionally, the degree of subcooling is expected to be small and limited to regions in the immediate vicinity of cooled surfaces, where condensation therefore occurs predominantly.
Motivated by the expected low condensation mass flux, this study employs a single-phase approach to investigate the influence of the release of latent heat during condensation on the fundamental structure of turbulent flow in isolation from other phenomena of multiphase flow.

This approach extends the study by Russo \emph{et al.}\cite{Russo2013}, which compares differentially heated turbulent channel flow in three cases: with phase transition and droplets in a two-way coupled multiphase simulation, with inert particles, and without either phase transition or particles.
Here, a fourth configuration is investigated, where phase transition effects are included, but no droplets are simulated.

The condensable species is modeled by a scalar concentration field which is transported through the system by convection and diffusion.
The simulation can be considered as multiphase only in so far as the phase transition process is modeled, changing the density of the fluid as a result of the sudden change in temperature and concentration.
Consequently, the coupling with the flow is realized not by direct exchange of momentum but via the change in fluid density and the resulting buoyant forces.

This simplified approach allows the use of fully resolved direct numerical simulations (DNS) to simulate flows including phase transition effects while incurring only a small penalty in computational cost compared to DNS of mixed convection without phase transition.
The use of DNS allows an undisturbed evaluation of the simplified condensation modeling approach since the solution of the flow equations is free from the influence of turbulence or subgrid models, thereby isolating the effects caused by the release of latent heat at the phase boundary.

Building upon the well-understood system of turbulent flow through a differentially heated vertical channel geometry\cite{Kasagi1997}, which is extended to include an additional active scalar field representing the vapor concentration\cite{Bahavar2019}, this approach allows the analysis of modifications of turbulence due to the phase transition within the framework of turbulent flow in the presence of buoyancy gradients.
\section{Governing equations}\label{section:2-Governing_equations}
The flow of the gas mixture is simulated by directly solving the incompressible Navier-Stokes equations for the carrier phase,
\begin{align}
\nabla\cdot \mathbf{u}=0,\label{eq:continuity}\\
\frac{\partial\mathbf{u}}{\partial t} + (\mathbf{u}\cdot\nabla )\mathbf{u} = -\frac{1}{\rho}\nabla p + \nu \nabla^2\mathbf{u}-\mathbf{B},\label{eq:momentum}
\end{align}
where $\mathbf{u}=(u_x,u_y,u_z)$ is the vector field describing the fluid velocity as a function of space and time, $\nu$ and $\rho$ represent the kinematic viscosity and density of the fluid, and $\nabla p$ is the pressure gradient driving the flow.
The evolution of the temperature $T$ and vapor concentration $c$, expressed in terms of molar fraction, are modeled with convection-diffusion equations,
\begin{align}
&\frac{\partial T}{\partial t} + \mathbf{u}\cdot\nabla T = \kappa\nabla^2 T + \frac{h_v}{c_p} f(T,c),\label{eq:cd_temperature}\\
&\frac{\partial c}{\partial t} + \mathbf{u}\cdot\nabla c = D\nabla^2 c - f(T,c),\label{eq:cd_concentration}
\end{align}
with the thermal diffusivity $\kappa$ and the mass diffusivity $D$.
The inclusion of these scalar fields into the simulation gives rise to buoyant forces $\mathbf{B}$ due to density changes induced by temperature and concentration differences within the system.
The buoyancy contributions are linearized within the framework of the Boussinesq approximation\cite{Gray1976},
\begin{align}
\mathbf{B} = \beta_T (T-T_{\mathit{ref}})\mathbf{g} + \beta_c (c-c_{\mathit{ref}})\mathbf{g}.\label{eq:buoyancy}
\end{align}
Here $\mathbf{g}$ is the gravitational acceleration and $\beta_T$ and $\beta_c$ refer to the expansion coefficients with respect to temperature and vapor concentration \cite{Bird2015},
\begin{align}
\beta_{T}=\frac{1}{\rho_{\mathit{ref}}}\left.\frac{\partial\rho}{\partial T}\right\rvert_{T_{\mathit{ref}}},\quad\beta_{c}=\frac{1}{\rho_{\mathit{ref}}}\left.\frac{\partial\rho}{\partial c}\right\rvert_{c_{\mathit{ref}}}.\label{eq:expansion_coeffs}
\end{align}
$T_{\mathit{ref}}$ and $c_{\mathit{ref}}$ are set reference values about which the density variation is calculated\cite{Barletta1999}, and $\rho_{\mathit{ref}}=\rho(T_{\mathit{ref}},c_{\mathit{ref}})$ is the reference density.

The coupling between temperature and concentration caused by phase transition effects is represented by the source term $f$ in equation (\ref{eq:cd_concentration}), which describes the loss of vapor concentration due to condensation.
Conversely, the temperature source term in equation (\ref{eq:cd_temperature}) differs by its sign and the ratio of latent heat of condensation to specific heat capacity of the fluid, $h_v/c_p$, thereby yielding the increase in temperature associated with the release of latent heat at the phase transition.
The mass flux $\dot{m}^\ast$ across the phase boundary during condensation is given by the simplified Hertz-Knudsen-Schrage equation\cite{Marek2001},
\begin{align}
\dot{m}^\ast = \frac{2\sigma_c}{2-\sigma_c}\sqrt{\frac{M}{2\pi RT}}(p_v - p_{\mathit{sat}}),\label{eq:HKS_mass_flux}
\end{align}
with the molar weight of the vapor $M$ and universal gas constant $R$.
$\sigma_c$ is the condensation accommodation coefficient, describing the probability of a vapor molecule remaining in the liquid phase after impinging on the physical boundary between the phases.
$p_v$ and $p_{\mathit{sat}}$ are the partial vapor pressure and saturation pressure, respectively, related to the vapor concentration via $c=p_v/p_{\mathit{total}}$.
In terms of the concentration, the rate of change in vapor content $f$ in a volume $V$ across a boundary of area $A$ is then given by
\begin{align}
f = \frac{2\sigma_c}{2-\sigma_c}\sqrt{\frac{RT}{2\pi M}}\frac{A}{V}(c - c_{\mathit{sat}}).\label{eq:HKS_vapor_source}
\end{align}
Vapor removed from the system by the source term is not considered further.
The liquid phase is not simulated in this single-phase approach.

The interdependence of the scalar transport equations via the phase change source terms in addition to the coupling between the equations for the scalar fields and the velocity field via convective transport on the one hand and buoyancy on the other completes the full coupling between the governing equations.
\section{Validation}
The validity of the numerical scheme for solving the governing equations in a simulation including scalar transport and buoyant forces along the streamwise direction is tested in a generic biperiodic channel geometry.

The finite volume DNS is performed using \mbox{OpenFOAM}, with second-order central differencing in space and explicit second-order accurate leapfrog-Euler time integration\cite{Kath2016}.
The projection method is used to determine a pressure field which corrects the initial solution for the velocity field by removing divergences\cite{Chorin1968}.

To solve the governing equations, the channel geometry is discretized into $N_x \times N_y \times N_z = 396 \times 180 \times 316$ hexahedral cells.
These cells are uniformly distributed in the streamwise ($x$) and spanwise ($z$) direction, with a resolution of $\Delta x^+=6$ and $\Delta z^+=3$, given in wall units $\nu/u_{\tau}$.
Along the wall-normal ($y$) direction, the cells are distributed following a hyperbolic tangent function, increasing the number of cells near the walls, where high gradients need to be resolved, while limiting the number of cells located in the bulk to reduce computational costs and at the same time minimizing adverse effects due to cell-to-cell stretching\cite{Vinokur1983}.
The resulting spacing ranges from $\Delta y^+=0.2$ to $3.7$.

To establish the accuracy of the DNS solver, a channel without the vapor concentration field including only the temperature is simulated.
Isothermal boundaries are applied at the walls, with temperatures set to $T_h$ at the heated and $T_c$ at the cooled wall.
No-slip and impermeability boundary conditions are imposed on the velocity field.

The flow is characterized by a bulk Reynolds number
\begin{align}
\mathrm{Re} = \frac{u_b \delta}{\nu} = 2280\label{eq:Re_validation}
\end{align}
with the bulk flow velocity $u_b$ and channel half-width $\delta$.
Buoyant forces are quantified by the Grashof number
\begin{align}
\mathrm{Gr} = \frac{g\delta^3\beta_T\Delta T}{\nu^2} = 120000\label{eq:Grashof_thermal_validation}
\end{align}
with the characteristic temperature difference \mbox{$\Delta T = T_h - T_c$.}
This results in a flow which is primarily determined by forced convection with a Richardson number of
\begin{align}
\mathrm{Ri}=\frac{\mathrm{Gr}}{\mathrm{Re}^2}=0.023.\label{eq:Richardson_thermal_validation}
\end{align}
The temperature diffusion is given by the Prandtl number $\mathrm{Pr}=\kappa/\nu=0.71$, corresponding to dry air.

\begin{figure}[t]
\centering
\includegraphics[scale=1]{fig_validation_ux-eps-converted-to.pdf}
\caption{Streamwise velocity profile obtained by DNS compared to reference data.}\label{fig:validation_ux}
\end{figure}
\begin{figure}[t]
\centering
\includegraphics[scale=1]{fig_validation_ux_rms-eps-converted-to.pdf}
\caption{Profile of streamwise velocity fluctuations normalized with $u_\tau$ compared to reference data.}\label{fig:validation_ux_rms}
\end{figure}
\begin{figure}[t]
\centering
\includegraphics[scale=1]{fig_validation_theta-eps-converted-to.pdf}
\caption{Temperature profile obtained by DNS compared to reference data.}\label{fig:validation_T}
\end{figure}
\begin{figure}[t]
\centering
\vspace*{-1.6mm}
\includegraphics[scale=1]{fig_validation_theta_rms-eps-converted-to.pdf}
\caption{Profile of temperature fluctuations obtained by DNS compared to reference data.}\label{fig:validation_T_rms}
\end{figure}

The flow parameters are chosen to allow comparison with an established reference DNS for differentially heated turbulent flow in a vertical channel geometry\cite{Kasagi1997}.
For statistical analysis, averaging is performed in time over an interval of $\Delta t^+ = 30$, with the dimensionless time unit $t^+=t u_{\tau}/\nu$, and along the homogeneous spatial directions.
Averaged quantities are denoted by $\langle\cdot\rangle$.

The nondimensional normalized temperature
\begin{align}
\theta = \frac{T-T_{\mathit{ref}}}{\Delta T}\label{eq:def_theta}
\end{align}
is introduced, with $T_{\mathit{ref}}=(T_h + T_c)/2$ used as the reference for this configuration.

Figures \ref{fig:validation_ux} and \ref{fig:validation_ux_rms} show the profile of the average streamwise velocity $u_x$ and the corresponding average fluctuations $u_x^{\mathit{rms}}$.
Very good agreement with the reference is found for the mean velocity at all points.
The fluctuations of the streamwise velocity are fully captured in the highly turbulent regions near the channel walls, but underestimated in the bulk region of the flow.

Analogously, Figures \ref{fig:validation_T} and \ref{fig:validation_T_rms} show the mean temperature and and its fluctuations.
Again, the agreement between the results obtained here and the reference is good for the mean field, with slight deficits in the magnitude of the fluctuations in the same regions as observed for the velocity.
This suggests that the error in the temperature fluctuations is a result of reduced mixing due to the underestimated turbulence intensity.
This discrepancy is expected due to the fundamental differences between the second-order finite volume approach employed in this study compared to the spectral method used for the reference case.

In addition to assessing solver accuracy, the biperiodic channel setup can be used to investigate the applicability of the Boussinesq approximation in the case of added solutal buoyancy.
To this end, uniform concentration boundary conditions are applied at the channel walls such that the buoyant forces caused by the resulting concentration gradients act in the same direction as the thermal buoyancy.
The solutal Grashof number is
\begin{align}
\mathrm{Gr}_c = \frac{g\delta^3\beta_c\Delta c}{\nu^2} = 248000\label{eq:Grashof_solutal_validation}
\end{align}
with $\Delta c$ describing the difference in concentration between the values at the walls.
The Richardson number quantifying the combined influence of solutal and thermal buoyancy is then $\mathrm{Ri}=0.07$.
The diffusive transport of the concentration field is quantified by the Schmidt number, $\mathrm{Sc}=D/\nu=0.48$.

In this configuration, the high vapor load required to cause high solutal buoyancy means that the assumption of constant fluid properties apart from the linearized density variation is not well justified\cite{Bahavar2019}.
Particularly, the coefficients for scalar transport are sensitive to the mixing ratio between carrier fluid and vapor, varying by more than 10\% across the range of $\Delta c$.
To investigate the inaccuracies introduced by applying the Boussinesq approximation regardless of these caveats, simulations were performed using an extended approximation including linearized changes of the thermal diffusivity $\kappa(c) = \beta_{c,\kappa}\,(c-c_{\mathit{ref}})$, analogously to the treatment of the density in equation (\ref{eq:buoyancy}).
The results obtained using this approach are compared to a simulation using the strict Boussinesq approximation, where the thermal diffusivity is fixed.

\begin{figure}[t]
\centering
\includegraphics[scale=1]{fig_boussinesq_ux-eps-converted-to.pdf}
\caption{Average streamwise velocity profile for channel flow with combined thermal and solutal buoyancy, comparison of strict and extended approximation.The black line shows the profile for the system without solutal buoyancy for comparison.}\label{fig:boussinesq_ux}
\end{figure}
\begin{figure}[t]
\centering
\includegraphics[scale=1]{fig_boussinesq_T-eps-converted-to.pdf}
\caption{Average temperature profile for channel flow with combined thermal and solutal buoyancy, comparison of strict and extended approximation.The black line shows the profile for the system without solutal buoyancy for comparison.}\label{fig:boussinesq_T}
\end{figure}

Figures \ref{fig:boussinesq_ux} and \ref{fig:boussinesq_T} show the resulting velocity and temperature profiles compared to the previous validation case without solutal buoyancy.
The effect of the increased Grashof number is clearly visible, resulting in a stronger asymmetry of the flow caused by aiding and opposing forces acting on the fluid and in steeper temperature gradients at the walls.
While no difference between the strict approximation and the extended formulation can be observed for the streamwise velocity, the influence of the varying thermal diffusivity is clearly visible in the temperature field.
The reference value was chosen as $\kappa_{\mathit{ref}}=\kappa(T_c)$, such that the differences are most pronounced directly at the heated boundary.
Nevertheless, the effects are smaller than the direct impact of the added buoyancy even at these high vapor loads.

Together, these simulations of biperiodic channel flow establish the viability of the chosen numerical approach and the range of applicability for the Boussinesq approximation employed in the governing equations.
They provide a solid starting point for the investigation of the effects of condensation in a comparable geometry.

On this basis, the active scalar approach to simulate flows including phase transition is compared to a subset of the simulations presented by Russo \emph{et al.}\cite{Russo2014}.
Again, a biperiodic channel geometry is considered.
Constant heat flux boundary conditions of $\pm 32\,\mathrm{W/m^2}$ are applied at the walls, such that the net heat flux into the channel is zero.
Flow with bulk Reynolds number $\mathrm{Re}=2333$, corresponding to a friction Reynolds number $\mathrm{Re_{\tau}}=150$, is simulated.
No gravity and consequently no buoyancy are considered.
At the start of the simulation, the fluid temperature is set to $T_0 = 293.15\,\mathrm{K}$, and the vapor concentration to saturation, $c_0 = c_{\mathit{sat}}(T_0)$.
The dimensionless flow parameters are set according to the values for humid air at such conditions, with Prandtl number $\mathrm{Pr}=0.74$ and Schmidt number $\mathrm{Sc}=0.63$.

Instead of an initial droplet distribution, a scalar field representing a continuous reservoir of liquid water dispersed across the whole channel volume is added .
This allows evaporation of liquid if undersaturation conditions are met, mirroring the treatment of condensation discussed previously.
After temperature and vapor concentration reached a statistically steady state, averaging was again performed in time and along the homogeneous directions.
\begin{figure}[htb]
\centering
\includegraphics[scale=1]{fig_yT_model_validation-eps-converted-to.pdf}
\caption{Average temperature profiles for four cases. The results from the approach taken in this work, including no particles, but with phase transition (PT), are compared to three cases from Russo \emph{et al.}\cite{Russo2014} covering the remaining combinations.}\label{fig:yT_model_validation}
\end{figure}

Russo \emph{et al.} performed simulations for cases with phase transition (PT) and without phase transition (nPT).
In both cases, particles were included in the investigation, liquid droplets for PT, and inert particles for nPT.
Additionally, a simulation without phase transition and particles was presented as a basis for comparison.

The approach taken in this work completes the possible combinations by including phase transition, but not simulating liquid droplets.

Figure \ref{fig:yT_model_validation} shows the average temperature profiles for the four different cases.
The difference between the case with no particles + nPT and with particles + PT illustrates the modification of the flow due to phase transitions obtained by a two-way coupled multiphase simulation.
These modifications are reflected to a large degree by the temperature profile obtained using the active scalar approach in the case without particles, PT.
The reduction of the temperature difference between the walls due to the latent heat released and absorbed in condensation and evaporation is captured, as well as the reduced temperature gradients.
Discrepancies exist primarily near the walls, where the droplet concentration is high due to turbophoresis\cite{Reeks1983} and the contribution of the droplets to the heat transfer is maximized.

Comparing the active scalar approach to the case with particles + nPT shows that including latent heat and neglecting droplets results a much closer approximation of the full multiphase effects than the converse approach.

Based on the comparison to the two-way coupled multiphase simulation, the simplified approach can be employed for the specific case considered in this study, with parameters reflecting flows at ambient conditions over cool surfaces.
Here, condensation will be primarily constrained to the surface, with condensate adhering to the surface. In this case, the influence of dispersed droplets following the flow is much smaller than in the periodic channel discussed above.
\section{Investigation setup}
\begin{figure}[htb]
\centering
\includegraphics[width=0.6\linewidth]{Bahavar_figure-domains-figure0-eps-converted-to.pdf}
\caption{Schematic overview of the coupled simulation domains with velocity mapping at the interface and temperature boundary conditions in the primary domain.}\label{fig:domains}
\end{figure}

The governing equations outlined and tested in the preceding sections allow the investigation of flows exposed to temperature gradients causing oversaturation conditions and, consequently, condensation.
A differentially cooled vertical channel geometry is considered in an inlet-outlet configuration to study the impact of condensation on the turbulent flow in such a setup.
The aperiodic setup limits the extent of the thermal boundary layer, ensuring that condensation occurs near the cooled surface and the error due to neglecting droplets in the flow is minimized.
Additionally, the progression of cooling and condensation along the channel length can be investigated in such a geometry.

Values for temperature $T_{\mathit{in}}$ and vapor concentration $c_{\mathit{in}}$ are prescribed at the inlet.
The channel walls are kept at constant temperature, with one wall at $T_h=T_{\mathit{in}}$, while the opposite wall is cooled with respect to the inlet temperature, $T_c<T_{\mathit{in}}$.
By additionally ensuring that the temperature of the cooled wall is below the dewpoint $T_{\mathit{dp}}$ of the fluid at the inlet, which is in turn below the inlet temperature, $T_c<T_{\mathit{dp}}(c_{\mathit{in}})<T_{\mathit{in}}$, a subcooled region can develop within the channel as a consequence of the cooling influence of the wall.
At the opposite wall, setting the temperature boundary to $T_{\mathit{in}}$ prevents the formation of a thermal boundary layer and, consequently, condensation.
Zero-gradient boundary conditions are enforced for the concentration field, letting vapor concentration at the walls evolve in tandem with the field inside the channel.

A precursor simulation is coupled to the system to provide a velocity inlet boundary condition consistent with fully developed turbulent flow\cite{Akselvoll1993,Lund1998,Bellec2017}.
The precursor simulation consists of a separate channel geometry with cyclic boundary conditions in streamwise direction without temperature and concentration fields.
To sustain the flow in the cyclic domain, a self-correcting uniform global pressure gradient is applied, ensuring a constant prescribed volume flux across the interface\cite{Kath2016}.
In addition to feeding back into the inlet of the cyclic domain, the velocity field at the outlet of the precursor simulation is mapped to the inlet of the primary domain.
This effectively creates an infinite-length inflow region where isothermal turbulence can fully develop before being continuously fed into the differentially cooled channel setup.

\begin{figure}[h]
\includegraphics[scale=1]{coupling-figure0}
\caption{Mapping of the velocity field from the precursor outlet back to its inlet (top) and forward towards the primary domain (bottom) visualized via continuous streaks of $u_x^\prime=2u_{\tau}$. The precursor region is replicated again in grey to illustrate the continuity of the streaks.}\label{fig:coupled_u}
\end{figure}

A schematic overview of the simulation domains is shown in Figure \ref{fig:domains}. The resulting flow setup is illustrated in Figure \ref{fig:coupled_u}.
A top view of isosurfaces of the fluctuating streamwise velocity $u_x^\prime=2u_{\tau}$ is drawn in both computational domains.
The top panel shows the continuity of the velocity fields between the cyclic inlet-outlet-planes of the precursor domain, which is repeated for reasons of visualization.
In the bottom panel, the undisturbed transfer of fluid velocity from the precursor to the primary domain is demonstrated in terms of the isosurfaces of the streamwise velocity.

At the channel walls, impermeability and no-slip boundary conditions are applied, while a standard outflow condition is set at the channel outlet.
The pressure boundaries are specified to complement the choices for the velocity field, with a fixed reference pressure set at the outlet and von-Neumann boundary conditions at the inlet.

Flow parameters are chosen to represent a binary mixture of dry air with trace amounts of water vapor, such as expected at standard atmospheric pressure at room temperature.
The diffusive transport of the scalar fields is characterized by the Prandtl number $\mathrm{Pr}=\kappa/\nu = 0.73$ for the temperature and the Schmidt number $\mathrm{Sc}=D/\nu=0.65$ for the vapor concentration reflecting humid air.
Gravity acts along the streamwise direction, causing aiding and opposing buoyant forces to act on the fluid.
Again, the separate contributions due to differences in temperature and vapor concentration can be quantified using the thermal Grashof number $\mathrm{Gr}_T = 38000$ and the solutal Grashof number $\mathrm{Gr}_c=1500$.
Here the characteristic temperature difference $\Delta T = T_{\mathit{in}} - T_c$ refers to the maximal span of temperatures possible within the system, and analogously $\Delta c = c_{\mathit{in}} - c_{\mathit{sat}}(T_c)$ to the corresponding range of concentration.

The bulk Reynolds number is $\mathrm{Re}=1994$, resulting in a friction Reynolds number of
\begin{align}
\mathrm{Re}_{\tau}=\frac{u_{\tau}\delta}{\nu}=135,\label{eq:Reynolds_tau}
\end{align}
where $u_{\tau}$ refers to the friction velocity.
Turbulent flow is expected for channel flow at these values of $\mathrm{Re}$.

The relation between the release of latent heat during condensation and the heating of the system as a consequence is expressed by the Jakob number,
\begin{align}
\mathrm{Ja}=\frac{c_p}{h_v}\Delta T = 0.012.
\end{align}

The combined Richardson number is
\begin{align}
\mathrm{Ri}=\frac{\mathrm{Gr}_T + \mathrm{Gr}_c}{\mathrm{Re}^2} = 0.01.\label{eq:Richardson}
\end{align}
The flow is therefore dominated by forced convection, with only small contributions due to the buoyant forces.
In particular, the Richardson number is much smaller than in the case for which the applicability of the strict Boussinesq approximation was examined ($\mathrm{Ri}=0.07$).
Therefore, using the strict approximation and considering all fluid properties apart from the density as constant is justified going forward.

Further, since the solutal Grashof number is smaller than the thermal Grashof number by an order of magnitude, the direct influence of the solutal buoyancy on the flow is negligible.
On the other hand, the release of latent heat means that a change in thermal buoyancy is leveraged by the concentration change due to condensation.
Given a concentration change $\Delta c^\prime$ due to condensation, the resulting solutal buoyant force is
\begin{align}
F_{\Delta c^\prime}=\beta_c \Delta c^\prime g,\label{eq:force_vapor_diff}
\end{align}
whereas the buoyant force due to the associated release of latent heat is
\begin{align}
\tilde{F}_{\Delta c^\prime}=-\beta_T \frac{h_v}{c_p}\Delta c^\prime g.\label{eq:force_temp_diff}
\end{align}
Given the parameters in this study, the leverage ratio between $\tilde{F}$ and $F$ is
\begin{align}
\frac{\tilde{F}}{F} = -\frac{\beta_T}{\beta_c}\frac{h_v}{c_p}=-22.99.\label{eq:buoyancy_leverage}
\end{align}
Since the condensation mass flux and the velocity field are connected only via the buoyancy term, this leverage ratio shows that the phase transition affects the flow primarily in terms of the associated release of latent heat, far more than by the change in vapor concentration itself.
Note that the negative sign of the leverage ratio signifies that the buoyant forces directly resulting from condensation are opposed to each other, while on the scale of the whole channel, a reduction in temperature leads to a reduction in vapor concentration such that the thermal and solutal buoyancy are aligned.

Two cases are considered in the following.
The primary case includes all fields and interactions as discussed above, in particular the modeling of phase transition (PT).
A reference case omits the vapor concentration field and does not take phase transitions into account (nPT).
By comparing both cases, the influence of condensation can be isolated.

The simulations are initialized using a perturbed laminar channel flow profile\cite{DeVilliers2006}, which is evolved in the absence of the scalar fields until turbulent flow has developed.
In the next phase, the temperature field is included in the simulation, and the system is simulated until the mean temperature within the thermal boundary layer converges.
At this point, the two branches are created.
In the first branch, calculations continue as before and temporal averaging is started to gather data for the reference case nPT, which is used as a baseline for comparison.
For the PT branch, the vapor concentration field and the phase transition model are added to the evolution, and temporal averaging begins after the mean concentration in the vapor boundary layer has converged.
Averages are calculated over intervals of $\Delta t^+=40$ for nPT and $\Delta t^+=36$ for PT.

Discretization with a resolution of $\Delta x^+=5.3$ and $\Delta z^+=2.7$, $\Delta y^+=0.2$ -- $3.3$ is applied consistently to both the primary and the precursor domain to avoid numerical disturbances at the interface.
The resulting mesh consists of $N_x \times N_y \times N_z = (480+320)\times 180 \times 316$ cells.

The upper limit for the discrete time step of the simulation $\Delta t$ is determined by the stability criterion for a second-order central differencing scheme in a system including diffusive scalar transport\cite{Shishkina2004}.
The resulting time step is $\Delta t^+ = 5\cdot 10^{-5}$.
For the combination of low vapor load and moderate subcooling present in the system considered in this study, the integration of the vapor source term derived from the Hertz-Knudsen-Schrage equation given in equation (\ref{eq:HKS_vapor_source}) results in the complete removal of excess vapor in every step of the simulation, since the time scale given by the condensation rate is much smaller than the simulation time step.
This behavior can be mimicked by setting the source term $f$ for equations (\ref{eq:cd_temperature}) and (\ref{eq:cd_concentration}) to
\begin{align}
f(T,c) = \begin{cases} \big( c-c_{\mathit{sat}}(T)\big)/\Delta t,\quad c>c_{\mathit{sat}}(T)\\ \phantom{(}0,\quad \text{else.}\end{cases}\label{eq:simplified_vapor_source}
\end{align}
In this simplified formulation, any vapor in excess of the local saturation value is removed instantaneously, yielding results indistinguishable from the original source term at reduced computational cost.
\section{Results}
\begin{figure}[htb]
\centering
\includegraphics[scale=1]{fig_source_pattern.png}
\caption{Normalized vapor concentration change due to condensation at the cooled wall. High values of $f$ correspond to hot spots of mass transfer across the phase boundary.}\label{fig:source_pattern}
\end{figure}
The DNS of the flow with condensation allows observation of the phase transition rates in the subcooled region within the thermal boundary layer at the cooled wall.
Figure \ref{fig:source_pattern} shows a heatmap of the normalized vapor source term
\begin{align}
\frac{f \Delta t}{\Delta c} = \frac{c-c_{\mathit{sat}}}{\Delta c},
\end{align}
evaluated on the cooled wall boundary patch.
A zone of high mass transfer is visible at $x=0$, where condensation occurs across the full width of the channel as the humid air first meets the cooled wall.
Further downstream, the pattern of condensation is organized into longitudinal, streak-like structures, reminiscent of the wall-layer streaks of fast and slow moving fluid, with a similar spacing of $\lambda_z^+\approx 100$ by visual estimation.\cite{Smith1983}

The inhomogeneity of the condensation pattern and the apparent footprint of coherent flow structures suggest that features of the turbulent flow directly govern the spatial distribution of phase transition events throughout the channel.
Fast-moving fluid originating from the bulk is transported towards the wall in sweeps.
Since it has not yet interacted with the cooled wall, the temperature and vapor concentration will be equal to the values prescribed at the channel inlet, small corrections due to diffusive transport notwithstanding.
At the wall, the fluid is slowed down and cooled by the interaction with the channel boundary, and the vapor concentration drops as a consequence of condensation.
Regions of high condensation rates therefore coincide with the impingement of sweeps onto the cooled surface.
Subsequently, the fluid is ejected from the wall and travels back towards the bulk, moving comparatively slower and carrying less vapor than before.

As a consequence of this vapor transport mechanism, fluid can be classified as pre- and post-interaction based on the instantaneous fluctuations of the scalar fields.
Analogous to equation (\ref{eq:def_theta}), a normalized concentration field is given by
\begin{align}
\xi = \frac{c-c_{\mathit{ref}}}{\Delta c}.\label{eq:def_xi}
\end{align}
The reference values for the given setup are $T_{\mathit{ref}}=T_{\mathit{in}}$ and $c_{\mathit{ref}}=c_{\mathit{in}}$.

Defining the thermal and solutal boundary layer as the set of locations $(x_b, y_b)$ where the cooled wall influences the average flow field due to turbulent mixing or diffusive transport means that $\langle\theta\rangle(x_b,y_b) < \theta_{\mathit{in}}$ and $\langle\xi\rangle(x_b,y_b) < \xi_{\mathit{in}}$ by construction.
Pre-interaction fluid will therefore cause positive fluctuations of the scalar fields, since it is characterized by $\theta_{\mathit{in}}$ and $\xi_{\mathit{in}}$ as outlined above.

Conversely, since this pre-interaction fluid contributes to the overall average, the value of the scalars in fluid that has already interacted with the cooled wall must be below this average value, and therefore cause negative fluctuations.
In combination, these observations yield conditions based on the instantaneous fluctuations of the scalar fields,
\begin{align}
\theta(x,y,z,t) - \langle \theta \rangle (x,y) \begin{cases} > \phantom{-}\epsilon_\theta \quad \text{warm fluid, } \theta_+,\\ < -\epsilon_\theta \quad \text{cold fluid, } \theta_- \end{cases}\label{eq:criterion-theta}
\end{align}
and
\begin{align}
\xi(x,y,z,t) - \langle \xi \rangle (x,y) \begin{cases} > \phantom{-}\epsilon_{\xi} \quad \text{humid fluid, } \xi_+,\\ < -\epsilon_{\xi} \quad \text{dry fluid, } \xi_-. \end{cases}\label{eq:criterion-c}
\end{align}
$\epsilon_{\theta,\,\xi}>0$ can be chosen to improve the signal-to-noise ratio of the criteria, the trade-off being loss of information about fluid closer to the local average than the threshold value.
Here, $\epsilon_{\theta,\,\xi} = 0.01$ are chosen as the limit for the classification.
Averages conditioned in such a way are denoted by $\langle \cdot \rvert \theta_{\pm}\rangle$ and $\langle \cdot \rvert \xi_{\pm}\rangle$, respectively.

Although averaging conditioned on the temperature fluctuations appears to be the natural choice for comparison---as they can be evaluated in both cases irrespective of the presence of water vapor and phase transition effects---it is mandatory to use the fluctuations of the concentration field in the PT case.
Temperature fluctuations cannot be used as the distinguishing criterion here since the injection of latent heat during condensation leads to a momentary increase in temperature, such that post-interaction fluid could possibly be classified as $\theta_+$ and sampled erroneously as pre-interaction.

The conditionally averaged velocity fields are then investigated to confirm the convective transport mechanism leading to the condensation pattern observed in Figure \ref{fig:source_pattern}.

Figure \ref{fig:conditional_ux_PT} shows profiles of the conditionally averaged streamwise velocity $\langle u_x \rvert \xi_{\pm}\rangle$ compared to the indiscriminately sampled average for PT.
The conditionally averaged values are only well-defined inside the boundary region, since the fluctuations of the scalar fields are zero outside the thermal and solutal boundary layers, resulting in truncated wall-normal velocity profiles.
\begin{figure}[hb]
\centering
\includegraphics[scale=1]{fig_conditional_ux_PT-eps-converted-to.pdf}
\caption{Total averaged streamwise velocity profiles {-}{-}{-} $\langle u_x \rangle$ compared to conditionally averaged profiles \textbf{\textcolor{humid}{---}} $\langle u_x \rvert \xi_+\rangle$ and \textbf{\textcolor{dry}{---}} $\langle u_x \rvert \xi_-\rangle$ at different positions along the channel length.}\label{fig:conditional_ux_PT}
\vspace*{2mm}
\centering
\includegraphics[scale=1]{fig_conditional_uy_PT-eps-converted-to.pdf}
\caption{Conditionally averaged profiles of the wall-normal velocity \textbf{\textcolor{humid}{---}} $\langle u_y \rvert \xi_+\rangle$ and \textbf{\textcolor{dry}{---}} $\langle u_y \rvert \xi_-\rangle$ at different positions along the channel length.}\label{fig:conditional_uy_PT}
\end{figure}

The ordering $\langle u_x \rvert \xi_-\rangle \leq \langle u_x \rangle \leq \langle u_x \rvert \xi_+\rangle$ can be observed at every position.
Additionally, the wall-normal component of the conditionally averaged velocity $\langle u_y \rvert \xi_{\pm}\rangle$ shows that $\xi_+$-tagged fluid moves towards the cooled wall, while $\xi_-$-tagged fluid exhibits the opposite trend, (Figure \ref{fig:conditional_uy_PT}).

This confirms the hypothesized transport mechanism.
Vapor is carried from the bulk towards the wall in fast-moving sweeps, where it is slowed down and the vapor concentration drops as a consequence of condensation.
After the interaction, the slowed and dried fluid is ejected from the wall and travels back towards the channel center line, causing the growth of the boundary layer as a function of the downstream position.

The evolution of the composition of the boundary layer can be observed in the relation between the sampled averages and the total average of the streamwise velocity.
At the most upstream position shown ($x/\delta=1\pi$), the overall average closely follows the profile of the pre-interaction fluid, suggesting that at this point, the majority of the fluid found in the boundary layer has not yet interacted with the wall.
In contrast, at the downstream positions ($x/\delta=5\pi$ and $6\pi$), a larger portion of the fluid has come into contact with the wall due to the longer residence time.
This is reflected by the overall average tending more towards the profile of $\langle u_x \rvert \xi_-\rangle$.

The analysis of the conditionally averaged velocity profiles illustrates the role of turbulence for the distribution of condensation at the cooled wall.
Convective transport of vapor from the bulk to the wall in turbulent sweeps is the determining factor of phase transition localization, leading to a pattern of condensation that acts as a footprint of the underlying turbulence, as shown in Figure \ref{fig:source_pattern}.

\begin{figure}[t]
\centering
\includegraphics[scale=1]{fig_isotherms_PT_nPT-eps-converted-to.pdf}
\caption{$\langle \theta \rangle = 0.75$ and $\langle \theta \rangle = 0.99$ isotherms compared between cases PT and nPT.}\label{fig:isotherms_PT_nPT}
\end{figure}
\begin{figure}[t]
\centering
\includegraphics[scale=1]{fig_ux_rms_PT_nPT-eps-converted-to.pdf}
\caption{Wall-normal profile of the fluctuations of the streamwise velocity component at the channel outlet compared between cases PT and nPT.}\label{fig:ux_rms_PT_nPT}
\end{figure}
The complete coupling implemented for the governing equations causes the influence of the turbulent flow on the occurrence of condensation to feed back to the overall flow.
The injection of heat into the wall-near regions during the phase transition directly affects the thermal boundary layer.
Figure \ref{fig:isotherms_PT_nPT} shows the growth of the thermal boundary layer as characterized by the position of the $\langle \theta \rangle = 0.99$ isotherm along the channel length, compared between cases PT and nPT.
Additionally, the modification of the temperature field in the wall-near region is illustrated by the $\langle \theta \rangle = 0.75$ isotherm.
The influence of the release of latent heat can be observed here.
Fluid with $\langle \theta \rangle<0.75$ extends further into the channel in the case without condensation, indicating that the direct cooling effect of the wall is partially counteracted by the phase transition.
However, the $\langle \theta \rangle = 0.99$ isotherm exhibits the opposite trend, with a thicker overall thermal boundary layer for PT.
This phenomenon is a consequence of the interaction of the localized temperature modification at the wall with the overall turbulent flow field via the buoyant force terms.
The increased average fluid temperature interferes with the damping effect of the aiding buoyant forces on the turbulence.
Figure \ref{fig:ux_rms_PT_nPT} shows the profile of the fluctuations of the streamwise velocity $u_{x\mathit{, rms}}^+$ at the channel outlet.
At this position along the channel, acceleration and deceleration effects of the buoyant forces are more pronounced than at upstream locations, since the residence time of the fluid in the differentially cooled system is maximized.
The well-established damping effect of the aiding buoyant forces acting at the cooled wall\cite{Kasagi1997, Wetzel2019} at $y/\delta=-1$ is visible in the asymmetric profile of the fluctuations, with a reduced peak height near the cooled wall.
Crucially, this damping effect is slightly reduced in the case with condensation, resulting in higher mixing rates and increased transport of cooled fluid towards the bulk, causing the thicker thermal boundary layer observed in Figure \ref{fig:isotherms_PT_nPT}.
This behavior results directly from the negative leverage between solutal and thermal buoyancy during condensation events given in equation (\ref{eq:buoyancy_leverage}).
While the loss in vapor concentration adds to the aiding force acting on the fluid and would therefore increase the damping of turbulence, the leveraged buoyant force caused by the associated release of latent heat results in the opposite behavior, amplifying the turbulence instead.

These modifications of the overall flow field observed in PT together with the connection between the spatial distribution of condensation events and the structure within the turbulent flow suggests that these structures themselves should be sensitive to phase transition in the flow.
To investigate this interaction, spectra of the streamwise velocity fluctuations conditioned on the fluctuations of the scalar fields,
\begin{align}
\langle u_x^\prime \rvert \xi_{\pm} \rangle = \left\langle\left. \sqrt{(u_x - \langle u_x\rangle )^2}\; \right\rvert\xi_{\pm} \right\rangle
\end{align}
and
\begin{align}
\langle u_x^\prime \rvert \theta_{\pm} \rangle = \left\langle\left. \sqrt{(u_x - \langle u_x\rangle )^2}\; \right\rvert\theta_{\pm} \right\rangle,
\end{align}
are extracted for PT and nPT, respectively.
Using the spanwise periodicity of the channel geometry, the energy spectra with respect to the spanwise wavenumber are calculated using fast Fourier transformations.
Since the global effect of buoyancy is small due to the short residence time of the fluid within the system and the limited vapor load and temperature difference (as expressed by the low Richardson number), the streamwise change in the mean velocity is small compared to the magnitude of the fluctuations, such that additional averaging can be performed along this axis.
For better comparability between the different cases, the spectra are normalized with the total energy integrated across all wavenumbers,
\begin{align}
E^\ast(\mathbf{k}) = E(\mathbf{k})\left( \int E(\mathbf{k})\,\mathrm{d}\mathbf{k}\right)^{-1}.
\end{align}

The nondimensional normalized pre-multiplied spectra of the streamwise velocity fluctuations $(k_z \delta)E^\ast_{xx}$ at a distance of $y^+=15$ from the cooled wall, sampled conditionally in accordance with the criteria to distinguish between pre- and post-interaction fluid (\ref{eq:criterion-theta}) and (\ref{eq:criterion-c}), are evaluated in the following.
Figure \ref{fig:kzExx_pre_PT_nPT} shows the spectra of pre-interaction fluid for PT and nPT.
Apart from the negligible effect on the global Grashof number, the vapor concentration field does not directly affect the flow field, only acting indirectly via the leveraged buoyancy during the phase transition.
Consequently, no difference between the two cases can be observed when sampling specifically for fluid that has not yet interacted with the subcooled region near the cooled wall.
In turn, this result underlines that the sampling criteria defined above are well-suited to distinguish between pre- and post-interaction fluid.

\begin{figure}[t]
\centering
\includegraphics[scale=1]{fig_kzExx_pre_PT_nPT-eps-converted-to.pdf}
\caption{Normalized pre-multiplied energy spectrum sampled from pre-interaction fluid \textbf{\textcolor{humid}{---}} $\xi_+$ and \textbf{\textcolor{warm}{---}} $\theta_+$ for PT and nPT.}\label{fig:kzExx_pre_PT_nPT}
\end{figure}
\begin{figure}[t]
\centering
\includegraphics[scale=1]{fig_kzExx_post_PT_nPT-eps-converted-to.pdf}
\caption{Normalized pre-multiplied energy spectrum sampled from post-interaction fluid \textbf{\textcolor{dry}{---}} $\xi_-$ and \textbf{\textcolor{cold}{---}} $\theta_-$ for PT and nPT.}\label{fig:kzExx_post_PT_nPT}
\end{figure}
\begin{figure}[htb]
\centering
\includegraphics[scale=1]{fig_correlations_ux_pre_post-eps-converted-to.pdf}
\caption{Correlation coefficient for the streamwise velocity fluctuations as a function of the spanwise separation $\Delta z^+$ compared between both cases. Only one curve is shown for pre-interaction fluid as the differences between PT and nPT are negligible.}\label{fig:uxux_pre_post_PT_nPZ}
\end{figure}

Figure \ref{fig:kzExx_post_PT_nPT} compares the same quantities after the interaction with the cooled wall.
In both PT and nPT, the spectrum of the fluid ejected from the wall is shifted to higher wavenumbers compared to the sweeps carrying pre-interaction fluid.
This shift is more pronounced in the case without phase transition with regards to both peak position and the high-wavenumber tail of the spectrum.
As with the averaged fluctuations of the streamwise velocity (Figure \ref{fig:ux_rms_PT_nPT}), the release of latent heat at or near the cooled wall during condensation inhibits the cooling of the fluid.
In turn, the turbulent flow structures are affected by the leveraged buoyant forces feeding back from the temperature to the velocity field.
Consequently, the shift towards larger wavenumbers caused by the interaction with the cooled wall without the vapor concentration field is opposed by the effects of the phase transition.
In terms of streak spacing, the difference between PT and nPT is clearly visible in the spanwise spatial correlation coefficient of the velocity fluctuations $\langle u_x^\prime u_x^\prime\rangle$, shown in Figure \ref{fig:uxux_pre_post_PT_nPZ}.
The correlation length $r$, defined as the separation $\Delta z$ from the maximum to the minimum of the correlation coefficient, is increased from $r_{\text{nPT}}^+=24.2$ to $r_{\text{PT}}^+=29.5$ for post-interaction fluid, a change of $\Delta r^+=5.3$.

These results show the impact of the mutual interplay between the release of latent heat during the phase transition and the turbulent flow.
Turbulent flow structures govern the transport of the condensable phase towards the cooled wall and determine the arrangement of the condensation sites.
In turn, the injection of latent heat into the fluid associated with the phase transition changes the turbulent flow structure along the channel by opposing the cooling influence of the channel wall.

Based on this knowledge about the leverage mechanism coupling condensation and the flow fields, the analysis of the behavior for flows with a different ratio of forced to natural convection reduces to the discussion of the effects of thermal buoyancy in such flows\cite{Wetzel2019}.
As $\mathrm{Ri}\to 0$, buoyant forces become negligible compared to the inertia of the fluid, and the effects of the phase transition explored here vanish.
In turn, for flows with a large contribution from natural convection, the opposing effect of the release of latent heat at the cooled wall will grow more important as well.
\section{Conclusions}
Direct numerical simulations of turbulent convective channel flow including contributions from thermal and solutal buoyancy were performed.
The aspect of latent heat release during phase transitions was included in the simulations by modeling condensation via appropriately chosen source terms for the concentration and temperature fields.
These fields were treated as active scalars, transported through the system by convection and diffusion.
This approach allowed the investigation of this aspect of phase transitions without a full multiphase simulation.

The setup of the system with an inlet-outlet simulation domain allowed for the analysis of the evolution of the flow as a function of the streamwise coordinate, showing that the residence time of the fluid within the region of influence of the cooled wall determines the degree to which phase transition effects affect the fluid on average.
Additionally, this configuration lends itself more readily to comparisons with experimental investigations of similar systems in the future.

Simulations of an identical system without the inclusion of the vapor field and phase transition modeling provided a baseline for comparisons against which the effect of the condensation on the flow could be evaluated.

The analysis of the mutual interplay between the turbulent flow and condensation revealed a clear bidirectional link between both aspects of the system.
The arrangement of hot spots of mass transfer across the phase boundary at the cooled wall reflects the structure of the alternating longitudinal streaks of high- and low-velocity fluid in the turbulent flow.
This condensation pattern results from the mechanism of convective vapor transport.
Fast-moving sweeps carry fluid with high vapor concentration towards the wall, thus creating oversaturated conditions and causing condensation events.

In turn, the turbulent flow itself is modified by the effects of the phase transition due to the coupling via the buoyant forces.
Here, the influence of buoyancy induced by the release of latent heat is opposed and far greater than the solutal contribution from the underlying change in vapor concentration during condensation.

Collectively, the modifications of the turbulent flow are characterized by the opposition to the cooling effect caused by the injection of heat into the subcooled region near the wall.
While the system is very efficient in removing the added thermal energy, such that an increase in average temperature is found only directly at the wall, the high spatial and temporal correlation between the associated generation of temperature spikes during condensation events and the underlying turbulent flow structure reduces the damping effect of the aiding buoyant force acting on the fluid.
This results in a slightly elevated peak turbulence intensity compared to the case without phase transition.
The small effect size is a consequence of the small overall importance of buoyancy compared to the inertia of the forced convection for the flow parameters investigated here.

Conditional sampling based on the instantaneous fluctuations of the scalar fields proves to be a valuable tool for distinguishing between pre- and post-interaction fluid.
Applied to the pre-multiplied spectra of the streamwise velocity fluctuations, they again show that the injection of latent heat during condensation serves to resist the changes exhibited by the system without phase transition, as the shift towards higher wavenumbers is considerably reduced.
This reflects a significant modification of the spanwise spacing of post-interaction streaks.

The treatment of condensation is limited to the effect of the release of latent heat.
In particular, the drop of the partial pressure of the vapor due to the phase transition\cite{Bukhvostova2014}, inducing compressibility effects, is not considered here.
Additionally, the deposition of condensate at the walls and the subsequent formation of droplets provides another mechanism for interaction with the turbulent flow.

Including these aspects of condensation is possible while at the same time preserving the single-phase approach, promising a simplified, yet accurate method of simulating flow with phase transition for systems similar to those presented in this investigation.
\section*{Acknowledgments}
The authors thank Annika K\"ohne for proofreading the paper.
\bibliographystyle{aipnum4-1}


\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.