Not so loosely bound rare gas atoms: finite-temperature vibrational fingerprints of neutral gold-cluster complexes

We present an experimental and theoretical study of the structure of small, neutral gold clusters—Au3, Au4 and Au7—‘tagged’ by krypton atoms. Infrared (IR) spectra of AuN·KrM complexes formed at 100 K are obtained via far-IR multiple photon dissociation in a molecular beam. The theoretical study is based on a statistical (canonical) sampling of the AuN·KrM complexes through ab initio molecular dynamics using density-functional theory in the generalized gradient approximation, explicitly corrected for long-range van-der-Waals (vdW) interactions. The choice of the functional is validated against higher-level first-principle methods. Thereby finite-temperature theoretical vibrational spectra are obtained that are compared with the experimental spectra. This enables us to identify which structures are present in the experimental molecular beam for a given cluster size. For Au2, Au3 and Au4, the predicted vibrational spectra of the Kr-complexed and pristine species differ. For Au7, the presence of Kr influences the vibrational spectra only marginally. This behavior is explained in terms of the formation of a weak chemical bond between Kr and small gold clusters that localizes the Kr atom at a defined adsorption site, whereas for bigger clusters the vdW interactions prevail and the Kr adatom is delocalized and orbits the gold cluster. In all cases, at temperatures as low as T = 100 K, vibrational spectra already display a notable anharmonicity and show, in comparison with harmonic spectra, different position of the peaks, different intensities and broadenings, and even the appearance of new peaks.


Introduction
As a bulk elemental metal, gold is a classic example of inertness [1]. However, at the nanoscale gold exhibits surprising chemical activity [2,3]. In fact, this property of nano-gold has already been used for commercial applications in offensive odor removal and gas sensors [4]. While most of the previous research efforts have been spent on deposited clusters, the experimental study of gas-phase clusters has the advantage of the reliable knowledge of the clusters' size (number of atoms) and charge. Such well-defined species represent an ideal situation for modeling and analyzing their properties by means of experiments and ab initio calculations. In the past years, the structure of small gold clusters has been studied in a series of works at various levels of theory [5][6][7][8][9][10]. Combined theoretical and experimental works on anions [11][12][13][14][15] and cations [14,16] yielded a consistent picture of the charged-cluster structures at all small sizes and, in particular, were able to identify the size at which three-dimensional (3D) structures become more stable than two-dimensional (2D) isomers. While theory has been equally applied to neutral gold clusters, their experimental characterization is more scarce and limited to the determination of ionization potentials (IPs) via electron impact (N = 1-22) [17] and optical absorption spectra [18].
Au clusters toward Kr. In section 4, we analyze Au 3 · Kr, Au 3 · Kr 2 and Au 4 · Kr 2 , where the unusually strong interaction between the clusters and Kr atom(s) is further detailed. In section 5, we finally re-examine the spectrum of Au 7 · Kr and show how the theoretical finite-temperature vibrational spectrum can explain finer details than the harmonic analysis previously performed in [19].

Experimental approach
The FIR-MPD experimental setup has been described elsewhere [19,27]; here we summarize the key aspects. Gold clusters are obtained by laser vaporization from a gold rod and then thermalized in a flow of He and Kr kept at T = 100 K. Thereby Kr atoms (one or two, rarely more) may adsorb on the pristine gold clusters. These Au N · Kr M complexes are subsequently investigated in a molecular beam, that is irradiated by a pulsed FIR beam from the Free Electron Laser for IR eXperiments (FELIX) [28]. The laser is tunable throughout the mid and far-IR (40-2300 cm −1 ). Subsequently, the neutral complexes are ionized by means of an F 2 -laser with an energy of 7.9 eV per photon and the ionized species are mass-analyzed in a time-of-flight mass spectrometer. When the FIR radiation is in resonance with an IR-active mode of a neutral complex, photons can be absorbed, the complex is heated and evaporation of the Kr ligand may follow. In this way, a depletion of the mass spectrometric signal of the gold-krypton complex results. Analyzing the frequency dependence of the depletion signal leads to the cluster-size specific IR spectra. Experimental IR intensities reported here are normalized for photon fluence rather than the laser intensity used in our previous studies. We have recently found that this gives better agreement with theoretical IR intensities if a wide spectral range is covered [29]. The observed vibrational bands of gold clusters are entirely in the FIR, namely at frequencies below 200 cm −1 . This is due to the large mass of the nuclei as well as the softness of the bonds.

