# Many-electron tunneling in atoms

## 1 Introduction

Many-electron ionization of atoms by laser field was firs observed by Suran and Zapesochny  in alkali-earth atoms (the review of that work as well as some earlier ones see in ). At present, such studies form one of main guidelines in physics of strong field interaction with atoms .

A number of theoretical models were proposed for interpretation of the gathered experimental data. Some models dealt with direct influence of laser radiation on the atomic electrons , the others consider highly-stripped ion formation due to nonelastic scattering of previously emitted electrons with the parent ion . These models allow to explain a number of observed features of the phenomenon . Nevertheless, there are some difficulties in theoretical description of highly-stripped ion formation in laser field which is not related to nonelastic collisions . Due to these difficulties, the above mechanisms cannot be properly used for explanation of the experiment.

At the same time, it is well-known fact that the single-charged ion formation by a laser field in tunnelling regime can be satisfactory described in terms of relatively simple formulae of the ADK theory . An empirical generalization of the ADK formulae for describing the highly-stripped ion formation was proposed in . So it would be reasonable to generalize the available theory of tunnelling in atoms to the case of non-sequentional multiple ionization of atom. Solution of this problem is the objective of the present work.

Obviously, the Josephson effect can be considered as a solid-state analogue of the considered phenomenon. Some considerations on difference between the one- and many-particle tunnelling are mentioned in reference . Comparison of these considerations with the results of the present work shows that the mentioned difference for tunnelling in atoms is not so trivial as it was described in .

## 2 Asymptotics of the many-electron wave function

Let us remembers some facts which make the main proposed concepts easer to understand. To describe optical transitions in complex atoms, Bates and Damgaard  modified the Slater method . Basically, the nodeless character of Slater orbitals was retained. Unlike the Slater method, the effective nuclear charge ceases to be a fitting parameter for valence electrons in atom since it coincide with the residual ion charge. But the effective principal quantum number is uniquely determined by the electron coupling energy. So the asymptotical region of electron motion is considered, where the atomic potential has Coulomb shape. High accuracy of oscillator strengths calculations  using the Bates–Damgaard method and its clear physical justification allows to use this method in calculations of other atomic characteristics determined by large electron-nucleus distances.

The tunnelling probability is also determined by large electron-nucleus distances where the energy of the electron interaction with the external field becomes comparable with the attractive energy of the residual ion. So the Bates–Damgaard method can be used for describing the tunnelling effect. Such a procedure was developed in recent work  for tunnelling calculation in Rydberg molecules. In that work some evaluations are presented for the applicability conditions of the method.

Let $`N`$ equivalent (i. e. belonging to the same atomic shell) electrons are removed from the atom via tunnelling. Then the asymptotic behaviour of the radial part of $`N`$-electron wavefunction in the Bates–Damgaard approximation is determined by the product of properly symmetrized one-electron function asymptotics:

$`\psi _{\nu lm}(𝒓)C_{\nu l}b^{3/2}\left({\displaystyle \frac{r}{b}}\right)^{\nu 1}\mathrm{exp}\left({\displaystyle \frac{r}{b}}\right)Y_{lm}\left({\displaystyle \frac{𝒓}{r}}\right),`$ (1)   
$`C_{\nu l}=(2\pi \nu )^{1/2}\left({\displaystyle \frac{2}{\nu }}\right)^\nu L(\epsilon ),L(\epsilon )=\left({\displaystyle \frac{1\epsilon }{1+\epsilon }}\right)^{\frac{1}{2}(l+1/2)}(1\epsilon ^2)^{\nu /2}.`$   
Here $`b=a\nu /Z`$, $`Z`$ is the residual ion charge, $`a=\mathrm{}^2/\mu e^2`$ is Bohr radius, $`\mu `$, $`e`$ are the mass of electron and the absolute value of its charge, $`\epsilon =(l+\frac{1}{2})/\nu `$. The $`C_{\nu l}`$ constant in (1) is determined in quasiclassical approximation not implying the condition $`l\nu `$, which was required in . It results in the arising of $`L(\epsilon )`$ function with $`L(\epsilon )1`$ at $`\epsilon 0`$. After passage to this limit the expression (1) for the $`C_{\nu l}`$ constant turns into the formula (11) of the reference  (with an inaccuracy corrected: the number $`\mathrm{e}=2.718\mathrm{}`$ should be omitted).

