Statistical Properties of Pairwise Distances between Leaves on a Random Yule Tree

A Yule tree is the result of a branching process with constant birth and death rates. Such a process serves as an instructive null model of many empirical systems, for instance, the evolution of species leading to a phylogenetic tree. However, often in phylogeny the only available information is the pairwise distances between a small fraction of extant species representing the leaves of the tree. In this article we study statistical properties of the pairwise distances in a Yule tree. Using a method based on a recursion, we derive an exact, analytic and compact formula for the expected number of pairs separated by a certain time distance. This number turns out to follow a increasing exponential function. This property of a Yule tree can serve as a simple test for empirical data to be well described by a Yule process. We further use this recursive method to calculate the expected number of the n-most closely related pairs of leaves and the number of cherries separated by a certain time distance. To make our results more useful for realistic scenarios, we explicitly take into account that the leaves of a tree may be incompletely sampled and derive a criterion for poorly sampled phylogenies. We show that our result can account for empirical data, using two families of birds species.


Introduction
The speciation process in evolution can be regarded as a branching process. One of the simplest stochastic models for a branching process is the so called Yule process [1,2]. In this model branches are assumed to split with a constant rate and both resulting branches will evolve independently in time. Starting from one branch, a tree will grow, such that the number of leaves on average increases exponentially in time. In a more general version of the Yule tree each branch can also die and get extinct with a constant rate.
Despite its simplicity, many phenomena in different fields of science have been successfully modeled using the Yule process [3,4]. Particular examples include statistical properties of the number of species in a genus [1], the number of members in protein and gene families [5,6] and phoneme frequencies in languages [7]. In stochastic modelling of biological evolution, the Yule process is often useful as an instructive null hypothesis [8][9][10][11], even when its assumptions are clearly violated.
As an illustrative example of the branching process we present the reconstructed phylogenetic tree of species in the Siilvidae family of birds in the left panel of Fig. 1. The basis of such a reconstructed tree is pairwise distances between individual species. The color-coded matrix of such distances for the species is shown in the right panel of Fig. 1. The statistical properties of such a matrix for a Yule tree is the focus of our article.
Statistical properties of Yule trees have been intensively studied and much is already known. One of the most useful results is the distribution of the number of leaves on a Yule tree [12]. This exact analytical result is widely exploited, in particular, for reconstruction of phylogenetic trees and for estimation of rates of speciation and extinction [10,11,13]. Other discrete properties have been studied in Refs. [14][15][16][17] as well as properties of the distribution of branch lengths [18,19].
Often the pairwise distances between all pairs of species in a group of species is the only available information useful for reconstruction of the evolutionary history of the group. For example, in phylogeny reconstruction, one can estimate the pairwise distance in time between two species (twice the time to their last common ancestor) using the molecular clock approach, together with morphological considerations and information about the fossil record [20]. Motivated by observations of mitochondrial DNA sequences with no recombination, the distribution of pairwise distances has been studied in Ref. [21] for a tree with discrete generations and a given number of leaves. In this study, the authors use a sort of mean-field approach, ignoring fluctuations in the number of leaves during the growth of the tree, to derive an approximate formula for the pairwise distances distribution on a tree.
Here we present a general method to derive the distribution of pairwise distances and other statistical properties on a continuous random Yule tree of a certain height with given birth and death rates. Using our method, we obtain exact, analytic, closed, non-recursive and compact formulas for the pairwise distance distribution, the distribution of distances to the closest neighbour, the distance distribution in so-called cherries, as well as a more general formula for the distribution distance to the n-th closest neighbour. One of the reconstructed trees for the Siilvidae family of species, taken from [28] (left) and its distance matrix (right). The tree includes only the branches which lead to survived and observed leaves.
Often, in biological context, one does not have an access to data about all existing species (i.e. leaves of a phylogenetic tree) [22]. Instead, species are incompletely sampled, or might have been subject to a recent massive extinction event [23]. As long as the extinction of species is random, both scenarios are equivalent on macroevolutionary timescales. In our study, we take the incomplete sampling explicitly into account, which allows us to make statements about the fraction of sampled species, using only the available data.
In the next section we will start with a formal definition of the Yule process and then derive the above mentioned distributions of pairwise distances. For illustrative purposes we also present numerical simulations perfectly matching our expectations. At the end of our article we apply our theoretical consideration to empirical data and analyze the speciation process in two families of birds for which data on speciation times and pairwise distances is available. One advantage of our approach is that we do not need to reconstruct a phylogenetic tree but can solely work with data on pairwise distances.

