[1]

[1]This document is the results of the research project PANDORA funded by the
National Institute of Nuclear Physics. S.T is also funded by the Bruno Kessler
Foundation.

[orcid=0000-0003-4010-8000] [1]

[orcid=0000-0001-9287-8651]

[2]

[cor1]Principal Corresponding author [cor2]Corresponding author

# Relativistic quantum theory and algorithms: a toolbox for modeling many-
fermion systems in different scenarios

Simone Taioli<EMAIL_ADDRESS>European Centre for Theoretical Studies in
Nuclear Physics and Related Areas (ECT*-FBK), Trento, Italy  Trento Institute
for Fundamental Physics and Applications (TIFPA-INFN), Trento, Italy  Peter
the Great St. Petersburg Polytechnic University, Russia  Stefano Simonucci
School of Science and Technology, University of Camerino, Italy  INFN, Sezione
di Perugia, Italy<EMAIL_ADDRESS>
###### Abstract

In this chapter we focus first on the theoretical methods and relevant
computational approaches to calculate the electronic structure of atoms,
molecules, and clusters containing heavy elements for which relativistic
effects become significant. In particular, we discuss the mean-field
approximation of the Dirac equation for many-electron systems, and its self-
consistent numerical solution by using either radial mesh or Gaussian basis
sets. The former technique is appropriate for spherical symmetric problems,
such as atoms, while the latter approach is better suited to study non-
spherical non-periodic polycentric systems, such as molecules and clusters. We
also outline the pseudopotential approximation in relativistic context to deal
with the electron-ion interaction in extended systems, where the unfavourable
computational scaling with system size makes it necessary. As test cases we
apply our theoretical and numerical schemes to the calculation of the
electronic structure i) of the gold atom, and ii) of the superatom W@Au12,
where the inclusion of spin-orbit effects is crucial to the accurate
understanding of the electronic properties. Furthermore, we describe the
extension of our relativistic approach to deal with nuclear reactions driven
by the weak force, such as the electron capture and $\beta$-decay, also at
finite temperature in astrophysical scenarios, using the Fermi–Dirac
statistics. The latter processes are indeed major drivers of the
nucleosynthesis of the elements in stars and, thus, their understanding is
crucial to model the chemical evolution of the Universe. Finally, we show the
application of our relativistic quantum mechanical framework to the assessment
of the elastic differential scattering cross section of electrons impinging on
molecular targets, notably liquid water. The latter process, together with
several inelastic scattering collisions by which secondary electrons deposit
their energy, represents a fundamental event of the chain of the physico-
chemical mechanisms initiated by the passage of fast ion beams through a bio-
medium. This technique is used in hadrontherapy for cancer cure.

###### keywords:

Relativistic quantum mechanics Computational methods Gaussian functions Real-
space grids Gold atom and W@Au12 $\beta$-decay Stars Cancer cure Electron-
molecule elastic scattering

##  1 Introduction

The relativistic extension of quantum mechanics to deal with many-particle
systems, such as nuclei, atoms, molecules and clusters, represents the most
fundamental theory of all molecular sciences [1]. Indeed, it combines the laws
of special relativity, to which natural systems obey, with the realm of
quantum mechanics.
Nonetheless, the description of the motion of valence electrons, which drive
the chemistry of materials, typically relies on a phenomenological approach
whereby the Schrödinger equation is solved with the possible addition of a
spin-orbit coupling term to the effective quantum-mechanical Hamiltonian for a
system of electrons. The latter term corrects the quantum-mechanical theory
for relativistic effects e.g. due to the presence of heavy elements or large
magnetic moments leading to significant spin-orbit interactions in the system.
However, this phenomenological approach, whose attractiveness lays essentially
in its affordable computational cost with respect to the fully relativistic
quantum electrodynamics methods, fails to deliver accurate results comparable
to experiments, most notably in heavy element compounds where electrons attain
sufficient speeds that the four-dimensional nature of the basic constituents
of matter cannot be neglected. In fact, the theory of light interacting with
many-electrons systems must adopt a rigorous gauge invariant and Lorentz
covariant form to describe both electromagnetic and fermionic degrees of
freedom and their interaction simultaneously [1].
This chapter aims to describe the theoretical and computational methods to
simulate the electronic structure of many-electron systems for which
relativistic effects are significant. Specifically, we present advanced
algorithms to find the numerical solution of the relativistic quantum-
mechanical Dirac equation using a mean field approximation. This theoretical
and computational framework, which is called the Dirac–Hartree–Fock (DHF)
method, is implemented either using real-space grids or by introducing
functional spaces of square-integrable ($L^{2}$) functions.
The former numerical method consists in the discretization in space by using a
grid of points to sample the (electron or nucleon) wavefunctions to a desired
accuracy. While grids can be used for problems with any number of dimensions,
they are typically preferred for dealing with spherical symmetric systems as
they reduce to the simple one-dimensional case upon factorizing out the
wavefunction angular dependence. The big advantages of real-space grid
calculations are their simplicity, versatility, and the possibility to
systematically improve the calculation accuracy by increasing the mesh density
as the error is proportional to a power law of the step size (defined as the
distance between neighboring points of the grid). Furthermore, boundary
conditions, particularly for solving the DHF equation for scattering states,
are easily imposed. Finally, the use of pseudopotentials to limit the explicit
treatment of the electron-electron and electron-ion interaction to outer
electrons only, which we will describe in the context of the Hermite Gaussian
Basis Functions [2, 3, 4] (HGBF), can be also efficiently implemented in grid-
space methods.
The grid-based approach can be considered a limiting case (for infinitesimally
small Gaussian widths) of the second numerical method that we describe here,
which relies on multi-centered HGBF – chosen as prototypes of $L^{2}$
functions – spanning the linear vector space $\cal G$ of bound and scattering
states [2, 3, 4, 5, 6]. Despite the solution cannot be systematically
improved, at variance with space grid methods (or plane-waves), due to Basis
Set Superimposition (or over-linear) errors [7], an advantage of using HGBF is
that it is possible to derive explicit analytical expressions of the integrals
necessary to carry out atomic and molecular calculations.
Applications will be focused on non periodic many-body systems interacting via
the Coulomb force, thus mainly atoms, molecules, and clusters. However, we
will show how our theoretical and computational methods are general enough to
deal with any (pseudo-) scalar (e.g. the pion field) and any (pseudo-) vector
potential (e.g. the Yukawa interaction between nucleons). In this regard, we
will use our approach also in problems typical of the nuclear physics, such as
the theory of beta decay of unstable nuclei induced by the electroweak force
(for which we refer the reader to Ref. [8]). This approach has been
implemented in a simulation package called DIRECT (DIRac Equation from
Camerino and Trento).
In particular, this fully relativistic quantum-mechanical framework will be
tested first to carry out the all-electron calculation of the electronic
structure of the gold ion (Au+). We compare the numerical results of some
electronic structure properties, such as the ionization potential (à la
Koopmans), by using either radial or HGBF basis sets. The chemistry of gold is
indeed shaped by strong relativistic effects, whereby the radius orbital of
the core electrons shrinks by more than 20%.
Furthermore, as a second test case we will discuss the W@Au12 superatomic
system, where the inclusion of spin-orbit effects is necessary to describe the
spectroscopic features. The 13-atom W@Au12 cluster has been experimentally
observed as an highly stable icosahedral cage of 12 gold atoms containing a
central tungsten, and characterised by photoelectron spectroscopy (PES) as a
closed-shell electron configuration with 18 valence electrons and an HOMO–LUMO
gap around 3 eV [9, 10]. This makes it chemically inert. Gold-based
nanoclusters of different size and shape have in general attracted increasing
interest, primarily due to their unique and tunable electronic and optical
properties. They can find application in the development of photovoltaic cells
or batteries, of biological and chemical sensors, or could replace toxic or
rare elements in biomedical devices, being biocompatible and chemically inert.
More recently, the use in radiation therapy for cancer cure of high-Z
materials, such as gold, platinum, gadolinium, and iron nanoparticles (NPs)
with a size that ranges between 1 to 100 nm, has been investigated and shown
to increase secondary electrons and free radical production [11, 12, 13, 14,
15, 16, 17, 18], eventually enhancing the Relative Biological Effectiveness
(RBE) of hadrontherapy. By selectively accumulating near the tumour region
gold NPs, the RBE can be potentially enhanced by 15–20% due to an increase in
both direct and indirect damage. However, while some of these materials can be
toxic for humans, gold NPs might offer a viable biocompatible alternative to
enhance localized dose and biodamage under photon [19] and proton irradiation
[20], once an adequate concentration of NPs in the tumour is achieved. In this
context, we will describe also the recent implementation of a fully
relativistic theory of the pseudopotential approach, which is a necessary
approximation to achieve an efficient scaling in the electronic structure
calculations of heavy-element compounds.
A second application of our fully relativistic framework, which goes beyond
the treatment of bound states, concerns the interpretation of beta-decay
spectra of heavy nuclei in both Earth conditions [8] and astrophysical [21]
scenarios.
The use of radioisotopes that decay via beta-emission is widespread in science
and medicine. On the one hand, $\beta$-decay processes are the dominant source
of background in liquid Xenon dark matter detectors [22] or in cerium-based
LaBr3-based crystals that are used to fabricate scintillators for applications
in gamma-ray spectroscopy and, more generally, to characterise ionizing
radiation [8]. $\beta$-emitting radiopharmaceuticals are widely used also in
clinical oncology to diagnose cancer and to deliver radiation therapy in
cancer treatment, or in the field of PET imaging [23].
On the other hand, stars have been known since decades to be the cauldrons
where heavy chemical elements are forged and matter becomes complex starting
from building blocks such as H, He and other very light elements ($A<8$)
produced by Big Bang nucleosynthesis (BBN) in the early stage of the Universe
[24, 25, 26]. This understanding relies on the discovery that some old stars
host isotopes, whose lifetime requires in situ formation. Theoretical models
and observations thus support the conclusion that reactions occurring within
stars are responsible for heavy elements formation. Nevertheless, the solution
to the global problem of studying the chemical evolution of the Universe, from
the evolutionary path after core He burning to the mechanisms driving core-
collapse supernovae and gamma ray bursts, is still missing and we are only now
obtaining first rough matches between models and observations [27]. Indeed,
while robust confirmation of the validity of the solar model can be derived by
the analysis of solar neutrinos, helio-seismology and improved experimental
nuclear physics data [28], a thorough understanding of the processes of
nucleosynthesis occurring at advanced stages of star evolution is uncertain
[29], even for masses low enough to experience only H- and He-burning phases
(Low and Intermediate Mass Stars, with $M<6-8$ solar masses). This could be
also the case of the cosmological Li abundance, which is still unexplained
within the framework of the Big Bang Standard Model, despite almost all its
predictions on isotopic abundances are confirmed [21, 30].
In fact, the problem of determining an ample selection of stellar elemental
abundances, the isotopic composition of meteorites and of the interstellar
medium, is cumbersome and its solution inherently needs an interdisciplinary
look. Indeed, this description must integrate a variety of approaches, ranging
from the description of nuclear reaction cross sections to the fine tuning of
evolutionary models. In particular, a thorough understanding of the
nucleosynthesis in stars must include the accurate determination of
$\beta$-decay and electron capture rates (other than slow and rapid neutron
capture processes) in astrophysical scenarios. Here we discuss the application
of our fully-relativistic quantum method for calculating weak interaction
processes in nuclear astrophysics context (typically occurring in post-main
sequence long-lived stages of stars). These methodological advances to
determine reaction rates will lead to a more quantitative-based knowledge of
the advanced evolutionary stages in massive stars and core-collapse
supernovae.
A last example to demonstrate the diverse applicability of our relativistic
quantum mechanical scheme will be the determination of the elastic mean free
path (EMFP) of electrons moving within a medium. This approach, based on
scattering theory, will be used in particular to assess the electron EMFP in
liquid water. The latter study is relevant for modelling the biodamage induced
by ion irradiation in materials of biological relevance [31, 32], such as the
human body tissue typically mimicked by liquid water [33]. Indeed, accurate ab
initio calculations of elastic events, together with the assessment of several
other chemical-physical mechanisms of electron energy loss (ionization,
excitation, dissociative electron attachment [34, 35]) initiated by ion
irradiation, make it possible to obtain reliable probabilities to assess
clustered DNA damage [33].
This chapter is organized as follows: section 2 outlines the theoretical
foundations of our relativistic quantum chemical method; section III deals
with the computational algorithms devised to find the numerical solution of
the DHF equation using either a space-grid or HBGF basis sets; sections IV, V,
and VI are devoted to the application of our method to the calculations of the
electronic properties of Au+ and W@Au12, to the simulation of $\beta-$decay of
31Si, 89Sr and 134Cs, and to the electron EMFP calculation in liquid water,
respectively.

##  2 Wavefunction calculation within the mean field approximation

###  2.1 The Dirac-Hartree-Fock method

The many-particle Dirac equation reads [1] (everywhere we use atomic units
otherwise stated):

$\left[\sum_{i}\left(c{\bf{\alpha}}_{i}\cdot\mathbf{p}_{i}+\beta_{i}mc^{2}+V_{i}\right)+\sum_{i<j}\left(1-\mathbf{\alpha}_{i}\cdot\mathbf{\alpha}_{j}\right)g_{ij}\right]\psi\left(\mathbf{r}_{1},\cdots\mathbf{r}_{N}\right)=E\psi\left(\mathbf{r}_{1},\cdots\mathbf{r}_{N}\right)$
(1)

where $\alpha_{i}$ and $\beta$ are the following Hermitian involutory 4
$\times$ 4 matrices

$\displaystyle\mathbf{\alpha}=\left(\begin{array}[]{cc}0&\mathbf{\sigma}\\\
\mathbf{\sigma}&0\end{array}\right),\qquad\beta=\left(\begin{array}[]{cc}I&0\\\
0&-I\end{array}\right)$ (6)

that mutually anticommute, $\sigma_{i}$ are the Pauli matrices, and $I$ is the
2$\times$2 identity matrix. In equation (1) $V_{i}$ is an external potential,
which may represent a general interaction between fermions.

When both scalar ($g_{S}$) and vector ($g_{V}$) potentials are in place, the
Dirac equation (1) can be written [8]:

$\left\\{\sum_{i}\left(c\mathbf{\alpha}_{i}\cdot\mathbf{p}_{i}+\beta_{i}mc^{2}+V_{i}\right)+\sum_{i<j}\left[\beta_{i}\beta_{j}g_{S,ij}+\left(1-\mathbf{\alpha}_{i}\cdot\mathbf{\alpha}_{j}\right)g_{V,ij}\right]\right\\}\psi\left(\mathbf{r}_{1},\cdots\mathbf{r}_{N}\right)=E\psi\left(\mathbf{r}_{1},\cdots\mathbf{r}_{N}\right)$
(7)

It is convenient to write the Dirac Hamiltonian in second quantization, as
fermionic anticommutation rules are therein automatically included in the
formalism. The Dirac equation (7) reads:

$\displaystyle H=\sum_{s_{1}s_{2}}\int
d\mathbf{r}\hat{\psi}_{s_{1}}^{+}(\mathbf{r})\left[-ic\mathbf{\alpha}_{s_{1}s_{2}}\cdot\mathbf{\nabla}+\beta_{s_{1}s_{2}}mc^{2}+\delta_{s_{1}s_{2}}V(\mathbf{r})\right]\hat{\psi}_{s_{2}}(\mathbf{r})+$
$\displaystyle\frac{1}{2}\sum_{s_{1}s_{2}s_{1}^{\prime}s_{2}^{\prime}}\int
d\mathbf{r}d\mathbf{r}^{\prime}\
\hat{\psi}_{s_{1}}^{+}(\mathbf{r})\hat{\psi}_{s_{1}^{\prime}}^{+}(\mathbf{r}^{\prime})\times\left[\beta_{s_{1}s_{2}}\beta_{s_{1}^{\prime}s_{2}^{\prime}}g_{S}\left(\mathbf{r},\mathbf{r}^{\prime}\right)+\left(\delta_{s_{1}s_{2}}\delta_{s_{1}^{\prime}s_{2}^{\prime}}-\mathbf{\alpha}_{s_{1}s_{2}}\cdot\mathbf{\alpha}_{s_{1}^{\prime}s_{2}^{\prime}}^{\prime}\right)g_{V}\left(\mathbf{r},\mathbf{r}^{\prime}\right)\right]\hat{\psi}_{s_{2}^{\prime}}(\mathbf{r}^{\prime})\hat{\psi}_{s_{2}}(\mathbf{r})$

where $s_{1},s_{2},s_{1}^{\prime},s_{2}^{\prime}$ label the bispinor upper and
lower components.
The Hartree-Fock (HF) approximation to the exact Dirac Hamiltonian (2.1) is
introduced by assuming that:

$\left\langle\hat{\psi}_{s_{1}}^{+}(\mathbf{r})\hat{\psi}_{s_{1}^{\prime}}^{+}(\mathbf{r}^{\prime})\hat{\psi}_{s_{2}^{\prime}}(\mathbf{r}^{\prime})\hat{\psi}_{s_{2}}(\mathbf{r})\right\rangle=\left\langle\hat{\psi}_{s_{1}}^{+}(\mathbf{r})\hat{\psi}_{s_{2}}(\mathbf{r})\right\rangle\left\langle\hat{\psi}_{s_{1}^{\prime}}^{+}(\mathbf{r}^{\prime})\hat{\psi}_{s_{2}^{\prime}}(\mathbf{r}^{\prime})\right\rangle-\left\langle\hat{\psi}_{s_{1}}^{+}(\mathbf{r})\hat{\psi}_{s_{2}^{\prime}}(\mathbf{r}^{\prime})\right\rangle\left\langle\hat{\psi}_{s_{1}^{\prime}}^{+}(\mathbf{r}^{\prime})\hat{\psi}_{s_{2}}(\mathbf{r})\right\rangle$
(9)

By defining the ($4\times 4$) density matrix

$\rho_{s_{2}^{\prime}s_{1}}\left(\mathbf{r}^{\prime},\mathbf{r}\right)=\left\langle\hat{\psi}_{s_{1}}^{+}(\mathbf{r})\hat{\psi}_{s_{2}^{\prime}}(\mathbf{r}^{\prime})\right\rangle$
(10)

the DHF Hamiltonian (2.1) in second quantization can be written:

$\displaystyle H_{DHF}=\sum_{s_{1}s_{2}}\int
d\mathbf{r}\hat{\psi}_{s_{1}}^{+}(\mathbf{r})\left\\{-\mathbf{\alpha}_{s_{1}s_{2}}\cdot\left[ic\mathbf{\nabla}+\mathbf{A}_{H}\left(\mathbf{r}\right)\right]+\beta_{s_{1}s_{2}}\left[mc^{2}+V_{HS}\left(\mathbf{r}\right)\right]+\delta_{s_{1}s_{2}}\left[V(\mathbf{r})+V_{H}(\mathbf{r})\right]\right\\}\hat{\psi}_{s_{2}}(\mathbf{r})$
$\displaystyle-\sum_{s_{1}^{\prime}s_{2}}\int
d\mathbf{r}d\mathbf{r}^{\prime}\hat{\psi}_{s_{1}^{\prime}}^{+}(\mathbf{r}^{\prime})\left[V_{FS,s_{1}^{\prime}s_{2}}\left(\mathbf{r}^{\prime},\mathbf{r}\right)-A_{F,s_{1}^{\prime}s_{2}}\left(\mathbf{r}^{\prime},\mathbf{r}\right)+V_{F,s_{1}^{\prime}s_{2}}\left(\mathbf{r}^{\prime},\mathbf{r}\right)\right]\hat{\psi}_{s_{2}}(\mathbf{r})$
(11)

where

