text:
# Analytic Evaluation of Four-Particle Integrals with Complex Parameters

## I Introduction

The problem of four particles with the Coulomb interaction plays an important role in atomic and nuclear physics. It forms a link between the three-body problem, which can be solved numerically with very high precision, and many-body problems, solutions of which are very approximate. Thus, profound studies of various four-particle systems can provide valuable insights into physics of systems with greater numbers of particles.

In addition to the methodological interest, the four-body problem has unquestionable practical significance. Positronium beams are extensively used in positronium-atom scattering experiments, but the positronium molecule, $`e^+e^{}e^+e^{}`$, has not been observed experimentally yet. All existing knowledge of its properties is based on numerical studies usu98 . Molecules and ions including $`\mu `$-meson have attracted much attention traditionally in connection with the problem of the muon catalyzed fusion. Calculations suggest fro97  that muonic molecules like $`p^+\mu ^{}p^+\mu ^{}`$ have higher nuclear reaction rates than the corresponding three-particle ions. These examples show that high-precision numerical solution of the four-body problem is essential for proper understanding of various physical phenomena.

The majority of four-particle systems are nonadiabatic, and cannot be treated within the adiabatic approximation. The only practical way to calculate their energy and properties is to use the variational approach, taking into account the correlated motion of all the particles. Basis functions of the Gaussian type, depending on six interparticle separations and several nonlinear parameters, have been extensively used for such calculations usu98 ; fro97 ; suz98 . An important advantage of the Gaussian functions is that all integrals can be easily evaluated. The nonlinear parameters are optimized stochastically suz98 : at each step of basis expansion, many functions with randomly generated parameters are examined, and the function, giving the largest decrease in energy, is added to the basis.

However, unlike real wavefunctions, the Gaussian functions do not decay exponentially, and do not satisfy the cusp condition. From this point of view, they are rather unphysical. As a result, convergence of the variational procedure is very slow, and many hundreds of basis functions must be used. A recent calculation of the positronium molecule usu98  involved 1600 Gaussian functions. It was suggested that further expansion of the basis was not practical because of increasing computation time and low probability of finding good parameters. Thus, more efficient basis functions are clearly required.

A method for analytic evaluation of four-particle integrals, proposed by Fromm and Hill hil87 , opened up possibilities of variational calculation of four-particle systems in a basis of exponential functions of interparticle separations. This method reduces computation of integrals, needed to determine matrix elements of a four-particle Hamiltonian, to evaluation of the dilogarithmic function lev81  of various arguments. Application of this method, however, is a very difficult problem. Because the dilogarithm is a multiple valued function, the entire algorithm cannot be used without an effective procedure of branch and singularity tracking.

This problem was initially solved by the authors for the case of real exponential parameters. The first calculations of the positronium molecule reb96 , and several mesic molecules zot98  in the exponential basis, depending on all six interparticle separations, have demonstrated high efficiency and great potential of this method. To the best of our knowledge, nobody else has done this yet fro01 .

Because one exponential function is as effective as eight Gaussians, a size of the basis can be reduced significantly. However, an amount of time, needed to compute one matrix element, is much larger than for the Gaussian basis. Thus, optimization of nonlinear parameters is the main difficulty. Deterministic optimization (gradient descent) gives excellent results for a relatively small number of exponential basis functions. Stochastic optimization (trial and error), used to expand the basis further, is inefficient due to a dramatic increase in computation time. This fact suggests that a possible alternative to an enormously large Gaussian basis is a relatively short basis of the most efficient and versatile functions with carefully optimized parameters.

A natural generalization of the exponential basis is the exponential-trigonometric basis, obtained by replacing real exponential parameters with complex ones reb90 . The exponential-trigonometric functions have been successfully employed in variational calculations of three-particle adiabatic systems zot94 . They are much more efficient, than the ordinary exponentials for two reasons. First, they contain twice as many nonlinear parameters, thus allowing better approximation of the wavefunction. Second, they exhibit nonmonotonic dependence on interparticle separations, being able to imitate sharp peaks in wavefunctions of adiabatic systems. The computation time increases only insignificantly in comparison with the case of real exponential parameters.

In order to use the exponential-trigonometric basis in the four-body problem, one has to evaluate the four-particle integrals with complex parameters. The problem of branch tracking in a general complex case is formidable. Every branch change for every multiple valued function has to be taken into account if correct values of the integrals are to be obtained. An original (and, inevitably, very nontrivial) procedure of numerical branch tracking has been developed by the authors. The first variational calculations of four-particle systems in the exponential-trigonometric basis proved extremely promising zot00 . They showed that one exponential-trigonometric function can replace seven exponential functions in calculation of $`e^+e^{}e^+e^{}`$, and several dozen exponentials in studies of adiabatic systems zot00 . Therefore, it presents a real alternative to both the exponential and Gaussian basis functions.

Even though results of the calculations involving the exponential and exponential-trigonometric functions have been published reb96 ; zot98 ; zot00 , details of the new method have not been reported yet. The purpose of the present paper is to fill this gap. We present a description of our algorithm that will enable a reader to implement it as a computer program.

The paper is organized as follows. Section II.A discusses what integrals are needed to compute matrix elements of a four-particle Hamiltonian, and how the number of them can be reduced. In Sec. II.B, principles of the original method by Fromm and Hill are outlined. Sec. II.C provides information about multiple valued functions used in the analysis. In Sec. II.D, a simplified procedure of branch tracking in the case of real parameters is described. Sec. II.E gives a detailed exposition of the method of branch tracking in the most general case, when all the parameters are complex. In Sec. II.F, a practical implementation of the branch tracking algorithm is described. The last Section presents our conclusions.

## II Description of the method

### II.1 Matrix elements of four-particle Hamiltonian

Let us consider a Hamiltonian of a four-particle system with the Coulomb interactions:

$$H=\frac{\mathrm{}^2}{2}\underset{j=1}{\overset{4}{}}\frac{\mathrm{\Delta }_j}{m_j}+\underset{j<k}{\overset{4}{}}\frac{q_jq_k}{r_{jk}}.$$
(1)
Here $`m_j`$ and $`q_j`$, $`j=1\mathrm{}4`$, are masses and charges, and $`r_{jk}=|𝐫_j𝐫_k|`$ are interparticle separations. Our purpose is to evaluate matrix elements of $`H`$ with exponential basis functions

$$\mathrm{\Phi }_b=\mathrm{exp}(\underset{j<k}{\overset{4}{}}b_{jk}r_{jk}),\mathrm{\Phi }_c=\mathrm{exp}(\underset{j<k}{\overset{4}{}}c_{jk}r_{jk}).$$
(2)
These functions depend on complex parameters $`\{b_{jk}\}`$ and $`\{c_{jk}\}`$. In what follows, the notation $`\{x_{jk}\}`$ will always refer to six quantities, $`x_{jk}`$, with $`j,k=1\mathrm{}4`$ and $`j<k`$, i.e. $`x_{12},x_{13},x_{14},x_{23},x_{24},x_{34}`$, assuming that $`x_{jk}=x_{kj}`$.

In order to compute matrix elements of the operator of kinetic energy in Eq. (1), one has to evaluate the following quantities:

$$\mathrm{\Phi }_b\mathrm{cos}\mathrm{\Theta }_{jkl}\mathrm{\Phi }_c=\mathrm{\Phi }_b\frac{r_{jk}^2+r_{jl}^2r_{kl}^2}{2r_{jk}r_{jl}}\mathrm{\Phi }_c,$$
(3)
where $`jk`$, $`kl`$, $`jl`$. The integrands in the last formula display linear and even quadratic dependence on certain interparticle separations. Therefore, in order to obtain matrix elements of the Hamiltonian, Eq. (1), one has to calculate a total of 43 integrals: one overlap integral, six integrals of the Coulomb interactions, and 36 integrals, given by Eq. (3).

It turns out, however, that it is possible to avoid computation of the integrals in Eq. (3). It has been shown by one of the authors that the matrix elements of the above Hamiltonian can be expressed in terms of the overlap integral and six Coulomb integrals only reb93 . Thus, one can write:

$$\mathrm{\Phi }_bH\mathrm{\Phi }_c=H_1H_2H_3.$$
(4)
The individual terms in Eq. (4) are given by the following expressions reb93 :

