  The design of accurate helium-solute interaction potentials for the
simulation of chemically complex molecules solvated in superfluid helium has
long been a cumbersome task due to the rather weak but strongly anisotropic
nature of the interactions. We show that this challenge can be met by using a
combination of an effective pair potential for the He-He interactions and a
flexible high-dimensional neural network potential (NNP) for describing the
complex interaction between helium and the solute in a pairwise additive
manner. This approach yields an excellent agreement with a mean absolute
deviation as small as 0.04 kJ/mol for the interaction energy between helium and
both, hydronium and Zundel cations compared to CCSD(T) reference calculations
with an energetically converged basis set. The construction and improvement of
the potential can be performed in a highly automated way, which opens the door
for applications to a variety of reactive molecules to study the effect of
solvation on the solute as well as the solute-induced structuring of the
solvent. Furthermore, we show that this NNP approach yields very convincing
agreement with the CCSD(T) reference for properties like many-body spatial and
radial distribution functions. This holds for the microsolvation of the
protonated water monomer and dimer by a few helium atoms up to their solvation
in bulk helium as obtained from path integral simulations at about 1 K.
