Helium atom in a box: I. doubly excited levels within the S-wave model

Darío M. Mitnik

Departmento de Física, FCEyN, Universidad de Buenos Aires, (C1428EGA) Buenos Aires, Argentina
and
Instituto de Astronomía y Física del Espacio (IAFE), Casilla de Correo 67, Sucursal 28, (C1428EGA) Buenos Aires, Argentina.

Abstract:

A complete non-perturbative close-coupling solution of the Helium atom in a box problem is presented by developing two numerical techniques. The first technique is the direct solution by diagonalization of the Hamiltonian, and the second is based on a constrained relaxation of the wave functions. A Feshbach projection-operator of the direct solutions to the bound-continuum subspaces allows a comparison of the low-lying autoionization levels obtained in both methods. Time-dependent propagation of these doubly excited wave functions is analyzed, allowing the calculation and the visualization of the autoionization process. In this work, results are presented for the S-wave model in which the electrons are restricted to spherical states and all angular correlations have been eliminated.

32.80.Dz, 31.70.Hq

Published in Physical Review A 70, 022703 (2004)


1. Introduction

Spectacular increases in computer power now provide opportunities to obtain numerical solutions to the three-body problem. A profound understanding of atomic few-body problems such as interference effects associated with resonances, correlation interactions between the charged particles, rearrangement processes, and others can be obtained by using fully quantal nonperturbative methods.

Among the fully quantum nonperturbative theories, the time-dependent close-coupling (TDCC) method has been successfully employed for calculations of electron-impact ionization  [1,2,3]. We have recently [4] studied the dielectronic capture into doubly excited resonances within the time-dependent framework. This study required the development of methods for generating accurate wave functions for doubly excited autoionizing states and time-dependent close-coupling calculations of the capture and the subsequent decay of autoionizing states. Since in future applications we wish to consider resonances in the time-dependent framework and to study processes that involve transitions to and among continuum states, we chose to discretize the wave functions and the action of operators which result from these procedures, working therefore in a numerical lattice, thus solving the problem of an ion inside a box.

Solving real atomic-physical problems with discrete numerical methods is an approximation that relies on strong logical grounds. It becomes a very good approximation to the physical bound states for all orbitals that fit well into the lattice (i.e., for such bound states in which the wave function at the boundaries approaches zero in a practical sense for numerical calculations). Use of discrete lattice functions allows the system to be described by the dynamics of always-square-integrable functions. Even if the states can be represented by fully analytical functions (such as the hypergeometrical functions for the continuum Coulomb waves), a proper calculation with these functions could involve a numerical work, and therefore, finite-lattice discretization. The method presented here has many other important advantages. It gives the exact solutions to the problem of an ion confined in a box. If the discretized allowed energy levels are positive, the eigenfunctions are true continuum solutions of the ion, not necessarily confined inside a box. These solutions form a set selected by a particular boundary condition, i.e., these are the particular true continuum wave functions with zero value at the boundaries of the box. All the solutions obtained by direct diagonalization are naturally orthogonal, and most importantly, the set is complete. Finally, the simplicity of the method is expected to lead to a greater understanding of the close-coupling formalism in general.

The problem of a spatially confined system has been a subject of interest in many branches of physics and chemistry since the early years of quantum mechanics. Nowadays, investigations on confined systems in physics (see, for example, [5]) have focused especially on the study of artificial atoms, also known as quantum dots (essentially a number of electrons confined in a potential well). However, an important theoretical motivation for the study of enclosed systems is to understand in detail the electron correlation effects on the properties of those systems.

In this work, we study how to obtain the exact wave functions in two-electron systems by the direct solution of the Schrödinger equation. Two methods are developed for this purpose. The first method consists of the direct diagonalization of the two-electron Hamiltonian on the radial grid. The second method is a constrained relaxation of the wave functions, until they relax on the successive doubly excited levels of He. The relaxation method was currently used in the TDCC method for calculation of the ground and first excited states of ions, and also for calculation of ground and low-lying excited states in Bose-Einstein condensates  [6]. However, a particular treatment is needed for calculation of doubly excited levels.

In order to explain the main features of the different methods, a detailed presentation of the theory is given for the spherically symmetric model [7,8] (also known as the Temkin-Poet model or the S-wave model). The remainder of the paper is the following. Section [3] shows the results for the first doubly excited wave functions obtained in both methods. Section [3b] shows how to use the Feshbach projection-operator formalism in order to compare the results obtained in both methods. In Sec. [3c] we propagate the first doubly excited autoionizing level, showing how the autoionization process evolves in time, and we also calculate the autoionization rate from this level by monitoring the autocorrelation function in time.

2. Theory

a. Direct solution of the Temkin-Poet He by diagonalization of the Hamiltonian

The Hamiltonian for a nonrelativistic spherically symmetric helium model may be written (in atomic units) as
\begin{displaymath} H(\vec{r_1},\vec{r_2}) = H(r_1,r_2) = -\frac{1}{2}\frac{\p... ...tial r_2^2} - \frac{Z}{r_1} -\frac{Z}{r_2} + \frac{1}{r_>} , \end{displaymath} (1)

