Instantaneous ionization rate as a functional derivative

The notion of the instantaneous ionization rate (IIR) is often employed in the literature for understanding the process of strong field ionization of atoms and molecules. This notion is based on the idea of the ionization event occurring at a given moment of time, which is difficult to reconcile with the conventional quantum mechanics. We describe an approach defining instantaneous ionization rate as a functional derivative of the total ionization probability. The definition is based on physical quantities, such as the total ionization probability and the waveform of an ionizing pulse, which are directly measurable. The definition is, therefore, unambiguous and does not suffer from gauge non-invariance. We compute IIR by numerically solving the time-dependent Schrödinger equation for the hydrogen atom in a strong laser field. In agreement with some previous results using attoclock methodology, the IIR we define does not show measurable delay in strong field tunnel ionization. There has been much recent experimental and theoretical interest in the physics of ionization, in particular the question of tunneling time. In this work, the authors derive a gauge invariant definition of the instantaneous ionization rate as a functional derivative of the total ionization probability.

T he notion of the instantaneous ionization rate (IIR) has proved extremely fruitful for understanding the physics of tunneling ionization, where the ionization regime characterized by small values of the Keldysh parameter γ ¼ ω=E 0 ffiffiffiffiffiffiffi 2 I j j p ≲ 1 1,2 . Here, ω, E 0 , and I are the frequency, field strength, and ionization potential of a target system expressed in atomic units.
Qualitatively, the fact that the IIR is a function which is sharply peaked near the local maxima of the electric field of the pulse can be used to pinpoint the most probable electron trajectory. This provides a basis for the design and interpretation of the results of well-known experimental techniques, such as attosecond streaking, or angular attosecond streaking [3][4][5] allowing one to follow electron dynamics at the attosecond level of precision. Recently, the temporarily localized ionization at the local maxima of the laser field has been used as a fast temporal gate to measure the laser field 6 .
Quantitatively, the notion of IIR underlies many successful simulations of tunneling ionization phenomena, relying on Classical Trajectory Monte Carlo method 7 or quantum trajectories 8,9 Monte-Carlo simulations. These methods become practically indispensable if the system in question is too complicated to allow an ab initio treatment based on the numerical solution of the time-dependent Schrödinger equation (TDSE). Even if a numerical solution of the TDSE is possible, these methods may provide a physical insight which is not obvious from the TDSE wave-function. Accurate quantitative calculations using these approaches, which agree well with the ab initio TDSE, have been reported in the literature 7,[10][11][12][13][14][15] . In these approaches, quantum-mechanical Keldysh-type theories 1,[16][17][18] are used to set up initial conditions for the subsequent electron motion [11][12][13] . IIR in the well-known Ammosov-Delone-Krainov (ADK) form 19,20 , or more refined Yudin-Ivanov IIR 21 provides a measure which allows us to assign probability to an ionization event occurring at a given moment of time inside the laser pulse duration.
The notion of the ionization event occurring at a given time and, correspondingly, the notion of IIR, is not free from ambiguity, however. An interpretation of the ionization event which is often used is based on the imaginary time method 18,22 . In this picture, an electron enters the tunneling barrier at some complex moment of time with complex velocity. Upon descending onto the real axis, the velocity and coordinates corresponding to the most probable electron trajectory become real, which can be interpreted as the electron's exit from under the barrier. This picture, however, cannot be taken unreservedly, since the path which descends onto the real axis is not unique and can be deformed, in principle, to cross the real time-axis at almost any given point 22 . An approach which allows us to define IIR from the solution of the TDSE has been described recently 23 . In this approach, the IIR is defined by projecting out contributions of the bound states from the solution of the TDSE. The authors found that the IIR thus defined lags behind the local maxima of the electric field, which suggests a nonzero tunneling time. A shortcoming of this definition of the IIR, however, is its nongauge-invariant character. The projection of the solution of the TDSE onto the subspace of the bound states performed during the interval of the pulse duration generally depends on the gauge used to describe atom-field interaction. Another approach which allows us to define IIR from the solution of the TDSE is based on the notion of the electron flux and was given in ref. 24 .
In the present work we describe a different approach to IIR, which is based on the notion of a functional derivative. This approach provides an unambiguous and gauge invariant definition of the IIR. We apply the IIR thus defined to the problem of the tunneling time for the process of the tunneling ionization.
Our definition does not show any appreciable tunneling delay in strong field ionization for the case of hydrogen.

