## Substrate stabilization of Jahn-Teller distortion in a single molecule

Moritz Frankerl, 1, 2, ∗ Laerte L. Patera, 3, 4 Thomas Frederiksen, 2, 5 Jascha Repp, 3, † and Andrea Donarini 1, ‡

- 1 Institute for Theoretical Physics, University of Regensburg, 93040 Regensburg, Germany 2 Donostia International Physics Center (DIPC), E-20018 Donostia-San Sebasti´n, Spain a 3 Institute for Experimental Physics, University of Regensburg, 93040 Regensburg, Germany 4 Department of Physical Chemistry, University of Innsbruck, 6020 Innsbruck, Austria 5 IKERBASQUE, Basque Foundation for Science, E-48013, Bilbao, Spain (Dated: August 2, 2024)

Charge-state transitions of a single Cu-phthalocyanine molecule adsorbed on an insulating layer of NaCl on Cu(111) are probed by means of alternate charging scanning tunneling microscopy. Real-space imaging of the electronic transitions reveals the Jahn-Teller distortion occurring upon formation of the first and second anionic charge states. The experimental findings are rationalized by a theoretical many-body model which highlights the crucial role played by the substrate. The latter enhances and stabilizes the intrinsic Jahn-Teller distortion of the negatively charged molecule hosting a degenerate pair of single-particle frontier orbitals. Consequently, two excess electrons are found to occupy, in the ground state, the same localized orbital, despite a larger Coulomb repulsion than the one for the competing delocalized electronic configuration. Control over the charging sequence by varying the applied bias voltage is also predicted.

Introduction Scanning tunneling microscopy (STM) is an excellent tool for investigating the electronic structure of single molecules on surfaces. By combining atomicscale imaging with spectroscopy of electronic states, STM enables a comprehensive understanding of molecularscale systems. Ultrathin insulating films of 2-3 monolayers (ML) have been used to separate the molecule from the metallic substrate and thereby inhibit their electronic hybridization, providing access to the (near-) pristine electronic structure of the molecule [1]. This enabled the observation of vibrational side bands [2], coherent electron-nuclear coupling [3], many-body correlations [4, 5], charge-state lifetimes [6], and coherent spin control [7] on a single molecule. Furthermore, this facilitated the demonstration of single-molecule fluorescence [8, 9] permitting vibronic spectroscopy with submolecular resolution [10], single-molecule electrofluorochromism [11], energy transfer in molecular dimers [12, 13], AntiKasha emissions [14], switching of the excitonic state of a single molecule [15], and bistability between π -diradical open- and closed-shell states [16].

By employing thicker insulating films ( &gt; 20 ML), electron exchange with the substrate can be completely inhibited, enabling the stabilization of multiple charge states in single-molecule investigations [17-19]. This approach requires atomic force microscopy (AFM) for the detection of electron exchange with the tip. It enabled probing and quantifying of single-electron transfer between molecules [20, 21], the reorganization energy upon single-molecule charging [22], excited molecular states [23], as well as their lifetime [24]. Further, based on this approach the imaging of charge-state transitions in single molecules with angstrom resolution [25-27] has been developed. This technique, dubbed single-electron alternate charging scanning tunneling microscopy (AC-STM), is based on driving a single electron to tunnel between the conductive AFM tip and a single molecule deposited on a non-conductive substrate. Charging events are controlled by short voltage pulses ( V ac ) of opposite polarity that are superimposed on the DC voltage ( V dc ) across the AFM junction, inducing alternating charging of the molecule. While several charge state transitions can be accessed by gating through V dc , the polarity of the voltage pulses of V ac provides selectivity in probing a specific direction of the charge-state transition [25].

Polar insulators strongly stabilize excess charges in adsorbed atoms and molecules [18, 19, 27-31], facilitating orbital localization through self-trapping of polarons [32, 33], revealing that the insulator does not merely electronically isolate the adsorbed molecule, but plays an active role in determining its electronic properties. Whereas the reorganization energy is a well studied manifestation of this effect, little is known about how the large polarizability of the environment influences other electronic properties, as, for example, the occupation sequence of the molecular orbitals.

Here we study the electronic states in a charged Cuphthalocyanine (CuPc) molecule as a model system that is known to exhibit a Jahn-Teller (JT) distortion. We employ AC-STM imaging to reveal the energetic sequence of orbital occupation for the JT distorted system. Injection of a second electron results in the occupation of the same formerly-degenerate orbital. Our theoretical analysis reveals that this sequence is not merely a molecular property but it is critically governed by the polarizable NaCl substrate, strongly enhancing the JT splitting. Moreover, we predict that the delocalized configuration is accessible at higher bias voltages where it dominates the electronic transition rate due to the state multiplicity.

Experiments We employ AC-STM [Fig. 1(a)], cf. SI for details, to image the spatial distribution of charge-state

FIG. 1. (a) Sketch of the experimental setup. On thick NaCl films ( &gt; 20 ML), tunneling from the CuPc molecule to the Cu substrate is inhibited, enabling charge-state manipulation and the visualization of electronic transitions via AC-STM. (b) Isosurfaces of the two degenerate LUMOs present in CuPc. (c-f) AC-STM images (oscillation amplitude A = 1 ˚) correA sponding to: (c) 0 → 1 -transition ( V dc = 1 0 V, . V ac = 0 75 . V peak-to-peak (V pp ), ∆ z = 3 7 ˚) (d) 1 . A -→ 0 transition ( V dc = 1 0 V, . V ac = 1 0 V . pp , ∆ = 4 4 ˚) (e) 1 z . A -→ 2 -transition ( V dc = 2 3 V, . V ac = 0 9 V . pp , ∆ = 4 8 ˚) (f) 2 z . A -→ 1 -transition ( V dc = 2 3 V, . V ac = 1 0 V . pp , ∆ z = 4 0 ˚ ) ∆ . A z are given with respect to an AFM set point of ∆ f = -1 5 Hz at . V dc = 0 V. Scale bar: 5 ˚. A

![image](https://ac.marcodsn.me/arxiv-markdown-images/2408.00478v1-figure-c3c7defd-4be2-46ba-ae66-d149aac70b04.jpg)

transitions for the neutral, anionic and dianionic CuPc species. The map of the 0 → 1 -transition [Fig. 1(c)] exhibits a 4-fold rotational symmetry, similar to the anionic resonance observed by STM on ultrathin NaCl layers (2 ML) [30]. Toggling of the molecule yields, on average, the superposition of two images rotated of about 32 degrees with respect to each other [26, 34]). As reported previously [26], for the reverse transition (1 -→ 0)[Fig. 1(d)] the symmetry reduces from 4- to 2-fold due to the JT distortion, which favors the occupation of one single LUMO [30]. In the dianionic ground state the second excess electron can either occupy the other (empty) LUMO or further populate the already occupied one. The actual preference results from a competition between JT distortion and Coulomb repulsion, the latter being stronger for two electrons occupying the same orbital. The AC-STM images of the 1 -→ 2 -and 2 -→ 1 -transitions, displayed respectively in Fig. 1(e) and (f), reveal that the second excess electron is occupying the same orbital as the first one.

Model We model the molecule in the AC-STM with the Hamiltonian

$$\hat { H } = \hat { H } _ { \text{mol} } + \hat { H } _ { \text{tip} } + \hat { H } _ { \text{tun} } + \hat { H } _ { \text{as} } + \hat { H } _ { \text{vib} } + \hat { H } _ { \text{e-ph} } \,, \ \ ( 1 )$$

where the minimal electronic description of the molecule (mol), metallic tip (tip) and tunnelling coupling (tun) is complemented by a small environmental asymmetry (as), the set of molecular and substrate vibrational modes (vib), and by their coupling to the electronic degrees of freedom (e-ph). We write the molecular Hamiltonian in the Fock space of the four frontier orbitals, which are the singly occupied molecular orbital (SOMO, S ), the highest occupied molecular orbital (HOMO, H ) and the two degenerate lowest unoccupied molecular orbitals (LUMOs, L xz and L yz ) [35-37]. The latter two are shown in Fig. 1(b). The single-particle energies, ε S = -12 0 eV, . ε H = -11 7 eV, . ε L xz /L yz = -10 7 eV, . are renormalized by a crystal field correction of ∆ = 2 44eV and . all symmetry-allowed interaction matrix elements V ijkl of the Coulomb interaction are retained. The tip is described as a non-interacting Fermi sea with harmonic lateral confinement [38, 39] ˆ H tip = ∑ kσ ( ε k + cV b ) ˆ c † kσ ˆ c kσ , where the tip states ε k are shifted by the applied bias V b subject to a lever arm c caused by the voltage drop in the insulating layer, ˆ c kσ destroys an electron in the tip with momentum k and spin σ . The tunneling is expressed as ˆ H tun = ∑ ikσ t tip ikσ ˆ c † kσ ˆ d iσ +Hc , only consider-. . ing tunneling processes between tip and molecule since molecule-substrate tunneling is inhibited. The tunneling amplitudes t tip ikσ are modelled to be proportional to the overlap between the molecule and the tip wave functions and thereby contain the dependence on tip position. The tip wave functions are composed by a Gaussian component stemming from the lateral confinement and an evanescent part along the direction perpendicular to the molecular plane [38, 39]. We incorporate an environmental asymmetry in our model via ˆ H as = -δ as (ˆ n xz -ˆ n yz ), with δ as = 4 k B T = 2meV.

The vibronic Hamiltonian is a collection of harmonic oscillators ˆ H vib = ∑ ζ ¯ hω ζ (ˆ a a † ζ ˆ ζ + 1 2 / ) containing molecular and salt modes. Based on symmetry considerations, the vibrational modes of the system fall into three categories, namely, (i) those that do not interact with charges in the molecule, (ii) those that affect both LUMOs the same way, and (iii) those that interact differently with each of the LUMOs. The C 4 v symmetry of the combined salt and molecular system facilitates this distinction of vibronic modes. We do not consider the modes falling into the first category. The modes constituting the second category lead to a stabilization of charge states in the molecule. These are the ones transforming as the A 1 irreducible representation and we call them symmetric (S) modes from now on since they do not reduce the symmetry of the system. They are only sensitive to the

FIG. 2. (a), (b) Examples of a A 1 and a B 1 molecular vibrational mode, respectively (Cu in brown, C in grey, N in blue, H in white). (c), (d) In-plane displacement of the salt atoms (Cl -in green, Na + in purple) corresponding to the A 1 ( B 1 ) mode, respectively.

