# Green’s function and LDOS for non-relativistic electron pair

Tomasz M. Rusin  email<EMAIL_ADDRESS>Institute of Physics, Polish Academy of
Sciences, Al. Lotników 32/46, 02-688 Warsaw, Poland

###### Abstract

The Coulomb Green’s function (GF) for non-relativistic charged particle in
field of attractive Coulomb force is extended to describe the interaction of
two non-relativistic electrons through repulsive Coulomb forces. Closed-form
expressions for the GF, in the absence of electron spins, are derived as one-
dimensional integrals. The results are then generalized to include electron
spins and account for the Pauli exclusion principle. This leads to a final GF
composed of two components, one even and the other odd with respect to
exchange particles, with closed-form expressions represented as one-
dimensional integrals. The Dyson equations for spin-independent potentials is
presented. The local density of states (LDOS) is calculated, which is a
combination of contributions from both even and odd GFs. This calculation
reveals the dependence of LDOS on inter-electron distance and energy. Separate
analysis of the impact of the Pauli exclusion principle is provided. An
examination of the pseudo-LDOS, arising from the two-body contribution to the
Green’s function, is undertaken. Complete suppression of the LDOS at $r=0$ is
ensured by this term, which exhibits a restricted spatial extent. The reasons
for the emergence of this pseudo-LDOS are elucidated.

##  I Introduction

In 1963, Hostler and Pratt introduced a closed-form solution for the Green’s
function (GF) of a non-relativistic particle in the presence of an attractive
Coulomb potential [1]. This GF is expressed in terms of Whittaker functions
with complex arguments. The derivation of this result was the culmination of
extensive efforts, which began with Meixner’s work in 1933 [2] and involved
multiple approaches [3, 4]. Hostler subsequently re-derived the Coulomb GF
using various methods, yielding several equivalent expressions [5, 6, 7].
Hameka pursued a different approach and obtained the Coulomb GF in an
alternative form [8, 9]. Schwinger calculated the Coulomb GF in momentum space
[10], while Blinder derived it in parabolic coordinates [11], extending
Hostler’s approach to cover repulsive Coulomb potentials. Swainson and Drake
further extended Hostler’s results to the relativistic case [12]. For a
comprehensive review of papers related to the Coulomb GF and its applications
in multi-photon process calculations, Maquety et al. provide an insightful
review [13]. Furthermore, recent reviews covering various aspects of the
Coulomb GF can be found in [14]. In Refs. [15, 16, 17, 18] the two-particle GF
for the Coulomb problem was derived through a convolution of two distinct one-
electron Coulomb GFs. These derived results were subsequently employed for the
computation of Sturmian matrix elements. It is noteworthy, however, that a
comprehensive and detailed examination of the two-electron GF was not carried
out within the context of the aforementioned works.

In this paper, we extend the findings of Holster and Pratt [1] and Blinder
[11] to the case of two non-relativistic electrons interacting through Coulomb
forces. This generalization is possible because on can separate the motion of
the electron pair into center-of-mass and relative motions. However, several
complexities arise. Firstly, the electrons repel each other, preventing the
formation of any bound states. Moreover, the electron pair can exist in either
a singlet state or one of three triplet states, which adds to the intricacy.
Additionally, the Pauli exclusion principle introduces extra limitations on
the pair’s wave function.

In contrast to the simpler Coulomb problem, the GF for the electron pair
depends on four variables instead of the usual two, making the calculations
more challenging. Most notably, the GF for the electron pair does not neatly
separate into center-of-mass and relative motion parts, which adds to the
complexity.

In the subsequent sections of this paper, we outline our approach to
overcoming these challenges and obtaining the GF and the local density of
states (LDOS) in terms of quadratures of special functions, particularly
Whittaker and Gamma functions. This work offers a solution to these issues.

Our calculation of the GF for an electron pair and the LDOS unfolds in three
steps. In the first step, we neglect the electron’s spin and temporarily
setting aside the Pauli exclusion principle. This phase involves the
generalization of the Coulomb GF to encompass a pair of particles governed by
the two-particle Schrodinger equation, influenced by the repulsive Coulomb
potential. As we progress to the second step, we take into account the
presence of electron spin and the presence of the Pauli exclusion principle.
This step involves deriving the even and odd components of the pair’s GF,
preserving the requisite symmetries concerning the exchange of particles. In
the final step, we calculate the trace of the imaginary parts of both the even
and odd components of the pair’s GF. This leads us to the calculation of the
local densities of states, corresponding to the singlet and triplet states of
an electron pair, respectively.

The paper is organized as follows. In Section II, we calculate the GF while
disregarding the presence of electron spins. We derive the general expression
for the two-particle GF, conduct a thorough analysis of its properties,
identify specific sets of arguments that lead to GF divergence, and compute
the local density of states in the limit of vanishing arguments. Section III
extends the results from the previous section to encompass the scenario of two
electrons, including their spins. This section also takes into account the
limitations imposed by the Pauli exclusion principle. Section IV is dedicated
to the derivation of the Dyson equation for spin-independent potentials. In
Section V, we calculate the LDOS as a function of inter-electron distance and
the pair’s energy. In this section also analyze the pseudo-LDOS term,
resulting from the Pauli exclusion principle, which leads to the complete
vanishing of the odd-part of the LDOS at $r=0$. Section VI engages in
discussions on several issues pertaining to the obtained results, as well as
potential possibilities for their experimental observation. The paper is
summarized in the Summary, followed by an appendix that offers additional
information and explanations to complement the main text.

##  II Green’s function of two electrons in absence of spin effects

###  II.1 General form of GF

Consider two charged particles labeled ’a’ and ’b,’ positioned at coordinates
$\mathbf{a}$ and $\mathbf{b}$, and both possessing a common mass $m_{e}$
equivalent to the electron mass. These particles carry charges of $Z_{a}|e|$
and $Z_{b}|e|$, where $|e|$ denotes the elementary charge. It’s important to
note that in this scenario, we are assuming spinless particles, and as a
consequence, the two-particle wave function is not constrained by the Pauli
exclusion principle. The Hamiltonian describing the system is

$\hat{H}_{p}=-\frac{\hbar^{2}}{2m_{e}}\nabla^{2}_{\mathbf{a}}-\frac{\hbar^{2}}{2m_{e}}\nabla^{2}_{\mathbf{b}}+\frac{Z_{a}Z_{b}}{4\pi\epsilon_{0}|\mathbf{a}-\mathbf{b}|},$
(1)

where $\epsilon_{0}$ represents the vacuum permittivity. Further in this paper
atomic units are introduced. In this paper, our focus is on a pair of two
electrons with charges $(Z_{a}=Z_{b}=-1)$, and they interact through the
repulsive Coulomb interaction. For $(Z_{a}=1,Z_{b}=-1)$, the Hamiltonian
$\hat{H}_{p}$ describes the positronium system. Moreover, in the limit where
$m_{a}$ approaches infinity while maintaining $(Z_{a}=1,Z_{b}=-1)$, the
Hamiltonian reduces to that of the hydrogen atom.

In center of mass $\mathbf{R}=(\mathbf{a}+\mathbf{b})/2$ and relative
$\mathbf{r}=\mathbf{a}-\mathbf{b}$ coordinates the Hamiltonian $\hat{H}_{p}$
separates

$\hat{H}_{p}=-c_{K}\nabla^{2}_{\mathbf{R}}+\left(-c_{k}\nabla^{2}_{\mathbf{r}}+\frac{1}{r}\right),$
(2)

where in atomic units $c_{K}=1/4$ and $c_{k}=1$. Note that for the hydrogen
atom $c_{k}=1/2$. The energy of the system sum of two terms:
$E_{K}=c_{K}K^{2}$ for the center-of-mass motion and $E_{k}$ of the relative
motion. In absence of external potential there is $E_{k}=c_{k}k^{2}$.

The Hamiltonian $\hat{H}_{p}$ does not possess any bound states, and its wave
functions are delocalized. We have

$\Phi(\mathbf{R},\mathbf{r})=N_{k}e^{i\mathbf{K}\mathbf{R}}\psi_{klm}(\mathbf{r}).$
(3)

Here, $N_{k}$ represents the normalization factor, $\mathbf{K}$ is the wave
vector for the center-of-mass motion, and $\psi_{klm}(\mathbf{r})$ denotes the
continuous states governed by the repulsive Coulomb Hamiltonian

$\psi_{klm}(\mathbf{r})=R_{kl}(r)Y_{l}^{m}(\theta,\phi).$ (4)

In the above expression, $Y_{l}^{m}(\theta,\phi)$ stands for the spherical
harmonics in standard notation. The parameters are defined as follows:
$k=\sqrt{E_{k}/c_{k}}$ represents the wave vector for relative motion, $l$
denotes the azimuthal quantum number, and the radial function $R_{kl}(r)$ is
given by [19]

$\displaystyle R_{kl}(r)$ $\displaystyle=$
$\displaystyle\frac{C_{kl}(2kr)^{l}e^{ikr}}{(2l+1)!}\
_{1}F_{1}(i/k+l+1,2l+2,-2ikr),\ \ \ $ (5) $\displaystyle C_{kl}$
$\displaystyle=$ $\displaystyle 2ke^{-\pi/2k}|\Gamma(l+1+i/k)|$ (6)
$\displaystyle=$ $\displaystyle\sqrt{\frac{8\pi
k}{e^{2\pi/k}-1}}\prod_{s=1}^{l}\sqrt{s^{2}+\frac{1}{k^{2}}},$

When $l=0$, the product in Eq. (6) simplifies to unity. The function
${}_{1}F_{1}(\alpha,\gamma,z)$ is the confluent hypergeometric function

$_{1}F_{1}(\alpha,\gamma,z)=1+\frac{\alpha}{\gamma}z+\frac{\alpha(\alpha+1)}{\gamma(\gamma+1)}z^{2}+\ldots.$
(7)

For a fixed value of $l$, the functions $R_{kl}(r)$ are normalized according
to the criterion [19]

$\int_{0}^{\infty}R_{kl}(r)R_{k^{\prime}l}(r)r^{2}dr=2\pi\delta(k-k^{\prime}).$
(8)

In the $(\mathbf{R},\mathbf{r})$ coordinates, the retarded (advanced) two-
particle GF is defined as

$\displaystyle
g^{\pm}(\mathbf{R}_{1},\mathbf{R}_{2},\mathbf{r}_{1},\mathbf{r}_{2};E)=$
$\displaystyle\int\frac{d^{3}\mathbf{K}}{(2\pi)^{3}}\int_{0}^{\infty}\\!\\!\\!dk\sum_{l=0}^{\infty}\sum_{m=-l}^{l}\frac{e^{i\mathbf{K}(\mathbf{R}_{1}-\mathbf{R}_{2})}\psi_{klm}(\mathbf{r}_{1})\psi_{klm}^{*}(\mathbf{r}_{2})}{(E-c_{K}K^{2})-c_{k}k^{2}\pm
i\eta}\ \ $ (9) $\displaystyle=\frac{1}{(2\pi)^{3}}\int
e^{i\mathbf{K}(\mathbf{R}_{1}-\mathbf{R}_{2})}\tilde{g_{c}}^{\pm}(\mathbf{r}_{1},\mathbf{r}_{2},\epsilon_{K})d^{3}\mathbf{K}.\
\ $ (10)

