Abdallah Ammar

Ugo Lorenzo

Arnaud Ancarani

Leclerc

A complex Gaussian

Keywords: complex Gaussian-type orbitals, continuum wavefunctions, photoionization, non-linear optimization, Gaussian integrals

come

INTRODUCTION

Ordinary Gaussian-type functions 1,2 are not especially renowned for their adequacy to the theoretical description of continuum states. Those are involved, for example, in dynamical processes where an electron escapes from the system. In a recent theoretical study about molecular photoionization, 3 one can read that "conventional Gaussian-or Slater-type orbitals are unsuitable for this purpose." Everyone can intuitively agree with that, given that an outgoing electron possesses a wavefunction with an oscillating and non-decreasing behaviour, whereas Gaussian-type functions do not feature any nodes and always go to zero, more or less rapidly depending on the exponent values. Such properties are intrinsically not compatible with the asymptotic behaviour of continuum states. This said, if only there were a chance that some Gaussian set manages to attenuate these apparent incompatibilities, this would be very interesting from a numerical point of view in many applications involving electronic detachment processes. Many Gaussian integrals could then be expressed in closed form and, in the molecular case, multicentric integrals which naturally arise in cross-section calculations could be drastically simplified by using the Gaussian product theorem. [START_REF] Shavitt | Methods in Computational Physics[END_REF][START_REF] Besalú | [END_REF] In this paper we focus on photoionization, which plays an important role in atomic and molecular physics and is strongly related to several other sub-domains of physics, including for example applications in astrophysics, 6,7 plasma science [8][9][10] or biology. 11,12 Photoionization has been studied since the early days of quantum mechanics, and the literature is quite vast. Nowadays, while a theoretical-experimental agreement can be generally observed for atoms, this is not always the case for molecular targets, especially in the near-threshold region where autoionizing resonant structures are difficult to describe, as pointed out for example in refs. 13,14. Such discrepancies are partly related to several theoretical difficulties coming from the multicentric nature of the target. Recent photoelectron spectroscopy experiments on molecules such as those presented in refs. 15-17 provide an ongoing and renewed interest, and a challenge for theoreticians. Indeed, the process provides an indirect tool to test the theoretical description of the target before and after the interaction, and thus our capacity of describing correlation and many-body effects for both bound and unbound electronic states. A number of methods have been proposed in the past, and some are cur-describing photoionization (see, for example, refs. 3,14,18). We are precisely interested in this difficult problem and we consider it as a good test case for a Gaussian representation of continuum states.

For bound states, the efficiency of Gaussian-type orbitals (GTOs) needs no further proof.

For a review about the use of Gaussian basis sets in molecular calculations, see for example ref. 2. However, in scattering processes, as for example photoionization of neutral species or photodetachment of anions, an efficient description of continuum states is also required.

What is meant by efficient is that the continuum wavefunction describing an incident or escaping charged particle, for example an outgoing electron, should be a good compromise between accuracy over a sufficiently large radial domain and numerical efficiency. In terms of accuracy, neither standard GTOs (with real exponents) nor Slater-type orbitals seem very adapted to the task of describing oscillating, non-decreasing continuum wavefunctions, as already stated above. However, as far as numerical efficiency is concerned, the use of Gaussian-type functions may lead to strong savings in terms of computational time. Among these two observations, the first (pessimistic) one about the inadequacy to deal with the continuum seems to have outshined the second (optimistic) one on numerical gain; indeed, in the literature there is a relatively low number of theoretical studies using Gaussian sets for describing continuum states, [19][20][21][22] sometimes within a hybrid basis set also including discrete variable representations 23 or B-splines. 24,25 In a recent work, we have developed an optimization method able to provide sets of complex Gaussians (i.e. Gaussian-type orbitals with a complex exponent) able to reproduce accurately Coulomb waves, and whose scattering applications in the benchmark case of atomic hydrogen demonstrated their reliability. 26 The complex Gaussian-type orbitals (cGTOs) are defined as

G(r) = e -αr 2 with α ∈ C, (1)
where Re(α) > 0 ensures that the functions are square integrable. cGTOs were first introduced in the context of resonance stabilization calculations. [27][28][29][30] More recently cGTOs were used in atomic photoionization calculations [START_REF] Matsuzaki | Highlights in Theoretical Chemistry[END_REF][START_REF] Matsuzaki | [END_REF] with an optimization method based on the Prony's algorithm used earlier by Huzinaga. 33 In ref. 26 we have focused our efforts on the numerical method used to perform cGTOs optimization in order to improve their numerical efficiency as a continuum representation. cGTOs oscillate with a spatial frequency increasing with the radial coordinate r, a feature that is not physically satisfactory. On the other hand, in comparison with real GTOs (rGTO), this oscillatory behavior due to the complex exponent leads to more stable expansions of physically sound continuum wavefunctions. We have also shown that optimized cGTO sets can actually be used to calculate accurate ionization cross-sections for atomic hydrogen, under photon or electron impact. It turns out that the necessary matrix elements can be evaluated analytically, making the approach numerically efficient.

Here we wish to adopt a similar strategy by presenting cGTO calculations of molecular photoionization differential cross-sections and asymmetric parameters, using a number of approximations: we use a monocentric description of the targets, considering only small molecules of type AH n which are known to reasonably bear such an approximation; we restrict ourselves to a one-active-electron description; we assume the molecular targets to be randomly oriented; we work in the dipolar approximation with linearly polarized light; the initial bound orbitals of the target are described in terms of Slater-type orbitals with parameters taken from Moccia. [34][35][36] The principal aim is not to develop a sophisticated chemical model or to get as close as possible to the available experimental results, but rather to prove that it is indeed possible to obtain converged results on differential photoionization crosssections for molecules using pure cGTO expansions of the radial continuum wavefunction.

Such a representation allows us to get closed-form expressions for the transition matrix integrals, the numerical evaluation of which becomes very fast. The monocentric (also often named Single Center) approximation is an important restriction, and in a near future we may definitely consider similar calculations with a multicentric description of the target.

The main idea of the present paper is to show that theoretical photoionization cross-sections keep the same level of accuracy within a given model, even when cGTO expansions are used to describe the outgoing electron.