A Yule tree with constant branching and extinction rates and incomplete sampling of leaves Definition of the Yule Tree
A Yule tree is defined as follows [1,2]. At time t = 0 there is one individual. As time progresses, this individual can branch and give birth to another individual. In an infinitesimally short time interval [t, t+dt], all individuals can give birth to another one, each with the probability λdt. The probability of an individual to die in the same time interval is μdt. We consider an ensemble of trees of age (height) T, referring to all existing individuals at this time as leaves. To make the model more realistic, we assume that due to incomplete sampling (or a short massive extinction event) just before the time T, each leaf is observed with a certain probability 0 σ 1. The described process is illustrated in Fig. 2. We assume that the incompleteness of the sampling is random and ignore possible biases due to different sampling schemes [24].

A Few Useful Results for Random Trees Generated by a Yule Process
Consider a Yule tree with birth rate λ and death rate μ, that have been grown for total time (height) T. In the case where all leaves are sampled (σ = 1), let P(MjT, σ = 1) be the probability that there are M leaves on a tree of age T. Following [25], we can then write the probability that no individual (M = 0) survives through to time T as For M > 0 we have We can derive corresponding equations also for the case where species are sampled incompletely. In this case, the probability that no species is observed is m 0 s 0 ð1 À sÞ mÀ0 PðmjT; s ¼ 1Þ ¼ e mT ðm À l þ slÞ À e lT ms e mT ðm À l þ slÞ À e lT ls ð3Þ Despite these complicated expressions, the average number of observed leaves in a tree of age T is simply given by and the average total number of pairs is However, the number of observed leaves is 2 (leaves 1 and 3) for the left subtree and 3 (leaves 5, 7 and 8) for the right one. The thick green line denotes the pairwise evolutionary distance between the two observed leaves 5 and 7. The horizontal dimension is meaningless. In this example for leaf 1 the first closest observed leaf is 3, the second (as well as the third and the fourth) is 5 (or 7 or 8). The tree has two observed cherry pairs: (1, 3) and (7,8).
The total length of all branches in a Yule tree is given by the integral: To derive a corresponding expression for a a tree reconstructed only from incompletely sampled leaves, we note that the average number of branches at time t with at least one observed descendant at time T is given by In the case where t = T, we have that hM(T, T)i = σhM(T)i. The average total branch length on the tree of length T excluding the branches which do not lead to an observed leaf is then given by Z T 0 hMðt; TÞi dt ¼ se TðlÀmÞ m À l þ sl ln ls þ ðl À sl À mÞe TðmÀlÞ l À m : ð9Þ In the limit of no extinction, μ ! 0, and exhaustive sampling, σ ! 1, Equation (9) is identical to Equation (7). We turn now to calculations of the statistical properties of pairwise distances, using the above formulas.

The Distribution of Pairwise Distances
In a biological context the available data often consist of the pairwise distances separating any pair in a group of species. Commonly these distances are used to reconstruct a phylogenetic tree representing the evolutionary history of a group of species. From such a tree one can then try to estimate rates of speciation and extinction [10,11]. Here we propose another approach of analysing such data on pairwise distances circumventing the reconstruction of a phylogenetic tree, provided that the pairwise distances between the leaves are properly estimated. Let N(tjT)dt be the average number of pairs of leaves on a tree of length (evolution time) T, separated by a time distance in the interval [t, t+dt], i.e. their last common ancestor lived in the time interval [T−t/2−dt/2, T−t/2]. Now consider the branching process as illustrated in Fig. 2. The first branching happened at time T 1 and the two resulting subtrees encompass, say, M 1 and M 2 leaves, respectively. In this situation one can derive the following recursion relation where the first part in the summation on the right hand side counts the pairs inside each of the two subtrees and the second one counts the pairs between them. The common multiplicative factor, e ÀmT 1 , expresses the probability that the first branch survives to the time T 1 (otherwise, N (tjT) = 0). The function I is the indicator function, defined by: & and δ(x) is the Dirac delta function. Averaging over M 1 , M 2 (using Equations (3,4) with time T −T 1 ) and then T 1 , which follows an exponential distribution with mean 1/λ, one obtains: In Laplace space one gets: where S is the Laplace conjugate variable of T. Solving and inverting the Laplace transform one finally gets the solution: for 0 t 2T and zero otherwise. Fascinatingly, this distribution is a simple exponential function in t. The distribution is cut off at t = 2T because in a tree of age T two leaves cannot be separated by a time larger than 2T. In Fig. 3(a) we show this distribution of pairwise distances for several parameter values together with results of numerical simulations, which match perfectly our theoretical expectations. This result, applied for trees of DNA sequences can account for statistics of exact sequence matches in genomes of eukaryotes [26]. One can also derive the same result (14) using the following simple arguments. Pairs, separated by a time in the interval [t, t+dt], branched at the time interval [T−t/2−dt/2, T−t/2]. The average number of branches in this interval is given by λe (λ−μ)(T−t/2) dt/2. The average number of observed pairs from a branch at this time is given by (σe (λ−μ)t/2 ) 2 . Multiplying the two factors one gets Equation (14). However, for other quantities, derived below, the recursive equation approach is more effective.