The expression (1) for $`C_{\nu l}`$ is obtained under $`\epsilon <1`$. For $`\epsilon >1`$, the quasiclassical approximation is not valid, so calculation of $`C_{\nu l}`$ requires numerical approaches (see, e. g., ).

The principal quantum number $`\nu `$ is determined by the electron coupling energy. Denoting the first, second etc. ionization potentials of the atom as $`E_1/e`$, $`E_2/e\mathrm{}`$, the principal quantum number of $`j`$-th removed electron is

$$\nu _j=\left(\frac{2aE_j}{Z^2e^2}\right)^{1/2}.$$

If the electron are equivalent and are *simultaneously* removed from the atom, then for all the electrons

$$\nu =\left(\frac{2aE_N}{NZ^2e^2}\right)^{1/2},$$
(2)   
where

$$E_N=\underset{j=1}{\overset{N}{}}E_j$$

is the coupling energy of $`N`$ electrons. Note that in framework of the considered model, the asymptotic behaviour of the bound electron wave function (1) *depends* on the number of the removed electrons. So a partial account is provided for many-electron effects in the initial state.

Now we consider $`N`$-electron ionization as removal of a $`N`$-electron “bundle” – a peculiar kind of quasiparticle of mass $`N\mu `$ and of charge $`Ne`$. In the region which determines the ionization process, we consider the distances between the electrons in the bundle to be much less than the separation between the atomic core and the center of bundle mass. Denoting the distance between the $`i`$-th and $`j`$-th electrons as $`𝒙_{ij}`$, and the position of the center of bundle mass as $`𝑹`$, we write the corresponding inequality:

$$x_{ij}R.$$
(3)   
Since the atom–laser radiation interaction is considered in dipole approximation, the influence of the field on $`N`$ individual electrons is completely equivalent to the influence of the field on a quasiparticle of charge $`Ne`$ which is located at the point $`𝑹`$. As for the interaction of this quasiparticle with the core Coulomb field, the correspondent error value is $`(x_{ij}/R)^2`$, which is small due to the accepted inequality (3).

For the mathematical description of the considered model, one should solve a problem which is analogous to that is occurred, e. g. in nuclear $`\alpha `$-decay theory. This problem is to construct the quasiparticle wave function $`\mathrm{\Psi }_{\{\nu lm\}}^{(N)}(𝑹,\{𝒙_i\})`$ at large distances from the residual system, using the one-particle wave functions of the system in the initial state. Symbols in the braces are sets of quantum numbers or coordinates of individual particles. To solve this problem we consider the asymptotics of the function $`\mathrm{\Psi }_{\{\nu lm\}}^{(N)}`$ at $`R\mathrm{}`$, which is a product of the one-electron function asymptotics (1). It is easy to see that the radial dependencies of the functions (1) bring the factor

$$\mathrm{exp}\left(\frac{NR}{b}\right)\left(\frac{R}{b}\right)^{N(\nu 1)}.$$

into the asymptotics of $`\mathrm{\Psi }_{\{\nu lm\}}^{(N)}`$. To obtain the angular dependence, the mean of the $`𝑹,\{𝒙_i\}`$ variables should be detalized. Since the problem has the axial symmetry for the linearly polarized field, the orbital moment projections of non-interacting electrons onto the polarization direction are conserved. So it is convenient to leave the azimuth angles $`\phi _i`$ the same that in the original spheric coordinate system centered in the atomic nucleus. The change of variables will effect only on the absolute values $`\{r_i\}`$ and polar angles $`\{\theta _i\}`$. At $`\theta 0`$, the behaviour of the Legendre polynomials involved in the spheric functions (1), is determined by

$$P_l^{|m|}(\mathrm{cos}\theta _i)(1)^{|m|}\frac{\mathrm{sin}^{|m|}\theta _i}{2^{|m|}|m|!}=(1)^{|m|}\frac{(r_i^2r_{iz}^2)^{|m|/2}}{2^{|m|}|m|!r_i^{|m|}}$$

Substituting here $`r_iR`$, $`r_{iz}R_z`$ and introducing the parabolic coordinates $`\xi =R+R_z`$, $`\eta =RR_z`$ for the center of the bundle mass, the asymptotics of the $`N`$-electron function at $`\xi \eta `$ can be written in the form