where $r_>$ denotes the larger of the two radii $r_1$ and $r_2$. This is the simplest model for two electrons interacting with each other and with a nucleus via long-ranged Coulomb forces. In this model, both electrons $\vec{r_1}$ and $\vec{r_2}$ are restricted to spherical states, and all angular correlations are eliminated. Therefore, the full six-dimensional problem is reduced to a two-radial dimensional problem and no further distinctions between the total wave functions and the radial wave functions will be made unless explicitly stated. However, this model retains most of the other features (and computational difficulties) associated with the full He calculation. Moreover, the S-wave model is quite a good approximation to the real helium for the bound $1sns$ configurations.

We first solve this problem by a direct diagonalization of the Hamiltonian,

\begin{displaymath} H(r_1,r_2) \Psi_q(r_1,r_2) = E_q \Psi_i(r_1,r_2) , \end{displaymath} (2)

where the functions $\Psi_q$ are represented by a combination of the natural basis vectors,
\begin{displaymath} \Psi_q(r_1,r_2) = \sum_{p=1}^P c^q_p \phi_p \end{displaymath} (3)

and
$\displaystyle \phi_1 = \left(\begin{tabular}{c}1\\ 0\\ 0\\ ...\\ 0 \end{tab... ...\phi_P = \left(\begin{tabular}{c}0\\ 0\\ 0\\ ...\\ 1 \end{tabular}\right) .$     (4)

The dimension $P= n \times m$, where $n$ and $m$ are the dimensions of the $r_1$ and $r_2$ coordenates, respectively. The natural basis represents also the discretization of the physical two-dimensional space $(r_1,r_2)$ on a uniform mesh in the following way: the first $n$ vectors ($\phi_1$ to $\phi_n$) represent the points
\begin{displaymath} \left\{(r_1 = \Delta r_1 , r_2 = \Delta r_2) \rightarrow \ph... ... \Delta_{r_1} , r_2 = \Delta r_2) \rightarrow \phi_n \right\} \end{displaymath} (5)

Next, the vectors ($\phi_{n+1}$ to $\phi_{2n}$) correspond to the points:
\begin{displaymath} \left\{(r_1 = \Delta r_1 , r_2 = 2 \Delta r_2) \rightarrow \... ...ta_{r_1} , r_2 = 2 \Delta r_2) \rightarrow \phi_{2n} \right\} \end{displaymath} (6)

and so on. Therefore, the point $(r_1,r_2)=(i \Delta r_1,j \Delta r_2)$ is represented by the vector $\phi_k$, where
\begin{displaymath} k=n(j-1)+i. \end{displaymath} (7)

By discretization of the derivatives with low-order finite differences, the matrix representation of the Hamiltonian in the natural basis becomes $(m \times m)$ blocks, each one of size $(n \times n)$, having the following structure :
\begin{displaymath} {\hat H} = \left( \begin{tabular}{ccccc} $\left[ \, \begi... ... \end{tabular}\right] $\par\\ \par\end{tabular} \right) \, \end{displaymath} (8)

where the diagonal elements $\hat{H}_{kk} \equiv h_{k}$ are:

\begin{displaymath} \hat{H}_{kk} \equiv h_{k} = \frac{1}{\Delta r_1^2} + \frac... ...Z}{j \Delta r_2} + \frac{1}{max(i \Delta r_1,j \Delta r_2)} , \end{displaymath} (9)


\begin{displaymath} \alpha = -\frac{1}{2} \frac{1}{\Delta r_1^2}\hspace{1cm} , \hspace{1cm} \beta = -\frac{1}{2} \frac{1}{\Delta r_2^2} , \end{displaymath} (10)

and $k$ is related to $i$ and $j$ through Eq. (7).

The complete set of wave functions corresponding to the S-wave model He are obtained from the direct diagonalization of the matrix $\hat{H}$ [Eq. (8)]. Standard diagonalization subroutines [9] produce a matrix $C$ with rank $(n \times m)$ where the column $q$ is the $q$ eigenvector of the Hamiltonian. In this case, the value of the eigenvector $\Psi_q$ at the $(r_1,r_2)$ point in the numerical grid is given by

\begin{displaymath} \Psi_q(i \Delta r_1,j \Delta r_2)= c_{kq}, \end{displaymath} (11)

where $c$ is a matrix element of the eigenvector matrix $C$.

In order to allow the diagonalization of large matrices, we wrote the computer programs to use on distributed-memory parallel computers. The Hamiltonian matrix is directly partitioned over the many processors, so memory requirements per processor are minimized and scalability in time is achieved. The procedure used in the parallelization of the codes is similar to that employed for the parallelization of the $R$-matrix package [10].

b. The constrained-relaxation method