Results
Total ionization probability can be considered as a functional P [E] of the waveform E(t), which maps the electric field of the pulse into a real number. For the regime of the tunneling ionization that we are interested in below, this functional is highly nonlinear and cannot be described in a closed form. A simplification is possible if we consider a waveform which can be represented as E(t) = E f (t) + δE(t), with fundamental field E f (t) and signal field δE(t). To be more specific, let us assume that the fundamental field is linearly-polarized (along z-axis, which we assume to be the axis of quantization) and is defined by the vector potential A f (t): with peak field strength E 0 , carrier frequency ω, and total duration T 1 = NT, where T = 2π/ω is an optical cycle (o.c.) corresponding to the carrier frequency ω, with N 2 N. We assume the signal field to vanish outside the interval (0, T 1 ) of the fundamental pulse duration. If the signal field δE(t) is sufficiently weak, we can write: where δP δE f ðtÞ is a functional derivative of the functional On the other hand, using the customary definition of IIR, we can write the probability of ionization driven by the combined field E(t) = E f (t) + δE(t) as: where W inst (E(t)) is the IIR, which by definition depends on the instantaneous value E(t) of the electric field. We introduced the notation P inst in Eq. (3) to emphasize that this expression pertains to the notion of the IIR. Eq. (3) gives a very particular case of the functional P[E]-clearly not every functional can be represented in this way. From Eq. (3), we obtain the change of the ionization probability P inst due to the presence of the weak signal field: Taking into account that δE(t) in Eqs. (2) and (4) is arbitrary (provided it is small and vanishes outside the interval (0, T 1 )), one can see that, if the treatment based on the notion of the IIR is justified (i.e., using Eq. (3) with IIR W inst (E(t)) depending only on the instantaneous value of the electric field gives an accurate value for the total ionization probability), one must have: We note that the quantity on the right-hand-side (r.h.s.) of (5) is an ordinary derivative which depends on time only through the instantaneous value of the electric field. The quantity on the left-hand-side (l.h.s.), on the other hand, is a functional derivative, which depends on time (through the time moment t at which differentiation is performed), and on the waveform E(t), i.e., on the complete history, past, and future, of the pulse. Suppose, for instance, that we represent the waveform (1) as a Fourier series: where Ω = 2π/T 1 . In general, the functional derivative on the l.h.s. of Eq. (5) depends on the whole set a m , while the ordinary derivative on the r.h.s. is a function of the particular combination of these coefficients. To express this differently, suppose that we fix the pulse shape in Eq. (1) and allow only the field amplitude E 0 to vary. The l.h.s. of Eq. (5) would then generally be a function of two variables E 0 and t, while the r.h.s. would depend only on a particular combination of E 0 and t, which gives the instantaneous field strength E(t). Exact equality in Eq. (5), therefore, cannot be achieved in general. This clearly demonstrates the approximate character of the notion of the IIR. The more accurate expression (2), based on the functional derivative of the ionization probability functional, depends not only on time, but on additional variables describing the waveform as well. For brevity, we will dub below the functional derivative in Eq. (2) an "exact ionization rate", though as one can see from Eq. (5), it rather corresponds to the derivative of the IIR with respect to the electric field.
Below, we will fix the functional form of the fundamental pulse shape in Eq. (1) and allow only the peak field strength E 0 of the fundamental field to vary. In Fig. 1 below, presenting results of the TDSE calculation, we show and analyze, not the functional derivative itself, but a proportional quantity, viz. the first variation of the ionization probability δP(E 0 ,τ), obtained if we substitute expression (17) for the signal field in Eq. (2).
In Fig. 1, we show the first variation δP(E 0 , τ), computed following the numerical procedure outlined above as a function of time τ and electric field strength E 0 for pulses (1) with total duration T 1 of one ( Fig. 1a) and two ( Fig. 1b) optical cycles.
In Fig. 2a, b, we show the functional derivative δP δEðtÞ obtained using the ADK expression 19,20 for the IIR. Since the ADK ionization probability assumes an instantaneous relation between the field strength and the probability of ionization at any given time, this naturally leads to a definition of the IIR following Eq. (4). The functional derivative in this case is just the derivative of the ADK ionization rate with respect to the instantaneous electric field. Plots showing results of the TDSE and ADK calculations show good qualitative agreement, with the functional derivative sharply peaked in the vicinity of the electric field maximum. We can expect such peaks to occur for every local maximum. The magnitudes of the secondary peaks, however, can be much smaller compared to the peak corresponding to the principal maximum of the field. This, in particular, is the case for the twocycle driving pulse shown in Fig. 1 (TDSE results) and Fig. 3a (ADK results). On a linear scale, the secondary peaks occurring in the vicinity of the secondary field maxima (at τ ≈ T/2 and τ ≈ 3T/ 2) are practically absent both in the TDSE and ADK cases. To see these peaks clearly we show in Fig. 3b the quantity δP δEðtÞ 1=3 . These properties are, of course, a direct consequence of the sharp dependence of the ionization rate on the electric field. This dependence is exponential in the case of the ADK IIR. The dependence on the electric field of the IIR that we define in the present work will be discussed below.
A more detailed picture of the IIR's emerges by looking at the contour plots (i.e., lines of equal elevation) of δP(E 0 , τ) in the E 0 , τ-plane. Of particular interest is, of course, the behavior of the IIR near the local maxima of the electric field where the IIR attains its largest values. We will concentrate, therefore, on the region of τ-values close to the highest maximum of the electric field strength.
We will consider below only pulses with a total duration of one optical cycle. The results for pulses with duration of 2 optical cycles have been found to be essentially the same as for the one optical cycle duration.
Since the ADK IIR is a function of the instantaneous electric field only, lines of constant elevation in this case are just the lines In ref. 23 it was found that the IIR defined in that work by projecting out bound states contributions from the TDSE wavefunction lags behind the electric field of the pulse, which was interpreted as a manifestation of the tunneling delay. Analogous result was reported in ref. 24 , where the authors monitored the probability current density at the (albeit adiabatic) exit point, as calculated by a one-dimensional TDSE solution. The point of view that the tunneling delay has a nonzero value has also been expressed in the works [25][26][27][28] , as opposed to papers advocating the view of tunneling as an instantaneous process [29][30][31] .
In the approach we pursue, a nonzero tunneling delay would lead to an asymmetry of the contour lines about the midpoint of the pulse, as can be seen from the following argument. The variation δP(E 0 , τ) is a function of E 0 and τ. In the vicinity of a local field maximum, we can introduce electric field and its derivative (E(τ), E′(τ)) as another pair of independent variables. The variation becomes a function δP(E(τ), E′(τ)). Expanding δPv(E(τ), E′(τ)) in the vicinity of the field maximum and keeping only first order terms in small E′(τ), one can then write: δPðEðτÞ; E′ðτÞÞ % δPðEðτÞ; 0Þ þ AðEðτÞÞE′ðτÞ with some A(E(τ)), Δ, and function f(x). The lag Δ in the latter equation can be interpreted as the tunneling delay. As can be seen from Eq. (7), in a small vicinity of the field maximum, the lines of constant elevation δP(E(τ), E′(τ)) = const locally assume the form E(τ − Δ) = const. Nonzero Δ, therefore, makes the contour lines asymmetric about the midpoint of the pulse. Figure 4a shows contour lines of δP(E 0 , τ) obtained from TDSE calculation for hydrogen with the Coulomb potential. The absence of any appreciable asymmetry in the Figure suggests that our definition of the IIR gives us an essentially zero time delay for the Coulomb potential.
We considered also the case of the Yukawa potentials VðrÞ ¼ ÀAe À r a =r with different screening parameters a. For every a we adjusted the value of A so that the resulting ionization potential for the ground state was always 0.5 atomic units (a.u.), corresponding to hydrogen. Results for the Yukawa potentials shown  to better reveal the peaks corresponding to the secondary field maxima in Fig. 4b, c again show an absence of any appreciable asymmetry of the contour lines, and consequently the time delay.
We used above the functional derivative δP δEðtÞ , which as Eq. (5) shows, provides a definition for the derivative of the IIR, dW inst dEðtÞ , rather than for the IIR itself. To obtain the IIR provided the functional derivative δP δEðtÞ is known, one can proceed as follows. Let us write the the functional P[E] (we remind readers that it is the total ionization probability considered as a functional of the laser field) as: WðEðtÞ; E′ðtÞ E ðnÞ ðtÞ; tÞdt; The IIR W inst (E(t)) in (3) depends, by definition, only on the instantaneous value of the electric field E(t). The function W in Eq. (8) is used to represent the exact functional P[E]. It may depend, therefore, not only on the electric field but on the (in principle infinite) number of its higher order derivatives and time, as reflected by the notation in Eq. (8). From the Eq. (8) where we used integration by parts to handle the electric field derivatives (the usual procedure which in the case of W in Eq. (8) depending only on E(t) and E′(t) gives the standard from of the Euler-Lagrange equations).
To advance further, we can proceed as follows. We can assume a trial form W trial (e), depending on the set e of variational parameters, for the function W. Using this trial function, we compute a trial functional derivative δP trial δEðtÞ as prescribed by the Eq. (9) for the field strengths E 0 we considered above and form a functional: In practice, of course, both E 0 and τ in Eq. (10) assume a discrete set of values for which we performed the TDSE calculations we described above, so the double integral in Eq. (10) is, in fact, a double sum. Minimizing expression (10) (we use the wellknown gradient decent method) we find the trial parameters e and hence a variational approximation to the exact IIR in Eq. (8).
We will consider the following trial IIR: with variational parameters e 1 -e 3 . The leading small-E term −2/3E(t) in the exponential reproduces, of course, the leading term in the ADK formula 19 , the coefficient e 1 with the logarithmic term would be −0.5, had we used the ADK expression. It is known, however 22 , that the Coulomb potential modifies considerably the pre-exponential factor in the ADK formula, we treat, therefore, e 1 as a fitting parameter.
With the fitting parameters given by the fitting procedure, we can compute the total ionization probability P trial (E) using Eq. (8) and judge thereby the quality of the variational approximations to the IIR. We illustrate this strategy in the case of a hydrogen atom with Coulomb potential driven by a single cycle pulse. Expression (10) was minimized using a set of τ-values in the interval of the pulse duration and peak field strengths in the interval between 0.055 and 0.07 a.u. for which we performed the TDSE calculations we described above.
The results are illustrated in Fig. 5, where we show the total ionization probability P(E) obtained from the TDSE calculation (Fig. 5a), IIR's obtained using the fitting procedure we described above (Fig. 5b), and the relative errors (P trial (E) − P(E))/P(E) which fits (a-c) in Eq. (12) give for the total ionization probability for the interval of the electric field strengths we consider (Fig. 5c).
It is noteworthy that inclusion of the first derivative in the expression for the trial function has no appreciable effect on the  quality of the fit. The value for the fitting parameter e 2 in the expression for W b trial , which our fitting procedure returns is very small (on the order of 10 −6 ). This result agrees with the absence of the tunneling delay for Coulomb potential which we found above by studying lines of constant elevation for δP(E 0 , τ). Indeed, were the term with the first derivative present with a nonnegligible coefficient in the expression W b for the IIR, we could rewrite it (keeping only first order terms in E′(t) as as we did above in deriving Eq. (7) for the first variation of the ionization probability) as: to obtain an ADK-type expression for the IIR without the first derivative but with a shifted argument E(t − Δ) of the electric field. This would imply a nonnegligible tunneling delay. We applied our method above in the case of the ionization in the tunneling regime (with Keldysh parameter γ ≲ 1/3 for the field parameters we considered). We should note that, at least formally, the procedure of calculating the functional derivative and retrieving the IIR from its values can be applied in the multiphoton or in the intermediate regimes as well. Indeed, by construction, the time-integral of the IIR in the Eq. (8) should give the total ionization probability regardless of the ionization regime. In the multiphoton regime, however, we cannot expect ADK-type trial expressions in Eq. (12) to be good trial guesses. In the multiphoton case and in the intermediate regime between the multiphoton and tunneling case, we can expect the role of the higher order time derivatives in Eq. (8) to be more important.
The theory and calculations we presented above were based on the delta-function signal pulse δE(t, τ) = αδ(t − τ) which has a DC-component (integral of the electric field along the interval of the pulse duration has nonzero value). We could use a more easily obtainable experimentally attosecond sine-pulse, with an electric field which can be approximated as a derivative of the deltafunction: δEðt; τÞ = α ∂ ∂t δðt À τÞ. For such a signal pulse, our Eq. (2) gives: Taking into account the correspondence between the functional derivative δP δE f ðtÞ and the derivative of the IIR we established in (5), we see that variation δP of the ionization probability in this case can be estimated as (provided the notion of the IIR depending only on the instantaneous value of the field is applicable): Instead of probing the first derivative of the IIR as we did above, we would be probing in this case its second derivative. The general behavior and the lines of the constant elevation in the case of the function given by the r.h.s. of the Eq. (13) are shown in Fig. 6a, b for the ADK IIR, and are quite different from those we have shown above in Figs. 1 and 4 for the case of the deltafunction pulse.
The reason for this difference is the factor E′(t) in Eq. (14) due to which the r.h.s of this equation vanishes near the field extrema. This, however, could be advantageous if we wish to obtain information about the IIR far from the field maxima. Indeed, for the variational approach based on Eq. (10), the time values close to the field maxima, where the quantity δP δEðtÞ has a large magnitude, contribute mostly to the variational functional. We can develop a similar variational procedure which would allow us to reconstruct the IIR from the quantity on the r.h.s. of Eq. (13). In this case, because of the factor which makes this quantity small near the field maxima, we could hope to determine the IIR more accurately for the time values farther away from the field extrema.

