**1. Introduction**

The properties of a two-electron atomic (helium-like) system with an infinitely massive nucleus of charge *Z* and nonrelativistic energy *E* are defined by the wave function (WF) <sup>Ψ</sup>(*<sup>r</sup>*1,*r*2,*r*12), where *r*1 and *r*2 are the electron–nucleus distances, and *r*12 is the distance between the electrons. The behavior of the ground state WF in the vicinity of the nucleus located at the origin is determined by the Fock expansion [1]

$$\Psi(r\_1, r\_2, r\_{12}) \equiv \Psi(R, a, \theta) = \sum\_{k=0}^{\infty} R^k \sum\_{p=0}^{\left[k/2\right]} \Psi\_{k,p}(a, \theta) \ln^p R,\tag{1}$$

where the hyperspherical coordinates *R*, *α* and *θ* are defined by the relations:

$$\theta = \sqrt{r\_1^2 + r\_{2'}^2} \quad \text{or} \quad \theta = 2 \arctan\left(\frac{r\_2}{r\_1}\right), \quad \theta = \arccos\left(\frac{r\_1^2 + r\_2^2 - r\_{12}^2}{2r\_1r\_2}\right). \tag{2}$$

The convergence of expansion (1) was proven in Ref. [2]. The angular Fock coefficients (AFC) *ψ<sup>k</sup>*,*<sup>p</sup>* satisfy the Fock recurrence relation (FRR)

$$\left[\Lambda^2 - k(k+4)\right] \psi\_{k,p}(\mathfrak{a}, \theta) = h\_{k,p}(\mathfrak{a}, \theta) \tag{3}$$

with the RHS of the form [3,4]:

*R* 
$$h\_{k,p} = 2(k+2)(p+1)\psi\_{k,p+1} + (p+1)(p+2)\psi\_{k,p+2} - 2V\psi\_{k-1,p} + 2E\psi\_{k-2,p}.\tag{4}$$

The dimensionless Coulomb potential representing the electron–electron and electron– nucleus interactions is

$$V \equiv \frac{R}{r\_{12}} - Z\left(\frac{R}{r\_1} + \frac{R}{r\_2}\right) = \frac{1}{\tilde{\xi}} - \frac{2Z\eta}{\sin \alpha'}\tag{5}$$

where we have introduced the important (in what follows) angular quantities:

$$
\zeta = \sqrt{1 - \sin a \cos \theta}, \quad \eta = \sqrt{1 + \sin a}. \tag{6}
$$

**Citation:** Liverts, E.Z.; Krivec, R. Fock Expansion for Two-Electron Atoms: High-Order Angular Coefficients. *Atoms* **2022**, *10*, 135. https://doi.org/10.3390/ atoms10040135

Academic Editors: Anatoli Kheifets, Gleb Gribakin and Vadim Ivanov

Received: 23 September 2022 Accepted: 2 November 2022 Published: 7 November 2022

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

**Copyright:** © 2022 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

The hyperspherical angular momentum operator, projected on *S* states, is defined as

$$
\Lambda^2 = -4 \left[ \frac{\partial^2}{\partial a^2} + 2 \cot a \frac{\partial}{\partial a} + \frac{1}{\sin^2 a} \left( \frac{\partial^2}{\partial \theta^2} + \cot \theta \frac{\partial}{\partial \theta} \right) \right]. \tag{7}
$$

It is clear that all circumnuclear features of the two-electron atoms (ions) are defined by the Fock expansion (1). There are a large number of methods for calculating the electronic structure of the two-electron atomic systems. An excellent review on this topic can be found in Refs. [3,5–8]. However, we know only one technique that correctly represents the WF <sup>Ψ</sup>(*<sup>r</sup>*1,*r*2,*r*12) near the nucleus. It is the so-called correlation function hyperspherical harmonic method (CFHHM) [9–11]. The expansion in hyperspherical harmonics (HHs) provides the correct representation of the AFCs. However, the HH expansion is known to converge very slowly. Although this method makes it possible to increase the convergence of the HH expansion, a sufficiently good accuracy requires a large HHs' basis size, which, in turn, creates grea<sup>t</sup> computational difficulties.

Thus, it would be extremely useful to develop a much simpler method for calculating the WF with correct behavior near the nucleus. In this regard, we would like to emphasize the following important peculiarities of the Fock expansion (FE). It follows from definition (1) that the FE can be split into individual power series (lines) associated with definite power of ln *R*. In other words, the FE can be represented in the form:

$$\begin{aligned} \Psi &= (\ln R)^0 \left( \psi\_{0,0} + R\psi\_{1,0} + R^2 \psi\_{2,0} + \dots \right) \\ &+ (\ln R)^1 R^2 \left( \Psi\_{2,1} + R\psi\_{3,1} + R^2 \psi\_{4,1} + \dots \right) \\ &+ (\ln R)^2 R^4 \left( \psi\_{4,2} + R\psi\_{5,2} + R^2 \psi\_{6,2} + \dots \right) \\ &+ (\ln R)^3 R^6 \left( \psi\_{6,3} + R\psi\_{7,3} + R^2 \psi\_{8,3} + \dots \right) \\ &+ (\ln R)^4 R^8 \left( \psi\_{8,4} + R\psi\_{9,4} + R^2 \psi\_{10,4} + \dots \right) + \dots \end{aligned} \tag{8}$$

It is seen that the leading term of each line represents the product (ln *R*)*<sup>k</sup>*/2 *<sup>R</sup><sup>k</sup>ψ<sup>k</sup>*,*k*/2(*<sup>α</sup>*, *θ*) with even *k*. The first AFCs ( *ψ*0,0 = 1) corresponding to *k* = 0, 2, 4 are well-known (see, e.g., [3,4]):

$$
\psi\_{1,0} = \frac{1}{2}\zeta - Z\eta\_{\prime} \tag{9}
$$

$$\psi\_{2,1} = -\frac{Z(\pi - 2)}{3\pi}(1 - \xi^2),\tag{10}$$

$$\psi\_{3,1} = \frac{Z(\pi - 2)}{36\pi} \left[ 6Z\eta (1 - \xi^2) + \xi (5\xi^2 - 6) \right],\tag{11}$$

$$\psi\_{4,2} = \frac{Z^2(\pi - 2)(5\pi - 14)}{540\sqrt{\pi}} \left[ \chi\_{40}(\mathfrak{a}, \theta) + \sqrt{2} \ \chi\_{42}(\mathfrak{a}, \theta) \right]. \tag{12}$$

The normalized HHs are

$$Y\_{40}(\mathfrak{a}, \theta) = \pi^{-3/2} (4 \cos^2 \mathfrak{a} - 1), \quad Y\_{42}(\mathfrak{a}, \theta) = 2 \sqrt{2} \pi^{-3/2} \sin^2 \mathfrak{a} P\_2(\cos \theta), \tag{13}$$

where the *Pn*(*x*) denote the Legendre polynomials.

In this paper, we present the theoretical calculations of the AFCs *ψ*5,2(*<sup>α</sup>*, *<sup>θ</sup>*), *ψ*6,3(*<sup>α</sup>*, *<sup>θ</sup>*), *ψ*7,3(*<sup>α</sup>*, *<sup>θ</sup>*), *ψ*8,4(*<sup>α</sup>*, *θ*) and *ψ*9,4(*<sup>α</sup>*, *θ*) included into the *k* = 4, *k* = 6 and *k* = 8 "lines" of the expansion (8), and also the AFC *ψ*10,5(*<sup>α</sup>*, *θ*) representing the leading term of the *k* = 10 "line". It is important to note that all mentioned angular coefficients represent the AFCs *ψ<sup>k</sup>*,*<sup>p</sup>* with the maximum possible *p* for a given *k*.

#### **2. Derivation of the Angular Fock Coefficient** *ψ***5,2(***<sup>α</sup>***,** *θ***)**

The FRR (3) and (4) for *k* = 5 and *p* = 2 reduce to the form

$$
\left(\Lambda^2 - 45\right)\psi\_{5,2}(\mathfrak{a}, \theta) = h\_{5,2}(\mathfrak{a}, \theta),
\tag{14}
$$

where

$$h\_{5,2}(\mathfrak{a}, \theta) = -2V\psi\_{4,2}(\mathfrak{a}, \theta). \tag{15}$$

Using Equations (5), (6), (12) and (13), it is convenient to represent the RHS of Equation (14) in the form

$$h\_{5,2}(\mathbf{a}, \theta) = -\frac{Z^2(\pi - 2)(5\pi - 14)}{270\sqrt{\pi}} \left[ 3\pi^{-3/2} (2h\_1 + h\_2) - 2Z(h\_3 + \sqrt{2}h\_4) \right],\tag{16}$$

where

$$h\_1 = \frac{(1 - \tilde{\xi}^2)^2}{\tilde{\xi}}, \quad h\_2 = \frac{\cos(2a)}{\tilde{\xi}}, \quad h\_3 = \frac{\eta \,\chi\_{40}(a, \theta)}{\sin a}, \quad h\_4 = \frac{\eta \,\chi\_{42}(a, \theta)}{\sin a}.\tag{17}$$

Accordingly, we obtain the solution of Equation (14) in the identical form

$$\psi\_{5,2}(\alpha,\theta) = -\frac{Z^2(\pi - 2)(5\pi - 14)}{270\sqrt{\pi}} \left[ 3\pi^{-3/2}(2f\_1 + f\_2) - 2Z(f\_3 + \sqrt{2}f\_4) \right],\tag{18}$$

where the AFC-components *fi* satisfy the individual Fock recurrence relations (IFRRs)

$$\left(\Lambda^2 - 45\right)f\_i = h\_i. \qquad \left(i = 1, 2, 3, 4\right) \tag{19}$$

We sequentially find solutions to each of the IFRRs (19) using various methods presented in Ref. [4].

*2.1. Solution of the IFRR* (Λ<sup>2</sup> − <sup>45</sup>)*f*3 = *ηY*40/ sin *α*

> Moving from simpler to more complex solutions, let us start with IFRR

$$\left(\Lambda^2 - 45\right)f\_3 = h\_3. \tag{20}$$

The RHS *h*3 ≡ *h*3(*α*) represents the function of only one angle variable *α*. It was shown [4] that the solution of the corresponding IFRR (20) reduces to the solution *g*(*ρ*) = *f*3(*α*) of the inhomogeneous differential equation

$$\frac{1}{2}(\rho^2+1)^2\mathbb{g}''(\rho) + 2\rho^{-1}(\rho^2+1)\mathbb{g}'(\rho) + 4\mathbb{5}\mathbb{g}(\rho) = -\mathbf{h}(\rho),\tag{21}$$

where *ρ* = tan(*α*/2), and

$$h(\rho) \equiv h\_3(a) = \frac{(4\cos^2 a - 1)\sqrt{1 + \sin a}}{\pi^{3/2} \sin a} = \frac{(\rho + 1)(3\rho^4 - 10\rho^2 + 3)}{2\pi^{3/2} \rho(\rho^2 + 1)^{3/2}}.\tag{22}$$

For convenience, we solve the Equation (21) with the RHS <sup>h</sup>(*ρ*) not containing the multiplier *π* −3/2. The final solution *f*3 will be multiplied by this factor.

Using the method of variation of parameters, one obtains [4] the particular solution of Equation (21) in the form

$$\log^{(p)}(\rho) = v \varepsilon 0(\rho) \int \frac{u\_{50}(\rho) \mathbf{h}(\rho) d\rho}{(\rho^2 + 1)^2 \mathcal{W}\_0(\rho)} - u \varepsilon\_{50}(\rho) \int \frac{v\_{50}(\rho) \mathbf{h}(\rho) d\rho}{(\rho^2 + 1)^2 \mathcal{W}\_0(\rho)}\tag{23}$$

where

$$\mathcal{W}\_0(\rho) = -(\rho^2 + 1)/\rho^2. \tag{24}$$

The independent solutions of the homogeneous equation associated with Equation (21) are [4]:

$$\mathfrak{u}\_{50}(\rho) = \frac{(\rho^2 + 1)^{9/2}}{\rho} \,\_2F\_1\left( 4, \frac{7}{2}; \frac{1}{2}; -\rho^2 \right) = \frac{1 - 7\rho^2(3 - 5\rho^2 + \rho^4)}{\rho(\rho^2 + 1)^{5/2}},\tag{25}$$

$$w\_{50}(\rho) = (\rho^2 + 1)^{9/2} \,\_2F\_1\left(4, \frac{9}{2}; \frac{3}{2}; -\rho^2\right) = \frac{1 - 35\rho^2 + 21\rho^4 - \rho^6}{7(\rho^2 + 1)^{5/2}},\tag{26}$$

where <sup>2</sup>*F*1(...) is the Gaussian hypergeometric function. The substitution of Equations (24)–(26) into the general representation (23) yields

$$g^{(p)}(\rho) = \frac{-7\rho\{\rho[5\rho(3\rho - 4)(3\rho + 5) - 24] + 23\} - 23}{420\rho(\rho^2 + 1)^{5/2}}.\tag{27}$$

The general solution of the inhomogeneous equation can be expressed as the sum of the general solution of the associated homogeneous (complementary) equation and the particular solution of the inhomogeneous equation, whence

$$\lg(\rho) = \lg^{(p)}(\rho) + c\_{\text{\tiny $\mu$ }}\mu\_{\text{\tiny $0}}(\rho) + c\_{\text{\tiny$ v $}}v\_{\text{\tiny$ 0}}(\rho),\tag{28}$$

where the coefficients *cu* and *cv* are currently undetermined. To choose these coefficients, it is necessary to determine the behavior of all independent solutions on the boundaries of the domain [0, <sup>∞</sup>]. We easily obtain:

$$\mathcal{g}^{(p)}(\rho) \underset{\rho \to 0}{=} -\frac{23}{420\rho} - \frac{23}{60} + \frac{451\rho}{840} + O(\rho^2),\tag{29}$$

$$\log^{(p)}(\rho) \underset{\rho \to \infty}{=} -\frac{3}{4\rho} - \frac{1}{4\rho^2} + \frac{85}{24\rho^3} + O(\rho^{-4}),\tag{30}$$

$$
\mu\_{50}(\rho) \underset{\rho \to 0}{=} \frac{1}{\rho} - \frac{47\rho}{2} + \mathcal{O}(\rho^3),
\tag{31}
$$

$$\mu\_{50}(\rho) \underset{\rho \to \infty}{=} -7 + \frac{105}{2\rho^2} + O(\rho^{-4}),\tag{32}$$

$$w\_{50}(\rho) \underset{\rho \to 0}{=} 1 - \frac{15\rho^2}{2} + O(\rho^4),\tag{33}$$

$$
v\_{50}(\rho) \underset{\rho \to \infty}{=} -\frac{\rho}{7} + \frac{47}{14\rho} + O(\rho^{-3}).\tag{34}$$

It is seen that the particular solution *g*(*p*)(*ρ*) is divergent at *ρ* = 0, whereas the solutions of the homogeneous equation associated with Equation (21) are divergent, at *ρ* = 0 and *ρ* = ∞ for *<sup>u</sup>*50(*ρ*) and *<sup>v</sup>*50(*ρ*), respectively. Thus, to avoid the divergence on the whole range of definition, one should set *cu* = 23/420 and *cv* = 0 in the general solution (28). Then, the final physical solution becomes

$$f\mathfrak{z}(a) = -\frac{(\rho+1)(23\rho^4 + 22\rho^3 - 122\rho^2 + 22\rho + 23)}{60\pi^{3/2}(\rho^2 + 1)^{5/2}} = $$

$$= -\frac{1}{60\pi^{3/2}}[11\sin a + 21\cos(2a) + 2]\sqrt{1 + \sin a}.\tag{35}$$

*2.2. Solution of the IFRR* (Λ<sup>2</sup> − <sup>45</sup>)*f*4 = *ηY*42/*sinα*

> It was shown in Ref. [4] that the solution *f*4 ≡ *f*4(*<sup>α</sup>*, *θ*) of the IFRR

$$\left(\Lambda^2 - 45\right)f\_4 = 2\pi^{-3/2}\eta\sqrt{2}\sin a P\_2(\cos \theta)\tag{36}$$

can be found in the form

$$f\_4 = 2\sqrt{2} \text{ } \pi^{-3/2} \sin^2 \alpha P\_2(\cos \theta) g\_4(\rho), \tag{37}$$

where the function *g*4(*ρ*) satisfies the equation

$$(\rho^2 + 1)^2 \mathbb{g}\_4^{\prime\prime}(\rho) + 2\rho^{-1}[1 + \rho^2 + 2(1 - \rho^4)]\mathbb{g}\_4^{\prime}(\rho) + 13\mathbb{g}\_4(\rho) = -\mathbf{h}\_4(\rho),\tag{38}$$

with

$$\mathbf{h}\_4(\rho) = \frac{\eta}{\sin \kappa} = \frac{(\rho + 1)\sqrt{\rho^2 + 1}}{2\rho}. \tag{39}$$

Using the method of variation of parameters, one obtains [4] the particular solution of Equation (38) in the form

$$\chi\_{4}^{(p)}(\rho) = v\_{52}(\rho) \int \frac{u\_{52}(\rho) \mathbf{h}\_{4}(\rho) d\rho}{(\rho^{2} + 1)^{2} \mathcal{W}\_{2}(\rho)} - u\_{52}(\rho) \int \frac{v\_{52}(\rho) \mathbf{h}\_{4}(\rho) d\rho}{(\rho^{2} + 1)^{2} \mathcal{W}\_{2}(\rho)},\tag{40}$$

where

$$\mathcal{W}\_2(\rho) = -\frac{5}{\rho} \left( \frac{\rho^2 + 1}{\rho} \right)^5. \tag{41}$$

The independent solutions of the homogeneous equation associated with Equation (38) are:

$$u\_{\mathbb{S}2}(\rho) = \frac{(\rho^2 + 1)^{13/2}}{\rho^5} \,\_2F\_1\left(4, \frac{3}{2}; -\frac{3}{2}; -\rho^2\right) = \frac{1 + 11\rho^2 + 99\rho^4 - 231\rho^6}{\rho^5\sqrt{\rho^2 + 1}},\tag{42}$$

$$w\_{52}(\rho) = (\rho^2 + 1)^{13/2} \,\_2F\_1\left(4, \frac{13}{2}; \frac{7}{2}; -\rho^2\right) = \frac{231 - 99\rho^2 - 11\rho^4 - \rho^6}{231\sqrt{\rho^2 + 1}}.\tag{43}$$

Thus, the particular solution (40) reduces to the form:

$$\mathcal{g}\_4^{(p)}(\rho) = -\frac{11(21\rho^5 + 9\rho^4 + \rho^2) + 1}{2772\rho^5\sqrt{\rho^2 + 1}}.\tag{44}$$

Considering the series expansions for *g*(*p*) 4 (*ρ*), *<sup>u</sup>*52(*ρ*) and *<sup>v</sup>*52(*ρ*) on the boundaries of the range of definition (*ρ* ∈ [0, ∞]), it can be shown that function

$$g\_4^{(p)}(\\\rho) + \frac{1}{2772}\mu\_{52}(\rho) = -\frac{\rho + 1}{12\sqrt{\rho^2 + 1}} = -\frac{1}{12}\sqrt{1 + \sin a} \tag{45}$$

represents the physical (finite) solution of Equation (38), and hence we finally obtain:

$$f\_4(u, \theta) = -\frac{\sqrt{2}}{6\pi^{3/2}} (\sin u)^2 \sqrt{1 + \sin u} \,\, P\_2(\cos \theta). \tag{46}$$

*2.3. Solution of the IFRR* (Λ<sup>2</sup> − <sup>45</sup>)*f*1 = (1 − *ξ*<sup>2</sup>)2/*ξ*

It is important to note that the RHS *h*1 (see Equation (17)) of the IFRR

$$\left(\Lambda^2 - 45\right)f\_1 = h\_1(\xi)\tag{47}$$

is a function of *ξ* (only) defined by Equation (6). For this case [4,12], the solution of Equation (47) coincides with the solution of the inhomogeneous differential equation

$$(\mathfrak{J}^2 - 2)f\_1^{\prime\prime}(\mathfrak{J}) + \mathfrak{J}^{-1}(5\mathfrak{J}^2 - 4)f\_1^{\prime}(\mathfrak{J}) - 45f\_1(\mathfrak{J}) = h\_1(\mathfrak{J}).\tag{48}$$

A particular solution of Equation (48) can be found by the method of variation of parameters in the form [12]

$$f\_1^{(p)}(\\\\\xi) = \frac{1}{7\sqrt{2}} \left[ \mathfrak{u}\_5(\xi) \int \upsilon\_5(\xi) w(\xi) d\xi - \upsilon\_5(\xi) \int \mathfrak{u}\_5(\xi) w(\xi) d\xi \right],\tag{49}$$

where

$$w(\xi) = h\_1(\xi)\xi^2\sqrt{2-\xi^2}.\tag{50}$$

The linearly independent solutions of the homogeneous equation associated with Equation (48) are defined by the relations

$$u\_5(\xi) = \frac{\mathbf{P}\_{13/2}^{1/2}(\xi/\sqrt{2})}{\xi \sqrt[4]{2-\xi^2}} = \frac{2^{1/4}(8\xi^6 - 28\xi^4 + 28\xi^2 - 7)}{\sqrt{\pi(2-\xi^2)}},\tag{51}$$

$$w\_5(\xi) = \frac{Q\_{13/2}^{1/2}(\xi/\sqrt{2})}{\xi \sqrt[3]{2-\xi^2}} = \frac{\sqrt{\pi}(-8\xi^6 + 20\xi^4 - 12\xi^2 + 1)}{2^{3/4}\xi}.\tag{52}$$

where *Pμν* (*x*) and *Qμν* (*x*) are the associated Legendre functions of the first and second kind, respectively. A substitution of the representations (50)–(52) into (49) yields the particular solution:

$$f\_1^{(p)}(\\\\\xi) = -\frac{1}{60}\xi(13\xi^4 - 30\xi^2 + 15). \tag{53}$$

It can be verified that the particular solution *f* (*p*) 1 (*ξ*) is finite on the whole range of definition (*ξ* ∈ [0, √2]), whereas the solutions of the homogeneous equation associated with Equation (48) are divergent at *ξ* = √2 (*α* = *π*/2, *θ* = *π*) and *ξ* = 0 (*α* = *π*/2, *θ* = 0) for *<sup>u</sup>*5(*ξ*) and *<sup>v</sup>*5(*ξ*), respectively. Thus, we can conclude that the final physical solution of the IFFR (47) coincides with the particular solution (53), whence

$$f\_1 = -\frac{1}{60}\overline{\xi}(13\xi^4 - 30\overline{\xi}^2 + 15) = -\frac{1}{60}\sqrt{1 - \sin a \cos \theta} [\sin a \cos \theta (4 + 13\sin a \cos \theta) - 2]. \tag{54}$$

*2.4. Solution of the IFRR* (Λ<sup>2</sup> − <sup>45</sup>)*f*2 = *cos*(<sup>2</sup>*α*)/*ξ*

> To solve the IFRR

$$\left(\Lambda^2 - 45\right)f\_2 = h\_{2\prime} \tag{55}$$

with the RHS *h*2 defined by Equation (17), first of all, it is necessary to recall Sack's representation [13] (see also [3,4]) for *ξν* with *ν* = −1:

$$\xi^{r-1} = \sum\_{l=0}^{\infty} P\_l(\cos \theta) \left(\frac{\sin \alpha}{2}\right)^l F\_l(\rho), \tag{56}$$

where

$$F\_l(\rho) = \,\_2F\_l\left(\frac{l}{2} + \frac{1}{4}, \frac{l}{2} + \frac{3}{4}; l + \frac{3}{2}; \frac{4\rho^2}{(\rho^2 + 1)^2}\right) = \begin{cases} \; \_l\mathfrak{F}\_l(\rho) & 0 \le \rho \le 1\\ \; \_l\mathfrak{F}\_l(1/\rho) & \rho \ge 1 \end{cases} \tag{57}$$

with

$$\mathfrak{F}\_l(\rho) = (\rho^2 + 1)^{l + \frac{1}{2}}.\tag{58}$$

This enables us to present the RHS of Equation (55) in the form

$$h\_2 \equiv \frac{\cos(2a)}{\overline{\xi}} = \sum\_{l=0}^{\infty} P\_l(\cos\theta) (\sin\alpha)^l \mathbf{h}\_l(\rho),\tag{59}$$

where

$$\mathbf{h}\_l(\rho) = 2^{-l} F\_l(\rho) \cos(2a) = 2^{-l} F\_l(\rho) \left[ 1 - \frac{8\rho^2}{(\rho^2 + 1)^2} \right]. \tag{60}$$

In turn, it was shown in Ref. [4] that in the case where the RHS is determined by Equation (59), the solution of the corresponding IFRR (55) can be found in the form

$$f\_2(a,\theta) = \sum\_{l=0}^{\infty} P\_l(\cos\theta) (\sin a)^l \sigma\_l(\rho),\tag{61}$$

where the function *<sup>σ</sup>l*(*ρ*) satisfies the inhomogeneous differential equation

$$[(\rho^2+1)^2\sigma\_l''(\rho)+2\rho^{-1}[1+\rho^2+l(1-\rho^4)]\sigma\_l'(\rho)+(5-2l)(9+2l)\sigma\_l(\rho)=-\mathbf{h}\_l(\rho). \tag{62}$$

The linearly independent solutions of the homogeneous equation associated with Equation (62) are:

$$\begin{split} u\_{5l}(\rho) &= \frac{(\rho^2+1)^{l+9/2}}{\rho^{2l+1}} \,\_2F\_1\left(4, \frac{7}{2}-l; \frac{1}{2}-l; -\rho^2\right) = \frac{(\rho^2+1)^{l-5/2}}{\rho^{2l+1}} \times \\ \left[\rho^6 \left(1+\frac{120}{2l-5}-\frac{120}{2l-3}+\frac{24}{2l-1}\right)+3\rho^4 \left(1+\frac{40}{2l-3}-\frac{24}{2l-1}\right)+3\rho^2 \left(1+\frac{8}{2l-1}\right)+1\right], \end{split} \tag{63}$$

$$v\_{5l}(\rho) = (\rho^2 + 1)^{l+9/2} \,\_2F\_1\left(4, \frac{9}{2} + l; \frac{3}{2} + l; -\rho^2\right) = (\rho^2 + 1)^{l-5/2} \times$$

$$\left[\rho^6 \left(1 - \frac{24}{2l+3} + \frac{120}{2l+5} - \frac{120}{2l+7}\right) + 3\rho^4 \left(1 + \frac{24}{2l+3} - \frac{40}{2l+5}\right) + 3\rho^2 \left(1 - \frac{8}{2l+3}\right) + 1\right]. \tag{64}$$

The method of variation of parameters enables us to obtain the particular solution of the inhomogeneous differential Equation (62) in the form

$$\sigma\_{l}^{(p)}(\rho) = v\_{5l}(\rho) \int \frac{u\_{5l}(\rho)\mathbf{h}\_{l}(\rho)d\rho}{(\rho^{2} + 1)^{2}\mathcal{W}\_{l}(\rho)} - u\_{5l}(\rho) \int \frac{v\_{5l}(\rho)\mathbf{h}\_{l}(\rho)d\rho}{(\rho^{2} + 1)^{2}\mathcal{W}\_{l}(\rho)}\tag{65}$$

where

$$\mathcal{W}\_l(\rho) = -\frac{2l+1}{\rho} \left( \frac{\rho^2 + 1}{\rho} \right)^{2l+1}. \tag{66}$$

Note that due to different representations for the function *Fl*(*ρ*) (see Equation (57)) at values of *ρ* less and greater than 1, we obtain special representations for a particular solution in these two regions:

$$\sigma\_l^{(0)}(\rho) = \frac{(\rho^2 + 1)^{l-5/2}}{2^{l+1}(2l-1)(2l-3)} \left[ (2l-3)\rho^4 - 4(l-2)\rho^2 - \frac{4l^2 + 4l - 27}{3(2l-5)} \right], \quad 0 \le \rho \le 1 \tag{67}$$

$$\sigma\_{l}^{(1)}(\rho) = \frac{\rho^{-2l-1}(\rho^2+1)^{l-5/2}}{2^{l+1}(2l+3)(2l+5)} \left[ \frac{4l^2+4l-27}{3(2l+7)} + 4(l+3)\rho^2 - (2l+5)\rho^4 \right]. \qquad \rho \ge 1 \quad \text{(68)}$$

It can be verified that both functions (67) and (68) have no singularities on their domains of definition. On the other hand, function *<sup>u</sup>*5*l*(*ρ*) is singular at *ρ* = 0, whereas *<sup>v</sup>*5*l*(*ρ*) is singular at *ρ* = ∞. This means that one should search the general solution of Equation (62) in the form:

$$
\sigma\_l(\rho) = \sigma\_l^{(0)}(\rho) + c\_{5l}^{(v)} v\_{5l}(\rho), \qquad \qquad 0 \le \rho \le 1 \tag{69}
$$

$$
\sigma\_l(\rho) = \sigma\_l^{(1)}(\rho) + c\_{\mathbb{S}l}^{(u)} u\_{\mathbb{S}l}(\rho). \qquad \qquad \rho \ge 1 \tag{70}
$$

Note that two coefficients *c*(*v*) 5*l* and *c*(*u*) 5*l* are presently undetermined. To calculate them, we need to find two equations relating these coefficients. The first equation is quite obvious. It follows from the condition that the representations (69) and (70) are coincident at the common point *ρ* = 1, that is

$$
\sigma\_l^{(0)}(1) + c\_{\mathfrak{S}l}^{(v)} v\_{\mathfrak{S}l}(1) = \sigma\_l^{(1)}(1) + c\_{\mathfrak{S}l}^{(u)} u\_{\mathfrak{S}l}(1). \tag{71}
$$

This relationship reduces to the first desired equation:

$$c\_{5l}^{(u)}(2l+3)(2l+5)(2l+7) = c\_{5l}^{(v)}(2l-5)(2l-3)(2l-1) + \frac{27-4l-4l^2}{3 \times 2^{l+1}}.\tag{72}$$

It can be verified that <sup>h</sup>*l*(0) = <sup>h</sup>*l*(∞) = 2−*l*. It follows from these relations that *σ*(0) *l* (0) + *c* (*v*) 5*l <sup>v</sup>*5*l*(0) = *σ*(1) *l* (∞) + *c* (*u*) 5*l <sup>u</sup>*5*<sup>l</sup>*(∞). It can be assumed that the last equation represents the second desired equation. However, this assumption turns out to be false, because it again leads to Equation (72).

We propose the following method to find the second desired equation. Recall that any suitable function of the angles *α* and *θ* may be expanded into HHs since they form a complete set:

$$f(\mathfrak{a}, \theta) = \sum\_{n=0(2)}^{\infty} \sum\_{l=0}^{n/2} f\_{n,l} \chi\_{n,l}(\mathfrak{a}, \theta)\_{l} \tag{73}$$

where (see, e.g., [3])

$$f\_{n,l} = \int f(\alpha, \theta) \, \chi\_{n,l}(\alpha, \theta) d\Omega \tag{74}$$

with

$$d\Omega = \pi^2 \sin^2 \alpha \sin \theta d\alpha d\theta. \qquad \text{or} \quad [0, \pi], \ \theta \in [0, \pi] \tag{75}$$

For the function *f*(*<sup>α</sup>*, *θ*) = *f*2(*<sup>α</sup>*, *θ*) represented by Equation (61), the expansion coefficient with *n* = 2*l* becomes

$$f\_{2\parallel l} = \pi^2 \int\_0^\pi \int\_0^\pi f\_2(a,\theta) \mathbf{Y}\_{2l,l}(a,\theta) \sin^2 a \sin\theta da d\theta = \frac{2\pi^2 N\_{2l,l}}{2l+1} \left[ \mathbf{K}\_0(l) + \mathbf{c}\_l^{(v)} \mathbf{K}\_0(l) + \mathbf{c}\_l^{(u)} \mathbf{K}\_0(l) \right],\tag{76}$$

where

$$\begin{split} K\_0(l) &= \int\_0^{\pi/2} (\sin a)^{2l+2} \sigma\_l^{(0)}(\rho) da + \int\_{\pi/2}^{\pi} (\sin a)^{2l+2} \sigma\_l^{(\infty)}(\rho) da = \\ &= \frac{\sqrt{2} (-24l^3 - 100l^2 + 198l + 249)}{(2l-5)(2l-3)(2l+3)(2l+5)(2l+7)(2l+9)}, \end{split} \tag{77}$$

$$K\_{\mathbb{P}}(l) = \int\_0^{\pi/2} (\sin a)^{2l+2} v\_{\mathbb{S}l}(\rho) d\mathfrak{a} = \frac{2^{l+3/2}(2l-1)}{(2l+5)(2l+9)'} \tag{78}$$

$$K\_{\mathfrak{u}}(l) = \int\_{\pi/2}^{\pi} (\sin \mathfrak{a})^{2l+2} u\_{5l}(\rho) d\mathfrak{a} = \frac{2^{l+3/2} (2l+3)(2l+7)}{(2l-5)(2l-3)(2l+9)}.\tag{79}$$

To derive the results (76)–(79), we use the representation *<sup>Y</sup>*2*l*,*<sup>l</sup>*(*<sup>α</sup>*, *θ*) = *<sup>N</sup>*2*l*,*<sup>l</sup>* sin*<sup>l</sup> αPl*(cos *θ*) for the particular case of the HHs, and the orthogonality property for the Legendre polynomials. It should be noted that the explicit form of the normalization constant *<sup>N</sup>*2*l*,*<sup>l</sup>* is not required.

On the other hand, expanding *f*2(*<sup>α</sup>*, *θ*) in HHs, and inserting this expansion into the LHS of the IFRR (55), we obtain

$$f\_1\left(\Lambda^2 - 45\right)f\_2(a,\theta) = \sum\_{n=0(2)}^{\infty} \sum\_{l=0}^{n/2} f\_{n,l} [n(n+4) - 45] Y\_{n,l}(a,\theta). \tag{80}$$

To derive the last equation, we use the fact that *Yn*,*<sup>l</sup>*(*<sup>α</sup>*, *θ*) is an eigenfunction of the operator Λ<sup>2</sup> with an eigenvalue equal to *n*(*n* + <sup>4</sup>), that is

$$
\Lambda^2 Y\_{n,l}(a,\theta) = n(n+4)Y\_{n,l}(a,\theta). \tag{81}
$$

The HH expansion of the RHS of Equation (55) is

$$h\_2(\mathfrak{a}, \theta) = \sum\_{n=0(2)}^{\infty} \sum\_{l=0}^{n/2} h\_{n,l} \chi\_{n,l}(\mathfrak{a}, \theta). \tag{82}$$

Hence,

$$f\_{2l,l} = \frac{h\_{2l,l}}{4l(l+2) - 45}.\tag{83}$$

Using again Sack's representation (56) and (57) and Equation (74), we obtain the expansion coefficient *<sup>h</sup>*2*l*,*<sup>l</sup>* in explicit form:

$$h\_{2l\downarrow} = \pi^2 \int\_0^\pi \int\_0^\pi h\_2(a,\theta) \mathbb{1}\_{2l\downarrow}(a,\theta) \sin^2 a \sin\theta da d\theta = -\frac{2^{7/2}\pi^2 N\_{2l\downarrow}(4l^2 + 24l + 19)}{(2l+1)(2l+3)(2l+5)(2l+7)}.\tag{84}$$

Thus, inserting (84) into the RHS of Equation (83) and equating the result to the RHS of Equation (76) we obtain the desired *second equation* in the form:

$$c\_{5l}^{(u)}(2l+3)(2l+5)(2l+7) = -c\_{5l}^{(v)}(2l-5)(2l-3)(2l-1) - 2^{-l-1}(2l+1).\tag{85}$$

Solving the system of two linear Equations (72) and (85) gives the desired coefficients:

$$\mathfrak{c}\_{5l}^{(u)} = -\frac{2^{-l-1}(l+4)(2l-3)}{3(2l+3)(2l+5)(2l+7)}, \quad \mathfrak{c}\_{5l}^{(v)} = \frac{2^{-l-1}(l-3)(2l+5)}{3(2l-5)(2l-3)(2l-1)}.\tag{86}$$

It should be noted that the method described above for calculating the coefficients *c*(*u*) 5*l* and *c*(*v*) 5*l* is very reliable, but quite complex. A much simpler method is based on the statement that the point *ρ* = 1 represents the match point for the functions defined by Equations (69) and (70). This means that not only these functions, but also their first (at least) derivatives must coincide at this point. Thus, the second required equation relating the coefficients *c*(*u*) 5*l* and *c*(*v*) 5*l*is:

$$\frac{d\sigma\_l^{(0)}(\rho)}{d\rho}\Big|\_{\rho=1} + c\_{5l}^{(v)}\frac{d\upsilon\_{5l}(\rho)}{d\rho}\Big|\_{\rho=1} = \frac{d\sigma\_l^{(1)}(\rho)}{d\rho}\Big|\_{\rho=1} + c\_{5l}^{(u)}\frac{d u\_{5l}(\rho)}{d\rho}\Big|\_{\rho=1}.\tag{87}$$

The solution of the system of two Equations (72) and (87) again gives the coefficients defined by Equation (86). Substituting these coefficients into the representations (69) and (70), we finally obtain:

$$f\_{\mathcal{Z}}(\alpha,\theta) = \frac{1}{6} \sum\_{l=0}^{\infty} \frac{\zeta\_l(\rho) P\_l(\cos\theta)}{(2l-1)(2l+3)},\tag{88}$$

where 
$$\mathcal{Z}\_{l}(\rho) = \begin{cases} \begin{array}{c} \chi\_{l}(\rho), \\ \chi\_{l}(1/\rho), \end{array} & 0 \le \rho \le 1 \\ \end{array} \tag{89}$$

with

$$\chi\_l(\rho) = \frac{\rho^l}{(\rho^2 + 1)^{5/2}} \left[ \frac{(l-3)(2l-1)\rho^6}{2l+7} + 9l\rho^4 - 9(l+1)\rho^2 - \frac{(l+4)(2l+3)}{2l-5} \right]. \tag{90}$$

It is clear that only the function (69) is required for calculating the function *<sup>χ</sup>l*(*ρ*). Thus, in fact, we need to calculate only one coefficient *c*(*v*) 5*l* to define this function. In this regard, it is important to emphasize that the representation (89) reflects the fact that the WF of a two-electron atomic system must preserve its parity when interchanging the electrons. For the singlet S-states (which include the ground state) this means that the AFC and/or its

component preserves its form (including the sign) under the transformation *α π* − *α*. For the AFC-component *f*2(*<sup>α</sup>*, *<sup>θ</sup>*), represented by the series expansion (61), this property corresponds (in terms of variable *ρ*) to the relationship:

$$
\sigma\_l^{(0)}(\rho^{-1}) + c\_{5l}^{(v)} v\_{5l}(\rho^{-1}) = \sigma\_l^{(1)}(\rho) + c\_{5l}^{(u)} u\_{5l}(\rho). \tag{91}
$$

The elimination of the RHSs between Equations (71) and (91) for *ρ* = 1 yields the identity, whereas the use of Equation (87) instead of Equation (71) yields the required equation:

$$\frac{d\sigma\_l^{(0)}(\rho)}{d\rho}\Big|\_{\rho=1} + c\_{5l}^{(v)}\frac{d\upsilon\_{\overline{l}l}(\rho)}{d\rho}\Big|\_{\rho=1} = \frac{d\sigma\_l^{(0)}(\rho^{-1})}{d\rho}\Big|\_{\rho=1} + c\_{5l}^{(v)}\frac{d\upsilon\_{\overline{l}l}(\rho^{-1})}{d\rho}\Big|\_{\rho=1}.\tag{92}$$

The solution of the last equation gives the coefficients *c*(*v*) 5*l* presented by Equation (86). Note that the coefficient *c*(*u*) 5*l*can then be calculated by the use of Equation (71) if needed.

 In the general case, we cannot sum the infinite series (88) to obtain the function *f*2(*<sup>α</sup>*, *θ*) in an explicit closed form. However, this can be done for some special angles *α* and/or *θ*. For example, it is worth noting that the angles *θ* = 0, *π* correspond to the collinear configuration [14] of the two-electron atomic system in question. For these cases we obtain

$$f\_2(\mathfrak{a},0) = \pm \frac{(\rho - 1)(12\rho^4 - 13\rho^3 - 88\rho^2 - 13\rho + 12)}{90(\rho^2 + 1)^{5/2}},\tag{93}$$

$$f\_2(\mathfrak{a}, \pi) = -\frac{(\rho + 1)(12\rho^4 + 13\rho^3 - 88\rho^2 + 13\rho + 12)}{90(\rho^2 + 1)^{5/2}}.\tag{94}$$

Sign "+" in Equation (93) corresponds to 0 ≤ *α* ≤ *π*/2 (0 ≤ *ρ* ≤ 1), whereas "−" to *π*/2 ≤ *α* ≤ *π* (*ρ* ≥ 1). The list of special *θ*-angles can be supplemented with an intermediate angle *θ* = *π*/2:

$$f\_2\left(\alpha, \frac{\pi}{2}\right) = -\frac{2(\rho^4 - 3\rho^2 + 1)}{15(\rho^2 + 1)^2}.\tag{95}$$

It is worth noting that for the important cases of the nucleus–electron and electron–electron coalescence, Equation (88), respectively, reduces to:

$$f\_2(0, \theta) = -\frac{2}{15}, \qquad f\_2\left(\frac{\pi}{2}, 0\right) = 0. \tag{96}$$

To derive the results (93)–(96) we used the relationships:

$$P\_{\rm tr}(0) = \sqrt{\pi} \Gamma^{-1} \left( \frac{1-n}{2} \right) \Gamma^{-1} \left( \frac{n}{2} + 1 \right), \quad P\_{\rm tr}(1) = 1, \quad P\_{\rm tr}(-1) = (-1)^n,\tag{97}$$

where <sup>Γ</sup>(*x*) is the gamma function.

#### **3. Derivation of the Angular Fock Coefficient** *ψ***6,3(***<sup>α</sup>***,** *θ***)**

We start this section by considering the FRR (3) and (4) for *k* = 6 and *p* = 2:

$$
\left(\Lambda^2 - 60\right)\psi\_{6,2} = 48\psi\_{6,3} - 2V\psi\_{5,2} + 2E\psi\_{4,2}.\tag{98}
$$

Next, let us expand each function in Equation (98) into HHs, using Equation (73). This gives

$$\psi\_{k,p} = \sum\_{n=0(2)}^{\infty} \sum\_{l=0}^{n/2} \mathfrak{c}\_{nl}^{(kp)} Y\_{n,l}(\mathfrak{a}, \theta), \tag{99}$$

with {*k*, *p*} = {6, <sup>3</sup>}, {6, <sup>2</sup>}, {4, 2} and

$$V\Psi\_{5,2} = \sum\_{n=0(2)}^{\infty} \sum\_{l=0}^{n/2} \mathbf{f}\_{nl} \mathbf{Y}\_{n,l}(\alpha, \theta),\tag{100}$$

where the dimensionless potential *V* is defined by Equation (5), whereas the expansion coefficient f*nl* can be calculated by the formula

$$\mathbf{f}\_{\rm nl} = \int V \psi\_{\rm 5,2} \mathbf{Y}\_{\rm n,l}(\boldsymbol{\alpha}, \boldsymbol{\theta}) d\boldsymbol{\Omega},\tag{101}$$

according to Equation (74).

It follows from Equation (4) that *hk*,*k*/2 = 0 for even *k*. Using additionally Equation (81), we can conclude that the AFC *ψ<sup>k</sup>*,*k*/2 (with even *k*) represents the linear combination of the HHs, *Yk*,*<sup>l</sup>*(*<sup>α</sup>*, *<sup>θ</sup>*). Hence,

$$
\varepsilon\_{nl}^{(63)} = 0 \quad for \quad n \neq 6,\tag{102}
$$

$$
\varepsilon\_{nl}^{(42)} = 0 \quad for \quad n \neq 4. \tag{103}
$$

Equating the coefficients for the HHs, *<sup>Y</sup>*6,*<sup>l</sup>*(*<sup>α</sup>*, *θ*) in both sides of Equation (98), we obtain:

$$0 = 48c\_{\text{\{d\}}}^{\text{(63)}} - 2\mathbf{f}\_{\text{\{d\}}}.\tag{104}$$

Hence, (using additionally Equation (101)),

$$\mathbf{c}\_{6l}^{(63)} = \frac{1}{24} \int V \psi\_{5,2} \mathbf{Y}\_{6,l}(\mathbf{a}, \theta) d\Omega. \tag{105}$$

Note that the LHS of Equation (104) equals zero, because (Δ<sup>2</sup> − <sup>60</sup>)*<sup>Y</sup>*6,*l* = 0 as follows from Equation (81).

Thus, according to Equations (99) and (102), the AFC *ψ*6,3(*<sup>α</sup>*, *θ*) represents a linear combination of four HHs, *<sup>Y</sup>*6,*<sup>l</sup>*(*<sup>α</sup>*, *θ*) with *l* = 0, 1, 2, 3. The contribution of each HH is determined by the coefficient *c*(63) 6*l* given by Equation (105). However, it is easy to prove that only the coefficients with odd values of *l* are nonzero for *ψ*6,3(*<sup>α</sup>*, *<sup>θ</sup>*). Indeed, this has already been mentioned in Section 2.4 that the WF of a two-electron atom/ion must preserve its parity when interchanging the electrons. For the singlet S-states this means that only the HHs, which preserve the sign under transformation *α π* − *α*, differ from zero in the expansion of the WF, and hence in the expansion of any AFC. In turn, it is easy to show that only *Yn*,*<sup>l</sup>*(*<sup>α</sup>*, *θ*) with even values of (*n*/2 − *l*) satisfy the above property. Hence, the AFC in question, becomes

$$
\psi\_{6,3}(\mathfrak{a},\theta) = a\_{61}\chi\_{6,1}(\mathfrak{a},\theta) + a\_{63}\chi\_{6,3}(\mathfrak{a},\theta),
\tag{106}
$$

where we denoted *<sup>a</sup>*6*l* ≡ *c*(63) 6*l* (*l* = 1, 3) for convenience and simplicity, and where the normalized HHs are

$$Y\_{6,1}(\mathbf{a},\theta) = \frac{2\left[\sin\alpha + 3\sin(3\alpha)\right]\cos\theta}{\pi^{3/2}\sqrt{5}}, \quad Y\_{6,3}(\mathbf{a},\theta) = \frac{8\sin^3\alpha P\_3(\cos\theta)}{\pi^{3/2}\sqrt{5}}.\tag{107}$$

Using Formula (105) and taking into account the representations (5) and (18) for the dimensionless potential *V* and the AFC *ψ*5,2, respectively, we can represent the desired coefficients in the form

$$a\_{\delta l} = -\frac{(\pi - 2)(5\pi - 14)}{6480} \left( I\_{l,4} Z^4 + I\_{l,3} Z^3 + I\_{l,2} Z^2 \right), \quad (l = 1, 3) \tag{108}$$

where

$$I\_{l,4} = 4\pi^{3/2} \int\_0^\pi \int\_0^\pi \left[ f\_3(a,\theta) + \sqrt{2} f\_4(a,\theta) \right] \eta \,\chi\_{6,l}(a,\theta) \sin a \sin \theta da d\theta,\tag{109}$$

$$I\_{l,3} = -2\int\_0^\pi \int\_0^\pi \left\{ \frac{3\eta[2f\_1(a,\theta) + f\_2(a,\theta)]}{\sin a} + \frac{\pi^{3/2}\left[f\_3(a,\theta) + \sqrt{2}f\_4(a,\theta)\right]}{\xi} \right\} \times \int\_{\omega \times (\pi - 0) + \sin^2 \theta \sin \theta d\omega d\theta} \times \quad (11.10)$$

<sup>×</sup>*Y*6,*<sup>l</sup>*(*<sup>α</sup>*, *θ*) sin<sup>2</sup> *α* sin *θdαdθ*, (110)

$$I\_{l,2} = 3\int\_0^\pi \int\_0^\pi [2f\_1(a,\theta) + f\_2(a,\theta)]\xi^{-1} \mathcal{Y}\_{6,l}(a,\theta)\sin^2 a \sin\theta da d\theta. \tag{111}$$

To calculate the integral (110), it is useful to separate the contributions which include the functions *f*1, *f*3 and *f*4, represented by the explicit closed expressions, and the function *f*2, represented by the infinite series (88). We obtain:

$$I\_{l,3} = I\_{l,3}^{(134)} - \Theta I\_{l,3}^{(2)},\tag{112}$$

where

$$I\_{l,3}^{(134)} = -2\int\_0^\pi \int\_0^\pi \left\{ \frac{6\eta f\_1(a,\theta)}{\sin a} + \frac{\pi^{3/2} \left[ f\_3(a,\theta) + \sqrt{2} f\_4(a,\theta) \right]}{\xi} \right\} Y\_{\theta,l}(a,\theta) \sin^2 a \sin\theta da d\theta,\tag{113}$$

$$I\_{l,3}^{(2)} = \int\_0^{\pi} \int\_0^{\pi} f\_2(\mathfrak{a}, \theta) \eta \, \mathcal{Y}\_{6,l}(\mathfrak{a}, \theta) \sin \mathfrak{a} \sin \theta \, dad\theta. \tag{114}$$

The integrals (113) can be taken in an explicit (closed) form that gives:

$$I\_{1,3}^{(134)} = \frac{3(45\pi - 122)}{35\pi^{3/2}\sqrt{5}}, \qquad I\_{3,3}^{(134)} = \frac{245\pi - 816}{70\pi^{3/2}\sqrt{5}}.\tag{115}$$

The problem of calculating the integrals (114) is that the corresponding integrands contain the function *f*2(*<sup>α</sup>*, *θ*) represented by the infinite series (88). Fortunately, using the orthogonality relationship for the Legendre polynomials, we can ge<sup>t</sup> these integrals also in explicit form. Changing the order of summation and integration, we easily obtain:

$$I\_{1,3}^{(2)} = \frac{\pi^{-3/2}}{3\sqrt{5}} \sum\_{l=0}^{\infty} \int\_0^{\pi} \left[ \sin a + 3\sin(3a) \right] \eta \left[ \frac{\zeta\_l(\rho)}{(2l-1)(2l+3)} \right] \sin a da \int\_0^{\pi} P\_l(\cos \theta) \cos \theta \sin \theta d\theta$$

$$= \frac{2\pi^{-3/2}}{45\sqrt{5}} \int\_0^{\pi} \left[ \sin a + 3\sin(3a) \right] \eta \left[ \zeta\_1(\rho) \sin a da \right] = \frac{7\pi + 22}{210\pi^{3/2}\sqrt{5}},\tag{116}$$

$$I\_{3,3}^{(2)} = \frac{4\pi^{-3/2}}{3\sqrt{5}} \sum\_{l=0}^{\infty} \int\_0^{\pi} \eta \left[ \frac{\zeta\_l(\rho)}{(2l-1)(2l+3)} \right] \sin^4 a \, da \int\_0^{\pi} P\_l(\cos \theta) P\_l(\cos \theta) \sin \theta \, d\theta$$

$$= \frac{8\pi^{-3/2}}{945\sqrt{5}} \int\_0^{\pi} \eta \, \zeta\_3(\rho) \sin^4 a \, da = \frac{3\pi - 32}{180\pi^{3/2}\sqrt{5}}.\tag{117}$$

Recall that *η* ≡ *η*(*α*) is defined by Equation (6) and *ρ* = tan(*α*/2).

It can be shown (using fairly long nontrivial derivations) that the integrals *Il*,<sup>2</sup> and *Il*,<sup>4</sup> vanish both for *l* = 1 and *l* = 3. This means that (according to the representations (106) and (108)) the AFC, *ψ*6,3(*<sup>α</sup>*, *θ*) is proportional to the third power of the nucleus charge *Z* (only), which is in full agreemen<sup>t</sup> with Formula (13) from Ref. [4].

Thus, combining the results of this section, we obtain the nonzero coefficients *<sup>a</sup>*6,*l* in the simple final form:

$$a\_{61} = -\frac{(\pi - 2)(5\pi - 14)(32\pi - 97)}{56700\pi^{3/2}\sqrt{5}}Z^3,\tag{118}$$

$$n\_{63} = -\frac{(\pi - 2)(5\pi - 14)(357\pi - 1112)}{680400\pi^{3/2}\sqrt{5}}Z^3. \tag{119}$$

#### **4. Derivation of the Angular Fock Coefficients** *ψ***7,3(***<sup>α</sup>***,** *θ***) and** *ψ***8,4(***<sup>α</sup>***,** *θ***)**

In Sections 2 and 3 we detailed the derivation of the AFCs *ψ*5,2(*<sup>α</sup>*, *θ*) and *ψ*6,3(*<sup>α</sup>*, *<sup>θ</sup>*), respectively. Therefore, for the AFCs *ψ*7,3(*<sup>α</sup>*, *θ*) and *ψ*8,4(*<sup>α</sup>*, *<sup>θ</sup>*), we give only abbreviated derivations, and include extended explanations only in case of significant differences.

$$4.1. \text{ The } A\text{FC } \psi \gamma \beta (\alpha, \theta) \text{ }$$

The FRR (3) and (4) for *k* = 7 and *p* = 3 reduces to the form

$$\left(\Lambda^2 - \mathcal{T}\mathcal{T}\right)\psi\_{\mathcal{T},3}(\mathfrak{a},\theta) = h\_{\mathcal{T},3}(\mathfrak{a},\theta),\tag{120}$$

where

$$h\_{7,3}(\mathfrak{a},\theta) = -2V\psi\_{6,3}(\mathfrak{a},\theta). \tag{121}$$

Using Equations (106), (107) and (5) the RHS of Equation (120) can be represented in the form:

$$h\_{7,3}(\mathbf{u}, \theta) = \frac{(\pi - 2)(5\pi - 14)Z^3}{340200\sqrt{5}\pi^{3/2}} \left\{ \frac{h\_1 + h\_2}{\sqrt{5}\pi^{3/2}} - 2Z \left[ 12(32\pi - 97)\mathbf{\hat{h}}\_3 + (357\pi - 1112)\mathbf{\hat{h}}\_4 \right] \right\},\tag{122}$$

where

$$\bar{h}\_1 = 20\bar{\xi}^{-1} \left[ 12(32\pi - 97)(1 - \bar{\xi}^2) + (357\pi - 1112)(1 - \bar{\xi}^2)^3 \right],\tag{123}$$

$$
\bar{h}\_2 = 60(688 - 255\pi) \mathfrak{J}^{-1} \sin^3 \alpha \cos \theta,\tag{124}
$$

$$
\bar{h}\_3 = \eta(\sin \alpha)^{-1} \mathcal{Y}\_{6,1}(a, \theta), \quad \bar{h}\_4 = \eta(\sin \alpha)^{-1} \mathcal{Y}\_{6,3}(a, \theta). \tag{125}
$$

Accordingly, the solution of the FRR (120) can be found in the form:

$$\psi\_{7,3}(\mathbf{a},\theta) = \frac{(\pi - 2)(5\pi - 14)Z^3}{340200\sqrt{5}\pi^{3/2}} \left\{ \frac{\vec{f}\_1 + \vec{f}\_2}{\sqrt{5}\pi^{3/2}} - 2Z \left[12(32\pi - 97)\vec{f}\_3 + (357\pi - 1112)\vec{f}\_4\right] \right\},\tag{126}$$

where the AFC components ¯ *fi* satisfy the IFRRs

$$(\Lambda^2 - 77)\vec{f}\_i = \vec{h}\_i. \qquad (i = 1, 2, 3, 4) \tag{127}$$

Note that the components ¯ *hi* of the RHS *h*7,3 of the FRR (120) for the AFC *ψ*7,3 are reasonably close to the components *hi* of the RHS *h*5,2 of the FRR (14). Therefore, we only briefly dwell on the conclusions of the corresponding results, as we noted earlier.

¯

It is seen from Equation (123) that the RHS *h*1 is a function of a single variable *ξ* defined by Equation (6). The solution of the corresponding IFRR was described in Section IV of Ref. [4] (see also Section II of Ref. [12]) and illustrated (among others) in Section 2.3 of the current article. Thus, following the technique mentioned above, we obtain:

$$\begin{split} \vec{f}\_{1} &= \left(\frac{41437\pi}{12} - \frac{74342}{7}\right) \mathfrak{j}^{\uparrow} + \left(36476 - \frac{35588\pi}{3}\right) \mathfrak{j}^{\downarrow} + \\ &\quad + \frac{5}{2} (4931\pi - 15156) \mathfrak{j}^{3} + 5(2276 - 741\pi) \mathfrak{j}. \end{split} \tag{128}$$

It can be verified that the RHSs ¯ *h*3 and ¯ *h*4 represent functions of the form *f*(*α*)*Pl*(cos *θ*) with *l* equals 1 and 3, respectively. The solution of the corresponding IFRR was described in Section V of Ref. [4] and illustrated in Sections 2.1 and 2.2 of the current article. This enables us to obtain:

$$\bar{f}\_3 = -\frac{\rho(1+\rho)(29+\rho\{16+\rho[\rho(16+29\rho)-114]\})\cos\theta}{9\sqrt{5}\ \pi^{3/2}(\rho^2+1)^{7/2}},\tag{129}$$

$$f\_4 = -\frac{\sin^3 \pi \sqrt{1 + \sin \pi}}{2\sqrt{5} \,\pi^{3/2}} P\_3(\cos \theta). \tag{130}$$

Recall that the *ρ* variable was defined previously in Section 2.1.

The RHS ¯ *h*2 represented by Equation (124) is slightly more complicated than *h*2 discussed in Section 2.4. In this regard, it would be useful to clarify two points.

First, using representation (56) for *ξ*−1, we can rewrite Equation (124) in the form:

$$
\hbar\_2 = 60(688 - 225\pi)\hbar\_2\tag{131}
$$

where

$$\tilde{h}\_2 = \sum\_{l=0}^{\infty} 2^{-l} (\sin a)^{l+3} F\_l(\rho) \cos \theta P\_l(\cos \theta), \tag{132}$$

and where *Fl*(*ρ*) is defined by Equations (57) and (58). In order to apply the solution of the corresponding IFRR by the method described in Section 2.4 (see also [4]), the *θ*-dependent *l*-componen<sup>t</sup> in the series expansion of ˜ *h*2 must be pure *Pl*(cos *<sup>θ</sup>*). To solve the problem, one could use the general formula representing the Clebsch–Gordan series for a product of two *spherical* harmonics. However, in our simple case, it is easier to use the recurrence relation for the Legendre polynomials

$$(l+1)P\_{l+1}(\mathbf{x}) - (2l+1)\mathbf{x}P\_l(\mathbf{x}) + lP\_{l-1}(\mathbf{x}) = 0,\tag{133}$$

which enables us to represent ˜ *h*2 in the desired form:

$$\ddot{h}\_2 = \sum\_{l=0}^{\infty} \ddot{h}\_l(\rho) (\sin \alpha)^l P\_l(\cos \theta),\tag{134}$$

where

$$\bar{\mathbf{h}}\_{l}(\rho) = \frac{l}{2^{l-1}(2l-1)} \sin^{2} a F\_{l-1}(\rho) + \frac{l+1}{2^{l+1}(2l+3)} \sin^{4} a F\_{l+1}(\rho). \tag{135}$$

The second point is related to the calculation of the coefficient

$$\tilde{\mathbf{h}}\_{2l,l} = \pi^2 \int\_0^\pi \int\_0^\pi \tilde{h}\_2(\alpha, \theta) Y\_{2l,l}(\alpha, \theta) \sin^2 \alpha \sin \theta da d\theta \tag{136}$$

in the HH expansion of ˜ *h*2 (see the corresponding Equation (84) for calculation of *ψ*5,2(*<sup>α</sup>*, *θ*)). Of course, we can use representation (134) and (135) and then apply the orthogonality condition for the Legendre polynomials. However, the simpler way is to use the original representation (132) taking into account that cos *θ* ≡ *<sup>P</sup>*1(cos *<sup>θ</sup>*). In this case, we can apply the well-known formula for the integral of three Legendre polynomials

$$\int\_{-1}^{1} P\_l(\mathbf{x}) P\_L(\mathbf{x}) P\_{l'}(\mathbf{x}) d\mathbf{x} = 2 \begin{pmatrix} l & L & l' \\ 0 & 0 & 0 \end{pmatrix}^2,\tag{137}$$

where the RHS represents twice the square of the Wigner 3-j symbol.

Thus, applying the methodologies outlined in Section 2.4, and given the above features, one obtains

$$\bar{f}\_2 = \frac{1}{48} \sum\_{l=0}^{\infty} \frac{\bar{\zeta}\_l(\rho) P\_l(\cos \theta)}{(2l-1)(2l+3)},\tag{138}$$

where

$$\mathcal{J}\_l(\rho) = \begin{cases} \bar{\chi}\_l(\rho)\_{\prime} & 0 \le \rho \le 1 \\ \bar{\chi}\_l(1/\rho)\_{\prime} & \rho \ge 1 \end{cases} \tag{139}$$

with

$$\begin{split} \dot{\chi}\_{l}(\rho) &= -\frac{\rho^{l}}{(\rho^{2}+1)^{7/2}} \Big\{ \frac{(32l^{2}+26l-25)\rho^{6}}{2l+5} \Big[ \frac{(2l-1)\rho^{2}}{2l+9} + 4 \right] + \\ &+ \frac{1}{2l-3} \Big[ \frac{6(84l^{2}+84l-95)\rho^{4}}{2l+5} - (32l^{2}+38l-19) \Big( \frac{2l+3}{2l-7} + 4\rho^{2} \Big) \Big] .\end{split} \tag{140}$$

Recall that the component ˘ *f*2 in the RHS of Equation (126) is equal to 60(688 − 225 *π*) ¯ *f*2 according to representation (131).

As in the case of the AFC *ψ*5,2(*<sup>α</sup>*, *<sup>θ</sup>*), there are combinations of special hyperspherical angles *α* and *θ* for which the component ¯ *f*2 ≡ ¯ *f*2(*<sup>α</sup>*, *θ*) of the AFC *ψ*7,3(*<sup>α</sup>*, *θ*) can be obtained in closed form. In particular, one obtains:

$$\bar{f}\_2(a,0) = \mp \frac{(\rho - 1)(95\rho^6 + 1166\rho^5 - 1879\rho^4 - 8844\rho^3 - 1879\rho^2 + 1166\rho + 95)}{5040(\rho^2 + 1)^{7/2}},\tag{141}$$

$$f\_2(a,\pi) = \frac{(\rho+1)(95\rho^6 - 1166\rho^5 - 1879\rho^4 + 8844\rho^3 - 1879\rho^2 - 1166\rho + 95)}{5040(\rho^2 + 1)^{7/2}},\tag{142}$$

$$\bar{f}\_2(\alpha, \frac{\pi}{2}) = \frac{19\rho^4 + 10\rho^2 + 19}{1008(\rho^2 + 1)^2}. \tag{143}$$

Sign "−" in Equation (141) corresponds to 0 ≤ *α* ≤ *π*/2 (0 ≤ *ρ* ≤ 1), whereas "+" corresponds to *π*/2 ≤ *α* ≤ *π* (*ρ* ≥ 1).

For the important cases of the nucleus–electron and electron–electron coalescence, representation (138)–(140) is simplified to:

$$
\bar{f}\_2(0, \theta) = \frac{19}{1008}, \qquad \bar{f}\_2\left(\frac{\pi}{2}, 0\right) = 0. \tag{144}
$$

*4.2. The AFC ψ*8,4(*<sup>α</sup>*, *θ*)

Having at our disposal the AFC *ψ*7,3 ≡ *ψ*7,3(*<sup>α</sup>*, *<sup>θ</sup>*), we can calculate the AFC *ψ*8,4 ≡ *ψ*8,4(*<sup>α</sup>*, *θ*) using the FRR (3) and (4) for *k* = 8 and *p* = 3:

$$
\Delta \left( \Lambda^2 - 96 \right) \psi\_{8,3} = 80 \psi\_{8,4} - 2V \psi\_{7,3} + 2E \psi\_{6,3}.\tag{145}
$$

It follows from Equation (81) and the FRR (3) and (4) for *k* = 8 and *p* = 4 that the AFC *ψ*8,4 is a linear combination of the HHs, *<sup>Y</sup>*8,*l* ≡ *<sup>Y</sup>*8,*<sup>l</sup>*(*<sup>α</sup>*, *<sup>θ</sup>*). Moreover, given that only *Yn*,*<sup>l</sup>*(*<sup>α</sup>*, *θ*) with even values of *n*/2 − *l* are suitable for singlet S-states, we obtain:

$$
\psi\_{8,4} = a\_{80}\chi\_{8,0} + a\_{82}\chi\_{8,2} + a\_{84}\chi\_{8,4}.\tag{146}
$$

For further derivations, it is advisable to represent the HHs in the form

$$\mathcal{Y}\_{8,l}(\mathfrak{a},\theta) = \mathcal{Y}\_{8l}(\mathfrak{a}) P\_l(\cos \theta), \tag{147}$$

where

$$y\_{80}(\alpha) = \pi^{-3/2} [2\cos(4\alpha) + 2\cos(2\alpha) + 1],\tag{148}$$

$$y\_{\text{82}}(a) = \frac{2}{\pi^{3/2}} \sqrt{\frac{10}{7}} \sin^2 a [4 \cos(2a) + 3],\tag{149}$$

$$y\_{84}(\alpha) = \frac{8}{\pi^{3/2}} \sqrt{\frac{2}{\mathcal{T}}} \sin^4 \alpha. \tag{150}$$

It was found in Section 3 that *ψ*6,3 ≡ *ψ*6,3(*<sup>α</sup>*, *θ*) is the linear combination of the HHs *<sup>Y</sup>*6,*<sup>l</sup>*(*<sup>α</sup>*, *<sup>θ</sup>*). Thus, expanding each function of Equation (145) in HHs, and equating the coefficients for *<sup>Y</sup>*8,*l*, we obtain (see the corresponding result (105) for *<sup>a</sup>*6*l*)

$$a\_{8l} = \frac{1}{40} \int V \psi\_{7,3} \chi\_{8\downarrow} d\Omega\_{\prime} \tag{151}$$

where the potential *V* is defined by Equation (5). When deriving the last equation, it was taken into account that (Δ<sup>2</sup> − <sup>96</sup>)*<sup>Y</sup>*8,*l* = 0, as follows from Equation (81).

A direct substitution of the representations (5), (126) and (147) into the RHS of Equation (151) yields:

$$a\_{8l} = \frac{(\pi - 2)(5\pi - 14)Z^3}{13608000\sqrt{5}\ \pi^{3/2}} \times$$

$$\chi \times \int \left(\frac{1}{\xi} - \frac{2Z\eta}{\sin a}\right) \left\{\frac{\check{f}\_1 + \check{f}\_2}{\sqrt{5}\pi^{3/2}} - 2Z\left[12(32\pi - 97)\check{f}\_3 + (357\pi - 1112)\check{f}\_4\right]\right\} Y\_{8J}(u, \theta) d\Omega. \tag{152}$$

It follows from Equation (13) of Ref. [4] that only the coefficients at *Z*<sup>4</sup> are nonzero on the RHS of the last equation. Hence, Equation (152) reduces to the form:

$$a\_{8l} = -\frac{(\pi - 2)(5\pi - 14)Z^4}{6804000\sqrt{\pi}} \left[ \frac{\mathcal{S}\_{1l} + \mathcal{S}\_{2l}}{\sqrt{5}\pi^{3/2}} + 12(32\pi - 97)\mathcal{S}\_{3l} + (357\pi - 1112)\mathcal{S}\_{4l} \right],\tag{153}$$

where

$$S\_{1l} = \pi^2 \int\_0^\pi \int\_0^\pi \check{f}\_1(\check{\xi}) \check{Y}\_{8,l}(a,\theta) \eta \sin a \sin \theta d\alpha d\theta,\tag{154}$$

$$S\_{2l} = \pi^2 \int\_0^\pi \int\_0^\pi \check{f}\_2(a,\theta) \chi\_{8J}(a,\theta) \eta \sin a \sin \theta da d\theta = $$

$$= \frac{5\pi^2 (688 - 225\pi)}{(2l - 1)(2l + 1)(2l + 3)} \int\_0^{\pi/2} \check{\chi}\_l(\rho) y\_{8l}(a) \eta \sin a da,\tag{155}$$

$$S\_{nl} = \pi^2 \int\_0^\pi \int\_0^\pi \check{f}\_n(a,\theta) \check{\chi}\_{8,l}(a,\theta) \check{\xi}^{-1} \sin^2 a \sin\theta da d\theta. \qquad (n=3,4) \tag{156}$$

The identifiers *ξ* and *η* are defined by Equation (6), whereas functions *<sup>χ</sup>*¯*l*(*ρ*) can be calculated by Formula (140). When deriving Equation (155), we applied the orthogonality condition for the Legendre polynomials. Fortunately, all integrals (154)–(156) can be taken in closed form. Thus, by collecting these results and substituting them into the RHS of Equation (153), we finally obtain the desired coefficients in the form:

$$a\_{8l} = \frac{Z^4(\pi - 2)(5\pi - 14)}{\pi^{5/2}} b\_{8l\prime} \tag{157}$$

with

$$b\_{80} = \frac{\pi (150339\pi - 927292) + 1430792}{19289340000}, \quad b\_{82} = \frac{\pi (751965\pi - 4654046) + 7200976}{1928934000\sqrt{70}},$$

$$b\_{84} = \frac{\pi (3190317\pi - 19828996) + 30802176}{25719120000\sqrt{14}}.\tag{158}$$

#### **5. Results and Discussions**

The angular Fock coefficients *ψ<sup>k</sup>*,*<sup>p</sup>* ≡ *ψ<sup>k</sup>*,*<sup>p</sup>*(*<sup>α</sup>*, *θ*) with the *maximum possible* value of subscript *p* were calculated on examples of the coefficients with 5 ≤ *k* ≤ 10. The results obtained in Sections 2–4 are summarized in Appendices A and B. The AFCs *ψ*9,4 and *ψ*10,5 are presented in Appendix C without derivations. To find the latter AFCs, the methods described in the main sections were used. The presented technique makes it possible to calculate such AFCs for any arbitrarily large *k*. These coefficients are leading in the logarithmic power series representing the Fock expansion (see Equation (8)). As such, they may be indispensable for the development of simple methods for calculating the helium-like electronic structure.

The proposed technique, as well as the final results, are quite complex. Therefore, both require verification. We are aware of two ways for the above-mentioned verification. The first one is to use the Green's function (GF) approach (see Ref. [1] and also Ref. [15], Section 4) which enables us to calculate (at least, numerically) the AFC (or its component) by the following integral representation:

$$\psi\_{k,p}(a,\theta) = \frac{1}{8\pi} \int\_0^\pi da' \sin^2 a' \int\_0^\pi d\theta' \sin\theta' \, h\_{k,p}(a',\theta') \int\_0^\pi \frac{\cos\left[\left(\frac{k}{2} + 1\right)\omega\right]}{\sin\omega} (1-\lambda) d\rho,\tag{159}$$

where *ω* is an angle defined by the relation

$$\cos\omega = \cos a \cos a' + \sin a \sin a' \left(\cos\theta \cos\theta' + \sin\theta \sin\theta' \cos\varphi\right), \tag{160}$$

whereas

$$
\lambda = \begin{cases} \ 0 \\ \ \omega/\pi \end{cases} \qquad \begin{array}{cc} k \ \operatorname{odd} \\ \ k \ \operatorname{even} \end{array} . \tag{161}
$$

For even *k* and maximum value of *p* = *k*/2, the RHS *hk*,*k*/2 of the FRR (3) equals zero. This implies that the GF formula (159) cannot be applied in this case. Hence, only the AFCs *ψ<sup>k</sup>*,*<sup>p</sup>* with odd values of *k* (and maximum *p*) can be verified with the GF method. Thus, numerically calculating (for various combinations of angles *α* and *θ*) the triple integrals (159) representing the AFCs *ψ*5,2(*<sup>α</sup>*, *<sup>θ</sup>*), *ψ*7,3(*<sup>α</sup>*, *θ*) and *ψ*9,4(*<sup>α</sup>*, *<sup>θ</sup>*), we verified that the representations obtained for them in Sections 2 and 4.1 and in Appendix C were correct.

The second verification method considered, covering all possible combinations of angles, being quite complex, is the only method known to us that correctly displays the WF near the nucleus. This is the CFHH method mentioned in the Introduction. It is based on decomposing the full WF into a form

$$\Psi^{\rm CFHH}(r\_1, r\_2, r\_{12}) = \exp[f(r\_1, r\_2, r\_{12})] \Phi^{\rm CFHH}(\mathbb{R}, \mathfrak{a}, \theta),\tag{162}$$

where the so-called correlation function *f* can be taken in a simple linear form

$$f(r\_1, r\_2, r\_{12}) = c\_1 r\_1 + c\_2 r\_2 + c\_{12} r\_{12}.\tag{163}$$

The so-called "cusp parametrization"

$$
\mathfrak{c}\_1 = \mathfrak{c}\_2 = -Z, \qquad \mathfrak{c}\_{12} = 1/2 \tag{164}
$$

is used as a rule. For a small enough hyperspherical radius *R*, the function Φ is represented as

$$\Phi^{\rm CFHHH}(\mathbb{R}, \mathfrak{a}, \theta) = \frac{1}{d\_{0,0}(\mathfrak{a}, \theta)} \sum\_{k=0}^{K} (2\kappa R)^k \sum\_{p=0}^{\left| k/2 \right|} d\_{k,p}(\mathfrak{a}, \theta) \ln^p(2\kappa R), \tag{165}$$

where *κ* = √−2*E*, and functions *dk*,*<sup>p</sup>*(*<sup>α</sup>*, *θ*) are expanded in *N* (basis size) HHs. It follows from representation (165) that the AFCs *ψ<sup>k</sup>*,*<sup>p</sup>*(*<sup>α</sup>*, *θ*) can be expressed in terms of the functions *dk*,*p*(*<sup>α</sup>*, *θ*) calculated by the CFHHM. For example, for the AFCs in question, one obtains:

$$
\psi\_{k,k/2}^{\text{CFHH}}(\mathfrak{a},\theta) = \frac{(2\kappa)^k d\_{k,k/2}(\mathfrak{a},\theta)}{d\_{0,0}(\mathfrak{a},\theta)}.\tag{166}
$$

We calculated all AFCs discussed in this article using CFHHM with *K* = 18 and *N* = 1600. The angles 0 ≤ *α* ≤ *π* and 0 ≤ *θ* ≤ *π* with step *π*/6 were considered. The relative difference |1 − *ψ<sup>k</sup>*,*<sup>p</sup>*(*<sup>α</sup>*, *θ*)/*ψ*CFHH *k*,*p* (*<sup>α</sup>*, *θ*)| was less than 10−<sup>7</sup> for all examined cases, including 1 ≤ *Z* ≤ 5. This indicates that all our theoretical calculations were correct.

**Author Contributions:** Conceptualization, E.Z.L.; Methodology, E.Z.L.; Software, E.Z.L. and R.K.; Validation, R.K.; Writing—original draft, E.Z.L. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research received no external funding.

**Conflicts of Interest:** The authors declare no conflict of interest.
