Absolute strong-field ionization probabilities of ultracold rubidium atoms

We report on precise measurements of absolute nonlinear ionization probabilities obtained by exposing optically trapped ultracold rubidium atoms to the field of an ultrashort laser pulse in the intensity range of $1 \times 10^{11}$ to $4 \times 10^{13}$ W/cm$^2$. The experimental data are in perfect agreement with ab-initio theory, based on solving the time-dependent Schr\"odinger equation without any free parameters. Ultracold targets allow to retrieve absolute probabilities since ionized atoms become apparent as a local vacancy imprinted into the target density, which is recorded simultaneously. We study the strong-field response of $^{87}$Rb atoms at two different wavelengths representing non-resonant and resonant processes in the demanding regime where the Keldysh parameter is close to unity.

Ultracold atoms serve as ideal systems for precise studies of light-matter interaction with a high degree of control over critical parameters. The creation of charged particles in form of ions and electrons out of these welldefined ensembles offers appealing possibilities for fundamental research as well as advanced applications. Hybrid atom-ion quantum systems benefit from long-range interactions introduced by the Coulomb potential [1,2] and are promising candidates for a quantum information platform [3,4]. Moreover, ions immersed in a neutral ultracold ensemble can form new states of matter such as mesoscopic molecular ions [5,6] and ultracold plasma [7][8][9].
Ultrashort laser pulses of durations in the femtosecond domain offer a mechanism for instantaneous creation of cold ions and electrons in ultracold environments with respect to the time-scales of the atomic motion in the trap. Experiments combining ultracold targets with ultrashort laser pulses are able to observe detailed space-charge dynamics [10,11] and provide ultracold electron bunches exhibiting a large degree of coherence for the use in highbrilliance accelerator sources [12,13]. The ionization is achieved at high peak intensities of the laser field either through absorption of multiple photons or a substantial modification of the Coulomb potential that leads to tunneling [14,15].
For all of the above mentioned examples of cold atom/ion and electron systems it is vital to understand the nature of the ionization process in detail since the subsequent dynamics strongly depend on the ionization regime. The high degree of controllability on the atomic side eventually has to be extended also to the process of * pwessels@physnet.uni-hamburg.de; http://www.cui.unihamburg.de/ creating electrons and ions on instantaneous time-scales with ultrashort laser pulses.
Traditionally, the ionization of atoms in strong light fields is categorized into intuitive regimes distinguished by the Keldysh adiabaticity parameter [16] relating the ionization potential E I to the average quiver energy of a free electron with mass m e and charge e. The latter is described by the ponderomotive potential U p generated by the driving laser field of angular frequency ω and intensity I. A perturbative multiphoton model (γ 1) is contrasted by an adiabatic tunnel ionization [17][18][19] regime (γ 1) that leads to over-the-barrier ionization for increasing laser intensities [14,15].
So far, the ionization process in strong laser fields has been studied mainly with rare gas atoms. For alkali atoms used in most cold atom experiments the transition regime separating multiphoton from tunnel ionization at a Keldysh parameter around γ ≈ 1 is reached at considerably lower intensities. Experiments accessing the strong-field regime with alkali atoms [20][21][22] have revealed that commonly used strong-field ionization models in form of a pure multiphoton-or tunnel ionization process can not be applied a priori which makes employing ab-initio calculations by solving the time-dependent Schrödinger equation (TDSE) inevitable.
An accurate comparison of theoretical models to experimental data sets is notoriously difficult since the ionization of atoms in strong laser fields depends nonlinearly on a variety of experimental parameters which are not always accessible or under control. Often, the intensity calibration is even performed using the theory under test, for example in advanced electron spectroscopy experiments [23,24].
Measuring absolute ionization probabilities is in general even more challenging due to the necessity of precise knowledge on the incoming photon-flux, target density as well as detector efficiencies. Ionization cross-sections can be accurately measured using magneto-optically trapped (MOT) atoms and observing the trap loss rates before and after introducing an ionizing laser beam [25][26][27]. This technique is very sensitive and can measure low ionization rates. However, it relies on calibrating the excited state fraction because of the intrinsic presence of MOT laser beams. Additionally, a focused laser beam always consists of a distribution of intensities, a phenomenon usually referred to as focal averaging that has to be taken into account.
We present an alternative methodology to overcome these obstacles by employing optically trapped, ultracold 87 Rb ensembles at nanokelvin temperatures. We observe local losses imprinted as vacancies into the simultaneously detectable atomic density profile which leads to a spatially resolved loss signal. This allows us to extract absolute ionization probabilities. We employ this technique to precisely measure strong-field losses in 87 Rb atoms and compare the measurements to the results obtained using ionization probabilities calculated by numerically solving the TDSE. We find an excellent agreement between experiment and theory without any free parameters. The sole outer-shell electron in alkali atoms establishes a reference system for studying strong-field ionization in the transitional γ ≈ 1 regime.

