Abdallah Ammar

Arnaud Leclerc

Lorenzo Ugo

L U Ancarani

Fitting continuum wavefunctions with complex Gaussians: Computation of ionization cross sections

Keywords: Continuum wavefunctions, Real Gaussians, Complex Gaussians, Non-linear optimization, Ionization

HAL is

Introduction

As early as 1950 Boys [1] emphasized the importance of Gaussians to simplify the calculation of multicenter integrals arising in the study of polyatomic molecular electronic states. One key property is the notorious Gaussian product theorem. For an overview of the use of Gaussian sets in molecular calculations see, e.g., reference [2] . The optimal set depends very much on the studied problem, e.g., system energy minimization, calculation of resonance states,. . . . Usually, to avoid a full nonlinear optimization, one may consider a particular choice of exponents such as eventempered sets [3,4] , well-tempered sets [5] , random-tempered sets [6] or polynomial expansions [7] . In the 60s many studies have been dedicated to the representation of bound states by a real Gaussian combination [8][9][10][11][12][13] . Despite their too fast decay at large distances, Gaussians manage to rather well reproduce molecular electronic orbitals over the physically relevant radial regions [14] . Another drawback is their mathematical inability to reproduce the correct cusp near the origin, an issue that is generally addressed numerically using a combination of Gaussians with large exponents [8] .

While Gaussians are pretty good for reproducing low-lying bound states, representing on large radial domains multinode states (highly excited or Rydberg states), or highly oscillating continuum states by nodeless Gaussians is numerically very challenging. In spite of the numerous potential applications (such as photoionization [15][16][17][18] , high harmonic generation [19][20][21] , ionization by particle impact [22] ) only few exploratory studies have been made [23][24][25][26] . Kaufmann et al [23] optimized a set of Gaussians to represent continuum functions with a reasonable accuracy up to 0.2 a.u. by the diagonalization of the attractive Coulomb Hamiltonian represented in a finite set of Gaussian-type functions. Nestmann and Peyerimhoff [24] proposed a least square approach to represent Bessel functions with real Gaussians. Faure et al [25] applied the same approach to fit Bessel and Coulomb functions. However, the fitting error for large energies (typically ≈ 3 Rydbergs) is visible to the naked eye. Fiori and Miraglia [26] optimized Gaussians to fit distortion functions by first removing the dominant fast oscillatory plane wave from the Coulomb function. An alternative approach to reproduce the oscillating behavior of continuum states is to add to Gaussian sets some supplementary functions, like B-splines [15] . We believe that expansions in complex Gaussians, that it to say Gaussians with complex exponents, offer an alternative and possibly more suitable way of representing continuum states.

Complex Gaussians have already been proven to be successful in resonance stabilization calculations, where they arise from a complex scaling transformation [27][28][29][30] . Such functions possess an intrinsically oscillatory behavior that may facilitate the task of representing continuum states over larger radial domains. Previous investigation in this line, include the work of Matsuzaki et al who used complex Slater orbitals [31] and complex Gaussians [32] in photoionization calculations. To represent the outgoing Coulomb function with complex Gaussians, they employed Pronys method used earlier by Huzinaga [8] . The main drawback of this method is that it cannot cover all the necessary radial ranges because the radial grid must satisfy a specific square root distribution. This is why in [32] a clear difference is visible between the exact continuum function and its complex fitting in the range 0 < r < 5 a.u.. Matsuzaki et al also optimized complex Gaussian sets to fit complex Slater functions [33] . The first aim of the present study is to explore further the capacity of complex Gaussians sets to represent continuum states by proposing and implementing a full nonlinear optimization method. In order to test its limitations and range of applicability we take, as an illustration, a set of regular radial Coulomb functions with different energies; the same can be done with any one-electron continuum function as, for example, a numerically generated distorted wave. We also highlight situations where real Gaussians become insufficiently accurate, and point out the advantage of using complex Gaussians instead.

The second purpose of this manuscript is to demonstrate the advantage of such Gaussian representation when evaluating, for example, ionization cross sections. The key idea is to propose an all-Gaussian approach which allows one to analytically evaluate all necessary integrals. It is very well known that using Gaussians to represent bound states renders most of the necessary integrals analytical, a feature that is especially important in the molecular multicenter case by applying the Gaussian product theorem. One question naturally arises: is it possible to exploit such analytical advantages also when continuum states are involved? If they are represented by standard real Gaussians, the answer is obviously positive. What about complex Gaussians? Is the Gaussian product theorem still applicable? It can be easily verified that when complex exponents are present, the product of two Gaussians with different centers will lead to a Gaussian centered at some complex position. Clearly the latter has no physical sense, but it is not a mathematical obstacle for performing multicenter integrals as shown by Kuang and Lin [34,35] . While our long-term goal is to deal with molecular systems, as a first step we consider here the one-center, atomic, case to put our all-Gaussian proposal on solid grounds. To this effect we consider some one-electron matrix elements ψ - k e (r)| O|φ i (r) corresponding to a transition from an initial state φ i (r) to a continuum state ψ - k e (r) with ejected electron's momentum k e ; the bound-continuum transition occurs via an ionization process represented by some operator O. To envisage an all-Gaussian integration approach, the main difficulty stands in a numerically robust representation of the continuum wavefunction by a set of Gaussian functions, and this up to a sufficiently large radial distance. We shall show that thereafter all integrations can be performed analytically with either real or complex sets. As indicated above, we consider here the one-center case and take hydrogen as a benchmark since its ionization cross sections are known exactly for the photon impact case and also, within the first Born approximation, for the electron impact case.

