Continuum Electronic States: The Tiresia Code

A multicenter (LCAO) B-spline basis is described in detail, and its capabilities concerning affording convergent solutions for electronic continuum states and wavepacket propagation are presented. It forms the core of the Tiresia code, which implements static-DFT and TDDFT hamiltonians, as well as single channel Dyson-DFT and Dyson-TDDFT descriptions to include correlation in the bound states. Together they afford accurate and computationally efficient descriptions of photoionization properties of complex systems, both in the single photon and strong field environments. A number of examples are provided.


Introduction
Quantum Chemistry (QC) approaches to the calculation of electronic bound states of complex molecules have reached a large degree of sophistication, and have given rise to a multitude of highly advanced computer codes. The treatment of continuum states to describe photoionization and electron scattering has significantly lagged behind, probably because of unfamiliarity with continuum states in the chemistry community, the unsuitability of bound state algorithms for the determination of approximate solutions, and the limitations of the ubiquitous Gaussian type orbital (GTO) basis sets in describing continuous states. A strong revival of interest in molecular continuum states, especially for ionization, has been spurred by the development of novel light sources, high intensity laser and Free Electron Lasers (FELs), capable of ultrashort pulses, together with the continuous development of Synchrotron Radiation (SR) sources, and the parallel development of ever more performant detectors for electrons and ions, which have produced a flurry of new experiments affording unprecedented detail in studies of light-matter interaction in molecules [1].
In the treatment of continuum states, two new problems appear. The first is the accurate solution of the one-particle Schrödinger equation for positive energies. The second is the treatment of many-electron systems, compatible with the asymptotic structure of the wavefunction. Initial studies on atoms naturally employed numerical radial grids. While this approach can be generalized to a one center expansion (OCE) of the wavefunction in the molecular case [2,3], it suffers from the very slow convergence of the partial wave expansion for heavy atoms far from the expansion center. Generalization to 3D grids is difficult, because a natural coordinate system is spherical around each nucleus, due to the Coulomb singularity and the relative cusps in the wavefunction. In the Continuum Multiple-Scattering Method (CMS) shape approximation of the potential allows efficient numerical integration [4]. Other adaptive grids have been suggested [5]. Cartesian grids are, however, in use, employing pseudopotentials to smooth the Coulomb singularities, especially for the treatment of time dependent phenomena [6]. Another difficulty is the implementation of asymptotic boundary conditions, which are most neatly expressed in a partial wave expansion. The obvious solution is suggested by the well established LCAO approach in molecules. The superposition of spherical basis sets centered on the various atoms effectively mimics the change of coordinates from spherical around one site to spherical around another, dealing very effectively with coulomb singularities. In the continuum case these have to be added to the OCE expansion, which is needed to describe the long range part of the wavefunction. This prompts the use of a basis set approach, requiring, however, new basis functions, capable of approximating the oscillatory behavior of continuum wavefunctions. Note however that the expansion can be limited to a finite spherical volume, large enough that all bound states of interest are negligible beyond this boundary, and the continuum states have reached their asymptotic form, so that they can be fitted to analytical forms (spherical Bessel or Coulomb functions), or to the solutions of a weak multipolar field by an OCE approach (Figure 1). Recent work has converged to the use of local polynomial functions as radial basis, finite elements [7] or B-splines [8][9][10], times spherical harmonics, at least for the long range OCE part, although GTOs have been used for small ranges and low electron kinetic energies [9]. For the short range, multicenter part, often GTOs are used, in various combinations [9,10].
The second issue is the treatment of the many electron problem. A single determinant description of the bound states is easily extended to the continuum, replacing a bound state orbital with a continuum one, obtained by solving a resultant one particle equation, like Frozen Hartree-Fock (HF) or Density Functional (DFT). If |Φ 0 denotes the HF or DFT ground-state of the N-electron system: then the set of degenerate solutions at energy E > 0 can be written as: where index i in Equation (2) labels the MO from which the electron is ejected. In HF theory, the approach is called Static Exchange (SE), while it has been termed SE-DFT in DFT theory.
Coupling among channels is implemented as Configuration Interaction Singles (CIS), also known as the Tamm-Dancoff approximation (TDA). Some initial state correlation can be included at the linear response level, giving the well known random phase approximation (RPA) or equivalent Time Dependent DFT (TDDFT), in the ab initio or DFT framework, respectively. This can be generalized to the use of correlated bound states, both initial and final, giving rise to the so called close-coupling form of the wavefunction In Equation (3), the first term comprises antisymmetrized products of N-1 particle bound eigenstates Ψ N−1 J , describing the final states of interest (target states), times continuum orbitals φ EI Jlm , and ensures the correct asymptotic limit. The second part describes short range correlations and autoionization resonances. Actually, some restrictions have to be included to avoid linear dependencies and ensure unique solutions [11]. Given its generality, it can be considered as the equivalent of CI in the continuum, and several variants can be designed depending on specific forms for the bound states included in the expansion. For instance, restricting the target states to single hole states relative to an initial closed shell state, and ignoring correlating functions, one recovers the CIS ansatz, or by including a single, correlated target state Ψ N−1 I , relative to the final state of interest, one has the Single Channel (SC) approximation.
At this level, correlation is fully included in the bound states, but neglected in the continuum. The close-coupling expansion can in principle be very accurate, but it is computationally expensive, and simpler approximations are often employed.

