Jing Li

Markus Holzmann

Ivan Duchemin

Xavier Blase

Valerio Olevano

Helium atom excitations by the GW and Bethe-Salpeter many-body formalism

Helium atom is the simplest many-body electronic system provided by nature. The exact solution to the Schrödinger equation is known for helium ground and excited states, and represents a workbench for any many-body methodology. Here, we check the ab initio many-body GW approximation and Bethe-Salpeter equation (BSE) against the exact solution for helium. Starting from Hartree-Fock, we show that GW and BSE yield impressively accurate results on excitation energies and oscillator strength, systematically improving time-dependent Hartree-Fock. These findings suggest that the accuracy of BSE and GW approximations is not significantly limited by self-interaction and self-screening problems even in this few electron limit. We further discuss our results in comparison to those obtained by time-dependent density-functional theory.

Introduction -

The solution of the Schrödinger equation in interacting many-body systems is a formidable problem in condensed matter and nuclear physics, as well as in quantum chemistry. Starting with the Hartree-Fock (HF) method as one of the first approaches to treat the quantum many-body problem, different theoretical formalisms have been developed over time [START_REF] Martin | Interacting Electrons[END_REF]: approaches relying on the many-body wave function, like quantum Monte Carlo or quantum chemistry methods, or on the electronic density, like densityfunctional theory (DFT) in its static or time-dependent (TDDFT) form or, finally, quantum field theoretical, Green function based methods. While exact in principle, all of these methods rely in practical calculations on approximations and recipes whose validity are difficult to judge.

