Apparent Reversal of Molecular Orbitals Reveals Entanglement

TL;DR

This study uses scanning tunneling microscopy to observe how molecular orbital sequences change with charge state, revealing many-body entanglement effects that challenge single-particle models.

Contribution

It demonstrates the reversal of molecular orbital order due to many-body entanglement in charged molecules, providing new insights into molecular electronic structure.

Findings

Orbital sequence reverses with charge state on different substrates.

Many-body theory explains the orbital reversal through entanglement.

Experimental results challenge single-particle interpretations.

Abstract

The frontier orbital sequence of individual dicyanovinyl-substituted oligothiophene molecules is studied by means of scanning tunneling microscopy. On NaCl/Cu(111) the molecules are neutral and the two lowest unoccupied molecular states are observed in the expected order of increasing energy. On NaCl/Cu(311), where the molecules are negatively charged, the sequence of two observed molecular orbitals is reversed, such that the one with one more nodal plane appears lower in energy. These experimental results, in open contradiction with a single-particle interpretation, are explained by a many-body theory predicting a strongly entangled doubly charged ground state.

Figures

Tables

Table 1. Many-body eigenenergies $E_{N,m}$ of ${\hat{H}}_{\mathrm{mol}}$, omitting the spin degrees of freedom
$m\backslash N$	0	1	2	3	4
0	0	${ϵ}_S$	$2{ϵ}_S+\Delta +U-\sqrt{{\Delta }^2+J^2}$	$3{ϵ}_S+\Delta +3U-J$	$4{ϵ}_S+2\Delta +6U-2J$
1		${ϵ}_S+\Delta$	$2{ϵ}_S+\Delta +U-J$	$3{ϵ}_S+2\Delta +3U-J$
2			$2{ϵ}_S+\Delta +U+J$
3			$2{ϵ}_S+\Delta +U+\sqrt{{\Delta }^2+J^2}$

Equations

$$\begin{array}{rl}{H}_{\mathrm{mol}}=& {ϵ}_{\mathrm{S}}{\hat{n}}_{\mathrm{S}}+{ϵ}_{\mathrm{AS}}{\hat{n}}_{\mathrm{AS}}+\frac{U}{2}\hat{N}(\hat{N}-1)\\ +& J\sum _{\sigma {\sigma }^{\prime }}{d}_{\mathrm{AS}\sigma }^{†}{d}_{\mathrm{S}{\sigma }^{\prime }}^{†}{d}_{\mathrm{AS}{\sigma }^{\prime }}{d}_{\mathrm{S}\sigma }\\ +& J\left(d_{\mathrm{AS}\uparrow }^{†}{d}_{\mathrm{AS}\downarrow }^{†}{d}_{\mathrm{S}\downarrow }{d}_{\mathrm{S}\uparrow }+d_{\mathrm{S}\uparrow }^{†}{d}_{\mathrm{S}\downarrow }^{†}{d}_{\mathrm{AS}\downarrow }{d}_{\mathrm{AS}\uparrow }\right),\end{array}$$

$$\begin{array}{rl}|2,0⟩& =\cos \theta d_{\mathrm{S}\uparrow }^{†}{d}_{\mathrm{S}\downarrow }^{†}|\Omega ⟩+\sin \theta d_{\mathrm{AS}\uparrow }^{†}{d}_{\mathrm{AS}\downarrow }^{†}|\Omega ⟩,\end{array}$$

$${\hat{H}}_{\mathrm{G}}={\hat{H}}_{\mathrm{mol}}-\delta {\hat{N}}^2+{\varphi }_0\hat{N}$$

$${\hat{H}}_{\mathrm{mol}}=$$

$$+$$

$$\begin{array}{rl}{U}_{\mathrm{S}}& =\frac{e^2}{4\pi {ϵ}_r{ϵ}_0}\int \mathrm{d}{\mathbf{r}}_1\int \mathrm{d}{\mathbf{r}}_2\frac{|{\psi }_{\mathrm{S}}({\mathbf{r}}_1){|}^2|{\psi }_{\mathrm{S}}({\mathbf{r}}_2){|}^2}{|{\mathbf{r}}_1-{\mathbf{r}}_2|},\\ U_{\mathrm{SAS}}& =\frac{e^2}{4\pi {ϵ}_r{ϵ}_0}\int \mathrm{d}{\mathbf{r}}_1\int \mathrm{d}{\mathbf{r}}_2\frac{|{\psi }_{\mathrm{S}}({\mathbf{r}}_1){|}^2|{\psi }_{\mathrm{AS}}({\mathbf{r}}_2){|}^2}{|{\mathbf{r}}_1-{\mathbf{r}}_2|},\\ U_{\mathrm{AS}}& =\frac{e^2}{4\pi {ϵ}_r{ϵ}_0}\int \mathrm{d}{\mathbf{r}}_1\int \mathrm{d}{\mathbf{r}}_2\frac{|{\psi }_{\mathrm{AS}}({\mathbf{r}}_1){|}^2|{\psi }_{\mathrm{AS}}({\mathbf{r}}_2){|}^2}{|{\mathbf{r}}_1-{\mathbf{r}}_2|},\end{array}$$