Discussion
To summarize, we described an approach which allows us to define the IIR as a functional derivative of the total ionization c The relative errors (P trial (E) − P(E))/P(E) which the fits give for the total ionization probability probability. This approach provides an unambiguous definition of the IIR. In particular, it is based on directly measurable quantities, such as the total ionization probability and the waveform of the pulse, which makes it gauge invariant. We applied the IIR thus defined to the problem of the tunneling time for the process of the tunneling ionization of systems with Coulomb and Yukawa potentials as examples. In agreement with some previous results using attoclock methodology (which assume the most probable electron trajectory to begin tunneling at the peak of the laser field), the IIR we define does not show any appreciable tunneling delay in strong field ionization for the case of hydrogen. Finally, we would like to note that we have pursued several goals in doing this work. One was to shed some light on the tunneling time problem. We tried to do that by giving a precise meaning to the concept of the IIR which, by itself, is an important and useful notion. As explained in the introductory section of the manuscript, the notion of the IIR is in the core of many successful approaches to simulations of strong field ionization. Use of a notion which has not been clearly and unambiguously defined for such purposes can lead to considerable confusion. The definition we give is mathematically clear, unambiguous, and useful.

Methods
To compute the exact ionization rate in Eq. (5) for a given moment of time τ and a given waveform E f (t), one can calculate numerically the modulation of the ionization probability δP using E f (t) as a fundamental field and employing a special form δE(t, τ) = αδ(t − τ) of the signal field, containing the Dirac delta-function. We did this by solving numerically the TDSE for a hydrogen atom in the presence of the pulse (1). We considered pulses with a low carrier frequency of ω = 0.02a.u. We need to work with low frequencies to stay within the framework of the adiabatic theory, which makes the notion of the IIR at least qualitatively applicable.
We report below results of the calculations with total pulse durations T 1 of one or two optical cycles respectively, and various peak field strengths, E 0 (chosen so as to remain in the tunneling regime of ionization).
The solution of the TDSE has been found by representing the wave function as a series in spherical harmonic functions, and discretizing the resulting system of radial equations on a grid with the step-size δr = 0.05a.u. in a box of the size R max = 700a.u., which was sufficient for the pulses of a short duration we consider below. More details about the numerical procedure used to solve the TDSE can be found in refs. 32,33 .
Total ionization probability, which we need for the practical implementation of the definition (5), is found by decomposing the wave function at the end of the pulse as: where ϕ ¼QΨ T 1 ð Þ, χ ¼Î ÀQ À Á Ψ T 1 ð Þ, andQ is the projection operator on the subspace of the bound states of the field-free atomic Hamiltonian. Total ionization probability P can be found then as the squared norm P ¼ χ k k 2 . The projection operatorQ is obtained by numerically computing (employing the same grid we used to solve the TDSE) the bound states nlm j i of the field-free Hamiltonian: where we use the fact that, for the geometry we employ, we need to consider only the states with zero angular momentum projection. We retain in Eq. (16) all the eigenvectors that we obtain for 0 ≤ l ≤ L b . We use L b = 12 in the calculations below after checking that results are well converged with respect to this parameter. We note that the definition of the ionization probability as the squared norm of χ is formally equivalent to the commonly-used definition in terms of the projection on the ingoing scattering states ϕ À k : P ¼ . Indeed, substituting the decomposition (15) into this equation, and using orthogonality of χ to the subspace of the bound states, we obtain again P ¼ χ k k 2 . In practice, the prescription based on the calculation of the norm of the part of the wave-function describing the ionized wave-packet in the coordinate space, as encapsulated in Eq. (15), is preferable. The reason for this is that the calculation of the total ionization probability by projecting Ψ(T 1 ) on the set of the scattering states implicitly presumes the strict orthogonality of the scattering and bound atomic states. We are interested below not in the ionization probability itself, but rather in its relatively small variations due to the weak signal field. Even small nonorthogonality of the scattering and bound states, unavoidable in numerical calculations, may lead to a significant loss of precision.
The definition of exact ionization rate in Eq. (2) is based on the physically observable quantities: the electric field of the pulse and modulation of the total ionization probability. It is clearly gauge invariant. It is immaterial, therefore, which gauge describing the atom-field interaction is used when solving the TDSE (provided, of course, that a sufficient level of numerical accuracy has been achieved). We used the length (L-) gauge. Convergence of the expansions in spherical harmonic functions that we use to represent the solution of the TDSE was achieved by including spherical harmonic functions of the rank up to L max = 70. In numerical calculations we have to regularize, of course, the expression δE(t, τ) = αδ (t − τ) for the signal field. We used the regularization: with ε = T/1000 (here T is an optical period (o.c.) corresponding to the carrier frequency we use), and α = 0.001a.u. The value of α is to be small enough to warrant omission of the higher order functional derivatives in the Taylor expansion for the ionization probability functional. That our choice of this parameter achieves this goal can be seen from the plot in Fig. 7, where the total ionization probability is shown as a function of the parameter α for the Yukawa potential VðrÞ ¼ ÀAe À r a=r with the screening parameter a = 20, and A = 1.05a.u. (the parameter A was adjusted so that the ionization potential equals that of the hydrogen atom), total duration of the fundamental pulse 1 o.c., its peak field strength E 0 = 0.061a.u. The parameters in Eq. (17) were τ = 0.5 o.c., ε = T/1000. One can see from Fig. 7 that for α ≲ 0.002 the ionization probability is very well approximated by a linear function of α. Least squares linear fit P(α) = P(0) + cα (here P 0 is the ionization probability in the absence of the signal field, c is a fitting parameter) performed using all data-points with α ≤ 0.002 allows us to quantify this statement. For α ≤ 0.002 the relative error of the linear fit is less than one percent. We may conclude, therefore, that for α = 0.001 the relative contributions of the terms containing higher order functional derivatives in the Taylor expansion of the ionization probability functional are of the order of 1%. To ascertain that the parameter ε is adequately chosen, we performed a series of calculations, using the same pulse and Yukawa potential parameters for different values of ε. Results, summarized in Table 1, show that the relative error introduced by the finite spread of the Gaussian in Eq. (17) is of the order of 10 −5 for the value of ε = T/1000 we use. We have ascertained thus that the combined relative error due to the omission of higher order functional derivatives and the finite width of the Gaussian regularization (17) is of the order of 1%. To conclude the discussion of the regularization procedure employed in the calculations, we note that for the Coulomb potential, the ionization probability is typically larger than that for the Yukawa potential, so the relative accuracy of the calculation in the case of the Coulomb potential (for the same signal pulse parameters) would be even higher.
It is worth noting that, even for ε = T/30, the relative error we get for the ionization probability is of the order of 1%. Full width at half maximum of the pulse (17) for this value of ε is about 300 attoseconds. Delta-function-like pulses of such duration have already been produced in the laboratory 34 , which makes possible experimental measurements of the IIR which we define in the present work.

Data availability
All relevant data are available from the authors upon request.