# Tunneling Ionization Rates from Arbitrary Potential Wells

## Abstract

We present a practical numerical technique for calculating tunneling ionization rates from arbitrary 1-D potential wells in the presence of a linear external potential by determining the widths of the resonances in the spectral density, $`\rho (E)`$, adiabatically connected to the field-free bound states. While this technique applies to more general external potentials, we focus on the ionization of electrons from atoms and molecules by DC electric fields, as this has an important and immediate impact on the understanding of the multiphoton ionization of molecules in strong laser fields.

While the behavior of atoms and molecules in strong laser fields has been studied for almost two decades, new phenomena continue to arise and a unified picture for either does not yet exist . However, simple models for tunneling ionization  have provided a crucial baseline for understanding atoms in strong fields, as deviations from this basic process in ion yields and electron energy spectra have pointed to new physics . Conversely, the lack of a general tunneling model for molecules has seriously impeded the interpretation of molecular data which clearly hints at even richer phenomena  than for atoms due to the extra degrees of freedom present in a molecule.

Ionization from even diatomic potential wells are more difficult to treat than atoms because the barrier to ionization can shift between the barrier separating the two wells and the outer barrier as a function of field strength and internuclear separation  (see Fig. 1). These situations are quite distinct and, indeed, this competition between barriers leads to interesting and complex behavior. Thus, it is critical for any tunneling ionization model to work for arbitrary potential wells so as to treat these different situations uniformly. Furthermore, a useful model for molecules in strong laser fields must be non-perturbative in both the rates and the Stark shifts. In this Letter we introduce a practical non-perturbative technique for computing tunneling ionization rates from an arbitrary 1-D potential well. After outlining the general method, we apply it to a single-well (‘atomic’) potential, yielding excellent agreement with existing theoretical analyses. Then we apply our method to a double-well (‘molecular’) potential, like that shown in Fig. 1, to illustrate the generality of the method and to show new features which arise due to the double-well structure. However, we emphasize that the method is exactly the same for any arbitrary shaped 1-D potential well.

The problem of determining ionization rates and Stark shifts non-perturbatively for an arbitrary potential in 1-D can be reduced to calculating the spectral density, $`\rho (E)`$. For potentials with strictly bound states $`\rho (E)`$ becomes a series of delta functions at the energy eigenvalues. Numerically integrating the Schrödinger equation (SE) determines these energies by finding the energy at which the wavefunction satisfies the given boundary conditions. This simple approach cannot be extended to the question of tunneling ionization because, by definition, the electron can escape to infinity. Given this fact, the electron’s state lies in the continuum; hence, all energies are allowed and simply integrating the SE at a particular energy provides little information, on its own. Rather, one must return to calculating $`\rho (E)`$ in general. In the presence of a linear potential, the delta functions in $`\rho (E)`$ for the bound states will evolve into Lorentzian-like resonances whose center positions reflect the Stark shift of the states and widths give their decay rates. As decay of the levels can only occur through ionization the width then directly gives the tunneling rate.

The limitation to 1-D potential wells may seem overly restrictive and requires some discussion. It is well known that ionization rates of atoms and molecules depend strongly on field strength . For this reason ionization occurs overwhelmingly in the direction of the electric field, essentially reducing the problem to 1-D. For diatomic molecules, the ionization rate is much greater when the field is along the internuclear axis  and, again, reduces the problem to 1-D.

While a full discussion of $`\rho (E)`$ is rarely presented in standard books on quantum mechanics, a general approach exists for calculating $`\rho (E)`$, based on the Weyl-Titchmarsh-Kodaira (WTK) spectral theorem . Remarkably, this approach is extremely well suited to numerical implementation. Essentially, the spectral theorem states that the “resolution of the identity” can be expressed as

$$\delta \left(xy\right)=\underset{i,j=1,2}{}_{\mathrm{}}^{\mathrm{}}𝑑E\varphi _i(x,E)\varphi _j(y,E)\rho _{ij}\left(E\right)$$
(1)   
where $`\varphi _{1,2}(x,E)`$ are two linearly independent eigenfunctions of the SE at energy $`E`$ and $`\rho _{ij}(E)`$ is the spectral density matrix. Before discussing how to calculate and interpret $`\rho _{ij}(E)`$, we must consider the question of normalization, as, in general, some of the states may lie in the continuum.

The usual concept of a finite square integral clearly cannot be used for the normalization of continuum wavefunctions. In the WTK approach, normalization is imposed locally by setting the Wronskian of the eigenfunctions equal to one which also ensures that $`\varphi _1`$ and $`\varphi _2`$ are linearly independent. Thus, we can choose the following initial conditions at any point $`x=c`$:

