Charge conjugation symmetry in the ﬁnite basis approximation of the Dirac equation

: 4-component relativistic atomic and molecular calculations are typically performed within the no-pair approximation where negative-energy solutions are discarded, hence the symmetry between electronic and positronic solutions is not considered. These states are however needed in QED calculations, where furthermore charge conjugation symmetry becomes an issue. In this work we shall discuss the realization of charge conjugation symmetry of the Dirac equation in a central ﬁeld within the ﬁnite basis approximation. Three schemes for basis set construction are considered: restricted, inverse and dual kinetic balance. We ﬁnd that charge conjugation symmetry can be realized within the restricted and inverse kinetic balance prescriptions, but only with a special form of basis functions that does not obey the right boundary conditions of the radial wavefunctions. The dual kinetic balance prescription is on the other hand compatible with charge conjugation symmetry without restricting the form of the radial basis functions. However, since charge conjugation relates solutions of opposite value of the quantum number κ , this requires the use of basis sets chosen according to total angular momentum j rather than orbital angular momentum (cid:96) . As a special case, we consider the free-particle Dirac equation, where the solutions of opposite sign of energy are related by charge conjugation symmetry. We note that there is additional symmetry in those solutions of the same value of κ come in pairs of opposite energy.


Introduction
Consider an electron of charge q = −e and mass m e , placed in an attractive Coulomb potential φ(r). Upon solving the time-independent Dirac equation, one gets a set of solutions ψ i associated with energy levels E i which forms the spectrum that is shown in a pictorial way in fig.(1.a). The charge conjugation operation [1] (C-operation) relates a particle to its anti-particle. The C-conjugated solution Cψ i , describes the solution of the same equation but with opposite charge (a positron), flipping the spectrum, as shown in fig.(1.c). In the free-particle case, φ(r) = 0, the two spectra, left and right, coalesce into the spectrum in fig.(1.b), that contains no bound solutions, and where the C-operation relates positive-and negative-energy solutions of the same equation, i.e. Cψ ±E i = ψ ∓E i . Note that since the free-particle equation does not require the specification of the charge, it describes equally well electrons and positrons. The Dirac equation is the starting point for 4-component relativistic atomic and molecular calculations. In the former case, the high symmetry of the problem allows the use of finite difference methods, whereas molecular applications generally call for the use of finite basis expansions. Early calculations using finite bases were flawed since the coupling of the large and small components was not respected. Spurious solutions appeared, and the calculations were converging to energy levels lower than it should be. It was observed by Schwarz and Wallmeier [2] as well as Grant [3] that in such calculations the matrix representation of the kinetic energy operator obtained in the non-relativistic limit of the Dirac equation did not match the Schrödinger one. It was realized that if the small components basis functions are generated from the large component ones by ϕ S i ∝ σ · pϕ L i , where σ are the Pauli spin matrices, then the non-relativistic limit of the kinetic energy operator goes directly to the Schrödinger one, and the spurious states disappear. This was further analyzed and formalized under the name of kinetic balance by Stanton and Havriliak [4] (see also Ref. 5). Calculations using this prescription were first done by Lee and McLean [6] (using unrestricted kinetic balance, see Section 2.3.3), and Ishikawa et al. [7].
Present-day 4-component relativistic atomic and molecular calculations are typically carried out within the no-pair approximation [8,9], where the electronic Hamiltonian is embedded by operators projecting out negative-energy orbitals, hence treating them as an orthogonal complement. However, going beyond the no-pair approximation and considering effects of quantum electrodynamics (QED), notably vacuum polarization and the self-energy of the electron, the negative-energy solutions take on physical reality and require a proper description. Charge conjugation symmetry also becomes an issue [10,11] and has to be considered when designing basis sets.
In the present work, we investigate the realization of charge conjugation symmetry, in short C-symmetry, of the one-electron Dirac equation within the finite basis approximation. Since basis functions are typically located at nuclear positions, we limit attention to the central field (spherically symmetric) problem. We shall consider three schemes for basis set construction: restricted kinetic balance [4,12], inverse kinetic balance [13], and dual kinetic balance [14]. As such our work bears some resemblance to the study of Sun et al. [13], but our focus will be on whether these schemes allow the realization of charge conjugation symmetry.