Here, the term $\epsilon_{K}$ is defined as:

$\epsilon_{K}=E-c_{K}K^{2}.$ (11)

Additionally, we introduce a small positive parameter $\eta>0$. It is
important to note that the factor $c_{k}$ is already accounted for in the
Coulomb GF, as outlined in Eq. (1.3) in Ref. [5]. The function
$\tilde{g_{c}}(\mathbf{r}_{1},\mathbf{r}_{2};E)$ represents the one-particle
GF for the Coulomb potential, as discussed in [1, 5]

$\displaystyle\tilde{g_{c}}^{+}(\mathbf{r}_{1},\mathbf{r}_{2};E)=-\frac{\Gamma(1-i\nu)}{4\pi|\mathbf{r}_{1}-\mathbf{r}_{2}|}\times$
$\displaystyle\left(W^{1/2}_{i\nu}(U)\frac{\partial}{\partial V}{\cal
M}^{1/2}_{i\nu}(V)-{\cal M}^{1/2}_{i\nu}(V)\frac{\partial}{\partial
U}W^{1/2}_{i\nu}(U)\right),$ (12)

where

$\displaystyle U$ $\displaystyle=$ $\displaystyle-
ik\big{(}r_{1}+r_{2}+|\mathbf{r}_{1}-\mathbf{r}_{2}|\big{)},$ (13)
$\displaystyle V$ $\displaystyle=$ $\displaystyle-
ik\big{(}r_{1}+r_{2}-|\mathbf{r}_{1}-\mathbf{r}_{2}|\big{)},$ (14)

$k=\sqrt{E/c_{k}}$ with ${\rm Im}\big{\\{}k\big{\\}}>0$ and $\nu=-1/k$. The
functions ${\cal M}_{\kappa}^{\mu}(z)$ and $W_{\kappa}^{\mu}(z)$ are the
Whittaker functions in notation used by Buchholtz [20, 21]. In Eq. (II.1), the
function $\tilde{g_{c}}$ represents the Coulomb GF that describes both
attractive ($\nu>0$) and repulsive ($\nu<0$) potentials [11]. When $\nu=0$,
the function $\tilde{g_{c}}$ simplifies to the GF of a free particle.

For specific applications, an alternative representation of the Coulomb GF
proves to be more practical. In this representation, the Coulomb GF is
expanded using partial waves [14]

$\tilde{g_{c}}({\bm{r}}_{1},{\bm{r}}_{2};E)=\sum_{l=0}^{\infty}\sum_{m=-l}^{l}g_{l}(r_{1},r_{2};E)Y_{l}^{m}(\Omega_{1})Y_{l}^{m*}(\Omega_{2}),$
(15)

where $\Omega=(\theta,\phi)$ denotes angular variables, and the radial Coulomb
GFs are defined as [5, 22]

$\displaystyle g_{l}(r_{1},r_{2};E)$ $\displaystyle=$
$\displaystyle-\Gamma(1+l-i\nu)\frac{i\nu}{r_{1}r_{2}}\times$ (16)
$\displaystyle\times$ $\displaystyle{\cal
M}^{l+1/2}_{i\nu}(\text{-}2ikr_{\scriptscriptstyle<})W^{l+1/2}_{i\nu}(\text{-}2ikr_{\scriptscriptstyle>}),$

where $r_{\scriptscriptstyle<}=\min(r_{1},r_{2})$ and
$r_{\scriptscriptstyle>}=\max(r_{1},r_{2})$. It is important to note that in
the subsequent sections of this paper, any form of
$\tilde{g_{c}}({\bm{r}}_{1},{\bm{r}}_{2};E)$ is permissible.

The Coulomb GF presented in Eq. (II.1) is an analytic function of energy in
the complex energy plane. It exhibits a branch cut along the positive real
axis, which corresponds to the continuous spectrum. In the case of negative
energies, the Coulomb GF becomes a real function.

For an attractive Coulomb potential, the GF in Eq. (II.1) possesses simple
poles at energies corresponding to the singularities of $\Gamma(1-i\nu)$,
specifically for $1-i\nu=0,-1,\ldots$. This results in the discrete spectrum
associated with the hydrogen atom. However, for a repulsive potential, the
Coulomb GF has no poles, and consequently, no bound states exist. This same
principle applies to the GF of an electron pair as described in Eq. (10),
where again, no poles are present, and therefore, no bound states are formed.

###  II.2 Properties of two-particle GF

When $\mathbf{R}_{1}$ is distinct from $\mathbf{R}_{2}$ and $\mathbf{r}_{1}$
is different from $\mathbf{r}_{2}$, the integral over $d^{3}\mathbf{K}$ in Eq.
(10) exhibit no singularities. The outcome of this integration depends on the
signs of $E$ and $\epsilon_{K}$, leading to three distinct scenarios: i) When
both $E$ and $\epsilon_{K}$ are greater than zero, the Coulomb GF described in
Eq. (10) has oscillatory behavior. ii) In the case of $E>0$ and
$\epsilon_{K}<0$, or iii) when both $E$ and $\epsilon_{K}$ are negative, the
Coulomb GF in Eq. (10) experiences exponential decay with increasing distance
between the particles. Specifically, for $E>0$

$g^{\pm}(\mathbf{R}_{1},\mathbf{R}_{2},\mathbf{r}_{1},\mathbf{r}_{2};E>0)=I^{\pm}+I^{0},$
(17)

where

$I^{\pm}=\int_{0}^{\sqrt{E/c_{K}}}\frac{K\sin(KR_{12})}{2\pi^{2}R_{12}}\tilde{g_{c}}^{\pm}(\mathbf{r}_{1},\mathbf{r}_{2},\epsilon_{K})dK,$
(18)

and

$I^{0}=\int_{\sqrt{E/c_{K}}}^{\infty}\frac{K\sin(KR_{12})}{2\pi^{2}R_{12}}\tilde{g_{c}}^{0}(\mathbf{r}_{1},\mathbf{r}_{2},-|\epsilon_{K}|)dK,$
(19)

and $R_{12}=|\mathbf{R}_{1}-\mathbf{R}_{2}|$. In the case of $E<0$, the two-
electron GF is a real valued function expressed by a single term

$\displaystyle
g(\mathbf{R}_{1},\mathbf{R}_{2},\mathbf{r}_{1},\mathbf{r}_{2};E<0)=$
$\displaystyle\int_{0}^{\infty}\frac{K\sin(KR_{12})}{2\pi^{2}R_{12}}\tilde{g_{c}}^{0}(\mathbf{r}_{1},\mathbf{r}_{2},\varepsilon_{K})dK,$
(20)

where

$\varepsilon_{K}=-|E|-c_{K}K^{2}<0,$ (21)

see Eq. (11). In the equations above, $\tilde{g_{c}}^{\pm}$ represent the
retarded and advanced Coulomb GFs respectively, while $\tilde{g_{c}}^{0}$
corresponds to the Coulomb GF associated with negative energies. It’s
important to note that the integrals involving $\tilde{g_{c}}^{0}$ do not make
any contribution to the density of states.

Subsequently in this paper, we will adopt a simplified notation for the two-
particle GF

$g^{\pm}(\mathbf{a}_{1},\mathbf{b}_{1},\mathbf{a}_{2},\mathbf{b}_{2};E)\equiv
g(a_{1}b_{1}a_{2}b_{2}),$ (22)

and

$\displaystyle g(a_{1}b_{1}a_{2}b_{2})$ $\displaystyle=$
$\displaystyle\frac{1}{(2\pi)^{3}}\int
e^{i\mathbf{K}(\mathbf{a}_{1}+\mathbf{b}_{1})/2}e^{-i\mathbf{K}(\mathbf{a}_{2}+\mathbf{b}_{2})/2}\times$
(23)
$\displaystyle\tilde{g_{c}}(\mathbf{a}_{1}-\mathbf{b}_{1},\mathbf{a}_{2}-\mathbf{b}_{2},\epsilon_{K})d^{3}\mathbf{K}.$

In this notation we do not make distinction between the retarder and advanced
GFs.

###  II.3 Divergences of $g(a_{1}b_{1}a_{2}b_{2})$

In the majority of cases, the integral in Eq. (23) converges. However, for
certain GF’s arguments, it diverges, leading to singular behavior of
$g(a_{1}b_{1}a_{2}b_{2})$. This integral diverges either for
$\mathbf{a}_{1}+\mathbf{b}_{1}=\mathbf{a}_{2}+\mathbf{b}_{2}$ or for
$\mathbf{a}_{1}-\mathbf{b}_{1}=\mathbf{a}_{2}-\mathbf{b}_{2}$. This occurs
either for vanishing exponent in Eq. (23) or for equal arguments of Coulomb
GF. The cases having zero, one, or two nonzero values among the vectors
$\mathbf{a}_{1}$, $\mathbf{b}_{1}$, $\mathbf{a}_{2}$, and $\mathbf{b}_{2}$ are
detailed in Table 1. Instances with three nonzero arguments of GF are omitted.