The Basis Set
Our current implementation uses B-splines for the OCE and also for the LCAO part. Bsplines (B is for basis) are functions defined over a set of points (knots) which divide a given interval, [0, R max ] in our case, in n int subintervals. Over each subinterval a B-spline is a different polynomial of a common fixed order n (degree n − 1), spanning n + 1 consecutive knots, and is zero outside. Different polynomial pieces are joined at the boundaries with a high degree of continuity, typically the maximum admissible, i.e., C (n−2) , continuity up to the derivative of order n − 2, but that can be relaxed employing coincident knots. In fact in the present implementation continuity is completely relaxed at the endpoints, 0 and R max , putting here n coincident knots, so that the first and last spline are nonzero there, while employing single knots, i.e., maximum continuity, at the interior points. This ensures a complete basis, entirely contained in the interval [0, R max ]. Actually the first spline is deleted, so as to enforce the boundary condition f (0) = 0. The same is done at R max when dealing with bound states, assumed to be negligible at the outer boundary. The last spline is instead retained in continuum calculations, as the solution naturally oscillates and is nonzero at the boundary. With this choice the dimension of the basis is n + n int − 1, so the B-spline order has a very limited impact on the dimension, but higher orders vastly increase accuracy, before degrading due to numerical truncation errors. Bsplines as well as their derivatives are easily evaluated at a given point, and integration by Gauss-Legendre quadrature is accurate down to machine precision. A full discussion including fortran routines is available in the monograph [12], and a detailed treatment of atomic and molecular applications can be found in [13].
In the present case the B-spline basis is completely defined by the polynomial order n (default value is n = 10) and the knot sequences that define intervals [0, R max ] and [0, R maxi ] of the OCE and LCAO expansions around each atom. In the OCE sequence it is important to have a knot at the radial position of each atom, as the singularity of the potential there degrades the accuracy of the numerical integration otherwise. As the primitive functions are the complete basis is specified by the maximum L values of the spherical harmonics for the OCE, L max and the LCAO part, L max,i . The OCE part is numerically very stable; in normal photoionization typical R max around 20-40 atomic units (au) are employed, but values up to several thousands have been employed for long wavepacket propagation without trouble. For photoionization a linear grid of knots is adequate, with spacing h being determined by the maximum electron momentum k max (electron kinetic energy KE = k 2 max /2). A choice h = 1/k max is generally very accurate. L max has to accommodate both convergence of the expansion in the molecular region, and the number of continuum partial waves carrying significant contribution at k max . For the latter, L max is roughly proportional to k max times the "size of the molecule", by the classical argument based of momentum times impact parameter. So large systems and high electron energies may require large L max , and small step h, making the basis very large. In the case of heavy atoms, steps have to be reduced close to the origin. If only bound states are needed, e.g., high Rydberg states, an exponential grid becomes more convenient. Then, with appropriate range, quantum numbers above 100 can be reached.
In contrast, LCAO R max,i are generally very small, of the order of about 1 au. As corresponding L max,i are the usual values of the maximum L in the free atom, or L + 1, the dimension of the LCAO basis is generally very small compared to OCE, which alone determines the computational cost of the calculation. The very short range of the LCAO expansion is dictated by the need to avoid overcompleteness, i.e., numerical linear dependence of the basis, which destroys numerical stability. Actually, the overcompleteness problem is even more severe with GTOs, which have a very large overlap. The ease of selecting the range of B-splines makes it easier to control overcompleteness. There is actually a bargain between LCAO radii and OCE L max , as with larger L max shorter R max,i have to be employed, and with low L max larger R max,i make up for convergence in the molecular region. Limiting the LCAO basis to non-overlapping ranges in the present implementation does not constitute a limitation and has the additional benefit of reducing the number of matrix elements to be evaluated, as products of LCAO functions on different centers are strictly zero, B i (r p )B j (r q ) = 0. The last LCAO splines (usually 3) are then deleted, so that the basis has overall continuous second derivative at all internal points. All two center integrals between OCE and LCAO functions are accurately evaluated by numerical integration within the LCAO spheres.
The basis is then fully adapted to point group symmetry, if present, which brings a great reduction in the dimension of the matrices and the computational effort for the OCE basis, and for the LCAO part. Here (λµ) are the labels for the irreducible representations, j ≡ lh counts the number of total angular independent components, with h counting the independent components relative to a given l value. For the LCAO functions, p counts the number of sets (shells) of symmetry equivalent atoms, and q runs over the atoms in this set. The final basis is then The final basis has a lot of flexibility, and can afford good convergence for a rather wide range of parameters. It requires, however, some care, as it is easy to overdo, which is signalled by the minimum eigenvalues of the overlap matrix S that become very low. It is then necessary either to restrict the range of LCAO functions or decrease the OCE L max . The upside is that it is possible to accomodate a large variety of situations.

