Article https://doi.org/10.1038/s41467-022-33653-7 Evidence of exciton-libron coupling in chirally adsorbed single molecules Received: 8 April 2022 Jiří Doležal 1,2 , Sofia Canola1, Prokop Hapala1, Rodrigo Cezar de Campos Ferreira 1, Pablo Merino 3,4 & Martin Švec 1,5 Accepted: 28 September 2022 Interplay between motion of nuclei and excitations has an important role in Check for updates molecular photophysics of natural and artificial structures. Here we provide a detailed analysis of coupling between quantized librational modes (librons) 1234567890():,; 1234567890():,; and charged excited states (trions) on single phthalocyanine dyes adsorbed on a surface. By means of tip-induced electroluminescence performed with a scanning probe microscope, we identify libronic signatures in spectra of chi- rally adsorbed phthalocyanines and find that these signatures are absent from spectra of symmetrically adsorbed species. We create a model of the libronic coupling based on the Franck-Condon principle to simulate the spectral fea- tures. Experimentally measured librational spectra match very well the theo- retically calculated librational eigenenergies and peak intensities (Franck- Condon factors). Moreover, the comparison reveals an unexpected depopu- lation channel for the zero libron of the excited state that can be effectively controlled by tuning the size of the nanocavity. Our results showcase the possibility of characterizing the dynamics of molecules by their low-energy molecular modes using µeV-resolved tip-enhanced spectroscopy. Coupling between excited electronic states and nuclear motion is an nanoscopic environment of the chromophores. Therefore, performing essential mechanism for conversion between optical, mechanical and experiments in well-controlled environments on the single-molecule chemical forms of energy in nanosystems. Such excitation-vibration level is a key to advancing our fundamental understanding of mole- coupling is relevant in biological processes such as photosynthesis1–3, cular librations as well as for the development of nanomachines and light-sensitive proteins in eyes4,5, in artificial molecular motors6,7 or nanodevices. organic solar cells8,9. Frustrated rotations (librations) represent a par- Recent progress in tip-enhanced single-molecule spectroscopy ticular type of vibration in which the molecule performs a torsional permits to overcome the limitations of traditional ensemble- oscillation when subjected to external stimuli and constraints that based spectroscopies and study neutral and charged excited restrict its orientation. Despite their efficient coupling to electronic states and their coupling to vibrations on single molecules. transitions, librations have largely eluded spectroscopic detection Scanning tunnelling microscope-induced electroluminescence10–18 because of naturally relatively small energy differences between their (STM-EL), photoluminescence19,20 (STM-PL) and tip-enhanced Raman quantised levels, making it difficult to derive any characteristic para- scattering21–23 (TERS) methodologies operating at cryogenic tempera- meters from the spectra, especially in large ensembles of molecules. tures can localise and amplify the interaction of electromagnetic Obtaining well-resolved spectra of librations directly is complicated radiation with a molecule located in the plasmonic nanocavity, formed because they become easily obscured by stochastic thermal motion, between the scanning probe and a metal sample by many orders of the effect of solvents and by a generally limited control over the magnitude. Using these approaches one can study individual 1 Institute of Physics, Czech Academy of Sciences, CZ16200 Praha 6, Czech Republic. 2Faculty of Mathematics and Physics, Charles University, CZ12116 Praha 2, Czech Republic. 3Catalan Institute of Nanoscience and Nanotechnology (ICN2), CSIC and BIST, Campus UAB, Bellaterra, E08193 Barcelona, Spain. 4Instituto de Ciencia de Materiales de Madrid; CSIC, E28049 Madrid, Spain. 5Institute of Organic Chemistry and Biochemistry, Czech Academy of Sciences, CZ16000 Praha 6, Czech Republic. e-mail: dolezalj@fzu.cz; svec@fzu.cz Nature Communications | (2022)13:6008 1 Article https://doi.org/10.1038/s41467-022-33653-7 photoactive molecules without the influence of stochastic thermal profile, determined from the intensity of Franck-Condon factors of the fluctuations or the presence of solvents and map the electron transi- transitions. tions in the optical near-field with submolecular resolution. This resolution is orders of magnitude higher than what can be achieved Results with in-solution spectroscopy and opens a new window to determine High-resolution STM-EL spectra of single phthalocyanine the mechanisms governing the photophysics of molecular systems. adsorbates Here we apply STM-EL spectroscopy to investigate coupling Zinc-phthalocyanine (ZnPc) and Magnesium-phthalocyanine (MgPc) between charged excited electronic states (trions) and quantised adsorb centred on the Cl− site of 2–5 monolayers (ML) of NaCl on Ag(111) librations (librons) of zinc, magnesium and free-base phthalocyanine and manifest a characteristic 16-lobe appearance in the occupied-state molecules (ZnPc, MgPc and H2Pc respectively), adsorbed on sodium STM images measured at −2.8 V (insets in Fig. 1b, c). This appearance is chloride (NaCl). Phthalocyanines are structurally similar to biological the result of averaging over two geometrically equivalent metastable fluorophores (e.g. chlorophyll), therefore their interaction with the chiral adsorption configurations, between which the molecule rapidly crystalline substrate provides a convenient controllable model for switches upon injection of electrons12,13,16,29,30. The motion between more complex interactions in vivo. Also, phthalocyanines on surfaces these two configurations (called switching16,29 or shuttling12,17) is repre- have been proposed as a model for molecular rotors and switches24–26 sented by the larger grey arrow in Fig. 1a. In contrast, free-base phtha- for their ability to jump between various adsorption geometries on the locyanine (H2Pc) adsorbs centred above the Na+ site and exhibits an surface upon electronic or mechanical excitation. We exploit the fact 8-lobe pattern in the STM image (inset in Fig. 1d). Here, the apparent that adsorbed fluorophores exposed to electric fields in the nanocavity symmetry is due to averaging over the H2Pc tautomers with different show propensity to charge and emit from excited trion states that configurations of the central two H atoms31. STM-EL spectra (see Fig. 1b, generally manifest substantially narrower lineshapes, compared to the c) acquired at bias voltages at −2.8 V on the molecular lobes of ZnPc and emission peaks of the neutral excited states27,28. Up to our knowledge, MgPc show distinct emission peaks, corresponding to neutral exciton this enables for the first time a high spectroscopic resolution suitable (Q) and of the cation trion (Q+). For H2Pc in the neutral state, the excited for studying the fine structure arising from the trion-libron coupling in states with electric transition dipoles oriented along the x and y axes of a molecule with high moment of inertia, which we rationalise using the the molecule (Qx and Qy) are not equivalent, which results in degeneracy Franck-Condon principle and a harmonic oscillator model based on lifting and two observable excitonic lines: Qx and Qy, the former having analysis of DFT and TD-DFT results. Using this approach we can pre- lower energy and higher intensity32. Interestingly, for H2Pc in the cisely extract parameters of the potential energy landscape of the cationic state only a single trion peak is detected. Based on quantum- systems in their respective ground and excited states, the libration chemical calculations (Supplementary Fig. 1 and Supplementary eigenenergies and estimate the probability distribution of the librons Table 1), we assign it to the Qy+ exciton, as the Qx+ is predicted to be of in the excited states. We can establish a general correlation between ~200 meV higher in energy (in the ref. 18, the analogous transition is adsorption configuration (chiral vs. non-chiral) and the spectral labelled as Xx+). a Fig. 1 | STM-EL characterisation of phthalocyanines and high-resolution spec- 1 nA. e–g Spectra measured at the same bias and tunnelling current with 400 μeV tra of their cations. a Scheme of the STM-EL measurement of the MgPc, H2Pc, ZnPc resolution on the ZnPc, MgPc and H2Pc cations, respectively, evidencing the fine molecules (from top to bottom, respectively), performing in-plane librations on the structure present in the first two cases. The scale is given relative to the peaks surface of NaCl/Ag. The red arrow denotes small-energy librations while the grey within the spectral manifold, which manifest lower intensity with respect to their arrow represents switching between two degenerate chiral equilibrium adsorption neighbours. The reference energies E00 are set to the assumed zero phonon lines in positions. b–d Overview STM-EL spectra of the ZnPc, MgPc and H2Pc at −2.8 V, each spectrum. All molecular spectra in b–g are normalised by dividing by their 100 pA, showing the neutral (Q, Qx, Qy) and cation (Q+, Qy+) emission fingerprints corresponding plasmonic background. Raw spectra are presented in Supplemen- plotted with blue solid lines. Grey-filled spectra on the background of each panel tary Fig. 2. Source data are provided as a Source Data file. are the responses of the nanocavities measured on a clean Ag(111) surface at 2.5 V, Nature Communications | (2022)13:6008 2 Article https://doi.org/10.1038/s41467-022-33653-7 The full width at half maximum (FWHM) of the ZnPc and MgPc vibronic transitions between excited and ground states in molecular neutral Q peaks is typically 8–20 meV (Supplementary Fig. 3 and ref. systems (Fig. 3a, b). Accordingly, the intensity Imn of an emission peak 33) depending on the exact NaCl thickness, tip-sample separation32 associated with the decay from the electronically excited librational and nanocavity plasmon-exciton matching34,35. Conversely, the Qx, Qy state m (ψ1m ) to the electronic ground librational state n (ψ0n ) can be linewidths of neutral H2Pc are narrower - as low as 4 meV20,35, but still approximated as: several orders of magnitude larger than a homogeneous broadening on a comparable system at 6 K36,37. The width of the ZnPc+ and MgPc+ 2 trion envelopes are in the range of 5–7 meV, however, high-resolution I mn / μ2 wm <ψ1m ∣ψ0n > ð2Þ spectra reveal a rich fine structure. When measured at 400 μeV reso- lution, the spectra manifest a manifold of narrow, nearly regularly where μ is the electronic transition dipole moment between the 2 spaced peaks (Fig. 1e, f), where typically one peak in the central part of ground and excited electronic states, and <ψ1m ∣ψ0n > is the squared the manifold shows less intensity than its neighbours. Surprisingly, modulus of the overlap integral between the wavefunctions ψ of such fine structure is absent in the high-resolution spectrum of the librational states, i.e. the Franck-Condon factors. wm describes the H2Pc+ trion (in Fig. 1g) which shows only a main line, accompanied by a probability of the system to be in the initial state m of energy ε1m at the second minor component, likely originating from the second tauto- moment of emission, and we model it by an exponential distribution, mer. Importantly, the spectral manifolds are appearing independently using an effective temperature Teff of the system: of the rotational switching motion of the molecules, as evidenced by the spectrum in Supplementary Fig. 4, obtained on MgPc on 4 ML   1 ε1 NaCl/Au(111), which is stabilised in a chiral adsorption geometry on a wm =  exp  m ð3Þ Z k B T eff step edge. P Theoretical model of the librations and fitting of the spectra where kB is the Boltzmann constant and Z = expðε1i =k B T eff Þ is the i The spectral fingerprints of the ZnPc and MgPc trions comprising + + partition function, whose sum runs over all libronic states i. multiple peaks indicate an efficient coupling between the molecular Next, we simulate the emission spectrum as the sum of all emis- libration transitions and electronic transitions, i.e. the energy differ- sion lines (up to m, n = 20) with energy shifts ε1m - ε0n and intensity Imn, ence generated by change of the quantised energy of libration modi- convolved with a dressing Gaussian (Lorentzian for H2Pc+, see Methods fies the energy of the exciton decay. To estimate the potential energy section) function to account for the additional spectral broadening, landscape that hosts the librations, we first perform calculations of the with a parametrically imposed full width at half maximum (γ). total energy E dependence on azimuthal angle ϕ for the planar mole- In order to yield an estimation of physically relevant values for cule cations in their doublet ground and first excited states (denoted each studied case we perform a qualitative comparison between the as D0, D1). The results of the E(ϕ) calculations are summarised in Fig. 2. measured spectra and the spectra simulated with the harmonic model, We find that ZnPc+ and MgPc+ in the ground and trion states have using different combinations of the potential stiffnesses k1, k0 and the double-well energy landscapes (Fig. 2a, b), separated by a barrier of angular displacement Δϕ0 (shown in Supplementary Fig. 5, for more ~200 meV (for the switching of the molecular geometry), with equili- details see Methods section). Thus, for Δϕ0 ≈ 0 the simulation resem- brium azimuthal angle value of ϕ0 ≈ ±15°, somewhat higher compared bles the spectrum of H2Pc+, and for nonzero Δϕ0 and a k1/k0 ratio to previously reported results29,38. The H2Pc+ energy landscape consists above 1.10 approximates the spectra observable in the case of ZnPc+ of a single well with the lowest energy configuration at ϕ0 = 45° and MgPc+. Although the envelope and energetic distribution of the (Fig. 2c). For the transitions D0 → D1 in ZnPc+ and MgPc+, parabolic fits simulated peaks are in a close agreement with the experiments, the of E(ϕ) around the equilibrium angle (Fig. 2d and Supplementary relative peak intensities for ZnPc+, MgPc+ cannot be perfectly repro- Table 2) quantify a shift in the equilibrium angle Δϕ0 of about 0.3° and duced for any combination of parameters using the initial state a change in the stiffness of a few percent. The nonzero Δϕ0 is a result of population probability defined in Eq. (3). In particular, the overall the different asymmetry of substrate electrostatic interaction acting intensity of the central peaks (including the zero-phonon line with between the NaCl substrate and the ground and excited states of the energy E00) generated by the m = n transitions is significantly lower chirally adsorbed molecules. H2Pc+ also shows a comparable change in than the one originating from m−n = ±1. We suggest that this can be the stiffness of the potential well upon excitation, but at the same time caused by an efficient depopulation channel for the lowest energy it does not rotate its equilibrium configuration (i.e. Δϕ0 = 0), due to its libron mode of the system in the excited state (i.e. the upper parabola symmetrical adsorption geometry with respect to NaCl. in Fig. 3a). To reflect this in the theoretical model we modified the Having learned from the calculations that at the lower energy limit distribution w’m in Eq. (3) by reducing the probability of the zero-level (well below the barriers) librating molecules can be treated as har- libron excited state, using a reducing factor A < 1 (w0 0 = A  w0 ) and monic torsional oscillators, we created a simple model for simulating fitted the experimental spectra (parameters optimised by fitting are the spectra. This model, based on the Franck–Condon principle, will summarised in Table 1). allow us to extract relevant quantities of each studied system by fitting. H2Pc represents a case in which the molecule does not undergo That is, using a harmonic molecule-surface interaction potential V(ϕ), any change of the equilibrium angle upon the excitation due to the defined by the stiffness according to the electronic state (k0 for ground symmetrical adsorption geometry (Fig. 2d). As a result, the m ≠ n state and k1 for excited state) and moment of inertia J, we can deter- transitions are mostly forbidden, except for the ones with non- mine the corresponding librational eigenenergies (ε0i , ε1i , respectively) negligible Franck-Condon factors resulting from the change in the and wavefunctions ψ0i , ψ1i separately in the ground and excited state as stiffness of the potential (k1/k0 ≠ 1). Conversely, due to non-negligible a function of ϕ, by solving numerically the stationary 1D Schrödinger Δϕ0, fitting of the rich spectral features on ZnPc+ and MgPc+ allows a equation precise determination of the k0, k1, Δϕ0 parameters (Fig. 3c, d) which " # are in accord with the values estimated from the calculations (see 2 Supplementary Table 2 for comparison). The effective temperature _2 ∂ V ðϕÞ  ψ ðϕÞ = εi ψi ðϕÞ ð1Þ above 50 K resulting from the fitting indicates that the excess energy of 2J ∂ϕ2 i inelastic tunnelling electrons of a few eV can excite high librational For the description of exciton-libron coupling we employ the states and create a transient initial state population above the zero- Franck-Condon principle in analogy with the description of vertical level libron state39. Nature Communications | (2022)13:6008 3 Article https://doi.org/10.1038/s41467-022-33653-7 Fig. 2 | TD-DFT evaluation of the potential energy landscape for rotation of and ½110 NaCl crystallographic direction, as shown in the inset of a). d Detailed phthalocyanines on NaCl. a ZnPc+, b MgPc+ and c H2Pc+ total energy (relative to comparison of the potential well minima of the three cationic chromophores the energy of the minima E0min) as a function of rotation by angle ϕ for the ground (computed as a difference with the energy of their respective minima E0min or E1min) and excited states. The computed energy is plotted with crosses and the corre- as a function of the shift in the equilibrium angle positions Δϕ0 between the ground sponding parabolic fits around the local minima with solid lines. The insets show and excited states. MgPc+ and H2Pc+ ground and excited state are vertically offset the schematic models of the respective ground state cations in their equilibrium by increments of 1 meV for clarity. Note the zero Δϕ0 for H2Pc+, imposed by the positions. The angle ϕ is defined as between the molecule x-axes (crossing two symmetry of the system. Source data are provided as a Source Data file. opposing isoindole groups along N–N atom direction) and the nearest of the [110] Fig. 3 | Scheme of the electron-libron coupling model and fitting of the spectra. H2Pc+. Measurement parameters were −2.8 V, 50 pA in c; −3 V, 60 pA in d; −3 V, a Scheme of the energy curves of ground and excited states illustrating the libronic 100 pA in e. The Franck-Condon factors are calculated including the modified vertical downward transitions. Excited-state population probability is considered as exponential distribution and are colour-coded according to the librational quan- exponential distribution or modified exponential distribution with suppressed tum number difference between the initial and final state, i.e. m–n (red - negative, probability of m = 0 state. b Scheme of the lowest-energy libronic transitions blue - positive, grey - null). The energies E00 obtained through the fitting are set as coupled to the electronic D1 → D0 transition with colours matching the peaks in the the reference in each spectrum. Parameters of the simulated spectra are listed in panels c–e. c–e Comparison of experimental (black dots) and simulated (solid blue Table 1. Source data are provided as a Source Data file. line) STM-EL fine spectra of Q+ peak of ZnPc+ and MgPc+ (c, d) and Qy1+ and Qy2+ of Nature Communications | (2022)13:6008 4 Article https://doi.org/10.1038/s41467-022-33653-7 Table 1 | Parameters of the simulated spectra in Fig. 3c–e k0 (meV/(°)2) k1 (meV/(°)2) Δϕ0(°) Teff (K) A γ (meV) E00 (eV) ZnPc+ 4 ML 1.64 1.84 0.603 62 0.5 0.64 1.5186 MgPc+ 4 ML 1.41 1.69 0.698 75 0.69 0.95 1.5147 H2Pc+ 3 ML (Q+y1) 1.57* 1.63* 0.001 50* 1* 0.91 1.3938 Parameter values for H2Pc denoted by an asterisk (*) were fixed; k0 and k1 were taken from the TD-DFT simulations and A was set to unity. For fitting of the H2Pc spectrum minor component, three additional free parameters (not listed in the table) characterising the emission peak Qy2+ of the second tautomer were used: relative intensity I = 0.29, FWHM of dressing function γ2 = 0.31 meV and mutual energy separation of the tautomers ΔE00 = 1.17 meV. a b c Fig. 4 | Evolution of libronic fingerprints with the gap size. a Experimental (black factor A as a function of Δz obtained from the fitting procedure of data in a. Average dots) and simulated (solid blue lines) STM-EL fine spectra of Q+ peak of MgPc+ as a values of fitting parameters: k0 = (1.42 ± 0.01) meV/(°)2, k1 = (1.69 ± 0.04) meV/(°)2, function of increasing tip-sample separation by Δz. V = −3 V, I ranges from 113 pA Δϕ0 = (0.69 ± 0.03)°, Teff = (82 ± 10) K. Source data are provided as a Source (bottom spectrum) to 19 pA (top spectrum). Acquisition time was 120 s for Δz = 0, Data file. 0.5, 1 Å and 360 s for Δz = 1.5 and 2 Å. b, c The energy E00, linewidth γ and reducing Nanocavity tuning of the initial zero-level libron state vibration excited states to the vibrational ground state of the trion is a probability likely explanation for this phenomenon. In order to shed light on the mechanisms leading to the modulation of At this point, one may wonder if librational peak manifolds could the peak around zero-phonon line in the spectra of the asymmetrically be also observed for the neutral excitons of the molecules. Based on adsorbed molecules, we measure how it is affected by opening/closing the theoretical calculations on neutral molecules (see Supplementary the nanocavity. Since it is known40 that the effective lifetime of the Fig. 6 and Supplementary Table 2), revealing very similar values of Δϕ0, excitation in the phthalocyanines is reduced by the confinement of the k0 and k1 as the cations (Fig. 2), we believe the manifolds are also optical density of states (Purcell factor), we are expecting a modulation present in the Q peak of neutral ZnPc and MgPc spectra, although of w0′ due to a variable radiative quenching by the nanocavity. The indistinguishable due to a naturally broad energy character of the resulting dependence of the MgPc+/4ML-NaCl spectra on the tip- electroluminescence emission (see Supplementary Fig. 3). This might sample distance is presented in Fig. 4a, along with the parameters be overcome with a resonant STM-PL that already demonstrated its A(Δz), the E00(Δz) and the overall line broadening γ in Fig. 4b, c that we resolution capability on Qx of H2Pc/4 ML NaCl (0.5 meV linewidth)20. determined by fitting of each spectrum individually. At the first glance, To conclude, we found the link between the observation of a significant increase of the central peak intensity with decreasing tip- libronic signatures in the electroluminescence spectra and the chiral sample distance Δz is apparent; there is also an overall redshift of the adsorption geometry of chromophores on NaCl. The molecules with entire spectra (decreasing E00) and peak broadening (increasing γ). chirally asymmetric adsorption configurations (Zn-, MgPc+) change the The latter two can be attributed to the known Lamb/Stark shift orientation upon excitation, which according to the Franck-Condon resulting from the coupling of the excited states to the nanocavity35. principle allows the transitions between different libration states. This However, the striking variation of the central peak intensity, with fac- gives rise to the observed manifolds of peaks in the spectra of chirally tor A changing from 0.40 to 0.93, is indicating a possible suppression adsorbed chromophores. In contrast, in the reference system (H2Pc+) of the deexcitation channel due to the faster radiative rate induced by with a mirror-symmetrical adsorption geometry, the adsorption compressing the nanocavity. At the moment, the mechanism respon- orientation remains unchanged upon excitation and, therefore, libra- sible for the decreasing intensity of the peak near the zero-phonon line tional sidebands are not arising. From the analysis of the experimental in ZnPc+ and MgPc+ is not fully clear to us. Nonetheless, a competition spectra, it follows that the process of excitation gives rise to a non- between the radiative and nonradiative decay rates from the zero- equilibrium initial libration states population, corroborating one of the Nature Communications | (2022)13:6008 5 Article https://doi.org/10.1038/s41467-022-33653-7 previously suggested mechanisms20 for the spectral broadening in the wavefunctions and energies, Franck-Condon factors are calculated STM-EL. Moreover, changes in the potential well stiffness of the using the parameters Teff and A. The Franck-Condon factors are plot- libration, associated with the excitation, leads to an additional peak ted in Fig. 3 as vertical bars at their corresponding energies. The entire broadening. All these effects have to be considered for a correct spectrum is calculated by summing all Franck-Condon factors, each interpretation of the STM-EL spectra of molecules in neutral and broadened by convolution with a broadening function of FWHM = γ. charged excited states. Finally, we have found experimental indication Finally, the spectrum energy range is offset by E00 and the simulated of a possible depopulation pathway predominantly affecting the zero- spectrum is resampled onto the experimental datapoint energy range. libration state of the trion. It can be effectively suppressed by overall The normalised sum of absolute differences between the simulated increase of the radiative decay rate by closing the STM-EL nanocavity. and experimental intensity values, corresponding by energy (I sim i and We anticipate that the newly emerging methodology of STM-PL20, I exp i respectively), called R-factor is calculated as enhanced with pump-probe capability could provide insight into the P P exp dynamics and physical origin of such deexcitation mechanism. ∣I exp i  cI sim i ∣ ∣I i ∣ i i R= P exp where c = P sim Methods ∣I i ∣ ∣I i ∣ i i Sample preparation and STM measurements All measurements were performed in ultrahigh vacuum low tempera- The value of R describes the agreement (R = 0 − perfect match). ture (at 7.5 K) STM/AFM microscope with base pressure below 5 × 10−10 The iterative optimisation procedure first calculates the Rðv0 Þ for mbar. The Ag(111) clean surface was prepared by standard cycles of an initial guess of the parameters – vector v0. Subsequently the para- sputtering and annealing. NaCl was evaporated from a source at 607 °C meters are repeatedly adjusted by introducing random adjustments on the surface kept at 120 °C during 3-5 min to obtain 2-4 ML thick (Δv). For each Δv that generates improvement of the R-factor such that islands. Once the sample was inserted into the microscope head and Rðvj Þ < Rðvj1 Þ (index j being the number of iteration), the change in the cooled, the phthalocyanine molecules (from Sigma Aldrich) were vector is retained, i.e. vj = vj1 + Δv, otherwise vj = vj1 is used. After deposited on it from a Ta crucible at 331 °C (H2Pc), 380 °C (MgPc) and several thousand iterations (n) we plot the simulated spectrum for vn. 415 °C (ZnPc). We used a PtIr tip, sharpened by a focused-ion beam The initial parameters (v0) were chosen from the calculations (k0, before inserting it into the scanner. Tips were cleaned and coated by k1, Δϕ0) or from initial guesses based on the typology of spectra for Ag (or eventually Au) material by applying voltage pulses and con- different parameters shown in Supplementary Fig. 5. trolled nanoscopic indentations into the clean substrate in order to To avoid underdetermination of the fit of H2Pc+ and to account for achieve a suitable near-infrared plasmonic response. its tautomerization-related dual-peak structure, we fixed parameters k0, k1,Teff and A as described in the caption of Table 1 and fitted STML measurements simultaneously the Qy1+ and Qy2+ peaks, using independent center Photons were collected by a ZEONEX® aspherical lens with 12 mm dia- energies E00(Qy1+) and E00(Qy2+), broadening factors γ1 and γ2 and with meter mounted on the scanner head 16.5 mm from the tunnelling a common Δϕ0. junction and directed out of the cryostat through a set of viewports. The The z-axis moment of inertia of the molecules was calculated from outgoing beam was refocused into an optical fibre bunch leading to the ground state cation optimised structure (B3LYP/6-31G*). The Andor Kymera 328i spectrograph equipped with a CCD sensor (Andor resulting Jz values are 113.7 mpnm2 for ZnPc+, 114.3 mpnm2 for MgPc+ Newton DU920P-BEX2-DD). Overview and high-resolution photon and 113.2 mpnm2 for H2Pc+ where mp is the proton mass. spectra were obtained using gratings of 150 and 1200 grooves/mm and Spectrum of ZnPc+ molecule (Fig. 3c) was fitted in the energy 100 µm wide slit, which provide spectral resolutions of 1.2 and 0.2 nm range (E00−5 meV, E00+10 meV) and spectra of MgPc+ molecules respectively. The best achievable energy resolution was 300, 400 and (Figs. 3d and 4a) were fitted in the energy range (E00 − 4 meV, E00 600 µeV for H2Pc Qy+, Q+ and Q of Mg/ZnPc peaks respectively. +10 meV). The spectrum of H2Pc+ (Fig. 3e) was fitted in the range (E00 ± 4 meV). (Gaussian for Zn- and MgPc+, Lorentzian for H2Pc+). The DFT and TD-DFT calculations fitting code with an example file is provided in ref. 44. The single-molecule calculations were run on H2Pc, ZnPc and MgPc in both their neutral and cationic state. The ground state molecular Reporting summary structure was optimised in vacuo with density functional method Further information on research design is available in the Nature (DFT) and ωB97X-D/6-31 G* level of theory41. The emission properties Research Reporting Summary linked to this article. are obtained by optimisation of the first (and second for H2Pc) excited states with TD-ωB97X-D/6-31 G* for neutral and TDA-ωB97X-D/6-31G* Data availability for cations42. To calculate the total energy as a function of the mole- Source data are provided with this paper. Additional data that support cular adsorption orientation, the electrostatic field of NaCl surface has the findings of this study are available from the corresponding authors been modelled as a slab of 3 layers of 6 × 6 point charges (+1 for Na+ upon request. and −1 for Cl-) with fixed position, at 3 Å of distance with the molecular plane. The optimised structure in vacuo of the neutral or cation Code availability molecule has been employed (B3LYP/6-31G*) and all coordinates are The fitting code with an example file is provided in Ref. 44. kept fixed except for the azimuthal angle ϕ, varying between 0 and 45° (step of 1°) for ZnPc or MgPc, and between 0 and 90° (step of 5°) for References H2Pc. All calculations were performed with the Gaussian16 package43. 1. Cao, J. et al. Quantum biology revisited. Sci. Adv. 6, eaaz4888 (2020). Spectroscopic fitting procedure 2. Tiwari, V., Peters, W. K. & Jonas, D. M. Electronic resonance with We used an iterative fitting procedure (described below) to minimise anticorrelated pigment vibrations drives photosynthetic energy the sum of absolute differences between the simulated spectra transfer outside the adiabatic framework. Proc. Natl Acad. Sci. USA intensities and experimental datapoints, by optimisation of the free 110, 1203–1208 (2013). parameters of the model, i.e. the vector v = (k0, k1, E00, Δϕ0, Teff, A, γ). 3. Panitchayangkoon, G. et al. Direct evidence of quantum transport in First, the Schrödinger equation in Eq. 1 is solved for the excited and photosynthetic light-harvesting complexes. Proc. Natl Acad. Sci. ground state using the parameters k0, k1, Δϕ0. In the next step, from USA 108, 20908–20912 (2011). Nature Communications | (2022)13:6008 6 Article https://doi.org/10.1038/s41467-022-33653-7 4. Schnedermann, C. et al. Evidence for a vibrational phase- phthalocyanine molecule adsorbed on an ultrathin NaCl film. Phys. dependent isotope effect on the photochemistry of vision. Nat. Rev. B 93, 165419 (2016). Chem. 10, 449–455 (2018). 30. Doppagne, B. et al. Vibronic spectroscopy with submolecular 5. Hahn, S. & Stock, G. Quantum-mechanical modeling of the fem- resolution from STM-induced electroluminescence. Phys. Rev. Lett. tosecond isomerization in rhodopsin. J. Phys. Chem. B 104, 118, 127401 (2017). 1146–1149 (2000). 31. Doppagne, B. et al. Single-molecule tautomerization tracking 6. Kottas, G. S., Clarke, L. I., Horinek, D. & Michl, J. Artificial molecular through space- and time-resolved fluorescence spectroscopy. Nat. rotors. Chem. Rev. 105, 1281–1376 (2005). Nanotechnol. 15, 207–211 (2020). 7. Conyard, J. et al. Ultrafast dynamics in the power stroke of a 32. Baeten, Y., Fron, E., Ruzié, C., Geerts, Y. H. & Van Der Auweraer, M. molecular rotary motor. Nat. Chem. 4, 547–551 (2012). Investigation of the Qx–Qy equilibrium in a metal‐free phthalocya- 8. Coropceanu, V., Chen, X.-K., Wang, T., Zheng, Z. & Brédas, J.-L. nine. ChemPhysChem 16, 992–3996 (2015). Charge-transfer electronic states in organic solar cells. Nat. Rev. 33. Luo, Y. et al. Electrically driven single-photon superradiance from Mater. 4, 689–707 (2019). molecular chains in a plasmonic nanocavity. Phys. Rev. Lett. 122, 9. Torres, A., Oliboni, R. S. & Rego, L. G. C. Vibronic and coherent 233901 (2019). effects on interfacial electron transfer dynamics. J. Phys. Chem. 34. Zhang, L. et al. Electrically driven single-photon emission from an Lett. 6, 4927–4935 (2015). isolated single molecule. Nat. Commun. 8, 580 (2017). 10. Kuhnke, K., Große, C., Merino, P. & Kern, K. Atomic-scale imaging 35. Rosławska, A. et al. Mapping lamb, stark, and Purcell effects at a and spectroscopy of electroluminescence at molecular interfaces. chromophore-picocavity junction with hyper-resolved fluores- Chem. Rev. 117, 5174–5222 (2017). cence microscopy. Phys. Rev. X 12, 011012 (2022). 11. Qiu, X. H., Nazin, G. V. & Ho, W. Vibrationally resolved fluorescence 36. Paulheim, A. et al. Inhomogeneous and homogeneous line broad- excited with submolecular precision. Science 299, ening of optical spectra of PTCDA molecules adsorbed at step 542–546 (2003). edges of alkali Halide surfaces. J. Phys. Chem. C 120, 12. Imada, H. et al. Real-space investigation of energy transfer in het- 11926–11937 (2016). erogeneous molecular dimers. Nature 538, 364–367 (2016). 37. Marquardt, C., Paulheim, A., Hochheim, M., Bredow, T. & Soko- 13. Zhang, Y. et al. Visualizing coherent intermolecular dipole–dipole lowski, M. Homogeneous and inhomogeneous line shape of the coupling in real space. Nature 531, 623–627 (2016). electronic excitation of a single molecule on a surface. Phys. Rev. B 14. Doppagne, B. et al. Electrofluorochromism at the single-molecule 104, 045415 (2021). level. Science 361, 251–255 (2018). 38. Patera, L. L., Queck, F., Scheuerer, P., Moll, N. & Repp, J. Accessing a 15. Merino, P., Große, C., Rosławska, A., Kuhnke, K. & Kern, K. Exciton charged intermediate state involved in the excitation of single dynamics of C60-based single-photon emitters explored by Han- molecules. Phys. Rev. Lett. 123, 016001 (2019). bury Brown-Twiss scanning tunnelling microscopy. Nat. Commun. 39. Martín-Jiménez, A. et al. Electronic temperature and two-electron 6, 8461 (2015). processes in overbias plasmonic emission from tunnel junctions. 16. Doležal, J. et al. Charge carrier injection electroluminescence with Nano Lett. 21, 7086–7092 (2021). CO-functionalized tips on single molecular emitters. Nano Lett. 19, 40. Doležal, J., Canola, S., Merino, P. & Švec, M. Exciton-trion conver- 8605–8611 (2019). sion dynamics in a single molecule. ACS Nano 15, 7694–7699 17. Hung, T.-C., Kiraly, B., Strik, J. H., Khajetoorians, A. A. & Wegner, D. (2021). Plasmon-driven motion of an individual molecule. Nano Lett. 21, 41. Chai, J.-D. & Head-Gordon, M. Long-range corrected hybrid density 5006–5012 (2021). functionals with damped atom-atom dispersion corrections. Phys. 18. Rai, V. et al. Boosting light emission from single hydrogen phthalo- Chem. Chem. Phys. 10, 6615–6620 (2008). cyanine molecules by charging. Nano Lett. 20, 7600–7605 (2020). 42. Martynov, A. G. et al. Methodological survey of simplified TD-DFT 19. Yang, B. et al. Sub-nanometre resolution in single-molecule pho- methods for fast and accurate interpretation of UV-Vis-NIR spectra toluminescence imaging. Nat. Photon. 14, 693–699 (2020). of phthalocyanines. ACS Omega 4, 7265–7284 (2019). 20. Imada, H. et al. Single-molecule laser nanospectroscopy with 43. Frisch, M. J. et al. Gaussian 16, Revision A.03 (Gaussian, Inc., 2016). micro-electron volt energy resolution. Science 373, 95–98 (2021). 44. Doležal, J., Švec, M., & Hapala, P. Franck–Condon model. Zenodo. 21. Zhang, R. et al. Chemical mapping of a single molecule by plasmon- https://doi.org/10.5281/zenodo.6726424 (2022). enhanced Raman scattering. Nature 498, 82–86 (2013). 22. Lee, J., Crampton, K. T., Tallarida, N. & Apkarian, V. A. Visualizing Acknowledgements vibrational normal modes of a single molecule with atomically S.C., R.C.C.F., M.Š. and J.D. acknowledge the Czech grant agency confined light. Nature 568, 78–82 (2019). funding no. 20-18741 S and the Charles University Grant Agency pro- 23. Jaculbia, R. B. et al. Single-molecule resonance Raman effect in a ject no. 910120. P.M. acknowledges grants EUR2021-122006, plasmonic nanocavity. Nat. Nanotechnol. 15, 105–110 (2020). RYC2020-029800-I and PID2021-125309OA-I00 funded by MCIN/AEI/ 24. Schaffert, J. et al. Imaging the dynamics of individually adsorbed 10.13039/501100011033 and European Union NextGener- molecules. Nat. Mater. 12, 223–227 (2013). ationEU/PRTR. 25. Peller, D. et al. Sub-cycle atomic-scale forces coherently control a single-molecule switch. Nature 585, 58–62 (2020). Author contributions 26. Kü gel, J., Klein, L., Leisegang, M. & Bode, M. Analyzing and tuning J.D., P.M. and M.Š. have conceived the experiment. J.D., R.C.C.F. and the energetic landscape of H2Pc tautomerization. J. Phys. Chem. C M.Š. have performed the experiments and preprocessed the data for 121, 28204–28210 (2017). analysis. TD-DFT calculations and their analysis were performed by S.C. 27. Kimura, K. et al. Selective triplet exciton formation in a single The model used for fitting was created by P.H., refined by J.D. and M.Š. molecule. Nature 570, 210–213 (2019). and the fitting was made by J.D. All authors have thoroughly discussed 28. Doležal, J. et al. Real space visualization of entangled excitonic the data and contributed to the creation of the manuscript including states in charged molecular assemblies. ACS Nano 15, figures and tables. 7694–7699 (2021). 29. Miwa, K., Imada, H., Kawahara, S. & Kim, Y. Effects of molecule- Competing interests insulator interaction on geometric property of a single The authors declare no competing interests. Nature Communications | (2022)13:6008 7 Article https://doi.org/10.1038/s41467-022-33653-7 Additional information Open Access This article is licensed under a Creative Commons Supplementary information The online version contains Attribution 4.0 International License, which permits use, sharing, supplementary material available at adaptation, distribution and reproduction in any medium or format, as https://doi.org/10.1038/s41467-022-33653-7. long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if Correspondence and requests for materials should be addressed to Jiří changes were made. The images or other third party material in this Doležal or Martin Švec. article are included in the article’s Creative Commons license, unless indicated otherwise in a credit line to the material. If material is not Peer review information Nature Communications thanks Hiroshi Imada, included in the article’s Creative Commons license and your intended Zhen-Chao Dong and the other, anonymous, reviewer(s) for their use is not permitted by statutory regulation or exceeds the permitted contribution to the peer review of this work. Peer reviewer reports are use, you will need to obtain permission directly from the copyright available. holder. To view a copy of this license, visit http://creativecommons.org/ licenses/by/4.0/. Reprints and permission information is available at http://www.nature.com/reprints © The Author(s) 2022 Publisher’s note Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations. Nature Communications | (2022)13:6008 8