$$J=\frac{e^2}{4\pi {ϵ}_r{ϵ}_0}\int \mathrm{d}{\mathbf{r}}_1\int \mathrm{d}{\mathbf{r}}_2\frac{{\psi }_{\mathrm{AS}}^{\ast }({\mathbf{r}}_1){\psi }_{\mathrm{S}}^{\ast }({\mathbf{r}}_2){\psi }_{\mathrm{AS}}({\mathbf{r}}_2){\psi }_{\mathrm{S}}({\mathbf{r}}_1)}{|{\mathbf{r}}_1-{\mathbf{r}}_2|},$$

$$e{V}_{\mathrm{b},\mathrm{res}}^{N,m\leftrightarrow N+1,n}=\frac{1}{{\alpha }_T}(E_{N+1,n}^{\mathrm{G}}-E_{N,m}^{\mathrm{G}}),$$

$$\begin{array}{rl}e{V}_{111}^{>0}& ={\alpha }_{\mathrm{T}}^{-1}\left({ϵ}_{\mathrm{S}}-\delta +{\varphi }_{111}\right),\\ e{V}_{111}^{>1}& ={\alpha }_{\mathrm{T}}^{-1}\left({ϵ}_{\mathrm{S}}+\Delta -\delta +{\varphi }_{111}\right).\end{array}$$

$$\begin{array}{rl}e{V}_{311}^{<0}& ={\alpha }_{\mathrm{T}}^{-1}\left({ϵ}_{\mathrm{S}}+\Delta +U-\sqrt{{\Delta }^2+J^2}-3\delta +{\varphi }_{311}\right),\\ e{V}_{311}^{>0}& ={\alpha }_{\mathrm{T}}^{-1}\left({ϵ}_{\mathrm{S}}+2U-J+\sqrt{{\Delta }^2+J^2}-5\delta +{\varphi }_{311}\right),\\ e{V}_{311}^{>1}& ={\alpha }_{\mathrm{T}}^{-1}\left({ϵ}_{\mathrm{S}}+\Delta +2U-J+\sqrt{{\Delta }^2+J^2}-5\delta +{\varphi }_{311}\right).\end{array}$$

$${\hat{H}}_{\eta =\mathrm{S},\mathrm{T}}=\sum _{\mathbf{k}\sigma }{ϵ}_{\eta \mathbf{k}}{\hat{c}}_{\eta \mathbf{k}\sigma }^{†}{\hat{c}}_{\eta \mathbf{k}\sigma }^{},$$

$${\hat{H}}_{\mathrm{tun}}=\sum _{\eta \mathbf{k}i\sigma }{t}_{\mathbf{k}i}^{\eta }{\hat{c}}_{\eta \mathbf{k}\sigma }^{†}{\hat{d}}_{i\sigma }^{}+h.c.,$$

$$\begin{array}{rl}{R}_{\sigma \eta }^{N,n\to N+1,m}& =\sum _{ij}{\Gamma }_{ji}^{\eta }(E_{N+1,m}^G-E_{N,n}^G)\\ & ⟨N+1,m|{\hat{d}}_{i\sigma }^{†}|N,n⟩⟨N,n|{\hat{d}}_{j\sigma }^{}|N+1,m⟩f^{+}(E_{N+1,m}^G-E_{N,n}^G-{\alpha }_{\eta }e{V}_{\mathrm{bias}},T)\\ R_{\sigma \eta }^{N,n\to N-1,m}& =\sum _{ij}{\Gamma }_{ij}^{\eta }(E_{N-1,m}^G-E_{N,n}^G)\\ & ⟨N-1,m|{\hat{d}}_{i\sigma }^{}|N,n⟩⟨N,n|{\hat{d}}_{j\sigma }^{†}|N-1,m⟩f^{-}(E_{N,n}^G-E_{N-1,m}^G-{\alpha }_{\eta }e{V}_{\mathrm{bias}},T),\end{array}$$

$${\dot{\rho }}_{\mathrm{red}}=\mathcal{L}[{\rho }_{\mathrm{red}}],$$

$$\mathcal{L}={\mathcal{L}}_{\mathrm{S}}+{\mathcal{L}}_{\mathrm{T}}+{\mathcal{L}}_{\mathrm{rel}}$$