The Galerkin Approach
The one particle Schrödinger equation: is converted to the generalized eigenvalue problem when the solution φ is expanded as a linear combination of basis functions In Equation (10) h and S denote the Hamiltonian and overlap matrices, respectively, of dimension equal to the number of basis functions employed in the expansion of φ. Bound state eigenvectors are obtained from a usual diagonalization. It is important to have accurate initial orbitals, in particular in the outer tails, which are often not well described by conventional basis sets, like GTO, because the transition matrix elements for the continuum prove sensitive to them. In the case of continuum eigenvectors several algorithms can be used. Note that in this case, because of the last nonzero spline at the end of the interval, the Hamiltonian matrix, h, is non-Hermitian. We employ the Galerkin approach, originally proposed in [14], which amounts to requiring orthogonality of the residual vector (h − E)φ to all basis functions. Defining the energy dependent matrix A(E) = h − ES, this amounts to solving: for the n o lowest eigenvalues, where n o is the number of open channels, i.e., the number of independent solutions. Actually, an alternative formulation [15,16], employing the hermitian product A T A is more stable and accurate. It turns out that, with a good basis, the gap between the lowest n o eigenvalues and the following ones is many orders of magnitude, and the eigenvectors are easily extracted with a block inverse iteration, which generally converges very fast. The resulting continuum eigenvectors are fitted to a linear combination of regular and irregular asymptotic functions and are then normalized to the K-matrix asymptotic form. Full details are in [15]. Note that, as the basis is close to complete within the sphere, any type of linear equation, either homogenous or inhomogeneous can be accurately solved in matrix form. We currently employ this technique to evaluate the Coulomb potential generated by a charge distribution solving the Poisson equation, to solve inhomogenous equations in TDDFT and so on.

The Static Exchange Approach
The Static Exchange (SE) approach is defined by the use of the same (frozen) orbitals for the ground state as for for the ionic states. The simplest approach is the DFT-SE scheme. We employ a fixed Kohn-Sham (KS), or DFT, single particle hamiltonian where V N is the nuclear attraction potential, V C is the classical Coulomb (Hartree) potential generated by the electronic density ρ(r) and V XC [ρ(r)] is a local exchange-correlation potential defined by the same density. The density is that of the ground state (GS), obtained by a previous self-consistent field (SCF) calculation with a separate QC program. The h and S matrices are computed in the B-spline basis, and the occupied orbitals obtained by diagonalization, defining the initial state in Equation (2). At each energy E > 0, n o independent degenerate solutions are obtained by inverse iteration: with the asymptotic behaviour (close to R max ): In Equation (15), the radial functions R El m lm (r) can be written as: where f l and g l are the regular and irregular radial solutions in the asymptotic region, and the matrices A and B are obtained by fitting the numerical solution to the asymptotic form at the boundary (we use the two last interior knots [14]). Multiplying the solution by A −1 gives the final K-matrix normalized form: From the final state Φ Eilm (Equation (2)) the transition dipole matrix elements are evaluated as: where γ is the dipole operator component. These, together with the K-matrix, are passed to a separate program that evaluates photoionization cross sections and angular distributions.

