Diego Sorbelli

Paola Belanzoni
email: paola.belanzoni@unipg.it

Trond Saue
email: trond.saue@irsamc.ups-tlse.fr

Leonardo Belpassi
email: leonardo.belpassi@cnr.it

energies on AuH 2 -using state-of-the-art relativistic calculations †

a l Na me Ground and excited electronic states of AuH 2 via detachment

Gold complexes have received increasing attention [1][2][3] since the late 1980s, when Haruta and his co-workers demonstrated their catalytic potentialities 4 (for an early history of catalysis by gold see Ref. 5). In this context, gold hydrides are very important species 6 , since they are believed to be key intermediates in gold catalysis [7][8][9] . Moreover, since gold is considered as a local maxi-mum of relativistic effects in the periodic table 10 , there has been a great attention to the physical and chemical properties of gold's simplest compounds, with AuH being often a target prototypical gold complex for investigating these effects [11][12][13][14][15][16] and benchmarking relativistic quantum chemistry methods. Unfortunately, gold dihydride (AuH 2 ) has not received the same amount of attention, mainly due to the difficulties in generating this open-shell species experimentally. In its ground state the AuH 2 molecule was predicted to have a bent structure 17,18 , in contrast to its closed-shell anion (AuH - 2 ), which is linear. A potential energy surface for Au reacting with the H 2 molecule has also been reported 19 . The first experimental characterizations of the openshell molecule AuH 2 were performed in a hydrogen solid matrix by Andrews et al. 20,21 using infrared spectroscopy. These seminal works, where a joint experimental and theoretical approach was employed, confirmed that the ground state structure of AuH 2 is bent and that the system is expected to be long-lived (compared to other coinage metal compounds) because of a large barrier to dissociation into Au + H 2 . This kind of information is important since the AuH 2 molecule (together with the other coinage metal dihydrides) may be of actual interest for the development of hydrogen storage devices 20 . These works gave useful information concerning the ground state of gold dihydride and its stability. Nevertheless, not much was reported concerning its low-lying excited states. A turning point is the more recent work by Liu et al., where AuH 2 has been characterized by means of photoelectron spectroscopy (PES) in its ground and low-lying excited states 22 . PES is a widely used experimental approach for investigating heavy element compounds [23][24][25][26][27][28] . The experiment was performed in the gas phase and consisted in removing an electron from the AuH - 2 anion in its ground state by using a laser beam. In this way, the AuH 2 molecule was generated in its ground and excited states and its formation was revealed by collecting the ejected electrons. Liu et al. obtained a detailed picture of the electronic structure of AuH 2 by combining the experimental PES information with quantum mechanical calculations. First of all, the Vertical Detachment Energies (VDE) for the ground and first five excited states of the neutral AuH 2 molecule were determined with high accuracy. Furthermore, because of the high resolution of the experimental spectrum, it was possible to analyze the vibrational progression and determine the Adiabatic Detachment Energy (ADE) for the ground state of AuH 2 . The latter contains information concerning the geometrical relaxation of the system from the Frank-Condon region (i.e. geometry of the anion) to the most stable nuclear configuration of ground state of AuH 2 . Despite its apparent simplicity, the VDE is challenging to reproduce by theoretical calculations, because it requires that the anionic and neutral species are described in a well-balanced manner. In order to reproduce the experimental VDEs and to assign the bands of the PES spectrum, it is mandatory to introduce relativity (including spin-orbit coupling) and electron correlation in the calculation at a high level of accuracy. For this reason, Liu et al. resorted to a multi-step theoretical approach based on a combination of ab initio methods. The approach included scalar-relativistic effects using pseudopotentials and the spin-orbit (SO) coupling as a perturbation. In particular, the SO splitting was calculated on the basis of the complete active space self-consistent field (CASSCF) method with the diagonal matrix elements replaced by the individual state energies evaluated at the unrestricted coupled cluster level with the full single and double excitations and estimating the triples contribution via perturbation theory [UCCSD(T)]. The CASSCF calculations were performed with an active space of 13 electrons in 8 valence orbitals (Au 5d, Au 6s and H 1s) denoted CASSCF(13e, 8o). The assignment of the experimental PES bands was facilitated by the fact the AuH - 2 anion is linear with D ∞h symmetry, thus the VDE signal (associated with ground and excited states of AuH 2 ) were labeled with the spectroscopic terms of centrosymmetric linear molecules. Moreover, by analyzing the vibrational progression of the ground state PES band, it was seen how this progression corresponds to a bending mode and thus it was once again confirmed (see also the work by Andrews et al. 21 ) that the AuH 2 ground state is bent. Calculations agreed with the previous findings in locating the equilibrium bond angle at about 129 degrees. The absence of any vibrational progression in the bands from B to E (i.e. bands concerning the second to fifth excited states) led the conclusion that the corresponding states maintain a linear geometry. A more intriguing issue came out when analyzing band A, i.e. of the first excited state. Indeed, for a laser energy close to the resonance, some anomalous weak vibrational bands in the spectrum at lower energies with respect to the VDE energy suggested that the structure of the first excited state could be characterized by a double-well resulting in a minimum energy bent structure, with a very small deviation from linearity. Nevertheless, none of the calculations performed so far found a bent structure for this electronic state. Liu et al. suggested that the application of more accurate theoretical methods to treat electron correlation and relativistic effects (i.e. including SO coupling in a variational manner) would shed light on this issue and find a more stringent link between the given interpretation of the experimental data and the theoretical findings.

In this work we have employed various theoretical methods which include relativistic effects and electron correlation with increasing accuracy for the evaluation of the electron detachment energies of AuH - 2 to reproduce the experimental PES and to give detailed information about the ground and excited states of the AuH 2 molecule. The role of relativistic effects, orbital relaxation and electron correlation has been analyzed in detail. In order to include both relativistic effects and electron correlation in a consistent manner, we have also applied the recently developed and current state-of-the-art relativistic approach for the prediction of the detachment energies, the relativistic ionization-potentialequation-of-motion coupled cluster approach 29 , based on the exact two-component (X2C) Dirac-Coulomb-Gaunt Hamiltonian using a molecular mean-field approach 30 (hereafter simply referred as "X2C-EOM-CCSD"). This method has been also employed for investigating the Adiabatic Potential Energy Surfaces (APESs) of the ground and excited electronic states of AuH 2 . This will shed light on the peculiarities of gold dihydride's electronic structure and reveal intriguing features of the ground state potential energy surface in order to build a bridge between the interpretation of the experimental spectra and the theoretical calculations.

Theory and Computational details

The key quantities theoretically used to interpret photoelectron spectra and probe the electronic structure of molecules are the Detachment Energies (DE), which can be defined as:

DE = E neutral -E anion (1)
Two different types of detachment energy can be distinguished. The Vertical Detachment Energy (VDE) is defined as the energy needed to eject an electron from the anion in its ground state without relaxing its internal degrees of freedom (the geometry of the neutral species is kept frozen at the ground state geometry of the anion). The Adiabatic Detachment Energy (ADE) is defined as the minimal energy required to remove an electron from the anion ground state to produce the neutral molecule. The ADE is accordingly calculated by considering the electronic states of the neutral molecule at their minimum nuclear configurations. This is clearly lower than the VDE value because here the system is allowed to relax its geometry from the anion structure (i.e. the nuclear arrangement it has in the Frank-Condon region) to its vibronic ground state. For a stringent comparison with experimental ADE one should add to the electronic contribution also the 0-0 vibrational transition correction which corresponds to the difference in the vibrational zero-point energy between the anion and the neutral molecule.

In this work, various theoretical methods have been employed for the prediction of the VDEs. The methods have been chosen according to their computational complexity and with the specific aim of emphasizing those contributions (such as relativistic effects, orbital relaxation and electron correlation) which are required to obtain accurate VDEs [START_REF] Cederbaum | Correlation Effects in the Ionization of Molecules: Breakdown of the Molecular Orbital Picture[END_REF][START_REF] Trofimov | [END_REF] . The simplest approach is to approximate the VDEs using Koopmans' theorem (KT) 33 . This requires a simple Hartree-Fock (HF) calculation on the anion and may be considered as a useful reference since it includes neither orbital relaxation nor electron correlation (usually, these two factors tend to cancel each other when dealing with valence ionization [START_REF] Pople | Ab Initio Molecular Orbital Theory[END_REF] ). The ∆SCF approach [START_REF] Bagus | [END_REF] has been applied to include the effect of the orbital relaxation (i.e. the response of the electronic structure to a hole created in the occupied orbitals [START_REF] Trofimov | [END_REF] ).In this case VDEs are calculated as in Eq.1 where the energy of the neutral molecule (ground and excited states) are computed by forcing the occupation of the orbitals and then by self-consistently solving the HF or the Kohn-Sham (KS) equation with this excited state population. In those cases in which this procedure properly converges, this approach remains computationally very cheap with the advantage that the orbital relaxation is included in the calculations. With the ∆SCF approach electron correlation can be neglected (in the case of ∆HF) or included through an exchangecorrelation functional at the Kohn-Sham level (∆KS). Reliable results for the VDEs can be obtained within the framework of the Ionization-Potential-Equation-of-Motion Coupled Cluster (EOM-IP-CCSD) 36,37 method, which is known to include in a consistent manner both orbital relaxation and electron correlation. The method has a computational cost comparable to that of the CCSD method for the ground state (scaling is approximately N 6 ) 36 . We refer the interested reader to a few pedagogical reviews which explain how the EOM (and EOM-IP) theory works and illustrate its applicability on typical target systems [38][START_REF] Watts | An Introduction to Equation-of-Motion and Linear-ResponseCoupled-Cluster Methods for Electronically Excited Statesof Molecules[END_REF][START_REF] Bartlett | [END_REF][41] .

