Complex Energies and Transition-Dipoles for the Uracil anion Shape-type Resonances from stabilization curves via Padé

Absorption of slow moving electrons by neutral ground state nucleobases have been known to produce resonance, metastable, states. There are indications that such metastable states may play a key-role in DNA/RNA damage. Therefore, herein, we present an ab-initio, non-Hermitian investigation of the resonance positions and decay rates of the low lying shape-type states of the uracil anion. In addition, we calculate the complex transition dipoles between these resonance states. We employ the resonance via Padé (RVP) method to calculate these complex properties from real stabilization curves by analytical dilation into the complex plane. This method has already been successfully applied to many small molecular systems and herein we present the first application of RVP to a medium-size system. The presented resonance energies are converged with respect to the size of the basis set and compared with previous theoretical works and experimental findings. Complex transition dipoles between the shape-type resonances are computed using the energy-converged basis set. The ability to calculate ab-initio energies and lifetimes of biologically relevant systems opens the door for studying reactions of such systems in which autoionization takes place. While the ability to also calculate their complex transition dipoles open the door for studying photo induced dynamics of such biological molecules.


Introduction
High-energy particles may lead to DNA damage, such process can be genotoxic causing mutagenic, recombinogenic as well as DNA single-and double-strand breaks. 1 Apparently, those processes are not produced by the primary high-energy quanta, but by secondary species such as free electrons with low energy. 2 Throughout the radiation, low energy electrons (LEE) are produced in very large quantities with an energy distribution lying below 10 eV; these electrons are more likely to induce significant amount of chemical damage within the cell. 3-6 LEE attached to molecules often lead to metastable anion state called resonances.
The resonance phenomenon is a state of the system that as time passes breaks into several subsystem. That is, even though the system has enough energy to break apart, it does not happen instantly but require some finite time (relate to the width/lifetime/decay of the resonance). In non-Hermitian quantum mechanics resonances are associated with complex eigenvalues: E = E r − i 2 Γ, where the real part relates to the energy position of the resonance, while the imaginary part relates to the decay rate of the metastable state. 7 Resonance states in biological molecules, like DNA and RNA nucleobases, may lead to breaking chemical bounds, which is often called dissociative electron attachment (DEA).
The outcome of DEA process highly depends on the energy and lifetime of the attached electron, i.e., it determines fragmentation or non-dissociative relaxation mechanism. [8][9][10] To get a better understanding of the mechanisms responsible for DNA damage by LEE via DEA requires knowledge of the interaction of LEE with the basic compounds, such as a nucleobase. 11 Studying the damage mechanism has been of great interest in the experimental and theoretical communities, however, showing various processes with very different energy values. [12][13][14][15][16][17][18] Therefore, it is desirable to have methods that are able to treat such systems in an ab-initio fashion. Moreover, in order to be able to study the interaction of resonances with light it is necessary to be able to calculate also other complex properties such as transition dipoles.
Recently we introduced and benchmarked an approach, known as, resonance via Padé (RVP), to calculate atomic and molecular resonances from standard (Hermitian) electronic structure packages combined with analytic dilation, via Padé, into the (non-Hermitian) complex plane. 19,20 This approach is based on the stabilization technique. [21][22][23] 30,31 Therefore, it is desirable to examine the performance of RVP in treating also larger chemical system. In this work we present the first application of RVP to a medium-size system, the uracil anion (with 59 electrons), while focusing on its shape-type resonances.
Uracil is the smallest nucleobase in RNA and it resembles the DNA base thymine, making it appealing for theoretical calculations. The resonance states are formed by attachment of an electron to one of the unoccupied virtual (π * and σ * ) orbitals of the neutral ground state.
The π * shape-type resonance have been observed experimentally and calculated theoretically, while the σ * have been only treated theoretically . 32 own advantages and disadvantages. 40 Nevertheless, RVP possess several clear advantages: as a method that relies on analytical continuation it is straightforward to test the analyticity of the employed data (by reconstructing the entire stabilization graph using only small part of the data). Further, the RVP fitting is less sensitive than in GPA to the specific chosen data set and/or to adding extra data points. 40 In addition, RVP allows calculations of other complex properties, such as the complex transition dipoles, 41 which are essential for studying light-matter interactions. Herein we also report the transition dipoles between the resonance states, these can be used, for example, to calculate the photoionization spectra and the Fano asymmetry parameter. 42 Thus, the information presented herein enables future investigation of uracil and the mechanisms involving DEA. Below, we shortly discuss a mechanism that may minimize radiation damage to biological systems, such as uracil, and illustrates the need for calculating complex energies as well as transition dipoles.

