Breakup of the aligned H 2 molecule by xuv laser pulses: A time-dependent treatment in prolate spheroidal coordinates
aa r X i v : . [ phy s i c s . a t o m - ph ] J a n submitted to Phys. Rev. A; 7 January 2011 Breakup of the aligned H molecule by xuv laser pulses: A time-dependent treatmentin prolate spheroidal coordinates Xiaoxu Guan , Klaus Bartschat , and Barry I. Schneider Department of Physics and Astronomy, Drake University, Des Moines, Iowa 50311, USA and Physics Division, National Science Foundation, Arlington, Virgina 22230, USA (Dated: June 29, 2018)We have carried out calculations of the triple-differential cross section for one-photon doubleionization of molecular hydrogen for a central photon energy of 75 eV, using a fully ab initio ,nonperturbative approach to solve the time-dependent Schr¨odinger equation in prolate spheroidalcoordinates. The spatial coordinates ξ and η are discretized in a finite-element discrete-variablerepresentation. The wave packet of the laser-driven two-electron system is propagated in timethrough an effective short iterative Lanczos method to simulate the double ionization of the hy-drogen molecule. For both symmetric and asymmetric energy sharing, the present results agreeto a satisfactory level with most earlier predictions for the absolute magnitude and the shape ofthe angular distributions. A notable exception, however, concerns the predictions of the recenttime-independent calculations based on the exterior complex scaling method in prolate spheroidalcoordinates [Phys. Rev. A , 023423 (2010)]. Extensive tests of the numerical implementationwere performed, including the effect of truncating the Neumann expansion for the dielectronic in-teraction on the description of the initial bound state and the predicted cross sections. We observethat the dominant escape mode of the two photoelectrons dramatically depends upon the energysharing. In the parallel geometry, when the ejected electrons are collected along the direction of thelaser polarization axis, back-to-back escape is the dominant channel for strongly asymmetric energysharing, while it is completely forbidden if the two electrons share the excess energy equally. PACS numbers: 33.80.-b, 33.80.Wz, 31.15.A-
I. INTRODUCTION
A measurement of the complete breakup of the atomichelium target by xuv radiation was achieved over 10 yearsago [1]. Since then, rapid developments in strong xuvlight sources and momentum imaging techniques havemade it possible to record all the reaction fragments,nuclei and electrons, in double photoionization of thesimplest two-electron hydrogen/deuterium molecule byone-photon absorption [2–5], and most recently also fortwo-photon absorption [6]. For the double ionizationof H by single-photon absorption, only randomly ori-ented molecules were investigated in earlier experiments(e.g. [7]). Using “fixed-in-space” molecules, more recentexperimental efforts include the measurements of energy-and angle-resolved differential cross sections by Weber etal. for either equal energy sharing [2, 3] or asymmetricenergy sharing [4], and by Gisselbrecht et al. [5] for equalenergy sharing at a photon energy of 76 eV. These ex-perimental studies were at least partially stimulated bythe goal of understanding the similarities and differencesbetween the hydrogen molecule and its atomic counter-part helium. However, all recorded fully differential crosssections to date suffer from some experimental uncertain-ties regarding the alignment angle of the molecule withrespect to the polarization vector of the laser and theemission angles of the photoelectrons.From a theoretical point of view, the hydrogenmolecule exhibits a significant complexity compared tohelium and, therefore, provides an enormous challengeto a fully ab initio description inherent in a multi- center, multi-electron system. The single-center con-vergent close-coupling method was used to model thedouble ionization of H by Kheifets and Bray [8, 9].Later McCurdy, Rescigno, Mart´ın and their collabora-tors [10–12] implemented a formulation based on time-independent exterior complex scaling (ECS) in spheri-cal coordinates, with the origin of the coordinate systemplaced at the center of the molecule, to treat the doublephotoionization at a photon energy of 75 eV. The radialcoordinates of the two electrons are measured from thecenter, and the radial parts of the wave function wereeither expanded in B -splines or using a finite-elementdiscrete-variable representation (FE-DVR).The time-dependent close-coupling (TDCC)method [13], again in spherical coordinates, wasalso extended to calculate the triple-differential crosssection (TDCS) for double photoionization of the H molecule. While the agreement between the publishedTDCSs from the ECS [11, 12] and TDCC [13] calcu-lations is basically acceptable, noticeable discrepanciesremained for a few particular geometries. In the parallelgeometry, for instance, where the molecular axis ζ ischosen along the laser polarization vector ǫ , the coplanarTDCS predictions from the ECS and TDCC calculationsdiffer by up to 40 percent when one of the electrons (wewill refer to it as the “fixed electron” below) is observedalong the direction perpendicular to the ǫ -axis. In someother cases, there exists a noticeable “wing” structurein the published TDCC predictions for equal energysharing. Additional TDCC calculations [14] suggestthat the agreement can be systematically improved,albeit the above-mentioned discrepancy still exists at asomewhat reduced level.Another independent approach [15] to this problem isthe very recent time-independent ECS treatment, for-mulated – as in the current work – in prolate spheroidalcoordinates. Quite surprisingly, the results of that calcu-lation differed from both the earlier ECS and also theTDCC predictions, both obtained in spherical coordi-nates. Specifically, the ECS prolate spheroidal calcula-tions showed differences from the earlier spherical coor-dinate calculation for the TDCSs, at a level of about20% depending on the details of the energy sharing. Aswill be demonstrated below, we have gone to consider-able lengths in an attempt to resolve these discrepancies.However, significant differences between the present re-sults and those of the ECS [15] still remain.Both the ECS [11, 12] and the TDCC [13] calcula-tions made some attempt to deal with the experimentaluncertainties in the scattering angles. Given the experi-mental uncertainties and the differences in the previouscalculations, however, it appeared worthwhile to inves-tigate the computational effort required to obtain accu-rate TDCSs before averaging over any experimental ac-ceptance angles. Consequently, the present calculationrepresents an independent implementation of the time-dependent FE-DVR approach in prolate spheroidal coor-dinates. As in the other approaches mentioned above, theinternuclear separation ( R ) was held fixed at its equilib-rium distance of 1 . molecule. The formulation of the Schr¨odinger equationin prolate spheroidal coordinates for diatomic moleculesis not new. The pioneering work of Bates, ¨Opik, andPoots [16] for the H +2 ion, which is exactly separable inprolate spheroidal coordinates, already revealed the ap-pealing features of the prolate system. In particular, theelectron-nuclear interaction is rendered benign in this co-ordinate system. A partial list of recent applications ofprolate spheroidal coordinates to diatomic molecules canbe found in [17–21].As has been demonstrated in a number of recent pub-lications, a grid-based approach provides a very ap-propriate description of laser-driven atomic and molec-ular physics when combined with an efficient time-propagation method such as the short iterative Lanczos(SIL) method [22, 23]. In the present work, we employthe FE-DVR/SIL approach in prolate spheroidal coordi-nates to study the correlated response of a two-electronmolecule in the double ionization process.The remainder of this manuscript is organized as fol-lows. After presenting the Hamiltonian of the hydrogenmolecule in Section II and providing some details aboutthe discretization of the system in an FE-DVR basis inSec. III, the solution of the two-center Poisson equationis presented in Sec. IV. This is followed by a descriptionof the procedure for extracting the cross sections of in-terest in Sec. V. The results are presented and discussedin Sec. VI, before we finish with a summary in Sec. VII. II. THE SCHR ¨ODINGER EQUATION INPROLATE SPHEROIDAL COORDINATES
The prolate spheroidal coordinates with the two fociseparated by a distance R are defined by ξ = r + r R , η = r − r R , (1)and the azimuthal angle ϕ . Here r and r are the dis-tances measured from the two nuclei, respectively. Thesecoordinates are specified in the ranges ξ ∈ [1 , + ∞ ), η ∈ [ − , +1], and ϕ ∈ [0 , π ]. The volume element is dV = ( R/ ( ξ − η ) dξdηdϕ . According to the asymp-totic behaviors as r , r → + ∞ , ξ and η approach 2 r/R and cos θ , respectively, where r and θ are the standardspherical coordinates. Consequently, ξ is the “quasi-radial” coordinate, while η is “quasi-angular”. TheHamiltonian of a single electron, H q , ( q = 1 , below), is written as H q = − R ( ξ q − η q ) " ∂∂ξ q ( ξ q − ∂∂ξ q + ∂∂η q (1 − η q ) ∂∂η q (2)+ 1( ξ q − ∂ ∂ϕ q + 1(1 − η q ) ∂ ∂ϕ q − ξ q R ( ξ q − η q ) . We solve the time-dependent Schr¨odinger equation(TDSE) of the laser-driven H molecule (with two elec-trons) in the dipole length gauge: i ∂∂t Ψ(1 , , t ) = h H + H + 1 r + E ( t ) · ( r + r ) i Ψ(1 , , t ) . (3)Here r q is the coordinate of the q -th electron measuredrelative to the center of the molecule and r = | r − r | is the interelectronic distance. Without loss of general-ity, we choose the molecular axis along the z axis, andthe plane formed by the molecular axis and the polar-ization vector as the xz plane. Generally, we can de-compose the polarization vector into its two components, ǫ = cos θ N e z + sin θ N e x , where θ N is the angle betweenthe ς and ǫ axes, and e z and e x are the unit vectors alongthe z and x axes, respectively. The dipole interaction istherefore given as E ( t ) · ( r + r ) = E ( t ) h ( z + z ) cos θ N +( x + x ) sin θ N i , (4)where the rectangular coordinates x and z and the pro-late spheroidal coordinates are related through x = R p ( ξ − − η ) cos ϕ, z = R ξη. (5)We expand the wave function for the H molecule in thebody-frame asΨ(1 , , t ) = X m m Π m m ( ξ , η , ξ , η , t ) × Φ m m ( ϕ , ϕ ) . (6)Here Φ m m ( ϕ , ϕ ) = e i ( m ϕ + m ϕ ) / (2 π ) is the angularfunction, where m and m denote the magnetic quantumnumbers of the two electrons along the molecular axis.Next, Π m m ( ξ , η , ξ , η , t ) is expanded in a prod-uct of normalized “radial” { f i ( ξ ) } and “angular” { g k ( η ) } DVR bases:Π m m ( ξ , η , ξ , η , t ) = X ijkℓ f i ( ξ ) f j ( ξ ) g k ( η ) g ℓ ( η ) C m m ijkℓ ( t ) . (7)Note that the basis is not symmetrized with respect tothe coordinates of the two electrons. Since we alwaysbegin the calculation with a properly symmetrized initialstate, however, that symmetry will be preserved in thecalculation.To discretize this partial differential equation, we em-ploy the FE-DVR approach for both the ξ and η vari-ables [21, 24]. If we normalize the DVR bases accordingto Z dξf i ( ξ ) f i ′ ( ξ ) = δ ii ′ and Z dηg k ( η ) g k ′ ( η ) = δ kk ′ , (8)respectively, the overall ( ξ, η ) DVR basis is not normal-ized with respect to the volume element in the prolatecoordinate system. This is corrected by defining the two-electron basis b m m ijkℓ (1 ,
2) = (cid:16) R (cid:17) q ( ξ i − η k )( ξ j − η ℓ ) (9) × f i ( ξ ) f j ( ξ ) g k ( η ) g ℓ ( η )Φ m m ( ϕ , ϕ ) , which satisfies the desired normalization Z Z dV dV b m m ∗ ijkℓ (1 , b m ′ m ′ i ′ j ′ k ′ ℓ ′ (1 ,
2) = δ ii ′ δ jj ′ δ kk ′ δ ℓℓ ′ δ m m ′ δ m m ′ , (10)to expand Π m m ( ξ , η , ξ , η , t ). Specifically, we haveΨ(1 , , t ) = X m m X ijkℓ b m m ijkℓ (1 , X m m ijkℓ ( t ) . (11)Introducing the normalized ( ξ, η ) DVR basis eliminatesthe complexities of matrix operations related to the over-lap matrix at each time, and hence makes the standardSIL algorithm directly applicable to study the temporalresponse of the molecule to laser pulses. III. THE FE-DVR BASIS
In our current implementation of the FE-DVR ap-proach for the time-dependent wave function in prolate spheroidal coordinates, we have chosen to work directlyin the FE-DVR basis. This differs from what is usuallydone for atoms in spherical coordinates, where spheri-cal harmonics are used for the angular variables and anFE-DVR for the radial coordinates. Consequently, theboundary conditions in prolate spheroidal coordinatesrequire some more discussion. Analyzing the asymp-totics reveals that in the region near the boundaries of ξ = 1 and η = ±
1, which correspond to the molecu-lar axis, the single-electron wave function behaves like( ξ − | m | / (1 − η ) | m | / . This indicates that the physi-cal wave function is finite for | m | = 0, whereas it goes tozero in the region close to the molecular axis for | m | 6 = 0.More importantly, the behavior of the wave function forodd | m | contains a square-root factor, giving a decidedlynonpolynomial behavior to the wave function that is im-possible to capture in a straightforward fashion using aDVR basis.The former problem is readily treated by using aGauss-Radau quadrature in the first DVR element for ξ ,where only the right-most point is constrained to lie onthe boundary between the first and second finite element.The volume element ensures that the integrand is wellbehaved near the end points and makes it unnecessaryto invoke a separate quadrature for different m values.For all the other elements, a Gauss-Lobatto quadratureis employed. This allows us to make the FE-DVR basiscontinuous everywhere and to satisfy the | m | -dependentboundary condition.To overcome the nonanalytic behavior of the basis forodd m , Bachau and collaborators [17] explicitly factoredout the ( ξ − | m | / (1 − η ) | m | / part before the wavefunction was expanded in terms of B -splines in the dis-cretization approach. We have adopted a similar ideain our FE-DVR treatment of the two-center problem toachieve much faster convergence, as was also done inRef. [21]. For the case of even | m | , no changes need to bemade to define the DVR basis, i.e., the normalized basisis written as f i ( ξ ) = 1 q ω iξ Y k = i ξ − ξ k ξ i − ξ k and g i ( η ) = 1 q ω iη Y k = i η − η k η i − η k . (12)For odd | m | , however, we define the DVR basis as f i ( ξ ) = 1 q ω iξ ( ξ − / ( ξ i − / Y k = i ξ − ξ k ξ i − ξ k (13)and g i ( η ) = 1 q ω iη (1 − η ) / (1 − η i ) / Y k = i η − η k η i − η k . (14)Here ω iξ and ω iη are the weight factors related to the DVRbases f i ( ξ ) and g i ( η ), respectively. The goal of using aunique set of mesh points, which are | m | -independent, todiscretize the ( ξ, η ) coordinates has now been achievedin this scheme. The same technique was employed inrecent calculations of one- and two-photon double ion-ization of H [15, 24]. In principle, it is also possibleto introduce the factors ( ξ − | m | / (1 − η ) | m | / intothe DVR bases to circumvent the difficulties related tothe nonanalytic behavior near the boundary. However,this results in an | m | -dependence of the DVR bases andquadrature points. This, in turn, leads to a number ofunnecessary complications in the practical implementa-tion of the computational methodology. One might arguethat an | m | -dependent discretization procedure could beuseful for a system in which the magnetic quantum num-ber m is conserved. An example is the H +2 ion in externalmagnetic fields along the molecular axis [25]. However,that is not the situation in the current calculation. IV. THE ELECTRON-ELECTRON COULOMBINTERACTION IN PROLATE SPHEROIDALCOORDINATES
Similar to the expansion of the electron-electron inter-action in terms of spherical coordinates, a counterpartexists in prolate spheroidal coordinates [26] through theNeumann expansion1 r = 1 a ∞ X l =0 l X m = − l ( − | m | (2 l + 1) (cid:18) ( l − | m | )!( l + | m | )! (cid:19) (15) × P | m | l ( ξ < ) Q | m | l ( ξ > ) P | m | l ( η ) P | m | l ( η ) e im ( ϕ − ϕ ) , where a ≡ R/
2. The two nuclei are located at ± R/ z axis and ξ > ( < ) = max(min)( ξ , ξ ). Boththe regular P | m | l ( ξ ) and irregular Q | m | l ( ξ ) Legendre func-tions [27], which are defined in the region (1 , + ∞ ), areinvolved in the expansion as the “radial” part, while the“angular” part is only related to P | m | l ( η ). Note thatwe chose to work in terms of an un-normalized “an-gular” basis rather than the usual spherical harmonics.The matrix elements of 1 /r in a traditional basis, forexample, a B -spline or Slater-type basis, can be com-puted through the well-known Mehler-Ruedenberg trans-formation [19, 28]. Due to the discontinuous derivativealong the line of ξ = ξ in the Neumann expansion,the straightforward computation of the matrix elementof 1 /r , using the value of this interaction potential atthe mesh points, is very slowly convergent. We seek amore robust representation of the 1 /r matrix whichretains both the underlying Gauss quadrature and theDVR property of all potentials being exactly diagonalwith respect to the highly localized DVR basis.In the following, we use the simplified notation | ijkℓm m i = | f i ( ξ ) f j ( ξ ) g k ( η ) g ℓ ( η )Φ m m ( ϕ , ϕ ) i to denote the basis. Essentially, we need the integral D ijkℓm m (cid:12)(cid:12)(cid:12) P | m | l ( ξ < ) Q | m | l ( ξ > ) P | m | l ( η ) P | m | l ( η ) (16) × e im ( ϕ − ϕ ) (cid:12)(cid:12)(cid:12) i ′ k ′ j ′ ℓ ′ m ′ m ′ E . After integrating over ϕ and ϕ , the matrix element of1 /r can be written as D ijkℓm m | r | i ′ j ′ k ′ ℓ ′ m ′ m ′ E (17)= 1 a ∞ X l = | m | ( − | m | (2 l + 1) (cid:18) ( l − | m | )!( l + | m | )! (cid:19) I i ′ j ′ k ′ ℓijkℓ ( l ) , where the selection rule m = m − m ′ = m ′ − m hasbeen used. Hence m is uniquely determined for a givenpair of angular partial waves. Above we introduced thereduced ( ξ, η ) integral I i ′ j ′ k ′ ℓ ′ ijkℓ ( l ) (18)= (cid:10) ijkℓ (cid:12)(cid:12) P | m | l ( ξ < ) Q | m | l ( ξ > ) P | m | l ( η ) P | m | l ( η ) (cid:12)(cid:12) i ′ j ′ k ′ ℓ ′ (cid:11) . Since | m | is fixed in the above equation, we omitted it in I i ′ j ′ k ′ ℓ ′ ijkℓ ( l ) and will do so in the related quantities belowas well. It is now worthwhile to define the two electrondensities [29]: ρ A ( ξ, η ) = f i ( ξ ) g k ( η ) f i ′ ( ξ ) g k ′ ( η ) ,ρ B ( ξ, η ) = f j ( ξ ) g ℓ ( η ) f j ′ ( ξ ) g ℓ ′ ( η ) . (19)After truncating the radial integral to the edge of thebox, ξ max , this yields I i ′ j ′ k ′ ℓ ′ ijkℓ ( l ) = Z Z dV ξ dV ξ ′ ρ B ( ξ, η ) P | m | l ( ξ < ) Q | m | l ( ξ > ) × P | m | l ( η ) P | m | l ( η ′ ) ρ A ( ξ ′ , η ′ )= Z ξ max dV ξ P | m | l ( η ) ρ B ( ξ, η ) U l ( ξ ) . (20)Here a convention for the volume element was made insuch a way that, for any function F ( ξ, η ), we define dV ξ F ( ξ, η ) ≡ dξa R +1 − ( ξ − η ) F ( ξ, η ) dη to simplify thenotation. Most importantly, the function U l ( ξ ) is definedby U l ( ξ ) = Q | m | l ( ξ ) Z ξ dV ξ ′ ρ A ( ξ ′ , η ′ ) P | m | l ( ξ ′ ) P | m | l ( η ′ ) (21)+ P | m | l ( ξ ) Z ξ max ξ dV ξ ′ ρ A ( ξ ′ , η ′ ) Q | m | l ( ξ ′ ) P | m | l ( η ′ ) . Instead of evaluating the above integrals directly, wesolve the differential equation satisfied by U l ( ξ ). As willbecome apparent later, this equation can be shown to bethe “radial” Poisson equation in the prolate spheroidalcoordinate system. The differential equations satisfiedby the Legendre functions P | m | l ( ξ ) and Q | m | l ( ξ ) suggeststhat we introduce the operator ∇ ξ = ddξ ( ξ − ddξ − l ( l + 1) − m ξ − l and m . This is the one-dimensional Laplacian operator in the ξ coordinate.After some algebra, we obtain ddξ U l ( ξ ) = dQ | m | l ( ξ ) dξ Z ξ dV t ρ A ( t, τ ) P | m | l ( t ) P | m | l ( τ ) (23)+ dP | m | l ( ξ ) dξ Z ξ max ξ dV t ρ A ( t, τ ) Q | m | l ( t ) P | m | l ( τ ) . and d dξ U l ( ξ ) = d Q | m | l ( ξ ) dξ Z ξ dV t ρ A ( t, τ ) P | m | l ( t ) P | m | l ( τ )+ d P | m | l ( ξ ) dξ Z ξ dV t ρ A ( t, τ ) Q | m | l ( t ) P | m | l ( τ )++ W ( P | m | l ( ξ ) , Q | m | l ( ξ )) a Z +1 − dτ ( ξ − η ) ρ A ( ξ, τ ) P | m | l ( τ ) . (24)Here the Wronskian of the Legendre functions P | m | l ( ξ )and Q | m | l ( ξ ) is given by W ( P | m | l ( ξ ) , Q | m | l ( ξ )) = ( − | m | (1 − ξ ) ( l + | m | )!( l − | m | )! . (25)Consequently, we obtain ∇ ξ U l ( ξ ) = ̺ ( ξ ) . (26)This is the second-order inhomogeneous Poisson equationsatisfied by U l ( ξ ) with the “source” term given by ̺ ( ξ ) =( − | m | +1 ( l + | m | )!( l − | m | )! a (27) × Z +1 − dη ′ ( ξ − η ′ ) ρ A ( ξ, η ′ ) P | m | l ( η ′ ) . After carrying out the integral over η ′ via Gauss quadra-ture, we recast the source term as ̺ ( ξ ) = δ kk ′ ( − | m | +1 ( l + | m | )!( l − | m | )! a f i ( ξ ) f i ′ ( ξ ) × ( ξ − η k ) P | m | l ( η k ) . (28)As one might expect, the inhomogeneous equation re-duces to the homogeneous one if k = k ′ . The solutionto the Poisson equations (26) and (28) can be uniquelydetermined by enforcing the boundary conditions U l (1) = P | m | l (1) Z ξ max dV ξ ′ ρ A ( ξ ′ , η ′ ) Q | m | l ( ξ ′ ) P | m | l ( η ′ )(29)at ξ = 1 and U l ( ξ max ) = Q | m | l ( ξ max ) Z ξ max dV ξ ′ ρ A ( ξ ′ , η ) P | m | l ( ξ ′ ) P | m | l ( η ′ )(30) at ξ = ξ max , respectively.There are two important points to realize, namely:First, on the right-hand boundary, the function U l ( ξ ) as-sumes a nonzero value, which is given by Eq. (30), for allpossible | m | values. Its asymptotic behavior relies on thefunction Q | m | l ( ξ ) at large ξ , which behaves like 1 /ξ l +1max .This indicates that it is nonzero generally, although itcould be small at the large ξ max values used in practicalcalculations. Second, on the left-hand boundary, the sit-uation depends on the value of | m | . U l ( ξ ) takes a nonzero value if | m | = 0, while it becomes zero if | m | 6 = 0.Following the philosophy employed to handle thespherical case [29], we first seek a solution, U l ( ξ ), to thePoisson equations (26)-(30) that satisfies the zero -valueboundary condition at ξ max by using exactly the same ξ mesh points as those for the wave functions. In otherwords, we have ∇ ξ U l ( ξ ) = ̺ ( ξ ) with U l (1) = U l (1) and U l ( ξ max ) = 0. After substituting the DVR expansion U (0) l ( ξ ) = P µ c µ f µ ( ξ ) of the solution into the differentialequation, we obtain a system of linear equations for theunknown coefficients { c µ } : X µ ′ c µ ′ T | m | µµ ′ =( − | m | ( l + | m | )!( l − | m | )! 1 q ω iξ a δ µi δ ii ′ δ kk ′ × ( ξ i − η k ) P | m | l ( η k ) . (31)The matrix T is defined by its elements T | m | µµ ′ = − Z ξ max dξf µ ( ξ ) (cid:20) ( ξ − d dξ + 2 ξ ddξ (32) − l ( l + 1) − m ξ − (cid:21) f µ ′ ( ξ ) . Therefore, the coefficient c µ can formally be written as c µ = [ T | m | ] − µi q ω iξ ( − | m | ( l + | m | )!( l −| m | )! a δ ii ′ δ kk ′ ( ξ i − η k ) P | m | l ( η k ) , (33)where [ T | m | ] − denotes the inverse of the matrix T | m | .In this case U l ( ξ ) fulfills the left-hand boundary condi-tion U l (1) = 0, and so does U l (1). Recall, however, thatits right-hand boundary condition differs from those of U l ( ξ max ). This suggests that the final answer to the func-tion U l ( ξ ) can be constructed as U l ( ξ ) = U l ( ξ ) + F l ( ξ ),i.e., we add the difference function F l ( ξ ) to U l ( ξ ). Thefunction F l ( ξ ) is also a solution to the homogeneousPoisson equation, subject to the boundary condition F l (1) = 0 and F l ( ξ max ) = U l ( ξ max ). After writing it as alinear combination of P | m | l ( ξ ) and Q | m | l ( ξ ), and imposingthe boundary conditions, F l ( ξ ) takes the form F l ( ξ ) = δ ii ′ δ kk ′ a ( ξ i − η k ) P | m | l ( ξ i ) P | m | l ( η k ) × Q | m | l ( ξ max ) P | m | l ( ξ max ) P | m | l ( ξ ) . (34)We finally arrive at U l ( ξ ) = ( − | m | q ω iξ ( l + | m | )!( l − | m | )! a δ ii ′ δ kk ′ ( ξ i − η k ) P | m | l ( η k ) × X µ [ T | m | ] − µi f µ ( ξ ) + δ ii ′ δ kk ′ ( ξ i − η k ) P | m | l ( ξ i ) × P | m | l ( η k ) Q | m | l ( ξ max ) P | m | l ( ξ max ) P | m | l ( ξ ) . (35)At this point, the DVR version of the solution U l ( ξ ) isready for all possible values of | m | , either | m | 6 = 0 or | m | = 0. Substituting Eq. (35) into Eq. (20) allows us toobtain the kernel integral, I i ′ j ′ k ′ ℓ ′ ijkℓ ( l ) = δ ii ′ δ jj ′ δ kk ′ δ ℓℓ ′ a ( ξ i − η k )( ξ j − η ℓ ) P | m | l ( η ℓ ) × " ( − | m | q ω iξ ω jξ ( l + | m | )!( l − | m | )! [ T | m | ] − ji P | m | l ( η k )++ P | m | l ( ξ i ) P | m | l ( ξ j ) P | m | l ( η k ) Q | m | l ( ξ max ) P | m | l ( ξ max ) . (36)The matrix element of 1 /r can finally be written as (cid:10) ijkℓm m (cid:12)(cid:12) r (cid:12)(cid:12) i ′ j ′ k ′ ℓ ′ m ′ m ′ (cid:11) = δ ii ′ δ jj ′ δ kk ′ δ ℓℓ ′ a ( ξ i − η k )( ξ j − η ℓ ) l max X l > | m | (2 l + 1) ( l − | m | )!( l + | m | )! P | m | l ( η k ) P | m | l ( η ℓ ) × " q ω iξ ω jξ [ T | m | ] − ji + ( − | m | ( l − | m | )!( l + | m | )! P | m | l ( ξ i ) P | m | l ( ξ j ) Q | m | l ( ξ max ) P | m | l ( ξ max ) , (37)where we truncated the l summation in the Neumannexpansion to l max . The above equation can be convertedto the normalized ( ξ, η ) bases with the help of Eq. (9).This results in a diagonal representation of the matrix el-ements of the electron-electron Coulomb interaction andthus considerably simplifies the FE-DVR discretizationprocedure. The above treatment of the 1 /r matrix wassuccessfully applied to the two-photon double ionizationof H [24]. The implementation of this representationwill be illustrated below. V. TIME EVOLUTION AND EXTRACTION OFCROSS SECTIONS
The time-dependent laser-driven electronic wavepacket in the hydrogen molecule is obtained by solvingthe TDSE on the ( ξ, η ) grid. Launched from the pre-viously determined ground state, the time evolution ofthe system is achieved by using our recently developedSIL method [30, 31]. The ground state is determinedby relaxing the system in imaginary time from an ini-tial guess of the wave function on the grid points. Ateach time step we only need to generate the values of thediscretized wave function on the selected grid points. Ifdesired, the information at arbitrary points within thespatial box can be obtained from the interpolation pro-cedure in terms of the DVR bases.A few remarks seem appropriate regarding the efficientimplementation of the SIL algorithm. The highest en-ergy, E max , which essentially depends on the smallestseparation between the ( ξ, η ) mesh points and also on the maximum values of | m | and | m | , determines the largesttime step ∆ t for the propagation in real time. Typically, E max is about 6 ,
000 atomic units (a.u.) in our calcu-lations. Although the chances of electrons populatingstates with such high energies are practically negligiblefor short time scales of the laser-molecule interaction, wegenerally require ∆ t . π/E max in order to resolve themost rapid oscillations in the time evolution. This meansthat at least a few time steps are needed during one pe-riod of 2 π/E max . We refer readers to Refs. [30–32] forfurther details and discussions behind the SIL method.In order to ensure that the double-ionization wavepacket is sufficiently far away from the nuclei, and alsothat the two photoelectrons are well separated, we al-low the system to evolve for a few more cycles in thefield-free Hamiltonian, i.e., after the laser pulse has diedoff. This is the wave packet we use to extract the phys-ical information. The ionization probabilities and thecorresponding cross sections are extracted by projectingthe time-dependent wave packet onto uncorrelated two-electron continuum states satisfying the standard incom-ing boundary conditions. The latter states of H are con-structed from the one-electron continuum state of the H +2 ion described in the following subsection. A. Continuum states of H +2 The field-free wave function Φ( ξ, η, ϕ ) of the one-electron molecular ion is completely separable in pro-late spheroidal coordinates. For a given, and conserved,magnetic quantum number m , the wave function takesthe form Φ( ξ, η, ϕ ) = T m ( ξ )Ξ m ( η )Φ m ( ϕ ), where the az-imuthal dependence, is given by Φ m ( ϕ ) ≡ e imϕ / √ π .The “radial” part T mq ( ξ ) and the “angular” part Ξ mq ( η )of the wave function satisfy the equations (cid:20) ∂∂ξ ( ξ − ∂∂ξ − m ( ξ −
1) + 2 Rξ + c ξ − A mq (cid:21) T mq ( ξ ) = 0(38)and (cid:20) ∂∂η (1 − η ) ∂∂η − m (1 − η ) − c η + A mq (cid:21) Ξ mq ( η ) = 0 , (39)respectively. Here c = kR/ k . In addi-tion, we need to introduce another quantum number q ,which denotes the number of nodes of Ξ m ( η ) in the re-gion η ∈ [ − , +1], to label the states, and finally theseparation constant A mq .When the angular function Ξ m ( η )Φ m ( ϕ ) is discretizedin terms of the relevant DVR bases, a few “spurious” so-lutions might be encountered. This is caused by the resid-ual errors associated with the Gauss quadratures. Conse-quently, we expand the angular function, or “spheroidalharmonics” function Y ℓm ( η, ϕ ) ≡ Ξ mq ( η )Φ mq ( ϕ ) with ℓ = | m | + q instead in terms of spherical harmonics. Thesefunctions are normalized according to Z +1 − dη Z π dϕ Y ∗ ℓm ( η, ϕ ) Y ℓ ′ m ′ ( η, ϕ ) = δ mm ′ δ ℓℓ ′ . (40)After obtaining the separation constant A mq by solv-ing Eq. (39), the “radial” function T mq ( ξ ) is again ex-panded in terms of the DVR bases. The last DVR pointat ξ = ξ max needs to be kept for the continuum state.Asymptotically, the radial function behaves like T mq ( ξ ) → ξR r π sin (cid:20) cξ + Rc ln(2 cξ ) − ℓπ mq ( k ) (cid:21) (41)as ξ → + ∞ . Here ∆ mq ( k ) is the two-center Coulombphase shift. The normalization factor on either the en-ergy or the momentum scale and the Coulomb phase shiftcan be determined by matching the numerical solutionof T mq ( ξ ) according to its asymptotic behavior given inEq. (41).The plane wave in prolate spheroidal coordinates canbe written as [33] e i k · r = 4 π X ℓm i ℓ Y ℓm ( η r , ϕ r ) Y ∗ ℓm ( η k , ϕ k ) R ( k ) ℓm ( ξ ) , (42)where R ( k ) ℓm ( ξ ) → / ( cξ ) sin (cid:2) cξ − ℓπ/ (cid:3) in the asymp-totic region. Note that η k and η r are related to thedirections of k and r in spherical coordinates through η k,r = cos θ k,r . The partial-wave expansion of the planewave e i k · r reminds us that the two-center Coulomb wave satisfying the incoming boundary condition can be ex-panded asΦ ( − ) k ( r ) = 1 k + ∞ X m = −∞ X ℓ > | m | i ℓ e − i ∆ mq ( k ) × Y ∗ ℓm ( k ) Y ℓm ( η r , ϕ r ) T ( k ) mq ( ξ ) . (43)This function is normalized in momentum space accord-ing to h Φ ( − ) k | Φ ( − ) k ′ i = δ ( k − k ′ ), provided the asymptoticsolution in Eq. (41) is satisfied.Uncorrelated two-electron continuum states with to-tal spin angular momentum S ( S = 0 in our case) cangenerally be constructed as follows:Φ ( − ) k k ( r , r ) = (44)1 √ h Φ ( − ) k ( r )Φ ( − ) k ( r ) + ( − S Φ ( − ) k ( r )Φ ( − ) k ( r ) i . With the help of Eq. (43), its partial-wave representationcan be written asΦ ( − ) k k ( r , r ) = (cid:16) R (cid:17) √ k k X ℓ m ℓ m i ℓ + ℓ × X ijkℓ b m m ijkℓ (1 , q ( ξ i − η k )( ξ j − η ℓ ) (45) × " e − i (cid:0) ∆ | m | ℓ ( k )+∆ | m | ℓ ( k ) (cid:1) Y ∗ ℓ m ( k ) Y ∗ ℓ m ( k ) C ℓ m ℓ m ijkℓ ( k , k ) + ( − S ( k ↔ k ) . Here we introduced C ℓ m ℓ m ijkℓ ( k , k ) = (46)˜ T ( k ) ℓ | m | ( ξ i ) ˜ T ( k ) ℓ | m | ( ξ j )˜Ξ ( k ) ℓ | m | ( η k )˜Ξ ( k ) ℓ | m | ( η ℓ ) , by representing the radial and angular parts on the ( ξ, η )grid points: T ( k ) ℓm ( ξ ) = X i f i ( ξ ) ˜ T ( k ) ℓm ( ξ i ) , (47)Ξ ( k ) ℓm ( η ) = X µ g µ ( η )˜Ξ ( k ) ℓm ( η µ ) . (48)Here the exchange symmetry C ℓ m ℓ m jiℓk ( k , k ) = C ℓ m ℓ m ijkℓ ( k , k ) (49)is satisfied. B. Extraction of double-ionization cross sections
It has been demonstrated [24, 30, 34] that using uncor-related two-electron continuum states is a good approxi-mation in a time-dependent propagation approach, pro-vided the two ejected electrons are well separated fromeach other. The probability amplitude of double ioniza-tion is then given by h Φ ( − ) k k | Ψ( t ) i = (50)1 k k X m ℓ m ℓ ( − i ) ℓ + ℓ e i (cid:0) ∆ | m | ℓ ( k )+∆ | m | ℓ ( k ) (cid:1) × Y ℓ m ( k ) Y ℓ m ( k ) F ℓ m ℓ m ( k , k ) , where F ℓ m ℓ m ( k , k ) = (51) √ X ijkℓ C ℓ m ℓ m ∗ ijkℓ ( k , k ) X m m ijkℓ ( t ) . Here the exchange symmetry F ℓ m ℓ m ( k , k ) = ( − S F ℓ m ℓ m ( k , k ) (52)is satisfied. We also see that the probability ampli-tude formulated in prolate spheroidal coordinates takesa similar form as for the atomic case in spherical coordi-nates. However, a subtle difference from the atomic caseis worth pointing out. Strictly speaking, the spheroidalharmonics involved in the probability amplitude gener-ally depend on the magnitude of the momenta k and k , in addition to their directions.The energy sharing of the two photoelectrons canbe specified by introducing the hyperangle α =tan − ( k /k ). This describes the double-ionization re-action with kinetic energies E = E exc cos α and E = E exc sin α for the two ionized electrons, respectively.Here E exc is the available excess energy above the double-ionization threshold. In the present work, specifically, E exc = 23 . k and ˆ k canbe extracted using the same formalism as in the corre-sponding He case [30, 34], i.e.,d σ d E dˆ k dˆ k = 1 k k cos α ωI T (1)eff Z k max d k ′ Z k ′ d k ′ × k ′ δ ( k ′ − k ′ tan α ) (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) X m ℓ m ℓ ( − i ) ℓ + ℓ e i (cid:0) ∆ | m | ℓ +∆ | m | ℓ (cid:1) × Y ℓ m ( k ′ , ˆ k ) Y ℓ m ( k ′ , ˆ k ) F ℓ m ℓ m ( k ′ , k ′ ) (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) . (53)Here ω and I are the central photon energy and thepeak intensity of the laser pulse, respectively, while T (1)eff denotes the effective interaction time between the tem-poral electric laser field and the electrons in the one-photon absorption process. For a laser pulse of time du-ration τ with a sine-squared envelope for the field ampli-tude, T (1)eff = (3 / τ . Note that T (1)eff corresponds to the special case of the generalized N -photon effective inter-action time T ( N )eff [35] for a one-photon reaction. Gener-alized cross sections for two-photon double ionization ofthe hydrogen molecule were extracted in the same for-malism [24].In the TDCC treatment [13], different strategies wereemployed to describe the linear one-photon and the non-linear two-photon double ionization processes of atomsand molecules. For the one-photon case, the cross sec-tions were obtained through the time derivative of thedouble-ionization probability, ∂P ( t ) /∂t . The laserfield does not need to be turned off in this case. On theother hand, a true laser pulse was used for the two-photoncase and an effective time, defined as the time integralunder a flat-top pulse with a smooth turn-on and turn-off, was introduced. In the present work, we employed aunified formulation through an effective interaction timefor both one-photon and multi-photon ionization in laserpulses.For the one-photon double ionization initialized fromthe lowest X Σ g state, the two ejected electrons can onlypopulate the final Σ u and Π u continuum states, withthe specifics depending on the relative orientation of themolecular axis and the laser polarization vector. Con-sequently only partial waves with ungerade parity [i.e.,( − ℓ + ℓ = −
1] need to be included in Eq. (53).
VI. RESULTSA. Preparation of the initial electronic X Σ g state For the nonsequential double-ionization process in-duced by one- or two-photon absorption, electronic cor-relation plays a dominant role, as the two photoelectronsmust share the available excess energy E exc . Double ion-ization by a single photon would not occur at all if thetwo-electron atom or molecule were approximated by anindependent-electron model. Therefore, the quality ofthe description of electron-electron correlation in a laser-driven system is crucially important for accurate resultsto be obtained. The Coulomb interaction between thetwo electrons has to be described in a consistent man-ner for both the initial bound state and the time-evolvedwave packet. Before we go any further, it is worth dis-cussing how we prepare the initial X Σ g state at theequilibrium distance of R = 1 . /r is uniquely determined by the angular bases. How-ever, this is not the case for the index l , if we choose todiscretize the coordinate η , rather then expanding thatpart of the wave function into spherical harmonics. Inpractice, the summation over l must be truncated at afinite value of l max . In principle, the higher-order ex-pansion terms always guarantee well-converged results.However, as mentioned earlier, we approximate the rel-evant η -integrals by using Gauss-Legendre quadrature.This is the price we have to pay for making the dielec-tronic Coulomb potential diagonal in the DVR bases. Asa consequence, we need to determine how the approxi-mation introduced in the η -integrals for the two-electronintegrals affects the results for the cross sections of inter-est.To answer this question, we first investigate the de-pendence of the energy obtained for the initial X Σ g state on the value of l max used in the Neumann expan-sion. Figure 1 shows the variation of the initial-stateelectronic energy of the hydrogen molecule with respectto l max , obtained with a ξ setup of ten elements in theregion of 1 < ξ .
82, five in the region 1 < ξ ξ .
82. Each element, inturn, contains five DVR points to further discretize theconfiguration space. Furthermore, we employ 9th-orderDVR points for η . For a given number of η mesh points( n η ) and | m | max = | m | max = | m | max , we observe thatthe resulting energy typically exhibits a plateau-like be-havior with increasing l max . For given n η , when l max is relatively small, the η -integral can be computed veryaccurately by using Gauss quadrature. However, when l max is too large, the numerical errors introduced fromthe Gauss quadrature cause the energy value to fluctu-ate. This occurs when l max approaches 2 n η and is shownby the grey stripes in Fig. 1. In this region of l max anunphysically low energy can be produced. Beyond thatpoint, the calculated energy increases to the next plateau.Ultimately, this is not too surprising, since any Gaussquadrature is only reasonably accurate up to a limitedpolynomial order of the integrand. Consequently, if wewant to keep more terms in the Neumann expansion,we have to increase n η correspondingly. This finding isfurther substantiated by the dependence of the energyfound for n η = 11 and 13. The plateaus are indeed ex-tended to the correspondingly larger values of 2 n η . Mostimportantly, the amplitude of the energy fluctuation issystematically reduced with increasing n η . The error inthe energy is lowered from 2 . × − to 1 . × − and finally 1 . × − a.u., when n η increases from 9to 11 and then 13 for | m | max = 4. We obtained the en-ergy at R = 1 . − . l max = 10, n η = 9, and | m | max = 4, resulting in a double-ionizationpotential of 51 .
394 eV. Keeping the other parameters un-changed, we obtained an energy of − . n η = 11. The benchmark energy in the literature is − . R [36], after we take outthe nucleus-nucleus interaction of 1 / . l max employed in the Neu-mann expansion in practical calculations, if we discretizethe coordinate η . However, this provides a way to ex-amine a potential sensitivity of the physical observablesof interest (here the differential cross sections) to theground-state wave functions generated by varying l max and other parameters. This will be further discussed be-low. (13 , , , E gs = − . , , , l max E l ec tr o n i ce n e r g y ( a u ) FIG. 1. (Color online) Energy of the lowest electronic X Σ g state at R = 1 . l max value used inthe Neumann expansion of 1 /r . The number of η points, n η , and the largest magnetic quantum number, | m | max , arelabeled as ( n η , | m | max ). The open symbols correspond to | m | max = 1, while the filled symbols are for | m | max = 4. Thebenchmark energy ( E gs ) from Ref. [36] is shown as well. B. Convergence of the TDCS
Before we present our results for the cross sections, letus take a closer look at the survival probability P surv = |h Ψ gs | Ψ( t ) i| (54)of the aligned H molecule in its ground state Ψ gs . Thisis shown in Fig. 2. For homonuclear molecules, the inde-pendent alignment angle θ N between the molecular axisand the polarization vector can be confined to the regionfrom 0 ◦ to 90 ◦ . In the xuv regime, we observe that thehydrogen molecule shows a larger probability of beingionized or excited (i.e., a lower probability of staying inthe initial state) at the end of the pulse in an alignedgeometry. This indicates that the perpendicular compo-nent of the temporal electric field exerts more influenceon the ionization process due to the larger dipole mo-mentum. Interestingly, at the earlier stages of the timeevolution (e.g., t . P surv . However, once the electric field has be-come sufficiently strong ( t & P surv for the tilted molecule.When the wave packet is driven back to the nuclearregion and therefore has a chance to recombine with theH +2 ion, a maximum in P surv appears. Not surprisingly,the parallel geometry always has the largest probabilityfor this to happen. Although the wave packet can alsobe scattered for the untilted molecule in the plane per-pendicular to the molecular axis, the probability is un-doubtedly larger if the laser electric field is perpendicular0to the molecular axis. A similar behavior of H +2 in xuvpulses was observed in Ref. [37]. θ N = 90 ◦ θ N = 60 ◦ θ N = 30 ◦ θ N = 0 ◦ ( a ) Time (au) Su r v i v a l p r o b a b ili t y b )Time (au) Su r v i v a l p r o b a b ili t y FIG. 2. (Color online) Survival probability of the hydrogenmolecule subjected to a sine-squared laser pulse with a peakintensity of 10 W/cm . The laser pulse lasts for 10 opticalcycles and the system is followed for a period of another 2 cy-cles of field-free propagation. The central photon energy ofthe laser pulse is 75 eV. For most calculations performed in this study, we ex-pose the hydrogen molecule to a laser pulse with a peakintensity of 10 W/cm . Looking at Fig. 2 we see thatthe depletion of the initial ground state can be safelyneglected for our typical interaction times. Even for θ N = 90 ◦ , P surv = 0 . W/cm mightseem very intense for most atomic and molecular targets.Here, however, we consider an xuv rather than an IRpulse. For an xuv pulse with central photon energy of75 eV, such laser fields definitely fall into the “weak-field”regime. The ponderomotive energy in the xuv regime ismuch smaller than the photon energy of interest.In this work, we are mainly interested in the triple- TABLE I. The discretization and expansion parameters of theH wave function in prolate spheroidal coordinates. Here ξ b stands for the border between the inner and outer regionsin the ξ coordinate, while ξ max is the size of the ξ box. Inaddition, n ξ denotes the number of ξ mesh points in eachelement. The numbers of ξ elements in the inner and outerregion are n inn and n out , respectively. These ξ parametersproduce the total number of ξ mesh points N ξ . The ξ grid Iand ξ grid II are used to examine the convergence of ourresults. ξ b n inn n out ξ max n ξ N ξ ξ grid I 5 5 67 150 5 288 ξ grid II 9 1 11 97 14 156 differential cross section, since it reveals the fine detailsof possible energy sharings and preferred directions of theejected electrons in the double-ionization process. Giventhe discrepancies between results from various calcula-tions found in the literature, we carried out comprehen-sive convergence tests for our predictions of the TDCSs.These tests are essentially divided into two groups. Thefirst group concerns the laser parameters, while the sec-ond one deals with the discretization and expansion pa-rameters. An example of two different parameter sets forthe ξ grid is given in Table I and will be further discussedbelow.In order to obtain a good handle on the sensitivity ofthe results to the various parameters and the resultinglevel of “convergence”, we try to only vary a single pa-rameter while keeping all others fixed if possible. Forthe dependence on the laser parameters, we use the ξ grid I combined with ( n η , | m | max , l max ) = (9 , , W/cm anda time scale of “10 + 2” optical cycles (o.c.) was used.Here “10 + 2” refers to a 10-cycle laser pulse with a sine-squared envelope for the field amplitude, followed by a2-cycle field-free propagation.Figures 3 and 4 show the convergence pattern of ourTDCS results for asymmetric energy sharing in the par-allel geometry ( θ N = 0 ◦ ). The energy sharing betweenelectron 1 (observed at the fixed angle θ ) and electron 2(observed at the variable angle θ ) is 20% : 80%. Onlythe electron that takes away 20% of the excess energy isrecorded at fixed positions either parallel or perpendicu-lar to the polarization axis.Since the laser pulse is explicitly involved in our time-dependent treatment, we first have to be sure that theextracted cross sections are essentially independent of thelaser intensity and the time scales. Only then are the cal-culations of cross sections meaningful. This also allowsus to compare the physical information extracted fromour time-dependent scenario to that obtained throughconventional time-independent treatments, which are ef-fectively equivalent to the weak-field approximation and“infinitely” long interaction times. Rather than comput-1 W/cm W/cm ( a ) θ = 0 ◦ θ N = 0 ◦ θ (deg) T D C S ( b / s r e V ) b ) θ = 90 ◦ θ N = 0 ◦ θ (deg) T D C S ( b / s r e V ) c ) θ = 0 ◦ θ N = 0 ◦ θ (deg) T D C S ( b / s r e V ) d ) θ = 90 ◦ θ N = 0 ◦ θ (deg) T D C S ( b / s r e V ) FIG. 3. (Color online) Convergence of the coplanar TDCSresults for the hydrogen molecule for asymmetric energy shar-ing with respect to the laser peak intensity and the time scale.The central photon energy is 75 eV. The slow reference elec-tron, observed at the fixed angle θ , takes away 20% of theavailable excess energy ( E = 4 . E exc ( E = 18 . b )-( d ) is 10 W/cm . The two columns show thecorresponding convergence of the TDCS for θ = 0 ◦ (left) and θ = 90 ◦ (right), respectively. 1 barn (b) = 10 − cm . ing the cross sections, it would be more appropriate toconsider ionization rates if the cross sections were foundto be sensitive to the laser parameters.In Fig. 3, we display the dependences of our TDCS re-sults upon the laser parameters. Note that the TDCSsextracted from I = 10 W/cm and 10 W/cm atfixed time evolution of “10 + 2” cycles are nearly identi-cal and agree with each other to better than the thick-ness of the line. When we turn to the dependence oftime scales at a fixed intensity of 10 W/cm , we usethe same pulse, but allow the system to freely evolve fora few additional cycles to extract the TDCS. This corre-sponds to the time scale of “10 + 4” o.c. Also, we mayincrease the laser-molecule interaction time, but extractthe TDCSs at the same cycles of field-free time evolu-tion after the pulse died off. This gives the scenario of“12 + 2” o.c. Since the total time durations are the same(14 o.c.), they allow us to examine the extracted TDCSsfrom different perspectives. The increased interactiontime yields a reduced bandwidth of the photon energy,while the longer field-free propagation ensures that thedouble-ionization wave packet is further away from thenuclear region [38]. The calculated TDCSs indeed showa slight, though in our opinion acceptable, sensitivity tothe time scales. Not surprisingly, the sensitivity is mostvisible for the smaller cross sections, when the two ejectedelectrons travel nearly parallel along the same direction(c.f. Fig. 3( c )).Having confidence in using the current sets of laserparameters, we now turn our attention to the scheme of n ξ = 7 n ξ = 5( a ) θ = 0 ◦ θ N = 0 ◦ θ (deg) T D C S ( b / s r e V ) b ) θ = 90 ◦ θ N = 0 ◦ θ (deg) T D C S ( b / s r e V ) ξ grid II ξ grid I( c ) θ = 0 ◦ θ N = 0 ◦ θ (deg) T D C S ( b / s r e V ) d ) θ = 90 ◦ θ N = 0 ◦ θ (deg) T D C S ( b / s r e V ) ξ max = 150 ξ max = 100( e ) θ = 0 ◦ θ N = 0 ◦ θ (deg) T D C S ( b / s r e V ) f ) θ = 90 ◦ θ N = 0 ◦ θ (deg) T D C S ( b / s r e V ) n η = 11 n η = 9( g ) θ = 0 ◦ θ N = 0 ◦ θ (deg) T D C S ( b / s r e V ) h ) θ = 90 ◦ θ N = 0 ◦ θ (deg) T D C S ( b / s r e V ) l max = 20 l max = 14 l max = 12 l max = 10( i ) θ = 0 ◦ θ N = 0 ◦ θ (deg) T D C S ( b / s r e V ) j ) θ = 90 ◦ θ N = 0 ◦ θ (deg) T D C S ( b / s r e V ) | m | max = 5 | m | max = 4( k ) θ = 0 ◦ θ N = 0 ◦ θ (deg) T D C S ( b / s r e V ) l ) θ = 90 ◦ θ N = 0 ◦ θ (deg) T D C S ( b / s r e V ) FIG. 4. (Color online) Same as Fig. 3, but for the convergenceof the TDCS results with respect to the discretization andexpansion parameters. See text for the details. spatial discretization ( n ξ , ξ max , n η ) and the convergenceof the expansion ( l max , | m | max ). The results are displayedin Fig. 4.For the discretization parameters, we obtain well-converged TDCSs by increasing n ξ from 5 to 7, n η from 9to 11, and extending the spatial box of ξ max from 100 to2 spherical TDCCspherical ECSprolate ECSpresent work ( a ) θ = 0 ◦ θ N = 90 ◦ θ (deg) T D C S ( b / s r e V ) b ) θ = 0 ◦ θ N = 60 ◦ θ (deg) T D C S ( b / s r e V ) c ) θ = 0 ◦ θ N = 30 ◦ θ (deg) T D C S ( b / s r e V ) ( d ) θ = 0 ◦ θ N = 0 ◦ θ (deg) T D C S ( b / s r e V ) FIG. 5. (Color online) Coplanar TDCS of the aligned hydrogen molecule for equal energy sharing ( E = E = 11 . θ = 0 ◦ with respect to the laser polarizationaxis. Also shown are the one-center spherical ECS results [11], the two-center prolate spheroidal results [15], and one-centerspherical TDCC results [14]. ξ mesh points: ξ grid I and ξ grid II (see Table I). Theprincipal motivation was to see whether or not we canreproduce the much lower TDCS values (by about 20%compared to the one-center spherical results) that wererecently obtained in an ECS calculation in two-center el-liptical coordinates by Tao et al. [15].We emphasize that these two grids in the “radial” ξ co-ordinate are completely different regarding both the dis-tribution of the elements and the number of grid pointsper element. In the ξ grid I, we divide the ξ space intotwo parts, an inner and an outer region with a border at ξ b = 5. We place a narrow span of elements in the innerregion, and then wider elements in the outer region. Incontrast to that, ξ grid II does not distinguish betweeninner and outer regions, i.e., the elements uniformly spanthe region from 1 to ξ max . The mesh setup in ξ grid II isthe same as that used in Ref. [15], except for the complexrotation. The ξ grid I has a much denser distribution ofmesh points than ξ grid II. Nevertheless, the extractedTDCSs from both sets of ξ grids are in excellent agree-ment with each other, even for the smallest cross sections. This strongly suggests that the results are well convergedat least with regard to the ξ grid. Both ξ sets are goodenough to capture the physics of interest. Differences atthe 20% level are unlikely to be caused by using differentsets of ξ meshes.Finally, we discuss the convergence of our results withrespect to the expansion parameters, | m | max and l max .As expected for a one-photon process, | m | max = 4 pro-duced well-converged results.Recall the discussion above regarding the ground state,especially how the truncated Neumann expansion of1 /r in our present FE-DVR implementation affects theinitial-state energy and therefore the quality of the wavefunction. For consistency, we use the same l max in thereal-time propagation and in the ground-state wave func-tion. As seen from Figs. 4( i ) and 4( j ), our truncated Neu-mann expansion has little effect on the calculated TDCSvalues. Well-converged results can be obtained even withan inappropriately large value of l max = 20, which yieldsa slightly higher energy of the ground state (c.f. Fig. 1).Overall, our detailed convergence tests only reveal avery weak sensitivity of the TDCS results to both the3 spherical TDCCspherical ECSprolate ECSpresent work ( a ) θ = 90 ◦ θ N = 90 ◦ θ (deg) T D C S ( b / s r e V ) b ) θ = 90 ◦ θ N = 60 ◦ θ (deg) T D C S ( b / s r e V ) c ) θ = 90 ◦ θ N = 30 ◦ θ (deg) T D C S ( b / s r e V ) ( d ) θ = 90 ◦ θ N = 0 ◦ θ (deg) T D C S ( b / s r e V ) FIG. 6. (Color online) Same as Fig. 5, except that the fixed electron is detected at the angle θ = 90 ◦ with respect to thelaser polarization axis. Since there was a plotting error in Fig. 3 of Tao et al. [15], we are comparing here with the propernumbers [39] from that calculation. time scales and the values of l max . Well-converged TDCSresults can be obtained by using either ξ grid I or ξ grid IIcombined with ( n η , | m | max , l max ) = (9 , , ξ grid I to discretize the two-electron wave packet and a “10 + 2” sine-squared laserpulse with a peak intensity of 10 W/cm . C. TDCSs for the aligned H molecule Figures 5, 6, and 7 display the coplanar TDCSs ofthe aligned hydrogen molecule at equal and asymmet-ric ( E : E = 20% : 80%) energy sharing. The twoelectrons are detected in the same (coplanar) plane de-fined by the ζ and ǫ axes. The angles θ , θ , and θ N are all measured with respect to the laser linear polariza-tion axis. We compare our TDCS predictions with thoseobtained in the time-independent one-center sphericalECS calculation [11], the time-independent two-centerspheroidal ECS model [15], and the time-dependent one-center spherical TDCC approach [14]. The TDCC num-bers were recently recalculated with a bigger box size and differ, in some cases substantially, from those pub-lished originally [13]. Except for the recent two-centerprolate spheroidal ECS results of Tao et al. [15, 39], theagreement between the other three sets of results is verysatisfactory. Once again, the largest relative differencesoccur when the cross sections are small (see Figs. 5( d )and 6( d )).Using spheroidal coordinates as well, as an illustrativeexample of their two-center ECS approach, Serov andJoulakian [40] recently presented the TDCS at the samephoton energy, but only for a single geometry of θ N = 20 ◦ and θ = 40 ◦ for asymmetric energy sharing of E : E =80% : 20%. Although not shown here, there is againgood agreement between their results, Vanroose et al. ’sone-center spherical ECS numbers [11], and our time-dependent FE-DVR predictions.It is also interesting to investigate the dominant es-cape modes for the various scenarios. These modes arestrongly dependent on how the electrons share the excessenergy. In an arbitrary geometry (0 ◦ θ N ◦ ), forexample, the back-to-back escape mode ( θ = 180 ◦ ) isforbidden for equal energy sharing. On the other hand, itbecomes the dominant mode for significantly asymmetric4 spherical ECSprolate ECSpresent work ( a ) θ = 90 ◦ θ N = 90 ◦ θ (deg) T D C S ( b / s r e V ) ( b ) θ = 90 ◦ θ N = 60 ◦ θ (deg) T D C S ( b / s r e V ) c ) θ = 90 ◦ θ N = 30 ◦ θ (deg) T D C S ( b / s r e V ) ( d ) θ = 90 ◦ θ N = 0 ◦ θ (deg) T D C S ( b / s r e V ) FIG. 7. (Color online) Coplanar TDCS of the aligned hydrogen molecule for asymmetric energy sharing. The electron detectedat the fixed angle θ = 90 ◦ takes away 20% of the available excess energy, while the second electron takes away 80% of E exc .The present time-dependent FE-DVR results are compared with those from time-independent one-center spherical ECS [11]and two-center prolate spheroidal ECS [15] calculations. k ( a ) k ( b ) k ( c ) k ( d ) FIG. 8. (Color online) Comparison of predicted relative coplanar TDCSs between H in the perpendicular (solid lines) andparallel (dashed lines) geometries, and He (chain lines) [42] for equal-energy sharing in polar coordinates. The polarizationaxis is taken along the horizontal direction. The photon energies for H and He are 75 eV and 99 eV, respectively. The fixedobservation angles for one of the electrons are 0 ◦ (a), 30 ◦ (b), 60 ◦ (c), and 90 ◦ (d) with respect to the laser polarization vector.Scaling factors were used to emphasize the shape comparison. energy sharing, including the 20%:80% scenario discussedin the present paper (see Fig. 3).These results can be understood from a symmetryanalysis [41]. Equal-energy sharing and back-to-backemission is equivalent to k = − k . When we con-sider the exchange and parity operations simultaneously in Eq. (45), we have Φ − k , − k = P ( − S Φ k , k . Here P = ± ungerade parity, we therefore must haveΦ − k , k ( r , r ) = 0 at any configuration of r and r . Al-though the magnitudes of the momenta k ′ and k ′ are5 spherical TDCCpresent work k out of plane ( a ) θ = 90 ◦ θ N = 90 ◦ θ (deg) T D C S ( b / s r e V ) k out of plane ( b ) θ = 90 ◦ θ N = 90 ◦ θ (deg) T D C S ( b / s r e V ) k out of plane ( c ) θ = 90 ◦ θ N = 0 ◦ θ (deg) T D C S ( b / s r e V ) FIG. 9. (Color online) Noncoplanar TDCS of the aligned hydrogen molecule for equal energy sharing. The present time-dependent FE-DVR results are compared with TDCC predictions [14]. not exactly conserved in the time-dependent picture, theionization events we collect must satisfy the condition k ′ = k ′ (because of the δ function in Eq. (53)) for theequal-energy sharing. This is the reason behind the for-bidden back-to-back ( θ = 180 ◦ ) escape mode for theequal energy sharing, as we observed in Figs. 5 and 6.Since the argument does not involve the relative align-ment angles, it is valid for all possible values of θ N . Onthe other hand, this is not the case when the excess en-ergy is not evenly distributed among the two electrons.Indeed, Figs. 3( a ) and 3( c ) show maxima in the back-to-back emission, thereby illustrating the dramatic changein the dominant escape mode.For equal-energy sharing in the parallel geometry( θ N = 0 ◦ ), the electron-electron Coulomb repulsion sug-gests that the TDCS should be dynamically small if thetwo electrons travel along the same direction. This isin agreement with the numerically small cross sections(not exactly zero, however) at θ = 0 ◦ or 360 ◦ seen inFig. 5( d ).Recall that the one-photon double-photoionizationprocess in helium [34, 42] shares the same property. Theback-to-back mode is forbidden for equal-energy sharing,and this can be explained by the above argument. Itis one of the similarities between the molecular hydrogenand the atomic helium targets for double photoionization.However, Figs. 5, 6, and 7 also reveal significant molecu-lar effects in the TDCS results. These are missing for the helium atom, not only in the shape of the angular distri-butions, but also in the magnitudes of the cross sections.Depending on the relative orientation (0 ◦ < θ N < ◦ ),there is interference between the Σ u and Π u symmetriesin H . A nice example of this effect was presented byReddish et al. [43]. Even without interference (i.e., for θ N = 0 ◦ or 90 ◦ ), the perpendicular geometry shows muchlarger magnitudes of the TDCS than the parallel geom-etry. Figure 8 shows the three cases of angular distribu-tions: H ζ ⊥ ǫ , H ζ k ǫ , and He at equal energy shar-ing. Interestingly, in most cases the angular distributionsof the perpendicular geometry resemble those of helium.The molecular effect can definitely not be ignored in theparallel geometry for θ = 0 ◦ (c.f. Fig. 8( a )). The for-ward escape mode of the second electron is dominant forthe H parallel geometry. In contrast, the backward modeis dominant for the H perpendicular geometry and alsofor helium.In Fig. 9, we show the TDCS for noncoplanar ge-ometries. Again, all angles are defined with respect tothe polarization vector. For the perpendicular geometry,Fig. 9( a ) depicts the escape modes for the configurationof k k ζ (the fixed electron) and at the same time k inthe plane perpendicular to the plane formed by ǫ and ζ .Figure 9( b ) shows the TDCS after exchanging the direc-tions of k and k in Fig. 9( a ). With the same directionsof k and k as in Fig. 9( b ), Fig. 9( c ) is for the case ofthe molecular axis orientated along the polarization vec-6tor. In the parallel case ( θ N = 0 ◦ ), we observe that anyescape modes of both electrons ejected in the directionperpendicular to ǫ are forbidden. This can be under-stood by analyzing the spheroidal harmonics in Eq. (53).In this case, only the Σ u states can be populated.Hence, only partial waves with ( m , m ) = ( − m, m )and ( − ℓ + ℓ can contribute to the cross sections, since Y ℓ m ( k ′ , ˆ k ) Y ℓ m ( k ′ , ˆ k ) at the angles of θ = θ = 90 ◦ vanish in spherical coordinates. Once again, the agree-ment between our FE-DVR noncoplanar TDCSs and therefined TDCC results [14] is excellent. VII. SUMMARY
We have presented calculations for one-photon dou-ble ionization of the hydrogen molecule at a photon en-ergy of 75 eV by solving the time-dependent Schr¨odingerequation in prolate spheroidal coordinates. The triple-differential cross sections were extracted through the pro-jection of the time-dependent wave packet onto uncorre-lated two-electron continuum states, a few cycles of field-free time evolution after the laser pulse died off.Exhaustive convergence studies of the TDCS resultswere performed with respect to a number of discretizationand expansion parameters, as well as the details of thelaser field. These tests provide a strong indication thatthe results for the triple-differential cross sections pre-sented here are well converged and numerically accurate.Excellent agreement was obtained between the current time-dependent results in prolate spheroidal coordinates,those obtained with the ECS approach in spherical co-ordinates [11] and, finally, larger TDCC calculations [14]than those published earlier [13].The present calculations do not confirm the signifi-cant reduction by about 20% in the TDCS results pre-dicted in recent ECS calculations in the two-center pro-late spheroidal coordinates [15]. Furthermore, our resultsdid not show the level of sensitivity to the description ofthe ground state that was also reported by Tao et al. [15].The detailed analysis reported in this study provides ahigh level of confidence in the present results. We hopethat they will be used as benchmarks for comparison infuture investigations. Tables of these results are availablein electronic format from the authors upon request.
ACKNOWLEDGMENTS
We thank Drs. T. N. Rescigno and J. Colgan for send-ing their results in numerical form and for helpful dis-cussions. This work was supported by the NSF undergrant PHY-0757755 (XG and KB) and generous super-computer resources through the NSF TeraGrid allocationaward TG-PHY090031 (Kraken at NICS, Oak Ridge Na-tional Laboratory) and by the Department of Energy al-location award MPH006 (Jaguar at NCCS, Oak RidgeNational Laboratory). Without these computational re-sources it would have been impossible to study the con-vergence properties of our cross section results on variousphysical and computational parameters. [1] H. Br¨auning, R. D¨orner, C. L. Cocke, M. H. Prior, B.Kr¨assig, A. S. Kheifets, I. Bray, A. Br¨auning-Demian, K.Carnes, S. Dreuil, V. Mergel, P. Richard, J. Ullrich, andH. Schmidt-B¨ocking, J. Phys. B , 5149 (1998).[2] Th. Weber, A. Czasch, O. Jagutzki, A. M¨uller, V. Mergel,A. Kheifets, J. Feagin, E. Rotenberg, G. Meigs, M. H.Prior, S. Daveau, A. L. Landers, C. L. Cocke, T. Osipov,H. Schmidt-B¨ocking, and R. D¨orner, Phys. Rev. Lett. , 437 (2004).[4] Th. Weber, PhD. Thesis, Universit¨at Frankfurt (unpub-lished, 2003).[5] M. Gisselbrecht, M. Lavoll´ee, A. Huetz, P. Bolognesi, L.Avaldi, D. P. Seccombe, and T. J. Reddish, Phys. Rev.Lett. , 153002 (2006).[6] Y. H. Jiang, A. Rudenko, E. Pl´esiat, L. Foucar, M.Kurka, K. U. K¨uhnel, Th. Ergler, J. F. P´erez-Torres, F.Mart´ın, O. Herrwerth, M. Lezius, M. F. Kling, J. Titze,T. Jahnke, R. D¨orner, J. L. Sanz-Vicario, M. Sch¨offler,J. van Tilborg, A. Belkacem, K. Ueda, T. J. M. Zouros,S. D¨usterer, R. Treusch, C. D. Schr¨oter, R. Moshammer, and J. Ullrich, Phys. Rev. A , 021401(R) (2010).[7] J. P. Wightman, S. Cvejanovi´c, and T. J. Reddish, J.Phys. B , 1753 (1998).[8] A. S. Kheifets, Phys. Rev. A , 022704 (2005).[9] A. S. Kheifets and I. Bray, Phys. Rev. A , 022703(2005).[10] W. Vanroose, F. Mart´ın, T. N. Rescigno, and C. W. Mc-Curdy, Phys. Rev. A , 050703(R) (2004).[11] W. Vanroose, D. A. Horner, F. Mart´ın, T. N. Rescigno,and C. W. McCurdy, Phys. Rev. A , 052702 (2006).[12] W. Vanroose, F. Mart´ın, T. N. Rescigno, and C. W. Mc-Curdy, Science , 1787 (2006).[13] J. Colgan, M. S. Pindzola, and F. Robicheaux, Phys. Rev.Lett. , 153001 (2007).[14] J. Colgan, private communication (2010).[15] L. Tao, C. W. McCurdy, and T. N. Rescigno, Phys. Rev.A , 023423 (2010).[16] D. R. Bates, U. ¨Opik, and G. Poots, Proc. Phys. Soc. A , 1113 (1953).[17] S. Barmaki, S. Laulan, H. Bachau, and M. Ghalim, J.Phys. B , 817 (2003).[18] S. Barmaki, H. Bachau, and M. Ghalim, Phys. Rev. A , 043403 (2004).[19] Y. V. Vanne and A. Saenz, J. Phys. B , 4101 (2004). [20] G. Lagmago Kamta and A. D. Bandrauk, Phys. Rev. A , 053407 (2005).[21] L. Tao, C. W. McCurdy, and T. N. Rescigno, Phys. Rev.A , 012719 (2009).[22] T. J. Park and J. C. Light, J. Chem. Phys. , 5870(1986).[23] D. J. Tannor, in Introduction to Quantum Mechanics, ATime-dependent Perspective , Chap. 11, p318 (UniversityScience Books, Sausalito, California 2007).[24] X. Guan, K. Bartschat, and B. I. Schneider, Phys. Rev.A , 041404(R) (2010).[25] X. Guan, B. Li, and K. T. Taylor, J. Phys. B , 3569(2003).[26] P. M. Morse and H. Feshbach, in Methods of TheoreticalPhysics , Parts I and II (McGraw-Hill, 1953).[27] I. A. Stegun, in
Handbook of Mathematical Functions withFormulas, Graphs, Mathematical Tables , Chap. 8, p331,eds. M. Abramowitz and I. A. Stegun (1964).[28] E. L. Mehler and K. Ruedenberg, J. Chem. Phys. ,2575 (1969).[29] C. W. McCurdy, M. Baertschy, and T. N. Rescigno, J.Phys. B , R137 (2004).[30] X. Guan, K. Bartschat, and B. I. Schneider, Phys. Rev.A , 043421 (2008).[31] X. Guan, C. J. Noble, O. Zatsarinny, K. Bartschat, andB. I. Schneider, Comp. Phys. Comm. , 2401 (2009).[32] K. T. Taylor, J. S. Parker, D. Dundas, K. J. Meharg,L. R. Moore, E. S. Smyth, and J. F. McCann, in Many- Particle Quantum Dynamics in Atomic and MolecularFragmentation , Chap. 9, p153, eds. J. Ullrich and V.Shevelko (Springer-Verlag, Berlin Heidelberg 2003).[33] C. Flammer, in
Spheroidal Wave Functions , (Dover Pub-lications, Inc., Mineola, New York 2005).[34] J. Colgan, M. S. Pindzola, and Robicheaux, J. Phys. B , L457 (2001).[35] L. A. A. Nikolopoulos and P. Lambropoulos, Phys. Rev.A , 063410 (2006).[36] J. Sims and S. Hagstrom, J. Chem. Phys. , 094101(2006).[37] S. X. Hu, L. A. Collins, and B. I. Schneider, Phys. Rev.A , 023426 (2009).[38] L. B. Madsen, L. A. A. Nikolopoulos, T. K. Kjeldsen,and J. Fern´andez, Phys. Rev. A , 063407 (2007).[39] T. N. Rescigno, private communication (2010).[40] V. V. Serov and B. B. Joulakian, Phys. Rev. A , 062713(2009).[41] F. Maulbetsch and J. S. Briggs, J. Phys. B , 551(1995).[42] X. Guan, unpublished (2008).[43] T. J. Reddish, J. Colgan, P. Bolognesi, L. Avaldi, M.Gisselbrecht, M. Lavoll´ee, M. S. Pindzola, and A. Huetz,Phys. Rev. Lett.100