The Dirac equation and C-symmetry
The relativistic behavior of an electron placed in an electromagnetic potential (φ, A) is described by the Dirac equation where are the Dirac matrices, anti-commuting amongst themselves. Dirac himself noted that if matrices α y and β are swapped, then the complex conjugate of a solution to eq.(1) will be the solution of the same equation, but with opposite charge [15]. Kramers coined the term charge conjugation for this symmetry linking particles and their anti-particles [1], and it has later been elevated to one of the three fundamental symmetries of Nature through the CPT-theorem [10,[16][17][18].
In the Dirac representation, the C-operator is given by where K 0 is the complex conjugation operator. The general form was investigated by Pauli [19]. For static potentials, the solution of eq.(1) has the form ψ(r, Since the action of the charge conjugation operation is Cĥ −e C −1 = −ĥ +e , we get the time-independent positronic equation asĥ +e Cψ(r) = −ECψ(r), with opposite sign of the energy. In passing, we note that the charge conjugation operator can be expressed as C = γ 5 βK, where K is the time-reversal operator.

Radial problem
We shall limit attention to the central-field case, with the vector potential A = 0, and a radial scalar potential φ(r). Solutions then have the general form where the imaginary number i is introduced to make both radial functions P κ and Q κ real. The Ω κ,m are 2-component complex eigenfunctions of theκ = −h − σ ·ˆ operator [20], and κ represent the corresponding eigenvalue. After separation of radial and angular variables, we obtain the radial Dirac equationĥ The C-operation shown in eq.(3), when applied to the spherical solution, eq.(6), gives [21] Cψ κ,m j = iS κ (−1) m j + 1 We observe that the κ and m j quantum numbers in the angular parts have switched sign, and that the radial components are swapped P κ Q κ .

Free-particle radial problem
Usually, the free-particle Dirac equation solutions are presented in the form of plane waves. However, we are interested in atomic (and molecular) calculations where we use spherical basis functions centered at the nuclear positions. It is therefore more appropriate to consider the free-particle solutions in spherical symmetry. By setting φ(r) = 0 in the radial Dirac equation eq.(7), the solution of this problem is found to be [21,22] where the large and the small components (upper and lower) are respectively given by where j are the spherical Bessel functions of the first kind, S x = x |x| is the sign function, and k(E) = E 2 − m 2 e c 4 /ch represent the wavenumber. These solutions are normalized to the delta function as they describe a continuum of solutions. By next applying the C-operator to the free-particle solution, we get the C-conjugated partner Finally, we see that it is possible to connect opposite energy solutions, in spherical symmetry, by the C-operation (as we expect from the trivial Dirac plane wave case) Details about the C-symmetry and the free electron in the spherical case can be found in [21][22][23]].

Finite basis approximation
Generally, the plan is to specify a finite number of basis functions, construct the matrix representation of the Dirac equation, then diagonalize it to get the set of eigenfunctions and eigenvectors. We start by introducing radial basis sets for the large and the small components, , which means that the radial functions P κ and Q κ are expanded as giving the matrix representation of the Dirac equation as with Preprints (www.preprints.org) | NOT PEER-REVIEWED | Posted: 31 May 2020 doi:10.20944/preprints202005.0492.v1 The matrix elements of the above submatrices are given by We note that the κ appearing in superscripts refers to the radial basis functions, whereas the κ appearing as a subscript is associated with the operator. We also note that the −e term appearing in the subscripts denotes the charge that appears in front of the scalar potential. Since we expect a (real) Hermitian matrix representation, the off-diagonal matrices should be related by the transpose operation; Using integration by parts, this implies that that the basis functions should vanish at the boundaries; 0 and ∞.

Gaussian type functions
We shall work with Gaussian type functions since they play a central role in quantum chemical calculations. The large and the small component radial Gaussian functions are given by with N X κ i the normalization constants. We choose the exponents γ p and γ q to reproduce the small r behavior of the radial functions in the case of a finite nucleus [21], that is This also corresponds to the small r behavior of the free-particle radial solutions eqs. (10,11). Furthermore, Sun et al., in calculations on Rn 85+ using the dual kinetic balance prescription (discussed later in Section 2.3.4), investigated the use of different integers powers γ p and γ q of r for large and small Gaussian-type functions [13] and concluded that optimal results, in particular avoiding variational collapse and divergent integrals, were obtained using the powers given in eq.(23).

C-symmetry in the finite basis approximation
We say that a basis set respects C-symmetry, and thus leading to C-symmetric matrix representation, if the C-conjugation of all the elements of the basis set belongs to the basis set itself For simplicity reasons we shall set the phase factor in eq.(8) to 1, since it does not contribute to expectation values. C is a map of C 4 → C 4 , the last condition is equivalent to say that the subspace Φ Preprints (www.preprints.org) | NOT PEER-REVIEWED | Posted: 31 May 2020 doi:10.20944/preprints202005.0492.v1 of C 4 consisting of basis functions {ϕ 1 , ϕ 2 , ..., ϕ n }, is preserved by the C-map. We have seen before that in the spherically symmetric case, the C-operation replaces κ → −κ, π L π S and m j → −m j , which means that the realization of the C-symmetry at the basis set level implies Under these conditions, we find that and using the last equations, we get the connection between eigenvalues and eigenvectors +e,−κ = − −e,+κ ,

Kinetic balance
Starting from the radial equation eq. (7), we get two coupled equations that relates the large and small radial components of the Dirac equation The exact couplings are energy and potential-dependent and therefore not appropriate for the construction of basis sets prior to the calculation of the energy. The energy-dependence can be eliminated by taking non-relativistic limit c → ∞ . It is sometimes stated that the expressions in square brackets go to one provided E ± eφ (r) << m e c 2 . However, the correct statement is rather that E ± eφ (r) should have a finite value as the limit is taken. For a point nucleus the scalar potential φ (r) is singular at r = 0, and so this condition is not satisfied. It can be restored by rather considering nuclei of finite charge distributions [24,25]. As it stands, the energy depends quadratically on the speed of light. This dependence can be eliminated by constant shifts, but implies taking different limits for the positive-and negative energy branches.
For the positive energy branch, we introduce the shifted energy E + = E − m e c 2 and from eq.(29) obtain For the negative energy branch, we introduce the shifted energy E − = E + mc 2 and from eq.(30) obtain Conventional atomic and molecular relativistic calculations in a finite basis focus on positive-energy solutions and so bases are constructed according to the prescription of kinetic balance which imposes the non-relativistic coupling eq.(31) at the level of individual basis functions, that is Preprints (www.preprints.org) | NOT PEER-REVIEWED | Posted: 31 May 2020 doi:10.20944/preprints202005.0492.v1 The numerical factor a S κ i , here set to half the reduced Compton wavelength in accordance with eq.(31), is arbitrary. For instance, if one did not introduce an imaginary i in the atomic spinor, eq.(6), our present choice would be multiplied with this factor. This particular choice of basis functions provides a proper representation of the kinetic energy operator in the non-relativistic limit [4], and therefore prevents the appearance of spurious states in the calculation. For calculations at finite values of the speed of light, it is necessary that the basis is sufficiently flexible so that the relativistic coupling can be realized through adjustment of basis expansion coefficients [26].
The one-to-one correspondence between large and small component basis functions can only be realized in a 2-component basis and is denoted restricted kinetic balance (RKB). The terminology was introduced by Dyall and Faegri [12] to contrast with the use of scalar basis functions, where the small component basis functions are taken as derivatives of the large component ones, in no particular linear combinations. This latter scheme is denoted unrestricted kinetic balance (UKB).
Restricted kinetic balance (eq.(33)) leads to the matrix eigenvalue equation, eq. (16) with elements given in eqs.(A1,A2) of the Appendix. The matrix representation of the Dirac equation in an RKB basis is that of the modified Dirac equation [27] (see also Ref. 28). From this Sun et al. [13] conclude that there is no 'modified' Dirac equation. However, this is formally incorrect, since the modified Dirac equation has an independent existence at the operator level.
A corresponding prescription that favors negative-energy solutions has been termed inverse kinetic balance [13] (IKB) and is based on eq.(32). In this prescription the small component basis functions are introduced first, then the large ones are generated using the last equation. This prescription leads to the eigenvalue equation whose elements are given in the Appendix eqs. (A3,A4). In order to respect the C-symmetry we impose the conditions of eq.(25) on the RKB prescription eq.(33), which leads the following equation for large component basis functions Its general solution is and the small components are then π S κ i = π L −κ i . The c i are arbitrary coefficients, j α (z) and y α (z) are spherical Bessel functions of the first and second kind respectively. If, on the other hand, we impose C-symmetry on the IKB prescription eq.(34), we obtain the same general solution, now with We also note that if we combine the RKB and IKB prescriptions, to describe both positive and negative energy solutions on the same footing, we again get the same general solution, The problem with this specific choice of functions is that the boundary conditions π L κ /S κ i (0) = 0 and π L κ /S κ i (∞) = 0, are not obeyed simultaneously. Therefore they are not useful for atomic and molecular calculations.

Dual kinetic balance
The kinetic balance prescription is widely employed in atomic and molecular calculations, but favors the positive-energy solutions. The dual kinetic balance prescription (DKB) ensures the right coupling between the large and the small components (in the non-relativistic limit) for both positive and negative energy solutions. It was introduced by Shabaev et al. [14] with the use of B-splines and tested by calculating the one-loop self-energy correction for a hydrogen-like ion. The radial function is expanded as where the first and second set of basis functions have the non-relativistic coupling of positive-and negative-energy solutions, respectively, as indicated by the [±] symbol on the coefficients. This particular expansion leads to a generalized eigenvalue problem whose elements are defined in eqs.(A9-A11) of the Appendix. Contrary to the case of RKB/IKB the conditions for C-symmetry, eq.(25), can now be imposed without putting constraints on the choice of basis functions. The two matrix representations associated with (+e, −κ) and (−e, +κ), become related by leading to the C-connection between the eigenvalues, and the eigenvectors +e,−κ = − −e,+κ , Note, however, that the condition in eq.(25) that ensures the C-symmetry, implies that one has to use the same exponents for both ±κ Gaussian type functions, as has also been pointed out by Dyall [29]. This corresponds, in the terminology of Dyall to j bases, where exponents are optimized for the total angular momentum j quantum number [30], contrary to conventional basis sets where functions are optimized for orbital angular momentum .

Computational details
To illustrate our findings we have written numerical codes using the Wolfram Mathematica program [31]. We built the matrix representations of the Dirac equation in the RKB, IKB and DKB schemes, using spherical Gaussian functions, eqs. (21,22), and a point nucleus.

Kinetic balance
We started by doing a simple free-particle calculation, φ(r) = 0, within the RKB scheme. Using spherical Gaussian functions, we specify ζ κ = {1, 2}, with κ = ±1 (s 1 2 -and p 1 2 -type functions). By solving the generalized eigenvalue problem for each κ-block, we get the eigenvalues κ reported in Table 1. At first glance, one gets the impression that C-symmetry is respected since eigenvalues come in pairs of opposite sign. However, the pairs occur for same κ and not opposite κ as predicted by C-symmetry. This is confirmed by inspection of the eigenvectors, as exemplified by showing the first two normalized eigenvectors of each κ block in Table 2. We see clearly that the expected connection (C-conjugation) between positive and negative energy solutions, does not hold here. In order to understand the reason, we set φ(r) to zero in the RKB matrix equation whose elements are given eqs.(A1-A2) of the Appendix. We then get the following equation The first and second lines of the last equation give and by combining these equations, we get We see that each eigenvalue λ κ corresponds to two values κ = ± 2m e c 2 λ κ + m 2 e c 4 . The corresponding eigenvectors are Although the eigenvalues exist in pairs, it is clear that the upper and lower components of two opposite energy solutions are not related by C-symmetry. In fact, as shown in Section 2.3.3, the RKB prescription does not generally respect C-symmetry.
Note that this pairing of energies can already be seen from the exact spherical free-particle solutions in eq.(9). Upon substitution E → −E and keeping in mind that E ∈ R \ [−m e c 2 , +m e c 2 ] we see that the solution of flipped energy sign can be expressed in terms of the original one Doing the same calculation using the IBK prescription, we get the sets of eigenvalues shown in Table 3. The first two eigenvectors of each spectrum are shown in Table 4. By comparing the eight eigenvectors we have chosen in Table 2 and Table 4, we see that positive and negative energy solutions that belongs to opposite κ-sign blocks are related by C-symmetry. Taking into account the condition in eq.(25), we see that RKB and IKB matrices eqs.(A1-A4) in the Appendix are indeed connected by C-symmetry, that is leading to the symmetry between RKB and IKB eigensystems IKB +e,−κ = − RKB −e,+κ , Since as we see, RKB and IKB are related by C-symmetry, this means that a combination of the two prescriptions would conserve the C-symmetry. And this is exactly what the DKB is about (eq.(37)).

Dual kinetic balance
We present two simple atomic calculations within the DKB prescription, where we used s For each calculation we get two sets of eigenvalues coming from each κ-block, shown in Table 5. Then we pick the first eigenvalues of each set and show their corresponding eigenvectors in Table 6.

Conclusion
We have investigated three basis set schemes for solving the Dirac equation in a central field: restricted, inverse and dual kinetic balance, and their compatibility with charge conjugation symmetry which connects solutions of opposite κ of the electronic and positronic problem. An interesting observation is that in the free-particle case, where the electronic and positronic problem coalesce, there is further symmetry such that pairs of eigenvalues of opposite sign also occur for same κ. We are not aware of any discussion of this feature in the literature.
Charge conjugation symmetry can be realized within restricted and inverse kinetic balance, but only using special functions which do not respect the boundary conditions of the radial Dirac solutions and which are not useful for atomic and molecular calculations. Dual kinetic balance, on the other hand, is compatible with charge conjugation symmetry for any type of radial basis function, provided j bases are used.
An alternative to dual kinetic balance, denoted dual atomic balance, has been proposed by Dyall [29]: In this scheme restricted and inverse kinetic balance is used separately for positive-and negative-energy solutions. This requires in principle two separate diagonalizations, followed by a final diagonalization in the dual basis. If one seeks to generate orbitals for use in QED calculations, then a possible simple alternative to the latter scheme is to first carry out a standard 4-component relativistic calculation within restricted kinetic balance and electronic charge q = −e and only retain the positive-energy solutions. Then a second calculation is carried out, again within restricted kinetic balance, retaining only positive-energy solutions and with the same potential, but now using the positronic charge q = +e. This scheme then has the intriguing property that the final set of orbitals is restricted to observable, positive-energy solutions only. We plan to study these schemes in future work.