Computational Details
The neutral geometry of uracil is optimized at the MP2/cc-pVTZ level of theory with a Cs symmetry (see Supporting Information for more details). 43,44 The energy positions and decay rates of the three lowest π * shape-type resonances of the uracil anion are calculated using the RVP technique. Equation-of-motion coupled-cluster with singles and doubles for electron affinities (EOM-EA-CCSD) 45 is used to calculate the stabilization graphs, where the singlet ground state of the neutral uracil serves as the reference state. All calculations were performed with the quantum chemistry package Q-Chem. 46 We perform a basis set convergence employing the following basis sets: Pople's 6-31+G, 6-311+G, and 6-311+G(2d,p), as well as Dunning's aug-cc-pVDZ, aug-cc-pVDZ+1s1p1d (in which we add one diffuse function of s, p and d angular momentum on each atom except for the hydrogen atoms) and aug-cc-pVTZ. Moreover, upon concluding that a triple-ζ basis set is essential, we investigate the effect of adding diffuse functions systematically on top of the cc-pVTZ basis set (while using aug-cc-pVTZ for the hydrogen atoms). We report six such cc-pVTZ sets, which are augmented with +1s1p, +2s2p, +3s3p, +2s2p1d, +3s2p1d and +2s2p2d. These additional diffuse functions were added in an even-tempered fashion (with respect to the value of the original most diffuse function with an even-temper parameter of 2). All the employed basis sets are presented in the Supporting Information (SI). Within RVP a Padé function is dilated into the complex plane based on data obtained from a real stabilization calculation. The data is taken from the stable (analytic) region of a branch, however several such branches may exist. The reported complex energies corresponds to the statistically best-behaved results, see Method Section for additional details. Fig. 1 presents such a stabilization graph (additional graphs, using different basis sets, are given in the SI). Fig. 1a . 2a) than for the energy calculations (Fig. 1c))