In section 2 we present the algorithm that performs the fitting. The numerical method, based on the approach of Nestmann and Peyerimhoff [24] , has been improved here on several aspects, by using a better optimization method and extending it to deal with complex Gaussians. We compare different fitting options with real or complex Gaussians and we point out the advantage of the latter.

In section 3 we illustrate the approach in two benchmark applications, the ionization of hydrogen by impact of either an electron or a photon. In both cases we compare the exact cross sections with those calculated with real or complex Gaussian fits of Coulomb continuum states. A brief conclusion is presented in section 4. Atomic units are used unless indicated otherwise.

Fitting with complex Gaussians

We wish to develop an efficient approach that fits complex Gaussians to represent a set of continuum functions arising, for example, in ionization calculations. We start with a detailed comparison between using real and complex Gaussians to highlight the potential benefits of the latter.

Fitting strategy

We aim to approximate a set of arbitrary functions f η (r), η = 1, . . . , η max by a linear combination of N Gaussians:

f η (r) ≈ f G η (r) = N ∑ i=1 [c i ] η exp(-α i r 2 ). (1)
To do so, in the case of real Gaussians, Nestmann and Peyerimhoff [24] proposed a least square approach which consists in minimizing, on some radial grid {r κ } κ=1,...,κ max , the function:

Ξ(α 1 , . . . , α N ) = ∑ η ∑ κ f η (r κ ) -f G η (r κ ) 2 ∑ κ ( f η (r κ )) 2 + D(α 1 , . . . , α N ). (2)
The Ξ function depends on N nonlinear parameters, the exponents {α i } i=1,...,N and η max × N linear parameters, the expansion coefficients {[c i ] η } i=1,...,N,η=1,...,η max . In [24] the standard Powell method [36] is used to optimize the exponents, while the linear coefficients are optimized by a standard least square method. Iterations are performed to alternate those two optimizations: after each variation of the exponents the coefficients are updated using least squares, and the process is repeated until convergence to a local minimum is reached. In eq. ( 2) a penalty function is added to avoid the convergence of two exponents to the same value. It is defined as:

D(α 1 , . . . , α N ) = N ∑ i=2 i-1 ∑ j=1 exp -g α i α j - α j α i , (3)
where g is a fixed parameter (generally g ≈ r κ max ).

Here, we generalize the approach of Nestmann and Peyerimhoff [24] for complex exponents

α i = ℜ(α i ) + iℑ(α i ), with ℜ(α i ) > 0.
The optimization function Ξ becomes:

Ξ(ℜ(α 1 ), . . . , ℜ(α N ), ℑ(α 1 ), . . . , ℑ(α N )) = ∑ η ∑ κ | f η (r κ ) -f G η (r κ )| 2 ∑ κ | f η (r κ )| 2 + D(ℜ(α 1 ), . . . , ℜ(α N )), (4)
and now depends on 2N non-linear real parameters {ℜ(α i ), ℑ(α i )} i=1,...,N so that Ξ is seen as a map from R 2N to R. The penalty function is the same as in eq. ( 3), applied only to the real part of the exponents. eq. ( 2) is a particular case of eq. ( 4) when the exponents {α i } and the coefficients {c i } are real. In order to minimize the fitting error Ξ, we choose to optimize the exponents {α i } by using the Bound Optimization BY Quadratic Approximation (BOBYQA) [37] , still alternating with a least square optimization of the coefficients {c i }. Both Powell and BOBYQA are gradient free methods and attempt to find a local minimum. Since Ξ has many local minima, the aim of the numerical optimization is to find a local minimum that gives a reasonable fitting accuracy. A critical issue in both methods is the choice of the initial values of the exponents. For BOBYQA, in addition to this, we have to fix the initial (∆ i ) and final (∆ f ) trust region radii where Ξ is approximated to a quadratic model. The optimization is stopped when the Euclidean dimension of the step is less or equal to ∆ f . On the other hand, the optimization with Powell is stopped when no further improvement is obtained after varying the exponents. A supplementary condition to stop the optimization may be the value of Ξ or the CPU time. The main difference between these two methods is that in the case of the standard Powell algorithm, we first determine the search directions and then find the optimal step along those directions, whereas using BOBYQA, we first set the step (by choosing the trust region) and then the directions are found in order to improve the quadratic model or minimize the objective function Ξ. For more details about the algorithms we refer the reader to Refs. [36] and [37] .

As an illustration, we consider a set of 6 regular Coulomb functions E : {F 1 (r)/k 1 , . . . , F 6 (r)/k 6 }, defined as [START_REF] Abramowitz | Handbook of Mathematical Functions[END_REF] :

F η (r) = F l,k η (r) = (2k η r) l+1 e πz 2k η Γ l + 1 -iz k η 2Γ (2l + 2) e ik η r 1 F 1 l + 1 - iz k η , 2l + 2; -2ik η r , (5)
with wavenumbers k η = 0.5 + 0.25(η -1) a.u., and angular momentum number l = 1. 1 F 1 is the Kummer confluent hypergeometric function. The real valued functions {F η (r)} are the exact solutions of the one particle Schrodinger equation with Coulomb potential -z/r and are strongly oscillating for large positive energies. For charge z = 1 we have the hydrogen continuum states while if we set z = 0 we obtain the spherical Bessel functions. The set E serves here as a test to compare real and complex Gaussian fittings, and will be used also in the cross section calculations of section 3.2. We apply the strategy presented above up to r κ max = 25 a.u. with a radial step 0.025 a.u. to fit E with either N = 30 real Gaussians (30 nonlinear real parameters to reproduce the real functions F η (r)) or N = 30 complex Gaussians (60 nonlinear real parameters to reproduce both the real functions F η (r) and the imaginary part which is 0 here), and we set g = 27. The number N is to be chosen sufficiently large as to reproduce the regular Coulomb functions in the considered range of energy and within the fitting box B. After several convergence tests on N by inspection of the reached Ξ opt , we found that N = 25 complex Gaussians could be judged as sufficient. In the physical application presented in section 3.2, however, integrations go up to 25 a.u. and in order to reduce cross section errors, we chose N = 30. For the sake of comparison, the same N = 30 is taken here also for the real Gaussian representation. Should one consider higher energies and/or larger radial domains, a convergence study should be envisaged possibly requiring a larger N.