The GW approximation [START_REF] Hedin | [END_REF][3][4][5][6] and Bethe-Salpeter equation (BSE) [7][8][9] are ab initio approaches to calculate the electronic structure within many-body Green function theory. The GW approximation of the self-energy was first applied to the homogeneous electron gas or jellium model [START_REF] Hedin | [END_REF] and later extended to real solids [3,5,6] and molecules [10][11][12]. The Bethe-Salpeter equation was originally developed to describe nuclear two-particle bound states [7] and later applied to calculate optical spectra in solids [8,[13][14][15] and excitations in molecules [16][17][18]. Judged by comparison with experiment [19], GW and BSE provide reasonably accurate results for electronic excitations. However, benchmarks against the experiment, as done so far for GW and BSE, are always affected by unaccounted effects not present in the theoretical description (non-Born-Oppenheimer, electronphonon, relativistic corrections, etc.), which can mask the real many-body performances of the approaches. Validation of the approximations against exact analytic solutions-or at least accurate, safe, and reliable numerical solutions-in realistic models or simple real systems is an unavoidable step for further improvement.

Neutral helium atom is the simplest real many-electron system in nature, at the extreme limit where "many electron" is reduced to only two electrons. It is the lightest atom presenting the complication of the many-body correlations, and, due to the full rotational symmetry of the ionic potential, it is simpler than the hydrogen molecule. This allows us to write the electronic wave functions on a basis of only three scalar coordinates, Hylleraas functions, whose parameters are then varied to obtain the ground or the excited states. Calculations can achieve extremely high accuracy, 10 -7 Ha on energies in the best cases [20], and their quantitative agreement with the experiment may be considered as one of the triumphs of quantum mechanics. The exact theoretical result makes helium, that is, almost a toy model and yet a real system, an ideal workbench for any many-body methodology. The purpose of this work is to check the validity of many-body GW and BSE calculations of the helium excitation spectrum against the exact results of the idealized nonrelativistic helium Hamiltonian [20]. This is a check of the core of the methodology, the solution to the many-body problem. Helium atom can represent a very severe test case for these many-body theories whose underlying approximations have been devised to describe the extreme opposite limit of many electrons, like bulk solids, where screening is much more important. The concept itself of screening is the way towards correlations beyond Hartree-Fock by the GW approximation. The validity of the BSE two-particle electron-hole propagator might be questioned in a system where the hole is dug in a Fermi sea of only two electrons. Important drawbacks of the GW and BSE approaches, like the self-interaction [21][22][23][24] or the self-screening problem [25,26], should directly manifest in helium atom calculations and strongly limit the accuracy of the results. On the other hand, the validation of GW and BSE on helium in atomic physics represents an important confirmation of this methodology; it quantifies performances and errors in condensed matter; it may face reluctance and skepticism within theoretical chemistry where GW and BSE were recently ex- ported; it may receive attention from nuclear physics for application to excitations in nuclei.

The high accuracy of experimental spectroscopy on helium have pushed theoretical calculations of higher-order effects, like nucleus finite mass recoil, relativistic fine structure, and QED radiative corrections [27][28][29]. Furthermore, the electron-hole (e-h) Bethe-Salpeter calculation of the present paper should not be confused with the electron-electron extended, three-body external potential, relativistic Bethe-Salpeter calculations [29][30][31] of higher-order corrections to the fine structure. Although an active research field, these corrections are not the purpose of the present work. Unless otherwise stated, we use atomic units in the following.

Calculation -The many-body perturbation theory methodology is a three step procedure. (1) First, the calculation of a starting, approximated electronic structure (both energies and wave functions) of a somewhat arbitrary Hamiltonian chosen as zero order. Today, the most common choice is DFT, e.g., in the local-density approximation (LDA) or beyond. Here, we preferred, rather, a return to the origins [3,4,8] and chose Hartree-Fock. This choice avoids hybridization with alternative theories, like DFT, so to allow, at the end, a clear comparison between GW -BSE and DFT-TDDFT results. In fact, the accuracy of DFT-TDDFT results on He critically re-lies on the first DFT step, and choice of the same initial step here too can raise doubts about the effective merit of the following GW and BSE steps. We then calculate the HF electronic structure energies and wave functions, ǫ HF i , φ HF i (r),

H H (r)φ HF i (r) + dr ′ Σ x (r, r ′ )φ HF i (r ′ ) = ǫ HF i φ HF i (r),
where

H H (r) = -∂ 2 r /2 + v ext (r) + v H (r)
is the Hartree Hamiltonian and Σ x is the Fock exchange operator. We use the full nuclear potential v ext (r) = -Z/r and perform an all-electron calculation to reduce sources of inaccuracy related to pseudopotential issues in our comparison to the exact result. (2) The starting electronic structure is then used to calculate the Green function G, the dynamically screened interaction, W , in the random-phase approximation, and the self-energy, Σ, that, in the GW approximation, is the convolution product of G and W :

Σ(r, r ′ , ω) = i 2π dω ′ G(r, r ′ , ω -ω ′ )W (r, r ′ , ω ′ ). ( 1
)
The GW electronic structure, ε i , ϕ i (r), is then calculated by solving the quasiparticle equation,

H H (r)ϕ i (r) + dr ′ Σ(r, r ′ , ω = ε i )ϕ i (r ′ ) = ε i ϕ i (r). (2)
In self-consistent (SC) GW , the new electronic structure is used to iterate the calculation of G, W , and Σ Eq. ( 1) and solve Eq. ( 2), until convergence is achieved. There are several ways to perform full or "quenched" SC-GW [32][33][34][35]. The question of which is better and even whether the SC form really improves on the spectral properties [32] is still open. Here, we want to check the standard procedure, so we refer to one iteration, nonself-consistent G 0 W 0 , or apply self-consistency only on the energies ε i , keeping unchanged the wave functions ϕ i (r). Indeed, the difference between HF and exact wave functions has been found to be barely discernible for He [36]. (3) The final step is the resolution of the Bethe-Salpeter equation on top of the GW electronic structure. This is equivalently done by solving for the eigenvalues E λ and the eigenvectors Ψ λ of an excitonic Hamiltonian,

H R H C -H C * -H R * Ψ λ = E λ Ψ λ , (3)
defining the resonant (R) and coupling (C) parts as

H R vc,v ′ c ′ = (ε c -ε v )δ vv ′ δ cc ′ + iΞ vc,v ′ c ′ , H C vc,c ′ v ′ = iΞ vc,c ′ v ′
, where v and c run over, respectively, the occupied and empty GW states ε i , ϕ i (r). Also, Ξ = δΣ/δG is the Bethe-Salpeter kernel that, for singlet and triplet states, can be written as [20] and experimental (EXP) result, and DFT-exact KS energies [42].

with w(r, r ′ ) = 1/|r -r ′ | being the bare Coulomb interaction. The excitonic eigenvalues E λ of the BSE equivalent equation ( 3) are the final excitation energies of the system, while the eigenvectors Ψ λ are the excitonic wave functions from which one can extract the excitation oscillator strength. In the Tamm-Dancoff approximation (TDA), the non-Hermitian coupling part is neglected,

H C = 0.
In addition to GW and BSE, we have also performed a time-dependent Hartree-Fock (TDHF) calculation, a common approximation in quantum chemistry, condensed matter, and frequently called RPA in the context of nuclear physics [START_REF] Ring | The Nuclear Many-Body Problem[END_REF] (not to be confused with the RPA approximation in condensed matter). TDHF is equivalent to a BSE calculation started from HF instead than GW , and taking W = w (the bare Coulombian) in the BSE kernel Eqs. ( 4) and (5).

We used a d-aug-cc-pV5Z [START_REF] Woon | [END_REF] correlation-consistent Gaussian basis set with angular momentum up to l = 5 and including a double set of diffuse orbitals. HF calculations were carried by the NWCHEM package [39], and GW and BSE by the Fiesta code [10,16,40] that integrates Eq. ( 1) by contour deformation and uses a Coulombfitting resolution of the identity (RI-V) with the associated auxiliary basis d-aug-cc-pV5Z-RI [41].

Results -In Table I we report our calculated He electron removal and addition HF and GW quasiparticle (QP) energies, in comparison to the exact energies [20] which coincide with the experiment within 10 -4 Ha. For the 1s removal energy, equal to minus the ionization potential (IP), the HF result presents an error of 0.3 eV which is reduced to only 0.1 eV in GW . There are no exact calculations of the 2s addition energy, but it is known from the experiment that the He electron affinity (EA) is negative, so that the addition of a third electron costs energy and three-electron states are unbound. Both HF and GW correctly present a spectrum of unbound states (positive energy) in the electron addition part of the spectrum. This is not the case of DFT Kohn-Sham (KS) energies, even the exact ones. Exact KS eigenvalues and the exact DFT exchange-correlation potential are available in He [42,[START_REF] Umrigar | Are unoccupied Kohn-Sham eigenvalues related to excitation energies?[END_REF] by inversion of the Kohn-Sham equation imposing the exact solution [20]. The last occupied exact KS eigenvalue is known [START_REF] Perdew | [END_REF][45][46] to coincide with minus the ionization potential, and this is the case here for the 1s. However, the rest of the KS spectrum has no physical interpretation, and here we see that even exact n S L HF G0W0 GW BSETDA BSE Exact TDHF 2 3 KS energies are qualitatively off. On the other hand, exact KS energy differences were found [42,[START_REF] Umrigar | Are unoccupied Kohn-Sham eigenvalues related to excitation energies?[END_REF] surprisingly closer to the neutral (e.g., optical) excitation spectrum than KS energies are to charged (e.g., EA, photoemission) excitations. In Fig. 1 and Table II, we report the first six helium excitation energies from HF and GW energy differences, corrected by the BSE calculation and compared to the exact result [20], which also coincides with the experiment within 10 -4 Ha, and to TDHF. We see that neutral excitations are qualitatively and quantitatively captured by BSE, which reaches a good overall agreement with the exact result. The very concepts of screening, exciton and BSE electron-hole propagators are questionable and pushed to their extreme limit of application in He, and yet BSE systematically improves upon TDHF. The error is reduced by a factor of 2 or better when passing from TDHF to BSE, i.e., introducing the screening. Although one might have expected severe self-interaction and selfscreening problems in He, the error is only 0.04 eV for the first 2 3 S excitation and degrades for higher excitations, up to 0.6 eV for 3 1 S. This degradation is merely a finite basis effect. Indeed, in Table III we show the convergence of selected states with respect to the size of the Gaussian basis sets, and, in particular, with respect to the presence of diffuse (augmented) orbitals. The ionization potential IP = -ǫ 1s is converged already at the level of a standard cc-pVTZ Gaussian basis [START_REF] Woon | [END_REF] (not shown). This is because the 1s state is highly localized and requires few Gaussians to be represented accurately. Higher states get more and more delocalized and, consequently, require larger (more diffuse) basis sets. The first 2 3 S neutral excitation is converged at an augmented basis (aug-cc-pV5Z), while the intermediate 2 1 P requires a double augmentation level (d-aug-cc-pV5Z). The results for 3 1 S are clearly less converged and accurate. Higher states are not converged even at this high level of double augmentation, and we do not report them. Anyway, states towards the continuum of hydrogenic He + (1s) plus a free electron are less interesting for the study of the many-body electron-electron interaction and would require better adapted bases, e.g., plane waves. The same also holds for double excitations [47][48][49], which are much more interesting, but, for He, lie deep in the continuum [50]. Our Gaussian calculation indicates an absolute accuracy of GW and BSE calculations of 0.1 eV for the physically relevant low-lying many-body excitation spectra, independent of finite basis or pseudopotential errors.

In the following, we discuss the comparison between BSE and TDDFT results. Very accurate TDDFT results are reported in Ref. [51] (see Fig. 1 and Tables I andII). However, these results fundamentally rely on the use of unapproximated, exact-DFT Kohn-Sham eigenvalues. Exact-DFT Kohn-Sham energy differences [42,[START_REF] Umrigar | Are unoccupied Kohn-Sham eigenvalues related to excitation energies?[END_REF]51] are in He already within 0.5 eV of the exact excitation energies. Therefore, the TDDFT kernel has the easier task of only splitting singlet and triplet states, which is already well done by the adiabatic LDA standard approximation. Exact Kohn-Sham eigenvalues are exceptionally available only for He and another few systems where the exact Kohn-Sham potential is known by reverse engineering from the exact solution. If both the DFT and TDDFT calculations are consistently done using approximations, e.g., LDA, as in the more standard procedure, the results drastically change. None of the excitations are bound in DFT+TDDFT LDA [52], and GGA does not improve. Because of the missing longrange 1/r decay in the exchange-correlation potential, there is no Rydberg series in LDA or GGA atoms. These standard approximations of DFT TDDFT are not even qualitatively validated on this workbench.

The physical mechanisms behind our HF-GW -BSE scheme are different. From Fig. 1, we see that HF as well as GW energy differences, contrary to exact-DFT KS, are substantially different from the exact neutral excitation energies. This is a general feature: GW quasiparticle energies physically represent electron removal or addition charged excitations, e.g., the IP and the EA or band plots measured in photoemission. Quasiparticle energy differences are physically distinct from neutral optical excitations because they miss the electron-hole interaction. The latter is introduced only at the level of the BSE. From the difference between GW and BSE levels in Fig. 1, one can see the importance of the e-h interaction. In contrast with TDDFT, the BSE kernel manages not only to split singlet and triplet states, but, more importantly, to bring GW levels, which are placed far in the continuum, down to the discrete bound spectrum region. Our calculation has not made use of any additional helium specific knowledge from the exact solution, like the use of exact-DFT KS energies, and relies only on standard, commonly accepted approximations. Specifically, it is derived on top of the GW approxima- tion, Ξ = δΣ GW /δG, it neglects a variational term in δW/δG, and its W in Eqs. ( 4) and ( 5) is taken to be static, W (ω = 0). Therefore, the present severe workbench on the helium limiting case represents a validation of an already broadly applied recipe.

In Table II we also quote the results of BSE in the TDA approximation. In He, TDA introduces an error not larger than other approximations in this methodology. However, we already know several cases [53][54][55] where TDA breaks down. In Table II we also report intermediate G 0 W 0 results. Their values are higher than the eigenvalues of self-consistent GW by only 8∼9 ×10 -4 Ha, and this difference roughly propagates through BSE to the final excitation energies. We confirm the conclusion of Ref. [10] that, in isolated systems, there is no appreciable improvement from going beyond G 0 W 0 when HF is used as a starting point.

Finally, our work also provides the possibility of quantifying the quality of wave functions, independently from energies, through excitation oscillator strengths f [20]. In BSE, f λ results from an expression implying a sum all electron-hole electric-dipole matrix elements between QP wave functions ϕ c |e -iqr |ϕ v times the exciton wave function Ψ vc λ e-h coefficients. In the inset of Fig. 1, we plot electron and hole distributions from the excitonic wave function Ψ 2 1 P (r h , r e ) for the 2 1 P state. In Table IV we report the oscillator strength f 1 1 S→2 1 P of the transition from the ground state 1 1 S to the 2 1 P excited state. Transitions to other states studied here are all forbidden by selection rules obeyed in our electric-dipole approximation BSE code. Starting from the HF value of 0.2009, we find for BSE 0.2763, in quantitative agreement with the exact value of 0.27616 [20] and improving the TDHF of 0.2916. The exact-DFT KS result is 0.3243 [56], and TDLDA is also around that value [52]. These values confirm the accuracy of our HF-GW -BSE approach and, specifically, the quality of exciton wave functions Ψ λ .

It is worth mentioning a last interesting finding of our BSE calculation. One expects that an nL excitation has, in general, a dominant contribution from the 1s → nl, l = L transition, like the graph in Fig. 1 suggests, plus an admixture of higher energy transitions introduced by the BSE kernel Ξ. We found that this is not the case for the triplet n 3 S excitations. The 1s → 2s transition contributes only for 42% to the 2 3 S excitation, while a 52% dominant contribution comes from 1s → 3s. The reverse holds for the highest energy 3 3 S excitation, which has a 58% contribution from 1s → 2s and only 38% from the 1s → 3s transition.

Conclusions -

We have benchmarked ab initio many-body GW and BSE against exact results of the helium atom. The standard approximations behind GW and BSE accurately capture the physics of this two-body only limiting case. Excitation energies are in good agreement with the exact spectra. The surprising agreement on the f 1 1 S→2 1 P oscillator strength indicates a good description of the wave functions.

FIG. 1 .

 1 FIG.1. He atom excitation spectrum. The 0 of the energy is set to the ground state 1 1 S. Inset: 2 1 P electron-averaged hole (blue) and hole-averaged electron (red) distributions, isosurfaces taken at 5 × 10 -4 a.u.

TABLE III .

 III Basis set convercenge of the ǫ1s = -IP, selected excitation energies E n S L , and the oscillator strength f of 2 1 P .

		cc-pV5Z aug-cc-pV5Z d-aug-cc-pV5Z Exact
	ǫ1s	-0.9066	-0.9076	-0.9075 -0.9037
	E 2 3 S	0.8538	0.7284	0.7271 0.7285
	E 2 1 P	1.3345	0.8684	0.7894 0.7799
	E 3 1 S	2.6041	1.1952	0.8637 0.8425
	f 1 1 S→2 1 P	1.7607	0.7272	0.2763 0.27616

  BSE Exact TDHF HF Exact-DFT f 1 1 S→2 1 P 0.2763 0.27616 0.2916 0.2009 0.3243 f1s→2p

TABLE IV .

 IV He first excitation oscillator strengths.

Ξ ij,i ′ j ′ = i ϕ j ϕ * j ′ |W |ϕ i ϕ * i ′ -2i ϕ * i ϕ j |w|ϕ * i ′ ϕ j ′ , (4) 3 Ξ ij,i ′ j ′ = i ϕ j ϕ * j ′ |W |ϕ i ϕ * i ′ ,(5)

Acknowledgements -We thank Peter Schuck for seminal discussions and continuous support. J. L. and X. B. acknowledge funding from the European Union Horizon 2020 Programme for Research and Innovation under Grant No. 646176.


Tonatiuh Rangel

Gian-Marco Rignanese

Valerio Olevano

Can molecular projected density-of-states (PDOS) be systematically used in electronic conductance analysis?

Keywords: nano-electronics, molecular electronics, quantum transport, DFT-Landauer, benzene-diamine, benzene-dithiol

Using benzene-diamine and benzene-dithiol molecular junctions as benchmarks, we investigate the widespread analysis of the quantum transport conductance G( ) in terms of the projected density of states (PDOS) onto molecular orbitals (MOs). We first consider two different methods for identifying the relevant MOs: 1) diagonalization of the Hamiltonian of the isolated molecule, and 2) diagonalization of a submatrix of the junction Hamiltonian constructed by considering only basis elements localized on the molecule. We find that these two methods can lead to substantially different MOs and hence PDOS. Furthermore, within Method 1, the PDOS can differ depending on the isolated molecule chosen to represent the molecular junction (e.g. benzene-dithiol or -dithiolate); and, within Method 2, the PDOS depends on the chosen basis set. We show that these differences can be critical when the PDOS is used to provide an physical interpretation of the conductance (especially, when it has small values as it happens typically at zero bias). In this work, we propose a new approach trying to reconcile the two traditional methods. Though some improvements are achieved, the main problems are still unsolved. Our results raise more general questions and doubts on a PDOS-based analysis of the conductance.

LUMO LUMO+1

?

Gas-Method 1

Junction-Method 2

INTRODUCTION

According to Moore's law, in a decade or so the downscaling of conventional silicon based electronics will achieve its ultimate nanoscale limits. Molecular electronics, or electronics at the nanoscale, is considered one of the most difficult technological challenges. The construction, measurement and understanding of electronic devices constituted by single molecules in between metal electrodes, is nowadays a major concern of fundamental research.

Today, different techniques are available to realize molecular junctions in laboratories, such as electromi-gration methods, mechanical strain and scanning tunneling microscopes to open small gaps in between gold leads which can host (with a small but non-negligible probability) single molecules from a wetting solution [1][2][3] . The complete characterization of such junctions (including the measurement of their current-voltage characteristics) is however still difficult to achieve. In order to obtain reliable single-molecule zero-bias conductances, it was suggested to resort to a statistically significant sample of tens of thousands measurements 4 . Thanks to this breakthrough work, it is now possible to quote the 0bias conductance of some molecular junctions such as benzene-diamine (BDA) and benzene-dithiol (BDT) in between gold leads. Nevertheless, important character-ization uncertainties still persist. For instance, in these experiments the junction geometry is not measured and hence unknown. Given these difficulties, resorting to theory could reveal a valid approach to understand and interpret the experimental observations.

The theoretical description of the electronic quantum transport in molecular junctions or nanostructures relies on established frameworks [START_REF] Datta | Electronic transport in mesoscopic systems[END_REF][START_REF] Ventra | Electrical transport in nanoscale systems[END_REF] like the Kubo-Greenwood [START_REF] Kubo | [END_REF][START_REF] Greenwood | Proc. Phys. Soc[END_REF] or the Landauer [START_REF] Landauer | [END_REF] formalisms, or the nonequilibrium Green's function theory [10][11][12] . In the last two decades, the combination of these formalisms with density-functional theory (DFT) or many-body perturbation (MBPT) theory allowed to establish ab initio approaches to quantum transport. The DFT-Landauer framework is one of the most popular. It has proven successful in calculating zero-bias conductances in good agreement with the experiment in some systems like the hydrogen molecule in between platinum wires 13 . In other systems, like organic molecule junctions, the DFT-Landauer estimate can be several orders of magnitude larger than the experiment 1,14 . Several solutions have been proposed to alleviate this discrepancy as by self-interaction corrections 15,16 , hybrid mixed Hartree-Fock approaches 17 , many-body model [18][19][20][21] or ab initio GW corrections 22,23 , arising in a not yet solved controversy [24][25][26][27][28][29][30][31][32] .

Besides calculating or measuring, a physical interpretation of the conductance is needed. At the end, we would like a complete picture of the mechanisms governing quantum transport in order to fully understand the behavior of the molecular junction as an electronic device. To this end, it is important to establish a relationship between the conductance and the electronic structure, for example determining the main ingredients influencing the absolute value of the zero-bias conductance A very common approach for providing such an interpretation proceeds as follows. A set of molecular orbitals (MOs) associated to the central molecule are identified and classified according to the energy levels, e.g. the highest occupied molecular orbital (HOMO), or the lowest unoccupied molecular orbital (LUMO), or the next one (LUMO+1), etc. Then, the total electronic density of states (DOS) is decomposed into the projected density of states (PDOS) associated to each different MO. Finally, by directly comparing the conductance profile G( ) with the various PDOS, one tries to establish a correspondence between conductance features and MOs. In particular, one tries to understand which MO has the largest influence on the zero-bias conductance.

The purpose of this work is to investigate how meaningful (or on the contrary misleading) this analysis is. How reliable are the interpretations that one can get? How pertinent is it to a correct understanding of the behavior of the system? We analyze two common benchmarks, the above mentioned molecular junctions of BDA and BDT in between gold leads in order to answer these questions and solve the problems evidenced in traditional methodologies. In particular, we propose a new method to identify MOs and the associated PDOS which clearly goes in this direction, though further work is still required. Though the findings of this work may seem quite theoretical at first sight, they will have an important impact in the experimental community. Indeed, theoretical analysis of quantum transport is often used for interpreting the measurements, predicting trends (for example, for the sign of the thermopower), for obtaining independent arguments, or checking the validity of the experimental work.

The paper is organized as follows: Sec. I introduces quantum transport ab initio theory, together with the definitions of all the relevant quantities and the two traditional methods to identify MOs and PDOS. In Sec. II and III, we present the results for the BDA and BDT molecular junctions, respectively. Sec. IV is devoted to the presentation of our new method and its results on BDT. Sec. V gives a critical discussion of the physical meaning of the interpretation provided by the traditional methods and our new one.

I. THEORY: MOLECULAR ORBITALS, PDOS AND CONDUCTANCE

In the DFT-Landauer framework, the molecular junction is modeled by a central region (C) connected to two semi-infinite leads (left L and right R). Its conductance G( ) as a function of the energy of the injected electrons is given by the Landauer formula:

G( ) = 2e 2 h M ( )T ( ) = 2e 2 h tr[Γ L ( )G r C ( )Γ R ( )G a C ( )].
M ( ) is the number of modes at a given energy . T ( ) is their transmittance. Γ L/R ( ) is the left/right leads injection rate. G r/a C ( ) is the retarded/advanced Green function for the central region. The quantities G r/a C ( ) and Γ L/R ( ) can be obtained from the DFT electronic structure [i.e. the energies n and wavefunctions φ n (r)] of the central region containing an "extended molecule" and of the leads (treated as infinite periodic solids), respectively. The central extended molecule actually consists of the molecule itself plus some layers actually belonging to the leads. The number of included layers (typically 3 or 4) should account for the relaxation of both the atomic and the electronic structures of the junction. The value assumed by G( ) at the Fermi energy F (which will be set to 0 in the following), G( = F =0), is an observable that can be directly measured in experiments and referred to as the zero-bias conductance.

The junction conductance depends on the nature and the shape of the leads, the geometric/atomic structure of the molecule-lead contact, and the molecule itself. Experiments and calculations very often only consider gold for the leads, so that these can be considered a constant ingredient. In contrast, the geometry of the moleculelead contact may vary quite a lot, but in many cases it is not known and furthermore it cannot be controlled experimentally. In practice, experiments only measure conductances averaged over the different possible geometries. In the end, the main factor influencing the junction conductance is the central molecule. Therefore, there is quite a lot of interest on how the conductance changes by varying the chemical composition or the atomic structure of the central molecule. Furthermore, by looking at the generic representation of the molecular junction, the central molecule appears as a "bottleneck" to the stream of electrons flowing from one lead to the other. For this reason, it is believed that the molecule itself and its electronic structure has a deep influence on the conductance.

The interpretation of the conductance profile G( ) or of the zero-bias conductance G(0) is often carried out by referring to the projected density-of-states onto molecular orbitals (see next Section). Traditionally, these are identified using two methods that will be detailed later.

A. Interpretation of the conductance by the PDOS Supposing that a set {m} of molecular orbitals with wavefunctions φ MO m (r) have been identified, the projected density of states, ρ m ( ), on the molecular orbital m is defined as

ρ m ( ) = n φ MO m |φ n δ( -n ), (1)
where n runs over all the states of the central extended molecule with wavefunction φ n (r) and energy n . In Eq. ( 1), the Dirac delta function is usually replaced by a Gaussian function with a given broadening.

As discussed above, G( ) is fundamentally determined by the electronic structure of the central extended molecule. In particular, the DOS ρ( )= n δ(n ) should play a major role. For instance, the conductance will be zero where the number of modes M ( )=0, and so will be the density-of-states. Hence, it is quite natural to interpret the conductance with the help of the DOS. More specifically, it has become very common to analyze G( ) in terms of the different partial molecular components which enter the full DOS, i.e. the PDOS on the various MOs [33][34][35][36][37] . Since the energy region of interest for the conductance is that one around the Fermi level, one usually takes into account the molecular orbitals around the fundamental gap, e.g. the highest occupied molecular orbital (HOMO), the lowest unoccupied molecular orbital (LUMO), and the next ones, the LUMO+1, LUMO+2, HOMO-1, etc.

The analysis of the conductance in terms of the PDOS is based on a one-to-one comparison of G( ) with ρ m ( ) for some chosen MOs. Whenever a peak in G( ) is in correspondence with a peak in a ρ m ( ), that molecular orbital m is said to "drive" the peak of conductance. The specific case of the 0-bias conductance is a bit particular. Indeed, very often, G(0) is quite small and the main conductance peaks are several eV away. The 0-bias conductance is actually interpreted as the tail of one of these peaks. But there is some ambiguity about the MO which will be said to drive G(0). Indeed, it can be chosen as: i) the MO corresponding to the peak closest to the Fermi level ( =0) 38,39 ; ii) the MO presenting the highest PDOS value at =0, no matter how far the PDOS maximum is from =0 23,40,41 .