$`\mathrm{\Psi }_{\{\nu lm\}}^{(N)}(𝑹,\{𝒙_i\})=B\varphi (\xi ,\eta )\chi (\{r_i,\theta _i\}){\displaystyle \underset{j=1}{\overset{N}{}}}{\displaystyle \frac{1}{\sqrt{2\pi }}}\mathrm{exp}(\mathrm{i}m_j\phi _j),`$   
$`B=a^{3/2}C_{\nu l}^N\left({\displaystyle \frac{Z}{\nu }}\right)^{3N/2}(2l+1)^N{\displaystyle \underset{j=1}{\overset{N}{}}}{\displaystyle \frac{(1)^{|m_j|}}{|m_j|!}}\left[{\displaystyle \frac{(l+|m_j|)!}{(l|m_j|)!}}\right]^{1/2},`$ (4)   
$`\varphi (\xi ,\eta )\mathrm{exp}\left[{\displaystyle \frac{N(\xi +\eta )}{2b}}\right]\left({\displaystyle \frac{\xi }{2b}}\right)^{N(\nu 1)}\left({\displaystyle \frac{\eta }{\xi }}\right)^{M/2},M={\displaystyle \underset{j=1}{\overset{N}{}}}|m_j|.`$   
Here $`\chi `$ is the normalized per unit wave function of the electron inner motion in the bundle. Note that there are only $`2(N1)`$ independent variables $`\{r_i,\theta _i\}`$ of $`2N`$. The function $`\varphi (\xi ,\eta )`$ describes the motion of the center of the bundle mass.

## 3 Tunnelling probability

The further calculation of the tunnelling probability is implemented according the standard technique , an account provided for that the electron bundle mass is $`N\mu `$ and its charge is $`Ne`$. Substituting the function $`\varphi (\xi ,\eta )`$ from (4) into the Schrödinger equation

$$\frac{\mathrm{d}}{\mathrm{d}\xi }\left(\xi \frac{\mathrm{d}\varphi }{\mathrm{d}\xi }\right)+\left(\beta \frac{E_NN\mu }{2\mathrm{}^2}\xi \right)\varphi =0,$$

describing the motion with respect to the parabolic $`\xi `$ coordinate at $`\xi \mathrm{}`$, we obtain the variables separation constant:

$$\beta =\frac{N}{b}\left[N(\nu 1)\frac{M1}{2}\right].$$
(5)   
The centrifugal potential is neglected since it vanishes rapidly at $`\xi \mathrm{}`$.

Now we consider the external field $`F(t)`$ to be slow-varying, and use quasiclassical approximation for the wave function $`\varphi _F(\xi ,\eta )`$ which describes the center of the bundle mass motion in the field. In the below-threshold domain

$`\varphi _F(\xi ,\eta )=\varkappa (\xi |p(\xi )|/\mathrm{})^{1/2}\mathrm{exp}\left({\displaystyle \frac{1}{\mathrm{}}}{\displaystyle _{\xi _1}^\xi }|p(\xi )|d\xi \right),`$ (6)   
$`p(\xi )=\mathrm{}\left({\displaystyle \frac{E_NN\mu }{2\mathrm{}^2}}+{\displaystyle \frac{\beta }{\xi }}+{\displaystyle \frac{1}{4\xi ^2}}+{\displaystyle \frac{N^2e\mu }{4\mathrm{}^2}}F\xi \right)^{1/2},`$ (7)   
where $`\xi _1`$ is the greater root of the equation $`p(\xi )=0`$. Comparing the expression (6) with the function $`\varphi (\xi ,\eta )`$ from (4) at the point $`\xi _0`$ lying in the region

$$\frac{2\mathrm{}^2\beta }{E_NN\mu }b\xi _0\frac{2E_N}{NeF}=\frac{eZ}{b\nu F},$$
(8)   
we obtain the $`\varkappa `$ value:

$$\varkappa (\eta ;\xi _0)\left(\frac{N\xi _0}{2b}\right)^{1/2}\mathrm{exp}\left(\frac{1}{\mathrm{}}_{\xi _0}^{\xi _1}|p(\xi )|d\xi \right)\varphi (\xi _0,\eta ).$$
(9)   
The condition of existence of the region (8) leads to following restriction to the external field:

$$FF_a\frac{eZ}{b^2\nu }=\frac{e}{a^2}\left(\frac{Z}{\nu }\right)^3,$$
(10)   
which differs from the condition arising in the one-electron tunnelling description only by the definition of the $`\nu `$ value. It should be noted that for $`\nu `$ essentially greater than 1 (what holds, e. g. for Rydberg states) the inequality (10) is changed by a stronger one:

$$F<\frac{Z^3e}{16\nu ^4a^2},$$
(11)   
which is deduced from the condition of existence of the potential barrier .

The formulae (6) and (9) determine the function $`\varphi _F(\xi ,\eta )`$ outside the barrier. With the account of inequality (8), its squared absolute value is 

$`|\varphi _F(\xi ,\eta )|^2`$ $`=`$ $`{\displaystyle \frac{\mathrm{}N\xi _0}{2b\xi p(\xi )}}\left({\displaystyle \frac{\xi _0}{2b}}\right)^{2N(\nu 1)}\left({\displaystyle \frac{\eta }{\xi _0}}\right)^M`$   
$`\times `$ $`\mathrm{exp}\left[{\displaystyle \frac{N\eta }{b}}{\displaystyle \frac{16\mathrm{}^2}{3N^2\mu eF}}\left({\displaystyle \frac{E_NN\mu }{2\mathrm{}^2}}\right)^{3/2}\beta \left({\displaystyle \frac{2\mathrm{}^2}{E_NN\mu }}\right)^{1/2}\mathrm{log}{\displaystyle \frac{NeF\xi _0}{8E_N}}\right].`$   
Using (2) and (5), it is easy to see that the dependence on the arbitrary parameter $`\xi _0`$ is actually disappeared in (3):

$$|\varphi _F(\xi ,\eta )|^2=\frac{\mathrm{}N(\eta /b)^M}{2^M\xi p(\xi )}\left(\frac{2F_a}{F}\right)^{2N(\nu 1)M+1}\mathrm{exp}\left(\frac{N\eta }{b}\frac{2NF_a}{3F}\right).$$
(13)   
The ionization probability is determined by the flux of probability density (13) through a plane perpendicular to $`z`$-axis :

$$W_{\nu l}^{(N)}(F)2\pi _0^{\mathrm{}}v_z|\varphi _F(\xi ,\eta )|^2\rho d\rho ,v_z=\frac{2p(\xi )}{N\mu },\rho =\sqrt{\xi \eta },\mathrm{d}\rho \sqrt{\frac{\xi }{\eta }}\mathrm{d}\eta .$$

Substituting here the formulae (4) and (13), we obtain: 

$`W_{\nu l}^{(N)}(F)`$ $`=`$ $`{\displaystyle \frac{\pi \mathrm{}}{a^2\mu }}{\displaystyle \frac{M!(2l+1)^NC_{\nu l}^{2N}}{2^{M2}N^{M+1}}}\left({\displaystyle \frac{Z}{\nu }}\right)^{3N1}{\displaystyle \underset{j=1}{\overset{N}{}}}{\displaystyle \frac{(l+|m_j|)!}{(|m_j|!)^2(l|m_j|)!}}`$ (14)   
$`\times `$ $`\left({\displaystyle \frac{2F_a}{F}}\right)^{2N(\nu 1)M+1}\mathrm{exp}\left({\displaystyle \frac{2NF_a}{3F}}\right).`$   
This formula determines the $`N`$-electron tunnelling probability in dc field within a factor accounting for the overlapping of wave functions of the electrons remaining in the atom, with the wavefunctions of the same electrons in the initial state. Obviously, this factor cannot exceed 1, and its more accurate evaluation can be performed only numerically. Note that the $`N`$ multiplier in the exponent in (14) in now ways gives an exhaustive account for the dependence of this exponent on $`N`$, as it was considered in . Due to the formulae (2), (10), this dependence is significantly more complicated and it is determined by the spectrum of the particular atom. We present below (figure 1) a numerical example illustrating this statement.

Now we consider that

$$F(t)=F_0\mathrm{cos}\omega t,$$
(15)   
where $`\omega `$ is the laser field frequency. It is a well-known fact that the tunnelling in a laser field is possible for small values of the Keldysh parameter 

$$\gamma =\frac{\sqrt{2\mu E_1}}{eF}\omega ,$$