$\displaystyle\mathbf{A}_{H}\left(\mathbf{r}\right)$ $\displaystyle=\int
d\mathbf{r}^{\prime}\sum_{s_{1}^{\prime}s_{2}^{\prime}}\left[\rho_{s_{2}^{\prime}s_{1}^{\prime}}\left(\mathbf{r}^{\prime},\mathbf{r}^{\prime}\right)\mathbf{\alpha}_{s_{1}^{\prime}s_{2}^{\prime}}\right]g_{V}\left(\mathbf{r},\mathbf{r}^{\prime}\right)$
(12) $\displaystyle V_{H}\left(\mathbf{r}\right)$ $\displaystyle=\int
d\mathbf{r}^{\prime}\sum_{s_{1}^{\prime}s_{2}^{\prime}}\left[\rho_{s_{2}^{\prime}s_{1}^{\prime}}\left(\mathbf{r}^{\prime},\mathbf{r}^{\prime}\right)\delta_{s_{1}^{\prime}s_{2}^{\prime}}\right]g_{V}\left(\mathbf{r},\mathbf{r}^{\prime}\right)$
(13) $\displaystyle V_{HS}\left(\mathbf{r}\right)$ $\displaystyle=\int
d\mathbf{r}^{\prime}\sum_{s_{1}^{\prime}s_{2}^{\prime}}\left[\rho_{s_{2}^{\prime}s_{1}^{\prime}}\left(\mathbf{r}^{\prime},\mathbf{r}^{\prime}\right)\beta_{s_{1}^{\prime}s_{2}^{\prime}}\right]g_{S}\left(\mathbf{r},\mathbf{r}^{\prime}\right)$
(14)

and

$\displaystyle
A_{F,s_{1}^{\prime}s_{2}}\left(\mathbf{r}^{\prime},\mathbf{r}\right)$
$\displaystyle=\sum_{s_{1}s_{2}^{\prime}}\mathbf{\alpha}_{s_{1}s_{2}}\cdot\mathbf{\alpha}_{s_{1}^{\prime}s_{2}^{\prime}}\rho_{s_{2}^{\prime}s_{1}}\left(\mathbf{r},\mathbf{r}^{\prime}\right)g_{V}\left(\mathbf{r},\mathbf{r}^{\prime}\right)$
(15) $\displaystyle
V_{F,s_{1}^{\prime}s_{2}}\left(\mathbf{r}^{\prime},\mathbf{r}\right)$
$\displaystyle=\sum_{s_{1}s_{2}^{\prime}}\delta_{s_{1}s_{2}}\delta_{s_{1}^{\prime}s_{2}^{\prime}}\rho_{s_{2}^{\prime}s_{1}}\left(\mathbf{r},\mathbf{r}^{\prime}\right)g_{V}\left(\mathbf{r},\mathbf{r}^{\prime}\right)$
(16) $\displaystyle
V_{FS,s_{1}^{\prime}s_{2}}\left(\mathbf{r}^{\prime},\mathbf{r}\right)$
$\displaystyle=\sum_{s_{1}s_{2}^{\prime}}\beta_{s_{1}s_{2}}\beta_{s_{1}^{\prime}s_{2}^{\prime}}\rho_{s_{2}^{\prime}s_{1}}\left(\mathbf{r},\mathbf{r}^{\prime}\right)g_{S}\left(\mathbf{r},\mathbf{r}^{\prime}\right)$
(17)

The sums over $s_{1}$ and $s_{2}^{\prime}$ in equations (15-17) must be
understood as matrix product. By writing the density matrix as a $2\times 2$
bloc matrix ($L~{}(S)$ labels the large (small) component of the Dirac
bispinor)

$\rho\left(\mathbf{r}^{\prime},\mathbf{r}\right)=\left(\begin{array}[]{cc}\rho_{LL}\left(\mathbf{r}^{\prime},\mathbf{r}\right)&\rho_{LS}\left(\mathbf{r}^{\prime},\mathbf{r}\right)\\\
\rho_{SL}\left(\mathbf{r}^{\prime},\mathbf{r}\right)&\rho_{SS}\left(\mathbf{r}^{\prime},\mathbf{r}\right)\end{array}\right)$
(18)

Equations (12-17) can be written explicitly:

$\displaystyle\mathbf{A}_{H}(\mathbf{r})$ $\displaystyle=\int
d\mathbf{r}^{\prime}\mathrm{Tr}\left[\left(\begin{array}[]{cc}\rho_{LL}\left(\mathbf{r}^{\prime},\mathbf{r}^{\prime}\right)&\rho_{LS}\left(\mathbf{r}^{\prime},\mathbf{r}^{\prime}\right)\\\
\rho_{SL}\left(\mathbf{r}^{\prime},\mathbf{r}^{\prime}\right)&\rho_{SS}\left(\mathbf{r}^{\prime},\mathbf{r}^{\prime}\right)\end{array}\right)\left(\begin{array}[]{cc}0&\mathbf{\sigma}\\\
\mathbf{\sigma}&0\end{array}\right)\right]g_{V}\left(\mathbf{r},\mathbf{r}^{\prime}\right)$
(22) $\displaystyle=\int
d\mathbf{r}^{\prime}\mathrm{Tr}\left\\{\left[\rho_{LS}\left(\mathbf{r}^{\prime},\mathbf{r}^{\prime}\right)+\rho_{SL}\left(\mathbf{r}^{\prime},\mathbf{r}^{\prime}\right)\right]\mathbf{\sigma}\right\\}g_{V}\left(\mathbf{r},\mathbf{r}^{\prime}\right)$
(23)

$\displaystyle V_{H}(\mathbf{r})$ $\displaystyle=\int
d\mathbf{r}^{\prime}\mathrm{Tr}\left[\left(\begin{array}[]{cc}\rho_{LL}\left(\mathbf{r}^{\prime},\mathbf{r}^{\prime}\right)&\rho_{LS}\left(\mathbf{r}^{\prime},\mathbf{r}^{\prime}\right)\\\
\rho_{SL}\left(\mathbf{r}^{\prime},\mathbf{r}^{\prime}\right)&\rho_{SS}\left(\mathbf{r}^{\prime},\mathbf{r}^{\prime}\right)\end{array}\right)\right]g_{V}\left(\mathbf{r},\mathbf{r}^{\prime}\right)$
(25) $\displaystyle=\int
d\mathbf{r}^{\prime}\mathrm{Tr}\left[\rho_{LL}\left(\mathbf{r}^{\prime},\mathbf{r}^{\prime}\right)+\rho_{SS}\left(\mathbf{r}^{\prime},\mathbf{r}^{\prime}\right)\right]g_{V}\left(\mathbf{r},\mathbf{r}^{\prime}\right)$
(26)

$\displaystyle V_{HS}\left(\mathbf{r}\right)$ $\displaystyle=\int
d\mathbf{r}^{\prime}\mathrm{Tr}\left[\left(\begin{array}[]{cc}\rho_{LL}\left(\mathbf{r}^{\prime},\mathbf{r}^{\prime}\right)&\rho_{LS}\left(\mathbf{r}^{\prime},\mathbf{r}^{\prime}\right)\\\
\rho_{SL}\left(\mathbf{r}^{\prime},\mathbf{r}^{\prime}\right)&\rho_{SS}\left(\mathbf{r}^{\prime},\mathbf{r}^{\prime}\right)\end{array}\right)\left(\begin{array}[]{cc}1&0\\\
0&-1\end{array}\right)\right]g_{S}\left(\mathbf{r},\mathbf{r}^{\prime}\right)$
(30) $\displaystyle=\int
d\mathbf{r}^{\prime}\mathrm{Tr}\left[\rho_{LL}\left(\mathbf{r}^{\prime},\mathbf{r}^{\prime}\right)-\rho_{SS}\left(\mathbf{r}^{\prime},\mathbf{r}^{\prime}\right)\right]g_{S}\left(\mathbf{r},\mathbf{r}^{\prime}\right)$
(31)

$\displaystyle A_{F}\left(\mathbf{r}^{\prime},\mathbf{r}\right)$
$\displaystyle=\left(\begin{array}[]{cc}0&\mathbf{\sigma}\\\
\mathbf{\sigma}&0\end{array}\right)\cdot\left(\begin{array}[]{cc}\rho_{LL}\left(\mathbf{r}^{\prime},\mathbf{r}\right)&\rho_{LS}\left(\mathbf{r}^{\prime},\mathbf{r}\right)\\\
\rho_{SL}\left(\mathbf{r}^{\prime},\mathbf{r}\right)&\rho_{SS}\left(\mathbf{r}^{\prime},\mathbf{r}\right)\end{array}\right)\cdot\left(\begin{array}[]{cc}0&\mathbf{\sigma}\\\
\mathbf{\sigma}&0\end{array}\right)g_{V}\left(\mathbf{r}^{\prime},\mathbf{r}\right)$
(38)
$\displaystyle=\left(\begin{array}[]{cc}\mathbf{\sigma}\cdot\rho_{SS}\left(\mathbf{r}^{\prime},\mathbf{r}\right)\cdot\mathbf{\sigma}&\mathbf{\sigma}\cdot\rho_{LS}\left(\mathbf{r}^{\prime},\mathbf{r}\right)\cdot\mathbf{\sigma}\\\
\mathbf{\sigma}\cdot\rho_{SL}\left(\mathbf{r}^{\prime},\mathbf{r}\right)\cdot\mathbf{\sigma}&\mathbf{\sigma}\cdot\rho_{LL}\left(\mathbf{r}^{\prime},\mathbf{r}\right)\cdot\mathbf{\sigma}\end{array}\right)g_{V}\left(\mathbf{r}^{\prime},\mathbf{r}\right)$
(41)

$V_{F}\left(\mathbf{r}^{\prime},\mathbf{r}\right)=\left(\begin{array}[]{cc}\rho_{LL}\left(\mathbf{r}^{\prime},\mathbf{r}\right)&\rho_{SL}\left(\mathbf{r}^{\prime},\mathbf{r}\right)\\\
\rho_{LS}\left(\mathbf{r}^{\prime},\mathbf{r}\right)&\rho_{SS}\left(\mathbf{r}^{\prime},\mathbf{r}\right)\end{array}\right)g_{V}\left(\mathbf{r}^{\prime},\mathbf{r}\right)$
(42)

$\displaystyle V_{FS}\left(\mathbf{r}^{\prime},\mathbf{r}\right)$
$\displaystyle=\left(\begin{array}[]{cc}1&0\\\
0&-1\end{array}\right)\cdot\left(\begin{array}[]{cc}\rho_{LL}\left(\mathbf{r}^{\prime},\mathbf{r}\right)&\rho_{LS}\left(\mathbf{r}^{\prime},\mathbf{r}\right)\\\
\rho_{SL}\left(\mathbf{r}^{\prime},\mathbf{r}\right)&\rho_{SS}\left(\mathbf{r}^{\prime},\mathbf{r}\right)\end{array}\right)\cdot\left(\begin{array}[]{cc}1&0\\\
0&-1\end{array}\right)g_{S}\left(\mathbf{r}^{\prime},\mathbf{r}\right)$ (49)
$\displaystyle=\left(\begin{array}[]{cc}\rho_{LL}\left(\mathbf{r}^{\prime},\mathbf{r}\right)&-\rho_{SL}\left(\mathbf{r}^{\prime},\mathbf{r}\right)\\\
-\rho_{LS}\left(\mathbf{r}^{\prime},\mathbf{r}\right)&\rho_{SS}\left(\mathbf{r}^{\prime},\mathbf{r}\right)\end{array}\right)g_{S}\left(\mathbf{r}^{\prime},\mathbf{r}\right)$
(52)

With these definitions in place, the most general Dirac equation (2.1)
including both scalar and vector potentials can be written in matrix form as
follows:

$\left(\begin{array}[]{ll}mc^{2}+W_{V}+W_{S}+\mathbf{A}_{P}\cdot\mathbf{\sigma}-E&-c\mathbf{\sigma}\cdot
i\mathbf{\nabla}-\mathbf{\sigma}\cdot\mathbf{A}+W_{PS}\\\
-c\mathbf{\sigma}\cdot
i\mathbf{\nabla}-\mathbf{\sigma}\cdot\mathbf{A}+W_{PS}&-mc^{2}+W_{V}+\mathbf{A}_{P}\cdot\mathbf{\sigma}-W_{S}-E\end{array}\right)\left(\begin{array}[]{l}\psi_{L}\\\
\psi_{S}\end{array}\right)=0$ (53)

where $W_{S},W_{PS}$ are generic scalar and pseudo-scalar potentials, while
$W_{V},\mathbf{A}_{P}$ are vector and pseudo-vector potentials. We notice that
the scalar and vector potentials sum up ($W_{S}+W_{V}$) in the head of the
Hamiltonian matrix; they act on the upper part of the spinor (the particle).
At variance, the scalar and vector potentials are subtracted ($W_{V}-W_{S}$)
in the tail of the Hamiltonian matrix; they act on the lower part of the
spinor (the anti-particle): this element thus accounts for the relativistic
corrections, that is the spin-orbit interaction.

###  2.2 Rotational invariant systems

The DHF equation can be written specifically to deal with spherical symmetric
systems for which both the electronic density $\rho(\mathbf{r})=\rho(r)$ and
the HF potential $V(\mathbf{r})=V(r)$ are rotational invariant and depend only
on the radial coordinate modulus.
Furthermore, one has:

$g_{V,S}\left(\mathbf{r},\mathbf{r}^{\prime}\right)=g_{V,S}\left(\left|\mathbf{r}-\mathbf{r}^{\prime}\right|\right)=\sum_{l=0}^{\infty}\frac{4\pi}{2l+1}g_{V,S,l}\left(r_{<},r_{>}\right)\sum_{m=-l}^{l}Y_{lm}^{*}\left(\vartheta,\varphi\right)Y_{lm}\left(\vartheta^{\prime},\varphi^{\prime}\right)$
(54)

where $r_{<}=\mathrm{min}(r,r^{\prime})$, $r_{>}=\mathrm{max}(r,r^{\prime})$,
and $Y_{lm}$ are spherical harmonics.
The solution to the DHF equation in spherical symmetry can be written as:

$\displaystyle H\psi_{\alpha,\kappa,m}(\mathbf{r})$ $\displaystyle=$
$\displaystyle E_{\alpha,\kappa}\psi_{\alpha,\kappa,m}(\mathbf{r})$ (55)
$\displaystyle\psi_{\alpha,\kappa,m}(\mathbf{r})$ $\displaystyle=$
$\displaystyle\left(\begin{array}[]{c}\psi_{\alpha,\kappa,L}\left(r\right)\chi_{\kappa,m}\left(\vartheta,\varphi\right)\\\
\psi_{\alpha,\kappa,S}\left(r\right)\chi_{-\kappa,m}\left(\vartheta,\varphi\right)\end{array}\right)$
(58)

where the angular eigenfunctions fulfill the following completeness relation

$\sum_{m}\frac{4\pi}{2l+1}\int\chi_{\kappa,m}^{+}\left(\vartheta^{\prime},\varphi^{\prime}\right)\chi_{\kappa,m}\left(\vartheta^{\prime},\varphi^{\prime}\right)Y_{l,0}^{*}\left(\vartheta^{\prime},\varphi^{\prime}\right)d\Omega^{\prime}=2\left|\kappa\right|\sqrt{4\pi}\delta_{l0}$
(59)

The electron density also simplifies to:

$\rho\left(\mathbf{r}^{\prime},\mathbf{r}\right)=\sum_{\alpha,\kappa,m}n\left(E_{\alpha,\kappa}\right)\left(\begin{array}[]{c}\psi_{\alpha,\kappa,L}\left(r^{\prime}\right)\chi_{\kappa,m}\left(\vartheta^{\prime},\varphi^{\prime}\right)\\\
\psi_{\alpha,\kappa,S}\left(r^{\prime}\right)\chi_{-\kappa,m}\left(\vartheta^{\prime},\varphi^{\prime}\right)\end{array}\right)\left(\begin{array}[]{cc}\psi_{\alpha,\kappa,L}^{*}\left(r\right)\chi_{\kappa,m}^{+}\left(\vartheta,\varphi\right)&\psi_{\alpha,\kappa,S}^{*}\left(r\right)\chi_{-\kappa,m}^{+}\left(\vartheta,\varphi\right)\end{array}\right)$
(60)

where

$n\left(E_{\alpha,\kappa}\right)=\frac{1}{e^{\frac{E_{\alpha,\kappa}-\mu}{K_{\mathrm{B}}T}}+1}$
(61)

is the Fermi-Dirac distribution function for the grand canonical ensemble
relative to a chemical potential $\mu$, $T$ is the temperature, and
$K_{\mathrm{B}}$ the Boltzmann constant.
The Hartree potential (see equation 25) reads:

$\displaystyle V_{H}\left(\mathbf{r}\right)=\int
d\mathbf{r}^{\prime}\mathrm{Tr}\left[\rho_{LL}\left(\mathbf{r}^{\prime},\mathbf{r}^{\prime}\right)+\rho_{SS}\left(\mathbf{r}^{\prime},\mathbf{r}^{\prime}\right)\right]g_{V}\left(\mathbf{r},\mathbf{r}^{\prime}\right)=$
$\displaystyle\sum_{\alpha,\kappa,m}n\left(E_{\alpha,\kappa}\right)\int
r^{\prime
2}dr^{\prime}\left|\psi_{\alpha,\kappa,L}\left(r^{\prime}\right)\right|^{2}\sum_{l=0}^{\infty}\frac{4\pi}{2l+1}g_{V,l}\left(r_{<},r_{>}\right)\times$
$\displaystyle\int\chi_{\kappa,m}^{+}\left(\vartheta^{\prime},\varphi^{\prime}\right)\chi_{\kappa,m}\left(\vartheta^{\prime},\varphi^{\prime}\right)Y_{l,0}^{*}\left(\vartheta^{\prime},\varphi^{\prime}\right)d\Omega^{\prime}Y_{l0}\left(\vartheta,\varphi\right)+$
$\displaystyle\sum_{\alpha,\kappa,m}n\left(E_{\alpha,\kappa}\right)\int
r^{\prime
2}dr^{\prime}\left|\psi_{\alpha,\kappa,S}\left(r^{\prime}\right)\right|^{2}\sum_{l=0}^{\infty}\frac{4\pi}{2l+1}g_{V,l}\left(r_{<},r_{>}\right)\times$
$\displaystyle\int\chi_{-\kappa,m}^{+}\left(\vartheta^{\prime},\varphi^{\prime}\right)\chi_{-\kappa,m}\left(\vartheta^{\prime},\varphi^{\prime}\right)Y_{l,0}^{*}\left(\vartheta^{\prime},\varphi^{\prime}\right)d\Omega^{\prime}Y_{l0}\left(\vartheta,\varphi\right)$
$\displaystyle=\int r^{\prime
2}dr^{\prime}g_{V,l}\left(r_{<},r_{>}\right)\sum_{\alpha,\kappa}n\left(E_{\alpha,\kappa}\right)\left[\left|\psi_{\alpha,\kappa,L}\left(r^{\prime}\right)\right|^{2}+\left|\psi_{\alpha,\kappa,S}\left(r^{\prime}\right)\right|^{2}\right]$
(62)

while $V_{HS}\left(\mathbf{r}\right)$ (see equation 30) turns out

$\displaystyle V_{HS}\left(\mathbf{r}\right)$ $\displaystyle=$
$\displaystyle\int
d\mathbf{r}^{\prime}\mathrm{Tr}\left[\rho_{LL}\left(\mathbf{r}^{\prime},\mathbf{r}^{\prime}\right)-\rho_{SS}\left(\mathbf{r}^{\prime},\mathbf{r}^{\prime}\right)\right]g_{S}\left(\mathbf{r},\mathbf{r}^{\prime}\right)$
(63) $\displaystyle=$ $\displaystyle\int r^{\prime
2}dr^{\prime}g_{S,l}\left(r_{<},r_{>}\right)\sum_{\alpha,\kappa}n\left(E_{\alpha,\kappa}\right)\left[\left|\psi_{\alpha,\kappa,L}\left(r^{\prime}\right)\right|^{2}-\left|\psi_{\alpha,\kappa,S}\left(r^{\prime}\right)\right|^{2}\right]$

