One-Electron Energy Spectra of Heavy Highly Charged Quasimolecules: Finite-Basis-Set Approach

Abstract

The generalized dual-kinetic-balance approach for axially symmetric systems is employed to solve the two-center Dirac problem. The spectra of one-electron homonuclear quasimolecules are calculated and compared with the previous calculations. The analysis of the monopole approximation with two different choices of the origin is performed. Special attention is paid to the lead and xenon dimers, Pb$${}^{82+}$$–Pb$${}^{82+}$$–e$${}^{-}$$ and Xe$${}^{54+}$$–Xe$${}^{54+}$$–e$${}^{-}$$, where the energies of the ground and several excited $$\sigma$$-states are presented in the wide range of internuclear distances. The developed method provides the quasicomplete finite basis set and allows for the construction of perturbation theory, including within the bound-state QED.

1. Introduction

Due to the critical phenomena of the bound-state quantum electrodynamics, such as spontaneous electron–positron pair production, quasimolecular systems emerging in ion–ion or ion–atom collisions have attracted much interest [1,2,3,4,5,6,7,8,9]. While collisions of highly charged ions with neutral atoms are presently available for experimental investigations, in particular at the GSI Helmholtz Center for Heavy Ion Research [10,11,12], the upcoming experiments at the GSI/FAIR [13], NICA [14], and HIAF [15] facilities might even allow the observation of the heavy ion–ion (up to U${}^{92+}$–U${}^{92+}$) collisions. The relativistic dynamics of the heavy-ion collisions has been investigated for decades by various methods; see, e.g., [7,8,9,16,17,18,19,20,21,22] and the references therein. The theoretical predictions of the quasimolecular spectra are also in demand for the analysis of the experimental data in these collisions.

Within the Born–Oppenheimer approximation, the one-electron problem is reduced to the Dirac equation with the Coulomb potential of two nuclei at a fixed internuclear distance D. This problem was investigated previously by a number of authors; see, e.g., [16,21,23,24,25,26,27,28,29,30,31,32,33,34,35,36,37,38]. The majority of these works relied on the partial-wave expansion of the two-center potential in the center-of-mass coordinate system. Alternative approaches include, e.g., the usage of the Cassini coordinates [34] and the atomic Dirac–Sturm basis-set expansion [36,37]. We consider the method based on the dual-kinetic-balanced finite-basis-set expansion [39] of the electron wave function for axially symmetric systems [40]. The results for the ground state of uranium dimers with one and two electrons were already presented in [41]. In the present work, we extend the one-electron calculations to the lowest excited $\sigma$-states and present the results for the one-electron dimers, Pb${}^{82+}$–Pb${}^{82+}$–e${}^{-}$ and Xe${}^{54+}$–Xe${}^{54+}$–e${}^{-}$. For the ground state, we demonstrate the accuracy of this method for the nuclear charge numbers Z from 1 to 100 at the so-called “chemical” distances, $D=2/Z$ a.u. We also investigate the difference between the two-center values and those obtained within the monopole approximation.

The relativistic units ($\hslash =1$, $c=1$, $m_e=1$) and the Heaviside charge unit ($\alpha =e^2/\left(4\pi \right),e<0$) are used throughout the paper.

2. Method

In heavy atomic systems, the parameter $\alpha Z$ ($\alpha$ is the fine-structure constant and Z is the nuclear charge number), which measures the coupling of electrons with nuclei, is not small. Therefore, the calculations for these systems should be performed within the fully relativistic approach, i.e., to all orders in $\alpha Z$. With this in mind, we start with the Dirac equation for the two-center potential,

(1)$\begin{array}{c}\left[\overrightarrow{\alpha }·\overrightarrow{p}+\beta +V\left(\overrightarrow{r}\right)\right]{\mathsf{\Psi }}_n\left(\overrightarrow{r}\right)=E_n{\mathsf{\Psi }}_n\left(\overrightarrow{r}\right),\end{array}$

(2)$\begin{array}{c}V\left(\overrightarrow{r}\right)=V_{\mathrm{nucl}}^{\mathrm{A}}\left(\right|\overrightarrow{r}-{\overrightarrow{R}}_1\left|\right)+V_{\mathrm{nucl}}^{\mathrm{B}}\left(\right|\overrightarrow{r}-{\overrightarrow{R}}_2\left|\right).\end{array}$