Group | $(\mathbf{a}_{1}\mathbf{b}_{1}\mathbf{a}_{2}\mathbf{b}_{2})$ | $\mathbf{r}_{1}$ | $\mathbf{r}_{2}$ | $\mathbf{R}_{1}\text{-}\mathbf{R}_{2}$ | Integrand in Eq. (10)
---|---|---|---|---|---
0 | $(\mathbf{0}\mathbf{0}\mathbf{0}\mathbf{0})$ | $\mathbf{0}$ | $\mathbf{0}$ | $\mathbf{0}$ | $\tilde{g_{c}}(\mathbf{0},\mathbf{0};\epsilon_{K})d^{3}\mathbf{K}$
1 | $(\mathbf{x}\mathbf{0}\mathbf{x}\mathbf{0})$ | $\mathbf{x}$ | $\mathbf{x}$ | $\mathbf{0}$ | $\tilde{g_{c}}(\mathbf{x},\mathbf{x};\epsilon_{K})d^{3}\mathbf{K}$
1 | $(\mathbf{x}\mathbf{0}\mathbf{0}\text{-}\mathbf{x})$ | $\mathbf{x}$ | $\mathbf{x}$ | $2\mathbf{x}$ | $\tilde{g_{c}}(\mathbf{x},\mathbf{x};\epsilon_{K})e^{i\mathbf{K}\mathbf{x}}d^{3}\mathbf{K}$
1 | $(\mathbf{x}\mathbf{0}\mathbf{0}\mathbf{x})$ | $\mathbf{x}$ | $\text{-}\mathbf{x}$ | $\mathbf{0}$ | $\tilde{g_{c}}(\mathbf{x},\text{-}\mathbf{x};\epsilon_{K})d^{3}\mathbf{K}$
2 | $(\mathbf{x}\mathbf{y}\mathbf{x}\mathbf{y})$ | $\mathbf{r}$ | $\mathbf{r}$ | $\mathbf{0}$ | $\tilde{g_{c}}(\mathbf{r},\mathbf{r};\epsilon_{K})d^{3}\mathbf{K}$
2 | $(\mathbf{x}\mathbf{y}\mathbf{y}\text{-}\mathbf{x})$ | $\mathbf{r}$ | $\mathbf{r}$ | $2\mathbf{x}$ | $\tilde{g_{c}}(\mathbf{r},\mathbf{r};\epsilon_{K})e^{i\mathbf{K}\mathbf{x}}d^{3}\mathbf{K}$
2 | $(\mathbf{x}\mathbf{y}\mathbf{y}\mathbf{x})$ | $\mathbf{r}$ | $\text{-}\mathbf{r}$ | $\mathbf{0}$ | $\tilde{g_{c}}(\mathbf{r},\text{-}\mathbf{r};\epsilon_{K})d^{3}\mathbf{K}$
Table 1: Three groups of arguments of GF for which the integrals in Eqs. (10)
and (23) is diverge. Group $0$ (all arguments vanish), group $1$ (one non-zero
argument $\mathbf{x}$), and group $2$ (two non-zero vectors, $\mathbf{x}$ and
$\mathbf{y}$). We define $\mathbf{r}=\mathbf{x}-\mathbf{y}$. Permutations
causing simultaneous sign changes of $\mathbf{r}_{1}$ and $\mathbf{r}_{2}$, as
well as $\mathbf{R}_{1}$ and $\mathbf{R}_{2}$, are excluded.

For GF arguments provided in Table 1, it is necessary to use the spectral
representation of the GF. We illustrate this method to calculate
$g^{\pm}(0000;E)\equiv g^{\pm}_{0}(E)$. From Eq. (9) we have

$\displaystyle g_{0}^{+}(E)=$
$\displaystyle\int\frac{d^{3}\mathbf{K}}{(2\pi)^{3}}\int_{0}^{\infty}\\!\\!\sum_{l=0}^{\infty}\sum_{m=-l}^{l}\frac{|\psi_{klm}(\mathbf{0})|^{2}dk}{E-c_{K}K^{2}-c_{k}k^{2}+i\eta}.$
(24)

At $\mathbf{r}=\mathbf{0}$, all functions $\psi_{klm}(\mathbf{r})$ described
in Eq. (24) become zero, except for those where $l$ and $m$ are both equal to
zero. There is

$f(k)\equiv|\psi_{k00}(0)|^{2}=\frac{8\pi k}{(2\pi)(4\pi)[\exp(2\pi/k)-1]},$
(25)

where $k=\sqrt{E_{k}/c_{k}}$, see Eq. (2). The factor $(4\pi)$ in the
denominator of Eq. (25) arises from the normalization of the spherical
function $Y_{0}^{0}$, as described in Ref. [23]. Similarly, the factor of
$(2\pi)$ comes from the normalization of the radial functions in Eq. (8). Then
we have

$g^{\pm}_{0}(E)=\int\frac{d^{3}\mathbf{K}}{(2\pi)^{3}}\int_{0}^{\infty}\\!\\!\frac{f(k)dk}{E-c_{K}K^{2}-c_{k}k^{2}\pm
i\eta}.$ (26)

Using the the Dirac identity: $1/(x\pm i\eta)={\cal P}(1/x)\mp i\pi\delta(x)$
to Eq. (26) and integrating over angular variables we have

$\displaystyle{\rm Im}\big{\\{}g_{0}^{+}(E)\big{\\}}$ $\displaystyle=$
$\displaystyle\frac{(-\pi)}{2\pi^{2}}\int_{0}^{\infty}\\!\\!\int_{0}^{\infty}\\!\\!f(k)K^{2}\times$
(27) $\displaystyle\times$ $\displaystyle\delta(E-c_{K}K^{2}-c_{k}k^{2})dKdk.$

Next we introduce the polar coordinates
$K=\frac{t}{\sqrt{c_{K}}}\cos(\alpha)$,
$k=\frac{t}{\sqrt{c_{k}}}\sin(\alpha)$, with $0\leq\alpha\leq\pi/2$ and the
volume element: $dKdk=tdtd\alpha/\sqrt{c_{K}c_{k}}$. Using the identity

$\delta(t^{2}-a^{2})=\frac{1}{|2a|}\big{[}\delta(t-a)+\delta(t+a)\big{]},$
(28)

we obtain

$\displaystyle{\rm
Im}\big{\\{}g_{0}^{+}(E>0)\big{\\}}=\frac{(-\pi)}{2\pi^{2}\sqrt{c_{K}^{3}c_{k}}}\times$
$\displaystyle\int_{0}^{\pi/2}\\!\\!\int_{0}^{\infty}\\!\\!\\!f[(t/\sqrt{c_{k}})\sin(\alpha)]t^{3}\cos(\alpha)^{2}\delta(E-t^{2})dtd\alpha$
$\displaystyle=\frac{(-\pi)E\Theta(E)}{4\pi^{2}\sqrt{c_{K}^{3}c_{k}}}\int_{0}^{\pi/2}f[\sqrt{E/c_{k}}\sin(\alpha)]\cos(\alpha)^{2}d\alpha,\
\ \ \ \ $ (29)

where $\Theta(E)$ is the step function. For negative energies in the second
line of Eq. (29) there is $\int_{0}^{\infty}\delta(-|E|-t^{2})dt=0$, and ${\rm
Im}\big{\\{}g_{0}^{+}(E<0)\big{\\}}$ vanishes.

The real part of $g_{0}^{+}(E)$ is Hilbert transform of ${\rm
Im}\big{\\{}g_{0}^{+}(E)\big{\\}}$, and it diverges for all energies. To
circumvent this divergence, a cutoff energy $W$ is introduced, beyond which
the density of states: $\rho_{0}(E)=(-1/\pi){\rm
Im}\big{\\{}g_{0}^{+}(E)\big{\\}}$ vanishes. In physical terms, this signifies
the finite width $W$ of the energy band. Consequently, we have

${\rm Re}\big{\\{}g_{0}^{+}(E)\big{\\}}=-\frac{1}{\pi}{\cal
P}\int_{0}^{W}\frac{{\rm
Im}\big{\\{}g_{0}^{+}(E^{\prime})\big{\\}}}{E-E^{\prime}}.$ (30)

The above integral exists for finite $W$ and it diverges for $W\to\infty$.

Figure 1: Panel (a), solid line: $\rho_{0}(E)$ for a pair of electrons in
presence of Coulomb potential, as given in Eq. (29). Panel (a), dotted line:
$\rho_{f0}(E)$ for a pair of electrons in the absence of Coulomb interaction,
as given in Eq. (80). Panel (b), solid line: $\rho_{c0}(E)$ for electron in a
repulsive Coulomb potential, as given in Eq. (78). Panel (b), dotted line:
$\rho_{e0}(E)$ for free electron, see Eq. (77). In the presence of the Coulomb
potential there is a non-zero overlap between pairs of electrons and a non-
negligible probability of finding an electron at the center of the Coulomb
potential. For two-electrons LDOS is in $r_{B}^{-6}$ units, while for one-
electron LDOS is in $r_{B}^{-3}$ units.

In Figure 1a, we present the local density of states denoted as
$\rho_{0}(E)=(-1/\pi){\rm Im}\big{\\{}g_{0}^{+}(E)\big{\\}}$, as per Eq. (29),
represented by the solid line. The results are displayed on a log-log scale.
The primary observation from Figure 1a is that the LDOS for an electron pair
subject to Coulomb forces does not vanish at $\mathbf{R}=\mathbf{r}=0$. In
other words, there is a non-zero overlap between the two electrons, a
manifestation of a purely quantum effect. This behavior is intriguing, as
classically, two electrons would repel each other and remain far apart. This
phenomenon, albeit somewhat mysterious, is also observable for a single
electron in a repulsive Coulomb potential.

In Figure 1b, we plot the local density of states $\rho_{c0}(E)$ for a single
electron in a repulsive potential, as given in Eq. (78). As Figure 1b
illustrates, $\rho_{c0}(E)$ also does not vanish for any energy. This implies
that there is a non-zero probability that the electron overlaps with the
center of the repulsive Coulomb potential, a phenomenon rooted solely in
quantum effects without a classical analogue.

The LDOS for an electron pair in Figure 1a and that for a single electron in a
repulsive potential in Figure 1b exhibit similar qualitative behavior. For low
energies, both LDOS profiles are negligibly small and gradually increase with
energy. At significantly high energy levels, a notable asymptotic behavior
becomes evident. The LDOS for the electron pair gradually converges towards
the limiting local density of states $\rho_{f0}(E)$, which represents an
electron pair in the absence of a Coulomb potential and is defined in Eq.
(80). This convergence is clearly illustrated in Figure 1a, where the dashed
line closely approaches the LDOS for the electron pair as energy increases.

Similarly, in this same high-energy limit, the density of states
$\rho_{e0}(E)$ for a single electron subjected to a repulsive potential
gradually approaches the density of states of a free electron $\rho_{e0}(E)$,
as governed by Eq. (77). This convergence can be observed in Figure 1b, where
the dashed line closely aligns with the LDOS for a single electron in the
presence of a repulsive potential as energy levels rise.

##  III Green’s Function for Two Electrons with Spin

In the preceding section, we conducted an analysis of the GF for a pair of
electrons interacting via Coulomb forces, while disregarding the influence of
electron spins and the constraints imposed by the Pauli exclusion principle.
In this section, we extend our study to encompass the GF of an electron pair,
accounting for their spin properties and the inherent limitations imposed by
the Pauli principle.

Let’s consider an arbitrary two-electron Hamiltonian, denoted as $\hat{H}$,
which is independent of electron spins, and let $|\rm N\rangle$ represent one
of its eigenstates. We can then decompose $|\rm N\rangle$ into two components

$|{\rm N}\rangle=|{\rm n}\rangle|\chi\rangle,$ (31)

Here, ${\rm n}$ encompasses all the quantum numbers that characterize state
$|{\rm n}\rangle$, while $|\chi\rangle$ signifies the state associated with
the electron spins.

The wave function $\langle\mathbf{a}\mathbf{b}\sigma_{a}\sigma_{b}|{\rm
N}\rangle$ separates on the spin-independent part $\Psi_{\rm
n}(\mathbf{a},\mathbf{b})$ and spin dependent function
$\chi(\sigma_{a},\sigma_{b})$, where
$\sigma_{a},\sigma_{b}\in\\{\uparrow,\downarrow\\}$ are electrons spins.
Because of the Pauli principle there is

$\Psi_{\rm n}(\mathbf{a},\mathbf{b})\chi(\sigma_{a},\sigma_{b})=-\Psi_{\rm
n}(\mathbf{b},\mathbf{a})\chi(\sigma_{b},\sigma_{a}).$ (32)