Moreover, the Fock term in the rotational invariant case becomes (see equation
42):

$V_{F}\left(\mathbf{r}^{\prime},\mathbf{r}\right)=\left(\begin{array}[]{cc}\rho_{LL}\left(\mathbf{r}^{\prime},\mathbf{r}\right)&\rho_{LS}\left(\mathbf{r}^{\prime},\mathbf{r}\right)\\\
\rho_{SL}\left(\mathbf{r}^{\prime},\mathbf{r}\right)&\rho_{SS}\left(\mathbf{r}^{\prime},\mathbf{r}\right)\end{array}\right)g_{V}\left(\mathbf{r},\mathbf{r}^{\prime}\right)$
(64)

In Equation (64) the diagonal term can be written explicitly:

$\displaystyle\int
d\mathbf{r}\rho_{LL}\left(\mathbf{r}^{\prime},\mathbf{r}\right)g_{V}\left(\mathbf{r},\mathbf{r}^{\prime}\right)\psi\left(r\right)\chi_{\kappa,m}\left(\vartheta,\varphi\right)=$
$\displaystyle\sum_{\alpha,\kappa^{\prime},m^{\prime}}n\left(E_{\alpha,\kappa^{\prime}}\right)\psi_{\alpha,\kappa^{\prime},L}\left(r^{\prime}\right)\chi_{\kappa^{\prime},m^{\prime}}\left(\vartheta^{\prime},\varphi^{\prime}\right)\int
r^{2}drd\Omega\psi_{\alpha,\kappa^{\prime},L}^{*}\left(r\right)\psi(r)\chi_{\kappa^{\prime},m^{\prime}}^{+}\left(\vartheta,\varphi\right)\chi_{\kappa,m}\left(\vartheta,\varphi\right)\times$
$\displaystyle\sum_{l=0}^{\infty}\frac{4\pi}{2l+1}f_{V,S,l}\left(r_{<}\right)h_{V,S,l}\left(r_{>}\right)Y_{l,m-m^{\prime}}^{*}\left(\vartheta,\varphi\right)Y_{l,m-m^{\prime}}\left(\vartheta^{\prime},\varphi^{\prime}\right)$
(65)

By defining

$\displaystyle\int\chi_{\kappa_{1},m_{1}}^{+}(\Omega)\chi_{\kappa_{2},m_{2}}(\Omega)Y_{lm}^{*}(\Omega)d\Omega=\int\chi_{\kappa_{2},m_{2}}^{+}(\Omega)\chi_{\kappa_{1},m_{1}}(\Omega)Y_{lm}(\Omega)d\Omega=\mathcal{F}_{m_{1},m_{2},m}^{\kappa_{1},\kappa_{2},l}=$
$\displaystyle\sqrt{\frac{(2l+1)(2j_{1}+1)(2j_{2}+1)}{4\pi}}\delta_{m,m_{2}-m_{1}}(-1)^{m_{2}+\frac{1}{2}}\left(\begin{array}[]{ccc}j_{2}&l&j_{1}\\\
-m_{2}&m&m_{1}\end{array}\right)\left(\begin{array}[]{ccc}j_{2}&l&j_{1}\\\
\frac{1}{2}&0&-\frac{1}{2}\end{array}\right)$ (70)

one can project equation (2.2) on $\chi_{\kappa^{\prime},m^{\prime}}$,
obtaining:

$\displaystyle\int
d\Omega^{\prime}\chi_{\kappa^{\prime},m^{\prime}}^{+}\left(\vartheta^{\prime},\varphi^{\prime}\right)\int
d\mathbf{r}\rho_{LL}\left(\mathbf{r}^{\prime},\mathbf{r}\right)g_{V}\left(\mathbf{r},\mathbf{r}^{\prime}\right)\psi\left(r\right)\chi_{\kappa,m}\left(\vartheta,\varphi\right)=$
$\displaystyle\sum_{\alpha,\kappa^{\prime\prime}}n\left(E_{\alpha,\kappa^{\prime\prime}}\right)\psi_{\alpha,\kappa^{\prime\prime},L}\left(r^{\prime}\right)\int
r^{2}dr\psi_{\alpha,\kappa^{\prime\prime},L}^{*}\left(r\right)\psi(r)\sum_{l=0}^{\infty}\frac{4\pi}{2l+1}f_{V,S,l}\left(r_{<}\right)h_{V,S,l}\left(r_{>}\right)\sum_{m^{\prime\prime}}^{\infty}\mathcal{F}_{m^{\prime\prime},m,m-m^{\prime\prime}}^{\kappa^{\prime\prime},\kappa,l}\mathcal{F}_{m^{\prime\prime},m^{\prime},m-m^{\prime\prime}}^{\kappa^{\prime\prime},\kappa^{\prime},l}$
$\displaystyle=\sum_{\alpha,\kappa^{\prime\prime}}n\left(E_{\alpha,\kappa^{\prime\prime}}\right)\psi_{\alpha,\kappa^{\prime\prime},L}\left(r^{\prime}\right)\int
r^{2}dr\psi_{\alpha,\kappa^{\prime\prime},L}^{*}\left(r\right)\psi(r)\sum_{l=0}^{\infty}f_{V,S,l}\left(r_{<}\right)h_{V,S,l}\left(r_{>}\right)\frac{\left|\kappa^{\prime\prime}\right|}{2l+1}\delta_{\kappa^{\prime},\kappa}$
(71)

Similar expression can be obtained for $\rho_{SS}$. The off-diagonal term in
equation (64) can be explicitly written:

$\displaystyle\int
d\mathbf{r}\rho_{LS}\left(\mathbf{r}^{\prime},\mathbf{r}\right)g_{V}\left(\mathbf{r},\mathbf{r}^{\prime}\right)\psi\left(r\right)\chi_{\kappa,m}\left(\vartheta,\varphi\right)=$
$\displaystyle\sum_{\alpha,\kappa^{\prime},m^{\prime}}n\left(E_{\alpha,\kappa^{\prime}}\right)\psi_{\alpha,\kappa^{\prime},L}\left(r^{\prime}\right)\chi_{\kappa^{\prime},m^{\prime}}\left(\vartheta^{\prime},\varphi^{\prime}\right)\int
r^{2}drd\Omega\psi_{\alpha,-\kappa^{\prime},S}^{*}\left(r\right)\psi(r)\chi_{-\kappa^{\prime},m^{\prime}}^{+}\left(\vartheta,\varphi\right)\chi_{\kappa,m}\left(\vartheta,\varphi\right)$
$\displaystyle\times\sum_{l=0}^{\infty}\frac{4\pi}{2l+1}f_{V,S,l}\left(r_{<}\right)h_{V,S,l}\left(r_{>}\right)Y_{l,m-m^{\prime}}^{*}\left(\vartheta,\varphi\right)Y_{l,m-m^{\prime}}\left(\vartheta^{\prime},\varphi^{\prime}\right)$
(72)

which can be projected on $\chi_{\kappa^{\prime},m^{\prime}}$, obtaining:

$\displaystyle\int
d\Omega^{\prime}\chi_{\kappa^{\prime},m^{\prime}}^{+}\left(\vartheta^{\prime},\varphi^{\prime}\right)\int
d\mathbf{r}\rho_{LS}\left(\mathbf{r}^{\prime},\mathbf{r}\right)g_{V}\left(\mathbf{r},\mathbf{r}^{\prime}\right)\psi\left(r\right)\chi_{\kappa,m}\left(\vartheta,\varphi\right)=$
$\displaystyle\sum_{\alpha,\kappa^{\prime\prime}}n\left(E_{\alpha,\kappa^{\prime\prime}}\right)\psi_{\alpha,\kappa^{\prime\prime},L}\left(r^{\prime}\right)\int
r^{2}dr\psi_{\alpha,-\kappa^{\prime\prime},S}^{*}\left(r\right)\psi(r)\sum_{l=0}^{\infty}\frac{4\pi}{2l+1}f_{V,S,l}\left(r_{<}\right)h_{V,S,l}\left(r_{>}\right)\sum_{m^{\prime\prime}}^{\infty}\mathcal{F}_{m^{\prime\prime},m,m-m^{\prime\prime}}^{-\kappa^{\prime\prime},\kappa,l}\mathcal{F}_{m^{\prime\prime},m^{\prime},m-m^{\prime\prime}}^{\kappa^{\prime\prime},\kappa^{\prime},l}$
$\displaystyle=\sum_{\alpha,\kappa^{\prime\prime}}n\left(E_{\alpha,\kappa^{\prime\prime}}\right)\psi_{\alpha,\kappa^{\prime\prime},L}\left(r^{\prime}\right)\int
r^{2}dr\psi_{\alpha,-\kappa^{\prime\prime},S}^{*}\left(r\right)\psi(r)\sum_{l=0}^{\infty}f_{V,S,l}\left(r_{<}\right)h_{V,S,l}\left(r_{>}\right)\frac{\left|\kappa^{\prime\prime}\right|}{2l+1}\delta_{\kappa^{\prime},-\kappa}$
(73)

Similar expression can be derived for $\rho_{SL}$. Finally, these analytical
formulae can be used to assess the Hartree and Fock potential energy terms
(see equations 2.2,64).

####  2.2.1 Coulomb and Yukawa potentials

An important rotational–invariant scattering potential for applications in
atomic and nuclear physics is the Yukawa potential [36]:

$V(r)=\frac{V_{0}e^{\mu r}}{\mu r}$ (74)

where the potential strength $V_{0}$ is independent of $r$, and $1/\mu$
represents the effective range of the interaction. Owing to the fact that $V$
goes to zero rapidly for $r>>1/\mu$, the Yukawa potential represents the
epitome of screened Coulomb interaction in both nuclei and atoms. Moreover, it
reduces to the bare Coulomb potential for $\mu\rightarrow 0$, provided that
the ratio $V_{0}/\mu=QQ^{\prime}e^{2}$ is constant ($Q,Q^{\prime}$ are the
nuclear charges of the interacting particles, e.g. the projectile and the
target). For electrons, the central point is the expansion of the bare Coulomb
operator $\frac{1}{\left|\mathbf{r}-\mathbf{r}^{\prime}\right|}$. It can be
performed by the following Laplace expansion in spherical polar coordinates
${\displaystyle(r,\theta,\varphi)}$:

$\frac{1}{\left|\mathbf{r}-\mathbf{r}^{\prime}\right|}=\sum_{l=0}^{\infty}\frac{4\pi}{2l+1}\frac{r_{<}^{l}}{r_{>}^{l+1}}\sum_{m=-l}^{l}Y_{lm}^{*}\left(\vartheta,\varphi\right)Y_{lm}\left(\vartheta^{\prime},\varphi^{\prime}\right)$
(75)

A similar expression can be obtained for the Yukawa potential (74):

$\frac{\mathrm{e}^{-\alpha\left|\mathbf{r}-\mathbf{r}^{\prime}\right|}}{\left|\mathbf{r}-\mathbf{r}^{\prime}\right|}=\sum_{l=0}^{\infty}\frac{4\pi}{2l+1}f_{\alpha,l}\left(r_{<}\right)h_{\alpha,l}\left(r_{>}\right)\sum_{m=-l}^{l}Y_{lm}^{*}\left(\vartheta,\varphi\right)Y_{lm}\left(\vartheta^{\prime},\varphi^{\prime}\right)$
(76)

where

$f_{\alpha,l}\left(r\right)=\frac{\left(2l+1\right)!!}{\alpha^{2l}}\left(\frac{1}{r}\frac{d}{dr}\right)^{l}\frac{\sinh\left(\alpha
r\right)}{\alpha r}$ (77)

and

$h_{\alpha,l}\left(r\right)=\frac{1}{\left(2l-1\right)!!}\left(\frac{1}{r}\frac{d}{dr}\right)^{l}\frac{\mathrm{e}^{-\alpha
r}}{r}$ (78)

are spherical Bessel functions with the following limiting properties:

$\lim_{\alpha\rightarrow
0}f_{\alpha,l}\left(r\right)=r^{l},\qquad\lim_{\alpha\rightarrow
0}h_{\alpha,l}\left(r\right)=\frac{1}{r^{l+1}}$ (79)

Furthermore, $f_{\alpha,l}\left(r\right)$ e $h_{\alpha,l}\left(r\right)$
satisfy the following recursive equations:

$\displaystyle
f_{\alpha,l-1}\left(r\right)-\frac{\alpha^{2}}{\left(2l+1\right)\left(2l+3\right)}f_{\alpha,l+1}\left(r\right)$
$\displaystyle=\frac{1}{r}f_{\alpha,l}\left(r\right)$ (80)
$\displaystyle-\frac{\alpha^{2}}{\left(2l-1\right)\left(2l+1\right)}h_{\alpha,l-1}\left(r\right)+h_{\alpha,l+1}\left(r\right)$
$\displaystyle=\frac{1}{r}h_{\alpha,l}\left(r\right)$ (81)

To solve the Dirac equation for a Yukawa potential, one needs to calculate
integrals such as:

$\int\limits_{0}^{\infty}dr\int\limits_{0}^{\infty}dr^{\prime}f_{\alpha}\left(r_{<}\right)h_{\alpha}\left(r_{>}\right)\xi_{1}\left(r\right)\xi_{2}\left(r^{\prime}\right)$
(82)

where $\xi(r)$ is the product of two wavefunctions. These integrals can be
decomposed as follows:

$\displaystyle I$
$\displaystyle=\int\limits_{0}^{\infty}dr\int\limits_{0}^{r}dr^{\prime}f_{\alpha}\left(r^{\prime}\right)h_{\alpha}\left(r\right)\xi_{1}\left(r\right)\xi_{2}\left(r^{\prime}\right)+\int\limits_{0}^{\infty}dr\int\limits_{r}^{\infty}dr^{\prime}f_{\alpha}\left(r\right)h_{\alpha}\left(r^{\prime}\right)\xi_{1}\left(r\right)\xi_{2}\left(r^{\prime}\right)$
(83)
$\displaystyle=\int\limits_{0}^{\infty}dr\left[\mathcal{F}_{f_{\alpha}\xi_{2}}\left(r\right)-\mathcal{F}_{f_{\alpha}\xi_{2}}\left(0\right)\right]h_{\alpha}\left(r\right)\xi_{1}\left(r\right)-\int\limits_{0}^{\infty}dr\left[\mathcal{F}_{h_{\alpha}\xi_{2}}\left(r\right)-\mathcal{F}_{h_{\alpha}\xi_{2}}\left(\infty\right)\right]h_{\alpha}\left(r\right)\xi_{1}\left(r\right)$

where $\mathcal{F}_{\xi}\left(r\right)$ is the antiderivative of $\xi(r)$. In
the next section we will show a numerical method to compute this integral in a
radial grid.

###  2.3 Computational methods for spherically symmetric systems

The general form of the Dirac equation ((2.1)) for a spherical symmetric
problem, such as particles interacting via mean-field or Yukawa potentials,
can be written as follows:

$\frac{d\psi}{dr}=F_{E}(\psi,r)$ (84)

where $\psi$ is a two-component vector. $F_{E}$, being a linear function of
both $\psi$ and of the energy $E$ (see Equation (53)), can be written as
$F_{E}(\psi,r)=A(r)\psi$ where $A(r)$ is a $2\times 2$ linear operator. The
difficulty to find the numerical solution of Equation (84) for many-fermion
systems is due to the fact that $F_{E}$ contains a non-local term, related to
the Fock exchange (see equation 64). To circumvent this issue, the non-local
Fock term is replaced by an expression derived using the local density
approximation (LDA) for the electron gas [37, 38]

$V_{ex}=\frac{3}{2}\alpha_{X}\Big{[}\frac{3}{\pi}\rho(r)\Big{]}^{1/3}$ (85)

where $\rho(r)$ is the electron density and $\alpha_{X}=3/2$.

According to (58) the solutions of the Dirac equation (53) in a spherically
symmetric potential can be factorised into the product of monodimensional
radial and angular functions [39]:

$\psi(\mathbf{r})=\psi_{\kappa,\mu}(\mathbf{r})=\left(\begin{array}[]{ll}\frac{u_{\kappa}(r)}{r}\
\chi_{\kappa,\mu}(\Omega)\\\ i\frac{v_{\kappa}(r)}{r}\
\chi_{-\kappa,\mu}(\Omega)\end{array}\right)$ (86)

where $\chi_{\kappa,\mu}$ are the tensor product of orbital and spin spherical
harmonics characterised by the Dirac quantum numbers $\kappa,\mu$,
respectively.

By using the factorization (86) and considering for the sake of simplicity a
system of many-fermions interacting via a mean-field potential, the Dirac
equation (53) reads:

$\frac{\partial}{\partial r}\left(\begin{array}[]{c}u_{\kappa}(r)\\\
v_{\kappa}(r)\end{array}\right)=\left(\begin{array}[]{cc}-\frac{\kappa}{r}&\frac{E}{c}-\frac{V(r)}{c}+c\\\
\frac{V(r)}{c}+c-\frac{E}{c}&\frac{\kappa}{r}\end{array}\right)\left(\begin{array}[]{c}u_{\kappa}(r)\\\
v_{\kappa}(r)\end{array}\right)$ (87)

where the total potential is given by ($Z$ is the atomic number):

$V(r)=-\frac{Z}{r}+\int\frac{\rho(r^{\prime})}{r_{>}}d^{3}r^{\prime}-V_{ex}(r)$
(88)

We notice that the Dirac equation (87) is valid for any central symmetric
potential. For example, by using the Wood-Saxon potential [40], which
typically models the nucleon-nucleon interaction, the calculation of the
nuclear wavefunctions can in principle be carried out [8].
Equations (87) with the local interaction (88) are solved self-consistently.
We stress that, as for the non-relativistic HF, the eigenvalues of the radial
equation (87) represent the electron binding energies.
To integrate numerically the equation (87), we introduce now a modified Runge-
Kutta method, which is known as the method of collocations. This method
proceeds by dividing the interval $[0,R]$ over which the solution is sought
for ($R$ can even tend to $\infty$) in a number of grid points
$[0=r_{0},\cdots r_{j},\cdots r_{N}=R]$. The latter can be determined as the
roots of the $n-$th Legendre polynomial. Each sub-interval $[r_{j-1},r_{j}]$
is further discretised in a number of $s_{j}$ points $\\{x_{j,l}\\}$, where
the derivative of the wavefunction appearing in the general Dirac equation
(84) can be approximated by

$\frac{d\psi}{dr}\simeq\sum_{l=1}^{s_{j}}q_{j,l}\phi_{j,l}^{\prime}(r)$ (89)

and $\psi$ is obtained by integrating the equation (89) as follows:

$\displaystyle\psi(r)$ $\displaystyle=$
$\displaystyle\psi(r_{j-1})+\int\limits_{r_{j-1}}^{r}\sum_{l=1}^{s_{j}}q_{j,l}\phi_{j,l}^{\prime}(x)dx$
(90) $\displaystyle=$
$\displaystyle\sum_{i=1}^{j-1}\sum_{l=1}^{s_{i}}q_{i,l}\int\limits_{r_{i-1}}^{r_{i}}\phi_{i,l}^{\prime}(x)dx+\sum_{l=1}^{s_{j}}q_{j,l}\int\limits_{r_{j-1}}^{r}\phi_{j,l}^{\prime}(x)dx$

imposing $\phi_{j,l}^{\prime}(x_{j,m})=\delta_{lm}$. By defining:

$\phi_{j,l}(r)=\int\limits_{r_{j-1}}^{r}\phi_{j,l}^{\prime}(x)dx,\qquad
a_{j,lm}=\int\limits_{r_{j-1}}^{x_{j,l}}\phi_{j,m}^{\prime}(x)dx,\qquad
b_{j,l}=\int\limits_{r_{j-1}}^{r_{j}}\phi_{j,l}^{\prime}(x)dx$ (91)

one obtains:

$\frac{d\psi}{dr}(x_{j,l})=q_{j,l}=F_{E}\left(\sum_{i=1}^{j-1}\sum_{l=1}^{s_{i}}b_{i,m}q_{i,m}+\sum_{l=1}^{s_{j}}a_{j,lm}q_{j,m},x_{j,l}\right)$
(92)

The Dirac equation (87) can now be recast as an eigenvalue problem

$c\gamma\left(\begin{array}[]{c}q_{1}\\\ q_{2}\\\ q_{3}\\\ \cdots\\\
q_{N}\end{array}\right)=\left(E-\beta
mc^{2}-V\right)\left(\begin{array}[]{ccccc}a_{1}&0&0&\cdots&0\\\
b_{1}&a_{2}&0&\cdots&0\\\ b_{1}&b_{2}&a_{3}&\cdots&0\\\
\cdots&\cdots&\cdots&\cdots&\cdots\\\
b_{1}&b_{2}&b_{3}&\cdots&a_{N}\end{array}\right)\left(\begin{array}[]{c}q_{1}\\\
q_{2}\\\ q_{3}\\\ \cdots\\\ q_{N}\end{array}\right)$ (93)

where $\gamma$ and $\beta$ are $2\times 2$ matrices. In compact form, the
equation (87) can be written:

$\left[c\gamma+\left(V+\beta mc^{2}\right)A-EA\right]q=0$ (94)

where

$A=\left(\begin{array}[]{ccccc}a_{1}&0&0&\cdots&0\\\ b_{1}&a_{2}&0&\cdots&0\\\
b_{1}&b_{2}&a_{3}&\cdots&0\\\ \cdots&\cdots&\cdots&\cdots&\cdots\\\
b_{1}&b_{2}&b_{3}&\cdots&a_{N}\end{array}\right)$ (95)

Finally, by transforming $A$ into a lower bidiagonal matrix, the equation (94)
reads:

$\left[c\gamma Q+\left(V+\beta mc^{2}\right)AQ-EAQ\right]q^{\prime}=0$ (96)

where $q^{\prime}=Q^{-1}q$, $\psi=Aq=AQq^{\prime}$, and

$Q=\left(\begin{array}[]{rrcrrr}1&0&\cdots&0&0&0\\\
-b_{2}^{-1}b_{1}&1&\cdots&0&0&0\\\ 0&-b_{3}^{-1}b_{2}&\cdots&1&0&0\\\
\cdots&\cdots&\cdots&\cdots&\cdots&\cdots\\\
0&0&\cdots&-b_{N-1}^{-1}b_{N-2}&1&0\\\ 0&0&\cdots&0&0&1\end{array}\right)$
(97) $AQ=\left(\begin{array}[]{rrcrrr}a_{1}&0&0&\cdots&0&0\\\
b_{1}-a_{2}b_{2}^{-1}b_{1}&a_{2}&0&\cdots&0&0\\\
0&b_{2}-a_{3}b_{3}^{-1}b_{2}&a_{3}&\cdots&0&0\\\
\cdots&\cdots&\cdots&\cdots&\cdots&\cdots\\\
0&0&\cdots&b_{N-2}-a_{N-1}b_{N-1}^{-1}b_{N-2}&a_{N-1}&0\\\
0&0&\cdots&0&b_{N-1}&a_{N}\end{array}\right)$ (98)

The tridiagonal form (96) is computationally efficient, particularly when we
are dealing with quasi-local potentials, as in that case the matrices $Q$,
$AQ$, $V$, and $VAQ$ are almost diagonal. The eigenvalues $E$ can be indeed
reckoned via the $LDU$ decomposition. The spatial extension of the region
$[0,R_{max}]$ where the differential equation (84) is solved depends weather
we are dealing with bound or continuum state problems. Bound state
wavefunctions typically die off after a few atomic radii ($R_{max}=25$ a.u.),
while continuum states extend to hundreds atomic units ($R_{max}\geq 100$
a.u.). In general, if one is to solve a continuum scattering state problem,
e.g. the simulation of electron emission from heavy atoms, the continuum
wavefunction asymptotic behaviour ($R_{max}\geq r\rightarrow\infty$) must be
matched via Coulomb functions [39], which represent the solution for a charged
particle in an ionic field:

$u\to\frac{\sqrt{\epsilon+c^{2}}}{\sqrt{\pi
p}}\mathrm{cos}(pr-\eta\ln(pr)+\delta)\\\
v\to-\frac{\sqrt{\epsilon-c^{2}}}{\sqrt{\pi
p}}\mathrm{sin}(pr-\eta\ln(pr)+\delta)$ (99)

where $\eta=Zm/p$ ($Z$=atomic number, $m$=electron mass=1), and
$\epsilon=\sqrt{p^{2}c^{2}+c^{4}}$ is the free electron kinetic energy. This
leads to the normalization condition of the continuum wavefunctions,
$\frac{\pi
p}{\sqrt{\epsilon+c^{2}}}\Big{(}u^{2}+\sqrt{\frac{\epsilon+c^{2}}{\epsilon-c^{2}}}v^{2}\Big{)}=1$.
To avoid a critical dependence of the wavefunction normalization on small
changes of the box size due to the $u$ and $v$ oscillating behaviour at the
borders of the integration region, we average over an entire oscillation
period of the wavefunction.
We point out that this numerical technique, based on the assessment of the
wavefunctions on a radial grid, scales almost linearly with the number of
points in the mesh. This approach is used below to calculate the electronic
structure of the Au+ ion.

##  3 Methods for polycentric, non-spherically symmetric, non-periodic systems

A second approach, which can be easily extended to polycentric molecular
systems, relies on the projected potential method [3, 4]. Within this
framework, the non-kinetic terms of equation (53) are projected into a finite
functional space $\cal{G}$ spanned by a number of $L^{2}$-functions
($\lambda$), which define a projector

$\pi=\sum_{\lambda}|\lambda\rangle\langle\lambda|$ (100)

where $g_{j}$ are $L^{2}$ basis functions. Inside that subspace the equation
of motion can be written:

$\left(H_{0}-E+\pi V\pi\right)\psi=0$ (101)

or, equivalently,

$\pi\psi=\pi\frac{1}{E-H_{0}}\pi V\pi\psi$ (102)

where

$H_{0}=\left(\begin{array}[]{cc}mc^{2}&-c\boldsymbol{\sigma}\cdot
i\boldsymbol{\nabla}\\\ -c\boldsymbol{\sigma}\cdot
i\boldsymbol{\nabla}&-mc^{2}\end{array}\right)$ (103)

and

$V=\pi\left(\begin{array}[]{cc}W_{V}+W_{S}+\mathbf{A}_{P}\cdot\boldsymbol{\sigma}&-\boldsymbol{\sigma}\cdot\mathbf{A}+W_{PS}\\\
-\boldsymbol{\sigma}\cdot\mathbf{A}+W_{PS}&W_{V}+\mathbf{A}_{P}\cdot\boldsymbol{\sigma}-W_{S}\end{array}\right)\pi$
(104)

The kinetic term (103) is not projected into the Hilbert subspace spanned by
the projector in order to recover the continuum for dealing with scattering
states. At variance, the potential terms in equation (104) are projected into
the functional space using the projector operator (100). For many-electron
systems the total Coulomb potential (88) upon projection reads:

$V_{\alpha}=\sum_{\lambda\mu\nu\tau}|\lambda\rangle
S^{-1}_{\lambda\mu}\langle\mu|V|\nu\rangle S^{-1}_{\nu\tau}\langle\tau|;\quad
S_{\lambda\mu}=\langle\lambda|\mu\rangle$ (105)

The projected potential (105) is well-suited for bound-state calculations,
where wavefunctions are appreciably confined into a finite region. To model
scattering states we require that the projected potential $V_{\alpha}$
correctly reproduces the effect of the long range part of the true potential
at least in a part of the asymptotic region, where $V_{\alpha}(r)\simeq V(r)$
(for $r$ large). This fact guarantees that the scattering wavefunction far
from the scattering center has the correct form [3, 4].
The solution of equation (101), or equivalently of equation (102), of course
delivers the eigenvalues of the projected Hamiltonian $H_{0}+\pi V\pi$ rather
than those of the complete Hamiltonian $H_{0}+V$. However, we notice that the
eigenvalues of the projected and complete Hamiltonian coincide provided that
the vectors $\psi$ and $V\psi$ belong to the projected functional subspace.
Equation (102) can then be solved if the matrix elements

$\pi\frac{1}{E^{2}-m^{2}c^{4}-c^{2}p^{2}}\pi,\quad\pi\frac{1}{E^{2}-m^{2}c^{4}-c^{2}p^{2}}\vec{p}\pi$
(106)

of the projected Dirac Hamiltonian are reckoned within the chosen functional
space. To find the elements of these Green’s operator we refer the interested
reader to Refs. [3, 4].

###  3.1 Gaussian functions

The functions that we use to project the equation (101) and to expand the
wavefunctions are HGBF. Gaussian functions (GF) are particularly suited for
representing both the bound and continuum orbitals [3] in multicentric
molecular systems, as analytic expressions of monoelectronic and bielectronic
integrals, necessary for the calculations of the spectral quantities, can be
explicitly derived. These integrals are the average values of different
operators evaluated using a mixed basis set of HGBF of every order and
centrature. The structure of the HGBF is the following:

$g(\mathbf{r})=g(u,v,w;a,\mathbf{R};\mathbf{r})=N{\frac{\partial^{u+v+w}}{\partial
X^{u}\partial Y^{v}\partial
Z^{w}}}(\frac{2\alpha}{\pi})^{3/4}exp[-\alpha(\mathbf{r}-\mathbf{R})^{2}]$
(107)

where $\mathbf{R}\equiv(X,Y,Z)$ gives the position where $g$ is centered,
$\alpha$ is a coefficient determining the GF width, the order of derivation
($u,v,w$) determines the symmetry type, and $N$ is a normalization factor

$N=[\alpha^{l}(2u-1)!!(2v-1)!!(2w-1)!!]^{-1/2}{,}\qquad l=u+v+w$

GFs build up a $L^{2}$ non-orthogonal basis set characterized by the
remarkable property that the product of GFs (represented by coefficients
$\alpha_{1},\alpha_{2}$ and centers $R_{1},R_{2}$, respectively) is a GF
itself:

$g\left(\alpha_{1},\mathbf{R}_{1};\mathbf{r}\right)g\left(\alpha_{2},\mathbf{R}_{2};\mathbf{r}\right)=g\left(\alpha_{1}+\alpha_{2},\frac{\alpha_{1}\mathbf{R}_{1}+\alpha_{2}\mathbf{R}_{2}}{\alpha_{1}+\alpha_{2}};\mathbf{r}\right)e^{-\frac{\alpha_{1}\alpha_{2}}{\alpha_{1}+\alpha_{2}}\left(\mathbf{R}_{1}-\mathbf{R}_{2}\right)}$
(108)

with normalization condition given by

$\left\langle
g_{1}|g_{2}\right\rangle=\left(\frac{2\alpha_{1}}{\pi}\right)^{\frac{3}{4}}\left(\frac{2\alpha_{2}}{\pi}\right)^{\frac{3}{4}}\left(\frac{\pi}{\alpha_{1}+\alpha_{2}}\right)^{\frac{3}{2}}e^{-\frac{\alpha_{1}\alpha_{2}}{\alpha_{1}+\alpha_{2}}\left(\mathbf{R}_{1}-\mathbf{R}_{2}\right)}=\left(\frac{2\sqrt{\alpha_{1}\alpha_{2}}}{\alpha_{1}+\alpha_{2}}\right)^{\frac{3}{2}}e^{-\frac{\alpha_{1}\alpha_{2}}{\alpha_{1}+\alpha_{2}}\left(\mathbf{R}_{1}-\mathbf{R}_{2}\right)}$

###  3.2 Mono-electronic integrals

GFs allow the analytical calculation of the projected Coulomb potential
(Hartree term):

$\int
d\mathbf{r}\frac{1}{r}g\left(\alpha,\mathbf{R};\mathbf{r}\right)=\left(\frac{\pi}{\alpha}\right)^{\frac{3}{2}}\frac{\mathrm{erf}\left(\sqrt{\alpha}R\right)}{R}$
(109)

This task can be achieved by exploting the fact that the Yukawa potential
reduces to the electron-electron Coulomb potential in the short wave-length
limit ($\mu\rightarrow 0$, see equation 74). Mono-electronic integrals for the
Yukawa potential can be written as follows:

$\displaystyle\small\left\langle
m\frac{e^{-\frac{r}{\xi}}}{r}\right\rangle=\int
d\mathbf{r}\frac{e^{-\frac{r}{\xi}}}{r}g\left(\alpha,\mathbf{R};\mathbf{r}\right)=2\pi\int\limits_{0}^{\infty}r^{2}dr\left[e^{-\alpha\left(r^{2}-2rR+R^{2}\right)}-e^{-\alpha\left(r^{2}+2rR+R^{2}\right)}\right]\frac{e^{-\frac{r}{\xi}}}{r}$
$\displaystyle=\frac{1}{2}\left(\frac{\pi}{\alpha}\right)^{\frac{3}{2}}\frac{\left[-\mathrm{erfc}\left(\frac{1}{2\sqrt{\alpha}\xi}+\sqrt{\alpha}R\right)e^{\left(\frac{1}{4\alpha\xi^{2}}+\frac{R}{\xi}\right)}+\mathrm{erfc}\left(\frac{1}{2\sqrt{\alpha}\xi}-\sqrt{\alpha}R\right)e^{\left(\frac{1}{4\alpha\xi^{2}}-\frac{R}{\xi}\right)}\right]}{R}$
$\displaystyle=\frac{1}{2}\left(\frac{\pi}{\alpha}\right)^{\frac{3}{2}}\frac{\left[-\mathrm{erfcx}\left(\frac{1}{2\sqrt{\alpha}\xi}+\sqrt{\alpha}R\right)+\mathrm{erfcx}\left(\frac{1}{2\sqrt{\alpha}\xi}-\sqrt{\alpha}R\right)\right]e^{-\alpha
R^{2}}}{R}$ (110)

where $\mathrm{erfcx}(x)=\mathrm{erfc}(x)e^{x^{2}}$ and

$\mathrm{erfc}(x)=\frac{2}{\sqrt{\pi}}\int\limits_{x}^{\infty}e^{-u^{2}}du,\qquad\mathrm{erfcx}(x)=\frac{2}{\sqrt{\pi}}\int\limits_{x}^{\infty}e^{x^{2}-u^{2}}du$

By Taylor expanding Equation (3.2) and observing that

$\displaystyle\frac{d}{dx}\mathrm{erfcx}(a+x)e^{-x^{2}}=-\frac{2}{\sqrt{\pi}}e^{-x^{2}}+2a\
\mathrm{erfcx}(a+x)e^{-x^{2}}$
$\displaystyle\frac{d}{dx}\mathrm{erfcx}(a-x)e^{-x^{2}}=\frac{2}{\sqrt{\pi}}e^{-x^{2}}-2a\
\mathrm{erfcx}(a-x)e^{-x^{2}}$ (111)

we obtain finally:

$\displaystyle\left\langle
m\frac{e^{-\frac{r}{\xi}}}{r}\right\rangle=\frac{2\pi^{\frac{3}{2}}}{\alpha}\Big{[}\frac{1}{\sqrt{\pi}}-a\
\mathrm{erfcx}(a)+\frac{1}{3}\left(-2a^{3}\mathrm{erfcx}(a)+\frac{2a^{2}-1}{\sqrt{\pi}}\right)x^{2}+$
$\displaystyle+\frac{2}{15}\left(-4a^{5}\mathrm{erfcx}(a)+\frac{4a^{4}-2a^{2}+3}{\sqrt{\pi}}\right)+\cdots\Big{]}$
(112)

where $a=\frac{1}{2\sqrt{\alpha}\xi}$ and $x=\sqrt{\alpha}R$.
Furthermore, one can also find analytically the integral of a Coulomb and
Gaussian potentials between Gaussian functions as:

$\left\langle
g_{1}|\frac{1}{\left|\vec{r}-\vec{C}\right|}|g_{2}\right\rangle=\left(\frac{2\sqrt{\alpha_{1}\alpha_{2}}}{\pi}\right)^{\frac{3}{2}}e^{-\frac{\alpha_{1}\alpha_{2}}{\alpha_{1}+\alpha_{2}}\left(\vec{R}_{1}-\vec{R}_{2}\right)^{2}}\left(\frac{\pi}{\alpha_{1}+\alpha_{2}}\right)^{\frac{3}{2}}\frac{\mathrm{erf}\left(\sqrt{\left(\alpha_{1}+\alpha_{2}\right)\left|\frac{\alpha_{1}\vec{R}_{1}+\alpha_{2}\vec{R}_{2}}{\alpha_{1}+\alpha_{2}}-\vec{C}\right|^{2}}\right)}{\left|\frac{\alpha_{1}\vec{R}_{1}+\alpha_{2}\vec{R}_{2}}{\alpha_{1}+\alpha_{2}}-\vec{C}\right|}$
(113)

and:

$\displaystyle\left\langle
g_{1}|e^{-\gamma\left(\vec{r}-\vec{C}\right)^{2}}|g_{2}\right\rangle=\left(\frac{2\sqrt{\alpha_{1}\alpha_{2}}}{\pi}\right)^{\frac{3}{2}}e^{-\frac{\alpha_{1}\alpha_{2}}{\alpha_{1}+\alpha_{2}}\left(\vec{R}_{1}-\vec{R}_{2}\right)^{2}}\int
d^{3}re^{-\left(\alpha_{1}+\alpha_{2}\right)\left(\vec{r}-\frac{\alpha_{1}\vec{R}_{1}+\alpha_{2}\vec{R}_{2}}{\alpha_{1}+\alpha_{2}}\right)^{2}}e^{-\gamma\left(\vec{r}-\vec{C}\right)^{2}}$
$\displaystyle=\left(\frac{2\sqrt{\alpha_{1}\alpha_{2}}}{\pi}\right)^{\frac{3}{2}}e^{-\frac{\alpha_{1}\alpha_{2}}{\alpha_{1}+\alpha_{2}}\left(\vec{R}_{1}-\vec{R}_{2}\right)^{2}}e^{-\frac{\left(\alpha_{1}+\alpha_{2}\right)\gamma}{\alpha_{1}+\alpha_{2}+\gamma}\left(\frac{\alpha_{1}\vec{R}_{1}+\alpha_{2}\vec{R}_{2}}{\alpha_{1}+\alpha_{2}}-\vec{C}\right)^{2}}\int
d^{3}re^{-\left(\alpha_{1}+\alpha_{2}+\gamma\right)\left(\vec{r}-\frac{\alpha_{1}\vec{R}_{1}+\alpha_{2}\vec{R}_{2}+\gamma\vec{C}}{\alpha_{1}+\alpha_{2}+\gamma}\right)^{2}}$
$\displaystyle=\left(\frac{2\sqrt{\alpha_{1}\alpha_{2}}}{\pi}\right)^{\frac{3}{2}}e^{-\frac{\alpha_{1}\alpha_{2}}{\alpha_{1}+\alpha_{2}}\left(\vec{R}_{1}-\vec{R}_{2}\right)^{2}}e^{-\frac{\left(\alpha_{1}+\alpha_{2}\right)\gamma}{\alpha_{1}+\alpha_{2}+\gamma}\left(\frac{\alpha_{1}\vec{R}_{1}+\alpha_{2}\vec{R}_{2}}{\alpha_{1}+\alpha_{2}}-\vec{C}\right)^{2}}\left(\frac{\pi}{\alpha_{1}+\alpha_{2}+\gamma}\right)^{\frac{3}{2}}$
(114)

###  3.3 Green’s function matrix elements.

The Green’s function projected on GFs in equation (102) reads:

$\left\langle g_{1}|\frac{1}{p_{0}^{2}+\nabla^{2}}|g_{2}\right\rangle$ (115)