The TDDFT Approach
We employ first-order linear response, which describes the screening of the external exciting field due to the response of the electron cloud. The fundamental equations governing linear response are where δρ is the first-order change in electron density induced by the field, χ is the electric susceptibility, δV is the first-order change in the potential induced by a density change δρ, and K is the linear kernel relating the two quantities. The final full potential (V SCF ) is the sum of the external perturbing potential (µ ext ) and the response potential (δV).
In the response framework, the effect of electron polarization gives rise to an effective potential, which is the sum of the external dipole field and the electronic response potential. The evaluation of the photoelectron matrix element then boils down to replacing the pure dipole field in Equation (18) by the effective potential V SCF : From Equation (19) one can obtain an equation for V SCF We then solve directly for the response potential V SCF as the prime dynamical variable, using a non-iterative numerical algorithm [17]. By representing V SCF as well as the operators K and χ in the B-spline basis, the integro-differential equation, Equation (21), is converted into a linear algebraic one The kernel K is the sum of the Coulomb potential and the linearized exchangecorrelation response The computationally expensive part is the calculation of the KS (non-interacting) linear susceptibility, which is energy dependent, and is computed for each photon energy by solving the inhomogeneous first order perturbative equations.
where Q in Equation (24) is a projector over the space of KS virtual orbitals and φ ±1 i are first order perturbations of the KS orbitals. The completeness of the basis ensures convergent solutions. Full details are available in Ref. [17].

The Correlated Single Channel Approach
One can improve the treatment of bound state correlation by employing fully correlated initial and final target bound states, computed by ab initio CI or equivalent approaches in a single channel framework, Equation (4). Here we label Ψ N I the initial state and Ψ N−1 F the target state. The SC approximation neglects correlation between the target and the continuum, between different channels (interchannel coupling) and autoionization resonances. However, important correlation effects, especially in the inner ionic states, are fully accounted for. Moreover, it deals equally well with multiconfigurational states, like ionization from excited initial states, and multielectron excitations in the final states.
In this approximation the many-particle transition moment reduces to: where and φ D FI and η D FI are called the Dyson orbital and the conjugate Dyson orbital, respectively, relative to the two bound states Ψ N−1 F and Ψ N I . So the full transition matrix element reduces to a single particle transition element between the Dyson and the continuum, plus an overlap term of the continuum with the conjugate Dyson [18]. The latter is, moreover, generally neglected, as it is expected to be very small (it is in fact zero in the SE or TDDFT approaches, as the continuum is rigorously orthogonal to the bound eigenstates of the same hamiltonian).
The continuum relative to a given target state can be computed variationally or can be separately approximated. In the present developments we use the same continuum as in the SE-DFT scheme, which is already quite accurate. We call it the Dyson-DFT approach [19]. In practice this amounts to substituting the Dyson orbital for the bound DFT orbital of the SE approach. Dyson orbitals can be computed by several ab initio approaches. A file containing the details of the ab initio data and Dyson orbitals is provided, and the latter, expressed in GTOs, are then reexpressed in the B-spline basis by projection, which is again very accurate, for the ease of computing transition moments with the final continuum. In the same spirit, one can couple the SC treatment, i.e., the Dyson orbital, to the TDDFT continuum, i.e., substituting in the computation of the transition matrix element the external perturbation µ ext γ , i.e., d γ with the perturbed one V SCF γ , giving: called the Dyson-TDDFT formulation [20]. This hybrid approach proves quite effective when the additional correlation included in TDDFT proves important. One important issue is the one-electron effective potential used in the continuum calculation. Our choice for the DFT V XC potential is the V LB94 potential from [21], which was designed to implement an asymptotic coulomb tail, which proves effective for photoionization [22]. Earlier a transition state density (i.e., half electron removed from the ionized orbital) V X α [23] or V VW N potentials [24] were employed. V LB94 proves at least comparable or better, and has the additional advantage that the same potential proves adequate to cover the entire ionization spectrum that is from outer valence to deep core states. The potential can be truncated to an analytical coulomb tail. This is necessary, for instance, with very large R max , when V LB94 becomes numerically unstable. Finally, the level of correlation for the treatment of bound states in the Dyson approaches is also important. Despite recent advances, highly correlated ab initio approaches are difficult to afford for relatively large molecules and good basis sets. This is an area where further work is certainly warranted.

