# A combined R-matrix eigenstate basis set and finite-differences propagation method for the time-dependent Schröddinger equation: the one-electron case

## Abstract

In this work we present the theoretical framework for the solution of the time-dependent Schrödinger equation (TDSE) of atomic and molecular systems under strong electromagnetic fields with the configuration space of the electron’s coordinates separated over two regions, that is regions and . In region the solution of the TDSE is obtained by an R-matrix basis set representation of the time-dependent wavefunction. In region a grid representation of the wavefunction is considered and propagation in space and time is obtained through the finite-differences method. It appears this is the first time a combination of basis set and grid methods has been put forward for tackling multi-region time-dependent problems. In both regions, a high-order explicit scheme is employed for the time propagation. While, in a purely hydrogenic system no approximation is involved due to this separation, in multi-electron systems the validity and the usefulness of the present method relies on the basic assumption of R-matrix theory, namely that beyond a certain distance (encompassing region ) a single ejected electron is distinguishable from the other electrons of the multi-electron system and evolves there (region II) effectively as a one-electron system. The method is developed in detail for single active electron systems and applied to the exemplar case of the hydrogen atom in an intense laser field.

## I Introduction

Exploration of the fundamental processes that occur when atomic and molecular systems are subject to extreme conditions is currently a major research area. Experimentally, such processes are realized by strong and/or short intense laser pulses radiating at infrared wavelengths Posthumus (2004); Zeidler et al. (2006) and have recently been utilised at a more practical level for reconstruction of nuclear probability distributions, visualisation of molecular orbitals, alignment of molecules as well as production of high-order harmonics which in turn are used for the generation of ultra short fields at the attosecond scale Niikura et al. (2002, 2003); Alnaser et al. (2004); Itatani et al. (2004); Alnaser et al. (2003); Rottke et al. (2002); Zeidler et al. (2006).

Theoretically, it is a huge task to treat the exact time-dependent (TD) response of a multi-electron system subject to a strong electromagnetic (EM) field by ab initio methods. In response to extensive experimental achievements using high-intensity Ti:Sapphire laser sources in the long wavelength regime, many theoretical studies employed the strong-field approximation where the influence of the Coulomb potential on the ejected electron wave function is neglected in favour of the external field. A more sophisticated approach that adopts the single-active-electron (SAE) approximation was also applied to the atomic case Kulander (1987). SAE models where one reduces the dimensionality of the multi-electron problem by freezing the most tightly bound electrons have proven to be very useful in cases where multiple electronic excitations are insignificant, and the SAE approximation is probably the most widely used approach when studying phenomena such as single ionization, above-threshold ionization (ATI) and high-harmonic generation (HHG).

For systems of only two electrons, such as the negative hydrogen ion, helium, molecular hydrogen, direct, ab-initio, solutions of the time-dependent Schrödinger equation (TDSE) have appeared in the early nineties (for a review see ref. Lambropoulos et al. (1998)). Since then, the computational power has increased steadily and as a result these methods have reached a high level of accuracy, efficiency and reliability, tackling successfully the very demanding theoretical problem, of single and double ionization of helium at 390 and/or 780 nm Parker et al. (2000, 2006).

Recently, the construction of FEL sources which deliver brilliant radiation in the soft- and (in the immediate future) hard X-ray regime have initiated new challenges in the field of atomic and molecular physics Wabnitz et al. (2005); Sorokin et al. (2006). However, in contrast to what occurs with conventional laser sources, more than a single electron at a time responds to short wavelength FEL light and X-ray FEL light will interact preferentially with the inner-most electrons, residing closer to the system’s core, rather than with the valence ones. An immediate consequence of the above property is that theories such as the SAE and models not taking into account interelectronic interactions at a sufficient level are inadequate to describe the processes involved. Moreover, high-order harmonic generation (HOHG) techniques are nowadays able to create pulses of subfemtosecond duration. Given that relaxation processes, such as Auger transitions, of the bound electrons are of the order of a femtosecond or less it can be concluded that the short time-variation of the EM field requires approaches where multi-electron dynamics can be reliably described.