The aforementioned condition can be satisfied under two distinct scenarios.
Firstly, when the state function $\chi(\sigma_{a},\sigma_{b})$ assumes the
singlet state form
$\frac{1}{\sqrt{2}}(|\uparrow\downarrow\rangle-|\downarrow\uparrow\rangle)$,
and the wave function $\Psi_{\rm n}(\mathbf{a},\mathbf{b})$ demonstrates even
symmetry with respect to the exchange of variables $\mathbf{a}$ and
$\mathbf{b}$. Alternatively, the condition is met when
$\chi(\sigma_{a},\sigma_{b})$ represents one of the triplet states:
$\frac{1}{\sqrt{2}}(|\uparrow\downarrow\rangle+|\downarrow\uparrow\rangle)$,
$|\uparrow\uparrow\rangle$, or $|\downarrow\downarrow\rangle$, and the wave
function $\Psi_{\rm n}(\mathbf{a},\mathbf{b})$ exhibits odd symmetry with
respect to $\mathbf{a}$ and $\mathbf{b}$.

Let $|\chi_{s}\rangle$ and $|\chi_{t}\rangle$ be singlet and triplet states,
respectively. Then the Green’s function $\hat{G}=(E-\hat{H})^{-1}$ is sum of
two terms

$\hat{G}^{\pm}=\sum_{\rm n}|\chi_{s}\rangle\langle\chi_{s}|\sum_{{\rm
n}_{e}}\frac{|{\rm n}_{e}\rangle\langle{\rm n}_{e}|}{E-H\pm
i\eta}+\sum_{\chi_{t}}|\chi_{t}\rangle\langle\chi_{t}|\sum_{{\rm
n}_{o}}\frac{|{\rm n}_{o}\rangle\langle{\rm n}_{o}|}{E-H\pm
i\eta}=\hat{G}^{\pm}_{e}+\hat{G}^{\pm}_{o},$ (33)

and the summation over $\chi_{t}$ extends across all three triplet states. In
the position representation we have

$\displaystyle G^{\pm}_{e}(a_{1}b_{1}a_{2}b_{2})$ $\displaystyle=$
$\displaystyle\Lambda_{s}\sum_{{\rm n}_{e}}\frac{\Psi_{{\rm
n}_{e}}(\mathbf{a}_{1},\mathbf{b}_{1})\Psi^{*}_{{\rm
n}_{e}}(\mathbf{a}_{2},\mathbf{b}_{2})}{E-E_{{\rm n}_{e}}\pm i\eta},$ (34)
$\displaystyle G^{\pm}_{o}(a_{1}b_{1}a_{2}b_{2})$ $\displaystyle=$
$\displaystyle\Lambda_{t}\sum_{{\rm n}_{o}}\frac{\Psi_{{\rm
n}_{o}}(\mathbf{a}_{1},\mathbf{b}_{1})\Psi^{*}_{{\rm
n}_{o}}(\mathbf{a}_{2},\mathbf{b}_{2})}{E-E_{{\rm n}_{o}}\pm i\eta},$ (35)

where $\Lambda_{s}=\chi_{s}^{\dagger}\chi_{s}$,
$\Lambda_{t}=\sum_{\chi_{t}}\chi_{t}^{\dagger}\chi_{t}$. By exploiting the
symmetry properties of the functions $\Psi_{{\rm
n}_{e}}(\mathbf{a}_{1},\mathbf{b}_{1})$ and $\Psi_{{\rm
n}_{o}}(\mathbf{a}_{1},\mathbf{b}_{1})$ with respect to the exchange of
coordinates $\mathbf{a}\leftrightarrow\mathbf{b}$ we find

$\displaystyle G_{e}(a_{1}b_{1}a_{2}b_{2})=$ $\displaystyle
G_{e}(b_{1}a_{1}a_{2}b_{2})=$ $\displaystyle G_{e}(a_{1}b_{1}b_{2}a_{2}),$
(36) $\displaystyle G_{o}(a_{1}b_{1}a_{2}b_{2})=$
$\displaystyle\text{-}G_{o}(b_{1}a_{1}a_{2}b_{2})=$
$\displaystyle\text{-}G_{o}(a_{1}b_{1}b_{2}a_{2}).$ (37)

As a consequence of Eq. (37) there is

$G_{o}(a_{1}b_{1}00)=G_{o}(00a_{2}b_{2})=G_{o}(0000)=0,$ (38)

i.e., the the odd part of GF vanishes for $\mathbf{R}=\mathbf{r}=\mathbf{0}$.
In Eqs. (34) and (35), the summation involves distinct functions, namely,
$\Psi_{{\rm n}_{e}}(\mathbf{a},\mathbf{b})$ for even states and $\Psi_{{\rm
n}_{o}}(\mathbf{a},\mathbf{b})$ for odd states. However, it is important to
note that the denominators in both cases are identical and equal to
$(E-c_{K}K^{2}-c_{k}k^{2}\pm i\eta)$. The functions $\Psi_{{\rm
n}_{e}}(\mathbf{a},\mathbf{b})$ and $\Psi_{{\rm
n}_{o}}(\mathbf{a},\mathbf{b})$ can be readily derived in the
$(\mathbf{R},\mathbf{r})$ coordinates, as outlined in Eq. (3)

$\displaystyle\Psi_{e}(\mathbf{R},\mathbf{r})$ $\displaystyle=$
$\displaystyle\frac{1}{2}e^{i\mathbf{K}\mathbf{R}}\left[\psi(\mathbf{r})+\psi(-\mathbf{r})\right],$
(39) $\displaystyle\Psi_{o}(\mathbf{R},\mathbf{r})$ $\displaystyle=$
$\displaystyle\frac{1}{2}e^{i\mathbf{K}\mathbf{R}}\left[\psi(\mathbf{r})-\psi(-\mathbf{r})\right].$
(40)

where $\psi(\mathbf{r})\equiv N_{k}\psi(\mathbf{r})_{klm}$ in Eq. (3).
Inserting Eqs. (39) and (40) into Eqs. (34) and (35) we find [see Eq. (10)]

$\displaystyle
g^{\pm}_{e/o}(\mathbf{R}_{1}\mathbf{r}_{1}\mathbf{R}_{2}\mathbf{r}_{2})=\frac{1}{4}\frac{1}{(2\pi)^{3}}\left\\{\begin{array}[]{c}\Lambda_{s}\\\
\Lambda_{t}\end{array}\right\\}\int
d^{3}\mathbf{K}\sum_{l=0}^{\infty}\sum_{m=-l}^{l}\int_{0}^{\infty}dk\frac{e^{\mathbf{K}(\mathbf{R}_{1}-\mathbf{R}_{2})}[\psi(\mathbf{r}_{1})\pm\psi(\text{-}\mathbf{r}_{1})][\psi^{*}(\mathbf{r}_{2})\pm\psi^{*}(\text{-}\mathbf{r}_{2})]}{E-c_{K}K^{2}-c_{k}k^{2}\pm
i\eta}$ (43)
$\displaystyle=\frac{1}{4}\frac{1}{(2\pi)^{3}}\left\\{\begin{array}[]{c}\Lambda_{s}\\\
\Lambda_{t}\end{array}\right\\}\int
e^{\mathbf{K}(\mathbf{R}_{1}-\mathbf{R}_{2})}\left[\tilde{g_{c}}^{\pm}(\mathbf{r}_{1},\mathbf{r}_{2},\epsilon_{K})\pm\tilde{g_{c}}^{\pm}(\mathbf{r}_{1},\text{-}\mathbf{r}_{2},\epsilon_{K})\pm\tilde{g_{c}}^{\pm}(\text{-}\mathbf{r}_{1},\mathbf{r}_{2},\epsilon_{K})+\tilde{g_{c}}^{\pm}(\text{-}\mathbf{r}_{1},\text{-}\mathbf{r}_{2},\epsilon_{K})\right]d^{3}\mathbf{K}$
(46) $\displaystyle=\frac{1}{2}\left\\{\begin{array}[]{c}\Lambda_{s}\\\
\Lambda_{t}\end{array}\right\\}\Big{[}g^{\pm}(\mathbf{R}_{1}\mathbf{r}_{1}\mathbf{R}_{2}\mathbf{r}_{2},\epsilon_{K})\pm
g^{\pm}(\mathbf{R}_{1}\mathbf{r}_{1}\mathbf{R}_{2}\text{-}\mathbf{r}_{2},\epsilon_{K})\Big{]}.$
(49)

Here, $\tilde{g_{c}}^{\pm}(\mathbf{r}_{1},\mathbf{r}_{2})$ represents the
Coulomb GF, as defined in Eq. (II.1), and we have utilized the property
$\tilde{g_{c}}^{\pm}(\text{-}\mathbf{r}_{1},\text{-}\mathbf{r}_{2})=\tilde{g_{c}}^{\pm}(\mathbf{r}_{1},\mathbf{r}_{2})$.
In this equation, the even function comprises the term $\Lambda_{s}$, whereas
the odd function incorporates $\Lambda_{o}$.

When the Coulomb GF is expanded in partial waves, see Eq. (15), the odd and
even GFs in Eq. (43) can be expressed in terms of spherical harmonic having
even and odd angular quantum numbers only

$\displaystyle
g^{\pm}_{e}(\mathbf{R}_{1}\mathbf{r}_{1}\mathbf{R}_{2}\mathbf{r}_{2})$
$\displaystyle=$
$\displaystyle\frac{\Lambda_{s}}{(2\pi)^{3}}\sum_{l=0}^{\infty}\sum_{m=-l}^{l}Y_{2l}^{m}(\Omega_{1})Y_{2l}^{m*}(\Omega_{2})\int
e^{\mathbf{K}(\mathbf{R}_{1}-\mathbf{R}_{2})}g^{\pm}_{2l}(r_{1},r_{2};\epsilon_{K})d^{3}\mathbf{K},$
(50) $\displaystyle
g^{\pm}_{o}(\mathbf{R}_{1}\mathbf{r}_{1}\mathbf{R}_{2}\mathbf{r}_{2})$
$\displaystyle=$
$\displaystyle\frac{\Lambda_{t}}{(2\pi)^{3}}\sum_{l=0}^{\infty}\sum_{m=-l}^{l}Y_{2l+1}^{m}(\Omega_{1})Y_{2l+1}^{m*}(\Omega_{2})\int
e^{\mathbf{K}(\mathbf{R}_{1}-\mathbf{R}_{2})}g^{\pm}_{2l+1}(r_{1},r_{2};\epsilon_{K})d^{3}\mathbf{K}.$
(51)

Returning to the $(\mathbf{a},\mathbf{b})$ coordinates, we obtain from Eq.
(43) the following expressions