$${\mathcal{L}}_{\mathrm{rel}}\left[\rho \right]=-\frac{1}{\tau }\left(\rho -\sum _N{\rho }^{\mathrm{th},N}\sum _l{\rho }_{ll}^N\right).$$

$${\rho }^{\mathrm{th},\mathrm{N}}=\sum _k\frac{{\mathrm{e}}^{-\beta E_{Nk}}}{{\sum }_l{\mathrm{e}}^{-\beta E_{Nl}}}|Nk⟩⟨Nk|,$$

$$⟨{\hat{I}}_{\eta }⟩={\mathrm{Tr}}_{\mathrm{mol}}\left(\hat{N}{\mathcal{L}}_{\eta }[{\rho }_{\mathrm{red}}^{\infty }]\right),$$

Full text

For the use of single molecules as devices, engineering and control of their intrinsic electronic properties is all-important. In this context, quantum effects such as electronic interference have recently shifted into the focus Cardamone et al. (2006); Donarini et al. (2009, 2010); Guédon et al. (2012); Vazquez et al. (2012); Ballmann et al. (2012); Xia et al. (2014). Most intriguing in this respect are electron correlation effects Zhao et al. (2005); Maruccio et al. (2007); Fernández-Torrente et al. (2008); Franke et al. (2011); Chiesa et al. (2013); Grothe et al. (2013); Ervasti et al. (2016), which are intrinsically strong in molecules due to their small size Begemann et al. (2008); Toroz et al. (2011, 2013); Schulz et al. (2015); Siegert et al. (2016).

In general, Coulomb charging energies strongly depend on the localization of electrons and hence on the spatial extent of the orbitals they occupy. Therefore the orbital sequence of a given molecule can reverse upon electron attachment or removal, if some of the frontier orbitals are strongly localized while others are not, like in e. g. phthalocyanines Liao and Scheiner (2001); Nguyen and Pachter (2003); Wu et al. (2006, 2008); Uhlmann et al. (2013). Coulomb interaction may also lead to much more complex manifestations such as quantum entanglement of delocalized molecular orbitals.

Here we show, that the energy spacing of the frontier orbitals in a single molecular wire of individual dicyanovinyl-substituted quinquethiophene (DCV5T) can be engineered to achieve near-degeneracy of the two lowest lying unoccupied molecular orbitals, leading to a strongly-entangled ground state of DCV5T2-. These orbitals are the lowest two of a set of particle-in-a-box-like states and differ only by one additional nodal plane across the center of the wire. Hence, according to the fundamental oscillation theorem of the Sturm-Liouville theory their sequence has to be set with increasing number of nodal planes, which is one of the basic principles of quantum mechanics Hückel (1931); Barford (2005). This is evidenced and visualized from scanning tunneling microscopy (STM) and spectroscopy (STS) of DCV5T on ultrathin insulating films. Upon lowering the substrate’s work function, the molecule becomes charged, leading to a reversal of the sequence of the two orbitals. The fundamental oscillation theorem seems strikingly violated since the state with one more nodal plane appears lower in energy. This contradiction can be solved, though, by considering intramolecular correlation leading to a strong entanglement in the ground state of DCV5T2-.

The experiments were carried out with a home-built combined STM/atomic force microscopy (AFM) using a qPlus sensor Giessibl (2000) operated in ultra-high vacuumat a temperature of $6$ K. Bias voltages are applied to the sample. All AFM data, d$I$/d$V$ spectra and maps were acquired in constant-height mode. Calculations of the orbitals and effective single particle electronic structure were performed within the density functional theory (DFT) as implemented in the SIESTA code Soler et al. (2002) and are based on the generalized gradient approximation (GGA-PBE). The many-body eigenstates are determined from a diagonalization of the many-body model Hamiltonian $H_{\rm mol}$, which is defined further below in the main text. Based on these, STM-image and spectra simulations were performed within a Liouville approach for the density matrix $\rho$. See Supplemental Material Sup for more details.

