Simple man model in the Heisenberg picture

We describe an approximate solution to the Heisenberg operator equations of motion for an atom in a laser field. The solution is based on a quantum generalization of the physical picture given by the well-known Simple Man Model (SMM). We provide justification of the plausibility of this generalization and test its validity by applying it for the calculation of the coordinate and velocity autocorrelation functions which, to our knowledge, have not been studied before in the context of the strong field ionization. Both our model and results of the ab initio numerical calculations show distinct types of correlations due to different types of electron's motion providing a useful insight into the strong field ionization dynamics.


Introduction
An atom exposed to a strong laser field can be ionized. The foundation of the quantum theory describing this process in strong fields has been laid out in the seminal paper by Keldysh [1] (also known as the strong field approximation or the SFA theory). The Keldysh theory introduces the well-known classification of the ionization phenomena based on the value of the Keldysh parameter γ = ω 2|ε 0 |/E (ω, E 0 and |ε 0 | are the frequency, field strength and ionization potential of the target system expressed in atomic units). The ionization regime corresponding to the values γ >> 1 is known as the multiphoton regime. The opposite limit γ 1 is known as the tunneling regime [2]. Depending on the ionization regimes, the ionization process is described in drastically different ways [1,2].
The tunneling regime is particularly interesting since many interesting and important phenomena occurring in this regime, such as the high harmonic generation (HHG), the attosecond pulse generation and the above threshold ionization (ATI), can be understood using fairly simple physical picture, the so-called simple man model (SMM) [2][3][4][5][6]. In the framework of this model the electron's motion after the ionization event is described using classical equations of motion for an electron in the presence of the laser field, neglecting the effect of the atomic potential. These equations can be easily solved leading to the testable predictions (such as, e.g., the maximum energy of the direct electrons in the ATI process, or the cutoff photon energy for the HHG process) which often agree remarkably well with results of the ab initio quantum simulations.
For the case of the electron's motion in the electromagnetic field in the absence of any other forces, not only the classical equations of motion but also their quantum * Electronic address: igorivanov@ibs.re.kr † Electronic address: kyungtaec@gist.ac.kr counterpart-the Heisenberg operator equations of motion for the operators of coordinate and momentum, can be solved fairly easily. Indeed, the solutions are practically identical. The question arises can this fact be somehow exploited? The great utility of the SMM, on one hand, and the fairly simple expressions for the quantum solutions to the Heisenberg equations of motion for the SMM physical settings suggest that we may try to somehow extend the SMM into the quantum domain.