$\displaystyle
g^{\pm}_{e}(a_{1}b_{1}a_{2}b_{2};E)=\frac{\Lambda_{ss}}{2(2\pi)^{3}}\int
e^{i\mathbf{K}(\mathbf{a}_{1}+\mathbf{b}_{1}-\mathbf{a}_{2}-\mathbf{b}_{2})/2}\big{[}\tilde{g_{c}}^{\pm}(\mathbf{a}_{1}-\mathbf{b}_{1},\mathbf{a}_{2}-\mathbf{b}_{2},\epsilon_{K})+\tilde{g_{c}}^{\pm}(\mathbf{a}_{1}-\mathbf{b}_{1},\mathbf{b}_{2}-\mathbf{a}_{2},\epsilon_{K})\big{]}d^{3}\mathbf{K},$
(52) $\displaystyle
g^{\pm}_{o}(a_{1}b_{1}a_{2}b_{2};E)=\frac{\Lambda_{tt}}{2(2\pi)^{3}}\int
e^{i\mathbf{K}(\mathbf{a}_{1}+\mathbf{b}_{1}-\mathbf{a}_{2}-\mathbf{b}_{2})/2}\big{[}\tilde{g_{c}}^{\pm}(\mathbf{a}_{1}-\mathbf{b}_{1},\mathbf{a}_{2}-\mathbf{b}_{2},\epsilon_{K})-\tilde{g_{c}}^{\pm}(\mathbf{a}_{1}-\mathbf{b}_{1},\mathbf{b}_{2}-\mathbf{a}_{2},\epsilon_{K})\big{]}d^{3}\mathbf{K}.$
(53)

Note that when dealing with non-vanishing exponents and distinct arguments for
the Coulomb GF, the integration over the angular variables of the vector
$\mathbf{K}$ is straightforward, as described in Eqs. (18)–(II.2). However, in
cases where the exponent’s argument becomes zero or the arguments of the
Coulomb GF are equal, as summarized in Table 1, the integrals in Eqs.
(43)–(53) diverge. In such instances, a different approach is required, which
is elaborated upon in Sections II and V.

##  IV Dyson equations for even and odd GFs

Let us examine the electron pair within the context of an external potential
$V(\mathbf{r})$ that does not depend on spin and a two-electron interaction
represented as $u(\mathbf{a},\mathbf{b})$. The system’s Hamiltonian is given
by

$\hat{H}=\hat{H}_{p}+V(\mathbf{a})+V(\mathbf{b})+u(\mathbf{a},\mathbf{b})\equiv\hat{H}_{p}+U(\mathbf{a},\mathbf{b}).$
(54)

Let $\hat{G}=\hat{G}_{e}+\hat{G}_{o}$ be the GF of the Hamiltonian in Eq. (54)
and $\hat{g}=\hat{g}_{e}+\hat{g}_{o}$ be the GF of the electron pair, see Eqs.
(52) and (53). Since $\hat{G}$ also separate into even and odd parts, see Eq.
(33), we have $\hat{G}_{e}=\Lambda_{s}\widetilde{G}_{e}$,
$\hat{G}_{t}=\Lambda_{t}\widetilde{G}_{o}$, and
$\hat{g}_{s}=\Lambda_{s}\widetilde{g}_{e}$, and
$\hat{g}_{t}=\Lambda_{t}\widetilde{g}_{o}$. Then the Dyson equation
$\hat{G}=\hat{g}+\hat{g}\hat{U}\hat{G}$ for the total GF reads

$\displaystyle\big{(}\Lambda_{s}\widetilde{G}_{e}+\Lambda_{t}\widetilde{G}_{o}\big{)}=$
$\displaystyle\big{(}\Lambda_{s}\widetilde{g}_{e}+\Lambda_{t}\widetilde{g}_{o}\big{)}+\big{(}\Lambda_{s}\widetilde{g}_{e}+\Lambda_{t}\widetilde{g}_{o}\big{)}\hat{U}\big{(}\Lambda_{s}\widetilde{G}_{e}+\Lambda_{t}\widetilde{G}_{o}\big{)}.\
\ \ $ (55)

In the presence of spin-independent potentials, it holds that
$\Lambda_{s}\hat{U}\Lambda_{t}=0$ due to the orthogonality between singlet and
triplet states. Consequently, Eq. (55) can be split into two separate Dyson
equations, one for the odd GFs and another for the even GFs

$\displaystyle\hat{G}_{e}$ $\displaystyle=$
$\displaystyle\hat{g}_{e}+\hat{g}_{e}\hat{U}\hat{G}_{e},$ (56)
$\displaystyle\hat{G}_{o}$ $\displaystyle=$
$\displaystyle\hat{g}_{o}+\hat{g}_{o}\hat{U}\hat{G}_{o}.$ (57)

By taking the matrix element of both sides of equations (56) and (57) between
the bra state $\langle\mathbf{a}_{1}\mathbf{b}_{1}|$ and the ket state
$|\mathbf{a}_{2}\mathbf{b}_{2}\rangle$, we obtain

$\displaystyle G_{e/o}(a_{1}b_{1}a_{2}b_{2})=g_{e/o}(a_{1}b_{1}a_{2}b_{2})$
$\displaystyle+\iint
g_{e/o}(a_{1}b_{1}a_{3}b_{3})U(a_{3}b_{3})G_{e/o}(a_{3}b_{3}a_{2}b_{2})d^{3}a_{3}d^{3}b_{3}.\
\ \ \ \ \ $ (58)

When the operator $\hat{U}$ depends on electron spins, such as when it
incorporates spin-orbit interactions, it becomes necessary to apply the
general formula as presented in Eq. (55).

##  V LDOS for electrons pair

For a system consisting of two electrons, the local density of states can be
derived from the GF as follows

$\varrho(\mathbf{a}_{1},\mathbf{b}_{1};E)=(-1/\pi)\mathrm{Im}\
\mathrm{Tr}\\{g(a_{1}b_{1}a_{1}b_{1};E)\\}.$ (59)

By taking limits $\mathbf{a}_{2}\to\mathbf{a}_{1}$,
$\mathbf{b}_{2}\to\mathbf{b}_{1}$ in Eqs. (52) and (53), setting
$\mathbf{r}=\mathbf{a}_{1}-\mathbf{a}_{2}$ and using Eq. (33) we have

$\varrho(\mathbf{r};E)={\cal S}_{s}\varrho_{e}(\mathbf{r};E)+{\cal
S}_{t}\varrho_{o}(\mathbf{r};E),$ (60)

where

$\displaystyle\varrho_{e}(\mathbf{r};E)$ $\displaystyle=$
$\displaystyle-\frac{1}{2}{\rm
Im}\Big{\\{}\frac{1}{2\pi^{3}}\int_{0}^{K_{m}}\tilde{g_{c}}^{+}(\mathbf{r},\mathbf{r};\epsilon_{K})K^{2}dK+\frac{1}{2\pi^{3}}\int_{0}^{K_{m}}\tilde{g_{c}}^{+}(\mathbf{r},-\mathbf{r};\epsilon_{K})K^{2}dK\Big{\\}}\equiv\frac{1}{2}\big{[}\varrho_{+}(\mathbf{r};E)+\varrho_{-}(\mathbf{r};E)\big{]},$
(61) $\displaystyle\varrho_{o}(\mathbf{r};E)$ $\displaystyle=$
$\displaystyle-\frac{1}{2}{\rm
Im}\Big{\\{}\frac{1}{2\pi^{3}}\int_{0}^{K_{m}}\tilde{g_{c}}^{+}(\mathbf{r},\mathbf{r};\epsilon_{K})K^{2}dK-\frac{1}{2\pi^{3}}\int_{0}^{K_{m}}\tilde{g_{c}}^{+}(\mathbf{r},-\mathbf{r};\epsilon_{K})K^{2}dK\Big{\\}}\equiv\frac{1}{2}\big{[}\varrho_{+}(\mathbf{r};E)-\varrho_{-}(\mathbf{r};E)\big{]},$
(62)

and ${\cal S}_{s}=\mathrm{Tr}{\Lambda_{ss}}=1$, ${\cal
S}_{t}=\mathrm{Tr}{\Lambda_{tt}}=3$, as detailed in Appendix. In the equations
above, we have employed the following definitions: $K_{m}=\sqrt{E/c_{K}}$ and
$\epsilon_{K}=E-c_{K}K^{2}$, as indicated in Eq. (11). The limits of
integration in Eqs. (61) and (62) arise from the condition that, for
$\epsilon_{K}<0$ the imaginary part of
$\tilde{g_{c}}^{+}(\mathbf{r}_{1},\mathbf{r}_{2};\epsilon_{K})$ vanishes, as
discussed in Section II and Ref. [5]. Finally we obtain

$\varrho(\mathbf{r};E)=2\varrho_{+}(\mathbf{r};E)-\varrho_{-}(\mathbf{r};E),$
(63)

For negative energies, both $\varrho_{+}(\mathbf{r};E)$ and
$\varrho_{-}(\mathbf{r};E)$ become null, as in such cases, $\epsilon_{K}\leq
0$, and the imaginary component of the Coulomb GF vanishes for all values of
$K$. However, for positive energies and finite values of $r$, the situation is
different. The quantity $\varrho_{-}(\mathbf{r};E)$ can be determined through
direct numerical integration [24] of the Coulomb GF for
$\mathbf{r}_{2}=-\mathbf{r}_{1}$

$\varrho_{-}(\mathbf{r};E)=-\frac{1}{2\pi^{3}}\int_{0}^{K_{m}}\left[\frac{\Gamma(1-i\nu^{\prime})}{8\pi
r}W^{1/2}_{i\nu}(-4ik^{\prime}r)\right]K^{2}dK,$ (64)

as detailed in Eqs. (II.1) and (13). Here,
$k^{\prime}=\sqrt{\epsilon_{K}/c_{k}}$, $\nu_{K}=-1/k^{\prime}$, see Eqs. (11)
and (II.1). In the above equation the Whittaker function is well-defined for
all $r>0$ and $k>0$.

To obtain $\tilde{g_{c}}^{+}(\mathbf{r},\mathbf{r};E>0)$ we compute
$\tilde{g_{c}}^{+}(\mathbf{r},\mathbf{r}_{2};E>0)$ in the limit as
$\mathbf{r}_{2}$ approaches $\mathbf{r}$. Specifically, we express
$\mathbf{r}_{2}$ as $\mathbf{r}+\bm{\delta}$, where $\bm{\delta}$ is a small
vector oriented in an arbitrary direction. It is important to note that we
assume both $\mathbf{r}$ and $\bm{\delta}$ to be non-vanishing, with the
condition that $\delta\ll r$. Then we have in Eq. (13)

$\displaystyle U$ $\displaystyle\simeq$ $\displaystyle-
ik\big{(}r+r+\delta\cos(\varphi)+\delta\big{)}\equiv Z+A+\Delta,$ (65)
$\displaystyle V$ $\displaystyle\simeq$ $\displaystyle-
ik\big{(}r+r+\delta\cos(\varphi)-\delta\big{)}\equiv Z+A-\Delta,$ (66)

