Photoionisation of Rubidium in strong laser fields

The photoionisation of rubidium in strong infra-red laser fields based on ab initio calculations was investigated. The bound and the continuum states are described with Slater orbitals and Coulomb wave packets, respectively. The bound state spectra were calculated with the variational method and we found it reproduced the experimental data within a few percent accuracy. Using the similar approach, ionisation of Rb was also successfully investigated. The effects of the shape and the parameters of the pulse to the photoionisation probabilities and the energy spectrum of the ionised electron are shown. These calculations may provide a valuable contribution at the design of laser and plasma based novel accelerators, the CERN AWAKE experiment.


Introduction
The study of the ultrafast dynamics of electrons in intense short light pulses is a hot topic nowadays [1]. Multiphoton photoionisation of one-, two-or many-electron atoms in intense laser field has been extensively studied by various non-perturbative quantum mechanical methods by different groups in the last decades [2,3,4,5,6]. One of them is the time-dependent close coupling (TDCC) method originally developed by Bray. The details can be found in the review of [7]. Another successful approach, the timedependent close coupling method on a two-dimensional finite lattice was introduced by M.S. Pindzola and F. Robicheaux [8]. Various R-Matrix matrix methods are also feasible to perform ab-initio photoionisation calculations [9,10].
As basis set expansions there are two popular ways. The first is the application of b-Splines by Bachau [11] and the other is the Sturmian basis set [12].
The problem of rubidium photoionizaton is almost ninety years old. The first experiment was performed by Lawrence [13] where the wavelength of the applied light was in the range of 220 to 313 nm. Multiphoton ionizaton of Rb atoms was also studied with the help of a tunable dyelaser over the wavelength of 460 to 650 nm by Collins [14]. Tamura and his group [15] were the first who measured the two-step selective photoionizatinon of Rb with Ti:sapphire solid state laser. Experimental determinaton of the photoelectron angular distribution of rubidium atoms in linearly and elliptically polarized lights were investigated by Wang and Elliott [16]. Courtade et al calculated two-photon ionization of cold Rb atoms with a near resonant intermediate state. [17].
Here we study the ionizaton of Rb atom problem with two different methods: full-fledged ab inito quantum mechanical method and classical trajectory Monte Carlo (CTMC) method. In the former one, computationally we realized a time-dependent coupled-channel method where the wave functions of the channels are constructed with Slater orbitals and regular Coulomb wave packets with equdistant finite widths in energy. This basis set was originally introduced to study heavy-ion and He atom collision in 2002 by Barna [18] later become sophisticated and culminated to describe the angular distribution in two-photon double ionisation of helium by intense attosecond soft-x-ray pulses [19].
In this work we apply both TDCC and CTMC method to calculate the photoionisation probability of Rb interacting with strong laser field. We analyse how the photoionisation probabilities depend on the pulse parameters. We also calculate the photoelectron energy spectrum at different laser intensities with the TDCC method.

Ab initio calculations
In our calculations we studied the photoionisation phenomena of rubidium via ab initio calculations. The rubidium atom was modelled with a single active electron model with a frozen atomic core. Considering the timedependent Schrödinger-equation (TDSE): The Hamiltonian operator has the form of whereĤ Rb andV I describe the free rubidium atom and the interaction with the external laser field, respectively. There are many approaches to model the free rubidium atom, or in more general, an alkali atom within the confines of the one active electron approach [20,21,22]. When choosing such a model potential, the most important criterion is if the potential reproduces the experimental energy spectrum of the free atom. The Hellmmann pseudopotential, i.e. the second term of (3) fulfils this criterion and also has the advantage that many matrix elements can be calculated fully analytically. It also provides a graphical point of view to our physical approach: the inner shell electrons shield the electric charge of the atomic core, hence the valence electron feels a reduced charge from the core: The shielding parameters depend on the structure of a specific alkali atom. Ref. [23] lists the parameters b and d for every alkali atoms. For rubidium, b = 4.5 and d = 1.09993. The solution of the time dependent Schrödinger equation can be expanded in terms of the eigenstates of the time-independent Schrödinger equation: We apply the following Ansatz for the TDSE: Inserting (5) into (1), using (4), we get: (6) Multiply the above equation by Φ k (r)| e iE k t . We get the following system of equations for the a j (t) coefficients: In (7), E kj := E k − E j and V kj := Φ k (r)|V I | Φ j (r) is the couplings matrix. The (highly) oscillatory term can be transformed out by introducingã Inserting (8) into (7), we get: In our calculations the initial conditions have been set up such that the valence electron is in the ground state: We should note, however, that this is not a general requirement. Any mixed state could be specified as well.
Integrating this system of ordinary differential equations on the time interval of the interaction, we get the wave function of the final state. Relevant physical quantities, e.g. the occupation probabilities of either the bound or the continuum states can be calculated from this wave function: with P k denoting the occupation probability of the k-th state. The energy spectrum is defined by with Φ l E (r) being an arbitrary continuum state with energy E and azimuthal quantum number l and Ψ (t = T, r) being the final state wave function i.e. the (numerical) solution of (9).