The Distribution of the Minimal-Distance to Other Leaves
Using the recursive method from the previous Section one can also compute other interesting quantities. For instance, in certain situations, the distance separating a leaf to its most closely relative may be estimated more precisely than its distance to other leaves in the tree. Thus, we might be interested in N 1 (tjT)dt-the average number of leaves on the tree of age T, separated by the time distance between t and t+dt from their most closely related leaf. Interestingly, calculating this quantity lets us make certain statements on the value of the sampling rate σ.
To calculate this distribution, we can again write a recursion relation, assuming that the first branching occurred at time T 1 . In this case one gets the distribution of the minimal distance time in the form where P(MjT) is the probability to observe M leaves after time T, as computed in Equations (3) and (4). In contrast to the recursion relation for the distribution of all pairwise distances, we count a branching point only if Equation (15).
Averaging Equation (15) over T 1 , one gets: The solution of this equation is given by for 0 t 2T and 0 otherwise. Results of numerical simulations perfectly match our theoretical expectations (see Fig. 3(b)). Interestingly, the function N 1 (tjT) from Equation (17) possesses a maximum only if and the position of the maximum is in the range [0, 2T]. This result is useful for a quick estimation of the data completeness. In particular, a maximum in the distribution of the minimal distance implies that the sampling of the considered tree is not complete and σ < 1/3. By similar arguments we can also derive expressions for the distributions of second minimal distances, N 2 (tjT) (see Appendix) and of the n-th minimal distance N n (tjT) (see Appendix) to other leaves. The latter quantity is computed to be for 0 t 2T and 0 otherwise. In Appendix we also calculate the distribution of distances in "cherries". Cherries are adjacent pairs of leaves, such that they are reciprocal closest neighbors to each other (see Fig. 2 for illustration of cherries): is in the range [0, 2T]. This result is useful for a quick estimation of the data completeness. In particular, a maximum in the distribution of the distance between cherries implies that the sampling of the considered tree is not complete and σ < 1/4. For illustration purposes we show the distributions for the second minimal distance in Fig. 3(c) and, for cherries, in Fig. 3(d).