$$H_1=\underset{j<k}{\overset{4}{}}\left[\frac{(m_j+m_k)}{2m_jm_k}a_{jk}+q_jq_k\right]\mathrm{\Phi }_a\frac{1}{r_{jk}}\mathrm{\Phi }_a;$$
(5)
$$H_2=\underset{j<k}{\overset{4}{}}\frac{(m_j+m_k)}{2m_jm_k}d_{jk}^{\mathrm{\hspace{0.17em}2}}\mathrm{\Phi }_a\mathrm{\Phi }_a;$$
(6)
$`H_3={\displaystyle \underset{j=1}{\overset{4}{}}}{\displaystyle \underset{\begin{array}{c}k<l\\ k,lj\end{array}}{\overset{4}{}}}{\displaystyle \frac{(a_{jk}s_{jk}+a_{jl}s_{jl}a_{jn}s_{jn})}{2m_ja_{jk}a_{jl}}}d_{jk}d_{jl}.`$ (9)
In these formulas, $`\mathrm{\Phi }_a`$ is a new function with parameters $`\{a_{jk}\}`$, defined as $`a_{jk}=(b_{jk}+c_{jk})/2`$:

$$\mathrm{\Phi }_a=\mathrm{exp}(\underset{j<k}{\overset{4}{}}a_{jk}r_{jk}).$$
(10)
The parameters $`\{d_{jk}\}`$ are defined as $`d_{jk}=(c_{jk}b_{jk})/2`$, and the quantities $`\{s_{jk}\}`$ are given by:

$$s_{jk}=\mathrm{\Phi }_a\frac{1}{r_{jk}}\mathrm{\Phi }_aa_{jk}\mathrm{\Phi }_a\mathrm{\Phi }_a.$$
(11)
The additional index $`n`$ in Eq. (7) is fixed by a condition $`nj,k,l`$.

Therefore, only seven integrals – the overlap integral and six Coulomb integrals, calculated with the function $`\mathrm{\Phi }_a`$, – are needed to determine the matrix elements of the Hamiltonian, Eq. (1). The above formulas are valid for both real and complex parameters. They are indispensable for any application of this method.

### II.2 Evaluation of four-particle integrals

In this paper, we generalize the method of analytic evaluation of four-particle integrals, proposed by Fromm and Hill hil87 , to include complex exponential parameters. First, we would like to recall basic ideas of this method. The following family of integrals is considered:

$$J(\{n_{jk}\},\{\alpha _{jk}\})=(\underset{j<k}{\overset{4}{}}r_{jk}^{n_{jk}1})\mathrm{exp}(\underset{j<k}{\overset{4}{}}\alpha _{jk}r_{jk})𝑑V.$$
(12)
Here, $`\{\alpha _{jk}\}`$ denotes a set of six exponential parameters, $`\alpha _{12},\alpha _{13},\alpha _{14},\alpha _{23},\alpha _{24},\alpha _{34}`$, and $`\{n_{jk}\}`$ is the corresponding set of non-negative integers. The integrand depends on six interparticle separations, $`\{r_{jk}\}`$. The integration is performed over 9-dimensional space of relative coordinates of four particles: $`dV=d^{\mathrm{\hspace{0.17em}3}}r_{12}d^{\mathrm{\hspace{0.17em}3}}r_{13}d^{\mathrm{\hspace{0.17em}3}}r_{14}`$.

An integral with all $`n_{jk}=0`$ is called “generating”:

$$I(\{\alpha _{jk}\})=(\underset{j<k}{\overset{4}{}}r_{jk}^1)\mathrm{exp}(\underset{j<k}{\overset{4}{}}\alpha _{jk}r_{jk})𝑑V.$$
(13)
All the integrals in Eq. (10) can be obtained from the generating integral, Eq. (11), by differentiation:

$$J(\{n_{jk}\},\{\alpha _{jk}\})=[\underset{j<k}{\overset{4}{}}(\frac{}{\alpha _{jk}})^{n_{jk}}]I(\{\alpha _{jk}\}).$$
(14)
The generating integral is given by the following formula:

$$I(\{\alpha _{jk}\})=\frac{16\pi ^3}{\sigma }[\underset{j=1}{\overset{4}{}}\underset{k=1}{\overset{4}{}}v(\gamma _k^{(j)}/\sigma )+\underset{j=2}{\overset{4}{}}u(\beta _1^{(1)}\beta _1^{(j)})].$$
(15)
The functions $`v(z)`$ and $`u(z)`$ are expressed in terms of the dilogarithmic function $`\mathrm{Li}_2(z)`$:

$$u(z)=\mathrm{Li}_2(z)\mathrm{Li}_2(1/z);$$
(16)
$`v(z)={\displaystyle \frac{1}{2}}\mathrm{Li}_2[(1z)/2]{\displaystyle \frac{1}{2}}\mathrm{Li}_2[(1+z)/2]`$ (17)
$`{\displaystyle \frac{1}{4}}\mathrm{ln}^2[(1z)/2]+{\displaystyle \frac{1}{4}}\mathrm{ln}^2[(1+z)/2].`$
In Eq. (13) for the generating integral, $`\gamma _k^{(j)}`$ are third-order polynomials in $`\alpha `$’s, defined in the following way:

$`\gamma _k^{(j)}=\mu _j^{(j)}\mu _k^{(j)}+\mu _l^{(j)}+\mu _m^{(j)},`$ (18)
$`\gamma _j^{(j)}=+\mu _1^{(j)}+\mu _2^{(j)}+\mu _3^{(j)}+\mu _4^{(j)},`$
where for each $`jk`$: $`lj,k`$; $`mj,k`$; $`lm`$. The polynomials $`\mu _k^{(j)}`$ are defined as follows:

$`\mu _k^{(j)}=\alpha _{lm}(\alpha _{jk}^2+\alpha _{kl}^2+\alpha _{km}^2),`$ (19)
$`\mu _j^{(j)}=2\alpha _{lm}\alpha _{kl}\alpha _{km},`$
with the same restrictions on values of $`j,k,l,`$ and $`m`$.

The function $`\sigma `$ is a square root of a sixth-order polynomial in $`\alpha `$’s: $`\sigma =\sqrt{s_1+s_2}`$. The quantity $`s_1`$ in this expression is given by

$$s_1=\underset{j=2}{\overset{4}{}}\alpha _{1j}^2\alpha _{lm}^2(\alpha _{1j}^2+\alpha _{lm}^2\alpha _{1l}^2\alpha _{1m}^2\alpha _{jl}^2\alpha _{jm}^2),$$
(20)
where for each $`j`$: $`l1,j`$; $`m1,j`$; $`lm`$. The quantity $`s_2`$ is determined as

$$s_2=\underset{j=1}{\overset{4}{}}\alpha _{jl}^2\alpha _{jm}^2\alpha _{jk}^2,$$
(21)
where for each $`j`$: $`l,m,kj`$; $`lm`$; $`mk`$; $`lk`$. Finally, $`\beta _k^{(j)}`$ is defined by the following expression:

$$\beta _k^{(j)}=(\sigma \gamma _k^{(j)})/(\sigma +\gamma _k^{(j)}).$$
(22)
In all these formulas, indices $`j,k,l,m`$ change from 1 to 4, and it is assumed that $`\alpha _{jk}=\alpha _{kj}`$ for each $`jk`$. If some indices are not defined uniquely, the formulas are symmetric under their permutations.

Eq. (13) is the main result of this method hil87 . It provides an analytic expression for the generating integral, Eq. (11). It was pointed out hil87  that there is no need to know an analytic dependence of the generating integral on the parameters $`\{\alpha _{jk}\}`$ to compute the family of integrals, Eq. (10). According to Eq. (12), all these integrals are derivatives of the generating integral. Special formulas can be used hil87  to calculate numerical values of derivatives of functions $`fg`$ and $`h(g)`$, if numerical values of derivatives of the functions $`f,g`$, and $`h`$ have already been computed. For example, derivatives of the term $`v(\gamma _k^{(j)}/\sigma )`$ in Eq. (13) can be obtained in the following way. First, derivatives of $`\sigma ^2`$ and $`\gamma _k^{(j)}`$ with respect to $`\{\alpha _{jk}\}`$ are calculated. Then derivatives of a function $`h(z)=z^{1/2}`$ are computed at $`z=\sigma ^2`$. After that, using a formula for derivatives of $`h(g)`$ with $`g=\sigma ^2`$, one finds derivatives of $`1/\sigma `$ with respect to $`\{\alpha _{jk}\}`$. Then, using a formula for derivatives of $`fg`$ with $`f=\gamma _k^{(j)}`$ and $`g=1/\sigma `$, one obtains derivatives of $`\gamma _k^{(j)}`$ with respect to $`\{\alpha _{jk}\}`$. After that, derivatives of a function $`h(z)=v(z)`$ at $`z=\gamma _k^{(j)}/\sigma `$ are calculated. Finally, using a formula for derivatives of $`h(g)`$ with $`g=\gamma _k^{(j)}/\sigma `$, one can find the derivatives of $`v(\gamma _k^{(j)}/\sigma )`$ with respect to $`\{\alpha _{jk}\}`$. Within this approach, all the integrals of Eq. (10) can be evaluated by means of an efficient recursive procedure, working with numbers only.