Multiphoton and Strong Field Processes
Let us consider for a moment the differential cross section for photoionization, which is expressed via a transition (dipole) matrix element between the initial and the final continuum (or bound) state which can also be seen as a projection of a "final wavepacket" or "prepared state" Φ DI onto a field free continuum (or bound) final state, eigenvector of the free hamiltonian. The same logic can be employed for wavepackets relative to more complex processes. One can consider separately 1.
Calculation of Φ DI , which is L 2 and depends on the particular excitation mechanism considered. 2.
Calculation of Ψ FElm , which is independent of the former, as eigenvector of the molecular hamiltonian already considered.
Projection on free hamiltonian eigenstates amount to wavepacket analysis, giving the corresponding transition probabilities.

1.
In the multiphoton domain, generally the lowest order perturbation theory is employed valid for a single photon field, and non-resonant intermediate states, but can be easily generalized, and can be evaluated recursively.

2.
For nonperturbative fields a standard approach is now the solution of the time dependent Schrödinger equation (TDSE) with a time dependent hamiltonian which includes the external field.
Both approaches have become rather standard for computational simulations of the relevant processes. As the time evolution is unitary, at all times the wavepacket is square integrable, and can be accurately described by the B-spline basis as long as its amplitude at the boundary remains negligible. This may require very long range bases, spanning hundreds or even thousands of atomic units [25], which can, however, be easily accomodated in the B-spline basis. Techniques like the use of complex absorbing potentials at the boundary can reduce the need of such extended ranges [26].
One typical approach is employing a spectral method, i.e., a preliminary full diagonalization of the free molecular hamiltonian, so that the resolvent can be expanded as and the amplitude (30) can be easily evaluated. Similarly, the TDSE equation, after expansion of the wavepacket is reduced to a well conditioned system of ordinary differential equations which can be solved by a variety of approaches.
In the static-DFT framework one is working with an independent particle hamiltonian Φ I = |φ I 1 · · · φ I N E I = ε I 1 + · · · + ε I N (37) The single particle h KS can be fully diagonalized even with large B-spline bases. Then all equations and matrix elements reduce to one-particle problems and large systems become affordable. While this is the most basic level, as it lacks correlation among the states, it is capable of giving a generally semiquantitative description [27]. Note that, as the coupling with the external electromagnetic field, h(t) = h KS + V(t), is a single particle operator, within this approximation the full hamiltonian remains of single particle type, and it is easy to show that the full N-particle TDSE is equivalent to the solution of individual one-particle TDSE for the initial orbitals This means that the single particle states φ i are propagated independently, without any restriction, and due to the unitarity of the evolution they remain orthonormal at all times. Moreover, as the exact solution of a model hamiltonian, the result is gauge independent, provided the basis is sufficiently accurate, as has been well verified numerically.
Final projection of the wavepacket Φ DI , either from multiphoton or time propagation, is straightforward, with the corresponding continuum functions already available in the same B-spline basis. In other formulations, techniques based on the evaluation of the electronic flux at the boundary can be employed [26].
It may be remarked that most algorithms reduce to straightforward operations of linear algebra that can be parallelized efficiently also with the use of optimized libraries and scale very well.

