Entanglement measures in doubly excited states and nondispersive wave packets in planar helium

We calculate the amount of entanglement for bound and doubly excited singlet states below the second, third and fourth ionization thresholds in planar helium. Furthermore, we investigate the dynamics of entanglement in nondispersive two-electron wave packets generated by the resonant coupling of frozen planet states of helium. We show that it is possible to obtain a two-electron nondispersive wave packet for which the time-average entanglement is greater than for each of the unperturbed coupled states.


Introduction
Entanglement is one of the most significant characteristics in the quantum mechanical description of multipartite systems [1] and plays an important role in understanding fundamental phenomena such as quantum correlations [2]. Motivated not only for the interest on the foundations of quantum mechanics but also for the possible technological applications in quantum information, quantum computing [3], quantum cryptography [4], and quantum teleportation [5] there has been an increasing activity in the entanglement properties in atomic systems. While most studies on entanglement in atoms have been focused on bound levels of simplified models [6][7][8][9][10][11], there is rather few information on the entanglement properties of realistic atomic systems. In this regard, the helium atom enjoys an important place. It is the simplest naturally available atomic species which contains the electron-electron interaction, and, therefore a natural candidate for the investigation of entanglement properties in atoms. Previous investigations have been restricted to singly excited states [9,[12][13][14]. More recently, entanglement measures have been calculated for doubly excited states (DES) below the second ionization threshold in helium [15,16]. All these investigations have been restricted to the amount of entanglement in atomic states of the unperturbed atom and, up to our knowledge, there is no information about the dynamics of the amount of entanglement in driven two-electron atoms.
The purpose of the present work is to calculate the amount of entanglement for singlet bound and doubly excited states below the second, third and fourth ionization thresholds in planar helium. Furthermore, we report the first calculations of the dynamics of entanglement in certain two-electron wave packets. We pay special attention to particular doubly excited states of helium called frozen planet states (FPS) which are associated with highly correlated classical configurations [17,18]. Under near-resonantly periodic driving, FPS transform into nondispersive two-electron wave packets (NDWP) [19][20][21][22]. That is quantum objects that propagate without dispersion along the classical periodic orbits of frozen planet configuration for long periods of time [22].
This paper is organized as follows. In section 2, we briefly discuss the spectrum of the unperturbed planar helium, the approach to investigate FPS under periodic driving, and the entanglement measure for our system of two identical fermions. In section 3 we present and discuss the results for the amount of entanglement in bound, doubly excited states and nondispersive wave packets. Finally, section 4 contains the summary, some conclusions and perspectives.

Unperturbed atom
In the center of mass system within the infinite nucleus mass approximation, the dynamics of a two-electron atom is governed by the following Hamiltonian given in atomic units (a.u.), where R(θ) is the complex rotation operator and θ the rotation angle. In particular, the product of probability amplitudes r r | , can be obtained within the single pole approximation [27]. Its expression reads , , (6) is written as , , which leads to the electronic density Expressions (10) and (11) seem to be independent of the rotation angle θ since the backward rotation r R q á - | ( )of the spatial coordinate compensates the rotation of the state i, Y ñ q | . However, due to the normalization condition (3) it holds e d r J r R 1, 12 where r p  is the fourth dimensional position vector in parabolic coordinates. As a consequence, going back to Cartesian coordinates, we obtain d r r R Re cos 3 . 13 Therefore, the normalization factor 1 cos 3q ( ) has to be included in the wave function r E f  ( ).

Driven helium atom
In the dipole approximation, length gauge and neglecting relativistic effects, the Hamiltonian for the helium atom in presence of an external electromagnetic field reads where H 0 is the unperturbed Hamiltonian (1) and the time-periodic field is linearly polarized along the x direction with amplitude F and frequency ω.
In our approach to solve the time-dependent Schrödinger equation (TDSE) we use a spectral method combined with Floquet theory [28,29]. The method is based on the solution of the TDSE in the atomic basis i L j ñ {| }, that is the set composed of eigenstates of the unperturbed Hamiltonian | , for a given value of the total angular momentum L. In this basis, and within the framework of Floquet theory, the matrix representation of (15) is given by where Φ i is the solution vector in the atomic basis, h 0 is the diagonal matrix containing the eigenvalues of H 0 , F the matrix representation of F(x 1 +x 2 ) and k is the Floquet quantum number which gives the amount of photons exchanged between the atom and the field.
The time evolution of a Floquet eigenstate i F ñ | is governed by the rotated time evolution operator [25,26] U t t e e e R R , , In this way, the dynamics of the probability amplitude product is given by

Unperturbed frozen planet states
The classical frozen planet configuration (FPC) is an asymmetric one-dimensional configuration where both electrons are located on the same side of the nucleus. In this dynamically stable configuration [30][31][32], the inner electron precesses on highly eccentric ellipses and the outer electron remains nearly 'frozen' around the equilibrium distance of the adiabatic potential obtained averaging the fast oscillations of the inner electron. The classical two-electron dynamics restricted to the collinear frozen planet configuration is regular. That can be observed in the phase space configuration presented in figure 1 (right), which is visualized within a Poincaré surface of section obtained by plotting the position x 1 and the momentum p x 1 of the outer electron every time the inner electron reaches the nucleus. The numerical computation of the classical dynamics is achieved by the previous regularization of the equations of motion by Kustaanheimo-Stiefel transformations [33,34].
Starting from the third series (N=3) of the helium spectrum, the regular phase space of the FPC supports frozen planet states. These are organized in series converging to single ionization thresholds. We denote by n N F  the FPS occupying the n F position in the series converging to the Nth ionization threshold.
To identify the FPS on the helium spectrum we take into account that these states have long lifetime compared to other resonances in the same energy regime. In addition, due to their localization properties along the FPC, the expectation value of the cosine of the angle θ 12 between the two electrons is close to unity ( cos 1 12 q á ñ  ). We verify these localization properties in phase space with the help of Husimi distributions (see [22] for technical details), for instance, figure 1 shows the phase space projections of the ground (middle) and first excited (right) FPS below the N=4 single ionization threshold (energies and decay rates for those states are identified in table A1). These quantum states are clearly localized along the periodic orbits of the frozen planet configuration in figure 1 (left).

Driven frozen planet states
In the presence of an external time-periodic electromagnetic field, the dynamics of the collinear FPC takes place in a five-dimensional phase space described by the positions and momenta of the electrons and by the phase ω t of the driving field. For near resonant driving ω≈ω I of amplitude F<F I (where ω I is the frequency of small oscillations of the outer electron around the equilibrium position and F I is the minimum static field necessary to ionize the FPC), the frozen planet dynamics is not affected in a sensitive way and the time scale separation between the fast Kepler oscillations of the inner electron and the slow motion of the outer electron makes possible to map the phase space structure onto a two-dimensional surface by a two-step Poincaré section method [19,35].
In figure 2 is depicted the two-step Poincaré surface of section obtained for field amplitude F=2×10 −5 a.u. and frequency ω=0.003 6 a.u., for three different field phases ω t. Under the action of the external electromagnetic field, the classical phase space separates into two regular regions: the intrinsic island (centred at 30 a.u.) and the resonance island which oscillates around the intrinsic island with the field frequency.
Theoretical studies [19][20][21][22] suggest that under near-resonantly periodic driving frozen planet states might transform into two-electron nondispersive wave packets, which propagate along the frozen planet trajectory without dispersion.
The solution of the time-dependent Schrödinger equation (15) for the description of such highly excited driven states is a non-trivial task [36]. Previous large scale computations have been unable to completely characterize such objects [22,37].
We have developed an efficient approach to study the dynamics of doubly excited states of helium under electromagnetic driving. The method, briefly described in section 2.4, has been used in [38] to characterize NDWP in helium. There, we have been found that the most important features of the nondispersive wave packets are a consequence of the one-photon resonant coupling between two frozen planet states of different angular momentum. In section 3.2, we shall use this result to identify NDWP in the Floquet spectrum of driven helium, below the N=3 and 4 ionization thresholds.  , the maximum of the probability density is localized along a periodic orbit with higher energy in the classical configuration. Figure 2. Classical phase space structure of the collinear driven frozen planet configuration for field parameters F=2×10 −5 a.u. and ω=0.003 6 a.u., and field phase ω t=0 (left), ω t=π/2 (center) and ω t=π (right). The intrinsic island remains basically unaffected by the field during the time evolution, while the resonance island oscillates around the intrinsic island with the same field frequency.

Entanglement measure
For a pure state Fñ | of two identical fermions, the Schmidt decomposition in an orthonormal basis of oneparticle states iñ {| }is given by [39,40] i where the Schmidt coefficients λ i satisfy 0 1 If the state Fñ | can be expressed as a single Slater determinant (only one Schmidt coefficient) the state has no entanglement [40][41][42]. In such a way, the decomposition (19) leads to a measure of the amount of entanglement exhibited by the two-fermion state, namely [40, 41] where Tr is the reduced density matrix (obtained by tracing the two-particle density matrix over one of the two particles).
The state of the two-electron system is described by wave functions of the form r r , , , )and χ(σ 1 , σ 2 ) correspond to the coordinate and spin wave functions, respectively. The corresponding density matrix is expressed as the product where r r r r r r d r , , .
The expressions for r r * f f ¢   ( ) ( )are obtained from (10) for unperturbed states and from (18) for driven FPS. For planar helium, the expression (24) corresponds to a 8-dimensional definite integral. We perform this calculation using the VEGAS algorithm included in the Cuba library 3 [43,44]. VEGAS is widely used for Monte Carlo multidimensional numerical integration and is primarily based on importance sampling, but it also does some stratified sampling as a variance-reduction technique. The algorithm iteratively builds up a piecewise constant weight function, represented on a rectangular grid, which is refined after each iteration. Further details can be found in [45].