B. Identification of the molecular orbitals

The molecular orbitals φ MO m (r) are the fundamental ingredient of the PDOS [see Eq. ( 1)]. As shown below, the approach chosen for identifying the MOs strongly affects the PDOS and the consequent interpretation of the conductance spectrum. Two main methods have been used so far in the literature for identifying MOs: Method 1: The φ MO m (r) are chosen to be the eigenfunctions of Hamiltonian of the uncontacted, gas phase, isolated molecule 41 . For consistency, they are usually determined using exactly the same supercell of the extended molecule, as in the molecular junction calculation, and removing the atoms of the leads.

Method 2:

The Hamiltonian of the extended molecule is first expressed on a real-space localized basis set. This can be achieved, for instance, using maximally localized Wannier functions (MLWFs) 42 . The φ MO m (r) are then chosen as the eigenfunctions of the submatrix constructed by considering only basis elements localized on the molecule 13 .

There is no obvious reason why the MOs identified using these two different procedures should coincide. Furthermore, it is not evident which method is preferred with respect to the assumed purpose, i.e. the analysis of the conductance. Method 1 coincides with the rigorous definition of MOs in the chemistry sense for the isolated molecule. However, the electronic structure of the extended molecule (taking into account charge transfer and other modifications induced by the contact between the molecule and the leads) is clearly much more important with respect to the conductance profile. So that Method 2 appears more relevant for the analysis of the conductance.

Note that choosing of one of these methods does not affect the conductance profile, provided that convergence is reached. What actually changes is rather the PDOS and hence the interpretation of the conductance in these terms.

C. Computational details

Our calculations are carried out within the DFT-Landauer framework. The exchange-correlation energy In contrast, the HOMO and LUMO in the junction present a clear bonding with the leads and thus slightly differ from the corresponding MOs in gas-phase. Hydrogen, carbon, nitrogen and sulfur atoms are represented by white, grey, brown and green spheres, respectively.

is approximated using the PBE functional 43 . We use ABINIT 44 for ground state calculations and WanT [START_REF] Ferretti | WanT code[END_REF]46 to construct Wannier functions and for conductance calculations. All the results presented here are obtained by well converged calculations, using the same convergence parameters as in Ref. 23, which are consistent and in agreement with the literature.

II. RESULTS FOR BENZENE-DIAMINE

A. BDA molecular orbitals