Given our intention to study multi-electron systems under intense EM ultrashort fields, there is considerable importance in the development of computationally tractable methods able to treat multi-electron systems with the least approximations possible. Such approaches have been developed in atomic and molecular physics studies, and include variants of time-dependent Hartree-Fock (TDHF) Kulander (1987). Though a vast number of theoretical efforts in the spirit of TDHF Kulander (1988); Pindzola et al. (1991); Krause et al. (1991); Pindzola et al. (1993, 1995) have appeared, even some extensions to include correlation between the electrons, the question of how much and under what conditions correlation beyond the Hartree-Fock model is important still remains unanswered. The underlying reason is the difficulties introduced by the nonlinear nature of the TDHF equations in combination with the fact that the single-configuration ansatz and the excitation process induced by the EM field are inconsistent. Improvements of the restricted Hartree-Fock ansatz and inclusion of exchange effects appear to be possible solutions to overcome such problems, although the applications so far are only in one-dimensional (1D) models Horbatsch et al. (1992); Dahlen and van Leeuwen (2001); Nguyen and Bandrauk (2006); Hu et al. (1998); Caillat et al. (2005).

An alternative ab-initio approach capable of treating multi-electron systems is R-matrix theory, with the basic formulation appeared first in the context of nuclear theory, and later on applied in the field of atomic physics (Burke and Seaton (1971); Burke and Robb (1975); Burke and Taylor (1975)). Traditionally, R-matrix theory is a theory where time is not involved in the study of the collision or photoionization processes. Variants of R-matrix theories and computational codes have been applied to an impressive number of systems, over the last 40 years Burke and Berrington (1993). With the advent of strong and/or short laser pulse technology an early application of R-matrix theory to multiphoton processes appeared in the form of a Floquet expansion of the driven time-dependent wavefunction Burke et al. (1991). Although able to treat the field non-perturbatively, the R-matrix Floquet approach cannot be considered as a fully TDSE solution methodology since it is only suited to laser pulses containing many cycles.

Similarly, the appearance of high power sources at the short wavelength regime has led a number of theoretical groups to develop TDSE approaches based on R-matrix theory (van der Hart et al. (2007a); Guan et al. (2007); van der Hart et al. (2007b)), with the first work to this end appearing some years ago Burke and Burke (1997). The basic assumption of R-matrix theory is very well suited to the physical situation of the photoionization process involved in light-matter interaction. Under strong radiation any system will ionize either multiply or singly. In the regime of single ionization the ejected electron, after some time, depending on its distance from the core, can be safely identified as distinguishable from the other electrons. In R-matrix theory this is taken into account through the division of configuration space into two regions where, in the inner region (region ), all interelectronic and exchange effects between all the electrons are treated, while in the outer region (region ) the ejected electron evolves effectively as a one-electron system under the influence of the residual core and the potential due to the remaining electrons. Thus in the outer region no matter what particular process has taken place the system wavefunction consists entirely of that of the wavefunction of the ejected electron.

The purpose of this work is two-fold. The first is to pursue development of a method which meets the above requirements for more complex systems than one- and two- electron systems and where atomic structure plays an important role in the processes. For this, a method based on R-matrix basis eigenstates appears to be tractable due to its success in describing such complex systems. Second, and equally important is the issue of efficiency and accuracy. It is inevitable that the demands of the calculations will make the study of such problems computationally very demanding. Finite-differences with high-order explicit time propagators Smyth et al. (1998) although difficult to use throughout configuration space in a direct extension to multi-electron systems, have proven to be very efficient and accurate in solving the TDSE for one- and two-electron systems. In fact the HELIUM code Smyth et al. (1998) using such methods to solve the TDSE fully for a two-electron atom exposed to intense laser fields is able to run with high efficiency in both computation and communication over many thousands of cores on the largest supercomputers presently available. This established efficiency makes their implementation for the outer region in our present approach a very reasonable one. In region , an R-matrix basis set is used to propagate the multi-electron wavefunction while in region amounting effectively to a one-electron problem, a finite-difference high-order propagation algorithm is used. Since to the best of our knowledge no such attempt has appeared, namely the propagation of the TDSE in a combined basis and grid representation of the TD wavefunction, we consider it essential to set out carefully in detail the basics of the method, free from complications arising from multi-electron considerations. Thus, we provide below the details of the method and its usefulness for one-electron systems and present results for the hydrogen system where accurate ab-initio methods, to compare with, are available to us.

The paper is organized as follows. In Sec. II we give an overview of the basic ideas and principles. Section III is the key section of this paper and there we set out in detail the theoretical formulation for a one electron system. In Sec. IV we apply the method to the hydrogen atom in an intense laser field which serves as an exemplar. We have relegated to appendices some of the more technical details. Finally we set out some conclusions and perspectives with regard to the new method in Sec. V. Atomic units are used () throughout.