Here, $\overrightarrow{r}$ and ${\overrightarrow{R}}_{1,2}$ are the coordinates of the electron and nuclei, respectively, $V_{\mathrm{nucl}}^{\mathrm{A},\mathrm{B}}\left(r\right)$ is the nuclear potential at the distance r generated by the nucleus with the charge $Z_{\mathrm{A},\mathrm{B}}$, and $\overrightarrow{\alpha }$ and $\beta$ are the standard $4\times 4$ Dirac matrices:

(3)$\beta =\left(\begin{array}{cc}I& 0\\ 0& -I\end{array}\right),\overrightarrow{\alpha }=\left(\begin{array}{cc}0& \overrightarrow{\sigma }\\ \overrightarrow{\sigma }& 0\end{array}\right),$

where $\overrightarrow{\sigma }$ is a vector of the Pauli matrices.

In the following, we consider the identical nuclei, i.e., $Z_{\mathrm{A}}=Z_{\mathrm{B}}$, with the Fermi model of the nuclear charge distribution:

(4)$V_{\mathrm{nucl}}\left(r\right)=-4\pi \alpha Z\underset{0}{\overset{\infty }{\int }}\frac{\rho \left(r^{\prime }\right)}{max(r,r^{\prime })}{r^{\prime }}^{2}d{r}^{\prime },\rho \left(r\right)=\frac{{\rho }_0}{1+exp(r-c)/a},$

where ${\rho }_0$ is the normalization constant, a is the skin thickness constant, and c is the half-density radius. For more details, see, e.g., [42].

The solution of Equation (1) is obtained within the dual-kinetic-balance (DKB) approach, which allows one to solve the problem of “spurious” states. These non-physical states occur in the spectrum upon solving the Dirac equation using the finite basis set [43]. Originally, this approach was implemented for spherically symmetric systems, such as atoms [39], using the finite basis set constructed from the B-splines [44,45]. Later, authors of [40] generalized it to the case of axially symmetric systems (A-DKB): they considered the atom in an external homogeneous field. This case was also considered within this method in [46,47,48] to evaluate the higher-order contributions to the Zeeman splitting in highly charged ions. In [41], we adapted the A-DKB method to diatomic systems, which also possess axial symmetry. Below, we provide a brief description of the calculation scheme.

The system under consideration is rotationally invariant with respect to the z-axis directed along the internuclear vector $\overrightarrow{D}={\overrightarrow{R}}_2-{\overrightarrow{R}}_1$. Therefore, the z-projection of the total angular momentum with the quantum number $m_J$ is conserved, and the electronic wave function can be written as,

(5)$\mathsf{\Psi }(r,\theta ,\phi )=\frac{1}{r}\left(\begin{array}{c}{G}_1(r,\theta )e^{i(m_J-\frac{1}{2})\phi }\\ G_2(r,\theta )e^{i(m_J+\frac{1}{2})\phi }\\ i{F}_1(r,\theta )e^{i(m_J-\frac{1}{2})\phi }\\ i{F}_2(r,\theta )e^{i(m_J+\frac{1}{2})\phi }\end{array}\right).$

The $(r,\theta )$-components of the wave function are represented using the finite-basis-set expansion:

(6)$\Phi (r,\theta )=\frac{1}{r}\left(\begin{array}{c}{G}_1(r,\theta )\\ G_2(r,\theta )\\ F_1(r,\theta )\\ F_2(r,\theta )\end{array}\right)\cong \sum _{u=1}^4\sum _{i_r=1}^{N_r}\sum _{i_{\theta }=1}^{N_{\theta }}{C}_{i_r,i_{\theta }}^u\mathsf{\Lambda }{B}_{i_r}\left(r\right)Q_{i_{\theta }}\left(\theta \right)e_u,$

where ${\left\{B_{i_r}\left(r\right)\right\}}_{i_r=1}^{N_r}$ are B-splines, ${\left\{Q_{i_{\theta }}\right\}}_{i_{\theta }=1}^{N_{\theta }}$ are Legendre polynomials of the argument $2\theta /\pi -1$, and ${\left\{e_u\right\}}_{u=1}^4$ are the standard four-component basis vectors:

(7)$e_1=\left(\begin{array}{c}1\\ 0\\ 0\\ 0\end{array}\right),e_2=\left(\begin{array}{c}0\\ 1\\ 0\\ 0\end{array}\right),e_3=\left(\begin{array}{c}0\\ 0\\ 1\\ 0\end{array}\right),e_4=\left(\begin{array}{c}0\\ 0\\ 0\\ 1\end{array}\right).$

The $\mathsf{\Lambda }$-matrix:

(8)$\begin{array}{c}\mathsf{\Lambda }=\left(\begin{array}{cc}I& -\frac{1}{2}{D}_{m_J}\\ -\frac{1}{2}{D}_{m_J}& I\end{array}\right),\end{array}$

(9)$\begin{array}{c}{D}_{m_J}=({\sigma }_{z}cos\theta +{\sigma }_{x}sin\theta )\left(\frac{\partial }{\partial r}-\frac{1}{r}\right)\\ +\frac{1}{r}({\sigma }_{x}cos\theta -{\sigma }_{z}sin\theta )\frac{\partial }{\partial \theta }+\frac{1}{rsin\theta }\left(i{m}_J{\sigma }_y+\frac{1}{2}{\sigma }_x\right),\end{array}$

imposes the dual-kinetic-balance conditions on the basis set. With the given form of $\mathsf{\Phi }$ and the finite basis set, one can find the corresponding Hamiltonian matrix $H_{ij}$ and the so-called basis function overlap matrix $S_{ij}$ (for more details, see [40]). The eigenvalues and eigenfunctions are found numerically by solving the generalized eigenvalue problem ${\sum }_j{H}_{ij}{C}_j={\sum }_{j}E{S}_{ij}{C}_j$ via the standard DSYGV subroutine from LAPACK. As a result, we obtain a quasicomplete finite set of wave functions and electronic energies for the two-center Dirac equation. Ground and several lowest excited states are reproduced with high accuracy while the higher excited states effectively represent the infinite remainder of the spectrum. The negative-energy continuum is also represented by the finite number of the negative-energy eigenvalues. This quasicomplete spectrum can be used to construct the Green function, which is needed for the perturbation theory calculations. So far, the leading non-trivial QED corrections—self-energy, vacuum polarization, and two-photon exchange—have been calculated [41] only within the so-called monopole approximation (see Equation (10) below). In order to evaluate these contributions on the basis of the two-center Dirac equation, the corresponding Green function has to be constructed, e.g., within the finite-basis-set approach presented in this paper.

3. Results

Relativistic calculations of the binding energies of heavy one-electron quasimolecules were presented, in particular, in [21,29,32,34,36,37]. Reference [37] provided nearly the most accurate up-to-date values for the very broad range of Z and taking into account the finite nuclear size. Therefore, we use these data for comparison with our results: see Table 1, where the ground-state energies are presented for $Z=1,\dots ,100$ at the so-called “chemical” distances, $D=2/Z$ a.u. We observe that the results are in good agreement, the relative deviation varying from $2\times {10}^{-6}$ for hydrogen to $5\times {10}^{-5}$ for $Z=100$. This deviation is consistent with our own estimation of the numerical uncertainty, which is evaluated by inspecting the convergence of the results with respect to the size of the basis set. In this calculation, up to $N_r=320$ B-splines and $N_{\theta }=54$ Legendre polynomials are used. For heavy nuclei, this number of basis functions ensures the uncertainty, which is comparable to or smaller than the uncertainty of the finite nuclear size effect at all internuclear distances from 0 to $2/Z$ a.u.

Next, we present the obtained one-electron spectra of the Pb${}^{82+}$–Pb${}^{82+}$–e${}^{-}$ and Xe${}^{54+}$–Xe${}^{54+}$–e${}^{-}$ quasimolecules in the wide range of the internuclear distances from a few tens of fermi up to the “chemical” distances.

In the present figures, only $\sigma$-states ($m_J=\pm \frac{1}{2}$) are shown. The precise quantum numbers are $m_J$ and parity, g (gerade) or u (ungerade). In addition, we determine the quantum numbers of the “merged atom”, i.e., the state of the system with the internuclear distance $D\to 0$, and put it to the left of the molecular term symbol, e.g., the ground state is $1s_{1/2}{\sigma }_g$.

In Figure 1, the energies of the ground ($n=1$) and first nine ($n=2,\dots ,10$) excited states of the Pb${}^{82+}$–Pb${}^{82+}$–e${}^{-}$ system as functions of the internuclear distance are shown. Here, n has no connection with the atomic principal quantum number; it simply enumerates the $\sigma$-states. To visually compare the data obtained with the ones by Soff et al., in Figure 2, we zoom the second plot from Figure 1 to match the scale of the corresponding figure from [16]. Although we cannot compare the numerical results, the plots for all the states under consideration appear to be in very good agreement: all the states are identified correctly, and all the crossings and avoided crossings appear at the same internuclear distances. The similar results for xenon, i.e., the energies of the ground ($n=1$) and first nine ($n=2,\dots ,10$) excited states of the Xe${}^{54+}$–Xe${}^{54+}$–e${}^{-}$ system are shown in Figure 3.