The molecular structure of DCV5T, shown in Fig. 1a, consists of a quinquethiophene (5T) backbone and a dicyanovinyl (DCV) moiety at each end. The delocalized electronic system of polythiophene and oligo-thiophene enables conductance of this material Feast et al. (1996); Yamada et al. (2008); Fitzner et al. (2011). The lowest unoccupied orbital of each of the thiophene rings couples electronically to its neighbors and forms a set of particle-in-a-box-like states Repp et al. (2010); Kislitsyn et al. (2016). The LUMO to LUMO+1 level spacing of the quinquethiophene (5T) backbone is approx. 0.7 eV Repp et al. (2010), which is in good agreement with the energy difference calculated for free 5T based on DFT, as shown in Fig. 1a, left. This DFT-based calculation also confirms the nature of the LUMO and LUMO+1 orbitals, both deriving from the single thiophene’s LUMOs and essentially differing only by one additional nodal plane across the center of the molecule. To enable the emergence of correlation and thus level reordering, we have to bring these two states closer to each other. This is achieved by substituting dicyanovinyl moieties with larger electron affinity at each end of the molecular wire. As the orbital density of the higher lying particle-in-a-box-like state, namely LUMO+1, has more weight at the ends of the molecule, it is more affected by this substitution than the lowest state, the LUMO. This is evidenced by corresponding calculations of DCV5T, for which the LUMO to LUMO+1 energy difference is reduced by more than a factor of two, see Fig. 1a, left. The increased size of DCV5T may also contribute to the reduced level spacing. For the rest of this work, we concentrate on the LUMO and LUMO+1 orbitals only. To avoid confusion, we refrain from labeling the orbitals according to their sequence but instead according to their symmetry with respect to the mirror plane perpendicular to the molecular axis, as symmetric (S) and antisymmetric (AS). Hence, the former LUMO and the LUMO+1 are the S and AS states, respectively.

To study the energetic alignment of the orbitals as well as their distribution in real space, we employ ultrathin NaCl insulating films to electronically decouple the molecules from the conductive substrate Repp et al. (2005). It has been previously shown that in these systems the work function can be changed by using different surface orientations of the underlying metal support Repp et al. (2005); Olsson et al. (2007); Swart et al. (2011). Importantly, this does not affect the (100)-terminated surface orientation of the NaCl film, such that the local chemical environment of the molecule remains the same, except for the change of the work function.

However, in the present case, this alone has a dramatic effect on the electronic structure of the molecular wires as is evidenced in Fig. 1b. There, the STM images are shown for voltages corresponding to the respective lowest lying molecular resonances at positive sample voltage for DCV5T adsorbed on NaCl/Cu(111) (top panel) and NaCl/Cu(311) (bottom panel). They both show a hot-dog like appearance of the overall orbital density as was observed and discussed previously Repp et al. (2010); Bogner et al. (2015). Importantly in the current context, however, the orbital density of DCV5T/NaCl/Cu(311) shows a clear depression at the center of the molecule, indicating a nodal plane, whereas DCV5T/NaCl/Cu(111) does not. Apparently, the energetically lowest lying state is not the same for the two cases, but S for DCV5T/NaCl/Cu(111) and AS in the case of DCV5T/NaCl/Cu(311). In contrast, STM images acquired at voltages well below the first resonance reflect the geometry of the molecule in both cases as wire-like protrusion (see insets of Fig. 1b).

We hence assume that the molecules are neutral on NaCl/Cu(111) and that the S state corresponds to the LUMO. According to the literature, changing the copper surface orientation from Cu(111) to Cu(311) results in a lowering of the work function by approximately $1$ eV  Gartland et al. (1972); Repp et al. (2005); Olsson et al. (2005). Hence, one may expect that the former LUMO, initially located $0.7$ eV above the Fermi level $E_{F}$ in the case of NaCl/Cu(111) will shift to below the Fermi level Swart et al. (2011); Uhlmann et al. (2013) for NaCl/Cu(311) such that the molecule becomes permanently charged.

To obtain a systematic understanding of the level alignment of the S and AS states of the molecule on both substrates, we acquired differential conductance (d$I$/d$V$) spectra and d$I$/d$V$-maps on DCV5T molecules. Typical spectra measured at the center and the side of the molecule are shown in Figs. 2a and b on NaCl/Cu(111) and NaCl/Cu(311), respectively. DCV5T exhibits two d$I$/d$V$ resonances at positive bias but none at negative voltages down to -2.5 V. According to the d$I$/d$V$ maps and consistent with the different intensities in the spectra acquired on and off center of the molecule, the S state at $\simeq 0.7$ V is lower in energy than the AS state occurring at $\simeq 1.1$ V. The energy difference of $\simeq 0.4$ eV is in rough agreement to our calculations (see Fig. 1a). As discussed above, in the case of NaCl/Cu(311), DCV5T exhibits the AS state as the lowest resonance at positive bias voltages, this time at $\simeq 0.9$ V. This is additionally evidenced by the constant-current STM image and the corresponding d$I$/d$V$ map in Fig. 2b. The S state is now located at higher voltages, namely at $\simeq 1.3$ V, as seen in the spectrum and the d$I$/d$V$ map. Obviously, the two states are reversed in their sequence. In this case, at negative bias voltages, a peak in d$I$/d$V$ indicates an occupied state in equilibrium, in stark contrast to DCV5T/NaCl/Cu(111) but in agreement with the assumption of the molecule being negatively charged. The constant-current image acquired at $-0.7$ V, corresponding to the first peak at negative bias, seems to be a superposition of both the S and AS states.

