t-SURFF: Fully Differential Two-Electron Photo-Emission Spectra

The time dependent surface flux (t-SURFF) method is extended to single and double ionization of two electron systems. Fully differential double emission spectra by strong pulses at extreme UV and infrared wave length are calculated using simulation volumes that only accommodate the effective range of the atomic binding potential and the quiver radius of free electrons in the external field. For a model system we find pronounced dependence of shake-up and non-sequential double ionization on phase and duration of the laser pulse. Extension to fully three-dimensional calculations is discussed.


Introduction
Differential double photo-electron spectra and the corresponding ionic recoil momentum spectra testify of dynamical correlation between the electrons. By sweeping extreme ultra violet (XUV) photon energies from below to above the threshold for single-photon double ionization of the He atom one probes correlation in initial and final states. In wave length in the infrared (IR) range, momentum distributions of recoil ions provide evidence for the importance of re-collision processes [1], where first one electron is ionized, which subsequently is re-directed by the oscillating laser field into a collision with its parent ion causing excitation and possibly detachment of the second electron. The early observation of unexpectedly enhanced double ionization of helium in IR fields [2] is now generally ascribed to this mechanism. Experimental data on strong field IR photo-ionization is also available for many other atomic and molecular systems and it was even proposed to use re-collision electron spectra for the analysis of structure and dynamics of molecules [3].
For the XUV wave-length, theoretical and experimental questions have matured, even if still under debate (see, e.g., [4] for a recent contribution to the debate with ample references to theory and experiment). At longer wave length, the large body of experimental data (see, e.g. [5,6,7,8,9] and references therein) and the somewhat smaller range of theoretical models largely based on classical or semi-classical methods (see, e.g., [10,11,12,13] ) are all plagued by the almost complete lack of reliable theoretical verification, with the notable exception of a few very large scale simulations of two-electron systems in strong fields [14,15,16], where, however, only in Ref. [15] the full two-electron dynamics is treated for laser wavelength of 800 nm. There are several reasons for this striking absence of complete ab initio simulations of ionization of twoelectron systems. Firstly, even single ionization is non-trivial to compute, if the external fields are non-perturbative. Roughly speaking, the effort for computing single ionization grows with the 4th power of the wave length λ 4 due to the growth of peak momenta ∝ λ, quiver amplitude of the electron motion in the field ∝ λ 2 and the growth of minimum pulse duration ∼ λ (see also the discussion in Ref. [17]). When the effect of the field is perturbative, the situation for single ionization relaxes somewhat, as one basically only needs to know the initial neutral state and the single-electron stationary scattering solutions in the energy range of interest. Although obtaining scattering solutions may be difficult, there is a well defined procedure and the whole technology of electronion scattering theory available to approach the problem. For double ionization, also in the perturbative regime the situation is more complex. The convenient partition into bound and singly ionized spectral eigenfunctions cannot be continued to above the double-ionization threshold: eigenfunctions above the double-ionization threshold will in general have, both, single-and double-ionization asymptotics. For distinguishing single from double ionization we therefore invariably need the solutions at large distances. Even without the need for asymptotic analysis, scattering with open double ionization channels is a challenging task.
Numerical simulations, in principle, can provide the full answer. The asymptotic analysis is usually done by propagating the wave function until after the end of the pulse and then analyzing its remote parts either in terms of momentum eigenfunctions, i.e. plane waves, or, when the tails of the ionic Coulomb potentials are not considered negligible, in terms of two-body Coulomb scattering solutions. Three-body scattering solutions, which would obviate a purely asymptotic analysis, unfortunately, are not accessible. The by far largest part of the computational effort in these simulations goes into following the solution to large distances until the pulse is over and analysis can begin. It may be irritating to think that a large effort is made to simulate dynamics that is known exactly in the case of finite range potentials or approximately at sufficient distances from all Coulomb centers: the free motion of one or two electrons in a laser field.
In a preceding publication [17] we have shown how, by absorbing the wavefunction and recording of flux before absorption, single particle photo-electron spectra can be computed using simulation volumes that only contain the relevant range of the atomic potential and the electronic quiver motion in the field. In this so-called time-dependent surface flux (t-SURFF) method, asymptotic information is accumulated during time propagation rather then drawn from the full wave function after the end of the pulse. The equivalence of both approaches was proven in [17] mathematically and by comparing numerical results with literature. Depending on the system parameters and accuracy requirements, absorption radii for typical strong field setups can be kept as small as 20 Bohr radii and as few as 90 linear discretization coefficients for the radial motion can give better that 1% accurate spectra over the complete spectral range. For Coulomb systems, due to the long interaction, larger radii of ∼ 100 Bohr and 200 to 300 discretization points may be needed.
In the present paper, we extend this approach to two-electron systems and multi-channel single ionization as well as double ionization spectra. A numerical demonstration is provided using a two times one-dimensional model system. As first physical results, we find phase and pulse-duration dependence of shake-up and doubleionization spectra, phenomena that are consistent with the spatial asymmetry of very short laser pulses and the re-collision mechanism for shake-up and non-sequential double ionization.