where $Z=-2ikr$, $\varphi$ is the angle between $\mathbf{r}$ and
$\bm{\delta}$, $A=-ik\delta\cos(\varphi)$, and $\Delta=-ik\delta$. On applying
the Taylor expansion to $W^{1/2}_{i\nu}(U)\equiv{\rm W}(U)$ and ${\cal
M}^{1/2}_{i\nu}(V)\equiv{\rm M}(V)$ ($n=0$) and their derivatives ($n=1$) we
have

$\frac{d^{n}}{dZ^{n}}{\rm W}(Z+A+\Delta)\simeq\frac{d^{n}}{dZ^{n}}{\rm
W}(Z)+(A+\Delta)\frac{d^{n+1}}{dZ^{n+1}}{\rm W}(Z),$ (67)
$\frac{d^{n}}{dZ^{n}}{\rm M}(Z+A-\Delta)\simeq\frac{d^{n}}{dZ^{n}}{\rm
M}(Z)+(A-\Delta)\frac{d^{n+1}}{dZ^{n+1}}{\rm M}(Z).$ (68)

When substituting Eqs. (67) and (68) into Eq. (II.1) while retaining terms at
the lowest order in both $A$ and $\Delta$, we arrive at

$\displaystyle\tilde{g_{c}}^{+}(\mathbf{r},\mathbf{r}+\bm{\delta};E)\simeq-\frac{\Gamma(1-i\nu)}{4\pi\delta}\Big{\\{}\Big{[}{\rm
W}(Z)\frac{d{\rm M}(Z)}{dZ}-{\rm M}(Z)\frac{d{\rm
W}(Z)}{dZ}\Big{]}+A\Big{[}{\rm W}(Z)\frac{d^{2}{\rm M}(Z)}{dZ^{2}}-{\rm
M}(Z)\frac{d^{2}{\rm{\rm W}}(Z)}{dZ^{2}}\Big{]}+$
$\displaystyle+\Delta\Big{[}{\rm M}(Z)\frac{d^{2}{\rm
W}(Z)}{dZ^{2}}-2\frac{d{\rm M}(Z)}{dZ}\frac{d{\rm W}(Z)}{dZ}+{\rm
W}(Z)\frac{d^{2}{\rm M}(Z)}{dZ^{2}}\Big{]}\Big{\\}}+\ldots.$ (69)

The first bracket in Eq. (V) represents the Wronskian of two Whittaker
functions [20, 21]

${\cal W}\left\\{W^{1/2}_{i\nu}(Z),{\cal
M}^{1/2}_{i\nu}(Z)\right\\}=\frac{\Gamma(1)}{\Gamma(1-i\nu)}.$ (70)

The term linear in $A$ vanishes since its is direct proportional to the
derivative of the Wronskian in Eq. (70) in respect of $Z$. In the last bracket
in Eq. (V) the second order derivatives $d^{2}{\rm W}(Z)/dZ^{2}$ and
$d^{2}{\rm M}(Z)/dZ^{2}$ are eliminated by employing the Whittaker equation

$\frac{d^{2}}{dZ^{2}}\left\\{\begin{array}[]{c}{\rm W}(Z)\\\ {\rm
M}(Z)\end{array}\right\\}+\left(-\frac{1}{4}+\frac{\kappa}{Z}\right)\left\\{\begin{array}[]{c}{\rm
W}(Z)\\\ {\rm M}(Z)\end{array}\right\\}=0,$ (71)

where $\kappa=i\nu$. Since $\Delta/\delta=-ik$ we obtain from Eq. (V)

$\tilde{g_{c}}^{+}(\mathbf{r},\mathbf{r}+\bm{\delta};E)\simeq-\frac{1}{4\pi\delta}+\frac{ik\Gamma(1-i\nu)}{4\pi}\Big{[}2\frac{d{\rm
M}(Z)}{dZ}\frac{d{\rm W}(Z)}{dZ}+2{\rm M}(Z){\rm
W}(Z)\Big{(}\frac{\kappa}{Z}-\frac{1}{4}\Big{)}\Big{]}+\ldots.$ (72)

In the limit as $\delta\to 0$, the divergence of
$\tilde{g_{c}}^{+}(\mathbf{r},\mathbf{r}+\bm{\delta};E)$ arises primarily from
the first term in Eq. (72), which is a real-valued quantity. For positive
energies, the second term in Eq. (72) is a complex number with a non-vanishing
imaginary component, thereby resulting in a finite LDOS in Eqs. (61) and (62).
The integration over this term is accomplished through standard numerical
techniques [24].

Figure 2: Local densities of states $\varrho(\mathbf{r};E)$, defined in Eq.
(63) as a functions of inter-electron distance $r$, calculated for three
distinct energy values. In Panel (c), the dotted line represents
$\rho_{f}(E)$, as given in Eq. (73), for $E=4$ Ha. LDOS is in $r_{B}^{-6}$
units.

In Figure 2, we have plotted the local densities of states, as defined in Eq.
(63) for three different energy levels $E$. In Figure 2a the energy is
significantly lower in comparison to the Hartree energy. The LDOS remains
generally low within the considered range of electron-electron distances $r$
and virtually approaches zero in the region of $0<r<3$rB, where rB represents
the Bohr radius. Beyond this region, the LDOS exhibits an almost parabolic
growth concerning distance, following the relationship:
$\varrho(E)\propto(r-3r_{B})^{2}$. Physically, this signifies that at lower
energies, the electrons are considerably distant from each other.

For an energy of $E=1$ Ha, as seen in Figure 2b, the LDOS displays a
transition-like behavior. It nearly diminishes for small electron-electron
distances and then gradually increases for larger values of $r$. At
approximately $r=4.5$rB, the LDOS curve shows an inflection point, suggesting
the onset of LDOS saturation. As energy increases, as shown in Figure 2c, we
observe a logistic-like dependence of the LDOS on $r$. It nearly vanishes as
$r$ approaches zero and then rapidly increases, seemingly reaching saturation
for $r>6$rB. The behavior of the LDOS for large values of $r$ indicates that
saturation occurs for LDOS values corresponding to those of a free electron
pair. According to Eq. (80) we have

$\rho_{f}(E)\approx 2\times\frac{E^{2}}{16\pi^{3}}.$ (73)

The factor of $2$ arises from the summation over two spin directions, as we
assume a common spin for both electrons in the case of spinless electrons.

Figure 3: Solid lines: Local densities of states $\varrho(\mathbf{r};E)$
defined in Eq. (63) as a functions of energy, calculated for four values of
inter-electron distance $r$. Dotted line: $\rho_{f}(E)$, as given in Eq. (73),
as a function of energy. LDOS is in $r_{B}^{-6}$ units. Note that for large
values of $r$, the local density of states $\varrho(\mathbf{r};E)$ approaches
$\rho_{f}(E)$.

To confirm the saturation level in Figure 2c, we calculate in Figure 3 the
LDOS at four values of the electron-electron distance $r$ as a function of
energy, represented by solid lines, and compare them with $\rho_{f}(E)$ in Eq.
(73), denoted by a dotted line. As seen in Figure 3, for high energies and
large electron-electron distances, the LDOS of the electron pair does not
depend on $r$ and tends to the uniform electron density of a pair of free
electrons as given in Eq. (73).

In conclusion, based on the results presented in Figures 2 and 3, we observe
the following quantitative behavior of the local density of states as a
function of energy and electron distance. For low energies and small $r$, the
LDOS is negligibly small but still nonzero. By increasing energy or electron
distance, the LDOS increases following a power-law trend up to an inflection
point. The position of this point depends on energy and shifts to lower values
of $r$ as $E$ increases. Ultimately, for large energies or distances, the LDOS
is uniform and corresponds to the LDOS of a pair of electrons in the absence
of the Coulomb potential.

The physical picture corresponding to the above dependence of LDOS on energy
and inter-electron distance is as follows. For every energy, there is an area
of low electron density around $r=0$, where the Coulomb repulsion between
electrons is strong, and the electrons do not overlap. The size of this area
decreases with energy. By increasing $r$, the Coulomb repulsion weakens, and
electrons start to overlap, leading to an increase in LDOS. Finally, at large
distances, the Coulomb repulsion between electrons vanishes, and they behave
as free particles.

Figure 4: Even part $\varrho_{e}(\mathbf{r};E)$ and odd part
$\varrho_{o}(\mathbf{r};E)$ of local density of states defined in Eqs. (61)
and (62), respectively, as functions of the inter-electron distance $r$. Both
quantities are in $r_{B}^{-6}$ units. Note that $\varrho_{o}(\mathbf{r};E)$
precisely vanishes at $r=0$ due to the Pauli exclusion principle, while
$\varrho_{e}(\mathbf{r};E)$ remains finite in this limit. For large values of
$r$, both quantities are in close proximity to each other.

It is important to observe that the LDOS shown in Figures 2 and 3 does not
approach zero as $r$ tends to $0$. However, for low energies, it is
exceptionally small and can be considered negligible.

Let us now examine the even and odd components of the local density of states,
as defined in Eqs. (61) and (62), corresponding to LDOS for singlet and
triplet states, respectively. In Figure 4, the local densities of states
$\varrho_{e}(r;E)$ and $\varrho_{o}(r;E)$ are plotted for $E=4$ Ha, focusing
on small inter-electron distances. The primary distinction between them lies
in their behavior as $r$ approaches $0$. In this limit, the odd component of
the LDOS precisely converges to zero, as implied by Eqs. (38) and (62).
Conversely, the even component of LDOS remains finite in this limit. As $r$
increases, a significant disparity between $\varrho_{e}(r;E)$ and
$\varrho_{o}(r;E)$ persists until $r>2.5$rB, beyond which the two quantities
become nearly identical.

Figure 5: Panel (a): Local density of states $\varrho_{e}(\mathbf{r};E)$ given
in Eq. (63) (solid line) and local density of states in the absence of spin
effects $\varrho_{s}(r;E)$ given in Eq. (74), as functions of the inter-
electron distance $r$. The difference between these quantities, denoted as
$\varrho_{-}(r;E)$ and defined in Eq. (64), is plotted in Panel (b). Note the
different scale in Panel (b). The term $\varrho_{-}(r;E)$ arises from the
Pauli principle, ensuring the vanishing of the odd part of LDOS at $r=0$ as
discussed in the text. LDOS is in $r_{B}^{-6}$ units.

The discrepancy between the local densities of states $\varrho_{e}(r;E)$ and
$\varrho_{o}(r;E)$ in Figure 4 arises as a consequence of many-body effects,
specifically, it is linked to the term $\varrho_{-}(r;E)$ in Eqs. (61)–(64).
For a more in-depth investigation of this phenomenon, Figure 5 facilitates a
comparison between $\varrho(r,E)$ and the LDOS in the absence of spin effects,
which is defined as

$\varrho_{s}(r;E)=2\varrho_{+}(r;E),$ (74)