Entanglement in bound and doubly excited states
We calculate the amount of entanglement for the lowest bound and doubly excited S 1 and P 1 states below the N=2, 3 and 4 ionization thresholds. Going beyond the 4th ionization threshold is a difficult task due to the huge number of integrand evaluations needed in (24) to obtain a converged Monte Carlo integration: the region in configuration space where the probability density is non-negligible scales in each direction as N 2 [22]. Therefore, the integration volume in the 8th dimensional integral (24) scales as N 16 . This is approximately the scaling factor of the number of integrand evaluations and also of the computation time. For instance, while the computation time for states below the 4th ionization threshold is three days, the expected one for states below the 5th ionization threshold is 35 times large. This also might explain why available full three-dimensional calculations are restricted to states below the second ionization threshold.
In table 1, we present the energies and the amount of entanglement for singlet bound states of planar helium. For bound S 1 states the entanglement increase with energy, whereas for bound P 1 states, the entanglement decrease with energy. These results are consistent with the entanglement behaviour observed in previous works in three-dimensional helium for low lying excited states [12,15]. In the case of DES, it can be observed that the amount of entanglement have different behaviour according to their value of cos 12 q á ñ, which, as we mentioned in section 2.3, is an important quantity to identify the frozen planet states which lead to the formation of NDWP. In figure 3 (left panels) we present the amount of entanglement as a function of the energy for doubly excited S 1 states. The results are organized in series that can  be related to the value of cos 12 q á ñ as we observe in the right panels of figure 3, where we show the expectation values of cos 12 q as a function of the energy. Similar behaviour is presented for doubly excited P 1 states in figure 4. The converged values of the entanglement measure are presented in the appendix.
In our results for DES below the N=2 ionization threshold, the amount of entanglement in the different series increases to reach a maximum value and then decreases, with exception of the P 1 states with  . This tendency is also exhibited by reported results in three-dimensional helium [15]. For DES below the N=3 ionization threshold, the maximum entanglement values correspond to FPS (marked with blue squares), while the minimum values correspond to the states associated with the classical eZe configuration (marked with red circles). In the case of DES below the N=4 ionization threshold, the FPS (blue squares) also have the maximum entanglement values (above 0.92), which confirms the high degree of correlations in FPS.
By considering only the results for frozen planet states we observed that: below each ionization threshold, entanglement values for P 1 states are greater than for S 1 states, and the amount of entanglement increases with N.
It has been shown that for cos 0.5 12  q á ñ - [46], doubly excited states gather into series labeled by F=N−K, where N and K are approximate quantum numbers from Herrick's classification [47][48][49]. Furthermore, for highly doubly excited states it is also possible to identify well defined series for the case of antisymmetric states ( cos 0.5 12  q á ñ ). From the relation K N cos 12 q  -á ñ [49,50], the classical eZe configuration is associated with K max =N−1 whereas that K min =−(N−1) is identified with the classical Zee configuration which corresponds to the frozen planet configuration. In this regard, our results suggest that the entanglement values could be used to classify this kind of highly doubly excited states. However, this assumption needs to be verified with entanglement calculations for higher excitations in planar and threedimensional helium.