The t-SURFF method for two-electron systems
We first briefly summarize t-SURFF for the single particle case. The basic requirement is that there is some radius R c beyond which the Hamiltonian reduces to an asymptotic one with known solutions. Let H(t) denote the time dependent Hamiltonian of our system and assume that there exists an exactly solvable H v (t) that at large distances agrees with H(t): (1) For a single particle in a laser field described in velocity gauge and a short range potential V ( r) ≡ 0 for | r| > R c , H v (t) is the Hamiltonian for the free motion in the laser field where Here and throughout we use atomic units = m e = e 2 = 1 unless indicated otherwise (m e and e denote electron mass and unit charge, respectively, the Bohr radius results as the atomic unit of length). The TDSE with H v (t) has the Volkov solutions i.e. plane waves times with time-dependence by the well-known Volkov phase Scattering describes the behavior of the time-dependent wave function Ψ( r, T ) at long times T and large distances | r| > R c . We choose T and R c large enough such that at T the pulse is over and the wave function has split into its bound and scattering parts with the properties The approximate sign in (7) refers to the tails of any bound state function, which decay exponentially with increasing R c . The approximate sign in (8) refers to the fact that electrons with very low energies k 2 /2 ∼ 0 may not have passed R c at time T . We only need to analyze Ψ s ( r, T ) in terms of asymptotic functions χ k ( r, T ). As Ψ s ( r, T ) vanishes inside the radius R c , we can multiply the full Ψ( r, T ) by the function write the emission amplitude for photo-electron momentum k as For obtaining a time integral over a surface flux, we write (10) as an integral of the time derivative and use the fact on the support of θ(R c ) both, χ k ( r, t) and Ψ( r, t) evolve by the same Hamiltonian H v (t). We obtain The commutator vanishes everywhere except on | r| = R c : the asymptotic information is obtained by integrating the time dependent flux through a surface at finite distance R c . For the further discussion of the single-particle case and numerical examples with short-range and Coulomb potentials we refer to Ref. [17]. Figure 1. Partitioning of coordinate space into bound (B), singly ionized (S, S) and doubly ionized D regions. Single and double photo electron spectra are obtained by integrating the flux across the boundaries between the regions.