## Ii Overview of the basic ideas and principles

As mentioned briefly in the introduction, the basic assumption of R-matrix theory for the outer-region wavefunction allows the derivation of a TDSE (in the outer region), where only one electron is involved reducing the dimensionality of the problem there to its minimum, namely to at most three, thus simplifying the computational problem considerably. To put this in a more quantitative fashion, let us recall the electron wavefunction beyond a certain distance, say (taken as the inner boundary of the outer region ) Burke and Berrington (1993):

(1) |

with , and . The are channel functions formed by coupling the target states of the residual atomic system , described by the Hamiltonian and the angular and spin quantum numbers of the ejected electron. The radial motion of the ejected electron (in the -channel) is described by the radial channel functions . The absence of the antisymmetrization operator is essential in the above expansion since it relies on the ejected electron and the remaining electrons occupying different portions of configuration space, thus making the ejected electron distinguishable from the others. Let us now consider the TDSE of the above system, in an external time-dependent radiation field. By writing the Hamiltonian for the field-free electron system as we end up with the following form for the TDSE:

(2) |

with denoting the interaction operator between the system and the external field, in the dipole approximation. Projection of the known channel states onto the TDSE and integration over and results in the following set of coupled partial differential equation for the radial motion in channels ,

(3) |

By properly ordering the radial channel functions into a column vector and the evolution operators and into a square matrix we may (in the outer region ) rewrite the TDSE of the ejected electron of any multi-electron system in the case of single ionization as,

(4) |

this equation having essentially the form of the one-electron TDSE. It is exactly this last equation, no matter how the inner region is treated, that allows us to utilize any propagation technique in the outer region of configuration space, which may have already been applied to one-electron ionization.

On the other hand, in the inner region an eigenstate representation of the TD wavefunction will result in a TDSE where only two dynamical quantities are needed to be provided for the forward propagation in time of the solution, namely eigenenergies and transition matrix elements between the system’s eigenstates. The key point in this case is that the whole information about the exact nature of the system described in the inner region, whether multi-electron or not, is contained in the values of the energies and the transition matrix elements together with the required selection rules for the transitions. Therefore, in a sense, without trying to oversimplify, one would expect the matching procedure between the two methods (inner region/basis representation - outer region/grid representation of the wavefunction) to hold regardless of the actual system being multi-electron or single-electron in nature. It is for this reason we believe the formulation in the present work should be readily extendable to complex multi-electron systems. The theoretical details and subsequent application will be more complicated, due to the multiplicity of ionizing channels for the ejected electron in such cases. In the following sections we will develop our approach for the one-electron atom case in detail thereby laying bare the basic concepts of what we believe a novel combination of basis set and finite-difference methods.

## Iii The theoretical framework

In this section we develop a theory for solving the TDSE using basis and grid representations in the inner and outer regions, respectively. The artificial R-matrix division of configuration space into two regions causes time-dependent boundary terms to appear in the corresponding TDSEs (in regions and ) which exactly account for the amount of probability current passing through the boundary during the interaction with the external field as well as after its turn off. Since the time-dependent wavefunction consists of two parts, a careful analysis is necessary in order to obtain the physical observables of interest such as bound and ionization probabilities as well as energy and angular information on the ejected electron.

In Sub-Sec. III.1 we present the calculation of the R-matrix eigenstates defined in region and derive the time-evolution equations for a wavefunction expanded over the R-matrix eigenstates of the field-free Hamiltonian. In Sub-Sec. III.2 we derive the finite-difference TDSE governing the radial motion of the ejected electron. In Sub-Sec. III.3 we summarize the calculational procedure for the forward in time propagation of the wavefunction. Finally, in Sub-Sec. III.4 we give the formal expressions for the calculation of experimental observables adapted to our methodology.

Before proceeding further we first define the inner and outer region as shown in Fig. 1. In region (defined as ) the TD wavefunction is expanded over the eigenstates of the Hamiltonian matrix representation in the interval . In region (defined as the interval ) the TD wavefunction is represented by its values [] at equidistant grid points .

### iii.1 R-matrix basis-set TDSE in the inner region

In the inner region we define the radial channel functions which are expanded over the R-matrix basis set ( defined in appendix B) as:

(5) |

with the bar at the top indicating that the channel function has been obtained by summing over the radial Hamiltonian eigenstates of the inner region. Note that we have ignored the dependence on the magnetic and spin quantum numbers. From the above definition and the TD wavefunction [Eq. (31)] we obtain in the inner region :