Unless otherwise stated, all the reported calculations have been carried out using the relativistic quantum chemistry code DIRAC [START_REF] Dirac | a relativistic ab initio electronic structure program, Release DIRAC18[END_REF][START_REF] Saue | [END_REF] (release 2018) which implements all the methods mentioned above in combination with accurate relativistic Hamiltonians. In particular, in order to investigate the impact of relativistic effects on VDEs we employed exact 2-component [START_REF] Aa | Douglas-Kroll the Easy Way[END_REF][START_REF] Kutzelnigg | [END_REF][46][47][48] (X2C) and 4-component (4c) Hamiltonians, with the latter also in its spinfree (SF) form 49 .

Dyall's Gaussian-type basis sets of different qualities, double-, triple-, and quadruple-zeta with polarization functions (dyall.dzp, dyall,tzp, dyall.qzp respectively) for both gold 50 and hydrogen 51 were employed for the KT, ∆HF and ∆KS calculations. For the ∆KS case the PBE 52 exchange-correlation functional was applied. The SCF convergence on the excited state of the AuH 2 molecule was obtained by overlap selection. [START_REF] Dirac | a relativistic ab initio electronic structure program, Release DIRAC18[END_REF][START_REF] Saue | [END_REF] Correlated EOM-IP-CCSD calculations are based on the recent relativistic formulation 29 and have been carried out using both Dyall valence (dyall.vdz, dyall.vtz, dyall.vqz) and core-valence (dyall.cvdz, dyall.cvtz, dyall.cvqz) basis sets 50,51 . For the former calculations the 5s, 5p, 5d and 6s electrons for Au are correlated, whereas for the latter the n=4 shell electrons of Au were included in the active space as well. The X2C-Dirac-Coulomb-Gaunt (DCG) Hamiltonian has been employed in the molecular mean field (mmf) approach 30 . This method (in the following simply referred as "X2C-EOM-CCSD") represents an excellent compromise between a good accuracy and a reasonable computational cost; the latter has been shown to yield results largely indistinguishable from their 4-component counterparts 29 . The X2C-EOM-CCSD approach has been used for the numerical optimization of the minimum energy structures of the AuH 2 molecule in its ground and (five) low-lying excited states and to investigate selected cuts of Adiabatic Potential Energy Surfaces (APESs). The calculation of the expectation values of the electric dipole moment for the ground state of AuH 2 was carried out using analytical gradients as implemented in Ref. 53. For calculating the VDEs, the geometry from Ref. 22

has been used for AuH

- 2 (r Au-H = 1.647 Å; α H-Au-H = 180 degrees).
We mention that some test calculations have been also applied including the relativistic effects through the approximate scalar (SR) and spinOrbit (SO) 2-component ZORA Hamiltonians [54][55][56] , for which the Amsterdam Density Functional (ADF) program package [57][58][START_REF] Baerends | ADF2014[END_REF] was used, together with Slater-type double-(DZ), triple-(TZ2P) and quadruple-zeta (QZ4P) basis sets. In particular this code has been used to estimate the zero-point-energy (ZPE) and the 0-0 vibrational correction to the ADE.

Results and Discussion

Vertical Detachment Energies (VDEs)

As mentioned, VDEs obtained with high accuracy in the experiment by Liu et al. 22 represent an ideal reference in order to establish a suitable theoretical protocol, based on first principles, for the quantitative estimation of these observables. Our results are summarized in Figure 1, where the VDEs are reported using Vertical Detachment Energy (eV) different theoretical methods and compared with the experimental data taken from Ref. 22. The numerical data are also reported in Table 1, together with the percent error with respect to the experimental data. All calculations have been carried out at the linear equilibrium geometry of the anion AuH - 2 , taken from Ref. 22 (r Au-H = 1.647 Å; α H-Au-H = 180 degrees). For the sake of simplicity, we keep the same notation used in Ref. 22 and label the VDEs using both the experimental bands (X, A, B, C, D and E in energy order) and the spectroscopic terms of the linear centrosymmetric molecules reported as 2Σ+1 Λ Ωu/g . It should be noted that the first two transitions (corresponding to the ground and first excited state of AuH 2 ) are challenging to assign unambiguously and represent a stringent test for our calculations (these two states are of different parity and are very close in energy, the experimental energy difference is of only 0.24 eV). Now, focusing on the results, we begin with the analysis of the VDEs values based on Koopmans' theorem (KT). As shown in Figure 1 and Table 1, the VDEs have been calculated both at the 4c level (4c-KT) and in a spinfree framework (SF-KT). Similar data can be obtained with the 2-component scalar and spinOrbit ZORA Hamiltonian, see Table S1 in Electronic Supplementary Information (ESI †). The results show that the inclusion of SO coupling is fundamental even for a qualitative description of these quantities. Indeed, concerning the higher-energy detachments (i.e. bands from B to E), when SO coupling is included, the experimental VDEs trend is nicely reproduced. The method is able to give the correct assignment of these electronic states even if the calculated VDEs are shifted at somewhat higher energy. In a SF framework, a completely different trend is observed. In this case, the C/D states (and the B/E states) become degenerate since they correspond to the doublet components of the 2 Π states (and 2 ∆ states). This clearly prevents the proper energy order and the correct assignment states. Although the 4c-KT method has only to be considered as a first approximation it is noteworthy that, the calculated energy splitting for the C/D and B/E states (numerical values are 1.05 eV and 1.86 eV, respectively) agrees surprisingly well with the experimental values (0.96 eV and 1.59 eV). The inclusion of SO coupling has a positive impact on the results for high energy states (i.e. bands from B to E), nevetheless, the situation is different at lower energies (X and A bands). The highest level calculations from Ref. 22 showed that the band X should be related to the electronic state of the ungerade symmetry, corresponding to the 2 Σ + 1/2u state. In our 4c-HF calculations, the AuH -

2

HOMO is an orbital of symmetry σ 1/2g , whereas the HOMO-1 has a σ 1/2u symmetry, thus the application of KT gives a ground state of the incorrect symmetry ( 2 Σ + 1/2g ). Indeed, there is a swap in energy between the HOMO and the HOMO-1 molecular orbitals of the anion. The 4c-KT method is clearly a too simple approach for an accurate estimation of DEs.

For the evaluation of the orbital relaxation effect, we use the 4c-∆HF approach. The experimental trend is qualitatively reproduced for bands B-E, although the calculated VDE values are now shifted at lower energies; nevertheless the 4c-∆HF method still fails in providing the correct energy order of the two lowestenergy states (bands X-A). In general, the error in the higher energy states is increased dramatically, probably because the error compensation (between the orbital relaxation and electron correlation contributions), typically expected in the KT framework, is now missing. The 4c-∆KS calculations show how including both orbital relaxation and electron correlation significantly improves the agreement with the experimental values. Unfortunately, the energy swap of the lower-energy states still remains and the method fails to reproduce the correct symmetry of the AuH 2 ground state. This finding suggests that a more accurate treatment of electron correlation is needed for reproducing these quantities. The results reported in Figure 1 and Table 1 show that the X2C-EOM-CCSD method works exceptionally well. First of all, the order of the electronic states found here fully supports the assignment reported in Ref. 22, where, as mentioned in the introduction, a multi-step theoretical recipe was specifically applied to include scalar relativistic effects, electronic correlation and spinorbit coupling in a balanced way. Our results unambiguously confirm that the ground state of AuH 2 at the linear centrosymmetric configuration (i.e. at the minimum energy of the anion) is the electronic state 2 Σ + 1/2u (ungerade symmetry). The computed VDEs are in very good agreement with the experimental data; indeed, as it can be seen, the percent errors for calculated VDEs are in the range 2.9-0.3% (Table 1). In addition to this, it must be pointed out that the results displayed in Table 1 can be consistently improved. Indeed, as shown in Table 2, one can increase both the quality of the basis set and the number of active electrons. It is evident from the data how increasing simultaneously both of these factors can lead to an even more stringent agreement with the experimental results. For instance, when electrons of the shell with n=4 of gold are correlated and a core-valence quadruple-zeta (cvqz) quality basis set is used, the error can be further reduced (the percent errors are in the range 0.1-1.4%, which corresponds to an absolute error below 0.06 eV). The error on the VDEs for the experimental reference values was estimated in the order of 0.02-0.03 eV 22 . The results are consistent (with only a slightly increase of the percent error), also if we apply the basis set extrapolation techniques (data labeled as v∞z and cv∞z in Table 2). In particular we have used the three-parameter exponential formula (see Eq. 20 in Ref. 32) already employed to estimate the ionization energies complete basis set limit [START_REF] Trofimov | [END_REF] .

Adiabatic Potential Energy Surfaces (APESs) of the ground and low-lying excited states of AuH 2