where $p_{0}\in\mathbb{C}$. The integral in Equation (115) can be more easily
calculated in momentum space, as the kinetic operator ($\nabla^{2}$) is
diagonal:

$\displaystyle\left\langle
g_{1}|\frac{1}{p_{0}^{2}+\nabla^{2}}|g_{2}\right\rangle=$
$\displaystyle\left(\frac{2\pi}{\sqrt{\alpha_{1}\alpha_{2}}}\right)^{\frac{3}{2}}\int
d\mathbf{p}\frac{e^{-\frac{\alpha_{1}+\alpha_{2}}{4\alpha_{1}\alpha_{2}}p^{2}+i\mathbf{p}\cdot\left(\mathbf{R}_{1}-\mathbf{R}_{2}\right)}}{p_{0}^{2}-p^{2}}=\left(\frac{2\pi}{\sqrt{\alpha_{1}\alpha_{2}}}\right)^{\frac{3}{2}}4\pi\int\limits_{0}^{\infty}\frac{p^{2}dp}{p\left|\mathbf{R}_{1}-\mathbf{R}_{2}\right|}\frac{e^{-\frac{\alpha_{1}+\alpha_{2}}{4\alpha_{1}\alpha_{2}}p^{2}}\sin\left(p\left|\mathbf{R}_{1}-\mathbf{R}_{2}\right|\right)}{p_{0}^{2}-p^{2}}=$
$\displaystyle\left(\frac{2\pi}{\sqrt{\alpha_{1}\alpha_{2}}}\right)^{\frac{3}{2}}2\pi\int\limits_{-\infty}^{\infty}\frac{pdp}{\left|\mathbf{R}_{1}-\mathbf{R}_{2}\right|}\frac{e^{-\frac{\alpha_{1}+\alpha_{2}}{4\alpha_{1}\alpha_{2}}p^{2}+ip\left|\mathbf{R}_{1}-\mathbf{R}_{2}\right|}}{p_{0}^{2}-p^{2}}$
(116)

These integrals can be reduced to integrals of the $\mathrm{erf}$ function in
the complex space by considering that $w(z)=e^{-z^{2}}\mathrm{erfc}(-iz)$ and
that

$\int\limits_{-\infty}^{\infty}\frac{e^{-t^{2}}}{z-t}dt=\left\\{\begin{array}[]{@{\kern
2.5pt}lL}\hfill-i\pi w(z)&ifIm(z)$>$0\\\ \hfill-i\pi
w(z)-2e^{-z^{2}}&ifIm(z)$<$0.\end{array}\right.$ (117)

###  3.4 The two-step recurrence relation of $G_{m}(x)$.

We observe that the integrals involving GFs of higher order can be obtained
through the differentiation of $s$-type GFs, because of the Leibnitz’s rule
that allows one to take the differentiation operator out of the integral sign.
In particular, the Yukawa potential matrix element can be calculated for
Gaussians of any order ($s,p,d$ symmetry) by a recurrence relation.
To demonstrate that, we define the following quantities:

$G_{m}(x)=\left(-\frac{1}{2x}\frac{d}{dx}\right)^{m}\left(\frac{\alpha}{\pi}\right)^{\frac{3}{2}}\left\langle
m\frac{e^{-\frac{r}{\xi}}}{r}\right\rangle=\left(-\frac{1}{2x}\frac{d}{dx}\right)^{m}\frac{\sqrt{\alpha}}{2}\frac{\left[-\mathrm{erfcx}\left(a+x\right)+\mathrm{erfcx}\left(a-x\right)\right]e^{-x^{2}}}{x}$
(118)

$I_{0}(a,x)=\left[-\mathrm{erfcx}\left(a+x\right)+\mathrm{erfcx}\left(a-x\right)\right]e^{-x^{2}}$
(119)
$I_{1}(a,x)=\left[\mathrm{erfcx}\left(a+x\right)+\mathrm{erfcx}\left(a-x\right)\right]e^{-x^{2}}$
(120)

The first element of the recurrence relation reads:

$\small G_{0}(x)=\frac{\sqrt{\alpha}}{2}\frac{I_{0}(a,x)}{x}$ (121)

while the second one can be obtained by observing that
$\frac{dI_{0}}{dx}=\frac{4}{\sqrt{\pi}}e^{-x^{2}}-2aI_{1}$ and
$\frac{dI_{1}}{dx}=-2aI_{0}$:

$\displaystyle G_{1}(x)$
$\displaystyle=\left(\frac{1}{2x}\frac{d}{dx}\right)G_{0}(x)=\frac{\sqrt{\alpha}}{4x^{2}}\left\\{-\frac{4}{\sqrt{\pi}}e^{-x^{2}}+2aI_{1}(a,x)+\frac{I_{0}(a,x)}{x}\right\\}$
(122)
$\displaystyle=\frac{1}{x^{2}}\left[\frac{1}{2}G_{0}(x)-\sqrt{\frac{\alpha}{\pi}}e^{-x^{2}}+\frac{\sqrt{\alpha}}{2}aI_{1}(a,x)\right]$

Finally, the recurrence relation can be generally written as:

$G_{m+2}(x)=\frac{1}{x^{2}}\left[\left(m+\frac{3}{2}\right)G_{m+1}(x)-\sqrt{\frac{\alpha}{\pi}}e^{-x^{2}}+a^{2}G_{m}(x)\right]$
(123)

###  3.5 Bi-electronic integrals

Bi-electronic integrals, which are the matrix elements of the screened Coulomb
potential between GFs centered in different positions, can be written:

$\displaystyle\left\langle
s_{1}s_{3}|\frac{e^{-\frac{\left|\mathbf{r}-\mathbf{r}^{\prime}\right|}{\xi}}}{\left|\mathbf{r}-\mathbf{r}^{\prime}\right|}|s_{2}s_{4}\right\rangle=$
$\displaystyle\left(\frac{2\alpha_{1}}{\pi}\right)^{\frac{3}{4}}\left(\frac{2\alpha_{2}}{\pi}\right)^{\frac{3}{4}}\left(\frac{2\alpha_{3}}{\pi}\right)^{\frac{3}{4}}\left(\frac{2\alpha_{4}}{\pi}\right)^{\frac{3}{4}}\int
d\mathbf{r}\int
d\mathbf{r}^{\prime}e^{-\alpha_{1}\left(\mathbf{r}-\mathbf{R}_{1}\right)^{2}}e^{-\alpha_{2}\left(\mathbf{r}-\mathbf{R}_{2}\right)^{2}}e^{-\alpha_{3}\left(\mathbf{r}^{\prime}-\mathbf{R}_{3}\right)^{2}}e^{-\alpha_{4}\left(\mathbf{r}^{\prime}-\mathbf{R}_{4}\right)^{2}}=$
$\displaystyle\left[\frac{2\left(\alpha_{1}\alpha_{2}\alpha_{3}\alpha_{4}\right)^{\frac{1}{4}}}{\pi}\right]^{3}e^{-\frac{\alpha_{1}\alpha_{2}}{\alpha_{1}+\alpha_{2}}\left(\mathbf{R}_{1}-\mathbf{R}_{2}\right)}e^{-\frac{\alpha_{3}\alpha_{4}}{\alpha_{3}+\alpha_{4}}\left(\mathbf{R}_{3}-\mathbf{R}_{4}\right)}\int
d\mathbf{r}\int
d\mathbf{r}^{\prime}e^{-\left(\alpha_{1}+\alpha_{2}\right)\left(\mathbf{r}-\mathbf{P}_{12}\right)^{2}}e^{-\left(\alpha_{3}+\alpha_{4}\right)\left(\mathbf{r}^{\prime}-\mathbf{P}_{34}\right)^{2}}\frac{e^{-\frac{\left|\mathbf{r}-\mathbf{r}^{\prime}\right|}{\xi}}}{\left|\mathbf{r}-\mathbf{r}^{\prime}\right|}$
(124)

where
$\mathbf{P}_{12}=\frac{\alpha_{1}\mathbf{R}_{1}+\alpha_{2}\mathbf{R}_{2}}{\alpha_{1}+\alpha_{2}}$
and
$\mathbf{P}_{34}=\frac{\alpha_{3}\mathbf{R}_{3}+\alpha_{4}\mathbf{R}_{4}}{\alpha_{3}+\alpha_{4}}$.
By defining $\mathbf{s}=\frac{\mathbf{r}+\mathbf{r}^{\prime}}{2}$ and
$\mathbf{q}=\mathbf{r}-\mathbf{r}^{\prime}$ one obtains:

$\displaystyle\int d\mathbf{r}\int
d\mathbf{r}^{\prime}e^{-\left(\alpha_{1}+\alpha_{2}\right)\left(\mathbf{r}-\mathbf{P}_{12}\right)^{2}}e^{-\left(\alpha_{3}+\alpha_{4}\right)\left(\mathbf{r}^{\prime}-\mathbf{P}_{34}\right)^{2}}\frac{e^{-\frac{\left|\mathbf{r}-\mathbf{r}^{\prime}\right|}{\xi}}}{\left|\mathbf{r}-\mathbf{r}^{\prime}\right|}=\int
d\mathbf{s}\int
d\mathbf{q}e^{-\left(\alpha_{1}+\alpha_{2}\right)\left(\mathbf{s}+\frac{\mathbf{q}}{2}-\mathbf{P}_{12}\right)^{2}}e^{-\left(\alpha_{3}+\alpha_{4}\right)\left(\mathbf{s}-\frac{\mathbf{q}}{2}-\mathbf{P}_{34}\right)^{2}}\frac{e^{-\frac{q}{\xi}}}{q}=$
$\displaystyle\left(\frac{\pi}{\alpha_{1}+\alpha_{2}+\alpha_{3}+\alpha_{4}}\right)^{\frac{3}{2}}\int
d\mathbf{q}e^{-\frac{\left(\alpha_{1}+\alpha_{2}\right)\left(\alpha_{3}+\alpha_{4}\right)}{\alpha_{1}+\alpha_{2}+\alpha_{3}+\alpha_{4}}\left[\mathbf{q}-\left(\mathbf{P}_{12}-\mathbf{P}_{34}\right)\right]^{2}}\frac{e^{-\frac{q}{\xi}}}{q}$
(125)

which can be reduced to an integral of the type reported in equation (3.2).
Further details on the calculation of mono- and bi-electronic integrals using
GFs of different orders (modified Gaussian functions) can be found in [41].
We point out that the self-consistent solution of the Dirac equation using the
HGFB numerical approach scales from cubically, owing to the matrix
diagonalization, to quartically, due to the calculation of the bi-electronic
integrals, as a function of the number of Gaussians in the basis set. This
approach is used below to calculate the electronic structure of the W@Au12
nanocluster.

##  4 Radial functions and Gaussian basis set at work: application to the
electronic structure calculation of heavy-element atoms and clusters

To show the potential of our relativistic approach, we calculate first the
electronic structure of the Au+ ion, by comparing the results obtained by both
the radial and Gaussian basis sets, and second of the W@Au12 superatomic
system, where the inclusion of relativistic effects is necessary to describe
accurately the electronic structure [42].

###  4.1 Au and W atoms from relativistic HF using radial grids and Gaussian
function basis sets

The calculations were carried out for Au+, as it is a closed-shell system
(neutral gold has indeed 79 electrons), using i) a contracted $ANO-R$ basis
set [43, 44, 45, 46] with (26$s$, 20$p$, 15$d$, 10$f$, 4$g$) GFs; ii) a
modified $ANO-R$ basis set obtained by uncontraction of the previous primitive
$ANO-R$ functions to (26$s$, 26$p$, 20$d$, 15$f$, 10$g$) GFs; iii) tempered
GFs, whose exponent follow a geometric progression; iv) $aug-cc-pVTZ-PP$ and
$cc-pVTZ-PP$ GF basis sets [47, 48], which are two-component relativistic
pseudo-potentials (PP) optimized to calculate ground and excited state
spectroscopic properties. We notice that the procedure to uncontract the
$ANO-R$ basis set is carried out to obtain GFs optimized also for the small
component of the Dirac spinor and not only for the large one (for which this
basis set is generated). Indeed, we have found out empirically that
uncontracted basis sets achieve better accuracy in fully relativistic
simulations. The calculations performed with both $ANO-R$ and tempered GF
basis sets are all-electron, while those labelled by $-PP$ are characterised
by semi-local PPs that mimic the ionic core-valence electrons interaction to
achieve high computational efficiency in DHF calculations. The number of
electrons for both Au and W included in the PP core is 60, and the remaining
electrons are treated explicitly. The semi-local pseudo-potential reads [48]:

$V_{PP}=\frac{Q}{r}+\sum_{l,j}V_{lj}(r){\cal{P}}_{lj}=\frac{Q}{r}+\sum_{ij}\sum_{k}B_{lj,k}e^{-\beta_{lj,k}r^{2}}$
(126)

where $\frac{Q}{r}$ represents the long-range part of the PP carrying the
charge $Q$, while $V_{lj}(r)$ is the short-range radial (spherical-symmetric)
potential characterised by the parameters $B_{lj,k},\beta_{lj,k}$ that are
obtained by fitting with ab initio all-electron atomic data. The projectors
${\cal{P}}_{lj}$ in equation (126) enforce the condition that each core-shell
characterised by a different combination of the angular-momentum quantum
numbers $l$ and $j=l\pm 1/2$ is included in the PP with different parameters.
Following equation (126), PP can be written as linear combination of spherical
symmetric GFs acting on functional spaces of different symmetry $g=(l,j)$
(e.g. $s_{1/2},p_{1/2},p_{3/2},d_{3/2},d_{5/2}\cdots$), that is:

$\tilde{V}_{g}\left(r\right)=\sum_{i}\tilde{c}_{i,g}e^{-\alpha_{i}r^{2}}$
(127)

We remind that in a fully relativistic picture, shells characterised by the
same orbital angular momentum $l$ differ in energy owing to the spin-orbit
interaction.
Since $l\cdot s=\frac{j(j+1)-l(l+1)-s(s+1)}{2}$, from a semi-relativistic
quantum formulation in spherical symmetry one has:

$\tilde{V}_{g}\left(r\right)=V_{l}\left(r\right)+\frac{1}{2m_{e}^{2}c^{2}}\frac{1}{r}\frac{dV_{so,l\pm
1}}{dr}\frac{\left(l\pm\frac{1}{2}\right)\left(l\pm\frac{1}{2}+1\right)-l(l+1)-\frac{1}{2}\left(\frac{1}{2}+1\right)}{2}$
(128)

which can be specified for the different shells ($l>0$) as follows:

$\displaystyle\tilde{V}_{l,l-1/2}\left(r\right)$
$\displaystyle=V_{l}\left(r\right)-\frac{\left(l+1\right)}{4m_{e}^{2}c^{2}}\frac{1}{r}\frac{dV_{so,l-1}}{dr}$
$\displaystyle\tilde{V}_{l,l+1/2}\left(r\right)$
$\displaystyle=V_{l}\left(r\right)+\frac{l}{4m_{e}^{2}c^{2}}\frac{1}{r}\frac{dV_{so,l+1}}{dr}$
$\displaystyle\tilde{V}_{l+1,l+1/2}\left(r\right)$
$\displaystyle=V_{l+1}\left(r\right)-\frac{\left(l+2\right)}{4m_{e}^{2}c^{2}}\frac{1}{r}\frac{dV_{so,l}}{dr}$
$\displaystyle\tilde{V}_{l-1,l-1/2}\left(r\right)$
$\displaystyle=V_{l-1}\left(r\right)+\frac{\left(l-1\right)}{4m_{e}^{2}c^{2}}\frac{1}{r}\frac{dV_{so,l}}{dr}$
(129)

By subtracting the second and first equation above, and by summing the third
and the fourth, one obtains the following relations:

$\displaystyle\frac{1}{4m_{e}^{2}c^{2}}\frac{1}{r}\frac{d\left[lV_{so,l+1}+\left(l+1\right)V_{so,l-1}\right]}{dr}=\tilde{V}_{l,l+1/2}\left(r\right)-\tilde{V}_{l,l-1/2}\left(r\right)$
$\displaystyle\left(l-1\right)V_{l+1}\left(r\right)+\left(l+2\right)V_{l-1}\left(r\right)=\left(l-1\right)\tilde{V}_{l+1,l+1/2}\left(r\right)+\left(l+2\right)\tilde{V}_{l-1,l-1/2}\left(r\right)$
(130)

Since for a given $l$, the scalar $\mathcal{V}_{s}$ and vector
$\mathcal{V}_{v}$ relativistic potentials are related to $V_{l}$ e $V_{so,l}$
by (see the upper and lower diagonal elements of the Dirac equation (53) in
matrix form):

$\displaystyle\mathcal{V}_{v,l}+\mathcal{V}_{s,l}$ $\displaystyle=V_{l}$
$\displaystyle\mathcal{V}_{v,l}-\mathcal{V}_{s,l}$ $\displaystyle=V_{so,l}$
(131)

one has:

$3V_{0}=3\tilde{V}_{s_{1/2}},\quad
V_{3}+4V_{1}=\tilde{V}_{f_{5/2}}+4\tilde{V}_{p_{3/2}},\quad
2V_{4}+5V_{2}=2\tilde{V}_{g_{7/2}}+5\tilde{V}_{d_{5/2}},\quad
6V_{3}=6\tilde{V}_{f_{7/2}},\quad 7V_{4}=7\tilde{V}_{g_{9/2}}$ (132)

Furthermore, by defining $\tilde{I}_{l,l\pm 1/2}=4m_{e}^{2}c^{2}\int
r\tilde{V}_{l,l\pm 1/2}dr$ we obtain:

$\displaystyle
V_{so,2}+2V_{so,0}=\tilde{I}_{p_{3/2}}-\tilde{I}_{p_{1/2}},\quad
2V_{so,3}+3V_{so,1}=\tilde{I}_{d_{5/2}}-\tilde{I}_{d_{3/2}},\quad
3V_{so,4}+4V_{so,2}=\tilde{I}_{f_{7/2}}-\tilde{I}_{f_{5/2}}$ $\displaystyle
5V_{so,3}=\tilde{I}_{g_{9/2}}-\tilde{I}_{g_{7/2}},\quad 6V_{so,4}=0$ (133)

and, finally,

$V_{so,0}=0,\quad
3V_{so,1}=\tilde{I}_{d_{5/2}}-\tilde{I}_{d_{3/2}}-\frac{2}{5}\left(\tilde{I}_{g_{9/2}}-\tilde{I}_{g_{7/2}}\right),\quad
4V_{so,2}=\tilde{I}_{f_{7/2}}-\tilde{I}_{f_{5/2}}$ (134)

