Joint probability calculation of the lateral velocity distribution in strong field ionization process

We describe an approach to the description of the time-development of the process of strong field ionization of atoms based on the calculation of the joint probability of occurrence of two events, event B being finding atom in the ionized state after the end of the laser pulse, event A being finding a particular value of a given physical observable at a moment of time inside the laser pulse duration. As an example of such an physical observable we consider lateral velocity component of the electron’s velocity. Our approach allows us to study time-evolution of the lateral velocity distribution for the ionized electron during the interval of the laser pulse duration. We present results of such a study for the cases of target atomic systems with short range Yukawa and Coulomb interactions.

We describe an approach to the description of the time-development of the process of strong field ionization of atoms based on the calculation of the joint probability of occurrence of two events, event B being finding atom in the ionized state after the end of the laser pulse, event A being finding a particular value of a given physical observable at a moment of time inside the laser pulse duration. As an example of such an physical observable we consider lateral velocity component of the electron's velocity. Our approach allows us to study time-evolution of the lateral velocity distribution for the ionized electron during the interval of the laser pulse duration. We present results of such a study for the cases of target atomic systems with short range Yukawa and Coulomb interactions.
The conventional framework for understanding the basic features of the tunneling ionization process, known as the Strong Field Approximation (SFA), was laid out by Keldysh 1 and developed and refined subsequently [2][3][4][5][6][7][8][9] . In this framework the tunneling regime of ionization is defined as the regime characterized by the small values of the so-called Keldysh parameter γ = ω √ 2|ε 0 |/E (here ω , E and |ε 0 | are the frequency, field strength and ionization potential of the target system expressed in atomic units).
A number of nonlinear phenomena, such as above-threshold ionization (ATI), high harmonic generation (HHG) and non-sequential double ionization (NSDI) may occur in this regime. The most essential features of these processes can often be understood with the help of the well-known simple man model (SMM) 6,10-13 , which represents electron's motion after the ionization event occurs as entirely classical. This fact underlies a variety of highly efficient semi-classical approaches [12][13][14][15][16][17][18][19] , in which the electron ionization event is treated quantum-mechanically, with electron's ionization occurring near the local peaks of the electric field, and the subsequent electron motion treated classically or semi-classically. These approaches include the well-known TIPIS model 12,16,20 , the quantum trajectory Monte Carlo model (QTMC) 19 , semi-classical two-step model 21 or Coulomb quantum orbit strong-field approximation (CQSFA) 22,23 .
Despite these semi-classical methods giving often impressive agreement with more rigorous calculations based on the solution of the time-dependent Schrödinger equation (TDSE), it is not quite clear what the very notion of the 'ionization event' means from the point of view of the conventional quantum mechanics (QM). Indeed, this notion is not easy to define even in classical physics. This problem manifests itself in the treatment of ionization in the standard QM framework in the following way. QM provides us with only a wave-function to describe evolution of the system in the laser field, and for the moments of time inside the interval of the laser pulse duration, when the wave-packet describing ionized electron has not left the atom yet, it is difficult to unambiguously single out the part of the wave-function describing ionized electron from the total wave-function of the system. At first glance such a separation of the total TDSE wave-function into the bound component and the ionized wave-packet arises quite naturally in the SFA or the Perelomov-Popov-Terent' ev (PPT) approaches [1][2][3][4]8,24 . We should note, however, that this splitting of the total wave-function in two components in these theories, although extremely useful and physically appealing, is not quite rigorous. That can been seen, for instance, by noting that this splitting is not gauge-invariant. It is this fact that is ultimately responsible for the lack of the gauge-invariance of the SFA or PPT approaches 9 . To see more clearly the origin of the problem we will recapitulate briefly, following the work 25 , the theoretical foundations of both approaches, which can done most easily by using the notion of the unitary time-evolution operator Û (t, τ ) which maps state of the system at time τ to the state at time t. For a system with the time-dependent Hamiltonian operator Ĥ (t) =Ĥ A (t) +B(t) . the integral Dyson equation for the evolution operator Û (t, τ ) can be written in two equivalent forms as: where Û 0 (t, τ ) = exp −iĤ atom (t − τ ) is the evolution operator describing evolution driven by the time independent field-free Hamiltonian Ĥ atom and Û F (t, τ ) is the so-called Volkov evolution operator with a simple known analytical form 9 which describes quantum evolution driven by the Volkov Hamiltonian Ĥ F (t).
Equations (3) are exact but are not very useful, being just as difficult to solve as the original TDSE with which we started. The widely used and important SFA and PPT approximations are obtained if we substitute Û F (t, τ ) for Û (t, τ ) under the integral in the first Dyson equation (1) and Û 0 (τ , 0) for Û (τ , 0) in the second equation (1): If we use the expression for the evolution operator in the first equation (4) to find the wave-function at the moment of time t = T 1 at the end of the laser pulse for a system which was initially in the field free state φ 0 , we obtain an approximate wave-function: � SFA (T 1 ) =Û SFA (T 1 , 0)φ 0 . Projecting this wave-function on a planewave state |k� and dropping the term �k|Û 0 (T 1 , 0)|φ 0 � , one obtains the well-known expression for the ionization amplitude used by Keldysh 1 , provided we use the length form for the interaction Hamiltonian Ĥ int (t) . If, instead, we use the velocity gauge to describe field-atom interaction, we obtain the expression for the ionization amplitude used in the well-known Strong Field Approximation (SFA) theory 2,3 . Using, on the other hand, the approximate evolution operator in the second equation (4) to evaluate the wave-function at the moment of time t = T 1 at the end of the laser pulse for a system which was initially in the field free state φ 0 , we obtain an approximate wavefunction � PPT (T 1 ) =Û PPT (T 1 , 0)φ 0 . Projecting this wave-function on a plane-wave state |k� and dropping the term �k|Û F (T 1 , 0)|φ 0 � which does not contribute to the probability current, one obtains the expression for the ionization amplitude used in the Perelomov-Popov-Terentiev (PPT) theory 4 .
We see, thus, that both in the PPT and SFA approaches we obtain expressions for the ionization amplitudes by omitting certain terms from the total wave-function. These terms, however, are different in the two approaches and, moreover, it is this omitting that breaks the gauge invariance which is of course present in the complete theory. Under a gauge transformation only the total wave-function of a system is transformed in a way ensuring gauge invariance of the final results for the observables. The parts of the wave-function obtained by splitting it into different components will not, in general, possess this property. One needs some, therefore, additional theoretical ingredients in the theory, which may allow to define the notion of the ionized part of the wave-function. An example of such a theoretical development is provided by the well-known back-propagation technique [26][27][28][29][30] . This procedure allows to single out only ionized electron and to follow its prehistory back in time by constructing the ionized wave-packet from the part of the coordinate wave-function localized far from the nucleus for times long after the end of the driving pulse.
In the present work we describe another procedure allowing to achieve such a separation of the total wavefunction describing the evolution of an atom in the laser field into the ionized and non-ionized components. We will use the notion of the conditional and joint probabilities for this purpose. In the Section Theory below  www.nature.com/scientificreports/ we describe our theoretical procedure. We will illustrate this procedure by applying it to a study of the evolution of the lateral velocity distribution (i.e. distribution of the velocity components in the directions orthogonal to the polarization vector of the driving field) of ionized electrons for times inside the laser pulse duration. These applications are described in the Section Results below. We conclude in the Section Conclusion making a brief summary of the results and future perspectives. Atomic units (a.u.) are used throughout the paper.