![image](https://ac.marcodsn.me/arxiv-markdown-images/2408.00478v1-figure-6dd7ff43-ab8e-4b44-8120-2be9dddebd59.jpg)

global charge on the LUMOs, thus

$$\hat { H } _ { e - \text{ph} } ^ { \text{S} } = \sum _ { \alpha \in \{ A _ { 1 } \} } g _ { \alpha } ( \hat { n } _ { x z } + \hat { n } _ { y z } ) \left ( \hat { a } _ { \alpha } ^ { \dagger } + \hat { a } _ { \alpha } \right ) \,. \quad ( 2 )$$

The modes in the third category affect the charge distribution within a specific charge state. We refer to them as antisymmetric (AS) modes since they lower the symmetry of the molecule. They belong to the B 1 and B 2 irreducible representation. However, due to energetic considerations (cf. SI) we completely neglect the B 2 modes. The antisymmetric ( B 1 ) modes couple to the imbalance between the populations of the L xz and L yz orbitals, i.e.

$$\hat { H } _ { e - \text{ph} } ^ { \text{AS} } = \sum _ { \beta \in \{ B _ { 1 } \} } g _ { \beta } ( \hat { n } _ { x z } - \hat { n } _ { y z } ) \left ( \hat { a } _ { \beta } ^ { \dagger } + \hat { a } _ { \beta } \right ) \,. \quad ( 3 ) \quad \underset { \text{mod}, \text{izati} } { \text{mod} }$$

We calculated the molecular mode energies as well as the corresponding electron-phonon couplings g α/β from first principles using Inelastica [40, 41], thus obtaining 14 symmetric ( A 1 ) and 14 antisymmetric ( B 1 ) modes. The coupled salt modes can be combined to obtain just a single symmetric ( A 1 ) and a single antisymmetric ( B 1 ) mode. We take the corresponding electron-phonon couplings as free parameters in our theory and set them to g salt A 1 = 68meV and g salt B 1 = 40meV, assuming for the energy the one of the transversal optical phonon mode ¯ hω TO = 20meV [42]. A microscopic derivation of this minimal model yielding couplings constants of comparable strength is given in the SI. We show a comparison between one symmetric molecular mode and the symmetric salt mode in Fig. 2(a) and (c). The antisymmetric case is displayed in Fig. 2(b) and (d).

The influence of the symmetric modes on the manybody spectrum is readily obtained by a canonical LangFirsov transformation, leading to a renormalization of the single-particle energies ˜ ϵ i = ϵ i -∑ α g 2 α / ¯ hω α as well as of the direct Coulomb interaction terms ˜ V ijji = V ijji -2 ∑ α g 2 α / ¯ hω α and i, j ∈ { L xz , L yz } . Thus, the energy renormalization induced by the symmetric salt mode explains the stabilization of charged molecules and atoms adsorbed on polar dielectrics due to the Franck-Condon blockade [43, 44]. For antisymmetric modes, the exchange and pair hopping between LUMOs shows an intricate interplay with the electron-phonon coupling. In this case, the Lang-Firsov transformation does not eliminate the coupling between the electronic and the mechanical degrees of freedom, which persists in the transformed exchange and pair-hopping terms. We identify, instead, a single JT active mode and resort for our analysis into a semi-classical adiabatic approximation. To this end, we introduce the canonical displacement operators

$$\hat { Q } _ { m } = \sum _ { \beta \in B _ { 1 } } \sqrt { \hbar { \omega } _ { \beta } } A _ { m \beta } \left ( \hat { a } _ { \beta } ^ { \dagger } + \hat { a } _ { \beta } \right ), \quad \ \ ( 4 )$$

where the first row of the transformation matrix reads A 1 β = 1 / g √ g 2 β / ¯ hω β and yields the JT active mode ˆ Q AS := ˆ Q 1 . The JT distortion is thus described by the operator

$$\hat { H } _ { J \text{T} } ( Q _ { \text{AS} } ) = \frac { Q _ { \text{AS} } ^ { 2 } } { 4 } + g \, ( \hat { n } _ { x z } - \hat { n } _ { y z } ) Q _ { \text{AS} } \quad ( 5 )$$

being Q AS the classical displacement of the JT active mode, with the effective coupling g = √ ∑ β g 2 β / ¯ hω β . On this basis, we simulate the experimentally investigated many-body eigenstates and energies by diagonalizing the effective Hamiltonian

$$\hat { H } _ { \text{eff} } = \hat { \tilde { H } } _ { \text{mol} } + \hat { H } _ { \text{as} } + \hat { H } _ { \text{JT} } ( Q _ { \text{AS} } ) \,, \quad \quad ( 6 )$$

along the coordinate Q AS of the JT-active antisymmetric mode, with ˆ ˜ H mol containing the Lang-Firsov renormalization due to the symmetric modes.

Electronic structure Some of the eigenenergies of ˆ H eff are the adiabatic potential energy surfaces (APES) plotted in Fig. 3. The lowest parabola in Fig. 3(d) is the APES of the neutral state. This neutral ground state defines the reference configuration ( Q AS = 0). The anionic states with an excess electron in one of the LUMO levels are orbitally degenerate for Q AS = 0. Their energy is lowered upon deforming along Q AS . Thus, the anionic ground state has one extra electron occupying the L xz ( L yz ) orbital for Q AS &lt; 0 ( Q AS &gt; 0). The small environmental asymmetry δ as determines here that the global minimum has the excess electron in L xz . The singlettriplet splitting due to SOMO-LUMO exchange is not resolved on this energy scale. So far, the energy spectrum

FIG. 3. Many-body states, energies and transitions. (a) Many-body states relevant for the transitions between neutral and anionic configurations. (b) 0 → 1 -at V b = 1 375 V, shows 4-fold symmetry. . All scale bars are 5 ˚ A and we use a Gaussian broadening of 4 ˚ . A (c) 1 -→ 0 at V b = 0 5 V, 4-fold symmetry is broken. . (d) Adiabatic potential energy surfaces (APES) as a function of the antisymmetric deformation coordinate ( Q AS ). The arrows represent the transitions to the corresponding ground states at the energy given by APES. (e) The 0 → 1 -transition involves a strong excitation of the symmetric modes ( Q S deformation coordinate) due to the large Huang-Rhys factor. The latter strongly suppresses the direct transition between ground states (black arrow). The blue arrow shows the transition from ground state to the vibronically excited state. The wavy downwards arrow indicates the relaxation dynamics towards the vibrational ground state. (f) Many-body states relevant for the transitions 1 -↔ 2 -. At a higher bias also the L yz orbital is populated, as indicated by the red arrow. (g) 1 -→ 2 -at V b = 2 75 V. . At the onset bias the second added electron occupies the same L xz orbital. The crosses belong to the positions reported in Fig. 4 (a). (h) 2 -→ 1 -at V b = 1 8 V. .

![image](https://ac.marcodsn.me/arxiv-markdown-images/2408.00478v1-figure-71463eab-8cc0-42da-b981-28e7115bdcf1.jpg)

̸

reflects the well-known JT effect. The richer structure of the dianionic APESs corresponds to two qualitatively different distributions of the excess electrons. For one excess electron in each of the L xz and L yz orbitals, the coupling to the JT mode vanishes. Thus, the lowest energy of these states (with several different spin multiplicities) is found at Q AS = 0 in the upper part of Fig. 3(d). On the other hand, states with both excess electrons in the same former LUMO strongly couple to the JT-mode. Their ground state is found at Q AS = 0, analogously to the anionic configuration, but with an even larger JT deformation and reorganization energy. Pair hopping between the real-valued orbitals induces the avoided crossing at Q AS = 0, characteristic of the pseudo JT effect [45].

Many-body transition rates A qualitative understanding of the measured charge-transitions is obtained by combining the electronic structure with few other constraints imposed by the experimental set-up. Specifically, the slow tip oscillations synchronized with the bias pulses allow for the system thermalization between consecutive tunnelling events and the large read-out tip-molecule distance yields single electron tunnelling. Moreover, the spectral broadening induced by the symmetric modes is much larger than the environmental asymmetry, which in turns exceeds the thermal broadening. Those premises imply, for example, that the neutral to anionic transition is 4-fold symmetric; it resembles the (incoherent) superposition of the L xz and L yz orbitals (with additional tog- gling), as it involves the spectrally unresolved states associated to both local minima with Q AS &gt; 0 and Q AS &lt; 0 [Fig. 1(c) and Fig. 3(b)]. In contrast, environmental asymmetry ensures, for the thermalized anion, the predominant occupation of the state with Q AS &lt; 0; consequently the anionic to neutral transition acquires the 2-fold rotational symmetry of the L xz orbital [Fig. 1(d) and Fig. 3(c)]. Mere energetic arguments would predict also for the transition between the anion and the dianion the participation of the states corresponding to both the deformed local minima in the dianionic APES. Single electron tunnelling, though, forbids the contribution of the state with Q AS &gt; 0 and double occupation of the L yz orbital, and ultimately yields a many-body transition resembling again L xz [Fig. 1(e) and Fig. 3(g)]. The same appearance is expected for the transition from the dianionic to the anionic ground state [Fig. 3(h)] although quantum fluctuations in the configuration of the dianion do not completely exclude small contributions from the L yz orbital, cf. Fig. 1(f) and the SI for the expression of the dianionic ground state as a function of the deformation coordinate.

Electron tunnelling implies displacements of both the symmetric and antisymmetric modes. Their interplay is visualized in Fig. 3(e). The large electron-phonon coupling to the symmetric modes strongly suppresses the transitions between ground states along the effective symmetric coordinate Q S [black arrow in Fig. 3(e)]. On

the contrary, w.r.t. the antisymmetric modes transitions between ground states are possible. This phenomenon, visualized for the neutral to anionic transition holds for all other considered transition between states with different charges.

The theoretical transition maps in Fig. 3 and Fig. 4 are calculated with Fermi's golden rule. The transition rate from an initial thermalized state towards a vibrationally excited state with neighbouring charge reads

$$R ^ { N \rightarrow N + 1 } & = \sum _ { \{ n _ { a } \}, f, i, j, \sigma } \prod _ { \alpha \in \{ A _ { 1 } \} } \Gamma _ { i \sigma, j \sigma } ( r _ { \text{tip} } ) P _ { \lambda _ { \alpha } } ( \{ n _ { a } \} ) f ^ { + } ( E _ { f } - E _ { g } - e e V _ { b } ) \langle N + 1, f | \hat { d } _ { i \sigma } ^ { \dagger } | N, g ) \langle N, g | \hat { d } _ { j \sigma } | N + 1, f ) \,, \, ( \text{f} ) \\ R ^ { N + 1 \rightarrow N } & = \sum _ { \{ n _ { a } \}, f, i, j, \sigma } \prod _ { \alpha \in \{ A _ { 1 } \} } \Gamma _ { i \sigma, j \sigma } ( r _ { \text{tip} } ) P _ { \lambda _ { \alpha } } ( \{ n _ { a } \} ) f ^ { - } ( E _ { g } - E _ { f } - e e V _ { b } ) \langle N + 1, g | \hat { d } _ { i \sigma } ^ { \dagger } | N, f ) \langle N, f | \hat { d } _ { j \sigma } | N + 1, g ) \,, \, ( \text{f} )$$

where the single-particle rate matrix is given by Γ iσ,jσ ( r tip ) = 2 π h / ¯ ∑ k t ∗ k,iσ t k,jσ δ E ( f -E g -ϵ k ). The symmetric modes' Franck-Condon factors yield the Poisson distributions P λ α ( { n α } ), which depend on the HuangRhys factors λ α = g α / ¯ hω α and { n α } is the set of the final state's vibrational excitations w.r.t. all modes α ∈ { A 1 } . Due to the large electron-phonon coupling λ salt ≈ 5, the salt contribution converges to a Gaussian and it is the main source of spectral broadening [46]. The Fermi function f + (and related f -= 1 -f + ) with the tip temperature T = 6K and chemical potential µ = -4 65 eV is . calculated at the difference between the molecular manybody energies shifted by the fraction c = 0 6 of sample . bias V b dropping between the molecule and the tip. The energy of the final state E f = E e f + ∑ α ¯ hω n α α includes the electronic potential energy E e f (the local minima of the APES in Fig. 3) and the excitation energies of the symmetric modes, while E g is the ground state energy of the initial state. Finally, Eqs. (7a) and (7b) includes the transition matrix elements between different charge states. To this end, we denote the initial state with | N,g ⟩ ( | N +1 , g ⟩ ) and the final state with | N +1 , f ⟩ ( | N,f ⟩ ). The sum over f extends to all possible states accessible from the designated initial state via a single electron tunnelling event. The tunneling rates in Eqs. (7a) and (7b) are depicted in Fig. 3 as a function of the tip position and show a strong similarity to the experimental results, with also a remarkable agreement between all the experimental and theoretical bias voltages for the different transitions.

We notice, in this respect, that the 1 -to 2 -transition denoted with a green arrow in Fig. 3 is calculated, in accordance to the experiment, at the onset of the anion to dianion transition ( V b = 2 75 V). . At higher biases we expect to open also transitions to the excited dianionic states, as sketched by the red arrows in Fig. 3(d) and (f). The bias dependence of the transition rate is shown in Fig. 4. We evaluate Eq. (7a) in the tip positions visualized by the black and teal crosses in the insets. They select tunnelling events involving a single orbital, respec-

FIG. 4. Bias dependence of the anionic to dianionic transition rate calculated in two different positions. The rates are given as a fraction of the maximum calculated value R max . The ratio between the rates varies with the bias, due to the opening of a further transition at higher bias [Fig. 2(f)]. The insets depict maps of the transition rate. The crosses mark the r 1 and r 2 positions used in the main figure. Star: Calculated at V b = 2 9 V. . Triangle: Map of the transition rate calculated at V b = 3 5 V. At this bias the transition involving . the L yz -orbital is dominant [Fig. 2(g)]. Scale bar is 5 ˚ . A

![image](https://ac.marcodsn.me/arxiv-markdown-images/2408.00478v1-figure-d8c349ae-28c7-49df-b395-69c450670092.jpg)

tively L xz for r 1 and L yz for r 2 . The form of the dianionic ground state determines the lower bias onset for the rate in position r 1 . At a higher bias voltages both rates saturate, with the one at r 2 being twice as large as the one at r 1 : the unoccupied orbital offers two spin channels as compered to the single one for an orbital already occupied by one electron. The predicted change of the full rate map is shown in the insets of Fig. 4. The comparison to Fig. 3(g) and Fig. 1 shows the emergence of the orthogonal LUMO ( L yz in our model) with an overall stronger signal.

Control over the landscape of the dianionic APES should also be obtained by changing the underlying salt.

It is experimentally known that RbI, NaCl, and LiF have increasing influence onto the spectral broadening of molecular electronic states [47, 48] and thus on the reorganization energy of the symmetric modes. We hypothesize that a similar trend also characterizes the antisymmetric modes, with a possible reversal of the orbital occupation sequence of CuPc on RbI with respect to the one on NaCl and LiF.

Conclusions We have shown how the electron-phonon coupling between a molecule and the underlying salt substrate determines the specific charge configuration within a molecule going beyond mere charge stabilization. We extend to this respect existing observations of charge localization in oligophenylene [27] and pentacene [31] to systems with JT distortion. The cooperative mechanism between salt substrate and molecule should generally appear in molecules and substrates sharing the same point group. This approach could pave the way to control charge configurations in molecules depending on the used substrate, with potential relevance in the study of intermolecular electron transfer, surface chemical reactivity and charge sensing. We also envisage the interplay between JT distortion of neighbouring molecules for the realization of cellular automata based on molecular arrays in which the local information is stored in the charge configuration of a single molecule.

Acknowledgement M.F., L.L.P., J.R. and A.D. acknowledge financial support from the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) within Project-ID 314695032 - SFB 1277 (project B02). J.R. gratefully acknowledges funding from the ERC Synergy Grant MolDAM (no. 951519) and the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) through RE2669/6-2. T.F. and M.F. acknowledge support from the Spanish MCIN/AEI/ 10.13039/501100011033 (PID2020-115406GB-I00), the Basque Department of Education (PIBA-2023-1-0021), and the Euskampus Transnational Common Laboratory QuantumChemPhys .