$$\begin{array}{cccc}\varphi _1\left(c\right)=1,& \varphi _1^{}\left(c\right)=0,& \varphi _2\left(c\right)=0,& \varphi _2^{}\left(c\right)=1.\end{array}$$
(2)   
For convenience, we generally choose $`x=0`$. The complete functions $`\varphi _i(x)`$ can now be found by numerically integrating the SE away from the origin using these initial conditions.

To construct the spectral density matrix $`\rho _{ij}(E)`$, we first find the linear combinations of $`\varphi _1`$ and $`\varphi _2`$,

$$\psi _\pm (x,\lambda )=\varphi _1(x,\lambda )+m_\pm \left(\lambda \right)\varphi _2,$$
(3)   
such that $`_{\mathrm{}}^0𝑑x\left|\psi _{}\right|^2`$ and $`_0^{\mathrm{}}𝑑x\left|\psi _+\right|^2`$ are both finite, where $`\lambda =E+i\epsilon `$ and $`\epsilon 0`$. In other words, when the energy has a positive imaginary part, $`\psi _\pm (x,\lambda )`$ must decay sufficiently rapidly as $`x\pm \mathrm{}`$ to have a finite square integral on their respective half-line. This defines the functions $`m_\pm \left(\lambda \right)`$. The spectral density matrix is then :

$$\rho _{ij}(E)=\frac{1}{\pi }\underset{\epsilon 0}{lim}\mathrm{Im}\frac{1}{m_{}m_+}\left[\begin{array}{cc}1& \frac{m_{}+m_+}{2}\\ \frac{m_{}+m_+}{2}& m_{}m_+\end{array}\right].$$
(4)   
If the asymptotic wavefunctions of the external potential are known analytically, as is often the case, all that is needed to evaluate $`m_\pm `$ is the logarithmic derivative of the correctly decaying asymptotic wavefunctions of just the external potential, $`\theta _\pm (x,\lambda )`$: $`f_\pm (x,\lambda )=\theta _\pm ^{}(x,\lambda )/\theta _\pm (x,\lambda )`$ . We can then solve for $`m_\pm \left(\lambda \right)`$:

$$m_\pm \left(\lambda \right)=\left(\frac{\varphi _1^{}f_\pm \varphi _1}{\varphi _2^{}f_\pm \varphi _2}\right).$$
(5)   
Note that the overall normalization of the asymptotic wavefunctions is not important, as the logarithmic derivative is sufficient to determine $`m_\pm \left(\lambda \right)`$.

This entire procedure is very straightforward to implement numerically: starting with the boundary conditions in Eq. (2), $`\varphi _1`$ and $`\varphi _2`$ are found by integrating the SE away from the origin with a standard Runge-Kutta technique at a particular complex energy, $`\lambda `$. At the same time, $`m_\pm `$ are evaluated from Eq. (5). Ideally, one would integrate out to $`\pm \mathrm{}`$ to evaluate $`m_\pm \left(\lambda \right)`$ exactly. Since this is not possible numerically, we simply evaluate $`m_\pm \left(\lambda \right)`$ as a function of $`x`$ far enough from the local potential to where it converges to a constant value. At this point, $`m_\pm \left(\lambda \right)`$ have been determined as well as $`\rho _{ij}(E)`$ from Eq. (4).

To finish the problem, we must interpret the matrix in Eq. (4). The spectral density matrix, $`\rho _{ij}(E)`$, is positive definite  and, thus, its eigenvalues $`\mathrm{\Lambda }_1`$ and $`\mathrm{\Lambda }_2`$ are positive and represent scalar spectral density functions. In general, there are two functions because of the possibility of a two-fold degeneracy (corresponding, for example, to left and right moving solutions) on the whole line. However, in the case of a linear external potential (Fig. 1), one direction is inaccessible to a free particle and no degeneracy exists. In this case, either $`\mathrm{\Lambda }_1`$ or $`\mathrm{\Lambda }_2`$ will be on the order of $`\epsilon `$ and is unimportant, leaving just one eigenvalue, resulting in a simple scalar spectral density. For problems on the half-line, $`x0`$, with a physical boundary condition at $`x=0`$, the spectral density matrix also reduces to a single spectral density function . Finally, we must consider taking the limit $`\epsilon 0`$. Again, numerically, this is simple to specify: if the state can tunnel, it will have a “natural” width $`\mathrm{\Gamma }`$ giving the tunneling rate and an “artificial” width determined by $`\epsilon `$. One then decreases $`\epsilon `$ until $`\mathrm{\Gamma }\epsilon `$. We have noticed that the numerical results are surprisingly robust to changes in the various numerical parameters. We attribute this to the local specification of the normalization and initial conditions for the numerical integration and the ease in determining the $`m_\pm `$ values.