Theory
We consider a hydrogen atom interacting with the laser pulse. In the Schrödinger picture the system evolves according to: whereĤ 0 =p/2 − 1/r, and we use velocity form H int (t) = A(t) ·p + A(t) 2 /2 for the interaction operator. Initial state of the system is Ψ(r, 0) = φ 0 (r)-the ground state of the hydrogen atom. Laser pulse is linearly polarized (along the z− axis) and is defined by the vector potential: with peak field strength E 0 , carrier frequency ω, and total duration T 1 = T , where T = 2π/ω is an optical cycle (o.c.) corresponding to the frequency ω. We use the dipole approximation to describe atom-field interaction, so the expression (2) for the vector potential does not contain the spatial variables. We will consider below the pulses with various peak field strengths E 0 . We use a short pulse of one o.c total duration, so that electric field of the pulse has a well-defined global maximum to facilitate study of the ionization dynamics. Fig. 1 shows arXiv:1910.05910v1 [quant-ph] 14 Oct 2019 the shape of the pulse (2) for E 0 = 0.0534 a.u. (intensity of 10 14 W/cm 2 ).
For the reader's convenience we recapitulate first some well-known facts and introduce few definitions which we will use below. Eq. (1) describes ionization process in the Schrödinger picture. In the Heisenberg picture the wave-function does not evolve in time while the operators evolve according to: whereQ stands for either coordinate or momentum operator, andÛ (t, 0) is the time-evolution operator. In the following we will reserve the subscript H for the operators in the Heisenberg picture, operators without subscripts will correspond to the Schroödinger picture. We will use also the time-dependent operatorsQ 0 (t) and Q V (t) which we define as: whereÛ 0 (t, 0) = exp −iĤ 0 t is the field-free atomic evolution operator, and is the so-called Volkov time-evolution operator [2,7] (hence the subscript 'V' in Eq. (4)). This time-evolution operator describes evolution driven by the Volkov Hamil-tonianĤ V (t) =T +Ĥ int (t) (T is the kinetic energy operator) of a free electron in the field of the pulse (2). We assume that bothQ 0 andQ V can be calculated numerically or analytically. Introducing a complete set of the eigenstates |n ofĤ 0 , one can write forQ 0 (t): where summations run over the spectrum (including the continuous spectrum) ofĤ 0 with eigenenergies ε n . As far asQ V (t) is concerned, the Heisenberg operator equations of motion can easily be solved for the case of the Volkov Hamiltonian, giving for the physically interesting cases of the coordinate and kinetic momentum operators the simple expressions: where r andp are the Schrödinger coordinate and momentum operators. In the atomic units system we use the kinetic momentum coincides with the velocity. We employed, therefore, the notationv V (t) in (7). We would like to remind also that we use the velocity gauge to describe atom-field interaction, hence presence of the vector potential in the second equation (7). Equations (7) look exactly like their classical counterparts used in the simple man model. This fact motivated us to try to incorporate them into a quantum version of the SMM, relying on the Heisenberg picture. The exact Heisenberg equations of motion for operatorsr H (t) and v H (t) following from the definitions (3): (hereV (r H (t)) = −1/|r H (t)| for hydrogen atom, [Â,B] =ÂB −BÂ denotes the commutator of operatorŝ A andB, andv H (t) =p H (t) + A(t) are, of course, too complicated to be solved exactly (or even numerically). We can try, however, to find an approximate solution using the time-dependent coordinate and momentum operators we introduced above in Eq. (6) and Eq. (7) and the physical insight we get from the classical SMM.
The predictive power of the SMM shows that the simple expressions for the classical electron's coordinate and velocity in the laser field which neglect all atomic interactions, can be used with success to explain many an ionization phenomena. It is natural, therefore, to use the quantum counterpart of these classical expressions given by the Eq. (7) in trying to construct an approximate solution to the Heisenberg equations of motion (8). In the limit of the vanishing field the solution should, of course, reduce to the field-free operatorsr 0 (t) andp 0 (t) obtained by substituting the Schrödinger coordinate and momentum operators, respectively, for the operatorQ in Eq. (6). This reasoning suggests the following tentative expression for the approximate solution to the Heisenberg equations of motion for the coordinate operator: On the right-hand side of this equationr 0 (t) is the field-free Heisenberg coordinate operator obtained using Eq. (6), r andp are time-independent Schrödinger operators, and we absorbed the c-number term in the first of the equations Eq. (7) into the c-number termr(t) in Eq. (9). That this term coincides with the expectation value of the coordinate can be easily seen from the Eq. (9) by noting that in the Heisenberg picture this expectation value is just the expectation value of the operatorr H (t) obtained using the wave-function of the ground state of atomic hydrogen. The expectation values ofr 0 (t), r and p in this state are all zero. The real function α(t) (it must be real to preserve hermicity of the coordinate operator) in Eq. (9) is yet unknown. Basing on the general structure of the Eq. (9) and the simple physical picture of ionization provided by the SMM, we may say that the first term on the r.h.s. of the Eq. (9) describes atomic and the second term describes ionized electrons. The function α(t) then must be related to the ionization probability. This assumption can be justified by the following reasoning.
Using the definition for the Heisenberg coordinate operator and the Dyson equation for the evolution operator U (t, 0), corresponding to the partitionĤ =Ĥ 0 +Ĥ int (t) of the total Hamiltonian: (10) one obtains the following approximate equation for the Heisenberg coordinate operator: where h.c. here and below stands for the Hermitian conjugate. In deriving Eq. (11) we neglected the term quadratic inĤ int (τ ) originating when the r.h.s of Eq. (10) is substituted into the equation defining time-dependent operators in the Heisenberg picture. Since the second term on the r.h.s. of the Dyson equation (10) is related to the ionization amplitude and is small for the not too strong fields which we consider below, the neglect of the terms quadratic inĤ int (τ ) is legitimate. From the Dyson equation written in a different form (which corresponds to the partitionĤ =Ĥ V +V atom of the total Hamiltonian): (12) where V atom is the atomic potential, one can see that replacing the operatorÛ 0 (0, t) on the r.h.s of the Eq. (11) withÛ V (0, t) would result in additional terms bilinear in H int (τ ) andV atom . Assuming again that these terms can be neglected comparing to the leading term, we can write: Consider the combination of the operatorŝ U V (0, t)rÛ (t, τ ) on the r.h.s of the Eq. (13). We can write for this combination: where we used multiplication properties (such aŝ is the Volkov coordinate operator introduced by Eq. (7). Furthermore, in deriving Eq. (14) we replaced in the third line the Volkov prop-agatorÛ V (0, t) with the total propagatorÛ (0, t). This is the replacement analogous to the one usually done in the derivation of the Keldysh expression for the ionization amplitude in the framework of the SFA (by replacinĝ U (0, t) on the r.h.s of the Eq. (10) with Volkov propaga-torÛ V (0, t) one can obtain the well-known expression for the Keldysh ionization amplitude). The justification of this replacement is again the assumption, which we have been using systematically, that the terms of the higher order (quadratic inĤ int (τ ) in this instance) can be neglected comparing to the leading term. From Eq. (13) and Eq. (14) we obtain: Introducing the complete set |n of the eigenstates of the field-free Hamiltonian, the operator expression under the integral sign in Eq. (15) can be written as: where Ψ n (τ ) =Û (τ, 0)φ n , φ m (τ ) =Û 0 (τ, 0)φ m = e −iεmτ φ m , φ m and ε m are, respectively, the eigenstates and eigenenergies of the field-free atomic Hamiltonian. The overlaps T nm = Ψ n (τ )|φ m (τ ) form the matrix of transition amplitudes between various field-free states of the atomic Hamiltonian. Substituting Eq. (16) into Eq. (15) we obtain: where the transition operatorT is given by: Eq. (17) is still pretty complicated. To advance further we have to make an assumption about the operator T in Eq. (18). This operator is clearly related to the transitions and ionization probabilities, becoming a zero operator in the absence of the electric field. Having in mind that we are developing a model based on the simple physical picture we will replace this operator with a c−number function which, taking into account the properties of the operatorT , can be considered proportional to the ionization probability. The Heisenberg operator r H (t) obtained in this way is then a sum of two terms, the field-free atomic time-dependent operatorr 0 (t) and the time-dependent Volkov coordinate operatorr V (t) entering the expression with the weight proportional to the ionization probability. Substituting the expression (7) for r V (t), we obtain: wherer 0 (t) is the field-free atomic time-dependent coordinate operator which can be obtained using Eq. (6), r and p are usual coordinate and momentum operators in the Schroödinger representation, Athe vector potential, and α(t) is a function which we assume to be proportional to the ionization probability. Taking into account that expectation value of the coordinate operator in the Heisenberg picture is just the matrix element φ 0 |r H (t)| φ 0 (where φ 0 is the initial state of the system which is the ground state of hydrogen in our case), and that expectation values of bothr 0 (t), r andp operators vanish in the ground state of hydrogen, we may absorb the c-number term in Eq. (19) into the expectation value, which gives us the Eq. (9) for the coordinate operator in the Heisenberg representation.
The expression for the the velocity operatorv H (t) follows from Eq. (9) by the time differentiation as indicated in the first of the Eq. (8). The canonical momentum can then be found asp H (t) =v H (t) − A(t). We note that while this procedure gives us the Hermitian Heisenberg coordinate and momentum operators, it does not preserve the correct commutation relations [p (n) n , This is a consequence of the fact that the transformation of the Schroödinger pair of operators r,p to the Heisenberg pairr H (t),p H (t) described by the Eq. (9) is not unitary. The lack of unitarity, however, should not be considered as a serious drawback. Many widely used approaches to the description of the ionization in strong fields, e.g., the SFA approach have a similar problem. The SFA is based on the Schroödinger picture and the non-unitarity of this method manifests itself as a non-unitary evolution of the state vector [8]. Consequently, the total sum of all the probabilities is not conserved and may not sum up to unity in the SFA. The great success and utility of the SFA show, however, this of unitarity of the method does not constitute a major impediment. We will study below some consequences and testable predictions we may derive from Eq. (9).

Results
An advantage which the Heisenberg picture offers, is the natural way in which the many-time correlation functions, containing detailed information about evolution of the system, can be introduced. We will be interested below in the two-time autocorrelation functions, defined as follows: where φ 0 is the ground state of the hydrogen atom, Q(t) is an operator in the Heisenberg picture andQ(t) is the expectation value ofQ(t). Since we consider only the ground state of hydrogen as the initial state in the present work, we will omit below φ 0 in the formulas, implicitly understanding that Â = φ 0 |Â|φ 0 for any operatorÂ. We will consider below as two examples the autocorrelation functions C zz (t 2 , t 1 ) and C vz,vz (t 2 , t 1 ) for the components of the coordinate and velocity operators in the direction of the laser field.
On one hand, we can compute the autocorrelation functions in an ab initio way without any approximations using the well-tested procedure [9] we use to solve the time-dependent Schrödinger equation (TDSE). Considering C zz (t 2 , t 1 ) as an example, using definition (3), and employing the well-known properties of the timeevolution operators, we may write: Calculation of this expression requires, thus, first propagating the TDSE starting with Ψ(0) = φ 0 -the ground state of hydrogen, on the interval (0, t 1 ), thus obtaining the state vector Ψ(t 1 ) at t = t 1 . Acting with the (Schrödinger) operator z on this vector we obtain the wave-function Ψ 1 (t 1 ) = zΨ(t 1 ). Ψ 1 (t 1 ) is further propagated on the interval (t 1 , t 2 ) yielding the wave-function Ψ 1 (t 2 ), from which we obtain Ψ 2 (t 2 ) = zΨ 1 (t 2 ). Finally, the autocorrelation function can be found by projecting Ψ 2 (t 2 ) on the state vector Ψ(t 2 ), obtained by solving the TDSE with the initial condition Ψ(0) = φ 0 on the interval (0, t 2 ). These calculations were performed using the numerical procedure [9] for the solution of the TDSE for hydrogen atom in presence of the laser field given by Eq. (2). Results of the ab initio TDSE calculations of the real and imaginary parts of C zz (t 1 , t 2 ) for different field strengths are shown in Fig. 2, Fig. 3.
The autocorrelation functions C zz (t 1 , t 2 ) shown in the Fig. 2, Fig. 3 show two distinctly different types of correlations. The correlation patterns for small fieldsthe diagonal stripes in the Fig. 2 Fig. 3 are reminiscent of those we would obtain for the field-free Heisenberg operators. The origin of the stripes is clear from the expression Eq. (6), from which we can obtain: where |0 is the hydrogen ground state. We can also explain the correlations pattern in the Fig. 2, Fig. 3 qualitatively using a simple classical picture of the atomic periodic motion. Electron's positions at the moments t and t + kT a , where T a is the period of the electron's orbital motion and k is an integer, are clearly correlated, analogously, the electron's positions at the moments t and t + kT a /2 are anticorrelated, hence the pattern of the maxima and minima of the autocorrelation function running diagonally. Strictly speaking, this simple classical explanation is applicable to the hydrogen atom with some reservations, since the different terms in the Eq. (6) are not equally spaced in energy. Consequently the frequencies in Eq. (22) are not integer multiples of some base frequency, as the simple classical picture we alluded to above implies. Had we considered instead of the hydrogen atom the harmonic oscillator with frequency Ω and equally spaced energy levels, Eq. (6) would result in the well-known relation for the coordinate operator in the Heisenberg pictureẑ H (t) = z cos Ωt +p z /Ω sin Ωt [10]. For suchẑ H (t) we would have obtained a pattern of equally spaced maxima and minima of the autocorrelation function running diagonally, for which the classical explanation would be perfectly adequate. Nevertheless, we shall use this line of arguments based on the picture of the electron's periodic motion for the purposes of the qualitative analysis even for hydrogen, having in mind that for this system this picture may have some limitations. In classical terms, thus, the pattern of correlations in the field-free case just reflects the periodic orbital motion of the electron. For non-zero external fields this field-free pattern of correlations becomes first perturbed (for E 0 = 0.04 a.u.) and then is almost completely superseded by a different pattern for E 0 = 0.06 a.u., which exhibits horizontal and vertical stripes rather than the diagonal ones, and which shows a good deal of correlated motion for both t 1 , t 2 near the end of the pulse. We shall see below that this new pattern of correlations induced by the field can be explained in considerable detail by our model based on the Eq. (9).
Using Eq. (9) and Eq. (20) we obtain for the autocorrelation function C zz (t 2 , t 1 ): where C 0 zz (t 2 , t 1 ) is the field-free autocorrelation function,ẑ 0 (t) is the z− component of the the field-free Heisenberg coordinate operator obtained using Eq. (6), z and p z are the usual time-independent Schrödinger operators of the z− components of coordinate and momentum. Calculations of the expectation values appearing in the second line of Eq. (23) can be done as follows. Considering the productẑ 0 (t 2 )p z as an example, we may write using the definition ofẑ 0 (t): ẑ 0 (t 2 )p z = e iεt2 zφ 0 |Û 0 (t 2 , 0)|p z φ 0 , whereÛ 0 (t 2 , 0) is the field-free atomic evolution operator. Calculation of this matrix element requires thus field-free propagation of the initial statep z φ 0 on the interval (0, t 2 ), which can be done without difficulties using the numerical procedure we use to solve the TDSE. Calculations of the expectation values in the third line of Eq. (23) does not pose any difficulties. There is one more ingredient which we have to provide to use Eq. (23), it is the function α(t). As we have noted above, the reasoning based on the physical picture of the ionization process suggests that α(t) should be related to the ionization probability P (t). We will use, therefore, the expression α(t) = cP (t) for the function α(t) in Eq. (9) and Eq. (23), where P (t) is the ionization probability, and c is a constant factor. We compute P (t) as Yudin-Ivanov instantaneous ionization rate (YI IIR) [11]. Our model, thus, has one free parameter c which we can vary to achieve better agreement between the ab initio C 0 zz (t 2 , t 1 ), and the C 0 zz (t 2 , t 1 ) we obtain from Eq. (23). A comparison of the results we obtain in this way, using our model and the ab initio TDSE results, is shown in Fig. 4, Fig. 5 . To reveal more detail we use a nonlinear scale in the Figures, presenting fractional power of Re(C zz (t 1 , t 2 )) and Im(C zz (t 1 , t 2 )). We can see from Fig.  4, Fig. 5 that both TDSE and our model exhibit a new type of correlations-the vertical and horizontal stripes. In the classical picture to which we alluded above, which associated the diagonal stripes in the pattern of correlations dominant for low fields with the correlations due to the periodic motion, this new type of correlations can be accounted for as follows. Correlations between coordinates of the electron at the moments of time t 1 and t 1 + kT a persist till the moment of time when either t 1 or t 1 + kT a gets equal to the time of ionization. After the occurrence of the ionization event, different types of correlations are introduced. In the Eq. (23) these types of correlations are described by the second and third lines of the equation. The second line describes the correlations between the electrons's coordinate along the fieldfree atomic trajectory and the electron's coordinate along the ionized trajectory. These correlations are responsible for the horizontal and vertical stripes in the correlations pattern in Fig. 4, Fig. 5. Employing the simplified classical picture we used above, we might say that this type of correlations arises because the z-coordinate of the ionized electron's wave-packet is correlated with the electron's z−coordinate either at the moment of the electron's birth t = t 1 , or the electron's coordinates at the moments of time t = t 1 − kT a with positive integer k, when electron is located at approximately the same spatial point in the course of its orbital motion. Similarly, the minima in the correlations pattern appear because the z-coordinate of the ionized electron's wave-packet is anticorrelated with the electron's z−coordinate at the moments of time t = t 1 − kT a /2 with positive integer k, when electron is located at approximately opposite spatial point of its orbit. The area of highly correlated motion seen for the real part of the correlation function in Fig. 4 in the region t 1 /T > 1/2, t 2 /T > 1/2 is due to the term in the third line of the Eq. (23), which describes essentially the propagation of correlations for the free electron motion, which grows with time as t 1 t 2 for large times. We see that all these features are present in both the ab initio TDSE and our model calculation based on the Eq. (23). Our model calculation reproduces these features quite accurately qualitatively, and even, as the Fig. 4 shows, quantitatively in the case of the real part of the autocorrelation function. Agreement between the imaginary parts C zz (t 1 , t 2 ) given by the ab initio TDSE and the model calculation based on Eq. (23) is less spectacular but, we believe, can still be considered pretty good given that we use only one adjustable parameter in our model calculations. That our model reproduces the real part of the correlation function better is, perhaps, not surprising, taking into account that in applying our Eq. (23) we employed only one real positive function α(t) proportional to the ionization probability, discarding thereby all information about phase. Going back to the real part of the autocorrelation function, we would like to draw attention to the fact that the transition be-tween the two types of correlations which we discussed above: the field-free correlations described by the term in the first line of the Eq. (23), and the correlations between the electrons's coordinate along the field-free atomic trajectory and the electron's coordinate along the ionized trajectory described by the terms in the second line of the Eq. (23), is quite distinct. Given that correlation function is, in principle, an experimentally measurable quantity [12], this offers an intriguing possibility of providing an experimental access to study the somewhat elusive notions such as the moment of the electron's birth into the continuum and the tunneling time, which have received considerable attention in the literature lately [13][14][15][16][17][18][19][20][21][22].
Having in our disposal the coordinate autocorrelation function C zz (t 1 , t 2 ), we can find other correlation functions. An example is shown in Fig. 6, Fig. 7 where we present the velocity autocorrelation function C vzvz (t 1 , t 2 ) obtained from C zz (t 1 , t 2 ) as: The velocity correlation pattern predicted by our model based on the Eq. (24), Eq. (23) agree qualitatively well with the ab initio TDSE results.
To summarize, we presented a simple tentative solution of the Heisenberg operator equations of motion for an atom exposed to electromagnetic field. We might call this solution a simple man model in the Heisenberg picture since our main equation Eq. (9) looks very much like its classical counterpart used in the SMM. We tested the veracity of this tentative expression by applying it to the calculation of the coordinate and velocity autocorrelation functions. As the ab initio TDSE calculations show, these functions are an interesting object of study and, to our knowledge, have not been studied before in the context of the strong field ionization. In particular, their exhibiting the distinct types of correlations due to different types of electron's motion may offer a useful insight into the physics of strong field ionization. We saw that the autocorrelation functions computed by fairly simple means using the main equation Eq. (9) agree well with the results of the TDSE calculations, thereby confirming the utility of our model.      (20). a Results of the TDSE calculation (Eq. (21)). b Results obtained using Eq. (23) with α(t) computed as ionization probability using Yudin-Ivanov ionization rate (YU IIR). For better visibility and to reveal more detail we plot the quantity (Re(Czz(t1, t2)) 1/5 .  (20). a Results of the TDSE calculation (Eq. (21)). b Results obtained using Eq. (23) with α(t) computed as ionization probability using Yudin-Ivanov ionization rate (YU IIR). The quantity (Im(Czz(t1, t2)) 1/5 is plotted.  21)). b Results given by the Eq. (23) with α(t) computed as ionization probability using Yudin-Ivanov ionization rate (YU IIR).