- ∗ moritz.frankerl@dipc.org
- † jascha.repp@ur.de
- ‡ andrea.donarini@ur.de
- [1] J. Repp, G. Meyer, S. M. Stojkovi´ c, A. Gourdon, and C. Joachim, Molecules on insulating films: Scanningtunneling microscopy imaging of individual molecular orbitals, Phys. Rev. Lett. 94 , 026803 (2005).
- [2] S. W. Wu, G. V. Nazin, X. Chen, X. H. Qiu, and W. Ho, Control of Relative Tunneling Rates in Single Molecule Bipolar Electron Transport, Phys. Rev. Lett. 93 , 236802 (2004).
- [3] J. Repp, P. Liljeroth, and G. Meyer, Coherent electronnuclear coupling in oligothiophene molecular wires, Nat. Phys. 6 , 975 (2010).
- [4] F. Schulz, M. Ij¨s, R. Drost, S. K. H¨m¨l¨inen, A. Harju, a a a a
- A. P. Seitsonen, and P. Liljeroth, Many-body transitions in a single molecule visualized by scanning tunnelling microscopy, Nat. Phys. 11 , 229 (2015).
- [5] P. Yu, N. Koci´ c, J. Repp, B. Siegert, and A. Donarini, Apparent reversal of molecular orbitals reveals entanglement, Phys. Rev. Lett. 119 , 056801 (2017).
- [6] K. Kaiser, L.-A. Lieske, J. Repp, and L. Gross, Chargestate lifetimes of single molecules on few monolayers of NaCl, Nat. Commun. 14 , 4988 (2023).
- [7] P. Willke, T. Bilgeri, X. Zhang, Y. Wang, C. Wolf, H. Aubin, A. Heinrich, and T. Choi, Coherent spin control of single molecules on a surface, ACS Nano 15 , 17959 (2021).
- [8] X. H. Qiu, G. V. Nazin, and W. Ho, Vibrationally resolved fluorescence excited with submolecular precision, Science 299 , 542 (2003).
- [9] L. Zhang, Y.-J. Yu, L.-G. Chen, Y. Luo, B. Yang, F.-F. Kong, G. Chen, Y. Zhang, Q. Zhang, Y. Luo, J.-L. Yang, Z.-C. Dong, and J. G. Hou, Electrically driven singlephoton emission from an isolated single molecule, Nat. Commun. 8 , 580 (2017).
- [10] B. Doppagne, M. C. Chong, E. Lorchat, S. Berciaud, M. Romeo, H. Bulou, A. Boeglin, F. Scheurer, and G. Schull, Vibronic spectroscopy with submolecular resolution from STM-induced electroluminescence, Phys. Rev. Lett. 118 , 127401 (2017).
- [11] B. Doppagne, M. C. Chong, H. Bulou, A. Boeglin, F. Scheurer, and G. Schull, Electrofluorochromism at the single-molecule level, Science 361 , 251 (2018).
- [12] H. Imada, K. Miwa, M. Imai-Imada, S. Kawahara, K. Kimura, and Y. Kim, Real-space investigation of energy transfer in heterogeneous molecular dimers, Nature 538 , 364 (2016).
- [13] S. Cao, A. Rosglyph[suppress]awska, l B. Doppagne, M. Romeo, M. F´ eron, F. Ch´ erioux, H. Bulou, F. Scheurer, and G. Schull, Energy funnelling within multichromophore architectures monitored with subnanometre resolution, Nat. Chem. 13 , 766 (2021).
- [14] H. Imada, M. Imai-Imada, X. Ouyang, A. Muranaka, and Y. Kim, Anti-Kasha emissions of single molecules in a plasmonic nanocavity, J. Chem. Phys. 157 , 104302 (2022).
- [15] J. Doleˇ zal, P. Mutombo, D. Nachtigallov´, a P. Jel´ ınek, P. Merino, and M. Svec, Mechano-optical switching of ˇ a single molecule with doublet emission, ACS Nano 14 , 8931 (2020).
- [16] S. Mishra, M. Vilas-Varela, L.-A. Lieske, R. Ortiz, S. Fatayer, I. Ronˇevi´, F. Albrecht, T. Frederiksen, D. Pe˜a, c c n and L. Gross, Bistability between π -diradical openshell and closed-shell states in indeno[1,2-a]fluorene, Nat. Chem. 16 , 755 (2024).
- [17] E. Bussmann and C. C. Williams, Single-electron tunneling force spectroscopy of an individual electronic state in a nonconducting surface, Appl. Phys. Lett. 88 , 263108 (2006).
- [18] L. Gross, F. Mohn, P. Liljeroth, J. Repp, F. J. Giessibl, and G. Meyer, Measuring the charge state of an adatom with noncontact atomic force microscopy, Science 324 , 1428 (2009).
- [19] W. Steurer, S. Fatayer, L. Gross, and G. Meyer, Probebased measurement of lateral single-electron transfer between individual molecules, Nat. Commun. 6 , 8353 (2015).
- [20] W. Steurer, J. Repp, L. Gross, I. Scivetti, M. Persson,

- and G. Meyer, Manipulation of the charge state of single au atoms on insulating multilayer films, Phys. Rev. Lett. 114 , 036801 (2015).
- [21] P. Scheuerer, L. L. Patera, and J. Repp, Manipulating and probing the distribution of excess electrons in an electrically isolated self-assembled molecular structure, Nano Lett. 20 , 1839 (2020).
- [22] S. Fatayer, B. Schuler, W. Steurer, I. Scivetti, J. Repp, L. Gross, M. Persson, and G. Meyer, Reorganization energy upon charging a single molecule on an insulator measured by atomic force microscopy, Nat. Nanotechnol. 13 , 376 (2018).
- [23] S. Fatayer, F. Albrecht, I. Tavernelli, M. Persson, N. Moll, and L. Gross, Probing molecular excited states by atomic force microscopy, Phys. Rev. Lett. 126 , 176801 (2021).
- [24] J. Peng, S. Sokolov, D. Hernang´ omez-P´rez, e F. Evers, L. Gross, J. M. Lupton, and J. Repp, Atomically resolved single-molecule triplet quenching, Science 373 , 452 (2021).
- [25] L. L. Patera, F. Queck, P. Scheuerer, and J. Repp, Mapping orbital changes upon electron transfer with tunnelling microscopy on insulators, Nature 566 , 245 (2019).
- [26] L. L. Patera, F. Queck, P. Scheuerer, N. Moll, and J. Repp, Accessing a charged intermediate state involved in the excitation of single molecules, Phys. Rev. Lett. 123 , 016001 (2019).
- [27] L. L. Patera, F. Queck, and J. Repp, Imaging charge localization in a conjugated oligophenylene, Phys. Rev. Lett. 125 , 176803 (2020).
- [28] J. Repp, G. Meyer, F. E. Olsson, and M. Persson, Controlling the charge state of individual gold adatoms, Science 305 , 493 (2004).
- [29] S. W. Wu, N. Ogawa, G. V. Nazin, and W. Ho, Conductance hysteresis and switching in a single-molecule junction, J. Phys. Chem. C 112 , 5241 (2008).
- [30] C. Uhlmann, I. Swart, and J. Repp, Controlling the orbital sequence in individual Cu-phthalocyanine molecules, Nano Lett. 13 , 777 (2013).
- [31] D. Hernang´ omez-P´rez, J. Schl¨r, D. A. Egger, L. L. Pae o tera, J. Repp, and F. Evers, Reorganization energy and polaronic effects of pentacene on NaCl films, Phys. Rev. B 102 , 115419 (2020).
- [32] B. Li, J. Zhao, K. Onda, K. D. Jordan, J. Yang, and H. Petek, Ultrafast interfacial proton-coupled electron transfer, Science 311 , 1436 (2006).
- [33] M. Setvin, C. Franchini, X. Hao, M. Schmid, A. Janotti, M. Kaltak, C. G. Van de Walle, G. Kresse, and U. Diebold, Direct view at excess electrons in TiO 2 rutile and anatase, Phys. Rev. Lett. 113 , 086402 (2014).
- molecule-insulator interaction on geometric property of a single phthalocyanine molecule adsorbed on an ultrathin NaCl film, Phys. Rev. B 93 , 165419 (2016).
- [35] B. Siegert, A. Donarini, and M. Grifoni, Effects of spinorbit coupling and many-body correlations in STM transport through copper phthalocyanine, Beilstein J. Nanotechnol. 6 , 2452 (2015).
- [36] B. Siegert, A. Donarini, and M. Grifoni, Nonequilibrium spin crossover in copper phthalocyanine, Phys. Rev. B 93 , 121406 (2016).
- [37] M. Frankerl and A. Donarini, Spin-orbit interaction induces charge beatings in a lightwave-STM -single molecule junction, Phys. Rev. B 103 , 085420 (2021).
- [38] S. Sobczyk, A. Donarini, and M. Grifoni, Theory of STM junctions for π -conjugated molecules on thin insulating films, Phys. Rev. B 85 , 205408 (2012).
- [39] A. Donarini, B. Siegert, S. Sobczyk, and M. Grifoni, Topographical fingerprints of many-body interference in STM junctions on thin insulating films, Phys. Rev. B 86 , 155451 (2012).
- [40] J. M. Soler, J. M. Artacho, J. D. Gale, A. Garc´ ıa, J. Junquera, P. Ordej´n, and D. S´nchez-Portal, The SIESTA o a method for ab initio orderN materials simulation, J. Phys.: Condens. Matter 14 , 2745 (2002).
- [41] T. Frederiksen, M. Paulsson, M. Brandbyge, and A.-P. Jauho, Inelastic transport theory from first principles: Methodology and application to nanoscale devices, Phys. Rev. B 75 , 205413 (2007).
- [42] G. Raunio and S. Rolandson, Lattice dynamics of NaCl, KCl, RbCl, and RbF, Phys. Rev. B 2 , 2098 (1970).
- [43] A. Mitra, I. Aleiner, and A. J. Millis, Phonon effects in molecular transistors: Quantal and classical treatment, Phys. Rev. B 69 , 245302 (2004).
- [44] J. Koch, F. von Oppen, and A. V. Andreev, Theory of the Franck-Condon blockade regime, Phys. Rev. B 74 , 205438 (2006).
- [45] I. B. Bersuker, The Jahn-Teller Effect (Cambridge University Press, 2006).
- [46] J. Repp, G. Meyer, S. Paavilainen, F. E. Olsson, and M. Persson, Scanning tunneling spectroscopy of Cl vacancies in NaCl films: Strong electron-phonon coupling in double-barrier tunneling junctions, Phys. Rev. Lett. 95 , 225503 (2005).
- [47] J. Repp and G. Meyer, Scanning tunneling spectroscopy of molecules on insulating films, CHIMIA 64 , 370 (2010).
- [48] N. Pavliˇek, c I. Swart, J. Niedenf¨hr, u G. Meyer, and J. Repp, Symmetry dependence of vibration-assisted tunneling, Phys. Rev. Lett. 110 , 136101 (2013).
- [34] K. Miwa, H. Imada, S. Kawahara, and Y. Kim, Effects of
- [49] F. J. Giessibl, Atomic resolution on Si(111)-(7 × 7) by noncontact atomic force microscopy with a force sensor based on a quartz tuning fork, Applied Physics Letters 76 , 1470 (2000).

## Experimental technique

Experiments were carried out with a low-temperature scanning tunneling and atomic force microscope equipped with a qPlus tuning fork [49] (resonance frequency f 0 ≈ 29 1 kHz, spring constant . k = 1 8 kN m . -1 , quality factor Q ≈ 3 × 10 ) in 4 ultrahigh vacuum ( p = 2 × 10 -10 mbar) and at a temperature of 6.2 K. As substrate we used Cu(111) single crystal covered with &gt; 20ML of NaCl. The molecules were deposited onto the sample kept at ≈ 7 K inside the microscope. Bias voltages are given as sample bias with respect to the tip. Positive constant-height offsets ∆ z correspond to a distance decrease with respect to the AFM ∆ f set point above the clean NaCl. For AC-STM

experiments, the ac voltage pulses were produced by an arbitrary waveform generator (Agilent 33522A) and fed to the microscope head through semirigid coaxial high frequency cables.

## Calculating molecular mode parameters

The dependence of the electronic wave functions of the molecule on the nuclear positions is obtained by expanding the electronic Hamiltonian around the reference configuration R 0 up to first order in the displacement coordinate q . By definition, the first derivative of the ground state energy for the reference electronic state vanishes. This is not the case, though, for the each component of the electronic Hamiltonian. In terms of the normal modes coordinates, we thus obtain