At this point, we already showed that the X2C-EOM-CCSD approach is a suitable and accurate theoretical tool for quantitatively probing the electronic structure of AuH 2 . Our principal aim here is to characterize in detail the ground and low-lying electronic excited states of AuH 2 , including their minimum energy nuclear configurations, and provide accurate ADEs. The latters are typically more difficult to extract from experiments than VDEs. Furthermore, we will analyze the first excited state of AuH 2 , which was suggested by Liu et al. 22 to possess a double well in the direction of the bending mode that would be consistent with a slightly bent structure, but that previous theoretical calculations were not able to predict. In their seminal work Liu et al. suggested that only those theoretical methods that included relativistic effects and electron correlation at a high level of accuracy could be useful for optimizing the geometry of this state, and therefore, for fully rationalizing the experimental data. The X2C-EOM-CCSD approach seems particularly suitable for answering this open question. Since an analytic gradient is not available in the current implementation of the method, we probe the potential energy surfaces resorting to single point energy calculations along selected cuts of the APES. We have preliminary optimized the Au-H bond length in a linear and centrosymmetric configuration (the results are reported in Table S2 in the ESI †). The results show that the optimal bond length is about 1.64 Å, almost independently of the specific electronic state (it ranges from 1.64 to 1.66 Å, for 2 Σ + 1/2u and 2 Π 3/2g , respectively). Thus, by freezing the bond length, we have performed a first scan in the bending direction (see Figure 2). Numerical data are reported in Table S3 of ESI †. In the ESI †we also report the results obtained using a small basis set (dyall.vdz) for comparison (see Figure S1 and Table S4 in the ESI †). Our results for the ground state matches completely with the experimental data and confirms previous calculations 17,21,22 : this state is found to be bent and the minimum bond angle has been calculated very close to 130 degrees (the system has a C 2v symmetry and using the spin-free notation we can label this electronic state as 2 B 2 ). For the second to fifth excited states, the photoelectron spectrum showed no vibrational progressions for the corresponding detachment bands and for this reason these states are expected to be linear (they maintain the linear structure as that of the anion and so they have a small nuclear relaxation). Our calculations confirm this picture and these states present a minimum at 180 degrees (the corresponding spectroscopic terms are those of the centrosymmetric linear configuration: 2 ∆ 5/2g , 2 Π 3/2g , 2 Π 1/2g and 2 ∆ 3/2g ). The interpretation given in the PES experiment suggested a slightly bent structure for the first excited state state, nevertheless, as can be seen clearly in Figure 2, this excited state of gold dihydride is found to be linear. Even scanning the adiabatic potential energy surface finely in the vicinity of the 180 degrees does not modify this picture: the minimum energy of the first excited state is found in its linear configuration. We mention that this result does not qualitatively depend on the quality of the basis set employed (for a comparison using different basis sets see Figure S2 and Tables S5-S6 in the ESI †). One possible pitfall of the above calculations Table 2 VDEs (in eV) for AuH - 2 calculated with the X2C-EOM-CCSD method by varying both the number of correlated electrons (i.e. only 5d and 6s electrons, all n=5 shell and 6s electrons or all n=4 and n=5 shells plus 6s electrons) and the basis set (valence and core-valence Dyall basis set of double-, triple-and quadruple-zeta quality). Complete basis set extrapolation data are also reported (see text for details). For each calculated value, the percent error with respect to the experimental reference is reported in parenthesis. Experimental data taken from Ref. 22 S3). Based on these findings, we can confidently say that the X2C-EOM-CCSD approach predicts a linear geometry for the first excited state of gold dihydride, which somehow contrasts with the conclusions derived by Liu et al. 22 based on the analysis of the PES spectra. We will return on possible reasons of this apparent inconsistency at the end of this section.

By performing the scan of the APES along both the bond length (maintaining the constraint that the two Au-H remain of the same length) and the bond angle at X2C-EOM-CCSD level, analogously to what we have done for the first excited state, we have found the equilibrium geometries and the ADEs for all electronic states associated with the experimental bands. The results are reported in Table 3. Concerning the AuH 2 ground state, the geometry found here matches quite well the one previously calculated by Liu et al. 22 with a bond length of 1.605 Å and a bond angle of 129.2 degrees. Moreover, the ADE value is very close to the one that was determined experimentally also reported in Table 3. We mention that our theoretical estimate of ADE does not include the 0-0 vibrational correction (the difference in the ZPE of the anion AuH - 2 and of AuH 2 at the respective equilibrium geometries). This contribution is relatively small an estimate at DFT level gives about 0.03 eV (see Table S7 in ESI †). Concerning the excited states, even though we have no experimental reference to validate our results, we find that ADEs have only a slightly smaller values compared to the VDEs values. This is consistent with the small geometrical rearrangements of AuH 2 in these electronic states with respect to the anion. The difference between VDE and ADE (∆ V DE-ADE ), can be considered as a stringent test for the ability of our calculations to describe the potential energy surface in a well-balanced and consistent manner. We find that the calculated ∆ V DE-ADE (0.687 eV) agrees extremely well with the experimental one (0.637 eV). The difference is 0.05 eV, which further demonstrates the consistency of the calculated APESs.

All the results reported above witness the accuracy of the X2C-EOM-CCSD computational protocol in describing the features of gold dihydride's electronic states, in close agreement with the experiment. The only discrepancy that still remains concerns the first excited state, which was expected to be characterized by a slightly bent equilibrium geometry. Returning to the experimental work of Liu et al. 22 , we have to observe that the bent structure for the first excited state was drawn on the basis of the presence of a fine structure PES spectrum that emerges at slightly lower energies with respect to the signal assigned to the VDE of the first excited state (band A, in Figure 3 and Table 1 of Ref. 22). These fine band structures are made of three or four vibrational states spaced from about 500 to 200 cm -1 with an increasing anharmonicity and were associated with some bound nuclear states of the bending mode. However, if this would be the case, the energetic range would be consistent with the presence of a double-well in the excited state with a barrier larger than 1600 cm -1 (about 0.2 eV), which is clearly far beyond the accuracy we expect from our level of theory and thus, it would have been easily revealed by our calculations. So, we have to conclude that, on the basis of our extensive theoretical study, the interpretation for which the first excited states should present a bent structure needs to be revised. Furthermore, as a matter of fact the vibrational progression, previously assigned to the first excited state of AuH 2 observed in the experiment occurs at an energy which is between the VDEs of the ground and the first excited state (3.67 and 3.90 eV, respectively) and has a spread of about 0.2 eV, which is really consistent with the small energy difference between these two quasi-degenerate states. We therefore suggest that the mentioned vibrational progression is not a sign of the presence of a bent structure of the first excited state, but may be more likely related to the dynamical effect of the nuclei, that arises from presence of quasi-degenerate adiabatic states (for which the Born-Oppenheimer approximation cannot be fulfilled) [START_REF] Köppel | Multimode Molecular Dynamics Beyond the Born-Oppenheimer Approximation[END_REF] . In such situations, the fine details of the PES spectra can only be extracted using appropriate nuclear dynamics coupling models based on the results of high-level electronic structure methods [START_REF] Trofimov | [END_REF][START_REF] Köppel | Multimode Molecular Dynamics Beyond the Born-Oppenheimer Approximation[END_REF] . Precisely because the non-adiabatic dynamics is expected to be very sensitive even to the fine details of the adiabatic potential energy surfaces, especially in the Frank-Condon region, it is highly desirable to investigate whether the presence of these quasi-degenerate states may give peculiar and unexpected features of their topology. In particular, in the following section we will demonstrate that the pseudo-Jahn-Teller effect [START_REF] Öpik | [END_REF]62 gives arise to a symmetry breaking in the adiabatic potential energy surface of the ground state, thus that, right in the Frank-Condon region (i.e. in the linear arrangement), the most stable configuration of AuH 2 is a non-centrosymmetric structure with unequal Au-H bond lengths.

Symmetry breaking and Pseudo-Jahn-Teller effect in the AuH 2 linear configuration

A number of linear triatomic open-shell molecules, such as ABA systems (with A and B from IA and VIIA groups respectively) 63 , ZnCl + 2 64,65 , CuCl 2 66 , and the more debated BNB 67 present a linear non-centrosymmetric structure. The distortion from the centrosymmetric symmetry of these linear molecules arises from the symmetry breaking (SB) consequence of the pseudo-Jahn-Teller effect (PJTE), where a significant vibronic coupling is present in near-degenerate electronic states [START_REF] Öpik | [END_REF]62 . In the ESI †, we briefly reviewed the main points of the theory behind the PJTE that may be useful for discussing our results. For further insights into PJTE theory, the reader may consult Refs. 61-63,68 and, for an even more specific analysis of PJTE in the linear X-Y-X systems, sections 3.1 and 3.2 of Ref. 62. For the PJTE treatment by means of quantum chemistry methods (including the EOM-IP-CCSD method) we refer the reader to a detailed discussion reported in Ref. 69. Other relevant works showing the specific application of the EOM-IP-CCSD method for describing PJTE, open-shell systems and quasidiabatic states are Refs. [START_REF] Stanton | [END_REF][71][72] .

The PJTE takes place only when the two (quasi-degenerate) electronic states have symmetry properties which are compatible with the symmetry of the nuclear displacement. In particular, the necessary condition in order to have the PJTE is that Γ = Γ × Γ " , where Γ is the irreducible representation of the nuclear displacement and Γ and Γ " are the irreducible representations to which the lower and higher energy coupled states belong in the high symmetry configuration. The other important factor to have a significant PJTE is electronic coupling parameter (∆), which is defined as half of the energy separation between the two coupled states in their high symmetry configuration. The closer these states are in energy (small ∆), the stronger the vibronic coupling is and an actual SB may occur. For gold dihydride in its linear configuration (D ∞h symmetry) we have that the ground ( 2 Σ + 1/2u ) and first excited state ( 2 Σ + 1/2g ) are very close in energy and are likely to undergo symmetry breaking along the coordinate of the asymmetric stretching (Σ + u symmetry). Note that the necessary condition in order to have the PJTE is fulfilled whether one uses the spin-free notation or the full picture including SO coupling (double group symmetry). In the spin-free framework one only considers spatial symmetry and the direct product of the involved electronic states (Σ + g × Σ + u ) gives Σ + u , which is the symmetry of the antisymmetric stretch. In the full picture the direct product of electronic states (E 1/2,g × E 1/2,u ) gives Σ + u + Σ - u + Π u which also includes the symmetry term of the antisymmetric stretch. For an interesting discussion about the impact of SO coupling on the selection rules for PJTE see Ref. 73. Since in our case the PJTE can take place by symmetry, the system may be stable in a configuration with unequal Au-H separation, and unstable (saddle point) in the centrosymmetric configuration. At first glance, this seems to be an exotic situation; however, as already mentioned above, it has been observed for a number of linear triatomic open-shell molecules of general formula X-Y-X [62][63][64][65][66][67] .