Level | Gaussian | Radial grid
---|---|---
| $HF$ | $DHF\,\,(ANO-R)$ |  $DHF\,\,(ANO-R)\,\,$mod. | Tempered |  $aug-cc-pVTZ-PP$ |  $cc-pVTZ-PP$ | $Loc.\,\,exch.$
$1s_{1/2}$ | $-2707.111$ | $-2960.001$ | $-2957.496$ | $-2957.543$ |  |  | $-2998.986$
$2s_{1/2}$ | $-457.459$ | $-532.437$ | $-529.119$ | $-529.566$ |  |  | $-534.299$
$2p_{1/2}$ |  21.2cm$-439.183$ | $-504.453$ | $-501.095$ | $-501.523$ |  |  | $-514.560$
$2p_{3/2}$ |  | $-443.371$ | $-440.020$ | $-440.443$ |  |  | $-445.810$
$3s_{1/2}$ | $-109.538$ | $-131.239$ | $-127.760$ | $-128.038$ |  |  | $-128.102$
$3p_{1/2}$ |  21.2cm$-100.900$ | $-119.602$ | $-116.176$ | $-116.417$ |  |  | $-118.838$
$3p_{3/2}$ |  | $-106.277$ | $-102.803$ | $-103.002$ |  |  | $-103.588$
$3d_{3/2}$ |  21.2cm$~{}-84.996$ | $-88.911$ | $-85.573$ | $-85.898$ |  |  | $-87.770$
$3d_{5/2}$ |  | $-86.498$ | $-82.958$ | $-83.272$ |  |  | $-84.430$
$4s_{1/2}$ | $-24.630$ | $-31.421$ | $-29.315$ | $-29.289$ |  |  | $-29.070$
$4p_{1/2}$ |  21.2cm$~{}-20.848$ | $-26.476$ | $-24.487$ | $-24.418$ |  |  | $-24.990$
$4p_{3/2}$ |  | $-23.285$ | $-21.383$ | $-21.312$ |  |  | $-21.292$
$4d_{3/2}$ |  21.2cm$~{}-13.933$ | $-15.721$ | $-14.032$ | $-14.039$ |  |  | $-14.459$
$4d_{5/2}$ |  | $-15.331$ | $-13.505$ | $-13.514$ |  |  | $-13.754$
$4f_{5/2}$ |  21.2cm$~{}-4.606$ | $-5.158$ | $-4.157$ | $-4.317$ |  |  | $-4.837$
$4f_{7/2}$ |  | $-1.036$ | $-4.069$ | $-4.224$ |  |  | $-4.677$
$5s_{1/2}$ | $-4.086$ | $-5.310$ | $-4.970$ | $-4.954$ | $-4.990$ | $-4.991$ | $-5.070$
$5p_{1/2}$ |  21.2cm$~{}-2.815$ | $-3.659$ | $-3.360$ | $-3.333$ | $-3.182$ | $-3.192$ | $-3.697$
$5p_{3/2}$ |  | $-3.075$ | $-2.897$ | $-2.879$ | $-2.898$ | $-2.899$ | $-3.040$
$5d_{3/2}$ |  21.2cm$~{}-0.801$ | $-0.897$ | $-0.791$ | $-0.791$ | $-0.795$ | $-0.792$ | $-1.120$
$5d_{5/2}$ |  | $-0.775$ | $-0.746$ | $-0.748$ | $-0.746$ | $-0.745$ | $-1.044$
Total energy |  | $-18941.867$ | $-18977.683$ | $-18961.916$ |  |  |
Table 1: Mean-field energy levels (Koopmans) of the Au+ ion obtained by a non-
relativistic approach ($HF$), a fully quantum mechanical relativistic approach
($DHF$) using $ANO-R$ contracted ($DHF\,\,(ANO-R)$), $ANO-R$ uncontracted
($DHF\,\,(ANO-R)\,\,$ mod.), and tempered Gaussian basis sets. Third and
second-last columns report the relativistic values of $\epsilon_{i}$ using the
pseudopotential method with the aug-cc-PVTZ-PP uncontracted and cc-PVTZ-PP
uncontracted Gaussian basis sets for the valence electrons [48]. Last column
reports the results obtained by using a radial grid with the local-exchange
interaction (see equation 85). Data in a.u.

The coefficients $\left({c}_{i,g},\alpha_{i}\right)$ of equation (127) (or
$\left(B_{lj,k},\beta_{lj,k}\right)$ of equation 126) are tabulated in Ref.
[48], and so the scalar and vector relativistic potentials (4.1) can be
assessed from equations (132-134). To show the impact of relativistic effects
on the energy levels (and ionization potential) we also carried out the mean-
field calculations using the non-relativistic approximation (standard HF from
the Schrödinger picture).

Figure 1: Left: relativistic $5d_{3/2}$ vs. non relativistic $5d$ valence
wavefunctions (bottom panel, in green the large part of the Dirac spinor, in
orange the small one) and densities (top panel, in orange color the DHF
results while in cyan the HF calculations) of Au+ represented in the radial
grid. The non-relativistic $5d$ density is similar to the relativistic
$5d_{3/2}$ density. Right: relativistic $6s_{1/2}$ vs. non relativistic $6s$
valence wavefunctions (bottom panel, in green the large part of the Dirac
spinor, in orange the small one) and densities (top panel, in orange color the
DHF results while in cyan the HF calculations) of Au+. The non relativistic
$6s$ density differs from the relativistic $6s_{1/2}$ density for a small
shift at larger distances, but shows similar oscillations. Figure 2: Optimized
molecular structure of the W@Au12 cage. Figure 3: a) Pseudo-atomic orbitals
energy-level diagram of the W@Au12 nanocage from DHF calculations. b) Virtual
state density at -12.7 eV. c) HOMO density. d) LUMO density. The coordinate
(0,0) from which the radial distance is measured corresponds to the position
of the W atom.

In Tab. 1 we report the values of the independent-particle eigenvalues
$\epsilon_{i}$ of the mean-field operator at different levels of
approximations, using either different Gaussian basis sets or the radial grid
method. We stress that within the Gaussian basis set approach, we included in
the relativistic Hamiltonian also the Gaunt term [1] for unretarded electron
interactions. The latter was at variance neglected in the radial grid
calculations. We notice that the DHF total energy is the lowest for the
$ANO-R$ modified basis set.

Level | Gaussian | Radial grid
---|---|---
| $HF$ | $DHF\,\,(ANO-R)$ |  $DHF\,\,(ANO-R)\,\,$mod. |  $aug-cc-pVTZ-PP$ | $Loc.\,\,exch.$
$1s_{1/2}$ | $-2362.726$ | $-2568.121$ | $-2563.117$ |  | $-2582.223$
$2s_{1/2}$ | $-393.356$ | $-453.526$ | $-447.916$ |  | $-450.206$
$2p_{1/2}$ |  21.2cm$-376.458$ | $-432.118$ | $-426.482$ |  | $-432.437$
$2p_{3/2}$ |  | $-382.822$ | $-377.199$ |  | $-381.707$
$3s_{1/2}$ | $-91.813$ | $-111.215$ | $-105.459$ |  | $-105.165$
$3p_{1/2}$ |  21.2cm$-83.950$ | $-102.003$ | $-96.163$ |  | $-96.916$
$3p_{3/2}$ |  | $-91.025$ | $-85.334$ |  | $-85.875$
$3d_{3/2}$ |  21.2cm$~{}-69.482$ | $-76.124$ | $-70.260$ |  | $-71.626$
$3d_{5/2}$ |  | $-73.575$ | $-67.943$ |  | $-69.174$
$4s_{1/2}$ | $-19.489$ | $-27.810$ | $-22.934$ |  | $-29.070$
$4p_{1/2}$ |  21.2cm$~{}-16.146$ | $-24.001$ | $-19.078$ |  | $-22.534$
$4p_{3/2}$ |  | $-21.207$ | $-16.535$ |  | $-16.421$
$4d_{3/2}$ |  21.2cm$~{}-10.078$ | $-15.032$ | $-10.284$ |  | $-10.479$
$4d_{5/2}$ |  | $-14.275$ | $-9.794$ |  | $-9.992$
$4f_{5/2}$ |  21.2cm$~{}-2.094$ | $-6.566$ | $-1.880$ |  | $-2.379$
$4f_{7/2}$ |  | $-3.078$ | $-1.784$ |  | $-2.281$
$5s_{1/2}$ | $-2.807$ | $/$ | $-3.433$ | $-3.259$ | $-3.515$
$5p_{1/2}$ |  21.2cm$~{}-1.746$ | $/$ | $-2.272$ | $-2.391$ | $-2.413$
$5p_{3/2}$ |  | $/$ | $-1.773$ | $-1.932$ | $-2.008$
$5d_{3/2}$ |  21.2cm$~{}-0.107$ | $/$ | $-0.350$ | $-0.378$ | $-0.528$
$5d_{5/2}$ |  | $/$ | $/$ | $/$ | $/$
Total energy | $-15286.935$ |  | $-16139.352$ |  |
Table 2: Mean-field energy levels (Koopmans) of the W atoms obtained by a non-
relativistic approach ($HF$), a fully quantum mechanical relativistic approach
($DHF$) using $ANO-R$ contracted ($DHF\,\,(ANO-R)$), $ANO-R$ uncontracted
($DHF\,\,(ANO-R)\,\,$ mod.) Gaussian basis sets. Third and second-last columns
report the relativistic values of $\epsilon_{i}$ using the pseudopotential
method with the aug-cc-PVTZ-PP uncontracted and cc-PVTZ-PP uncontracted
Gaussian basis sets for the valence electrons [48]. Last column reports the
results obtained by using a radial grid with the local-exchange interaction
(see equation 85). Data in a.u.

Furthermore, in Fig. (1) we plot the wavefuntions (bottom panels) and the
densities (top panels) for the $5d_{3/2}$ (left) and $6s_{1/2}$ (right)
states, respectively, obtained by using a radial mesh of 2000 points. By
comparing the non-relativistic HF calculations (cyan lines) with the DHF
results (orange lines), we notice that the non-relativistic $5d$ density is
similar to the relativistic $5d_{3/2}$ density.
A similar calculation was performed in the case of tungsten, which is the
foreign atom hosted in the gold cage. In Tab. 2 we report the single-level
energies in a.u. using either different Gaussian basis sets or the radial
mesh.

###  4.2 The W@Au12 superatomic system

The atomic positions of the W@Au12 nanocluster have been optimized by using
the VASP program suite [49, 50]. To find the equilibrium structure of W@Au12
we used the PBE-PAW pseudo-potential [51, 52, 53] with a plane-wave cut-off
equal to 250 eV to treat the ion-valence Coulomb interaction, and a Gaussian
smearing of the partial occupancies with a width $\sigma=0.05$ eV. The
simulation supercell is a cubic box with 15 Å side, which is large enough to
avoid spurious interactions between periodic images (long-range van der Waals
forces are not present in this calculation, which makes it equivalent to an
infinite box with an affordable computational cost). The Brillouin zone was
sampled only at the $\Gamma$ point. The equilibrium geometry of the W@Au12
nanocage is represented in Figure (2). W@Au12 presents a icosahedral molecular
geometry, having W at the center of the cage with W to Au equilibrium bond
distance equal to 2.8 Å and Au-Au 2.95 Å.
The stability of the W@Au12 nanocage, and in general of gold clusters, is
improved by encapsulating hetero-atoms owing to increased intermolecular
metallophilic interactions. Several theoretical and experimental studies have
been carried out on gold nano-cages hosting foreign atoms [9, 10, 42, 54, 55,
56]. The existence of the icosahedral W@Au12 cluster was first theoretically
predicted [57] and soon after successfully synthesized under gas-phase
conditions [9]. Indeed the Au12-cage itself is unstable, thus its stability
was attributed to the aurophilic attraction, in turns generated by the strong
relativistic effects in gold complexes, and the 18-valence electron rule (1
electron from each of the twelve 6$s^{1}$ orbitals of Au, and 6 unpaired
electrons from the $6s^{2}5d^{4}$ shells of W). Indeed, the electronic
structure of this compact nanocage is characterised by a closed-shell
configuration. In our DHF calculation 19 valence electrons for each Au atom
and 14 electrons for the W atom, populating the $n$=5,6 orbitals, have been
explicitly included in the model Hamiltonian, while the remaining ones (all
orbitals up to $n$=4) are treated via the frozen-core pseudo-potentials [48],
as outlined in the previous section. The choice of a small-core PP, where all
5$spd$ orbitals of Au hosting 19 electrons (14 for W) with a similar $\langle
r\rangle$ are considered explicitly as valence space, increases the accuracy
of the calculations due to a better treatment of the electron correlation
using a reduced “non-chemical” core [48]. To perform the DHF calculations we
used $aug-cc-PVTZ-PP$ Gaussian basis sets, including triple-zeta plus
polarization functions ($PTZ$), for the valence orbitals of all the atoms. In
Fig. 3a) we report the energy-levels diagram and, furthermore, a plot of the
density (square of the wavefunctions) for several levels, including a virtual
state (Fig. 3b), the HOMO (Fig. 3c), and the LUMO (Fig. 3d)). We notice that
the HOMO has a character mainly derived by the Au$6s$ and $5d$ (two major
lobes) orbitals with a small contribution from W (central lobe), while the
LUMO is composed primarily of W$5d$ orbitals (central lobe). The virtual state
is a mixture of Au$ns$ and $nd$ orbitals delocalized all over the cage. The
HOMO-LUMO gap from our DHF calculations is around 6 eV, which is an
overestimation of the value of 3 eV obtained from previous relativistic
density functional calculations using the generalized-gradient Perdew-Wang
exchange-correlation functional [9] and B3LYP/LANL1DZ [58] with a similar
value for the electron affinity. It is well known that the values of these
observables are typically overstimated using the HF approximation.

##  5 Application to the $\beta$-decay of heavy atoms

###  5.1 General theory of $\beta$-decay

$\beta$-decays occur in unstable nuclei, whereby an excess of neutrons
(protons) leads to the conversion into a proton (neutron) with the emission of
an electron (positron) and an electron antineutrino (neutrino). A typical
$\beta^{-}$-decay process is the following:

${}^{A}_{Z}X_{N}\to{}^{A}_{Z+1}X^{{}^{\prime}}_{N-1}+e^{-}+\bar{\nu}_{e}$
(135)

where the energy released upon decay ($Q$-value), defined as the difference
between the initial and final nuclear mass energies, determines whether the
transition can (or cannot) occur. In particular, a positive $Q$-value means
that the transition is energetically allowed. Other constraints are related to
the conservation of the electric charge, lepton number and baryon number of
the initial parent nucleus (${}^{A}_{Z}X_{N}$) and of the reaction final
products (${}^{A}_{Z+1}X^{{}^{\prime}}_{N-1}+e^{-}+\bar{\nu}_{e}$). In the
traditional approach, beta-decay spectra of allowed and forbidden unique
transitions are interpreted by using the following expression [59]:

$\frac{dN}{dW}\propto pWq^{2}F(Z,W)C(W)$ (136)

where

  * •

$pWq^{2}$ is a phase-space factor that accounts for both the relativistic
total kinetic energy ($W$, including thus the electron rest mass) and the
momentum sharing between the $\beta$-electron ($p$) and (anti-)neutrino ($q$);

  * •

a Fermi function $F(Z,W)$ that accounts for the static corrections due to the
Coulomb field of the nucleus;

  * •

a shape factor $C(W)$ that accounts for the coupling between the nuclear and
lepton dynamics. Typically, this factor reads:

$C(W)=(2L^{\prime}-1)!\sum_{k=1}^{L^{\prime}}\lambda_{k}\frac{p^{2(k-1)}q^{2(L^{\prime}-k)}}{(2k-1)![2(L^{\prime}-k)+1]!}$
(137)

where
$\lambda_{k}=\frac{(\alpha_{-k}^{2}+\alpha_{k}^{2})}{\alpha_{-1}^{2}+\alpha_{1}^{2}}$,
$\alpha_{k}$s are the Coulomb amplitudes of the electron wave functions, and
$L^{\prime}$ is related to the order of forbiddance (e.g. $L^{\prime}=1$ if
$\Delta J=0$ for an allowed transition or $L^{\prime}=\Delta J$ for any
($L^{\prime}-1$)th forbidden unique transition). $C(W)$ is assumed independent
of the nuclear structure details of the decaying atom [60]. However, nuclear
structure effects cannot be neglected when dealing with forbidden non-unique
transitions [61], and there is no such a simple relation for $C(W)$ in that
case as given by Equation (137).

To simplify Equation (137) one generally assume all $\lambda_{k}=1$. Despite
convincing theoretical arguments [62], a rigorous approach to model forbidden
non-unique transitions is still missing [63]. Our relativistic approach can
help to unravel the physics of beta-decay of these transitions.
In our model, the $\beta$-decay rate from an initial state $|i>$ to a final
state $|f>$ can be described by using the Fermi golden rule as:

$N_{i\rightarrow f}=2\pi\frac{\left|\left\langle
i|V|f\right\rangle\right|^{2}}{\left\langle i|i\right\rangle\left\langle
f|f\right\rangle}\delta(E_{i}-E_{f})$ (138)

where

$\hat{V}=\int
d{\bf{r}}\frac{G_{F}}{\sqrt{2}}\left[\hat{\bar{\psi}}_{p}({\bf{r}})\gamma^{\mu}\left(1-x\gamma^{5}\right)\hat{\psi_{n}}(\bf{r})\right]\left[\hat{\bar{\psi}}_{e}(\bf{r})\gamma_{\mu}\left(1-\gamma^{5}\right)\hat{\psi_{\nu}}(\bf{r})\right]=\int
d{\bf{r}}V(\bf{r})$ (139)

is the phenomenological weak interaction of the Standard Model of particles
[64] with $x=1.26\pm 0.02$ and $G_{F}=1.16637\times 10^{-5}GeV^{-2}$,
$\hat{\psi_{p}}(\bf{r})$, $\hat{\psi_{n}}(\bf{r})$, $\hat{\psi_{e}}(\bf{r})$,
$\hat{\psi_{\nu}}(\bf{r})$ are the field operators that annhilate a proton, a
neutron, and electron, and a neutrino, respectively, and
$\gamma_{\mu},\gamma^{\mu}$ are the Dirac matrices
($\hat{\bar{\psi}}=\hat{\psi}^{+}\gamma_{0}$) [1]. In particular, we are
interested in reckoning the probability per unit time that the atomic system
decays from a statistical mixture of initial states
$\hat{\rho}_{i}=p_{i}|i><i|$ to a mixture of final states
$\hat{\rho}_{f}=p_{f}|f><f|$, that is:

$N_{i\rightarrow
f}=2\pi\mathrm{Tr}(\hat{\rho}_{i}\hat{V}P_{f}\hat{V})\delta(E_{i}-E_{f})+h.c.$
(140)

where $P_{f}=\sum_{f}|f><f|$ is the projector onto the final states. In the
typical approximation to the general theory of $\beta$-decay the initial and
final states can be written as:

$\displaystyle|i>=|h_{i}>\otimes|e_{i}>$ (141)
$\displaystyle|f>=|h_{f}>\otimes|e_{f}>\otimes|\bar{\nu}_{f}>$ (142)

where $|h_{i,f}>$ are initial and final multi-nucleon states, $|e_{i,f}>$ are
initial and final multi-electron states, and $|\bar{\nu}_{f}>$ is the final
anti-neutrino state. We point out that the initial and final multi-nucleon
states and the initial multi-electron states are characterised by a discrete
spectrum, while the final multi-electron state, which describes the
$\beta$-emission, is of course a continuum state that can be written as linear
combination of external product $|e_{f}>=\sum_{j}I_{j,f}\wedge|\eta_{j,f}>$
where $|\eta_{j,f}>$ describes the $\beta$-electron continuum wavefunction.
Using this, the evaluation of the transition operator

$V_{fi}=\int d^{3}r\left\langle
f|\frac{G_{F}}{\sqrt{2}}\left[\bar{\psi}_{p}(\vec{r})\gamma^{\mu}\left(1-x\gamma^{5}\right)\psi_{n}(\vec{r})\right]\left[\bar{\psi}_{e}(\vec{r})\gamma_{\mu}\left(1-\gamma^{5}\right)\psi_{\nu}(\vec{r})\right]|i\right\rangle$
(143)

can be simplified to:

$\displaystyle V_{fi}$ $\displaystyle\simeq$ $\displaystyle\int
d^{3}r\frac{G_{F}}{\sqrt{2}}\left\langle
f_{h}|\left[\bar{\psi}_{p}(\vec{r})\gamma^{\mu}\left(1-x\gamma^{5}\right)\psi_{n}(\vec{r})\right]|i_{h}\right\rangle\left\langle
f_{l}|\left[\bar{\psi}_{e}(\vec{r})\gamma_{\mu}\left(1-\gamma^{5}\right)\psi_{\nu}(\vec{r})\right]|i_{l}\right\rangle$
(144) $\displaystyle=$ $\displaystyle\int
d^{3}r\frac{G_{F}}{\sqrt{2}}J_{i\rightarrow f}^{H,\mu}(\vec{r})J_{i\rightarrow
f,\mu}^{L}(\vec{r})$

that is, it can be factorized in the product of leptonic

