**2. Calculation Technique**

The simplicity of the WFs under consideration implies that the form

$$f\_N(r) = \sum\_{k=1}^{N} \mathbb{C}\_k \exp(-b\_k r) \tag{1}$$

represents the sum of a few single exponentials, whereas the compactness means that their number 3 ≤ *N* ≤ 6 in Equation (1), unlike the basis sizes mentioned in the introduction. The relevant accuracy will be discussed later. It is seen that the RHS of Equation (1) includes *N* linear parameters *Ck* and *N* nonlinear parameters *bk* with *k* = 1, 2, . . . , *N*.

The *collinear* arrangemen<sup>t</sup> of the particles consisting of the nucleus and two electrons can be described by a single scalar parameter *λ* as follows [1]:

$$r\_1 = r, \qquad r\_2 = |\lambda|r, \qquad r\_{12} = (1 - \lambda)r,\tag{2}$$

where *λ* ∈ [−1, 1], and *r* is the distance between the nucleus and the electron most distant from it. Clearly *λ* = 0 corresponds to the *electron–nucleus* coalescence, and *λ* = 1 to the *electron–electron* coalescence. The boundary value *λ* = −1 corresponds to the collinear **e-n-e** configuration with the same distances of both electrons from the nucleus. In general, 0 < *λ* ≤ 1 corresponds to the collinear arrangemen<sup>t</sup> of the form **n-e-e** where both electrons are on the same side of the nucleus. Accordingly, −1 ≤ *λ* < 0 corresponds to the collinear arrangemen<sup>t</sup> of the form **e-n-e** where the electrons are on the opposite sides of the nucleus. The absolute value |*λ*| measures the ratio of the distances of the electrons from the nucleus.

Thus, for the particles with collinear arrangemen<sup>t</sup> we can introduce the *collinear* WF of the form

$$\Phi(r,\lambda) \equiv \Psi(r, |\lambda| r, (1-\lambda)r) / \Psi(0,0,0). \tag{3}$$

It should be emphasized that, e.g., the PLM WF with *collinear* configuration reduces to the form

$$\Phi\_{PLM}(r,\lambda) = \exp(-\delta\_\lambda r) \sum\_{p=0}^{\Omega} c\_p(\lambda) r^p \,. \tag{4}$$

where Ω = 25 for the current (standard) consideration, as it was mentioned earlier.

We can give an example of the physical problem where the collinear WF of the form (4) cannot be applied, but the quite accurate WF of the form (1) is required instead. In Refs. [2,3], the mechanism of photoionization in the two-electron atoms is investigated. Calculations of various differential characteristics (cross sections) of ionization are based on computation of the triple integral of the form

$$\int d^3 \mathbf{r} \, e^{i\mathbf{q}\mathbf{r}} \,\_1F\_1(i\mathbf{j}\_1^\mathbf{r}, 1, i\mathbf{p}\_1r - i\mathbf{p}\_1\mathbf{r}) \,\_1F\_1(i\mathbf{j}\_2^\mathbf{r}, 1, i\mathbf{p}\_2r - i\mathbf{p}\_2\mathbf{r}) \Phi(r, 1),\tag{5}$$