So we have used the X2C-EOM-CCSD approach for investigating the APES along the asymmetric stretching coordinate. The method is expected to be accurate to describe the PJTE 38,[START_REF] Krylov | Reviews in Computational Chemistry[END_REF][START_REF] Stanton | [END_REF] . The CI-like form of the EOM excitation operator enables access to multi-configurational open-shell wave functions, including degenerate and nearly degenerate states, as well as interacting states of different character 38 . The notorius issue of coupledcluster theory to treat conical intersections [74][75][76] between states of the same symmetry is not expected to play a role here. The numerical examples by Kohn et al. 74 showed that this problem arises only very close to the intersection (within an energy range of ±0.02 eV). This is clearly not the case in point here, in that the involved electronic states have different symmetry for the high symmetry nuclear configuration and their energy gap is above 0.15 eV. For sake of completeness, we mention that the problem of EOM with treating conical intersections between states of the same symmetry has recently found a solution within the framework of similarity-constrained theory 76 . The selected cut of potential energy surface has been worked out by scanning several distances between the two terminal hydrogen atoms (R H-H ) and by varying the single Au-H bond lengths asymmetrically (the asymmetric displacement with respect to the centrosymmetric Surprisingly, the results support our hypothesis: the ground state APES of gold dihydride is not centrosymmetric around its linear configuration (as demonstrated by the lower sheet in Figure 3, showing a double-well APES for the ground state), whereas the first excited state is linear and centrosymmetric. We find the minimum equilibrium structure to be at R H-H = 3.3 Å with the minimum energy asymmetric displacement being ∆R = 0.06 Å. For a more quantitative perspective on the energies, we also performed calculations with larger basis sets (i.e. vtz, vqz, with n=5 and 6s electrons correlated) and by correlating more electrons (i.e. cvtz, cvqz basis, with n=4, n=5 and 6s electrons correlated) by keeping fixed the overall R H-H distance at 3.3 Å. As shown in Table 4, a SB is always observed, with the minimum energy asymmetric displacement being ∆R=0.06 Å. Both the PJTE stabilization energy separating the minimum energy asymmetric structure from the centrosymmetric one (E PJT E ) and the electronic coupling parameter ∆ change, becoming smaller and greater respectively; indeed, the most accurate calculations we performed (i.e. with a cvqz basis set) yield a E PJT E of 31 cm -1 , which is small, yet is perfectly compatible with previous data referred to other systems in which the same phenomenon takes place 62 . One of the most important and straightforward parameters for the evalua-tion of the PJTE (see Eqs. 3-6 in the ESI †) and SB is the electron coupling parameter, ∆, which can be directly evaluated from the APES. We see, in Table 4, how actually this value is related with the depth of the double-well: when, with smaller basis-sets, we have a greater PJTE stabilization energy, we observe a smaller ∆ and vice versa with larger basis set. One may argue that this symmetry breaking could be an artefact of the calculation; nevertheless, we have means to point out that we avoid artificial SB. First of all, as already mentioned, the X2C-EOM-CCSD method, which is based on EOM-IP-CCSD, is an accurate approach since it recovers most of the necessary correlation needed for treating these phenomena 37,38,65,[START_REF] Krylov | Reviews in Computational Chemistry[END_REF][START_REF] Stanton | [END_REF] .

Furthermore, as commonly used in the literature, in order to demonstrate that the SB is a real effect and not an artefact of the calculations is to ascertain that wavefunction modifies its character continuously as the geometrical parameters are modified along a symmetry breaking mode. Thus, for the linear AuH 2 molecule, we calculate the expectation values of the dipole moment for the ground state at different asymmetric displacement (∆R). A real SB should feature a zero dipole moment when ∆R=0 and a continuous curve for other values of ∆R 65,67 (reflecting a wavefunction which varies continuously as a consequence of the physical mixing of the coupled electronic states). This criterion is fully respected by our calculations: the dipole moment expectation value is zero at ∆R=0 and as the geometry undergoes distortion the corresponding dipole varies with continuity as shown in the previous cited examples 65,67 . The numerical data are reported in the ESI †, see Table S8.

We have shown that the PJTE is very sensitive to the energy separation of the electronic states at their highest symmetric configuration (∆ parameter). Since AuH 2 is an heavy-element containing molecule (to our knowledge the first that displays a SB in the asymmetric stretching direction due to PJTE), it is interesting to investigate if the distortion in the asymmetric direction can be affected by spin-orbit coupling. The answer to this question can be given by scanning the APESs with the EOM-IP-CCSD approach, using a spinfree Hamiltonian in order to switch off spinorbit interaction. The spinfree APESs are reported in Figure 4. Figure 4 shows how, by using a spinfree Hamiltonian in the calculations, there is no evidence of a SB in the asymmetric stretching direction. Indeed, the APES for the ground state results to be very flat, but both the ground and first excited state's APESs show a minimum at ∆R=0, i.e. the minimum configuration is the centrosymmetric one when spin-orbit coupling is neglected. This finding is surprisingly in that the symmetry analysis just reported above shows that PJTE is also allowed also within the spinfree framework. A more quantitative perspective can be obtained by considering the electronic coupling parameter (and therefore the separation of the two states at ∆R=0, 2∆). Indeed, as shown in Table 4, where relativistic effects (including spin-orbit coupling) are considered, the separation between these ground and first excited state in their linear nuclear configuration (when using the same computational setup, i.e. dyall.vtz basis set) is less than 1000 cm -1 (2∆ = 932 cm -1 ). In the case of spinfree calculations instead, switching off the spin-orbit interaction, the separation is much larger and almost the double (1724 cm -1 ). This is consis- tent with the PJTE theory, where to larger ∆ should correspond weaker couplings between states.

Conclusions

A very detailed photoelectron spectroscopy (PES) experiment was previously reported, in which the electronic structure of AuH 2 was extensively probed. Detachment energies (DEs) were determined with high accuracy and, coupled with ab initio calculations, the ground state of AuH 2 was found to be bent. The first excited state was also supposed to be slightly bent due to the presence in the spectrum of a small vibrational progression which was associated with a bending mode, although, to this day, every theoretical approach tested yielded a linear structure. In this work, based on this highly accurate experiment as reference, we have applied several methods for reproducing the experimental DEs of AuH 2 , by taking into account the factors that may influence calculations, such as relativistic effects, orbital relaxation and electron correlation. The X2C-EOM-CCSD method, which includes the electronic correlation at the EOM-IP-CCSD level and relativistic effects based on the X2C-Dirac-Coulomb-Gaunt Hamiltonian using a molecular mean-field approach is found to be highly accurate in both reproducing the correct order of the experimental PES bands and the DEs values. The method has been used for describing the Adiabatic Potential Energy Surfaces of gold dihydride's ground and a few low-lying excited states. For the ground state, for which detailed information concerning its structure and DEs was available, the X2C-EOM-CCSD approach yields a very accurate adiabatic potential energy surface, whose features tightly match the experimental findings. Indeed, both the bending angle and the Adiabatic Detachment Energies (ADEs) were reproduced with high accuracy (typical error is of the order 0.03-0.06 eV). All the other excited state (including the first) are found to be linear. The picture that emerges from our theoretical investigation is intriguing and gives a new perspective on the features of this peculiar gold complex.

In particular, on the basis of our X2C-EOM-CCSD calculations we have demonstrated that the interpretation which suggested for a slightly bent structure of the first excited state of AuH 2 , needs to be revised. The first excited state in AuH 2 is found to be linear and, in this nuclear configuration (which also corresponds to the equilibrium geometry of the anion, AuH - 2 and characterizes the vertical detachment process), it is quasi-degenerate with the ground electronic state. Surprisingly, these states are found to be coupled by the pseudo-Jahn-Teller effect (PJTE) which induces a symmetry breaking along the asymmetric stretching coordinate, for which the centrosymmetric nuclear configuration becomes a saddle point and asymmetric configuration (two Au-H bonds with different lengths) is the most stable. Even though similar triatomic open-shell molecules have been reported to be affected by this phenomenon, this is, to our knowledge, the first reported example of an heavy-metal-containing molecule affected by PJTE in the direction of the asymmetric stretching. Since AuH 2 is an heavy metal containing compound, the role of the spin-orbit coupling has been analyzed. Its effect is to reduce the energy gap between these quasi-degenerate states and enhancing the PJTE. Interestingly, thus, we have shown how the symmetry breaking is observed only when the calculations accurately account for spinorbit coupling. Our findings strongly suggest that the mentioned vibrational progression observed in the reference experiment is not a sign of the presence of a bent structure of the first excited state, but, more likely, may be related to the dynamical effect of the nuclei, due to the presence of quasi-degenerate adiabatic states (for which the Born-Oppenheimer approximation may not be fulfilled) in the Frank-Condon region associated with the photodetachment process. It is clear that in order to confirm the last hypothesis a full treatment of the nuclear degrees of freedom would be required. The latter, together with the surprising symmetry breaking reported here paves the way for further theoretical and experimental investigation on the peculiar electronic structure of AuH 2 , which is a simple molecular system only in appearance.

Fig. 1

 1 Fig.1Gray and green diamonds correspond to VDEs estimated using Koopmans' theorem (KT) at the 4-component (4c-KT) and spin-free (SF-KT) levels, respectively. Black and orange crosses correspond to the 4c-∆KS and 4c-∆HF calculations, respectively. Basis set for all the latter is dyall.tzp. Blue dots represent X2C-EOM-CCSD calculations using the dyall.vtz basis set. Red bars represent experimental values from Ref.22.

Fig. 2

 2 Fig.2X2C-EOM-CCSD calculated cuts of the APESs in the bending direction for AuH 2 (the Au-H bond length frozen at 1.64 Å) in its ground and first five excited states. Basis set: dyall.vtz. The energy has been shifted relatively to the minimum energy. The electronic states are labeled with spectroscopic terms of the centrosymmetric linear configuration corresponding with the bending angle of 180 degree.

Fig. 3

 3 Fig.33D plot of the ground (lower sheet) and first excited (upper sheet) state APESs of AuH 2 in the direction of the asymmetric stretching. Calculations have been performed at the X2C-EOM-CCSD level. Basis set: dyall.vdz. The energy has been shifted relatively to the minimum.

Fig. 4

 4 Fig.43D plot of the ground (lower sheet) and first excited (upper sheet) state APESs of AuH 2 in the direction of the asymmetric stretching. Calculations have been performed at the EOM-IP-CCSD level by switching off the spin-orbit coupling effect with a spinfree Hamiltonian. Dyall.vtz basis set was used. The energy has been shifted relatively to the minimum.

Table 1

 1 VDEs (in eV) for AuH - 2 calculated at different levels of theory. Basis set for KT and ∆HF(KS) calculations is dyall.tzp, whereas for X2C-EOM-CCSD calculations is dyall.vtz. For each calculated value, the percent error with respect to the experimental reference is reported in parenthesis. a Experimental VDEs (Exp.) are taken from Ref.22.

	State/Exp. band	SF-KT	4c-KT	4c-∆HF	4c-∆KS	X2C-EOM-CCSD	Exp. a
	2 Σ + 1/2u /X 2 Σ + 1/2g /A	3.982 (8.3) 3.676 (5.8)	4.021 (9.3) 3.599 (7.8)	3.386 (7.9) 2.299 (41.1)	4.013 (9.1) 3.848 (1.4)	3.667 (0.3) 3.790 (2.9)	3.678 3.904
	2 ∆ 5/2g /B	6.398 (38.1)	5.699 (22.9)	2.900 (37.4)	4.850 (4.6)	4.731 (2.1)	4.635
	2 Π 3/2g /C	6.640 (38.8)	5.878 (22.8)	3.021 (36.9)	4.823 (0.8)	4.825 (0.8)	4.785
	2 Π 1/2g /D	6.640 (15.4)	6.933 (20.5)	4.128 (28.3)	5.832 (1.5)	5.765 (0.3)	5.745
	2 ∆ 3/2g /E	6.398 (2.9)	7.566 (21.5)	4.477 (28.0)	6.450 (3.7)	6.275 (0.8)	6.22

  (the numerical values are also reported in Table1of this work for easy reference).

	State/exp. band					Active electrons and basis set
			5d+6s			5s5p5d + 6s			4s4p4d5s5p5d +6s
		vdz	vtz	vqz	v∞z	vdz	vtz	vqz	v∞z	cvdz	cvtz	cvqz	cv∞z
	2 Σ + 1/2u /X	3.454 3.667 3.714 3.727	3.458	3.673 3.725 3.742	3.462	3.678 3.727 3.741
		(6.1)	(0.3)	(1.0)	(1.3)	(6.0)	(0.1)	(1.3)	(1.7)	(5.9)	(0.0)	(1.3)	(1.7)
	2 Σ + 1/2g /A	3.574 3.790 3.873 3.925	3.462	3.734 3.839 3.905	3.478	3.748 3.849 3.909
		(3.8)	(2.9)	(0.8)	(0.5)	(11.3)	(4.3)	(1.7)	(0.2)	(10.9)	(4.0)	(1.4)	(0.1)
	2 ∆ 5/2g /B	4.457 4.731 4.858 4.968	4.126	4.496 4.667 4.814	4.142	4.513 4.667 4.814
		(3.8)	(2.1)	(4.1)	(7.2)	(11.0)	(3.0)	(0.7)	(3.9)	(10.6)	(2.6)	(0.9)	(3.9)
	2 Π 3/2g /C	4.551 4.825 4.950 5.055	4.238	4.595 4.771 4.942	4.254	4.619 4.781 4.910
		(4.9)	(0.8)	(3.4)	(3.8)	(11.4)	(4.0)	(0.3)	(3.3)	(11.1)	(3.6)	(0.1)	(2.6)
	2 Π 1/2g /D	5.500 5.765 5.884 5.981	5.259	5.593 5.764 5.943	5.274	5.621 5.773 5.891
		(4.3)	(0.3)	(2.4)	(4.1)	(8.5)	(2.6)	(0.3)	(3.4)	(8.2)	(2.1)	(0.5)	(2.5)
	2 ∆ 3/2g /E	5.999 6.275 6.404 6.521	5.731	6.097 6.276 6.450	5.746	6.118 6.286 6.396
		(3.5)	(0.8)	(2.9)	(4.8)	(7.9)	(1.9)	(1.0)	(3.7)	(7.6)	(1.6)	(1.1)	(2.8)

may be the fact we used a frozen Au-H bond length (fixed it at 1.64 Å). In order to relax also this last constraint, a larger 3D scan of the first excited state has been performed, where bond length (Au-H) and bond angle (H-Au-H) have been varied simultaneously. The results are shown in the ESI †(see Figure

Table 3

 3 X2C-EOM-CCSD calculated equilibrium geometries (the Au-H bond length in Å) and Adiabatic Detachment Energies (ADEs) (in eV) for the ground and first five excited states of gold dihydride. Dyall.vtz basis set was used. The equilibrium bond lengths (Å) and angles (degrees) have been determined through polynomial fits. a Experimental data taken from Ref. 22 State/PES band Bond angle Bond length Calculated ADE Exp. ADE a

		(degrees)	(Å)	(eV)	(eV)
	2 B 2 /X 2 Σ + 1/2g /A	129.8 180.0	1.59 1.63	2.98 3.72	3.030 -
	2 ∆ 5/2g /B	180.0	1.64	4.49	-
	2 Π 3/2g /C	180.0	1.65	4.60	-
	2 Π 1/2g /D	180.0	1.66	5.60	-
	2 ∆ 3/2g /E	180.0	1.64	6.10	-

Table 4

 4 PJTE stabilization energy separating the centrosymmetric configuration and the distorted one (E PJT E ), electronic coupling parameter (∆) and asymmetric displacement at the minimum energy configuration (∆R) for gold dihydride ground state calculated at the X2C-EOM-CCSD level with different Dyall valence (dyall.vdz, dyall.vtz, dyall.vqz) and corevalence (dyall.cvtz, dyall.cvqz) basis sets.

	Basis set	vdz	vtz	vqz	cvtz	cvqz
	E PJT E (cm -1 )	268	126	34	150	31
	∆ (cm -1 )	23	221	466	186	472
	∆R (Å)	0.06	0.06	0.06	0.06	0.06
	structure is indicated as ∆R).

1-11 J o u r n a l N a me , [ y e a r ] , [ v o l . ] ,

[1][2][3][4][5][6][7][8][9][10][11]

J o u r n a l N a me , [ y e a r ] , [ v o l . ] ,

Acknowledgements

The Ministero Istruzione dell'Università e della Ricerca (MIUR) and the University of Perugia are acknowledged for the financial support through the program "Dipartimenti di Eccellenza 2018-2022" (grant AMIS).

Notes and references

Conflicts of interest

There are no conflicts to declare.


Ayaki Sunaga

Trond Saue

Towards highly accurate calculations of parity violation in chiral molecules: relativistic coupled-cluster theory including QED-effects

Keywords: Relativistic quantum chemistry, parity violation, radiative corrections, molecular property, coupled cluster

Parity-violating energies EPV of the H 2 X 2 (X = O, S, Se, Te, Po) molecules are reported, calculated as analytical expectation values at the relativistic coupled-cluster singles-and-doubles (CCSD) level using property-optimized basis sets. Radiative corrections to the EPV was investigated using effective QED-potentials and found to reach a maximal value of 2.38% for H 2 Po 2 . However, this result depends on the choice of effective self-energy potential and may indicate limitations to their domain of validity.

Introduction

In early 2012 the journal Nature had a News Feature article on "tough science" listing five experiments "as hard as finding the Higgs" [1]. One experiment in the list attempts to provide the first observation of parity violation (PV), associated with the weak force, in molecular systems using high-resolution laser spectroscopy [2]. It is an example of a spectroscopic test of fundamental physics in the low-energy regime, using atoms and molecules, in vivid contrast to the Large Hadron Collider, 27 km long and operating at several TeV [3,4].

It was Lee and Yang who in 1956 pointed out that there was insufficient experimental evidence for the conservation of parity in processes involving the weak interaction and proposed experiments to test possible parity violation [5], soon to be confirmed by experiments [6][7][8][9]. The processes that were studied are mediated by charged vector bosons, in contrast to the neutral photon of electromagnetic interactions, and are therefore incompatible with stable atoms and molecules. The situation changed with the development of electroweak theory [10][11][12] and in particular the prediction of the existence of a neutral partner, Z 0 , of the charged W ± bosons. Early studies about the feasibility of observing parity violation in stable atoms were rather pessimistic [13,14],

CONTACT T. Saue. Email: trond.saue@irsamc.ups-tlse.fr but limited attention to the hydrogen atom. M. A. and C. Bouchiat noted that PV has a strong scaling with nuclear charge and therefore suggested to look for such effects in heavy atoms [15,16] and, indeed, the first observation of atomic PV was in the form of optical activity in a vapor of bismuth atoms [17,18]. Since then PV has been observed in several heavy atoms: thallium [19][20][21][22][23], cesium [24][25][26][27][28][29], lead [30][31][32], bismuth [33] and ytterbium [34,35]. For a recent review on atomic parity violation, see Ref. 36.

For chiral molecules parity violation induces a finite, albeit minute energy difference between enantiomers [37][38][39], making them formally diastereomers, an intriguing fact in view of the observed chiral preferences in biological systems [40][41][42], but yet to be observed. The experiment under construction in Paris will search for the signature of parity violation in the form of a frequency shift in the vibrational spectra of leftand right-handed molecules. It aims for a measurement precision below 0.1 Hz and is based on the method of the Doppler-free two-photon Ramsey fringes in a supersonic molecular beam [2,43]. Recent progress is reported in Ref. 44, whereas information about other experiments and molecular (as well as atomic) PV in general is found in an excellent review by Berger and Stohner [45].

The role of theory at this point is to guide experiment. We have developed a computational protocol for the rapid screening of candidate molecules, based on 2-component relativistic density functional theory (DFT) [46]. However, calculations of higher accuracy are needed, in the short run to calibrate the current computational protocol for molecular PV calculations, in the longer run to confront a possibly successful experiment. In this context coupled cluster theory stands out as the method of choice [47][48][49][50]. In 2000 Thyssen et al. reported 4-component wavefunction-based correlated calculations of the parity-violating energy of the H 2 O 2 and H 2 S 2 molecules using a finitefield approach at the level of 2nd-order Møller-Plesset perturbation theory (MP2), coupled-cluster singles-and-doubles (CCSD), CCSD with approximate triples correction (CCSD(T)) as well as configuration-interaction singles-and doubles (CISD) [51]. The H 2 X 2 -series of molecules are not good candidates for a PV experiment due to the almost free rotation around the X-X bond, but have been much used for benchmarking and analysis purposes due to their simple structures (see for instance Refs. [52][53][54][55][56][57][58][59][60]). In 2005 van Stralen et al. reported the first implementation of analytic one-electron properties at the 4-component relativistic MP2 level and reported PV energies for H 2 X 2 , (X=O,S,Se,Te) [61]. More recently, in 2016, these analytic calculations were extended to the CCSD level by Shee et al. [62].

In the present work we repeat the calculations of Shee et al. [62], but with a number of improvements: i) The coupled cluster calculations have been carried out using a new coupled cluster code, ExaCorr, implemented for massively parallel GPU-accelerated computing architectures [63]. ii) We have extended the Dyall-type all-electron basis sets with tight exponents optimized for PV calculations. iii) Using H 2 Se 2 , as detailed in the Supplemental Material, we have investigated the effect of truncation of the occupied and virtual spaces on the PV energy. iv) We have extended the calculations to include H 2 Po 2 . v) Most importantly, we have to our knowledge carried out the first study of the effect of quantum electrodynamics (QED) on molecular PV energies.