Description of the bound and the continuum states
As stated above, the wavefunction of the valence electron is expanded on the basis of the eigenfunctions of the free Hamiltonian operator. It depends on the problem which basis is a sufficiently good choice. We describe the bound and the continuum states of the valence electron with Slater-type orbitals (13) and Coulomb wavepackets (15), respectively [18]: with n, l, m being the principal, azimuthal and magnetic quantum numbers, respectively, and κ the screening constant, which is, in our case, a variational parameter that specifies the energy of a bound state. The Slater orbitals form a normed, but not orthogonal basis, with the normalization factor The Coulomb wave packets discretize the continuum via integrating the Coulomb wave functions (17) on a finite energy (momentum) interval. Therefore, they form the probability amplitude of an electron being in the state with its energy lying between E − ∆E and E + ∆E. The energy and the width of such a state is given by the integration limits, i.e. k and ∆k.
with k and ∆k being the center and the width of the covered momentum range, l and m the azimuthal and magnetic quantum numbers, respectively, andZ the charge of the ion,Z = 1 for a singly ionized atom. The Coulomb wave packets form an orthonormal basis such that if their corresponding energy ranges do not overlap, then their overlap integral is zero. The normalization factor reads: Finally, the Coulomb wave function has the form of For further details about the Coulomb wavefunctions, consult [24]. Sometimes, as it is in our case, it is more natural to specify the Coulomb wave packets with their energy range instead of the corresponding momentum range. There is a simple connection between the parameters of the momentum and the corresponding energy range: We chose the parameters of the Coulomb wave packets such that they subdivide the total energy range equidistantly.

Final formula for the energy spectra
Determining the energy spectrum of the bound states leads to a variational problem. In general, the basis functions describing the bound states should contain at least one free parameter. In our case, every Slater orbital has a free parameter, the κ screening constant, which will be determined via minimizing the energy functional corresponding to our physical system. By doing so, one gets a generalized eigenvalue problem: with and being the Hamiltonian and the overlap matrices and E the energy of a bound state, which is a generalized eigenvalue and c denotes a generalized eigenvector. Here ψ j can refer either to a Slater-function or to a Coulomb wave packet.
Having M states, an energy eigenstate of the free Hamiltonian operator with energy E can be expanded on the chosen basis, the expansion coefficients being the components of the corresponding generalized eigenvector:

Interaction with the external laser field
In the configuration-interaction approach, introduced at the beginning of this section, every interaction with the external environment is incorporated into the couplings matrix. Its exact shape depends on the kind of the interaction, in our case, on the shape of the external laser field. Considering electromagnetic interaction, it is easy to show that the couplings matrix is proportional to the dipole matrix: In length gauge, the interaction operator has the form of Since in the present study we are investigating the interaction of a single atom and the external laser field, it is clear that the distances characterizing an atom are much smaller than the wavelength of the laser field. Therefore, the dipole approximation is valid. The electric field has the form of with ε being the polarization vector, E 0 the amplitude and f (t) an arbitrary function that is feasible for modelling a laser pulse. The dipole matrix elements are therefore: Let d denote the dipole matrices corresponding the basis functions: d pq := e ψ q (r) |rε ε ε| ψ p (r) .
Using this notation, we get a compact formula for the dipole matrix elements: We chose the electric field such that it is polarized into the z direction and has a sine-square envelope. The most common choice for the envelope is a Gaussian function, however, the sine-square envelope has a slight advantage at numerical calculations: such a pulse results in a function with compact support.
with T being the pulse duration and ω L being the initial laser frequency.

Classical trajectory Monte Carlo method
A classical trajectory Monte-Carlo method is used to calculate the ionisation probabilities of Rb in intense laser fields. The Rb atom is characterised as two body system with an active electron and the remaining core. The interaction between the active electron and the target core is described either by a simple Coulomb potential as in the conventional CTMC method (model 1) or with the a model potential (model 2). In both versions of the present CTMC approach, Newton's classical non-relativistic equations of motions for a two-body system are solved numerically for a statistically large number of trajectories for given initial parameters. The equations of motion were integrated with respect to time as an independent variable by the standard Runge-Kutta method. Atomic units were used throughout the calculations. In this calculation, the total number of recorded trajectories was 1 · 10 5 .