At this point, we can appreciate importance of Eqs. (4)-(9). In order to obtain the matrix elements of the Hamiltonian, Eq. (1), we have to compute the mixed derivatives given by Eq. (12) up to the sixth order only, i.e. for $`n_{jk}=0,1`$, where $`j,k=1\mathrm{}4`$, and $`j<k`$. This means that, at each step of the recursive procedure, we calculate $`2^6=64`$ derivatives. If we tried to evaluate the integrals of Eq. (3) directly, it would be necessary to compute the mixed derivatives in Eq. (12) up to 18th order, i.e. for $`n_{jk}=0\mathrm{}3`$. The number of derivatives, calculated at each step, would increase quadratically. An amount of time, required to carry out the entire recursive procedure, would be enormous. Therefore, the original method by Fromm and Hill hil87 , used by itself, does not make high precision calculations of four-particle systems possible. Only in conjunction with the method reb93  for reducing the number of integrals can it produce valuable results.

### II.3 The multiple valued functions

The main difficulty in using Eq. (13) for the generating integral is the fact that the functions in this formula are multiple valued. Indeed, the functions $`u(z)`$ and $`v(z)`$, given by Eqs. (14) and (15), are expressed in terms of the dilogarithmic function $`\mathrm{Li}_2(z)`$. The dilogarithm is defined as follows lev81 :

$$\mathrm{Li}_2(z)=_0^z\frac{\mathrm{ln}(1\zeta )}{\zeta }𝑑\zeta .$$
(23)
This function is analytic inside the unit circle in the complex plane:

$$\mathrm{Li}_2(z)=\underset{n=1}{\overset{\mathrm{}}{}}\frac{z^n}{n^2},|z|<1.$$
(24)
Its values outside the unit circle can be determined using a relation lev81 :

$$\mathrm{Li}_2(z)=\frac{\pi ^2}{6}\frac{1}{2}\mathrm{ln}^2(z)\mathrm{Li}_2(1/z).$$
(25)
In the immediate vicinity of the unit circle, where convergence of the series in Eq. (22) is slow, the following relations can be used to shift the argument of $`\mathrm{Li}_2(z)`$:

$$\mathrm{Li}_2(z)=\frac{\pi ^2}{6}\mathrm{ln}(z)\mathrm{ln}(1z)\mathrm{Li}_2(1z).$$
(26)
$$\mathrm{Li}_2(z)=\frac{1}{2}\mathrm{Li}_2(z^2)\mathrm{Li}_2(z).$$
(27)
Presence of the logarithm in Eqs. (23) and (24) clearly indicates that the function $`\mathrm{Li}_2(z)`$ is, in general, multiple valued. In order to specify its principal branch we need to fix the principal branch of the logarithm. The complex logarithm has branch points at 0 and $`\mathrm{}`$. We choose its branch cut to run along the negative real axis and define the principal branch as follows:

$$\mathrm{ln}(z)=\mathrm{ln}|z|+i\mathrm{arg}z,\pi <\mathrm{arg}z<\pi .$$
(28)
This choice determines branch cuts and fixes the principal branch for the dilogarithm, and the functions $`u(z)`$ and $`v(z)`$.

The function $`\mathrm{Li}_2(z)`$ has branch points at 1 and $`\mathrm{}`$; its branch cut runs from 1 to $`\mathrm{}`$ along the positive real axis. The function $`u(z)`$ has branch points at 0, 1, and $`\mathrm{}`$; its branch cut goes from 0 to $`\mathrm{}`$ along the positive real axis. The function $`v(z)`$ has branch points at 1, -1, and $`\mathrm{}`$; its branch cuts run from $`\mathrm{}`$ to -1 along the negative real axis, and from 1 to $`\mathrm{}`$ along the positive real axis.

Fig. 1 exhibits the branch points and cuts for these multiple valued functions.

It is important to note that the function $`\sigma `$, which is present in Eq. (13) and defined using Eqs. (18) and (19), is also a multiple valued function. The complex square root has branch points at 0 and $`\mathrm{}`$. We choose its branch cut to run along the positive real axis and define the principal branch as follows:

$$\sqrt{z}=\sqrt{|z|}\mathrm{exp}(\frac{i}{2}\mathrm{arg}z),0<\mathrm{arg}z<2\pi .$$
(29)
It can be seen from the definition of the generating integral, Eq. (11), that it is a continuous function of parameters $`\{\alpha _{jk}\}`$ for all values of these parameters satisfying the following conditions:

$`\alpha _{12}+\alpha _{13}+\alpha _{14}>0,`$ $`\alpha _{12}+\alpha _{23}+\alpha _{24}>0,`$ (30)
$`\alpha _{13}+\alpha _{23}+\alpha _{34}>0,`$ $`\alpha _{14}+\alpha _{24}+\alpha _{34}>0.`$
$`\alpha _{12}+\alpha _{13}+\alpha _{24}+\alpha _{34}>0,`$
$`\alpha _{12}+\alpha _{14}+\alpha _{23}+\alpha _{34}>0,`$ (31)
$`\alpha _{13}+\alpha _{14}+\alpha _{23}+\alpha _{24}>0.`$
These conditions mean, physically, that the wavefunction of a system of four particles decreases exponentially when any of the interparticle separations become infinitely large. If the parameters $`\{\alpha _{jk}\}`$ are complex, the above inequalities must be satisfied by their real parts.

The continuity of the generating integral, Eq. (11), implies that the right hand side of Eq. (13) is also a continuous function of $`\{\alpha _{jk}\}`$. This fact has two important consequences.

First, the multiple valued functions $`u(z)`$, $`v(z)`$, and $`\sigma (z)`$ in Eq. (13) remain continuous while their branches change. As a point of interest moves in 12-dimensional space of six complex parameters, $`\{\alpha _{jk}\}`$, the arguments of these functions move freely in the complex plane, and their branches change repeatedly. However, a computer can evaluate only the principal branch of the logarithm, given by Eq. (26), and the principal branch of the square root, given by Eq. (27). Therefore, only the principal branches of the functions $`\mathrm{Li}_2(z)`$, $`u(z)`$, $`v(z)`$, and $`\sigma (z)`$, defined in the complex plane with the branch cuts, can be calculated directly. Thus, a special procedure of branch tracking is necessary to restore continuity of these functions every time their arguments cross the branch cuts.

Second, all singularities, which different terms in Eq. (13) can have, cancel mutually. These singularities arise when $`\sigma =0`$, and when any of the following equalities are satisfied:

$`\alpha _{jl}+\alpha _{jm}+\alpha _{jn}=0,`$
$`\alpha _{jl}\alpha _{jm}+\alpha _{jn}=0,`$ (32)
$`\alpha _{jl}+\alpha _{jm}\alpha _{jn}=0,`$
where for each $`j=1\mathrm{}4`$: $`l,m,nj`$; $`lm`$; $`mn`$; $`ln`$. These singularities are unphysical, and should have no effect on the value of the generating integral. As a point under consideration moves in the space of the parameters $`\{\alpha _{jk}\}`$, the arguments of the functions $`u(z)`$ and $`v(z)`$ can frequently appear in the vicinity of the singular (branch) points. As a result, the values of these functions can exhibit considerable change, even if the parameters $`\{\alpha _{jk}\}`$ change only slightly. Therefore, a special procedure for dealing with the singularities is needed in order to carry out explicit cancellation of all singular terms.