(6) |

The time evolution of the TD wavefunction is now entirely contained in the coefficients . The time evolution of the is determined by the TDSE. However in writing the TDSE we must take care that the Hamiltonian and dipole operators which act on are Hermitian over region (where is only defined). The Hermitian inner-region Hamiltonian is given by and the dipole operator by , where the Bloch surface terms and are set out in Eqs. (38) and (41) respectively. In these circumstances the TDSE over region is written:

(7) |

with . This equation is a key one to the method. The second term on the right hand side compensates for the Bloch terms introduced to make and Hermitian. Note that it makes a contribution only at and brings into play there [Eq. (31)], a wavefunction form which we have defined throughout both regions. This term is central to any time propagation scheme in region because it connects the wavefunction form specific to that region (which may be multi-electronic in a more general formulation) with a wavefunction form that at represents a single electron and which in calculations is obtained from region . We obtain from Eq.(7) the evolution equations for the coefficients by projection over the states :

The quantity is defined as:

(8) |

where is the time-dependent field potential in the Coulomb gauge (see appendix A) and is an angular factor given in Eq. (33c). If the coefficient vector is structured as the inner-region TDSE in matrix notation is as follows:

(9) |

The amplitudes have been defined as in appendix B. The matrix has the block-triangular form of Eq. (36) with the block-diagonal matrices and the lower- and-upper block matrices having matrix elements as:

(10) | ||||

(11) |

where are matrix elements defined in Eq. (42).

### iii.2 Finite-difference TDSE in the outer region

In the external region a grid representation of the TD wavefuction is adopted:

(12) |

with . The time-dependence of the wavefunction is represented by the values of the radial channel functions on a equidistant discretized grid, with . The grid is defined such that and . Furthermore by constructing the vector from the values of the radial channels at the grid points, we obtain a vector of length structured as . The FD representation of the TDSE takes the form:

(13) |

In the FD representation of the time-dependent Hamiltonian the entries and are square matrices of order . The explicit form of these operators depends on the approximation chosen for the derivatives. In the present case, the first and the second derivative of a function are approximated with a 5-point central difference scheme as follows:

(14) |

with chosen so that polynomials of order 4 are differentiated exactly. Given the above, the finite-difference approximation of the diagonal operator in the FD Hamiltonian is,

(15) |

The velocity form of the non-diagonal operator is given by:

The FD form of the TDSE in Eq. (13) is sufficient to propagate the wavefunction in time provided it vanishes at both ends of the spatial grid at all times. This is certainly the case when the FD grid has its innermost point at the origin. In contrast, in the present case, vanishing boundary conditions occur only at the far end of the grid (). More specifically, we assume that for all and this forms the set of boundary conditions imposed on the wavefunction at the far boundary.

Thus some further consideration of the differential operators involved in the FD representation of the TDSE is necessary and we shall shortly see that non-zero function values at an inner boundary bring about contributions from functions values at points below the inner boundary point to the propagation.

We begin by appreciating that since the FD method is a local method, the evaluation of function derivatives at any point relies on function values at neighbouring points, and which of these come into play depends on the approximation chosen for the derivatives, as mentioned earlier. In the FD method the operators are also discretized in a similar way to the functions, i.e. as . The action of a non-derivative operator on a function is trivial, since at the th grid point, but the same is no longer true when operators contain derivatives. Then the rule of differentiation should be given. The central characteristic of the differential operators in the FD method is that values of the wavefunctions at neighbouring points are involved in the calculation of the derivative function. It is then obvious that since the diagonal operators in the finite-difference TDSE [Eq. (15)] involve the second-order differential operator (due to the kinetic term) the complete determination of the requires knowledge of the at points since these enter the determination of second-order derivatives at points and according to Eq. (14). If the propagation is done in the velocity gauge a similar conclusion is reached by considering the non-diagonal operators . The modified form of the TDSE corresponding to a non-vanishing solution on the inner boundary is then

(16) | |||||

where

(17a) | ||||

(17b) | ||||

(17c) |

and are given by,

The elements are given by Eq.(33c) but when the term with is missing and when the term with is also missing from the corresponding equations. The bar on the emphasizes that these radial function values have been evaluated by use of the R-matrix basis set expansion form of the wavefunction in region .