Entanglement in nondispersive wave packets
Here, we focus our attention on the entanglement properties of nondispersive wave packets produced by the one-photon coupling of two frozen planet states. This is achieved by solving (16) in the atomic basis containing only the states a S N . In that basis, the matrix representation of the Floquet Hamiltonian is given by , and E α is the complex energy of the resonance state añ | . In the Floquet spectrum obtained after diagonalization of (26), we identify nondispersive wave packets below the N=3 and 4 ionization thresholds. For N=3 we use the field parameters F=4×10 −5 a.u. and ω=0.006 6 a.u.. The wave packet below the N=3 ionization threshold has energy E=−0.342 610 1 a.u. and decay rate Γ=8×10 −7 a.u. For the case of N=4 the field parameters are F=2×10 −5 a.u. and ω=0.003 6 a.u., and the energy and decay rate wave of the packet are E=−0.175 769 8 a.u. and Γ=4.5×10 −6 a.u., respectively.  Figure 5 displays the phase space projection of the N=4 singlet wave packet, which can be compared with the classical phase space structure presented in figure 2, obtained for the same driving field parameters. Clearly, the wave packet motion in figure 5 is associated with the dynamics of the resonance island in the classical phase space. Figure 6 shows the amount of entanglement calculated as a function of the driving field phase for the the two wave packets. For comparison, we also show the amount of entanglement for the FPS involved in the formation of the NDWP and the time-average entanglement  over one period T For the wave packet below the third ionization threshold ( figure 6 (top)) the time-average entanglement  was approximately the mean value between the amount of entanglement of the S 1 2 3  and P 1 1 3  frozen planet states, while for the wave packet below the fourth ionization threshold (figure 6 (bottom)), the time-average entanglement was greater than the values for the initial FPS. In both cases N=3 and N=4, we observe that the maximum values for the wave packet entanglement were for the driving field phases ω t=π/2 and ω t=3π/2.