In the rest of the paper, we first summarize the photoionization model and parameters in section 1. We introduce the fitting method that we have used to find optimal cGTO expansions of the (monocentric) continuum functions on a definite domain, adapted to molecular applications. Then, integrals needed for the photoionization cross-section calculations are detailed. Finally, in section 2 come the results on two illustrative examples, the Ne-like molecules NH 3 and H 2 O, followed by some concluding remarks. Unless otherwise stated, we use Hartree atomic units throughout the paper.

1 GAUSSIAN-SLATER MOLECULAR PHOTOIONIZATION MODEL

Photoionization parameters

During the photoionization process, one electron of the target is ejected from an initial bound orbital φ i (r) to a continuum state ψ - ke (r) with wavevector k e (kinetic energy E ke = k 2 e /2 and solid angle k e ). We restrict ourselves to a simple one-active-electron model with electronic coordinate r, all other electrons remaining unaffected. We work in the dipole approximation and the photon of energy E γ is linearly polarized along a direction ˆ . The energy conservation reads

E γ = E ke + V ion , (2)
where V ion is the ionization energy associated with the initial state.

The photoelectron angular distribution (sometimes quoted with the acronym PAD in the literature) is given by the differential cross-section, which should be averaged over all the possible molecular orientations. In a time-independent framework, it is given by: [START_REF] Bethe | Quantum mechanics of one-and two-electron atoms[END_REF][START_REF] Drake | Springer Handbook of Atomic, Molecular, and Optical Physics[END_REF] dσ

d k e = N i 4π 2 k e E γ c d R T ike ( R) 2 (3)
where R stands for the molecular orientation with respect to the laboratory frame (Euler angles), N i denotes the initial number of electrons in the concerned molecular orbital and c is the speed of light. Although not explicitly written, the cross-section (3) is also differential with respect to the photoelectron energy. The key quantities to be calculated are the dipole transition moment elements T ike . In length and velocity gauges, they are respectively defined by

T (L ) ike = ψ - ke |ˆ • r|φ i (4)
or

T (V ) ike = 1 ıE γ ψ - ke |ˆ • p|φ i (5)
with p = -i∇. It is worth reminding that, while calculated cross-sections should be in principle gauge-invariant, using an approximate Hamiltonian or non-exact wavefunctions may lead to substantial differences between length and velocity gauges results. [START_REF] Drake | Springer Handbook of Atomic, Molecular, and Optical Physics[END_REF] In the case of a linearly polarized photon, the differential cross-section for randomly oriented molecules can also be expressed using an energy-dependent asymmetry parameter β, through: [START_REF] Chandra | [END_REF]40 dσ