We are here interested in the electron self-energy (SE) and the vacuum polarization (VP), both at the origin of the Lamb shift [64]. QED, as expressed through the scattering matrix formalism, can describe atoms and molecules to amazing precision, but is in practice limited to few-electron systems, in particular due to the slow convergence of electron correlation within this formalism [65][66][67][68][69][70][71]. In recent years there has been an increased interest in the use of effective QED potentials [72][73][START_REF] Artemyev | Handbook of Relativistic Quantum Chemistry[END_REF][START_REF] Pašteka | [END_REF]. An early example of such a potential is the Uehling potential for vacuum polarization [76]. The electron self-energy is more complicated to represent by a potential due to its delocal nature, but a number of such potentials are now available [77][78][79]. QED-effects have also been included in the fitting of relativistic effective core potentials [80,81]. We have implemented effective QED potentials in the DIRAC code for general relativistic molecular calculations [82]. A full account of our implementation will be given elsewhere [83]; very recently, a similar implementation has been reported by Skripnikov [84]. QED-effects on parity violation have previously been studied in atoms and then with focus on transition amplitudes [85][86][87][88][89][90][91][92][93]. We believe that the present work is the first to study radiative corrections to molecular PV energies.

Our paper is organized as follows: In Section 2.1 we explain how effective QED potentials are included in our calculations. Section 2.2 gives a brief introduction to parity violation in chiral molecules. Computational details are given in in Section 3. In Section 4 we present and discuss our results, before concluding in Section 5. We employ SI-based atomic units throughout this paper.