Ab initio potential energies and forces, and their validation.
The theoretical results presented in this work were obtained using the FHI-aims [30] program package for an accurate all-electron description based on numeric atom-centered basis functions. Where not differently specified, for our analysis we employed (collinear) spin polarized DFT at the PBE [31] generalized gradient approximation (GGA) level, corrected for long range vdW interactions via the Tkatchenko-Scheffler (TS) scheme, i.e. a sum over C 6 [n]/R 6 tails, and C 6 coefficients derived from the self-consistent electron density n and reference values for the free atoms [32]. This functional will be referred throughout the paper as PBE + vdW. We used 'tight' integration grid settings and accurate 'tier 2' basis sets [30]. The scaled zeroth-order regular approximation to the Dirac equation (ZORA) (scalar) relativistic correction [33] was employed for the static calculations. However, the computational cost of evaluating forces with the latter method is prohibitive; for the MD runs, where forces need to be evaluate at each timestep, the 'atomic ZORA' scalar relativistic correction [30] was used. This scheme has been shown to provide remarkably good binding energy and bond distance for Au 2 in comparison to benchmarks methods [30]. We performed a test on the larger clusters analyzed here and we find that the 'atomic ZORA' and scaled ZORA yield binding energies that differ less than 0.02 eV atom −1 for all functionals. When the difference in binding energies between isomers, of the same size N of Au N , is examined, the two methods agree within 0.01 eV atom −1 . Harmonic vibrational frequencies and intensities were computed from finite differences of the analytic forces. The binding energy of Kr to the gold dimer, equilibrium geometry, static electric dipole moment, and harmonic spectrum of Au 2 · Kr, as calculated with PBE + vdW, were compared to a hierarchy of electronic structure methods, from the local-density approximation (LDA) functional, through the hybrid functional PBE0 and the double-hybrid XYG3 [34], to RPA + rSE and rPT2 applied on both PBE and PBE0 orbitals [35]; furthermore, also MP2 and CCSD(T) values were calculated. All methods beyond GGA, except CCSD(T), were calculated with FHIaims and we used really tight settings and tier 4 basis set. CCSD(T) values are calculated with Gaussian03 (revision D.01) [36] and aug-cc-pVTZ-PP basis set of Figgen et al [37] and Peterson and co-workers [38,39]. For some of the vdW complexes (see below), the interactions between the cluster and the RG atom were calculated at the MP2 level including the recently introduced correction to the dispersion interactions, MP2+ vdW [26]; MP2+ vdW energies have been shown to be in excellent agreement with CCSD(T) calculations for systems bonded by dispersion forces. A very good agreement between PBE + vdW and MP2 + vdW is found, which confirms the reliability of our calculated results reported and discussed below.
As explained in section 'Experimental approach', we can only detect a particular species and record its IR spectrum if its IP is lower than the energy of the F 2 -laser used to ionize the species in the molecular beam. For this reason, we have evaluated the vertical IP (vIP) of the species here analyzed, in two different ways (vide infra in table 2). (i) Energy difference between the (relaxed) neutral cluster and the cationic cluster (in the neutral cluster geometry), both evaluated at the PBE + vdW level. (ii) By evaluating the one-shot perturbative single-particle excitation (G 0 W 0 [40]), starting from PBE orbitals of the neutral cluster. The latter evaluation was performed with FHI-aims, with 'really-tight' settings and 'tier 4' basis set. For our systems, the difference of predicted vIPs between the two methods is within 0.3 eV (vide infra in table 2).

Statistical mechanics.
The IR spectra beyond the harmonic approximation of the clusters were calculated by performing Born-Oppenheimer MD simulations in the canonical ensemble at the experimental temperature (see next session for the definition of temperature) and extracting the Fourier transform of the dipole-dipole autocorrelation function from the trajectories. Thus, the IR intensities are computed via where M(t) is the total electric dipole of the cluster at time t, β = 1/k B T and the angular brackets indicate an ensemble average in the canonical thermodynamical ensemble. We assume the system as ergodic: this means that a time average performed on a long thermostatted trajectory is equivalent to an ensemble average in the N V T ensemble. A trajectory is judged 'long enough' when the vibrational spectrum calculated for the whole trajectory does not change any more. The scalar product in the integral is averaged by selecting several times t = 0 along the same trajectory. The interval between two subsequent t = 0 is chosen to be longer than the time for the decay of the dipole-dipole autocorrelation function from one to the long-time average. This is because at short times the correlation between dipole moments (as for any other property of the system) is nearly 1 (the vector has still a similar modulus and direction). At the time at which the scalar product reaches the long-time average, the memory of the initial time is lost and thus a new dipole can be used as initial one for the statistical average. The factor βω 2 in front of the integral is the result of the product of the classical pre-factor βω(1 − exp(−βhω)) and the quantum correction factor ω/(1 − exp(−βhω)) [41,42]. The classical factor results from the assumption of Boltzmann statistics for the ensemble of oscillators, while the quantum factor corrects for the so-called detailed balance, which reflects in an asymmetry of the peaks in the spectrum. The quantum correction is not uniquely defined, but the one we applied was shown to be the most accurate [43,44] when comparing theoretical and experimental spectra.
In the literature, autocorrelation functions are normally calculated from simulations in the microcanonical ensemble (N V E, i.e. constant number of particles N , constant volume V , constant energy E) and then referred to the average temperature of the run (such simulations are typically pre-equilibrated with a thermostat in order to impose the target temperature). In our case, though, the small number of DoF required a thermostat during the sampling of the correlation function. The reason for this is that, when the DoF are few, the distribution of the kinetic energy in a NVE ensemble departs from the distribution of the canonical ensemble (constant number of particles, constant volume, constant temperature, NVT) at the same average temperature (the latter distribution is nothing else than the Boltzmann distribution). The NVT distribution has a thick tail at large kinetic energies [45], while the NVE distribution is a Gaussian function around the average temperature (both distribution have 2/(3N ) relative variance). When the number of DoF is large (rigorously, at the thermodynamic limit), the two distributions converge to the same shape and a simpler NVE simulation, after thermalizing the system at the desired temperature, would be a good approximation of the rigorous NVT sampling.
There are two ways to overcome this problem, either averaging the correlation function over an ensemble of NVE trajectories, where the initial states (coordinates and velocities) are extracted from a canonically distributed set at the target temperature, or using a thermostat that does not perturb the dynamics. The first solution is computationally very demanding and the second requires a not trivial implementation. Since a thermostat always acts on the velocities, it is difficult to design one that does not destroy the dynamical correlations. Recently, Bussi et al introduced a stochastic thermostat that fulfills this requirement [45]. We tested the thermostat by calculating spectra via equation (1) at very low temperatures. The results reproduced the harmonic spectra impressively well. Furthermore, we observe that the finitetemperature spectrum is practically independent of the only tuning parameter that the thermostat has, which can be interpreted as a relaxation time, over a wide range of its values.
In the figures where we compare theoretical and experimental spectra (figures 3-6) we have shifted the theoretical spectra in order to align the frequencies of the peaks to the frequencies at which the experimental peaks occur (note that a rigid shift is sufficient to align all the peaks). The necessity of such a shift is due to force inaccuracies caused by the approximate exchangecorrelation functional and the finiteness of the basis set, but also they are caused by the finiteness of the MD timestep and the granularity of the mesh onto which the basis functions are projected. The sensitivity of the theoretical spectra toward the above mentioned settings is debated in the appendix. All the theoretical spectra shown in this paper are obtained with the same setting and we found that a rigid frequency (blue-)shift of 8 cm −1 optimized the matching for all cases. In figures 3-6, we also report the harmonic spectra, where the frequencies of the peaks (bars) were scaled by a factor 1.05. In this way the peaks are approximately aligned to the experimental and finite-T theoretical spectra, in order to help the visual comparison. For converged vibrational spectra, MD runs of at least 100 ps were needed, and we used a timestep of 10 fs. A stable integration of the equations of motion with such an unusually large timestep is allowed by the low value of the highest frequency phonon (∼200 cm −1 ) in our systems.