This discussion demonstrates that the method of hil87  is impossible to use without an effective algorithm for numerical branch and singularity tracking.

### II.4 Branch tracking in the real case

Before discussing a general algorithm of branch tracking, it is beneficial to consider a particular case, when all the exponential parameters, $`\{\alpha _{jk}\}`$, are real numbers. Let us introduce the following parametrization:

$$\alpha _{jk}(p)=(\alpha _{jk}1)p+1,0p1.$$
(33)
As the real parameter $`p`$ changes from 0 to 1, the corresponding point in 6-dimensional space moves from (1,1,1,1,1,1) to $`(\alpha _{12},\alpha _{13},\alpha _{14},\alpha _{23},\alpha _{24},\alpha _{34})`$. If the parameters $`\{\alpha _{jk}\}`$ satisfy the conditions of Eqs. (28) and (29), the parameters $`\{\alpha _{jk}(p)\}`$ will satisfy these conditions for any $`p`$ between 0 and 1. Therefore, the generating integral, given by Eq. (13), must be a continuous function of $`p`$. It is known hil87  that Eq. (13) with the functions $`u(z)`$, $`v(z)`$, and $`\sigma (z)`$, represented by their principal branches, yields the correct value for the generating integral at the reference point (1,1,1,1,1,1). If this value changes continuously, as the parameter $`p`$ goes from 0 to 1, one can be sure that the generating integral will be computed correctly at the final point $`(\alpha _{12},\alpha _{13},\alpha _{14},\alpha _{23},\alpha _{24},\alpha _{34})`$. Therefore, continuity of this integral is a criterion of the correct branch tracking.

Let us define a function $`S(p)`$ as the sum in the square brackets of Eq. (13) when the parameters $`\{\alpha _{jk}(p)\}`$ are used instead of $`\{\alpha _{jk}\}`$:

$$S(p)=\underset{j=1}{\overset{4}{}}\underset{k=1}{\overset{4}{}}v(\gamma _k^{(j)}/\sigma )+\underset{j=2}{\overset{4}{}}u(\beta _1^{(1)}\beta _1^{(j)}).$$
(34)
Our purpose is to ensure that this function is continuous along the path from $`p=0`$ to $`p=1`$.

First, we consider a case when $`\sigma ^2(p)>0`$. The function $`\sigma (p)`$ is real, and all the arguments of the functions $`u(z)`$ and $`v(z)`$ in Eq. (32) are real as well. It will be shown in Sec. II.E that only imaginary parts of these functions exhibit discontinuities, when their arguments cross the branch cuts. Because the generating integral is real, the imaginary parts of the functions $`u(z)`$ and $`v(z)`$ in Eq. (32) must cancel anyway. Therefore, discontinuities in the real part of $`S(p)`$ may appear near the singular points of the functions $`u(z)`$ and $`v(z)`$ only. The singularities of different terms in Eq. (32) should cancel one another. However, because of possible branch changes, complete cancellation may not happen. The formulas of Sec. II.E suggest that, near the singular points of $`u(z)`$ and $`v(z)`$, the real part of the function $`S(p)`$ can undergo changes by $`m\pi ^2`$, where $`m`$ is some integer. Thus, the function $`S(p)`$ can have finite discontinuities, which are integer multiples of $`\pi ^2`$.

From now on, the branch tracking is only a technical problem. To solve it, it is necessary to find all values of the parameter $`p`$ between 0 and 1, which correspond to singular points. They include zeros of the sixth-order polynomial $`\sigma ^2(p)`$, and values of $`p`$, at which the parameters $`\{\alpha _{jk}(p)\}`$ satisfy any of the conditions of Eq. (30). Let us denote the resulting set of numbers as $`\{p_j\}`$, $`j=1\mathrm{}n`$. The correction function, needed to remove discontinuities of the function $`S(p)`$, is given by the following expression:

$$C(p)=\pi ^2\underset{p_j<p}{}\mathrm{Nint}[(S(p_j+ϵ)S(p_jϵ))/\pi ^2].$$
(35)
Here, the function $`\mathrm{Nint}[x]`$ returns an integer number, nearest to the real number $`x`$. The value of $`ϵ`$ in actual calculations was set to $`10^2`$. The correct value of the generating integral can now be determined from the formula:

$$I=\frac{16\pi ^3}{\sigma }[\underset{j=1}{\overset{4}{}}\underset{k=1}{\overset{4}{}}v(\gamma _k^{(j)}/\sigma )+\underset{j=2}{\overset{4}{}}u(\beta _1^{(1)}\beta _1^{(j)})+C(1)].$$
(36)
Therefore, in the case of the real parameters $`\{\alpha _{jk}\}`$, the procedure of branch tracking can be implemented without a detailed numerical analysis of behavior of the multiple valued functions. All we need to do is to calculate the function $`S(p)`$ twice for each singular point $`p_j`$, encountered along the path from $`p=0`$ to $`p=1`$, and subtract discontinuities, proportional to $`\pi ^2`$. The time, needed to determine the correction $`C(1)`$ in Eq. (34), is shorter than the time, required to carry out the recursive procedure for the family of integrals. It does not increase the overall computation time significantly.

The case of $`\sigma ^2(p)<0`$ is also straightforward. The quantity $`\sigma `$ is now imaginary. The function $`S(p)`$ is imaginary as well, thus giving a real value of the generating integral. $`\mathrm{Im}(S(p))`$ can be expressed in terms of Clausen’s function $`\mathrm{Cl}_2(\theta )`$, which is a real function of a real argument hil87 ; lev81 . Eqs. (33) and (34) are valid also in this case, if $`S(p)`$ is replaced by $`\mathrm{Im}(S(p))`$, and $`\sigma (p)`$ is replaced by $`\mathrm{Im}(\sigma (p))`$. Therefore, in both cases ($`\sigma ^2>0`$ and $`\sigma ^2<0)`$ the entire algorithm for analytic evaluation of the four-particle integrals can be presented in the real form without any use of complex numbers.

This simplified method of branch tracking has been successfully employed in variational calculations of four-particle systems reb96 ; zot98 . Therefore, the described method of branch tracking in the case of real exponential parameters is both theoretically correct and practically reliable.

### II.5 Branch tracking in the complex case

Let us now describe a method of branch tracking in a general case, when the exponential parameters, $`\{\alpha _{jk}\}`$, are complex numbers. It is assumed that their real parts satisfy Eqs. (28) and (29). We use the same parametrization as before, but with a complex parameter $`p`$ :

$$\alpha _{jk}(p)=(\alpha _{jk}1)p+1,0\mathrm{Re}(p)1.$$
(37)
As $`p`$ moves in the complex plane from 0 to 1, the corresponding point in 12-dimensional space of six complex parameters moves from (1,1,1,1,1,1) to $`(\alpha _{12},\alpha _{13},\alpha _{14},\alpha _{23},\alpha _{24},\alpha _{34})`$. The generating integral, Eq. (13), must be a continuous function of $`p`$. Moreover, its value, computed at the final point, $`\{\alpha _{jk}\}`$, should not depend on a choice of the path from $`p=0`$ to $`p=1`$. However, an optimal choice of this path can facilitate branch tracking considerably.

Fig. 2 exhibits three examples of paths in the complex $`p`$ plane. In case a), there are no singular points on or near the real axis between 0 and 1. The path is simply a straight lime segment between these points. In case b), there is one point, $`p_{\mathrm{\hspace{0.17em}1}}`$, at which different terms in Eq. (13) exhibit singular behavior. The path is the same as before, except for a small semicircle in the vicinity of this point. In case c), there are two singular points, $`p_{\mathrm{\hspace{0.17em}1}}`$ and $`p_{\mathrm{\hspace{0.17em}2}}`$, near the real axis between 0 and 1. The path is more complicated, as shown in the figure. In general, only those singular points in the $`p`$ plane, which are close to the real interval between 0 and 1, are of interest. The path should be carefully defined in the vicinity of every such point to allow precise analysis of behavior of all arguments of the multiple valued functions. The values of $`p`$, at which singularities may arise, can be found from the polynomial equation $`\sigma ^2(p)=0`$, and from twelve linear equations, contained in Eq. (30).

In order to obtain correction functions for the function $`u(z)`$, defined by Eq. (14), we have to consider behavior of this function near its branch points 0, 1, and $`\mathrm{}`$:

$`u(z0)={\displaystyle \frac{1}{2}}\mathrm{ln}^2(z)+u_{(0)}(z);`$
$`u(z1)=2\mathrm{ln}(z)\mathrm{ln}(1z)+u_{(1)}(z);`$ (38)
$`u(z\mathrm{})={\displaystyle \frac{1}{2}}\mathrm{ln}^2(z)+u_{(\mathrm{})}(z).`$
In these formulas, the functions with subscripts (0), (1), and $`(\mathrm{})`$ are functions, analytic in the vicinities of 0, 1, and $`\mathrm{}`$, respectively.

Let us introduce the following notations. A complex function $`z(p)`$ will represent any of the arguments, $`\beta _1^{(1)}\beta _1^{(j)}`$, of the function $`u(z)`$ in Eq. (13). It depends on $`p`$ through the parameters $`\alpha _{jk}(p)`$, given by Eq. (35). Let $`\{p_j\},j=1\mathrm{}N,`$ denote values of the parameter $`p`$, for which $`z(p_j)`$ are singular points 0, 1, $`\mathrm{}`$, or any points, where $`z(p)`$ crosses the real axis. It is assumed that $`0\mathrm{Re}(p_j)1`$ for each $`j=1\mathrm{}N`$, and $`\mathrm{Re}(p_j)<\mathrm{Re}(p_{j+1})`$. Each point $`p_j`$ will be characterized by an index $`n_j`$, and either integer $`m_j`$, or real $`\delta _j`$. The index $`n_j=1\mathrm{}6`$ specifies a type of singular behavior, as explained below. The number $`m_j`$ provides information about direction, in which the real axis is crossed by $`z(p)`$. We set $`m_j=+1`$, if the axis is crossed from below (i.e. $``$), and $`m_j=1`$, if it is crossed from above (i.e. $``$). The real quantity $`\delta _j`$ is equal to a change in $`\mathrm{arg}(z(p)z(p_j))`$, when $`z(p)`$ moves in the vicinity of a singular point $`z(p_j)`$. If $`z(p_j)=\mathrm{}`$, the quantity $`\delta _j`$ denotes a change in $`\mathrm{arg}(z(p))`$. These notations will allow us to present the algorithm of branch tracking as a series of formulas.

Five correction functions, $`u_{n_j}(z,j)`$, are needed to restore continuity of the computed function $`u(z)`$.

If $`z(p)`$ crosses the branch cut $`]\mathrm{\hspace{0.17em}1},+\mathrm{}[`$ at $`z(p_j)`$,
let $`n_j=1`$, and

$$u_1(z,j)=+2\pi ^22m_j\pi i[\mathrm{ln}(z)+U_j].$$
(39)
If $`z(p)`$ crosses the branch cut $`]\mathrm{\hspace{0.17em}0},1[`$ at $`z(p_j)`$,
let $`n_j=2`$, and

$$u_2(z,j)=2\pi ^2+2m_j\pi i[\mathrm{ln}(z)+U_j].$$
(40)
If $`z(p)`$ moves near the singular point $`z(p_j)=1`$,
let $`n_j=3`$, and

$$u_3(z,j)=2i\delta _j[\mathrm{ln}(z)+\widehat{U}_j].$$
(41)
If $`z(p)`$ moves near the singular point $`z(p_j)=0`$,
let $`n_j=4`$, and

$$u_4(z,j)=\delta _j^2/2i\delta _j[\mathrm{ln}(z)+U_j].$$
(42)
If $`z(p)`$ moves near the singular point $`z(p_j)=\mathrm{}`$,
let $`n_j=5`$, and

$$u_5(z,j)=+\delta _j^2/2+i\delta _j[\mathrm{ln}(z)+U_j].$$
(43)
If $`z(p)`$ crosses the branch cut $`]\mathrm{},0[`$ of the
function $`\mathrm{ln}(z)`$ at $`z(p_j)`$, let $`n_j=6`$.

The logarithms in these formulas are multiple valued functions themselves. Their branches can also change, and they can exhibit singular behavior, while an argument $`z(p)`$ moves further in the complex plane. Because only the principle branch of the logarithm is calculated by a computer, the additional terms, $`U_j`$ and $`\widehat{U}_j`$ are included to correct values of these functions. These terms are given by the following formulas:

$$U_j=+\underset{\begin{array}{c}k>j\\ n_k=1,2\end{array}}{\overset{N}{}}2m_k\pi i\underset{\begin{array}{c}k>j\\ n_k=4,5\end{array}}{\overset{N}{}}i\delta _k;$$
(44)
$$\widehat{U}_j=\underset{\begin{array}{c}k>j\\ n_k=6\end{array}}{\overset{N}{}}2m_k\pi i\underset{\begin{array}{c}k>j\\ n_k=4,5\end{array}}{\overset{N}{}}i\delta _k.$$
(45)
The condition $`n_k=1,2`$ in these formulas means that we have to sum up only those indices $`m_k`$, that correspond to situations, when $`z(p)`$ crosses the branch cuts $`]\mathrm{\hspace{0.17em}1},+\mathrm{}[`$ and $`]\mathrm{\hspace{0.17em}0},1[`$. The condition $`n_k=4,5`$ limits the summation of $`\delta _k`$ to those cases, when $`z(p)`$ moves near the singular points 0 and $`\mathrm{}`$. If $`n_k=6`$, we consider only situations, when $`z(p)`$ crosses the real axes in the interval $`]\mathrm{},0[`$.

Thus, each singular or crossing point $`z(p_j)`$, encountered by the argument $`z(p)`$ of the function $`u(z)`$, gives rise to a correction function, $`u_{n_j}(z,j)`$, required to make $`u(z)`$ continuous. However, the structure of this correction function at the end of the path, $`p=1`$, will depend on behavior of $`z(p)`$ near all the following singular and crossing points, $`z(p_k),j<kN`$. The resulting correction function, $`u_c(z)`$, obtained after passing all the points $`z(p_j),j=1\mathrm{}N`$, is given by the following expression:

$$u_c(z)=\underset{j=1}{\overset{N}{}}u_{n_j}(z,j).$$
(46)
In order to see, how these correction functions operate, consider values of the principal branch of $`u(z)`$ at the edges of the branch cut $`]\mathrm{\hspace{0.17em}1},+\mathrm{}[`$:

$$u(x\pm iϵ)=\frac{\pi ^2}{3}2\mathrm{L}\mathrm{i}_2(1/x)\frac{1}{2}\mathrm{ln}^2(x)\pm i\pi \mathrm{ln}(x).$$
(47)
In this formula, $`x>1`$ is real, and $`ϵ+0`$. Imagine that the branch cut is crossed from below ($``$). Then $`m_j=+1`$, and the value of the correction function $`u_1(z,j)`$, defined by Eq. (37), at the point $`z=x+iϵ`$ is equal to

$$u_1(x+iϵ,j)=2\pi i\mathrm{ln}(x).$$

If the branch cut is crossed from above ($``$), $`m_j=1`$, and the value of this correction function at the point $`z=xiϵ`$ is equal to

$$u_1(xiϵ,j)=+2\pi i\mathrm{ln}(x).$$

Thus, the correction function $`u_1(z,j)`$, added after the branch cut is crossed, eliminates the finite discontinuity of the principal branch, Eq. (45), of the function $`u(z)`$ along $`]\mathrm{\hspace{0.17em}1},+\mathrm{}[`$. The correction function $`u_2(z,j)`$, defined by Eq. (38), acts in a similar way at $`]\mathrm{\hspace{0.17em}0},1[`$.

Imagine now that the argument of $`u(z)`$ goes around the singular point at $`\mathrm{}`$, starting from $`z=x+iϵ`$, and coming back to $`z=xiϵ`$, without crossing the branch cut along the positive real axis. The value of $`u(z)`$ exhibits a singular change (from Eq. (45)) by

$$\mathrm{\Delta }u=2\pi i\mathrm{ln}(x).$$

In this case, $`\delta _j=2\pi `$, and the value of the correction function $`u_5(z,j)`$, defined by Eq. (41), at the point $`z=xiϵ`$ is equal to

$$u_5(xiϵ,j)=+2\pi i\mathrm{ln}(x).$$

