Quantification of finite-temperature effects on adsorption geometries of $\pi$-conjugated molecules

The adsorption structure of the molecular switch azobenzene on Ag(111) is investigated by a combination of normal incidence x-ray standing waves and dispersion-corrected density functional theory. The inclusion of non-local collective substrate response (screening) in the dispersion correction improves the description of dense monolayers of azobenzene, which exhibit a substantial torsion of the molecule. Nevertheless, for a quantitative agreement with experiment explicit consideration of the effect of vibrational mode anharmonicity on the adsorption geometry is crucial.

Precise experimentally determined structures of large organic adsorbates are indispensable -for the detailed understanding of their wide-ranged functionalities, but also for the benchmarking of ab initio electronic structure calculations [1][2][3]. For large molecules with polarizable π-electron systems, van der Waals (vdW) interactions are substantial and may critically influence the adsorption geometry [4][5][6][7]. Accounting for these interactions in ab initio calculations remains a challenge, and different approaches to this problem at varying degrees of accuracy are currently explored [8][9][10][11][12][13][14]. Due to the system sizes inherent to large molecular adsorbates, efficient semi-empirical dispersion correction schemes to density-functional theory (SEDC-DFT) are particularly promising [13]. However, their approximate nature makes them even more dependent on reliable experimental benchmark structures. This holds in particular for adsorption at metal surfaces, where the non-local collective substrate response (many-body electronic screening) requires advancements beyond the traditional pairwise summation of vdW interactions in these schemes [14,15].
With SEDC-DFT now striving for the approximate inclusion of the collective substrate response [14], the accuracy increases to approximately 0.1Å for the predicted adsorption heights [14,[16][17][18]. At this level of accuracy, a new issue arises: Experiments for structure determination are often carried out close to room temperature, while in SEDC-DFT the ground state (at 0 K) is normally calculated. The complex internal vibrational structure of large organic adsorbates which may sensitively influence the experimental time-averaged geometry is thus neglected. In this Letter we show that the inclusion of such thermal expansion effects into SEDC-DFT is indeed necessary to reach quantitative agreement between experiment and theory. Hence, benchmarking at the current level of sophistication requires the careful analysis of finite-temperature effects. Otherwise misleading conclusions with respect to the SEDC-DFT accuracy might be obtained.
Our experiments have been carried out on azobenzene (AB, cf. Fig. 1a) adsorbed at Ag(111), by the normal incidence x-ray standing wave technique (NIXSW). AB is a widelyused molecular switch [20]. Investigations of its substrate interaction are driven by the challenge to preserve the switching functionality in the presence of a surface. In this context, the knowledge of the adsorption structure is essential. NIXSW is an established method to determine the adsorption geometry (in particular adsorption heights) of large organic adsorbates [21,22]. The AB/Ag(111) system has been studied by NIXSW before and the results were compared to the SEDC-DFT approaches of the time to conclude on the importance of (then untreated) electronic screening effects [15]. Using a most recent SEDC-DFT revision that approximately includes non-local collective substrate response we here confirm this proposition. Also, accounting in our refined analysis for coverage dependence, we nevertheless show that the crucial missing link to achieve quantitative agreement with experiment lies not on the electronic structure level, but in hitherto generally neglected finite-temperature effects.
NIXSW experiments were carried out at ESRF, beamline ID32, under ultra-high vacuum conditions (≈ 5 × 10 −10 mbar) [15]. The Ag(111) surface was cleaned by several cycles of Ar + ion sputtering and annealing at 820 K. Multilayers of AB were deposited from an effusion cell held at room temperature onto the atomically ordered Ag(111) crystal at 220 K. AB monolayers were then prepared by desorption from multilayers, by heating with a rate of 1 K/s until the multilayer desorption peak had decayed and before the monolayer peak was observed [23]. This guarantees coverages close to one monolayer. The AB fragment mass of 77 amu (C 6 H + 5 ) was monitored on-line with a quadrupole mass spectrometer. The Ag crystal was kept at 210 K during the NIXSW experiments to prevent desorption.
Vibrations are expected to influence the average geometry of the adsorbate (via vibra- tional mode anharmonicity) and to broaden the distribution of atoms around their average positions (via vibrational dynamics) [21]. While this will affect both the coherent position P c and the coherent fraction F c of the NIXSW signal, prevalent (harmonic) Debye-Waller theory only considers temperature effects on F c [21,24]. P c defines the average adsorption height of a species, while F c quantifies the corresponding height distribution. A coherent fraction of 1 means that all photoemitters of a certain species have precisely the same adsorption height above the relevant family of Bragg planes, while a coherent fraction of 0 indicates a homogeneous distribution of the photoemitters throughout the Bragg spacing. In general, F c < 1 due to unavoidable structural disorder [25], adsorbate and substrate thermal vibrations [26]. However, the coherent fractions of different chemical species also contain information about the internal geometry of the adsorbate that has so far been left aside in most NIXSW studies. Here we recover this information by including differences between the F c of different species into our analysis.
In the present case of AB/Ag(111), NIXSW provides coherent positions P C c = 0.27±0.02, [27]. Therefore, while the respective coherent positions are identical within the errors, the coherent fraction of C is 29% smaller than the one of N. In our refined structure determi-  nation, we ascribe this difference to the internal geometry of AB, assuming that F C i c and F N i c , the coherent fractions of individual C and N atoms, are equal (and smaller than 1 due to disorder) [26]. To solve the AB structure, two internal degrees of freedom are considered: the tilt angle ω (Fig. 1a, [28]) and the torsion angle β (Fig. 1b), defined as dihedral angles CNNC and CCNN [19], respectively. It is impossible to explain the ratio F C c /F N c in a model in which ω is the only internal degree of freedom of the molecule, because any distortion along ω that would lead to a decrease of F C c would at the same time result in an increase of the coherent position P C c which is related to the average adsorption height of the carbon atoms. Hence, an additional degree of freedom must be considered to explain the measured NIXSW structure parameters. A plausible choice is the torsion angle β, because for small angles ω a finite β would broaden the carbon distribution essentially without changing the average carbon height (Fig. 2 magenta curve). Note that this broadening could in principle be due to a static distortion of the molecule and/or due to its vibrational dynamics. However, for a purely dynamical reduction of the average coherent fraction F C c by 29 % an unreasonably large C vibrational amplitude of the order ±0.30Å (with fixed N atoms) would be required. Therefore, we will first consider a static distortion before coming back to a possible dynamical contribution.
and constructing the molecular geometry such that the measured values for P C c , P N c , F C c , F N c are obtained, we find an adsorption geometry with d N−Ag of 2.97 ± 0.05Å, a tilt angle ω of −0.7 • and a torsion angle β of 17.7 • from our NIXSW data (cf . Table I) [29,30].
A torsion angle β of more than 17 • is difficult to rationalize for a single molecule adsorbed on the surface without neighbors. Yet, all calculations so far [15,31,32] have been carried out for single molecules (while our experiment is performed on a condensed layer, see above). We therefore need to analyze the coverage-and packing-dependence of the adsorption geometry of AB/Ag(111) theoretically. While our previous SEDC-DFT calculations [15] for this system were carried out at the level of the dispersion-correction scheme by Tkatchenko and Scheffler (TS) [12], we now employ the more recent vdW surf scheme [14], which accounts for non-local collective substrate response via renormalization of the dispersion coefficients on the basis of Lifshitz-Zaremba-Kohn theory. Details of the calculations can be found in the supplement [26]. We determine the optimized adsorption geometries for a range of different surface unitcells [26], with one AB per (6 × 7) cell representing the low-coverage (LC) limit and two AB molecules in a (2 × 5) cell leading to the highest considered molecular surface density (cf. Fig. 3).
The PBE+vdW surf results compiled in Fig. 3 show that the adsorption geometry indeed varies substantially with increasing molecular surface density. While in the LC limit the adsorbed molecule is essentially flat (Fig. 3b), the tilt and torsion angles ω and β increase with the packing density (Fig. 3c-d). As a consequence of the internal distortion of the molecule, the vertical adsorption height of the azo-bridge also increases [26]; this tendency continues beyond the critical coverage of 1.56 AB · nm −2 at which the now nearly upright molecules start to flatten out again within the increasingly dense overlayer (Fig. 3d). The flattening of the upright molecule implies an increasingly asymmetric position of the azobridge, i.e. different vertical adsorption heights of the two nitrogen atoms, with a consequent lifting from the surface [26].
With most of the adsorption energy of the flat AB molecule in the LC limit coming from dispersive interactions with the substrate, the binding energy per molecule naturally decreases in the distorted high density structures (cf. supplementary Table I [26]). Due to the denser packing, the adsorption energy per surface area E ads /area nevertheless increases and reaches a maximum at 1.17 AB · nm −2 (Fi.g 3a). The intermolecular vdW interactions in the high density phases further increase E ads /area, which reaches a second maximum at a density of 1.87 AB · nm −2 . Our calculations therefore predict the existence of two optimum packing densities, a phase A (Fig. 3c) at 1.17 AB · nm −2 and a phase B (Fig. 3d) at 1.87 AB · nm −2 . Qualitatively similar findings are obtained with the TS scheme. In detail, however, there are decisive differences. For example, TS fails to predict the maximum of E ads /area corresponding to phase A, cf. Fig. 3a.
To decide which structure -if any of the above -we have in our experiment, we take the calculated ground state geometries and compare them to experiment (Table I). Phase LC can be ruled out, both because of its small torsion angle, and because of our sample preparation procedure which yields dense layers. The average adsorption height of N atoms in phase B is 4.39Å, which is inconsistent with the experimental value of 2.97Å ormodulo a Bragg spacing -5.32Å. We therefore conclude that our NIXSW experiment has been carried out on a structure similar to phase A. This conclusion is consistent with the expectation that neither of the two dispersion-correction schemes will work reliably at the packing density of phase B that is close to the one of the molecular crystal, in which even the  vdW surf scheme will overestimate lateral interactions [33], because higher-order many-body terms are neglected [34]. This neglect will contribute to a spurious stabilization of phase B in the calculation.
In Table I  has a large impact on the predicted adsorption height d N−Ag . It tends to improve the TS prediction, although the height is still not perfect, and the calculated ω is too large.
We will now show that the predictions of the vdW surf theory can be substantially improved toward a quantitative agreement with experiment, if the effect of finite temperature is taken into account. In particular, anharmonic contributions to the vibrational motion may modify the time-averaged geometries that are experimentally observable [21]. We demonstrate this by explicitly calculating the harmonic vibrations for the adsorbed AB molecule at the optimum density of 1.17 AB · nm −2 , both at PBE+vdW surf and PBE+TS levels. We then map out the anharmonic regimes of these modes at energies around the experimentally employed 210 K, by distorting the molecule along the corresponding vibrational eigenvectors. Next, we fit a Morse potential [35] to these data points for every harmonic mode and integrate the motion in the Morse potentials analytically to obtain the shifts of the average positions at 210 K relative to the harmonic minima. Summing these shifts over all vibrational modes, we finally construct an anharmonically corrected geometry for the adsorbed AB molecule [26].
With this procedure we arrive at the following finite-temperature To check the self-consistency of the finite-temperature geometry, we have simulated NIXSW results on its basis, with the aim to evaluate the influence of vibrational excitations on the coherence of the NIXSW signal (cf. supplementary Section III [26] for details). Note that our experimental values for ω and β in Table I [26]. ), hence closer to experiment. The nearly equal reduction of F C c and F N c due to thermal vibrations confirms a posteriori that in deriving the experimental structure we can interpret the different experimental F c s as being due to static distortion, and hence the experimental geometry parameters in Table I are the correct reference point for the calculated finite-temperature geometry.
In conclusion, we have analyzed the structure of the archetypal molecular switch azobenzene on the Ag(111) surface. We find that the inclusion of collective substrate response into SEDC-DFT correction schemes is absolutely essential for a correct description of the adsorption geometry. However, since the vdW surf scheme leads to a reduction of both the dispersion coefficients and the van der Waals radii, it may -counterintuitively -decrease the adsorption height compared to SEDC-DFT without collective substrate response. This is clearly observed for the LC phase. However, we have identified two effects which increase the adsorption height again. First, this is the dense molecular packing and the associated molecular distortion, which increase d N−Ag by 0.20Å. Second, the anharmonicity of molecular vibrations raises d N−Ag by another 0.17Å. The remaining disagreements between experiment and theory notwithstanding, there are three important findings which can be generalized: Firstly, information regarding the molecular conformation beyond the average positions of certain chemical species can be retrieved from the coherence of the NIXSW signal with suitable simulations. Secondly, thermal expansion due to the anharmonicity of molecular vibrations not captured in Debye-Waller theory [21] may contribute 0.1 − 0.2Å to the adsorption height. This must be taken into account in future benchmarks of high-level ab initio theory against NIXSW, either by carrying out the experiments at low temperature or by inclusion of finite-temperature vibrational effects into the calculation, as sketched in the present paper. And thirdly, our observation that all three geometry parameters d N−Ag , ω, and β develop into the correct direction if anharmonic effects are included proves that the PBE+vdW surf scheme captures the essential physics of both chemical and dispersion interactions and is therefore a good starting point as a ground state calculation for the adsorption of large π-conjugated molecules.