In this method, the energies and wave functions are calculated by relaxation of an initial wave function $\Phi$ in a fictitious imaginary time $\tau=it$  [11]. That means a transformation of the time-dependent Schrödinger equation into a diffusion equation,
\begin{displaymath} \frac{\partial \Phi(r_1,r_2,\tau)}{\partial \tau} = -H \Phi(r_1,r_2,\tau). \end{displaymath} (12)

The solution of this equation is given by:
\begin{displaymath} \Phi(r_1,r_2,\tau) = e^{-H \tau} \Phi(r_1,r_2,0) \end{displaymath} (13)

Expanding the solution in terms of the time-independent energy-eigenvector basis,
\begin{displaymath} \Phi(r_1,r_2,\tau) = \sum_{q=1}^\infty a_q \Psi_q(r_1,r_2) ... ...sum_{q=2}^\infty a_q \Psi_q(r_1,r_2)e^{-(E_q - E_1) \tau} \}, \end{displaymath} (14)

where $\Psi_1$ is the wave function of the lowest level having the same symmetry as $\Phi$, and $E_1$ is its energy, as in Eq. (2).

Since $(E_q - E_1) > 0$, the net result from this imaginary time propagation is the enhancement of those components of the solution with smaller eigenvalues of $H$ relative to those with larger eigenvalues. At the limit $\tau \rightarrow \infty$, $\Phi \rightarrow \Psi_1$. Thus, after many iterations (renormalizing continuously the wave function), only the lowest level eigenvalue (i.e., the ground state, or the first metastable level, according to the parity of the initial function $\Phi$) survives from the relaxation. Higher eigenvectors can be calculated by imposing constraints at the iteration of the relaxation which requires the state to be orthogonal, thus preventing its collapse to lower levels. This procedure is a well-known method which has been employed to determine energies of atomic and molecular systems [12], as well as the ground-state energy of a quartic oscillator [13]. In particular, for the case of He, Kulander et al. [14] integrated Hartree-Fock time-dependent equations in imaginary time, to obtain the $1s^2 ~ ^1$S ground state, and Bottcher et al. [15] calculated the energies of the $1s^2 ~ ^1$S, and the $1s2s ~ ^1$S and $1s2p~ ^1$P singly excited levels, by relaxing coupled Hylleras-type functions with a damping kinetic operator.

The relaxation method, however, will work well for the low-lying single-excited levels (e.g., $1snl$), but not for the doubly excited levels. The reason for that is that the orthogonalization procedure becomes unstable, even for a relatively low number of wave functions. In general, the relaxed wave function will collapse in any of the $1snl$ bound levels, or any of the $1skl$ continuum levels, without any warranty that this procedure will ends in a doubly excited level. doubly excited autoionizing states are calculated by imposing additional constraints at the iteration of the relaxation, projecting out also the one-electron component of the lower wave functions. Therefore, we do not allow the function to have any $1s$ component in any of the radial coordenates. The relaxation of coupled Hylleras-type functions with a damping kinetic operator procedure was used by Schultz et al. [16] for the calculation of the $2s^2$ autoionizing level of He. In this paper, the authors outlined the basic theoretical method without providing broad details concerning the computationally procedures. However, we found that this technique requires particular care in order to obtain convergent results and deserves detailed explanations of the computational algorithm employed. Based on our experience, the suggested receipt for calculating (for example) the $\Phi_{2s^2}$ wave function (corresponding to the $2s^2$ level) is the following

i)
Beforehand we calculate the $\Phi_{1s^2}(r_1,r_2)$ and $\Phi_{\gamma}(r_1,r_2)$ ($\gamma=1snl$) wave functions by the traditional relaxation method, which does not involve any other constraints than the mutual orthogonality of the wave functions.
ii)
We start the relaxation of $2s^2$ by constructing a symmetric product of two one-electron wave functions. Denoting $\phi_{2s}$ the $2s$ level of He$^{+}$ initially (i.e., at $t=0$),
\begin{displaymath} \Phi_{2s^2}(r_1,r_2,0) = \phi_{2s}(r_1)\phi_{2s}(r_2). \end{displaymath} (15)

iii)
We let the wave function $\Phi$ relax a few $n_\tau$ time steps (typically $n_\tau$=50-100, $\Delta_t=0.01$ a.u.), and then, at this time $\tau$, we project out the $1s$ state from the $r_1$ component of the wave function using
\begin{displaymath} \Phi_{2s^2}(r_1,r_2,\tau) \leftarrow \Phi_{2s^2}(r_1,r_2,\... ...t_0^\infty dr_1 ~\Phi_{1s^2(r_1,r_2)} \Phi_{1s^2}(r_1,r_2)} . \end{displaymath} (16)

This ensures that the imaginary-time-propagated wave function will not fall into any $1s(r_1)nl(r_2)$ state during the relaxation. In particular, one can easily check that
\begin{displaymath} \int_0^\infty dr_1 ~\Phi_{2s^2}(r_1,r_2,\tau) \Phi_{1s^2}(r_1,r_2) = 0 \end{displaymath} (17)