Measuring absolute local losses
In order to achieve nanokelvin temperatures, 87 Rb captured and laser-cooled in a 2D/3D MOT is further cooled by transferring the atoms into a hybrid trap [28] which consists of a magnetic quadrupole field and a crossed optical dipole trap at 1064 nm wavelength. This realizes a well defined sample of ground state 87 Rb atoms at final temperatures around 100 nK close to quantum degeneracy (for details see Methods section). As presented in Fig. 1, a single femtosecond laser pulse delivered by a pulse-picker is focused onto the trapped atoms where it locally ionizes 87 Rb atoms in the cloud. The remaining atomic density is mapped by resonant absorption imaging at the 87 Rb D 2 transition (780 nm wavelength) by a 50 µs long imaging pulse. The focal intensity distribution is imaged and the pulse energy is recorded simultaneously to accurately characterize the intensity in the focal volume for calibrated peak intensities.
This experiment is performed at two different wavelengths representing non-resonant and resonant ionization processes. The non-resonant case is realized at a wavelength of 511 nm with pulses of 215 +20 −15 fs FWHM duration while for the resonant case pulses of 1022 nm wavelength are used at a pulse duration of 300 +30 −23 fs FWHM. In both cases the femtosecond laser beam is focused to a waist of w 0 ≈ 10 µm using a lens of f = 200 mm focal length (for exact numbers and details see Methods section).
Data evaluation and comparison to theory Figure 2a shows an in-situ absorption image of the atomic density exposed to a single femtosecond laser pulse of I 0 = 6.9 +0.7 −0.8 × 10 12 W/cm 2 peak intensity. The ionization losses lead to a vacancy imprinted into the atomic density profile. In contrast to typical MOT temperatures of approximately 100 µK, here the motion of atoms is negligible during imaging so that we have access to the local losses in the focal intensity distribution and simultaneously measure the line-of-sight integrated twodimensional density profile of the neutral atoms. Most other experiments are limited to an integrated signal over the focal volume as sole observable and can not access the local density profile. In order to quantify the fraction of lost atoms, the difference of two Gaussian functions is fitted to the pixel row sum of the image in Fig. 2c where A 0 determines an offset while A c , x c and σ c correspond to amplitude, central position and width of the atomic cloud and A v , x v and σ v describe the same parameters for the vacancy, respectively. The loss fraction Quantifying ionization losses in an ultracold atomic ensemble. a Measured and b calculated in-situ density distribution of ultracold 87 Rb atoms after applying a focused femtosecond laser pulse of 511 nm wavelength and 6.9 +0.7 −0.8 × 10 12 W/cm 2 peak intensity. c The row sum profiles for the measured (gray crosses) and simulated (red dashed line) density images are used to quantify the measured and simulated fraction of ionized atoms by applying a double Gaussian fit (blue dotted line for measured profile).
relating the number of ionized atoms N i to the total initial number of atoms N 0 can be derived by the fit parameters as shown in Eq. (3).
In order to compare the measurements to theoretically obtained ionization probabilities, the predicted lossfraction due to photoionization is calculated. First, the probability P of an atom to be ionized after the laser pulse of peak intensity I 0 is derived from numerical TDSE calculations. Using these results, the distribution of peak intensities within the laser focus I 0 (x, y, z) is mapped to spatially resolved ionization probabilities P (x, y, z) which are superimposed with the 3D density map of the atomic sample for each grid unit (for details see Methods section). This results in a simulated 3D density distribution of the remaining, non-ionized atoms that is projected onto the image plane of the camera and convoluted with the imaging resolution as shown in Fig. 2b. The resulting simulated 2D projection is evaluated in the same manner as the absorption image. A remarkable agreement is obtained between the calculated and measured row sum as depicted in Fig. 2c. The simulated loss fraction f sim l can be extracted from the fit to the computed profile and compared to the experimentally obtained value.
Non-resonant two-photon ionization at 511 nm The measured fraction of lost atoms f l according to Eq. (3) with respect to the peak intensity I 0 is collected in Fig. 3b for 511 nm wavelength. The corresponding number of generated ions N i is computed by multiplying the fraction of lost atoms with the initial number of atoms in the cloud and is depicted on the right axis. The increasing measured loss fraction (blue dots) originates from higher ionization probabilities P at higher peak intensities presented in Fig. 3a. Once the ionization probability is close to unity in the focal center (here at 5 × 10 12 W/cm 2 ), losses still increase since the volume of high ionization probabilities expands over the atomic cloud.
The calculated loss fraction f sim l expected by strongfield ionization probabilities obtained from ab-initio TDSE calculations is shown as red dashed line in Fig. 3b and perfectly reproduces the experimental result. It should be emphasized that this curve neither has free parameters to fit to the data points, nor is it scaled or shifted so that the absolute ionization yield of the focal intensity distribution can be directly compared to the measurement because we image local losses in a simultaneously observable atomic density profile.
Although the TDSE method uses no other approximation than the single active electron approximation, it is still helpful to compare our measurements to widely used strong-field ionization models to gain intuition regarding the applicability of these models as well as the sensitivity of our measurements to deviations from the TDSE ionization probabilities. For these reasons we now compare our datasets to common strong-field ionization models.
In a multiphoton ionization process at 511 nm wavelength, the rubidium atom has to absorb two photons. For reference, the expected ionization probability described by lowest order perturbation theory using a multiphoton ionization (MPI) cross-section of σ 2 = 1.5 × 10 −49 cm 4 s at the given wavelength [29] is included in Fig. 3a as green solid line. This model agrees surprisingly well with the ab-initio results over the whole intensity range even though the ionization probability as well as the Keldysh parameter are close to unity.
In contrast to rare gas atoms, the electric field that classically suppresses the Coulomb barrier of the atom to free the valence electron is exceeded at much lower intensities in alkali atoms due to the low binding energy. For rubidium this over-the-barrier (OBI) intensity is reached at I OBI = 1.22 × 10 12 W/cm 2 (vertical black dashed line). Previous experiments [20,22] and theoretical studies [30,31] have already mentioned the peculiar situation that in alkali atoms the OBI intensity is reached well within the multiphoton regime, in particular for rubidium at a Keldysh parameter of γ = 8.4 for an optical wavelength of 511 nm (compare Fig. 3a). Assuming that the ionization probability is unity once the OBI intensity is exceeded [32] and is given by the multiphoton ionization cross section otherwise, we can test for contributions of over-the-barrier ionization processes. The expected loss fraction in this case significantly overestimates the observed loss fraction (purple lines in Fig. 3).
Although the Keldysh parameter is still bigger than (but close to) unity, we can also test the tunnel-ionization Rb. a Keldysh parameter γ and calculated ionization probabilities P after a single laser pulse for different ionization models and b measured fraction of lost atoms f l (blue dots) for laser pulses of 511 nm wavelength. The nominal over-the-barrier intensity IOBI is depicted as vertical black dashed line. The data is compared to the calculated loss fraction f sim l using ab-initio TDSE calculations (red) and a twophoton ionization cross section of σ2 = 1.5 × 10 −49 cm 4 s [29] (green). The red dotted lines indicate the I 2 0 scaling of a two-photon process. The additional assumption in the calculation that the electron is free once IOBI is exceeded leads to an overestimation of the losses (purple lines) just as including tunnel-ionization losses using pure ADK theory (yellow lines). The inset in b presents the same data in linear instead of log/log scaling. model by calculating ionization probabilities obtained by the Ammosov-Delone-Krainov (ADK) theory [18] that is strictly valid only for Keldysh parameters much smaller than unity (γ 1). Details on this calculations are given in the Supplemental Material. Including this loss channel into the calculation (yellow lines in Fig. 3) still clearly overestimates the observed ionization yield. Regarding the sensitivity of our method it should be noted that the small intensity range from 1.5 to 4 × 10 12 W/cm 2 where the ADK ionization probability exceeds the MPI proba-bility, leads to an experimentally expected loss fraction that lies significantly outside the error bars of the measurements so that our technique is sensitive to even subtle details in the ionization probability curve.
We observe considerable deviations from all conventional adiabatic models which is expected from a Keldysh parameter still bigger than unity. This means that the modification of the Coulomb barrier of the atom cannot be regarded as adiabatic which also renders the concepts of an over-the-barrier ionization intensity not applicable in this case. Consequently, one has to rely on calculating the ionization process with TDSE methods. Our results underline the dominance of multiphoton ionization in alkali atoms consistent with previous experiments [20][21][22] and theory [30,31].  scaling is plotted along with the theoretical results as the most compatible scaling law with the calculated ionization probability curve.
The experimental data points in Fig. 4b also agree extremely well with the expected loss fraction from the TDSE calculations in view of the fact that there is no free parameter and no arbitrary scaling involved. Compared to the resonant measurement the focal intensity distribution here was not perfectly Gaussian. Figure 7 depicts the measured foci for both wavelengths obtained by the focal diagnostics (see also Fig. 1). While the focus of the 511 nm laser pulse is well represented by a single 2D Gaussian fit, the focus image at 1022 nm wavelength features additional weak intensity regions next to the central main intensity peak. In the calculations of the loss fraction with the TDSE ionization probabilities we have accounted for this by approximating the measured intensity distribution by a superposition of three Gaussian profiles. The reconstruction is discussed in the Methods section and also presented in Fig. 7. The remaining minor discrepancy between calculated and measured loss fraction is believed to be resolved by a more accurate intensity model. The larger experimental vertical error bars compared to the 511 nm results are due to less statistics collected in this measurement.
The resonant I 2 0 scaling at 1022 nm can be explained using multiphoton transition pathways and energy levels of 87 Rb depicted in Fig. 5. The unperturbed energy levels are experimental NIST values [33]. This figure indicates that the apparent high ionization probability can be induced by the 4 2 D and 5 2 F states. The 4 2 D resonance is energetically separated from the 5 2 S ground state by 2.40 eV which is close to the 2.43 eV energy difference corresponding to two 1022 nm photons but just outside the spectral bandwidth of the laser pulse. The 5 2 F resonance, however, is separated from the ground state by 3.63 eV so that it is accessible using three photons which add up to an energy of 3.64 eV with 0.02 eV bandwidth (see Methods for details). By exposing the atom to the femtosecond laser beam, the atomic energy levels are modified by a timedependent AC Stark shift during the laser pulse so that bound states can shift into or out of resonance. In order to answer the question whether the resonant process is a 2+1+1 photon process involving both, the 4 2 D and 5 2 F resonances or rather a 3 + 1 photon process via the 5 2 F state we present further results from the TDSE calcula- tions regarding the population of bound states after the laser pulse in Fig. 6a. By projecting the electron wave packet to the 5 2 S ground state and to the 4 2 D or 5 2 F excited state, respectively, we can trace the final population of these states with respect to the intensity. Additionally, we show the sum over all bound states of the atom. This survival fraction is used to derive the ionization probability by subtracting it from unity. It is apparent that the depopulation of the ground state after the pulse results in a final population of the 5 2 F state and ionization whereas the contribution of the 4 2 D state is negligible. The oscillations of the population of the ground and 5 2 F states presumably are due to Rabi-type oscillations in a complex system of three bound states and the continuum.
Because bound states are clearly populated and we experimentally detect remaining, non-ionized atoms on the 87 Rb D 2 transition, one could argue that the experiment only detects non-ionized atoms if they are in the 5 2 S ground state. However, if solely the 5 2 F excited state is populated in addition to the ground state, a single photon of the 780 nm detection pulse can also be absorbed on the 5 2 F to continuum transition which counts as a non-ionized atom in the absorption imaging. This is why the ionization probability here is derived from the sum over all bound states in Fig. 6a. Figure 6b shows the population of bound states during a laser pulse of I 0 = 7 × 10 10 W/cm 2 peak intensity. Al-though the 4 2 D resonance is negligibly populated after the pulse, a transient population during the pulse can be observed which enhances the ionization probability. The plot shows an oscillatory population transfer with a frequency of approximately 7 THz which modulates much faster oscillations not resolved in the figure. A Fourier analysis of the fast oscillations shows that the 4 2 D and 5 2 F populations follow a signal corresponding to mainly twice the laser frequency. This is an indication for a twophoton population transfer, thus explaining the I 2 0 scaling. Similar theoretical studies have also reported such oscillations of final states with respect to the intensity and time for alkali atoms [30,31,34].

