Oscillator strengths of electronic excitations with response theory using phase including natural orbital functionals

The key characteristics of electronic excitations of many-electron systems, the excitation energies ω α and the oscillator strengths f α , can be obtained from linear response theory. In one-electron models and within the adiabatic approximation, the zeros of the inverse response matrix, which occur at the excitation energies, can be obtained from a simple diagonalization. Particular cases are the eigenvalue equations of time-dependent density functional theory (TDDFT), time-dependent density matrix functional theory, and the recently developed phase-including natural orbital (PINO) functional theory. In this paper, an expression for the oscillator strengths f α of the electronic excitations is derived within adiabatic response PINO theory. The f α are expressed through the eigenvectors of the PINO inverse response matrix and the dipole integrals. They are calculated with the phase-including natural orbital functional for two-electron systems adapted from the work of L¨owdin and Shull on two-electron systems (the phase-including Löwdin-Shull functional). The PINO calculations reproduce the reference f α values for all considered excitations and bond distances R of the prototype molecules H 2 and HeH + very well (perfectly, if the correct choice of the phases in the functional is made


I. INTRODUCTION
The excitation energies ω α and the oscillator strengths f α are the key characteristics of electronic excitations in manyelectron systems. The latter quantity provides a measure of the intensity of the corresponding excitation in the electronic spectrum. Linear response theory applied to an independent electron model system provides an efficient framework for the evaluation of these quantities. The response of a one-electron property, such as the electron density or the one-particle density matrix, δ g(ω), to an external perturbation δv ext (r, ω) with frequency ω is evaluated as follows: Here, the vector V collects the matrix elements of the perturbation in the chosen one-electron basis, V pq = φ p |δv ext |φ q . χ −1 (ω) is the inverse response function corresponding to the chosen perturbation and the electronic property. It can be expressed in the case of real perturbations and real property functions in terms of a matrix [g 0 ](ω) with the latter being considered in response theory as a functional of the ground-state function g 0 . With (1) and (2), the excitation energies ω α are obtained as the frequencies at which a response exists even in the absence of an external per-turbing field (ω's where χ (ω) diverges, which are the zeros of χ −1 (ω)). These frequencies are the roots ω α of the following matrix equations: Virtually, all applications of the theory are made within the adiabatic approximation, in which the frequency-dependent (ω) is replaced with a frequency-independent static matrix, (ω) ≈ 0 , so that all ω α and F α of the adiabatic eigenvalue problem can be simultaneously obtained with a single diagonalization of 0 . The eigenvectors F α of Eq. (4) are used to evaluate the oscillator strengths f α (See Sec. II where the derivation will be given).
In particular, the adiabatic response approach of timedependent density functional theory (TDDFT) 1 has become the method of choice for the calculation of electronic excitations in large molecules. In this approach, the electron density ρ(r, ω) is chosen as the one-electron function g(ω) in Eq. (1). The squared differences ( a − i ) 2 between the energies of the virtual ψ a (r) and occupied ψ i (r) Kohn-Sham (KS) orbitals of DFT, constitute the leading diagonal contribution to 0 of Eqs. (2)- (4). Defining the diagonal matrix E ia,j b = δ ij δ ab ( a − i ), the matrix 0 can be written as 0 = E 2 In order to produce a better approximation to the true excitations ω α , the orbital energy differences are shifted and coupled with the coupling matrix K constructed with the Hartree-exchange-correlation (Hxc) kernel f Hxc (r 1 , r 2 ).
In spite of its general success, adiabatic response TDDFT also exhibits serious failures. Using the response of the KS orbitals to obtain the density response, it is found that δρ can be expanded only with the virtual-occupied pairs of the ground-state KS orbitals as i∈occ a∈virt δρ ia ψ i ψ * a . Because of the simple structure of the density response δρ(ω), which includes only "transitions" ψ i → ψ a from occupied to virtual KS orbitals, adiabatic TDDFT lacks double excitations. 2,3 As was shown in Refs. 2 and 4-8, double excitations could be, in principle, recovered with a frequency dependent Hxc kernel f Hxc (r 1 , r 2 , ω) with a rather involved frequency dependence. Furthermore, because of the features of the KS orbital spectrum, conventional adiabatic TDDFT yields vanishing excitation energies in the case of bond breaking [9][10][11][12] and too low excitation energies in the case of long-range charge transfer (CT). 13,14 As was established in Refs. 15 and 16, in order to provide the correct excitation energies, the Hxc kernel should diverge, which classifies the bond breaking and CT cases as problematic for adiabatic TDDFT.
These problems of adiabatic TDDFT are illustrated with Figures 1-4, which show the excitation energies ω TDDFT and oscillator strength f TDDFT for the lowest two 1 + u excitations in the H 2 molecule, and the lowest two 1 + excitations in the HeH + molecule. The comparison is between standard TDDFT calculations (BP86, basis set aug-cc-pVTZ) with the reference quantities ω FCI and f FCI of the full configuration interaction (FCI) method. The failure of TDDFT for these excitations (and not only for these) have been spelled out in Refs. 11, 12, and 17. We note that the (larger) basis used here shows considerable difference with the previous (smaller) one, notably for the 2 1  the TDDFT curve. This is a consequence of the use of an augcc-pVTZ basis rather than the nonaugmented cc-pVTZ one of Ref. 12. Further extension to aug-cc-pVQZ has little effect. It is surprising that, whereas the TDDFT excitation energy is totally wrong for the first excited state of 1 + u symmetry of H 2 , the oscillator strength seems to be still qualitatively right. As we will demonstrate below, this is fortuitous, being caused by a cancellation of errors.
The problems of TDDFT were addressed with the development of time-dependent density matrix functional theory (TDDMFT). 12,17,18 In this theory, the first-order reduced density matrix (1RDM) γ (r, r , ω) is chosen as the FIG one-electron function g(ω) of Eq. (1). The 1RDM response in Eq. (1) is represented as the vector δγ (ω) = (δγ R (ω), δn(ω)) T , where δγ R (ω) collects the responses of the real part of the off-diagonal matrix elements δγ R pq (ω), while δn(ω) is the response of the diagonal elements n p = γ pp . The matrix elements of the perturbed 1RDM are expressed in the natural orbital (NO) basis. The natural orbitals are defined as the eigenfunctions of the 1RDM and the corresponding eigenvalues are called the occupation numbers.
The adiabatic approximation (4) is also applied in TDDMFT, 19 but then it encounters a serious problem with the proper evaluation of the response of the occupation numbers δn(ω). 17,[19][20][21][22] It can be proven 17,21,22 that an adiabatic TD-DMFT with a functional 0 [γ 0 ] produces an incorrect zero response of the NO occupations, δn(ω) = 0. This makes it impossible to describe diagonal double excitations of the type (φ i ) 2 → (φ a ) 2 , which entail a response in the occupation numbers (diagonal elements of the 1RDM). This deficiency also holds for the conventionally used functionals that depend on the NOs and NO occupations, 0 [{φ 0 i }, {n 0 i }], when they do not depend, just as the 1RDM itself, on the phases of the NOs φ 0 i . Furthermore, it has been shown that the adiabatic TD-DMFT results are quite poor if no phase dependence, as introduced in the PINO theory, is allowed in the functional. One can write the exact ground state energy of the two-electron system with exchange integrals ij|ji , so that the energy is not dependent on the phases of the NOs. This defines the socalled Löwdin-Shull density-matrix functional (DMLS). 21 In spite of the fact that the DMLS functional is exact for the ground state energy, and is a true density matrix functional, the adiabatic TDDMFT calculations with this functional, for e.g., the 1 + u H 2 excitation energies yield a completely spurious excitation spectrum, see Figure 5. 23 This problem of adiabatic TDDMFT was addressed recently with the development of phase-including NO (PINO) functional theory, 21, 24, 25 which goes beyond the TDDMFT proper. In the response PINO theory, the same response δγ (ω) = (δγ R (ω), δn(ω)) T is considered in (1) as in TD-DMFT. However, more variables are inserted into the representation of (ω) in (2) and (3). These additional variables assume the form of the phases of the NOs, so that 0 is considered in adiabatic PINO theory as a functional 0 [{ π 0 i }, {n 0 i }] of orbitals and their phases, i.e., the phase-including NOs or PINOs π 0 i and their occupations n 0 i . A proof of principle has been given for prototype two-electron systems that adiabatic response PINO theory is able to resolve the above mentioned problematic cases of adiabatic TDDFT in calculation of excitation energies ω α . 26 In this paper, the expression for the oscillator strengths f α of electronic excitations is derived within the adiabatic response PINO theory and f α is calculated for the prototype two-electron molecules HeH + and H 2 . In Sec. II, a full account of the response equations (1) and (2), as well as the eigenvalue problem (4) of adiabatic PINO theory is presented. The adiabatic PINO matrix 0 [{ π 0 i }, {n 0 i }] is obtained by taking (functional) derivatives of the ground state functional for the electron-electron interaction energy with respect to the PINOs π 0 i and their occupations n 0 i . The dynamic dipole polarizability tensor α(ω) is expressed through the response function χ(ω) of (2). The eigendecomposition of χ (ω) is employed to represent α(ω) with the sum-over-states formula. From the residues of this formula, an expression is derived, which relates the oscillator strengths f α directly to the eigenvectors of 0 . In Sec. III, the results of the PINO calculations of f α are presented for the lowest excitations along the dissociating bond coordinate R of the prototype two-electron molecules HeH + and H 2 . The response matrices are obtained from the electronelectron interaction energy in terms of NOs, occupation numbers and orbital phases, which are available from the work of Löwdin and Shull. 27 This is an explicit PINO functional (the phase-including Löwdin-Shull functional, PILS) for the electron-electron interaction energy. The PINO calculations reproduce the reference f α values for all considered excitations in HeH + and H 2 if the phase information in PILS is completely correct. We comment on deviations for the case there are errors in the phase choices in Sec. III. Importantly, even with a severely reduced size of the response matrix 0 (ω), the oscillator strengths f α remain very accurate, as we have found earlier for the excitation energies. 26 The shape of the f α (R) curves is rationalized with a decomposition analysis of the corresponding transition densities and dipole integrals in the basis of the KS molecular orbitals (MOs). Such an analysis also explains why the TDDFT f (1 1 + u )(R) does not seem to fail so spectacularly as the excitation energy does. In Sec. IV, the conclusions are drawn.

II. PINO OSCILLATOR STRENGTHS
In adiabatic PINO theory, the matrix response Eqs. (1) can be written in the ground-state PINO basis as 21,24 ( Here, the vector X(ω) contains the response of the density matrix where δγ R (ω) is the real part of the change in unique offdiagonal density matrix elements and δn(ω) is the change in the diagonal elements (occupation numbers). The right-hand side of (5) contains a modified perturbing potential (a denotes axis x, y, or z) which is partitioned into its off-diagonal (OD: v a pq , p > q) and diagonal (D: v a pp (ω)) parts. We consider a simple dipole perturbation along the axis a, so v a pq = π 0 p |r a | π 0 q , so that V a is expressed through the dipole integrals. The inverse response matrix χ −1 (ω) in (5) has the form (2) where the adiabatic PINO 0 is expressed as The matrices A + in Eqs. The PINO response matrix χ (ω) is obtained with the inversion of the symmetric matrix in the left-hand side of (5) through the eigendecomposition where F α are the eigenvectors of the adiabatic 0 of (8). Insertion of (9) in (5) with the subsequent inversion of the left-hand-side matrix yields the response matrix With Eq. (10), the dynamic dipole polarizability tensor α ab (ω) can be obtained as minus the dipole-dipole response function Insertion of (10)) in (11) produces the sum-over-states formula for the polarizability The physical interpretation of this sum-over-states representation of α ab (ω) is that its poles are the excitation energies, while the residues yield the corresponding oscillator strengths. Specifically, for the diagonal component of (12) one has For each excitation, the average of the residues of (13) over three dipole directions yields the desired equation for the cor- The product V T a √ A + F α in this equation can be related to the transition dipole moment M τ a (α), which enters the expression for f α , see below. In order to make this relation, we consider the real part of the spectral Lehman representation of the response function where γ R (α) is the real part of transition 1RDM From the comparison of (15) with (10) follows the relation so that, finally, we have the expression of M τ a (α) through the dipole integral vector V a and the transition 1RDM γ R (α). With (17), the master Eq. (14) for f(α) is formally converted to the conventional formula In practice, we obtain the eigenvectors F α from the diagonalization of 0 . The transition density matrix elements, which we use to interpret the corresponding excited state 26 then follow directly from Eq. (17).

III. CALCULATION OF OSCILLATOR STRENGTHS WITH ADIABATIC PINO THEORY
In this section, the results of the adiabatic PINO calculations of the oscillator strengths f α of (14) are presented.
The curves of f α as a function of the bond distance are calculated for the lowest two 1 + u excitations in the H 2 molecule and the lowest two 1 + excitations in the HeH + molecule. At first, the ground-state calculations are carried out with the , which is the energy expression for the electron-electron interaction energy of the two-electron wavefunction in NO basis according to the analysis of Löwdin and Shull (PILS) 27 for closed-shell singlet systems, This expression follows immediately from the spatial part of the CI wavefunction in NO basis, (A completely occupied NO has an occupation number of 2 in our notation.) The prefactor f i is either +1 or −1. We consider the PINOs i can generate the +1 or −1 prefactor with θ 0 i = 0 and π /2, respectively. In principle, the phases can also be optimized, but here we have simply chosen the phases according to the known rules that they usually obey. Specifically, for simple two-electron systems (atoms, molecules at R e ) the phase of the first (heavily occupied) NO is conventionally chosen θ 0 i = 0, f i = +1 (fixing the one free overall phase factor in the wavefunction), and for all the "virtual" (weakly occupied) NOs θ 0 i = π/2, f i = −1. So, the standard DMFT-LS expression for the ground-state electron-electron interaction energy is recovered where f 1 = 1 and f i = −1, i > 1. We also use this Ansatz for the H 2 molecule for R H-H < 5.0 Bohr. However, this sign convention gives errors in the considered excitation energies of up to 2 mHartree when R H-H > 5.0 Bohr. Even larger errors are observed for excitations not considered in this article. This stems from the fact that the used sign convention is erroneous in those cases. From accurate ground state CI calculations, it is known that the signs of the weakly occupied second σ u NO, 2σ u , and first weakly occupied π g NO, 1π g , are +1 when R H-H > 5 Bohr. 28 We use this sign convention for the R H-H > 5 Bohr calculations. These sign changes correct the deficient excitation energy behavior. The minimum of the total energy, which uses W of (21), is found in an iterative process. Each iteration consists of two steps: in the first step, the real NOs φ 0 i (r) are optimized while keeping the occupations fixed and in the second step the occupations n 0 i are optimized. The optimization of the NOs φ 0 i (r) at a given set of occupation numbers uses two strategies: it is initially carried out by (iterative) diagonalization of the effective Fockian for the NOs of Ref. 29. Upon approach of convergence, we switch to a direct optimization scheme using gradients of the NO coefficients. 30,31 The optimization of the occupation numbers at a given set of NOs is carried out with the gradients of the energy with respect to the occupation numbers. At the next stage, when NOs and occupations have been determined, the adiabatic PINO response calculations are performed according to Refs. 21 and 24. The PINO response matrices A + and D (see the Appendix for their definitions) are obtained with differentiation of the functional (20) with respect to the PINOs π 0 i , π 0 * i and their occupations n 0 i . The integrals occurring in these matrices are obtained with the GAMESS-UK package. 32 With these matrices, the eigenvalue problem (4) is solved for the excitation energies ω α and the oscillator strengths f α are calculated from the formula (14). Reference values for the f α are obtained from the response coupled cluster singles and doubles (CCSD) calculations with the DALTON package. 33 All calculations are performed in the aug-cc-pVTZ basis. 34 Note, that full PINO calculations have a higher computational cost than adiabatic TDDFT because the dimension of the response matrices is of the order n(n + 1)/2 if n is the number of basis functions (the TDDFT size is n occ n virt ). However, as was shown in Ref. 26, one can strongly reduce the size of the total response matrix 0 to roughly the same size as the TDDFT response matrix and still retain a good accuracy. To assess the performance of these restricted approaches in calculations of f α , the so-called R0, R1, and R2 restricted variants are used in this paper. These variants are characterized by restrictions on the length of the response vector δγ ia that is used (i.e., restriction on the dimension of the response matrix). We denote the N/2 NOs with the highest occupations as the "strongly occupied" ones, or just the occupied orbitals. Their occupation numbers are all >1, indeed most are close to 2 and only a few will have, upon dissociation, occupations that approach 1. The remaining NOs with occupations <1, usually close to 0 are called the "weakly occupied" orbitals, or simply the "unoccupied orbitals" or "virtuals." The (strongly) occupied orbital with the lowest occupation usually corresponds to the HOMO in the Hartree-Fock model and is denoted the "highest occupied NO" or HONO, the "unoccupied" NO with the highest occupation is denoted the "lowest unoccupied NO" or LUNO. If a response vector obtained for a certain excitation energy ω α has a single large δγ ia element, it corresponds to a singly excited state i → a, see Ref. 26. If two diagonal elements δγ ii and δγ aa are large, we are dealing with a diagonal double excitation from configuration (φ 0 i ) 2 → (φ 0 a ) 2 . All elements related to the occupation number (diagonal density matrix) response δγ pp = δn 0 p , p = 1, . . . , n are always retained in the response vector, so double excitations can be obtained, which is not the case in adiabatic TDDFT. In the variants R0-R2, the off-diagonal elements of δγ are restricted as follows. In R0, one only retains off-diagonal elements δγ ia , where i indexes one of the N/2 strongly occupied NOs. The index a always runs over all NOs. The R1 variant extends the range of the i index to the lowest "unoccupied" natural orbital (HONO+1=LUNO). Both off-diagonal double excitations, to a configuration (φ 0 LU NO ) 1 (φ 0 a ) 1 , and single excitation out of the doubly occupied (1σ u ) 2 (which configuration is present in the ground state wavefunction at elongated distances) into a virtual φ 0 a , are represented with an offdiagonal δγ LUNO,a element. 26 Finally, the R2 variant extends the range of the i index one step further to include LUNO+1. Figure 6 displays the excitation energy curves ω α (R) calculated for the two lowest excited 1 + u states of the H 2 molecule. As was already established in Ref. 26, the R1 and R2 variants reproduce the reference curves fairly well, while the R0 variant has a much larger deviation. The composition of the PINO transition 1RDM γ R (α) in a suitable basis describes the orbital nature of the excitation. Tables I and II compare the expansions of γ R (α) in the NO and KS basis sets for two lowest 1 + u excitations of H 2 at the equilibrium bond distance (R e = 1.4 Bohr). The KS MOs are obtained from the corresponding ground-state calculation with the BP86 DFT functional.

A. Hydrogen molecule
This comparison reveals a striking difference between the KS and NO γ R (α) representations. In particular, in the KS representation the leading term of γ R is about an order of magnitude smaller, while other terms are negligible. In its turn, the leading term of γ R for the second transition, γ R (2 1 + u ), is γ R 1σ g 2σ u (2 1 + u ). This means, that the 1 1 + u excitation can be interpreted as the single KS orbital transition 1σ g → 1σ u , while the 2 1 + u excitation is, mainly, the 1σ g → 2σ u orbital transition. In contrast, the NO representation does not provide such a simple orbital interpretation of the considered excitations. Indeed, for both 1 1 + u and 2 1 + u many NO elements of γ R (α) have a comparable magnitude. The revealed difference between the KS and NO representations is understandable. In the former case, the KS orbitals ψ i are ordered according to their orbital energies i . Then, one can expect that the lowest excitations are produced, predominantly, with the transitions ψ i → ψ a between the frontier orbitals with the smallest energy differences a − i . Indeed, the latter quantity serves as a reliable zeroorder adiabatic TDDFT estimate of the energies of the lowest excitations in compact systems. It is a fortunate property of the KS MO basis, and the KS orbital energies for the virtual (and occupied) orbitals, that excitations can, in contrast to the HF case, often be described as predominantly single orbital transitions, with the orbital energy difference as a very good first approximation to the excitation energy. 16,35 In contrast, the NOs φ i are ordered according to their occupations n i and this ordering provides no information about the relative "energy" of the ordered NOs. The occupation number does not provide information on the possible involvement of an NO in an excitation: a (very) low occupation does not imply at all that the NO is less likely to be occupied in a low-lying excitation. This is corroborated by the results of   Tables I and II for the lowest 1 1 + u and 2 1 + u excitations. This feature of the NO representation makes its use for the purpose of the orbital analysis of excited states rather inconvenient. Our calculations do not find a confirmation of the result that an excited state can be described as a simple change in occupation number of largely unmodified NOs, as proved to be the case for a low-lying excited state in the model calculations of Ref. 37. On the other hand, the KS representation, as observed above, appears to be favorable for an orbital description of excitations, so that only the latter orbital representation of the PINO quantities will be used in the rest of this section. Figure 7 presents oscillator strength curves f α (R) for the 1 1 + u and 2 1 + u states of H 2 . As in the above mentioned case of excitation energies, the R0 variant produces substantial errors for the calculated f α (R), while R1 and R2 perform very well. The only exception is the region R > 7.5 Bohr for the higher 2 1 + u state, where all the restricted variants exhibit appreciable deviations. The remarkable feature of Fig. 7 is a very different shape of the curves f α (R) for the two lowest excited states. The curve f (1 1 + u )(R) has a bell-like shape with the maximum near 3 Bohr. In contrast, the curve f (2 1 + u )(R) exhibits two maxima near 0.5 and 7.5 Bohr and it passes through a near-zero minimum near 4.0 Bohr (see Fig. 7).
To interpret the shape of the f α (R) curves in Fig. 7, the factors V z ia and γ R ia (α) of the transition dipole moment M τ ia (α) = V z ia γ R ia (α) are displayed for the contributing orbital transitions. Table IV presents M τ (1 1 + u ) and its leading components M τ ia (1 1 + u ), V z ia , and γ R ia (1 1 + u ) calculated with R2 at various bond distances of H 2 . The element γ R 1σ g 1σ u (1 1 + u ) steadily decreases with R and it vanishes at R = 7 Bohr. Eventually, this leads to the decreasing oscillator strength f (1 1 + u )(R) for larger R, which is also true for the reference CCSD curve (see Fig. 7). Yet, for all R considered, the 1σ g 1σ u configuration is the main term of the FCI wave function for the excited 1 1 + u state, 1 disappears is a direct manifestation of the strong non-dynamical electron correlation in the ground 0 (1 1 + g ) state of the stretched H 2 . The point is that non-dynamical correlation dictates a strong mixing of the configurations 1σ 2 g and 1σ 2 u in 0 , 0 (1 1 + g ) ≈ c g |1σ 2 g | − c u |1σ 2 u |, c g ∼ c u . As a result, the corresponding terms of the opposite signs compensate each other in the transition Thus, we have the seemingly paradoxical situation that the transition 1RDM misses the "single excitation" element γ R 1σ g 1σ u while the excited state has this singly-excited character. One may call this an accidentally forbidden (due to the strong non-dynamical correlation in the ground state) 1σ g → 1σ u transition. This trend derived from the character of the CI wavefunction is faithfully reproduced by the transition 1RDM γ R 1σ g 1σ u of the PINO response theory. The shape of the oscillator strength curve (Fig. 7) can now be understood as follows. For the 1σ g → 1σ u transition, it is easy to see, using the fact that these orbitals are approximately (1s a ± 1s b )/ √ 2S ± 2, that the dipole integral V z 1σ g 1σ u = dr z 1σ g (r)1σ u (r) becomes proportional to the bond length R. Therefore, as long as non-dynamical cor-relation is not strong, which is the case till ca. 3 Bohr, the transition dipole M τ (α) (the square of which enters the oscillator strength f(α)) increases because of the increase of the V z 1σ g 1σ u factor. At longer distances, however, the 1σ g → 1σ u transition becomes accidentally forbidden due to the strong nondynamical correlation in the ground state and γ R 1σ g 1σ u (1 1 + u ) vanishes by virtue of Eq. (22). This combination of the increasing V z 1σ g 1σ u and decreasing γ R 1σ g 1σ u (1 1 + u ) fully explains the observed bell-like shape of the oscillator strength f (1 1 + u )(R) (see Fig. 7). We note that at larger distances also other excitations enter the 1 1 + u state, and affect the oscillator strength, cf. the 1σ g 2σ u , 2σ g 1σ u and notably 3σ g 1σ u contributions at 7 Bohr. Table V presents M τ (2 1 + u ) and its leading components M τ ia (2 1 + u ), V z ia , and γ R ia (2 1 + u ) (evidently, V z ia does not depend on the excited state, so that the corresponding columns of Tables IV and V are identical). This state has predominant 1σ g → 2σ u character at all distances, although at very long distances (cf. 7 Bohr in the Table) considerable admixture of other transitions (notably, 1σ u → 2σ g ) occurs. However, unlike the increasing dipole integral V z 1σ g 1σ u mentioned above, the dipole integral V z 1σ g 2σ u passes through a near-zero minimum at R = 3 Bohr (see Table V). This is caused by the atomic orbital character of the 2σ u . For instance, when this orbital becomes predominantly 2s a − 2s b the dipole integral becomes proportional to R 1s a |2s b , but the two-center overlap 1s a |2s b goes to zero. At larger R values, exemplified by the 7 Bohr entry in the Table, the character of the excited state becomes strongly mixed and several contributions to the tran-  sition dipole moment arise, which make it deviate from the near zero value at ca. 4 Bohr. This explains the observed shape of the reference curve f FCI (2 1 + u )(R), which passes through a near-zero minimum at R = 4 Bohr, while the PINO R1 and R2 curves faithfully reproduce this shape (see Fig. 7).

B. Failure of adiabatic TDDFT
We can now explain why the adiabatic TDDFT oscillator strength for the lowest excited state, 1 1 + u , in Fig. 2, although not quantitatively accurate, appears to give qualitatively the correct behavior as function of the distance, in spite of the completely wrong TDDFT excitation energy. [9][10][11][12] This, indeed, is caused by a fortuitous cancelation of errors. The influence of non-dynamical correlation on the shape of the KS 1σ g and 1σ u orbitals participating in the leading 1σ g → 1σ u transition does not affect the above mentioned proportionality to R of the dipole integral V z 1σ g 1σ u and this leads to a much too large transition dipole M τ TDDFT (1 1 + u ). However, the excitation energy ω TDDFT (1 1 + u ) goes (erroneously) to zero. Eventually, this last factor (cf. Eq. (19)) makes the oscillator strength decline to zero. We conclude that this decline in the oscillator strength beyond 3 Bohr is not caused, as it should, by the change in character of the ground state due to nondynamic correlation at long distance, but happens because the excitation energy tends erroneously to zero. Figure 8 displays the excitation energy curves ω α (R) calculated for the two lowest excited 2 1 + and 3 1 + states of the HeH + molecule. In this case, already R0 reproduces the shapes of both reference curves rather closely, while the R1 and R2 curves virtually coincide with the reference ones (see Fig. 8). Judging from the transition 1RDM elements γ R ia (2 1 + ) in the KS basis presented in Table VI, the 2 1 + excitation can be interpreted at all R considered as, mainly, the orbital transition 1σ → 2σ , which corresponds to the CT excitation 1s(He)→ 1s(H). In its turn, judging from the elements γ R ia (3 1 + ) of Table VII, the 3 1 + excitation can be interpreted as, mainly, the orbital transition 1σ → 3σ , which also corresponds to a CT excitation, this time the 1s(He) → 2s(H) excitation. Figure 9 displays oscillator strength curves f α (R) for the 2 1 + and 3 1 + excited states of HeH + . In this case, only R0 displays a visible error of the curve f(2 1 + )(R) in the interval R < 4 Bohr, while in all other cases the restricted PINO curves virtually coincide with the reference ones. Interestingly, the f α (R) curves for the two lowest states of HeH + have a similar shape as those for the two lowest states of H 2 (Compare Figures 7 and 9). In particular, the curve f(2 1 + )(R) has the bell-like shape with the maximum near the equilibrium bond distance R e = 1.46 Bohr. In its turn, f (3 1 + )(R) passes through a near-zero minimum near R = 4 Bohr and it exhibits a maximum near 5.5 Bohr. Both f α (R) vanish at longer R (Fig. 9).

C. The HeH + molecule
To interpret the shape of the f(2 1 + )(R) curve, the corresponding total transition dipole moment M τ (2 1 + ), the leading contributions M τ ia (2 1 + ) as well as the factors V z ia and γ R ia (2 1 + ) are presented in Table VI. The leading transition 1RDM element γ R 1σ 2σ (2 1 + ) does not change much with R, so that the shape of f(2 1 + )(R) is, apparently, determined by the corresponding dipole integral V z 1σ 2σ and the excitation energy ω(2 1 + ), the latter influences f(2 1 + )(R) through (19). Then, the combination of the decreasing ω(2 1 + ) (Fig. 8) and the increasing (for R ≤ 2.5 Bohr) V z 1σ 2σ (See Table VI) produces the maximum of f(2 1 + )(R) near R e = 1.46 Bohr (Fig. 9) (Table VII), thus causing the resultant f(3 1 + )(R) also pass through a near-zero minimum near this point (Fig. 9). The dip in the oscillator strength at R = 4 Bohr is caused by configuration mixing: the 3 1 + state acquires majority 1σ → 4σ character here, and its contribution to the transition dipole is counteracted by opposite sign contributions from the remaining 1σ → 2σ , 1σ → 3σ , and 1σ → 5σ excitation character. At large R V z 1σ 3σ vanishes due to the vanishing overlap of the orbitals 1σ ≈ 1s(He) and 3σ ≈ 2p z (H) localized on the different atoms, so that the resultant f(3 1 + )(R) also vanishes ( Fig. 9).

IV. CONCLUSIONS
In this paper, the previous development of the TDDMFT and PINO adiabatic response theories for electronic excitation energies ω α is completed with the derivation of the expression for the corresponding oscillator strengths f(α). The latter quantity and the related transition dipole moment M τ (α) are expressed through the eigenvector F α and the A + part of the coupling matrix of the PINO eigenvalue equations for ω α as well as through the dipole integral vector V a .
The main emphasis of this paper is on the quality of the oscillator strength curves f α (R) for various bond distances R, an important subject, which has received surprisingly little attention in the literature. It is well-known that conventional (adiabatic) TDDFT yields wrong excitation energies upon bond elongation. 9,11,12 In the present paper, the problem that TDDFT has with reliably estimating the oscillator strengths f α (R) for dissociating bonds is highlighted. For the prototype case of the 1 1 + u excitation in the stretched H 2 molecule, the qualitative resemblance of f TDDFT (1 1 + u )(R) to the reference curve is achieved only because of a compensation of unphysical errors in various components of f TDDFT (1 1 + u ). The highlighted problem of TDDFT is addressed with the present application of the PINO theory to the calculation of f(α). The extension of adiabatic TDDMFT to the phaseincluding adiabatic PINO leads to the exact 1RDM response theory for two-electron systems. Then, full adiabatic PINO with the ground-state PILS functional (in which the PINO phases are properly chosen) produces accurate f α (R) curves for the lowest 1 + u excitations in the prototype "two-electronbond" molecule H 2 and for the lowest 1 + CT excitations in HeH + . Already, the restricted R1 and R2 PINO variants with the strongly reduced size of the eigenvalue problem produce high quality f α (R) curves in these cases.
An interesting result of this paper is the very different shape in both H 2 and HeH + of the f α (R) curves calculated for the lowest excitation (1 1 + u and 2 1 + , respectively) and the second lowest excitation (2 1 + u and 3 1 + , respectively). For both molecules, the f α (R) curves of the lowest excitations have a bell-like shape, while those for the second lowest excitations pass through a near-zero minimum near R = 4 Bohr. This difference has been traced to the (changes in the) composition of the ground state and the excited states along the dissociation coordinate.
The interpretation of the nature of excitations is based on the orbital-pair composition of the transition 1RDM γ R (α), while the established trends of the f α (R) curves have been analyzed with an orbital decomposition of the transition dipole moment M τ (α). The corresponding orbital analysis is carried out in the representation of the KS MOs, since each lowest excitation in question can be conveniently interpreted (for the majority of the considered bond distances) just as a transition between a single pair of frontier KS MOs. The NO representation lacks this convenient feature and quite a few NO transitions contribute significantly to a given lowest excitation.
Further application of the PINO theory to the calculation of ω α and f(α) for larger molecules requires the development of an approximate N-electron functional, which, just as the two-electron PILS functional of the present paper, would provide a good quality of ω α (R) and f α (R) with the restricted PINO calculations. This work is in progress.