Thus, the correction function $`u_5(z,j)`$, added after $`z`$ has moved near the singular point at $`\mathrm{}`$, eliminates the singular contribution to the value of the function $`u(z)`$. The correction functions $`u_3(z,j)`$ and $`u_4(z,j)`$, given by Eqs. (39) and (40), produce similar results for the other singular points.

If, in the above examples, the argument of $`u(z)`$ first crosses the branch cut, and then moves around the singular point, the correction function $`u_1(z,j)`$ has to be modified by adding nonzero $`U_j`$ to the logarithm according to Eq. (37).

The same principles of branch tracking apply to the function $`v(z)`$, defined by Eq. (15). First, we consider behavior of this function near its branch points 1, -1, and $`\mathrm{}`$:

$`v(z1)={\displaystyle \frac{1}{4}}\mathrm{ln}^2[(1z)/(1+z)]+v_{(1)}(z);`$
$`v(z1)={\displaystyle \frac{1}{4}}\mathrm{ln}^2[(1+z)/(1z)]+v_{(1)}(z);`$ (48)
$`v(z\mathrm{})={\displaystyle \frac{1}{2}}\mathrm{ln}(z^2/4)\mathrm{ln}[(z+1)/(z1)]+v_{(\mathrm{})}(z).`$
In these formulas, the functions with subscripts (1), (-1), and $`(\mathrm{})`$ are functions, analytic in the vicinities of 1, -1, and $`\mathrm{}`$, respectively.

Let us again consider a complex function $`z(p)`$, which can represent each of the arguments $`\gamma _k^{(j)}/\sigma `$ of the function $`v(z)`$ in Eq. (13). Let $`\{p_j\},j=1\mathrm{}N`$, denote values of the parameter $`p`$, such that $`z(p_j)`$ are singular points 1, -1, $`\mathrm{}`$, or $`z(p)`$ crosses the real axis. It is assumed that their real parts form an increasing set of numbers between 0 and 1. As before, each point $`p_j`$ is characterized by an index $`n_j`$, and either $`m_j`$ or $`\delta _j`$. Values of $`n_j`$ will be assigned below, and meanings of $`m_j`$ and $`\delta _j`$ remain the same.

Five correction functions, $`v_{n_j}(z,j)`$, are used to make the computed function $`v(z)`$ continuous.

If $`z(p)`$ crosses the branch cut $`]\mathrm{\hspace{0.17em}1},+\mathrm{}[`$ at $`z(p_j)`$,
let $`n_j=1`$, and

$$v_1(z,j)=+\pi ^2+m_j\pi i\{\mathrm{ln}[(1+z)/(1z)]+V_j\}.$$
(49)
If $`z(p)`$ crosses the branch cut $`]\mathrm{},1[`$ at $`z(p_j)`$,
let $`n_j=2`$, and

$$v_2(z,j)=\pi ^2m_j\pi i\{\mathrm{ln}[(1+z)/(1z)]+V_j\}.$$
(50)
If $`z(p)`$ moves near the singular point $`z(p_j)=1`$,
let $`n_j=3`$, and

$$v_3(z,j)=+\frac{\delta _j^2}{4}i\frac{\delta _j}{2}\{\mathrm{ln}[(1+z)/(1z)]+V_j\}.$$
(51)
If $`z(p)`$ moves near the singular point $`z(p_j)=1`$,
let $`n_j=4`$, and

$$v_4(z,j)=\frac{\delta _j^2}{4}i\frac{\delta _j}{2}\{\mathrm{ln}[(1+z)/(1z)]+V_j\}.$$
(52)
If $`z(p)`$ moves near the singular point $`z(p_j)=\mathrm{}`$,
let $`n_j=5`$, and

$$v_5(z,j)=i\delta _j\{\mathrm{ln}[(z+1)/(z1)]+\widehat{V}_j\}.$$
(53)
If $`z(p)`$ crosses the branch cut $`]1,1[`$ of the
function $`\mathrm{ln}[(z+1)/(z1)]`$ at $`z(p_j)`$, let $`n_j=6`$.

The additional terms, $`V_j`$ and $`\widehat{V}_j`$, necessary to correct behavior of the logarithms, are given by the following formulas:

$$V_j=\underset{\begin{array}{c}k>j\\ n_k=1,2\end{array}}{\overset{N}{}}2m_k\pi i+\underset{\begin{array}{c}k>j\\ n_k=3\end{array}}{\overset{N}{}}i\delta _k\underset{\begin{array}{c}k>j\\ n_k=4\end{array}}{\overset{N}{}}i\delta _k;$$
(54)
$$\widehat{V}_j=+\underset{\begin{array}{c}k>j\\ n_k=6\end{array}}{\overset{N}{}}2m_k\pi i+\underset{\begin{array}{c}k>j\\ n_k=3\end{array}}{\overset{N}{}}i\delta _k\underset{\begin{array}{c}k>j\\ n_k=4\end{array}}{\overset{N}{}}i\delta _k.$$
(55)
As in the previous case, a correction function $`v_{n_j}(z,j)`$ has to be added to the function $`v(z)`$ every time its argument $`z(p)`$ passes a singular or crossing point $`z(p_j)`$. This way, the calculated function $`v(z)`$ can be made continuous. However, the form of this correction function at the end of the path depends on behavior of $`z(p)`$ near all the points $`z(p_k)`$, following $`z(p_j)`$. The resulting correction function, $`v_c(z)`$, is the following:

$$v_c(z)=\underset{j=1}{\overset{N}{}}v_{n_j}(z,j).$$
(56)
Let us now briefly discuss the effect of using these correction functions. Consider values of $`v(z)`$ at the edges of the branch cut $`]\mathrm{\hspace{0.17em}1},+\mathrm{}[`$:

$`v(x\pm iϵ)={\displaystyle \frac{1}{2}}\mathrm{Li}_2[2/(1+x)]{\displaystyle \frac{1}{2}}\mathrm{Li}_2[2/(1x)]+`$
$`+{\displaystyle \frac{1}{2}}\mathrm{ln}^2[2/(1+x)]{\displaystyle \frac{1}{2}}\mathrm{ln}^2[2/(x1)]\pm `$ (57)
$`\pm i{\displaystyle \frac{\pi }{2}}\mathrm{ln}[(x1)/(x+1)].`$
Here, $`x>1`$ is real, and $`ϵ+0`$. Imagine that the branch cut is crossed from below ($``$). Then $`m_j=+1`$, and a value of the correction function, $`v_1(z,j)`$, defined by Eq. (47), at a point $`z=x+iϵ`$ is equal to

$$v_1(x+iϵ,j)=i\pi \mathrm{ln}[(x1)/(x+1)].$$

If the branch cut is crossed from above ($``$), then $`m_j=1`$, and a value of this correction function at $`z=xiϵ`$ is equal to

$$v_1(xiϵ,j)=+i\pi \mathrm{ln}[(x1)/(x+1)].$$

Therefore, the function $`v_1(z,j)`$, added to the function $`v(z)`$ after the branch cut is crossed, removes the discontinuity of the principal branch along $`]\mathrm{\hspace{0.17em}1},+\mathrm{}[`$. The correction function $`v_2(z,j)`$, given by Eq. (48), makes $`v(z)`$ continuous at $`]\mathrm{},1[`$.

Imagine now that the argument of $`v(z)`$ moves around the singular point +1, starting from $`z=x+iϵ`$ and returning to $`z=xiϵ`$, without crossing the branch cut. The value of $`v(z)`$ undergoes a change (from Eq. (55)) by

$$\mathrm{\Delta }v=i\pi \mathrm{ln}[(x1)/(x+1)].$$

Because $`\delta _j=2\pi `$, a value of the correction function $`v_3(z,j)`$, defined by Eq. (49), at $`z=xiϵ`$ is equal to

$$v_3(xiϵ,j)=+i\pi \mathrm{ln}[(x1)/(x+1)].$$

Thus, by adding the correction function $`v_3(z,j)`$, it is possible to eliminate the singular contribution to the value of $`v(z)`$, when $`z`$ goes around the singular point +1. The correction functions $`v_4(z,j)`$ and $`v_5(z,j)`$, given by Eqs. (50) and (51), produce the same results for the other two singular points.

If, in the above examples, the argument of $`v(z)`$ first crosses the branch cut, and then moves around the singular point, the correction function $`v_1(z,j)`$ should be modified by adding nonzero $`V_j`$ according to Eq. (47).