Eq. (16) is the second (and last!) key equation of the method. It does for region what Eq. (7) above did for region . The communication with the solution in region is provided through the terms involving radial function evaluations at two FD points in region immediately inside the boundary with region . Although our detailed exposition above has centred around one-electron wavefunctions throughout both regions, it is clear how the concept embodied in Eq. (16) can be extended to handle a region that is multi-electron in character. The crucial requirement of such a multi-electron inner region is that it must collapse to one-electron character within a few FD points of its outer boundary at . Since in multi-electron R-matrix calculations anyway the inner region must be one-electron in nature by , our additional requirement provides no great extra overhead.

### iii.3 Calculational procedure

Having set out the form of the TDSE in the two regions [Eq. (9)] and [Eq. (16)] we now present briefly the computational procedure involved in the propagation of the wavefunction through one time-step from time to time .

#### Outer region: calculation of :

Assuming at time the wavefunction is known throughout the inner and outer regions we first consider the
outer region TDSE [Eq. (16)]. Although there is a wide variety of methods in the literature we have chosen to employ the standard Taylor propagator as prescribed in Eq. (44). The evaluation of the Taylor series terms requires the quantities which bring into play values of the partial waves
evaluated *in the internal region*
at time [Eq. (17)]. These inner-region partial wave values are formed using Eq. (5).

#### Inner region: calculation of :

In a similar way as done for the outer region, the propagation of the coefficients from time through one time-step to gain their values
at time is now based on the inner-region TDSE in the form of Eq. (9) and the Taylor expansion Eq. (44). For this evaluation
knowledge of the quantity at time is required. The latter quantity includes the
outer-region partial wave and its derivative evaluated *on the boundary *.
Having calculated the coefficients we can immediately form the wave function
according to Eq. (6).

By this stage the wavefunction is known at time throughout regions and and we can proceed further in time by repeating the above procedure for successive time-steps .

### iii.4 Observables within the dual representation

In this section we develop the necessary formulation for the calculation of observables given the different representation used of the time-dependent wavefunction in the inner and outer region (regions and respectively). These representations are given by Eq. (6) and Eq. (12), respectively. Any spatially dependent observable represented by the operator is calculated through the standard formula, which in our case separates into two pieces. To link with the standard experimental setups we assume that any calculation of the observables is performed for times where the external field has vanished. In the following formulas, taking the pulse duration as , we assume the projection time such that . To obtain the population in an eigenstate of the physical system at time , we use the projection operator with the result:

(18) |

with given by Eq. (5) and denoting radial integrations over the inner and outer regions, respectively. Complete information about the final state (ignoring spin variables) is possible by recalling the partial wave expansion of a continuum electron with asymptotic momentum , namely:

(19) |

where defines the direction of the photoelectron with respect to the polarization axis (quantization axis), is normalized on the energy scale and the amplitudes are chosen so that the wavefunction fulfils incoming spherical wave boundary conditions. In the present case, where the ionizing target is hydrogen and (in the following again we drop the dependence) we have with the long-range Coulomb phase shift analytically known Messiah (1999). Therefore the desired angular distribution is obtained through the projection operator which gives:

with the volume element in momentum space. Integration of the above formula over the kinetic energies () results in the photoelectron angular distribution (PAD),

(20) |

while integration over the photoelectron ejection angles provides the angle-integrated photoelectron energy distribution (PES),

(21) |

Finally, further integration over the photoelectron kinetic energies of the last equation results in the total ionization probability (yield) at time as:

(22) |

At this point we have completed the present theoretical formulation leading to the calculation of the most important experimental observables following the interaction of an electromagnetic field with a one-electron atomic target in the dipole approximation.

## Iv Illustrative application to hydrogen

In the present section we apply our approach to the case of ionization of the hydrogen atom by a strong EM field. The reasons we have chosen hydrogen are as follows: (a) it represents the simplest among the atomic systems having just one electron participating in the ionization process, thus being free from complications that may arise from inter-electronic effects in the case of multi-electron systems, (b) angular momentum considerations are reduced to the minimum level where a simple partial wave expansion is adequate to represent the TD wavefunction throughout the electron’s configuration space and (c) last but not least very reliable methods treating one-electron systems Nikolopoulos and Maragakis (2001); Nikolopoulos et al. (2007); Meharg et al. (2005) are at our disposal for a systematic study of the reliability and accuracy issues surrounding the present method. In the present application we have chosen an explicit type time-propagator based on a Taylor expansion [Eq. (43)]. In all the calculations the order of the propagator was and the time step a.u. .