Summary and conclusions
In this work, we have calculated the amount of entanglement of bound and doubly excited states below the N=2, 3 and 4 ionization thresholds for planar helium. Additionally, we computed the entanglement of nondispersive two-electron wave packets produced by frozen planet states driven by an external periodic field.
We found that for bound S 1 states the entanglement increase with energy while for bound P 1 states the entanglement decrease with energy. For DES below the N=2 ionization threshold, the entanglement values are organized in series where the entanglement reaches a maximum before decreasing. These results are consistent with that observed in three-dimensional calculations [15]. In both cases N=3 and N=4, the entanglement values for FPS are higher than for the other series of states (above 0.92 for N=4 and 0.89 for N=3) and we also found that the amount of entanglement of P 1 FPS is greater than for S 1 FPS for each series. We expect that this behaviour could be verified on three-dimensional computations. On the other hand, it would be interesting to extend the calculations to highly doubly excited states and higher angular momenta, and for those states verify the general trend of entanglement to increase with the strength of the inter-particle interaction and its nonvanishing value in the limit of zero interaction [51].
In the case of driven FPS, we found that the time-average entanglement for NDWP below the N=4 ionization threshold is greater than the entanglement of the unperturbed FPS involved in its formation. Although only fixed driving field parameters were used here, the NDWP produced by driven FPS maintain their behaviour under driving field amplitude and frequency variations [38]. Thus, investigation of the effect of the field parameters variation on the amount of entanglement in NDWP is feasible.

Appendix. Entanglement data for DES
This appendix contains the entanglement values for the doubly excited S 1 and P 1 states below the N=2, 3 and 4 ionization thresholds, discussed in section 3.1.
In tables A1 and A2 we present the energies, decay rates and the amount of entanglement for doubly excited S 1 and P 1 states, respectively. Additionally, we include the expectation value of cos 12 q . In the second column of each table we identify the lowest FPS of the series converging to the third and fourth ionization thresholds.  Table A2. Energies and entanglement of the lowest P 1 doubly excited states below the Nth single ionization threshold. Decay rates and the expectation value of cos 12 q are also shown.