In Fig. 1, we show the molecular orbitals of BDA calculated with Methods 1 and 2. They are analogous to those found previously for instance using Method 1 41 . While the HOMO-1 molecular orbitals are very similar, the HOMO show non-negligible differences: the bonding character with the leads is more important when using Method 2, as indicated by the more pronounced lobes on the N atoms that point towards the gold adatoms.

We observe a close similarity between the LUMO from Method 1 and the LUMO+1 from Method 2, like if there were a change in the ordering of the corresponding eigenvalues between the two methods. Notice that the energy difference between the LUMO and the LUMO+1 is ∼0.5 eV, so enough to exclude their degeneracy. Vice versa, the LUMO+1 from Method 1 resembles the LUMO from Method 2 but there are some small differences: the bonding character with the leads is again more pronounced when using Method 2. In fact, the corresponding density arises from a MLWF basis element which is localized on the gold-amino bond and not clearly identifiable as purely belonging to gold or to the molecule. In this MO important differences are also found for the lobes on the benzene-ring: in Method 1, the lobes are mainly on the opposite C atoms along the molecule long axis; whereas, in Method 2, they are on the C atoms close to the Au adatom.

These differences will induce non-negligible differences in the PDOS analysis, as will see in the next section.

B. PDOS and interpretation of the conductance

In Fig. 2(c), we show the conductance of BDA calculated in the Landauer-DFT framework using the PBE approximation. And, as it is usually done in literature for providing a physical interpretation of the conductance, we also present the PDOS as calculated using Methods 1 [Fig. 2(a)] and 2 [Fig. 2(b)]. The position and height of the main features are in very good agreement with previous work 41 .

The two PDOS look quite similar but with differences that can be associated to the already discussed discrepancies between MOs. In particular, we observe the change in the ordering between the LUMO and the LUMO+1 from Method 1 to 2. The PDOS onto non-hybridized MOs (HOMO-1 and LUMO/LUMO+1 in Method 1/2) look similar, whereas the PDOS onto the HOMO and LUMO+1/LUMO in Method 1/2 present differences, as expected from the MO plots. Finally, the PDOS onto HOMO-2 seems to have more weight in Method 1 than in Method 2.

When interpreting the conductance profile, one can associate the small conductance peak at ∼1.5 eV with the intense LUMO and LUMO+1 PDOS peaks observed respectively in Methods 1 and 2. The conductance struc- The insets show a zoom on the PDOS zoom around the Fermi energy region. Notice that, in the inset of (b), the PDOS is presented in logarithmic scale.

ture arising at energies > 0 eV with maximum at 2 eV could be correlated to the other unoccupied molecular orbital (LUMO+1 of BDA-gas alias LUMO of BDAjunction), as well as the LUMO+2. The peak in the conductance at ∼-1.5 eV could be related to the HOMO PDOS peak at ∼-1.4 eV, and so also the structure from -2 down to -3.8 eV. The HOMO-1 and its PDOS peak at ∼-3.6 eV does not reflect in the conductance. However, when performing a one-to-one comparison of the conductance with the total PDOS on the various MOs (Fig. 3), the relationship does not look that direct, even qualitatively.

We finally discuss the interpretation of the 0-bias conductance. Following one possible interpretation scheme very common in the literature, the zero-bias conductance appears on the tail of the conductance peak at -1.5 eV (HOMO), though the smallest peak at +1.5 eV (asso- ciated to the PDOS onto the LUMO/LUMO+1 from Method 1/2) is equally distant. According to this interpretation, the zero-bias conductance is driven by the HOMO, though a contribution from the LUMO from Method 1 (alias the LUMO+1 from Method 2) is expected.

-4 -2 0 2 ε -ε F [eV] 0 0.2 0.4 0.6 Conductance [G 0 ] BDA -4 -2 0 2 ε -ε F [eV]
These conclusions are contrasted by another approach which rather looks at the absolute values of the PDOS at the Fermi energy (see Fig. 2 insets showing zooms on the Fermi energy regions). According to this scheme, the other unoccupied MO (the LUMO+1 from Method 1, alias the LUMO from Method 2) drives the zero-bias conductance. In fact, both methods agree on the fact that this MO (labeled differently) presents the largest PDOS value at the Fermi energy. Nevertheless, its corresponding PDOS value at 0 eV disagrees by one order of magnitude: from 0.1 in Method 2 to 0.02 in Method 1. The next MO presenting an important PDOS value at the Fermi energy is the HOMO-2 from Method 1, with a value even not much smaller than the LUMO+1, implicating that the HOMO-2 has a certain weight on the zero-bias conductance. However, this is the HOMO from Method 2 with a marked gap (from 0.1 to 0.01). Both methods agree about the HOMO PDOS absolute value (∼0.01) at = 0, probably by mere coincidence given the disagreements mentioned above.

Summarizing, when interpreting the BDA zero-bias conductance we are in front of 3 problems: i) arbitrariness in the labeling of MOs (the LUMO in Method 1 becomes the LUMO+1 in Method 2, and vice versa); ii) dependence on the method to identify the MOs; iii) dependence on the interpreting approach. Hence PDOS analyses of G( = 0) are affected by some ambiguity.

III. RESULTS FOR BENZENE-DITHIOL

A. BDT molecular orbitals

We now consider a more complex case: the benzene-dithiol (BDT) gold junction.

Experimen-tal and theoretical works concluded that the BDTgold junction can be stable in several different atomic structures/geometries 16,[47][48][49][50][51] . To account for different hybridizations and bonding motifs, three geometries are studied here: the sulfur atom of the benzene-dithiol molecule can adsorb to an extra gold adatom without loosing the bound hydrogen atom (BDT-n); the benzenedithiol molecule can loose the hydrogen, thus becoming benzene-dithiolate, and bind its sulfur atom to an extra gold adatom in a pyramid structure (BDT-p); or the benzene-dithiolate can bind to 3 equidistant gold atoms in the hollow structure (BDT-h). These geometries are shown in Fig. 4. In Fig. 4, we show the MOs of BDT calculated with different methods. For Method 2, we show the molecular orbitals obtained for the 3 different junction geometries: BDT-n, BDT-p and BDT-h. They are very similar to those obtained previously 13 , especially given the differences in the considered geometries. In Ref. 52, an alternative set of MOs are shown for BDT-h, obtained within Method 2 by considering only the localized orbitals on the benzene molecule (excluding the S atoms). For Method 1, we depict both the cases of benzene-dithiol and benzenedithiolate. The latter might better represent the BDT molecule in the BDT-p and BDT-h junctions where it looses a hydrogen atom before binding. But this is not so straightforward: besides the effective chemical composition of the molecule in the junction, other chemical/physical effects, e.g. saturation of bonds, transfer of charge, may be considered [24][25][26][27] .

We start by analyzing the MOs from Method 1. The MOs for the dithiol and dithiolate molecules present a few similarities. The LUMO+1 are similar in shape. The HOMO of the dithiol molecule resembles to the LUMO of the dithiolate molecule, with an exchange of the ordering as already seen in BDA (see previous Section). Nevertheless, other MOs strongly differ. So, identification of MOs using Method 1 strongly depends on the molecule (dithiol vs. dithiolate).

Now we analyze the MOs obtained with Method 2. We focus on BDT-p, the junction in which the interpretation of conductance using the PDOS is the most critical of all the cases considered here, as will be seen later. The LUMO+1 from Method 2 looks very similar to the LUMO+1 from Method 1 for both the dithiol and dithiolate molecules, though with differences on the sulfur atom. The LUMO from Method 2 corresponds to the LUMO from Method 1 for the dithiol molecule, but it has no correspondence to any MO from Method 1 for the dithiolate molecule. On the other hand, the HOMO from Method 2 is similar to the HOMO from Method 1 for the dithiolate molecule, but it differs from all MOs from Method 1 for the dithiol molecule. Finally, the MOs which look closer to the HOMO-1 from Method 2 are the HOMO from Method 1 for the dithiol molecule and the LUMO from Method 1 for the dithiolate molecule.

From the above discussion, it appears that no one-toone correspondence can be established between the MOs obtained with the two methods nor between the MOs from Method 1 both for the dithiol and dithiolate isolated molecules. The BDT-p MOs from Method 2 seem at halfway between the MOs from Method 1 for the dithiolate and dithiol molecules.

ε -ε F [eV] ε -ε F [eV] ε -ε F [eV] 0.2 -0.2 0 0.2 -0. (c) (d)

B. PDOS and interpretation of the conductance

We now move to the analysis of the most critical case (among the examples investigated here) regarding the interpretation of the conductance in terms of the PDOS: benzene-dithiol in the pyramid geometry (BDT-p). In Fig. 5 we present the Landauer-DFT conductance of BDT-p. On top, we present also 3 different PDOS calculated following Method 1 (gas phase) and Method 2 (junction), for the former both the dithiol and dithiolate molecules are considered.

Without entering into all details, it is clear that the PDOS strongly depends on the method used to calcu-FIG. 6. Scheme representing the integration volume (shadowed area passing through the two S atoms of the BDT-h junction) used for our Method 3. late it, reflecting previously seen differences in the MOs. For instance, the zero-bias conductance seems dominated by the HOMO from Method 1 for the dithiol molecule, since the PDOS onto the HOMO is the closest to the Fermi level and it also provides the highest contribution at that level (see the inset), with minor contribution from the HOMO-1. When using Method 1 for the dithiolate molecule, the zero-bias conductance seems equally driven by the HOMO and LUMO, with also some contribution from the HOMO-1 and the HOMO-2. Finally, using Method 2, the HOMO-1, the HOMO and the HOMO-2 (in decreasing order) are the most important contributions at zero-bias . Though some discrepancies can be ascribed to simple re-labeling of the same MO, one cannot pass over more important differences among the methods.

In conclusion, we could not find a rigorous definition of the MOs and associated PDOS for the BDT-p case when using the traditional methods. As a consequence, the PDOS interpretation of the conductance does not rely on stable grounds.

IV. NEW METHOD FOR IDENTIFYING MOLECULAR ORBITALS

A. Charged isolated molecules

In order to reconcile the two main methods found in literature, that is, smooth their differences and solve the difficulties, we here propose a new approach which is based on an evolution of Method 1.

Method 3:

The φ MO m (r) are chosen as the eigenfunctions of the Hamiltonian of the uncontacted, gas phase, isolated molecule, to which some charge is added accounting for metal-molecule chargetransfer. The same supercell is used as in the contacted molecule junction calculation, but removing the atoms of the leads. The added charge is calculated from a three steps procedure: i) the density ρ(r) of the complete junction is computed;

ii) the density ρ (r) of the molecule is also calculated using the same geometry and simulation box as in the junction; iii) the added charge is given by integrating ρ(r)-ρ (r) over the volume spanned by the molecule. For BDT-h, this volume is given by the region between two planes perpendicular to the S-S axis and passing through the two S atoms (see Fig. 6).

