Total photoionization cross section of planar helium: scaling laws and collision orbits

The total photoionization cross section of planar helium has been calculated up to the single ionization threshold I22 of triple P states. The cross section shows chaotic fluctuations as the energy E approaches the double ionization threshold E = 0. By analyzing the fluctuating part of the cross section, we show that its amplitude decreases as ∣ E ∣ μ for E → 0 − as predicted in Byun et al (2007 Phys. Rev. Lett. 98 113001). The Fourier transform of the fluctuating part reveals peaks at the classical actions of closed triple collision orbits. Furthermore, the relative height of the peaks is consistent with the semiclassical predictions. Our findings underline that the fluctuating part of the photoionization cross section can be described by classical triple collision orbits in the semiclassical limit. These orbits all lie in the collinear eZe subspace, demonstrating that the fluctuations are dominated by the dynamics of this low dimensional phase space.


Introduction
Since Madden and Codling observed the auto-ionization of helium from the ground state to the doubly excited state in 1963 [1], two-electron atoms have served as a prototype for understanding electron-electron correlations in atomic systems. To describe two-electron atoms fully quantum mechanically, many efficient numerical methods, such as exterior complex scaling [2,3], convergent closed coupling methods [4,5], time-dependent close-coupling methods [6][7][8][9][10] or advanced hyperspherical R-matrix techniques using semiclassical outgoing waves [11][12][13], have been applied.
In recent years, the photoionization cross section has been fully explored experimentally up to the N≈15 ionization threshold by virtue of the experimental advances in synchrotron radiation sources and detector technology [14,15]. However, only the spectrum of the bound states and low-lying doubly excited resonances is well understood, see for example [11] for an overview. Furthermore, as one approaches the double ionization threshold E=0 from below, the photoionization cross section shows highly chaotic behavior in the region of overlapping Rydberg series [11]. In [16,17], based on the extended semiclassical closed orbit theory, it is shown that this behavior is intricately linked to the complexity of the underlying classical dynamics of three-body Coulomb problem helium. In particular, the paper predicts two semiclassical results: (i) the Fourier spectrum of the fluctuating part of the photoionization cross section shows distinct peaks at the positions of the actions of so-called closed triple collision orbits (CTCO), that is, classical trajectories of coupled two-electron dynamics starting and ending in nonregularizable triple collision; (ii) the amplitude of the fluctuations in the total and partial cross sections for single electron photoionization in twoelectron atoms decays algebraically in terms of a threshold law m | | E as  E 0, where the exponent μ is given as Re 100  9  4  1  2  4  9  4  1  1 with Z, the charge of the nucleus. The exponent is related to the stability exponents of the CTCOs and differs from Wannier's exponent [18] dominating double ionization processes. The prediction implies that the photoionization cross section can be described semiclassically in terms of classical triple collision orbits all lying in the collinear eZe space. By numerically calculating the total photoionization cross section (TPICS) of collinear eZe helium up to the 50th single ionization threshold, I 50 , the semiclassical predictions have been confirmed in [16,17]. For fully three-dimensional helium, however, a verification of these semiclassical results based on numerical computations or experiments has been missing so far. By analysing the experimental data in [14] on partial photoionization cross sections (PPICS) for helium up to N≈13, it could be shown in [19] that the fluctuations in the PPICS are semiclassically dominated by the contributions from closed triple collision orbits. The threshold laws could, however, not be extracted due to lack of data. In [15], Jiang et al have measured the total photoionization cross section up to single ionization threshold I 15 of He + and computed the spectrum up to I 17 numerically by using the complex rotation method [20]. It was shown that the spectra are dominated by the principal Rydberg series; however, semiclassical results could not been extracted due to lack of resolution.
For planar helium, the fluctuating parts of TPICS have been calculated up to the 20th single ionization threshold of triplet P states using the complex rotation method [21]. The authors of [21] argued that the fluctuations in the low energy region are mainly due to a dominant series of resonances associated with an approximate quantum number F=N− K, related to the collinear eZe dynamics. But as the energy increases, new series start to contribute significantly to the cross section, and it is suggested that the scaling law starts to become invalid. This is in contrast to findings in [16,17], where it is predicted that planar helium should also follow the scaling law and other semiclassical results stated above.
Based on ab initio numerical calculation for the total photoionization cross section of planar helium and by employing a quantum map approach [16,17], we will show in this paper that the above mentioned semiclassical results are indeed valid. The paper is organized as follows. In section 2, we will briefly describe the semiclassical theory and the numerical implementation of the ab initio method. In section 3, the numerical result of the photoionization cross section of the planar helium up to I 22 is presented. By analysing the numerical data, the semiclassical results are confirmed for the first time in higher dimensions. Section 4 concludes the paper. Unless stated otherwise, atomic units are used throughout this paper.