Theory

Relativistic Hamiltonian with effective QED potentials

Our calculations are based on the Dirac-Coulomb Hamiltonian

ĤDC = V nn + i ĥ(i) + 1 2 i =j 1 r ij (1)
Here, the first and third terms represent the nucleus-nucleus and electron-electron repulsion, respectively. The one-electron part is given by the Dirac Hamiltonian ĥ = βm e c 2 + c (α

• p) + V en ; V en = -eφ n (2)
in the scalar potential φ n of fixed nuclei and will be extended by effective QEDpotentials V effQED ĥ → ĥ + V effQED ;

V effQED = A V VP A + V SE A (3)
where V VP A and V SE A refer to the contributions from the vacuum polarization and the electron self-energy, respectively, associated with nucleus A. We employ the Uehling potential [76] for the vacuum polarization term V VP . The Uehling potential for a spherically symmetric extended nucleus can be expressed as follows [94]

V VP A (r) = - 2αλ e 3 Z A e 2 r ∞ 0 dr r ρ A r K 0 2 λe r -r -K 0 2 λe r + r . (4)
where

K 0 (x) = ∞ 1 dt exp [-xt] 1 t 3 + 1 2t 5 t 2 -1 1/2 , (5)
α the fine structure constant, λe the reduced Compton wavelength and ρ A is the normalized charge density of nucleus A.

For the self-energy term V SE , we have employed two kinds of effective potentials. The first one is the Gaussian model potential proposed by Pyykkö and Zhao (PZ) [77],

V SE A → V PZ A (r) = B(Z A ) exp -β(Z A )r 2 , (6)
where the parameters B and β are second-order polynomials in nuclear charge with coefficients obtained by fitting to SE data (2s 1/2 SE energy shift of hydrogen-like ions [95,96] and relative SE shift of hyperfine transition energies of lithium-like ions [97]) in the range 29 ≤ Z ≤ 83.

The second effective SE potential was proposed by Flambaum and Ginges (FG) [78]. It is formulated in the language of form factors expressing the effective interaction of an electron with the electromagnetic field, radiative corrections included [START_REF] Schweber | An Introduction to Relativistic Quantum Field Theory[END_REF][START_REF] Berestetskii | Quantum Electrodynamics (Course of Theoretical Physics[END_REF]. The potential is written as a sum of three contributions

V SE A → V FG A (r) = V MAG A + V HF A + V LF A , (7)
where the first term is associated with the magnetic form factor

V MAG A (r) = αλ e 4π iβα • ∇ φ n;A (r) ∞ 1 dt 1 t 2 √ t 2 -1 e -2tr λe -1 (8)
and contains the gradient of the scalar potential φ n;A associated with nucleus A. The contribution from the electric form factor has been split into a high-frequency (HF) part

V HF A (r) = -A (Z A , r) α π φ n;A (r) ∞ 1 dt 1 √ t 2 -1 1 - 1 2t 2 × ln t 2 -1 + 4ln 1 Z A α + 0.5 - 3 2 + 1 t 2 e -2tr λe , (9)
obtained from theory after introducing a fixed low-frequency cutoff, and a lowfrequency (LF) part

V LF A (r) = -B (Z A ) Z 4 A α 5 m e c 2 e -ZAr . (10)
for which an approximate form has been chosen. Parameters A (Z, r) and B (Z) have been fitted to SE shifts of hydrogen-like ions [START_REF] Mohr | [END_REF]101]. These effective QED-potentials have been implemented by grafting routines from the numerical atomic code GRASP [102] onto the DFT-grid of the DIRAC code [82,83,103]. Since the QED effects are principally generated in the vicinity of nuclei, on the order of reduced Compton wavelengths [START_REF] Artemyev | Handbook of Relativistic Quantum Chemistry[END_REF], the numerical integration is only carried out on the nucleus associated with the QED-potential and with a radial cutoff. A problematic term in this respect is the low-frequency part of the Flambaum-Ginges potential since it by construction has the range of a hydrogen-like 1s orbital and for light elements may therefore overlap with other atoms in a molecule. Therefore, by default in the current implementation, QED-potentials are only activated for centers with Z > 18.

The potentials can also be used at the eXact 2-Component relativistic (X2C) level [104]. In the present work we employ the eXact 2-Component (X2C) molecularmean field (mmf) Hamiltonian [105]. In this approach an initial HF calculation is car-ried out at the 4-component relativistic level. The converged Fock matrix is then transformed to 2-component form and employed in the ensuing correlated calculations.

Parity-violating energy shift