CONCLUSION AND OUTLOOK
In summary, by cooling the atomic system considerably below MOT temperatures to approximately 100 nK in an optical far-off resonance trap we image local losses imprinted into the measurable atomic density profile. We have experimentally determined absolute ionization yields in ultracold atomic ensembles and have extracted absolute ionization probabilities for non-resonant and resonant ionization processes at two different wavelengths. We provide accurate data sets for precise tests and compare our experimental observation, without any arbitrary scaling, to ionization probabilities obtained by numerically solving the TDSE.
Our methodology allows for a pixel-wise evaluation of the losses for truly overcoming focal-averaging and measuring ionization processes at one single intensity. This would provide a simultaneous measurement at several intensities in the focal intensity distribution and will be highly beneficial for sources where parameters such as pulse duration and pulse energy fluctuate. Moreover, this would give a more direct comparison to theory where single intensities are used for calculations. Additionally, by varying the wavelength of the detection pulse, sensitive tracing of intermediate state population after the laser pulse is possible. In principle, the aforementioned methods can also be applied in a time-resolved manner by employing pump-probe techniques.
In contrast to magneto-optical traps, atoms trapped optically by far-off resonant laser light have the advantage of a vanishing excited state fractions and the absence of electric or magnetic fields. This is beneficial for electron or ion spectroscopy after ionization [35] which can reveal further insight into strong-field physics in alkali atoms. Our experimental techniques can also be used for absolute calibration of detector efficiencies in the rapidly growing number of experiments using cold atoms with charged particle detection [11][12][13][36][37][38].
Recent publications discuss that ponderomotive shifts for alkali atoms substantially lower the ionization probability compared to direct ADK results due to the high polarizability of alkali atoms [39][40][41]. When using longer pulsed laser wavelengths where the adiabatic approxima-tion is genuinely valid, these effects have to be taken into account and can be tested with a high sensitivity with our method. Note, however, that ADK does not take into account intermediate resonances which become important for longer wavelengths. Moreover, the population of bound states after the laser pulse and the critical influence of an AC Stark shift can impose coherent control schemes for excitations in ultracold samples [42].
Particularly interesting experiments include the investigation of the strong-field ionization of ultracold atoms before and after the Bose-Einstein condensation (BEC) phase transition. Previous results have been inconclusive whether the ionization rate increases or decreases in this case [43,44].
Finally, the presented measurements reveal the strongfield response of alkali atoms by accessing a transitional regime where the Keldysh parameter and the ionization probability are close to unity. The single outer-shell electron and the low ionization potential combined with a high polarizability compared to commonly used rare gas atoms make ultracold alkali atoms ideal model systems for studying strong-field physics.