where **p***j* (*j* = 1, 2) are the momenta of photoelectrons, **q** is the recoil momentum, *ξj* = *<sup>Z</sup>*/*pj*, *i* is the imaginary unit, and <sup>1</sup>*F*1(...) is the confluent hypergeometric function of the first kind. The most important for our consideration is the fact that integral (5) contains the *collinear* WF <sup>Φ</sup>(*<sup>r</sup>*, 1) describing the case of the electron–electron coalescence (*λ* = 1) in the helium-like atom/ion with the nucleus charge *Z*. It is clear that the numerical computation

of the triple integral (5) is not impossible, but rather a difficult problem, especially for building the relevant graphs. Fortunately, already in 1954 [18], the explicit expression for the triple integral which is very close to integral (5) was derived. In fact, integral (5) can be calculated by simple differentiation (with respect to a parameter) of the explicit form for the integral mentioned above, but only under condition that the WF, <sup>Φ</sup>(*<sup>r</sup>*, 1) is represented by a single exponential of the form exp(−*br*) (with positive parameter *b*, of course).

According to the Fock expansion [19,20] (see also [21,22]), we have:

$$\Psi(r\_1, r\_2, r\_{12}) / \Psi(0, 0, 0) \underset{R \to 0}{=} 1 - Z(r\_1 + r\_2) + \frac{1}{2}r\_{12} - Z\left(\frac{\pi - 2}{3\pi}\right) \left(R^2 - r\_{12}^2\right) \ln R + O(R^2),\tag{6}$$

where *R* = (*r*21 + *r*22)1/2 is the hyperspherical radius. Using Equation (6) and the collinear conditions (2), we obtain the Fock expansion for the collinear WF in the form:

$$\Phi(r,\lambda) \underset{r \to 0}{=} 1 + \eta\_{\lambda}r + \zeta\_{\lambda}r^{2}\ln r + \zeta\_{\lambda}r^{2} + \dots \tag{7}$$

where

$$
\eta\_{\lambda} = -Z(1 + |\lambda|) + \frac{1 - \lambda}{2},
\tag{8}
$$

$$\mathcal{Z}\_{\lambda} = -\frac{2Z\lambda(\pi - 2)}{3\pi},\tag{9}$$

and the general form of the coefficient *ξλ* being rather complicated will be discussed later. The necessity of the equivalent behavior of the model WF, (1) and the variational WF, <sup>Φ</sup>(*<sup>r</sup>*, *λ*) near the nucleus (*r* → 0) results in the following two coupled equations for 2*N* parameters **C***N* ≡ {*<sup>C</sup>*1, *C*2,..., *CN*} and **b***N* ≡ {*b*1, *b*2,..., *bN*} of the model WF:

$$\sum\_{k=1}^{N} \mathbb{C}\_{k} = 1,$$

$$\sum\_{k=1}^{N} \mathbb{C}\_{k} b\_{k} = Z(1 + |\lambda|) + \frac{\lambda - 1}{2}. \tag{11}$$

Equation (10) follows from the condition Φ(0, *λ*) = 1, whereas Equation (11) is obtained by equating the linear (in *r*) coefficients of the power series expansion of the model WF (1) and the Fock expansion (7).

As it was mentioned above, to obtain the fully defined model WF of the form (1) one needs to determine 2*N* coefficients. To solve the problem with given Equations (10) and (11), we need to find extra 2(*N* − 1) coupled equations for parameters of the exponential form (1). To this end, we propose to use the definite integral properties of the collinear WF (3).

A number of numerical results presenting expectation values of Dirac-delta functions *δ*(**<sup>r</sup>**1)-, *δ*(**<sup>r</sup>**12)- ≡ *δ*(**<sup>r</sup>**1 − **<sup>r</sup>**2)- and *δ*(**<sup>r</sup>**1)*δ*(**<sup>r</sup>**2)- for the helium-like atoms can be found in the proper scientific literature (see, e.g., [15,17,23] and references therein). It was shown [1] that expectation values mentioned above represent the particular cases of the more general expectation value

$$
\langle \delta(\mathbf{r}\_1 - \lambda \mathbf{r}\_2) \rangle = 4\pi \langle \delta(\mathbf{r}\_1) \delta(\mathbf{r}\_2) \rangle \int\_0^\infty |\Phi(r, \lambda)|^2 r^2 dr,\tag{12}
$$

where

$$
\langle \delta(\mathbf{r}\_1) \delta(\mathbf{r}\_2) \rangle = \Psi^2(0,0,0) / \int \psi^2(\mathbf{r}\_1,\mathbf{r}\_2) d^3 \mathbf{r}\_1 d^3 \mathbf{r}\_2 \tag{13}
$$

is a square of the normalized WF taken at the nucleus. It is seen that the expectation value (12) is fully defined by the *collinear* WF, <sup>Φ</sup>(*<sup>r</sup>*, *<sup>λ</sup>*).

We propose to use the integrals of the form

$$S\_{\mathfrak{n}} = \int\_0^\infty |\Phi(r, \lambda)|^2 r^n dr, \qquad (n = 0, 1, 2, \dots) \tag{14}$$

for deriving 2(*N* − 1) extra coupled equations required, in its turn, for determining 2*N* coefficients defining the model WF, (1). Replacing <sup>Φ</sup>(*<sup>r</sup>*, *λ*) in the RHS of Equation (14) by the model WF (1) and using the closed form of the corresponding integral, one obtains *n* equation of the form

$$S\_{ll} = n! \sum\_{j=1}^{N} \sum\_{k=1}^{N} \frac{\mathbb{C}\_{j}\mathbb{C}\_{k}}{(b\_{j} + b\_{k})^{n}},\tag{15}$$

where, in fact, **C***N* ≡ **<sup>C</sup>***N*(*<sup>Z</sup>*, *<sup>λ</sup>*), **b***N* ≡ **b***N*(*<sup>Z</sup>*, *λ*) are the coefficients we are requested, whereas the integrals *Sn* ≡ *Sn*(*<sup>Z</sup>*, *λ*) can be computed using, for example, the PLM WFs according to definition (14). The technique proposed, in fact, represents a variant of the "Method of Moments" (see, e.g., [24]) supplemented by the boundary conditions (10) and (11). 

The problem is that it is necessary to select a set (sample) of integers 5*n*1, *n*2,..., *<sup>n</sup>*2(*<sup>N</sup>*−<sup>1</sup>)6describing Equations (14) and (15) for each triple of numbers (*<sup>Z</sup>*, *λ*, *<sup>N</sup>*). Those selected samples are presented in Tables 1–4, along with the corresponding parameters of the model WFs.

**Table 1.** Parameters of the model WFs *f*3(*r*).


**Table 2.** Parameters of the model WFs *f*4(*r*).


5

9.54371206725

0.878994353474

11.8356042739

0.109401348385


**Table 2.** *Cont.*

22.5015053896

0.0108088103065

91.8749584959

0.000795487835141

 2, 3, 4, 5, 6, 7

> 0.91


**Table 3.** Parameters of the model WFs *f*5(*r*).

**Table 4.** Parameters of the model WFs *f*6(*r*) for the negative ion of hydrogen (*Z* = 1).


To solve the set of Equations (10), (11) and 2(*N* − 1) nonlinear equations of the form (15) we apply, as the first step, the built-in function **NSolve[**... **]** of the Wolfram *Mathematica*. The additional conditions (inequalities) **b***N* > 0 are used. The program **NSolve** generates all possible solutions. However, only one of them represents the nodeless solution that

corresponds to the ground-state WF. We have computed and presented the parameters of the model WFs for 3 ≤ *N* ≤ 6. It was mentioned above that the **NSolve** is used only at the first step. The reason is that this program works normally (with no problems) only for *N* ≤ 3, that is for number of equations 2*N* ≤ 6. Even for 2*N* = 6 computer freezes for a few second capturing 100 % of CPU time, and then normal operation is restored. However, for 2*N* = 8 *Mathematica* (through **NSolve**) takes all CPU time, and computer freezes for an indefinite time. This is happened for any settings of *Mathematica*, e.g., for any settings in "Parallel Kernel Configuration". We checked that this problem persists in different computers and for different version of *Mathematica* (9, 10.3, 11.0, 12.1). Therefore, to solve the relevant set of nonlinear equations for the number of exponentials *N* > 3 we employed the built-in (*Mathematica*) program **FindRoot[**... **]**. Unlike **NSolve** this program generates only one solution (if it exists, of course) starting its search from some initial values **C**(0) *N* , **b**(0) *N* for which we take the values **C***N*−1, **b***N*−<sup>1</sup> of the corresponding calculation on the *N* − 1 exponentials. The conditions of the positive exponents and the WF nodeless are certainly preserved.

To estimate the accuracy of the model WF we employ the following integral representation

$$R\_N = \int\_0^\infty r |f\_N(r) - \Phi(r, \lambda)| dr \left( \int\_0^\infty r |\Phi(r, \lambda)| dr \right)^{-1}. \tag{16}$$

Note that the function *r*Φ is more indicative than Φ, at least, for the ground state.