A quantum map approach to the total photoionization cross section
In the Poincaré map approach introduced in [17], one starts from writing the photoionization cross section σ(E) in terms of the retarded Green function G(E), Here, f 0 is the wave function of the initial state and  = + · ( ) D r r 1 2 is the dipole operator with  , the polarization of the incoming photon and r i , the position of the i-th electron. The direction of the polarization is taken along the z-axis.
Partitioning the whole configuration space into inner and outer regions by introducing a dividing surface Σ, we can express the cross section σ(E) in terms of the local scattering matrices r, s and the scattering solutions of the half-space Hamiltonians related to the inner and outer regions as where Σ is taken here at a fixed hyperradius R=R 0 with = + R r r satisfying outgoing boundary conditions and H I is the Hamiltonian for the inner space. The subscript o denotes open channels for the scattering system described by H I . The merit of this formula is that one can divide σ(E) into a smooth background signal σ bg and a fluctuating part of the photoionization cross section. The latter contains information about doubly excited states and thus most of the resonance states making up the chaotic fluctuations in the cross section for E → 0 − . Starting from equation (3), one can now derive the semiclassical results described in the introduction as detailed in [17].

Quantum mechanical treatment in hyperspherical coordinates
The Hamiltonian of the helium atom with infinite nucleus mass fixed at the origin is written as where r i and = i p , 1, 2 i denote the positions and the conjugate momenta of electron i respectively. Z is the nuclear charge. We will consider in the following the planar helium atom restricted to the x-y plane, that is, both electrons move in the z=0 plane. This is an invariant subspace of the full classical dynamics for initial conditions and is considered here as a lower-dimensional problem quantum mechanically exhibiting the proliferation of Rydberg series typical for the full helium problem.
In hyperspherical coordinates (R, α, θ 1 , 2 1 and θ 1 and θ 2 , the polar angles with respect to the x-axis, the Schrödinger equation is written as 3 2 . Here, Ȳ is a solution of the Schödinger equation )H E 0 and H R is the adiabatic Hamiltonian for angle variables Ω=(α, θ 1 , θ 2 ). The adiabatic Hamiltonian H R is written as with azimuthal angular momentum operators of the electrons 1 and 2, where θ=θ 2 −θ 1 . The total angular momentum is where Φ n (Ω; R) are the eigenfunctions of the adiabatic Hamiltonian, that is, ; . 1 0