Ab Initio Developments
The B-spline basis maintain excellent performance in an ab initio setting [28]. Actually, most algorithms employed remain unchanged, the main difference being the choice of the hamiltonian, or implicitly the structure of the many electron basis, and the evaluation of matrix elements. An obvious choice is the close-coupling approach, but other structures like multiconfigurational SCF can be equally employed.
Clearly, a prerequisite is the availability of two electron integrals, relative to the Coulomb interaction, in the LCAO basis relative to the electron densities ρ ik (1) = φ i (1)φ k (1) and similarly for ρ jl (2). The one-center radial two electron integrals have been long available [29,30], and the angular part is analytical, as is well known in atomic theory. Actually, for not too large systems, that may be a reasonable compromise, as it is partly used in some current well known codes [2,3]. The full evaluation of LCAO B-spline two electron integrals, essentially by numerical approaches, is demanding but by no means unfeasible. One straightforward approach is full numerical integration. Due to the limited range of functions on atomic centres, all integrals will have one density (particle 1, ρ(1)) restricted to an atomic sphere, whose Coulomb potential may be easily evaluated in multipolar form and integrated over the second particle density ρ (2), which in many cases is also restricted to an atomic sphere. A second approach is to solve the Poisson equation for the first density. As the basis is very accurate, so is the computed potential. It requires implementation of multipolar boundary conditions at R max , which is easily implemented with B-splines. After the evaluation of the the expansion of the density in the basis, ρ µ1 = χ µ |ρ 1 (1) the corresponding potential is obtained by solving the linear system V 1 is then integrated with the second density ρ 2 (42) which reduces to the last equation if ρ µ2 = χ µ |ρ 2 is also available. This approach is particularly effective if not all integrals over primitive χ µ functions are required. In a close coupling expansion only integrals of the type are required and the relative smooth ρ 1 densities are accurately approximated in the Bspline basis. Moreover, this avoids the need for integral transformation from primitive to the MO basis. A general routine for these two electron integrals has been recently completed and will be separately reported with full details [31]. It is used in the B-spline ADC(2) code developed by Ruberti [32], mainly devoted to strong field and attosecond processes.
The second main point is the evaluation of the matrix element in the chosen N-particle basis, including dealing with nonorthonormality. In the case of close-coupling expansions, explicit formulas entail up to three particle transition density matrices between the bound states {Ψ N−1 I , Ψ N K } [7,9]. Some implementations have been discussed in the literature, and represent generally a computationally expensive task, although special algorithms may afford significant speedup. Explicit formulas have been derived for the lowest level, CIS wavefunctions, and have been implemented at the OCE level [28]. Implementation with RASCI wavefunctions is currently studied. These will be available in a future release of Tiresia. Moreover, separate treatment of Auger and resonant Auger decay [33] via golden rule formulation is currently underway.