$$\hat { H } _ { \text{e-ph} } = \sum _ { l, \alpha } \left ( \frac { \partial \hat { H } _ { \text{e} } } { \partial q _ { l, \alpha } } \right ) _ { q _ { l, \alpha } = 0 } \sum _ { \zeta } v _ { l, \alpha } ^ { \zeta } \sqrt { \frac { \hbar } { 2 M _ { l } \omega _ { \zeta } } } \left ( \hat { a } _ { \zeta } ^ { \dagger } + \hat { a } _ { \zeta } \right ),$$

with M l the mass of the atom l and v ζ the displacement vector of the normal mode ζ with frequency ω ζ . The molecular Hamiltonian ˆ H mol presented in the main text corresponds to ˆ H e ( r R , 0 ). The matrix element describing the electron phonon coupling strength is given by

$$g _ { i j } ^ { \zeta } = \sum _ { l \alpha } \langle i \, | \, \frac { \partial \hat { H } _ { e } } { \partial q _ { l \alpha } } \, | \, j \rangle \mathfrak { v } _ { l \alpha } ^ { \zeta } \sqrt { \frac { \hbar { \, } } { 2 M _ { l } \omega _ { \zeta } } }.$$

These matrix elements can be evaluated by the INELASTICA python package [41]. To evaluate them we first obtained the molecular geometry with the SIESTA DFT implementation [40]. We employed a Single-Zeta-Pseudo (SZP) basis set and used the GGA-PBE exchange functional. For details about the calculation of the mode energies ¯ hω ζ and electron phonon couplings g ζ ij we refer to the paper by Frederiksen et al. [41].

## Separation of symmetric and antisymmetric modes

In this section, we aim to demonstrate the contrasting effects of symmetric and antisymmetric modes. We find that the symmetric modes solely lead to an energy renormalization depending on the total number of charges on the molecule. Conversely, the coupling to the antisymmetric modes varies across different orbital occupations, thereby exerting an influence on the charge distribution of the molecular ground state.

Our discussion begins with the component of the Hamiltonian that encompasses the vibronic degrees of freedom as well as the electron-phonon coupling

$$\hat { H } _ { \text{vib} } + \hat { H } _ { e - \text{ph} } = \sum _ { \zeta } \hbar { \omega } _ { \zeta } \left ( \hat { a } _ { \zeta } ^ { \dagger } \hat { a } _ { \zeta } + \frac { 1 } { 2 } \right ) + \sum _ { \zeta \sigma } \sum _ { i, j = L \pm } g _ { i j } ^ { \zeta } \hat { d } _ { i \sigma } ^ { \dagger } \hat { d } _ { j \sigma } \left ( \hat { a } _ { \zeta } ^ { \dagger } + \hat { a } _ { \zeta } \right ).$$

The sum is taken over all contributing modes ζ , the spin σ , and the degenerate LUMOs i and j , as we assume that any excess charge on the molecule to be confined solely to the LUMOs, thus restricting the electron-phonon coupling to these two orbitals. The electron phonon coupling strength is proportional to the matrix elements with respect to the derivative of the electronic Hamiltonian, see Eq. (9),

$$g _ { i j } ^ { \zeta } \in \langle i \, | \, \frac { \partial \hat { H } _ { e } } { \partial q _ { l \alpha } } \, | \, j \rangle.$$

The derivative ∂H ˆ e ∂q lα is an irreducible tensor operator which transforms in the same way as the respective normal mode ζ . The wave functions i and j transform according to the E g irreducible representation of the D 4 h point group. Therefore, for the matrix element in Eq. (11) to be nonzero we can only use modes ζ such that their irreducible representation Γ ζ is contained in the symmetric product

$$[ E _ { g } \otimes E _ { g } ] _ { + } = A _ { 1 g } \oplus B _ { 1 g } \oplus B _ { 2 g }.$$

The modes transforming according to the A 1 g irreducible representation couple to the global excess charge on the LUMOs. We call these the symmetric modes. The associated molecular distortion preserves the symmetry of the

system. Therefore, they are not contributing to the JT effect. It is important to consider them, as they give a renormalization to the eigenenergies of the molecule depending on the charge state and they lead to a broadening in the transition rates.

On the other hand, the B 1 g and B 2 g modes couple to the unbalance between the occupation of the LUMOs. In particular, this imbalance concerns, for the B 1 g modes, the real valued basis L xz and L yz shown in Fig. 5. Contrary to the totally symmetric A 1 g modes, these modes can lower the symmetry of the molecule. Thus, they are responsible for the JT distortion.

## Energy renormalization due to symmetric modes

In this section, we delve into the impact of the symmetric modes. All relevant parameters are shown in Table I. The symmetric modes couple diagonally to the LUMOs, i.e. i = j in Eq. (10),

$$\hat { H } _ { \text{vib} } ^ { A } + \hat { H } _ { e - \text{ph} } ^ { A } = \sum _ { \alpha \in \{ A _ { 1 g } \} } \hbar { \omega } _ { \alpha } \left ( \hat { a } _ { \alpha } ^ { \dagger } \hat { a } _ { \alpha } + \frac { 1 } { 2 } \right ) + \sum _ { \alpha \in \{ A _ { 1 g } \} } g _ { \alpha } \sum _ { i \sigma } \hat { d } _ { i \sigma } ^ { \dagger } \hat { d } _ { i \sigma } \left ( \hat { a } _ { \alpha } ^ { \dagger } + \hat { a } _ { \alpha } \right ).$$

All modes couple equally to the occupation of the L xz and L yz orbitals and thus we omit the subscript specifying the LUMO from the coupling strength g α . The sum over α encompasses all modes transforming according to the A 1 g irreducible representation and g α is the coupling strength to the LUMOs of mode α . We apply the Lang-Firsov transformation to the full Hamiltonian

$$\tilde { \hat { H } } = e ^ { \hat { S } } \left ( \hat { H } _ { \text{mol} } + \hat { H } _ { \text{tip} } + \hat { H } _ { \text{tun} } + \hat { H } _ { \text{as} } + \hat { H } _ { \text{vib} } + \hat { H } _ { e - \text{ph} } \right ) e ^ { - \hat { S } } \,,$$

where the operator

$$\hat { S } = \sum _ { \alpha \in \{ A _ { 1 g } \} } \sum _ { i \sigma } \lambda _ { \alpha } \hat { n } _ { i \sigma } \left ( \hat { a } _ { \alpha } ^ { \dagger } - \hat { a } _ { \alpha } \right ).$$

with λ α = g α ¯ hω α . This transforms the annihilation and creation operators on the molecule as

$$\tilde { \hat { d } } _ { i \sigma } = \hat { d } _ { i \sigma } \hat { X }, \quad \tilde { \hat { d } } _ { i \sigma } ^ { \dagger } = \hat { d } _ { i \sigma } ^ { \dagger } \hat { X } ^ { \dagger }.$$

with ˆ = X e -∑ α λ α ( ˆ b † α -ˆ b α ) . The phonon operators are shifted by

$$\tilde { \hat { a } } _ { \alpha } = \hat { a } _ { \alpha } - \lambda _ { \alpha } \sum _ { i \sigma } \hat { n } _ { i \sigma }, \quad \tilde { \hat { a } } _ { \alpha } ^ { \dagger } = \hat { a } _ { \alpha } ^ { \dagger } - \lambda _ { \alpha } \sum _ { i \sigma } \hat { n } _ { i \sigma }.$$

The transformed molecular Hamiltonian thus obtains a renormalization of the single-particle LUMO energies ˜ ϵ i = ϵ i -∑ α ∈ A 1 g g 2 α ¯ hω α . and of the direct Coulomb interaction between the LUMOs ˜ V ijji = V ijji -2 ∑ α ∈ A 1 g g 2 α ¯ hω α , with i, j = L xz , L yz . Furthermore, the symmetric electron-phonon coupling term vanishes, leaving only the harmonic vibrations in the Hamiltonian which describes the effect of the symmetric modes. The renormalization of the singleparticle LUMO energies and the corresponding Coulomb interaction terms leads to a stabilization of charged states on the molecule. The Lang-Firsov transformation of the tunneling Hamiltonian yields the ˆ X operator introduced in Eq. (16). Its matrix elements in the harmonic oscillator eigenbasis are the Franck-Condon factors giving the Poisson distribution in the rates of Eqs. (7) of the main text. Ultimately, the ˆ X operator and thus the coupling to the A 1 g modes is responsible for the broadening of the calculated transition rates.

## Treatment of antisymmetric modes

The antisymmetric modes are more sensitive to the actual distribution of the charge between the LUMOs, rather then to their global occupation. The electron vibron coupling Hamiltonian, written for the real basis depicted in Fig. 5, reads

$$\hat { H } _ { e - p h } ^ { B } = \sum _ { \beta \in \{ B _ { 1, g } \} } g _ { \beta } \left ( \hat { n } _ { x z } - \hat { n } _ { y z } \right ) \left ( \hat { a } _ { \beta } ^ { \dagger } + \hat { a } _ { \beta } \right ) + \sum _ { \beta \in \{ B _ { 2, g } \} } g _ { \beta } \left ( \hat { d } _ { x z } ^ { \dagger } \hat { d } _ { y z } + \hat { d } _ { y z } ^ { \dagger } \hat { d } _ { x z } \right ) \left ( \hat { a } _ { \beta } ^ { \dagger } + \hat { a } _ { \beta } \right ),$$

TABLE I. Reorganization energy, dimensionless coupling λ and energy quantum of the A 1 g , B 1 g and B 2 g molecular modes. The reorganization energy is obtained as ∆ E reorg = ¯ hωλ 2 , with λ = g/hω. ¯ The vibrational energy and the reorganization energy are given in meV.

| A 1 g     | modes   | modes   | B 1 g modes   | B 1 g modes   | B 1 g modes   | B 2 g modes   | B 2 g modes   | B 2 g modes   |
|-----------|---------|---------|---------------|---------------|---------------|---------------|---------------|---------------|
| ∆ E reorg | λ       | ¯ hω    | ∆ E reorg     | λ             | ¯ hω          | ∆ E reorg     | λ             | ¯ hω          |
| 61.0      | 1.435   | 29.6    | 28.0          | 0.389         | 185.0         | 11.8          | 0.292         | 148.0         |
| 25.5      | 0.384   | 173.3   | 9.6           | 0.235         | 172.7         | 3.0           | 0.239         | 53.7          |
| 21.0      | 0.522   | 76.7    | 6.1           | 0.269         | 84.9          | 2.9           | 0.130         | 172.5         |
| 16.2      | 0.406   | 98.4    | 1.1           | 0.085         | 159.1         | 2.3           | 0.132         | 134.1         |
| 8.3       | 0.222   | 167.3   | 0.838         | 0.078         | 138.2         | 1.4           | 0.108         | 120.4         |
| 2.6       | 0.126   | 164.1   | 0.78          | 0.080         | 120.3         | 1.2           | 0.108         | 106.2         |
| 2.2       | 0.109   | 186.5   | 0.38          | 0.047         | 169.3         | 1.0           | 0.298         | 11.2          |
| 1.3       | 0.080   | 195.8   | 0.34          | 0.056         | 110.5         | 0.7           | 0.061         | 178.2         |
| 0.5       | 0.059   | 135.2   | 0.29          | 0.12          | 19.4          | 0.4           | 0.048         | 197.0         |
| 0.5       | 0.062   | 120.5   | 0.10          | 0.034         | 90.9          | 0.2           | 0.093         | 25.4          |
| 0.02      | 0.012   | 110.5   | 0.046         | 0.027         | 62.7          | 0.03          | 0.020         | 76.7          |
| 0.007     | 0.005   | 359.7   | 0.051         | 0.0038        | 359.7         | 0.009         | 0.005         | 358.6         |
| 0.006     | 0.010   | 67.0    | 0.0024        | 0.0035        | 195.4         | 0.005         | 0.006         | 117.1         |
| 0.004     | 0.003   | 358.9   | 0.0001        | 0.0006        | 359.0         | 0.002         | 0.003         | 359.04        |

where β encompasses all modes transforming according to either B 1 g or B 2 g . This Hamiltonian gives rise to the JahnTeller (JT) effect, as each deformation favours the occupation of a determined LUMO thus lifting their degeneracy. The more complex form of the electron-phonon coupling combined with the structure of the exchange Hamiltonian makes the Lang-Firsov transformation unuseful. The latter does not decouple in this case the electronic from the vibronic degrees of freedom. Therefore, we treat the antisymmetric displacements classically and estimate the ground state confirmation by searching for the minima of adiabatic potential energy surfaces (APES). To this end, we introduce canonical displacements and momenta as follows

$$\hat { x } _ { \beta } = \sqrt { \hbar { \omega } _ { \beta } } \left ( \hat { a } _ { \beta } ^ { \dagger } + \hat { a } _ { \beta } \right ), \quad \hat { p } _ { \beta } = \frac { i \hbar { \ } } { 2 \sqrt { \hbar { \omega } _ { \beta } } } \left ( \hat { a } _ { \beta } ^ { \dagger } - \hat { a } _ { \beta } \right ).$$

This transformation allows us to express the Hamiltonian associated to the modes with the B 1 g and B 2 g symmetry as