Single-ionization into ionic ground and excited state channels
The simplest extension of the single electron method is for computing single ionization into ground and excited state ionic channels. The complication is that in presence of a strong field the ionic state can differ from the field free ionic state due to polarization and one must make sure to count flux passing the surface into the correct ionic channel. In this section we consider only the lowest ionic states that remain bound and do not contribute to double ionization.
For decomposing coordinate space into bound and asymptotic regions we define on both coordinates r 1 and r 2 projector functions θ 1 ( r 1 , R c ) and θ 2 ( r 2 , R c ), respectively, as in the single particle case (9). Again picking sufficiently large times T and a sufficiently large surface radius R c we can partition the wave function Ψ(T ) into its bound, singly ionized, and doubly ionized parts (see Figure 1) Here and in the following we suppress the arguments r 1 , r 2 and R c of θ 1 and θ 2 . Note that the assignment of singly and doubly ionized character to the different regions is asymptotically exact: the error can be made arbitrarily small for any specific solution Ψ(T ) by choosing sufficiently large T and R c . The single ionization Hamiltonian agrees with the exact Hamiltonian on the support of θ 1 (1−θ 2 ) (region S in Figure 1). The corresponding Hamiltonian on S is obtained by particle exchange. Channel solutions have the form Here φ c ( r 2 , t) solves the ionic TDSE Rather than an initial we use a final condition that φ c ( r 2 , t) for t > T should be the desired final ion state (we need to solve the TDSE backward in time). The channel solution on the support S of the particle exchanged projector the probability density for finding at time T an electron with momentum k in ionic channel c is The factor 2 arises from adding the two identical exchange symmetric contributions. For converting this integral into a time integral over a surface we make the simplifying assumption that the ionic solution never leaves the bound area This is the precise version of the assumption that the ionic states considered does not get further ionized. The approximate sign refers to the fact that again there is always an exponential tail reaching to arbitrary distances and further that interaction with any pulse will lead to a small amount of ionization. Using (20) we neglect the flux between the singly ionized regions S, S and the doubly ionized region D at all times. Then, using the same procedure as for a single particle we obtain where ψ c ( r 1 , t) is the channel projected wave function For evaluating the integral (21) we only need to know values and radial derivatives of ψ c on the surface | r| = R c . For the computation it means solving the full two-electron problem up to time T and up to radius R c . Beyond R c one can absorb all amplitudes. In addition, for each channel c, we need to solve one single electron problem up to the same time and radius (which is usually a much simpler task).

Double ionization spectra
When double ionization occurs, flux passes from the bound region B through the singly ionized regions S or S into the doubly ionized region D. Our naming of the areas is suggestive but does not imply any bias as to an actual state that the electrons occupy within any of these areas. It is unsubstantial for the present discussion whether some intermediate ionic bound state is occupied by one of the electrons in S or S (sequential ionization) or whether both electrons must be considered unbound. The sole purpose of the partitioning is to have well-defined surfaces outside the ranges of the respective potentials where we will integrate fluxes.
For double ionization there is one obvious limitation of the discussion so far: on the line | r 1 − r 2 | = a the electron-electron interaction is constant and not negligible for small a. This problem is not related to the long-range Coulomb potential, rather it is a general problem for any multi-particle breakup, which is why break-up processes are more complex than single particle scattering. Within the framework of the present approach the problem can be solved for short-range electron-electron repulsion without making approximations, which will be discussed below. A pragmatic solution has been effectively employed in many earlier publications, which is to neglect electron-electron repulsion at large distances from the nucleus. This is what using any projection onto products of single electron states implies, be it Coulomb scattering waves or plane waves (both approaches are discussed, for example, in [18]). Sensitivity to this approximation can be tested by varying the distance from the nucleus where the projection starts.
As the pragmatic solution was found to work well in many cases, we make this approximation explicit in the present paper by smoothly turning off all potentials including the electron-electron interaction before the surface radius R c . In that case we can always use the free (Volkov) Hamiltonian H v (t) beyond R c .
By our assumptions, in the region | r 1 | ≥ R c , the Hamiltonian is identical to H s (t), Eq. (14), which motivates the ansatz with the Volkov solutions χ k ( r 1 , t) on coordinate r 1 and an expansion into a timeindependent, complete, but otherwise arbitrary set of functions ξ n ( r 2 ) on r 2 . Using orthogonal projection onto the expansion functions χ k 1 and ξ n , the coefficients b ′ ( k 1 , n, t) are obtained as The integral (24) over q θ accounts for the fact that the plane waves are not δ-orthonormal when the integration is restricted by θ 1 : the inverse overlap is defined by For notational brevity, we assume here that the basis functions ξ n are orthonormal. The time-derivative of the b( k 1 , n, t) is By the double right bracket we emphasize that integration is over both coordinates r 1 and r 2 . Inserting the representation (23) into the first term we obtain an inhomogeneous equation for the b( k 1 , n, t): where we have used (26). The inhomogeneity is the flux through the surface | r 1 | = R c . Initial conditions are b( k 1 , n) ≡ 0, i.e. no electrons outside R c . We write the double ionization amplitude and use one more time the conversion to integrals over surface flux The two terms B and B are related by exchange symmetry For computing B, we only need to know θ 1 Ψ(t), for which we insert the representation (23) to obtain The inverse overlap q θ ( k 1 , k ′ 1 ), (26), cancels with the overlap integrals and never needs to be evaluated explicitly.
B( k 1 , k 2 , t) is the contribution to the double ionization spectrum passing at time t through surface | r 2 | = R c from region S into D. In S, the first electron is already detached and has a fixed canonical momentum k 1 that is carried into the double-ionized region D. B( k 1 , k 2 , t) is the alternate contribution going through region S.