The factor of $2$ in Eq. (74) appears due to summation over two spins, see Eq.
(73). As shown in Figure 5a, for $E=4$ Ha, the deviation between
$\varrho(r;E)$ and $\varrho_{s}(r;E)$ is primarily localized at small values
of inter-electron separation $r$. In Figure 5b, the quantity
$\varrho_{-}(r;E)$, predominantly concentrates at small $r$ distances and
exhibits a decaying and oscillatory behavior, diminishing beyond $r>3r_{B}$.
It is worth noting that $\varrho_{-}(r;E)$ cannot be considered a genuine
density of states as it lacks positive definiteness.

This quantity should be interpreted as a many-body contribution to the overall
density of states, $\varrho(r;E)$, arising from electron correlations. Its
emergence can be attributed to the Pauli exclusion principle, which imposes
specific spatial symmetries on electron wave functions concerning the exchange
of particles. The presence of $\varrho_{-}(r;E)$ ensures that, for triplet
states, both the wave function and the LDOS precisely vanish at
$\mathbf{R}=\mathbf{r}=\mathbf{0}$, indicating that two electrons in any
triplet state cannot occupy the same spatial point.

It is noteworthy that the Pauli exclusion principle demonstrates a
considerable potency in comparison to the Coulombic interaction between
electrons. The latter allows for a minor yet non-negligible degree of electron
overlap at the point where the positions of the electrons
$\mathbf{R}=\mathbf{r}=\mathbf{0}$ coincide, as visually depicted in Figure 1.
The necessity for the term $\varrho_{-}(r;E)$ to exist is a direct consequence
of the presence of a non-zero density of states at $r=0$. In the event that
$\varrho_{+}(r;E)$ were to vanish at $r=0$ in Eq. (62), the Pauli exclusion
principle would be intrinsically satisfied by the $\varrho_{+}(r;E)$ term
alone, rendering any corrections at $r=0$ unnecessary. However, owing to the
finite value of $\varrho_{+}(r;E)$ for $r\to 0$, the supplementary term
$\varrho_{-}(r;E)$ becomes indispensable to ensure the complete suppression of
electron overlap at $r=0$.

The Pauli exclusion principle has a limited range of influence, as shown in
Figure 5b. The term $\varrho_{-}(r;E)$, which reflects the Pauli exclusion
effect, is predominantly concentrated at small inter-electron distances $r$
and diminishes rapidly beyond a certain point. Consequently, for large values
of $r$, the density of states $\varrho_{t}(r;E)$ converges to
$\varrho_{s}(r;E),$ as depicted in Figure 5a. Similarly, at larger $r,$
there’s no distinction between $\varrho_{+}(r;E)$ and $\varrho_{-}(r;E)$, as
demonstrated in Figure 4.

##  VI Discussion

In the previous sections, various issues related to the obtained results have
already been discussed. In this section, we shift our focus to other aspects
of the considered problem.

In our method, we analytically evaluate the double summation over $l$ and $m$
and integrate over $k$ in Eq. (9) using the Coulomb GF. This yields the GF
expressed as an integral over
$\tilde{g_{c}}(\mathbf{r}_{1},\mathbf{r}_{2};\epsilon_{K})$ with energy
$\epsilon_{K}$ depending on the wave vector $\mathbf{K}$. Alternatively, it’s
possible to reverse the order of integration by considering the variable
$\mathbf{K}$ first. This alternative approach leads to the same GF results as
those presented in Eq. (10), albeit in a more tedious way.

The results in this paper pertain to three-dimensional space but can be
extended to one-dimensional systems as the closed-form Coulomb GF is known for
$1D$ [2]. Generalization is also possible for systems with dimensions
$D=5,7,\ldots$ due to a discovered relationship between Coulomb GF in
different dimensions [7]

$\tilde{g_{c}}(x,y;E;D+2)=-\frac{1}{2\pi y}\frac{\partial}{\partial
y}\tilde{g_{c}}(x,y;E;D).$ (75)

This relation applies for integer values of $D\geq 1$. However, there is no
closed-form Coulomb GF for $2D$ systems. Thus, our method extends to systems
with dimensions $D=2,4,6,\ldots$ only when employing the partial-wave
expansion of the Coulomb GF, as shown in Eq .(15). The two-dimensional
equivalent of Eq. (15) can be found in Ref. [7].

The LDOS calculations in Figures 2–5 used an analytical expression for the
Coulomb GF, as shown in Eq. (II.1). However, alternative representations of
the Coulomb GF are also applicable. These include the partial waves expansion
as indicated in Eq. (15), the momentum representation developed by Schwinger
[10], or representations involving hyperbolic functions, as referenced in Ref.
[13]. It is worth noting that deriving an explicit closed-form expression for
the GF for two particles seems unattainable, and in any case, one typically
arrives at the GF in the form of a definite integral.

Regarding the numerical aspects of the calculation, the Whittaker functions in
Eq. (II.1) are computed using series expansions for small arguments [21], and
for large arguments, they are obtained from corresponding asymptotic
expansions, as described in the reference [21].

Additionally, the guidelines provided in Ref. [25] for calculating confluent
hypergeometric functions have been taken into consideration. Accuracy was
ensured by monitoring the Wronskian of the two Whittaker functions in Eq. (70)
for each parameter $k$, $\nu$, and $r$ in Eq. (II.1). For fixed $\nu$ this
Wronskian value should remain constant, independent of the arguments of the
Whittaker functions, providing a dependable metric for assessing accuracy.

The method presented in this paper is applicable, with some modifications, to
two-particle molecules subject to attractive Coulomb potentials such as
excitons, positronium, and muonium, among others. The key distinctions in the
GF calculations for these systems, compared to the results in this paper, are
as follows: i) Existence of bound states due to attractive Coulomb
interactions in these systems. ii) Absence of the Pauli exclusion principle,
because of the opposite charge signs in these systems. iii) Transition from a
negative parameter $\nu$ in the Coulomb GF in Eq. (II.1) for repulsive Coulomb
potentials to a positive parameter $\nu$ for attractive potentials, leading to
different behavior of the Whittaker functions within the Coulomb GF. In
addition to these primary differences, other system-specific factors may come
into play, such as varying electron and hole effective masses in excitons or
finite recombination lifetimes in positronium. Nevertheless, in principle, our
approach remains applicable to these systems as well.

The method outlined in our paper is specifically designed for two-particle
systems. It excels at separating the motion of two electrons into center-of-
mass and relative motion (as detailed in Eq. (2)), making it universally
applicable to all two-particle systems with interactions depending on
$|\mathbf{r}_{1}-\mathbf{r}_{2}|$. However, this approach cannot be extended
to systems with three or more particles. The fundamental reason is that the
motion separation, which forms the basis of our approach, is infeasible in
systems involving three or more particles. Therefore, our method is not well-
suited for computing the GF in multi-particle systems.

When either $\mathbf{R}_{1}=\mathbf{R}_{2}$ or
$\mathbf{r}_{1}=\mathbf{r}_{2}$, we calculate the real part of the GF as a
Hilbert transform of the imaginary part with a cutoff energy $W$, see Eq.
(30). In the case of a pair of non-relativistic electrons in a vacuum, there
isn’t a natural or intrinsic cutoff energy. One potential choice for a cutoff
is $W\approx m_{e}c^{2}=0.511$ MeV, which exceeds the characteristic energy of
a pair, typically on the order of $27.21$ eV, by four orders of magnitude. For
electrons confined in a quantum dot with barrier height $U_{0}$, choosing
$W\approx U_{0}$ is a reasonable cutoff. However, the selection of a cutoff
energy is model-specific and should be rigorously justified based on the
system’s characteristics. Cutoff determination lacks a universal method; it
hinges on the unique attributes of the model under investigation.

For experimental observation of LDOS in Figures 2–5, one viable approach is
measuring the LDOS of an electron pair within quantum dots. Figure 2
demonstrates that the dot’s size should exceed the effective Bohr radius
$r_{B}^{*}$ of electrons within the dot, given by [26]

$r_{B}^{*}=r_{B}\kappa\frac{m_{e}}{\mu^{*}},$ (76)

where $\kappa$ represents the material’s dielectric constant, and $\mu^{*}$ is
the reduces effective mass of electron pair. Typically, $r_{B}^{*}$ falls in
the range of approximately $50-100$Å. In such a system, the presence of dot
barriers may minimally affect pair motion. The LDOS can be observed using
techniques like scanning tunneling microscopy, scanning tunneling
potentiometry, scanning tunneling spectroscopy, or tip-enhanced Raman
spectroscopy, see e.g. [27, 28]. In addition to the total LDOS, it is feasible
to measure the LDOS for singlet and triplet states of electrons within the
quantum dot. To control the spin states of an electron pair, various methods
can be employed, such as electron tunneling between electrodes, optoelectronic
and magneto-optical techniques, or injecting electrons in specific spin states
from an external source.

Regardless of the experimental method used, the system allows for the
measurement of the following effects: i) The overall shape of the LDOS
concerning energy and inter-electron distance, as shown in Figure 2. ii) The
free-particle limit of the LDOS at sufficiently high energies, demonstrated in
Figures 2 and 3. iii) The distinction between the LDOS for singlet and triplet
states of the electron pair, as seen in Figure 4. iv) A complete absence of
the LDOS at the origin, i.e., when the inter-electron distance is $r=0$. These
effects remain robust even in the presence of confinement resulting from the
quantum dot structure.

##  VII Summary

In this paper, we extended the Coulomb Green’s function to a system involving
two electrons interacting through repulsive Coulomb forces. We derived closed-
form expressions for the GF, which are represented by one-dimensional
integrals, as shown in Eqs. (18) to (II.2). It’s important to emphasize that
these equations are valid for systems where electron spins are not considered.
The obtained GF has no poles, and no bound states exist.

For positive energies, the obtained GF comprises a complex oscillatory term
with a non-vanishing imaginary component and a real term that decays
exponentially with inter-electron distance. For negative energies, the GF is
real and also decays exponentially with inter-electron distance.

For certain combinations of GF arguments, specifically when
$\mathbf{R}_{1}=\mathbf{R}_{2}$ or $\mathbf{r}_{1}=\mathbf{r}_{2}$, the
integrals for the real parts of the GF diverge. However, the imaginary parts
of the GF remain finite, resulting in a finite local density of states. We
examined these cases in detail, including situations where
$\mathbf{R}_{1}=\mathbf{R}_{2}=\mathbf{r}_{1}=\mathbf{r}_{2}=\mathbf{0}$, as
illustrated in Figure 1. It has been discovered that for any pair energy, the
LDOS at this point remains finite, indicating a non-zero overlap of electrons
wave functions. This phenomenon lacks a classical counterpart. We also
considered the scenario where $\mathbf{r}_{1}=\mathbf{r}_{2}\neq\mathbf{0}$,
as described in Eq. (72).