The rationale behind our new method is to modify the electronic structure of the gas-phase isolated molecule with the purpose to account for the lead-molecule charge transfer. Thus, the isolated molecule is placed into an environment closer to that one of the molecular-junction. Previous studies 24,26 have already underlined the importance of the lead-molecule charge transfer and the significance of its role in transport properties of molecular junctions. Here, it constitutes the basis for the construction of a new method of analysis.

B. Application of the new method to BDT-h

We apply our new method to the case of BDT-h (hollow geometry), which presents contradictory results using standard methods, as explained later. According to our recipe, the extra charge to be added to BDT-thiolate to simulate the environment of the BDT-h junction was found to be ∼0.5 e -. However we observe that the modifications of the MOs are slightly affected by the precise value of the added charge, apart when the charge crosses integer values, ρ=0, 2, . . ., of the electronic unit charge e -, at the onset of the occupation of new levels. The MOs found with this procedure are shown at the bottom of Fig. 4. Remarkably, these MOs look now much more similar to the MOs found with Method 2 for BDT-h, as it can be clearly seen. Furthermore, they present marked differences with the original Method 1 for dithiolate, and in some cases are even closer to Method 1 for the dithiol molecule.

Fig. 7 shows the PDOS for the BDT-h junction calculated with the traditional methods (Method 1 for dithiol and dithiolate isolated molecules and Method 2 for a selected set of MLWFs for the BDT-h junction) and our new Method 3. We focus on the PDOS around =3 eV where traditional methods present the most important differences. In that energy region, our new Method 3 provides an evident improvement. Coming from a dithiolate isolated molecule, the PDOS from Method 3 is closer to the one from Method 1 for the dithiol molecule than for the dithiolate molecule, thus bridging the gap between the dithiol and dithiolate molecules. Moreover, in this same energy region, when considering the relative height between the PDOS peaks of LUMO and LUMO+1, Method 3 evidently bridges the gap between Method 1 for the dithiol and dithiolate molecules and Method 2. We can probably conclude the same also for the =0 eV region, though restricting the discussion to the PDOS of the HOMO. One can observe the evolution of the PDOS peak of the HOMO at the Fermi energy from Method 1 for the dithiol molecule, from Method 3 and from Method 2. We can say that Method 3 is somehow successful in reconciling the traditional Methods 1 and 2. However, we do not notice any other evident improvement. We still find clear differences among the PDOS when zooming on the =0 eV region (not shown). The MO ordering problem continues to exist: the PDOS peak at ∼-2.5 eV from Method 3 is attributed to yet another MO, the HOMO-4. The same ambiguous attribution remains for the PDOS of the intermediate HOMO orbitals.

ε -ε F [eV] 0 0.5 Conductance [G 0 ] 0 (e) (a) (b) (c) (d)
ε -ε F [eV] ε -ε F [eV] ε -ε F [eV] ε -ε F [eV] -0.2 0.2 0 -0.2 0.2 0 -0.2 0.2 0 -0.
We have tried Method 3 also on the more complex case of BDT-p. The MOs from Method 3 (not shown) do not resemble to those from Method 2, and consequently we get no satisfactory results on the PDOS. BDT-p continues to be an unsatisfactory case also for Method 3. This is so probably because the metal-molecule charge transfer is not the only, or the main, parameter affecting the electronic structure of BDT-p, due to a may be higher metal-molecule coupling and hybridization.

In conclusion, Method 3 provides encouraging partial satisfactory results, in particular in reconciling the two traditional methods as in BDT-h. However, this is not general and not all problems are solved. The metalmolecule charge transfer is not the only mechanism at play. One should probably take into account also the metal-molecule hybridization. This is not an easy task if the purpose is to keep the picture of an isolated molecule.

V. DISCUSSION: FURTHER CONSIDERATIONS ON THE PDOS ANALYSIS

As discussed in the previous section, Method 3 aims at overcoming the drawbacks related to the identification of MOs using Method 1. Instead, one could have explored the possibility to improve upon Method 2. However, as we argue hereafter, this path appears to us less physically-grounded. It actually opens even more fundamental questions on the implicit hypotheses at the basis of the interpretation of the conductance based on the PDOS, and raises further doubts on the validity of the whole procedure.

A. Dependence of MOs and PDOS from the choice of Wannier Functions basis set

At first sight, Method 2 (for which MOs originate from the junction) would seem more meaningful for studying the conductance. However, it presents a severe drawback for which it seems very difficult to find a solution. There is a certain arbitrariness in the criterion establishing the spatial limits of a molecule and thus the basis elements that will be considered as being "localized on the molecule". For instance, there can be MLWFs localized on the molecule-lead bonds as we have pointed out for BDA. It is then quite arbitrary to say whether they are localized on the molecule or on the leads. This choice clearly affects the resulting submatrix, as well as the number and the shape of the MOs found after its diagonalization.

Intuitively, these basis elements should have an important effect on the junction conductance, so that it makes a lot of sense to keep them when generating the MOs. Coming back to the case of BDA, the most important PDOS at the Fermi energy was precisely the one associated to the MO presenting the major localization on the molecule-bond MLWF (i.e. the LUMO). If we had discarded the latter from those "localized on the molecule", we would have excluded this important MO from the analysis of the zero-bias conductance. It is actually reassuring that this MO also appeared when using Method 1, though labeled LUMO+1 due to the already discussed inverted ordering (see Fig. 1) and it was also the most important PDOS at F . But, at the same time, it shows that the exclusion of some MLWFs based on their localization may lead to very different interpretations starting from Method 1 or Method 2.

A strategy to circumvent this drawback is to select a different set of Wannier functions (WFs), or any other localized basis set with elements presenting a well-defined localization (on the molecule or on the leads). For instance, atom-centered basis sets would resolve this ambiguity, such as symmetry-adapted WFs 53 , WFs obtained from LCAO projections 54 , or LCAO basis sets. Furthermore, it is well-known that, in some cases, the Marzari-Vanderbilt 42 algorithm can lead to different sets of WFs. For instance, silicon bulk presents at least 3 different sets of WFs with a similar degree of localization (as measured by the spread S). When starting the Marzari-Vanderbilt algorithm from a random initial guess, there is a high probability to fall down into the global minimum (S=2.56 Å2 ) for which the lowest 8 MLWFs are of the sp 3 -backward kind [Fig. 8(c)] which do not correspond to the real chemical orbitals. It is obviously possible to obtain the 8 sp 3 -forward WFs [Fig. 8(b)] which correspond to the physical chemical sp 3 orbitals, but at slightly higher local minimum (S=2.95 Å2 ). Finally, the set of WFs with 4 bonding orbitals on one Si atom and 4 anti-bonding orbitals on the other atom [Fig. 8(a)] has a relatively large spread (S=5.09 Å2 ). However, when performing the search of the MLWFs for the 4 valence states only, the minimum spread is obtained for a set containing the 4 bonding orbitals.

The previous discussion points to a possible ambiguity in Method 2 for identifying the MOs and hence in using the corresponding PDOS to interpret the conductance. For a single junction, one may find several sets of WFs. The one presenting the minimum spread (the most localized) does not necessarily correspond to the real physical situation, and this cannot be known a priori. The calculated conductance must and does not depend on the chosen basis set, provided the basis is complete and at convergence. On the other hand, the submatrix of the junction Hamiltonian does depend on the chosen basis set. So do its eigenfunctions (which define the MOs) and the resultant PDOS. Consequently, the physical interpretation of the conductance by the PDOS does depend on the chosen WF or other basis set. A basis-dependent interpretation method is questionable.

Starting from this point, we are led to ask even more fundamental questions: Is the conductance really related to a MO, or a PDOS, or to some MOs and a total PDOS? Before answering these questions, let us try to answer a question even further upstream.

B. Is the conductance directly related to the full DOS?

The conductance G( ) is certainly directly related to the electronic structure of the junction, i.e. to both the electronic energies n and wavefunctions φ n (r) of the extended molecule. Hence, there should also be a relationship to the total DOS, ρ( ), though somehow indirect and not one-to-one. For instance, wherever ρ( )=0 (there are no states available at that energy), the conductance G( ) must be also zero. The reverse is not true: the conductance can be zero at energies where the total DOS is finite. This can happen at energies associated to strongly localized wavefunctions, with zero spatial overlap among them, for example core states. There can also be other factors beyond localization altering the direct relationship between G( ) and ρ( ). For instance, not all delocalized wavefunctions are good conducting channels 55 . As a result, direct conclusions cannot be drawn from the inspection of the DOS only.

C. Is the conductance related to some kind of PDOS?

Whether the conductance is directly related to some kind of PDOS, be it onto a given MO or onto some MOs or even the total PDOS, is actually less obvious to answer than for the full DOS. And so is the physical interpretation of the conductance based on such quantities.

Taking the example of BDT-h (Fig. 7), one can see that the conductance profile is qualitatively related to a total PDOS including the MOs which are close to the Fermi energy. Nevertheless, it is not possible to observe a quantitative relationship between the conductance value and the total PDOS height. This is more evident in the case of BDT-p (Fig. 5), one cannot explain why the conductance is larger at -1 eV than at 3.5. At below -1 eV, the agreement worsens even qualitatively. In the case of BDA (see Fig. 2), the relationship between the conductance and the total PDOS is even less evident.

This work has made it clear that the conductance analysis depends on a suitable choice of the MOs. For this reason, the interpretation of the conductance in terms of the PDOS is quite questionable. We should first provide an answer to the following fundamental questions: i) which set of MOs do physically represent the molecule in the junction? ii) given the lead-molecule hybridization, are MOs obtained from an isolated molecule (i.e. from Methods 1 or 3) meaningful for analyzing a metalmolecule junction? iii) are MOs obtained by diagonalizing a submatrix of the Hamiltonian (Method 2) physical, given the fact that they depend on the choice of basis set?

MOs identified as the eigenvectors of the gas phase, isolated Hamiltonian (Methods 1 and 3) have a physical meaning. But, this is only true for the isolated molecule not necessarily for the junction. For the latter, the eigenfunctions of the isolated molecule are nothing but yet another basis set (just like the atomic orbitals for a solid). Furthermore, the actual choice of the molecule may not be unique (e.g., dithiol or dithiolate). As for Method 2, an interpretation which depends on the chosen basis set (WFs, LCAO, Gaussians or Wavelets) cannot be considered physical.