$$\hat { H } _ { \text{vib} } ^ { B } + \hat { H } _ { e - \text{ph} } ^ { B } = \sum _ { \beta \in \{ B _ { 1 / 2 _ { \alpha } } \} } \left ( \omega _ { \hat { \beta } } ^ { 2 } \hat { \beta } _ { \beta } ^ { 2 } + \frac { 1 } { 4 } \hat { x } _ { \beta } ^ { 2 } \right ) + \sum _ { \beta \in \{ B _ { 1 \alpha } \} } \frac { g _ { \beta } } { \sqrt { \hbar { \omega } _ { \beta } } } \left ( \hat { n } _ { x z } - \hat { n } _ { y z } \right ) \hat { x } _ { \beta } + \sum _ { \beta \in \{ B _ { 2 \alpha } \} } \frac { g _ { \beta } } { \sqrt { \hbar { \omega } _ { \beta } } } \left ( \hat { d } _ { x z } ^ { \dagger } \hat { d } _ { y z } + \hat { d } _ { y z } ^ { \dagger } \hat { d } _ { x z } \right ) \hat { x } _ { \beta }. \ ( 2 0 )$$

We further simplify the problem by ignoring the kinetic terms and by turning ˆ x β into the classical coordinate x β . At this point we should diagonalize the effective Hamiltonian

$$\hat { H } _ { \text{eff} } ( \{ x _ { \beta } \} ) = \tilde { \hat { H } } _ { \text{mol} } + \hat { H } _ { \text{as} } + \sum _ { \beta } \frac { x _ { \beta } ^ { 2 } } { 4 } + \hat { H } _ { \text{e-ph} } ^ { B } ( \{ x _ { \beta } \} )$$

and look for the ground state energy and configuration for each charge sector. Given the presence of 14 B 1 g and 14 B 2 g modes, determining the energy minimum of this potentials proves impractical. A first screening on the modes can be performed using their reorganization energy ∆ E reorg ≡ ¯ hωλ 2 . As can be seen from Table I, the B 1 g modes provide a higher reorganization energy than the B 2 g modes. Therefore, we neglect the B 2 g modes and focus our attention solely on the B 1 g modes.

## Many-body eigenstates

A crucial role for the understanding of the experiments is played by the effective Hamiltonian in Eq. (21) from which the low energy many-body eigenstates of CuPc in its neutral, anionic and dianionic configuration can be extracted. We restrict ourselves to the Fock space generated from the single particle basis of four frontier orbitals: i.e. the SOMO ( S ), HOMO ( H ) and the two degenerate LUMOs. In particular, we will only consider, in the following, the neutral anionic and dianionic states with respectively 0,1 and 2 electrons in the LUMOs. This further simplification does not