Computational remarks
The substantial gain of the method is that, rather than computing the full solution in region D and then analyzing it, we only need to integrate the flux through the surface separating S from D. In the direction parallel to that surface the wave function is represented in terms of the free solutions χ k 1 ( r 1 , t), where we do not need to expand the wave function completely, but can restrict propagation to the momenta k 1 that we are interested in. For each k 1 we need to solve equation (28), which is an ionic TDSE with an additional source term (the flux entering S from the bound region B). Although the derivation may appear complex, the implementation of the procedure can be done by the following simple algorithm • Set up b( k 1 , n, t) on a k 1 grid of the desired density and initialize to 0.
• Get time grid points t i with a time step that resolves the fastest amplitude oscillations for the desired energy range ∆t < 2π/E max • In a loop through all times t i , get the surface terms from file or from a simultaneous calculation of |Ψ(t i ) and use (28) to advance to b( k 1 , n, t i ).
• Add the contribution from region S to the spectral amplitude The contribution s from region S is obtained from s by particle exchange. The total spectral amplitude is the sum of both contributions The asymptotic value -the spectral amplitude -is attained at times T when the flux through the surfaces becomes negligible. The approach can only be successful, when absorption does not significantly distort the solution at the integration surfaces. The "infinite range exterior complex scaling" (irECS) absorber was shown in Refs. [17,19] to provide traceless absorption over a very wide energy range at low computational cost. It outperforms standard complex absorbing potentials by several orders of magnitude in accuracy. If needed, irECS can be pushed to full machine precision using not more than 20 discretization coefficients per coordinate in the absorbing region | r| > A 0 . Absorption can begin at any A 0 ≥ R c . More mathematical and numerical detail and ample numerical examples using irECS can be found in Refs. [17,19].
In general discretization errors for two electrons are similar to those for a singleelectron system with the same ionization potential. This has the practical advantage that a suitable discretization can be determined from the single-electron problem. Only a few consistency checks on specific two-electron observables need to be performed for the computationally heavier two-electron calculation.
2.3.1. General single ionization spectra Once we have obtained the b( k 1 , n, T ) we can reconstruct the wave function for | r 1 | > R c . In particular, the amplitude for single ionization spectra for any ionic state φ c is where a (c) n are the expansion coefficients of φ c with respect to the basis functions ξ n : φ c ( r) = n ξ n ( r)a (c) n .
This approach is not limited to non-ionizing φ c , as it analyzes ionic population at time T after the end of the pulse. For non-ionizing φ c , it is an alternative to the single-ionization procedure above. Note that here we need to solve the inhomogeneous ionic problem (28) for each photo-electron momentum k. The total spectrum is a linear combination of the individual contributions from each n. Where it is applicable, the advantage of the singleionization procedure of section 2.1 twofold: firstly, for each final φ c one needs to solve only one ionic TDSE and compute values and derivatives of the channel surface function ψ c . The complete spectrum can be obtained by time-integration with little numerical effort. Secondly, as one directly obtains the channel spectrum without intermediated decomposition and final resummation of the wave function, results are in general more robust numerically.