for every $r_2$.

iv)
We project out the $1s$ state from the $r_2$ component of the wave function using:
\begin{displaymath} \Phi_{2s^2}(r_1,r_2,\tau) \leftarrow \Phi_{2s^2}(r_1,r_2,\... ...t_0^\infty dr_2 ~\Phi_{1s^2}(r_1,r_2) \Phi_{1s^2}(r_1,r_2)} . \end{displaymath} (18)

Combining the last two operations is similar to performing a Feshbach projection [17] $P_q$ to the propagated wave function, where
\begin{displaymath} P_q = \vert\omega_q(r_1)\rangle\langle\omega_q(r_1)\vert + \vert\omega_q(r_2)\rangle\langle\omega_q(r_2)\vert , \end{displaymath} (19)

the $\omega_q$ are bound eigenstates of the one-electron (He$^+$) Hamiltonian, and enforcing
\begin{displaymath} P_{1s}\Phi(r_1,r_2,\tau)=0 . \end{displaymath} (20)

v)
The next step is to keep the total wave function orthogonal with all the low-lying levels previously calculated following a two-dimensional Gramm-Schmidt procedure,
\begin{displaymath} \Phi_{2s^2}(r_1,r_2,\tau) \leftarrow \Phi_{2s^2}(r_1,r_2,\t... ...y dr_1 dr_2 ~\Phi_{\gamma}(r_1,r_2) \Phi_{\gamma}(r_1,r_2)} , \end{displaymath} (21)

where the $\Phi_{\gamma}(r_1,r_2)$ are all the finally relaxed $1snl$ wave functions, already calculated with the standard relaxation procedure. The procedure outlined below should be repeated every $\tau$ time steps, until convergence is achieved (it requires a time of about 10-50 a.u. for the different wave functions).

vi)
At this point, the one-electron projection method developed here produces a wave function that is not spatially symmetric, even after applying the two-dimensional orthogonalization. Thus, at every $\tau$ step, the order of the one-electron projections [Eqs. (16) and (18)] should be alternated.

The computer codes which implement this method are also adapted to run on parallel computers. In this case, the wave functions are partitioned over the many processors in such a way that the communications between the processors are minimized and performed at every time step only for the partitioned domain borders. This parallelization scheme is a standard procedure for many of the TDCC different works (see, for example [18]).

3. Results

a. wave functions and energies

In principle, the methods outlined here are exact, and we can obtain solutions with arbitrary precision. However, our intention in the present work has not been to obtain the best energies and wave functions for the helium atom. Instead, we are interested in presenting a complete solution to the problem which could be used to understand the nature and physical significance of many-body interactions in confined atomic systems. We are interested also in the calculation of such atomic processes which are strongly dependent on these interactions. In this work, only results for the S-wave model are presented, which contains many of the electron-electron correlations, but it is simple enough to illustrate the main features and the time evolution of the solutions. Our numerical results can be improved, and we will show better close-coupling results in forthcoming work, together with the calculation of different atomic processes.

We have computed first the ground state of the S-wave model He with different numerical grids in order to check the sensitivity and convergence of the calculations and the compatibility of the two different numerical methods. Table I shows that the energy of the $1s^2$ level, obtained with a numerical grid having a mesh spacing $\Delta r = 0.2$ is $\epsilon_{1s^2}=-2.759$ a.u., compared with the best value available in the literature of $-2.87903$ a.u. (with 14 other digits that are not relevant in our comparisons) obtained by Goldman [19]. It is important to notice that for this numerical lattice, the one-electron energy of the He$^+$ $1s$ level is $\epsilon_{1s}=-1.926$ a.u., compared to the exact value of $-2.000$ au. We found an excellent agreement between the results obtained with the constrained relaxation method and with the direct diagonalization method. We performed a better calculation, with a mesh size of $\Delta r = 0.15$ a.u., where $\epsilon_{1s}=-1.9569$ a.u., and obtained a ground-state energy of $\epsilon_{1s^2}=-2.809$ a.u. with both theoretical methods. Better energies can be generated by decreasing the mesh step size and increasing the number of points. Convergence is demonstrated by using a grid with $\Delta r_1 = \Delta r_2 = 0.075$ a.u., where $\epsilon_{1s}=-1.9889$ a.u. and $\epsilon_{1s^2}=-2.861$ a.u., and with $\Delta r = 0.01$, in which $\epsilon_{1s}=-1.9998$ a.u. and $\epsilon_{1s^2}=-2.8787$ a.u. Table I also shows the energies for the first excited levels $1s2s$ and $1s3s$, for both the singlet and triplet terms. The results are in good agreement with the results given by Draeger, up to a 2% error.


Table I: Calculated energies for the first levels in the S-wave model He, in atomic units. CR means the constrained relaxation method, and DD means the direct diagonalization method.