Furthermore, in Table 2 and Table 3, we compare the ground-state binding energies obtained within our approach for the two-center (TC) potential with those for the widely used monopole approximation (MA), where only the spherically symmetric part of the two-center potential is considered:

(10)$V_{\mathrm{mono}}\left(r\right)=\frac{1}{2}\underset{0}{\overset{\pi }{\int }}d\theta sin\theta [V_{\mathrm{nucl}}^A\left(\right|\overrightarrow{r}-{\overrightarrow{R}}_1\left|\right)+V_{\mathrm{nucl}}^B\left(\right|\overrightarrow{r}-{\overrightarrow{R}}_2\left|\right)].$

The corresponding results are obtained within the DKB approach [39]. Within the MA, the potential and all the results depend on where to place the origin of the coordinate system (c.s.). At the same time, for the TC potential, the results should be identical within the numerical error bars. We compare two different placements of the c.s. origin: $\left(1\right)$ at the center of mass of the nuclei and $\left(2\right)$ at the center of one of the nuclei; see Figure 4.

The agreement between TC(1) and TC(2) within the anticipated numerical uncertainty serves as a non-trivial self-test of the method, since the basis-set expansion (6) is essentially different for the two cases. In fact, due to the lower symmetry of the second c.s., the uncertainty of the TC(2) values is much larger and completely determines the difference between TC(1) and TC(2). The differences between the TC(1) and MA(1) results, presented in the second-to-last column, can be interpreted as the inaccuracy of the MA. In the last column, the differences between the MA(1) and MA(2) results are given, a kind of “inherent inconsistency” of the MA.

As one can see from these data, $\left|\mathrm{MA}\right(1)-\mathrm{MA}(2\left)\right|$ is roughly comparable to $\left|\mathrm{TC}\right(1)-\mathrm{MA}(1\left)\right|$. In other words, in case the TC values are unknown, one can estimate the magnitude of $\left|\mathrm{TC}\right(1)-\mathrm{MA}(1\left)\right|$ from $\left|\mathrm{MA}\right(1)-\mathrm{MA}(2\left)\right|$, except for the regions where the latter is anomalously small due to the sign change. This observation can be used to quantify the inaccuracy of the MA results for the various quantities evaluated within the MA, but not yet within the rigorous TC approach, which is generally much more demanding. In particular, this is the case of the self-energy, vacuum-polarization, and two-photon-exchange contributions calculated for heavy quasimolecules in the MA within the DKB approach, e.g., in [41].

4. Summary

In this work, the two-center Dirac equation is solved within the dual-kinetic-balance method [39,40]. The energies of the ground and several excited $\sigma$-states in such heavy diatomic systems as Pb${}^{82+}$–Pb${}^{82+}$–e${}^{-}$ and Xe${}^{54+}$–Xe${}^{54+}$–e${}^{-}$ are plotted as a function of the internuclear distance D. The ground-state energies at the “chemical” distances ($D=2/Z$ a.u.) are presented for one-electron dimers with $Z=1,\dots ,100$. The obtained data are compared with the available previous calculations, and a good agreement is observed. The comparison of the results for different origin placement of the coordinate system is used as a self-test of the method. The values obtained within the monopole approximation are also presented. It is shown that their dependence on the origin placement can serve to estimate the deviation from the two-center results, at least by the order of the magnitude.

The developed method, in addition to the energies and wave functions of the ground and lowest excited states, provides the quasicomplete finite spectrum. The Green function computed on the basis of this spectrum gives access, in particular, to the evaluation of the Feynman diagrams within the bound-state QED.

Author Contributions

A.A.K., D.A.G., V.M.S., and G.P. have contributed equally to this work. All authors have read and agreed to the published version of the manuscript.

Funding

The work was supported by the Foundation for the Advancement of Theoretical Physics and Mathematics “BASIS”, by the Russian Foundation for Basic Research (Grant Number 19-02-00974), by TU Dresden (DAAD Programm Ostpartnerschaften), and by G-RISC.

Institutional Review Board Statement

Not applicable.

Informed Consent Statement

Not applicable.

Data Availability Statement

Not applicable.

Acknowledgments

Valuable discussions with Timur Isaev, Ilia Maltsev, Aleksei Malyshev, Leonid Skripnikov, and Ilya Tupitsyn are gratefully acknowledged.

Conflicts of Interest

The authors declare no conflict of interest.

Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Figures and Tables

Figure 1. Electronic terms of the one-electron Pb82+–Pb82+ quasimolecule.