The one-dimensional two-electron model system
As even with the simplifications by t-SURFF the solution of the full three-dimensional two-electron problem remains a very large scale computational task, we use for demonstration purposes the standard one-dimensional two-electron model Hamiltonian For making all potentials strictly finite range, we have chosen a "truncation radius" C p and a truncation function M(x) that is ≡ 1 up to |x α | = C p − 5 and goes differentiably smoothly to 0 at |x α | = C p . Where not indicated otherwise we use C p = 20. The screening factor of 1/2 in the ionic potential was chosen for esthetic reasons, as then the exact ionic ground state (without potential truncation) is E 0 = −2. The first excited state occurs at E 1 = −0.93, substantially stronger bound than the first excited He + energy of −1/2. With the electron-electron screening of 0.3 one obtains an ionization potential of 0.88, which we consider a fair approximation of the actual He ionization potential of 0.90. For all computations we use a finite elements discretization where orders between 8 and 20 where used for convergence studies. For the present purposes we found order 8 or 10 sufficient. A deeper discussion of the finite element method for irECS calculations can be found in [19]. The irECS radius A 0 which marks the beginning of absorption was varied between 20 and 100 atomic units. Results do not depend on A 0 on the level of accuracy shown.
In all calculations below we use cos 2 pulse shapes defined in terms of the vector potential where T F W HM is the full-width-half-maximum of the envelope of the vector potential, ω is the laser central photon energy, and φ CEO the carrier-envelope offset phase. The peak field amplitude E 0 is related to the pulse peak intensity I by E 0 = √ 2I. By construction, the electric field E(t) = − d dt A(t) has no zero frequency component, which is important when studying effects of φ CEO with very short pulses. For a "cosine pulse" φ CEO = 0 the peak of the electric field approximately coincides with the maximum of the pulse envelope with minor deviations due to the derivative of the envelope.

Two photon double ionization in the extreme ultraviolet
In Ref. [4] perturbative two-photon double ionization of He by short XUV pulses with photon energies between the two-and single-photon ionization thresholds was investigated and a remarkable universal description of the process was found. The findings were confirmed by numerical solutions of the fully three-dimensional twoelectron TDSE. Pulse durations of T = 4.5 f s and photon energies ω between 42 and 80 eV were used. The perturbative regime at ω = 42 eV extends up to intensities near I = 10 16 W/cm 2 , where one first observes sizable deviations from the perturbative scaling ∝ I at the single-photon single-ionization peak and ∝ I 2 for shake-up and double ionization.
As an example, Figure 2 shows the two-electron spectrum for photon energy ω = 70 eV at peak intensity 10 16 W/cm 2 . Three two-electron at energies |E 1 |+|E 2 | = nω+E 0 are clearly visible. With ground state energy E 0 = −2.88 au ≈ −78 eV the first peak appears for n=2-photon ionization. The ridges parallel to the energy axes are artefacts of the method: there is significant shake-up and the exponential tail of the excited ionic states reaches into region D. Note, however, that with the choice of R c = 40 for the present calculation, the ridges are suppressed by at least a factor 10 −6 relative to the surrounding signal. Even stronger suppression can be obtained by further increasing R c or, alternatively and with slower convergence rate, propagating to longer times. If the ionic states are known accurately, their effect can be completely removed by projection (for a detailed discussion of the procedure see Ref. [17]).
In Figure 3 we show line-outs of the two-photon two-electron energy ridges for photon energies between 42 and 80 eV. The curves are converged within the resolution of the plots with respect to T , R c and spatial discretization. Although the general characteristics are determined by the single-excitation resonances as in [4], for the present system we do not reproduce the universal behavior found there for the angleintegrated spectra. The likely reason are shake-up interferences that also modify the the angular distributions in the 3D case [4].
Even at the present benign, nearly perturbative parameters, photo-electrons with   energies of about E ϕ ∼ 3au ≈ 80 eV travel to a distance r max ∼ 2T F W HM 2E ϕ /m e ≈ 900 au during the duration of the pulse. To correctly represent the corresponding momenta one needs a grid spacing of ∆r 2π/ √ 2E max ≈ 2.5, which results in at least about 400 radial discretization points in a numerical calculation. In practice, finer grid spacings need to be used with a correspondingly larger number of discretization points. In addition, for asymptotic analysis of the wave function, one may need to propagate further for a certain period of time after the pulse, leading to further increase in the required box radius r max proportional to propagation time.
For comparison, with t-SURFF we obtain converged results using only 49 discretization coefficients on the positive half axis [0, ∞). The total number of discretization points per coordinate is twice as large because of the negative halfaxis, which would correspond to a different angular direction in a three-dimensional calculation. 32 out of the 49 points were used on the interval [0, 20], i.e. an effective grid spacing of ∼ 20/32 = 0.625, which sets a theoretical limit for the maximum photoelectron energy of 50 au ≈ 1300 eV . The remaining 17 points were used for irECS absorption. The t-SURFF simulation volume radius R c is independent of pulse duration.