In order to test this method we considered tunneling ionization of a hydrogen atom by a DC electric field:

$$V(x)=\frac{1}{\sqrt{x^2+a^2}}Fx,$$
(6)   
where $`F`$ is the field strength and $`a=\sqrt{2}`$ for a hydrogen atom. This 1-D form of the Coulomb potential  has proven to be an excellent model for strong field calculations because it preserves the long range $`1/r`$ dependence and has an infinite Rydberg-like series of bound states while removing the catastrophic singularity that occurs for 1-D. Furthermore, it has been used recently in direct calculations of the time-dependent Schrodinger equation for molecules . Fig. 2 shows the spectral function at a moderate field strength. The first two low lying states in the Rydberg series can be seen, while the higher states have broadened and merged together. We stress that the blurring of the higher states is a real physical effect, due to the lowering of the potential barrier by the external field, and is not a limitation of the method. In a weaker field more of the excited states would be seen as well-defined peaks. Fig. 3 shows a typical unbound wavefunction at the peak of a resonance. At this point, most of the wavefunction is ”inside” the potential well, although it is still a continuum state. Fig. 4 shows the DC Stark shift as well as the first term in the second order perturbation theory expansion for the 1-D potential. Finally, Fig. 5 shows the tunneling ionization rate of the ground state as a function of field strength. In order to compare to a 3-D potential, we performed an angular average of the 1-D result:

$$\mathrm{\Gamma }(F)=\frac{1}{2\pi }𝑑\mathrm{\Omega }\mathrm{\Gamma }\left(F\mathrm{cos}\theta \right).$$
(7)   
As seen in Fig. 5, this averaged result agrees quite well with the standard 3-D static atomic tunneling model of ADK . However, the slight discrepancy at the higher intensity is intriguing. The ADK model depends on the ionization potential and the field strength, but does not take into account the DC Stark shift. While the shifts are small, ionization rates depend exponentially on ionization potential and the effect of the Stark shift on the ionization rate will increase with intensity. We recalculated the ADK rate using our calculation of the DC Stark shift, and, as seen in Fig. 5 the agreement between our angularly averaged model and the ADK result including the DC Stark shift is excellent.

The ’molecular’ or double-well potential has a much richer spectrum in an external field than an atom, and we cannot give a comprehensive treatment, here. However, as an example, we show, in Fig. 6, the effect of the separation of the double-well potential on ionization rates for both the ground and first excited states. The potential consists of two 1-D Coulomb potentials (Eq. 6) separated by a distance R, and the external field, as in Fig. 1. Using exactly the same method of calculating the spectral function and finding the widths of the peaks associated with the ground and first excited states, we find that the ionization rate is strongly dependent on the well separation, becoming greater at larger separations. Moreover, the qualitative shape of the ionization rate vs. intensity curve depends on the separation with apparently no saturation setting in for the excited state at the larger separation. In fact, these are just the conditions for electron localization and enhanced ionization predicted in Ref. . These ionization rate curves can provide a framework to analyze experimental ion yield data from molecules.

In summary, we have described a practical numerical technique for calculating spectral density functions for arbitrary 1-D potentials on the whole line from which we can extract level shifts and tunneling rates. This will be of great value for interpreting data from atoms and molecules in strong laser fields.

We thank S. A. Fulling and W. P. Reinhardt for helpful discussions. We would like to acknowledge support from the NSF under grant PHY-9502935. G. N. G. was also supported through funding as a Cottrell Scholar of Research Corporation. K. J. B. was supported by the Research Experience for Undergraduates (REU) program of the NSF. G. D. acknowledges support from the DOE under grant DE-FG02-92ER40716.00.

Figure Captions

1. Molecular potential wells showing the importance of (a) the outer barrier and (b) the inner barrier.
2. Spectral density function for a model atom in an external field.
3. Potential in an external field with wavefunction on resonance.
4. Comparison of the 2nd order Stark shift calculation with perturbation theory, in atomic units.
5. Comparison of the angularly-averaged WTK tunneling rate model with a 3-D calculation (ADK) from Ref. .
6. Ionization rate of the ground and first excited state of a double-well potential for two different separations.