### iv.1 Initial state calculation

We start by calculating the radial function by numerically solving the radial SE for [Eq. (39)] within the inner region. The initial state, made up of and in the inner and outer region, respectively, is then calculated by means of an imaginary time propagation of the field-free versions of Eqs. (9) and (16) with initial conditions:

It is important to emphasize here, that the R-matrix eigenstates do not actually represent the eigenstates of the system, instead they only serve as a complete basis for the representation of the physical state exclusively in the interval . The B-splines basis used consisted of basis functions of order . The knot-sequence is chosen to be linear with discretization step equal to that of the outer region spatial step a.u.. In Fig. 2 we plot the squared amplitude of the radial part of the R-matrix basis for , (black curve) and the state as converged after the imaginary-time field-free propagation. Black-solid and red-dashed curves represent the inner () and outer () region values. The inner boundary has been set at a.u. while the outer boundary at a.u. The R-matrix eigenstate has zero-derivative on the boundary (due to the chosen boundary condition ) while the imaginary time propagation has converged to the state with non-vanishing derivative on the boundary as actually is the case for the ground state of hydrogen. Our initial state has the following form:

(23) |

### iv.2 Real time propagation

Having obtained an accurate initial state [Eq. (23)], through imaginary time propagation, we proceed to the propagation of the TDSE in the outer/inner region in the presence of an external EM field. The EM field chosen was linearly polarized along the -axis with vector potential:

(24) |

where is the field frequency and the pulse duration ( being the number of cycles contained in the pulse and the field period).

To test our approach we calculated the bound state population using three different methods. First, within the present method (BS/FD), we obtained the population in the region by simply summing (at a sufficiently long time ) over all R-matrix eigenstates as:

(25) |

The second method consisted of the standard finite-difference (FD) approach over both regions and i.e. over the whole range (with the same spatial and time step) thereby invoking no division of the electron’s configuration space. Formally, within the present method this is equivalent to setting . We calculated the ionization probability as:

(26) |

where was chosen sufficiently large so that all the outgoing components of the ionized wavepacket were able to travel beyond the chosen distance . With no absorbing potential present we always have for the bound state probability . When an absorbing potential is present then the bound state probability is obtained as:

(27) |

Finally, we performed calculations using a standard basis set (BS) to span the whole range with again no division of configuration space. This is formally equivalent to setting . We obtained the bound state population by summing only over the bound part of the spectrum:

(28) |

In Fig. 3 we show the bound state population of hydrogen as a function of time when the atom is irradiated by a pulse of central frequency a.u. ( eV), total duration of 10 cycles ( a.u.) and peak intensity W/cm. For the BS calculation the bound state population is calculated at the end of the pulse (10 cycles) and no further field-free propagation of the wavefunction is required since the population distribution remains unchanged. In the case of the FD and BS/FD calculations the propagation is extended for a further 10 cycles (field-free propagation) after the end of the pulse until a sufficiently large part of the wavepacket has passed the artificial boundary at a.u.. Given the photon frequency, the hydrogen ionization potential and the rather modest field intensity, we expect the dominant partial wave in the outer region to be the partial wave with the electron’s kinetic energy peaked around a.u. (8.16 eV). By assuming an outgoing wavepacket with central energy of a.u. (thus of velocity a.u.) we can estimate that 10 cycles of field-free propagation is sufficient for our purposes. In connection with this latter point note that this wavepacket travels a distance of 14.5 a.u. in approximately 2.5 field-cycles. This is why the FD and BS/FD bound state populations exhibit a time-delay compared to the BS bound state population. The maximum angular momenta allowed was . Results remained practically unchanged against further increase in angular momentum. We have performed similar calculations with peak intensities W/cm and found similar results with analogous agreement between the BS/FD, FD and BS bound state populations.