where $`E_1`$ is the coupling energy of one electron. Following the technique developed in  for “particle” of mass $`N\mu `$ and charge $`Ne`$, it is easy to see that the $`N`$-electron tunnelling is possible for small values of the parameter

$$\gamma _N=\frac{\sqrt{2\mu E_N/N}}{eF}\omega .$$
(16)   
Since the coupling energy is increasing for each subsequent electron, $`N`$-electron tunnelling requires field values lower than $`N`$-electron tunnelling cascade.

Substituting (15) into (14), we average the result over the time interval $`t[\pi /2\omega ,\pi /2\omega ]`$ <sup>1</sup><sup>1</sup>1The values $`t[\pi /2\omega ,3\pi /2\omega ]`$ leads to $`F(t)<0`$ and the tunnelling takes place in the direction of negative $`z`$ semiaxis.. Due to the inequality (10), the integral arising here can be calculated using the saddle-point method. Under the condition (11) fulfilled, the saddle point is $`t=0`$, and the final formula is:

$`W_{\nu l}^{(N)}(F_0)`$ $`=`$ $`{\displaystyle \frac{\sqrt{3\pi }\mathrm{}}{a^2\mu }}{\displaystyle \frac{M!(2l+1)^NC_{\nu l}^{2N}}{2^{M3/2}N^{M+3/2}}}\left({\displaystyle \frac{Z}{\nu }}\right)^{3N1}{\displaystyle \underset{j=1}{\overset{N}{}}}{\displaystyle \frac{(l+|m_j|)!}{(|m_j|!)^2(l|m_j|)!}}`$ (17)   
$`\times `$ $`\left({\displaystyle \frac{2F_a}{F_0}}\right)^{2N(\nu 1)M+1/2}\mathrm{exp}\left({\displaystyle \frac{2NF_a}{3F_0}}\right).`$   
Remember that the exponent dependence on $`N`$ in (17) is not reduced to the factor $`N`$ which is written explicitly.

## 4 Numerical examples

Unfortunately, the obtained formulae cannot be immediately related to an experiment, because, along with the direct $`N`$-fold ions formation, there are a number of cascade processes as well as other ionization mechanisms due to nonelastic collisions of electrons and ions . For the relation of the theory with an experiment, the correspondent kinetic equations are to be solved, that should be a subject for another work. So only some illustrative examples are considered in this section.

The figure 1 presents the relation of probabilities of 3-fold ions formation in the noble gases resulted from two 2-cascade processes: $`\mathrm{A}\mathrm{A}^+\mathrm{A}^{3+}`$ and $`\mathrm{A}\mathrm{A}^{2+}\mathrm{A}^{3+}`$. These probabilities are denoted as $`W(1;2)`$ and $`W(2;1)`$ correspondingly. They have similar dependence on the laser pulse duration. As it is seen, the relation $`W(1;2)/W(2;1)`$ is not equal to 1, as it is follows from the results of reference .

The following result seems to be curious. The 2-electron tunnelling probabilities for neutral atoms can be greater than the one-electron tunnelling probabilities in correspondent singly charged ions. E. g., for Ar atom the 2-electron tunnelling probability exceeds the 1-electron process probability for $`\mathrm{Ar}^+`$ ion at the intensities $`I>10^{14.88}\mathrm{W}/\mathrm{cm}^2`$. The same result takes place for Kr at $`I>10^{14.76}\mathrm{W}/\mathrm{cm}^2`$, for Xe at $`I>10^{14.34}\mathrm{W}/\mathrm{cm}^2`$. At the same time, for light noble gases atoms He and Ne, the probabilities of one-electron tunnelling in singly charged ions are approximately by two orders greater than the probabilities of two-electron process in the correspondent neutral atoms at $`I10^{15}\mathrm{W}/\mathrm{cm}^2`$. These facts shows wide range of experimental situations arising in multiphoton tunnelling effect.

This work was stimulated by the report . The author is grateful to Professor W. Sandner for the interest to the work, and to WE–Heraeus-Stiftung for the offered opportunity to participate in the seminar work. I also express my deep gratitude to Professor N. B. Delone and to the participants of his seminar in IOF RAN for helpful discussion. This work was partially supported by Russian Foundation for Basic Researches (grant no. 97-02-18035).