Definition of temperature: classical versus quantum statistics for nuclei
In Born-Oppenheimer MD the nuclei are propagated as classical (point) particles. As a consequence, the population of their vibrational modes in the canonical ensemble obeys classical (Boltzmann) statistics. However, nuclei are quantum particles and also the population of the vibrational modes is quantized. Strictly speaking, one cannot define a joint {positions, momenta} phase-space probability distribution for a quantum system. For the simple case of a harmonic oscillator, however, the semiclassical Wigner distribution associated with thermal populations of the vibrational states is a Gaussian distribution of position and momentum, with fluctuations that depend parametrically on frequency and temperature [46]. For each normal mode and for a given 'classical temperature' T , one could then define an effective 'quantum temperature' as the temperatureT at which the quantum oscillator would exhibit the same fluctuations as a classical oscillator of the same frequency and at temperature T . For a set of non-interacting harmonic oscillators, by equating the widths of the classical and quantum harmonic oscillators momenta distributions, one obtains [47] where ν i are the (harmonic) vibrational frequencies of the cluster under consideration, T is the classical temperature (the one by the thermostat adopted for our simulations), N is the number of DoF andT is the quantum temperature, which we adopt as an estimate of the 'real' temperature for the equilibrated system. We note that the two temperatures converge for large T . The discrepancy between the two temperatures can be intuitively understood in terms of zero point energy: a classical system has to use some temperature in order to give kinetic energy to vibrational modes, while for a quantum system these modes are already active at T = 0 K. Thus the classical temperature has to be higher in order to give the same kinetic energy to the vibrational modes. Interestingly, this mapping also defines a lowest classical temperature, which is, by taking the limitT → 0: N k B T = i hν i 2 . This is the (classical) temperature needed to activate all the zero point vibrations. In the rest of the paper, for each cluster we will give both the classical temperature T at which the thermostat was set and the estimated quantum temperatureT .
The thermalization of clusters in sources similar as used here to prepare the Au N · Kr M complexes has been characterized before and allows the conclusion that under our conditions equilibration to the source temperature is achieved [48,49]. Nevertheless, it has to be noted that the experimental FIR-MPD spectra may not come from an exactly canonically distributed population. This is due to the fact that in the molecular beam cluster complexes belonging to the hotter tail of the canonical distribution may spontaneously dissociate and thus not contribute to the depletion spectrum. The experiment would then be sampling only the colder part of the full distribution.