In Fig. 4 the values of are plotted as calculated with the present BS/FD and the standard FD method at a.u.. In region (within 14.5 a.u. of the nucleus) the partial wave function () was obtained from Eq. (5). In region (from 14.5 a.u out to 174 a.u.) the values of come directly from the propagation of the outer-region TDSE [Eq. (16)]. Similarly for the FD calculation we obtained by solving Eq. (13) over the whole range . The figure displays excellent agreement between such results from the present (BS/FD) method and the standard FD method. We have chosen to plot only the partial wave since this is the dominant outgoing channel with all other partial wave channels being an order of magnitude lower. This observation simplifies the analysis of the physics involved in the process. We briefly elaborate on this plot. The peak probability for the travelling wavepacket appears around a.u. with a much smaller secondary peak inside region . In an energy representation of the wavepacket, the large peak is associated with the continuum states contribution while the second peak is related to the bound states contribution. Whereas the bound contribution is trapped in the inner region the outgoing component (corresponding to the continuum spectrum) travels a distance of about a.u. which is rather close to the maximum of the wavepacket probability in the plot. We have allowed 15 cycles of travelling time for the wavepacket since significant ionization only takes place around the maximum of the applied pulse which occurs at approximately 5 cycles after the turn-on.

In Fig. 5 the bound state population of hydrogen is shown after exposure to an EM field of central frequency a.u. ( eV), total duration 10 cycles and peak intensity W/cm. Since the photon energy is comparable to the energy gap ( eV) between the ground and first excited states () an appreciable population in these excited states appears at the end of the pulse. At the end of the pulse we obtained from the BS calculation a value for the ground state probability ; a value for the total population in all the excited states ( ) and thus a total bound state probability of [Eq. (28)]. The bound state probability as a function of time is shown in the figure (blue line). We have also performed a FD calculation (with no absorbing boundary present) and calculated the probability within the region a.u. using Eq. (26) and . We chose a box with a.u. to prevent reflection of the wavepacket at the outer boundary over the time interval of interest. In this case the calculated bound state probability for the FD method is given by the green curve (empty cycles). Next, we applied the present method (BS/FD) for a.u. and a.u. To maintain the same accuracy in the calculations in the inner region we increased the number of B-splines basis members to . To compare with the BS and FD calculations we obtained the various probabilities as follows: Black curve (filled squares) in the figure was calculated using Eq. (25) which includes a summation only over those R-matrix eigenstates that have negative energies such that [equivalent to Eq. 28)]. This curve follows closely the bound state probability calculated using the BS method. If in Eq. (25) we include all R-matrix states then the probability enclosed in region is given by the red curve (filled circles) and matches perfectly with the bound state probability from the FD calculation. Similar evaluation through Eq. (26) and with the BS/FD method results in practically the same curve and verifies the equality of results obtained using Eq. (25) and Eq. (26). In other words any increase/decrease of probability within region is matched by an equal decrease/increase of probability within region .

Finally in Fig. 6 we have calculated the photoelectron energy spectrum up to about 20 eV kinetic energy of the ejected electron. In the hydrogen case, although the analytical solutions for the bound and the continuous spectrum are available, the numerical calculation of the eigenstates proves more advantageous for the evaluation of the necessary integrals. In expression (21) we may use an asymptotic expansion for the Coulomb functions Burgess (1963); Madsen et al. (2007) provided that (a) the evaluation is performed at times where the outgoing part of the electron wave packet has travelled sufficiently far away from the residual system (b) the projection operator is constructed either from Coulomb wavefunctions or plane waves depending on the chosen projection time () and (c) the inner-region contribution is ignored since it is only the bound part of the wavepacket that still remains there as time grows. The results can easily be checked by tracing their convergence in time. A detailed discussion of this approach, very well suited to our approach, can be found in ref. Madsen et al. (2007). The solid black curve represents the result of a BS calculation while the dashed red curve the result of the present calculation. Had we used a larger box and a finer mesh for the outer region we would be able to calculate even higher (in energy) the corresponding PES for a full comparison with the BS calculation, but this is not the purpose of the present work.

## V Conclusion and perspectives

In conclusion a new *ab initio* time-dependent method for the treatment of the single-electron ionization of atomic
and molecular systems under an external electromagnetic field has been set out. It has been developed in detail for systems that are single-electron throughout
and applied to the simplest case, namely the hydrogen
atom. The method is based on the division of the configuration space of the ejected electron into two regions and .
In region (which may be multi-electronic) the time-dependent wavefunction is expanded on the basis of R-matrix eigenstates and propagated
through the time evolution of the expansion coefficients.
In region a grid-representation of the time-dependent
wavefunction is adopted and a finite-difference technique is employed for the representation of the operators. In both
regions the chosen time-propagator in illustrative calculations is a high-order explicit Taylor propagator. The key point in the present method is
the time-dependent matching conditions that the inner (region ) and outer (region ) wavefunctions
should simultaneously satisfy at each time step. Although these matching conditions have been developed here
for an explicit time-propagator, the methodology can also be applied for implicit time propagators. The present work
represents an important step towards the implementation of such a methodology in multi-electron systems
(atomic and molecular) where the full advantage of the R-matrix technique can be taken into account.
The straightforward extension of the present approach to the case of a truly multi-electron system is discussed in
Sec. II and is currently the subject of our work. In addition to our fundamental interest in gaining
an ab-initio description of multi-electron systems under strong laser fields the present work is mainly motivated
by the development of sources of short wavelength laser light residing well into the VUV or soft XUV regime (HOHG/FEL sources).
In contrast to long wavelength laser light, the light from such sources tends to interact directly with more than just a single electron
and is able to probe directly the innermost electrons of multi-electrons systems
thus making the development of new suitable theoretical methods a necessary and formidable task.