It is worth emphasizing that once the optimal set of N exponents {α i } will be found, they can be employed to represent with a reasonable accuracy any other function F η (r) within the considered energy range k η ∈ [0.5; 1.75]. Indeed, simply performing a linear least squares method will provide the corresponding optimal coefficients [c i ] η .

Limitations of real Gaussians

In this subsection we focus on the use of real Gaussians. We first show the efficiency of the quadratic method BOBYQA [37] and then highlight the limitations of real Gaussians in representing continuum functions.

Comparison between BOBYQA and Powell

We wish all the initial exponents α i to increase slowly and consistently within an interval α 1 = a and α N = b. From our numerical experience we found that the distribution

α i+1 α i = α N α 1 1 N-1 leads
to satisfactory results and is obtained by picking up the initial exponents as:

ln(α i ) = 1 N -1 [(N -i) ln(a) + (i -1) ln(b)] . ( 6
)
The value of a should be chosen small enough to reach the end of the fitting box: e -ar 2 max ∼ 1.

For the optimization with BOBYQA, we set a = 10 1.

From our numerical experience, BOBYQA is faster and more efficient in the present context. This may not be true for all optimizations but generally BOBYQA turns out to be at least as efficient as Powell method or better especially for a large number of Gaussians.

Deviations at large distances

Looking at the fitted functions (not shown), both optimizations lead to a fitting quality which is very good inside the fitting box B = [0, 25]. However, the error increases very quickly for r > 25.

This is due to the fact that the price to pay for a good optimization within B is the presence of small exponents α i ∼ 10 -4 and very large associated coefficients c i ∼ 10 10 . The presence of such diffused Gaussians with very important amplitudes at large distance may not be a serious problem in the calculation of matrix elements; indeed, the integration over continuum functions is usually accompanied by a decreasing radial exponential (factor coming from the Hermitian product with bound states) and the fitting error beyond some physical distance will not affect the numerical calculation. This is why we choose to examine the effective functions

F e f f η (β , r) = F η (r)e -β r r 2 (7)
instead of F η (r) itself. The exponential corresponds to a bound decreasing factor where β is a positive number and r 2 comes from the integration volume element in spherical coordinates. The deviation at large distances between the fitting and the original function F e f f η is obviously very sensitive to the value of β .

As an example, we examine in Figure 1 the fitting of F e f f 2 (β , r) for two cases: β = 1 and β = 0.5, corresponding for example to the hydrogen 1s and 2s exponents, respectively. We clearly see that these fittings will cause trouble in the case β = 0.5 since the deviation after r = 25 will jeopardize radial integrals involving F e f f 2 (β , r). However when β = 1 the bound state cancels this deviation. For higher energies (η = 3, . . . , 6) the same problem arises (not shown). The faster the functions oscillate the more important the deviation. For η = 1, on the other hand, the decreasing term e -β r cancels this deviation for both β = 0.5 and β = 1.

Using reduced bounds (RB) for the search of the exponents

We try in this section to soften the errors coming from the diffused part of the fitting. To do so with BOBYQA, the optimization lower bound is modified to α min = 0.01, so that α i does not fall below 0.01. The optimal set of exponents in this reduced bound (RB) case is shown in the fourth column of Table 1. They are overall of the same order of magnitude as those obtained without constraint but the lower bound forbids the smallest, possibly troublesome, exponents. Figure 2 shows the resulting improvement at large distances of F e f f 2 (RB BOBYQA) in the case β = 0.5.

While the overall fitting is better, the accuracy is slightly worse in the fitting box B because imposing a constraint on the lower bound of α i reduces the overall flexibility. Even in this RB approach the coefficients {c i } η remain large. For example, the second column of Table 2 shows the magnitude of the coefficients [c i ] 5 for η = 5 (k 5 = 1.5). All coefficients are larger than 10 6 and most coefficients are of order 10 10 or more. When the continuum functions are substituted by the Gaussian combinations, this ill-conditioning generates a numerical error that is related to the limited machine precision. = F 2 (r)e -0.5r r 2 with its fitting using real Gaussians and BOBYQA without constraint on the exponents (rG BOBYQA) or after reducing the bounds (rG BOBYQA RB). The absolute error on the fitting is plotted in the bottom panel.

Complex Gaussians

In section 2.2 we have shown that representing oscillating functions with real Gaussians requires very small exponents {α i } and implies very large coefficients {c i }. This causes a large deviation out of the fitting box that may turn out to be troublesome in a given application. A partial solution could be to reduce the exponents bounds but the presence of very large values of {c i } does not disappear. We will now explore the ability of complex Gaussians to soften these problems.

In order to optimize the set E of functions (5) defined in section 2.1, we pick the initial complex exponents as:

     ln(ℜ(α i )) = 1 N -1 [(N -i) ln(a) + (i -1) ln(b)] ℑ(α i ) = 0 (8)
with a = 10 -4 , b = 100, and we fix the following research bounds:

   10 -4 ≤ ℜ(α i ) ≤ 1000 -0.1 ≤ ℑ(α i ) ≤ 0.1. (9)