The experimentally observed reversal of the orbital sequence is in striking disagreement with the fundamental oscillation theorem. To understand this apparent orbital reversal we go beyond the single particle picture and invoke the role of electronic correlations. In the double-barrier tunneling junction geometry employed here, the resonances in d$I$/d$V$ are associated with a temporary change of electron number on the molecule. In this terms the two peaks of DCV5T/NaCl/Cu(111) at positive bias are DCV5T$\leftrightarrow$DCV5T- transitions (See Fig. 3), and, in the same spirit, the ones of DCV5T/NaCl/Cu(311) at positive and at negative bias should be interpreted as DCV5T${}^{2-}\leftrightarrow$DCV5T3- and DCV5T${}^{2-}\leftrightarrow$DCV5T*-* transitions, respectively.

Both the topographical and the spectroscopic data presented so far suggest that the electronic transport through DCV5T involves, in the present bias and work function ranges, only the symmetric (S) and the antisymmetric (AS) orbitals. We concentrate on them and freeze the occupation of the other lower (higher) energy orbitals to 2 (0). In terms of these S and AS frontier orbitals we write the minimal interacting Hamiltonian for the isolated molecule,

[TABLE]

where $d^{\dagger}_{\mathrm{S(AS)}\sigma}$ creates an electron with spin $\sigma$ in the symmetric (antisymmetric) orbital, $\hat{n}_{\mathrm{i}}$ counts the number of electron in the orbital with ${\mathrm{i}}={\mathrm{S,AS}}$ and $\hat{N}$ represents the total number of electrons occupying the two frontier orbitals. The interaction parameters $U=1.4$ eV and $J=0.75$ eV are obtained from the DFT orbitals by direct calculation of the associated Coulomb integrals and assuming a dielectric constant $\epsilon_{r}=2$ which accounts for the screening introduced by the underlying frozen orbitals Ryndyk et al. (2013); Siegert et al. (2016). As expected from their similar (de-)localization, the Coulomb integrals of the S and AS states are almost identical 111For the Coulomb integrals we obtain $U_{\mathrm{S-S}}=1.37\,$eV, $U_{\mathrm{AS-AS}}=1.43\,$eV, $U_{\mathrm{S-AS}}=1.37\,$eV. Besides a constant interaction charging energy $U$, the model defined in Eq. (1) contains exchange interaction and pair-hopping terms, both proportional to $J$, which are responsible for the electronic correlation. The electrostatic interaction with the substrate is known to stabilize charges on atoms and molecules Kaasbjerg and Flensberg (2008, 2011); Olsson et al. (2007) due to image charge and polaron formation. We account for this stabilization with the additional Hamiltonian $H_{\mathrm{mol-env}}=-\delta\hat{N}^{2}$. The orbital energies $\epsilon_{\mathrm{S}}=-3.1\,$eV and $\epsilon_{\mathrm{AS}}=-2.8\,$eV as well as the image-charge renormalization $\delta=0.43\,$eV are obtained from the experimental resonances of the neutral molecule and previous experimental results on other molecules Sup

Many-body interaction manifests itself most strikingly for the ground state DCV5T2-, which will therefore be discussed at first. Consider the two many-body states, in which the two extra electrons both occupy either the S or the AS state: They differ in energy by the energy $2\Delta$, where $\Delta=\epsilon_{\mathrm{AS}}-\epsilon_{\mathrm{S}}$ is the single-particle level spacing between the S and the AS state. These two many-body states interact via pair-hopping of strength $J$, leading to a level repulsion. As long as $\Delta\gg J$, this effect is negligible. In DCV5T, though, the single-particle level spacing $\Delta$ is small compared to the pair-hopping $J$, leading to an entangled ground state of DCV5T2- as

[TABLE]

with $\quad\theta=\frac{\arctan\left(J/\Delta\right)}{2}$ and where $|\Omega\rangle$ is the ground state of neutral DCV5T. Note that here, as $J/\Delta\approx 2.6$, this state shows more than 30% contribution from both constituent states, is strongly entangled, and therefore it can not be approximated by a single Slater determinant. The first excited state of DCV5T2- is a triplet with one electron in the S and one in the AS orbital at about $54\,$meV above the ground state, as shown in Fig. 3.

The level repulsion in DCV5T2- mentioned above leads to a significant reduction of the ground state energy by roughly 0.5 eV. This effect enhances the stability of the doubly charged molecule to the disadvantage of DCV5T*-*, which has just a single extra electron and therefore does not feature many-body effects.

Within the framework of the many-body theory, as sketched in Fig. 3, the apparent orbital reversal between Fig. 2a and Fig. 2b is naturally explained. To this end, as mentioned above, tunneling events in the STM experiments have to be considered as transitions between the many-body states of different charges $N$ (see arrows in Fig. 3). The spatial fingerprints of the transitions and hence their appearance in STM images is given by the orbital occupation difference between the two many-body states and is indicated by the labels S and AS in Fig. 3.