Figure 2. Electronic terms of the one-electron Pb82+–Pb82+ quasimolecule. Energies of states with n=6,…,9 (scaled).

Figure 3. Electronic terms of the one-electron Xe54+–Xe54+ quasimolecule.

Figure 4. Two different coordinate systems considered. Left: (1) the origin is at the center-of-mass of the system; right: (2) the origin is at the center of one of the nuclei.

Table 1

Comparison of the ground-state energies (in a.u.) of one-electron quasimolecular systems with $Z=1,\dots ,100$ at the internuclear distance $D=2/Z$ a.u.

Z	This Work	Dirac–Sturm [37]
1	−1.1026433	−1.102641581032
2	−4.4106607	−4.410654714140
10	−110.33722	−110.3371741499
20	−442.23969	−442.2392996469
30	−998.4194	−998.4214646525
40	−1783.5479	−1783.563450815
50	−2804.5304	−2804.571434254
60	−4070.971	−4071.036267926
70	−5595.889	−5595.926978290
80	−7397.003	−7397.028800116
90	−9498.452	−9498.588788490
92	−9957.567	−9957.775519122
100	−11,935.89	−11,936.41770218

Table 2

Ground-state binding energy $E_{1{\sigma }_g}$ (in eV) of the Pb${}^{82+}$–Pb${}^{82+}$–e${}^{-}$ quasimolecule for the two-center potential (TC) and for the monopole-approximation potential (MA), with the coordinate system origin at the center of mass of the nuclei $\left(1\right)$ and at the one of the nuclear centers $\left(2\right)$.

D, fm	$\mathbf{TC}\left(1\right)$	$\mathbf{TC}\left(2\right)$	$\mathbf{MA}\left(1\right)$	$\mathbf{MA}\left(2\right)$	$\left|\mathbf{TC}\right(1)-\mathbf{MA}(1\left)\right|$	$\left|\mathbf{MA}\right(1)-\mathbf{MA}(2\left)\right|$
40	−646,254	−646,254	−637,032	−598,564	9222	38,468
50	−614,504	−614,504	−604,643	−568,188	9861	36,455
80	−550,575	−550,575	−539,861	−506,742	10,714	33,119
100	−521,373	−521,373	−510,350	−478,423	11,023	31,927
200	−433,348	−433,347	−421,146	−392,345	12,202	28,801
250	−405,450	−405,450	−392,687	−365,185	12,763	27,502
500	−319,773	−319,769	−304,337	−283,510	15,436	20,827
$1/Z\left[\mathrm{a}.\mathrm{u}.\right]\simeq 645$	−289,068	−289,067	−272,212	−255,389	16,856	16,823
700	−279,462	−279,464	−262,095	−246,756	17,367	15,339
1000	−238,887	−238,873	−218,905	−211,937	19,982	6968
$2/Z\left[\mathrm{a}.\mathrm{u}.\right]\simeq 1291$	−212,020	−212,003	−189,652	−190,174	22,368	522

Table 3

Ground-state binding energy $E_{1{\sigma }_g}$ (in eV) of the Xe${}^{54+}$–Xe${}^{54+}$–e${}^{-}$ quasimolecule. The notations are the same as in Table 2.

D, fm	$\mathbf{TC}\left(1\right)$	$\mathbf{TC}\left(2\right)$	$\mathbf{MA}\left(1\right)$	$\mathbf{MA}\left(2\right)$	$\left|\mathbf{TC}\right(1)-\mathbf{MA}(1\left)\right|$	$\left|\mathbf{MA}\right(1)-\mathbf{MA}(2\left)\right|$
40	−192,031	−192,031	−191,860	−190,033	171	1827
50	−190,845	−190,845	−190,607	−188,314	238	2293
80	−187,217	−187,216	−186,775	−183,199	442	3576
100	−184,805	−184,805	−184,228	−179,895	577	4333
200	−173,425	−173,425	−172,190	−165,031	1235	7159
250	−168,242	−168,242	−166,695	−158,621	1547	8074
500	−146,803	−146,802	−143,860	−133,919	2943	9941
700	−133,710	−133,710	−129,828	−120,118	3882	9710
$1/Z\left[\mathrm{a}.\mathrm{u}.\right]\simeq 980$	−119,414	−119,413	−114,421	−105,971	4993	8450
1000	−118,529	−118,528	−113,464	−105,233	5065	8231
$2/Z\left[\mathrm{a}.\mathrm{u}.\right]\simeq 1960$	−89,269	−89,276	−81,433	−79,488	7836	1945