Molecular Photoionization
The most common application is the description of single photon molecular photoionization. This is interesting per se, as it is an important probe of molecular electronic structure, or as an initial step of strong field or ultrafast phenomena, as well as a probe in time dependent studies. From the transition amplitudes and K-matrix the cross section observables σ and β, which define the differential cross section for randomly oriented molecules, and linear polarization can be obtained. P 2 (cos θ) is the second Legendre polynomial and θ is the angle between light polarization and photoelectron momentum. More complex angular distribution, for partially aligned or oriented molecules, nondipole contributions, etc, can also be obtained. Photoionization of C 60 has been investigated several times [34][35][36][37]. Illustrative results from recent studies are reported in Figure 2. A comparison of DFT and TDDFT results show a very similar shape of the total cross section, with a significant increase of the absolute value in TDDFT (left upper panel). Note that the latter incorporates oscillator strength due to discrete excitations above threshold that do not appear in the DFT cross section but show up as autoionization resonances in TDDFT. Actually the fixed nuclei cross section shows an amazing array of very sharp structures, arising both from shape and autoionization resonances. These are washed away by nuclear motion, and we have actually broadened the computed profiles for a better comparison with the experiment, which is reproduced quite accurately. To show, however, the high resolution attainable, a very fine energy scan of the total DFT core cross section is reported in the central upper panel for the first 20 eV KE. Considering that only shape resonances can appear at the DFT level, the presence of these sharp structures is very surprising. Clearly a specific feature of C 60 , it hints at the presence of very localized quasi bound states in the one-particle continuum. The absence of dramatic effects in going from DFT to TDDFT results appears to disprove the claim of giant plasmonic effects in C 60 , which derives from results obtained with jellium models, which smear the carbon cores over a spherical shell. Indeed, seen on a large energy scale, total valence and core ionization of C 60 follows closely the sum of 60 free-atom cross sections [37]. TDDFT becomes very close to DFT some 150 eV above threshold. Indeed the jellium models, lacking the hard cores which are responsible for the high energy cross section, give results that decay much too fast with increasing KE. Still, a lot of structure survives up to pretty high energies, as strong oscillations in the cross section. These are best seen in individual ionizations, like the highest occupied molecular orbital (HOMO) 6H u channel (left lower panel). The phase of oscillation is opposite for states of opposite parity, so that their cross section ratio enhances the oscillations (right lower panel). These have been detected since the first photoemission studies in C 60 . The DFT results match closely the experimental data available up to about 300 eV [38]. Further results suggest that the oscillations persist undamped up to at least 600 eV, after which they become irregularly damped [37]. However, nuclear motion neglected in this study may induce some earlier damping. Oscillations in the cross sections are indeed quite ubiquitous, caused by diffraction effects due to multicenter emission of the electron wave [39][40][41][42]. Further structures at high energies are associated with non-dipole effects [43].
Correlation effects in bound states due to ionization have been noticed and explained long ago [44,45], especially as concerns ionization energies and the appearance of satellite states. One subtle manifestation is a modification of the orbital structure in the ion. In most small molecules the orbital structure is quite rigid, and ground state orbitals, either HF or DFT, change very little in the ion. However, in larger systems, a so called hole-mixing effect, associated with a rotation of the GS orbitals in the ion, can be relevant [19,46]. An example is in low symmetry molecules, where many MOs can remix without symmetry constrains. Another typical situation is in transition metal compounds, where the mixing between metal d and ligand orbitals can change dramatically by electron removal. Whie often signalled by large changes in ionization energies with respect to Koopman's values, the orbital change, which is embodied in the corresponding Dyson orbital, can be detected only by transition properties, like photoionization observables or electron momentum profiles [47]. A recent study on a prototypical system, bis-allyl nickel Ni(C 3 H 5 ) 2 [19], shows the large changes in σ and β parameters from those afforded by HF or DFT orbitals ( Figure 3). Moreover, it shows the value of photoionization observables concerning aiding in the difficult problem of assigning the spectrum in strongly correlated systems, e.g., the appearance of the strong 3p → 3d autoionization resonance in ionizations involving metal 3d participation. A new parameter β 1 appears in the angular distribution of photoelectrons in photoionization of chiral molecules with circularly polarized light: where m r is +1 or −1 for left and right circular polarization. This can be measured and is an important tool for the study of molecular chirality. It arises purely from an electric dipole transition, and the effect is orders of magnitude stronger than circular dichroism, which is very difficult to measure for diluted samples in gas phase. DFT calculations have often afforded a fair reproduction of the experimental results [48]. However, the effect is largest close to the threshold, and dies after about 50 eV KE. Unfortunately the threshold region is the most difficult to describe, and is very sensitive to the specific DFT potential employed, as illustrated in the case of camphor HOMO ionization in Figure 4. The left upper panel taken from [49] shows results from a multiple scattering X α calculation and present LB94 results, compared to experiments. Other panels show results obtained with additional choices of the potential, employing GS or transition state densities. A further problem is the proper treatment of correlation effects, e.g., in the case of chiral transition metal molecules, which are still pretty large for higher level ab initio approaches. A case is photoionization of Cobalt tris-acetylacetonate. The photoelectron spectrum, with superimposed KS eigenvalues from the LB94 potential, and IEs calculated with the outer valence greens function (OVGF) approach [44], a usually pretty accurate treatment for ordinary organic molecules, is presented in Figure 5. While HOMO ionization is reasonably described by DFT, the following ones are clearly inadequate. As a consequence the computed β 1 parameter is in fair agreement with the experiment for HOMO ionization, but pretty off for the following ones [50,51] (Figure 6).   All results considered are evaluated at a fixed nuclear geometry. To include nuclear motion effects it is necessary to evaluate transition dipoles at various geometries and integrate over vibrational wavefunctions [40,41]. Only one-dimensional problems have been considered, but it is easy to generalize to multidimensional models currently available. Vibronic coupling effects have been amply observed [40][41][42]. More demanding is the simulation of pump-probe femtosecond experiments via time resolved photoelectron spectra (TRPES) that require evaluation of photoionization observables for hundreds or more nuclear configurations. For this reason up to now only rather crude approximations, like plane or coulomb waves, have been considered, and computational efficiency becomes of paramount importance. A recent simulation of pyrazine employing a surface-hopping method and the Tiresia code has been presented [52].
A comparison of experimental and calculated time resolved spectra is presented in Figure 7. The agreement is quite satisfactory, some deficiency at the lowest KE is probably due to inaccuracy of the cross section close to threshold. Up to now laser ionization has been employed almost exclusively in this field, and photon energy available is quite low. It is expected that such experimental limitations will be soon overcome. In this respect the simulation may even suggest an optimal photon range for maximum sensitivity to the dynamical aspects investigated. Moreover, additional information may be gained if molecular orientation is not random, but prealigned with an additional laser pulse, so that a richer angular distribution is achieved.