Single ionization and shake-up by an IR pulse
Possibly the most elementary correlated process is shake-up: after forcefully removing a single electron, the remaining electrons rearrange not exclusively into the ionic ground state, but a fraction goes into excited states. The exact distribution of excited states strongly depends on system parameters. Shake-up by very strong, short infrared pulses has been observed recently using the "attosecond transient absorption spectroscopy" [20], which can monitor the time-evolution of an excited ionic state after IR ionization.
Within the present simplified model we can only demonstrate basic qualitative features of shake-up in IR ionization. We have studied dependence on the carrierenvelope offset phase φ CEO and on pulse duration. We use an intensity of 2×10 14 W/cm 2 and ω = 0.057au (corresponding to wave length 800 nm). Converged results were obtained using the same discretization parameters as in the XUV case discussed above.
In Figure 4 photo electron spectra for the ground and first excited states calculated by formula (21, smooth line) and, alternatively, by Eq. (37, coarser energy grid) for a single cycle cosine pulse: T = 2π/ω and φ CEO = 0. The spectra show some generic short pulse features: pronounced asymmetry and absence of individual photo-electron peaks at energies, where emission occurs only during a single laser half-cycle. Emission to the left is lower, as the field amplitude to the left is smaller.
Clearly, the spectral structure is highly sensitive to φ CEO . However, also the total shake-up yield is strongly φ CEO dependent. Table 1 list the yields in the ground and first excited ionic states for different phases and pulse durations. While shake-up is strongly suppressed for a single cycle cosine pulse, it increases to a sizeable ∼ 10% of the ground state ionic channel for φ CEO = 3π/4. We tentatively ascribe this fact to a recollision mechanism for shake-up. The fact that longer pulses with 2 and 3 optical    Table 1. Shake-up as a function of φ CEO . Total yield in the ionic ground Y 0 and first excited Y 1 states and ratio Y 0 /Y 1 .
cycles in general have shake-up fractions comparable to the φ CEO = 3π/4 is consistent with this hypothesis. We abstain from a more detailed analysis of the phenomenon for the present model, as due to the absence of transverse wave-packet spreading, recollision effects are greatly exaggerated in one dimensional models.
3.4. Double ionization spectra generated by an IR pulse Figure 5 gives an overview of doubly-differential photo emission spectra obtained for different φ CEO and pulse durations up to three optical cycles T F W HM = 6π/ω au ≈ 7.5 f s. The calculations were performed using a slightly more accurate discretization with 60 points on the half-axis [0, ∞).
As to be expected, we see pronounced asymmetries and strong effects of φ CEO for the shortest pulses. While the single-cycle cosine pulse essentially only shows weak and unidirectional emission of both electrons, backward emission and significantly higher yields arise at φ CEO = 3π/4. Most likely, the responsible mechanism is re-collision. However, due to the limitations of the one-dimensional model, we defer a more profound discussion of the effects to a future fully three-dimensional calculation.
3.4.1. Dependence on the potential cutoff By truncating all potentials at C p = 20 au, the Volkov Hamiltonian H v (t) is exact outside R c and all errors of t-SURFF are either due to discretization or absorption. As a pragmatic approach such a truncation is very appealing, but the truncation error must be controlled. For our simple model, we made a series of calculations for truncation radii up to C p = 100. Figure 6 shows the single photo-emission spectrum for the first excited ionic channel, a diagonal lineout E 1 = E 2 of the double-emission spectrum, and two lineouts where one energy energy is fixed at E 1 = 0.1 au(2.7 eV ) and E 1 = 0.5 au(13.6 eV ). The single electron spectra are hardly affected by the truncation. The strongest effect, as to be expected, is along the line E 1 = E 2 . With C p = 20 qualitative differences in the two-electron plots can be observed, while good qualitative agreement can be found for C p = 50. The values C p where results are acceptable must be chosen depending on the system and on accuracy requirements. In our case, the on-diagonal spectra are still incorrect up to C p = 50.  In contrast, convergence for larger momentum differences appears quite acceptable at C p ≈ 50.
Alternatively to increasing C p , one may use a different arrangement of surfaces, where an additional surface is placed at | r 1 + r 2 | = R c (see Figure 7). In scaled centerof-mass and interparticle coordinates r ± := ( r 1 + r 2 )/ √ 2 the asymptotic Hamiltonian H D (t) in region D separates into the center-of-mass motion with a Volkov solution and field-free relative motion in the repulsive potential of the electrons where ∇ + and ∆ − denote the Nabla and Laplace operators in coordinates r + and r − , respectively. Thus the asymptotic problem reduces to field-free potential scattering and free motion in the laser field. The contribution going from region B directly into D can be conveniently expanded in this form. We leave the technical discussion and numerical demonstration of this procedure for future work.