![image](https://ac.marcodsn.me/arxiv-markdown-images/2408.00478v1-figure-b1155751-beb6-4e83-8a2b-f3ab13f7f5b3.jpg)

FIG. 5. Frontier molecular orbitals of CuPc. We show the SOMO ( S ), HOMO ( H ) and the LUMOs in their real valued representation L xz and L yz .

![image](https://ac.marcodsn.me/arxiv-markdown-images/2408.00478v1-figure-0db8baa6-04f2-4351-a1da-6101e6023c2f.jpg)

![image](https://ac.marcodsn.me/arxiv-markdown-images/2408.00478v1-figure-225e9a7d-c53c-4bd3-b27b-93666a69ffbf.jpg)

![image](https://ac.marcodsn.me/arxiv-markdown-images/2408.00478v1-figure-57e13b82-a07f-4cf1-a4ce-61e6318c782a.jpg)

TABLE II. Coulomb integrals between the frontier orbitals. We indicate with U i = V iiii the Hubbard energy of the orbital i , with U ij = V iijj and J ex ij = V ijji the direct Coulomb integral and the exchange between the orbitals i and j . Finally, J p ij = V ijij and ˜ J p ijk = V ijik refer to the pair hopping and the split pair hopping, respectively. Taken from Ref. [35].

| U S               | 11.352 eV J ex HL = - ˜ J p H+ - 548   | meV     |
|-------------------|----------------------------------------|---------|
| U H               | 1.752 eV J ex + -                      | 258 meV |
| U L = U ± = U + - | 1.808 eV J P + -                       | 168 meV |
| U SH              | 1.777 eV J ex SL = - ˜ J p S+ -        | 9 meV   |
| U SL              | 1.993 eV J ex SH = J p SH              | 2 meV   |
| U HL              | 1.758 eV                               |         |

affect significantly the low energy eigenstates of the neutral and anionic molecule, as obtained by comparison with numerical diagonalization performed in the full Fock space associated to the four frontier orbitals. In the dianionic case this approximation excludes a pair of doublets with triple occupation of the LUMOs. The latter are though not relevant to explain the experimental results, as these states cannot be accessed with a single electron tunnelling event from the anionic ground state.

As the LUMOs are degenerate, they allow for different representations. We show in Fig. 5 the frontier orbitals with the LUMOs in their real valued representation ( L xz , L yz ). However, in the discussion of the many-body eigenstates, we will also use their rotationally invariant representation which is given by [36]

$$| L + \rangle = \frac { 1 } { \sqrt { 2 } } \left ( | L _ { x z } \rangle + i \, | L _ { y z } \rangle \right ), \ \ | L - \rangle = \frac { 1 } { \sqrt { 2 } } \left ( | L _ { x z } \rangle - i \, | L _ { y z } \rangle \right ).$$

The Hamiltonian we discuss is ˆ H eff ( x 1 ) = ˜ ˆ H mol + x / 2 1 4 + ˆ H B e -ph ( x 1 ), where, to start, we restrict ourselves to the B 1 g mode with the highest reorganization energy. The renormalization due to the symmetric modes only introduces the same shift to all eigenenergies of a given particle number. In practice, we have

$$\hat { H } _ { \text{eff} } ( x _ { 1 } ) = \sum _ { i } ( \tilde { \varepsilon } _ { i } + \Delta ) \hat { n } _ { i } + \frac { 1 } { 2 } \sum _ { i j k l } \sum _ { \sigma \sigma ^ { \prime } } \tilde { V } _ { i j k l } \hat { d } _ { \omega } ^ { l } \hat { d } _ { k \sigma ^ { \prime } } \hat { d } _ { l \sigma ^ { \prime } } \hat { d } _ { j \sigma } + \frac { x _ { 1 } ^ { 2 } } { 4 } + \sum _ { \sigma } \frac { g _ { 1 } } { \sqrt { h \omega _ { 1 } } } \left ( \hat { d } _ { L + \sigma } ^ { l } \hat { d } _ { L - \sigma } + \hat { d } _ { L - \sigma } ^ { l } \hat { d } _ { L + \sigma } \right ) x _ { 1 } \,.$$

The summation indices i , j , k and l run over all four frontier orbitals, i.e. i j k l , , , = S H L , , + , L -. The model comprises the renormalized single particle energies ˜ ε i = ε i -∑ α g 2 α /hω ¯ α with ε S = -12 0 eV, . ε H = -11 7 eV, . ε L + /L -= -10 7 eV, . also shifted by a crystal field correction of ∆ = 2 44eV. . The Coulomb parameters concerning the direct interaction are also renormalized by the reorganization energy of the symmetric modes, i.e. ˜ V ijkl = V ijkl -2 ∑ α g 2 α /hω ¯ α . Exchange and pair hopping contributions are instead not influenced by the LangFirsov transformation. The bare Coulomb parameters, calculated by direct integration from the frontier orbitals [36], are given in Table II. The energy of the considered mode is ¯ hω 1 = 185meV and the coupling strength g 1 = 72meV.

FIG. 6. Adiabatic potential energy surfaces along the deformation coordinate belonging to the B 1 g mode with the highest reorganization energy. The renormalization of the single particle and direct Coulomb interaction due to the A 1 g modes is taken into account and the chemical potential µ = -4 65 eV. (a) The potential energy of the neutral state is a parabola with . its minimum in the reference configuration ( x 1 = 0), as the electron-phonon coupling does not contribute due to the empty LUMOs. (b) The anionic potential energy is not minimal in the reference configuration. Either the triplet | 4 , T S z 0 , xz ⟩ or the triplet | 4 , T S z 0 , yz ⟩ are the ground state depending on the sign of the deformation coordinate in equilibrium.

![image](https://ac.marcodsn.me/arxiv-markdown-images/2408.00478v1-figure-c4caa00a-99da-455e-b2e7-01999b242e13.jpg)

## The neutral molecule

We start our discussion of the eigenstates form the neutral molecule. By diagonalizing Eq. (23) in the reference position R 0 , i.e. x 1 = 0, we find a spin degenerate doublet

$$\left | 3, D _ { 0 } ^ { \uparrow } \right \rangle = \hat { d } _ { S \uparrow } ^ { \dagger } \left | \Omega \right \rangle, \quad \left | 3, D _ { 0 } ^ { \downarrow } \right \rangle = \hat { d } _ { S \downarrow } ^ { \dagger } \left | \Omega \right \rangle,$$

where | Ω ⟩ is the reference state with doubly occupied HOMO but empty SOMO and LUMOs. We observe the characteristic unpaired spin in the SOMO, due to its localization and correspondingly high charging energy. The first-excited neutral energy level lies about 0 8 eV above the ground state [36] and thus we do not consider it in . our calculations. The lowest eigenvalue of the Hamiltonian Eq. (23) is depicted in Fig. 6 (a) as a function of the displacement coordinate x 1 given in units of √ ¯ hω 1 and represents the adiabatic potential energy surface (APES). As the LUMOs remain unoccupied in the neutral ground state, the electron-phonon coupling does not contribute to the potential energy and we only see the harmonic contribution due to the first antisymmetric mode.

## The anionic molecule

The effective Fock space for the CuPc anion comprises the eight states given by ˆ d † Sσ ˆ d Liσ ′ | Ω ⟩ with σ, σ ′ = ↑ , ↓ and i = ± , which are degenerate in the reference configuration as far as it concerns the direct Coulomb interaction. The exchange coupling between the electrons occupying the SOMO and the LUMO splits this eightfold degeneracy. The eigenstates in the reference configuration can thus be organized as two sets of triplets and singlets

$$\L e c \, \text{cong} \, \text{quantation} \, \text{can} \, \text{thus} \, \text{ be organized as two sets of triplets and singlets} \\ | 4, T _ { 0 } ^ { + 1 }, i \rangle & = \hat { d } _ { S \uparrow } ^ { \dagger } \hat { d } _ { L i \uparrow } ^ { \dagger } \left | \Omega \right \rangle, \\ | 4, T _ { 0 } ^ { 0 }, i \rangle & = \frac { 1 } { \sqrt { 2 } } \left ( \hat { d } _ { S \uparrow } ^ { \dagger } \hat { d } _ { L i \downarrow } ^ { \dagger } - \hat { d } _ { L i \uparrow } ^ { \dagger } \hat { d } _ { S \downarrow } ^ { \dagger } \right ) | \Omega \rangle \,, \\ | 4, T _ { 0 } ^ { - 1 }, i \rangle & = \hat { d } _ { S \downarrow } ^ { \dagger } \hat { d } _ { L i \downarrow } ^ { \dagger } \left | \Omega \right \rangle, \\ | 4, S _ { 1 }, i \rangle & = \frac { 1 } { \sqrt { 2 } } \left ( \hat { d } _ { S \uparrow } ^ { \dagger } \hat { d } _ { L i \downarrow } ^ { \dagger } + \hat { d } _ { L i \uparrow } ^ { \dagger } \hat { d } _ { S \downarrow } ^ { \dagger } \right ) | \Omega \rangle \,, \\ \, \text{d} \, \text{d} \, \text{egenerate ground states separated from the singlets by} \, 2 J _ { s r } \, = \, 1 \, & \text{meV. These two sets.}$$

with the triplet orbitally degenerate ground states separated from the singlets by 2 J SL = 18 meV. These two sets, labeled by i = ± , arise due to the degenerate LUMOs and are mixed by the coupling to the B 1 g vibrational mode.

̸

First, we discuss the two triplet states with S z = +1. The eigenstates of the Hamiltonian Eq. (23) when x 1 = 0 read

$$| 4, T _ { 0 } ^ { + 1 }, x z \rangle = \frac { 1 } { \sqrt { 2 } } \left ( | 4, T _ { 0 } ^ { + 1 }, + \rangle + | 4, T _ { 0 } ^ { + 1 }, - \rangle \right ), \ \ | 4, T _ { 0 } ^ { + 1 }, y z \rangle = \frac { i } { \sqrt { 2 } } \left ( | 4, T _ { 0 } ^ { + 1 }, + \rangle - | 4, T _ { 0 } ^ { + 1 }, - \rangle \right ).$$

We notice that these states are also eigenstates of Eq. (23) if x 1 = 0, thus giving a single set of eigenstates at all deformations, although the corresponding eigenenergies vary as function of x 1 , as shown in Fig. 6. From the definition of the rotationally invariant representation given in Eq. (25) we can express these anionic eigenstates states in the real valued representation as

$$| 4, T _ { 0 } ^ { + 1 }, x z \rangle = \hat { d } _ { S \uparrow } ^ { \dagger } \hat { d } _ { L _ { x x } \uparrow } ^ { \dagger } \, | \Omega \rangle \,, \ \ | 4, T _ { 0 } ^ { + 1 }, y z \rangle = \hat { d } _ { S \uparrow } ^ { \dagger } \hat { d } _ { L _ { y x } \uparrow } ^ { \dagger } \, | \Omega \rangle \,.$$

The same procedure can be applied to the other 4 triplet states. We obtain

$$& | 4, T _ { 0 } ^ { 0 }, i z \rangle = \frac { 1 } { \sqrt { 2 } } \left ( \hat { d } _ { S _ { 1 } ^ { \uparrow } } ^ { \downarrow } \hat { d } _ { L _ { i z } \downarrow } ^ { \uparrow } - \hat { d } _ { L _ { i z } \uparrow } ^ { \downarrow } \hat { d } _ { S _ { 1 } \downarrow } ^ { \uparrow } \right ) | \Omega \rangle \,, \\ & | 4, T _ { 0 } ^ { - 1 }, i z \rangle = \hat { d } _ { S _ { 1 } \downarrow } ^ { \dagger } \hat { d } _ { L _ { i z } \downarrow } ^ { \uparrow } | \Omega \rangle \,,$$

with i = x/y . Analogously, for the singlet states we obtain, in the real valued basis,

$$| 4, S _ { 1 }, i z \rangle = \frac { 1 } { \sqrt { 2 } } \left ( \hat { d } ^ { \dagger } _ { S \uparrow } \hat { d } ^ { \dagger } _ { L _ { i z \downarrow } } + \hat { d } ^ { \dagger } _ { L _ { i z \uparrow } } \hat { d } ^ { \dagger } _ { S \downarrow } \right ) | \Omega \rangle \,,$$

again with i = x/y . In each of the two dimensional subspaces with a given spin S and S z the Jahn-Teller Hamiltonian in Eq. (23) acquires the canonical form

$$H _ { J T } ^ { 1 - } = \frac { x _ { 1 } ^ { 2 } } { 4 } \mathbf 1 _ { 2 } + \frac { g _ { 1 } } { \sqrt { \hbar { \omega } _ { 1 } } } x _ { 1 } \sigma _ { x }$$

typical of a E ⊗ b 1 Jahn-Teller problem [45]. The sign of the displacement along the x 1 coordinate determines whether the triplet-singlet set with the added electron in the L xz or in the L yz orbital experience a decrease or increase in energy. This leads to a preference of one set of triplet-singlet states over the other depending on the direction of distortion. The APES depicted in Fig. 6 (b) illustrates this phenomenon. The plot shows four parabolas representing the potential energy of both triplet-singlet sets.

## The dianionic molecule

The effective Fock space for the dianionic molecule consists of twelve states, which, in the rotationally invariant basis, are represented as

$$\hat { d } _ { S \sigma } ^ { \dagger } \hat { d } _ { L i \sigma ^ { \prime } } ^ { \dagger } \hat { d } _ { L j \sigma ^ { \prime \prime } } ^ { \dagger } \left | \Omega \right \rangle,$$

where i, j = ± , and σ ′ = ¯ σ ′′ if i = j due to the Pauli exclusion principle. These states are degenerate for what concerns the direct Coulomb interaction (terms involving U in ˆ H mol ) with an energy E 5 = ˜ ε S +2(˜ ε H + ˜ ) + 5∆ + ε L ˜ U H + ˜ U L + 2( ˜ U SH + ˜ U SL ) + 4 ˜ U HL . By considering the exchange and pair hopping interactions (terms involving J in ˆ H mol ), this degeneracy is lifted. We first discuss the spectrum and the eigenstates in the reference configuration x 1 = 0. The eigenenergies are given with respect to the degenerate energy. The scheme of the energy levels is given in Fig. 7 (a). The ground state is a spin quadruplet

$$\{ \alpha \cdot \, \text{$1$} \, \text{$1$} \, \text{$1$} \, \text{$1$} \, \text{$1$} \, \text{$1$} \, \text{$1$} \, \text{$1$} \, \text{$1$} \, \text{$1$} \, \text{$1$} \, \text{$1$} \, \text{$1$} \, \text{$1$} \, \text{$1$} \, \text{$1$} \, \text{$1$} \, \text{$1$} \, \text{$1$} \, \text{$1$} \, \text{$1$} \, \text{$1$} \, \text{$1$} \} \, \text{$1$} \, \text{$1$} \, \text{$1$} \, \text{$1$} \, \text{$1$} \, \text{$1$} \, \text{$1$} \, \text{$1$} \, \text{$1$} \, \text{$1$} \, \text{$1$} \, \text{$1$} \, \text{$1$} \, \text{$1$} \, \text{$1$} \, \text{$1$} \, \text{$1$} \, \text{$1$} \, \text{$1$} \, \text{$1$} \, \text{$1$} \, \text{$1$} \\ & | 5, Q _ { \delta } ^ { \frac { 3 } { 2 } } \rangle = \frac { 1 } { \sqrt { 3 } } \left ( \hat { d } _ { S _ { 1 } } ^ { \dagger } \hat { d } _ { L + \downarrow } \hat { d } _ { L - \uparrow } \hat { d } _ { S _ { 1 } } ^ { \dagger } - \hat { d } _ { S _ { 1 } } ^ { \dagger } \hat { d } _ { L - \uparrow } \hat { d } _ { L + \downarrow } + \hat { d } _ { S _ { 1 } } ^ { \dagger } \hat { d } _ { L - \downarrow } \hat { d } _ { L + \uparrow } \hat { d } _ { L - \downarrow } \hat { d } _ { L - \downarrow } + \hat { d } _ { L - \uparrow } \hat { d } _ { S _ { 1 } } ^ { \dagger } \hat { d } _ { L + \downarrow } \hat { d } _ { L + \downarrow } \right ) | \Omega \rangle \,, \\ & | 5, Q _ { 0 } ^ { - \frac { 3 } { 2 } } \rangle = \hat { d } _ { S _ { 1 } } ^ { \dagger } \hat { d } _ { L + \downarrow } \hat { d } _ { L - \downarrow } \hat { d } _ { L + \downarrow } ^ { \dagger } \hat { d } _ { L - \downarrow } \hat { d } _ { L + \downarrow } ^ { \dagger } | \Omega \rangle \,.$$

![image](https://ac.marcodsn.me/arxiv-markdown-images/2408.00478v1-figure-3cd4796d-9711-4118-858b-8203af0fbe2d.jpg)

(b)

FIG. 7. (a) Scheme of the exchange interaction's influence on the dianionic spectrum at x 1 = 0. The 12-fold degeneracy gets lifted by including the exchange interaction. The states | 5 , Q S z 0 ⟩ , | 5 , D σ 1 ⟩ and | 5 , D σ 2 ⟩ correspond to the three low lying parabolas in panel (b). On the contrary, the states | 5 , D σ 3 ⟩ and | 5 , D σ 4 ⟩ are affected by the JT distortion and exhibit an anti-crossing in (b). (b) APES of the dianionic molecule. There are 3 parabolas which have their minimum in the x 1 = 0 reference configuration. They correspond to the 8 states with the two added electrons distributed between the L xz and L yz orbitals. The two lines exhibiting an anti-crossing are describing the energy of the states with the electrons located in one of the two orbitals. One mode is not enough to lower their energies enough to become the ground state.

![image](https://ac.marcodsn.me/arxiv-markdown-images/2408.00478v1-figure-70e7e471-a512-4bde-95e8-82cc5a727686.jpg)

The energy of this quadruplet is E Q 0 = E 5 -J ex + --2 J ex SL . The first excited states form a spin degenerate doublet, with energy E D 1 = E 5 -J ex + -+ J ex SL . The corresponding eigenstates are given by

$$| 5, D _ { 1 } ^ { \uparrow } \rangle & = \frac { 1 } { \sqrt { 6 } } \left ( 2 \, \hat { d } _ { L + \uparrow } ^ { \downarrow } \hat { d } _ { L - \uparrow } ^ { \downarrow } \hat { d } _ { S \downarrow } ^ { \dagger } + \hat { d } _ { S \uparrow } ^ { \dagger } \hat { d } _ { L - \uparrow } ^ { \dagger } \hat { d } _ { L + \downarrow } ^ { \dagger } - \hat { d } _ { S \uparrow } ^ { \dagger } \hat { d } _ { L + \uparrow } ^ { \dagger } \hat { d } _ { L - \downarrow } ^ { \dagger } \right ) | \Omega \rangle \,, \\ | 5, D _ { 1 } ^ { \downarrow } \rangle & = \frac { 1 } { \sqrt { 6 } } \left ( 2 \, \hat { d } _ { S \uparrow } ^ { \downarrow } \hat { d } _ { L + \downarrow } \hat { d } _ { L - \downarrow } ^ { \downarrow } + \hat { d } _ { L + \uparrow } ^ { \dagger } \hat { d } _ { S \uparrow } ^ { \dagger } \hat { d } _ { L - \downarrow } ^ { \dagger } - \hat { d } _ { L - \uparrow } ^ { \dagger } \hat { d } _ { S \downarrow } ^ { \dagger } \hat { d } _ { L + \downarrow } ^ { \dagger } \right ) | \Omega \rangle \,.$$

The following three energy levels are also spin doublets, with the following eigenenergies

$$E _ { D _ { 2 } } = E _ { 5 } - J _ { + - } ^ { \text{p} } - J _ { \text{SL} } ^ { \text{ex} } \,, \ \ E _ { D _ { 3 } } = E _ { 5 } + J _ { + - } ^ { \text{p} } - J _ { \text{SL} } ^ { \text{ex} } \,, \ \ E _ { D _ { 4 } } = E _ { 5 } + J _ { + - } ^ { \text{ex} } - J _ { \text{SL} } ^ { \text{ex} } \,.$$

The corresponding eigenstates are,

$$\mu _ { 1 } \mu _ { 2 / 3 } & = \frac { 1 } { \sqrt { 2 } } \left ( \hat { d } _ { S \uparrow } ^ { \dagger } \hat { d } _ { L + \uparrow } ^ { \dagger } \hat { d } _ { L + \downarrow } ^ { \dagger } \mp \hat { d } _ { S \uparrow } ^ { \dagger } \hat { d } _ { L - \uparrow } ^ { \dagger } \hat { d } _ { L - \downarrow } ^ { \dagger } \right ) | \Omega \rangle \,, \\ & | 5, D _ { 2 / 3 } ^ { \downarrow } \rangle & = \frac { 1 } { \sqrt { 2 } } \left ( \hat { d } _ { L + \uparrow } ^ { \dagger } \hat { d } _ { S \downarrow } ^ { \dagger } \hat { d } _ { L + \downarrow } ^ { \dagger } \mp \hat { d } _ { L - \uparrow } ^ { \dagger } \hat { d } _ { S \downarrow } ^ { \dagger } \hat { d } _ { L - \downarrow } ^ { \dagger } \right ) | \Omega \rangle \,,$$

and, respectively,

$$| 5, D _ { 4 } ^ { \uparrow } \rangle & = \frac { 1 } { \sqrt { 2 } } \left ( \hat { d } _ { S \uparrow } ^ { \dagger } \hat { d } _ { L + \uparrow } ^ { \dagger } \hat { d } _ { L - \downarrow } ^ { \dagger } + \hat { d } _ { S \uparrow } ^ { \dagger } \hat { d } _ { L - \uparrow } ^ { \dagger } \hat { d } _ { L + \downarrow } ^ { \dagger } \right ) | \Omega \rangle \,, \\ | 5, D _ { 4 } ^ { \downarrow } \rangle & = \frac { 1 } { \sqrt { 2 } } \left ( \hat { d } _ { L + \uparrow } ^ { \dagger } \hat { d } _ { S \downarrow } ^ { \dagger } \hat { d } _ { L - \downarrow } ^ { \dagger } + \hat { d } _ { L - \uparrow } ^ { \dagger } \hat { d } _ { S \downarrow } ^ { \dagger } \hat { d } _ { L + \downarrow } ^ { \dagger } \right ) | \Omega \rangle \,. \\ \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \,.$$

The first eight states corresponding to the lower energies are not affected by the electron vibron coupling to the modes of symmetry B 1 g . This is evident in Fig. 7, where we observe three parabolas - one is fourfold and the other two are twofold degenerate - with a minimum at x 1 = 0. This antisymmetric mode couples instead the eigenstates | 5 , D σ 3 ⟩ and | 5 , D σ 4 ⟩ . For each of the two spin components, the Jahn-Teller Hamiltonian projected on the corresponding two dimensional basis reads

$$H _ { J T } ^ { 2 - } ( x _ { 1 } ) & = \frac { x _ { 1 } ^ { 2 } } { 4 } \mathbf 1 _ { 2 } - \frac { J _ { + - } ^ { \text{ex} } - J _ { + - } ^ { \text{p} } } { 2 } \sigma _ { z } + 2 \frac { g _ { 1 } } { \sqrt { \hbar { \omega } _ { 1 } } } x _ { 1 } \sigma _ { x } \,.$$

Due to the finite energy splitting in the reference configuration ( x 1 = 0) this Hamiltonian describes the pseudo-JahnTeller effect [45]. The existence of localized minima corresponding to deformed configurations is only guaranteed for a

sufficiently large coupling, i.e. g &gt; √ ¯ hω 1 ( J ex + --J p + -) / 4, as can be deduced by analysing the two highest eigenvalues of the dianionic effective Hamiltonian:

$$E _ { D _ { 3 / 4 } } ( x _ { 1 } ) = E _ { 5 } - J _ { \text{SL} } ^ { \text{ex} } + \frac { J _ { + - } ^ { \text{ex} } + J _ { + - } ^ { p } } { 2 } + \frac { x _ { 1 } ^ { 2 } } { 4 } \mp \sqrt { \left ( \frac { J _ { + - } ^ { \text{ex} } - J _ { + - } ^ { p } } { 2 } \right ) ^ { 2 } + \frac { 4 g _ { 1 } ^ { 2 } } { \hbar { \omega } _ { 1 } } x _ { 1 } ^ { 2 } } \,.$$

Differently from the anionic case, here also the eigenstates depend on the deformation coordinate x 1 . In particular we obtain

$$| 5, D _ { 3 } ^ { \sigma } ( x _ { 1 } ) \rangle & = \cos \theta ( x _ { 1 } ) \, | 5, D _ { 3 } ^ { \sigma } \rangle - \text{sgn} \, x _ { 1 } \sin \theta ( x _ { 1 } ) \, | 5, D _ { 4 } ^ { \sigma } \rangle \, \\ | 5, D _ { 4 } ^ { \sigma } ( x _ { 1 } ) \rangle & = \text{sgn} \, x _ { 1 } \sin \theta ( x _ { 1 } ) \, | 5, D _ { 3 } ^ { \sigma } \rangle + \cos \theta ( x _ { 1 } ) \, | 5, D _ { 4 } ^ { \sigma } \rangle$$

where the mixing angle is defined as

$$\theta ( x _ { 1 } ) = \frac { 1 } { 2 } \arccos \left ( 1 + \frac { 1 6 g _ { 1 } ^ { 2 } x _ { 1 } ^ { 2 } } { \hbar { \omega } _ { 1 } ( J _ { + - } ^ { \text{ex} } - J _ { + - } ^ { p } ) ^ { 2 } } \right ) ^ { - \frac { 1 } { 2 } } \,.$$

Within the single mode approximation considered here we obtain a mixing angle at equilibrium θ eq ≈ 0 24 . π . The larger is the electron phonon coupling g with respect to the vibrational energy ¯ hω 1 and the energy splitting J ex + --J p + -of the reference configuration, the larger is the absolute value of the equilibrium configuration. Thus, the mixing angle converges to θ = π/ 4. In view of the enhancement of the coupling introduced by the other molecular modes and, above all, by the insulating substrate, we expand the equilibrium eigenstate in the vicinity of this limiting case. The lowest energy state having the equilibrium configuration with x 1 eq , &lt; 0 reads

$$\cdots \quad \omega \cdots \quad \omega \cdots \quad \omega \cdots \quad \cdots \\ | 5, D _ { 3 } ^ { \sigma } ( x _ { 1, \epsilon } ) ) & = \cos \left ( \theta _ { \epsilon } - \frac { \pi } { 4 } \right ) \frac { 1 } { \sqrt { 2 } } \left ( | 5, D _ { 3 } ^ { \sigma } \right ) + | 5, D _ { 4 } ^ { \sigma } \right ) \right ) - \sin \left ( \theta _ { \epsilon } - \frac { \pi } { 4 } \right ) \frac { 1 } { \sqrt { 2 } } \left ( | 5, D _ { 3 } ^ { \sigma } \right ) - | 5, D _ { 4 } ^ { \sigma } \right ) \right ) \\ & = \text{sgn} \, \sigma \, \hat { d } _ { S \sigma } ^ { \dagger } \left [ \cos \left ( \theta _ { \epsilon } - \frac { \pi } { 4 } \right ) \, \hat { d } _ { L _ { x x } } ^ { \dagger } \hat { d } _ { L _ { x x } \downarrow } ^ { \dagger } + \sin \left ( \theta _ { \epsilon } - \frac { \pi } { 4 } \right ) \, \hat { d } _ { L _ { y x } \uparrow } ^ { \dagger } \hat { d } _ { L _ { y x } \downarrow } ^ { \dagger } \right ] | \Omega \right \rangle \\ & = \cos \left ( \theta _ { \epsilon } - \frac { \pi } { 4 } \right ) | 5, D _ { x x } ^ { \sigma } \right \rangle + \sin \left ( \theta _ { \epsilon } - \frac { \pi } { 4 } \right ) | 5, D _ { y y } ^ { \sigma } \right \rangle \,. \\ \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \, \end{array}$$

This state shows a predominant double occupation of the real valued L xz orbital ( 5 | , D σ xx ⟩ ) with a small admixture of the state with a doubly occupied L yz orbital ( 5 | , D σ yy ⟩ ). The roles of L xz and L yz are exchanged if x 1 eq , &gt; 0.

The experimental results we showed in Fig. 1 of the main text suggest that the dianionic ground state has both electrons in the same real valued LUMO. This is clearly not the case in our theoretical model if we consider one mode since the ground state is in the reference configuration where the lowest lying states are the quadruplets which have one electron in each LUMO. We can estimate the reorganization energy required to achieve the desired dianionic ground state by first considering that E D 3 -E Q 0 = J ex + -+ J p + -+ J ex SL = 435 meV. The total reorganization energy scales with the square of the electron-phonon coupling. Thus, the one associated to the B 1 g modes scales with the square of the electronic unbalance in the real valued LUMOs. It is thus necessary to achieve a reorganization energy per electron due to B 1 g modes of least 109 meV for the state | 5 , D σ 3 ( x 1 eq , ) ⟩ to become the dianionic ground state.

## Combination of multiple modes with the same symmetry

In this section, our focus is on combining multiple modes that belong to the same irreducible representation into a single effective mode. This allows us to analyze the potential energy of the system with contributions from all modes of the same irreducible representation. The combined system of molecule and salt has C 4 v symmetry. The B 1 g modes we calculated for the molecule are all in plane modes and thus they transform as B 1 modes with respect to the C 4 v point group. Therefore, we can not only combine the molecular modes among themselves but also include salt contributions to them. We begin this procedure by introducing the transformations of the displacements and momenta given by

$$\hat { Q } _ { n } = \sum _ { \beta \in \{ B _ { 1 } \} } A _ { n \beta } \hat { x } _ { \beta }, \ \hat { P } _ { m } = \sum _ { \beta \in \{ B _ { 1 } \} } B _ { m \beta } \hat { p } _ { \beta }.$$

Imposing canonical commutation relations on these new coordinates yields

$$\left [ \hat { Q } _ { n }, \hat { P } _ { m } \right ] = \sum _ { \beta \beta ^ { \prime } } A _ { n \beta } B _ { m \beta ^ { \prime } } \left [ \hat { x } _ { \beta }, \hat { p } _ { \beta ^ { \prime } } \right ] = i \hbar { \sum } _ { \beta \beta ^ { \prime } } A _ { n \beta } B _ { m \beta ^ { \prime } } \delta _ { \beta \beta ^ { \prime } } = i \hbar { \sum } _ { \alpha } A _ { n \beta } B _ { m \beta } \stackrel {! } { = } i \hbar { \delta } _ { n m }.$$

We further assume A to be orthogonal and thus deduce B = A . Moreover, we fix the first row of the transformation matrix to be

$$A _ { 1 \beta } = \frac { 1 } { g } \sqrt { \frac { g _ { \beta } ^ { 2 } } { \hbar { \omega } _ { \beta } } },$$

with g = √ ∑ β hω g 2 β ¯ β . This change of coordinate allows us rewrite the Jahn-Teller Hamiltonian as

$$\hat { H } _ { J T } ( \{ Q _ { n } \} ) = \frac { 1 } { 4 } \sum _ { n } Q _ { n } ^ { 2 } + g Q _ { 1 } \left ( \hat { n } _ { x z } - \hat { n } _ { y z } \right ) \,,$$

in which the electron-phonon coupling only concerns one effective JT mode, with effective coupling strength g . However, this transformation comes at the cost of a non-diagonal kinetic part in the Hamiltonian. We are still free to construct the matrix transformation A such that all the Jahn-Teller inactive modes remain orthogonal among themselves. We can thus consider them as phonon bath for the JT active mode Q 1 and we neglect them in our further discussions. This is justified as far as we are not interested in the dynamics of the JT distortion, but rather on the new equilibrium configurations. √

Only considering molecular contributions to the effective coupling strength yields a value of g = 218 meV corresponding to a reorganization energy of ∆ E B 1 = 48meV. From our previous discussion in ?? we know that this reorganization energy is not enough for the states with two electrons in one LUMO to become the new ground states. We can also deduce this from the APES in Fig. 8 (a), which is plotted along the effective JT coordinate Q AS := Q 1 in units of characteristic quantity

$$L _ { c } ^ { A S } = \sqrt { \frac { 1 } { g } \sqrt { \sum _ { \beta } \hbar { \omega } _ { \beta } g _ { \beta } ^ { 2 } } } \,,$$

where β runs over all considered modes. We notice that an effective Huang-Rhys factor can be defined, starting from Eq. (45), for the JT active mode Q AS as half of the ratio between the deformed equilibrium configuration Q AS min , of the anion and the characteristic scaling constant L AS c . In other terms, it reads

$$\lambda _ { \text{eff} } ^ { \text{AS} } = \frac { g } { L _ { c } ^ { \text{AS} } } = \left [ \frac { \left ( \sum _ { \beta } \frac { g _ { \beta } ^ { 2 } } { \hbar { \omega } _ { \beta } } \right ) ^ { 3 } } { \sum _ { \beta ^ { \prime } } \hbar { \omega } _ { \beta ^ { \prime } } g _ { \beta ^ { \prime } } ^ { 2 } } \right ] ^ { 1 / 4 }.$$

The latter reduces to the common expression λ = g β / hω (¯ β ) when the sums are restricted to a single mode. The deformed minima in the dianionic configuration are closer to the dianionic minimum energy but this still is not enough to explain the double charging of the same LUMO. However, if we consider also the contributions from the salt modes with the parameters given in the main text, we obtain g = 359 √ meV, yielding a reorganization energy of ∆ E B 1 = 129meV. We show the corresponding APES along the Q AS displacement coordinate in Fig. 8 (b). For the neutral and anionic charge states we observe no qualitative changes. For the dianionic molecule we see instead a qualitative change, as | 5 , D σ xx ⟩ becomes the new global ground state. The next highest lying states are the | 5 , D σ yy ⟩ states. These two pairs of states are split by 2 δ as , due to the asymmetry which we included in our Hamiltonian. Otherwise, analogously to the anionic pairs of triplet and singlet states, they would be degenerate. The energy splitting between the dianionic ground states | 5 , D σ xx ⟩ and the quadruplet states is, with our choice of parameters, 70 meV. This analysis confirms that the salt substrate in this system has an effect on the charge distribution within a single molecule going beyond mere stabilization of extra charges. We notice, moreover that, even with the contribution of the salt, λ eff ≈ 1 thus allowing for direct transitions between vibrational ground state with different charges. The classical limit of non-adiabatic vertical transition does not apply to the JT active mode. Analogous considerations can be reserved also to the symmetric ( A 1 ) modes, yielding a significant difference in the effective Huang-Rhys factors

$$\lambda _ { \text{eff} } ^ { S } = 2. 3 6 \,, \quad \lambda _ { \text{eff} } ^ { A S } = 1. 1 1 \,.$$

FIG. 8. APES of the neutral, anionic and dianionic molecule along the displacement coordinate of the active JT mode. (a) Considering all molecular modes is not enough to obtain the deformed states as ground states. (b) By considering all molecular modes combined with a contribution from the salt substrate the states | 5 , D σ xx ⟩ become the ground state of the dianionic molecule.

![image](https://ac.marcodsn.me/arxiv-markdown-images/2408.00478v1-figure-f51afa75-1254-4bbc-b29d-f0dfccbf09a8.jpg)

## Microscopic model of the salt electron-phonon coupling

We present in this section a simple microscopic model to estimate the electron-phonon coupling between an excess charge on the molecule and surface phonons of the underlying salt. We consider the bare Coulomb interaction between the excess electrons in the molecule and the salt ions. Comparison with experimental results indicate that this method seems to overestimate the values for the coupling.

## Linear vibronic coupling values

We consider vibration of the ions in all three space directions, thus accounting for both the longitudinal and the transversal phonon modes in the salt. For simplicity, we limit ourselves to an Einstein model of independent harmonic oscillators with the energies ¯ hω L = 32meV for the vibrations parallel to the surface and ¯ hω T = 20meV for the ones perpendicular to it. The linear electron-phonon couplings are calculated from the change of the Coulomb interaction between the LUMOs and the salt ions upon displacing of the latter. Specifically we evaluate

$$[ g _ { l } ^ { \alpha } ] _ { i j } = \frac { q _ { l } e ^ { 2 } } { 4 \pi \epsilon _ { 0 } } \left ( \frac { \partial } { \partial u _ { l \alpha } } \int \, d r \frac { \phi _ { i } ^ { * } ( r ) \phi _ { j } ( r ) \Delta _ { 0 l } ^ { \alpha } } { | r - R _ { l } ^ { 0 } - u _ { l } | } \right ) _ { \{ u _ { l \alpha } \} = 0 }.$$

Here l labels an atom in the salt, i, j = L ± specify one of the degenerate LUMOs and α = x, y or z represents the spatial direction of the displacement. Furthermore, we have q i = ± 1 depending on the atom being a chloride (+1) or sodium (-1) atom, e is the elementary charge and ϵ 0 the vacuum permittivity. The atom's equilibrium position is specified by R 0 l , u l is the displacement coordinate of the considered atom and we have introduced the zero point fluctuations ∆ x/y 0 l = √ ¯ h 2 m ω l L and ∆ z 0 l = √ ¯ h 2 m ω l T . The central copper atom of the molecule, which we define to be at the origin, is centered on a chloride ion and we take 20 layers of salt into account. We checked that the integral in Eq. (49) is converged in the z direction (perpendicular to the salt surface) when taking into account this number of layers. In the x and y direction we consider a square with an area of 400 ˚ A , which again suffices for the convergence of 2 the electron phonon couplings. Moreover, we choose a distance of 2 8 ˚ between salt and molecule, which is consistent . A

with values found in the literature [22, 31]. By calculating the derivative in Eq. (49), we obtain the expression

$$[ g _ { l } ^ { \alpha } ] _ { i j } = \frac { q _ { l } e ^ { 2 } } { 4 \pi \epsilon _ { 0 } } \int \, d x d y d z \frac { \phi _ { i } ^ { * } ( x, y, z ) \phi _ { j } ( x, y, z ) ( \alpha - R _ { l \alpha } ^ { 0 } ) \Delta _ { 0 l } ^ { \alpha } } { \left [ ( x - R _ { l x } ^ { 0 } ) ^ { 2 } + ( y - R _ { l y } ^ { 0 } ) ^ { 2 } + ( z - R _ { l z } ^ { 0 } ) ^ { 2 } \right ] ^ { \frac { 3 } { 2 } } } \,,$$

which represents the starting point for the numerical calculation of the coupling constants.

## Numerical evaluation of the integral

The numerical integration of Eq. (50) faces the problem of a divergent integrand when the integration variable coincides with one of the equilibrium positions of the atoms in the salt r = R l . The integral is though not diverging. To better understand this issue, we expand the product of the wave functions around the diverging point as

$$\phi _ { i } ^ { * } ( r ) \phi _ { j } ( r ) = & \phi _ { i } ^ { * } ( R _ { i } ^ { 0 } ) \phi _ { j } ( R _ { i } ^ { 0 } ) + \sum _ { \alpha } \left ( \frac { \partial } { \partial r _ { \alpha } } \phi _ { i } ^ { * } ( r ) \phi _ { j } ( r ) \right ) _ { r = R _ { i } ^ { 0 } } \left ( r _ { \alpha } - R _ { i \alpha } ^ { 0 } \right ) + \\ & \frac { 1 } { 2 } \sum _ { \alpha } \sum _ { \beta } \left ( \frac { \partial ^ { 2 } } { \partial r _ { \alpha } \partial r _ { \beta } } \phi _ { i } ^ { * } ( r ) \phi _ { j } ( r ) \right ) _ { r = R _ { i } ^ { 0 } } \left ( r _ { \alpha } - R _ { i \alpha } ^ { 0 } \right ) \left ( r _ { \beta } - R _ { i \beta } ^ { 0 } \right ) + \dots,$$

where α, β = x, y or z . We now evaluate the integral Eq. (50) for the first few terms of this expansion in the vicinity of a diverging point r = R 0 l . Furthermore, we drop all constants since we are only interested in the behavior of the integral when the integrand is diverging and without loss of generality set α = z . By shifting to cylindrical coordinates with ζ = z -R lz and ρ = √ ( x -R lx ) 2 +( y -R ly ) 2 and integrating over a small volume around the diverging point, we obtain, for the constant term in Eq. (51)

$$\int _ { 0 } ^ { 2 \pi } \, d \phi \int _ { 0 } ^ { R } \, d \rho \int _ { - \epsilon } ^ { \epsilon } \, d \zeta \frac { \rho \zeta } { ( \rho ^ { 2 } + \zeta ^ { 2 } ) ^ { \frac { 3 } { 2 } } } = 0.$$

For the linear component of Eq. (51), we distinguish two cases. The first case is for α = z , in which it is convenient to choose spherical coordinates with r = √ ρ 2 + ζ 2 and θ = arctan ( ρ ζ ) . We thus express the integral in the vicinity of the diverging point in these coordinates and obtain

$$- \int _ { 0 } ^ { 2 \pi } \, d \phi \int _ { 0 } ^ { \pi } \, d \theta \int _ { 0 } ^ { R } \, d r \, r \cos ^ { 2 } ( \theta ) \sin ( \theta ),$$

̸

which is not diverging and vanishes in the limit R → 0. The second case we need to discuss is when α = z . We can choose α = x as the case with α = y is completely analogous. Now, we rotate our coordinate system and obtain ζ ′ = y -R ly and ρ ′ = √ ( x -R lx ) 2 +( z -R lz ) 2 . With this transformation, the local integral is proportional to

$$\int _ { 0 } ^ { 2 \pi } \, d \phi \int _ { 0 } ^ { R } \, d \rho ^ { \prime } \int _ { - \epsilon } ^ { \epsilon } \, d \zeta ^ { \prime } \frac { \rho ^ { \prime 3 } \cos ( \phi ) \sin ( \phi ) } { ( \rho ^ { \prime 2 } + \zeta ^ { \prime 2 } ) ^ { \frac { 3 } { 2 } } } = 0.$$

Analogously, one can prove that also the integrals involving higher orders expansions will be finite and vanish in the limit of infinitesimal integration volume around the point with diverging integrand. Therefore, in the numerical evaluation we can set the integrand at the diverging point to 0 and can just use standard numerical integration schemes. The electron-phonon coupling relative to each vibrational mode of a salt ion is a Hermitian 2x2 matrix which can be expanded in terms of the Pauli matrices. We obtain

$$g _ { l } ^ { \alpha } = g _ { l, A _ { 1 } } ^ { \alpha } \mathbf 1 _ { 2 } + g _ { l, B _ { 1 } } ^ { \alpha } \sigma ^ { x } + g _ { l, B _ { 2 } } ^ { \alpha } \sigma ^ { y }.$$

as, due to the form of the LUMO orbitals, the σ z component vanishes exactly. We labelled the different electronvibron coupling components according to the representations of the corresponding irreducible tensor operator with respect to the C 4 v point symmetry group. In fact,

$$\sum _ { \sigma } \sum _ { i j } \hat { d } _ { L i \sigma } ^ { \dagger } \hat { \delta } _ { i j } \hat { d } _ { L j \sigma } \in A _ { 1 } \,, \ \sum _ { \sigma } \sum _ { i j } \hat { d } _ { L i \sigma } ^ { \dagger } \sigma _ { i j } ^ { x } \hat { d } _ { L j \sigma } \in B _ { 1 } \,, \ \text{ and } \ \sum _ { \sigma } \sum _ { i j } \hat { d } _ { L i \sigma } \sigma _ { i j } ^ { y } \hat { d } _ { L j \sigma } \in B _ { 2 } \,.$$

So far we have shown how to calculate the electron-phonon couplings for each individual atom in the salt. We are, however, more interested in the collective behavior of the substrate. Thus, despite of the simplified Einstein model, we look for the modes with finite coupling to the excess charge in the molecular LUMOS. Analogously to the scheme proposed in Sec. we thus combine all modes with the same symmetry. We start with collective transversal modes, whose coordinates read

$$x _ { \Gamma } ^ { T } = \sum _ { l } \frac { g _ { l, \Gamma } ^ { z } } { \sqrt { \sum _ { m } \left ( g _ { m, \Gamma } ^ { z } \right ) ^ { 2 } } } \tilde { u } _ { l z },$$

being Γ = A , B , 1 1 and B 2 the corresponding irreducible representation and ˜ u lz = u lz ∆ z 0 l the displacement of each individual atom in z direction normalized by the respective zero point fluctuation. Analogously, we define for the longitudinal modes the collective mode coordinates

$$x _ { \Gamma } ^ { L } = \sum _ { l } \frac { g _ { l, \Gamma } ^ { x } \tilde { u } _ { l x } + g _ { l, \Gamma } ^ { y } \tilde { u } _ { l y } } { \sqrt { \sum _ { m } \left ( g _ { m, \Gamma } ^ { x } \right ) ^ { 2 } + \left ( g _ { m, \Gamma } ^ { y } \right ) ^ { 2 } } } \,.$$

These coordinate transformations identify six collective salt modes with a finite electron phonon coupling. All the other modes are not influenced by an excess charge on the CuPc molecule and are treated as a dissipative phonon bath. The coupled modes exhibit respectively the transversal and longitudinal coupling strengths g T Γ = √ ∑ ( l g z l, Γ ) 2

and g L Γ = √ ∑ ( l g x l, Γ ) 2 + ( g y l, Γ ) 2 , with Γ = A , B 1 1 and B 2 . Starting from the numerical calculation of the integrals in (50), one obtains the following coupling strengths for the collective modes

$$g _ { A _ { 1 } } ^ { T } & = 1 1 6 \, \text{meV}, \quad g _ { B _ { 1 } } ^ { T } & = 2 0. 6 \, \text{meV}, \quad g _ { B _ { 2 } } ^ { T } & = 2. 5 \, \text{meV} \\ g _ { A _ { 1 } } ^ { L } & = 9 2. 7 \, \text{meV}, \quad g _ { B _ { 1 } } ^ { L } & = 1 1. 1 \, \text{meV}, \quad g _ { B _ { 2 } } ^ { L } & = 1 1. 7 \, \text{meV}.$$

Similarly to what we observed with the molecular modes, also in the salt the contribution to the electron-phonon coupling of the B 1 mode dominates over the one of the B 2 modes. Therefore, we concentrate our discussion solely on the A 1 and B 1 modes. They yield a total reorganization energy

$$\Delta E _ { r e o r g } ^ { \text{salt} } = \frac { \left ( g _ { A _ { 1 } } ^ { T } \right ) ^ { 2 } + \left ( g _ { B _ { 1 } } ^ { T } \right ) ^ { 2 } } { \hbar { \omega } _ { T } } + \frac { \left ( g _ { A _ { 1 } } ^ { L } \right ) ^ { 2 } + \left ( g _ { B _ { 1 } } ^ { L } \right ) ^ { 2 } } { \hbar { \omega } _ { L } } = 9 6 6 \, \text{meV}.$$

These values of the reorganization energy due to the salt deformation can be taken as an estimate for their order of magnitude. Their quantitative accuracy is though problematic.

On the one hand, by following the conclusions from the experiments of Fatayer et al. [22] we expect that the total reorganization energy upon charging should be in the range of 300 meV ≤ ∆ E tot reorg ≤ 500 meV. As the molecular component is ∆ E mol reorg = 187meV, a salt reorganization energy in the range 113 meV ≤ ∆ E salt reorg ≤ 313 meV is to be expected. This leads us to the conclusion that our model overestimates the total reorganization energy stemming from the salt.

On the other hand, in the molecular reference condition the energy difference between the dianionic ground state (both electrons in different LUMOs) and the localized configuration is ≈ 200 meV. Thus the minimum reorganization energy stemming from the antisymmetric modes, in order to obtain the localized configuration as the dianionic state with lowest energy, must be ≈ 100 meV.

We assume that these problems can be solved by a more realistic treatment of the salt and salt-molecule interaction which includes the dispersion relation of the phonons and a screening of the bare Coulomb interaction between the molecule and the underlying ions. While keeping the form of the Hamiltonian Eq. (23), we have considered the effective couplings as free parameters and explored their effects on the many body transitions. We have chosen them to be ∆ E salt A 1 = 231 2 meV and ∆ . E salt B 1 = 80meV. Combined with the molecular contributions this choice of parameters gives a total reorganization energy of ∆ E A 1 +∆ E B 1 = 500meV. Thus not overestimating the total reorganization energy while having the localized dianionic configuration as the ground state.
Moritz Frankerl Laerte L. Patera Thomas Frederiksen Jascha Repp Andrea Donarini
Substrate stabilization of Jahn--Teller distortion in a single molecule
Charge-state transitions of a single Cu-phthalocyanine molecule adsorbed on
an insulating layer of NaCl on Cu(111) are probed by means of alternate
charging scanning tunneling microscopy. Real-space imaging of the electronic
transitions reveals the Jahn--Teller distortion occurring upon formation of the
first and second anionic charge states. The experimental findings are
rationalized by a theoretical many-body model which highlights the crucial role
played by the substrate. The latter enhances and stabilizes the intrinsic
Jahn--Teller distortion of the negatively charged molecule hosting a degenerate
pair of single-particle frontier orbitals. Consequently, two excess electrons
are found to occupy, in the ground state, the same localized orbital, despite a
larger Coulomb repulsion than the one for the competing delocalized electronic
configuration. Control over the charging sequence by varying the applied bias
voltage is also predicted.