We believe that a completely different direction should be taken in order to provide an answer to these questions. What matters for a physical interpretation of the conductance is the full electronic structure of the extended molecule (containing also some layers of the leads). Considering the extended molecule system needed to converge the conductance, which typically contains of the order of 10 2 gold and 10 1 molecule atoms, one can realize that the molecule does not even have such an important weight on the determination of the electronic structure of the junction. Following these arguments, we can give the indication that a meaningful procedure to provide a physical interpretation of a junction conductance should rely on the wavefunctions and energies directly identified for the extended molecule electronic structure. Thus, in order to provide a physical interpretation of the conductance, we believe that the local density-of-states (LDOS), a quantity independent from the basis set and directly built on the extended molecule wavefunctions and energies, is the most meaningful. Actually, we have already presented an application which uses the LDOS for the interpretation of the quantum transport conductance 23 .

Regarding an interpretation of the molecular junction conductance rooted on the molecular PDOS, this work first tried to reconcile the two traditional methods (Methods 1 and 2) by introducing a new one (Method 3). Some success was achieved in this direction, but we cannot consider the problem to be solved. Further work is clearly needed. However, our considerations led us to doubt that a fully satisfactory solution exists along this direction.

CONCLUSIONS

Taking as examples two reference molecular junctions (benzene-diamine and benzene-dithiol between gold leads), we have investigated the interpretation of the conductance based on the projected density of states (PDOS) onto molecular orbitals. The latter are usually identified following two procedures: diagonalization of Hamiltonian of the gas-phase isolated molecule (Method 1); and diagonalization of a submatrix of the junction Hamiltonian constructed by considering only basis elements localized on the molecule (Method 2). We have shown that these two methods can lead to substantially different MOs and hence PDOS. Furthermore, within Method 1, the PDOS depends on the isolated molecule chosen to represents the junction (e.g. with or without dangling bonds); and, within Method 2 the PDOS depends on the chosen basis set. As a consequence, the analysis of the conductance based on the PDOS can lead to different, if not contrasting, conclusions. This is particularly true for the analysis of the zero-bias conductance which can be found to be driven by, e.g., the LUMO in one method and the HOMO in another. To go beyond these drawbacks, we proposed an alternative method (Method 3) as an improvement to Method 1. This new method somehow reconciles Methods 1 and 2, but still presents problems which point to more fundamental questions. An analysis of the conductance based on the PDOS seems not to rely on well established roots due to the arbitrariness in the identification of MOs. Our proposal provided some indications toward possible solutions to the problem of interpreting the molecular junction conductance.

FIG. 1 .

 1 FIG.1. Electronic density isosurfaces (red) of the HOMO-1, HOMO, LUMO and LUMO+1 molecular orbitals of BDA as obtained with the two traditional methods (see text). The ordering of the LUMO and LUMO+1 is inverted in the two methods. The localized MOs [HOMO-1 and LUMO (gas) or LUMO+1 (junction)] look remarkably similar for both methods. In contrast, the HOMO and LUMO in the junction present a clear bonding with the leads and thus slightly differ from the corresponding MOs in gas-phase. Hydrogen, carbon, nitrogen and sulfur atoms are represented by white, grey, brown and green spheres, respectively.

FIG. 2 .

 2 FIG. 2. Projected density of states (a,b) and conductance (c) of benzene-diamine (BDA). The PDOS for the different molecular orbitals (from HOMO-2 to LUMO+1) have been obtained with (a) Method 1 and (b) Method 2 (see text).The insets show a zoom on the PDOS zoom around the Fermi energy region. Notice that, in the inset of (b), the PDOS is presented in logarithmic scale.

2 FIG. 3 .

 23 FIG.3. Total molecular PDOS (red line) and conductance (black line) of BDA. The total molecular PDOS is the sum of the PDOS onto the MOs from HOMO-2 to LUMO+2 as obtained from Method 2. Note that the PDOS onto LUMO+2 is not shown in Fig.2.

FIG. 4 .

 4 FIG.4. Electronic density isosurfaces (red) of the HOMO-1, HOMO, LUMO and LUMO+1 molecular orbitals of BDT as obtained with the two traditional methods as well as with the new method (see text). For Method 1, the dithiol and dithiolate molecules are considered. For Method 2, the different molecular junction geometries (BDT-n and BDT-p and BDT-h) are examined. For Method 3, a charge of +0.5 e -was added to the dithiolate molecule in order to account for the transfer of charge to the molecule from gold atoms in the BDT-h junction. The resulting orbitals are very similar to those obtained with Method 2 for BDT-h. Hydrogen, carbon, nitrogen and sulfur atoms are represented by white, grey, brown and green spheres, respectively.

FIG. 5 .

 5 FIG. 5. Projected density of states (a,b,c) and conductance (d) of benzene-dithiol in the pyramid geometry (BDT-p). The PDOS for the different molecular orbitals (from HOMO-3 to LUMO+1) have been obtained with (a) Method 1 based on the dithiol molecule, (b) Method 1 based on the dithiolate molecule, and (c) Method 2 (see text). The insets show a zoom on the PDOS zoom around the Fermi energy region.

3 FIG. 7 .

 37 FIG. 7. Projected density of states (a,b,c,d) and conductance (e) of benzene-dithiol in the hollow geometry (BDT-h). The PDOS for the different molecular orbitals (from HOMO-3 to LUMO+1) have been obtained with (a) Method 1 based on the dithiol molecule, (b) Method 1 based on the dithiolate molecule, (c) Method 2, and (d) Method 3 (see text). The insets show a zoom on the Fermi energy region.

  FIG. 8.Using the Marzari-Vanderbilt algorithm, three different sets of Wannier functions (WFs) with comparable spread can be obtained for bulk silicon. While a) bond-ing+antibonding and b) sp 3 -forward are the most "physical" WFs, though not the most localized ones, c) sp 3 -backwards are the maximally localized WFs.

ACKNOWLEDGMENTS

We thank Pierre Darancet, Jeff Neaton and Xavier Blase for useful discussions. GMR acknowledges the F.R.S.-FNRS for financial support.

Computational resources have been provided by the supercomputing facilities of the Université catholique de Louvain (CISM/UCL), by the Consortium des Équipements de Calcul Intensif en Fédération Wallonie Bruxelles (C ÉCI), and by the French GENCI supercomputing center (Project i2012096-655).


Z Yong

P E Trevisanutto

L Chiodo

I Santoso

A R Barman

F Della Terentjevs

Sala

T C Asmara

S Dhar

A Kotlov

A Terentjevs

F Della Sala

V Olevano

M Rübhausen

T Venkatesan

A Rusydi

Emerging giant

come

Titanium dioxide (TiO2) has rich physical properties with potential implications in both fundamental physics and new applications. Up-to-date, the main focus of applied research is to tune its optical properties, which is usually done via doping and/or nano-engineering. However, understanding the role of d-electrons in materials and possible functionalization of d-electron properties are still major challenges. Herewith, within a combination of an innovative experimental technique, high energy optical conductivity, and of the state-of-the-art ab initio electronic structure calculations, we report an emerging, novel resonant exciton in the deep ultraviolet region of the optical response. The resonant exciton evolves upon low concentration Ta-substitution in anatase TiO2 films. It is surprisingly robust and related to strong electron-electron and electron-hole interactions. The dand f -orbitals localization, due to Ta-substitution, plays an unexpected role, activating strong electronic correlations and dominating the optical response under photoexcitation. Our results shed light on a new optical phenomenon in anatase TiO2 films and on the possibility of tuning electronic properties by Ta substitution.

I. INTRODUCTION

Doped or defective titanium dioxide (TiO 2 ) exhibits rich physical phenomena in electronic transports and optical properties [START_REF] Dulub | Electron-Induced Oxygen Desorption from the TiO2(011)-21 Surface Leads to Self-Organized Vacancies[END_REF][START_REF] Crossland | Mesoporous TiO2 single crystals delivering enhanced mobility and optoelectronic device performance[END_REF][START_REF] Richter | Exciton-like trap states limit electron mobility in TiO2 nanotubes[END_REF][START_REF] O'regan | A low-cost, high-efficiency solar cell based on dye-sensitized colloidal TiO2 films[END_REF] . TiO 2 is opaque in the visible sun light whereas it is very efficient in absorbing ultraviolet (UV) light rendering it interesting especially for photocatalysis applications [START_REF] Chung | All-solid-state dye-sensitized solar cells with high efficiency[END_REF] . The first step of photoexcitation is the formation of electron-hole pair quasiparticles (excitons), which may either recombine or decay into free charges. Eventually the free charges react with molecules on the surface enhancing photocatalytic effects and formation of reactive free radicals [START_REF] Mcfarland | A photovoltaic device structure based on internal electron emission[END_REF] . Excitons, and their spatial behaviour, play therefore a key role for both fundamental physics as well as for applications, but the precise nature and behaviour of excitons in TiO 2 based materials remains unclear in some respects.

The many-body electron-electron (e-e) and electronhole (e-h) interactions determine the physical properties of excitons, with different contributions depending on the system and on the considered energy range. Excitons usually occur below direct bandgap in semiconductor and insulator materials, but they may involve higher energy bands in the case of strong electronic correlation. With the recent development of supercomputing and ab initio calculations [START_REF] Albrecht | Ab Initio Calculation of Excitonic Effects in the Optical Spectra of Semiconductors[END_REF][START_REF] Benedict | Optical Absorption of Insulators and the Electron-Hole Interaction: An Ab Initio Calculation[END_REF][START_REF] Rohlfing | Electron-Hole Excitations in Semiconductors and Insulators[END_REF] , theoretical studies have shown that when both e-e and e-h interactions are strongly coupled, they yield to a new type of optical phenomenon, the so-called high-energy resonant excitonic effect. In fact, resonant excitons have been predicted [START_REF] Yang | Excitonic Effects on the Optical Response of Graphene and Bilayer Graphene[END_REF][START_REF] Trevisanutto | Ab initio high-energy excitonic effects in graphite and graphene[END_REF] and later observed [START_REF] Mak | Seeing Many-Body Effects in Single-and Few-Layer Graphene: Observation of Two-Dimensional Saddle-Point Excitons[END_REF][START_REF] Santoso | Observation of room-temperature highenergy resonant excitonic effects in graphene[END_REF] in two-dimensional graphene. Unlike excitons in conventional semiconductors, the resonant excitons can occur at arXiv:1605.06887v1 [cond-mat.str-el] 23 May 2016 energies even well above the corresponding optical band gap of the material, and they can be probed directly using high-energy optical conductivity [START_REF] Santoso | Observation of room-temperature highenergy resonant excitonic effects in graphene[END_REF] . A detailed understanding of the role of e-e and e-h interactions in TiO 2 based materials remains elusive, and resonant excitons have not been observed in the material, mainly because both experimental and theoretical studies at high-energy optical conductivity are challenging and limited in number.