Beyond the Averages
Above results are average expectations. For instance, in The Distribution of Pairwise Distances Section we derive N(tjT), defined as the average density number of pairs, separated by a certain time distance t, on a tree of length T. The average is over many realizations, say S many, of the Yule trees with a given set of parameters λ, μ, σ and T. Namely, where N s (tjT) is the density number of pairs separated by a time distance in the interval [t, t +dt] in an individual sample tree number s. In reality one often possesses information only about one specific tree s = 1, i.e. N 1 (tjT). Therefore, we are interested not only in the derived averages of N(tjT), N n (tjT), N Λ (tjT) etc. but also their distributions in finite time intervals. The last becomes especially important in the maximum likelihood fitting and model testing. In the discussion below we refer to the distribution of the number of pairs separated by a certain time, N 1 (tjT). However, the same arguments can be applied to other quantities, like the n-th minimal distance or the distance in cherries, which we mention above.
Consider an infinitesimal (in practice very small) interval, [t, t+dt], such that N(tjT)dt ( 1. The number of pairs N 1 (tjT)dt in this interval is distributed with the mean N(tjT)dt. However, in the considered small bin limit, the mean does not represent well the typical value because the distribution of N 1 (tjT)dt is not well peaked but possesses a very small probability of having any positive value, while probability of having zero is almost one (see Appendix).

Pairs separated by the time in the interval [t, t+dt] branched at the time interval [T−t/2−dt/ 2, T−t/2]. The probability to have a branch in this interval is given by λe (λ−μ)(T−t/2) dt/2.
Given that there is a branching point in this interval it can lead to different number of leaves. The probability that no observed pairs survive from this branching is given by 1−[1−P(0jt/2)] 2 , where P(MjT) is the probability to observe M leaves on a tree of age T and is given in Equations  (3, 4). Therefore, the probability that there are no observed pairs separated by the time in the interval [t, t+dt] is given by In sum, in the small bin limit it is convenient to break the full distribution in two distributions: One comprising only the peak at zero and a second representing all samples with N 1 (tjT) dt 6 ¼ 0. The total average can be broken as follow: HereÑ ðt j TÞ is the average of N 1 (tjT) over the tree realizations with N 1 (tjT) > 0. It can be computed to be: whereSðtÞ ¼ P S s¼1 1 À d N s ðtjTÞ;0 h i is the number of samples with N 1 (tjT) > 0. Since, 1−Pr (N 1 (tjT)dt = 0) ( 1, the value of N(tjT)dt is not representative of the expected empirical average of N 1 (tjT)dt for finite S and, in particular, S = 1. However, the value ofÑ ðt j TÞ, derived above (see Equation (27)), is representative of the expected empirical average of positive values of N s (tjT)dt. We illustrate this in Fig. 4 Constrains on the sampling fraction One can easily see that all the derived above results do not depend explicitly on the parameters λ, μ and σ, but only on their combinations: λ−μ and σλ. Therefore, one cannot estimate the sampling fraction, σ, based on fitting the empirical data to the derived formulas (see examples in the next Section). The same loss of information in reconstructed trees was reported, based on an analysis of the density of bifurcation times in the reconstructed tree [27]. However, the information about the values of λ, μ and, most intriguingly, σ is not lost completely. For instance, observing a maximum in the distribution of the minimal distances one can deduce that σ < 1/3 (see Equation (18)). Observing a maximum in the distribution of the distances between cherries one can deduce that σ < 1/4 (see Equation (23)). It is of an interest to construct other distributions which, possessing a maximum, provide information about the value of the sampling fraction, σ.
Consider an average density of pairs of leaves with the following property. Given that the first (second) leaf of the pair has a nearest neighbor at a distance (if a leaf is alone in the tree we define the distance to its nearest neighbor as twice the height of the tree) t 1 (t 2 ) the quantity min(t 1 , t 2 ) is given by t. We denote this density by N min2 (tjT). The recursive equation for this quantity is given for a given time of first bifurcation, T 1 by After average over T 1 the solution is given by N min2 ðtjTÞ ¼ 2l 2 s 3 ðl À mÞ 4 e TðlÀmÞ 3l À m 2le tðlÀmÞ þ ðl À mÞe À 1 2 tðlþmÞ þ ðm À 3lÞe 1 2 tðlÀmÞ e TðlÀmÞ ls À l þ m À lse This function possesses a maximum only if Therefore, observing a maximum in the distribution of the minimal distance to the closest neighbors between two leaves one can deduce that σ < 1/5. Using our recursive method one can calculate different distributions (say, the minimal distance to the closest neighbor among three leaves etc.) which, exhibiting a maximum, provide direct information about an upper limit on the sampling fraction.

Comparison of the derived results to empirical data
In this Section we demonstrate the relevance of the obtained analytic formulas to empirical data, studying the pairwise distances between species in families of the evolutionary tree. For comparison with the derived results we choose N(tjT), N n (tjT) with n = 1, 2, 3, 4 and N Λ (tjT). The results are presented in Fig. 5 for the Siilvidae family of birds (see one of the reconstructed trees for this family and its distance matrix in Fig. 1) and for the Tyrannidae family of birds in  (14) and (27) 6. For every family we analyze Bayesian sampling of 1000 trees downloaded from the database [28]. Namely, we collect pairwise distances, n-minimal distances and distances between cherries of all 1000 trees and plot the histograms of these distances (with the y-axis divided by 1000) in Figs. 5 and 6. We fit all the points in a figure using the iterative reweighted least squares algorithm [29] in Matlab. Unfortunately, the explicit dependencies on λ and μ in Equations (14,20,21) are insufficient to estimate all parameters. Instead one can estimate from the fit only the effective growth rate, λ−μ and λσ. The value of σ can be obtained assuming a certain ratio μ/λ. In the captions of Figs. 5 and 6 we present the obtained estimates for σ for different assumptions about the ratio μ/λ. Over all, the fits to empirical data look satisfactory and result in a reasonable set of parameters, which roughly agree with the ones given in [28]. This indicates that certain statistical properties of speciation can be well captured by a simple Yule process. However, in some cases, deviations can be observed. For example, for the Sylviidae family the pairwise distances distribution deviates from the prediction for t > 30 Myr, while for the Tyrannidae family we observe a clear deviation for distances around 55 Myr in all our estimates. This indicates a massive radiation event in the considered family of birds around 27.5 Myr ago, as already reported in [28], or other violation of the Yule process assumptions.
Interestingly, we can state that the Sylviidae family of birds is currently not well sampled. In fact, the estimator for the upper limit of the sampling fraction σ is 30% (see Fig. 5).

Summary and concluding remarks
In this paper we present a novel method to calculate statistical properties of Yule trees. The method is based on a recursive equations which can be solved using the Laplace transform. We demonstrate the strength of our method deriving formulas for (i) average number of pairs separated by a certain time (Equation (14)), (ii) the number of most closely related pairs separated by a certain time (Equation (17)), (iii) the number of next-most closely related pairs separated by a certain time (Equation (33)), (iv) the number of n-most closely related pairs separated by a certain time (Equation (20)) and (v) the number of cherries separated by a certain time (Equation (21)).
Our results can be compared to empirical data using only the information about pairwise distances between leaves of a considered tree. We assume that the estimation of the pairwise distances is precise enough. If the distances are estimated using genetic divergence, this assume that the molecular clock reflect adequately the real time distance. If this holds the reconstruction of the tree structure is not required. This is a particular strength of our method because the reconstruction of such trees for a large number of leaves is sometimes problematic. In such cases one often considered a posterior distribution of trees which is generated by Bayesian sampling [30,31]. Such a distribution of trees can still be easily analyzed using our method, based on recursive equations. Analyzing such ensembles of trees we use only their distance matrices.
We demonstrate the relevance of our results to statistical properties of pairwise evolutionary time distances between biological species. We find that in some cases the speciation process is well described by the Yule model. Significant deviations from the derived distributions are expected to be indicative for massive extinction or radiation events. In the case where the assumptions of the Yule process are justified, we expect our results to be useful for estimation of the incompleteness of the data sampling, i.e. the fraction of observed leaves out of all existing leaves, σ. However, similarly to the method developed in Ref. [11], all the derived results depend only on three parameters: λ−μ, λσ and σe (λ−μ)T . Therefore, even knowing those three parameters one cannot estimate the values of the four unknown parameters: the rates λ, μ, the height of the tree, T and the sampling fraction, σ, without an additional assumption about one of these parameters, for instance the fraction μ/λ. After estimation of (λ−μ) and (λσ) one can get an upper bound for the sampling fraction in the form (note that μ ! 0) s ðslÞ ðl À mÞ : If the death rate is known to be much smaller than the birth rate, 0 μ ( λ, the upper bound is expected to be a good estimate for σ. If it is known that the sampling is perfect, σ = 1, one can estimate both the birth and the death rate. However, in contrast to Ref. [11], the method presented here does not require the reconstruction of the tree, but is solely based on statistical properties of pairwise distances between the leaves of the tree. In the general case, one can get an upper limit for the sampling fraction and a lower limit for the birth rate by setting μ/λ = 0. These bounds are expected to be useful for analysis of exponentially growing trees. Such trees can appear in phylogeny when analyzing the evolution of taxa, but also in population genetics, for instance, when considering an exponentially growing sub-population under the influence of a positive selection. To simulate Yule process for the generation of phylogenetic trees we use a Kinetic Monte Carlo algorithm. For a given birth rate λ, death rate μ, and sampling fraction σ, the system is initiated with one "alive" lineage M = 1 at time t = 0. The system is then iteratively propagated to the time t = T. In each iterative step one alive lineage is chosen at random and either either split into two alive lineages (with probability λ/(λ+μ)) or killed (with probability μ/(λ+μ)). In each step the time is incremented by an amount Δt that is exponentially distributed with mean 1/(M (λ+μ)), where M is the number of alive lineages. After the time t = T has been reached, alive lineage are kept in the set of sampled leaves with probability σ.
During the whole simulation the complete tree-especially information about all branching points and branching times-are kept in memory. This way the distribution of pairwise distances or other quantities described in the text can easily be computed. To obtain the mean of such distributions we usually generated at least 10 6 trees and computed the averages.

Second-minimal-distance distribution
Let N 2 (tjT)dt be the average number of leaves on the tree of length T, separated by the time distance t from their second-most closely related leaf. Then, if the first branching occurs at time T 1 and the two resulting subtrees possess M 1 and M 2 leaves, respectively, one gets the distribution of the minimal distance time in a form After average over T 1 and solving the resulting equation one obtains for 0 t 2T. Similarly, one can obtain any third-minimal distance distribution forth-etc. The general formula for the n-minimal-distance distribution is calculated in the following.

n-minimal-distance distribution
Let N n (tjT)dt be the average number of leaves on the tree of length T, separated by the time distance t from their n-most closely related leaf. This notation means that 1-most closely related leaf is the closest one, 2-most closely related leaf is the second-most closest one etc. Then, if the first branching happens at time T 1 and the two resulting subtrees possess M 1 and M 2 leaves, respectively, one gets the distribution of the minimal distance time in a form N n ðtjTÞ ¼ 2N n ðtjT À T 1 Þe ÀmT 1 þ2 nPðnjt=2ÞP > ð0jt=2Þ þ ðn À 1ÞPðn À 1jt=2ÞP > ð1jt=2Þ þ ::: þ Pð1jt=2ÞP > ðn À 1jt=2Þ ½ Âd t À 2 T À T 1 ð Þ ð Þ I 0 t 2T ð Þ e ÀmT 1 ¼ 2N n ðtjT À T 1 Þ þ 2d t À 2 T À T 1 ð Þ ð Þ I 0 t 2T ð Þ X n k¼1 k Pðkjt=2ÞP > ðn À kjt=2Þ " # e ÀmT 1 : Here P > ðkjTÞ ¼ s kþ1 ðm À lÞl k e TðmÀlÞ À 1 ½ k e Tl l À e Tm m ð Þ k e Tl e Tm ðm À l þ slÞ À e Tl ls ½ kþ1 l À e TðmÀlÞ m ½ k ð35Þ is the probability to observe more than k leaves on a tree of age T and P(njT) is given in Equations (3,4). After average over T 1 and solving the resulting equation one obtains

Cherries-distance distribution
A cherry is a pair of adjacent tips on a tree (see Fig. 2). Let N Λ (tjT)dt be the average number of cherry pairs on the tree of length T, separated by the time distance t. Then, if the first branch splits at time T 1 and the two resulting subtrees possess M 1 and M 2 leaves, respectively, one gets the distribution in the form After average over T 1 and solving the resulting equation one obtains The distribution of N 1 (t|T)dt In this Appendix we derive the distribution of N 1 (tjT)dt. Consider an infinitesimal (in practice very small) interval, [t, t+dt], such that N(tjT)dt ( 1. The number of pairs N 1 (tjT)dt in this interval is distributed with the mean N(tjT)dt. The full distribution can be derived using the following arguments. Pairs, separated by the time in the interval [t, t+dt], branched at the time interval [T−t/2−dt/2, T−t/2]. The probability to have a branch in this interval is given by λe (λ−μ)(T−t/2) dt/2. Given that there is a branching point in this interval it can lead to different number of leaves and, therefore, pairs separated by the time in the interval [t, t+dt]. The probability that no observed pairs survive from this branching is given by 1−[1−P(0jt/2)] 2 , where P(njT) is the probability to observe n leaves on a tree of age T and is given in Equations (3, 4). The probability that there are no observed pairs separated by the time in the interval [t, t+dt] is given by Equation (25). The probability that there are n > 0 observed pairs separated by the time in the interval [t, t+dt] is given by Pr N 1 ðtjTÞdt ¼ n ð Þ¼ le ðlÀmÞðTÀt=2Þ dt=2 X n n 1 ;n 2 ¼1 Pðn 1 jt=2ÞPðn 2 jt=2Þd n 1 n 2 ;n ¼ le ðlÀmÞðTÀt=2Þ dt=2 X n 1 jn Pðn 1 jt=2ÞPðn=n 1 jt=2Þ: The last sum runs over all divisors of n, including 1 and n. One can see the comparison of Equations (25) and (39) to numerical results in Fig. 7.