Localized bonding of Kr: Au 2 · Kr and Au 2 · Kr 2
In [50,51], the cationic gold atom is found to form a strong bond, suggested to be covalent on the basis of orbitals-population analysis, with the heavier rare gases. The CCSD(T) binding energies between Au+ and Ar, Kr and Xe are 0.29, 0.51 and 0.91 eV, respectively [50]. While the neutral gold atom only forms a weakly bonded vdW dimer with Kr (as well as Ar and Xe) [52], we find that Au 2 · Ar, Au 2 · Kr, and Au 2 · Xe are linear molecules where the RG-Au 2 interaction at equilibrium is unexpectedly strong. The PBE + vdW bonding energies of the RG with the gold dimer are 0.11, 0.22 and 0.43 eV for Ar, Kr and Xe, respectively. The values for Ne and He are 0.02 and 0.01 eV, i.e. there is practically only a vdW interaction between the two lighter RGs and the gold dimer (the equilibrium geometry is in these cases an isosceles triangle with Au 2 as the short basis). We find a similar trend for the other coinage metals, Cu and Ag. Consistently with the behavior of the surfaces of these coinage metals, Ag proves to be overall less binding, with a maximum of 0.13 eV for Ag 2 Xe, while Cu 2 has the interaction energies roughly halved when compared to the corresponding Au 2 · RG molecule. The detailed analysis of this unusual bonding between the dimer and rare gases will be presented elsewhere [53]. We note in passing that when Au is treated non-relativistically 9 , Kr would exhibit a negligible bonding, namely purely vdW. Here we focus on the vibrational properties of the Au N · Kr and Au N · Kr 2 complexes. Recently, in [54] a study was presented of the equilibrium distance and binding energy of Xe, Kr and Rn to Au, Ag and Cu small clusters, calculated at the CAM-B3LYP level 10 . Our equilibrium geometries and energies for Kr adsorbed on Au 2 , Au 3 and Au 4 qualitatively agree with the results there presented (cf figure 1). In particular, Kr was found to adsorb at distances between 2.7 and 2.9 Å to one Au atom in Au 2 , Au 3 and Au 4 with binding energies of 0.1-0.2 eV and a Au-Kr stretching frequency around 70 cm −1 . However, on one issue our results disagree: while in [54] a very small increase of the Au-Au distance upon adsorption of RG was found, we find a small decrease (see table 1). This is of no relevance for this paper. However, we note in passing that we carefully tested our results and the physical mechanism of the decrease of the Au-Au distance is in fact interesting. It will be discussed in a separate paper [53]. 9 Also Kr is consistently treated non-relativistically, but this is less crucial. When compared to relativistically corrected Au 2 , the PBE + vdW NR gold dimer has a bond distance significantly increased from 2.51 to 2.77 Å and a bonding energy changing from −2.38 to −1.57 eV. 10 Validation at M06-2X, MP2 and CCSD(T) level was performed only for dimers bound to Xe; MP2 binding energies were calculated for all studied system, but on CAM-B3LYP geometries.  [36] and aug-cc-pVTZ-PP basis set of Figgen et al [37] and Peterson and co-workers [38,39]. The binding energy of the gold dimer is E b (Au 2 ) = E(Au 2 ) − 2E(Au). The adsorption energy of the Kr atom(s) onto the Au 2 is is the total energy of the relaxed system.
The occurrence of an interaction between RG atoms and metal clusters, so strong that the vibrational spectrum of the pristine cluster is perturbed, was also observed by Gehrke et al [22] for charged Co clusters and Ar. In that case, the interaction was explained in terms of electrostatic interactions between the static charge at the metal cluster and the induced dipole at the RG. In our case, for neutral systems, this electrostatic explanation cannot be invoked.
The static dipole of the Au N · Kr molecule is non-zero (see table 1) and vibrations become IR active, with two marked lines in the harmonic spectra: the Au-Au stretch at 185 cm −1 and a Kr-Au 2 stretch at 82 cm −1 . Note that, while the higher frequency eigenmode is still recognizable as an Au-Au stretch with just a small blue-shift (cf figure 2), and Kr just makes it IR active by breaking the symmetry of the molecule, the lower frequency line is entirely due to the presence of Kr.
The bonding predicted at the PBE + vdW level is confirmed at higher level of calculation, as shown in  The harmonic frequencies are labeled corresponding to the eigenmodes of Au 2 · Kr 2 as illustrated on the right. ν 4,5 and ν 6,7 are doubly degenerate modes, respectively. Thin lines show the main character of the bands observed in the finite-temperature spectra, for which the thermostat was set to T = 23 K (T = 0 K) and T = 103 K (T = 100 K). These spectra were neither shifted nor scaled.
level calculations exhibit similar Au-Kr bond distances, static dipole moments, harmonic vibrational frequencies and harmonic IR intensities suggests that the nature of the bond is also the same in the different approaches.
In the case of Au 2 · Kr 2 , the vibrational properties of the dimer are also strongly modified by the adsorption of Kr. Like Au 2 , the Au 2 · Kr 2 molecule is linear and inversion symmetric. Thus, fewer modes than for Au 2 · Kr are IR active. Au 2 · Kr 2 clearly shows anharmonic features already at relatively low temperature (see figure 2). At T = 23.5 K (T = 0 K) the IR spectrum simulated by MD mimics closely the harmonic spectrum (this is an indication of the reliability of the MD settings, and in particular of the thermostat). At T = 103 K (T = 100 K), however, clear differences become apparent: the band related to the antisymmetric Kr versus Au 2 stretching is red-shifted to 57 cm −1 as compared to the harmonic frequency of 69 cm −1 and a new satellite peak appears at 48 cm −1 . This is due to the interaction between the mentioned antisymmetric stretching mode with the symmetric one (harmonic frequency 62 cm −1 ), which is IR inactive in the harmonic approximation. During the MD simulation of both Au 2 · Kr and Au 2 · Kr 2 , the Kr atoms were found to stay localized at their bonding site. In the following, we will label this adsorption sites as 'chemisorption sites', in order to distinguish them from the pure vdW adsorption sites (vide infra). The actual nature of these adsorptions, which is largely covalent but also involves a complex charge polarization and redistribution, is discussed in [53].
Both Au 2 · Kr and Au 2 · Kr 2 , however, are not detected in the FIR-MPD experiment because their IPs are significantly higher than the energy of the ionizing UV laser. The theoretical values of their vIP are 9.0 eV for Au 2 · Kr and 8.7 eV for Au 2 · Kr 2 (see table 2), i.e. well above the photon energy of the UV laser (7.9 eV) used in the experiment.  figure 1. E b (equation (3)) is the total binding energy, E b (equation (4)) is the interaction energy between the Kr atom(s) and the relaxed gold cluster, E vdW (equation (5)) is the vdW part of the interaction between the Kr atom(s) and the gold cluster. For bare clusters, this number (reported between brackets) is the total intra-cluster vdW interaction. For Au 3 , (a) labels the acute-angled and (o) the obtuse-angled isomer. For Au 4 , (rh) means rhombus and (Y) Y-shaped isomer. For Au 7 , (1) is the second isomer from the top in figure 6 while (2) is the upper one in the same figure. vIPs are evaluated as energy difference of two single-point calculations with PBE + vdW for the neutral and cationic cluster (column marked with SCF) and via G 0 W 0 [40], on PBE orbitals. Experimental values of the IPs for Au 2 , Au 3 , Au 4 and Au 7 are 9.5, 7.5, 8.6 and 7.8 eV, respectively [17].
The 'adsorption' energy of the Kr atom(s) onto the cluster is defined as where E(Au N ) is the total energy of the relaxed Au N cluster, Au N · Kr M is the relaxed adsorbate system and E(Kr) the total energy of a single Kr atom. Furthermore, we report the vdW interaction energy between the Kr atom(s) and the gold clusters, calculated as the overall vdW correction minus the vdW correction within the bare cluster: Due to the fact that some vdW interaction is present among the atoms of the gold cluster (and the larger the cluster, the larger is the intra-cluster vdW interaction), with the latter definition we single out the part of vdW interaction between the gold cluster and the adsorbed Kr atom(s). When this value is small in comparison to E b , the interaction between the Kr atom(s) and the gold cluster has some covalent character [53].