We report in this paper on optical studies of TiO 2 doped at different concentrations of Tantalum, via optical conductivity measurements and ab initio Time Dependent Density Functional Theory (TDDFT) calculations. We have observed resonant excitonic effects in the deepultraviolet (DUV) in anatase Ta x Ti 1-x O 2 films, with only a small amount of Ta-substitution. A series of unusual phenomena arise, in particular the spectral-weight transfer from high towards low energies, and the emergence of an intense resonant exciton at ∼ 6 eV. Based on our theoretical calculations, we relate these effects to a peculiar manifestation of strong e-e and e-h interactions. The paper is organized as follows: in section II, experimental and theoretical-computational used techniques are described. In section III, optical spectra, both measured and calculated, are described. In section IV, the main conclusions are drawn.

II. MATERIALS AND METHODS

Details of samples preparation and characterization, optical conductivity measurements and theoretical calculations are described in this section.

A. Experimental techniques

The optical conductivity was obtained using a combination of spectroscopic ellipsometry (0.5 -5.6 eV) and UV-VUV reflectivity (3.7 -35 eV) measurements 14,15 . The spectroscopic ellipsometry measurements were performed in the spectral range between 0.5 and 5.6 eV by using an SE 850 ellipsometer at room temperature. Three different incidence angles of 60 • , 70 • and 80 • from the sample normal were used, and the incident light was 45 • linearly polarized from the plane of incidence. For reflectivity measurements in the high-energy range between 3.7 and 35 eV, we used the SUPERLUMI beamline at the DORIS storage ring of HASYLAB (DESY). The incoming photon was incident at an angle of 17.5 • from the sample normal with linear polarization parallel to the sample surface. The sample chamber was outfitted with a gold mesh to measure the incident photon flux after the slit of the monochromator. The measurements were performed in ultra-high vacuum environment (chamber pressure of 5 × 10 -10 mbar) at room temperature. Before these measurements, the samples were heated up to 400 K in ultra-high vacuum to ensure that there were no additional adsorbate layers on the surface of the samples. The obtained UV-VUV reflectivity data were calibrated by comparing it with the luminescence yield of sodium salicylate (NaC 7 H 5 O 3 ) and the gold mesh current. These as-measured UV-VUV reflectivity data were further normalized by using the self-normalized reflectivity extracted from spectroscopic ellipsometry 16 .

B. Experimental samples and preparations

Ta 2 O 5 and TiO 2 powders with high-purity (99.999%) were ground for several hours before sintering in a furnace at 1000 • C in air for 20 h. Subsequently, target pellets were made and sintered at 1100 • C in air for 24 h. Anatase Ta x Ti 1-x O 2 epitaxial thin films (with x=0, 0.018 and 0.038) of thickness 280 nm were deposited on high-quality (001) LaAlO 3 substrates by Pulsed Laser Deposition, using a 248 nm Lambda Physik excimer laser with an energy density of 1.8 J cm -2 and a repetition rate of 210 Hz. Depositions were performed for 0.51 h in a stable oxygen partial pressure of 1 × 10 -5 Torr while the substrate temperature was maintained at 750 • C. The chemical and structural properties of the samples were studied by X-ray Photoelectron Spectroscopy, electrical transport measurements, Rutherford backscattering spectrometry (RBS)/channelling, X-ray diffraction (XRD) and time-of-flight secondary-ion mass spectrometry (TOF-SIMS) as reported elsewhere 17,18 . Ion channelling measurements indicated near perfect substitutional Ta atoms in Ti sites.

C. Theoretical and computational methods

All ground state electronic calculations are carried out by using Density Functional Theory (DFT) based on Quantum Espresso 28 and Abinit 29 codes, with the Perdew-Burke-Ernzerhof (PBE)-GGA approximation for the exchange-correlation functional 30 . Norm conserving pseudopotentials in Troullier-Martins scheme [START_REF] Troullier | Efficient Pseudopotentials for Plane-Wave Calculations .2. Operators for Fast Iterative Diagonalization[END_REF] are used, and semicore electrons are included in Ti and Ta pseudopotentials. The cutoff energy for the expansion of plane-wave basis is up to 170 Ry 22 . For pristine anatase TiO 2 , we used a 12×12×8 Monkhorst-Pack kpoint mesh sampling the Brillouin zone. For anatase Ta-TiO 2 bulk we used supercells with 48 atoms and a 4×4×4 Monkhorst-Pack k-point mesh grid. One Ti atom is replaced with one Ta atom (substitutional doping, modelling a 6.5 % Ta-doping in TiO 2 bulk almost equivalent to the experimental doped sample (3.8 %)).

The excited state calculations have been performed within two approaches: solving the Bethe Salpeter Equation (BSE, which implicitly includes both e-h and e-e interactions), and applying the jellium with gap model (JGM) kernel [START_REF] Trevisanutto | Optical spectra of solids obtained by time-dependent density functional theory with the jelliumwith-gap-model exchange-correlation kernel[END_REF] within the Time Dependent Density Functional Theory (TDDFT). The latter method includes the e-e and e-h interactions maintaining the com-putational feasibility for such compelling calculations. The complex dielectric function has also been evaluated at the random phase approximation (RPA) level, with electrons and holes treated as independent particles, without correlation. The Bethe-Salpeter Equation has been solved using Yambo code 32 and EXC 33 . The screening dieletric matrix has been evaluated by using the static inverse dielectric function, with cutoffs of 21 Ry for the correlation (exchange) part, and unoccupied states are summed over 176 empty states. In BSE calculations defined hereafter as Low Resolution (LR), 28 occupied bands and 52 empty bands are included in the diagonalization, to describe the region above 5 eV, on a k-point grid of 4×4×2. For the High Resolution (HR) BSE calculation, used to describe in more detail the adsorption threshold, 8 occupied and 8 empty bands are included in the diagonalization, on a 12×12×12 k-point grid. Haydock recursive approach for diagonalization is used, with threshold accuracy of -0.02. The DP-EXC code 34 is used for the TDDFT calculations. In pristine TiO 2 bulk, 200 bands are included for the RPA and JGM-TDDFT calculations.

III. RESULTS AND DISCUSSION

Here, we provide some general information on pure TiO 2 , useful in the following discussion on optical conductivity. The electronic ground state structure of pristine and doped TiO 2 , based on DFT results for total Density of States (DOS) and Partial DOS (PDOS), is shown in Fig. 1. The DFT-PBE band gap of pristine TiO 2 is ∼2.20 eV. The valence band (Fig. 1-(a)) mainly consists of O 2p orbitals slightly hybridized with Ti 3d orbitals. The conduction band is comprised by Ti 3d orbitals with a small hybridized amount O 2p orbitals.

In Ta-TiO 2 (Fig. 1-(b)), the 3d Ta orbitals fill the bottom of conduction band, making the system metallic, and they are hybridized with the adjacent O 2p orbitals up to 8 eV in the conduction band. The proper inclusion of correlation removes this spurious metallicity described by DFT, as shown for Nb-doped rutile [START_REF] Zhang | All-electron GW calculation of rutile TiO2 with and without Nb impurities[END_REF] . Main features of electronic band structure (Fig. 2-(a) ) and optical absorption spectrum of anatase TiO 2 have been studied for a long time, and they have been deeply revised and reanalized recently 36 .

Next, we focus on the large spectral changes induced by Ta-doping on the optical response , even for small amounts of Ta-substitution. Fig. 3 shows the optical conductivity of TiO 2 and Ta-doped TiO 2 films, at increasing Ta-doping and on a broad energy range, up to ∼35 eV. Ta x Ti 1-x O 2 films are measured for x = 0, 0.018 and 0.038. For x = 0, the pure TiO 2 sample, we observe a first sharp optical excitation at 3.48 eV (P1, see inset Fig. 3), followed by bulk resonances from 3.85 eV to 4.6 eV (P2) 36 . A well-defined large peak, at ∼6.12 eV (P3), is a newly observed intense bulk resonance in pure TiO 2 . It is followed by broad and multiple structures up to ∼35 eV.

Upon Ta-substitution (x= 0.018) we observe an emerging new giant peak at 6.0 eV (E2, Fig. 3), three times more intense than the P3 peak of pure anatase. At 6.8 eV there is an intense shoulder (E3) of the giant peak E2. The first optical excitation is also affected by Tadoping, as it occurs at higher energy ( 3.75 eV, P1 ) and reduced in intensity with respect to the pure sample. Furthermore, upon substitutional doping, the spectrum shows a significant reduction of the spectral-weight in a broad energy range ( from ∼8 eV to ∼20 eV) and a slight spectral-weight gain, singular at even higher energy (from ∼20 eV to ∼35.0 eV). For higher Ta-concentration (x=0.038), the E2 at 6.0 eV peak shows further enhanced intensity, without any significant change in the remaining structures with respect to lower Ta-concentration. To summarize, we have, upon Ta-doping: (i) an anomalous spectral-weight transfer from energies as high as 35 eV towards the 6 eV region; (ii) the emergence of a novel resonant exciton E2 at 6 eV; and (iii) the strong modification of TiO 2 optical conductivity, with an augmented optical bandgap. The optical conductivity here measured in such a broad energy range results crucial to investigate the nature of E2. Based on the optical f -sum rule, we find that the total spectral-weight (up to 35 eV) is nearly conserved for all three investigated doping ratios (Fig. 3). This directly implies that the oscillator strength at 6.0 eV is coming from spectral-weight transfers of the higher bands, i.e. from 8 to 20 eV. Such a collective spectral-weight transfer is a fingerprint of strong electronic correlations [START_REF] Eskes | Anomalous transfer of spectral weight in doped strongly correlated systems[END_REF][START_REF] Ohta | Evolution of the in-gap state in high-Tc cuprates[END_REF][START_REF] Phillips | Colloquium : Identifying the propagating charge modes in doped Mott insulators[END_REF] . Our theoretical analysis (see below) shows that the E3 peak in Ta x Ti 1-x O 2 has similar origin as the P3 peak at 6.12 eV in undoped TiO 2 , while E2 can be associated to an evolution of the bulk resonance P2. Further, based on our theoretical calculations, we could investigate the role of e-e and e-h interactions, clarifying the nature of the giant exciton E2 and confirming that observed optical behaviour is due to a manifestation of strong e-e and e-h interactions.