It is important to note that the function $`\sigma (z)`$ is also a multiple valued function. Its principal branch, defined by Eq. (27), changes sign each time the argument $`z`$ crosses the branch cut along the positive real axis. If this happens $`N`$ times while the parameter $`p`$ changes from 0 to 1, the corrected value $`\sigma _c(z)`$ of this function at $`p=1`$ is equal to

$$\sigma _c(z)=(1)^N\sigma (z),$$
(58)
where $`\sigma (z)`$ is the value of the principal branch of the complex square root.

We are now in a position to write a corrected expression for the generating integral:

$$I=\frac{16\pi ^3}{\sigma _c}[\underset{j=1}{\overset{4}{}}\underset{k=1}{\overset{4}{}}(v+v_c)(\gamma _k^{(j)}/\sigma _c)+\underset{j=2}{\overset{4}{}}(u+u_c)(\beta _1^{(1)}\beta _1^{(j)})]$$
(59)
In this formula, each of 16 terms with the function $`v(z)`$ contains its own correction function $`v_c(z)`$, which is characterized by its own $`N`$ and sets of numbers $`\{p_l\}`$, $`\{n_l\}`$, $`\{m_l\}`$, and $`\{\delta _l\}`$. The same is true for each of the three terms with the function $`u(z)`$.

Eq. (57) is profoundly different from the original formula, Eq. (13), for the generating integral. In Eq. (13), the functions $`u(z)`$, $`v(z)`$, and $`\sigma (z)`$ are expressed in terms of the multiple valued logarithm and square root. When Eq. (57) is used, it is assumed, on the contrary, that all the logarithms and square roots are represented by their principal branches, and, therefore, can be readily evaluated by a computer. The multiple valued nature of the functions $`u(z)`$, $`v(z)`$, and $`\sigma (z)`$ is taken into account explicitly by means of the additional correction terms and factors. Also, singular contributions from different terms in Eq. (13) are expected to cancel each other to yield a correct value for the generating integral, Eq. (11). When we use Eq. (57), all the singular contributions from different terms are cancelled explicitly and separately, so that each function $`(v+v_c)(\gamma _k^{(j)}/\sigma _c)`$ or $`(u+u_c)(\beta _1^{(1)}\beta _1^{(j)})`$ is continuous along the path from (1,1,1,1,1,1) to $`(\alpha _{12},\alpha _{13},\alpha _{14},\alpha _{23},\alpha _{24},\alpha _{34})`$. As a result, the generating integral, obtained from Eq. (57), is a continuous function of the complex parameters $`\{\alpha _{jk}\}`$. Thus, the problem of branch tracking is successfully solved in the most general case.

Let us consider an example of branch tracking for the function $`v(z)`$ with the argument $`z(p)`$. A sample path in the complex $`p`$ plane is displayed in Fig. 3a, and the corresponding path from $`z(0)=A`$ to $`z(1)=B`$ in the complex $`z`$ plane is shown in Fig. 3b. There are five points of interest: $`z(p_{\mathrm{\hspace{0.17em}1}})=1`$, $`z(p_{\mathrm{\hspace{0.17em}3}})=\mathrm{}`$, and $`z(p_{\mathrm{\hspace{0.17em}5}})=1`$ are singular points of the function $`v(z)`$; $`z(p_{\mathrm{\hspace{0.17em}2}})`$ and $`z(p_{\mathrm{\hspace{0.17em}4}})`$ are points, where the argument $`z(p)`$ crosses the branch cuts. According to the chosen classification: $`n_1=3,\delta _1=\pi `$; $`n_2=1,m_2=+1`$; $`n_3=5,\delta _3=\pi `$; $`n_4=2,m_4=1`$; $`n_5=4,\delta _5=\pi `$. Using Eqs. (47)-(53), one can easily obtain expressions for the correction functions:

$`v_3(z,1)={\displaystyle \frac{\pi ^2}{4}}i{\displaystyle \frac{\pi }{2}}\{\mathrm{ln}[(1+z)/(1z)]i\pi \};`$
$`v_1(z,2)=\pi ^2+i\pi \{\mathrm{ln}[(1+z)/(1z)]+i\pi \};`$
$`v_5(z,3)=0i\pi \{\mathrm{ln}[(z+1)/(z1)]i\pi \};`$
$`v_2(z,4)=\pi ^2+i\pi \{\mathrm{ln}[(1+z)/(1z)]i\pi \};`$
$`v_4(z,5)={\displaystyle \frac{\pi ^2}{4}}i{\displaystyle \frac{\pi }{2}}\{\mathrm{ln}[(1+z)/(1z)]+0\}.`$
The resulting correction function, $`v_c(z)`$, is

$`v_c(z)={\displaystyle \frac{3\pi ^2}{2}}+i\pi \mathrm{ln}[(1+z)/(1z)]`$
$`i\pi \mathrm{ln}[(z+1)/(z1)].`$
One can see that this function is different from zero, even if $`A=B`$, i.e. the contour is closed. This is not surprising. Even though $`v(z)`$ is represented by its principal branch, the sum $`v(z)+v_c(z)`$ is still a multiple valued function. Its value generally undergoes a finite change, if its argument traverses a closed loop, encircling branch points. Consider a value of $`v_c(z)`$ at $`z=xiϵ`$, where $`1<x<1`$, and $`ϵ+0`$. In this case, the logarithms cancel, and

$$v_c(xiϵ)=\pi ^2/2.$$

Contributions of this type from different terms in Eq. (57) produce an additional constant $`m\pi ^2`$, needed to correct a value of the generating integral in the case of real parameters $`\{\alpha _{jk}\}`$. Thus, Eq. (34) is a particular case of Eq. (57).

This example demonstrates that the branch tracking in the general case requires a comprehensive numerical analysis of behavior of all the arguments in Eq. (13).

### II.6 Numerical procedure and results

Practical implementation of the method, described in the previous sections, is inevitably a very complicated task. Detailed information about the recursive procedure, needed to compute the family of integrals, Eq. (10), can be found in hil87 . Here we describe only the procedure for numerical branch tracking.

First, the set of points, $`\{p_k\}`$, at which different terms in Eq. (13) can exhibit singular behavior, is determined. This is done by solving the sixth-order equation $`\sigma ^2=0`$, and linear equations of Eq. (30), with the parametrization according to Eq. (35). Only those values of $`p`$, that lie in or near the real interval $`]\mathrm{\hspace{0.17em}0},1[`$, are included in the set $`\{p_k\}`$. Then, a path from 0 to 1 in the complex $`p`$ plane is chosen. Fig. 2 gives an idea of this. The whole path is shifted downward by a small imaginary quantity $`iϵ`$ to avoid possible ambiguities, when arguments of the functions $`u(z)`$ and $`v(z)`$ are real. All the arguments in Eq. (57) are computed at the final point of this path, $`p=1iϵ`$. In actual calculations, $`ϵ`$ was set to $`10^{28}`$. This did not affect values of the integrals, but was enough to shift the arguments from the real axis.

The path in the complex $`p`$ plane is divided into small intervals, as shown in Fig. 4. The intervals $`]P_l,P_{l+1}[`$, into which the linear segments between the singular points, $`\{p_k\}`$, are divided, have a typical length of $`10^2`$. Each small semicircle beneath a singular point, $`p_k`$, has a radius $`r=10^9`$, and divided into six parts. The corresponding boundary points are

$$P_l=p_k+r\mathrm{exp}[i(\pi l/6)]iϵ,l=0\mathrm{}6.$$
(60)
In order to obtain full information about behavior of different arguments in Eq. (13), all these arguments should be computed at all the points $`P_l`$ along the path. The quantities $`\gamma _k^{(j)}`$, $`\beta _k^{(j)}`$, and $`\sigma ^2`$, given by Eqs. (16)-(20), are simple functions of $`\{\alpha _{jk}(p)\}`$, so this calculation can be performed almost immediately. The values of each argument are analyzed, and the numbers $`N`$, $`n_j`$, $`m_j`$, and $`\delta _j`$, $`j=1\mathrm{}N`$, needed to apply the formulas of Sec. II.E, are determined. This procedure works as follows. To find out, if an argument crosses the real axis, the imaginary parts of its values, computed at points $`P_l`$ and $`P_{l+1}`$, are compared. If they have opposite signs, dichotomy is used to reduce the interval, and determine, where the real axis is crossed, and in which direction. This is also done for the intervals along each small semicircle, but without the dichotomy. In this manner, all the crossing points can be found and analyzed. At each interval, $`]P_l,P_{l+1}[`$, this analysis is carried out for $`\sigma ^2`$ first. If $`\sigma ^2`$ crosses the positive real axis, when $`p`$ is between $`P_l`$ and $`P_{l+1}`$, the sign of $`\sigma `$ should be changed. Then, all $`\gamma _k^{(j)}`$, computed at $`P_{l+1}`$, should be multiplied by (-1), and all $`\beta _k^{(j)}`$ – inverted. The quantity $`\delta _j`$ is determined as follows:

$$\delta _j=\underset{l=0}{\overset{5}{}}[\mathrm{arg}(z(P_{l+1}))\mathrm{arg}(z(P_l))].$$
(61)
Here, the points $`P_l`$ are specified by Eq. (58), and $`z`$ can stand for any of the arguments $`\gamma _k^{(j)}/\sigma `$ and $`\beta _1^{(1)}\beta _1^{(j)}`$ as functions of $`\{\alpha _{jk}(p)\}`$. The right hand side of this formula is presented as a sum, because as $`p`$ changes from $`P_{\mathrm{\hspace{0.17em}0}}`$ to $`P_{\mathrm{\hspace{0.17em}6}}`$, the argument $`z(p)`$ may go around the singular point, $`z(p_k)`$ several times. Thus, behavior of each argument in the vicinity of each singular point can be analyzed. The described procedure of numerical branch tracking provides all the information, necessary for successful use of Eq. (57).

The entire algorithm for analytic evaluation of the four-particle integrals with complex parameters was tested in four different ways.

First, real parameters $`\{\alpha _{jk}\}`$ were used, and the results, obtained using the complex algorithm of Sec. II.E, were compared with results, provided by the method of Sec. II.D for the real case. The real parts of the computed integrals, Eq. (10), were invariably in excellent agreement. The imaginary parts, given by the new method, were at least 20 orders of magnitude smaller than the real parts, and could be considered negligible. Therefore, the complex algorithm works correctly for any acceptable real parameters.

Second, the parameters $`\{\alpha _{jk}\}`$ were multiplied by an arbitrarily chosen complex number $`\lambda `$. A resulting integral, Eq. (10), with a particular set $`\{n_{jk}\}`$ must be equal to its original value, multiplied by $`f=1/(\lambda )^K`$, where $`K=n_{12}+n_{13}+n_{14}+n_{23}+n_{24}+n_{34}+3`$. Various $`\lambda `$’s were used, and values of the integrals, calculated directly, were compared with the rescaled original values. Remarkable agreement was observed in all these cases. Note that different values of $`\lambda `$ correspond to different paths in the space of parameters according to Eq. (35).

Third, if two exponential parameters, $`\alpha _{13}`$ and $`\alpha _{24}`$, are equal to zero, the six Coulomb integrals and one overlap integral, needed to determine matrix elements of the Hamiltonian according to Eqs. (4)-(9), can be obtained analytically in terms of rational functions and logarithms. Values of these integrals, calculated with various sets of complex parameters, $`\alpha _{12}`$, $`\alpha _{14}`$, $`\alpha _{23}`$, $`\alpha _{34}`$, were compared with the same integrals, computed using the new method. They were always in complete agreement.

Fourth, different paths in the complex $`p`$ plane were chosen. They included singularities, located not only near the real axis, but also further away. The results did not depend on the choice of the path. This fact suggests that the described method of numerical branch tracking is stable and reliable. Of course, the path in actual calculations should be as simple as possible, provided that all nearby singularities are carefully taken into account.

Table I displays values of the integrals for three different sets of parameters $`\{\alpha _{jk}\}`$, used to test the computer program. Many other sets of parameters were also considered. All the integrals were calculated using the general algorithm for numerical branch tracking, described in Sec. II.E. Implementation of this algorithm requires quadruple precision. The program computes a family of 64 integrals, Eq. (10), with two possible values for every index: $`n_{jk}=0,1`$. Only seven integrals, necessary to obtain matrix elements of the Hamiltonian according to Eqs. (4)-(9), are presented in Table I for each set of parameters.

Our results demonstrate that the developed algorithm allows precise evaluation of the four-particle integrals with arbitrary complex parameters, provided that the integrals themselves converge.

The described method makes it possible to use the highly versatile exponential-trigonometric basis functions in variational calculations of four-particle Coulomb systems. In order to illustrate efficiency of the new basis, we would like to mention some results, obtained previously zot00  for the following systems: $`e^+e^{}e^+e^{}`$, $`p^+\mu ^{}p^+\mu ^{}`$, $`\mu ^+e^{}\mu ^+e^{}`$, and $`p^+e^{}p^+e^{}`$. The calculations were performed using one exponential-trigonometric basis function:

$$\mathrm{\Psi }=\widehat{S}\mathrm{exp}(\underset{j<k}{\overset{4}{}}A_{jk}r_{jk})\mathrm{sin}(\underset{j<k}{\overset{4}{}}B_{jk}r_{jk}+C).$$
(62)
This function includes 12 real nonlinear parameters, $`\{A_{jk}\}`$ and $`\{B_{jk}\}`$, and one linear parameter, $`\mathrm{tan}(C)`$. It can be considered a linear combination of two exponential functions, Eq. (2), with the complex parameters $`A_{jk}\pm iB_{jk}`$. The operator $`\widehat{S}`$ ensures that this function has correct symmetry with respect to permutations of particles.

All integrals, necessary to determine matrix elements of the Hamiltonian, Eq. (1), with the function $`\mathrm{\Psi }`$, were computed according to the method, described in this paper. The nonlinear parameters were subjected to careful gradient optimization. For more details about this calculation, see zot00 .

Table II exhibits values of the ground-state energy, $`E`$, for $`e^+e^{}e^+e^{}`$, $`p^+\mu ^{}p^+\mu ^{}`$, $`\mu ^+e^{}\mu ^+e^{}`$, and $`p^+e^{}p^+e^{}`$, determined using the variational method with the trial function $`\mathrm{\Psi }`$. The table also displays the most accurate energy values, $`E_0`$, available for these systems usu98 ; fro97 ; reb97 ; kin99 . One can see that the relative errors are 0.2%, 0.7%, 2.4%, and 3.6%, respectively. The results for two adiabatic systems, $`\mu ^+e^{}\mu ^+e^{}`$ and $`p^+e^{}p^+e^{}`$, with very low mass ratios, $`m/M`$, are very impressive. Neither Gaussian, nor exponential functions are even nearly as efficient zot00 . Thus, a single symmetrized exponential-trigonometric basis function, Eq. (60), provides a remarkable accuracy in variational calculations of various four-particle systems.

## III Conclusion

The method for analytic evaluation of four-particle integrals with complex parameters, described in this paper, can be regarded as both further theoretical development and practical implementation of the original method by Fromm and Hill hil87 . Validity of this method is not limited to the case of real parameters. Moreover, because the integrals are expressed in terms of multiple valued complex functions, it is more natural to consider a general case when all the parameters are complex. The original formula, Eq. (13), for the generating integral can be used only in the immediate vicinity of the standard reference point where all the parameters are equal to 1. The procedure of numerical branch tracking, proposed in this paper, allows computation of the integrals at any other point in the space of six complex parameters, by taking into account all branch changes along the path. The simplified method of branch tracking for real parameters is also discussed.

The new method makes possible high precision variational solution of the Coulomb four-body problem in the basis of exponential-trigonometric functions. The first calculations have shown high efficiency of this basis zot00 . They have also demonstrated correctness of the branch tracking algorithm, described in this paper.

However, if the full potential of the exponential-trigonometric basis is to be revealed, an efficient procedure for selecting optimal values of the nonlinear parameters is necessary. Ideally, all the parameters should be chosen *a priori*, and all matrix elements – computed only once. Such a procedure has been developed by the authors for the case of adiabatic three-particle systems zot94 . All nonlinear parameters of the exponential-trigonometric functions had been chosen *before* the computation, which yielded 10 correct significant figures for the ground state energy of $`H_2^+`$ zot94 . We believe that the exponential-trigonometric basis can provide similar precision in calculations of four-particle systems.