In the subsequent phase, we further generalize the results to include the
spins of electrons and account for the influence of the Pauli exclusion
principle. In this context, we discovered that the GF is composed of two
terms, each characterized as either odd or even concerning the exchange of
particles, as described in Eq. (33). We derived closed-form expressions for
the even and odd GF as sums and differences of GFs with the appropriate
arguments, as outlined in Eqs. (43) to (53).

We also discovered that for spin-independent potentials, the Dyson equation
separates into two distinct equations, one for even and the other for odd GFs,
respectively. After obtaining the GF for an electron pair in the presence of
spins, we proceeded to calculate the LDOS for the system. The LDOS is a sum of
contributions from the even and odd parts of the GF, as described in Eq. (60).
In Figures 2 to 5, we computed the LDOS as a function of inter-electron
distance and the pair’s energy. Additionally, we separately calculated the odd
and even contributions to the LDOS, highlighting the significance of the Pauli
exclusion principle.

We further investigated the pseudo-local density of states, denoted as
$\varrho_{-}(\mathbf{r};E)$, which signifies the many-body contribution to the
GF and guarantees the total suppression of the local density of states at
$r=0$. The necessity of including this term arises from the non-zero local
density of states at $r=0$, as depicted in Eq. (62) and Figure 4. This term
exhibits a relatively limited spatial extent and diminishes as the inter-
electron distances increase.

We hope our paper facilitates a deeper understanding of the non-relativistic
electrons interacting through repulsive Coulomb forces and underscores the
significance of the Pauli exclusion principle in few-electron systems.

## Appendix A

The retarded GF for a free particle with positive energies is

$g^{+}_{e}(\mathbf{r}_{1},\mathbf{r}_{2};E)=-\frac{1}{4\pi
r}\exp(i\sqrt{E}r),$ (77)

where $r=|\mathbf{r}_{1}-\mathbf{r}_{2}|$. In the limit $r\to 0$ there is
$\rho_{e0}(E)=\sqrt{E}/(4\pi^{2})$. For the Coulomb GF corresponding to a
repulsive potential, we obtain, see Eq. (26),

$\rho_{c0}(E)=-\frac{1}{\pi}\int_{0}^{\infty}\frac{f(k)dk}{E-c_{k}k^{2}+i\eta}=\frac{1}{2\pi}\frac{f(\sqrt{E})}{\sqrt{E}},$
(78)

where $c_{k}=1$ and we have utilized Eq. (28).

We will now proceed to calculate the LDOS for a system of two electrons in
absence of Coulomb interaction. At the specific point
$\mathbf{R}=\mathbf{r}=\mathbf{0}$, the LDOS for a noninteracting electron
pair, denoted as $\rho_{f0}(E)$, is given by

$\rho_{f0}(E)=\frac{-1}{\pi}\frac{1}{4\pi^{4}}\int_{0}^{\infty}\int_{0}^{\infty}\frac{k_{a}^{2}k_{b}^{2}dk_{a}dk_{b}}{E-\frac{1}{2}k_{a}^{2}-\frac{1}{2}k_{b}^{2}+i\eta}.$
(79)

By introducing polar coordinates, where $(k_{a},k_{b})\to(t,\alpha)$ and
applying Eq. (28), we obtain

$\displaystyle\rho_{f0}(E)$ $\displaystyle=$
$\displaystyle\frac{1}{2\pi^{4}}\int_{0}^{\pi/2}\int_{0}^{\infty}t^{5}\cos^{2}(\alpha)\sin^{2}(\alpha)^{2}\delta(2E-t^{2})dt$
(80) $\displaystyle=$ $\displaystyle\frac{E^{2}}{16\pi^{3}}\Theta(E),$

where the step function $\Theta(E)$ ensures that the integral is non-zero only
for $E>0$, and vanishes for $E<0$.

For arbitrary positions $\mathbf{R}$ and $\mathbf{r}$, the LDOS for a pair of
non-interacting electrons also can be obtained in a closed form. The retarded
GF is

$g^{+}_{f}(\mathbf{R},\mathbf{r};E)=\frac{1}{(2\pi)^{6}}\int\frac{e^{i\mathbf{K}\mathbf{R}}e^{i{\mathbf{k}}\mathbf{r}}d^{3}\mathbf{K}d^{3}{\mathbf{k}}}{E-c_{K}K^{2}-c_{k}k^{2}+i\eta}.$
(81)

By integrating over $d^{3}{\mathbf{k}}$, we arrive at the expression

$g^{+}_{f}(\mathbf{R},\mathbf{r};E)=\frac{1}{(2\pi)^{3}}\int
e^{i\mathbf{K}\mathbf{R}}g^{+}_{1e}(\mathbf{r},\epsilon_{K})d^{3}\mathbf{K},$
(82)

where $\epsilon_{K}=E-c_{K}K^{2}$, and $g^{+}_{1e}(\mathbf{r},E)$ is given in
Eq. (77). For $E<0$, the imaginary part of $g^{+}_{0}(\mathbf{r};E)$ vanishes,
leading to a reduction of the LDOS for negative energies. For $E>0$ we have

$\displaystyle{\rm Im}\big{\\{}g^{+}_{f}(\mathbf{R},\mathbf{r};E>0)\big{\\}}=$
$\displaystyle=-\frac{1}{2\pi^{2}}\int_{0}^{\infty}\left[\frac{\sin(r\epsilon_{K})}{4\pi
r}\right]\left(\frac{K\sin(KR)}{R}\right)dK.$ (83)

To evaluate the above integral, we use the identity (2.5.25.1) in Ref. [29]

$\displaystyle\int_{0}^{a}\sin(r\sqrt{a^{2}-K^{2}})\cos(qK)dK=$
$\displaystyle=\frac{\pi}{2}\frac{ar}{\sqrt{q^{2}+r^{2}}}J_{1}(a\sqrt{q^{2}+r^{2}}),$
(84)

By differentiating both sides of Eq. (A) with respect to $dq$ and introducing
$E^{\prime}=E/c_{K}=4E$, we derive the following from Eq. (A).

$\displaystyle\varrho_{f}(\mathbf{R},\mathbf{r};E)$ $\displaystyle=$
$\displaystyle\frac{1}{16\pi^{3}}\left[-\frac{E^{\prime}J_{0}\big{(}t\sqrt{E^{\prime}}\big{)}}{2t^{2}}+\frac{E^{\prime}J_{2}\big{(}t\sqrt{E^{\prime}}\big{)}}{2t^{2}}\right.+$
(85)
$\displaystyle\left.+\frac{\sqrt{E^{\prime}}J_{1}\big{(}t\sqrt{E^{\prime}}\big{)}}{t^{3}}\right],$

where $t=\sqrt{R^{2}+r^{2}}$, and $E>0$. It is worth noting that the limits
$\mathbf{R}\to\mathbf{0}$ and $\mathbf{r}\to\mathbf{0}$ yield the LDOS as
given in Eq. (80).

To calculate the traces over $\Lambda_{ss}$ and $\Lambda_{tt}$ in Eqs. (61)
and (62), we consider a singlet state $|s\rangle$ for an electron pair and a
vector of triplet states
${\mathbf{t}}=[|t_{1}\rangle,|t_{2}\rangle,|t_{3}\rangle]$ for the pair. In
this context, the trace of $\Lambda_{s}=|s\rangle\langle s|$ equals unity.
Meanwhile, the matrix $\Lambda_{t}={\mathbf{t}}\cdot{\mathbf{t}}^{\dagger}$ is
a $3\times 3$ identity matrix, and its trace equals three.

## References

  * [1] C. Hostler and R. H. Pratt, Phys. Rev. Lett. 10, 469 (1963).
  * [2] I. Meixner, Math. Z. 36, 677 (1933).
  * [3] R. A. Mapleton, J. Math. Phys. 2, 478 (1961).
  * [4] K. Mano, J. Math. Phys. 4, 522 (1963).
  * [5] L. Hostler, J. Math. Phys. 5, 591 (1964).
  * [6] L. Hostler, J. Math. Phys. 8, 642 (1967).
  * [7] L. Hostler, J. Math. Phys. 11, 2966 (1970).
  * [8] H. F. Hameka, J. Chem. Phys. 47, 2728 (1967).
  * [9] H. F. Hameka, J. Chem. Phys. 48, 4810 (1968).
  * [10] J. Schwinger, J. Math. Phys. 5, 1606 (1964).
  * [11] S. M. Blinder, J. Math. Phys 22, 306 (1981).
  * [12] R. A. Swainson and G. W. F. Drake, J. Phys. A: Math. Gen. 24, 79 (1991).
  * [13] A. Maquety, V. Veniard, and T. A. Marianz, J. Phys. B: At. Mol. Opt. Phys. 31, 3743 (1998).
  * [14] Robert N. Hill, in Springer Handbook of Atomic, Molecular, and Optical Physics, ed. G. W. F. Drake, (New York: Springer, 2006).
  * [15] Z. Papp, C.-Y. Hu, Z. T. Hlousek, B. Konya, and S. L. Yakolev, Phys. Rev. A 63, 062721 (2001).
  * [16] Z. Papp, J. Darai, C.-Y. Hu, Z. T. Hlousek, B. Konya, and S. L. Yakolev, Phys. Rev. A 65, 032725 (2002).
  * [17] R. Shakeshaft, Phys. Rev. A 70, 042704 (2004).
  * [18] D. Feranchuk and V. V. Triguk Phys. Lett. A 375, 2550 (2011).
  * [19] L. D. Landau and E. M. Lifshitz Quantum Mechancs, Non-relativistic Theory (3rd ed.) (Oxford: Pergamon Press 1977), Section 36, Secion 33.
  * [20] H. Buchholz, The Confluent Hypergeometric Function (Springer, New York, 1969). Note a diffecence between functions $M^{\mu}_{\kappa}(z)$ and ${\cal M}^{\mu}_{\kappa}(z)$ used in this work.
  * [21] http://dlmf.nist.gov/13.14; http://dlmf.nist.gov/13.19, (2023).
  * [22] A. Zon, N. L. Manakov, and L. I. Rapoport, Sov. Phys. JETP 28, 480 (1969).
  * [23] https://en.wikipedia.org/wiki/Table_of_spherical_harmonics (2023).
  * [24] W. H. Press, S. A. Teukolsky, W. T. Vetterling, and B. P. Flannery Numerical Numerical Recipes: The Art of Scientific Computing (3rd ed.), (New York: Cambridge 2007), Section 4.
  * [25] J. W. Pearson, S. Olver, and M. A. Porter, Numer. Algor. 74, 821 (2017).
  * [26] https://en.wikipedia.org/wiki/Quantum_dot, (2023).
  * [27] S. Schintke and W. D. Schneider, J. Phys.: Cond. Matt. 16, R49 (2004).
  * [28] J. Kroger, N. Neel, and L. Limot, J. Phys.: Cond. Matt. 20, 223001 (2008).
  * [29] A. P. Prudnikov, Yu. A.Brychkov, and O. I. Marichev, Integrals and Series, Elementary Functions, Vol. 1 (Taylor and Francis: London 1998), Eq. (2.5.25.1).