Strong Field and Ultrafast Processes
The program has been used rather extensively to describe and analyze wavepackets in strong field ionization and attosecond pump-probe experiments [53,54] (see also the program by Ruberti [32] based on the same B-spline basis). One of the first applications was to study the angular dependence of strong field ionization probability, and to highlight the role played by ionization of deeper levels besides HOMO [55]. Some results of a study of water [56] are reported in Figure 8. It illustrates the dependence of ionization yield on the laser intensity (at λ = 800 nm). The inset shows how the simple MO-ADK approach that was often used, based on the exponential dependence on ionization energy, overestimates the HOMO contribution except at the highest intensities. Moreover, the different angular dependence of the ionization yield on the angle with the laser polarization, for different orbitals, may even switch the dominant contribution. These results are obtained at the DFT level, and it has been recently shown that important modifications are induced by interchannel coupling [57].
The full energy resolved photoelectron spectrum relative to HOMO ionization, with laser polarization along Z is reported in Figure 9 for different laser intensities. It is interesting to see as, starting with a regular series of peaks corresponding to individual photon absorption, the spectrum becomes chaotic at the highest intensities.  Finally, we present a calculation of HHG emission in CO 2 at different field intensities ( Figure 10). It was experimentally found that a minimum in the HHG emission appears, shifting to higher harmonics as intensity is increased. The DFT results reproduce the trend, although they underestimate the harmonic order of the minimum. A more refined calculation, which included interchannel coupling at essentially the CIS level, is able to obtain complete agreement [58]. It is clear that a lot of information and insight can be obtained already at the DFT level, as might be expected from the great success of oneelectron models to explain most of the strong field phenomena. For quantitative agreement, however, some level of correlation has to be included, even at the lowest level, like CIS of linear response TDDFT. This is important since higher-level correlated approaches become often prohibitively expensive in larger systems.

Conclusions
The structure of the Tiresia code for describing continuum states as well as electronic wavepackets has been presented. The multicenter (LCAO) B-spline basis set has been illustrated in detail. It affords: • A dense set within a finite range (a sphere) of arbitrary length. It can approach completeness and therefore converge to the required solutions. • A complete control of the overlap matrix, hence numerical linear independence and stability. • Accurate solutions of homogeneous and inhomogeneous equations within the range, with proper boundary conditions, which are easily implemented.
In particular for photoionization studies, it affords convergent solutions of the multichannel problem, at selected energies from threshold up to very high KE (10 keV have been reached) and arbitrary energy resolution, for ionizations from valence to deep core, and complex molecular systems. In the present static-DFT and TDDFT formulations, as well as Dyson-DFT and Dyson-TDDFT, it provides a rather accurate description of a vast phenomenology. The algorithms implemented are simple, mostly based on linear algebra, and easily parallelized, leading to efficient computation that makes the treatment of complex systems, and calculations of many geometries easily affordable. Important steps towards the implementation of fully ab initio close coupling formulations have been performed, and further work is in progress. The same basis and hamiltonians allow calculation of multiphoton and strong field processes, in particular time propagation of wavepackets in external electromagnetic fields, and their analysis by projection on free field states of the same hamiltonian. A number of examples from past results and current work are provided for illustration.