Au 3 · Kr and Au 3 · Kr 2
Au 3 has two (meta)stable isomers, (a) an obtuse-angled isosceles triangle 11 , with an obtuse angle of about 140 • and (b) an almost equilateral triangle 12 . According to PBE + vdW, the latter is 0.12 eV less stable than the former. The linear isomer is only a saddle point for the neutral Au 3 . PBE0 + vdW also finds the obtuse-angled isomer more stable, but only by 0.04 eV (after relaxing both structures with PBE0 + vdW). For higher level methods, however, the most stable isomer is the acute-angled. Namely, for XYG3 by 0.11 eV and for RPA + rSE@PBE by 0.05 eV. Also in [10], GGA functionals predict the obtuse-angled isomer to be more stable and higher level functionals the acute-angled. However, all functionals underestimate the formation energy of the gold trimer when compared to the experimental value (3.80 ± 0.13 eV) [57] and the best agreements comes from GGA functionals. Thus, the accurate relative energetics between pristine isomers has to be regarded as still an open issue. According to PBE + vdW, the binding of one or two Kr atoms to these two isomers of Au 3 brings the two Au 3 · Kr structures to approximately the same energy (see table 2). For PBE0 + vdW, the Au 3 · Kr structure with acute-angled Au 3 is 0.06 eV more stable than the other and the Au 3 · Kr 2 structure with acute-angled Au 3 0.10 eV more stable. With XYG3, these values become 0.20 and 0.24 eV, while for RPA + rSE@PBE 0.13 and 0.23 eV. For heavier RG atoms, in particular for Xe, we find that the acute-angled isomer becomes even (slightly) more stable than the obtuse-angled one. The final structure of Au 3 · Kr is an isosceles Au 3 triangle with the two equal-length bonds of 2.6 Å (i.e. 0.02 Å contracted with respect to the bare isomer) and 11 There has been some confusion in literature, when referring to this 'obtuse-angled' triangle. In [6,7] an 'obtuseangled' triangle is mentioned, but it turns out that its internal angle is not bigger than 90 • , it is just bigger than 60 • . It looks like the terms 'obtuse' and 'acute' have been referred to the equilateral reference, rather than the usual right angle. The first clear reference to the 'obtuse' triangle as the one with internal angle about 140 • wide we find in [8]. 12 Actually, when spin-orbit coupling is not considered, as in this paper, two nearly equilateral isosceles triangles are identified as local minima, with the internal angle between the two equal length bonds of 66 • and 56 • , respectively. These are determined by a Jahn-Teller distortion of the perfectly equilateral triangle. By considering spin-orbit coupling [10], the degeneracy that causes the Jahn-Teller distortion is removed and only the perfectly equilateral triangle is found. In this paper we consider only one obtuse-angled triangular isomer. This is justified by the fact that, when one Kr is chemisorbed to the acute-angled triangle, only a 64 • triangle is found; when two Kr are chemisorbed, only a 56 • triangle is stable. What the FIR-MPD experiment sees is thus only one acute-angled isosceles triangular species at a given number of adsorbed Kr, while the isolated perfect equilateral Au 3 is not observable in a FIR-MPD experiment. Furthermore, when Kr is adsorbed to the acute-angled isomer, by binding to one of the Au atoms, the final structure is the same, whether the bonding Au atom is initially at the 66 • corner or one of the two 57 • . the angle between them 64 • wide (i.e. 2 • smaller than in the bare cluster); Kr is bonded to the gold atom at the 64 • vertex, with a Au-Kr distance of 2.74 Å. When Kr is adsorbed to the obtuse-angled isomer, it is chemisorbed to one of the two one-fold coordinated Au atoms. The Kr-Au bond length, 2.94 Å, is much longer than in the acute-angled case; it follows that the Au 3 -Kr interaction energy is about half than for the acute-angled isomer. As a consequence, by adding one Kr atom, e.g. by letting Kr approach toward the central (two-fold coordinated) Au atom, the obtuse-angled Au 3 isomerizes into the acute-angled isomer. This would be an unusual example of a RG-induced isomerization of a metal cluster, but similar to the observations for Cu 3 RG [58].
Both Au 3 isomers can also bind two Kr atoms, with binding energy slightly smaller than double the binding energy of one Kr (see table 2). A third Kr atom on the acute-angled isomer of Au 3 is only vdW bound. This Au 3 · Kr 3 has the geometry of the acute-angled Au 3 · Kr 2 , with the third Kr in plane, but only vdW-bonded at 4.0 Å from the third Au atom.
For both Au 3 isomers, there are also other equilibrium positions for Kr. For the lowest energy structure besides the chemisorbed sites, i.e. the complex with the acute-angled triangle, Kr is out-of-plane and its trace is on the center of mass of the gold cluster, the interaction energy is −0.05 eV, which is about 25% of the bonding energy upon chemisorption. As it is easy to predict, finite-temperature MD simulations find Kr localized at the bonding site(s).
While the IR spectrum of bare Au 3 is dominated by one intense mode at 95 cm −1 , related to the antisymmetric stretching of the acute-angled triangle, many peaks appear in the IR spectrum when one or two Kr atoms are adsorbed. Some of these peaks are associated with eigenmodes that correspond to IR-inactive modes of the pristine cluster, whereas the change of symmetry provoked by Kr adsorption makes them visible; some are new modes involving Kr as well.
The finite-temperature theoretical spectrum of Au 3 · Kr (figure 3) agrees well with the experimental spectrum when only the acute-angled isomer is considered 13 . In fact, the obtuseangled isomer would have a peak at ∼130 cm −1 , which is not present in the experimental spectrum. It turns out that, even if the energetics of the obtuse-angled isomer is close to that of the acute-angled, the vIP of the obtuse-angled isomer is calculated to be close to the energy of the UV laser (7.9 eV), and thus it may not be (or not efficiently) ionized, while the vIP of the acute-angled isomer is far below the photon energy. At low wavenumbers (50-70 cm −1 ), though, the experimental spectrum does not show the band predicted by theory. The reason for this behavior is that also Au 3 · Kr 2 absorbs at those frequencies (see below). However, at such low photon energies many Au 3 · Kr 2 appear to loose only one Kr in the photodissociation and the formation of Au 3 · Kr compensates for the dissociated fraction.
Provided that the relative energy of the obtuse-angled di-krypton complex is close to the energy of the acute-angled, the calculated vIPs suggest that the di-krypton complexes of both isomers can be ionized at 7.9 eV and could contribute to the experimental spectrum of Au 3 · Kr 2 , although the obtuse-angled may be slightly less efficiently ionized as it has the higher vIP. Indeed, a superposition of the two theoretical spectra shows a remarkable agreement with the experimental one ( figure 4). For instance, the broad band around 50-60 cm −1 is reproduced quite well and shows a similar substructure. The subpeaks are an anharmonic feature, since only one peak per isomer is found in the harmonic spectrum in that region. Inclusion of the obtuse-angled isomer can explain the extension of this band toward lower frequencies.