Ultracold rubidium ensembles
Rubidium atoms are thermally evaporated from a dispenser into a vacuum system where the 87 Rb isotope is captured in a 2D MOT and transferred into a 3D MOT by a near-resonant pushing laser beam. After the 12 s long MOT loading phase, the atoms are cooled in a 4 ms expansion during an optical molasses phase. Subsequently, the atoms are transferred into a hybrid trap [28]. This trap consists of a magnetic quadrupole field generated by a pair of coils in anti-Helmholz configuration and a far-detuned crossed optical dipole trap at 1064 nm wavelength. The attractive potential generated by the dipole trap laser beams prevents Majorana losses in the node of the magnetic field. The temperature of the trapped atoms is reduced by forced radio-frequency (rf) evaporation. Subsequently, the magnetic field is switched off and the atoms are evaporatively cooled by reducing the laser intensity of the crossed dipole trap. The atoms are being held in the trap for 0.5 s before the ensemble is exposed to a femtosecond laser pulse.
After exposure to the femtosecond laser pulse, the remaining atomic density is mapped onto a CCD camera with a spatial resolution of (3 ± 0.5) µm by resonant absorption imaging (Fig. 1). An advantage of our setup is the availability of considerably lower temperatures compared to experiments with MOT targets so that we can directly image the ionization vacancy imprinted into the atomic density in real space. During the exposure time of the 50 µs long imaging light pulse the motion of atoms at nanokelvin temperatures is negligible because a typical distance covered is on the order of a few hundred nanometers which is far below the imaging resolution. The detection is performed at the D 2 transition of 87 Rb (780 nm wavelength). In the data processing, an exposureand saturation correction are applied to account for laser intensity fluctuations and population of excited states during imaging. The optical densities are kept considerably below the maximum detectable optical density before saturation correction ODmax = 3.7 given by the dynamic range and the noise level of the 12-bit CCD camera.
The parameters of both ultracold ensembles for non-resonant and resonant ionization are summarized with the imaging parame-  Table I.
The laser emits linearly-polarized pulses of 300 +30 −23 fs FWHM duration of the Gaussian temporal intensity distribution at a repetition rate of 100 kHz. Isolated laser pulses can be extracted using a pulse-picker driven by a pockels cell in the regenerative amplifier of the laser system. The pulses are focused onto the ultracold atoms by a f = 200 mm achromatic lens where the axis of the beam is tilted by 13.7 • with respect to the camera axis. In case of the four-photon ionization the pulses are directly transported to the focusing lens.
To study ionization yields at 511 nm wavelength the pulses are converted by second-harmonic generation (SHG) in a Beta-Barium-Borate (BBO) crystal resulting in pulses with a central wavelength of 511.4 nm and 1.75 nm FWHM bandwidth. The pulse duration of the Gaussian temporal intensity distribution in this case is 215 +20 −15 fs FWHM. Both pulse durations are measured in an intensity autocorrelation.
In case of the non-resonant measurement at 511 nm the pulse is applied in-situ in the optical dipole trap while for the resonant measurement at 1022 nm the pulse is applied after 2 ms time-offlight (TOF) to avoid ionization losses induced by leakage of highrepetition rate pulses from the oscillator cavity through the pulse picker. Since these pulses are below the conversion threshold for the SHG, leakage is not an issue in the 511 nm measurements.
The focus for both wavelengths is characterized using a focal diagnostics camera (compare Fig. 1) with 2.2 × 2.2 µm 2 pixel size. The pulse energy Ep is recorded simultaneously for each measurement to accurately determine the intensity in the focal volume for calibrated peak intensities I 0 . Figure 7 depicts the focal intensity distribution together with a 2D Gaussian fit for both wavelengths. While the 511 nm focus is very well represented by a single Gaussian intensity distribution, the 1022 nm focus deviates from this ideal model by two additional weak features on the bottom left corner of the image. In order to calculate the loss fraction expected by theory, we reconstruct the intensity distribution using a two-dimensional fit.  Table II.
The 3D intensity distribution of the laser pulse in space and time in this general case is given by where the sum accounts for the superposition of up to three elliptical Gaussian profiles k ∈ {I, II, III} offset along the y axis by y k 0 . For each profile k and each coordinate i ∈ {x, y} the waist is given by w k i (z) = w k 0i 1 + (z/z k Ri ) 2 with the Rayleigh length z k Ri = πw k 0i 2 /λ. The peak intensity includes the peak power P 0 = Ep/( √ 2πτ ) with pulse energy Ep and rms pulse duration τ = τ FWHM / 2 2 ln (2) . The amplitude of the main peak is normalized to A I = 1 while the amplitudes for the second and third beam (A II and A III ) are set to zero in the 511 nm case and determined by the fit in the 1022 nm case. Integration over the spatiotemporal intensity profile returns the pulse energy Ep.
The 2D fit is performed using Eq. (4) in the focal plane at z = 0. A global rotation of the coordinate system around the z axis by an angle αz is included into the fit. The resulting fit parameters with confidence intervals are summarized in Table II and the reconstruction is plotted in the bottom right display of Fig. 7.