In Fig. 4 (a), we show the optical conductivity of the pure TiO 2 calculated using the three above mentioned theoretical methods. The GW-BSE calculations, taking into account both the e-e and e-h interactions, give rise to bound and resonant excitons or other excitonic effects along with spectral-weight transfers. The comparison between the GW-BSE and GW-RPA results gives therefore a direct measure of the excitonic nature of a resonance. The GW-RPA calculation fails in reproducing the P2 bulk resonance and the structures near the absorption edges, while both GW-BSE and JGM-TDDFT, with some differences in their detail, are able to describe the P1 and P2 peaks. This confirms that the e-e and e-h interactions are significant and important, not only for doped anatase, but even for pure TiO 2 , in agreement with previous results [START_REF] Chiodo | Self-energy and excitonic effects in the electronic and optical properties of TiO2 crystalline phases[END_REF][START_REF] Kang | Quasiparticle and optical properties of rutile and anatase TiO2[END_REF] . We plot here the optical conductivity, but we note that our GW-BSE result (HR) for optical absorption (not shown) is comparable with previous calculations for TiO 2 dielectric function [START_REF] Mak | Seeing Many-Body Effects in Single-and Few-Layer Graphene: Observation of Two-Dimensional Saddle-Point Excitons[END_REF][START_REF] Santoso | Observation of room-temperature highenergy resonant excitonic effects in graphene[END_REF][START_REF] Eskes | Anomalous transfer of spectral weight in doped strongly correlated systems[END_REF]36 .

From Ref. [START_REF] Chiodo | Self-energy and excitonic effects in the electronic and optical properties of TiO2 crystalline phases[END_REF]36 , we know that the P1 peak (experimentally at ∼3.48 eV) related to bound exciton whereas P2 (at ∼4.6 eV) comes from a bulk resonance.

The LR and HR GW-BSE data allow to properly align and identify the JGM-TDDFT spectral features with respect to experimental data, having as reference the P2 peak. The JGM-TDDFT and HR-GW-BSE coincide in intensity and energy for the peak P2, whereas, for the peak P3 and higher energy features, LR-GW-BSE calculations are in good agreement with the JGM-TDDFT. Upon Ta-substitution, the solution of the GW-BSE becomes computationally cumbersome. We turn therefore to JGM-TDDFT, which is equally reliable, as just shown in the case of pure TiO 2 , but computationally feasible also for large supercells.

We focus on optical features in the region of the ∼6.0 eV, and we use TDDFT to qualitatively study the relationship between E g (and therefore the screening properties of the material), and resonant excitonic effects in Ta x Ti 1-x O 2 . In Fig. 4(b), we show JGM-TDDFT results for increasing band-gap values, E g = 3.2, 3.5, 3.7 eV for Ta x Ti 1-x O 2 . The strong correlation mimicked by E g is reflected in the optical response, in particular in the behaviour of the peak at 6.0 eV. Peak P2 undergoes a redshift of almost 1 eV, and at the same time its intensity increases. Other features in the spectrum (as P3) undergo a similar shift, but no intensity changes are observed other than for P2.

In Fig. 5-(a), we compare the experimental findings with the theoretical calculations. Even tough the results differ in intensity, both P2 and P3 are present. Nevertheless, when the JGM-TDDFT with E g = 3.7 eV is compared with the Ta-TiO 2 optical conductivity (Fig. 5b), the theoretical calculations qualitatively suggest that peak P2 is evolving in the E2 exciton at 6.0 eV whereas the P3 peak at ∼6.1 eV is transforming in the shoulder E3 at 6.8 eV (Fig 5(b)). This seems to be counterintuitive looking at only the experimental optical conductivity results for pure and Ta-substituted TiO 2 , but becomes clear when the proper alignment and assignment of optical features are performed. From the current results, it seems that E2, evolving from P2, is indeed a resonant exciton emerging from an electron-hole continuum which exists at higher energy bands, well beyond a continuum spectrum.

Our JGM-TDDFT calculations display an interplay between E g and the resonant excitons, i.e. larger E g reflects in an enhancement of the resonant excitonic effects. Furthermore, the JGM-TDDFT calculations support the following scenario: the resonant exciton E2 at ∼6 eV in experimental spectra can be related to a modification of the electronic structure under Ta-substitution, leading also to the opening of the bandgap. This result is in contrast to a conventional picture, where Ta-substitution would lead to a simple electron doping and metallization of TiO 2 . In fact, our findings imply that Ta-substitution in Ta x Ti 1-x O 2 does not act as a conventional dopant, but plays instead an unusual role in enhancing strong electronic correlations. A behavior partaking some similarities with these results has been recently reported for a magnetic-doped TiO 2 system. Theoretical investigation of magnetic Cr-doped TiO 2 [START_REF] Da Pieve | Origin of Magnetism and Quasiparticles Properties in Cr-Doped TiO2[END_REF] shows that upon Crdoping the electronic properties undergo a transformation, and the (initially Charge Transfer insulator) system becomes a strongly correlated Mott-Hubbard crystal. We observe here an optical response which is consequent to an analogous effect in the electronic structure: the excitonic strength resulting in an enhancement of the absorption peak is occurred by the band gap opening. Upon electron doping via Ta-substitution, a possible scenario about the increasing electronic correlation is related to the Ti d-d and Ti d -O p orbital repulsions, but a more detailed analysis is left to frameworks with a better treatment of strongly correlated interactions. The scenario for Ta-doped TiO 2 optics presents conceptual similarities to strong correlated materials, as cuprates like doped La 2-x Sr x CuO 4 [START_REF] Das | Strong correlation effects and optical conductivity in electron-doped cuprates[END_REF][START_REF] Weber | Optical weights and waterfalls in doped charge-transfer insulators: A local density approximation and dynamical mean-field theory study of La2xSrxCuO4[END_REF] . In this respect, TiO 2 is widely considered an intermediate oxide between the Charge Transfer insulator and Mott-Hubbard regimes [START_REF] Zaanen | Band gaps and electronic structure of transition-metal compounds[END_REF][START_REF] Bocquet | Electronic structure of early 3dtransition-metal oxides by analysis of the 2p core-level photoemission spectra[END_REF] . We observe some analogies between our Ta-doped semiconductor spectra, and the optical behaviour of undoped Mott insulators (as cuprates), where intense optical absorption in the DUV are due to transitions from the lower to the upper Hubbard band. In the case of Mott-insulators, the change of E g (or the Mott-gap 23 ) as function of doping gives a signature of the e-e correlation. By using the Dynamical Mean Field Theory (DMFT), it has been shown that different percentages of doping enable a change of phase. The increase of d-states modifies their electronic Density Of States (DOS), and the pseudo-gap material becomes insulating. This results in a strong enhancement of peaks intensity in the optical conductivity. A similar behavior may be revealed in the present case where the inclusion of Ta d-and f-orbitals seems to have a role in changing the TiO 2 physics increasing both e-e and e-h correlations.

IV. CONCLUSIONS

In conclusion, we have presented the emergence of an intense resonant exciton induced by Ta-substitution in anatase TiO 2 . This result is of primarily importance for possible industrial applications. We argue that in these experimental findings tunable e-e and e-h correlations play a key role in the observed resonant excitons in Ta x Ti 1-x O 2 system, and can be used in a model for resonant excitonic effects. Further works will be devoted to improve our qualitative description in a more quantitative agreement with the experimental results.

It is then also important for future theoretical study to explicitly include the on-site Coulomb repulsion in the optical spectra calculations.

FIG. 1 .

 1 FIG. 1. (a) Total DOS of anatase TiO2 from DFT calculations. Black line is referred to the total DOS. Red line is the Oxygen partial DOS. (b) DOS of 6.5 % Ta-substituted anatase TiO2. The blue line is the partial DOS of Ta atom multiplied by a factor 5 to make more clear the peaks position.

FIG. 2 .

 2 FIG. 2. (a) Anatase band structure (black line), along the high-symmetry directions Z-Γ-X of Brillouin Zone. A scissor of 1.2 eV has been applied on top of PBE-DFT band structure. Ta-derived states, extrapolated from the PDOS analysis, has been superimposed on the pristine TiO2 band structure (blue dashed line). We can rule out a direct effect of Ta-doping on optical properties, as optical transitions involving Ta-states neither coincide in energy with peaks observed in optical conductivity, neither the PDOS associated to the low Ta-doping here considered could generate the intense optical features we observed. (b) A graphical representation of the long-range correlation effects of the low Ta-doping in TiO2 anatase crystal. (c) A simplified cartoon representation of the effects of Ti -Ta substitution that turns on the on-site Coulomb repulsion involving d-Ti and p-O states. The onsite repulsive interaction induces a bandgap opening, giving optical transitions above the optical gap.

FIG. 3 .

 3 FIG.3. Room temperature measurements of the real part of the optical conductivity for pure TiO2 (solid black line), 1.8% (blue short dotted line), and 3.8% Ta-substituted TiO2 (red dashed-dotted line). The polarization vector is perpendicular to the [001] direction. Inset: details of the real part of the optical conductivity in the VIS and low UV regions.

FIG. 4 .

 4 FIG. 4. (a) The calculated optical conductivity ([001] inplane polarization) for anatase TiO2 in the GW-RPA approximation (orange solid line), solving the Bethe-Salpeter equation in low resolution (LR-BSE, dark red dashed-dotted line)and High Resolution BSE (HR-BSE, dark red solid line), and within TDDFT using the JGM kernel with Eg = 3.2 eV -experimental indirect bandgap of TiO2(black solid line). The HR-BSE calculations include more stringent convergence parameters but with less included conduction bands. The LR-BSE calculations are performed to describe high energy features (beyond 6 eV). Both of them are performed to validate the JGM-TDDFT calculations. (b) Optical conductivity ([001] in-plane polarization) for anatase Tax Ti1-xO2 calculated by TD-DFT with JGM kernel and Eg of 3.2 eV (solid black line), 3.5 eV (dotted purple line), and 3.7 eV (dasheddotted red line). Increasing Eg results in an enhanced peak at 6 eV.

FIG. 5 .

 5 FIG. 5. (a) Comparison between experimental (grey ball lines) and BSE (dot-dashed dark red) and JGM-TDDFT (Eg =3.2 eV, black solid line) Optical conductivity calculations for pure anatase TiO2. A rigid shift of -0.6 eV has been applied to align the theoretical spectra to experimental data. (b) Experimental optical conductivity (pink balls) and JGM-TDDFT calculations (solid red line) for Ta-substituted ( 3.8 % TiO2 with Eg=3.7 eV). A rigid shift of -0.6 eV has been applied to align the theoretical spectra to experimental data.

V. ACKNOWLEDGMENTS

This work is supported by Singapore National Research Foundation under its Competitive Research Funding (NRF-CRP 8-2011-06), MOE-AcRF Tier-2 (MOE2015-T2-1-099), and FRC. We acknowledge the CSE-NUS computing centre, Centre for Advanced 2D Materials and Graphene Research Centre for providing facilities for our numerical calculations. We also acknowledge the National Research Foundation, Prime Ministers Office, Singapore, under its Medium Sized Centre Programme and Competitive Research Funding (R-144-000-295-281). L.C. acknowledges M. Lauricella and J. Sofo for useful discussions. Authors acknowledge F. Da Pieve for a critical reading of the manuscript. Z.Y., P.E.T. and L.C. contributed equally.