Au 4 · Kr 2
Similarly to Au 3 , Au 4 has two low energy isomers, a rhombus and a Y-shaped cluster, with a difference in energy of 0.02 eV, which is further reduced by the adsorption of one or two Kr atoms. Au 4 · Kr is not ionized in the experiment, and indeed the calculated vIPs (table 2) are consistent with this observation. The calculated vIPs of the two Au 4 · Kr 2 isomers suggest that the rhombus isomer is more efficiently ionized (vIP = 7.5 eV), as the vIP of the Y-shaped isomer is with 7.7 eV already rather close to the photon energy. The FIR-MPD spectrum for Au 4 · Kr 2 (figure 5) is well reproduced by the theoretical finite-temperature spectrum of the rhombic isomer. Inclusion of a fraction of the Y-shaped isomer could explain a further broadening of the low frequency peak, however, there are no signs of the other, though less intense, bands predicted for this isomer.
If the Kr atom is placed above the plane of the rhombic Au 4 , a vdW complex is formed with E vdW (Au 4 -Kr) = −0.08 eV, which is less than half the energy of the localized bonding. For comparison, we calculated the MP2 + vdW interaction energy for this complex 14 , which is with −0.07 eV in very good agreement with the PBE + vdW value. The vdW complex Au 4 · Kr 2 has the two Kr atoms symmetrically above and below the plane of Au 4 . The interaction energy is still −0.08 eV per Kr, again less than half than the bonding energy of the bonded Au 4 · Kr 2 . Even in this case at finite temperature the Kr atoms are practically always found at the bonding site(s).  In conclusion, theory predicts that Au 4 · Kr 2 is present as a mixture of two low lying isomers, which are nearly equally present. The experiment clearly identifies the rhombic isomer, but there is no compelling verification of the Y-shaped structure, which may be explained by its low ionization probability.