Time-dependent Schrödinger equation
Ionization of rubidium is investigated theoretically by numerically solving the Time Dependent Schrödinger Equation (TDSE) within the single active electron approximation. In brief, TDSE is treated in spherical coordinates within the split-propagation paradigm. Angular dynamics is described within the basis of 20 spherical σ harmonics. The propagation step over the non-uniform grid of the radial coordinate is performed with the Crank-Nicolson algorithm.
The l-dependent interaction of the electron with the residual positive rubidium ion is parametrized to reproduce energies of the 5s ground state and six lowest excited states of each s, p, d, and f series. The average over the corresponding NIST tabulated multiplet energies is taken as experimental energies of the rubidium states. The eigenenergies of the used potentials reproduce the experimental values with 10 −4 accuracy.
The region of the wave packet propagation is 5200 a.u. Initially, only the 5s state is populated and the intensity of the Gauss shaped excitation pulse of 215 fs FWHM duration for 511 nm and 300 fs for 1022 nm is very small at this moment. An absorbing potential at the boundary of the radial mesh is used for suppression of wave packet reflection from this boundary. The propagation time is substantially larger than the pulse duration and therefore a quite noticeable decrease of the wave packet population at r ≤ 400 a.u. is obtained.
The results of the propagation are the time-dependent populations of the involved excited states. As the full probability of the system survival after ionization, we took the full wave packet population inside the sphere r ≤ 400 a.u. The part of the wave packet beyond this sphere is very small and contains contributions from slowly receding electrons and high Rydberg states. An analysis of the results reveals the role of resonant states in the case of a laser pulse of 1022 nm wavelength. Applying a laser pulse of 511 nm wavelength does not excite resonances.