When on NaCl(2ML)/Cu(111), the DCV5T molecule is in its neutral ground state, see green panel in Fig. 3.A sufficiently large positive sample bias triggers transitions to the singly charged DCV5T-: The S and AS transitions subsequently become energetically available in the expected order of the corresponding single-particle states. A fast tunnelling of the extra electron to the substrate restores the initial condition enabling a steady-state current.

When on NaCl(2ML)/Cu(311) the molecule is doubly charged and in the entangled ground state described by Eq. (2), see Fig. 3. At sufficiently high positive sample bias the transitions to DCV5T3- are opening, enabling electron tunnelling from the tip to the molecule. The topography of these transitions is again obtained by comparing the 2 and the 3 (excess) electron states of DCV5T (cf. Fig. 3). The transition to the 3 particle ground state occurs by the population of the AS state and it involves the first component of the entangled 2 electron ground state only. The second component cannot contribute to this transition, which is bound to involve only a single electron tunneling event. Correspondingly, at a larger bias the first excited 3 particle state becomes accessible, via a transition involving the second component of the 2 particle ground state only. This transition has a characteristic S state topography. Hence, although the electronic structure of the 3 electron states does follow the Aufbau principle, the entanglement of the 2 particle ground state leads to the apparent reversal of the orbital sequence.

As described in the SI, in addition to the many-body spectrum we calculated the full dynamics of subsequent tunneling processes for all relevant situations, resulting in the calculated d$I$/d$V$ characteristics, constant current maps and constant height d$I$/d$V$ maps for a DCV5T single molecule junction presented in Fig. 4. A qualitative agreement with the experimental results of Fig. 2 can be observed both for the relative strength of the spectral peaks and the dI/dV maps. The above discussed apparent orbital reversal is fully consistent with the calculations.

The experimental data of DCV5T on the Cu(311) substrate at negative bias also show a non-standard feature. The d$I$/d$V$ map at resonance resembles a superposition of the S and AS orbital, see Fig. 2b. The effect is also reproduced in the theoretical simulations presented in Fig. 4. This can be rationalized in terms of a non-equilibrium dynamics associated to a population inversion recently predicted by some of the authors Siegert et al. (2016).

In conclusion, we showed that a reduction of the single-particle level spacing of two frontier orbitals enables the manifestation of strong electron-correlation effects in single molecules. Here, the single-particle level spacing engineered by dicyanovinyl-substitution is leading to an apparent reversal of orbital sequence and a strongly-entangled ground state of DCV5T2-. The many body description of the electronic transport is capable to reconcile the experimental observations of the orbital reversal with the fundamental oscillation theorem of quantum mechanics and shows how to achieve quantum entanglement of frontier orbitals in molecules.

Appendix A Parametrization of the many-body Hamiltonian

The grandcanonical many-body Hamiltonian used in this work,

[TABLE]

where

[TABLE]

is characterized by six parameters: i.e. the single particle energies $\epsilon_{\mathrm{S}}$ and $\epsilon_{\mathrm{AS}}$ of the frontier orbitals, the direct interaction and exchange integrals $U$ and $J$, the image charge and polaron renormalization energy $\delta$, the substrate work function $\phi_{0}$. Some of these parameters ($U$ and $J$) are obtained from first principle calculations, others ($\epsilon_{\mathrm{S}}$ and $\epsilon_{\mathrm{AS}}$ and $\delta$) are fitted to the present experimental data or taken from the literature ($\phi_{0}$) Olsson et al. (2005).

The direct interaction parameter results from a simplification of a more general model in which all possible combinations of density-density interaction including the symmetric and the antisymmetric orbital are taken into account. In the most general case one should consider the three parameters:

[TABLE]

where the screening introduced by $\epsilon_{r}=2$ is justified, even for an isolated molecule, by the presence of a polarizable core electrons. We have performed the integrals in Eq. (5) using the DFT wave functions plotted in Fig. 1 of the main text, with the help of a Montecarlo method. The very similar numerical results ($U_{\rm S}=1.37e$V, $U_{\rm AS}=1.43e$V, and $U_{\rm SAS}=1.37e$V), together with the partial convergence of the Montecarlo method when applied to the DFT molecular orbitals suggested us to simplify the model to a single parameter $U=1.4e$V. The robustness of the level ordening in the many-body spectrum with respect to variations of the interaction parameters within the estimated error given by the Montecarlo method has been tested. The exchange parameter $J$ has been similarly calculated from the formula:

[TABLE]