\begin{ruledtabular} \begin{tabular}{cccccc} ~ & \multicolumn{2}{c}{$(\Delta_r =... ...-1.984 & -1.984 & -2.016 & -2.016 & -2.06079 \\ \end{tabular}\end{ruledtabular}


Figure 1 shows the probability densities of the first wave functions obtained with the direct diagonalization method, for the numerical grid with a mesh spacing $\Delta r_1 = \Delta r_2 = 0.15$ a.u., and a number of $n=m=150$ points. We show only one set of figures because the wave functions obtained from both methods, for these levels, are indistinguishable. We used the same numerical grid to calculate the energies of the He$^+$ $1s$ level and found five bound orbitals ($1s$ to $5s$), therefore we expect the plots of the figures to be similar to those obtained without a confining box.

We have computed the energy of many doubly excited states by direct diagonalization, and after the proper identification of the states, we compared them with the results obtained with the relaxation method. The numerical grid having a mesh spacing $\Delta r_1 = \Delta r_2 = 0.15$ and a number of $n=m=150$ points has 57 negative eigenenergies. Many of their eigenfunctions are bound functions of the form $1sns$, and most of them are single-continuum wave functions of the form $1sks$. We also found in this group of functions with negative energy doubly excited levels corresponding to $2s^2$, $2sns$, and $3s^2$. Table II shows the results of the calculated energies for the first doubly excited levels in the S-wave model He. Very good agreement is obtained for the $2s3s$ terms (both singlet and triplet) between the constrained relaxation and the direct diagonalization results.


Table II: Calculated energies for the first doubly excited levels in the S-wave model He, in atomic units. CR means the constrained relaxation method, and DD means the direct diagonalization method.

\begin{ruledtabular} \begin{tabular}{cccc} ~ & \multicolumn{2}{c}{$(\Delta_r = 0... ...188 \\ $3s^2$\ & -0.365 & -0.320 & -0.32142 \\ \end{tabular}\end{ruledtabular}


The probability density of these wave functions, obtained with the relaxation method, is shown in Figure  2. Among the 57 negative-energy levels obtained with the diagonalization, we can easily identify the $1snl$ and $1skl$ functions. It is also very simple to find the $2s^2$ wave function by direct inspection between the many continuum functions. This is shown in Fig. 3, where the $22^{nd}$ through $25^{th}$ wave functions (with eigenvalues $-0.844$, $-0.817$, $-0.718$, and $-0.612$ a.u., respectively) are plotted. However, for the $2s^2$ level, the functions calculated with both theoretical methods are different. In the next section, we will establish the relationship between them.

b. Feshbach projection-operator (PQ) formalism

Figure 4 shows, in a detailed scale, the first doubly excited wave function $\Phi_{2s^2}$ obtained with the constrained damped relaxation method (a), and $\Psi_{2s^2}$, obtained with the direct diagonalization method (b). It is noticeable that both functions are different because the relaxed function $\Phi_{2s^2}$ is a ``pure" bound function which does not contain any interaction with the continuum. Following Fano [21], we will observe that both functions are related by

\begin{displaymath} \Psi_{2s^2} = a \Phi_{2s^2} + \int dk ~ b_{k} \Phi_{1skl}, \end{displaymath} (22)

where $\Phi_{1skl}$ represents the ``pure" single-electron continuum functions. The coefficients $a$ and $b$ are functions of the energy. While the functions $\Psi$ are eigenfunctions of the Hamiltonian, the $\Phi$ functions are not. In this section, we will discuss the relationship between these functions.

In order to test the results obtained for the $2s^2$ wave functions with both methods, we employ a Feshbach projection formalism, but in a way inverse to the procedure used by Fano [21]. In general, the problem is how to obtain the true eigenvector $\Psi$ of the total Hamiltonian, starting from the restricted functions $\Phi$, and combining them as expressed in Eq.  (22). Here, we will obtain the restricted $\Phi_{2s^2}$ wave function from the eigenfunction of the total Hamiltonian (the wave function $\Psi_{2s^2}$). That is equivalent to projecting the eigenfunction $\Psi_{2s^2}$ onto the subspace $Q$, where $Q=1-P$. The procedure applied in order to obtain the projected wave $Q\vert\Psi_{2s^2}\rangle$ is to calculate

$\displaystyle Q \vert \Psi_{2s^2}(r_1,r_2)\rangle$ $\textstyle =$ $\displaystyle \vert \Psi_{2s^2}(r_1,r_2)\rangle -$ (23)
    $\displaystyle \vert\phi_{1s}(r_1)\rangle\langle\phi_{1s}(r_1)\vert\Psi_{2s^2}(r_1,r_2) \rangle -$  
    $\displaystyle \vert\phi_{1s}(r_2)\rangle\langle\phi_{1s}(r_2)\vert\Psi_{2s^2}(r_1,r_2) \rangle +$  
    $\displaystyle \vert\phi_{1s}(r_1)\phi_{1s}(r_2)\rangle \langle\phi_{1s}(r_1)\phi_{1s}(r_2)\vert\Psi_{2s^2}(r_1,r_2)\rangle.$  

If the numerical procedures are correct, the new wave function $Q\vert\Psi_{2s^2}\rangle$ must be very similar to the $\Phi_{2s^2}$ wave function obtained from the relaxation method. Graphic results of the probability densities for the functions $Q\vert\Psi_{2s^2}\rangle$ and $P\vert\Psi_{2s^2}\rangle$ are shown in Fig. 5. It is clear that the $Q\vert\Psi_{2s^2}\rangle$ function is very similar to the $\Phi_{2s^2}$ function displayed in Fig.  2 (first figure). On the other hand, the $P\vert\Psi_{2s^2}\rangle$ function has the form of a $1sks$ function, like the second of the continuum functions displayed in Figure 3. The overlap between the $\Psi_{2s^2}$ function (obtained by diagonalization) and the $\Phi_{2s^2}$ function (obtained by relaxation) is $\langle\Phi_{2s^2}\vert\Psi_{2s^2}\rangle= -0.990$, which accounts for all the $1skl$ components presented in the eigenvector $\Psi_{2s^2}$. The overlap of the relaxed function with the new projected function is $\langle\Phi_{2s^2}\vert Q \Psi_{2s^2}\rangle= -0.998$, an indication of the great similarity between the relaxed and the Feshbach projected functions. It is also consistent with the fact that the overlap $\langle\phi_{1s}(r1)\phi_{1s}(r2)\vert\Phi_{2s^2}(r_1,r_2)\rangle=0.002$, meaning that the projection mechanism in the relaxation did not eliminate completely the $1s$ components.

c. Time-dependent propagation

As we propagate the Schrödinger equation in imaginary time [Eq.(13)], we approach the (real) energy of the doubly excited level. We need to propagate it further in real time in order to obtain the imaginary part of the energies, i.e., the lifetime of the level. Figure 6 shows the evolution of the probability density of the wave function $\Phi_{2s^2}(r_1,r_2,t)$ under the time propagation with the Schrödinger equation,

\begin{displaymath} \Phi_{2s^2}(r_1,r_2,t) = e^{-iH(r_1,r_2)t} \Phi_{2s^2}(r_1,r_2,t_0) , \end{displaymath} (24)

where $\Phi_{2s^2}(r_1,r_2,t_0)$ is the wave function $\Phi_{2s^2}$ obtained using the constrained relaxation method.

The different frames of the figure show snapshots at early times from the beginning of the propagation. The velocity of the autoionizing electron can be obtained approximately from the energy difference between the $2s^2$ level and the He$^+$ ground state, i.e., $k = \sqrt{2 \epsilon_k}$, where $\epsilon_k= \epsilon_{2s^2} - \epsilon_{1s} = -0.71 - (-2) = 1.3$ a.u. is the energy of the free electron. That gives a velocity of 1.6 a.u. for the free autoionized electron. Taking into account the size of the box (22.5 a.u.), this free electron will rebound from the box borders at an approximate time of 14 a.u. Snapshots of the propagation of the $\Phi_{2s^2}(r_1,r_2,t)$ wave function are shown at the times 3.6, 7.3, 10.8, and 14.4 a.u.. The pictures show how the initial pure $2s^2$ wave function acquires a continuum feature, as it evolved in time. At the beginning, a slow deformation of the wave is developed, with the probability dropping through the wings of the wave. Then, the characteristic pattern of a $1sks$ continuum wave is formed at the axes, advancing rapidly to the box borders. After this time, the wave function bounces back and we found oscillations which under specific conditions can return the propagated function to roughly the original shape at particular times. We are aware of the possibility of introducing a mask function in order to absorb the wave at the borders preventing the rebound, but this revival of the original wave function can have a particular further interest. The full movie for the propagation can be downloaded at the author's personal webpage:

Propagation of the wavefunction - Movie in gif format (15 Mb).

Propagation of the probability density - Movie in gif format (11 Mb).

Movie in mpeg format (193 Kb).

The dynamical behavior of autoionization in two-electron systems has been studied previously by monitoring the decay of the autoionizing state $\Phi$ in time [16]. This requires the computation of the autocorrelation amplitude defined by

\begin{displaymath} A_{\Phi}(t) = \vert \langle\Phi(r_1,r_2,t) \vert \Phi(r_1,r_2,t_0)\rangle \vert, \end{displaymath} (25)

where $\Phi(r_1,r_2,t)$ represents the initial relaxed wave function evolved in time. This method has been used to study autoionization in a one-dimensional two-electron model [22], a two-dimensional two-electron model  [16], and in an S-wave model [16,23]. The result for $A_{2s^2}(t)$ is shown in Fig. 7. In this case, in order to allow the study of the propagation for longer times, we extended the numerical mesh to 75 a.u. (i.e., we used 500 points, in place of 150). An exponential fit of the form exp( $-\frac{\Gamma}{2}t)$ to the data (shown in the same figure by a dash curve) yields a value of $\Gamma_{2s^2}=2.7\times10^{-3}$ a.u., compared to the width determined by Draeger et al. [20] of $3.24\times10^{-3}$ a.u.

We found many other interesting features propagating our wave functions in time. First, we propagated the diagonalized eigenvectors $\Psi_i$, constructing the functions $\Psi_i(r_1,r_2,t)$. We found that the autocorrelation amplitudes $A_{\Psi_i}(t)$ for these functions are unity, even for a very long time. This is a good test for the numerical accuracy of the propagation method. We checked that the real and imaginary part of the autocorrelation function oscillate with a period corresponding to $2\pi/\epsilon_i$, where $\epsilon_i$ is the energy of the level $i$. We also tested the overlap of the propagated relaxed wave function $\Phi_{2s^2}(r_1,r_2,t)$ with the diagonal eigenvectors $\Psi_i$. These overlapping values also remained constant for a very long time, even after the wave function bounces back many times against the boundaries. Projection of the propagated relaxed wave function $\Phi_{2s^2}(r_1,r_2,t)$ with product functions $1sks$, on the other hand, does depend on time. After an initial predominance of low energies, the overlap profile shows a steep peak centered at $k$ close to 1.6 a.u. (as is expected).

One can notice the presence of wiggles in the autocorrelation amplitude at the initial steps of the time evolution (see Fig. 7.). In order to understand the origin of these wiggles, we repeated the calculations with a smaller number of mesh points, in this case 75 points separated by $\Delta_r= 1$ a.u., and the results are plotted in Fig. 8. Since the eigenfunctions $\Psi_i$ obtained from the diagonalization process constitute a complete group, we can expand the relaxed $\Phi_{2s^2}$ wave function as

\begin{displaymath} \Phi_{2s^2} = \sum_i^P c_i \Psi_i . \end{displaymath} (26)

Thus, the autocorrelation amplitude can be writen
$\displaystyle A_{\Phi_{2s^2}}(t)$ $\textstyle =$ $\displaystyle \vert \langle\Phi_{2s^2}(r_1,r_2,t) \vert \Phi_{2s^2}(r_1,r_2,t_0... ..._1,r_2)t}\Phi_{2s^2}(r_1,r_2,t_0) \vert \Phi_{2s^2}(r_1,r_2,t_0)\rangle \vert =$  
  $\textstyle =$ $\displaystyle \vert \langle \sum_j^P c_j e^{-iE_jt}\Psi_j(r_1,r_2) \vert \sum_k... ...i_k(r_1,r_2) \rangle \vert = \vert \sum_j^P \vert c_j\vert^2 e^{-iE_jt} \vert =$  
  $\textstyle =$ $\displaystyle \sqrt{ \sum_j^P \vert c_j\vert^4 + 2\sum_{j>k}^P \vert c_j\vert^2\vert c_k\vert^2\cos{(E_k - E_j)}t },$ (27)

where $P$ is the dimension of the space defined in Eq. (4).

We calculated the final expression in Eq. (26), cutting off the number P of terms in the sum to a different number of terms p; these results are shown in Fig. 8. We found that the $c_i$ coefficients in the expansion (26) have a very sharply peaked distribution around the $\Psi_i$ states located close to the $\Phi_{2s^2}$ function. However, a large number of terms in the expansion (27) are needed in order to reproduce an exponential behavior similar to the autocorrelation amplitude of the propagated wave function. The relatively large size of the box allows us to study the propagation of the wave function until long times (about 150 a.u.). At times beyond 150 a.u., an important part of the propagating wave function reaches the boundaries, bouncing back to the origin. Therefore, at long times the autocorrelation amplitude is no longer a decreasing exponential; it has an oscillatory behavior- it shows interferences between the propagating and bouncing back portion of the wave and other irregularities. Surprisingly, the expansion (27) reproduces perfectly the autocorrelation amplitude of the propagated wave function including the oscillations and the irregularities at long times. It is also interesting to notice how the frequency of the oscillations is related to the size of the box. For bigger boxes the propagating wave function takes more time to reach the boundaries, so the frequency of the oscillations is smaller. On the other hand, for bigger boxes the eigenvalue spectra are more dense. Therefore, the energy differences $(E_k-E_j)$ in the leading terms of the sum (27) are smaller, tracking exactly the oscillations in the autocorrelation amplitude.

As we pointed out above, our result may improve by using a finer radial mesh. For a summary of our findings, Table III shows our results for the calculated energies E $_{2s^2}= \epsilon_{2s^2} + i \Gamma_{2s^2}/2$ for the $2s^2$ level in He, and $\epsilon_{1s}$, the energy of the $1s$ level in He$^+$, for different numerical grids.

Table III: Calculated energies (E $_{2s^2}= \epsilon_{2s^2} + i \Gamma_{2s^2}/2$) for the $2s^2$ level in the S-wave model He, and $\epsilon_{1s}$, the energy of the $1s$ level in He$^+$, in atomic units.

\begin{ruledtabular} \begin{tabular}{lcc} ~&He E$_{2s^2}$\ & ~~~~He$^+$\ $\epsil... ...aeger:94} & $-0.7227 - i~0.00162$\ & ~~~$-2.0000$\end{tabular}\end{ruledtabular}


4. Conclusions

We have developed an exact numerical method for calculating the full spectrum of the helium atom in a box. It consists of a direct diagonalization of the Hamiltonian in a numerical grid. Since we have a particular interest in the doubly excited autoionizing levels, we developed an additional method, consisting of a constrained relaxation of the wave function. Calculations were performed for the S-wave model He, in which only spherical states are taken into account. Even though our goal is not to present the best wave functions but to obtain a good insight into the physical correlations, both methods show very good agreement among them and with other existing calculations. A Feshbach projection of the $\Psi_{2s^2}$ eigenvector is absolutely consistent with the relaxed $\Phi_{2s^2}$ wave function, proving the reliability of the numerical procedure. Time propagation of the relaxed wave function allows the calculation of the autoionization rate by monitoring the autocorrelation amplitude in time. Furthermore, computer animations were produced to visualize the dynamics of the autoionization process. Work is in progress to present better numerical results, including total close-coupling results, and for the calculations of other processes involving transitions among the different continuum states.

5. Acknowledgments

We would like to thank Professor Jorge Miraglia for many fruitful discussions. Computational work was carried out at the National Energy Research Supercomputer Center in Oakland, CA. These computational resources are supported through a grant from Scientific Discovery through Advanced Computinc (SciDAC) (U.S. Department of Energy), administered through Auburn University. \end{acknowledgments}">

6. Bibliography

1
M.S. Pindzola and D.R. Schultz, Phys. Rev. A 53, 1525 (1996).

2
M.S. Pindzola and F. Robicheaux, Phys. Rev. A 54, 2142 (1996).

3
M.S. Pindzola and F. Robicheaux, Phys. Rev. A 55, 4617 (1997).

4
D.M. Mitnik, D.C. Griffin, and M.S. Pindzola, Phys. Rev. Lett. 88, 173004 (2002).

5
W. Jaskolski, Phys. Rep. 271, 1 (1996); R.C. Ashoori, Nature (London) 379, 413 (1996); S. Bednarek, B. Szafran, and J. Adamowski, Phys. Rev. B 59, 13036 (1999); L.S. Costa, F.V. Prudente, P.H. Acioli, J.J. Soares Neto, and J.D.M Vianna, J. Phys. B 32, 2461 (1999); J.P. Connerade, V.K. Dolmatov, and P.A. Lakshmi, ibid. 33, 251 (2000).

6
D.L. Feder, M.S. Pindzola, L.A. Collins, B.I. Schneider, and C.W. Clark, Phys. Rev. A 62, 053606 (2000)

7
A. Temkin, Phys. Rev. 126, 130 (1962).

8
R. Poet, J. Phys. B 11, 3081 (1978).

9
For example, the real-symmetric diagonalization routine DSYEVX from the LAPACK package, or its parallel version PDSYEVX.

10
D.M. Mitnik, M.S. Pindzola, D.C. Griffin, and N.R. Badnell, J. Phys. B 32, L479 (1999).

11
S.E. Koonin, Computational Physics (Benjamin-Cummings, Menlo Park, CA, 1986), p. 172.

12
J.B. Anderson, J. Chem. Phys. 63, 1499 (1975).

13
P. Balcou, A. L'Huillier, and D. Escande, Phys. Rev. A 53, 3456 (1996).

14
K.C. Kulander, K.R. Sandhya Devi, and S.E. Koonin, Phys. Rev. A 25, 2968 (1982).

15
C. Bottcher, D.R. Schultz, and D.H. Madison, Phys. Rev. A 49, 1714 (1994).

16
D.R. Schultz, C. Bottcher, D.H. Madison, J.L. Peacher, G. Buffington, M.S. Pindzola, T.W. Gorczyca, P. Gavras, and D.C. Griffin, Phys. Rev. A 50, 1348 (1994).

17
H. Feshbach, Ann. Phys. (N.Y.) 19, 287 (1962).

18
M.S. Pindzola, D. Mitnik, and F. Robicheaux, Phys. Rev. A 62, 062718 (2000).

19
S.P. Goldman, Phys. Rev. Lett. 78, 2325 (1997).

20
M. Draeger, G. Handke, W. Ihra, and H. Friedrich, Phys. Rev. A 50, 3793 (1994).

21
U. Fano, Phys. Rev 124, 1866 (1961).

22
S.L. Haan, R. Grobe, and J.H. Eberly, Phys. Rev. A 50, 378 (1994).

23
W. Ihra, M. Draeger, G. Handke, and H. Friedrich, Phys. Rev. A 52, 3752 (1995).

24
URL: http://www.df.uba.ar/users/dmitnik

7. Figures

Figure 1.

Figure 2.

Figure 3.

Figure 4.

Figure 5.

Figure 6.

Figure 7.

Figure 8.