Calculation of the loss fraction
Using ultracold atoms, a local 2D density distribution of the atoms is accessible via the absorption imaging technique. In order to calculate the fraction of lost atoms and compare it to different ionization models, a double Gaussian function according to Eq. (2) is fitted to the row sum of the absorption image where A 0 determines an offset while Ac, xc and σc correspond to amplitude, central position and width of the atomic cloud and Av, xv and σv describe the same parameters for the vacancy, respectively. The fraction of lost atoms can be derived by Eq. (3). For an example the reader is referred to Fig. 2.
In order to calculate the expected loss fraction, the full 3D initial atomic density is modeled by a Gaussian cloud consisting of N 0 atoms distributed within σa axial rms width and σr radial rms width, all given in Table I.
From the 3D intensity distribution in Eq. (4) the ionization probability P (x, y, z) after the laser pulse is obtained for each grid unit. In case of the TDSE results the ionization probability presented in Fig. 3a and Fig. 4a is used. For the MPI, OBI and ADK models the ionization probability is calculated as demonstrated in the Supplemental Material. This ionization probability map is now superimposed with the initial 3D atomic density map in Eq. (6) and N i = P N 0 atoms at a specific grid unit are removed from the initial 3D density leading to the final density after ionization. Subsequently, this 3D density distribution ρ f (x , y , z ) is projected along the imaging axis in order to obtain a simulated 2D density similar to the absorption image. The small angle between the ionizing laser beam and the camera axis of 13.7 • has been taken into account. Therefore, the 2D coordinate system is rotated by this angle with respect to the 3D coordinate system. The 2D density distribution ρ 2D (x, y) is convoluted with the optical resolution of the imaging system and evaluated in the same manner as the absorption image. The row sum ρ 1D (x) is calculated and fitted to Eq. (2) which results in a simulated loss fraction f sim l extracted from the fit according to Eq. (3).
A simulated absorption image compared to the corresponding experimental image and the projected row sum is presented in Fig. 2. The result for the simulated loss fractions with respect to the intensity is given in Fig. 3b and Fig. 4b. All calculations were performed on a three dimensional grid (x × y × z) of 500 × 120 × 120 grid units each 1 × 1 × 1 µm 3 in size so that the total simulated area spans 500 × 120 × 120 µm 3 . able online at: https://doi.org/10.1038/s42005-018-0032-5