Theoretical studies of parity violation in atoms and molecules are based on an effective Hamiltonian describing the weak interaction between electrons and nucleons mediated by the exchange of virtual Z 0 bosons. A detailed derivation and discussion of this Hamiltonian, starting at the level of quarks, is for instance found in Ref. [START_REF] Greiner | Gauge Theory of Weak Interactions[END_REF] (the reader may also consult Refs. 53,59,[START_REF] Berger | Relativistic Electronic Structure Theory. Part 2. Applications[END_REF][START_REF] Halzen | Quarks & Leptons[END_REF][START_REF] Hung | [END_REF][110]. The Z 0 boson is neutral, like its electroweak partner, the photon. On the other hand, it has mass -91.1876 (21) GeV/c [111], which is close to that of 98 Mo, the most abundant isotope of molybdenum [59]. This confers a short-range nature to the interaction. The effective Hamiltonian accordingly has the form of a contact interaction between two fermion 4-currents, ĥeff ∼ j µ;e (r, t)

Z i j p µ;i (r, t) + N i j n µ;i (r, t) d 3 r (11)
in line with the form suggested by Fermi in 1934 [112]. The fermion currents of electromagnetic interaction have the form ψγ µ ψ = ψ † (I 4 , α)ψ, which is termed a vector current since the space-like part of the current transforms as a vector under space inversion. The fermion currents associated with the weak interaction are of mixed vector and axial vector (V-A) nature [113][114][115], the latter on the form ψ † (γ 5 , Σ)ψ, where appears

γ 5 = 0 2 I 2 I 2 0 2 . ( 12
)
For the nucleon currents a non-relativistic approximation is used, setting the small components of the 4-spinors to zero, which in practice means that terms involving the odd operators γ 5 and α vanish whereas even operators are reduced to their 2component counterparts, that is, I 4 , Σ → I 2 , σ. Parity-violation arises from the V e -A p,n and A e -V p,n cross-terms. The former coupling leads to a nuclear spin-dependent effective Hamiltonian, which has so far been ignored in the calculation of parity-violation energies, but has been employed in studies of parity-violation in the NMR spectra of chiral molecules [116][117][118][119][120][121]. Here we limit attention to the latter coupling, which leads to a nuclear spin-independent effective Hamiltonian on the form

ĥPV eff = G F 2 √ 2 A Q W,A i γ 5 i ρ A (r i ). ( 13
)
where

Q W,A = -N A + Z A 1 -4sin 2 θ W , (14)
is the weak charge [16] of nucleus A given in terms of the numbers N A and Z A of neutrons and protons, respectively, as well as the weak mixing angle θ W [122,123]. The current recommended value1 from the Particle Data Group (PDG) [START_REF] Zyla | [END_REF] is sin2 θ W = 0.23122(4). The presence of the Fermi coupling constant G F = 2.22255 × 10 -14 E h a 3 0 in the effective Hamiltonian (13) shows that the interaction is truly weak. The operator can therefore not be added to the electronic Hamiltonian in an energy calculation, since its contribution would be masked by numerical noise, rather it will be obtained perturbatively as an expectation value. The expectation value with respect to the electronic wavefunction Ψ can in turn be expressed as a sum over atomic contributions,

E PV = A E PV A = G F 2 √ 2 A Q A W M PV A ; M PV A = i Ψ γ 5 i ρ A (r i ) Ψ . ( 15
)
The reduced PV matrix element M PV A is of highly local nature, not only because of the presence of the normalized nuclear charge density ρ A of each center A, but also because the γ 5 matrix (12) couples large and small components, the latter being of highly local atomic nature for positive-energy orbitals.

The parity-violating energy E PV is strictly zero in the absence of spin-orbit coupling and its meaningful calculation in a non-relativistic framework hence requires double perturbation theory [39,[52][53][54][127][128][129]. The non-relativistic form of the effective Hamiltonian ( 13) and spin-orbit interaction then contributes factors Z 3

A and Z 2 B , respectively, giving an overall Z 5 scaling. Single-center contributions (A = B) vanish in a non-relativistic framework [128], but not necessarily in a relativistic setting [59].

In a previous study, we analyzed the parity-violating energy of the H 2 X 2 -series of molecules at the Hartree-Fock(HF) and Kohn-Sham(KS) levels [59]. The E PV expectation value was decomposed in terms of pre-calculated atomic orbitals. It was found that a contribution E PV A is completely dominated by contributions from the same center A. Inserting the form of relativistic atomic orbitals

ψ(r) = R L (r)χ κ,mj (θ, φ) iR S (r)χ -κ,mj (θ, φ) , (16)
the expectation value E PV A can be expressed in terms of products of radial and angular contributions. From the angular part of it was found that only pairs of atomic orbitals with opposite values of κ and the same value of m j can contribute. From the small r behaviour of the radial functions it was found that for point charge nuclei only (s 1/2 , p 1/2 ) -pairs of atomic orbitals can contribute, whereas for extended nuclei contributions with |κ| > 1 are allowed. Finally, for the H 2 X 2 -series of molecules we found that the reduced PV matrix element M X PV can be remarkably well approximated by the mixing of valence ns 1/2 and np 1/2 orbitals on the same center X, that is

M X PV ≈ M X PV;nn × 2Re i c(ns X 1/2 ) * i c(np X 1/2 ) i , (17)
where

M X PV;mn = ms X 1/2 |γ 5 ρ X |np X 1/2 = R L;X ms1/2 |ρ X |R S;X np1/2 r + R S;X ms1/2 |ρ X |R L;X np1/2 r . (18)
In ( 17) c(ψ) i is the coefficient of atomic orbital ψ in the expansion of molecular orbital i. The quantity 2Re i c(ns X 1/2 ) * i c(np X 1/2 ) i thereby determines the mixing of atomic orbitals ns 1/2 and np 1/2 contributing to the PV energy, but it should be noted that this particular form refers to a specific phase convention for atomic orbitals. Further details are found in Ref. 59.

Computational details

Unless otherwise stated all calculations have been carried out using the DIRAC program [82] (revisions 9e12f6d89, 13d3f86a3, 9533dd59e, 09012db84, and ccb40d212). We employed a Gaussian model for the nuclear charge distribution in all calculations [130]. Molecular geometries were taken from Table I of Ref. 55 and were optimized at the MP2 level using quasirelativistic pseudopotentials for the heavy elements. In all calculations the dihedral angle was fixed at 45 • of the P enantiomer. For the weak mixing angle, we have used the value sin 2 θ W = 0.2319 in line with previous works [53,55] and corresponding to the PDG 1994 recommended value [131].

Various flavours of Dyall basis sets [132][133][134] have been employed for all elements of the target molecules. These basis sets have been extended by tight exponents for better description of the PV operator. These optimizations were carried out using the simplex method [START_REF] Press | Numerical recipes in Fortran[END_REF], minimizing the target function

∆M PV (%) = 100 × 1 N nm abs M num PV,nm -M bas PV,nm M num PV,nm , (19)
where M num PV,nm and M bas PV,nm refer to reduced PV matrix elements, (18), calculated numerically or in the current basis set, respectively. For this, we employed a development version of the atomic GRASP code [102] allowing the use of Gaussian basis sets, in this case the Dyall.nzp (n = 2, 3, 4) basis sets for SCF calculations.

For the MP2 and CCSD calculations a computational protocol was developed based on a calibration study of the parity-violating energy E P V of the H 2 Se 2 molecule, described in the Supplemental Material. All calculations employed the molecular-meanfield X2C Hamiltonian ( 2 DC M ) [105] based on the Dirac-Coulomb Hamiltonian. For technical reasons, MP2 and CCSD analytical expectation values were calculated using the RELCCSD [62,[START_REF] Visscher | [END_REF] and ExaCorr [63] modules, respectively. The expectation values were calculated using the Lagrangian-based analytical energy derivative technique described in Ref. 62 and references therein. For the CCSD calculations, we used the dyall.aae3z basis sets [132][133][134], extended with tight functions, as described above. We correlated (n -1)s(n -1)p(n -1)dnsnp electrons, and truncated the virtual space at 100 E h . For MP2 the calibration study indicated that only valence orbitals need to be correlated, so we used the the dyall.acv3z basis.

For comparison we carried out DFT calculations using the extended dyall.aae3z basis sets and functionals from different rungs of the Jacob's ladder of density functional approximations [START_REF] Perdew | Density Functional Theory and Its Applications to Materials[END_REF]: LDA (SVWN5) [START_REF] Slater | [END_REF]139], BLYP [140][141][142], and B3LYP [143][144][145], which is Becke's series, and PBE [146,147] and PBE0 [148,149], which is Perdew's series. We also used the long-range correlated functional CAMB3LYP [150].

Results and discussion.

For convenience, we have collected our HF, MP2 and CC results in Table 6 of the Supplemental Material together with previous results [51,61,62] calculated at the same geometries [55]. Our outputs are available via the zenodo repository [START_REF] Sunaga | Towards highly accurate calculations of parity violation in chiral molecules: relativistic coupled-cluster theory including QED-effects: Dataset[END_REF].

Optimization of tight exponents

In Table 1 we show the average relative error ∆M PV (19) of reduced PV matrix elements for various basis sets with respect to the reference numerical calculation. For Se, we see a reduction from 19.8% to 3.2% when going from dyall.2zp to dyall.4zp unmodified basis sets. At the dyall.3zp level the relative error is reduced for the heavier elements since these basis sets already contain quite tight exponents. Adding tight exponents, however, the error can be brought below 1% at all basis set levels. Specifically, we added 2p for O, 1s3p for S and 1p for Te at the dyall.3zp level. For the Se atom, we added 1s3p for dyall.2zp, 2p for dyall.3zp, and 1p for dyall.4zp, respectively. The values of the optimized exponents are given in the Supplemental Material. For Po no exponents were added since the relative error was already well below 1%.

Comparison of methods

Parity-violating energies of the H 2 X 2 molecules, calculated at the 45 • dihedral angle of the P enantiomer using different methods, are collected in Table 2. Using the CCSD results as reference, we find that the mean relative error (MRE) of the HF method, corresponding to the relative correlation effect, is -2.8%, so rather modest, in line with previous studies [51,62]. On the other hand, the standard deviation ∆ std approaches ten percent, showing quite a spread. MP2 has a slightly smaller MRE of opposite sign (2.0 %) and similar spread. No density functional outperforms HF; their average MRE and ∆ std are -11.1 % and 7.0%, respectively. The best performing functional is PBE0 (-5.7± 5.0 %), the worst CAMB3LYP (-17.2±4.3 %). These observations are as such disappointing since, as already mentioned, the computational protocol developed for the scan of candidate molecules for the Paris experiment is based on DFT. On the other hand, the protocol calculates a vibrational shift, corresponding to the PV energy difference of two vibrational levels, rather than an absolute electronic PV energy as here. Further study is therefore required to judge the performance of the computational protocol.

Radiative corrections to PV energies

In Table 3 we show QED-effects on parity-violating energies for the heavier members of the H 2 X 2 -series, calculated at the CCSD level. We used the Uehling term (U) for the inclusion of vacuum polarization (VP) and the Pyykkö-Zhao (PZ) or Flambaum-Ginges (FG) potentials for electron self-energy (SE). Fitting our data to a simple power law ∼ Z γ X with respect to the nuclear charge of the heavy atom X, we find that the uncorrected PV energy has exponent 7.0, whereas the radiative corrections show stronger scaling, with γ equal to 8.1 and 7.7, for ∆(U+PZ), and ∆(U+FG), respectively. In passing we note that this trend may not hold for superheavy systems, since SE and VP terms tend to have opposite sign, as known for valence atomic properties (see for instance Refs. 73,[START_REF] Pyykkö | [END_REF].

From Table 3 we see that the radiative correction may reach a few percent of the PV energy, but only if we look at the results obtained with the FG potential. The discrepancy between the results obtained with the PZ and FG potentials could signal that one or both are being used beyond their domain of validity, since they have been principally designed to reproduce QED energy shifts.

Conclusions

In this article we have given improved values for the parity violating energy of H 2 X 2 molecules (X = O, S, Se, Te, Po). Improvements include: i) Property-optimized basis sets, which go beyond conventional energy-optimized basis sets. ii) We have used a new relativistic coupled cluster code, ExaCorr, geared towards massively parallel calculations [63], allowing extensive calibration studies in order to elaborate a CCSD computational protocol for PV energies and to include H 2 Po 2 among the target molecules iii) Inclusion of QED-effects through effective potentials.

Using our CCSD numbers as reference, we find the performance of DFT somewhat disappointing. This is possibly worrisome, since DFT is the basis of the computational protocol elaborated for the scan of candidate molecules for the high-resolution laser spectroscopy experiment which is under development in Paris and which aims to pro-vide the first observation of parity violation in molecules as a shift in the vibrational spectra of enantiomers of chiral molecules. We shall pursue our investigation of the reliability of this protocol by calculating at the coupled cluster level the vibrational shift, involving energy differences, of candidate molecules. On the positive side, it has been noted that HF, here better performing, tends to predict larger PV vibrational shifts than DFT [46].

A major novelty of this work is the first calculation of molecular PV energies E PV including QED effects. For the heaviest molecules in the H 2 X 2 -series the radiative corrections are on the order of a few percent. However, we do note that this depends on what effective self-energy potential is used. This may indicate limitations to their applicability, an issue that encourages us to develop more rigid approaches to incorporate QED effects.

Further improvements would include going beyond the zeroth-order term, that is, the instantaneous Coulomb term, of the fully relativistic electron-electron interaction [153]. The contribution from the Breit term, or more precisely spin-otherorbit interaction, has been estimated by Berger within the Breit-Pauli framework [58]. The Breit effects for H 2 X 2 (X = O, S, Se, Te, Po) are 12%, 4.7%, 1.9%, 1.2%, 0.7%, respectively and hence monotonically decreasing along the series, contrary to the radiative corrections.

Another improvement would be to increase the excitation level beyond that of CCSD. An indication of the role of higher-order correlation effects comes from the finite-field study of Thyssen et al [51]: The relative difference between CCSD(T) and CCSD was found to be -1.0% and -2.5% for H 2 O 2 and H 2 S 2 , respectively. This suggests that these effects may be at the percentage level and increase down the series. This clearly calls for further study.

Table 1 .

 1 Average relative error ∆M PV (%) in reduced PV matrix elements in different basis sets (with and without added tight exponents) with respect to the reference numerical values.

		O	S		Se		Te	Po
	basis	dyall.3zp dyall.3zp dyall.2zp dyall.3zp dyall.4zp dyall.3zp dyall.3zp
	without	11.4	11.4	19.8	8.5	3.2	3.4	0.1
	with	0.9	0.5	0.9	0.6	0.2	0.1

Table 2 .

 2 Summary of the values of E PV (E h ), calculated at the 45 • dihedral angle of the P enatniomer of H 2 X 2 molecules at various levels of theory. "err" is the percentage error with respect to CCSD. H 2 O 2 H 2 S 2 H 2 Se 2 H 2 Te 2 H 2 Po 2

	Method	E PV	err (%)	E PV	err (%)	E PV	err (%)	E PV	err (%)	E PV	err (%)
	CCSD	-5.917E-19		-2.166E-17		-2.463E-15		-3.536E-14		-1.337E-12
	MP2	-6.525E-19	-10.3 -2.067E-17	4.6 -2.191E-15	11.1 -3.236E-14	8.5 -1.390E-12	-3.9
	HF	-6.808E-19	-15.0 -2.160E-17	0.3 -2.301E-15	6.6 -3.403E-14	3.8 -1.467E-12	-9.7
	LDA	-6.581E-19	-11.2 -2.659E-17	-22.7 -2.910E-15	-18.2 -4.164E-14	-17.7 -1.308E-12	2.2
	BLYP	-6.303E-19	-6.5 -2.663E-17	-22.9 -2.864E-15	-16.3 -4.095E-14	-15.8 -1.324E-12	1.0
	B3LYP	-6.515E-19	-10.1 -2.576E-17	-18.9 -2.796E-15	-13.5 -4.032E-14	-14.0 -1.372E-12	-2.6
	CAMB3LYP -6.627E-19	-12.0 -2.621E-17	-21.0 -2.941E-15	-19.4 -4.263E-14	-20.5 -1.515E-12	-13.3
	PBE	-6.269E-19	-5.9 -2.498E-17	-15.3 -2.653E-15	-7.7 -3.814E-14	-7.8 -1.261E-12	5.7
	PBE0	-6.487E-19	-9.6 -2.404E-17	-11.0 -2.570E-15	-4.4 -3.729E-14	-5.4 -1.313E-12	1.8

Table 3 .

 3 QED effects on the parity-violating energy E PV (E h ), calculated at the 45 • dihedral angle of the P enantiomer of H 2 X 2 molecules. The radiative corrections have been included using the Uehling potential (U) for vacuum polarization and either the Pyykkö-Zhao (PZ) or Flambaum-Ginges (FG) effective potentials for electron self-energy. Absolute and relative differences are also shown.

		non-QED	U+PZ	U+FG	∆(U+PZ)(E h ) ∆(U+FG)(E h ) ∆(U+PZ)(%) ∆(U+FG)(%)
	H 2 Se 2	-2.463E-15 -2.458E-15 -2.434E-15	4.6E-18	2.9E-17	-0.19	-1.18
	H 2 Te 2	-3.536E-14 -3.525E-14 -3.471E-14	1.2E-16	6.6E-16	-0.33	-1.89
	H 2 Po 2	-1.337E-12 -1.331E-12 -1.306E-12	6.8E-15	3.1E-14	-0.51	-2.38

The value depends on the renormalization scheme employed. The PDG employs the modified minimal subtraction scheme, whereas the conceptually simpler, but as such less accurate[124] CODATA value

sin 2 θ W = 1 -(m W /m Z ) 2 = 0.22290(30) [START_REF] Tiesinga | The 2018 CO-DATA Recommended Values of the Fundamental Physical Constants[END_REF], given explicitly in terms of the masses of the Z 0 and W ± bosons is obtained from the on-shell scheme.

Acknowledgement(s)

This paper is dedicated to John F. Stanton, a driving force of coupled cluster theory for 35 years ! We would like to thank David Sanchez(Toulouse) and Johann Pototschnig(Vienna) for a lot of technical support. This work is part of the project molQED funded by the Agence Nationale de la Recherche (ANR, France). This research used resources of the Oak Ridge Leadership Computing Facility, which is a DOE Office of Science User Facility supported under Contract DE-AC05-00OR22725. We are thankful to CALMIP (Calcul en Midi-Pyrénées) and LCPQ (Laboratoire de Chimie et Physique Quantiques) clusters for the additional computation time. We are also thankful to the supercomputer of ACCMS (Kyoto University) and Research Institute for Information Technology, Kyushu University (General Projects) for test calculations. A.S. acknowledges financial support from the Japan Society for the Promotion of Science (JSPS) KAKENHI Grant No. 17J02767 and 20K22553, and JSPS Overseas Challenge Program for Young Researchers, Grant No. 201880193.