which yields $J=0.75e$V. Since the molecular orbitals are real functions, Eq. (6) also gives the pair hopping integral.

The remaining parameters ($\epsilon_{\mathrm{S}}$ and $\epsilon_{\mathrm{AS}}$ and $\delta$) are fitted to experimental data from the d$I$/d$V$ measurements. The position of the resonant peaks in the differential conductance curves corresponding to the transition $|N,m\rangle\leftrightarrow|N+1,n\rangle$ is given by the relation:

[TABLE]

where $\alpha_{\mathrm{T}}$ ,the voltage drop at the tip, is defined by $\mu_{T}=\mu_{0}+\alpha_{\mathrm{T}}eV_{\rm b}$, and for the the grand canonical energies $E^{\mathrm{G}}_{N,m}$ it holds $E^{\mathrm{G}}_{N,m}=E_{N,m}-\delta N^{2}+\phi_{0}N$, where $E_{N,m}$ is the eigenvalue of the $m$-th excited $N$ particle state of $\hat{H}_{\mathrm{mol}}$ (analytical expressions for $E_{N,m}$ are given in Table 1. In the experimental data (Fig. 3 in the main text) we identify 5 resonance biases $V_{\rm b,res}$ from which, besides $\epsilon_{\mathrm{S}}$ and $\epsilon_{\mathrm{AS}}$ and $\delta$, also the parameter $\alpha_{T}$ can be extracted.

We assign the resonances seen in the measurements on Cu(111), $\phi_{111}=4\,e$V from low to high bias to specific transitions between 0 and 1 particle states: $V_{111}^{>0}=V_{\mathrm{b,res}}^{0\leftrightarrow 10}=0.7$V and $V_{111}^{>1}=V_{\mathrm{b,res}}^{0\leftrightarrow 11}=1.1$V, yielding

[TABLE]

where $\Delta=\epsilon_{\mathrm{AS}}-\epsilon_{\mathrm{S}}$. Analogously, we assign the resonances in the d$I$/d$V$ measurements on Cu(311), $\phi_{311}=3\,$eV. The only negative bias resonance is associated to a transition between 1 and 2 particle states: $V_{311}^{<0}=V_{\mathrm{b,res}}^{1\leftrightarrow 20}=-0.7$V. The ones at positive bias involve instead 2 and 3 particle states: $V_{311}^{>0}=V_{\mathrm{b,res}}^{2\leftrightarrow 30}=0.9$V, and $V_{311}^{>1}=V_{\mathrm{b,res}}^{2\leftrightarrow 31}=1.3$V, we get

[TABLE]

where $e$ is the elementary charge taken with positive sign. It is now straightforward to determine the bias drop and the parameters of the Hamiltonian: $\alpha_{\mathrm{T}}=0.70$, $\delta=0.43\,\mathrm{eV}$, $\epsilon_{\mathrm{S}}=-3.08\,\mathrm{eV}$, and $\epsilon_{\mathrm{AS}}=-2.8\,\mathrm{eV}$.

Appendix B Dynamics and transport

The transport characteristics for the STM single molecule junction with thin insulating film, are obtained following the approach already introduced by some of the authors Sobczyk et al. (2012); Donarini et al. (2012); Siegert et al. (2013) in earlier works. We summarize here only the main steps of the calculation. The junction is described by the Hamiltonian $\hat{H}_{\mathrm{}}=\hat{H}_{\mathrm{G}}+\hat{H}_{\mathrm{S}}+\hat{H}_{\mathrm{T}}+\hat{H}_{\mathrm{tun}}$, where, beside the grand canonical Hamiltonian $\hat{H}_{\mathrm{G}}$ for the molecule, $\hat{H}_{\mathrm{S}}$ and $\hat{H}_{\mathrm{T}}$ correspond to substrate (S) and tip (T), respectively and $\hat{H}_{\mathrm{tun}}$ contains the tunnelling dynamics. The tip and the substrate are treated as noninteracting electronic leads:

[TABLE]

where $\hat{\operatorname{c}}_{\eta\mathbf{k}\sigma}^{{\dagger}}$ creates an electron in lead $\eta$ with spin $\sigma$ and momentum $\mathbf{k}$. The tunneling Hamiltonian $\hat{H}_{\mathrm{tun}}$ is given by

[TABLE]

and it contains the tunneling matrix elements $t^{\eta}_{\mathbf{k}i}$, which are obtained by calculating the overlap between the lead wavefunctions $\ket{\eta\mathbf{k}}$ and the molecular orbitals $\ket{i}$ Sobczyk et al. (2012). The latter are the starting point for the calculation of the single particle tunnelling rate matrices $\Gamma^{\eta}_{ij}(E)=\frac{2\pi}{\hbar}\sum_{\mathbf{k}}\,t^{\eta}_{i\mathbf{k}}\left(t^{\eta}_{j\mathbf{k}}\right)^{*}\delta\left(\epsilon_{\eta\mathbf{k}}-E\right)$ and, eventually, of the many-body rates:

[TABLE]

where $f^{+}(E)=(1+\exp(\beta E))^{-1}$ is the Fermi distribution with $\beta=(k_{\mathrm{B}}T)^{-1}$ and $f^{-}(E)=1-f^{+}(E)$. Eq. (12) clearly show how each manybody rate is in general the superposition of several molecular orbitals, whose population is changed by the creation (annihilation) operator $\hat{\operatorname{d}}_{i\sigma}^{{\dagger}}$ ($\hat{\operatorname{d}}_{i\sigma}^{\phantom{{\dagger}}}$).

The system dynamics is calculated by means of the generalized master equation,

[TABLE]

for the reduced density operator Darau et al. (2009); Sobczyk et al. (2012) $\rho_{\mathrm{red}}=\operatorname{Tr}_{\mathrm{S,T}}\left(\rho\right)$. The Liouvillian superoperator in Eq. (13)

[TABLE]

contains the terms $\mathcal{L}_{\mathrm{S}}$ and $\mathcal{L}_{\mathrm{T}}$ describing tunneling from and to the substrate and the tip, respectively. These superoperators are combinationsSobczyk et al. (2012) of the many body rates in Eq. (12). To account for relaxation processes independent from the electron tunnelling, similarly to Ref. Koch and von Oppen (2005), we included the term $\mathcal{L}_{\mathrm{rel}}$:

[TABLE]

$\mathcal{L}_{\mathrm{rel}}$ combines different relaxation processes associated e.g. with the phonon emission or with particle-hole excitation in the substrate within the relaxation time approximation. This term induces the relaxation of each N-particle subblock of $\rho$ towards its (canonical) thermal distribution:

[TABLE]

with $\beta=\left(k_{\mathrm{B}}T\right)^{-1}$. The speed of the process is set by the relaxation time $\tau$. Since $\mathcal{L}_{\mathrm{rel}}$ acts separately on each $N$-particle subblock, it conserves the particle number on the molecule and thus does not contribute directly to the transport. For the calculation of the long time dynamics, we are interested in the stationary solution $\rho_{\mathrm{red}}^{\infty}$ for which $\dot{\rho}_{\mathrm{red}}^{\infty}=\mathcal{L}[\rho_{\mathrm{red}}^{\infty}]=0$. Eventually, the stationary current through the system is evaluated as

[TABLE]

being $\hat{I}_{\eta}=\hat{N}\mathcal{L_{\eta}}$ the current operator for the lead $\eta$.

Appendix C Level alignment

In previous studies of molecules on insulating films it was observed that, due to the electronic decoupling by the film, the molecular levels are roughly aligned with the vacuum level. From an electrochemical characterization the electron affinity of DCV5T in solution was determined to be at $-3.73$ eV relative to the vacuum level Fitzner et al. (2011). The polarizability of the solution lowers the electron affinity level, such that here the LUMO transport level can be expected at some tenths of an eV higher in energy. Considering the work function of NaCl/Cu(111) of about $4\,$eV Repp et al. (2005); Bennewitz et al. (1999), this expectation in good agreement with the experimentally observed position of the S state for this system.

Appendix D Kelvin probe force spectroscopy measurements

We performed Kelvin probe force spectroscopy (KPFS) measurements along the molecules for both substrates, as is shown in Fig. 5. From a fit to the parabolic shape of the frequency shift $\Delta f(V)$ as a function of sample voltage $V$, the local contact potential difference (LCPD) between tip and sampleNonnenmacher (1991); Sadewasser et al. (2009); Mohn et al. (2012); Steurer et al. (2015) is extracted. Next to the molecules, on the clean NaCl films, the LCPD differs by slightly more than $1$ eV for the two systems providing a rough estimate of the work function difference for the two systems 222For the two KPFS measurements of the two different systems the tip apex condition was not necessarily identical. in accordance with literature values Gartland et al. (1972); Repp et al. (2005); Olsson et al. (2005).

Since local surface charges and dipoles affect the LCPD above adsorbates, the latter should qualitatively reflect the charge state Gross et al. (2009), the electron affinity Schuler et al. (2014), and the charge distribution Neff and Rahe (2015); Mohn et al. (2012); Albrecht et al. (2015). The decrease of about 20 meV in LCPD over the molecule in the case of DCV5T/NaCl/Cu(111) we assume to be due to the large electron affinity of DCV5T. On the NaCl/Cu(311) substrate the observed increase of LCPD is consistent with an anionic state of DCV5T/NaCl/Cu(311) Ikeda et al. (2008).