Discussion and conclusions
With the reduction of the simulation volume afforded by t-SURFF, fully dimensional calculations of double photo-ionization for a broad range of wave-length and intensities has come into reach. Scaling to three dimensions depends on the laser frequency: in the perturbative regime, the number of angular momenta required may remain as low as 3 [18], which, in linear polarization, scales the problem size by ∼ 3 3 = 27. In that case the computational effort is less than one order of magnitude larger than the calculations presented in section 3.2, if one considers that the positive and negative half-axes of the one-dimensional model lead to a 4-fold scaling of the problem compared to a radial problem. At infrared wave length, the complete angular phase-space becomes activated. At 800 nm wave length and intensities of ∼ 2 × 10 14 W/cm 2 ≈ 30 angular momenta are needed for full convergence of electron spectra. For two-electron spectra, this entails an increase in problem size by ∼ 30 3 /4 ≈ 7000 compared to the present calculation, if we keep the same standards of accuracy. Considering that a typical solution for our problems takes 2 hours on a single CPU, one sees that 3D calculation are quiet feasible on a moderate size parallel computer. The fact that in the 3D case the fraction of the phase space, where electron repulsion is important, is much smaller than in 1-d, may possibly be exploited for further reduction of the problem size. The development of t-SURFF is motivated by the photo-ionization problem. What is special about that problem is that in dipole approximation the time-dependent field affects the whole wave function, including the asymptotic region and therefore final momenta cannot be determined before the end of the pulse, unless one uses knowledge about the time-dependence of the asymptotic solution. However, also in situations without time-dependence of the asymptotic Hamiltonian, the method may turn out to be useful for obtaining fully differential momentum spectra. All it requires is a reliable absorption method and knowledge of the solution in the asymptotic region. Among the candidates for further application are reactive scattering and chemical break-up processes.