Theory
Joint probability in Quantum Mechanics. We begin by presenting a short summary of the theory of the conditional and joint probabilities in QM following work 31 . Suppose we have a system described by the state vector | � , and two observables with the corresponding Hermitian operators Â and B . The probability to find the system upon measurement in an eigenstate of B belonging to a spectral interval B of the operator B , can be computed, as is well known, as follows: where Q (� B ) is the Hermitian projection operator which can be expressed in terms of the eigenstates | � of the operator B as: . Analogous formula, where we will have to use a projection operator Q (� A ) , which can be similarly expressed in terms of the eigenstates of the operator Â , can be written, of course, for the probability P(� A ) to find the system in an eigenstate of Â belonging to a spectral interval A of the operator Â . Suppose now that we are interested in the joint probability P(� A &� B ) of finding observable A in the spectral interval A , and observable B in the spectral interval B . One might try to define such a joint probability by an expression: This definition would make perfect sense if the quantum projection operators Q (� A ) and Q (� B ) were classical quantities or at least quantum commuting operators. Unfortunately, this is not necessarily always the case. If the operators Â and B do not commute, the corresponding projection operators Q (� A ) and Q (� B ) do not commute either, the operator product Q (� A )Q(� B ) in Eq. (6) becomes non-Hermitian and the joint probability defined by Eq. (6) is generally complex-valued. This is, of course, another manifestation of the well-known difficulty that one encounters when trying to assign a meaning to the joint probability distributions of the observables described by non-commuting operators using the quasi-probability distributions, such as Wigner and Husimi distributions 32-35 which may not necessarily be strictly positive. With Eq. (6) we have the same problem in another disguise, it can give, as we see, complex-valued joint probability distributions. It can be argued 31,36,37 that such quasi-probability distributions, which are not strictly positive, can be incorporated in the framework of the QM, and that they appear naturally in the description of the classically forbidden processes such as the tunneling process 36 . Negative probabilities are, however, difficult to interpret, and many physicists are somewhat reluctant in accepting such a notion.
It is the tunneling process due to atomic ionization in the external electric field that interests us in the present work. We will see below that in some circumstances Eq. (6) gives us joint probability distributions which are perfectly legitimate in an ordinary probabilistic sense, and which may provide useful information about the development of the tunneling ionization process. Joint probabilities for the strong field ionization process. We will be studying below ionization of a single-electron atom in a strong electromagnetic field. The evolution of the system in time is described by the three-dimensional TDSE: with H atom =p 2 2 + V (r) -field free atomic Hamiltonian, and Ĥ int (t) -interaction Hamiltonian describing atomfield interaction. We use the length gauge for the latter operator: where E(t) is electric field of the pulse. We assume that pulse is linearly polarized along z-axis. The pulse is defined in terms of the vector potential: where T 1 is the total pulse duration, which we choose to be one optical cycle (o.c.): T 1 = 2π/ω , for the majority of the calculations reported below. We use a single cycle pulse for purely computational reasons as the calculations become rather time-consuming for longer pulses, our approach can be generally applied for long pulses as well. We will present below results for different values of the pulse base frequency ω , peak field strength E 0 and different atomic potentials V(r). www.nature.com/scientificreports/ To apply Eq. (6) and the notion of the joint probability distribution in QM to the process of the tunneling ionization we should first specify the operators Â and B in this equation. We note first that it is by no means necessary that both projective measurements in Eq. (6) are performed at the same moment of time. We may well assume that one of the projective measurements (let's say the measurement described by the projection operator Q (� B ) ) is performed at the time t 2 , while the measurement described by the operator Q (� A ) is performed at the time t 1 , with t 2 > t 1 . To modify Eq. (6) accordingly it is convenient to use the Heisenberg picture of the QM, in which the state vector � = φ 0 is independent of time (here φ 0 is the initial state of the system), while the projection operators Q (� A ) and Q (� B ) evolve in time. With the use of the Heisenberg picture the generalization of the Eq. (6) for the case of different times is straightforward: with where Û (t, 0) is the evolution operator, driving quantum evolution of the system, and |φ 0 � is the initial atomic state (which we assume to be the ground state of the field-free atomic Hamiltonian H atom . Definition Eq. (10) of the joint probability distribution is quite similar to the definition of the two-time correlation function [38][39][40][41] , describing correlations between two observables with corresponding quantum mechanical operators Q (� A ) and Q (� B ) at times t 1 and t 2 , respectively. Correlation functions can be employed 42,43 to study correlations between different observables in the strong field ionization process. Such a study may provide useful information about presence or absence of the correlations between observables at different moments of time.
On the other hand, the two-time correlation functions, being generally complex-valued objects, do not have direct physical interpretation. The approach we describe in the present work differs from the approach based on the two-time correlation functions analysis in one important aspect. We assume that the quantum operators Q H (� A , t 1 ) and Q H (� B , t 2 ) in Eq. (10) are are not arbitrary operators, but quantum mechanical projection operators. This implies that both is again a positive definite Hermitian operator with positive real expectation value. It is easy to see, in fact, that P is then a projection operator, satisfying P 2 =P . If Q H (� A , t 1 ) and Q H (� B , t 2 ) commute we can, therefore, assign direct physical meaning of the joint probability to the quantity defined by the Eq. (10). We will now specify the observables A and B and the projection operators Q (� A ) and Q (� B ) in Eq. (11).
Choice of the observable B. We will assume from now on that in Eq. (10) the later moment of time t 2 = T 1 , the moment of time when the laser pulse (9) is switched off. For the projection operator Q (� B ) we will use the projection operator projecting the state vector on the continuous part of the energy spectrum of the atomic Hamiltonian H atom ,and by the spectral interval B we will understand the whole range of positive electron energies of the continuous spectrum of the field-free atomic Hamiltonian. In other words, Q(� B ) =Q c , where Q c is the projection operator on the continuous spectrum of the field-free atomic Hamiltonian, so that for any state vector | �: where the sum on the r.h.s includes all bound states of the field-free atomic Hamiltonian H atom . It is worth discussing briefly what we have achieved by choosing this particular form of the projection operator Q (� B ) . Observable B with this choice of the projection operator provides basically the answer to the question: will the electron be ionized or not, it might be pictured as an observable having the value one if electron is ionized, and zero otherwise. The joint probability distribution (10) (provided it gives a legitimate probability distribution, we will touch on this question later) gives us therefore a joint distribution of the observable A at the time t 1 inside the laser pulse duration upon the condition that electron will be ionized in the future. By using the joint probability approach and our choice of the observable B we are able, therefore, to separate ionized and non-ionized electrons. As we noted in the Introduction such a separation is by no means trivial and requires some additional theoretical notions, such as the notions based on different ionization criteria employed in the back-propagation technique [26][27][28][29][30] . We can look at the Eq. (10) as yet another way to achieve this separation and study prehistory of the ionized electron based on the notion of joint probability distribution and the choice of the observable B we made above.
Choice of the observable A. As far as the observable A in Eq. (10) is concerned, we can choose, in principle, any observable which will provide useful information about the ionization process. Our only limitation is that the choice of A must be such that the joint distribution (10) be at least approximately legitimate probability distribution. What it means in practice is that we must choose A so that the imaginary part of the distribution (10) be small comparing to its real part for all values of A of interest, i.e. for the bulk of the distribution where the joint probability has non-negligible values. If this condition is satisfied, the projection operators in Eq. (10) will (approximately) commute, and Eq. (10) will automatically give a joint distribution which will be positive at least for the bulk of the distribution. We will, of course, still have non-physical negative probabilities on the edges of www.nature.com/scientificreports/ the distribution, where the joint probability is small. That is inevitable unless the projection operators Q (� A , t 1 ) and Q (� B , t 2 ) in Eq. (10) strictly commute, but that is presumably not very important.
That was the strategy we pursued in choosing the observable A. It turns out that, while it is not generally possible to choose a physically interesting observable A so that imaginary part of the joint distribution (10) be always small, it is possible to achieve this goal for at least some values of time t 1 inside the pulse. That can still provide us with valuable information about development of the ionization process and evolution of the observables characterizing ionized electrons inside the interval of the pulse duration .
An observable A which we will study below in detail is the lateral component of electron's velocity, i.e component of the velocity perpendicular to the polarization vector of the laser pulse ( z−direction for the geometry we use in Eq. (9)). It will be convenient in the following to use a cylindrical coordinate system (ρ, φ, z) in the vector space of electron's velocities v with the longitudinal axis along the z− direction. Let us consider a set of regions k of the electron's velocity space, where integer k = 0, 1 . . . , such that the components of the vector v in the cylindrical coordinate system we introduced above satisfy for any integer k: In other words, k is the space between two cylinders of infinite height and radii of kd and (k + 1)d . In the calculations below we use d = 0.1 a.u. Parameter d affects the overall 'resolution' in the velocity space which we may hope to achieve in the framework of our approach. We should choose is so as to not to loose any important fine details of the lateral velocity distributions. We will justify the choice of the value d = 0.1 a.u. that we made below. Let χ � k (v) be the characteristic function of the We define now the action of the projection operator Q k on a state vector | � as follows. We first perform Fourier transform of the coordinate wave-function �(r) obtaining momentum space wave-function � (v) . We define now Q k | � as an inverse Fourier transform of � (v)χ � k (v) . With this definition � |Q k | � is clearly the probability to detect electron's velocity in the region k of the velocity space we defined above. The regions k with k = 0, 1 . . . in the velocity space cover all the space and do not intersect, they satisfy, therefore, Q iQj = δ i jQ i . We now identify Q k with Q (� A ) where, for every integer k it is understood that spectral region A coincides with the region k in the velocity space. In physical terms, distribution (10) in which we substitute projection operator Q k for observable A gives us a probability of finding electron's velocity anywhere in the region k of the velocity space at the moment of time t 1 provided that electron is found to be ionized at the end of the laser pulse.
Calculation of the joint probability. Following the discussion presented in the two preceding subsections, we rewrite the formula (10) for the joint probability using more compact notation: where in the first line we used the Heisenberg representation for the projection operators (as in Eq. (10)), while to derive the second line we used the transformation equations (11) and the property |�(t)� =Û(t, 0)|φ 0 � of the evolution operator, so that |�(t)� is the solution of the TDSE obtained at the moment t from the initial (ground) atomic state |φ 0 � . It is the quantity on the r.h.s of the second line of Eq. (13) that was actually computed in our calculations. We did it as follows. The TDSE (7) is propagated forward in time on the interval of the laser pulse duration (0, T 1 ) starting with an initial atomic state |φ 0 � , thus obtaining the state vector |�(T 1 )� describing atomic system at the end of the laser pulse. Next, we apply the operator Q c to the vector |�(T 1 )� and propagate back in time both the resulting vector and the original vector |�(T 1 )� , obtaining two time-dependent state-vectors: |� 1 (t)� =Û(t, T 1 )Q c |�(T 1 )� and |�(t)� =Û(t, T 1 )|�(T 1 )� . For any given moment t inside the interval of the pulse duration we can now compute the joint probability defined by the Eq. (13) as a matrix element �Q i �(t)|� 1 (t)� (we used here Hermicity of the operator Q i ). This calculation requires multiple solutions of the TDSE for different regions k defining the projection operators Q k . The 3D TDSE was solved numerically using the procedure we tested and described in detail in [44][45][46] . The procedure relies on representing the coordinate wave-function as a series of spherical harmonics with quantization axis along the pulse polarization direction. Spherical harmonics with orders up to L max = 70 were used. The radial variable is treated by discretizing the TDSE on a grid with the step-size δr = 0.05 a.u. in a box of the size R max = 200 a.u. Necessary checks were performed to ensure that for these values of the parameters L max and R max convergence of the calculations has been achieved. The solution of the 3D TDSE was propagated both forward and backward in time using the matrix iteration method 47 . We did calculations for two atomic systems: system with the field-free dynamics governed by the short range (SR) Yukawa potential V (r) = − 1.9083e −r r and the hydrogen atom with Coulomb potential V (r) = − 1 r . Both these systems have the same ionization potential I p = 0.5 a.u. and we use their ground s-states as initial states |φ 0 � in the calculations below.
The joint probability distribution in Eq. (13) depends on the time moment t and the region k in the velocity space. We remind that this region was defined above as the volume in the velocity space between two cylinders of radii kd and (k + 1)d where k = 0, 1, . . . : � k = (kd ≤ v ρ < (k + 1)d; −∞ < v z < +∞) , and we used d = 0.1 a.u.
By introducing a variable v ⊥ = (k + 1/2)d the joint probability distribution in Eq. (13) can be considered as a (discretized) function of the lateral velocity v ⊥ and time t computed on a lateral velocity grid v ⊥ = (k + 1/2)d , k = 0, 1, . . . . We will use this fact to simplify the notation yet a bit more. We define now a function: www.nature.com/scientificreports/ which tells whether the imaginary part of the joint probability can be considered as small for a given moment of time t. The choice of the particular form of the function G(v ⊥ , t) is not unique, of course, and even is not very important. We only need a convenient measure allowing us to gauge the relative magnitudes of the real and imaginary parts of the joint probability defined in Eq. (13) and to judge to what extent the joint probability is a real quantity. The function defined in Eq. (14) is just the simplest choice allowing to achieve this purpose. As we discussed above, in general, we cannot expect that for a given t G(v ⊥ , t) vanishes for all v ⊥ altogether. We can, however, have G(v ⊥ , t) small for the bulk of the v ⊥ -distribution, i.e., in the region of lateral velocities where Eq. (13) has non-negligible values. In such situations we can use Eq. (13) to define a physically legitimate probability distribution for lateral velocities v ⊥ . Our strategy, therefore, is to look at the regions in the (t, v ⊥ )-plane where we have G(v ⊥ , t) ≪ 1 . If such regions in the (t, v ⊥ )-plane can be found for time t inside the laser pulse duration we will obtain a means of calculating a physically sensible probability distribution for lateral velocities v ⊥ . We can compare the distributions thus obtained with the theoretical predictions obtained in the framework of the SFA. The well-known SFA expression for the lateral velocity distribution reads 8,9 : where I is the ionization potential, E 0 -electric field strength, and v 2 ⊥ = v 2 x + v 2 y for the geometry we employ. We do not specify the constant A in this expression as we will be interested below in the distribution shape which is described by the exponential function. Expression (15) gives us the final velocity distribution at the detector when the laser pulse is gone, but it can also be used as a plausible expression for the lateral velocity distribution within interval of the pulse duration, as is done e.g., in the semi-classical simulations 12,16,20,21 . The rational behind that is that a linearly polarized electric field does not affect the distribution in the lateral direction during the electron's motion subsequent to the ionization event. The width of the lateral velocity SFA distribution (15) can, therefore, be used as a guide in choosing the value for the parameter d we used to foliate the velocity space. For the field strengths of the order of E 0 ≈ 0.07 we are interested in, the full width at half maximum (FWHM) of the distribution (15) is approximately 0.45 a.u., so our choice of d = 0.1 a.u. allows us to foliate the velocity space into the layers thin enough so as not to loose important details of the velocity distribution. The accuracy could be increased by using a smaller value of the parameter d at the expense of additional computing time, but as we shall see below, with the foliation parameter d = 0.1 a.u. the lateral velocity distribution obtained in the framework of our approach agrees pretty well with the SFA expression (15) for the short-range interaction, confirming the overall consistence and accuracy of the procedure.
To be able to compare the SFA distribution (15) to the joint distribution (13) we should take into account the fact that (13) refers to the joint probability of finding electron's velocity in the region k : � k = (kd ≤ v ρ < (k + 1)d; −∞ < v z < +∞) between the two cylindrical surfaces, while expression Eq. (15) refers to the Cartesian volume element in the velocity space. To take into account this difference of the volume elements we can observe that the volume k is proportional to v k,⊥ , where v k,⊥ = (k + 1/2)d . We must therefore divide the joint probability in Eq. (13) by v k,⊥ obtaining the distribution: which refers to the Cartesian volume element in the velocity space and can, therefore, be compared with the results predicted by the SFA formula (15). We will perform such a comparison below, by fitting the results we will obtain using Eq. (16) with analytical formulas. We will employ two types of fit, one for the SR Yukawa potential and another for Coulomb potential. The SR interactions fall into the domain of the validity of the standard SFA 8,9 . For the SR interaction, therefore, we will fit our results using the fitting expression based on the SFA formula (15): where A and β are used as fitting parameters. Parametrization in terms of the parameter β is convenient as for both Yukawa and Coulomb systems we have I = 0.5 a.u. Therefore, for the lateral velocity distributions having shapes similar to the ones predicted by the SFA formula (15), parameter β should be approximately equal to the peak field strength E 0 of the laser pulse.
For the Coulomb potential situation is a bit more complex. Due to the well-known effect of the Coulomb focusing 48 , lateral velocity distribution acquires a cusp at zero transverse velocity. Such a cusp-like structure can be described using an expression 16,44,49 .
where A, α and β are used as fitting parameters. Depending on the value of the non-integer parameter α the distribution function in (18) exhibits at v ⊥ = 0 a discontinuity in first or higher order derivatives, thus serving as a model of the cusp-like behavior.
The lateral velocity distributions that we compute using Eq. (16) pertain to the moments of time insider the interval of the pulse duration. We will compare them also to the asymptotic (obtained in the limit of large times) velocity distributions P as (v ⊥ ) , which we compute using the standard prescription, by projecting the Such a comparison is of interest, especially for the Coulomb potential, as it shows how cusp in the lateral velocity distribution develops in time. It is known that the Coulomb atomic system needs to evolve for some time for the cusp to develop 50 . Comparing asymptotic lateral velocity distribution (19) to the distributions (16) we will be able to have a glimpse at how this development of the cusp actually occurs.

Results and discussion
Yukawa potential. We begin by presenting the results we obtain for the Yukawa potential and field parameters ω = 0.03 a.u., E 0 = 0.07 , which places us relatively deep in the tunneling regime. In the framework of the strategy we outlined above, we will analyze first the function defined in Eq. (14) which will tell us at what particular (if any) times inside the laser pulse the notion of the joint probability can be sensibly used. This function is shown in Fig. 1a. The plot shows a rather complicated pattern, important for us are the "good" regions in the (t, v ⊥ )-plane, where G(t, v ⊥ ) is small. We need, of course, some quantitative criterion of "smallness". As such, we choose the condition G(t, v ⊥ ) ≤ 0.1 , which, as we will try to show, is sufficiently small threshold value for the notion of the joint probability distribution to make sense. One can see from Fig. 1 that there are areas (shown in black) for times inside the laser pulse duration, where this criterion is satisfied for the large enough intervals of lateral velocities. In these areas, as we argued above, the joint probability distribution is a sensible notion, and we can compute the lateral velocity distributions by taking cuts of the Cartesian distribution P(v ⊥ , t) we defined above in Eq. (16) along the lines t = const. The distribution P(v ⊥ , t) is shown in Fig. 1b, and its cuts P(v ⊥ , t i ) taken along the lines t = 0.548T and t = 0.7T are shown in Fig. 2. As an inspection of the Fig. 1 shows, for both these cuts the areas where G(t, v ⊥ ) is small extend sufficiently far in v ⊥ -direction so that Cartesian lateral velocity distributions (16) could be computed for the intervals of v ⊥ containing the bulk of the velocity distribution. The results shown in Fig. 2 for the two moments of time inside the laser pulse duration show very nice agreement with the Gaussian fit (17) based on the SFA expression (15). The prediction given by the SFA for the value of the parameter β in Eq. (17) is β = 0.07 a.u. for the field strength E 0 we use. The values we obtain for β by fitting the computed lateral velocity distributions are: β = 0.0729 a.u. for t = 0.548T and β = 0.0710 a.u. for t = 0.7T . In Fig. 3 we show evolution of the Cartesian lateral velocity distributions for the three values of time inside the interval of the laser pulse duration (indicated by the dotted white lines in Fig. 1) and compare these distributions to the asymptotic distribution P as (v ⊥ ) defined in Eq. (19). As one can see from Fig. 3a, all these distributions, both the ones computed for times inside the interval of the pulse duration and the asymptotic distribution (19) closely agree and shape of the distributions evolves little with time. This fact is illustrated also in Fig. 3b, where we show evolution of the parameter β in time. The pattern of the evolution we see in Fig. 3b agrees completely with the picture provided by the SFA, where electron emerges into the continuum at the moment of the maximum field strength, with the lateral velocity distribution given by Eq. (15), which does not change subsequently. This latter fact is due to the short range character of the potential which is assumed in the SFA, and the fact that in the absence of any long range forces the electric field of the pulse cannot affect the lateral components of electron's velocity.  www.nature.com/scientificreports/ Coulomb potential. Results that we obtain for the Coulomb potential are somewhat different. This is to be expected as we have in this case the strong long range Coulomb force. We use the same field parameters: ω = 0.03 a.u., E 0 = 0.07 as above. As in the case of the short range interaction, we proceed by calculating the function defined in Eq. (14) to find the values of t which could be used for calculating lateral velocity distributions for the times inside the laser pulse duration. This function is shown in Fig. 4a. For the Coulomb potential the region in the (t, v ⊥ )-plane where G(t, v ⊥ ) is small, proves larger than in the case of the Yukawa potential, allowing continuous scan of the velocity distribution for all times t after the midpoint of the pulse. The distribution P(v ⊥ , t) is shown in Fig. 4b, and the cuts P(v ⊥ , t k ) taken along the lines t = 0.548T and t = T are shown in Fig. 5. We also show results of the fitting procedure applied to the calculated distributions. The fitting procedure that we employ is based on the expression (18) which is non-analytical at v ⊥ = 0 , thereby taking into account presence of a cusp 48 for the Coulomb case 48 .
As one can see, the fits based on (18) reproduce fairly well the calculated lateral velocity distributions. One can also note that the distribution is more sharply peaked at the end of the pulse than at the moment t = 0.548T near the midpoint of the pulse. This evolution of the lateral velocity distribution is presented in more detail in Fig. 6 for different moments of time inside the laser pulse duration. One can observe from Fig. 6a that the non-analytical character of the velocity distribution and its departure from the Gaussian become progressively more pronounced as we approach the end of the laser pulse. This evolution of the lateral velocity distribution does not stop there, of course, long range Coulomb force keeps distorting the velocity distribution long after the pulse is gone. The asymptotic lateral velocity distribution obtained using Eq. (19) is shown in Fig. 6b. One can see that it is much sharper yet than the distribution we obtain for T = 1 . This evolution of the later velocity distributions we see in the Coulomb case agrees with the observation made in 50 , where it was noted that the cusp formation requires large (strictly speaking infinite) time. At any finite moment of time, the lateral velocity distribution remains an analytical function of v ⊥ which however, becomes, as time progresses, increasingly sharper in a vicinity of the point v ⊥ = 0 50 , approaching thus the cusp-like behavior present in the asymptotic (i.e. obtained for the infinite time) distribution shown in Fig. 6b. This process of the cusp formation is illustrated in Fig. 7, where we show evolution of the fitting parameters in Eq. (18) with time. One can see that parameter α in Eq. (18) whose role in  www.nature.com/scientificreports/ this expression is to mimic the cusp-like behavior, progressively decreases from the nearly Gaussian value α ≈ 2 , at times close to the instant of ionization, to the value α ≈ 1.35 at t = T . Such a low value of α is responsible for the rather sharp character of the lateral velocity distribution we observe for t = T in Fig. 6a. Indeed, the value α ≈ 1.35 at t = T implies that the first derivative of the lateral velocity distribution at v ⊥ = 0 is continuous and finite, while the second derivative is infinite. The asymptotic infinite time distribution shown in Fig. 6b is sharper yet, we obtain the value α = 0.55 if we apply the fit based on Eq. (18) to the distribution we calculate using Eq. (19). That means that for the asymptotic distribution it is the first derivative which is infinite at v ⊥ = 0 . These observations illustrate again the statement made in 50 to which we referred above, the cusp in the lateral velocity distribution becomes fully formed in the limit of infinite time.
Field dependence of the lateral velocity distributions. We present in this section the results we obtain for the lateral velocity distribution for different driving pulse parameters. The aim of performing these calculations was primarily to convince ourselves that the procedure we devised for the calculation of the lateral velocity distributions can be applied for a wide range of the laser field parameters. For the calculations presented below we use the frequency ω = 0.057 a.u., and we vary the peak electric field strength. We follow the same strategy we used above and begin with presenting the analysis of the function G(t, v ⊥ ) defined in Eq. (14). This function is shown in Fig. 8 for both Yukawa and Coulomb potentials. Choosing the values of t where the function G(t, v ⊥ ) is small for the wide enough intervals of v ⊥ we can, as we did above, calculate the distributions defined in Eq. (16) for the range of the lateral velocities in which the distribution is predominantly concentrated.  The results of these calculations of the lateral velocity distributions for different driving field strengths are shown in Fig. 9. Also shown are the results of the fitting procedures based on Eqs. (17) and (18) for the Yukawa and Coulomb potentials, respectively. The fitting expressions (17) and (18) reproduce calculated distributions fairly well. The values of the fitting parameters we obtain are indicated in Fig. 9. For the short range Yukawa potential the values of the fitting parameter β are of primary interest. As one can see, these values indeed satisfy the relation β ≈ E 0 , which one would expect basing on the SFA formula (15). We see again that in the case of the short range interaction, the shapes and widths of the lateral distributions we obtain for the time moments www.nature.com/scientificreports/ inside the laser pulse duration, coincide with the asymptotic distribution (15), describing velocity distribution at the detector. This agrees fully with the SFA-based scenario we mentioned above, according to which the lateral velocity distribution does not evolve with time once electron is ionized because of the short range character of the atomic potential assumed in the standard form of the SFA. For the Coulomb potential the fitting parameters are less informative. As we mentioned above, the long range Coulomb forces exert considerable effect on the lateral velocity distributions so that the system must evolve for a long time for the distribution to approach its limiting asymptotic form. We illustrate the character of this approach for the driving laser frequency ω = 0.057 a.u. in Fig. 10a, where we show time-evolution of the fitting parameter β in Eq. (18). The asymptotic distribution (19) and its fit (18) are shown in Fig. 10b. One can see that after the ionization event β starts decreasing, so that the lateral velocity distribution becomes progressively sharper at v ⊥ = 0 . At the end of the pulse β ≈ 1 , which is still far from the value β ≈ 0.59 which our fitting procedure (18) gives for the asymptotic distribution shown in www.nature.com/scientificreports/ Fig. 10b. We see thus the same behavior of the cusp formation that we saw in Fig. 7 for the driving laser frequency ω = 0.03 a.u. The cusp is not fully formed at the end of the pulse, the parameter characterizing the sharpness of the cusp needs more time yet to reach its asymptotic value.

Conclusion
We presented a study of the lateral velocity distribution and its development in time during the interval of the laser pulse duration. Our approach is based on the application of the notion of the joint probability. Though this notion is not very well defined in QM and may sometimes lead to unphysical (e.g. complex) probability distributions, one may find situations when such a joint probability becomes a meaningful concept. This opens a way to tackle the questions which are somewhat difficult to address, in particular, to study the characteristics of the ionized electron for times inside the interval of the laser pulse duration. As we mentioned above, for the moments of time inside the interval of the laser pulse duration, when the wave-packet describing ionized electron is not fully formed and is still partly inside the atom, it is not easy to unambiguously single out the part of the wave-function describing ionized electron from the total wave-function of the system. In the framework of our approach this ambiguity is resolved by calculating joint probability of two events, event A being detection of the lateral electron's velocity in a given volume of the electron's velocity space and event B being electron's detection by a detector placed far away from the atom. We saw, that though the joint probability thus defined cannot be used for all moments of time inside the laser pulse duration, since operators corresponding to the observables A and B generally do not commute, one can find the instances when these operators commute approximately, allowing us to define sensible joint probability distribution of the events A and B.
By following this strategy we were able to track the time-evolution of the lateral velocity distributions on the time interval between the moment of ionization and the end of the laser pulse. We saw that in the case of the short range Yukawa potential, this evolution reproduces the standard SFA scenario. The shape of the lateral velocity distributions for times inside the laser pulse is reproduced pretty accurately by the SFA formula (15), and it changes very little after the moment of ionization. This is, of course, a consequence of the short range character of the Yukawa potential, which has little influence of the motion of the ionized electron while electric field alone cannot alter velocities in lateral directions. We observe quite different behavior in the case of the Coulomb potential, which is expected since long range Coulomb force can, strictly speaking, never be neglected. Since the average of the electric field over an optical cycle is zero, the long range Coulomb force may affect electron's motion to a considerable degree even when electron is far away from the ionic core. This effects manifests itself, in particular, in the process of the time-development of the cusp in the lateral velocity distributions in the Coulomb case. We saw that the lateral velocity distributions we compute using Eq. (16) become progressively more sharply peaked near the point v ⊥ = 0 , as illustrated in Figs. 7 and 10, thus showing the development of a cusp. Nevertheless, even at the end of the pulse this development is far from being complete, the values of the parameter β characterizing "sharpness" of the cusp are still far from the value we retrieve from the asymptotic (infinite time) distribution (19).
In the present work we used as objects of study the simple one-electron targets, the hydrogen atom and its counterpart described by the short-range Yukawa potential with the same ionization potential, in the field of a linearly polarized electromagnetic wave. The theoretical framework we used can be applied for more complicated targets as well, such as heavier atomic species or multiply charge ions with larger ionization potential. We would www.nature.com/scientificreports/ expect qualitatively similar results for these targets. In the absence of the long range Coulomb force the formation of the lateral velocity distributions, as we saw, is very well captured by the SFA, which does not rely on any information about atomic target at all, except the value of the ionization potential. We believe, therefore, that for the targets with higher ionization potential, described by the short range forces, we should obtain results similar to the results we presented above for the Yukawa SR potential, with the shape of the lateral velocity distribution reproduced fairly well for times inside the laser pulse by the SFA formula (15) and changing insignificantly after the moment of ionization. For the targets with the ionic potential exhibiting long range Coulomb tail (e.g. the multicharged ions) we would expect stronger yet effect of the Coulomb field on the cusp formation. That should give us the overall picture of the cusp formation qualitatively similar to the case of hydrogen we studied above, resulting in a considerable evolution of the cusp parameter in the lateral velocity distribution, like the one shown in Fig. 10. For the systems with stronger Coulomb forces the particular values of the parameters characterizing the cusp can, of course, be completely different from those shown in Fig. 10.
We can apply this technique to other field geometries as well, in particular to the case of the circularly or elliptically polarized electromagnetic fields. The case of the driving field polarization different from the linear one is of particular interest since for elliptically polarized field the tunneling electronic wave packet possesses an initial transverse momentum due to the non-adiabatic effects 51,52 . Calculations for such field geometry are considerably more time consuming because of the necessity of performing multiple solutions of the TDSE for the much more computationally demanding case of the non-linear polarization. We plan, nevertheless, to perform such calculations in the future to study these highly interesting non-adiabatic effects.
We also note that our approach is by no means limited to the study of the lateral velocity distributions only. By choosing other observables to represent event A at times inside the interval of the laser pulse duration and applying the strategy we described above, we can study evolution of the corresponding distributions. As an illustration of this statement we consider briefly some results we obtain for another choice of the event A. Let us define A as an event consisting in finding the electron in a region around the point r 0 in space at the moment of time t. More specifically, we define the projection operator Q r 0 representing this observable as: Q r 0 = |φ r 0 (r)��φ r 0 (r)| . For |φ r 0 (r)� we use the Gaussian form: φ r 0 (r) = Ne −a(r−r 0 ) 2 ) , where N is the normalization factor and we use the value a = 2 ln 2 for the parameter a. This parameter defines the 'resolution' with which we can scrutinize the coordinate space, and it is approximately one atomic unit of length for the choice of the parameter a we made. In the Eq. (13) we can replace now the projection operator Q i with the projection operator Q r 0 , and the resulting expression will give the joint probability P(r 0 ) of finding electron ionized after the end of the pulse and finding it in a location nearby r 0 at the moment t. In other words, the joint probability thus calculated gives us a distribution of the ionized electron coordinates at times inside the laser pulse. We proceed with the calculations in exactly the same way as we did before, obtaining the results shown in Fig. 11 for the case of the short range Yukawa potential. We use the same pulse parameters as above and report the results for the peak field strengths E 0 = 0.1 a.u. To simplify the notation we drop the subscripts and we will write simply r instead of r 0 in the formulas below which should not cause confusion.
We show in Fig. 11a the function G(z, t) calculated according to Eq. (14), with projection operator Q k replaced with the operator Q r for r = (0, 0, z) . We remind that in our geometry z-direction is the polarization direction, so by choosing this direction in space we can study the coordinate distribution for the ionized electron along the laser polarization direction. One can see from Fig. 11a that for sufficiently large z and all t G(z, t) is small. Joint distributions computed according to the Eq. (13) give us, therefore, reliable distributions for the ionized electron z− coordinate for all t and at least some interval of z-values. This can also be seen from Fig. 11b,c where we show imaginary and real parts of the joint distributions P(z) computed according to Eq. (13). One can see that for the values of z of interest (near the peaks of the real parts of P(z)) imaginary part of P(z) is at least an order of magnitude smaller than the real part. We obtain thus the meaningful ionized electron's coordinate distributions for the most important and interesting intervals of z, where the bulk of a distribution is concentrated. As one can see from Fig. 11c electron ionized at the moment near the field maximum at t = 0.5T is located at z ≈ 7 a.u. This is not very different from a simple estimate for the tunnel exit location based on the energy conservation formula for an electron moving in the combined field of the ionic core and the static electric field with the amplitude E 0 (the so-called Field Direction Model (FDM) 18 ), which, for the Yukawa potential we consider and E 0 = 0.1 a.u. is z ≈ 5 a.u.
Using this approach one can obtain information about other nonlinear phenomena as well. The rescattering process, for instance, is at the core of the HHG process and is responsible for the formation of the high energy part of the ATI spectra 11 . Consider, for instance, the ATI process. By choosing event B to be detection of the electron in the ionized state and event A to be finding the electron's coordinate in a certain region of space (as we did above studying the ionized electron's coordinate distributions), we can single out the part of the wavefunction describing the ionized wave-packet, and study rescattering process in detail by following spatial and temporal evolution of the ionized wave-packet.