Absolute strong-field ionization probabilities of ultracold rubidium atoms

Understanding strong-field ionization requires a quantitative comparison between experimental data and theoretical models which is notoriously difficult to achieve. Optically trapped ultracold atoms allow to extract absolute nonlinear ionization probabilities by imaging the atomic density after exposure to the field of an ultrashort laser pulse. We report on such precise measurements for rubidium in the intensity range of 1 × 1011 – 4 × 1013 W cm−2. The experimental data are in perfect agreement with ab-initio theory, based on solving the time-dependent Schrödinger equation without any free parameters. We investigate the strong-field response of 87Rb 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. The authors report on absolute strong-field ionization probabilities of rubidium atoms exposed to femtosecond laser pulses and show that Ab-initio theory is in perfect agreement with the data at Keldysh parameters near unity.

U ltracold 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 well-defined 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 timescales 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 high-brilliance 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 creating electrons and ions on instantaneous timescales 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 γ ¼ ffiffiffiffiffi ffi 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. For small laser intensities, photoionization is considered as a perturbative process described by the absorption of multiple photons γ ) 1 ð Þ. At larger intensities, the laser field substantially modifies the Coulomb potential so that the electron can tunnel into the continuum γ ( 1 ð Þin an adiabatic picture [17][18][19] . Beyond a critical intensity, the barrier formed by Coulomb potential and light field is suppressed completely so that the electron can classically escape, which is referred to as over-the-barrier ionization 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 due to the small ionization potential E I . 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 and especially the measurement of absolute ionization probabilities is notoriously difficult since the ionization of atoms in strong laser fields depends on a variety of experimental settings (target density, excited state fraction, incoming photon-flux, detector efficiencies, and relative spatial alignment of target and laser pulse). These parameters are not always accessible or under control although the measured quantity depends nonlinearly on some of these factors. Often, the intensity calibration is even performed using the theory under test, for example in advanced electron spectroscopy experiments 23,24 . 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.
Ionization cross sections can be accurately measured using magneto-optically trapped (MOT) atoms by observing the trap loss rates with and without an ionizing laser beam 25,26 . This technique is very sensitive and gives access to low ionization rates. However, the excited state fraction needs to be calibrated because a MOT inherently relies on repeated absorption and emission cycles to perform cooling and trapping forces. In order to measure ionization losses from the ground state, the ionizing and trapping lasers have to be switched in an alternating fashion to first allow a decay back into the ground state and, after ionization, a fluorescence measurement with near-resonant laser light [27][28][29] .
We present an alternative methodology to overcome the above mentioned obstacles by employing optically trapped, ultracold 87 Rb ensembles at nanokelvin temperatures. We image the lineof-sight integrated two-dimensional (2D) density of the remaining atoms after ionization using absorption imaging. Here, regions of locally reduced atomic densities relate to 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.

Results
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 30 which consists of a magnetic quadrupole field and a crossed optical dipole trap at 1064 nm wavelength. This realizes a welldefined 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 as described in the Methods section. This experiment is performed at two different wavelengths corresponding to non-resonant and resonant ionization regimes. While a resonant ionization process gives rise to the population of intermediate bound state, the non-resonant case involves purely virtual states to reach the continuum. The latter case is realized at a wavelength of 511 nm with pulses of 215 þ20 À15 fs full width at half maximum (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 ≈ 100 μK, here the motion of atoms is negligible during imaging so that we have access to local atomic losses within the focal intensity distribution. 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 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 relating the number of ionized atoms N i to the total initial number of atoms N 0 is derived from the fit parameters as shown in Eq. (3). In order to compare the measurements to theoretically obtained ionization probabilities, the predicted loss fraction 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. Calculated ionization probabilites in Fig. 3a are presented along with the measured fraction of lost atoms f l according to Eq. (3) with respect to the peak intensity I 0 in Fig. 3b for 511 nm wavelength. The corresponding number of generated electrons or ions N i after the pulse 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. This quantity in connection with the excess energy of the electrons and the focal volume is important when considering the creation of ultracold electron bunches, ultracold plasma, or when investigating space-charge dynamics. In this work, collisions of charged particles with neutral atoms are assumed to be negligible due to the diluteness of the atomic sample and the low recollision probabilities of electrons accelerated during the laser pulse.
The increasing measured loss fraction (blue dots) originates from higher ionization probabilities P at higher peak intensities presented in Fig. 3a. The ionization probability can not exceed unity which is reached at 5 × 10 12 W cm −2 . In the spatial distribution of intensities in the focus, losses still increase beyond this intensity due to higher ionization probabilities in the wings of the intensity profile.
The calculated loss fraction f sim l expected by strong-field 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.
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 data sets 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 ionization (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 31,32 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 33 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 for the presented data points is never below 1.5, we can also test the tunnel-ionization 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 Supplementary Note 1. 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 probability, 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. This is especially evident from the data presented in linear scaling in Fig. 3c.
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 31,32 .
Resonant multiphoton ionization at 1022 nm. Figure 4 covers the resonant ionization at smaller photon energies using pulses of 1022 nm wavelength. Compared to the measurement at 511 nm in Fig. 3, significant losses become evident at lower intensities. The calculated ionization probabilities obtained by TDSE calculations are presented in Fig. 4a. In the absence of resonances, a multiphoton ionization process at this wavelength involves four photons and would scale with I 4 0 . For comparison, an I 2 0 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 non-resonant measurement, the focal intensity distribution here was not perfectly Gaussian. 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 scaling of a two-photon process. The additional assumption in the calculation that the electron is free once I OBI is exceeded leads to an overestimation of the losses (purple lines) just as including tunnel-ionization losses using pure Ammosov-Delone-Krainov (ADK) theory (yellow lines). c Same data in linear instead of log/log scaling. Vertical error bars represent the standard deviation of the individual loss measurements. Horizontal error bars consist of the standard deviation of the intensity in the respective bin and include systematic errors in the intensity measurement TDSE ionization probabilities, we have accounted for this by approximating the measured intensity distribution by a superposition of three Gaussian profiles. The measured foci for both wavelengths obtained by the focal diagnostics (see also Fig. 1) and the reconstruction are discussed in detail in the Methods section. 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 34 . 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 time-dependent 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 calculations 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 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 nonionized 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. Although 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 two-photon 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 31,32,35 .

Discussion
In summary, 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. By cooling the atomic system considerably below MOT temperatures to ≈ 100 nK in an optical far-off resonance trap, we image local losses imprinted into the measurable atomic density profile. 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 36 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][37][38][39] .
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 [40][41][42] . When using longer pulsed laser wavelengths where the adiabatic approximation 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 43 .
Particularly interesting experiments include the investigation of the strong-field ionization of ultracold atoms before and after Bose-Einstein condensation (BEC). Previous results have been inconclusive whether the ionization rate increases or decreases in this case 44,45 .
Finally, the presented measurements reveal the strong-field 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.

Methods
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  300 fs full width at half maximum (FWHM) duration is shown with respect to the peak intensity I 0 (red solid line). Population of the 5 2 S ground state (blue) and 5 2 F resonant state (yellow) shows that after the laser pulse the population is mainly transferred from the ground state into the continuum and into the 5 2 F state. The final 5 2 F population oscillates with the peak intensity. b Time-dependent population of the 5 2 S ground state (blue) as well as the 5 2 F (yellow) and 4 2 D (purple) excited states during the laser pulse of I 0 = 7 × 10 10 W cm −2 peak intensity (gray dashed line) 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 30 . 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 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 charge-coupled device (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 exposure-and 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 OD max = 3.7 given by the dynamic range and the noise level of the 12bit CCD camera.
The parameters of both ultracold ensembles for non-resonant and resonant ionization are summarized with the imaging parameters in Table 1.
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°w ith 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 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-of-flight (TOF) to avoid ionization losses induced by leakage of high-repetition rate pulses from the oscillator cavity through the pulsepicker. Since these pulses are below the conversion threshold for the SHG, leakage is not an issue in the 511 nm measurements.
The pulse energy E p is recorded by a silicon photodiode for each measurement to quantify the spatiotemporally integrated intensity. In order to calibrate the photodiode signal, the peak voltage U 0 is averaged over at least 20 pulses released by the pulse-picker at a repetition rate of 0.33 Hz for a constant laser amplifier output power. The low repetition rate ensures exceeding the recovery time of the carriers in the photodiode between two pulses so that saturation is avoided. The average power P at a high repetition rate f is measured with a thermal powermeter so that the pulse energy E p can be derived by E p = P/f and related to the peak voltage U 0 . This is done for the whole pulse energy range so that an interpolated calibration curve relates the measured peak voltage U 0 to the experimentally applied pulse energy E p . The error bars in Figs. 3 and 4 include systematic as well as statistical errors of the calibration procedure.
The focus for both wavelengths is characterized using a focal diagnostics camera (compare Fig. 1) with 2.2 × 2.2 μm 2 pixel size to obtain calibrated peak intensities I 0 using the measured pulse energy E p and pulse duration τ. 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 2D fit.
The 3D intensity distribution of the laser pulse in space and time in this general case is given by Iðx; y; z; tÞ ¼ I 0 exp À t 2 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 with the Rayleigh length z k Ri = πw k 0i 2 =λ. The peak intensity includes the peak power P 0 ¼ E p = ffiffiffiffiffi 2π p τ À Á with pulse energy E p and rms pulse duration τ ¼ τ FWHM = 2 ffiffiffiffiffiffiffiffiffiffiffiffiffi 2lnð2Þ p À Á . 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 E p .
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 2 and the reconstruction is plotted in the bottom displays of Fig. 7.
The difference of the two intensity distributions originates from the intermediate second-harmonic generation. First, the beam path changes and additional optical elements are used for the frequency doubling into 511 nm photons. Secondly, due to the inherent nonlinearity of the SHG process, low intensity regions are not as efficiently converted as regions of high intensity so that the beam profile is spatially filtered by the SHG. Third, both beams have a different wavelength, divergence and beam diameter so that the transport over the same optical elements to the ultracold atoms changes and spherical aberrations increase for the larger 1022 nm beam.   Table 2 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. (for conversion between atomic units and SI units see Supplementary Note 2). 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.
In the present paper only direct multiphoton ionization is addressed. There is a number of secondary ionization pathways which are not considered here. Freed electrons may either directly ionize secondary atoms or excite atoms which undergo subsequent photoionization. An analysis of the energy necessary for such processes shows that the corresponding primary electrons have to absorb more than five 1022 nm photons. Additionally, two excited atoms in the 5 2 F state can undergo a Penning-Auger-type process, when one of the atoms deexcites into the ground, 4 2 D or 5 2 P states and another atom is ionized. One can presume that those processes give a small contribution. This is strongly supported by the very good correspondence of the experimental data and the computational results. To analyze these processes, the spectra of ejected electrons have to be both measured and calculated. This is a challenging problem for further investigation.
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 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 fraction of lost atoms is derived from the fit parameters according to 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 ρ 0 ðx; y; zÞ ¼ N 0 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 1.
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 Figs. 3a and 4a is used. For the MPI, OBI, and ADK models the ionization probability is calculated as demonstrated in Supplementary Note 1. This ionization probability map is now superimposed with the initial 3D atomic density map in Eq. (6) and N i = PN 0 atoms at a specific grid unit are removed from the initial 3D density leading to the final density ρ f ðx; y; zÞ ¼ ρ 0 ðx; y; zÞ 1 À Pðx; y; zÞ ½ ð7Þ 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 Figs. 3b and 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 .
Data availability. The data that support the findings of this study are available from the corresponding author upon reasonable request.
Code availability. Any custom computer code used in this study is available from the corresponding author upon reasonable request.  Two-dimensional fit parameters with confidence intervals of the measured focal intensity distributions for 511 nm and 1022 nm wavelength as presented in Fig. 7