A quantum simulation of dissociative ionization of H + 2 in full dimensionality with time dependent surface flux method
AA quantum simulation of dissociative ionization of H +2 in full dimensionality withtime-dependent surface flux method Jinzhen Zhu ∗ Physics Department, Ludwig Maximilians Universit¨at, D-80333 Munich, Germany
The dissociative ionization of H +2 in a linearly polarized, 400 nm laser pulse is simulated by solvinga three-particle time-dependent Schr¨odinger equation in full dimensionality without using any datafrom quantum chemistry computation. The joint energy spectrum (JES) is computed using a time-dependent surface flux (tSurff) method, the details of which are given. The calculated groundenergy is -0.597 atomic units and internuclear distance is 1.997 atomic units if the kinetic energyterm of protons is excluded, consistent with the reported precise values from quantum chemistrycomputation. If the kinetic term of the protons is included, the ground energy is -0.592 atomic unitswith an internuclear distance 2.05 atomic units. Energy sharing is observed in JES and we findpeak of the JES with respect to nuclear kinetic energy release (KER) is within 2 ∼ PACS numbers: 32.80.-t,32.80.Rm,32.80.Fb
I. INTRODUCTION
Understanding the three-body Coulomb interactionproblem is an on-going challenge in attosecond physics.The typical candidates for investigation include Heliumatom and H +2 molecule. In attosecond experiments, ashort, intense laser pulse is introduced as a probe for themeasurements. Various mechanisms were proposed inthe recent decades to describe the dissociation and disso-ciative ionization of H +2 , including bond softening [1],thecharge-resonance enhanced ionization (CREI) [2], bondhardening [3], above threshold dissociation (ATD) [4, 5],high-order-harmonic generation (HHG) [6] and abovethreshold explosion [7]. One may find a summary of theabove mechanisms in theoretical and experimental inves-tigations of H +2 in literature [8, 9]. Experimental studieson the H +2 ion exposed to circular and linearly polarizedpulses for angular and energy distributions of electronswere reported recently [10–13].In theory, the joint energy spectra (JES) of the kineticenergy release (KER) for one electron and two protonsof the H +2 ion are predominant observables that showhow energy distributes around the fragments, where theJES is represented by the KER of two electrons for dou-ble ionization (DI) [14–16]. In theory, the JES computa-tions for double ionization in full dimensionality was veryscarce for laser pulses with wavelengths beyond the XUVregime ( ≥
400 nm) because the computational consump-tion scales dramatically with the wavelength and inten-sity of the laser field [16]. With tSurff method, which wasfirst introduced in Ref. [17], full dimensional simulationof the JES for double ionization was available with mod-erate computational resources for 800 nm [16] and 400nm [18] laser pulses. The tSurff method was also suc- ∗ [email protected] cessfully applied to the dissociative ionization of the H +2 ion [19, 20] in a two-dimensional (2D) model, where theenergy sharing of the photons and electron is observed inJES.The dissociative ionization of the H +2 ion has been sim-ulated by many groups [10, 14, 21–26]. However, they areall in reduced dimensionality. Quantum simulation in fulldimensionality is not available yet. Although the corre-lation among the fragments could be observed in the 2Dmodel, the peaks of the JES with total nuclear KER arealways above 10 eV. This is far from experimental observ-ables [11–13], which are usually below 5 eV. The tRecXcode, which successfully implements the tSurff methodin full dimensionality, has been applied successfully inthe simulations of the double ionization of Helium [16]and the single ionization of polyelectron molecules [27–31]. The dissociative ionization of the H +2 ion has notbeen computed using the tRecX code from before, evenin reduced dimensionality.In this paper, we will introduce simulations of the dis-sociative ionization of the H +2 ion by solving the time-dependent Schr¨odinger equation (TDSE) in full dimen-sionality based on the tRecX code. We will first presentthe computational method for scattering amplitudes withtSurff methods, from which the JES can be obtained.Then we will introduce the specific numerical recipes forthe H +2 ion based on the existing discretization methodsof tRecX code. With such numerical implementations,the ab initio calculation of field free ground energy of theHamiltonian is available. Finally we will present resultsof dissociative ionization in a 400 nm laser pulse, the JES,and projected energy spectrum on the azimuth angle. II. METHODS
In this paper, atomic units with specifying (cid:126) = e = m e = 4 π(cid:15) ≡ a r X i v : . [ phy s i c s . a t o m - ph ] S e p mass of two protons is chosen to be the origin. Insteadof using the vector between two protons (cid:126)R as an coordi-nate [14, 19, 20], we specify the coordinates of the protonsand electrons as (cid:126)r , − (cid:126)r and (cid:126)r . We denote M = 1836atomic units as the mass of the proton. A. Hamiltonian
The total Hamiltonian can be represented by sumof the electron-proton interaction H EP and two tensorproducts, written as H = H B = H (+) ⊗ + ⊗ H ( − ) − H EP , (1)where the tensor products are formed by the identity op-erator multiplied by the Hamiltonian for two protons( H (+) ), or that for the electron ( H ( − ) ). H B is called theHamiltonian in the B region and will be detailed later.With the coordinate transformation used in Ref. [32],which is also illustrated in Appendix A for our specificcase, the single operator for the electron is H ( − ) = − ∆2 m − i β (cid:126)A ( t ) · (cid:126) (cid:79) , (2)and the Hamiltonian for protons can be written as H (+) = − ∆4 M + 12 r , (3)where we introduce reduced mass m = M M +1 ≈ β = MM ≈ (cid:126)A ( t ) is the vectorpotential. The Hamiltonian of the electron-proton inter-action can be written as H EP = 1 | (cid:126)r + (cid:126)r | + 1 | (cid:126)r − (cid:126)r | . (4) B. tSurff for dissociative ionization
The tSurff method is applied here for the dissocia-tive ionizations, which was successfully applied to thepolyelectron molecules and to the double emission of Heatom [16, 28–31]. In this section, we will follow a similarprocedure as is done in Ref. [16].According to the approximations of tSurff method, be-yond a sufficient large tSurff radius R (+ / − ) c , the interac-tions of protons and electrons can be neglected, with thecorresponding Hamiltonians being H (+) V = − ∆4 M for twoprotons and H ( − ) V = − ∆2 m − i β (cid:126)A ( t ) · (cid:126) (cid:79) for the electron.The scattered states of the two protons, which satisfyi ∂ t χ (cid:126)k ( (cid:126)r ) = H (+) V χ (cid:126)k ( (cid:126)r ), are χ (cid:126)k ( (cid:126)r ) = 1(2 π ) / exp( − i (cid:90) tt k M dτ ) exp( − i (cid:126)k (cid:126)r ) , (5) FIG. 1. The regions of dissociative ionization time propaga-tion. The B stands for bound region, D for dissociation regionwhere the two protons are out of R (+) c but electron not ion-ized and stays inside. I represents the ionization region whereelectron is out-of-box R ( − ) c but two protons are still inside R (+) c . DI stands for the dissociative ionization region whereboth the electron and the protons are out of R (+ / − ) c . R (+ / − ) c are the tSurff radii for r = | (cid:126)r | or r = | (cid:126)r | . and those of the electron, which satisfies i ∂ t χ (cid:126)k ( (cid:126)r ) = H ( − ) V χ (cid:126)k ( (cid:126)r ), are χ (cid:126)k ( (cid:126)r ) = 1(2 π ) / exp( − i (cid:90) tt k m − i β (cid:126)A ( τ ) · (cid:126) (cid:79) dτ ) exp( − i (cid:126)k (cid:126)r ) , (6)where we assume the laser field starts at t and (cid:126)k / denote the momenta of the protons or the electron.Based on the tSurff radius R (+ / − ) c , we may splitthe dissociative ionization into four regions namely B, I, D, DI , shown in figure 1, where bound region B preserves the full Hamiltonian in Eq. (1), D, I are timepropagations by single particles with the Hamiltonian H D ( (cid:126)r , t ) = H ( − ) V ( (cid:126)r , t ) = − ∆2 m − i β (cid:126)A ( t ) · (cid:126) (cid:79) (7)and H I ( (cid:126)r , t ) = − ∆4 M + 12 r , (8)and DI is an integration process. The treatment was firstintroduced in the double ionization of Helium in Ref. [15]and then applied in a 2D simulation of the H +2 ion inRef. [19]. Without considering the low-energy free elec-trons that stay inside the box after time propagation, wemay write ψ B ( (cid:126)r , (cid:126)r , t ) ≈ , r ≥ R (+) c , or r ≥ R ( − ) c ψ D ( (cid:126)r , (cid:126)r , t ) ≈ , r < R (+) c , or r ≥ R ( − ) c ψ I ( (cid:126)r , (cid:126)r , t ) ≈ , r ≥ R (+) c , or r < R ( − ) c ψ DI ( (cid:126)r , (cid:126)r , t ) ≈ , r < R (+) c , or r < R ( − ) c (9)We assume that for a sufficiently long propagation time T , the scattering ansatz of electron and protons disentan-gle. By introducing the step functionΘ / ( R c ) = (cid:40) , r / < R (+ / − ) c , r / ≥ R (+ / − ) c , (10)the unbound spectra can be written as P ( (cid:126)k , (cid:126)k ) = P ( φ , θ , k , φ , θ , k ) = (cid:12)(cid:12)(cid:12) b ( (cid:126)k , (cid:126)k , T ) (cid:12)(cid:12)(cid:12) . (11) b ( (cid:126)k , (cid:126)k , T ) is the scattering amplitudes and can be writ-ten as b ( (cid:126)k , (cid:126)k , T ) = (cid:104) χ (cid:126)k ⊗ χ (cid:126)k | Θ ( R c )Θ ( R c ) | ψ ( (cid:126)r , (cid:126)r , t ) (cid:105) = (cid:90) T −∞ [ F ( (cid:126)k , (cid:126)k , t ) + ¯ F ( (cid:126)k , (cid:126)k , t )] dt (12)with two sources written as F ( (cid:126)k , (cid:126)k , t ) = (cid:104) χ (cid:126)k ( (cid:126)r , t ) (cid:12)(cid:12)(cid:12) [ H ( − ) V ( (cid:126)r , t ) , Θ ( R c )] (cid:12)(cid:12)(cid:12) ϕ (cid:126)k ( (cid:126)r , t ) (cid:105) (13)and¯ F ( (cid:126)k , (cid:126)k , t ) = (cid:104) χ (cid:126)k ( (cid:126)r , t ) (cid:12)(cid:12)(cid:12) [ H (+) V ( (cid:126)r , t ) , Θ ( R c )] (cid:12)(cid:12)(cid:12) ϕ (cid:126)k ( (cid:126)r , t ) (cid:105) . (14)The single particle wavefunctions ϕ (cid:126)k ( (cid:126)r , t ) and ϕ (cid:126)k ( (cid:126)r , t )satisfyi ddt ϕ (cid:126)k ( (cid:126)r , t ) = H D ( (cid:126)r , t ) ϕ (cid:126)k ( (cid:126)r , t ) − C (cid:126)k ( (cid:126)r , t ) (15)andi ddt ϕ (cid:126)k ( (cid:126)r , t ) = H I ( (cid:126)r , t ) ϕ (cid:126)k ( (cid:126)r , t ) − C (cid:126)k ( (cid:126)r , t ) . (16)The sources are the overlaps of the two-electron wave-function on the Volkov solutions shown by C (cid:126)k ( (cid:126)r , t ) = (cid:90) d (cid:126)r χ (cid:126)k ( (cid:126)r , t )[ H (+) V ( (cid:126)r , t ) , Θ ( R c )] ψ ( (cid:126)r , (cid:126)r , t )(17)and C (cid:126)k ( (cid:126)r , t ) = (cid:90) d (cid:126)r χ (cid:126)k ( (cid:126)r , t )[ H ( − ) V ( (cid:126)r , t ) , Θ ( R c )] ψ ( (cid:126)r , (cid:126)r , t ) , (18)with initial values being 0, where · · · means complexconjugate. The two tSurff radii could be set to equiv-alent R (+) c = R ( − ) c , because all Coulomb interactions areneglected when either the protons or electron is out ofthe tSurff radius. According to our previous researches,the spectrum computation is independent of the R c ifall Coulomb terms are removed and the wavefunction ispropagated long enough after the pulse [15, 16]. ThetSurff for double emission of two particles was firstly in-troduced in Ref. [15]. The above derivations are very similar to what was reported in Ref. [16] of double emis-sion of Helium, where the only differences are constantsbefore different operators, say ∆ , (cid:126) (cid:79) and r . Thus, detailedformulas are omitted here and the interested readers canrefer to Ref. [15, 16].The computation for photoelectron spectrum includesfour steps, similar to the one used in Ref. [16], detailedas 1. Solve full 6D TDSE with the Hamiltonian in the B region, given in Eq. (1), and write the time-dependent surface values in the disk.2. Evolve the single-particle wave packets in the D region by Eq. (13) with surface values given in the B region time propagation.3. Evolve the single-particle wave packets in I regionby Eq. (14) with the surface values given in the B region time propagation.4. Integrate the fluxes calculated from surface valueswritten in the D and I regions’ time propagationby Eq. (12). III. NUMERICAL IMPLEMENTATIONS
The numerical methods here are similar to what wasdetailed in Ref. [16, 18]. In fact, the code in this paperis developed based on the double ionization frameworkof the tRecX code used in the reference [16, 18]. Thus,we will focus on the electron-protons interaction whichwas not mentioned about before and only list relevantdiscretization methods in this paper.
A. Discretization and basis functions
The 6D wavefunction ψ is represented by the productof spherical harmonics for angular momentum and radialfunctions as ψ ( (cid:126)r , (cid:126)r , t ) = ψ ( r , θ , φ , r , θ , φ , t )= (cid:88) m ,l ,m ,l Y m l ( θ , φ ) Y m l ( θ , φ ) R m ,m ,l ,l ( r , r , t ) , (19)where Y m l ( θ , φ ) and Y m l ( θ , φ ) are the spherical har-monics of the two electrons and the radial function isrepresented by the finite-element discretization variablerepresentation (FE-DVR) method as R m ,m ,l ,l ( r , r , t ) = (cid:88) n ,n R n ,n m ,m ,l ,l ( r , r , t ) R n ,n m ,m ,l ,l ( r , r , t ) = (cid:88) p ,p f ( n ) p ( r ) f ( n ) p ( r ) 1 r r c m ,m ,l ,l n ,n ,p ,p ( t ) (20)where f ( n / ) p / ( r / ) are p / th basis functions on n / thelement, and the time-dependency of the three parti-cles are included in the radial functions and coefficients c m ,m ,l ,l n ,n ,p ,p ( t ), as is used in Ref. [15, 16]. The infinite-range exterior complex scaling (irECS) method is utilizedas an absorber [33]. The tSurff expression for computingspectra of such discretization can be found in Ref. [16]. B. Electron-protons interaction
The first part of electron-protons interaction can bewritten in a multi-pole expansion as1 | (cid:126)r − (cid:126)r | = 1 (cid:112) r + r − r r cos γ = 1 r > (cid:112) h − h cos γ = ∞ (cid:88) l =0 h l r > P l (cos γ ) , (21)where r > = max( r , r ) , r < = min( r , r ) , h = r < r > , γ is the angle between (cid:126)r , (cid:126)r and P l (cos γ ) are Legendrepolynomials. Similarly, we have1 | (cid:126)r + (cid:126)r | = 1 r > (cid:112) h + 2 h cos γ = ∞ (cid:88) l =0 ( − l h l r > P l (cos γ ) . (22)And the summation goes as1 | (cid:126)r + (cid:126)r | + 1 | (cid:126)r − (cid:126)r | = 2 ∞ (cid:88) l =0 h l r > P l (cos γ ) l mod 2 = 0 . (23)where l mod 2 = 0 means l is even. With the Legendrepolynomials expanded by spherical harmonics Y ml ( θ , φ )and Y m ∗ l ( θ , φ ) , we have H EP = 2 ∞ (cid:88) l =0 l (cid:88) m = − l π l + 1 r l< r l +1 > Y ml ( θ , φ ) Y m ∗ l ( θ , φ ) l mod 2 = 0 . (24)The matrix elements of electron-protons are (cid:104) ψ ( n (cid:48) ,n (cid:48) ) m (cid:48) ,m (cid:48) ,l (cid:48) ,l (cid:48) | | (cid:126)r − (cid:126)r | + 1 | (cid:126)r + (cid:126)r | | ψ ( n ,n ) m ,m ,l ,l (cid:105) =2 (cid:88) λµ π λ + 1 (cid:104) Y m (cid:48) l (cid:48) Y µλ | Y m l (cid:105)(cid:104) Y m (cid:48) l (cid:48) | Y µλ Y m l (cid:105)(cid:104) R n (cid:48) ,n (cid:48) m (cid:48) ,m (cid:48) ,l (cid:48) ,l (cid:48) | r λ< r λ +1 > | R n ,n m ,m ,l ,l (cid:105) , λ mod 2 = 0 , (25)which could be obtained by dropping the odd λ termsand doubling the even λ terms in the standard multi-poleexpansion for electron electron interactions from Ref. [16] as (cid:104) ψ ( n (cid:48) n (cid:48) ) m (cid:48) m (cid:48) l (cid:48) l (cid:48) | | (cid:126)r − (cid:126)r | | ψ ( n n ) m ,m ,l ,l (cid:105) = (cid:88) λµ π λ + 1 (cid:104) Y m (cid:48) l (cid:48) Y µλ | Y m l (cid:105)(cid:104) Y m (cid:48) l (cid:48) | Y µλ Y m l (cid:105)(cid:104) R n (cid:48) n (cid:48) m (cid:48) ,m (cid:48) ,l (cid:48) ,l (cid:48) | r λ< r λ +1 > | R n n m ,l ,m ,l (cid:105) . (26)Here ψ ( n ,n ) m ,m ,l ,l = Y m l ( θ , φ ) Y m l ( θ , φ ) R n ,n m ,m ,l ,l ( r , r , t ) . (27)Therein, the matrix for electron-protons interactioncould be obtained by the numerical recipes used inRef. [16, 34] with limited changes. Numerically, we find λ does not need to go to infinity and a maximum valueof 8 already suffices our simulations. IV. NUMERICAL RESULTS
A numerical convergence study shows, unlike the 6Ddouble emission of He, where m / = 0 , ≤ l / ≤ ≤ m / ≤ ≤ l / ≤ +2 ion. The R (+) c = R ( − ) c =12 . R ( − ) c does not change the quality of the spectrum butintroduces longer propagation time for low-energy par-ticles to fly out. R (+) c = 12 . R = 25 atomic units as used inRef. [19]. According to the convergence study in ap-pendix B, R (+) c = R ( − ) c = 12 . R ( − ) c = 12 . . × W/cm laser pulse. Ref. [35] shows theCoulomb explosion at large distances contributes to lowenergy fragments of protons, which are highly correlatedto the resonances between two H +2 eigenstates. Thus weneed a large simulation box to include the all possibleeigenenergies. In Ref. [35], the maximum internucleardistance for the low energy fragments is 11 atomic units.The molecular eigenenergies are nearly invariant withinternuclear distance far below our 2 R (+) c = 25 atomicunits here. Any potential dynamic dissociation quench-ing effect (DDQ) is also included in B region for the H +2 is in dissociative limit with internuclear distance over 12atomic units [36]. The wavefunction is propagated longenough after the pulse for to include the unbound stateslow kinetic energies.If the kinetic energy of protons is included, the fieldfree ground energy value is E = − .
592 atomic units andthe internuclear distance is 2.05 atomic units. With thekinetic energy of protons excluded, the ground eigenen-ergy is -0.597 atomic units, three digits exact to groundenergy from quantum chemistry calculations in Ref. [37],where the internuclear distance is fixed. The internucleardistance is 1.997 atomic units, three digits exact to thatfrom the precise computations in Ref. [38].
A. Laser pulses
The dipole field of a laser pulse with peak intensity I = E (atomic units) and linear polarization in z -directionis defined as E z ( t ) = − ∂ t A z ( t ), phase φ CEP = 0 with A z ( t ) = E ω a ( t ) sin( ωt + φ CEP ) . (28)A pulse with λ = 400 nm is given with intensities8 . × W/cm close to 2D computation in refer-ence [19] and 5 . × W/cm close to experimentalconditions in Ref. [12]. We choose a ( t ) = [cos( t/T )] as a realistic envelope. Pulse durations are specified asFWHM=5 opt.cyc. w.r.t. intensity. To compare withthe published results, a 400 nm, sin envelope laser pulseat 8 . × W/cm from Ref. [19] and a 791 nm, cos laser pulse at 7 . × W/cm used in Ref. [39] are alsoapplied. B. Joint energy spectra
The JES of the two dissociative protons and the elec-tron is obtained by integrating Eq. (11) over angular co-ordinates as σ ( E N , E e ) = (cid:90) dφ (cid:90) dφ (cid:90) dθ sin θ (cid:90) dθ sin θ P ( φ , θ , (cid:112) M × E N , φ , θ , (cid:112) m × E e ) , (29)where E N , E e are kinetic energies of two protons and anelectron, respectively. σ ( E N , E e ) is presented in Fig. 2 (a,b). The tilt lines with formula E N + E e = N ω + E − U p with ponderomotive energy U p = A m specify the energysharing of N photons for both the computations from8 . × W/cm and 5 . × W/cm , indicating corre-lated emissions of the electron and protons, which is alsoobserved in the 2D computations [19, 20]. The yields areintense around nuclear KER from 2 eV to 4 eV in the cos envelope pulse, consistent with the experimental valuesreported in Ref. [12, 13]. The peak of JES for dissocia-tive ionization is for lower nuclear KER than that (3-4eV) from Coulomb explosion from ground eigenstate ofthe H +2 ion, which property is also close to experimen-tal observables [39]. The Coulomb explosion JES is ob-tained with the same method as dissociative ionizationexcept that H EP is removed from B region Hamiltonianas H ( CS ) B = H (+) ⊗ + ⊗ H ( − ) , but the initial stateis still obtained from Hamiltonian H B in Eq. (1). Wefind that the contribution from time-propagation in sub-region D → DI (see Eq. (13)) is small, as the numerical FIG. 2. Log-scale JES log σ ( E N , E e ) represented by totalenergy of two protons E N and that of an electron E e . Linearpolarized, 400 nm, with (a) cos envelope with FWHM=5opt.cyc. at 8 . × W/cm and (b) cos envelope withFWHM=5 opt.cyc. pulses at 5 . × W/cm is appliedto the H +2 ion. The dashed lines represent the energy sharingbetween the protons and electron with formula E N + E e = Nω + E − U p , where ω is the photon energy. (c) JES fromCoulomb explosion simulation from the ground eigenstate ofthe H +2 ion. (d) Log-scale error log( δ ( σ )) of two spectrafrom cos envelope laser pulse at 8 . × W/cm with andwithout the contribution from D → DI (from Eq. (13)) by δ ( σ ) = 2 | σ (cid:48) ( E N ,E e ) − σ ( E N ,E e ) || σ (cid:48) ( E N ,E e )+ σ ( E N ,E e ) | . σ ( E N , E e ) of (a) and (b) arenormalized by dividing the maximum value. error of JES δ ( σ ) of σ computed from I → DI , and σ (cid:48) computed from two subregions ( I → DI and D → DI ),is always below 1% the main contribution of the JES(2 < E N < eV ), see Fig. 2 (d). This numerical propertyis also observed in two-dimensional (2D) simulations [19].This is because the electrons are much faster than pro-tons and the H +2 ion tends to release first.Here we compare our computations with the publishedresults. First, the JES with a 400 nm, sin envelopelaser pulse at 8 . × W/cm as used in the 2D simu-lations in Ref. [19] is computed, see Fig. 3 (a), comparedto the results from 2D simulations in Fig. 3 (b). Oneclearly sees the JES is most considerable with nuclearKER around 2-4 eV in our computation but is with nu-clear KER above 10 eV in the 2D simulation. We alsoattach the JES with a linearly polarized, 791 nm laserpulse at 7 . × W/cm as used in Ref. [39], whereJES is considerable with nuclear KER around 3 eV. Forthe 791 nm computation, R ( − ) c = 15 atomic units is ap-plied which is slightly above the quiver radius of theelectron. For a direct comparison, we integrate the JESaround the electron KER and obtain the photoelectronspectrum with respect to the nuclear KER in Fig. 3 (d),where the experimental data from Ref. [39] is also at-tached. The peak of the spectrum around 1.42 eV in our FIG. 3. Log-scale JES log σ ( E N , E e ) represented by totalenergy of two protons E N and that of an electron E e . JESof H +2 in a linearly polarized, 400 nm laser pulse at 8 . × W/cm as is used in Ref. [19] by (a) 6D computation com-pared to (b) the normalized 2D data from Ref. [19], and (c) ina linearly polarized, 791 nm laser pulse at 7 . × W/cm asused in Ref. [39]. (d) Blue dots: The normalized single JES iscreated from σ S ( E N, / ) = (cid:82) σ (2 E N, / , E e ) √ E e dE e , where E N, / depicts the kinetic energy of a proton. Green triangles:The normalized signals extracted from Ref. [39] with respectto the kinetic energy of a proton. computation is in good agreement with the experimentaldata, and the position of a minor peak around 1.7 eVin our computation also matches the experimental obser-vation. The observation that JES is most considerablewith nuclear KER around 2-4 eV can also be found inthe Coulomb explosion computation shown in Fig. 2 (c).In other experiments, the distribution of emitted protonspeaks at nuclear KER=4 eV for a 780 nm laser pulse at6 × W/cm [11], Ref. [12, 13] reported the nuclearKER is most probable around 3 eV for two protons for400 nm laser pulses. These observables at different exper-imental conditions show nuclear KERs are around 2 ∼ R = 2 r = 2 atomic units,which gives nuclear KER 0.5 atomic units. We now con-sider the wave packet dispersion in different coordinates.For 1D simulation on protons with a symmetric Gaussianwave packet, the center of the wave packet moves withthe velocity of the classical particle, thus most probablenuclear KER is 0.5 atomic units. In 3D simulation onprotons, the wave packet keeps expanding perpendicularto the polarization direction and is not symmetric withany r >
0, see Fig. 6. A numerical study with a simpleunphysical toy model in Sec. C shows for the Coulomb explosion in a spherical coordinate, the most probable nu-clear KER can be shifted to 1/4-1/3 of the total kineticenergy and the right half curve of the integrated JES isless steep, see Fig. 6 (a). The longer tail in higher nuclearKER for integrated JES of H +2 is also observed both incomputation and experiments, but not in 1D simulationon protons, see Fig. 3 (d) and Fig. 7. The existing2D simulations for the dissociative ionization put correc-tions to the electron-proton interaction with a softeningparameter to give the correct ground energy of electrons H +2 [14, 19]. However, the pure Coulomb repulsion ofthe two protons 1 /R ( R is the internuclear distance) isincluded without a softening parameter. We would liketo point out that, for the 2D simulation, for consistencyof the correction of Coulomb interaction of the electron,the Coulomb repulsion term of the two protons may alsoneed a softening parameter, whose value needs furtherinvestigations. C. Angular distribution
The projected energy distribution on the azimuth angleof the electron and the protons is calculated by integrat-ing the 6D scattering amplitudes as p N ( θ , E ) = (cid:90) d (cid:126)k (cid:90) dφ | (¯ (cid:126)k , (cid:126)k , T ) | ,(cid:126)k = [ φ , θ , (cid:112) M × E ] T (30)for protons, and p e ( θ , E ) = (cid:90) d (cid:126)k (cid:90) dφ | (¯ (cid:126)k , (cid:126)k , T ) | ,(cid:126)k = [ φ , θ , (cid:112) m × E ] T (31)for electron, where E and E are kinetic energies for anindividual proton and electron.As is observed in Fig. 4, the probability distributions ofelectron and protons reach the highest value in the polar-ization direction, which is consistent with the experimen-tal observations for linearly polarized laser pulses [11, 13].The probability of the dissociative protons is most con-siderable with 1 ≤ E ≤ eV , higher than the E < eV for dissociative channels reported in Ref. [11, 39], butin the range of their Coulomb explosion channel, wherethe laser wavelength is 800 nm. For higher intensity8 . × W/cm , the angular distribution of releasedprotons and electron extends more in the polarizationdirection. For distribution of protons, tiny yields around3 eV in radial coordinates indicates the Coulomb explo-sion channel, close to what is observed in experiments,however, for different laser pulses [39]. V. CONCLUSION AND DISCUSSIONS
We simulate the dissociative ionization of the H +2 ion infull dimensionality and obtain the ground energy same as FIG. 4. The log-scale probability distribution of (left col-umn) protons by log p N ( θ , E ) and (right column) protonsby log p e ( θ , E ), θ , ∈ [0 , π ]. The plot is symmetrized by p N (2 π − θ , E ) = p N (2 π − θ , E ) and p e (2 π − θ , E ) = p e (2 π − θ , E ) The upper row is computed from laser pulseat intensity 8 . × W/cm and lower row represents the5 . × W/cm . The values of the radial coordinates E / are represented in eV. The polarization direction is along thehorizontal axis and the direction electric field is labeled ateach sub-figure with an arrow above ” E ( t )”. The values areall normalized by dividing the maximum value. the quantum chemistry methods. Using tSurff methods,we obtained the JES where energy sharing is observed,which indicates a correlation between proton and elec-trons. The JES peaked at E N from 2 eV to 4 eV, whichis different from the previous 2D simulations, but is con-sistent with the experimental data. The projected energydistribution on angles shows that the electron and pro-tons tend to dissociate in the direction of polarization ofthe laser pulse.The simulation of the single emission spectrumshowing dissociation channels, is however not possi-ble yet. The difficulty lies mainly in constructing theinternuclear-distance-dependent electronic ansatz of H with a given ionic state in a single emission TDSE on (cid:126)r , which might be solved by reading the energy surfacesfrom quantum chemistry calculations or another tRecXcalculation. We leave work for future studies. ACKNOWLEDGMENTS
J.Z. was supported by the DFG Priority Programme1840, QUTIF. We are grateful for fruitful discussionswith Dr. Lun Yue from Louisiana State University, Dr.Xiaochun Gong from East China Normal University, andProf. Dr. Armin Scrinzi from Ludwig Maximilians Uni-versity.
Appendix A: Coordinate transformation
We use sub-indices ”a”, ”b” and ”e” represent the twoprotons and the electron of an arbitrary coordinate. Thesub-indices ”0”, ”1” and ”2” represent the center of thetwo protons, the relative position of an proton to the cen-ter and the electron in our transformed coordinate, re-spectively. Suppose originally the coordinates of the twoprotons and electrons are represented by vectors (cid:126)x a , (cid:126)x b and (cid:126)x e of an arbitrary origin, respectively. The new co-ordinates (cid:126)r and (cid:126)r satisfies (cid:126)r = (cid:126)x a + (cid:126)x b (cid:126)r = (cid:126)x a − (cid:126)x b (cid:126)r = (cid:126)x e − (cid:126)x a + (cid:126)x b , (A1)where (cid:126)r is the coordinate of the center of the two pro-tons. The Laplacian of the two protons (cid:79) a , (cid:79) b and elec-tron (cid:79) e are (cid:79) a = (cid:79) (cid:79) (cid:79) (cid:126) (cid:79) · (cid:126) (cid:79) − (cid:126) (cid:79) · (cid:126) (cid:79) − (cid:126) (cid:79) · (cid:126) (cid:79) (cid:79) b = (cid:79) (cid:79) (cid:79) − (cid:126) (cid:79) · (cid:126) (cid:79) (cid:126) (cid:79) · (cid:126) (cid:79) − (cid:126) (cid:79) · (cid:126) (cid:79) (cid:79) e = (cid:79) . (A2)Thus the kinetic energy of the system can be representedby −
12 ( (cid:79) a M + (cid:79) b M + (cid:79) e − (cid:79) M − (cid:79) M − (cid:79) M + (cid:126) (cid:79) · (cid:126) (cid:79) M + (cid:79) ≈ − (cid:79) M − (cid:79) m , (A3)where m = M M , and ” ≈ ” means the motion of the (cid:126)r is neglected. The interaction of the two protons with thelaser pulse can be written asi M (cid:126)A · ( (cid:126) (cid:79) a + (cid:126) (cid:79) b ) = i (cid:126)A · M ( (cid:126) (cid:79) − (cid:126) (cid:79) ) ≈ − i M (cid:126)A · (cid:126) (cid:79) , (A4)with which the total interaction with the laser field canbe written as − i (cid:126)A · ( (cid:126) (cid:79) + 1 M (cid:126) (cid:79) ) = − i β (cid:126)A · (cid:126) (cid:79) , (A5) FIG. 5. The (blue crosses) error δ (+) ( σ ) for R (+) c by Eq. B1with fixed R ( − ) c = 12 . δ ( − ) ( σ ) for computations for R ( − ) c with fixed R (+) c = 12 . . × W/cm isapplied. For details of the laser pulse, refer to Sec. IV A. Theparameters for angular momenta are L max = 7, M max = 5. where β = M +1 M . Appendix B: Convergence study
The errors are computed by the difference of JES fromtwo subsequent calculations σ ( E N , E e ) and σ (cid:48) ( E N , E e )with respect to R (+) c and R ( − ) c δ ( σ ) = max E N ,E e | σ ( E N , E e ) − σ (cid:48) ( E N , E e ) || σ ( E N , E e ) + σ (cid:48) ( E N , E e ) | , (B1)as used in previously in Fig. 2 (c). As depicted in Fig. 5,the JES is converged at R (+) c = R ( − ) c = 12 . Appendix C: Coulomb explosion with dispersion ofwave packets
We here illustrate the dispersion of the wavefunc-tion with r in spherical coordinate by a simple but un-physical model: the Coulomb explosion of a hydrogenatom. It starts with the ground eigenstate of a hydro-gen atom with electronic wavefunction ψ ( (cid:126)r ) = ψ ( r ) = √ π exp( − r ). The initial kinetic energy of the electron is0.5 atomic units and the potential energy is -1 atomicunits. Then the charge of the nucleus suddenly changesfrom +1 to -1, the kinetic energy remains constant butthe Coulomb potential reverses its sign. Thus the systemexplodes and the electron finally becomes a free particle.So what is the most possible final kinetic energy of thefree electron?In classical mechanics, the answer is 1.5 atomic unitsbecause of energy conservation. However, our quantumsimulation by tSurff gives an answer around 0.5 atomicunits, see Fig. 6 (a). We can also see that the distribu-tion of the spectrum has long tails to high energy region.The integration of σ ( E ) gives the total energy 1.5 atomicunits, which is consistent with the energy conservation,see Fig. 6 (a).The numerical simulations show the nuclear wave pack-ets also expand in space during the time propagation, seeFig. 6 (b), where the probability of the protons duringtime is computed by p wfN ( φ , θ , r , R , R , t ) = (cid:90) dφ (cid:90) sin θ dθ (cid:90) R R r dr | ψ ( (cid:126)r , (cid:126)r , t ) | . (C1)We split the radial coordinates into inner region r , r ∈ [ R , R ] = [0 ,
5] and outer region r , r ∈ [ R , R ] =[5 , [1] P. H. Bucksbaum, A. Zavriyev, H. G. Muller, and D. W.Schumacher, Phys. Rev. Lett. , 1883 (1990).[2] T. Zuo and A. D. Bandrauk, Phys. Rev. A , R2511(1995).[3] G. Yao and S. I. Chu, Phys. Rev. A , 485 (1993).[4] A. Giusti-Suzor, F. H. Mies, L. F. DiMauro, E. Charron,and B. Yang, J. Phys. B At. Mol. Opt. Phys. , 309(1995).[5] G. Jolicard and O. Atabek, Phys. Rev. A , 5845 (1992).[6] T. Zuo, S. Chelkowski, and A. D. Bandrauk, Phys. Rev.A , 3837 (1993).[7] B. D. Esry, A. M. Sayler, P. Q. Wang, K. D. Carnes, andI. Ben-Itzhak, Phys. Rev. Lett. , 013003 (2006). [8] J. H. Posthumus, Reports Prog. Phys. , 623 (2004).[9] A. Giusti-Suzor, X. He, O. Atabek, and F. H. Mies,Phys. Rev. Lett. , 515 (1990).[10] M. Odenweller, N. Takemoto, A. Vredenborg, K. Cole,K. Pahl, J. Titze, L. P. H. Schmidt, T. Jahnke, R. D¨orner,and A. Becker, Phys. Rev. Lett. , 143004 (2011).[11] M. Odenweller, J. Lower, K. Pahl, M. Sch¨utt, J. Wu,K. Cole, A. Vredenborg, L. P. Schmidt, N. Neumann,J. Titze, T. Jahnke, M. Meckel, M. Kunitski, T. Haver-meier, S. Voss, M. Sch¨offler, H. Sann, J. Voigtsberger,H. Schmidt-B¨ocking, and R. D¨orner, Phys. Rev. A - At.Mol. Opt. Phys. , 013424 (2014). FIG. 6. (a) The energy spectrum σ ( E ) = (cid:82) dθ (cid:82) dφk | b ( (cid:126)k, T ) | , E = k , where b ( (cid:126)k, T ) are the sin-gle electron scattering amplitudes. The spectrum is com-puted from advancing a hydrogen electronic ground statein the Coulomb potential + Zr . Z=0 means no externalCoulomb potential. The kinetic energy by (cid:82) ∞ σ ( E ) √ EdE = Z + 0 . p wfN ( φ , θ , r , , , t ) , r ∈ [0 ,
5] of the inner shell andlog p wfN ( φ , θ , r , , , t ) , r ∈ [5 ,
10] of the outer shell. Thetimes are selected where the | E z ( t ) | reaches a maximum, andmeasured in atomic units, except the 465 a.u. where the elec-tric field strength is weak. The two shells are split by thered, dashed circle, and are both normalized of each circle.The polarization direction is along the horizontal axis andthe direction electric field is labeled at each sub-figure withan arrow above ” E ( t )”. The absolute value of the outer shellis several orders smaller than the inner shell.[12] J. Wu, M. Kunitski, M. Pitzer, F. Trinter, L. P. H.Schmidt, T. Jahnke, M. Magrakvelidze, C. B. Madsen,L. B. Madsen, U. Thumm, and R. D¨orner, Phys. Rev.Lett. , 023002 (2013).[13] X. Gong, P. He, Q. Song, Q. Ji, K. Lin, W. Zhang, P. Lu,H. Pan, J. Ding, H. Zeng, F. He, and J. Wu, Optica ,643 (2016). FIG. 7. The normalized single JES is created from (a) 2D and(b) 6D σ S ( E N ) = (cid:82) σ (2 E N , E e ) √ E e dE e , where E N depictsthe kinetic energy of a proton. A laser pulse used in Ref. [19]is applied here. The 2D data is taken from Ref. [19].[14] C. B. Madsen, F. Anis, L. B. Madsen, and B. D. Esry,Phys. Rev. Lett. , 163003 (2012).[15] A. Scrinzi, New J. Phys. , 085008 (2012).[16] A. Zielinski, V. P. Majety, and A. Scrinzi, Phys. Rev. A , 023406 (2016).[17] L. Tao and A. Scrinzi, New J. Phys. , 013021 (2012).[18] J. Zhu and A. Scrinzi, Phys. Rev. A , 063407 (2020).[19] L. Yue and L. B. Madsen, Phys. Rev. A , 063420(2013).[20] L. Yue and L. B. Madsen, Phys. Rev. A , 063408(2014).[21] G. L. V. Steeg, K. Bartschat, and I. Bray, J. Phys. BAt. Mol. Opt. Phys. , 3325 (2003).[22] W. Qu, Z. Chen, Z. Xu, and C. H. Keitel, Phys. Rev. A , 013402 (2002).[23] R. E. Silva, F. Catoire, P. Rivi`ere, H. Bachau, andF. Mart´ın, Phys. Rev. Lett. , 113001 (2013).[24] N. Takemoto and A. Becker, Phys. Rev. Lett. ,203004 (2010).[25] B. Feuerstein and U. Thumm, Phys. Rev. A - At. Mol.Opt. Phys. , 043405 (2003).[26] K. C. Kulander, F. H. Mies, and K. J. Schafer, Phys.Rev. A , 2562 (1996).[27] V. P. Majety and A. Scrinzi, Phys. Rev. Lett. ,103002 (2015).[28] V. P. Majety, A. Zielinski, and A. Scrinzi, J. Phys. B , 025601 (2015).[29] V. Majety and A. Scrinzi, Photonics , 93 (2015).[30] V. P. Majety, A. Zielinski, and A. Scrinzi, New J. Phys. , 63002 (2015).[31] V. P. Majety and A. Scrinzi, J. Phys. B , 245603(2015).[32] J. R. Hiskes, Phys. Rev. , 1207 (1961).[33] A. Scrinzi, Phys. Rev. A , 053845 (2010).[34] C. W. McCurdy, M. Baertschy, and T. N. Rescigno, J.Phys. B , R137 (2004).[35] S. Chelkowski, T. Zuo, O. Atabek, and A. D. Bandrauk,Phys. Rev. A , 2977 (1995).[36] F. Chˆateauneuf, T. T. Nguyen-Dang, N. Ouellet, andO. Atabek, J. Chem. Phys. , 3974 (1998).[37] D. Bressanini, M. Mella, and G. Morosi, Chem. Phys.Lett. , 370 (1997).[38] L. J. Schaad and W. V. Hicks, J. Chem. Phys. , 851(1970).[39] D. Paviˇci´c, A. Kiess, T. W. H¨ansch, and H. Figger, Phys.Rev. Lett.94