$J_{i\rightarrow
f,\mu}^{L}(\vec{r})=\psi_{f,e}^{+}(\vec{r})\gamma_{0}\gamma_{\mu}\left(1-\gamma^{5}\right)\psi_{i,\nu}(\vec{r})$
(145)

and hadronic

$J_{i\rightarrow
f,\mu}^{H}(\vec{r})=\psi_{f,p}^{+}(\vec{r})\gamma_{0}\gamma^{\mu}\left(1-x\gamma^{5}\right)\psi_{i,n}(\vec{r})$
(146)

currents. This is ultimately due to the large rest mass of the $W$ boson,
which is the vector that mediates the weak interaction. We can safely assume
that electrons and neutrinos are not coupled and thus the leptonic current can
be further factorised into the independent product of the neutrino and
electron field operators (or wavefunctions). Furthermore, also the hadronic
current can be factorised into the product of neutron and proton wavefunctions
in systems where the mean-field approximation is expected to work rather well
[8] (such as in the case of the odd-even nuclei 63Ni, 129I, 241Pu). Odd-odd
nuclei, where many-body effects are expected to play an important role (such
as in the case of 36Cl, 210Bi, and 138La), are more complex to model and need
a more correlated approach to the nucleon-nucleon binding energy [8]. Indeed,
the separability of the hadronic current requires basically that only one
nucleon is involved in the decay (one neutron to one proton) and it is
independent of the remaining nucleons (which build the “core”, typically
approximated by a linear combination of angular momentum wavefunctions).
Therefore, in our approach the hadronic current is only partially correlated
via the presence of a core of nucleons that has the only goal of recovering
the total angular momentum of the parent reactant and of the final daughter
nucleus. However, this does not represent an intrinsic limit of our method, as
the hadronic current can be assessed via more accurate first-principles
approaches [65].

###  5.2 Ab-initio calculation of the hadronic and leptonic currents

In this work we consider a $\beta$-decaying system composed by an isolated
nucleus surrounded by atomic electrons (which can be partially stripped,
depending on temperature). Thus, we can safely assume to deal with a spherical
symmetric problem. Here we show how to reckon the fermion wavefunctions
appearing in the expression of the leptonic and hadronic currents, equations
(145) and (146), respectively.

####  5.2.1 Leptonic current

##### Neutrino wavefunction

The neutrino mass is negligible, therefore we can assume to deal with a free-
particle, whose energy-normalized wave functions

$\frac{u_{\kappa}(r)}{r}=k_{\nu}\sqrt{\frac{1}{\pi
c}}j_{l_{\kappa}}(k_{\nu}r),\qquad\frac{v_{\kappa}(r)}{r}=k_{\nu}S_{\kappa}\sqrt{\frac{1}{\pi
c}}j_{l_{-\kappa}}(k_{\nu}r)$ (147)

can be found by solving the free-particle Dirac equation (equation (87) with
$V=0$). For an emitted antineutrino ($\beta^{-}$ process), the energy-
normalized solution reads:

$\psi_{i,\nu}({\bf r})\equiv\psi_{\kappa_{\nu},m_{\nu}}^{-}({\bf
r})=\frac{1}{r}\left(\begin{array}[]{c}k_{\nu}\sqrt{\frac{1}{\pi
c}}j_{l_{\kappa_{\nu}}}(k_{\nu}r)\chi_{\kappa_{\nu},m_{\nu}}(\Omega)\\\
-iS_{\kappa_{\nu}}k_{\nu}\sqrt{\frac{1}{\pi
c}}j_{l_{-\kappa_{\nu}}}(k_{\nu}r)\chi_{-\kappa_{\nu},m_{\nu}}(\Omega)\end{array}\right)$
(148)

##### Electron (corrected) wavefunction

In a spherical symmetric atomic system the spinor wavefunctions $\psi$ can be
labelled by means of two azimutal quantum numbers $\kappa$ and $m$ as follows:

$\psi_{\kappa,m}({\bf
r})=\frac{1}{r}\left(\begin{array}[]{c}u_{\kappa}(r)\chi_{\kappa,m}(\Omega)\\\
iv_{\kappa}(r)\chi_{-\kappa,m}(\Omega)\end{array}\right)$ (149)

where $u_{\kappa}(r)$ and $v_{\kappa}(r)$ are radial functions and
$\chi_{\kappa,m}(\Omega)$ are the spin spherical harmonics. The solution of
the Dirac equation for the atomic electrons is thus simplified by the
spherical symmetry but made it cumbersome by the interaction with the other
electrons (in both the initial and in the final states). Therefore, several
level of approximations are in place.

  1. 1.

The electron wavefunction is found by solving self-consistently the DHF
equation in spherical symmetry (see equation 87) for electrons interacting via
a mean-field potential (see equation 88).

  2. 2.

Asymptotically, when the potential is negligible, they will behave as free
particles:

$\frac{u_{\kappa}(r)}{r}\simeq c_{u}\frac{\sin\left(kr-
l_{\kappa}\frac{\pi}{2}+\varphi\right)}{kr},\qquad\frac{v_{\kappa}(r)}{r}\simeq
c_{v}\frac{\sin\left(x-l_{-\kappa}\frac{\pi}{2}+\varphi\right)}{kr}$ (150)

with normalization condition given by

$\frac{\pi}{2k^{2}}\left(\left|c_{u}\right|^{2}+\left|c_{v}\right|^{2}\right)\delta\left(k-k^{\prime}\right)$
(151)

which requires that

$\left(\left|c_{u}\right|^{2}+\left|c_{v}\right|^{2}\right)=\frac{2k^{2}}{\pi}$
(152)

and

$c_{u}=k\sqrt{\frac{E+mc^{2}}{\pi E}}A\qquad,\qquad
c_{v}=kS_{\kappa}\sqrt{\frac{E-mc^{2}}{\pi E}}A$ (153)

  3. 3.

Inclusion of the non-orthogonality between the bound initial
$\varphi_{f_{i}}^{b}$ and final $\varphi_{f_{f}}^{b}$ orbitals, which are
obtained by solving the DHF equation (87) with two different atomic numbers.
The continuum electron wavefunction $\psi_{f,e}(\bf{r})$ is then calculated in
the field produced by the nucleus and the surrounding electrons and thus
modified by their presence. Finally, taking into account the Pauli’s principle
the bare Coulomb continuum wavefunction $\psi_{f}^{c}(\bf{r})$ corrected for
the non-orthogonality between initial and final states can be written as a
Slater determinant:

$\psi_{f,e}({\bf
r})=\left|\begin{array}[]{ccccc}\langle\varphi_{f,1}^{b}\big{|}\varphi_{i,1}^{b}\rangle&\langle\varphi_{f,1}^{b}\big{|}\varphi_{i,2}^{b}\rangle&\ldots&\langle\varphi_{f,1}^{b}\big{|}\varphi_{i,N}^{b}\rangle&\varphi_{f,1}^{b}(\bf{r})\\\
\langle\varphi_{f,2}^{b}\big{|}\varphi_{i,1}^{b}\rangle&\langle\varphi_{f,2}^{b}\big{|}\varphi_{i,2}^{b}\rangle&\ldots&\langle\varphi_{f,2}^{b}\big{|}\varphi_{i,N}^{b}\rangle&\varphi_{f,2}^{b}(\bf{r})\\\
\vdots&\vdots&\ddots&\vdots&\vdots\\\
\langle\varphi_{f,N}^{b}\big{|}\varphi_{i,1}^{b}\rangle&\langle\varphi_{f,N}^{b}\big{|}\varphi_{i,2}^{b}\rangle&\ldots&\langle\varphi_{f,N}^{b}\big{|}\varphi_{i,N}^{b}\rangle&\varphi_{f,N}^{b}(\bf{r})\\\
\langle\psi_{f}^{c}\big{|}\varphi_{i,1}^{b}\rangle&\langle\psi_{f}^{c}\big{|}\varphi_{i,2}^{b}\rangle&\ldots&\langle\psi_{f}^{c}\big{|}\varphi_{i,N}^{b}\rangle&\psi_{f}^{c}(\bf{r})\end{array}\right|$
(154)

  4. 4.

Inclusion of the shake-up and shake-off states of the final ion. Indeed, while
we expect that at room conditions the (both nuclear and electronic) ground-to-
ground $\beta$-decay gives a major contribution to the rate, to obtain the
total rate one has to include all possible electronic (and, in principle, also
nuclear) excitations that the final ion can undergo. Typically, the inclusion
of shake-up and shake-off has a sizable effect on the $\beta$-spectrum only at
low energy.

For a more detailed description of the calculation of the leptonic matrix
elements on a real-space grid we refer the reader to Ref. [8].

####  5.2.2 Hadronic current

##### Neutron and proton wavefunction calculations

While we expect that the effective shape of the $\beta$ spectra is almost
independent of the behaviour of the hadronic current, however explicit
numerical solutions of the Dirac equations for both protons and neutrons must
be found in order to asses the transition matrix elements (144). Assuming that
the recoil energy of the final nucleus is small with respect to their rest
mass, we can solve the Dirac equation (87), with a semi-empirical scalar and
vector relativistic Wood-Saxon potential to describe the nucleon-nucleon
interaction [40]. This potential represents essentially a screened Coulomb
interaction as follows:

$V_{C}(r)=-V_{C}\left[1+\exp{\left(\frac{r-R}{a}\right)}\right]^{-1}$ (155)

where $V_{C}$ represents the isospin-dependent total strength:

$V_{C}=V_{0}\left(1\pm\chi\frac{N-Z}{A}\right)$ (156)

and the upper sign is for protons and the lower one for neutrons, while $\chi$
is a constant to be optimized case-by-case. In equation (155) $R=R_{0}A^{1/3}$
($R_{0}$ is a constant, $A$ is the mass number) and $a$=constant are the size
of nuclear potential and diffuseness of the surface, respectively. However,
the nuclear force is also spin-dependent. Thus, one has to include the nuclear
spin-orbit contribution, which can be modelled by:

$\tilde{V}_{SO}(r)=\tilde{V}_{SO}\left[1+\exp{\left(\frac{r-R_{SO}}{a_{SO}}\right)}\right]^{-1}$
(157)

where $\tilde{V}_{SO}$ is the strength, while $R_{SO}=R_{0,SO}A^{1/3}$ and
$a_{SO}$=constant are the radius and the diffuseness of the spin-orbit term
(typically $a=a_{SO}$). The spin-orbit interaction strength is determined by
$\tilde{V}_{SO}=\lambda V_{C}$, where $\lambda$ is a proportionality constant.
We point out that the parameters $V_{0}$, $\chi$, $\lambda$, $a=a_{SO}$,
$R_{0},R_{0,SO}$ must be optimized by minimizing the difference between
experimental and calculated orbital energies. In particular, the values that
have been used in our simulations are the following $V_{0}=52.06~{}$ MeV,
$\chi=0.639$, $R_{0}=1.260~{}$ fm, $R_{0,SO}=1.160~{}$ fm, $\lambda=24.1$, and
$a=a_{SO}=0.662~{}$ fm.
For a more detailed description of the calculation of the hadronic current
term in a real-space radial grid we refer the reader to Ref. [8].

###  5.3 Test cases: the $\beta$-decay spectra of 31Si and 89Sr

As a test case of our relativistic approach to $\beta$-decay of heavy-mass
nuclei, we show in Fig. 4 a comparison between our model and the experimental
measurements recorded for the allowed $\beta$ transition of ${}^{31}_{14}$Si17
(left panel) and for the first forbidden unique decay of ${}^{89}_{38}$Sr51
(right panel). We stress that the agreement between experimental data [66] and
calculations is rather good.

Figure 4: Left panel: allowed $\beta$ spectrum of 31Si. Experimental data
(empty squares) from Ref. [66]. Right panel: First forbidden unique $\beta$
spectrum of 89Sr. Experimental data (empty squares) from Ref. [66].

The ${}^{31}_{14}$Si17 (odd neutron, even proton) nucleus $\beta^{-}$-decays
100% as follows [67]:

$^{31}_{14}\mathrm{Si}_{17}\to{}^{31}_{15}\mathrm{P}_{16}+e^{-}+\bar{\nu}_{e}$
(158)

with a ground state to ground (GS-GS) state $Q$-value (it first decays to the
first excited state of P) equal to 1491.50 keV and half-life of 157.36 minutes
at room conditions. The parent-nucleus ground state is characterized by the
quantum number $3/2^{+}$ ($J=3/2$ is the total angular momentum and $\pi=+$
the parity of the wavefunction). Protons are in a closed-shell configuration,
while one neutron is uncoupled: a mean-field single-particle approximation can
describe rather accurately its wavefunction. In the final state, before
undergoing a $\gamma$ transition to the ground state, the
${}^{31}_{15}\mathrm{P}_{16}$ nucleus is still in the $3/2^{+}$ configuration.
Thus, parity and total angular momentum are both conserved: this is an allowed
transition [8].

The ${}^{89}_{38}$Sr51 (odd neutron, even proton) nucleus $\beta^{-}$-decays
100% as follows [67]:

$^{89}_{38}\mathrm{Sr}_{51}\to{}^{89}_{39}\mathrm{Y}_{50}+e^{-}+\bar{\nu}_{e}$
(159)

with a ground state to ground state $Q$-value (it first decays to the first
excited state of Y) equal to 1500.9 keV and half-life of 50.563 days at room
conditions. The parent nucleus ground state is characterized by the quantum
numbers $J=5/2^{+}$. Protons are in a closed-shell configuration, while one
neutron is uncoupled: again mean-field single-particle approximation can
describe rather accurately its wavefunction. In the final state, before
undergoing a $\gamma$ transition to the ground state, the
${}^{89}_{39}\mathrm{Y}_{50}$ nucleus is in the $9/2^{+}$ configuration. Thus,
parity is conserved while the total angular momentum changes by two units:
this is first forbidden unique transition [8].

###  5.4 $\beta$-decay in astrophysical scenarios: the case of 134Cs

Our understanding of the cosmic chemical abundances is tightly bound to the
advancements in measuring or modelling the rates at which nuclear reactions,
and in particular $\beta$-decays, occur in stellar environment. However,
experimental set-up and procedures to reproduce the (ionized) plasma density
and temperature in astrophysical scenarios aiming to measure $\beta$-decay
rates are still in its infancy [68, 69]. Therefore, to model the isotopic
abundances of chemical elements in stars we must rely on theoretical modelling
[21, 70].
It is known that the abundance of Ba in evolved (AGB) stars is determined by a
chain of processes that include the slow neutron-capture process occurring in
the stable isotope 133Cs, whose presence is in turn affected by the
concentration of Xe in the solar system [71]. Here, we discuss specifically
the following $\beta$-decay:

${}_{55}^{134}\mathrm{Cs}~{}\rightarrow_{56}^{134}\mathrm{Ba}+e^{-}+\bar{\nu}$

${}^{134}_{55}$Cs decays to several nuclear excited states of
${}^{134}_{56}$Ba characterised by different total angular momentum. A
transition scheme is reported in the left panel of Fig. 5 [67].

Figure 5: Left panel: decay scheme of ${}^{134}_{55}$Cs. Data from [72]. Right
panel: half-life (years) in $log~{}x-y$ scale of ${}^{134}_{55}$Cs obtained by
using DHF with FD probability distribution for electrons (cyan line), without
including the decay to bound states (yellow line, no Fermi-Dirac temperature,
orbitals are occupied in their HF ground state). Brown line: Takahashi-Yokoi
(TY) calculations [73, 74] based on systematics [75]. Nuclear excited state
decays are included in all simulations. The electron density is 1026
electrons/cm3.

The GS-GS $Q$-value is 2058.7 keV and in Earth’s conditions the $\beta$-decay
occurs in 2.0652 years to the stable nucleus of ${}^{134}_{56}$Ba.
Additionally to the $4^{+}$ GS, several nuclear excited states of
${}^{134}_{55}$Cs can be identified, the lowest energy being $5^{+}$ (11 keV
above the GS) and possibly $3^{+}$ (60 keV above the GS, unsafe attribution).
Given the high temperature in stellar environment and assuming a degeneracy-
weighted Boltzmann distribution $p_{i}=e^{(-E_{i}/K_{\mathrm{B}}T)}$, where
$E_{i}$ is the energy of the nuclear level, and $T$ the temperature of the
stellar plasma ($1$ keV $\approx$ 11.6 MK), Cs nuclear excited states may be
populated and contribute to the decay. The excited state level population can
modify dramatically the isotopic half-life (typically reducing it by several
orders of magnitude).
The calculation of the $\beta$ decay has been carried out then by solving the
DHF equations for both the electron liquid (see section 5.2.1) and the nucleus
(see section 5.2.2), using the following approximations:

  * •

Our rate has been renormalized at all temperatures by a constant factor, which
has been obtained so as to recover the experimental $log(ft)$ in our room
temperature simulations [73, 74, 75].

  * •

We assumed that the neutron in the parent Cs nucleus weak decays from the
$2d_{3/2}$ shell to the $1g_{7/2}$ shell of daughter Ba (according to the
nuclear shell model).

  * •

The electronic levels of the Cs atom have not been re-optimized at each
temperature and are populated according the Fermi-Dirac (FD) distribution for
fermions
$F(T,\mu)=\frac{1}{1+e^{(\epsilon_{i}-\mu_{e^{-}})/(K_{\mathrm{B}}T)}}$, where
the energies $\epsilon_{i}$ of the levels $i$ are obtained by the self-
consistent solution of the DHF equation for electrons.

  * •

The chemical potentials of electrons and positrons are calculated as ideal
Fermi gases in a box using the relativistic energy-momentum dispersion
$E^{2}=c^{2}p^{2}+m_{e}c^{4}$ by the following relation:

$n_{p}=\int_{0}^{\infty}dp~{}p^{2}/\pi^{2}\times(F((c\times\sqrt{(}p^{2}+c^{2})-\mu_{e})/K_{\mathrm{B}}T)-F((c\times\sqrt{(}p^{2}+c^{2})+\mu_{e})/K_{\mathrm{B}}T))$
(160)

At variance, protons are non-relativistic particles. The density of electrons
in the charge-neutral plasma is taken 1026 electrons/cm3.

Following the decay scheme reported in the left panel of Fig. 5, we included
several possible decay paths from the $4^{+},5^{+},3^{+}$ nuclear states of
Cs, that is Cs$(4^{+})\rightarrow\mathrm{Ba}(4^{+},3^{+},4^{+})$;
$\mathrm{Cs}(5^{+})\rightarrow\mathrm{Ba}(4^{+},3^{+},4^{+})$, and
$\mathrm{Cs}(3^{+})\rightarrow\mathrm{Ba}(4^{+},3^{+},4^{+},2^{+},2^{+})$. In
the right panel of Fig. 5 we show the cumulative half-life, taking into
account both electron and nuclear excitations, using a FD and a degeneracy-
weighted Boltzmann distribution for electrons and nuclei, respectively.
We notice that the half-life of the ${}^{134}_{55}$Cs isotope is dramatically
affected by the presence of the three nuclear excited states of Cs (see the
fast decrease of the cyan line reported in the right panel of Fig. 5 above 10
keV, 1 keV=11.6 MK) as well as by the electronic temperature (compare the
yellow and cyan line shown in the right panel of Fig. 5 without and with
electronic temperature, respectively), which may lead to partial or total
ionization. We observe that the nuclear excited state dynamics, which takes
into account the population of fast-decaying nuclear excited states, is the
most relevant of the two concurrent effects, as it can decrease the rate by a
factor of 15 at 100 keV (1 GK) to 23 at 1000 keV with reference to room
temperature conditions. This effect is of course in place at higher
temperatures ($>$10 keV) according to the high nuclear excitation energies of
Cs (11 and 60 keV, respectively).

##  6 Electron elastic cross section in liquid water