R n n n
Substitution of the expansion equation (9) into the Schrödinger equation (5) produces a set of N CH ×N CH coupledchannel equations about F(R), where N CH is the number of channels included in equation (9). The coupled-channel equations are solved using the SVD generalized log-derivative propagation method, called L-matrix propagation [17,[22][23][24], which is highly suitable for parallel computation. In the limit  ¥ R , we can ignore the channel-channel coupling and we used the WKB form of the channel equation for F n (R) described in [25] as boundary condition in this limit. The equations for F n (R) reduce then to single channel equations of the form for each channel n. The asymptotic form of potentials U n are given by n n 2 where a n and b n are channel dependent constants. The dividing surface Σ is taken at hyperradius R 0 =20. Note that if R 0 is greater than the support of the initial wave function f 0 , the TPICS results become independent of R 0 . By solving the Schödinger equation for the half-space Hamiltonians, the scattering matrices s, r, and the dipole transition amplitude d, necessary for equation (3), are obtained.
For the Hamiltonian (4), there are five symmetry operators, namely Π x , Π y , Π, P 12 , L z [26]. Here, Π x and Π y are parity operators about the x and y-axis respectively, Π is a point reflection (  -r r r r , 1 1 2 2 ), L z coincides here with the total angular momentum operator, and P 12 is the particle exchange operator. The spectrum of H are classified by taking the mutual commuting operators P 12 , Π x , and L z .
For a given total angular momentum l the adiabatic function Φ(Ω; R) can be expanded as where Φ m (α; R) is a function of α at fixed R. By rearranging the basis, we can write the function Φ l (Ω; R) in terms of fixed quantum numbers also for Π x and P 12 . For example, the triplet states with Π x =1, P 12 =1 and l=1 or l=0 can be written as where A m (α; R) and B m (α; R) are functions of α at a given R. We expand A m (α; R) and B m (α; R) in terms of the B-spline basis function [27]. In this basis, the adiabatic Hamiltonian H R is represented by a banded matrix of large size, e.g., 9572×1055 for triple state in the range I 19 to I 20 . We calculate the adiabatic potential U n (R) and the channel functions Φ n (Ω; R) up to N CH by using the ARPACK package [28].

Numerical results
The initial state is taken as the lowest lying triplet bound state, with angular momentum l=0 and Π x =+1 as done in [21]. In order to get the initial state, we use the matching method introduced in [ which is very close to the value given in [21]. Here nbps is the number of B-spline basis functions for a M , max denotes the number of basis functions for θ, and h is the step size in the R direction.
The final state Y + ( ) E f is taken as the scattering state satisfying outgoing boundary conditions with l=1, Π x =+1, and P 12 =1 (triplet state). By employing the mixed L-and S-matrix propagation method [17], we obtain the scattering matrix r for the outer region ( > R R 0 ) and Y + ( ) E f in the inner region (R<R 0 ). The outgoing boundary condition for r is applied at a sufficiently large = ¥ R R . The dipole transition amplitude d is obtained using the J-integral algorithm [29] for parallel computation. The parameter values used are varied with the energy region considered as shown in 200, and N CH =200. The value is in good agreement with [31]. Note that the state is different from the initial ground state mentioned above, of which the angular momentum is not l=1 but l=0.
In figure 1, the total photoionization cross section (TPICS) is presented as a function of N for the energy region of I 1 to I 15 , where I N and N are related by We can see some peaks in the smooth part σ bg at low energies with N<4, which are related to resonant states localized in the region R<R 0 . However, σ bg is smooth for N>4 and the high lying resonances are indeed all contained in the fluctuating part of σ as expected. This is because the support of the resonant wave functions expands far beyond the boundary R=R 0 and their contribution to σ bg becomes negligible. σ bg converges to a constant value as N increases. σ fl , on the other hand, is highly fluctuating for N>4 with decreasing amplitude of fluctuation as N increases. We note a slight overall increase of the average value of σ fl which is a remanent of the smooth part contributing to the R>R 0 region. Such behavior has also been seen in the cross section σ fl of eZe collinear helium [17].
In figure 2, the TPICS and its fluctuating part are plotted from I 10 to I 22 . Because of the non-zero average for σ fl , a smooth part is subtracted directly from total the cross section σ(E). This smooth part is obtained by fitting the function f (N)=a+b/N 2 to σ(E) and subtracting it to obtain s ( ) E fl as shown in figure 2(b). We note again that the amplitude of the fluctuating signal decreases as  | | E 0 from below. Hereafter, σ fl will mean the fluctuation part obtained by the above fitting method.
In order to confirm the semiclassical prediction (i) stated in section 1, we calculate the Fourier transform (FT) of the signal in the from , is the scaling exponent given in equation (1) with μ=1.306 for Z=2, and E 0 is the energy of the initial state. The factor m -| | E is introduced to compensate for the decrease in amplitude. In figure 3, the Fourier transform of F(z) with respect to z is plotted as a function of the scaled variable = | | S S E ; note that the variable S can be identified with a classical action. When comparing the position and height of the peaks in the FT of the total cross sections (obtained from a full quantum calculation) with the semiclassical results deduced from the actions and stabilities of closed triple collision orbits (CTCO) shown as blue circles, one finds excellent agreement. The position of the circles represents the scaled classical action of a specific CTCO whose trajectory is shown as an inset above the peak in figure 3. The heights are obtained by considering the stability matrix of the orbit as described in [17]. In the figure, the relative heights are plotted. For the range below p » ( ) S 2 8.6 we can see that there is a one to one correspondence between the peaks of the quantum results and the positions of the classical actions of CTCOs. For larger actions, the correspondence deteriorates due to the finite width of the quantum signal σ(E) and thus finite resolution of  the FT. We expect that further peaks will be resolved if the TPICS is obtained for higher N such as has been demonstrated in the lower-dimensional eZe collinear helium in [16,17]; here TPICS data were obtained up to N=50 and peaks were well resolved up to p » ( ) S 2 15. The FT result for planar helium is almost identical to one for eZe collinear helium, except for the resolution.
A validation of the scaling law, the second semiclassical prediction in section 1, is obtained by extracting a decay exponent μ from the strongly fluctuating signal shown in figure 2. This is a somewhat delicate task and is done here by fitting a decay law directly to the signal σ fl after suitable smoothing. In order to avoid artifacts due to the variation in curvature in the original signal σ fl (N) when using large averaging intervals, we average the absolute value of σ fl (N) over fixed intervals of size ΔN, that is, we consider The results are shown in figure 4 on a log-log scale, here for ΔN=0.1, 0.5, 1.0. The smoothed data show a linear behavior on the log-log scale indicating a single-valued decay law. We extract the slope of the curve in figure 4 using a linear regression model with least-square fitting, that is, with fit parameters a and b. The parameter b is compared with the theoretically predicted decay exponent μ, see inset in figure 4. We would expect that b is independent of N if N is sufficiently large, thus b should be independent of the fitting range. We take two regions N= [10,16] and N= [16,22] from the whole region N= [10,22], and then determined the slope of the averaged s | | fl for each regions separately. The results for the fit parameter b for each of these intervals are shown in the inset of figure 4 for different averaging intervals ΔN. The data coincide with the expected value within ±6%, that is, we obtain an exponent of the order μ=1.3±0.08.

Conclusions
By using an ab initio numerical method, the total photoionization cross section has been calculated for planar helium up to I 22 for the first time. The analysis of the fluctuating part of the TPICS confirms the semiclassical predictions in [16]. This is further strong evidence that the fluctuation is completely dominated by the classical triple collision orbits all lying in the collinear eZe subspace embedded in the three-dimensional space. Although an analysis of the full 3D helium problem is still out of reach both theoretically and experimentally due to the limitations in resolving the full spectrum for higher energies up to N≈20, this gives further confidence in expecting that semiclassical predictions are also valid in the case of the full helium problem as well as other two-electron atoms.