The authors gratefully acknowledge discussions with Dr. Michael Lysaght, Dr. Hugo Van der Hart and Prof. P. G. Burke. The present work has been supported through a Marie Curie Intra-European Fellowship under contract TDRMX-040766 awarded to LAAN and by the UK Engineering and Physical Sciences Research Council.

## Appendix A TDSE of single-active-electron atomic systems over a spherical harmonic basis

The field-free SAE Hamiltonian reads,

(29) |

with the potential equal to for a purely hydrogenic system (of atomic number). Alternatively could be constructed as a model or Hartree-Fock potential. The TDSE of the system in an external time-dependent radiation field, is written as:

(30) |

with the system wavefunction and the interaction operator between the system and the external field, in the dipole approximation. In our present numerical implementation we choose a spherical coordinate system for the active electron. We represent the angular variables in a basis of spherical harmonics and write the wave function as,

(31) |

where the spin-variables of the wavefunctions are ignored. In an actual calculation we must truncate the spherical harmonics expansion at some maximum value . In the remaining formulas we abbreviate the truncated double summation by .

The time-propagation of the wave function proceeds in spherical coordinates as follows. Substituting Eq. (31) in Eq. (30) and projecting onto the spherical harmonic basis we obtain the following coupled differential equations for the radial channel functions as,

(32) |

For the special case of linearly polarized light along the -axis and in dipole approximation the radial time-evolution operators are given by,

(33a) | ||||

(33b) | ||||

(33c) |

with and the radial dipole operator. The time-dependent radial dipole operator is given as,

(34) |

in the velocity form. The quantity represents the field potential in the Coulomb gauge. Within the present context the interaction operator couples atomic states of equal magnetic quantum number, hence we drop the dependence on in the subsequent formulation.

By properly arranging the radial channel functions according to their angular momentum label we form the radial vector wavefunction . In this case the matrix representation of the TDSE [Eq. (32)] is written as,

(35) |

where and

(36) |

## Appendix B R-matrix eigenstates in the inner region

### b.1 Hamiltonian operator in the inner region

In the inner region the radial wavefunctions are expanded over the eigenstates of the radial Hamiltonian:

(37) |

with the radial Bloch operator,

(38) |

and given by Eq. (29). The eigenstates of the R-matrix Hamiltonian operator are uniquely determined if we set the boundary conditions needed to be fullfiled at the boundaries and . In the present case on physical considerations we take all solutions to vanish at the origin while at the solutions take non-vanishing values. This choice makes the radial R-matrix operator Hermitian over the inner region . Therefore for each value of the angular momentum we solve the following eigenvalue problem:

(39) |

where is an integer labelling the eigenstate. The above eigenvalue differential equation is transformed to solving a matrix diagonalization problem by employing a B-spline basis set of size , order de Boor (1978) for the representation of the solutions in region :

(40) |

In the expansion the first B-spline [] is excluded in order to conform to the boundary condition at the origin . Note that by definition of the B-splines the amplitude of the eigenstates on the boundary [] is simply the last coefficient in the expansion, namely, . All required integrals are evaluated, with the Gaussian quadrature rule, to machine accuracy.

For each partial wave the solutions constitute an orthonormal basis with members,

with real eigenvalues .

### b.2 Dipole operator in the inner region

While the velocity form of the radial dipole operator includes a first-order derivative term [Eq. (34)] which taken together with the non-vanishing values of the eigenstates at the boundary makes it non-hermitian.We can make this operator Hermitian by adding the dipole Bloch operator for the first-order derivative in a similar way as done for the field-free Hamiltonian . Thus if we define the dipole velocity operator in region as:

(41) |

we find for the radial velocity operator: