text:
# Quantum three-body system in 𝐷 dimensions

## I Introduction

From the very early stage of the progress in quantum mechanics in the real three-dimensional world, it has been pointed out that the essence of these theories would be easily understandable if their mathematics is constructed in the non-relativistic hyperspace worlds . The mathematical tools for generalization of the orbital angular momentum in an arbitrary $`D`$-dimensional space have been presented . Recently, the $`D`$-dimensional Coulombic and the harmonic oscillator problems in a two-body system have been studied in some detail by many authors .

Exact solutions played very important roles in the development of physics. The exact solutions of the Schrödinger equation in the real three dimensional space for a hydrogen atom and for a harmonic oscillator were important technical achievements in quantum mechanics , which provided strong evidence in favor of the theory being correct, at least as far as atomic physics is concerned. The next simplest atom is the helium atom, for which the Schrödinger equation cannot be solved analytically, but only numerically . In the numerical calculation, one of the main difficulties is how to separate the global rotational degrees of freedom.

In our previous paper  we separated completely the global rotational degrees of freedom in the Schrödinger equation for an $`N`$-body system in the real three-dimensional space from the internal ones. We have determined a complete set of $`(2l+1)`$ independent base functions for a given total orbital angular momentum $`l`$, which are the homogeneous polynomials in the components of coordinate vectors and do not contain the Euler angles explicitly. Any function with the given angular momentum $`l`$ in the system can be expanded with respect to the base functions, where the coefficients are the functions of the internal variables, called the generalized radial functions. The generalized radial equations satisfied by the functions are established explicitly . For the typical three-body system in the real three dimensional space , such as a helium atom  and a positronium negative ion , the generalized radial equations  have been solved numerically with high precision.

With the interest of higher dimensional field theory recently, we have a try to generalize the study of the $`D`$-dimensional two-body system to the $`D`$-dimensional three-body system. The purpose of this paper is, for a three-body system in an arbitrary $`D`$-dimensional space, to find a complete set of independent base functions with any given total orbital angular momentum and to reduce the Schrödinger equation with a spherically symmetric potential $`V`$ to the generalized radial equations, where only three internal variables are involved. Any function with the given angular momentum in the system can be expanded with respect to the base functions. It provides a possibility to calculate numerically the energy levels of the three-body system in $`D`$-dimensions with high precision.

From the viewpoint of mathematics, the separation of the global rotational degrees of freedom from the internal ones is a typical application of group theory to physics. The properties of the independent base functions for a given total orbital angular momentum would be more clear if they are constructed in arbitrary $`D`$-dimensional space than that in the real three dimensional space. The total orbital angular momentum for a three-body system in a $`D`$-dimensional space is described by an irreducible representation denoted by a Young diagram with one or two rows. For the real three-dimensional space, the rotational symmetry group is SO(3) group, and its only irreducible representations denoted by the Young diagrams with two rows are $`[l,1]`$, which are equivalent to the representations denoted by the one-row Young diagrams $`[l,0]`$, respectively. This is the reason why the angular momentum can be described by only one quantum number $`l`$ for the real three-dimensional space.

This paper is organized as follows. After separating the motion of the center of mass by the Jacobi coordinate vectors in Sec. 2, we review in Sec. 3 the generalization of the orbital angular momentum operators and the properties of the spherical harmonics  and the harmonic polynomials  for a two-body system in $`D`$ dimensions. In Sec. 4 we will define the generalized harmonic polynomials for a three-body system in $`D`$ dimensions and prove that they constitute a complete set of independent base functions for a given total orbital angular momentum in the system. The generalized radial functions are defined and the generalized radial equations are derived in Sec. 5. Some conclusions will be given in Sec. 6.

## II Schrödinger equation in $`D`$ dimensions

For a quantum $`N`$-body system in an arbitrary $`D`$-dimensional space, we denote the position vectors and the masses of $`N`$ particles by $`𝐫_k`$ and by $`m_k`$, $`k=1,2,\mathrm{},N`$, respectively. $`M=_km_k`$ is the total mass. The Schrödinger equation for the $`N`$-body system with a pair potential $`V`$, depending upon the distance of each pair of particles, $`|𝐫_j𝐫_k|`$, is

$$\frac{1}{2}\underset{k=1}{\overset{N}{}}m_k^1_{𝐫_k}^2\mathrm{\Psi }+V\mathrm{\Psi }=E\mathrm{\Psi },$$
$`(1)`$
where $`_{𝐫_k}^2`$ is the Laplace operator with respect to the position vector $`𝐫_k`$. For simplicity, the natural units $`\mathrm{}=c=1`$ are employed throughout this paper. The total orbital angular momentum operators $`L_{ab}`$ in $`D`$ dimensions are defined as

$$L_{ab}=L_{ba}=i\underset{k=1}{\overset{N}{}}\left\{r_{ka}\frac{}{r_{kb}}r_{kb}\frac{}{r_{ka}}\right\},a,b=1,2,\mathrm{}D,$$
$`(2)`$
where $`r_{ka}`$ denotes the $`a`$th component of the position vector $`𝐫_k`$.

Now, we replace the position vectors $`𝐫_k`$ by the Jacobi coordinate vectors $`𝐑_j`$:

$$𝐑_0=M^{1/2}\underset{k=1}{\overset{N}{}}m_k𝐫_k,𝐑_j=\left(\frac{m_{j+1}M_j}{M_{j+1}}\right)^{1/2}\left(𝐫_{j+1}\underset{k=1}{\overset{j}{}}\frac{m_k𝐫_k}{M_j}\right),$$

$$1j(N1),M_j=\underset{k=1}{\overset{j}{}}m_k,M_N=M,$$
$`(3)`$
where $`𝐑_0`$ describes the position of the center of mass, $`𝐑_1`$ describes the mass-weighted separation from the second particle to the first particle. $`𝐑_2`$ describes the mass-weighted separation from the third particle to the center of mass of the first two particles, and so on. An additional factor $`\sqrt{M}`$ is included in $`𝐑_j`$ for convenience. The mass-weighted factors in front of the formulas for $`𝐑_j`$ are determined by the condition

$$\underset{k=1}{\overset{N}{}}m_k𝐫_k^2=\underset{j=0}{\overset{N1}{}}𝐑_j^2,$$

One may determine the factors one by one from the following schemes. In the center-of-mass frame, if the first $`j`$ particles coincide with each other and the last $`(Nj1)`$ particles are located at the origin, the factor in front of $`𝐑_j`$ is determined by

$$𝐫_1=𝐫_2=\mathrm{}=𝐫_j=m_{j+1}𝐫_{j+1}/M_j,\underset{k=1}{\overset{j+1}{}}m_k𝐫_k^2=𝐑_j^2.$$
$`(4)`$
A straightforward calculation by replacement of variables shows that the Laplace operator in Eq. (1) and the total orbital angular momentum operator $`L_{ab}`$ in Eq. (2) are directly expressed in $`𝐑_j`$:

$$\begin{array}{c}^2=\underset{k=1}{\overset{N}{}}m_k^1_{𝐫_k}^2=\underset{j=0}{\overset{N1}{}}_{𝐑_j}^2,\hfill \\ L_{ab}=i\underset{j=0}{\overset{N1}{}}\left\{R_{ja}\frac{}{R_{jb}}R_{jb}\frac{}{R_{ja}}\right\}.\hfill \end{array}$$
$`(5)`$
In the center-of-mass frame, $`𝐑_0=0`$. The Laplace operator (5) obviously has the symmetry of the O$`(NDD)`$ group with respect to $`(N1)D`$ components of $`(N1)`$ Jacobi coordinate vectors. The O$`(NDD)`$ group contains a subgroup SO($`D)\times `$O$`(N1)`$, where SO$`(D)`$ is the rotation group in the $`D`$-dimensional space. The space inversion and the different definitions for the Jacobi coordinate vectors in the so-called Jacobi tree  can be obtained by O$`(N1)`$ transformations. For the system of identical particles, the permutation group among particles is also a subgroup of the O$`(N1)`$ group .

It is easy to obtain the inverse transformation of Eq. (3):

$$𝐫_j=\left[\frac{M_{j1}}{m_jM_j}\right]^{1/2}𝐑_{j1}\underset{k=j}{\overset{N1}{}}\left[\frac{m_{k+1}}{M_kM_{k+1}}\right]^{1/2}𝐑_k+M^{1/2}𝐑_0,$$

$$𝐫_j𝐫_k=\left[\frac{M_j}{m_jM_{j1}}\right]^{1/2}𝐑_{j1}+\underset{i=k}{\overset{j2}{}}\left[\frac{m_{i+1}}{M_iM_{i+1}}\right]^{1/2}𝐑_i\left[\frac{M_{k1}}{m_kM_k}\right]^{1/2}𝐑_{k1}.$$
$`(6)`$
Thus, the potential $`V`$ is a function of $`𝐑_j𝐑_k`$ and is rotationally invariant.

## III Harmonic polynomials in $`D`$ dimensions

In the center-of-mass frame, $`𝐑_0=0`$. Hence, for a two-body system there is only one Jacobi coordinate vector $`𝐑_1`$, which will be denoted by $`𝐱`$ for simplicity:

$$\begin{array}{c}𝐱=\left(\frac{m_1m_2}{m_1+m_2}\right)^{1/2}\left\{𝐫_2𝐫_1\right\},\\ ^2=_𝐱^2,L_{ab}=i\{x_a\frac{}{x_b}x_b\frac{}{x_a}\},\end{array}$$
$`(7)`$
Louck  introduced the hyperspherical coordinates

$$\begin{array}{c}x_1=r\mathrm{cos}\theta _1\mathrm{sin}\theta _2\mathrm{}\mathrm{sin}\theta _{D1},\hfill \\ x_2=r\mathrm{sin}\theta _1\mathrm{sin}\theta _2\mathrm{}\mathrm{sin}\theta _{D1},\hfill \\ x_k=r\mathrm{cos}\theta _{k1}\mathrm{sin}\theta _k\mathrm{}\mathrm{sin}\theta _{D1},3kD1,\hfill \\ x_D=r\mathrm{cos}\theta _{D1}.\hfill \end{array}$$
$`(8)`$
The spherical harmonics $`Y_{l_{D2},\mathrm{},l_1}^l`$ in $`D`$ dimensions  are the simultaneous eigenfunctions of the commutant operators $`𝐋_k^2`$:

$$\begin{array}{c}𝐋_1^2=\frac{^2}{\theta _1^2},𝐋_k^2=\left\{\frac{1}{\mathrm{sin}^{k1}\theta _k}\frac{}{\theta _k}\mathrm{sin}^{k1}\theta _k\frac{}{\theta _k}\frac{𝐋_{k1}^2}{\mathrm{sin}^2\theta _k}\right\},\hfill \end{array}$$
$`(9)`$
$$\begin{array}{c}𝐋_1^2Y_{l_{D2},\mathrm{},l_1}^l(\theta _1\mathrm{}\theta _{D1})=l_1^2Y_{l_{D2},\mathrm{},l_1}^l(\theta _1\mathrm{}\theta _{D1}),\hfill \\ 𝐋_k^2Y_{l_{D2},\mathrm{},l_1}^l(\theta _1\mathrm{}\theta _{D1})=l_k(l_k+k1)Y_{l_{D2},\mathrm{},l_1}^l(\theta _1\mathrm{}\theta _{D1}),\hfill \\ ll_{D1}=0,1,\mathrm{},l_k=0,1,\mathrm{},l_{k+1},l_1=l_2,l_2+1,\mathrm{},l_21,l_2,\hfill \end{array}$$
$`(10)`$
where $`𝐋^2𝐋_{D1}^2`$, $`0r<\mathrm{}`$, $`\pi \theta _1\pi `$, $`0\theta _k\pi `$, and $`2kD1`$. The volume element of the configuration space is

$$\underset{j=1}{\overset{D}{}}dx_j=r^{D1}dr\underset{j=1}{\overset{D1}{}}\left(\mathrm{sin}\theta _j\right)^{j1}d\theta _j.$$
$`(11)`$
Through a direct calculation by replacement of variables, one obtains

$$_𝐱^2=\frac{1}{r^{D1}}\frac{}{r}r^{D1}\frac{}{r}\frac{𝐋^2}{r^2},$$
$`(12)`$
Due to the spherical symmetry, the wave function can be expressed as

$$\psi _{l_{D2},\mathrm{},l_1}^l(𝐱)=\varphi _l(r)Y_{l_{D2},\mathrm{},l_1}^l(\theta _1\mathrm{}\theta _{D1}),$$
$`(13)`$
and the $`D`$-dimensional Schrödinger equation (1) for a two-body system in the center-of-mass frame reduces to the radial equation

$$\frac{1}{r^{D1}}\frac{}{r}r^{D1}\frac{}{r}\varphi _l(r)\frac{l(l+D2)}{r^2}\varphi _l(r)=2\left[EV(r)\right]\varphi _l(r).$$
$`(14)`$
Bander and Itzykson  introduced the harmonic polynomials in $`D`$ dimensions

$$𝒴_{l_{D2},\mathrm{},l_1}^l(𝐱)=r^lY_{l_{D2},\mathrm{},l_1}^l(\theta _1\mathrm{}\theta _{D1})r^lY_{l_{D2},\mathrm{},l_1}^l(\widehat{𝐱}),$$
$`(15)`$
to avoid the angular functions $`Y_{l_{D2},\mathrm{},l_1}^l(\theta _1\mathrm{}\theta _{D1})`$. $`𝒴_{l_{D2},\mathrm{},l_1}^l(𝐱)`$ is a homogeneous polynomial of degree $`l`$ in the components of $`𝐱`$ and satisfies the Laplace equation

$$_𝐱^2𝒴_{l_{D2},\mathrm{},l_1}^l(𝐱)=0.$$
$`(16)`$
The number of linearly independent homogeneous polynomials of degree $`l`$ in $`D`$ components of $`𝐱`$ is $`N(l)=(l+D1)!/l!(D1)!`$. The Laplace equation (16) gives $`N(l2)=(l+D3)!/(l2)!(D1)!`$ constraints. Hence, the number of the harmonic polynomials $`𝒴_{l_{D2},\mathrm{},l_1}^l(𝐱)`$ of degree $`l`$ as well as the number of the spherical harmonics $`Y_{l_{D2},\mathrm{},l_1}^l(\widehat{𝐱})`$ in $`D`$ dimensions is

$$N(l)N(l2)=\frac{(2l+D2)(l+D3)!}{l!(D2)!}=d_D([l,0,\mathrm{},0]).$$
$`(17)`$
$`d_D([l,0,\mathrm{},0])`$ is the dimension of the irreducible representation of SO$`(D)`$ denoted by the one-row Young diagram $`[l,0,\mathrm{},0]`$. $`[l,0,\mathrm{},0]`$ describes the symmetric traceless tensor representation. In fact, any polynomial in the components of one vector $`𝐱`$ has to belong to a symmetric representation.

Due to the spherical symmetry, one only needs to write the explicit form of the highest weight state

$$𝒴_{l,\mathrm{},l}^l(𝐱)=N_l(x_1+ix_2)^l.$$
$`(18)`$
where $`N_l`$ denotes the normalization factor. The partners of $`𝒴_{l,\mathrm{},l}^l(𝐱)`$ can be simply generated by rotation. Now, the solution to the Schrödinger equation in the center-of-mass frame can be re-expressed as

$$\psi _{l,\mathrm{},l}^l(𝐱)=R_l(r)𝒴_{l,\mathrm{},l}^l(𝐱),$$
$`(19)`$
and the radial equation is easy to be derived:

$$\frac{1}{r^{D1}}\frac{}{r}r^{D1}\frac{}{r}R_l(r)+\frac{2l}{r}\frac{}{r}R_l(r)=2\left[EV(r)\right]R_l(r).$$
$`(20)`$
Recall $`R_l(r)=r^l\varphi _l(r)`$. Eq. (20) coincides with Eq. (14) but the angle variables do not appear explicitly in calculation.

The number (17) of the harmonic polynomials $`𝒴_{l_{D2},\mathrm{},l_1}^l(𝐱)`$ of degree $`l`$ can be understood from another viewpoint. After removing those homogeneous polynomials in the form $`r^2f(𝐱)`$, where $`f(𝐱)`$ is a homogeneous polynomial of degree $`(l2)`$, Eq. (17) shows the number of the remaining linearly independent homogeneous polynomials of degree $`l`$ in the components of $`𝐱`$. Therefore, the harmonic polynomials $`𝒴_{l_{D2},\mathrm{},l_1}^l(𝐱)`$ construct a complete set of linearly independent base functions for the homogeneous polynomials of degree $`l`$ in the components of $`𝐱`$, excluded those in the form of $`r^2f(𝐱)`$.

## IV Three-body system in $`D`$-dimensions

For a three-body system, in the center-of-mass frame there are two Jacobi coordinate vectors $`𝐑_1`$ and $`𝐑_2`$, which will be denoted by $`𝐱`$ and $`𝐲`$, respectively:

$$\begin{array}{c}𝐱=\left[\frac{m_1m_2}{m_1+m_2}\right]^{1/2}\left\{𝐫_2𝐫_1\right\},𝐲=\left[\frac{(m_1+m_2)m_3}{m_1+m_2+m_3}\right]^{1/2}\left\{𝐫_3\frac{m_1𝐫_1+m_2𝐫_2}{m_1+m_2}\right\},\hfill \end{array}$$

$$\begin{array}{c}^2=_𝐱^2+_𝐲^2,\hfill \\ L_{ab}=L_{ab}^{(x)}+L_{ab}^{(y)}=i\left\{x_a\frac{}{x_b}x_b\frac{}{x_a}\right\}i\left\{y_a\frac{}{y_b}y_b\frac{}{y_a}\right\},\hfill \end{array}$$
$`(21)`$
The Schrödinger equation (1) reduces to

$$\begin{array}{c}\{_𝐱^2+_𝐲^2\}\mathrm{\Psi }(𝐱,𝐲)=2\{EV(\xi _1,\xi _2,\xi _3)\}\mathrm{\Psi }(𝐱,𝐲),\\ \xi _1=𝐱𝐱,\xi _2=𝐲𝐲,\xi _3=𝐱𝐲.\end{array}$$
$`(22)`$
where $`\xi _j`$ are the internal variables. Since Eq. (22) is rotational invariant, the total orbital angular momentum is conserved. The wave function $`\mathrm{\Psi }(𝐱,𝐲)`$ with the given total angular momentum has to belong to an irreducible representation of SO($`D`$). In the traditional method, one calculates the wave function by the Clebsch-Gordan coefficients:

$$\begin{array}{c}\underset{l_{D2},\mathrm{},l_1l_{D2}^{},\mathrm{},l_1^{}}{}𝒴_{l_{D2},\mathrm{},l_1}^l(𝐱)𝒴_{l_{D2}^{},\mathrm{},l_1^{}}^l^{}(𝐲)l,l_{D2},\mathrm{},l_1;l^{},l_{D2}^{},\mathrm{},l_1|L,M.\hfill \end{array}$$
$`(23)`$
As usual, $`𝒴_{l_{D2},\mathrm{},l_1}^l(𝐱)`$ and $`𝒴_{l_{D2}^{},\mathrm{},l_1^{}}^l^{}(𝐲)`$ are called the partial angular momentum states, and their combination is called the total angular momentum state, which is a homogeneous polynomial of degrees $`l`$ and $`l^{}`$ in the components of $`𝐱`$ and $`𝐲`$, respectively.

There are three problems. First, what kinds of representations (or total angular momentum $`L`$) appear in the Clebsch-Gordan series for decomposition of the direct product of two representations denoted by one-row Young diagrams $`[l,0,\mathrm{},0]`$ and $`[l^{},0,\mathrm{},0]`$? This problem has been solved in group theory by the Littlewood-Richardson rule and traceless conditions. A new character is that the representations denoted by two-row Young diagrams appear in the Clebsch-Gordan series for a three-body system when $`D>3`$. Those representations denoted by the Young diagrams with more than two rows could not appear because there are only two Jacobi coordinate vectors. For simplicity we denote a one-row or two-row Young diagram by $`[\mu ,\nu ][\mu ,\nu ,0,\mathrm{},0]`$. Hence, we have the Clebsch-Gordan series:

$$[l,0][l^{},0]\underset{s=0}{\overset{n}{}}\underset{t=0}{\overset{ns}{}}[l+l^{}s2t,s],$$
$`(24)`$
where $`n`$ is the minimum between $`l`$ and $`l^{}`$. The representations with $`t=0`$ are calculated by the Littlewood-Richardson rule, and the remaining are calculated by the traceless conditions. The dimension of a representation denoted by a two-row Young diagram is

$$\begin{array}{cc}\hfill d_D([\mu ,\nu ])& =(D+2\mu 2)(D+\mu +\nu 3)(D+2\nu 4)(\mu \nu +1)\hfill \\ & \times \frac{(D+\mu 4)!(D+\nu 5)!}{(\mu +1)!\nu !(D2)!(D4)!}.\hfill \end{array}$$
$`(25)`$
When $`D=4`$, the representation denoted by a two-row Young diagram reduces to a direct sum of a selfdual representation $`[(S)\mu ,\nu ]`$ and an antiselfdual one $`[(A)\mu ,\nu ]`$. Their dimensions are equal to half of $`d_4([\mu ,\nu ])`$ given in Eq. (25). When $`D=3`$, due to the traceless condition, the only representations with the two-row Young diagrams are representations $`[\mu ,1]`$, which are equivalent to that with the one-row Young diagrams $`[\mu ,0]`$, respectively. Eq. (25) still holds for $`D=3`$. Second, how to calculate the Clebsch-Gordan coefficients? The calculation must be very complicated. We will avoid the difficulty by the method of determining the highest weight states directly. Third, how many base functions are independent for a given total orbital angular momentum such that any wave function with the same angular momentum can be expanded with respect to the base functions where the coefficients are the functions of the internal variables. We are going to solve the last two problems by group theory.

Let us sketch some necessary knowledge of group theory. From the representation theory of Lie groups , the Lie algebras of the SO(2$`n`$+1) group and the SO($`2n`$) group are $`B_n`$ and $`D_n`$, respectively. Their Chevalley bases with the subscript $`j`$, $`1jn1`$, are same:

$$\begin{array}{c}H_j=L_{(2j1)(2j)}L_{(2j+1)(2j+2)},\hfill \\ E_j=\left(L_{(2j)(2j+1)}iL_{(2j1)(2j+1)}iL_{(2j)(2j+2)}L_{(2j1)(2j+2)}\right)/2,\hfill \\ F_j=\left(L_{(2j)(2j+1)}+iL_{(2j1)(2j+1)}+iL_{(2j)(2j+2)}L_{(2j1)(2j+2)}\right)/2.\hfill \end{array}$$
$`(26a)`$
But, the bases with the subscript $`n`$ are different:

$$\begin{array}{c}H_n=2L_{(2n1)(2n)},\hfill \\ E_n=L_{(2n)(2n+1)}iL_{(2n1)(2n+1)},\hfill \\ F_n=L_{(2n)(2n+1)}+iL_{(2n1)(2n+1)},\hfill \end{array}$$
$`(26b)`$
for SO($`2n+1`$), and

$$\begin{array}{c}H_n=L_{(2n3)(2n2)}+L_{(2n1)(2n)},\hfill \\ E_n=\left(L_{(2n2)(2n1)}iL_{(2n3)(2n1)}+iL_{(2n2)(2n)}+L_{(2n3)(2n)}\right)/2,\hfill \\ F_n=\left(L_{(2n2)(2n1)}+iL_{(2n3)(2n1)}iL_{(2n2)(2n)}+L_{(2n3)(2n)}\right)/2,\hfill \end{array}$$
$`(26b)`$
for SO($`2n`$). $`H_k`$ span the Cartan subalgebra, and their eigenvalues are the components of a weight vector $`𝐦=(m_1,\mathrm{},m_n)`$:

$$H_k|𝐦=m_k|𝐦,1kn.$$
$`(27)`$
If the eigenstates for a given weight $`𝐦`$ are degeneracy, this weight is called a multiple weight, otherwise a simple one. $`E_k`$ are called the raising operators and $`F_k`$ the lowering ones. For an irreducible representation denoted by a Young diagram $`[\mu _1,\mu _2,\mathrm{}]`$ of SO$`(D)`$, $`\mu _j\mu _{j+1}`$, there is a highest weight $`𝐌=(M_1,M_2,\mathrm{})`$, which must be simple:

$$\begin{array}{cc}M_j=\mu _j\mu _{j+1},\hfill & 1jn2,\hfill \\ M_{n1}=\mu _{n1}\mu _n,M_n=2\mu _n,\hfill & \mathrm{for}\mathrm{SO}(2n+1),\hfill \\ M_{n1}=\mu _{n1}\mu _n,M_n=\mu _{n1}+\mu _n,\hfill & \mathrm{for}\mathrm{selfdual}\mathrm{representation}\mathrm{in}\mathrm{SO}(2n),\hfill \\ M_{n1}=\mu _{n1}+\mu _n,M_n=\mu _{n1}\mu _n,\hfill & \mathrm{for}\mathrm{antiselfdual}\mathrm{representation}\mathrm{in}\mathrm{SO}(2n).\hfill \end{array}$$
$`(28)`$
We are not interested here in the spinor representations where $`M_n`$ is odd for SO($`2n+1`$) and $`M_{n1}+M_n`$ is odd for SO($`2n`$). For a given irreducible representation $`[\mu _1,\mu _2,\mathrm{}]`$ of SO$`(D)`$, we only need to consider the highest weight state $`|𝐌`$, which satisfies

$$H_k|𝐌=M_k|𝐌,E_k|𝐌=0,1kn,$$
$`(29)`$
because its partners can be calculated by the lowering operators $`F_k`$. In this paper the highest weight state will simply be called the wave functions with the given angular momentum $`[\mu ,\nu ]`$ for simplicity.

Now, we return to our problems. Recalling the Clebsch-Gordan series in Eq. (24), we can rewrite Eq. (23) for the highest weight state $`𝐌`$:

$$\begin{array}{c}𝒴_𝐌^{l,l^{},s,t}(𝐱,𝐲)=\underset{𝐦}{}𝒴_𝐦^l(𝐱)𝒴_{𝐌𝐦}^l^{}(𝐲)l,𝐦,l^{},(𝐌𝐦)|[(l+l^{}s2t),s],𝐌,\hfill \end{array}$$
$`(30)`$
where the subscripts of the harmonic polynomials are changed to the weights for simplicity. $`𝒴_𝐌^{l,l^{},s,t}(𝐱,𝐲)`$ is the highest weight state of the representation $`[(l+l^{}s2t),s]`$. It is a homogeneous polynomial of degrees $`l`$ and $`l^{}`$ in the components of $`𝐱`$ and $`𝐲`$, respectively. Generally speaking, some $`𝒴_𝐌^{l,l^{},s,t}(𝐱,𝐲)`$ may be expressed as a sum where each term is a product of an internal variable $`\xi _j`$ and a homogeneous polynomial $`f(𝐱,𝐲)`$ of lower degree (see p. 042108-5 in ). Since $`𝒴_𝐌^{l,l^{},s,t}(𝐱,𝐲)`$ will be used as a base function for the wave function with a given angular momentum and the combinative coefficient is the function of the internal variables, in this meaning, the base function in the form of $`\xi _jf(𝐱,𝐲)`$ is not independent, and we should find out the independent and complete base functions for any given angular momentum $`[\mu ,\nu ]`$. In the following we are going to prove $`𝒴_𝐌^{q,(\mu +\nu q),\nu ,0}(𝐱,𝐲)`$ and their partners, where $`l=q`$, $`l^{}=\mu +\nu q`$, $`s=\nu `$, $`t=0`$, and $`\nu q\mu `$, constitute a complete set of independent base functions for the total orbital angular momentum $`[\mu ,\nu ]`$. In other words, those total angular momentum states $`𝒴_𝐌^{l,l^{},s,t}(𝐱,𝐲)`$ with $`t>0`$ are not independent, where the sum of the partial angular momentum quantum number $`l`$ and $`l^{}`$ is larger than $`\mu +\nu `$ for the total angular momentum $`[\mu ,\nu ]`$.

The highest weight for the representation $`[\mu ,\nu ]`$ is $`𝐌=(\mu \nu ,\nu ,0,\mathrm{},0)`$. Removing the normalization factor in $`𝒴_𝐌^{q,(\mu +\nu q),\nu ,0}(𝐱,𝐲)`$, which is irrelevant here, we can determine the explicit form for $`𝒴_𝐌^{q,(\mu +\nu q),\nu ,0}(𝐱,𝐲)`$ according to its orders in the components of $`𝐱`$ and $`𝐲`$ and the property of the highest weight state (29), and denote it by the generalized harmonic polynomial $`Q_q^{\mu \nu }(𝐱,𝐲)`$:

$$\begin{array}{cc}\hfill Q_q^{\mu \nu }(𝐱,𝐲)& =\frac{X_{12}^{q\nu }Y_{12}^{\mu q}}{(q\nu )!(\mu q)!}\left(X_{12}Y_{34}Y_{12}X_{34}\right)^\nu \hfill \\ & 𝒴_𝐌^{q,(\mu +\nu q),\nu ,0}(𝐱,𝐲),0\nu q\mu ,\hfill \\ \hfill X_{12}=x_1+ix_2,& X_{34}=x_3+ix_4,Y_{12}=y_1+iy_2,Y_{34}=y_3+iy_4.\hfill \end{array}$$
$`(31)`$
The formula for $`Q_q^{\mu \nu }(𝐱,𝐲)`$ holds for $`D=3`$ ($`x_4=y_4=0`$, $`\nu =0`$ or $`1`$)  and $`D>4`$. When $`D=4`$ we denote the highest weight states by $`Q_q^{(S)\mu \nu }(𝐱,𝐲)`$ and $`Q_q^{(A)\mu \nu }(𝐱,𝐲)`$ for the selfdual representations and the antiselfdual representations, respectively:

$$\begin{array}{c}Q_q^{(S)\mu \nu }(𝐱,𝐲)=\frac{X_{12}^{q\nu }Y_{12}^{\mu q}}{(q\nu )!(\mu q)!}\left(X_{12}Y_{34}Y_{12}X_{34}\right)^\nu \\ Q_q^{(A)\mu \nu }(𝐱,𝐲)=\frac{X_{12}^{q\nu }Y_{12}^{\mu q}}{(q\nu )!(\mu q)!}\left(X_{12}Y_{34}^{}Y_{12}X_{34}^{}\right)^\nu \\ X_{34}^{}=x_3ix_4,Y_{34}^{}=y_3iy_4.\end{array}$$
$`(32)`$
The generalized harmonic polynomial $`Q_q^{\mu \nu }(𝐱,𝐲)`$ is a homogeneous polynomial of degrees $`q`$ and $`(\mu +\nu q)`$ in the components of $`𝐱`$ and $`𝐲`$, respectively. It is a simultaneous eigenfunction of $`_𝐱^2`$, $`_𝐲^2`$, $`_𝐱_𝐲`$, and the total angular momentum operator $`𝐋^2`$ \[see Eq. (9)\],

$$\begin{array}{c}_𝐱^2Q_q^{\mu \nu }(𝐱,𝐲)=_𝐲^2Q_q^{\mu \nu }(𝐱,𝐲)=_𝐱_𝐲Q_q^{\mu \nu }(𝐱,𝐲)=0,\hfill \\ 𝐋^2Q_q^{\mu \nu }(𝐱,𝐲)=C_2([\mu ,\nu ])Q_q^{\mu \nu }(𝐱,𝐲),\hfill \\ C_2([\mu ,\nu ])=\mu (\mu +D2)+\nu (\nu +D4),\hfill \end{array}$$
$`(33)`$
where $`C_2([\mu ,\nu ])`$ is the Casimir calculated by a general formula (see (1.131) in Ref. ). The parity of $`Q_q^{\mu \nu }(𝐱,𝐲)`$ is obviously $`(1)^{\mu +\nu }`$.

It is evident that $`Q_q^{\mu \nu }(𝐱,𝐲)`$ do not contain a function of the internal variables as a factor, neither do their partners due to the rotational symmetry. Therefore, $`Q_q^{\mu \nu }(𝐱,𝐲)`$ are independent base functions for the given angular momentum described by $`[\mu ,\nu ]`$. Now, we are going to prove that $`(\mu \nu +1)`$ base functions $`Q_q^{\mu \nu }(𝐱,𝐲)`$ where $`\nu q\mu `$ are complete for the angular momentum $`[\mu ,\nu ]`$. That is, $`Q_q^{\mu (l\mu )}(𝐱,𝐲)`$ with $`0l\mu q\mu `$ and their partners construct a complete set of linearly independent base functions for the homogeneous polynomials of degree $`l`$ in the components of $`𝐱`$ and $`𝐲`$, excluded those in the forms of $`\xi _jf(𝐱,𝐲)`$, where $`f(𝐱,𝐲)`$ is a homogeneous polynomial of degree $`(l2)`$.

The number of linearly independent homogeneous polynomials of degree $`l`$ in the components of $`𝐱`$ and $`𝐲`$ is

$$M_D(l)=\left(\begin{array}{c}l+2D1\\ 2D1\end{array}\right).$$

After removing those polynomials in the form $`\xi _jf(𝐱,𝐲)`$, the number $`M(l)`$ reduces to $`K(l)`$:

$$\begin{array}{cc}\hfill K_D(l)& =M_D(l)3M_D(l2)+3M_D(l4)M_D(l6)\hfill \\ & =4(l+D3)\left[2l(l+2D6)+(D2)(2D5)\right]\frac{(l+2D7)!}{l!(2D4)!},\hfill \end{array}$$
$`(34)`$
when $`l+2D7`$, which only excludes one case of $`l=0`$ and $`D=3`$, where $`K_3(0)=1`$.

On the other hand, the number of $`Q_q^{\mu (l\mu )}(𝐱,𝐲)`$ with $`0l\mu q\mu `$ and their partners can be calculated directly from Eq. (25):

$$\underset{l/2\mu l}{}(2\mu l+1)d_D([\mu ,(l\mu )])=K_D(l).$$
$`(35)`$
Eqs. (34) and (35) are checked by Mathematica. Thus, we have proved that $`(\mu \nu +1)`$ base functions $`Q_q^{\mu \nu }(𝐱,𝐲)`$ where $`0\nu q\mu `$ are independent and complete for the angular momentum $`[\mu ,\nu ]`$. Any function with the angular momentum $`[\mu ,\nu ]`$ in the system can be expanded with respect to the base functions $`Q_q^{\mu \nu }(𝐱,𝐲)`$, where the coefficients are functions of internal variables.

From Eq. (30), for a given total orbital angular momentum $`[\mu ,\nu ]`$ there are infinite number of wave functions $`𝒴_𝐌^{(q+t),(\mu +\nu +tq),\nu ,t}(𝐱,𝐲)`$ combined from different partial angular momentum states. Now, what we have proved is that only a finite number of partial angular momentum states ($`t=0`$) are involved in the complete set of independent base functions $`Q_q^{\mu \nu }(𝐱,𝐲)`$ for a given total orbital angular momentum $`[\mu ,\nu ]`$.

## V Generalized radial equations

In the preceding section we proved that any function with angular momentum $`[\mu ,\nu ]`$ in the quantum three-body system of $`D`$ dimensions can be expanded with respect to the base functions $`Q_q^{\mu \nu }(𝐱,𝐲)`$

$$\mathrm{\Psi }_𝐌^{[\mu ,\nu ]}(𝐱,𝐲)=\underset{q=\nu }{\overset{\mu }{}}\psi _q^{\mu \nu }(\xi _1,\xi _2,\xi _3)Q_q^{\mu \nu }(𝐱,𝐲),$$
$`(36)`$
where the coefficients $`\psi _q^{\mu \nu }(\xi _1,\xi _2,\xi _3)`$ are called the generalized radial functions. When substituting Eq. (36) into the Schrödinger equation (22), the main calculation in the derivation is to apply the Laplace operator (21) to the function $`\mathrm{\Psi }_𝐌^{[\mu ,\nu ]}(𝐱,𝐲)`$. The calculation consists of three parts. The first is to apply the Laplace operator to the generalized radial functions $`\psi _q^{\mu \nu }(\xi _1,\xi _2,\xi _3)`$, which can be calculated by replacement of variables:

$$\begin{array}{c}^2\psi _q^{\mu \nu }(\xi _1,\xi _2,\xi _3)=\{4\xi _1_{\xi _1}^2+4\xi _2_{\xi _2}^2+2D(_{\xi _1}+_{\xi _2})+(\xi _1+\xi _2)_{\xi _3}^2\hfill \\ +4\xi _3(_{\xi _1}+_{\xi _2})_{\xi _3}\}\psi ^{\mu \nu }_q(\xi _1,\xi _2,\xi _3),\hfill \end{array}$$
$`(37)`$
where $`_\xi `$ denotes $`/\xi `$ and so on. The second is to apply it to the generalized harmonic polynomials $`Q_q^{\mu \nu }(𝐱,𝐲)`$. This part is vanishing because $`Q_q^{\mu \nu }(𝐱,𝐲)`$ satisfies the Laplace equation. The last is the mixed application

$$\begin{array}{c}2\{\left(_{\xi _1}\psi _q^{\mu \nu }\right)2𝐱+\left(_{\xi _3}\psi _q^{\mu \nu }\right)𝐲\}_𝐱Q_q^{\mu \nu }+2\{\left(_{\xi _2}\psi _q^{\mu \nu }\right)2𝐲+\left(_{\xi _3}\psi _q^{\mu \nu }\right)𝐱\}_𝐲Q_q^{\mu \nu }.\hfill \end{array}$$
$`(38)`$
From the definition (31) for $`Q_q^{\mu \nu }(𝐱,𝐲)`$ we have

$$\begin{array}{cc}𝐱_𝐱Q_q^{\mu \nu }=qQ_q^{\mu \nu },\hfill & 𝐲_𝐲Q_q^{\mu \nu }=(\mu +\nu q)Q_q^{\mu \nu }\hfill \\ 𝐲_𝐱Q_q^{\mu \nu }=(\mu q+1)Q_{q1}^{\mu \nu },\hfill & 𝐱_𝐲Q_q^{\mu \nu }=(q\nu +1)Q_{q+1}^{\mu \nu }.\hfill \end{array}$$
$`(39)`$
Hence, we obtain the generalized radial equation, satisfied by the $`(\mu \nu +1)`$ functions $`\psi _q^{\mu \nu }(\xi _1,\xi _2,\xi _3)`$:

$$\begin{array}{c}^2\psi _q^{\mu \nu }+4q_{\xi _1}\psi _q^{\mu \nu }+4(\mu +\nu q)_{\xi _2}\psi _q^{\mu \nu }+2(\mu q)_{\xi _3}\psi _{q+1}^{\mu \nu }+2(q\nu )_{\xi _3}\psi _{q1}^{\mu \nu }\\ =2\left(EV\right)\psi _q^{\mu \nu },\end{array}$$
$`(40)`$
where $`^2\psi _q^{\mu \nu }`$ is given in Eq. (37). Only three invariant variables $`\xi _1`$, $`\xi _2`$ and $`\xi _3`$ are involved both in the equations and in the functions. When $`D=4`$, Eq. (40) holds for the generalized radial functions either in $`[(S)\mu ,\nu ]`$ or in $`[(A)\mu ,\nu ]`$, because two representations incorporate to one irreducible representation of the O(4) group when the space inversion is considered. When $`D=3`$, the equations for the functions in $`[\mu ,0]`$ and in $`[\mu ,1]`$ are different although two representations $`[\mu ,0]`$ and $`[\mu ,1]`$ are equivalent, because the functions have different parity.

At last, we discuss rotational variables and the volume element of the configuration space. We fix the body-fixed frame such that $`𝐱`$ is parallel with its $`D`$th axis, and $`𝐲`$ is located in its $`(D1)D`$ hyperplane with a non-negative $`(D1)`$th component. That is, in the body-fixed frame, the nonvanishing components of two Jacobi coordinate vectors $`𝐱^{}`$ and $`𝐲^{}`$ are

$$x_D^{}=\xi _1^{1/2},y_{D1}^{}=\left(\xi _2\xi _3^2/\xi _1\right)^{1/2},y_D^{}=\xi _3\xi _1^{1/2}.$$
$`(41)`$
Let $`R=R^{(1)}R^{(2)}`$SO($`D`$) rotate the center-of-mass frame to the body-fixed frame:

$$\begin{array}{c}R^{(1)}=R_{12}(\theta _1)R_{31}(\theta _2)R_{43}(\theta _3)R_{54}(\theta _4)\mathrm{}R_{D(D1)}(\theta _{D1}),\hfill \\ R^{(2)}=R_{12}(\phi _1)R_{31}(\phi _2)R_{43}(\phi _3)R_{54}(\phi _4)\mathrm{}R_{(D1)(D2)}(\phi _{D2}),\hfill \\ R𝐱^{}=𝐱,R𝐲^{}=𝐲,\hfill \end{array}$$
$`(42)`$
where, for example, $`R_{12}(\theta )`$ is a rotation on the hyperplane with the first and the second axes through $`\theta `$ angle:

$$R_{12}(\theta )=\left(\begin{array}{ccc}\mathrm{cos}\theta & \mathrm{sin}\theta & 0\\ \mathrm{sin}\theta & \mathrm{cos}\theta & 0\\ 0& 0& \mathrm{𝟏}_{D2}\end{array}\right).$$

$`(D1)`$ $`\theta _j`$ and $`(D2)`$ $`\phi _k`$ are the rotational variables, called the generalized Euler angles. Through a straightforward calculation, we obtain

$$\begin{array}{c}x_1+ix_2=\xi _1^{1/2}e^{i\theta _1}\underset{a=2}{\overset{D1}{}}\mathrm{sin}\theta _a,x_3+ix_4=\xi _1^{1/2}\left(\mathrm{cos}\theta _2\mathrm{sin}\theta _3+i\mathrm{cos}\theta _3\right)\underset{a=4}{\overset{D1}{}}\mathrm{sin}\theta _a,\hfill \\ y_1+iy_2=\xi _3\xi _1^{1/2}e^{i\theta _1}\underset{a=2}{\overset{D1}{}}\mathrm{sin}\theta _a+\left(\xi _2\xi _3^2/\xi _1\right)^{1/2}e^{i\theta _1}\hfill \\ \times \left\{i\underset{a=1}{\overset{D2}{}}\mathrm{sin}\phi _a+\underset{a=1}{\overset{D2}{}}\mathrm{cos}\theta _{a+1}\mathrm{cos}\phi _a\left(\underset{b=2}{\overset{a}{}}\mathrm{sin}\theta _b\right)\left(\underset{c=a+1}{\overset{D2}{}}\mathrm{sin}\phi _c\right)\right\},\hfill \end{array}$$

$$\begin{array}{c}y_3+iy_4=\xi _3\xi _1^{1/2}\left(\mathrm{cos}\theta _2\mathrm{sin}\theta _3+i\mathrm{cos}\theta _3\right)\underset{a=4}{\overset{D1}{}}\mathrm{sin}\theta _a+\left(\xi _2\xi _3^2/\xi _1\right)^{1/2}\hfill \\ \times \{\mathrm{cos}\phi _1\mathrm{sin}\theta _2\underset{a=2}{\overset{D2}{}}\mathrm{sin}\phi _a+(\mathrm{cos}\theta _2\mathrm{cos}\theta _3i\mathrm{sin}\theta _3)\mathrm{cos}\phi _2\underset{a=3}{\overset{D2}{}}\mathrm{sin}\phi _a\hfill \\ +(\mathrm{cos}\theta _2\mathrm{sin}\theta _3+i\mathrm{cos}\theta _3)\underset{a=3}{\overset{D2}{}}\mathrm{cos}\theta _{a+1}\mathrm{cos}\phi _a\left(\underset{b=4}{\overset{a}{}}\mathrm{sin}\theta _b\right)\left(\underset{c=a+1}{\overset{D2}{}}\mathrm{sin}\phi _c\right)\},\hfill \end{array}$$
$`(43)`$
where $`_{a=b+1}^bF_a=1`$. The volume element of the configuration space is

$$\underset{j=1}{\overset{D}{}}dx_jdy_j=\frac{1}{4}\left(\xi _1\xi _2\xi _3^2\right)^{(D3)/2}d\xi _1d\xi _2d\xi _3\underset{j=1}{\overset{D1}{}}\left(\mathrm{sin}\theta _j\right)^{j1}d\theta _j\underset{k=1}{\overset{D2}{}}\left(\mathrm{sin}\phi _k\right)^{k1}d\phi _k.$$
$`(44)`$
## VI Conclusions

After separating the motion of center of mass, we have defined the homogeneous polynomial $`Q_q^{\mu \nu }(𝐱,𝐲)`$ of degree $`q`$ and $`(\mu +\nu q)`$ in the components of the Jacobi coordinate vectors $`𝐱`$ and $`𝐲`$, respectively. $`Q_q^{\mu \nu }(𝐱,𝐲)`$ is a solution of the Laplace equation. We have proved that $`(\mu \nu +1)`$ generalized harmonic polynomials $`Q_q^{\mu \nu }(𝐱,𝐲)`$ constitute a complete set of independent base functions for the total orbital angular momentum $`[\mu ,\nu ]`$. Any wave function with the given angular momentum in the system can be expanded with respect to the base functions, where the coefficients are the functions of the internal variables, called the generalized radial functions. The three-body Schrödinger equation with a spherically symmetric potential $`V`$ in $`D`$ dimensions reduces to the generalized radial equations satisfied by the generalized radial functions. Only three internal variables are involved in the functions and equations. The number of both the functions and the equations for the given angular momentum $`[\mu ,\nu ]`$ is finite and equal to $`\left(\mu \nu +1\right)`$. Only a finite number of partial angular momentum states are involved in constructing the generalized harmonic polynomials $`Q_q^{\mu \nu }(𝐱,𝐲)`$, and the contributions from the remaining partial angular momentum states have been incorporated into those from the generalized radial functions.

The generalization of this method to a quantum $`N`$-body system in $`D`$-dimensions is straightforward. The difficulty is how to obtain the unified forms for the generalized harmonic polynomials, because it needs $`D1`$ vectors to determine the body-fixed frame and there are $`N1`$ Jacobi coordinate vectors. The cases with $`N<D`$ are very different to the cases with $`ND`$. We will study this problem elsewhere.

ACKNOWLEDGMENTSThis work is supported by the National Natural Science Foundation of China.


format: text

text:
# 1 . Introduction

## 1 . Introduction

As estimated by Mao and Paczyński (1991) several percent of all microlensing events in our Galaxy should be caused by binaries acting as lenses. The first such case, the event OGLE-7 (Udalski et al. 1994a) was analyzed also by Dominik (1999), who showed the existence of several different binary lens models of the event.

The amount of information one may obtain from a binary lens light curve depends strongly on the frequency of observations and the kind of the event itself. The events with well sampled caustic crossing may be used to study the effect of limb darkening in the atmospheres of lensed stars (e.g., Albrow et al. 1999b) or to set limits on the source velocity relative to the lens (e.g., Albrow et al. 1999a). In some cases with caustic crossings and a pronounced effect of Earth motion parallax, the mass of the lens can be found (An et al. 2001). Such finding is not possible for single gravitational lenses and rather uncommon for doubles. The typical caustic crossing takes several hours, and the observations should be taken every ten or twenty minutes to be useful for detailed modeling. The standard frequency of observations used by the gravitational lens survey teams (including OGLE-II – Udalski et al. 2000) is one or two per night – insufficient to treat the caustic crossings as separate issue. Some properties of binary lenses can be found even without well sampled caustic crossings.

The only systematic study to date of the binary lens population was made for the MACHO data (Alcock et al. 2000). It gives the light curves and binary lens models for 21 events observed in 1993–1998. Many of the events were known to be caused by binary lenses early enough to make the follow up observations with the dense sampling of the caustic crossing. In four cases the relative proper motion of the lens was measured with the accuracy sufficient to study the likelihood of the lens masses. The distribution of the binary mass ratio is also given.

Woźniak et al. (2001) present a list of 520 gravitational lens candidates among OGLE-II stars based on the new data processing using DIA photometry (Alard and Lupton 1998, Alard 2000). The database includes observations of seasons 1997–1999. Most of the candidates have been found algorithmically. The list is supplemented by candidates selected by visual inspection of the light curves and also by some events which were known as candidates in earlier study (Udalski et al. 2000). The sample of the gravitational lens candidates belongs to a larger sample of transient sources. The transient sample (Woźniak 2000, Woźniak et al. 2001) has been selected in fully algorithmic manner from all OGLE-II light curves and serves as the basis of the present study.

Our goal is to find the candidates for binary lens events, find their models and make some survey of the binary lens population properties. In the next Section we describe our selection method for finding both single and binary lens candidates in the transient sample of Woźniak et al. (2001). In Section 3 we describe the procedure of fitting binary lens models to the data. The results are described in Section 4, while the graphical material is shown in Appendix. Discussion follows in Section 5.

## 2 . Choice of Candidates

Woźniak (2000) and Woźniak et al. (2001) describe the construction of the transient sample based on OGLE-II data for seasons 1997–1999 processed using the DIA photometry. The sample contains 4424 sources. The transients are selected algorithmically, as sources which have (almost) constant flux of radiation and undergo at most two brightening episodes during observation span. We follow Woźniak in defining a threshold value for physical flux variation. We assume that the variability of the source is limited to a sequence containing not more than 50% of all observations. The remaining data points serve to define a constant flux base. The choice of the “variable” part is made in a way minimizing the scatter of the remaining observations (the “base”). We calculate $`\sigma _{\mathrm{base}}`$ – a rms deviation of flux measurements from its mean value $`F_{\mathrm{base}}`$. Following Woźniak (2000) we adjust the calculated parameter taking into account possibly non-Gaussian distribution of flux values around the mean. Finally, we obtain the threshold parameter $`\delta =\mathrm{max}\{\sigma _{\mathrm{base}},\sigma _{\mathrm{ph}}\}`$, where $`\sigma _{\mathrm{ph}}`$ is the mean photometric error estimated for the points belonging to the constant flux database. Woźniak et al. (2001) define a brightening episode as 3 consecutive points at least $`5\delta `$ above the base flux or 4 points at least $`4\delta `$ above.

The light curves of binary lens events usually have two maxima and it is necessary to have at least five measurements to distinguish them. For this reason and having in mind some safety margin, we require the transients to have at least 7 points which are $`5\delta `$ above the base flux. Some experience with Woźniak et al. (2001) database suggests an extra criterion for selecting the interesting events, which helps removing DIA photometry artifacts from the investigated light curves. In further investigation we include only the light curves with at least 7 points which are $`5\delta `$ above the base flux with at least one point among them at $`15\delta `$ level.

Using our criteria we automatically search the transient database getting 390 light curves. We visually inspect all the selected cases. For each light curve we try a single lens fit taking into account blending. In the majority of cases comparison of the best single lens fit with the data shows a good agreement, in the sense of visual impression and the $`\chi ^2`$ value. There are cases of weak transients, which can be formally fitted with single lens model based on $`\chi ^2`$ value, which are nevertheless rejected after visual inspection, usually the changes in flux look as random with a period of increased mean flux value. There are also cases which under visual inspection show a very good agreement with a single lens model light curve, but have high $`\chi ^2`$ value. Usually these are the strongest lensing events. We treat them as single lens events not investigating further the reasons for low formal quality of the fits. The second most common category among selected transients are artifacts of DIA photometry. They show irregular behavior dominated by noise. In many cases series of sources close to each other show very similar light curves of this kind. In a few cases the light curves suggest periodic variability. We exclude the possibility that there are binary lenses belonging to this category. The binary lens candidates are chosen mostly on the basis of visual inspection. Some show the characteristic “U-shape” of the caustic crossing. Some have multiple flux maxima. The other show asymmetry. All visually selected candidates have high values of $`\chi ^2`$ for a single lens model, usually exceeding 3 per degree of freedom. There are also several cases of good quality light curves, which do not resemble a single or binary lens model variability. We put them into “other objects” category.

Among 390 transient selected we find 215 single lenses and 141 unclassifiable cases of poor photometry. We treat 27 of the remaining 34 light curves as possible binary lens candidates, and 7 – as “other objects”.

## 3 . Fitting Binary Lens Models

The models of the two point mass lens were investigated by many authors (Schneider and Weiss 1986, Mao and DiStefano 1995, DiStefano and Mao 1996, Dominik 1998, to mention few). The methods applicable for extended sources have recently been described by Mao and Loeb (2001). While we use mostly the point source approximation, we extensively employ their efficient numerical schemes for calculating the binary lens caustic structure and source magnification.

We fit binary lens models using the $`\chi ^2`$ minimization method for the light curves. It is convenient to model the flux at the time $`t_i`$ as:

$$F(t_i)=F_b+F_s\times A(t_i)F_{\mathrm{base}}+F_s\times \left(A(t_i)1\right)$$
$`(1)`$
where $`F_b`$ is the blended flux (from the source close neighbors and possibly the lens), $`F_s`$ – the flux of the source which is lensed, and $`A(t_i)`$ is the lens magnification. The combination $`F_{\mathrm{base}}F_b+F_s`$ is the total flux long before or long after the event. The DIA photometry only gives the flux differences, so one can use $`\stackrel{~}{F}_iF_iF_{\mathrm{base}}`$ as “measured” quantities. Using this notation one has for the $`\chi ^2`$ (Woźniak 2000):

$$\chi ^2=\underset{i=1}{\overset{N}{}}\frac{\left((A_i1)F_s\stackrel{~}{F}_i\right)^2}{\sigma _i^2}$$
$`(2)`$
where $`\sigma _i`$ are the errors of the flux measurement taken from the DIA photometry. For self-consistency we check the value of $`\chi ^2`$ calculated for the points defining the constant flux base.

Binary lens models are possibly and typically non unique (Dominik 1999). The presence of caustics and cusps in the lens theory (Schneider, Ehlers and Falco 1992, Blandford and Narayan 1992) makes the $`\chi ^2`$ dependence on the model parameters complicated and discontinuous for point sources. This suggests that the search for the $`\chi ^2`$ global minimum may be difficult because of the likely existence of many local ones. For this reason we combine the scan of the parameter space with the standard minimization techniques in our search.

The light curves we investigate do not have well sampled caustic crossings, so we use the point source model in the majority of calculations. For the same reason we are also not able to use the strategy of Albrow et al. (1999c) for finding binary lens models. We also neglect the parallax and/or binary rotation effects, which further simplifies the models.

The binary system consists of two masses $`m_1`$ and $`m_2`$, where by convention $`m_1m_2`$. The Einstein radius of the binary lens is defined as:

$$r_\mathrm{E}=\sqrt{\frac{4G(m_1+m_2)}{c^2}\frac{d_{\mathrm{OL}}d_{\mathrm{LS}}}{d_{\mathrm{OS}}}}$$
$`(3)`$
where $`G`$ is the constant of gravity, $`c`$ is the speed of light, $`d_{\mathrm{OL}}`$ is the observer–lens distance, $`d_{\mathrm{OS}}`$ is the observer–source distance, and $`d_{\mathrm{LS}}`$ is the distance between the lens and the source. The Einstein radius serves as a length unit and the Einstein time: $`t_\mathrm{E}=r_\mathrm{E}/v_{}`$, where $`v_{}`$ is the lens velocity relative to the line joining the observer with the source, serves as a time unit. The passage of the source in the lens background is defined by seven parameters: $`qm_1/m_2`$ ($`0<q1`$) – the binary mass ratio, $`d`$ – binary separation expressed in $`r_\mathrm{E}`$ units, $`\beta `$ – the angle between the source trajectory as projected onto the sky and the projection of the binary axis, $`b`$ – the impact parameter, $`t_0`$ – the time of closest approach of the source to the binary center of mass, $`t_E`$ – the Einstein time, and $`F_s`$ the (not magnified) source flux. For all other parameters fixed, $`F_s`$ can be obtained from the (linear) equation $`\chi ^2/F_s=0`$ with the obvious limitation $`0F_sF_{\mathrm{base}}`$. Thus we are left with the six dimensional parameter space to look for minima. We begin with a scan of the parameter space using a logarithmic grid of points in $`(q,d)`$ plane ($`10^3q1`$, $`0.1d10`$) and allowing for continuous variation of the other four parameters (inclusion of $`q`$ and $`d`$, which change the caustic structure of the binary, in the automatic routine finding $`\chi ^2`$ minimum may cause problems with its convergence – see An et al. (2001) Jaroszyński and Mao (2001)). The other dimensions are searched using AMOEBA routine (Press et al. 1992). We use many starting points in $`(\beta ,b,t_0,t_E)`$ subspace. The expected value for the Einstein time $`t_E`$ is based on the event duration. The expected time of the closest approach to the binary center of mass $`t_0`$ should be close to the time of the maximal observed flux. The direction of the source motion relative to the binary axis $`\beta `$ is chosen at random. For each choice of $`\beta `$ there is a defined range of impact parameter $`b`$ giving rise to a “binary character” of the light curve. Practically it means that the source should either cross a caustic or pass not too far from it.

Using the point source approximation we find the best binary lens models for every $`(q,d)`$ on the grid. With the fitted value of the Einstein time $`t_\mathrm{E}`$ and the typical velocity of a source relative to the observer–lens line $`v`$, we are able to estimate the linear size of the Einstein radius as projected into the source plane $`\stackrel{~}{r}_\mathrm{E}vt_\mathrm{E}`$. The fitted deblended source flux allows for a crude estimate of the star radius $`R_s`$, and corresponding dimensionless source radius $`r_sR_s/\stackrel{~}{r}_\mathrm{E}`$. For self consistency we repeat the fitting procedure at every grid point using the extended source model of the estimated size and starting from the solution obtained in the point source approximation. Finally we search the regions of $`(q,d)`$ plane with the lowest values of $`\chi ^2`$ on a finer grid.

## 4 . Results

We have applied the fitting procedure to the 27 selected events. We also visually compare the top quality models with the data, checking whether the fitted light curves are “natural”, not having superficial peaks in places not covered by observations. In eighteen cases the binary lens fits are formally good and natural. They are also substantially better than single lens fits applied to the same data (the binary fits give the $`\chi ^2`$ values lower at least by 50, and typically by several hundreds or more as compared with single lens fits.) There are three candidates with formally good fits which are rejected after the visual inspection of the models. The other six are rejected because of the high $`\chi ^2`$ value. The plots showing light curves and our best fits for the rejected candidates are included in Appendix.

Nine of the events are on the Udalski et al. (2000) list of microlensing events, and five of them were detected in real time and included in the alerts of the Early Warning System (Udalski et al. 1994b). One of the events (SC21\_6195) has been included in Alcock et al. (2000) study of binary lenses as candidate event 97-BLG-d2. Not all events marked as possible binary lens candidates by Udalski et al. (2000) are on our list. We give the comparison of the two lists in the Section 5.

### 4.1 The Binary Lens Models

Of the eighteen binary lens candidates, twelve show caustic crossing, five at least one close cusp approach with sharp light maximum, and one has two shallow maxima explained by passages at greater distances from cusps.

In Table 1 we give the names, positions and observed magnitudes in $`I`$ band of the objects long before or long after the event.

In Table 2 we give the parameters for the best models of each of the events. In some cases we include also another solution if it belongs to a distinct $`\chi ^2`$ minimum and lies inside the confidence region. Our best fits are also shown on the plots in Appendix. For each binary there are three separate plots. The first shows the source trajectory as projected onto the lens plane with caustic structure and binary components included. The model light curves and observed fluxes are shown in the second plot. The third diagram shows the $`\chi ^2`$ confidence regions in the $`\mathrm{lg}q`$$`\mathrm{lg}d`$ plane.

The source radius $`r_s`$ in Einstein units, which we give in the last column is not a fitted parameter. It is only estimated and used for the self consistency check.

We concentrate mostly on the physical parameters of the binary ($`q`$ and $`d`$). Examining the plots one can see that not all events can be well constrained and for some the $`1\sigma `$ confidence regions are rather large and irregular.

### 4.2 Distribution of Mass Ratios and Separations

There are 12 caustic crossing binary lens events among the total of 18. The conditional probability of a binary lens to cause a caustic crossing event is well defined (Mao and Paczyński 1991). It is proportional to the angle averaged width of caustic measured in the direction perpendicular to the source path and can be calculated for any binary lens as $`w(\mathrm{lg}q,\mathrm{lg}d)`$, where we use the dependence on logarithms as more convenient. The most efficient lenses have $`q1`$ and $`d1`$.

The rate of observed caustic crossing events caused by lenses with given parameters $`P_{\mathrm{BCC}}(\mathrm{lg}q,\mathrm{lg}d)`$ also depends on the distribution of binary stars parameters $`P_{\mathrm{bin}}(\mathrm{lg}q,\mathrm{lg}d)`$. Following Mao and Paczyński (1991) we assume, that the uniform distribution of binary systems in the logarithm of orbital period, which holds over six decades (Abt 1983) translates into uniform distribution in the logarithm of semimajor axis, and consequently in $`\mathrm{lg}d`$. Under this assumption the distribution of parameters among observed events is given by:

$$P_{\mathrm{BCC}}(\mathrm{lg}q,\mathrm{lg}d)P_{\mathrm{bin}}(\mathrm{lg}q)w(\mathrm{lg}q,\mathrm{lg}d)$$
$`(4)`$
and, in principle, may be used to constrain $`P_{\mathrm{bin}}(\mathrm{lg}q)`$. The binary systems database has been investigated by Trimble (1990), who shows that the mass ratio distribution is uniform in $`\mathrm{lg}q`$. This dependence is not valid for small $`q`$ where the binaries are less abundant. We consider distributions of the form:

$$P_{\mathrm{bin}}(\mathrm{lg}q)\{\begin{array}{cc}q/q_0\hfill & qq_0\hfill \\ 1\hfill & qq_0\hfill \end{array}$$
$`(5)`$
using $`0.03q_00.3`$ which defines the location of the power law break. Using the above formulae one can find the integral probability distribution for the mass ratio and separation among the binaries causing BCC events. It is not obvious how to model the probability distribution for the cases without caustic crossings. We assume that our calculation, which includes only the caustic crossing events, may be applied to a sample including both caustic crossing and cusp approach events. The comparison of the theoretical predictions with the distribution of the binary lens sample, excluding the event 1998-BUL-28, which shows no sharp feature in the light curve. is shown in Fig. 1.

Our binary sample consists of the best models for each event. There are two cases with extreme mass ratios: SC20\_1793 and SC20\_3525. Both are cusp approach events. The first of the cases has also many other fits within the $`1\sigma `$ confidence limits with $`q>0.1`$. The other case is a more serious candidate for an extreme mass ratio binary lens, but models with $`q>0.1`$ are present within $`2\sigma `$ confidence limits.

The mass ratio and separation distributions of our sample (excluding possible extreme cases) is in rough agreement with theoretical expectations. The number of cases is low, so it is too early to make definite conclusions regarding the distribution of parameters among the binary stars.

### 4.3 Distribution of Einstein Radius Crossing Time

The Einstein radius crossing time is a physical parameter depending on the distances, lens mass and relative velocity, which can be found both for single and binary lenses. We show the histogram of this quantity for our sample of binary lens candidates in Fig. 2.

The statistics is poor, but the over-representation of the long lasting events, as compared with our single lens sample is apparent. The binary lenses should be twice as massive on average, which translates to the $`\sqrt{2}`$ factor in Einstein time. This is, however, not sufficient to explain the difference between the distributions.

### 4.4 The Blending Parameter $`f`$

This parameter, showing the fraction of the base flux belonging to the lensed star (while the rest is blended) can be fitted with greater confidence than for single lenses. We force this parameter to belong to the region $`0.01f0.99`$. Only one of our models lies on its upper boundary. (Sometimes one can obtain unphysical solutions with $`f>1`$ and low $`\chi ^2`$ value which represent blends of negative fluxes. Other rather unlikely models can have $`f0`$ and very high amplification – usually related to crossing of a small caustic. None of our models falls into the second category.) We show a histogram of $`f`$ values in Fig. 3.

For better comparison of single and binary lens models we show their positions on $`\mathrm{lg}(t_\mathrm{E})`$$`f`$ diagram. We use the convention employed also by Woźniak et al. (2001), which preserves single lens fits with unphysical values of blending parameter ($`f>1`$) implying the negative flux of light from the blend. (If forced to $`f1`$, these models cluster on the boundary of permitted region.) Another unphysical region on the diagram consists of points with $`f0`$ and high $`t_\mathrm{E}`$. The presence of single lens models in the unphysical regions of the diagram is explained by Woźniak and Paczyński (1997) and is related to the degeneracy among fitted parameters in certain situations.

In Fig. 4 we show the single and binary lens models belonging to our sample. The median points for the two kinds of lenses are different: the binaries have systematically longer duration (73<sup>d</sup>versus 30<sup>d</sup>for singles) and are more blended ($`f_{\mathrm{med}}=0.36`$ versus 0.62 for singles). It is interesting to point out that for all of the 520 single lenses of Woźniak et al. (2001) the median duration is even shorter (23<sup>d</sup>) and the blending less important ($`f_{\mathrm{med}}=0.88`$). The simple rejection of all models with $`f>1`$ gives almost equal values of median blending parameters ($`f_{\mathrm{med}}0.37`$) for binaries, singles from our sample and singles from Woźniak et al. (2001). The median Einstein times become longer for $`f<1`$ subsamples of single lens models, but still shorter than for binaries (42<sup>d</sup>/33<sup>d</sup>versus 73<sup>d</sup>). There is no simple criterion for removing models with $`f0`$ and we do not pursue this issue.

### 4.5 Contribution of Binary Lenses to Microlensing Optical Depth

Gravitational microlensing is detected because of the temporary presence of an extra flux of light caused by source amplification. The detailed modeling of the detection process is a separate issue, but it is clear that both the duration of the event and the strength of the lens play here important roles. The customary definition of lensing event (Paczyński 1996) requires the object to be amplified at least 1.34 times, which corresponds to the source position inside the Einstein ring. Thus the contribution of a single event to the optical depth is proportional to the time it spends above the threshold amplification. For a blended source the observed amplification is not so simply related to the lens position. In DIA photometry it is rather the flux difference, which is important in discovering lensing events. We estimate the relative contribution of binary and single events to the microlensing optical depth comparing the extra light related to the two kinds of phenomena. The extra light corresponding to a particular event is given as:

$$a_i=F_\mathrm{s}(A_i(t)1)dt$$
$`(6)`$
where index $`i`$ enumerates the events. Comparing the whole extra light of our binary lenses with the extra light of single lenses we obtain:

$$\frac{_{\mathrm{bin}}a_i}{_{\mathrm{sin}}a_i}=0.03$$
$`(7)`$
which suggests a 3% contribution of binary lenses to optical depth. The base fluxes of objects belonging to our sample which undergo single lensing are systematically higher than for objects undergoing binary lensing. This may be the result of our selection criteria. The effect diminishes the importance of binary lenses. It is possible that sources amplified by binary lenses are on average fainter as compared to sources undergoing lensing by single masses for a wide class of selection criteria. We check also other possibilities using our sample of single and binary lenses but changing the population of sources. This should be treated as a thought experiment which serves only to estimate the likely systematic errors of the above result – it is not guaranteed that for a different population of sources the same population of lenses would be discovered.

In the first experiment we assume that the base flux ($`F_{\mathrm{base}}`$) of all objects is the same. We keep unchanged the other parameters of the lens models (including $`f`$). Calculating the ratio of extra light contributed by binary and single lenses we obtain the value 0.11. In the second experiment we assume that the source fluxes ($`F_\mathrm{s}`$) are the same for all objects. This gives 0.05 for the calculated ratio.

The preliminary analysis of our sample shows that the recognizable binary lenses make about 8% of all lensing events in excellent agreement with Mao and Paczyński (1991). Estimating the contribution to optical depth by binaries is a more difficult issue. Using three different approaches we obtain the values between 3% and 11%.

## 5 Discussion

Our fitting procedure is based on two important simplifications. First we use mostly a point source model while calculating the source magnification. We make some self consistency checks using also extended source models, but we are not able to fit its size. The effects of source size are important during caustic crossing and/or cusp approaches. Some of the models would have probably been changed if the data were better sampled. Only the event SC36\_4030 (1999-BUL-11), for which there are many observations of the second caustic crossing might be used to constrain the source size. The observations near the caustic are, however, marked as not to be trusted by DIA photometry pipeline and we do not include them. The second simplification is the fact that we neglect all effects of the Earth motion, rotation of the binary, and possible acceleration of the source. In fact we have attempted fits including the effects of parallax/acceleration, but finally we have decided not to use them. Formally for long lasting events the Earth motion should be included. In the case of single lens events, the Einstein time is simply related to the duration of strong amplification. In the case of a binary this is not necessarily so: the passage through or near a small caustic, which implies strong lensing, may last much shorter than the Einstein time. The observations of only weekly magnified source are not discriminating between the models. Thus, even for formally “long” events it is not easy to measure the effects of changes in the relative velocity of the source. Introducing extra parameters describing the relative acceleration of the source may serve as to lower the $`\chi ^2`$ value appreciably (at least for some events), but the extra parameters cannot be well constrained.

The estimates of the caustic size which we use for self consistency checks may also be used to place limits on possible lens mass and luminosity. In some cases the resulting lens brightness is too high. We do not consider this a serious problem, since the relative velocity of the lens may be substantially different from the “typical” value used in the estimates.

For two events (SC20\_1793 and SC20\_3525) the formally best binary lens models have extreme mass ratios ($`q<0.02`$), which may suggest they contain brown dwarfs or Jupiter-like objects. The same events can be, however, modeled using binary lens models with much higher mass ratio, which cannot be rejected with high confidence. Some parts of the parameter confidence region also reach low mass ratio values ($`q10^2`$) for SC35\_2526, but this is a case of very weakly constrained solutions. Statistically acceptable models may belong to a wide region of mass ratios. Thus we cannot conclude a presence of a planetary system among binary lens candidates.

The distribution of Einstein radius crossing times among the models requires further investigation. The occurrence of long lasting events is certainly more likely in our sample as compared to single lens sample. The same effect but to lesser extent is also present in Alcock et al. (2000) sample of binary lens candidates. Since the recognition of binary events requires longer baselines, there is certainly some selection effect involved. The essential part of the light curves – caustic crossings or cusp approach (approaches) may last only a small fraction of the Einstein time. That puts some lower limit on the inter-caustic time, depending on the frequency of observations. With observations made once a day we are biased toward long lasting events.

### 5.1 Comparison with Previous Study

We are using DIA processed observations of OGLE-II in the direction of Galactic bulge made in seasons 1997–1999 (Woźniak et al. 2001). The same observational material, but processed differently, has been analyzed by Udalski et al. (2000), who explicitly mark the binary lens candidates in their Table 2. We list all their candidates in our Table 3, and show their status according to the present study.

Four of the Udalski et al. (2000) binary lens candidates (see Table 3) are not in the transient database of Woźniak et al. (2001), and another one has not passed our selection criteria. Two of the events were also on the list of our 27 candidates, but have been rejected after the attempt to fit them with binary lens model. In five cases we believe the events are due to single lens rather than binary. (Three cases are strong: the DIA photometry does not show the irregularities of the light curves near maximum, which are present on the plots in the Udalski et al. (2000) catalog. Two other cases are examples of noisy light curves and we treat them as poor quality single lens cases.) In nine cases we agree that the events are due to binary lenses, but in two of them we think there is no caustic crossing.

The comparison of our sample with the list of 520 microlensing events of Woźniak et al. (2001) is also due, because in both cases the same database has been used. Our approach differs since we apply stronger selection criteria to have a better chance of distinguish the binary lens candidates. The possible binary lens events are not marked in the other study, so a direct comparison is not possible. Of our original 27 binary lens candidates only 3 are missing on the Woźniak et al. (2001) list. After attempting to model the 27 candidate events as binary lenses, we reject 9 of them, including the three cases mentioned above. According to our study, the rejected events are not due to microlensing, and this includes 6 cases from the Woźniak et al. (2001) list; at the same time the list includes all 18 of our final binary lens candidates. The single lenses have not been studied extensively. Two events, which we classify as single microlensing are missing in the Woźniak et al. (2001) sample – these are rather poor cases of microlensed light curves. Two events, which we mark as “other objects” belong to the list of Woźniak et al. (2001).

The event SC\_21 6195 has been considered also by Alcock et al. (2000) (their symbol: 97-BLG-d2.) Both models are qualitatively similar, representing intermediate separation binaries and source trajectory approaching one cusp and then crossing the caustic twice near another cusp. Model parameters differ from 7% (impact parameter) to 40% (mass ratio). The parameter defining blending can not be compared directly because of different filters used by OGLE and MACHO teams.

Acknowledgements. We thank Bohdan Paczyński for many helpful discussions. Special thanks are due to Shude Mao for the permission of using his binary lens modeling software, to Przemek Woźniak for the permission to use DIA photometry transient database before its publication, and to OGLE team for their generous early release of data. This work was supported in part by the Polish KBN grant 2-P03D-013-16 and the NSF grant AST 98-20314.

## REFERENCES

* Abt, H.A. 1983, Ann. Rev. Astron. Astrophys., 21, 343.
* Alard, C. 2000, Astron. Astrophys. Suppl. Ser., 144, 363.
* Alard, C., and Lupton, R.H. 1998, Astrophys. J., 503, 325.
* Albrow, M.D., et al.  1999a, Astrophys. J., 512, 672.
* Albrow, M.D., et al.  1999b, Astrophys. J., 522, 1011.
* Albrow, M.D., et al.  1999c, Astrophys. J., 522, 1022.
* Alcock, C., et al.  2000, Astrophys. J., 541, 270.
* An, J.H., et al.  2001, Astrophys. J..submitted, astro-ph/0110095
* Blandford, R.D., and Narayan, R. 1992, Ann. Rev. Astron. Astrophys., 30, 311.
* DiStefano, R., and Mao, S. 1996, Astrophys. J., 457, 93.
* Dominik, M. 1998, Astron. Astrophys., 333, L79.
* Dominik, M. 1999, Astron. Astrophys., 341, 943.
* Jaroszyński, M., and Mao, S. 2001, MNRAS, 325, 1546.
* Mao, S., and DiStefano, R. 1995, Astrophys. J., 440, 22.
* Mao, S., and Loeb, A. 2001, Astrophys. J. Letters, 547, L97.
* Mao, S., and Paczyński, B. 1991, Astrophys. J. Letters, 374, L37.
* Paczyński, B. 1996, Ann. Rev. Astron. Astrophys., 34, 419.
* Press W.H., Teukolsky S.A., Vetterling W.T., and Flannery B.P. 1992, in: “Numerical Recipes in Fortran” (NY: CUP).
* Schneider, P., Ehlers, J., and Falco, E.E. 1992, “Gravitational Lenses”, Springer, Berlin.
* Schneider, P., and Weiss, A. 1986, Astron. Astrophys., 164, 237.
* Trimble, V. 1990, MNRAS, 242, 79.
* Udalski, A. et al.  1994a, Astrophys. J. Letters, 436, L103.
* Udalski, A., Szymański, M., Kaluzny, J., Kubiak, M., Mateo, M., Krzemiński, W., and Paczyński, B. 1994b, Acta Astron., 44, 227.
* Udalski, A., Żebruń, K., Szymański, M., Kubiak, M., Pietrzyński, G., Soszyński, I., and Woźniak, P. 2000, Acta Astron., 50, 1.
* Woźniak, P.R. 2000, Acta Astron., 50, 421.
* Woźniak, P.R., and Paczyński, B. 1997, Astrophys. J., 487, 55.
* Woźniak, P.R., Udalski, A., Szymański, M., Kubiak, M., Pietrzyński, G., Soszyński, I., and Żebruń, K. 2001, Acta Astron., 51, 175.


text:
# Shortest paths and load scaling in scale-free trees

## I Introduction

Recently, many examples have been found of systems whose innate topology is not homogenous and can rather be described in terms of a scale-free, random structure. Examples range from the Internet to cellular metabolism networks. The interest of the physics community in this field stems from the fact that the behavior of many systems on such networks or graphs changes drastically and often attains characteristics close-to but not quite like the mean-field limit.

A scale-free graph consists of a set of nodes or vertices $`V`$ and bonds or edges $`E`$ connecting the vertices to a structure. The essential measure of the scales or lack thereof is the connectivity or degree distribution of the nodes $`V`$: the probability of any node to have $`k`$ edges (one may distinguish between directed and undirected graphs; in the former case the in-going and out-going pdf’s can differ). If this probability $`P_k`$ follows a power-law behavior, a structure arises that does not have any intrinsic scale. The Internet is an example of such a $`P_kk^\gamma `$, and several models have been designed that fit the same description. Later on enhanced models have been devised to capture the characteristics of more elaborate phenomena, like the tendency of clustering tunable\_clustering ; clustered\_scale\_free ; develop\_decay ; heritable\_connectivity .

The models lead to evolving graphs which grow continuously in time by the addition of new nodes, with only a limited number of notable exceptions where the scale-free graph is generated by means of a Monte-Carlo algorithm Krzywicki1 . The degree distribution $`P_k`$ and average connectivity become stationary in the thermodynamic limit, save for the tail of the distribution which is subject to finite-size cutoff effects generic\_scale ; scaling\_evolving . A practically minded question in the same spirit is the growth mechanism of the Internet internet .

It is a common feature of growing networks that they spontaneously develop degree-degree correlations between adjacent nodes stat\_mech ; evolution . This is a manifestation of the preferential attachment BA\_original  principle, where more connected nodes are to attract a larger proportion of new links as the network grows. One recent study hints that the correlation between neighboring node connectivities is the mechanism behind the logarithmic scaling of the network diameter or the average shortest distance between two randomly chosen vertices, with respect to the system size Krzywicki2 . The support for the argument is empirical evidence from simulation results of a broader class of scale-free graph ensembles, where a power law growth of the diameter has been indeed identified. The question of viability of logarithmic scaling in real-world networks is particularly essential, since it has an impact on efficiency and percolation issues (communication over the Internet, spreading phenomena epidemic\_spreading ; epidemic\_dynamics , community structures community\_structure ; local\_events ).

Until recently, less attention has been paid to the probability distribution of shortest path lengths or sometimes referred to as chemical distances in scale-free graphs, possibly owing to the fact that it has been implicitly assumed that the average diameter is an adequate measure of distance properties in the networks. The particular form of the distribution function may have bearings on the performance of search algorithms in scale-free graphs path\_finding . On the other hand, the distribution of shortest paths has been analytically calculated for the small-world model, employing the underlying lattice structure exact\_small\_world  and arriving at a Gaussian-like distribution for large system sizes. Likewise, a model for deterministic scale-free graphs has been proposed and analyzed lately pseudofractal , where a Gaussian is again obtained in the asymptotic limit.

In this paper we focus on a subset of scale-free graphs described by the Barabási-Albert model BA\_original  that in addition are loopless rooted trees from a topological point of view, i.e. the $`m=1`$ case where one connects new nodes by only one link to the existing structure. By removing the redundancy of interconnecting loops it is possible to consider the distance properties on a mean-field level, and also to analyze “load” or “betweenness” universal\_load ; scientific\_coll , the number of any shortest paths passing through vertices. The essential fact here is that the hub of the tree, e.g. the node with the highest connectivity for simplicity, transmits connections between all the branches emanating from it. We show that a stochastic branching process rooted in the preferential attachment rule gives rise to the logarithmic scaling of the diameter and that the pdf of the minimal paths approaches a Gaussian.

Since throughout this text we are interested in tree structures, it is useful to overview their basic features. In the context of random networks, one often resorts to Ref. arbitrary\_degree  and the derivation therein, which suggests that the diameter of graphs grows logarithmically. Although the calculations there are performed for random graphs containing loops, the result obtained closely resembles that for balanced Cayley trees with uniform coordination numbers (except for the coordination number of the central node, which is different). According to this, the number of nodes separated from node zero by $`k`$ nonrecurring steps goes as $`z^k`$ where $`z`$ is the coordination number for the Cayley tree. It then follows simply from the sum of a geometric series that both the longest distance between nodes and the average distance should behave as $`\overline{l}\mathrm{ln}N`$.

It is obvious that trees have unique shortest paths between any two nodes in the sense that without traversing the same edge twice it is not possible to find an alternate minimal route (unlike in unweighted graphs with loops, where there is usually more than one minimal path). We can then define one of the nodes as the root of the tree and unambiguously arrange all the other nodes into layers depending on their minimal distance to the root. Finding the shortest path between two chosen nodes is nothing but identifying the deepest common node along the paths leading from the root to the source and target vertices and then connecting the two nodes via this common fork. Notice that the choice of the root here is slightly arbitrary; one would prefer to use balanced trees.

We study scale-free Barabási-Albert trees BA\_original , starting with a single vertex. In each time step then we add a new vertex with only one outgoing edge. The other end of the edge is connected to one of the nodes already present in the system with a connection probability proportional to the connectivity or degree of a particular node. All edges are thought of as bidirectional and having the same weight, namely 1. As a slight modification to the original model, the connection symmetry of the first two nodes is broken by introducing a “virtual” edge to the very first vertex which only gives preference for this node over the second one when it comes to the subsequent addition of further nodes. This way, we can automatically identify the most connected node in the network and call it its root. To have a balanced tree one needs every subtree of the root to have the same number of nodes in the configurational average. This is only attained when the root is the most connected node, since the BA model ensures that the order of the nodes in terms of connectivity does not change in the course of addition of new nodes and is fully determined by the time of introduction of a node.

## II Mean-field approach

To begin with, we will investigate the shortest path distribution in a mean-field model of a tree network, between the root of the tree and all the other nodes. This argument extends also to general graphs in the case that the new nodes added (e.g. $`m>1`$ Barabási-Albert networks) do not cause a significant amount of shortcuts between already existing nodes.

Let us consider a uniform branching process for each of the layers in the tree so that every node on a certain layer has the same number of offsprings to produce the next layer beneath; it shall amount to $`b(l)`$ for layer $`l`$ for short. This way the original stochastic model is approximated by a deterministic graph pseudofractal ; BRV . The number of nodes $`n(l)`$ with a separation $`l`$ from the root is then $`n(l)=n(l1)b(l1)`$ with the condition that $`n(0)=1`$. The actual form of $`b(l)`$ can be obtained by making use of the preferential attachment rule for BA networks. According to this, the probability that a newly introduced node will connect to any given set of nodes is proportional to the cumulative connectivity of the set in question. Thus, the number of nodes on layer $`l+1`$ changes according to the following rate equation, due to the addition of a new node:

$`{\displaystyle \frac{n(l+1)}{N}}={\displaystyle \frac{1}{2N}}(b(l)+1)n(l),`$ (1)
since the right hand side describes the attachment probability to layer $`l`$, where $`N`$ is the number of nodes in the system and $`2N`$ is the normalization factor for the connectivity. Writing $`n(l+1)=n(l)b(l)`$, expanding the derivation, and dividing by $`n(l)b(l)`$ give

$`{\displaystyle \frac{1}{b(l)}}{\displaystyle \frac{b(l)}{N}}={\displaystyle \frac{1}{2N}}\left({\displaystyle \frac{1}{b(l)}}{\displaystyle \frac{1}{b(l1)}}\right).`$ (2)
If we substitute $`B(N,l)=1/b(l)`$ by explicitly indicating the size dependence on $`N`$ and assume that $`B(N,l)`$ is a slowly varying function,

$`{\displaystyle \frac{\mathrm{ln}B(N,l)}{N}}{\displaystyle \frac{1}{2N}}{\displaystyle \frac{B(N,l)}{l}}.`$ (3)
It is straightforward to expect a solution in the decomposed form of $`B(N,l)=B_l(l)/B_N(N)`$:

$`{\displaystyle \frac{\mathrm{ln}B_N(N)}{N}}{\displaystyle \frac{1}{2N}}{\displaystyle \frac{1}{B_N(N)}}{\displaystyle \frac{B_l(l)}{l}},`$ (4)
and since the left hand side is a function of only $`N`$, $`B_l(l)=2\alpha l`$ with a constant $`\alpha `$. Finally, we get $`B_N(N)=\alpha \mathrm{ln}N`$ and

$`b(l)={\displaystyle \frac{1}{2}}{\displaystyle \frac{\mathrm{ln}N}{l}}.`$ (5)
This relation does not apply to the root ($`l=0`$) for obvious reasons. Eq. (5) also implies that the number of nodes with a given distance to the root $`n(l)`$ keeps growing with $`l`$ until $`b(l)`$ drops below $`1`$ and then starts to decrease, as the bottom of the tree is approached ($`n1`$). This is in strong contrast to the formal prediction of a constant branching for any random graph arbitrary\_degree , which would result in a monotonous exponential growth, as would be the case in usual Cayley trees.

## III Distance scaling

Using the recursion relation $`n(l)=n(l1)b(l1)`$ for the number of nodes on a given level, we can give an estimate for the shortest path distribution function with the source of the paths at the root of the tree. Instead of Eq. (5) we take now the more general form of $`b(l)=(A/l)^\lambda `$ and approximate the sum with an integral in the following expression, which reads

$`n(l)=n(0){\displaystyle \underset{i=0}{\overset{l1}{}}}b(i)=b(0)\mathrm{exp}\left[\mathrm{ln}{\displaystyle \underset{i=1}{\overset{l1}{}}}b(i)\right]`$ (6)
$`b(0)A^{\lambda (l1)}\mathrm{exp}\left(\lambda {\displaystyle _1^{l1}}\mathrm{ln}xdx\right)=`$
$`={\displaystyle \frac{b(0)}{e^\lambda }}\left({\displaystyle \frac{Ae}{l1}}\right)^{\lambda (l1)}.`$
The result above for $`n(l)`$ approaches a non-normalized Gaussian in the large-$`N`$ limit as $`A\mathrm{ln}N`$ goes to infinity, which can be seen from Fig. 1 where a $`C`$ corresponding to a very large network has been used. In order to draw further conclusions, we will determine the parameters of the Gaussian that give a best fit to $`n(l)`$. For the sake of simplicity, let us now consider the function of the form $`f(x)=(C/x)^{\lambda x}`$:

$`f(x)=\left({\displaystyle \frac{C}{x}}\right)^{\lambda x}R\mathrm{exp}\left[{\displaystyle \frac{(x\mu )^2}{2\sigma ^2}}\right].`$ (7)
We first match the extremal point of $`f`$ to the mean of the Gaussian, resulting in $`\mu =C/e`$. The maximum value is thus $`R=\mathrm{exp}(C\lambda /e)`$; the standard deviation $`\sigma `$ can be obtained by the requirement that the derivative functions of $`f`$ and the Gaussian be the same in the vicinity of $`\mu `$ up to first order, giving $`\sigma =\sqrt{C/(\lambda e)}`$. Using the parameters acquired this way, we can find a very good approximation to $`f(x)`$, which is almost identical to that of a least-square fit.

Furthermore, additional information can be gained if we look into the normalization conditions for $`n(l)`$. Trivially, the sum of $`n(l)`$ over all layers should return the total number of nodes in the system, $`N`$. Again, we approximate the sum with an integral:

$`N`$ $`=`$ $`{\displaystyle \underset{l=1}{\overset{\mathrm{}}{}}}n(l){\displaystyle _1^{\mathrm{}}}n(x)𝑑x={\displaystyle \frac{b(0)}{e^\lambda }}{\displaystyle _0^{\mathrm{}}}\left({\displaystyle \frac{Ae}{x}}\right)^{\lambda x}𝑑x`$ (8)
$``$ $`{\displaystyle \frac{b(0)}{e^\lambda }}{\displaystyle _0^{\mathrm{}}}e^{A\lambda }\mathrm{exp}\left[{\displaystyle \frac{(xA)^2}{2A/\lambda }}\right]𝑑x`$
$``$ $`{\displaystyle \frac{b(0)}{e^\lambda }}e^{A\lambda }\sqrt{{\displaystyle \frac{2\pi A}{\lambda }}},`$
where we assumed that $`A`$ is large enough so that we can neglect the correction of the error function to the Gaussian integral. We should also be aware that $`l`$ has a finite cutoff because of the bounded depth of the tree—yet, the quickly vanishing $`n(l)`$ makes it possible to take $`l`$ to infinity. Finally, $`f(x)1`$ as $`x0`$ so that the integrand is bounded everywhere.

Recall now that the degree of a node in BA networks grows with the power of $`N`$, $`b(0)N^\beta `$ stat\_mech . Apart from $`b(0)`$, the only term on the right hand side of Eq. (8) that may contribute to the overall linear growth in $`N`$ is $`e^{A\lambda }`$, which increases much faster than $`\sqrt{A}`$, so the latter can be taken a constant. The consistency condition with the left hand side requires that $`e^{A\lambda }N^{1\beta }`$ should hold, and thus

$`A={\displaystyle \frac{1\beta }{\lambda }}\mathrm{ln}N+\text{const}.`$ (9)
Disregarding the constant, we end up with a very similar but more general expression than that of Eq. (5) for $`b(l)`$:

$`b(l)=\left({\displaystyle \frac{\upsilon \mathrm{ln}N}{l}}\right)^\lambda ,\text{with}\beta +\upsilon \lambda =1.`$ (10)
This implies that if a scale-free tree is characterized by a branching process decaying as a power law as a function of the distance from a suitable root node with the highest connectivity, the relation (10) should necessarily be satisfied. Not surprisingly, it is true in the case of BA trees, where $`\beta =1/2`$ and according to Eq. (5), $`\upsilon =1/2`$ and $`\lambda =1`$. One should note that in the process of constructing the mapping we rely on the fact that the number of nodes in a layer depends only on the average branching ratio $`b(l)`$; the fluctuations in the degrees of nodes is omitted. For this reason the degree distribution exponent $`3`$ is not present in the tree representation.

The node-to-node distances in the mean-field model are calculated as follows. We traverse each node of the tree and enumerate the routes with certain lengths that start at or go through this node and have both of their ends in the subtree of the node. Practically speaking, we can think of this node as the root for its subtree and perform the same calculations as we would do for the ’global root’ of the tree. If we by $`n^{(s)}(l)`$ denote the number of possible paths going out to the subtree of a node on level $`s`$ that have length $`l`$ and one end fixed at the node on level $`s`$,

$`n^{(s)}(l)={\displaystyle \underset{i=0}{\overset{l1}{}}}b(s+i).`$ (11)
Let now $`r^{(s)}(l)`$ be the number of all routes that go through or end at a particular node on level $`s`$ and have a length of $`l`$,

$`r^{(s)}(l)=`$
$`=n^{(s)}(l)+\mathrm{\Theta }(b(s)1){\displaystyle \underset{i=1}{\overset{l1}{}}}\left({\displaystyle \genfrac{}{}{0pt}{}{b(s)}{2}}\right){\displaystyle \frac{n^{(s)}(i)}{b(s)}}{\displaystyle \frac{n^{(s)}(li)}{b(s)}}=`$
$`=n^{(s)}(l)+\mathrm{\Theta }(b(s)1){\displaystyle \frac{b(s)1}{2b(s)}}{\displaystyle \underset{i=1}{\overset{l1}{}}}n^{(s)}(i)n^{(s)}(li).`$ (12)
The second term in the sum has contribution to $`r^{(s)}(l)`$ only when there are branches left going out from a node, in average when $`b(s)1`$. The number of paths with a specific length in the whole system is therefore

$`r(l)={\displaystyle \underset{s=0}{\overset{L}{}}}r^{(s)}(l)n(s),`$ (13)
where $`n(s)`$ is defined as earlier in Eq. (6), the number of nodes on a given level $`s`$.

The Barabási-Albert model allows for more rigorous derivations of the relation for $`n(l)`$. The mathematics community often refers to the tree interpretation of the model as *recursive trees*, and thus exact results have been obtained for both the distance distribution and the diameter of the trees pittel\_rec\_trees ; poisson\_approx\_trees ; distr\_recursive\_trees . Bollobás and Riordan give a general proof for the diameter scaling of scale-free BA graphs bollobas\_riordan . The mapping to Cayley trees also resembles the work of Krapivsky and Redner, who arrive at a closed recursive analytical form for $`n(l)`$, in a more general context than that of scale-free trees organization\_growing\_nets . It also resembles Cayley models of Internet traceroutes fractal\_internet  by Caldarelli and co-workers.

## IV Comparison with simulations

Numerical simulations of BA scale-free trees fully confirm the inferences drawn in the preceding section. Most important of all, the average number of branches per node on a given level is shown in Fig. 2. The numerical parameters of the power-law fit conform with the mean-field values: the exponent of the decay is almost exactly $`1`$, and the prefactor of the logarithm with $`0.43`$ is also close to that of the predicted $`1/2`$. It is also worth to note that if we rescale the distance variable by the logarithm of the system size, we can attain a data collapse with a very good accuracy. This means that for BA trees in practice $`b(l)`$ can be approximated as

$`b(l)=\{\begin{array}{cc}\frac{0.43\mathrm{ln}N}{l}\hfill & \text{if }lL(N)\hfill \\ 0\hfill & \text{otherwise}\hfill \end{array}`$ (16)
From the inset of Fig. 2 it is also apparent that the cutoff $`L(N)`$ is a little over the value of $`\mathrm{ln}N`$, by a factor of about $`1.3`$. On the other hand, the drop of $`b(l)`$ at $`L(N)`$ is measured to be either an exponential or a power-law with a very large exponent. The mean-field prediction for the maximum of the shortest path length, $`L(N)`$, can be obtained by equating $`n(L)=1`$ in Eq. (6) and using the Gaussian approximation of Eq. (7). The solution up to first order in $`\mathrm{ln}N`$ is that $`L(N)\frac{1+\sqrt{2}}{2}\mathrm{ln}N`$, which again is in reasonable agreement with the mean-field argument.

The derived quantities $`n(l)/N`$ and the node-to-node distance distribution is shown in Fig. 3 for two distinct cases:

* the root-to-node and node-to-node shortest path distribution is measured in an ensemble of random BA trees using simulations. Instead of every possible pair, the node-to-node distances are measured only between a large but finite number of randomly selected vertex pairs, for practical reasons;
* both distribution functions are estimated by utilizing the mean-field tree mapping, using the asymptotic form of Eq. (16) for $`b(l)`$ with a cutoff at $`L=1.2\mathrm{ln}N`$.

It is to be seen that a very good correspondence is found between the root-node distribution functions, but the overall two-point pdf’s are sensibly close as well.

While it has been relatively easy to derive analytical results for the root-node distances in the mean-field trees, Eq. (13) and the quantities it is constructed of turn out to be too complex to handle without numerics. The formulas (11)–(13) above are used to calculate the approximate values of the node-to-node path length distribution in the mean-field trees using the expression of Eq. (16) instead of the analytical form of Eq. (5), so as to better represent the random BA trees. It is reassuring that the generic form of the node-to-node distance pdf also follows a $`(C^{}/x)^x`$ function, only with a different $`C^{}`$ constant from that of $`C`$ for the root-node distances \[Eq. (6)–(7)\]; see Fig. 3. The diameter of the trees relative to the logarithm of the system size can be seen on Fig. 4. $`l\mathrm{ln}N`$ or in other words twice the mean of root-to-node distances. This is somewhat expected as the main contribution to the node-to-node paths arises from passing through the root, for large graphs. It leads to a convolution-type distribution (from the two ’legs’). It can easily be seen that the diameter cannot exceed twice the depth of the tree, which gives rise to a logarithmic growth in any case.

## V On the “load” on trees

On a hierarchical structure the total number of minimum paths going through a node (the “load”) can be divided into two contributions. First, those paths that connect nodes in separate sub-branches of the node to each other, and, second, those that connect the nodes belonging to the branches to the rest of the tree. Call $`d(l)`$ the number of the descendants of a node on level $`l`$. In other words, $`d(l)`$ is the size of the subtree for the node. Then the load can be written simply as

$$\mathrm{\Lambda }(l)=\left(\genfrac{}{}{0pt}{}{b(l)}{2}\right)\left[\frac{d(l)}{b(l)}\right]^2+d(l)\times (Nd(l)),$$
(17)
where the last term counts the connections towards the hub.

For the particular example we are concerned with, it is easy to see that the latter term dominates ($`Nd(l)`$) and moreover that a good approximation is given by just simply $`\mathrm{\Lambda }(l)d(l)N`$. Thus one may investigate the dependence of the load on the level (or depth) of the tree, $`l`$. For $`d(l)`$ in the mean-field picture one has that

$`d(l)={\displaystyle \frac{_{i=l+1}^Ln(i)}{n(l)}},`$ (18)
and for the layer immediately below

$`d(l+1)`$ $`=`$ $`{\displaystyle \frac{_{i=l+2}^Ln(i)}{n(l+1)}}={\displaystyle \frac{d(l)n(l)n(l+1)}{n(l+1)}}=`$ (19)
$`=`$ $`{\displaystyle \frac{d(l)}{b(l)}}1{\displaystyle \frac{d(l)}{b(l)}},`$
where we also used the recursion relation for $`n(l+1)`$. Finally, the load changes for the layer underneath as

$`\mathrm{\Lambda }(l+1)=d(l+1)N={\displaystyle \frac{d(l)}{b(l)}}N={\displaystyle \frac{\mathrm{\Lambda }(l)}{b(l)}}.`$ (20)
Since the load $`\mathrm{\Lambda }(l)`$ is the same for each of the nodes on a particular level $`l`$, the distance–load distribution is directly given by the normalized $`n`$$`\mathrm{\Lambda }`$ function, thus hiding the implicit dependence on $`l`$. Considering that

$`\mathrm{\Lambda }(l+1)`$ $`=`$ $`{\displaystyle \frac{\mathrm{\Lambda }(l)}{b(l)}}`$
$`n(l+1)`$ $`=`$ $`n(l)b(l),`$ (21)
$`\mathrm{\Lambda }(l+1)n(l+1)=\mathrm{\Lambda }(l)n(l)=\text{const.}`$ and therefore

$`n={\displaystyle \frac{\text{const.}}{\mathrm{\Lambda }}}.`$ (22)
We then expect to see that the load is inversely proportional to the number of nodes on the levels, which is indeed the case according to Fig. 5. The same result holds for normal Cayley trees from Eq. (17).

Note that we have to use mean-field trees which would correspond to random networks with a large number of nodes so that the number of levels is of the order of 10. For the load distribution we consider only levels for which $`b(l)1`$ because otherwise subtrees do not exist in the average sense. It is surprising that the load distribution exponent does not depend on the actual form of $`b(l)`$, being universally $`1`$ \[Eq. (21)\]. Indeed, *the exponent of the distance–load pdf is independent of the choice of the node that all the distances are taken relative to*.

Another common way of defining the importance of the nodes in terms of shortest paths passing through them is the one called betweenness, favorable for its algorithmic feasibility and simplicity. Newman presents a breadth first search algorithm for efficient calculation of the betweenness of nodes on random graphs scientific\_coll . The only notable difference to Eq. (17) comes from the fact that the betweenness also accounts for paths that originate from the nodes themselves, which nevertheless amounts only to a constant system size.

We will calculate the betweenness on the trees, now focusing as a goal on the probability distribution of the load. An estimation can be given for a node by considering the contributions to it, and by separating the network to a descendants part with $`d`$ nodes in the branches and all the rest with $`N1d`$ nodes. The node being the source, we have $`N`$ shortest paths to any other node; if the source is among the descendants, we have $`d(N1d+1)`$ ones going through; if the source is any other node from the network, we have $`(N1d)(d+1)`$. A fourth contribution, coming from paths passing through the node but having both ends in the descendant tree, has been neglected. They add up to an estimated betweenness $`K`$ of

$`K(d)2N1+2(N1)d2d^2.`$ (23)
Here it is to be seen that for small $`d`$’s the linear term dominates, just as in our previous load calculation; Fig. 6 justifies our estimations.

The betweenness probability distribution $`P(K)`$ taken over all nodes in the network can then be concluded to asymptotically follow a power law decay with a universal exponent of $`2`$. This is since $`K`$ is linear in the number of descendants $`d`$ and moreover that the pdf of $`d`$ scales universally with an exponent of $`2`$ for supercritical trees supercritical\_PL . Strictly speaking, the conclusions here are only true for the supercritical part of the tree, i.e. where $`b(l)>1`$. The subcritical leaves of the tree have an increasingly smaller number of descendants, though, which drop exponentially with each new layer, and it can be verified that the descendant pdf decay exponent is indeed above 2 if only this part of the tree is considered. Nevertheless, Fig. 7 shows that both the descendant pdf and the load pdf are accurately described by inverse square functions. A scaling of the load distribution has been experimentally found on other scale-free networks as well universal\_load , only with a slightly different universal exponent of about $`2.2`$.

A further, practically more far-reaching observation is that the average betweenness as measured as a function of the locally known node degree grows as a power-law of the degree with an exponent of about $`1.8`$ (Fig. 8). A mean-field approach can be used to estimate for the exponent, though, if we consider that the preferential attachment principle for large degrees gives rise to a descendant–degree scaling of $`dk^{\frac{1}{\beta }}`$ ($`\beta =1/2`$), which is the inverted relation for the time evolution of the degree of a parent node stat\_mech . In this particular case, time is measured as the size of the node’s subtree. A substitution of the latter into the linear load equation would suggest an exponent of $`2`$; the deviation from it may come from the rather restricted range of the degree that the relatively small system sizes allow.

## VI Conclusions

In this paper we have mapped scale-free Barabási-Albert trees to a deterministic model of a rooted tree with a uniform branching process on each layer of the tree. This idea resembles work on the Internet structure fractal\_internet  and on the structure of branched cracks, where an inverse relation of the branching to distance has been observed bouchaud .

Simulations show that the distribution of the number of branches on one particular layer of the tree follows a power-law function, but it turns out to be a good approximation to describe the branching only by its mean, $`b(l)`$. In the simple case of BA trees it can be shown by means of this mapping that the diameter of the networks is bounded by the logarithm of the network size and the asymptotic form of the distance distribution functions follows immediately. In other words, we can examine the slow convergence of this function to the limiting Gaussian form for infinite system sizes. Given an effective description in terms of a tree plus a branching process, further information can be found, e.g. one may consider the scaling of the number of shortest-distance paths (load, or betweenness). Non-uniform critical trees could perhaps be constructed in a self-organized fashion, as is possible for the statistically uniform case SO\_branching .

One should note the close relation of the Cayley representation to minimal spanning trees (MST) on scale-free (random) networks; for the $`m=1`$-networks these two coincide. This makes it an interesting prospect to study the load and distance properties of MST’s in other scale-free networks.

We thank Bálint Tóth for useful discussions.

This work has been supported by the Academy of Finland’s Centre of Excellence Programme and by OTKA T029985.