The trust regions are ∆ i = 0.01 and ∆ f = 10 -6 .

The optimal set of complex exponents obtained for the set E are shown in the fifth column of Table 1. We recall that the complex Gaussian expansion optimizes a set of real-valued radial Coulomb functions. Exponents do not appear in complex conjugate pairs, thus necessarily requiring complex coefficients to build up a real function. We verified that the imaginary part resulting from the complex combinations of the optimal complex Gaussians is indeed negligible. Concerning the real part, Figure 3 shows two functions of the set, F 2 (r) and F 6 (r), with their fitting and the corresponding absolute errors, using: (i) 30 real Gaussians and BOBYQA with constraints on {α i } (rG BOBYQA RB) or (ii) 30 complex Gaussians and BOBYQA (cG BOBYQA). It also shows the corresponding effective functions F 7)) with β = 0.5 are compared to their cG BOBYQA representation, and this on a radial scale extended to 50 a.u..

rG

Illustrative applications to ionization problems

For illustration purposes, we consider hereafter a one-electron description. Computing ionization cross sections involves the calculation of transition matrix elements T ik e = ψ - k e (r)| O|φ i (r) where φ i (r) represents the initial (bound) wavefunction and ψ - k e (r) represents the final (continuum) wavefunctions of the ejected electron (with momentum k e ). In order to keep the present investigation free of extra numerical uncertainties and easily reproducible, we choose as continuum state the analytical Coulomb function

ψ - k e (r) = N(a)
e ik e r (2π)

3 2 1 F 1 (-ia, 1; -i(k e r + k e r)) , (10)
where N(a) = e πa 2 Γ(1 + ia) with the Sommerfeld parameter a = z/k e and z the charge seen by the ejected electron. O is the transition operator that connects the initial to final states: 4π q 2 e iq•r in the case of particle impact (q is the momentum transfer vector) and -ε • r for photoionization in length gauge (ε is the polarization vector). In what follows, we will show that if the radial parts of both ψ - k e and φ i are expanded in Gaussians, the calculation of the transition matrix elements becomes analytical for both processes. As mentioned in the introduction the ultimate goal is to implement such an all-Gaussian approach to treat scattering from polyatomic molecules. Here, in order to illustrate the feasibility and the numerical robustness, we consider first an atomic case with an initial wavefunction given by:

φ i (r) = R n i l i (r)Y m i l i (r), (11)
where n i ,l i ,m i are the usual quantum numbers. For the numerical illustration, we shall take as benchmark the hydrogen atom (z = 1) for which exact cross sections are available and serve as a solid benchmark.

Hydrogen ionization by electron impact

We consider the ionization of a hydrogen atom by electron impact: e -+H → H + +2e -. In the first Born approximation the colliding electron is described by a plane wave before (momentum k i ) and after the collision (momentum k s ), while the wavefunction of the ejected electron is the Coulomb function ψ - k e (r) of eq. ( 10). The cross section calculation involves the transition matrix element

T ik e = 4π q 2 F ik e (q), (12)
where q = k i -k s is the momentum transfer vector and

F ik e (q) = ψ - k e |e iqr |φ i (13)
is the atomic form factor with φ i the initial wavefunction (11).

The standard way to separate angular and radial variables is to use a partial wave expansion of the whole continuum wavefunction (10) over the spherical harmonics Y m * l (r)Y m l ( ke ):

ψ - k e (r) = 2 π ∑ l,m i l e iδ l F l,k e (r) kr Y m * l (r)Y m l ( ke ), (14)
where F l,k e (r) is the regular radial Coulomb function ( 5) and δ l = arg Γ(l + 1 + zi k e ) is the Coulomb phase shift. Alternatively, one can extract from the wavefunction (10) the highly oscillating behavior of the plane wave e ik e r and expand in partial waves only the distortion factor represented by the confluent hypergeometric function. The hydrogen continuum state can then be written as

ψ - k e (r) =
e ik e r (2π)

3 2 ∑ l,m D l,k e (r)Y m * l (r)Y m l ( ke ), (15)
where the radial factors [START_REF] Spielberger | [END_REF] D l,k e (r) = 4πe

π 2k e Γ 1 + i k e Γ l -i k e (-i) l Γ -i k e Γ(l + 1)(2l + 1)!! (k e r) l 1 F 1 l - i k e , 2l + 2; -2ik e r (16)
are complex functions. Fiori and Miraglia [26] proposed to follow this second path and found that taking l max = 8 provides a sufficient accuracy in the calculation of hydrogen ionization by proton impact in a given energy range. If both D l,k e (r) and the radial part of φ i are represented by Gaussian combinations, then the transition matrix integral becomes analytical.

We should first shed light on the behavior of the radial D l,k e (r) functions which, supposedly without strong oscillations [26] , should be easier to represent through a Gaussian set. Actually they are smooth only for low l. Functions with large values of l manifest non negligible oscillations as illustrated by Figure 4 which shows the behavior of D 1,k e (r) and D 7,k e (r) for k e = 1.75. While the small l case can be easily represented with Gaussians, the oscillations of the large l case result to be very difficult to reproduce. This l dependence is even more pronounced for larger k e values. On the other hand, it is known that the first partial terms are generally those contributing the most to the continuum wavefunction. Consequently, the efficiency of the Gaussian fitting approach depends on the physical application under consideration and on the number of partial waves needed to achieve an overall acceptable accuracy.

To compare real and complex Gaussians, we use the strategy presented in section 2.1 for 9 different sets of distortion functions, D l : {D l,k e1 , D l,k e2 , . . . , D l,k e5 } for l = 0, 1, . . . , 8. In each set we consider 5 functions corresponding to k eη = 0.25 + 0.75(η -1). The fitting is performed up to r = 20 and we set the parameter g = 18 in eq. ( 3). Concerning the fitting with real Gaussians, for a given l we optimize {ℜ D l,k e1 , ℜ D l,k e2 , . . . , ℜ D l,k e5 } using 20 real Gaussians and {ℑ D l,k e1 , ℑ D l,k e2 , . . . , ℑ D l,k e5 } with 20 other real Gaussians. On the other hand, in the complex case we use 20 complex Gaussians to fit the set D l . In order to somehow take into account the r l behavior of D l,k e (r) close to r = 0, we add a r γ l factor in the Gaussian expansion,

D l,k e (r) = r γ l ∑ s [c s ] l,k e exp -[α s ] l r 2 ( 17
)
where γ 0 = 0 and γ l = 1 for l > 0. We use γ l = 1 instead of γ l = l as a compromise that gives the vanishing behavior close to r = 0 for l > 0 while avoiding error amplifications at large radial distances for large values of l. Figure 5 shows the hydrogen continuum state at two energies k e2 = 1 and k e4 = 2.5, calculated using eq. ( 15) up to l max = 8, for an angle (r, ke ) = 0, and with D l,k e (r) represented by either real or complex Gaussians using eq. ( 17). The corresponding error on the fitting is also shown. For small energies the errors corresponding to the complex fitting are smaller. For large energies this advantage is lost since the partial waves with high order l 3 can not be reproduced accurately either by real or by complex Gaussians. Nevertheless the continuum wavefunction is overall well reproduced with both options because the first few termswhich contribute the most to ψ - k eare sufficiently well fitted. In summary, with expansion (15) the 15) with l up to l max = 8 and comparison with that calculated by either real Gaussians (rG) or complex Gaussians (cG), for k e2 = 1.00 a.u. (upper panels) and k e4 = 2.50 a.u. (bottom panels), and (r, ke ) = 0 in both cases. The corresponding absolute error on the fitting is shown in the bottom sub-panels. fitting difficulties due to oscillations are not completely removed as suggested in [26] and we should remain aware of this weakness depending on the application. In the numerical illustration presented hereafter, the difficulty appears only in partial wave terms that do not contribute substantially and therefore do not affect the overall cross section calculation.

Let us now turn to the calculation of the atomic form factor (13) employing Gaussian expansion (17). The bound radial function of eq. ( 11) is also represented with either real or complex Gaussians:

R n i l i (r) = ∑ t b t e -β t r 2 . ( 18
)
Using expansions ( 15), ( 17) and ( 18), the form factor becomes:

F ik e (q) = 1 (2π) 3 2 ∑ l,m Y m * l ( ke ) ∑ s,t [c s ] * l,k e b t dr r γ l e iQr e -([α s ] * l +β t )r 2 Y m i l i (r)Y m l (r), (19)
where Q = qk e is the momentum of the hydrogen ion after the collision. To separate angular and radial variables we use the Rayleigh expansion:

e iQr = 4π ∑ λ ,µ i λ j λ (Qr)Y µ * λ ( Q)Y µ λ (r), (20)
where j λ are the spherical Bessel functions. Therefore, in eq. ( 19) we have an integral over 3 spherical harmonics multiplied by the following radial integral,

I rad = ∞ 0 dr r 2+γ l e -([α s ] * l +β t )r 2 j λ (Qr) (21)
which can be calculated analytically (eq. 6.6.31 of [START_REF] Gradshteyn | Table of Integrals, Series, and Products[END_REF] ):

I rad = √ π 4 Q 2 λ Γ λ +γ l +3 2 Γ λ + 3 2 [α s ] * l + β t - λ +γ l +3 2 1 F 1 λ + γ l + 3 2 , λ + 3 2 ; -Q 2 4 [α s ] * l + β t . (22)
Finally the form factor (19) can be written as

F ik e (q) = 2l i + 1 32π ∑ l √ 2l + 1 l+l i ∑ λ =|l-l i | iQ 2 λ Γ λ +γ l +3 2 Γ λ + 3 2 2λ + 1S rad l,λ S ang l,λ , (23)
where

S rad l,λ = ∑ s,t [c s ] * l,k e b t [α s ] * l + β t - λ +γ l +3 2 1 F 1 λ + γ l + 3 2 , λ + 3 2 ; -Q 2 4 [α s ] * l + β t (24)
and

S ang l,λ =   l l i λ 0 0 0   l ∑ m=-l   l l i λ m m i -m -m i   Y m * l ( ke )Y -(m+m i ) * λ ( Q) (25)
where

  j 1 j 2 j 3 m 1 m 2 m 3 
 denote the 3 j Wigner coefficients. Assuming that the hydrogen target is in its ground state 1s, it is possible to compare to the exact atomic factor [START_REF] Mcdowell | Introduction to the Theory of Ion-Atom collisions[END_REF] :

F 1sk e (q) = 2 √ 2e π 2k e Γ 1 -i k e π(1 + Q 2 ) 2 U i k e 1 - i k e + k e + i k e U , (26)
where

U = q 2 -(k e + i) 2 1 + Q 2 . (27)
The triple differential cross section (TDCS) is defined as:

d 3 σ dΩ s dΩ e dE e = 1 4π 2 k s k e k i |T 1sk e (q)| 2 , (28)
with E e = k 2 e /2 the energy of the ejected electron. Ω e and Ω s are the solid angles for the ejected and the scattered electron, respectively. We have calculated the TDCS for coplanar geometry at

E i = k 2
i /2 = 250 eV and scattering angle of 3 • for k e = 0.25, 1.00 and 1.75, with either real or complex Gaussian expansions. Both results, shown in Figure 6, are very close to the exact values given by eq. ( 26), and display a similar accuracy. As discussed before, this can be related to the fact that the first partial functions D l,k e (r) (making the largest contribution to the TDCS) are relatively smooth and can be well reproduced by both real and complex Gaussians with an equal precision.

The very good reproduction of the exact cross section is also due to the presence of a fast decaying initial state (e -r ) so that any error from outside the fitting box will not affect the correct evaluation of the matrix element (see discussion in section 2).

Other scattering quantities can be tested. For example, by integrating over the solid angle of 26), and dotted/dashed curves represent the TDCS calculated with either real Gaussians (rG) or complex Gaussians (cG). The corresponding absolute error is shown in the bottom panels.

the scattered electron, one can define a doubly differential cross section (DDCS)

The integration of the TDCS (28) is to be performed numerically from either the closed form of the form factor (Eq. ( 26) and ( 27)) or its partial expansion (Eq. ( 23), ( 24) and ( 25)) up to l max terms that make use of the Gaussian representation. As an illustration, in figure (7) we show the results as a function of the ejected energy E e = k 2 e /2, for an incident energy E i = 250 eV, a coplanar geometry and the ejected angle fixed at θ e = 0. We clearly observe that with l max = 8 in the continuum state representation, the expected cross section is recovered (with an error of less than 1 percent).

Photoionization

We now consider a hydrogen atom in the initial state φ i illuminated by a photon of energy E γ = k 2 e /2+V ion , where V ion is the energy needed to ionize the target. This photon interacts with the atom leading to an ion and a photo-electron of energy E k e = k 2 e /2. The photoionization cross section is The doubly differential cross section (DDCS) with respect to the energy of the ejected electron, for an incident energy E i = 250 eV, a coplanar geometry and ejected angle θ e = 0. The "exact" result (solid line) is obtained by integrating the TDCS (28) expressed using the closed form (Eq. ( 26) and ( 27)), while the cG results (symbols) are obtained by using the partial wave expansion (Eq. ( 23), ( 24) and ( 25)) and the complex Gaussian representation with l max = 2 or 8. defined as [START_REF] Burke | R-Matrix Theory of Atomic Collisions[END_REF] σ = 4π 2 k e E γ c

1 2l i + 1 ∑ m i dΩ k e |T ik e | 2 , ( 30
)
where c is the speed of light in vacuum and Ω k e is the solid angle of the ejected photo-electron. In the length gauge and the dipole approximation the transition matrix element is given by

T ik e = ψ - k e (r)| -ε • r|φ i (r) , (31)
where ε defines the polarization direction. Here we expand the hydrogen continuum wavefunction (10) in the standard way; expansion (14) presents the advantage that selection rules reduce the summation over l to just l = l i -1 and l = l i + 1. In contrast, when using expansion (15), converged photoionization calculations require many partial terms in order to cover large physical distances. For example we found that for a 1s hydrogen target, about 10 partial terms D l,k e (r) are needed to achieve a relative error of the order 0.01, and one has to consider ≈ 30 partial terms to reach a similar accuracy if the hydrogen is initially in the radially more extended 2s state. With expansion ( 14), only two partial terms appear in both 1s and 2s cases.

The calculation of transition matrix elements (31) involves two integrals:

J ang = Y m * l (r)Y 0 1 (r)Y m i l i (r)dΩ r , (32)
and

J rad l = ∞ 0 F l,k e (r) R n i l i (r) r 2 dr. (33)
The angular integral J ang imposes the selection rules l = l i ± 1 and m = m i . The radial function R n i l i (r) is fitted by real Gaussians as in eq. ( 18), and the functions F l,k e (r) are written as

F l i ±1,k e (r) = ∑ s [c s ] l i ±1,k e exp(-[α s ] l i ±1 r 2 ). (34)
with the different Gaussian parameters optimized to fit the E set, given in Table 1. Using these Gaussian representations, the radial integral becomes ). The hydrogen atom is initially in the 1s state (upper panels), and in the 2s state (bottom panels). Black crosses, green circles and blue squares refer respectively to cross sections computed with real Gaussians with Powell, real Gaussians with BOBYQA by imposing reduced bounds, and complex Gaussians with BOBYQA. On the right panels, the histograms show the corresponding relative errors on the cross section computed with Gaussians.

J rad l i ±1 = √ π 4 ∑ s,t [c s ] l i ±1,k e b t ([α s ] l i ±1 + β t ) -3 2 . ( 35
)
After summation over the magnetic numbers m i , the cross section can be simply expressed as:

σ = 8π E γ 3 (2l i + 1) k e c l i J rad l i -1 2
+ (l i + 1) J rad l i +1 2 (36)
For a hydrogen target initially in a ns state the cross section calculation involves only the l = 1 partial function (5) which is expanded using the different Gaussian sets optimized to fit the E set given in section 2.

Figure 8 shows the photoionization cross section as a function of the photo-electron wavenumber k e for initial states 1s or 2s, calculated with real Gaussians optimized using Powell (rG POW-ELL), real Gaussians optimized using BOBYQA with reduced bounds (rG BOBYQA RB), com-plex Gaussians optimized using BOBYQA (cG BOBYQA) and the exact results given by [START_REF] Harriman | [END_REF] :

σ 1s = 2 5 π 2 e -4 arctan(k e ) k e 3cE 4 γ 1 -e -2 π k e , (37)
and

σ 2s = π 2 8 + 3 E γ e -4 arctan(2k e ) k e 6cE 4 γ 1 -e -2 π k e . ( 38
)
The corresponding relative errors are plotted in a histogram. The results are consistent with the quality of the fitting detailed in section 2. Concerning H(1s), there is no important difference between the different Gaussian basis sets except for the largest energy, k e = 1.75 where the rG POWELL method fails. As explained in section 2, this is due to large deviations outside the fitting box. For H(2s) the calculation with rG POWELL completely fails except for the small energy case k e = 0.50. The calculation with rG BOBYQA RB gives a small relative error at low energies, but fails above k e = 1.25. Only the cG BOBYQA choice achieves a very good accuracy in both 1s and 2s cases for all energies. These photoionisation cross section calculations demonstrate that care must be taken when using Gaussian expansions for continuum states, especially at higher energies and if they appear in matrix elements whose integrands extend on a large radial domain. For these delicate cases, complex Gaussians are clearly superior to their real counterpart.

Conclusions

In this work we applied a full nonlinear optimization to represent oscillating functions with both real and complex Gaussians by using a least square approach in combination with a quadratic approximation method. We have shown that real Gaussians may be sufficient to correctly span the space of functions within the fitting box, however they behave erroneously outside with sometimes very large errors due to diverging coefficients. Real Gaussian representation of continuum radial functions can nevertheless be used in ionization integrals if the decay factor associated with the bound state decreases fast enough to compensate this weakness. This may not be the case for highly excited states. Imposing a constraint on the Gaussian exponents reduces the error magnitude at small energies but it also reduces the fitting accuracy in the box and does not resolve the problem of ill-conditioned coefficients. On the other hand, complex Gaussians have an intrinsically oscillating behavior and are clearly more appropriate to fit bound functions that spread over large distances (highly excited or Rydberg states) and scattering functions.

As an illustrative application, we have used the regular Coulomb wavefunctions fitted by Gaussians in hydrogen ionization problems. The accuracy of cross sections reproduced with complex Gaussians is as good as that obtained with real Gaussians at small energies while it is clearly superior at large energies, where real Gaussians fail completely. The application also allowed us to show the advantage provided by an all-Gaussian approach: the matrix elements can be evaluated analytically once all involved functions, bound and continuum, are expanded in real or complex Gaussians. In the present work we considered pure Coulomb continuum functions but our numerical strategy can be applied in the same manner to any one-center distorted wavefunctions generated analytically or numerically. We have already successfully tested the complex Gaussian representation of positive energy generalized Sturmian functions [44] ; with an adequately chosen asymptotic behavior, the latter can be used as an efficient basis to expand any distorted wave and thus to study collision problems.

Extension to more complicated situations, including many-electron atoms and those molecules that can be described successfully with one-center expansions, do not present any additional technical difficulties. Going beyond scattering from a central potential, electron-electron integrals can be treated within a multipolar approach, and dealing with the ensuing integrals is part of our current investigations. Work is also ongoing to extend our proposal to study ionization processes with molecular targets, and thus to deal with multicenter integrals. For molecules the initial electronic wavefunctions are often already calculated in Gaussian bases. The representation of continuum states by Gaussians, even with complex exponents, should provide a way to simplify some of the multidimensional numerical integrations needed in scattering calculations. The use of the Gaussian product theorem will lead to some meaningless complex centers. These, however, do not lead to mathematical obstacles that cannot be dealt with as indicated in the work of Kuang and Lin [34,35] .

- 6 ,

 6 b = 1. Two research bounds are defined α min = 10 -6 and α max = 10. The initial trust region is ∆ i = 0.01 and the final one is ∆ f = 10 -6 . For Powell optimization there are no constraints on {α i } except being strictly positive, and slightly different initialization parameters are selected: a = 10 -4 and b = 10. The time taken to perform the optimization of the set E with Powell is ≈ 26 times that needed with BOBYQA. The final minimum value of the error (Ξ -D) found with Powell is 0.18 × 10 -3 and 0.02 × 10 -3 with BOBYQA. The optimal sets of real exponents obtained by BOBYQA or Powell are displayed in the second and third columns of Table

Figure 1 :

 1 Figure 1: Effective function F e f f 2 (β , r) = F 2 (r)e -β r r 2 (l = 1 and k 2 = 0.75 a.u.) and its fitting with real Gaussians using Powell (rG POWELL) or BOBYQA (rG BOBYQA) optimizations, for β = 0.5 (left upper panel) and β = 1 (right upper panel). The absolute errors are plotted in the corresponding bottom panels.

Figure 2 :

 2 Figure 2: Upper panel: F e f f 2

η

  (β , r), with no deviation due to the diffused Gaussians outside the box B. One can clearly see that complex Gaussians reproduce the sample Coulomb functions with a better accuracy, and this becomes particularly evident for large energy values. The faster the functions oscillate, the more gainful complex Gaussians become. All expansion coefficients for the F e f f 5 case, shown in the third column of Table2, remain moderate in contrast to the real Gaussian fitting.

Figure 3 :

 3 Figure 3: In panels (a) and (b) the sample Coulomb function F η (r), with l = 1, is plotted for η = 2 and 6 respectively with the fitting by 30 real Gaussians using BOBYQA with reduced bounds (rG BOBYQA RB) or by 30 complex Gaussians (cG BOBYQA). The associated panels show the corresponding absolute errors within the fitting box. In the bottom panels (c) and (d) the effective function F e f f η (β , r) (see eq. (7)) with β = 0.5 are compared to their cG BOBYQA representation, and this on a radial scale extended to 50 a.u..

Figure 4 :

 4 Figure 4: Real and imaginary parts of the D l,k e (r) functions defined by eq. (16) for k e = 1.75 and l = 1 (upper panel) or l = 7 (lower panel).

Figure 5 :

 5 Figure 5: Real (left panels) and imaginary (right panels) part of the hydrogen continuum states ψ -k e (r) calculated using eq. (15) with l up to l max = 8 and comparison with that calculated by either real Gaussians (rG) or complex Gaussians (cG), for k e2 = 1.00 a.u. (upper panels) and k e4 = 2.50 a.u. (bottom panels), and (r, ke ) = 0 in both cases. The corresponding absolute error on the fitting is shown in the bottom sub-panels.

Figure 6 :

 6 Figure6: Upper panels: Triple differential cross section (TDCS) for the ionization of a hydrogen atom in the ground state by a 250 eV electron, in coplanar geometry, with a scattering angle of 3 • and for ejected electron momentum k e = 0.25, 1.00, or 1.75 a.u.: Solid curves represent the exact analytical TDCS computed with eq. (26), and dotted/dashed curves represent the TDCS calculated with either real Gaussians (rG) or complex Gaussians (cG). The corresponding absolute error is shown in the bottom panels.

d 2 σ dE e dΩ e = d 3 σ

 3 dΩ s dΩ e dE e dΩ s .

Figure 7 :

 7 Figure7: The doubly differential cross section (DDCS) with respect to the energy of the ejected electron, for an incident energy E i = 250 eV, a coplanar geometry and ejected angle θ e = 0. The "exact" result (solid line) is obtained by integrating the TDCS (28) expressed using the closed form (Eq. (26) and (27)), while the cG results (symbols) are obtained by using the partial wave expansion (Eq. (23), (24) and (25)) and the complex Gaussian representation with l max = 2 or 8.

Figure 8 :

 8 Figure8: Hydrogen photo-ionization cross section (in Mb) in terms of photo-electron wavenumber (in a.u.). The hydrogen atom is initially in the 1s state (upper panels), and in the 2s state (bottom panels). Black crosses, green circles and blue squares refer respectively to cross sections computed with real Gaussians with Powell, real Gaussians with BOBYQA by imposing reduced bounds, and complex Gaussians with BOBYQA. On the right panels, the histograms show the corresponding relative errors on the cross section computed with Gaussians.

Table 1 :

 1 Optimal exponents obtained after fitting the set E of Coulomb functions defined by eq. (5). rG POWELL means that the standard Powell method is used with real Gaussians, rG BOBYQA means that BOBYQA is used with real Gaussians, rG BOBYQA RB means that BOBYQA is used with reduced bounds of real exponents, and cG BOBYQA means that BOBYQA is used with complex Gaussians and coefficients.

		POWELL rG BOBYQA rG BOBYQA RB	cG BOBYQA
	α 1	0.0000092	0.0000128	0.0122973	0.0001000 + 0.0151526 j
	α 2	0.0000912	0.0000244	0.0145577	0.0001418 -0.0171753 j
	α 3	0.0009435	0.0005575	0.0170202	0.0013875 -0.0021493 j
	α 4	0.0019943	0.0016484	0.0197536	0.0019035 + 0.0216135 j
	α 5	0.0027596	0.0034884	0.0228103	0.0025827 + 0.0298661 j
	α 6	0.0044246	0.0084723	0.0262398	0.0034908 -0.0241050 j
	α 7	0.0057793	0.0111337	0.0301176	0.0046754 -0.0338441 j
	α 8	0.0072420	0.0136592	0.0344988	0.0061876 + 0.0378398 j
	α 9	0.0089783	0.0163089	0.0394645	0.0156781 + 0.0504851 j
	α 10	0.0114303	0.0194778	0.0451376	0.0196315 -0.0395359 j
	α 11	0.0141988	0.0226689	0.0516205	0.0246020 -0.0352098 j
	α 12	0.0178271	0.0267680	0.0591117	0.0308382 -0.0302620 j
	α 13	0.0221354	0.0311299	0.0677411	0.0386508 -0.0436192 j
	α 14	0.0272471	0.0358749	0.0778144	0.0484105 + 0.0624525 j
	α 15	0.0341612	0.0412564	0.0896070	0.0618203 + 0.0513264 j
	α 16	0.0436030	0.0473579	0.1036178	0.0816826 -0.0552154 j
	α 17	0.0622928	0.0543047	0.1205385	0.1936246 + 0.0097332 j
	α 18	0.1058245	0.0616580	0.1407606	0.2915597 + 0.0024838 j
	α 19	0.1361998	0.0696498	0.1648189	0.4971406 -0.0113366 j
	α 20	0.1729779	0.0779709	0.2204008	0.8532849 + 0.0200975 j
	α 21	0.2157985	0.0871103	0.2668776	1.3797714 + 0.0062079 j
	α 22	0.2755926	0.0970341	0.3164855	2.1958583 -0.0311758 j
	α 23	0.3476818	0.1087089	0.3754007	3.5672935 -0.0087858 j
	α 24	0.4425831	0.1228380	0.4395371	5.7419028 + 0.0012293 j
	α 25	0.6059363	0.1384401	0.5106744	9.2535835 -0.0212019 j
	α 26	0.7983113	0.1591001	0.5897495	14.868739 -0.0267870 j
	α 27	1.3120641	0.3614618	0.6853581	23.954872 -0.0406721 j
	α 28	1.6592168	0.6148974	0.7913713	38.575864 + 0.0168896 j
	α 29	2.0833228	1.0857781	0.9151002	62.088037 -0.0059998 j
	α 30	2.6339627	1.5257161	1.0660043	99.986651 + 0.0164035 j