Results and Discussions
Complex Energies -Convergence with respect to the size of the one-electron basis set Herein, we present the computed complex energies (positions and widths) of the first three low-lying shape type resonance states of the uracil anion. We check for convergence of the complex energies with respect to the size of the basis set, where spacial attention is given to the role of the diffuse functions. The results in Table 1 are sub-divided into four groups or panels. In the first panel, we use different Pople's basis sets. A comparison between the first three bases clearly highlight the importance of using a triple-ζ (TZ) basis set. The resonance energy position for all the three shape-type resonances decrease on going from the 6-31+G to the 6-311+G basis sets, whereas the effect of additional polarization functions [6-311+G(2d,p)] is even more pronounced.
The importance of polarization functions in the Pople basis sets, naturally led towards the highly polarised Dunning basis sets. On the second panel, we now compare Dunning's correlation-consistent augmented-DZ and TZ basis sets. We observe that even at the double ζ level (aug-cc-pVDZ) we obtain very similar results to the largest Pople basis set. Of note, the decay rate for the 3π * resonance state had a more abrupt increase to ∼0.4 eV as compared to ∼0.2 eV in the Pople basis sets. Additional diffuse basis functions play a key role in the evaluation of the widths (-2ImE) as can be seen from comparing the aug-cc-pVDZ results with aug-cc-pVDZ+1s1p1d. On going from aug-cc-pVDZ to aug-cc-pVTZ the energy positions (ReE = E r ) are also clearly affected. Therefore, in the third and forth panels we employ the cc-pVTZ basis set and systematically augment it with increasing number of diffuse functions.
Notice that the energy position differences, in the 3 rd and 4 th panels, are within 0.1 eV, thus the energy position is converged with respect to the additional diffuse basis functions.
As for the width, we see in the 3 rd panel that adding only +1s1p basis functions is not sufficient and that adding +2s2p is essential. Although adding +3s3p does make an effect (on the 2π * width), the effect is small. Therefore, in the 4 th panel we examine the effect of adding d-functions on top of the +2s2p diffuse functions. Clearly, adding only one d-type function (+2s2p vs. +2s2p1d) as well as augmentation via +3s2p1d have a very small effect, while adding two d-diffuse functions (+2s2p vs. +2s2p2d) does have some effect. Therefore, we conclude that the largest cc-pVTZ+2s2p2d is the optimal basis set. Table 2 Figure 2: (a) Energy stabilization plot for uracil anion, at the EOM-EA-CCSD/cc-pVTZ+2s2p2d (H: aug-cc-pVTZ) level. The black squares represent the sets of data points that are common between the two stabilised areas (the overlap). (b) Transition dipole stabilisation plot for 2π * ↔ 3π * . The black data points corresponds to a stable part on the plot, which has the same α-range as the overlap area in the energy graph. These points are taken as inputs for the analytical dilation in the Padé approximate.

Complex Transition Dipoles between the Uracil Anion Resonances
Using the converged basis set, cc-pVTZ+2s2p2d, we also calculated the complex transition dipoles between the three resonance states of the uracil anion. Employing the RVP technique to compute complex transition dipoles have been successfully benchmarked in Ref 41. Herein, we present the first application of RVP to property calculation of a medium-size molecular system of biological interest. The complex transition dipoles are an essential in order to study light-matter interactions within the non-Hermitian formalism. Figure 2 illustrates the RVP procedure for calculating complex transition dipoles for the 2π * ↔ 3π * case, i.e., between the second and third shape resonance states. Figure 2a depicts the energy stabilisation plot for these shape-type states. The overlap between the two stable regions is highlighted in black. The overlap represents an area, in the parameters space, in which the transition dipole is also stable, as expected form any physical property.
That is, the overlap region in Figure 2a is the same as the stable regions in Figure 2b.
The complex transition dipoles are calculated in a similar manner to the procedure for the complex energies. The data points marked in black in Figure 2b serve as input points for the Padé fitting. This is followed by analytically dilation of the Padé function into the complex plane and search for clusters of stationary points, which corresponds to the complex dipoles. The outcome of this procedure, i.e., the complex transition dipoles between the three resonance states, are given in Table 3. The real part dominant the three transition dipoles, where the imaginary part corresponds to about 1% of it or less. The transition dipoles are needed for describing the interaction of the uracil anion with light. Table 3: Complex transition dipole moments, real and imaginary parts (in a.u.) between the three lowest shape-type resonances of uracil anion obtained in this work using RVP. Basis set: cc-pVTZ+2s2p2d Reµ Imµ Reµ Imµ Reµ Imµ 1π * ↔ 2π * 1π * ↔ 3π * 2π * ↔ 3π * 5.089e-01 -3.599e-03 8.782e-01 -6.017e-03 8.204e-01 -1.628e-02 A simple example that illustrates the need for complex transition dipoles is to be able to transfer the uracil anion from one resonance state to another in order to minimize DNA damage via dissociative electron attachment. According to Matsika and coworkers the 3π * state, unlike 1π * and 2π * , is reactive with respect to CO elimination. 10 Therefore it is desirable to eliminate the attachment of an electron to the 3π * state, i.e., we wish to transfer the electron from 3π * to the 1π * or 2π * states. However, simply using a laser to induce a transition from one resonance to another dose not solve the problem. Since a typical laser pulse is of low intensity and their envelop supports many optical cycles, the oscillations between the two resonances are unavoidable. Nevertheless, by using chirped laser pulses an asymmetric switch can be obtained, which enables the transfer of 3π * to 2π * or 1π * without the backward transfer to 3π * . The asymmetric switch is obtained as the laser encircles a special point in the frequency-intensity 2D space. The special point, know as an exceptional point, is a point in this 2D space in which two resonances coalesce (i.e non-Hermitian degeneracy). 7 Explaining this mechanism is out of the scope of this work, for a general theoretical representation see Ref 51 and for an experimental demonstration of such an asymmetric switch see Ref 52. Notwithstanding, we perceive that in order to design an asymmetric switch one requires the complex energies (Table 2) and transition dipoles (Table 3).