d k e = σ(k e ) 4π [1 + βP 2 (cos (θ)] (6)
where

σ(k e ) = dσ d k e d k e (7)
is the integrated cross-section (or partial ionization cross section, sometimes quoted with the acronym PICS in the literature), P 2 (x) = 1 2 (3x 2 -1) is the second Legendre polynomial, and θ = (ˆ , k e ) is the scattering angle in the laboratory frame. Compared to the integrated cross-section, the asymmetry parameter β regulates the photoelectron angular distribution and contains thus more detailed information; as a consequence, it is also more sensitive to the quality of the involved wavefunctions.

A typical photoionization calculation thus needs sufficiently accurate wavefunctions for the initial and continuum states to ensure that integrals (4) or (5) yield reasonable results.

Model wavefunctions

The initial molecular orbitals of the target are taken from Moccia (ref. 35 for NH 3 and ref. 36 for H 2 O), who used a self-consistent-field calculation. He proposed a monocentric expansion in terms of N i Slater-type orbitals with real coefficients and real spherical harmonics, centered on the heaviest atom. We transform the data to work with complex spherical harmonics Y m l (r) and complex coefficients C ij , and use the following form,

φ i (r) = N i j=1 C ij r n j -1 e -ζ j r Y m j l j (r). ( 8
)
with normalization φ i |φ i = 1.

For the outgoing electron, we shall use a standard single center partial wave expansion, 41

ψ - ke (r) = 2 π l,m (ı) l e -ıδ l u l,ke (r) k e r Y m l (r)Y m * l ( k e ) (9)
where δ l denotes the phase shift for a given angular momentum l. The radial part of the wavefunction u l,ke (r) satisfies the following equation,

- 1 2
d 2 dr 2 + l(l + 1) 2r 2 + U mol (r) u l,ke (r) = k 2 e 2 u l,ke (r) (10)
for a central potential U mol (r). This one-dimensional differential equation can be numerically solved by using, for example, the RADIAL code. [START_REF] Salvat | [END_REF] As a very first approximation, we assume that the potential felt by the photoelectron is Coulombic, U mol (r) -z/r, with a point charge z = 1, so that the continuum wavefunction is a pure Coulomb wave. In this case the phase shift δ l = arg (Γ(l + 1 + ıη)) with the Sommerfeld parameter η = -z/k e , and the radial functions are the (real) regular Coulomb

functions u l,ke (r) =F l (η, k e r) =(2k e r) l+1 e -πη 2 |Γ (l + 1 + ıη)| 2Γ (2l + 2) e ıker × 1 F 1 (l + 1 + ıη, 2l + 2; -2ık e r) , (11)
where Γ is the Gamma function and 1 F 1 is the Kummer confluent hypergeometric function. [START_REF] Buchholz | The confluent hypergeometric function: with special emphasis on its applications[END_REF] Should one be interested in studying the photodetachment of anions, the remaining core is neutral and the asymptotic charge of U mol (r) is z = 0. In this case, the Coulomb wave should be replaced by a plane wave, and the radial function ( 11) is then a spherical Bessel function.

In reality, the main motivation for using a Coulomb wavefunction is methodological: within such an approximation, the transition matrix integrals are fully analytical, allowing for a rapid check of the convergence.

In subsequent calculations, the photoelectron is assumed to feel a distorted radial potential U mol (r), obtained as the angular average of the anisotropic potential calculated from the direct term in the static exchange approximation. [START_REF] Granados-Castro | Electron Correlation in Molecules ab initio Beyond Gaussian Quantum Chemistry[END_REF][START_REF] Fernández-Menchero | [END_REF][45][START_REF] Castro | Application of Generalized Sturmian Basis Functions to Molecular Systems[END_REF] At short distances, this potential behaves as -Z center /r, where Z center is the approximate charge number at the central nucleus, and features a z = 1 Coulomb asymptotic behavior at large distances (again, in the case of photodetachment of anions, the asymptotic charge would be z = 0). This model potential, which takes somehow into account the charge distribution of the core, can clearly be improved as to include exchange and many-electron effects. As we shall see in section 2, the potential U mol (r) not only provides a reasonable model that allows to capture the correct physics, it serves here also to illustrate that our approach can be applied to any central potential.

The fact that we use a central potential U mol (r) (pure Coulomb or distorted) thus putting aside the multicentric nature of the target, is clearly a strong approximation. It would not make much sense if we were interested for example in photoionization for a fixed molecular orientation. Since in the present study we are looking at AH n molecules with a heavy central atom, and the molecules are randomly oriented as in the experiments, an angular average of the multicentric potential is not an heresy.

Complex Gaussian representation of the radial continuum wavefunction

In this subsection we give an overview of the cGTO representation method introduced in ref.

26, as a generalization of previous works on rGTOs. 20,21 We also provide a table of optimal exponents used later in section 2.

Fitting approach

We consider a set of arbitrary functions f p (r), p = 1, . . . , p max that we would like to approximate by linear combinations of N cGTOs:

f p (r) ≈ f G p (r) = N i=1 [c i ] p exp(-α i r 2 ), (12)
where the exponents are complex, α i = Re(α i ) + iIm(α i ), with Re(α i ) > 0. A minimization is performed for the objective function Ξ defined as

Ξ (Re(α 1 ), . . . , Re(α N ), Im(α 1 ), . . . , Im(α N )) = p κ |f p (r κ ) -f G p (r κ )| 2 κ |f p (r κ )| 2 + D(Re(α 1 ), . . . , Re(α N )), (13)
over some given radial grid {r κ } κ=1,...,κmax . The Ξ function depends on 2N non-linear real parameters, {Re(α i ), Im(α i )} i=1,...,N (the exponents), and N × p max linear parameters {[c i ] p } i=1,...,N,p=1,...,pmax (the expansion coefficients). In eq. ( 13), D is a penalty function whose aim is 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 , ( 14
)
where g is a fixed parameter (in general g ≈ r κmax ).

We start the optimization with some reasonable set of exponents {α i }. Based on our numerical experience, the real parts are picked between two research bounds Re(α 1 ) = a and Re(α N ) = b and follow the distribution

Re(α i+1 ) Re(α i ) = b a 1 N -1 . ( 15
)
This choice ensures that most exponents start with a small real part. The imaginary parts are initially set to zero for all α i .

The fitting error Ξ is then minimized following a two step iterative algorithm: (i) a least square optimization gives an approximation for the coefficients {c i } and (ii) the exponents {α i } are optimized by using the Bound Optimization BY Quadratic Approximation (BOBYQA). [START_REF] Powell | The BOBYQA algorithm for bound constrained optimization without derivatives[END_REF] We iterate over steps (i) and (ii) until some reasonable convergence is reached.

The BOBYQA algorithm requires research bounds that allow us to constrain the domain for the optimal exponents.

Optimal exponents and accuracy

In preparation for molecular applications, basis sets are required for several partial waves, l = 0, 1, 2, 3, 4 being sufficient for our purpose. An optimization has been performed with a set of p max = 6 regular Coulomb functions with z = 1 defined as

F l : {F ν (r) = F l (η, k ν r)} ν=1,...,6 (16)
on a momentum grid k ν = 0.5 + 0.25(ν -1), ν = 1, . . . , 6. Each Coulomb set F l is fitted in a radial box r ∈ [0; 25] using N = 30 cGTOs, i.e., for each value of l there are 60 non-linear real parameters to be optimized along with 30 × 6 linear coefficients.

The initial value for the real parts of the exponents are selected between a = 10 -4 and b = 10 2 and research bounds are chosen as to constrain Re(α i ) ∈ [10 -4 , 10 3 ] and Im(α i ) ∈

[-0.1, 0.1]. Table 1 reports the optimal exponents for l = 0, 1, 2, 3, 4 ordered according to their real parts. The number of digits is necessary as to ensure accuracy. The fitting error for An optimization reaching (Ξ -D) < 10 -4 for one value of l typically requires 4 -5 hours using 8 CPU at 2.66 GHz. Optimizations have been made only for Coulomb functions (and not distorted waves). However we expect that optimal exponent sets can also be used for distorted radial functions in the same energy range without significant loss of accuracy. This will be later confirmed in the result section 2.

1.4 Closed-form expression of the cross-section using Gaussian-Slater integrals

General expressions

We now turn to the explicit calculation of the photoionization differential cross-section and asymmetry parameter. Let us start with a few considerations about the molecular orientation. The incident photon polarization is defined in the laboratory frame, and so is the photoelectron wavevector. Spherical harmonics transform from the molecular frame (from now on, unprimed variables) to the laboratory frame (from now on, primed variables) using Wigner rotation matrices, 48

Y m l r = l m=-l D l mm R Y m l (r) . (17)
For an arbitrary molecular orientation, the dipole operator in (4) can be written as (for convenience we replace the modulus r = r hereafter)

T (L ) ≡ ˆ • r = r cos (θ ) = 2 π 3 rY 0 1 ( r ) = 1 µ=-1 D 1 µ0 R T (L ) µ , (18)
with

T (L ) µ = 2 π 3 rY µ 1 (r). (19)
A similar rotated expression holds for the gradient operator in velocity form (5),

T (V ) ≡ 1 ıE γ ˆ • p = 1 µ=-1 D 1 µ0 R T (V ) µ , (20)
with

T (V ) µ = - 1 E γ ∇ µ , (21)
∇ µ being the spherical tensor components of the gradient operator (see section 5 of ref. 48).

We should also take into account the inverse rotation matrix between the molecular and the laboratory frame since the photoelectron wavevector is measured in the laboratory frame while φ i in eq. ( 8) corresponds to a given molecular orientation. Thus the spherical harmonic in eq. ( 9) is written as

Y m * l ( k e ) = l m=-l D l m m R Y m *
l k e with k e the wavevector direction in the laboratory frame.

Following Chandra [START_REF] Chandra | [END_REF] and in the case of a linearly polarized photon, the explicit formula of (3) after averaging over all possible Euler angles becomes

dσ (G ) k e d k e = L A (G ) L (k e ) Y 0 L k e ( 22
)
where

A (G ) L (k e ) ≡ 8πE γ k e c l 1 ,m 1 ,µ 1 l 2 ,m 2 ,µ 2 (-ı) l 1 -l 2 e ı(δ l 1 -δ l 2 ) M (G ) l 1 ,m 1 ,µ 1 M (G ) l 2 ,m 2 ,µ 2 * × (-1) m 1 -µ 1 (2l 1 + 1)(2l 2 + 1)(2L + 1) 4π   l 1 l 2 L 0 0 0     1 1 L 0 0 0   ×   1 1 L µ 1 -µ 2 -µ 1 + µ 2     l 1 l 2 L -m 1 m 2 µ 1 -µ 2   , (23)
and

M (G ) l 1 ,m 1 ,µ 1 stand for the integrals M (G ) l,m,µ ≡ u l,ke (r) r Y m l (r) * T (G ) µ φ i (r)dr. (24)
In (23), the superscript G refers to the gauge choice, the phase shifts δ l 1 and δ l 2 come from the partial wave expansion (9), and the standard notation is used for Wigner 3j symbols. [START_REF] Edmonds | Angular momentum in quantum mechanics[END_REF] The summation in (22) actually involves only two non-zero terms, L = 0 and L = 2, and is thus equivalent to expression (6) written in terms of the asymmetry parameter. Using the above notations, the total cross-section and asymmetry parameter can be respectively calculated as

σ (G ) (k e ) = A (G ) 0 (k e )Y 0 0 k e + A (G ) 2 (k e )Y 0 2 k e d k e = √ 4π A (G ) 0 (k e ) = 8πE γ 3 k e c l,m,µ M (G ) l,m,µ 2 , (25)
and

β (G ) (k e ) = √ 20π A (G ) 2 (k e ) σ (G ) (k e ) = √ 5
A (G ) 2 (k e ) A (G ) 0 (k e ) , (26)
with

A (G )
2 (k e ) given by ( 23). The numerical calculation of ( 25) and ( 26) is greatly facilitated by the fact that integrals (24) becomes analytical by using cGTOs expansions for the continuum radial function u l,ke (r).

Transition integrals in length or velocity gauges

In the length gauge, we first insert the transition operator (19) in integral (24), and then make use of the initial state expansion (8). The angular part yields the Gaunt coefficients [START_REF] Edmonds | Angular momentum in quantum mechanics[END_REF][START_REF] Gaunt | [END_REF] denoted by

l f m f |L M |l i m i ≡ Y m f * l f (r) Y M L (r) Y m i l i (r) dr.
We thus obtain

M (L ) l,m,µ = 2 π 3 N i j=1 C ij l, m|1, µ|l j , m j I (2) j,l (27)
where the radial integral

I (g) j,l = ∞ 0 (u l,ke (r)) * e -ζ j r r n j +g-1 dr, (28)
is defined here for a given integer g. Note that the sum over j in (27) refers to different contributions in the initial state expansion (8); each value of j is associated with a pair (n j , ζ j ). The Gaunt coefficient imposes the restrictions l = l j ± 1 and m = µ + m j , selecting therefore some specific terms in this expansion.

The formulation is a little bit more complicated in velocity gauge. First, the gradient operator ( 21) is inserted in (24), leading to

M (V ) l,m,µ = - 1 E γ N i j=1 C ij ∞ 0 r 2 dr (u l,ke (r)) * r l, m ∇ µ r n j -1 e -ζ j r l j , m j . (29)
Replacing the gradient formula given in section 5.7 of ref. 48, the angular integral in (29) becomes

l m|∇ µ r n j -1 e -ζ j r |l j m j = Π m µ m j l 1 l j ∆ l l j ∂ r + b l l j r r n j -1 e -ζ j r = Π m µ m j l 1 l j ∆ l l j -ζ j r + n j -1 + b l l j r n j -2 e -ζ j r (30)
with angular coefficients

Π m µ m j l 1 l j = (-1) m   l 1 l j -m µ m j     l 1 l j 0 0 0   -1 , (31)
∆ l l j =          l j + 1 (2l j + 1)(2l j + 3) if l = l j + 1, l j (2l j -1)(2l j + 1) if l = l j -1, (32)
and

b l l j =    -l j if l = l j + 1, l j + 1 if l = l j -1. (33)
Using again the notation (28), the dipole integral becomes

M (V ) l,m,µ = - 1 E γ N i j=1 C ij Π m µ m j l 1 l j ∆ l l j -ζ j I (1)
j,l + n j -1 + b l l j I (0) j,l . (34)
Now let us assume that the continuum radial function, as defined in (10) or (11), is expanded in an optimal cGTO set of the form (12), i.e.,

(u l,ke (r)) * ≈ N s=1 [c s ] l,ke e -[αs] l r 2 . ( 35
)
Then radial integrals (28) can be written in closed form as

I (g) j,l = N s=1 ∞ 0 [c s ] l,ke e -[αs] l r 2 e -ζ j r r n j +g-1 dr = (n j + g -1)! 2 n j +g N s=1 [c s ] l,ke [α s ] n j +g 2 l U n j + g 2 , 1 2 ; ζ 2 j 4 [α s ] l , (36)
where U (•, •; •) is the Tricomi function. [START_REF] Bateman | Higher transcendental functions[END_REF][START_REF] Gradshteyn | Table of integrals, series, and products[END_REF] It can be numerically evaluated, for example, using the Python library mpmath. [START_REF] Johansson | mpmath: a Python library for arbitrary-precision floating-point arithmetic[END_REF] Expression (36) depends on the optimal complex exponents {α s } and expansion coefficients {c s }. Note that if the radial continuum function u l,ke (r) is real, the complex conjugation in (28) has no effect. However, in the case its cGTOs expansion is not perfectly real, we approximate the integral to its real part, I (g) j,l

Re(I (g) j,l ). In the particular case of a pure Coulomb continuum function given by ( 11) with η = -1 ke , it is even possible to compute integrals (28) in closed form without using cGTO expansions of the radial function. Indeed, we have

I (g)Coul j,l = ∞ 0 F l (η, k e r) e -ζ j r r n j +g-1 dr = k l+1 e e π 2ke
2 l Γ l + 1 -ı ke Γ (2l + 2) Γ (n j + g + l + 1) (ζ j -ık e ) n j +g+l+1 × 2 F 1 l + 1 + ıη, n j + g + l + 1, 2l + 2; -2ık e ζ j -ık e (37)
where 2 F 1 is the Gauss hypergeometric function (see formula 13.10.3 of ref. 53). This analytical result can be easily evaluated. [START_REF] Johansson | mpmath: a Python library for arbitrary-precision floating-point arithmetic[END_REF] It can be used as a benchmark to validate the evaluation of ( 36) that makes use of the cGTO representation.

Before presenting our results, we should emphasize that in ref. 54 the photoionization of atomic hydrogen was considered and used as a testbed for the analytical evaluation of matrix elements. Indeed, the cross-section from the exact initial states n i l i m i is known in closed form, and formulae ( 36) and (37) could be tested successfully.

RESULTS

We shall now present partial cross-sections σ(k e ) and asymmetry parameters β(k e ) for NH 3 and H 2 O photoionization. The initial state is described with Moccia's monocentric expansion as explained in subsection 1.2 and the final continuum state is represented with cGTO optimal sets as introduced in subsection 1.3. In each case, assuming that the cGTO set has been previously optimized, the numerical evaluation of closed-form expressions given by eqs. ( 25), ( 26), ( 27) and (34) for one value of σ(k e ) and β(k e ) in the two gauges takes about 2 minutes of computational time using one CPU at 2.67 GHz. For each molecule, we first validate the cGTO approach using internal comparisons for one selected valence orbital, i.e., we compare the results obtained using cGTO representations for the radial continuum function with results given by the original wavefunction. We do this for the outgoing electron described by either a pure Coulomb wave or a distorted wave, and in both velocity and length gauges. Then, we present the distorted wave results for cross-sections and asymmetry parameters for some other orbitals, along with comparisons with experimental and other theoretical results from the literature.

NH 3 photoionization parameters

The ground state electronic configuration of NH 3 is 1a 2 1 2a 2 1 1e 4 3a 2 1 1 A 1 . Fig. 2 shows the partial cross-section and asymmetry parameter for the lowest energy ionization (3a 1 valence orbital), obtained using either optimal cGTO representation for the radial continuum wavefunction, or the original wavefunction. Results are shown for length and velocity gauges, and using either a pure Coulomb wave or a distorted wave as explained in subsection 1.2.

In the Coulomb wave case, we compare the results obtained using cGTOs and Gaussian integrals (eq. ( 36)) with those obtained by direct Coulomb integrals (eq. ( 37)). The curves being almost undistinguishable, the quality of the cGTO representation is clearly sufficient to reproduce the exact results, in both length and velocity gauges. In the case of a distorted radial function, the comparison of cGTO results with those obtained by numerical integration of the original distorted continuum wave, leads to the same conclusion. The results from the two continuum models, Coulomb wave vs distorted wave, are expectedly different, illustrating the strong sensitivity of the results to the physical model used to describe the outgoing electron wavefunction. Besides, the observed gauge dependence should not surprise since both initial and final states are not exact eigenfunctions of the Hamiltonian. The important point here is that, from a methodological point of view, the Gaussian approach is validated in both cases, and in both gauges. In other words, using cGTOs does not jeopardize the quality of photoionization parameters within the limits of the selected model.

To put our results in perspective with respect to the literature, Fig. 3, 4 and 5 show the partial cross-sections and asymmetry parameters for ionization from orbitals having the lowest ionization energies, i.e. from orbitals 3a 1 , 1e and 2a 1 , respectively. Results are compared with other theoretical results, some of them using more sophisticated models, [START_REF] Novikovskiy | [END_REF][START_REF] Stener | [END_REF] and with experimental points. 56,57 Note that experimental results are typically associated with uncer- Results using cGTOs integrals (cG) are compared with results from exact integrals (th) in the Coulomb wave case and with results from numerical integrals (num) in the distorted wave case. Calculations are performed in both length (L) and velocity (V) gauges.

tainties of about 3 -5% (not shown in the figures). We get an overall reasonable agreement for the partial cross-sections, given the simplicity of the monocentric model. The difference between length and velocity gauge results is sometimes quite large but the experimental cross-sections are always framed between the two results. A meaningful comparison with other theoretical cross-sections should be done within the same gauge, i.e. in velocity gauge for ref. 14, while in length gauge for ref. 55. Results for the asymmetry parameter, which is related to the photoelectron angular distribution and thus is a more sensitive quantity, are good for the outer valence orbital 3a 1 but larger deviations are visible in the 1e case. For the 2a 1 orbital, the model yields less satisfactory results. We thus observe two general trends:

(i) our results for outer orbitals give an overall better agreement with experiments and other calculations, for both the cross section and the asymmetry parameter and (ii) cross section results are in general better than the asymmetry results. We think that these variations are mainly due to the simplicity of the static-exchange potential (used in the one-active-electron model for the description of the final state) which is not sufficiently realistic for more internal target orbitals. Moccia's wavefunctions describing the initial bound orbitals could also be suspected to be responsible for such variations. However, this can be ruled out by the fact that the results of ref. 14, using a similar level of approximation (a single-center Hartree-Fock wavefunction) for the initial state but a different model for the final state, differ from ours.

Once again, the aim of the present study is not to defeat previous results from the literature in terms of physical description but rather to prove the feasibility of fast cGTO calculations in the context of molecular photoionization processes. In this respect, the present results are fully satisfactory.

H 2 O photoionization parameters

We now turn to the results for water. The ground state electronic structure of H 2 O is

1a 2 1 2a 2 1 1b 2 2 3a 2 1 1b 2 1 1 A 1 .
An internal comparison between the partial cross-section and asymmetry parameter (for the outer valence orbital 1b 1 ) obtained using the original functions or their cGTO representation, is shown in Fig. 6. As previously observed for ammonia, the agreement between cGTO calculations and results from either theoretical integrals (in the Coulomb wave approximation) or numerical integrals (in the distorted wave case) is very Formulas have been here derived within the framework of the monocentric approximation.

Transition integrals involving cGTOs have been written in closed form which can be easily evaluated at little numerical cost. The preliminary optimization giving optimal exponents for the Gaussian sets within a given energy range is costly but it had only to be performed once for each orbital quantum number. We have provided a table of such exponents that can be used for other applications involving continuum Coulomb or distorted states with energy up to 45 eV; the linear coefficients can be easily obtained through linear optimization. Our complex Gaussian approach has been validated here by studying photoionization of two molecules of type AH n , namely ammonia and water. Practically no difference is seen between a conventional calculation and a closed form approach through the proposed cGTO representation.

A similar conclusion may be reached by studying the photodetachment of anions. In this case, the Coulomb wave with charge z = 1 are to be replaced by plane waves, and the cGTO exponents have to be calculated for a set of spherical Bessel functions of first kind instead of the regular Coulomb functions, 54 i.e. formula (11) with z = 0.

The present investigation focused on a cGTO representation of the continuum states. For the initial bound state φ i (r) we took a single center expansion over Slater-type orbitals. In the case where GTO were used instead, the involved matrix elements are also all analytical.

Actually, in such an all Gaussian approach, the radial integrals turn out to be even simpler because the integrands involves only Gaussian type exponentials and powers of the radial coordinate. The resulting closed form expressions have been presented in ref. 26, and have been numerically validated with an rGTO representation of the exact states of the hydrogen atom 26 and of approximate states of small atoms. [START_REF] Ammar | Représentation des états du continuum par des gaussiennes complexes : application aux processus d'ionisation atomiques et moléculaires[END_REF] We are currently working on the natural extension of the proposed complex Gaussian strategy to represent continuum states in two directions. The first one is to use cGTO to deal with other applications such as (e, 2e) molecular ionization under electron impact, and the study of High Harmonic Generation. [61][62][63][64] The second direction, even more interesting, is related to the well-known Gaussian product theorem leading to drastic simplifications in the calculation of multicentric Gaussian integrals. This well-known property remains valid in the presence of complex exponents 65,66 and this motivates us to develop a generalization of the cGTO transition integrals, considered here in the monocentric approximation, to a more realistic multicentric model. This should be made by taking into account the multicentric nature of the molecular target, in both initial (bound) and final (continuum) states. In a preliminary investigation, [START_REF] Ammar | Représentation des états du continuum par des gaussiennes complexes : application aux processus d'ionisation atomiques et moléculaires[END_REF]67 we have already been able to show that the analytical character of the calculations, which is related to using cGTO representations of the continuum, holds with a multicentric test-case initial state of Gaussian type. Whether the latter is described in spherical or cartesian coordinates, the matrix elements evaluation involves more laborious mathematics. Remains to apply the formulation to photoionization of a realistic molecular model. As a second step, a similar multicentric modification in the final state should be considered. This would ultimately drastically speed up numerical calculation for processes implying electronic continua in larger molecules. as defined in eq. ( 16), for l = 0, 1, 2, 3, 4. The optimal exponents for l = 1 have been already published in ref. 26 and were obtained with slightly different numerical parameters in the optimization code.

Figure 1 :

 1 Figure 1: Top panel : regular Coulomb function with wavenumber k e = 1 a.u. (thin continuous line), superimposed with its cGTO fit (Real part: dashed line; null imaginary part: dotted line). The bottom panel shows the corresponding absolute errors.

NH 3 (3a 1 )Figure 2 :

 312 Photon energy (eV)

Figure 3 :Figure 4 :Figure 5 :H 2 O (1b 1 )Figure 6 :

 345216 Photon energy (eV)

Figure 7 :

 7 Figure 7: Partial cross-section σ(k e ) (top panel) and asymmetry parameter β(k e ) (bottom panel) as a function of the photon energy (in eV), for orbital 1b 1 of H 2 O. Present results using cGTOs integrals (cG) with a distorted continuum wave are compared with results from other theoretical methods (TD-DFT by Stener et al 55 , single-center method of Novikovskiy et al 14 and EOM-CCSD Dyson/B-spline DFT calculations of Moitra et al 3 ), and experimental points (Banna et al 60 , Truesdale et al 59 and Tan et al 58 ). Our calculations are performed in both length (L) and velocity (V) gauges.

Figure 9 :

 9 Figure 9: Same as Fig. 7 for orbital 1b 2 of H 2 O.

Figure 10 :

 10 Figure 10: Same as Fig. 7 for orbital 2a 1 of H 2 O.

  Figure 8: Same as Fig. 7 for orbital 3a 1 of H 2 O.

						H 2 O (3a 1 )
	Cross Section (Mb)	0 1 2 3 4 5 6 7 8 9 11 10
		1.3
		1.0
	β	0.7 0.4			cG-(L) cG-(V)	Moitra et al Banna et al
		14 -0.2 0.1	19	24	29 Pho)on ener y (eV) 34 39 44 TD-DFT Novikovskiy et al	49 Truesdale et al 54 59 Tan et al
						H 2 O (1b 2 )
	Cross Section (Mb)	0 1 2 3 4 5 6 7 8 9 10 11 12			cG-(L) cG-(V) TD-DFT Novikovskiy et al	Moitra et al Banna et al Truesdale et al Tan et al
		1.2
		0.9
	β	0.6
		0.3
		0.0
	3 CONCLUSION 18 22 26 30 34 38 42 46 50 54 58 62 Photon energy (eV) -0.3
	Studying molecular photoionization processes requires a good description of the outgoing
	electron. At first sight, Gaussian representations do not seem very appropriate for highly