AUTHOR CONTRIBUTION
The experiment was coordinated by J.S., M.D., K.S. and P.W. Construction and set up of the experiment was done by P.W., B.R., T.K. and J.S.. Measurements and data analysis were performed by P.W. with support from B.R., T.K. and J.S. The theoretical calculations were performed by A.K.K. and N.M.K. The manuscript was written by P.W. with discussions involving all co-authors.

Competing financial interests
The authors declare no competing financial interests.
Supplemental Material: Absolute strong-field ionization probabilities of ultracold rubidium atoms UNITS All equations are given in SI units where ε 0 = 1/(µ 0 c 2 ) represents the electric constant defined by the speed of light in vacuum c and the magnetic constant µ 0 . For conversion between atomic units (a.u.) and SI units these factors can be used: where E h is the Hartree energy, e the elementary charge and a 0 the Bohr radius. A conversion between the peak electric field E SI 0 and peak intensity I SI 0 in SI and atomic units (F a.u. (S3) A conversion to the peak electric field strength in atomic units F 0 from a given peak intensity I 0 in SI units can be achieved by using the relation

MODEL IONIZATION PROBABILITIES
In order to quantify an ionization process, the ionization rate w(t) is calculated as a measure of how many atoms per time-interval are ionized. For photoionization this rate depends on the instantaneous intensity I of the light field which is time-dependent for a laser pulse. It is described by an oscillating field term and an intensity envelope I(t) ⇒ w(I) = w(I(t)) ≡ w(t). When combining two ionization models, the total ionization rate w(t) = max {w 1 (t), w 2 (t)} (S5) is determined by using the dominating ionization rate of the two models. Once the ionization rate w(t) is known from the models presented below, the ionization probability P (t), i.e. the fraction of ionized atoms with respect to the total initial number of atoms N i /N 0 , is determined by the following differential equation: Solving the equation leads to the solution function that saturates at zero and unity for low and high ionization rates, respectively. Equation (S7) refers to the cumulated ionization probability at a specific time t in the intensity envelope. In order to calculate the ionization probability after a single laser pulse P = lim t→∞ P (t) the limit for t → ∞ has to be calculated. The ionization probabilities for rubidium using the models discussed below are presented in Fig. 3a and Fig. 4a in the main manuscript.

Multiphoton ionization
In a multiphoton process the photon energy E ph =hω is lower than the binding energy E I of the atoms so that the atom has to absorb multiple photons in order to release the bound electron into the continuum. The number of absorbed photons m = E I /E ph refers to the order of the m-photon process and the ionization rate The cross section depends on the wavelength and for a twophoton process in rubidium at 511 nm wavelength it has been measured and calculated to σ 2 = 1.5 × 10 −49 cm 4 s [29]. At 1022 nm only a calculated five-photon cross section σ 5 = 5×10 −139 cm 10 s 4 has been reported in the literature [45] that describes an ATI process. For reference, a calculated four-photon cross section for a slightly different wavelength at 1064 nm (Nd:YAG laser) σ N 4 = 1.32 × 10 −107 cm 8 s 3 [46] is also available in the literature.

Tunnel ionization
The tunnel rate of a bound electron inside an atom penetrating the barrier formed by the superposition of Coulomb potential and light field can be calculated using the Ammosov-Delone-Krainov (ADK) model [17,18,39,47] and the time-averaged rate for a linearly polarized laser pulse is obtained by w ADK (I) = |C n * 0 | 2 6 πR(I) E Ī h ×R 2n * −1 (I) exp − R(I) 3 (S10) where R(I) = 4E 3 2 Ī he ε 0 cme I (S11) and |C n * 0 | 2 = 2 2n * n * Γ(n * + 1)Γ(n * ) (S12) using the gamma function Γ and the effective principal quantum number n * = Z E I,H E I (S13) given by the charge state Z of the ion and the ionization potential E I related to the hydrogen ionization potential E I,H = 13.6 eV. Note that this is already the time-averaged result for an oscillating electric field of a laser pulse so that the total ionization rate can be obtained by integrating over the intensity-envelope of the laser pulse. If an integration over the instantaneous intensity I(t) = (cε 0 /2) |E(t)| 2 of the field oscillations is performed, the prefactor 6/[πR(I)] has to be omitted in Eq. (S10).