Cancer treatment based on fast ion beams – the so-called hadrontherapy – is
much more efficient than conventional radiotherapy that uses photons as energy
carriers [32]. Indeed, in the former approach a substantial fraction of the
ion initial kinetic energy is deposited in a region near the end of their
trajectories – the so-called Bragg peak – while in the latter a majority of
the radiation dose is delivered by photons near the entrance and then smoothly
distributed along the path inside the body, irradiating also the healthy
tissues and causing undesired biodamage especially in the therapy of deep-
seated tumors [31].
In hadrontherapy the energy lost by the ions in their way within the medium
typically results in the emission of secondary electrons ejected via
ionizations of the constituents. Most of these secondary electrons are
produced with very low energies (below 100 eV) and specific angular pattern
[76, 77]. The assessment of the scattering cross sections for the different
collisional processes, such as excitation, ionization or elastic scattering
[78], initiated by the ion beam inside the living tissue is of paramount
importance to determine the generated secondary electron cascade, the
concentration of energy deposition and the damaging events in biotargets.
Typically human tissue is well mimicked by liquid water, which is its main
constituent. Here we show how our relativistic approach based on HGBF can help
our understanding of the elastic scattering in this medium.
Elastic scattering of electrons by atomic or molecular centers [79] can be
dealt either with the relativistic Mott theory [80] or with the direct self-
consistent solution of the Dirac equation.
The former approach solves the Dirac equation typically using either a semi-
empirical electron-electron interaction central potential or a best fit of
data from Hartree-Fock (HF) simulations [81]. Then, the differential elastic
scattering cross section – defined as the ratio of the number of particles
elastically scattered into the solid angle $d\Omega$ per unit time divided by
the incident flux – is obtained for an unpolarized electron beam by the
following expression [81]:

$\frac{d\sigma_{\mathrm{el}}}{d\Omega}\,=\,|f(\theta)|^{2}+|g(\theta)|^{2}$
(161)

where $f(\theta),g(\theta)$ are the direct and spin-flip scattering
amplitudes, and $\theta$ represents the scattering angle between the incident
and the emitted electron. Finally, by integrating over all solid angles one
obtains the total elastic scattering cross section:

$\sigma_{\mathrm{el}}(T)=\int_{\Omega}\frac{d\sigma_{\mathrm{el}}(T,\theta)}{d\Omega}{d\Omega}.$
(162)

and from the latter observable the elastic mean free path reads
$\lambda_{el}=\frac{1}{{\cal{N}}\sigma_{el}}$, where $\cal{N}$ is the target
atomic or molecular number density. Equation (161) is valid to model the
elastic scattering of an unpolarized beam of electrons impinging on a single
water molecule, which is the typical of the gas phase. To take into account
the presence of randomly distributed water molecules in the liquid phase
extensions have been proposed [81].
The second approach computes the differential and total elastic cross sections
by directly solving the DHF equation (53) or (equation (101) if HGBF are used)
for a single molecule or a cluster of water molecules. Water is a tri-atomic
system, thus the numerical solution will be searched by using our fully
relativistic quantum mechanical method with HGBF basis sets centered on the
atomic nuclei of the molecule(s). We project out the potential interaction
only, so as to recover the elastic continuum (see equations 103 and 104). In
particular, the electronic wavefunctions and the mean-field electron-electron
interaction potential were expanded using the aug-cc-pVTZ Gaussian functions
basis set optimized for both hydrogen and oxygen atoms [44]. Mono– and bi–
electronic molecular integrals of the bare Coulomb and exchange interaction
are computed at each SCF cycle as outlined in sections 3.2-3.5.
Once the scattering stationary states $|\psi_{\bf{k}}^{+}(E)\rangle$ of the
projected Dirac Hamiltonian (101) are obtained at energy $E$, the differential
elastic cross section for solid angle unit is finally reckoned via the Fermi
Golden rule as follows:

$\frac{\mbox{d}\sigma_{el}}{\mbox{d}\Omega}=\frac{m^{2}}{4\pi^{2}}|\langle\phi_{k\hat{n}}|T^{+}(E)|\phi_{\bf{k}}\rangle|^{2}=\frac{m^{2}}{4\pi^{2}}|\langle\phi_{k\hat{n}}|V|\psi_{\bf{k}}^{+}\rangle|^{2}$
(163)

where $\phi_{k\hat{n}}$ is the impinging electron wavefunction represented by
a plane-wave with momentum $k$ ($E=\frac{k^{2}}{2m}$) along the direction of
motion $\hat{n}$, $\phi_{\bf{k}}$ is the outgoing free plane wave elastically
scattered in the direction $\bf{k}$, and $T^{+}(E)$ is the on-shell $T$-matrix
defined by $T|\phi_{\bf{k}}=V|\psi_{\bf{k}}^{+}\rangle$. The scattering
wavefunction $\psi_{\textbf{k}}^{+}(E)$ is characterised by the so-called
outgoing $(+)$ wave boundary conditions, that is the eigensolutions of the
Dirac equation are matched asymptotically by Coulomb functions [3, 4, 8]
defining the normalization condition. When the electron kinetic energy is
large, the equation (163) can be simplified by using the first-order Born
approximation, whereby the $T$-operator expansion in powers of the molecular
relativistic potential $T=V+VG_{0}V+VG_{0}VG_{0}V+....$ ($V$ is obtained via
the self-consistent solution of the Dirac equation,
$G_{0}=\frac{1}{E-H_{0}-i\epsilon}$ is the free Green’s function, see equation
103) consists only of the first term, that is $T=V$. Adopting the first-order
Born approximation, the formula (163) is equivalent to the Mott cross section
(161). To assess the differential and total elastic cross section for
electrons moving within liquid water one should also consider the presence of
neighboring molecules. Thus, we trimmed a cluster of 6 water molecules after
carrying out molecular dynamics (MD) simulations with several thousands
molecules reproducing the experimental water density at room conditions
($\rho=1$ g/cm3), using the empirical TIP3P force-field [82] implemented in
the LAMMPS package [83]. The presence of a number of water molecules takes
into account the possibility of multiple scattering effects resulting in
possible destructive interference between the escaping waves.

Figure 6: Left panel: differential elastic scattering cross section (DESCS) of
an electron beam impinging with kinetic energy of 20 eV on liquid water as a
function of the scattering angle. Calculations from the DHF approach for a
cluster of 6 water molecules (see equation 163) is reported in solid red line
[78]. Calculations from the Mott approach for a single water molecule (see
equation 161) is reported in solid black line [78]. Symbols represent
experimental data recorded on water vapour from: Cho [84] (squares), Danjo
[85] (circles), Johnstone [86] (triangles) and Katase [87] (diamonds). Right
panel: total elastic scattering cross section (ESCS) as a function of the
incoming electron beam kinetic energy. Red line reports our ab-initio data for
a cluster of 6 water molecules, while black line is obtained by the Mott
theory with one water molecule. Experimental data from Refs. [87, 88].

In the left panel of Fig. 6 we report the differential elastic cross sections
(DESCS) of an electron beam incident on liquid water with kinetic energy equal
to 20 eV. Red and black solid lines correspond to calculations based on
equation (163) with a water molecule cluster and on equation (161) with a
single water molecule, respectively [78], while symbols represent experimental
data [84, 85, 86, 87]. In equation (161) the atomic potential consists of a
superposition of Yukawa functions, whose parameters were set according to a
best fit of data from Hartree-Fock simulations [89], while exchange effects
were taken into account via the Furness and McCarthy formula [90]. We notice
that particularly at small scattering angles the ab-initio calculations
performed by solving numerically the Dirac equation for a cluster of molecule
is in better agreement with experimental data (even though the latter
measurements are performed on water vapour). In the right panel we report a
comparison between the total elastic cross section obtained by the equation
(162) using the fully relativistic approach (red line) for the 6 water
molecule cluster, the Mott approach on single water molecule (black line), and
the experimental data. We notice that the condensed phase nature of liquid
water emerges as a significant deviation from the single water molecule
particularly at low energy, where the elastic cross section assessed on the
cluster is appreciably reduced with respect to the single molecule, stressing
the importance of including multiple scattering in the modelling.

##  7 Conclusions

This chapter is aimed at describing the relativistic theory of many-fermion
systems, which stems from the merging of the two most fundamental theory of
modern physics, that is quantum mechanics and the special theory of
relativity. In particular, we focused on the recent theoretical and
computational advances of a method for solving the DHF equation for
interacting relativistic fermions. This approach can make use of both radial
and Gaussian basis sets, which allows one to study any nuclear, atomic, and
molecular system with controllable and systematically improvable accuracy.
The attractiveness of the HGBF method is due to the the possibility of writing
analytical expressions of the integrals needed for computing the matrix
elements of the Coulomb potential for both bound and scattering states.
Furthermore, this approach can be also easily applied to multicentric systems,
such as molecules and clusters. We showed that HGFB sets, typically created
for bound states of quantum systems within the Schrödinger picture (thus, in
the limit of an infinitely large speed of light), should be modified by
uncontraction procedures to account properly for the four-dimensional nature
of electrons and nucleons, which ultimately leads to the spin-orbit
relativistic effects.
On the other hand, the spherical symmetry, typically found in atomic and
nuclear systems, allows one to formulate the many-fermion problem in the
radial coordinate only (in absence of symmetry-destroying fields), resulting
in the solution of a mono-dimensional Dirac equation. The choice between these
two computational schemes basically depends on the symmetry of the problem
under investigation and on the accuracy needed, keeping in mind that the
radial basis approach scales almost linearly with the number of mesh points
while the HGBF approach grows in between a cubic (due to matrix
diagonalization) and quartic (due to the bi-electronic integrals calculation)
manner as a function of the number of basis functions. In this regard, we
discussed a few applications of our method to both atomic and molecular
systems based on gold, where relativistic affects dictate the electronic
structure properties.
The computational effort increasing with system size is faced with the
implementation of a relativistic pseudopential approach, where the interaction
between ion cores and outer electrons is treated by replacing the core
electrons by a weaker pseudopotential that acts on a set of pseudo
wavefunctions with the aim to suppress the strong oscillations in the Coulomb
potential close to the nucleus that require several Gaussian functions to be
accurately described at an increased computational cost.
Furthermore, we showed the extension of our relativistic approach based on
these two complementary numerical solutions to deal with i) the elastic
continuum for describing the electron-molecule scattering; ii) the
$\beta$-decay emission from heavy elements in both Earth and astrophysical
scenarios by including the electron temperature via the FD distribution and
the nuclear temperature via the Boltzmann population.
We point out that in all the case studies discussed in this chapter the
description and the accuracy achievable with our numerical methods rely
essentially on the interplay between the so-called relativistic effects,
connected to the Dirac four-dimensional representation with respect to the
one-dimensional nature of the Schrödinger quantum mechanics, and the electron-
electron correlation.
Finally, desirable and prospective developments of our theoretical and
computational schemes will be sought along two different lines. On the
theoretical side, we plan to move from the solution of electronic structure
problems to the study of electron spectroscopy, most notably to deliver a
relativistic theory of the Auger effect in heavy-elements compounds within the
framework of the formal theory of scattering. Owing to the paramount
importance of the many-body interaction in Auger and resonant decays, this
goal requires to devise and develop tailor-made methods for treating dynamic
correlation effects beyond the HF approximation. Furthermore, to predict
quantitatively the spectral lineshapes one needs to account on the same level
of theory not only for the intrinsic dynamical properties of the perturbed
system, but also for the feature of the incident beam, be it x-ray photon or
electron beams. A further notable application of our approach, which will
require serious methodological development, will be in the field of nuclear
astrophysics, where nucleosynthesis models still show severe discrepancies
with the observations, possibly due to inaccurate (if not wrong) nuclear input
data (e.g. the cosmological abundance of 7Li). These experiments along with
our theoretical and computational tools will also be useful to support the
conclusions drawn so far by using Standard Model of particles or to search for
a new physics beyond the Standard Model (BSM).

##  8 Acknowledgments

S.S. and S.T. acknowledge funding from the Pandora project provided by the
National Institute of Nuclear Physics. S.T. acknowledges fruitful discussions
with F. Stella.

## References

  * [1] M. Rehher, A. Wolf, Relativistic Quantum Chemistry: The Fundamental Theory of Molecular Science, Wiley-VCH Verlag GmbH & Co. KGaA, 2009.  doi:10.1002/9783527627486.
  * [2] S. Taioli, S. Simonucci, M. Dapor, Surprises: when ab initio meets statistics in extended systems, Computational Science & Discovery 2 (1) (2009) 015002.
  * [3] S. Taioli, S. Simonucci, L. Calliari, M. Dapor, Electron spectroscopies and inelastic processes in nanoclusters and solids: theory and experiment, Physics Reports 493 (5) (2010) 237–319.
  * [4] S. Taioli, S. Simonucci, A computational perspective on multichannel scattering theory with applications to physical and nuclear chemistry, in: Annual Reports in Computational Chemistry, Vol. 11, Elsevier, 2015, pp. 191–310.
  * [5] S. Taioli, S. Simonucci, L. Calliari, M. Filippi, M. Dapor, Mixed ab initio quantum mechanical and monte carlo calculations of secondary emission from sio 2 nanoclusters, Physical Review B 79 (8) (2009) 085432.
  * [6] S. Taioli, S. Simonucci, The resonant and normal auger spectra of ozone, Symmetry 13 (3) (2021) 516.
  * [7] D. Tzeli, A. A. Tsekouras, Mind the basis set superposition error, Chemical Physics Letters 496 (1) (2010) 42–45.  doi:https://doi.org/10.1016/j.cplett.2010.07.053.
URL https://www.sciencedirect.com/science/article/pii/S0009261410009851

  * [8] T. Morresi, S. Taioli, S. Simonucci, Nuclear beta decay: Relativistic theory and ab initio simulations of electroweak decay spectra in medium-heavy nuclei and of atomic and molecular electronic structure, Advanced Theory and Simulations 1 (11) (2018) 1870030.
  * [9] X. Li, B. Kiran, J. Li, H.-J. Zhai, L.-S. Wang, Experimental observation and confirmation of icosahedral w@au12 and mo@au12 molecules, Angewandte Chemie International Edition 41 (24) (2002) 4786–4789.  arXiv:https://onlinelibrary.wiley.com/doi/pdf/10.1002/anie.200290048, doi:https://doi.org/10.1002/anie.200290048.
URL https://onlinelibrary.wiley.com/doi/abs/10.1002/anie.200290048

  * [10] L. Chen, Y. Gao, Y. Cheng, H. Li, Z. Wang, Z. Li, R.-Q. Zhang, Nonresonant chemical mechanism in surface-enhanced raman scattering of pyridine on m@au12 clusters, Nanoscale 8 (2016) 4086–4093.  doi:10.1039/C5NR07246H.
URL http://dx.doi.org/10.1039/C5NR07246H

  * [11] S. Jain, D. Hirst, J. O’Sullivan, Gold nanoparticles as novel agents for cancer therapy, Br. J. Radiol. 85 (2012) 101–13.
  * [12] E. Porcel, S. Liehn, H. Remita, N. Usami, K. Kobayashi, Y. Furusawa, C. L. Sech, S. Lacombe, Platinum nanoparticles: a promising material for future cancer therapy?, Nanotechnology 21 (8) (2010) 085103.  doi:10.1088/0957-4484/21/8/085103.
URL https://doi.org/10.1088/0957-4484/21/8/085103

  * [13] L. Sancey, F. Lux, S. Kotb, S. Roux, S. Dufort, A. Bianchi, Y. Crémillieux, P. Fries, J.-L. Coll, C. Rodriguez-Lafrasse, M. Janier, M. Dutreix, M. Barberi-Heyob, F. Boschetti, F. Denat, C. Louis, E. Porcel, S. Lacombe, G. Le Duc, E. Deutsch, J.-L. Perfettini, A. Detappe, C. Verry, R. Berbeco, K. T. Butterworth, S. J. McMahon, K. M. Prise, P. Perriat, O. Tillement, The use of theranostic gadolinium-based nanoprobes to improve radiotherapy efficacy, The British Journal of Radiology 87 (1041) (2014) 20140134, pMID: 24990037.  arXiv:https://doi.org/10.1259/bjr.20140134, doi:10.1259/bjr.20140134.
URL https://doi.org/10.1259/bjr.20140134

  * [14] J. C. Polf, L. F. Bronk, W. H. P. Driessen, W. Arap, R. Pasqualini, M. Gillin, Enhanced relative biological effectiveness of proton radiotherapy in tumor cells with internalized gold nanoparticles, Applied Physics Letters 98 (19) (2011) 193702.  arXiv:https://doi.org/10.1063/1.3589914, doi:10.1063/1.3589914.
URL https://doi.org/10.1063/1.3589914

  * [15] J.-K. Kim, S.-J. Seo, H.-T. Kim, K.-H. Kim, M.-H. Chung, K.-R. Kim, S.-J. Ye, Enhanced proton treatment in mouse tumors through proton irradiated nanoradiator effects on metallic nanoparticles, Physics in Medicine and Biology 57 (24) (2012) 8309–8323.  doi:10.1088/0031-9155/57/24/8309.
URL https://doi.org/10.1088/0031-9155/57/24/8309

  * [16] Z. Kuncic, S. Lacombe, Nanoparticle radio-enhancement: principles, progress and application to cancer treatment, Physics in Medicine & Biology 63 (2) (2018) 02TR01.  doi:10.1088/1361-6560/aa99ce.
URL https://doi.org/10.1088/1361-6560/aa99ce

  * [17] J. Schuemann, A. F. Bagley, R. Berbeco, K. Bromma, K. T. Butterworth, H. L. Byrne, B. D. Chithrani, S. H. Cho, J. R. Cook, V. Favaudon, Y. H. Gholami, E. Gargioni, J. F. Hainfeld, F. Hespeels, A.-C. Heuskin, U. M. Ibeh, Z. Kuncic, S. Kunjachan, S. Lacombe, S. Lucas, F. Lux, S. McMahon, D. Nevozhay, W. Ngwa, J. D. Payne, S. Penninckx, E. Porcel, K. M. Prise, H. Rabus, S. M. Ridwan, B. Rudek, L. Sanche, B. Singh, H. M. Smilowitz, K. V. Sokolov, S. Sridhar, Y. Stanishevskiy, W. Sung, O. Tillement, N. Virani, W. Yantasee, S. Krishnan, Roadmap for metal nanoparticles in radiation therapy: current status, translational challenges, and future directions, Physics in Medicine & Biology 65 (21) (2020) 21RM02.  doi:10.1088/1361-6560/ab9159.
URL https://doi.org/10.1088/1361-6560/ab9159

  * [18] I. Kempson, Mechanisms of nanoparticle radiosensitization, WIREs Nanomedicine and Nanobiotechnology 13 (1) (2021) e1656.  arXiv:https://onlinelibrary.wiley.com/doi/pdf/10.1002/wnan.1656, doi:https://doi.org/10.1002/wnan.1656.
URL https://onlinelibrary.wiley.com/doi/abs/10.1002/wnan.1656

  * [19] A. Briggs, S. Corde, S. Oktaria, R. Brown, A. Rosenfeld, M. Lerch, K. Konstantinov, M. Tehei, Cerium oxide nanoparticles: influence of the high-z component revealed on radioresistant 9l cell survival under x-ray irradiation, Nanomedicine: Nanotechnology, Biology and Medicine 9 (7) (2013) 1098–1105.  doi:https://doi.org/10.1016/j.nano.2013.02.008.
URL https://www.sciencedirect.com/science/article/pii/S1549963413000749

  * [20] R. Brown, M. Tehei, S. Oktaria, A. Briggs, C. Stewart, K. Konstantinov, A. Rosenfeld, S. Corde, M. Lerch, High-z nanostructured ceramics in radiotherapy: First evidence of ta2o5-induced dose enhancement on radioresistant cancer cells in an mv photon field, Particle & Particle Systems Characterization 31 (4) (2014) 500–505.  arXiv:https://onlinelibrary.wiley.com/doi/pdf/10.1002/ppsc.201300276, doi:https://doi.org/10.1002/ppsc.201300276.
URL https://onlinelibrary.wiley.com/doi/abs/10.1002/ppsc.201300276