Model 1
The CTMC method is applicable to hydrogen-like target atoms (i.e., A (z−1)+ ions, in general) in a natural manner. An extension to this picture requires an effective charge Z eff of the target nucleus as seen by the active electron: The effective charge of 2.2 and corresponding binding energy of 0.154 a.u. were used for Rb(5s).

Model 2
The potential of the Rb ion is represented by a central model potential developed by Green [20] which was based on Hartree-Fock calculations: where Z is the nuclear charge and Using the energy minimization, Garvey et al. [22] obtained the following parameters for Rb: H = 4.494 a.u. and d = 0.777 a.u. The initialization parameters of Rb were selected as described by Reinhold and Falcon [25] developed for non-Coulombic systems. The potential from Eq. (31) is used to describe the interaction between the projectile and the target electron with the Rb ion core. The initial state of the target is characterized by a micro-canonical ensemble, which is constrained to an initial binding energy of 0.154 a.u., at a relatively large distance from the collision center, choosing the initial parameters randomly. The distance between the projectile and the target was large enough that the interaction with the target was negligible.

Results
In our calculations we modeled the rubidium atom by including the first 35 bound states. According to our the theoretical model we approximated the experimental energy levels by solving the corresponding variational problem that ended up in a generalized eigenvalue problem. By finding the proper values of the screening constants, the generalized eigenvalues of the Hamiltonian matrix corresponding to the bound states provide a sufficiently good approximation of the experimental energy levels. We managed to solve this optimization problem by using genetic algorithm (see e.g. Ref. [26]). For the corresponding linear algebraic calculations we used the Armadillo library [27]. The calculated energy levels of the bound states, compared with the experimental are listed in tables 1, 2, 3 and 4, respectively. The experimental data has been taken from [28]. E exp and E calc denote the experimental and the calculated energy eigenvalues, respectively.
Then, we defined an energy range in the continuum. We chose 0 − 15 eV and subdivided it into 400 parts. The Coulomb packets with different azimuthal quantum numbers have been constructed to all energy levels accordingly, the azimuthal quantum numbers lying in the range l = 0 − 3 with m = 0. Since we included 35 bound states and there are 400 energy levels for every azimuthal quantum number in the continuum, we have a total of 1635 states in our model. For the first glance the energy resolution of the continuum spectra may seem to be too fine. However, during our calculations we saw that this is the resolution where the spectra converge to the same peak structures even at high intensities. That is, less fine resolution would have resulted in spectra with high numerical noise, whereas higher resolution would have been redundant. At this point we mention that choosing the proper energy resolution is a crucial point in similar calculations as the runtime of a single simulation grows rapidly as the intensity of the applied laser field grows. On the other hand, the runtime scales with N 2 with N denoting the total number of states. At the highest applied intensity, using 20 processor cores took for almost two months.  Table 3. Energy spectrum of the bound states for l = 2. The experimental values have been taken from Ref. [28].
u.] 4f -0.0314329 -0.031225 5f -0.0201073 -0.0199783 6f -0.0139554 -0.0138734 7f -0.0102476 -0.0101937 8f -0.00784234 -0.00780359 Table 4. Energy spectrum of the bound states for l = 3. The experimental values have been taken from Ref. [28]. We also took care of solving the coupled channel equations with the least numerical noise as possible. Since the pulse length is quite long, we decided to not to apply the common fourth order Runge-Kutta method when integrating our ODE system. Instead, we applied the Bulirsch-Stoer method with an eighth order embedded Runge-Kutta method as controller method. We set both the absolute and relative error tolerances to 10 −8 . This was sufficientand also required-to preserve the unity norm of the wave function, i.e. to not to numerically violate the unitarity. In our calculations the difference of the norm of the wave function from unity was not higher than 10 −5 .
After theses preparations, the spectrum and other quantities have been calculated such that the intensities lie in the range 10 12 W/cm 2 ≤ I < 10 14 W/cm 2 with λ = 800 nm wavelength and τ = 120 fs pulse duration. The parameters have been chosen such that they fit to the CERN-AWAKE experiment [29]. Note that the Keldyshparameter runs from 0.5916 to 5.916 if the intensities lie between 10 14 W · cm −2 and 10 12 W · cm −2 .
According to Morales et. al, [30], stabilisation of the total ionisation probability is expected (see Fig. 1). This means that as the laser intensity increases, the total ionisation probability does not grow steadily, but instead, it either saturates at a given level, (Fig. 1), or reaches its maximum at a critical point, then drops slightly and remains constant at higher intensities (Fig. 2). The critical point in the latter case is also referred as the bottom of the "Death Valley". Fig. 2 shows the photoionisation probabilities at different intensities with λ = 800 nm and τ = 120 fs. The blue line corresponds to the ab initio calculations, the orange and the green to CTMC simulations with a classical Coulomb potential with Slater screening and Garvey model potential, respectively. According to the classical suppression [31,32], the CTMC simulations predict lower ionisation probabilities at lower intensities, but they also predict that the ionisation probability is going to be saturated as the laser intensity is growing. However, the dis- crepancy between the classical and quantum mechanical calculations is greater than expected. This difference could be corrected by fine-tuning the model potential used in the CTMC simulations. The ab initio calculations suggest that the ionisation probability has been already saturated even at low intensities. This is a pleasant result since at the intensities planned to be applied in the CERN-AWAKE experiment practically fully ionisation occurs. This is a necessary condition for the experiment to work properly.
On Fig. 2 one can see that the classical suppression is larger than expected. Our idea is that this difference could be decreased by choosing a better model potential in the CTMC simulations When should also note that in our case the stabilisation occurs at a much lower intensity than for potassium, even though the ionisation potentials of rubidium and potassium are quite similar (4.177 eV and 4.341 eV, respectively). The reason behind this that, in our model, the length of the laser pulse (τ = 120 fs) is almost the double of the one used in Ref.
[30] (τ = 65 fs). Therefore, the photon absorption is more "efficient" both in the bound states and the continuum states, resulting in a lower saturation intensity.
Photoelectron energy spectra have been calculated as well, see Fig. 3. Note that at relatively low intensities there is a strong peak near 0 eV. This means that at 10 12 W/cm 2 most of the electrons have relative low kinetic energy, that suggests that the plasma applied in the CERN-AWAKE experiment can have uniform electron density. At higher intensities none of the peaks is dominant, but all the spectra have a mutual characteristic. The difference between the abscissae of two neighbouring peaks is approximately 1.55 eV, i.e. one photon energy at λ = 800 eV wavelength. This is a fingerprint of above-threshold ionisation (ATI) at these intensities. The width of these peaks originates from the non-resonant excitation during the ionisation process. Ref. [33] provides an excellent review about this phenomenon.