Orbiting Kr: the case of Au 7 · Kr
For clusters bigger than Au 4 , the chemisorption sites for Kr are still present (always single coordinated adsorptions to perimetral two-or three-fold coordinated Au atoms), but, starting from Au 5 , the bonding energy is weakened to ∼ − 0.1 eV in the most favorable geometry. Since with increasing size the number of Au-Kr pairs that show a significant vdW attraction increases, the total vdW interaction between the gold cluster and the Kr atom(s) increases and grows comparable to the bonding energy to a specific site. For planar clusters, this interaction is maximized when Kr is out of plane, and the trace of the Kr position onto the cluster plane lies near the center of mass of the cluster. We find that at size 5 the interaction of Kr sitting at the best bonding site starts to compete energetically with the (vdW) interaction energy of the purely vdW bound complex. In facts, for Au 5 the purely vdW bound Kr yields a E vdW equal to −0.09 eV. Interestingly, a E vdW of about −0.1 eV is also the strongest interaction we found for a Au N · Kr complex. We tested up to Au 20 (which is a perfect tetrahedron, with four triangular {111} surfaces each made of ten atoms [19]) for which we find a E vdW (Kr above the center of one of the faces) of −0.11 eV. The fact that the vdW interaction between Kr and the cluster saturates with the cluster size is due to geometrical reasons, but also  to the polarizability of the Au atoms in the cluster, which, at least in the size interval that we have probed, decreases with the cluster size. The comparably strong vdW binding of Kr to the planar fcc sites is well in line with the experimental findings for anionic Au clusters, where Ar binding has been used to discriminate between 2D and 3D structural isomers [13]. Au 7 is the first cluster size larger than Au 4 for which we have a clearly structured experimental spectrum to compare to (Au 5 · Kr shows only relative week features in the spectrum, and no Au 6 · Kr M is ionized in the experiment). The IR spectrum of Au 7 had been the subject of analysis before and its structure has been identified as planar edge-capped triangle [19], which is here confirmed as the clear global minimum.
We find that Kr can bind, within the plane of the Au 7 cluster, to single Au atoms, similarly to what was described before [19]. The strongest binding ( E b = −0.10 eV) is in the two geometries shown in the top part of figure 6. However, the vdW complex, with Kr above the Au 7 plane, depicted as third from the top in figure 6 has an interaction energy of −0.09 eV, i.e. comparable to the bound case. For this vdW complex, too, we have checked the accuracy of the PBE + vdW interaction energy with the MP2 + vdW approach. The latter predicts an interaction of −0.09 eV as well.
The experimental FIR-MPD spectrum of Au 7 · Kr is already well reproduced by the harmonic spectrum of the bare cluster and in particular the peak positions are in excellent agreement (see figure 6). Nonetheless, the relative intensities of the peaks, in particular the more pronounced ones at 165, 186 and 201 cm −1 , do not match. When the harmonic spectra of the complexes with chemisorbed Kr are considered, it is found that the peak positions do not shift significantly. But the relative intensities do change, in such a way that for one of the in-plane IR intensity / arb. u. Figure 6. Theoretical harmonic IR spectra of Au 7 , Au 7 · Kr and their calculated finite-temperature IR spectra at T = 100 K (T = 96 K) compared to the experimental FIR spectrum of Au 7 · Kr (lower panel). The lower right structure depicts the isosurface enclosing the region where Kr is found 80% of the time during a 0.5 ns MD run at T = 100 K, when forming a vdW complex with Au 7 . Between the atoms surrounded by a square the weakest bond in Au 7 is formed (inset lower panel, see text for details). The Au-Kr distances are in Å.
binding sites an optimal matching with the experiment is found (see figure 6), as already noted in [19]. Due to the competing energetics between the bonded and the vdW complex, though, the picture suggested by MD is slightly different.
In a MD run at the experimental temperature of 100 K, the Kr atom, even when prepared in an initial position at one of the bonding sites, soon starts to orbit around the planar cluster, with a preference for the 'polar' regions (if the planar Au 7 is regarded as the equatorial plane of the approximate sphere onto which Kr slides). The lower right structure in figure 6 shows the isosurface that encloses the region in which Kr spends 80% of its time during a 0.5 ns long MD run. In practice, the simulation box is divided in small cubes and for each cube the average Kr-density is evaluated as the (normalized) number of times the Kr nucleus is found in the cube during the MD sampling. The enhanced density at the polar regions can be interpreted as a clear preference for forming the vdW complex, despite its energetic quasi-degeneracy with the localized bonding situations. This is easily understood on entropic grounds. The vdW complex offers a large number of energetically degenerate levels, as shown by the extension of the isosurface shown in figure 6. In contrast, when Kr is localized at an adsorption site, the system visits only a small number of configurational states, just because of the localization! Since entropy is a measure of the number of the states accessible to the system at a given temperature, the vdW complex has a larger (configurational) entropy and thus a lower free energy compared to the localized bonding case.
The finite-temperature spectrum obtained from a MD simulation of Au 7 · Kr is shown in figure 6 (second lowest panel, continuous trace). The relative peak heights are well reproduced. Furthermore, we find that the finite-temperature spectrum of the bare cluster (same panel, dashed line) is very similar to the spectrum of Au 7 · Kr. In particular the correct peak height ratio for the three peaks at higher frequency is found. This suggests that the Kr atom is not significantly affecting the spectrum, but that the differences to the harmonic spectrum of Au 7 are related to an intrinsic behavior of Au 7 . The marked broadening of the highest frequency peak in the spectra of Au 7 and Au 7 · Kr is indeed an anharmonic feature. Analysis of the MD trajectory reveals that the internal bond of the inner rhombus (i.e. the bond between the highlighted atoms in figure 6) is the weakest, in the sense that the variance of its length is about twice than the average variance of the other bond lengths. Elongation of this bond implies shortening of the distance between the other two atoms belonging to the inner rhombus of Au 7 . Furthermore, a trajectory at higher temperature reveals that this isomer undergoes a fluxional transformation by swapping the role of the atoms arranged in the inner rhombus, i.e. the long and short diagonal interchange and an isomer with identical topology, but with scrambled atoms is formed. This feature will be analyzed in detail in a subsequent publication, by comparing the fluxional behavior of this cluster with other similar behaviors of larger gold clusters. At this point, we just note that the anharmonic broadening of the highest frequency peak in Au 7 · Kr (as well as in pristine Au 7 ) is related to this fluxional behavior. The theoretical understanding of Au 7 (·Kr) vibrational spectrum is not in contradiction with the analysis reported in [25]. In fact, in the work of Mancera and Benoit the analysis was carried out for small distortions of the global minimum structure (the same as ours) at T = 0 K. Only including larger distortion as sampled in a canonical MD trajectory, the anharmonic features of Au 7 become visible.

Conclusions
We report the FIR-MPD vibrational spectra of small Au N · Kr M complexes and provide their assignment by simulating finite-temperature spectra via DFT (with vdW tail correction) molecular dynamics. This approach led us to the identification of the structural information of the considered species. For the MD simulations, we have used the PBE + vdW functional, but some static properties of the Au N · Kr were compared to higher level methods, for validation of our approach. In particular we have tested the validity of the TS scheme [32] for the vdW tail correction of the PBE functional for the system considered here against the MP2 + vdW [26] results. We always found a remarkably good agreement between the adopted method and the higher levels one.
Similarly to what was found for Ar adsorbed on cationic Co clusters [22] or Y doped gold clusters complexed with Xe [20,21], krypton is not a mere tag and does perturb the vibrational spectra of small neutral gold clusters, sometimes even largely. While a somewhat strong electrostatic interaction between a charged cluster or molecule and a RG atom is predictable, an interaction between a RG and a neutral gold cluster, strong enough to visibly perturb the vibrational spectrum of the pristine cluster, may not be expected.
When Kr binds to Au 2 , Au 3 and Au 4 , it is found localized at its binding sites also at the experimental temperature of 100 K and vibrational spectra are greatly affected by Kr adsorption, i.e. Kr does not act as a mere messenger for the detection. Therefore the interpretation of the vibrational spectra has to consider the whole Au N · Kr M complex. Theory predicts that Au 3 · Kr, Au 3 · Kr 2 and Au 4 · Kr 2 appear in two different isomers each. However, in all these cases the experimental spectra can be explained also without assuming the presence of the higher energy isomer. The latter may be invisible due to low population or because of poor ionization yield.
In Au 7 · Kr, the Kr atom is not bound to a single site but delocalized with a clear preference for being above/below the Au 7 plane. We also find that anharmonic effects leading to unusual broadening of the peaks or new peaks due to mode interactions are present already at comparably low temperatures (T ∼100 K). In the case of Au 7 · Kr, the inclusion of anharmonicity results in the best agreement with the experimental spectrum. Nevertheless, the comparison with the harmonic spectra demonstrates that these can be sufficient to establish the metal cluster structure as long as one realizes intensities may be perturbed.
The role of the Kr atom(s), i.e. localized versus orbiting around the cluster, and its influence on the measured vibrational spectra, is fully revealed only by the statistical sampling of the canonical ensemble, due to the fact that the system is investigated at finite, albeit relatively low, temperatures.
We are currently extending the calculation of finite-temperature vibrational spectra via molecular dynamics to larger Au clusters, where at some sizes fluxional behaviors (i.e. relatively frequent structural interchanges between neighboring isomers) must be taken into account for a full understanding of the FIR-MPD spectra. The circles on the x-axes mark the experimental peaks (atT = 100 K). Along the sets of diamonds, the temperature is changed fromT = 0 K (harmonic analysis) to T = 100 K (T = 96 K). Along the sets of squares the grid is changed from 'light' to 'really tight 974', i.e. including the extended angular grid [30]. Along the sets of downward triangles the basis set is changed independently for the two elements, from 'tier1 − gh' for Au and 'tier1 − f ' for Kr to 'tier2' for both elements. Along the upward triangles the timestep is changed from 10 to 1 ps. The labels for the basis set are 'T1' for 'tier1', 'T2' for 'tier2', while 'tier1 − gh' means that a g and a h atomic basis function are missing from 'tier1'. On the right the status of the settings that are not changed along the sets are specified. In that case, when the basis set is marked with just 'T2', it means 'tier2' for both Au and Kr. associated to the lowest-frequency of these modes appears in the experimental spectrum (see figure 4) as divided into three sub-peaks (we assign the fourth sub-peak at the lowest frequency to the obtuse-angled isomer of this cluster). For the analysis reported in figure A.1, we have reported the sub-peak at highest frequency, i.e. at 63 cm −1 . The position of the three experimental peaks (the other two are at 88 and 190 cm −1 ) are marked on the x-axes of the plot as filled circles.
This particular isomer was chosen because it has IR visible modes near the extremes and in the middle of the (experimentally accessible) frequency interval spanned by the gold clusters. There are also modes at frequencies lower than 50 cm −1 , but they fall below the experimental window. Right above the experimental points, the positions of the theoretical peaks at various temperatures, all the other settings being fixed, appear. One can note that the peaks are generally red-shifted at increasing temperature, but, more strikingly, the change of frequency as a function of temperature is not the same at all frequencies: lower frequency modes shift more than higher frequency ones. This is one of the effects of anharmonicity. One can note that in the harmonic approximation (i.e. at T = 0 K) the lowest frequency mode is closer to the experimental value than the highest frequency mode to its experimental counterpart. As a consequence, in literature an empirical scaling factor for aligning the theoretical and experimental numbers is often employed . AtT = 100 K the three peaks are then red-shifted by different amounts from their positions calculated within the harmonic approximation, in such a way that a (rigid) blue-shift would bring them in approximate correspondence with the experimental values. The necessity of a rigid shift rather than a scaling factor for theoretical finite-temperature spectra has been already noted by other authors (see e.g. [42] and references therein). Here we make the systematic observation that (a) when the density of the integration grid is increased (set of squares in figure A.1), (b) when the size of the basis set is increased (set of downwards triangles) or (c) when the MD timestep is decreased (set of upward triangles), i.e. whenever the accuracy of the evaluation of the potential-energy surface improves, the position of the peaks are always blueshifted, i.e. toward the experimental values. These shifts are also approximately rigid (except when passing from the least accurate basis set or grid to the next step) across the frequency window. It has to be noted that the sensitivity of the peak positions toward the settings is quite small, unless when passing from the coarsest 'light' integration grid to the next level ('tight'), or when passing from the very small basis set ('tier1' for Kr, without the f function) to the next level.