Summary
We study the three lowest shape-type resonances of the uracil anion using the resonance via Padé (RVP) approach. RVP was already successfully employed in studying atomic and small molecular resonances, and herein we apply this approach to a medium-size chemical system with a biological interest. The presented results are converged with respect to the oneelectron basis set and compared to other theoretical and experimental values. In addition, we present the calculated complex transition dipoles between the first resonance states and the second and third ones. These properties are essential in studying light-matter interaction within the uracil anion, which may lead to progress in understanding the mechanism and even minimize DNA/RNA damage.

Methods
The resonances energies are calculated with the RVP method. RVP employs Hermitian electronic structure calculations and move into the Non-Hermitian regime via the Padé approximant. The method is based on analytic continuation of stabilization graphs. We start by generating the stabilization graphs, which we computed on the real axis using standard Hermitian electronic structure formalism. [21][22][23] To generate such graphs, we scale a finite Gaussian basis set by a real factor (α) , i.e., we divide the exponent of the most diffuse basis functions with α. For α < 1 the space spanned by the basis set compress and for α > 1 it expands. The plot of the eigenvalues (energies) as function of the real factor is known as a stabilization graph. Continuum, resonance and bound states behave differently upon scaling due to the different nature of their associated wave-functions. The continuum states are associated with a delocalized wave-function, whereas the resonance and bound states are more localized in the interaction region. 7 Consequently, the resonance and bound eigenvalues would not depend strongly on the scaling parameter. Contrary, continuum eigenvalues will strongly depend on the scaling factor. For the resonance eigenvalues, unlike for the bound ones, we expect crossing attempts by continuum eigenvalues, also known as, avoided crossings.
To calculate the resonances positions as well as its decay rate, we are using analytical continuation for the real axis (stabilization graphs) to the complex plane. We used the Padé approximant, where an energy function, E(η) (where η = αe iθ with θ = 0) is fitted to a ratio between two polynomials P (η) and Q(η) E(η) = P (η) Q(η) . (1) In practice, we use the Schlessinger point method 53 to generate a numerical expression to E(η). The data for the fitting were only taken from the stable part of the stabilization graphs, which is know to be an analytic region. 19 α and θ are the real stretching and rotating (into the complex plane) scaling parameters. Next, we substitute a complex parameter η = αe iθ into the fitted energy function and the resonances are identified as stationary points (SPs) in the complex α and θ plane, i.e., they satisfy: Graphically, these stationary points are associated with cusps in the α-and θ-trajectories. A θ-trajectory is generated by fixing α and varying θ over a range of values, the opposite is true for the α trajectory. Convergence is achieved when those two trajectories from cusps that meet. 54 Finally, we checked the stability of the complex resonance energies with respect to small modifications of the Padé input data set using a statistical approach. 55 We choose many input data sets for the Padé fitting and we are looking for clusters in the complex energy space, these clusters are associated with resonance plotted as a function the Padé input. That is, the average of the clusters with the smallest deviation is the reported result.
Notice that those chosen data sets successfully reproduced the original stabilization curve; we consider this as a numerical proof for the analiticity of the selected stable region. 19,56 Data Availability Statement The data that support the findings of this study are available from the corresponding author upon reasonable request.