Conclusion
In our work photoionisation of rubidium with strong laser pulses has been studied. We applied two different techniques: time-dependent close coupling and classical trajectory Monte Carlo method. In the former one, the continuum states have been described with Slater-type orbitals and the continuum states have been approximated with Coulomb wave packets. First we reconstructed the eigenfunctions of the free Hamiltonian operator on our basis set and found that the generalized energy eigenvalues of the Hamiltonian operator are in a good agreement with the experimental values. Then we investigated the dependence of the photoionisation probability and the photoelectron energy spectrum on the laser intensity. We chose the laser parameters such that they are similar to the corresponding ones in the CERN-AWAKE experiment. Since this parameter range is very specific, we could not directly compare our calculations with experimental data. However, our calculations reproduced the expected saturation of the photoionisation probability as a function of laser intensity.The ab initio calculations predicted that staturation occurs at moderate intensities, near 10 12 W/cm 2 , whilst CTMC simulation predicted this phenomena only near 10 14 W/cm 2 . Both methods agree qualitatively with each other and also with the presented simulation of Morales et al. [30]. We also calculated the photoelectron energy spectrum at different laser intensities. The ATI peaks are clearly visible and it can also be seen that the average kinetic energy of the continuum electron grows monotonically with the laser intensity.
For a better understanding of the physics of photoionisation of Rb further work is in progress. We plan to calculate photoelectron angular distributions and to investigate the additional effects regarding the variation of additional laser parameters like pulse duration or wavelength. We also plan to study the phenomena related to non-zero chirp.

Acknowledgement
The ELI-ALPS project (GINOP-2.3.6-15-2015-00001) is supported by the European Union and co-financed by the European Regional Development Fund.One of us K.

Statement
This article is a basis of Author M.A. Pocsai's Ph. D Thesis. He is responsible for preforming and evaluating the ab-inito calculations, including implementing the corresponding software, written in C++11. He is primarily responsible for writing the manuscript. Author I.F. Barna is the supervisor of M.A. Pocsai. He also participated in writing and reviewing the manuscript. Author K. Tőkési is primarily responsible for performing and interpreting the CTMC calculations that have been compared with the abinitio results. Furthermore, he participated in writing and reviewing the manuscript.