Table 1 :

 1 Optimal cGTO exponents {α i } l obtained after fitting the sets of Coulomb functions

	i	{α i } l=0	{α i } l=1	{α i } l=2	{α i } l=3	{α i } l=4
	1	0.00035555 + 0.00901614ı	0.0001000 + 0.0151526ı	0.00100000 + 0.00845442ı	0.00014763 + 0.00587055ı	0.00039019 -0.01494676ı
	2	0.00053211 + 0.02065728ı	0.0001418 -0.0171753ı	0.00121886 + 0.01803494ı	0.00109586 + 0.02784779ı	0.00262171 -0.03255562ı
	3	0.00080009 -0.01912061ı	0.0013875 -0.0021493ı	0.00148684 + 0.01813637ı	0.00148998 -0.01168882ı	0.00339622 + 0.01408910ı
	4	0.00117085 -0.01123799ı	0.0019035 + 0.0216135ı	0.00181381 + 0.02949728ı	0.00205346 + 0.01943320ı	0.00443045 -0.02409112ı
	5	0.00171276 + 0.02161259ı	0.0025827 + 0.0298661ı	0.00220490 + 0.02447032ı	0.00280227 + 0.03734991ı	0.00574632 + 0.02828016ı
	6	0.00249700 + 0.03158584ı	0.0034908 -0.0241050ı	0.00267806 -0.01024824ı	0.00358933 + 0.01238024ı	0.00744649 + 0.02342353ı
	7	0.00351122 -0.03766107ı	0.0046754 -0.0338441ı	0.00325391 -0.01961417ı	0.00459090 -0.02321602ı	0.00970895 + 0.03880653ı
	8	0.00494586 -0.02649727ı	0.0061876 + 0.0378398ı	0.00395324 -0.02477460ı	0.00588757 -0.02219062ı	0.01223182 + 0.02067056ı
	9	0.00695348 -0.02692584ı	0.0156781 + 0.0504851ı	0.00480634 -0.01801510ı	0.00757308 -0.03698798ı	0.01537803 -0.03985771ı
		0.00974499 -0.04400129ı	0.0196315 -0.0395359ı	0.00585038 -0.02149702ı	0.00967950 -0.02915329ı	0.01953215 -0.03489981ı
		0.01341500 + 0.03611724ı	0.0246020 -0.0352098ı	0.00712456 -0.02735801ı	0.01239829 + 0.04894840ı	0.02499666 + 0.05443031ı
		0.01853987 + 0.04014546ı	0.0308382 -0.0302620ı	0.00867609 -0.03810200ı	0.01699053 -0.05180981ı	0.03507324 -0.05430569ı
		0.02554489 -0.05549019ı	0.0386508 -0.0436192ı	0.01058491 + 0.03742598ı	0.03292671 + 0.06650370ı	0.04435803 -0.02396048ı
		0.03437502 + 0.04391841ı	0.0484105 + 0.0624525ı	0.01305188 + 0.03757777ı	0.04217102 -0.07408910ı	0.05617605 + 0.02052618ı
		0.04626602 + 0.05104620ı	0.0618203 + 0.0513264ı	0.01646720 -0.05003223ı	0.06196807 + 0.03814310ı	0.07128007 + 0.02106443ı
		0.06218273 -0.06475675ı	0.0816826 -0.0552154ı	0.02076651 + 0.05258585ı	0.19225626 + 0.03964722ı	0.16136565 + 0.02940124ı
		0.08554072 -0.02426909ı	0.1936246 + 0.0097332ı	0.03514922 -0.06965770ı	0.30147472 + 0.00811515ı	0.39894781 + 0.02631795ı
		0.11861245 + 0.07471136ı	0.2915597 + 0.0024838ı	0.04424413 + 0.07302765ı	0.51617609 + 0.04399335ı	0.55207468 + 0.07285449ı
		0.16691697 + 0.08944899ı	0.4971406 -0.0113366ı	0.09707468 -0.10000000ı	0.82002155 + 0.05919349ı	0.86945734 -0.01791204ı
		0.23869247 -0.09999999ı	0.8532849 + 0.0200975ı	0.12043971 + 0.08910818ı	1.19395577 + 0.01279526ı	1.12048249 + 0.08602790ı
		1.19614463 -0.02970423ı	1.3797714 + 0.0062079ı	2.10573048 -0.03921720ı	1.48783709 + 0.02469442ı	1.41353953 + 0.04522398ı
		5.33188880 -0.01441299ı	2.1958583 -0.0311758ı	3.94194435 + 0.02888234ı	2.19692967 -0.08785385ı	2.14872735 -0.01352875ı
		7.71182758 -0.09811990ı	3.5672935 -0.0087858ı	6.10200559 -0.01548483ı	3.54466390 + 0.02398460ı	3.60699097 -0.03483872ı
		10.9962323 -0.09305650ı	5.7419028 + 0.0012293ı	9.14648581 -0.00307356ı	5.54622769 + 0.00914500ı	5.64732513 + 0.01125528ı
		15.5241861 -0.09730549ı	9.2535835 -0.0212019ı	13.7297058 -0.03200391ı	9.24463017 -0.00535522ı	9.19300771 + 0.04174570ı
		21.7519013 -0.09157168ı	14.868739 -0.0267870ı	20.4351559 -0.00138094ı	14.8708853 -0.02233793ı	14.8869767 -0.03464065ı
		30.2989535 -0.07925670ı	23.954872 -0.0406721ı	30.3651118 + 0.02650005ı	23.9508076 + 0.08258257ı	23.9066166 + 0.00726798ı
		42.0163884 -0.09927635ı	38.575864 + 0.0168896ı	45.2156191 -0.00722570ı	38.5477254 + 0.09948084ı	38.5167312 -0.01803004ı
		58.1483705 -0.08991001ı	62.088037 -0.0059998ı	67.2295184 + 0.01284461ı	62.0907276 -0.03652093ı	62.0736229 -0.06204804ı
		80.7346781 -0.08528403ı	99.986651 + 0.0164035ı	99.9856251 -0.06227840ı	100.021625 + 0.00412096ı	99.9767615 -0.03846239ı

ACKNOWLEDGMENTS

We acknowledge Torsha Moitra, Sonia Coriani and Nikolay M. Novikovskiy for sending us their raw data associated with refs 3 and 14. The PMMS (Pôle Messin de Modélisation et de Simulation) is also acknowledged for providing us with computer time.

DATA AVAILABILITY STATEMENT

The data that support the findings of this study are available from the corresponding author upon reasonable request.


Though much critiqued and criticised, the generally ill-defined notion of the "failed state" continues to dot and enthral academic and even more so public debates. Over the years, it has become a habit to diagnose an increasing number of internationally recognised states marked by serious internal conflicts and flagging delivery of public goods as "failing", "failed", or "collapsed". More often than not, the label serves to legitimate external intervention, including projects of reconstruction and "state re-building" to contain or combat domestic disorder supposed to threaten the rest of the world, be it in the form of migration or terrorism [START_REF] Woodward | The Ideology of Failed States: Why Intervention Fails[END_REF]). In the Middle East and North Africa (MENA) area alone, the dubious honour of this title has been repeatedly awarded to Iraq, Syria, Lebanon, Yemen, Sudan, and Libya.

The present paper does not seek to repeat the countless attempts to question popular definitions of "state failure", or to refine them, or to endorse the often fruitful search for alternatives such as "limited statehood" (Risse 2011). Rather the aim is to start from a simple question that is rarely asked: if states in the MENA area and beyond are actually failing or have failed for protracted periods of time-in some cases several decades-then why have they not completely disappeared as actors, coloured patches on maps, names, and rhetorical references? For example, in spite of some expectations, Iraq since the overthrow of Saddam Husayn has not broken up into a Shi'i state in the south, a Sunni one in the centre, and a Kurdish one in the north-east. No doubt, the Kurdish area is in many respects autonomous, and demands for full independence have been strong; however, the area has also maintained ties with the rest of Iraq, such as the use of the dinar as legal tender, and after the 2017 referendum it remained part of the state.

The first part of the seemingly obvious answer to the question of why ostensibly failed states have failed to disappear is that borders once drawn by governments, often the major powers of historical imperialism, and endorsed by official treaties and the "international community" need not necessarily reflect the will of those who live in the lands and geographies concerned. Numerous borders thus reflect the interests of external actors and related path dependencies which continue to be defended, enforced, and reproduced by these very actors. In the case of Iraqi Kurdistan, there was indeed no shortage of statements and other pressures from around the world discouraging the referendum and a vote in favour of independence. The second part of the seemingly obvious answer is that the external defenders of such "artificial" states provide them with the means that are necessary to avoid complete disintegration and to ensure minimal performance. For instance, the Iraqi monarchy was propped up by the British air force which dropped bombs on recalcitrant "tribes", while the current Assad government in Syria is under military transfusion from Russia, Hizb Allah, and Iran, the latter occasionally shipping crude oil to the country.

In a sense, the "obvious" answer dovetails with Robert H. [START_REF] Jackson | Quasi-States: Sovereignty, International Relations, and the Third World[END_REF] famous distinction between "states" and "quasi-states": many of the "new" states in the areas formerly occupied by European powers were created by these very powers rather than by local actors, forces, and alliances. In that sense their origin and trajectory differed substantially from that of historically consolidated states such as Egypt, Tunisia, and Moroccowhich does not always prevent observers from wondering whether the latter will be able to resist "failure" in the light of the severe economic, ecological, and attendant social and societal challenges they are facing. Due to the support of their creators, the "new" states managed to conquer a place on the map independently of their ability to effectively monopolise the means of coercion, to control the territory, and to pursue coherent policies. In other words, they obtained at some stage "negative sovereignty", defined as the absence of-all too visible-external interference, but they failed to acquire "positive sovereignty", defined as the ability to decisively shape and reshape matters at home. As a result, they remained quasi-states rather than states as such. Though too neat to be true and not always confirmed by developments on the ground, the distinction nonetheless pinpoints two different aspects of the leeway states have as actors that need not necessarily go hand in hand. States and their borders may be, or remain, internationally recognised, even when their governments lose control over much of the territory, when their presidents more or less temporarily morph into glorified mayors of the capital city, or when they leave the country altogether. It should be added that, in practice, the international recognition of a state and its borders need not prevent those offering such recognition from supporting opposition or rebel forces within that state. However, the staying power of states on the map and indeed on the ground cannot simply be reduced to an issue of international recognition or to a complete fiction. Even the states that at the outset were largely built by foreign powers consisted not only of a territory and a population but also of a political regime. A set of institutions, including organisations, the political regime impacted the lives of the inhabitants and influenced their activities, not least their political activities and choices, through constraints and incentives. Material and symbolical values were distributed by means of legislation and other forms of regulations or simply by brute, physical force. The repetition of government action and specific policies, for instance in the area of education or indoctrination, as well as the ensuing expectations and path dependencies, added "soft" or "infrastructural" power to the "hard" or "despotic" power of coercion (see [START_REF] Mann | The Autonomous Power of the State: Its Origins, Mechanisms, and Results[END_REF][START_REF] Mann | The Sources of Social Power[END_REF]. This applies even to countries where the reach of the central government was highly uneven and shaped the lives of people far less than, say, primary schools in the French Third Republic. In a nutshell, institutions, even imposed from the outside, gained weight and weighed down on the individuals and groups that had to face them. Recalling the triptych of exit, voice, and loyalty [START_REF] Hirschman | Exit, Voice, and Loyalty: Responses to Decline in Firms, Organizations, and States[END_REF], they led inhabitants not only to retreat into apathy or to forcefully stage protests, but also to play by the imposed rules and accept them. As long as some rules were observed by some, the state was-however precariously-kept alive and managed to resist-however imperfectlythose who challenged its existence.

As a matter of course, the "new" states such as Iraq, Syria, or Lebanon were heavily contested when they were established as dependent political entities and even when they became formally independent members of an emerging international community epitomised by the League of Nations and the United Nations. Many inhabitants of what became independent Syria, Lebanon, and Iraq considered the borders established by the victors of the First World War as imposed and "artificial". Although not necessarily mourning the predecessor state, the Ottoman Empire, they developed and promoted projects of Greater Syrian, Fertile Crescent, or Arab unity. Like more or less explicitly secessionist projects, these schemes, with important variations over time, have heavily influenced the politics of the countries ever since. No less than the borders of state territory and the criteria governing citizenship or the right of abode, the forms of government in these places have been contested, resulting in coups and the succession of rather different, if mainly authoritarian, political regimes. However, especially from the late 1950s and early 1960s when governments came to power and for a variety of reasons emphasized the delivery of public goods, growing parts of the population developed increasing stakes in these states and considered their own welfare as closely tied to the states' survival. Public-sector employment, be it in the administration, armed forces, or state-owned companies, as well as related social benefits established bonds of loyalty to the state. These distributive measures were partly the result of a mix of coercion and cooptation, as manifest in land reforms and the expropriation of private capital. Partly, such redistribution was the result of the very existence of the stateshowever "artificially" created. In actual fact, the latter was not only a legal construct endowed with negative sovereignty but, by virtue of such sovereignty, was also the recipient and domestic distributor of various types of external resources, including development aid, military support, and budget transfers from wealthier neighbours seeking to stabilise or instrumentalise their poorer counterparts. No less importantly, thanks to its power to sign treaties, the state was the doorkeeper of access to trade regimes and other international economic arrangements. Throughout the decades, forces which contested polices, governments, and the states themselves faced other forces which defended these very policies, governments, or states, either selectively or as a package. The question is not whether the forces of contestation were stronger or weaker than the ties of loyalty; the important observation is that in many periods and at many junctures there have been sufficient domestic defenders of states who kept them alive, even when, in the eyes of others, they had lost their legitimacy or usefulness.

In a world subdivided into states, the very existence of a state contributes to an extent to consolidate it, even against such odds as civil wars and the emergence of powerful challengers who dominate large swaths of its land, as did for instance the Islamic State. This is not to say that states are eternal or "functional" by definition, only that they do command resources that help them to resist.