text:
# Folding Pathways of Prion and Doppel

## I Introduction

Neurodegenerative diseases causing transmissible spongiform encephalopathies (TSE) are the subject of intense research. They affect humans and animals and include scrapie in sheep, mad-cow disease in bovines and Creutzfeldt-Jacob disease in humans. Most of them arise sporadically, the rest of the instances are inherited or transmitted by inoculation through the dietary or infected tissues.

The belief is that the normal, and benign, form of the prion protein $`PrP^c`$, which is found in various organs of vertebrates, including the brain, can undergo a post translational process leading to a conformational change of its native state . This new form is designated $`PrP^{sc}`$ and contrarily to $`PrP^c`$ it is insoluble. The NMR investigation of the three-dimensional structure of recombinant $`PrP^c`$ of various species has revealed that the N-terminus is unstructured whereas the C-terminus consists of three $`\alpha `$-helices and two short $`\beta `$-strands.

On the other hand the structure of the malign form, $`PrP^{sc}`$, is mostly unknown. However, spectroscopic investigations, mainly based on circular dichroism, showed that $`PrP^c`$ contains about 40% $`\alpha `$-helices and little $`\beta `$-sheet, whereas$`PrP^{sc}`$ is composed of about 30% $`\alpha `$-helices and 45% $`\beta `$-sheet .

Unfolding and re-folding experiments suggest that the conversion between the $`\alpha `$-rich and $`\beta `$-rich form occurs through a complete unfolding of the protein. Experiments and clinical tests suggest two main processes for the onset and spreading of the prion infection. In the first case, the spreading of the scrapie form in healthy tissues contaminated with $`PrP^{sc}`$ reveals that the presence of $`PrP^{sc}`$ helps further conversion from normal to scrapie form, i.e. $`PrP^{sc}`$ acts as a template for the restructuration of $`PrP^c`$. In the second case, the genetic influence on the propensity of cellular $`PrP^c`$ to rearrange into $`PrP^{sc}`$ is revealed by the discovery of, at least, 20 single point mutations in the PrP gene in Humans which favor the spontaneous onset of the disease .

Within the picture presented in  the role of the mutations in the arousal of the disease could be caused by several factors ranging from the lowering of a free energy barrier in the conversion process to the increase in the oligomerization rate. One of the goals of the present study is to understand which are the sites where mutations are expected to have the major impact in the stabilization/destabilization of the $`PrP^c`$ structure and/or its folding process.

Recently, some experiments on transgenic mice lacking the $`PrP`$ gene ($`PrnP^{0/0}`$) have revealed the onset of neurodegeneration (Purkinje-cell death)  completely different from usual prion disease. This disease was ultimately traced back to a prion homologous gene, $`Prnd`$, (not expressed in normal subjects), which encoded a protein named doppel, $`Dpl`$ which has almost the same topology of $`PrP^c`$, despite the low sequence identity (25%).

A striking difference between $`Dpl`$ and $`PrP^{sc}`$, is that the former affects the central nervous system while retaining its native structure, i.e. without a need to convert to a scrapie-like conformation. Various hypothesis have been formulated  to explain these different behaviours.

The theoretical framework we elaborated and adopted in our study helps to shed light on the putative causes for such remarkably different behavior of doppel and prion. Indeed, the different folding behavior and stability of prion and doppel are related, at least in part, to the small, but important, differences in the topology of the native state which affect in an amplified fashion, the folding routes and ultimately facilitating structural rearrangements of prion.

The model we adopt builds on the importance of the native state topology in steering the folding process, that is in bringing into contact pairs of amino acids that are found in interaction in the native state. In the past few years an increasing amount of experimental  and theoretical  evidence (such as all-atom molecular dynamics (MD) simulations in implicit solvent) have confirmed this view .

A further confirmation of this view is provided by the remarkable accord of the key folding stages predicted by topology-based models and the available experimental results . Indeed, the success of the topological picture itself helps to explain why the folding process is not too sensitive to the detailed chemical composition of most residues in a protein. The dependence of the folding process on the detailed chemistry is much more subtle than that given by the native state topology.

Usually the folding mechanism is affected only when mutations occurs in a small set of key residues. Within our theoretical framework those key residues take part to contacts that are crucial for the folding process .

The establishment of the key contacts leads to a rapid formation of further interactions. Interestingly, the folding bottlenecks can be identified by just knowing the topology of the native conformation . Thus, the topology itself also dictates the impact of chemistry on the folding process.

The purpose of the present work is three-fold:

i) determination of the folding bottlenecks for the cellular prion and the identification of the key amino acids taking part to the corresponding crucial contacts;

ii) making connection between the set of such key residues with those that are known to be associated with harmful $`PrP^{sc}`$ mutations. As argued before, a random mutation on key positions will usually result in a disruption of the folding process. Only fine-tuned mutations can lead to a wild-type-like native state (as for viral enzymes mutating under drug attack ) or in another viable structure (as postulated by Prusiner for the prion );

iii) furthermore, we will focus our attention on how the topological differences between the native states of $`Dpl`$ and $`PrP^c`$ may have impact on the folding process thus aiding or avoiding the formation of misfolded conformers. It is important to stress that our study is based on the native-state structures of $`PrP^c`$ and $`Dpl`$. Hence, although we can confidently identify the crucial folding residues, we cannot confirm explicitly that their mutation results in a different native state form.

## II Materials and Methods

### A The model for prion

Our model builds on a schematic representation of the three-dimensional structure of the protein, where the amino acids are replaced by effective centroids identified with the $`C_\alpha `$ atoms. The bias towards the native state is introduced through a Go-like  energy scoring function that rewards the formation of native contacts between the centroids. The list of native interactions was compiled from the knowledge of the coordinates of atoms of the human prion protein (PDB code 1qlx). A pair of amino acids is considered in interaction if any pair of their heavy atoms, $`i`$ and $`j`$, have a native separation smaller than the distance, $`1.244(R_i+R_j)`$, where the point of inflection of the Van der Waals interaction occurs ($`R_i`$ denotes the Van der Waals radius of atom $`i`$). The values of $`R_i`$ are taken from Ref. . In this way one defines a symmetric matrix, known as contact map, $`\mathrm{\Delta }`$, whose entries, $`\mathrm{\Delta }_{ij}`$ are equal to 1 if the $`i`$-th and $`j`$-th centroids interact, and zero otherwise.

The energy function for our system includes terms that are routinely used in standard molecular dynamics (MD) simulations on biopolymers . It is composed by a ”bonded” and a ”non-bonded” term: $`E=V_B+V_{NB}`$. The former accounts for the constraints such as the peptide bond length and Ramachandran angles bias acting on the aminoacids at a local level. Its explicit form is:

$`V_B`$ $`=`$ $`gV_p+hV_a+kV_d,`$ (1)
$`V_p`$ $`=`$ $`\sigma {\displaystyle \underset{i=1}{\overset{N1}{}}}\left(r_{i,i+1}r_{i,i+1}^{(n)}\right)^2,`$ (2)
$`V_a`$ $`=`$ $`\sigma {\displaystyle \underset{i=1}{\overset{N2}{}}}\left(\theta _{i,i+1,i+2}\theta _{i,i+1,i+2}^{(n)}\right)^2,`$ (3)
$`V_d`$ $`=`$ $`\sigma {\displaystyle \underset{i=1}{\overset{N3}{}}}\left(1\mathrm{cos}\left(\tau _{i,i+1,i+2,i+3}\tau _{i,i+1,i+2,i+3}^{(n)}\right)\right),`$ (4)
where $`r_{i,j}`$ is the distance between residue $`i`$ and $`j`$, $`\theta _{i,j,k}`$ is the angle with the $`j`$-th amino acid as vertex and the $`i`$-th and the $`k`$-th as edges and $`\tau _{i,j,k,l}`$ is the dihedral generated by the $`i`$-th, the $`j`$-th, the $`k`$-th and the $`l`$-th amino acid. The $`n`$ as superscript denotes the native state value. $`\sigma `$ represents a suitable scale factor to fix the temperature scale (set to $`10`$ in our simulations) which is given in dimensionless units.

The minimum of expression (1) is precisely attained in correspondence of the native state. However the cooperative character of the folding process  can be captured only by introducing an explicit bias towards the formation of native interactions . This is accomplished through the second term, $`V_{NB}`$, which weights the interaction of any pair of non-consecutive centroids with a Van der Waals-like potential:

$$V_{NB}=\sigma \underset{i<j3}{\overset{N}{}}V_{ij}=\sigma \underset{i<j3}{\overset{N}{}}\left(5\left(\frac{r_{ij}}{r_{ij}^{(n)}}\right)^{12}\mathrm{\hspace{0.33em}\hspace{0.33em}6}\left(\frac{r_{ij}}{r_{ij}^{(n)}}\right)^{10}\mathrm{\Delta }_{ij}\right).$$
(5)
This choice rewards native interactions and disfavors non-native ones. The form of the energy scoring function ensures that the global minimum is attained in correspondence of the native state, regardless of the precise values for the (positive) couplings parameters $`g`$, $`h`$, $`k`$. The values used in the present study are: $`g=50`$, $`h=5`$ and $`k=0.3`$; we have checked that the results are robust to their variation within a certain range .

The dimensionless unit of time in our MD simulations is $`510^3`$ and conformations were sampled every $`500`$ time steps to allow a sufficient uncorrelation at temperatures below the folding transition. The multiple-histogram technique  has been applied to reweight quantities such as the average internal energy collected in more than 50 equilibrated runs at different temperatures in the range 2-8. The optimal reweighting allows to calculate thermodynamic quantities, such as the average internal energy and the specific heat for an arbitrary temperature, $`T`$.

As visible in Fig. 1A, the specific heat exhibits a single peak which signals the folding transition . In addition to this overall characterization of the equilibrium thermodynamics, it is possible to monitor how the individual interactions that are present in the native state are formed as the folding progresses. Indeed, although all native contacts are energetically favored in the same way, their entropic cost of formation may significantly vary according to their locality, burial/exposure etc. .

The crucial contacts for the folding process are identified according to the method given in Refs. , which singles out the contacts giving the largest contribution to the overall specific heat at the folding transition temperature, $`T_F`$.

$$c_{ij}=\frac{dV_{ij}}{dT}=\left(V_{NB}V_{ij}V_{NB}V_{ij}\right)/kT^2$$
(6)
The physical interpretation of this procedure is that these special contacts are precisely those that, having a significant and sudden formation at $`T_F`$, act as bottlenecks for the folding process. $`c_{i,j}`$ measures the sensitivity of the average energy of the system to perturbation occurring at the contact $`(i,j)`$.

In our studies, a long run carried out at $`T_F`$ of $`210^7`$ time steps, allowed a detailed analysis for all contacts. In addition to the identification of crucial contacts, it is useful to introduce a similar measure for the individual amino acids through the single residue specific heat (SRSH) resulting from the contribution of all contacts to which a given site, $`i`$, takes part to:

$$c_i=\underset{j}{}\mathrm{\Delta }_{ij}c_{ij}$$
(7)
SRSH’s have been used to identify folding and unfolding phases and to pinpoint residues that play a relevant role in the process . The crucial folding steps can thus be identified with two (related) criteria: namely either through the contacts or the sites with the largest specific heat of formation.

### B Folding pathways

Another issue that we have addressed in our study is to establish the presence of alternative pathways taking from unfolded conformations to the native state. Again, this is accomplished by examining in detail both the thermodynamic and the kinetic relevant features extracted from the model. The comparison of dynamic trajectories connecting unfolded and native-like conformations has revealed the presence of several folding events that have to occur for the protein to fold. The main ones that emerged from our analysis, and that are particularly apt for characterizing the folding process, correspond to the formation of contacts involving three sets of contiguous residues. The first set, $`Bh`$, from residue 125 to residue 172 comprises all N-terminal residues ($`\beta `$-sheet included) before the helix 2 (defined according to the indications in the PDB file). The second set, $`h1`$, from residue 173 to 194 comprises residues of helix 2, and finally the third set, $`h2`$, from residue 200 to 228 is composed of residues from helix 3.

The four observed events are:

* Event $`Bh`$: formation of contacts between residues within $`Bh`$,
* Event $`Bh`$-$`h1`$: formation of contacts between $`Bh`$ and $`h1`$,
* Event $`Bh`$-$`h2`$: formation of contacts between $`Bh`$ and $`h2`$,
* Event $`h1`$-$`h2`$: formation of contacts between $`h1`$ and $`h2`$.

For convenience of presentation, in order to describe the four events on a similar footing, we focus on the top 10 contacts (in terms of specific heat contributions $`c_{ij}`$) in each interacting group. It is interesting to note that, among this top contacts, those occurring within an $`\alpha `$ helix almost never occur. This can be explained since local contacts, among others, are easily formed, even at high temperature, and hence do not constitute a bottleneck for overcoming the configurational barrier to the native state.

We computed the fraction of formed contacts, $`Q_i`$ ($`i`$ = $`Bh`$, $`Bh`$-$`h1`$, $`Bh`$-$`h2`$, $`h1`$-$`h2`$), within each of the 4 groups along the dynamical trajectories of $`210^7`$ steps at $`T_F`$. The natural binning size for each of the $`Q_i`$ was $`1/10`$. From this, we obtained a vivid picture of the dynamical effects of the free energy landscape at the folding transition (see Fig. 5 and the Results section).

We also carried out 100 (“quenched”) folding runs at $`T=\left(1ϵ\right)T_F`$, where $`ϵ=1/30`$, slightly below $`T_F`$, starting from completely unfolded conformations (thermalised at high temperature runs). In all cases, the native state was reached in less than $`210^5`$ steps. The average $`Q_i`$’s as a function of the elapsed time using all $`100`$ trajectories was also calculated (Fig. 5 ).

### C Doppel protein

We carried out a parallel analysis for the doppel protein, $`Dpl`$, whose NMR-resolved structure was retrieved from the PDB (PDB code 1i17, model 1). For this protein too, we built the contact map considering all-atom distances, and obtained both the overall specific heat (Fig. 1B ) and the one associated to individual contacts using the same methods described above. A crucial difference in the folding events of the two proteins was revealed by monitoring the four events, analogous to those of the previous section.

For $`Dpl`$ the sets of amino acids structurally homologous of the prion ones are as follows: the set $`Bh`$ corresponds to residues from 1 to 50 before helix 3; $`h1`$ to residues from 51 to 76 (comprising helix 3 and 4) and the set $`h2`$ to residues from 78 to 96 (comprising helix 5). The numeration of the helices is consistent with the secondary structure indications present in the PDB file. Free energy profiles and time dependent averages of the reaction coordinate are shown in Fig. 4 and 5 respectively.

## III Results

### A Determination of key residues

The specific heat curves as a function of the temperature for $`PrP^c`$ and $`Dpl`$ are shown in Figures 1A and and 1B, respectively. The single peak in each of them allows to identify unambiguously the folding transition in these systems.

From the long run of $`210^7`$ time steps performed at $`T_F`$, the SRSH’s for prion have been computed according to eq.(7). We have ranked the residues according to their SRSH (see Table 1). It is natural to compare the top key sites from our analysis with those whose mutation is known to favor the emergence of $`PrP^{sc}`$ (highlighted in Table 2). In the top 10 amino acids of our list we found 4 of the known key mutating sites; whereas in the top 30 amino acids their number increases to 12. The statistical relevance of this match can be obtained from a combinatorics calculation of the probability to have at least the observed number of matches by pure chance. For the top 10 amino acids the probability to get a better results by chance is $`6.8\%`$ whereas for the top 30 amino acids that probability becomes $`0.02\%`$. Thus, the criterion based on the SRSH allows to identify confidently amino acids important for the folding process of the prion protein.

Their structural location and role is made apparent in Fig. 2 where the prion contact map has been colored according to the value of the SRSH. As visible, crucial contacts connect secondary-structure elements thus forming the protein tertiary structure. In particular they connect the $`\beta `$sheet-helix-1 region to the helix-2 and helix-3 ($`Bh`$-$`h1`$ and $`Bh`$-$`h2`$ sets, respectively) and helix-2 - helix-3 ($`h1`$-$`h2`$ set). The parts of the protein contributing more to the specific heat, within our model, are located at the C-terminal strand of the $`\beta `$-sheet, and at the two C-terminal helices, especially helix-3 (see Fig. 3).

### B Folding pathways

The analysis of the long run of $`210^7`$ time steps carried out at $`T_F`$ has led to the construction of an effective free energy landscape. It is based on the reaction coordinates $`Q_i`$ as defined in the Methods section. The free energy as a function of pairs $`(Q_i,Q_j)`$ for all six combinations of $`i`$ and $`j`$ have been calculated as

$$F(Q_i,Q_j)=\mathrm{ln}\left(P(Q_i,Q_j)\right)$$
(8)
where $`P(Q_i,Q_j)`$ is the fraction of times the dynamical trajectory spends in the bin around $`(Q_i,Q_j)`$. The contour plots of $`F`$ are shown in Figures 4A and 4B for $`PrP^c`$ and $`Dpl`$ respectively. The two energy minima around $`Q=0`$ and in the higher $`Q`$ region correspond to the unfolded and (partially) folded conformations respectively. For prion the plots involving set $`h1`$-$`h2`$ show the presence of multiple pathways connecting these two minima. For example, the contour plot in the bidimensional space spanned by $`Q_{h1h2}`$ and $`Q_{Bhh1}`$ presents two possible successions of events: either the $`Bh`$-$`h1`$ set is formed followed by formation of $`h1`$-$`h2`$ set or vice versa (see Fig. 4A). The same observation can be done for the other plots involving set $`h1`$-$`h2`$. Such pathway ambiguity is absent for the case of doppel (Fig. 4B), where trajectories can be grouped essentially in a single set. Indeed, for doppel, on average, the formation of the $`h1`$-$`h2`$ contacts follows the other sets. This result is further confirmed by the analysis of the time-dependent average of the reaction coordinates $`Q_i`$’s as obtained from the $`100`$ non-equilibrium folding runs. Figures 5A and 5B show the results for $`PrP^c`$ and $`Dpl`$ respectively. For $`Dpl`$ the order in the formation of the contacts is the same detected by the the equilibrium run: set $`Bh`$ is formed before set $`Bh`$-$`h1`$ that is followed by $`Bh`$-$`h2`$ eventually followed by $`h1`$-$`h2`$ interhelical contacts. On the other hand, such definite succession of events is not observed for prion.

## IV Discussion

The goal of our numerical study was to characterize the folding process of the prion protein elucidating the structural mechanisms responsible for the propensity of $`PrP^c`$ to mutate into the harmful scrapie form, $`PrP^{sc}`$.

At the heart of our analysis is a topology-based energy function that, by rewarding the formation of native-like interactions among the residues, allows to mimic the progressive build-up of native structure in a folding process. Remarkably, the progress towards the native state does not advance smoothly, but occurs through the overcoming of configurational barriers, which are responsible for the appearance of a marked peak in the specific heat curve. Crucial contacts (and sites) for the folding process have thus been identified as those giving the largest contributions to the specific heat peak (Fig. 2).

The top sites isolated by our procedure contain a highly significant fraction of residues that are known to enhance the misfolding propensity of $`PrP^c`$. This indicates that there is little freedom in choosing the chemical identity of the amino acids at the sites which have a special structural role during the folding process.

The analysis presented here, based on topological arguments, provides a physically-appealing interpretation for the connection between certain site mutations and the arousal of the prion disease. Conversely, the key mutating sites not captured by our analysis are presumably crucial because they affect some aspects of the process leading to the scrapie form that cannot be accounted for by the present method. For example, the alternative mechanism leading to misconformation is the template binding of $`PrP^c`$ to $`PrP^{sc}`$ involving the interconversion of the former into the latter one. Consistently with this interpretation, amino acids that do not present high SRSH and that are known to be involved in the development of the disease, specifically MET 129, GLU200 and GLU 219, (i.e. mostly charged residues with a large exposed surface area).

Despite the substantial topological similarities of $`Dpl`$ and $`PrP^c`$, the detailed analysis of the folding trajectories of both proteins have revealed significant differences of the folding pathways. Strikingly, for $`PrP^c`$ the folding pathways can be classified in two main distinct routes, whereas for $`Dpl`$ only one of these is essentially present. For the latter, the series of folding steps is more distinctly marked than for the former. Such noticeably different behavior can be ascribed to the different set of contacts stabilizing the assembly of the $`\alpha `$-helices and $`\beta `$-sheets in the two proteins. It is interesting to note that the contact order  of $`PrP^c`$ is 0.2056 and is lower than that of $`Dpl`$ which is 0.2298, even though the former has somewhat higher number of contacts which is, in turn, reflected by the shift of the folding transition peak towards higher temperatures (see Fig. 1). This difference could explain the different average folding time measured in our models for the two proteins. Indeed, according to Fig. 5, the sets of contacts in prion are all formed within the first 25000 time steps while contacts in doppel contacts can take 4-5 times longer to form. This difference is of the same order of magnitude of that predicted by Plaxco  on the basis of the contact order of the two proteins. Although the present picture would certainly be influenced by the introduction of amino acid specific interactions, it is appealing to connect the different folding routes to the appearance of misfolded conformers. According to this view, the latter would result from following the folding route of $`PrP^c`$ where the inter C-terminal helices contacts (set $`h1`$-$`h2`$) are formed while the contacts with the $`\beta `$-sheet (set $`Bh`$-$`h1`$ and $`Bh`$-$`h2`$) are not yet formed (Figures 4A and 5A). Indeed, along this route, which is alternative to the one where the two events are interchanged as it occurs in the single route of doppel (Fig. 4B and 5B), the N-terminal part of the protein is free to rearrange its structure before reaching the native-state. This scenario is consistent with several recent models (see e.g. Ref. ) that identify the N-terminal part of the protein as the one that undergoes the conformational change that leads the scrapie form of prion.

Figure captions

1. Figure 1. Temperature dependent specific heat for the prion (A) and doppel (B) models. The folding transition temperature is identified at the point where the specific heat peak occurs.
2. Figure 2. Contact map of $`PrP^c`$ (PDB code 1qlx) where contacts are highlighted according to their specific heat $`c_{ij}`$ of equation (6). Red spots correspond to contacts with high $`c_{ij}`$, while blue ones correspond low SRSH. The relevant sets of contacts (as defined in the text) are framed within the rectangles; (a) set $`Bh`$, (b) set $`Bh`$-$`h1`$, (c) set $`Bh`$-$`h2`$, (d) set $`h1`$-$`h2`$. Notice that red spots are prevalent in all sets but $`Bh`$.
3. Figure 3. The three-dimensional structure of prion protein where amino acids have been colored according to the value of their SRSH. Red spots are located mainly on the C-terminal strand of the $`\beta `$sheet, on helix-2 and helix-3.
4. Figure 4. Effective free energy landscape for prion (A) and doppel (B) model. Upper and right squares represent the contact maps where the representative set of contacts defining the reaction coordinates are highlighted (see Method section). Contour plots are shown for each of the six pairs of reaction coordinates $`(Q_i,Q_j)`$. White (black) spots represent low (high) values of the free energy. The minimum close to (0,0) corresponds to the unfolded state whereas the minimum closer to (1,1) corresponds to the folded state. In most of the contour plots for prion (A) more than one class of favorable paths from the unfolded to the folded conformation is possible whereas for doppel (B) only one class of paths is the most favorable. Following these pathways allows to assign a preferential order to the events. Thus, on average, the sequence of contact formation for doppel is $`Bh`$, $`Bh`$-$`h1`$, $`Bh`$-$`h2`$ and $`h1`$-$`h2`$ (see also Fig. 5). For prion, two class of equivalent pathways are present: in one case the $`h1`$-$`h2`$ contacts are formed before the other sets of contacts; in the other case $`h1`$-$`h2`$ contacts follow the formation of all the other contacts as for $`Dpl`$.
5. Figure 5. Average time dependence of reactions coordinates $`Q_i`$’s for prion (A) and doppel (B) computed on 100 folding runs starting from completely unfolded conformations. (A) In case of prion the curves corresponding to the different sets (red for $`Bh`$, green for $`Bh`$-$`h1`$, blue for $`Bh`$-$`h2`$ and magenta for $`h1`$-$`h2`$) present intersections. Indeed the average over the different folding pathways results in a non well-defined sequence of folding events (B) The folding events in the case of doppel protein have a well defined sequence: $`Bh`$ set (red) is the first to form; then $`Bh`$-$`h1`$ (green); then $`Bh`$-$`h2`$ (blue); finally $`h1`$-$`h2`$ (magenta).


format: text

text:
# 1 Area as the affine parameter

## 1 Area as the affine parameter

Suppose $`N`$ is a 2dim spacetime, i.e. a surface with a metric tensor of signature $`(1,1)`$. It is fairly easy to find the isotropic geodesics of $`N`$ – we just integrate the isotropic directions. It is perhaps more interesting to find affine parameters for these geodesics: one can use the area of a strip between two infinitesimally close integral curves:

To prove this fact we use the following useful characterization of geodesics in a pseudo-Riemannian manifold $`M`$ (for the purpose of teaching of general relativity, it appears as a convenient definition): If we choose coordinates near a curve $`\gamma `$ so that the metric tensor is a constant plus $`O(r^2)`$, where $`r`$ is (say) the Euclidean distance to $`\gamma `$, then $`\gamma `$ is a geodesic iff it is a straight line. Now our claim that area can be used as the affine parameter is clear, since it is surely so for constant metric tensor.

We shall use our result in the next section to prove Birkhoff theorem, but now, as a digression, we mention some other elementary applications. We can use it conveniently to check isotropic geodesic completeness of 2dim spacetimes. Consider, for example, the Eddington-Finkelstein metrics:

$$ds^2=(1\frac{2m}{r})du^2+2dudr.$$

We prove that the geodesic $`\gamma `$ on the following picture is incomplete to the left; we do not use the exact form of the metrics, only the fact that it is invariant w.r.t. horizontal translations and that the geodesics converge on the left (on the picture, $`u`$ is the horizontal coordinate):

Indeed, the green area between $`\gamma `$ and $`\overline{\gamma }`$ is equal to the red area: to see it, just take the green triangle and translate it to the right. The red area is finite, q.e.d.

As another example, consider a $`2n`$-gon with isotropic sides. Widen each of its sides to an infinitesimally narrow strip between two isotropic curves, and compute the expression

$$\mu =\frac{A_1A_3\mathrm{}A_{2n1}}{A_2A_4\mathrm{}A_{2n}},$$

where $`A_i`$ is the area at the $`i`$’th corner, as on the picture:

$`\mu `$ is clearly independent of the choice of the widening of the sides: if we widen the $`i`$’th side in a different way, $`A_i`$ and $`A_{i+1}`$ get multiplied by the same number, so that $`\mu `$ doesn’t change. We leave it as an excercise to the reader to prove that $`\mu `$ is actually the result of the parallel transport along the polygon, so that it is equal to the exponential of $`\pm `$ the integral of the curvature inside the polygon.

## 2 Birkhoff theorem

Suppose $`M`$ is a 4dim spacetime on which $`SO(3)`$ acts by isometries, so that all the orbits are spheres. Birkhoff theorem states that under some assumption on the Ricci tensor, there is a 1-parameter group of isometries of $`M`$, commuting with $`SO(3)`$.

As one easily sees, $`M`$ is of the form $`M=N\times S^2`$, where $`N`$ is a surface, and

$$ds_M^2=ds_N^2+r^2ds_{S^2}^2,$$

where $`r>0`$ is a function on $`N`$ and $`ds_{M,N,S^2}^2`$ are the metrics on $`M`$, $`N`$ and the unit sphere. Let $`R`$ and Ric be the Riemann and Ricci tensor on $`M`$, and $`\omega _N`$ the area form on $`N`$. Finally, let $`X_r`$ be the Hamiltonian vector field on $`N`$ generated by $`r`$ with symplectic form $`\omega _N`$, i.e.

$$dr=\omega _N(,X_r).$$

Birkhoff theorem: If $`\text{Ric}(v,v)=0`$ for any isotropic $`v`$ tangent to $`N`$, then $`X_r`$ is a Killing vector field on $`M`$.

The proof is split into two lemmas:

Lemma 1: If $`N`$ is a 2dim spacetime then a vector field $`w`$ on $`N`$ is conformal iff for any isotropic geodesic $`\gamma `$, $`\omega _N(\dot{\gamma },w)`$ is constant along $`\gamma `$.

Lemma 2: Under the same assumptions as in Birkhoff theorem, if $`\gamma `$ is any isotropic geodesic in $`N`$ and $`s`$ is its affine parameter, then $`dr/ds=\text{const.}`$ along $`\gamma `$.

Proof of the theorem: We have to prove that the flow of $`X_r`$ preserves $`r`$ and $`ds_N^2`$. First we show that $`X_r`$ is conformal on $`N`$. It follows from the lemmas, since

$$\text{const.}=dr/ds=\dot{\gamma },dr=\omega _N(\dot{\gamma },X_r).$$

It remains to check that $`X_r`$ preserves $`\omega _N`$ and $`r`$, but it certainly does, since it is a Hamiltonian vector field q.e.d.

Proof of Lemma 1: This follows immediately from our result in §1. Indeed, $`w`$ is conformal iff its flow transports isotropic curves to isotropic curves. As we noticed, the area between $`\gamma `$ and an infinitesimally close curve $`\overline{\gamma }`$ is a linear function of $`s`$ iff $`\overline{\gamma }`$ is isotropic q.e.d.

Proof of Lemma 2: If $`\gamma `$ is an isotropic geodesic in $`N`$ and $`PS^2`$ then $`\gamma \times \{P\}`$ is a geodesic in $`M`$. It follows from the $`O(2)`$-symmetry of isometries of $`S^2`$ preserving $`P`$: a geodesic is uniquely determined by its velocity at a point, so that if the velocity is $`O(2)`$ invariant, the geodesic must be pointwise $`O(2)`$-invariant, i.e. it lies in $`N\times \{P\}`$.

If $`Q`$ is a different point in $`S^2`$ then $`\gamma \times \{Q\}`$ is also a geodesic; therefore, if we take $`Q`$ infinitesimally close to $`P`$ we get that any constant vector $`aT_PS^2`$ satisfies the Jacobi (geodesic deviation) equation

$$\ddot{a}+R(a,\dot{\gamma })\dot{\gamma }=0.$$

On the other hand, the parallel transport along $`\gamma \times \{P\}`$ must be $`O(2)`$-equivariant, i.e. it acts as a multiple of the identity on $`T_PS^2`$. It also preserves lengths, so that $`a/r`$ is parallel, and

$$\ddot{a}=(ra/r)\ddot{}=\ddot{r}a/r.$$

To prove that $`\ddot{r}=0`$ it remains to show that if $`v`$ is an isotropic vector tangent to $`N`$ and $`a`$ is a vector tangent to $`S^2`$, then $`R(a,v)v=0`$. By definition, $`\text{Ric}(v,v)`$ is the trace of the linear map $`A`$ defined by $`A(w)=R(w,v)v`$. From $`O(2)`$ symmetry we see that $`A(a)=\lambda a`$ for some number $`\lambda `$ (independent of $`a`$); on the other hand, $`A`$ restricted to $`TN`$ is nilpotent, since $`A(v)=0`$ and $`A(w)`$ is orthogonal to $`v`$ for any $`w`$, hence it is a multiple of $`v`$. Therefore $`0=\text{Tr}A=2\lambda `$, i.e $`A(a)=0`$, q.e.d.


text:
# Detection of non-Gaussian Signatures in the VIRMOS-DESCART Lensing SurveyBased on observations obtained at the Canada-France-Hawaii Telescope (CFHT) which is operated by the National Research Council of Canada (NRCC), the Institut des Sciences de l’Univers (INSU) of the Centre National de la Recherche Scientifique (CNRS) and the University of Hawaii (UH).

## 1 Introduction

The large-scale structure of the Universe are expected to form from the gravitational growth of initial density perturbations obeying Gaussian statistics. As the Universe expands and the perturbations grow, non-Gaussian features are expected to emerge in the density field due to gravitational dynamics. These features can be characterized with perturbation theory calculations, which allow to compute for instance the skewness, third moment of the local density probability distribution function (Peebles 1980, Fry 1984, Bernardeau 1992). The reduced skewness of the density field has been found to be quite insensitive to the variance and the cosmological parameters, $`\mathrm{\Omega }_m`$ (Juszkiewicz et al. 1992, Bernardeau 1994). In contrast, weak lensing surveys are sensitive to $`\mathrm{\Omega }_m`$ since they trace the integrated mass along the line-of-sight which is roughly proportional to the density parameter of the Universe. Weak lensing by large scale structures has been measured by several teams as a coherent distortion field of distant galaxies over large angular distances (Bacon et al. 2000 and 2002; Hämmerle et al.; 2002; Hoekstra et al. 2002; Kaiser et al. 2000; Maoli et al. 2001; Réfrégier et al 2002; Rhodes et al. 2001; Van Waerbeke et al. 2000, 2001 and 2002; Wittman et al. 2000). The projected mass density reconstructed from the distortion field ($`i.e.`$ the convergence field) can be used for non-Gaussian signatures searches, as shown in Bernardeau, van Waerbeke & Mellier (1997, hereafter BvWM97). This work and further studies (Jain & Seljak 1997; van Waerbeke, Bernardeau & Mellier 1999) have shown that the non-Gaussian properties of the convergence field can be used as a probe of the cosmological density parameter, with a weak dependence on the cosmological constant $`\mathrm{\Omega }_\mathrm{\Lambda }`$, provided that the redshift distribution of the sources is known.

However, a straight application of these theoretical considerations to real data sets turned out to be arduous. The convergence field has to be recovered from a mass reconstruction process which uses a continuous shear field obtained from a smoothed map of the discrete galaxy ellipticities. Unfortunately, survey topologies are generally complex and are alterated by masked areas due to light scattering, bright stars, comet-like reflections, asteroid/airplane tracks, very bright galaxies, etc…. Since the mask sizes cover a range of scales from a few arc-seconds to two degrees and are strongly anisotropic (for instance bright stars preferentially saturate along CCD columns), results of mass reconstruction in such data sets are not yet reliable. An alternative approach is the aperture mass applied to cosmic shear (Schneider et al. 1998), which allows the measurement of the skewness from the distortion field directly, bypassing the mass reconstruction process. So far, our attempts for measuring the skewness of the aperture mass lead to very noisy and un-significant results. Recently, Bernardeau, van Waerbeke & Mellier (2002, hereafter BvWM02) have proposed a new method using some specific patterns in the shear three-point function. This method has also the advantage to bypass the mass reconstruction process. Despite the complicated shape of the three-point correlation pattern, BvWM02 uncovered a specific geometrical property and demonstrated it can be used for the measurement of non-Gaussian features. Their detection strategy based on this method has been found to be robust, usable in patchy catalogs, and quite insensitive to the topology of the survey. In the following we apply this method to the VIRMOS-DESCART weak lensing survey done at the Canada-France-Hawaii Telescope.

## 2 Optimized 3-point correlation function applied to VIRMOS-DESCART data

### 2.1 Method

Let us consider a triplet of galaxies at locations $`x_1`$, $`x_2`$ and $`x^{}`$, and their shear estimates $`\gamma (x_1)`$, $`\gamma (x_2)`$ and $`\gamma (x^{})`$. BvWM02 introduced a 2D angular vector field (representing 2 components of the general shear 3-point correlation functions),

$$\xi _3(x^{})=(\gamma (x_1)\gamma (x_2))\gamma (x^{})$$
(1)
viewed as a function of $`x^{}`$ for a fixed $`|x_2x_1|`$ distance. At fixed separation $`x_{12}`$, the components of $`\xi _3`$ are expected to scale like the square of the shear 2-point correlation function, $`\xi _2(x_{12})=\gamma (x_1)\gamma (x_2)`$. Its sensitivity to the cosmological density parameter ($`\mathrm{\Omega }`$) is expected to be similar to the convergence skewness ($`\mathrm{\Omega }^{0.8}`$), with an additional dependence on the slope of the mass power spectrum. The way $`\xi _3(x^{})`$ varies as a function of $`x^{}`$ is in general complicated although asymptotic properties can be obtained analytically. However $`\xi _3(x^{})`$ has been found to be rather uniform, and perpendicular to $`x_{12}`$, over an elliptic area that covers the segment joining $`x_1`$ to $`x_2`$. This central pattern turns out to be robust against different cosmologies and smoothing scales, with an amplitude which can be related to the cosmological parameters. The results obtained by BvWM02 in synthetic catalogs suggested that present-day cosmic shear surveys were already large for a secure detection.

We consider the geometrical average,

$`\overline{\xi _3}(|x_1x_2|)={\displaystyle _{\mathrm{Ell}.}}{\displaystyle \frac{\mathrm{d}^2x^{}}{V_{\mathrm{Ell}.}}}\xi _3(x^{})`$ (2)
which corresponds to the average three-point function inside an elliptic area, $`V_{\mathrm{Ell}}`$, encompassing the points $`x_1`$ and $`x_2`$ (they are actually the foci of the ellipse). In this domain $`\xi _3(x^{})`$ is expected to vary weakly and its average value $`\overline{\xi _3}`$ is a vector quantity orthogonal to the $`x_{12}`$ direction111mathematically it means that the only non-zero component of $`\overline{\xi _3}`$ is along $`x_1x_2`$, but has a negative value. that depends only on the $`x_1`$-$`x_2`$ pair separation. An optimum selection of pair points, where close pairs and highly elliptical galaxies are rejected, turns out to provide an optimal information on non-Gaussian features in simulated catalogues. BvWM02 checked that different sources of noise produced by the intrinsic ellipticity distribution of galaxies, by realistic galaxy shape measurements and PSF corrections and by masking effects do not spoil the result. Final tests and validations were carried out with simulated catalogues containing as many galaxies as real data, with ellipticity distribution, PSF anisotropy and masking templates similar to the VIRMOS-DESCART sample. In all configurations the global signal to noise remains higher than 5 for scales between 30 arc-seconds to 5 arc-minutes.

### 2.2 The VIRMOS-DESCART sample

The VIRMOS-DESCART sample used in this work is part of the DESCART cosmic shear programme222http://terapix.iap.fr/Descart which uses the VIRMOS photometric/imaging survey333http://www.astrsp-mrs.fr for wide field cosmic shear. It covers 11.5 deg2 of CFH12K images spread over four uncorrelated fields. All data were obtained in I-band up to a limiting magnitude of $`I_{AB}=24.5`$ (within 5 arc-second aperture, 5-$`\sigma `$), which is consistent with a mean source redshift of $`z0.9`$ of the VIRMOS-DESCART sample (van Waerbeke et al. 2002). The data contain all observations used by van Waerbeke et al. (2001) plus new fields obtained in September and November 2000. All the images were processed as described in van Waerbeke et al (2001) at the TERAPIX data center444http://terapix.iap.fr. From an initial detection we build up a sample containing 1.6 millions objects. After masking and all galaxy selection processes, it reduces to 580,000 galaxies, covering an effective area of 8.5 deg2. Close-pairs with angular separation smaller than 10 arc-seconds can produce systematics and are rejected (see van Waerbeke et al. 2000).
Following standard notation, the shear is defined as the mean source ellipticities

$$\gamma _1=e\mathrm{cos}(2\theta ),\gamma _2=e\mathrm{sin}(2\theta )$$
(3)
where $`\theta `$ is the angle between the major axis of the source galaxy and the $`x`$ axis, and $`e=(ab)/(a+b)`$ is determined by the major axis length $`a`$ and minor axis length $`b`$. It is computed according to the rules and weighting schemes given in Pen, van Waerbeke & Mellier (2002) and the PSF anisotropy corrections used in van Waerbeke et al. (2002). In particular, we consider in this Letter the decomposition into $`E`$ and $`B`$ modes described in Pen et al. (2002) which is used for residual systematics checks for the 2-point statistics.

The statistical estimators for the binned two and three point functions are given by

$`\xi _2(d_{ij})`$ $`=`$ $`{\displaystyle \frac{\underset{ij}{}w_iw_j(e_ie_j)}{_{ij}w_iw_j}}`$ (4)
$`\overline{\xi }_3^t(d_{ij})`$ $`=`$ $`{\displaystyle \frac{\underset{ijk}{}w_iw_jw_k(e_ie_j)e_k^{(ij)}}{_{ijk}w_iw_jw_k}}`$ (5)
where $`e_k^{(ij)}`$ is opposite to the component of the ellipticity of galaxy $`k`$ along the ($`x_jx_i`$)-direction. The summations are made for pairs or triplets such that $`d_{ij}=|x_ix_j|`$ is in the chosen bin. $`x_k`$ lies within the ellipse defined by $`|x_kx_i|+|x_kx_j|>1.1|x_ix_j|`$ and $`w_i`$ are weights associated to each galaxy according to the scheme discussed in BvWM02.

### 2.3 The VIRMOS-DESCART 2 and 3-point correlation functions

Figure 1 shows the estimated 2-point (top) and reduced 3-point (bottom) correlation functions ($`\xi _3`$ in units of $`\xi _2^2`$). Measurements have been made in regularly spaced bins of width 400 pixels (e.g. 1.3 arcmin). For comparison, the thick dotted lines shows the corresponding quantity for an open CDM model, the thick dashed line for a $`\tau `$CDM and the dot-dasled line for a $`\mathrm{\Lambda }`$CDM model. The source redshift is one, very close to our mean source redshift of $`0.9`$. In this plot the error bars are the measurement errors and do not include the cosmic variance (see BvWM02). The dot-dot-dashed line represents the signal corrected from the residual systematics (the $`B`$-mode contribution to the two-point correlation function subtracted from the $`B`$ mode contribution as discussed in Pen et al. 2002), while the solid line shows the total ($`E+B`$) measured signal. The closeness of the two curves reveals the small amount of residual systematics still present in the 2-point correlation function.
A quality assessment of the 3-point function measurements can be done by studying the effect of PSF correction. It is inferred from the star shapes (they ought to be round) following the standard procedure described initially in Kaiser, Squires & Broadhurst (1995). It can be thought as a star shear field that has to be subtracted off the measured galaxy shapes. One can then measure the 2- and 3-point correlation functions of the star field and compare them to the PSF corrected data set. This is presented on Figure 2. The solid line in the top and bottom plots respectively show the star two and three point functions, which are compared to the same quantities measured on the corrected galaxies (dashed lines and dashed-dot line for the corrected $`E`$-mode 2-point function). The star two-point correlation function is significantly larger than the galaxy signal, but this is known not to be a problem for the two-point function: Erben et al. (2001) and Bacon et al. (2001) have shown that the PSF correction can account for star anisotropies as large as $`10\%`$ (the stars r.m.s. anisotropy is $`6\%`$ in our data) to a precision better than than one percent. This is also demonstrated in BvWM02 using cosmic shear simulated data. The small amplitude of the $`B`$-mode is another proof of the robustness of the correction.

BvWM02 have also shown that the correction scheme works for the shear three-point function. We see that the situation for the three-point function (Figure 2, bottom panel) is far more interesting than the two-point function: the correction is smaller (that is in principle more robust), and the angular dependence of the star three point function is totally different from the galaxy three-point function. If our three-point function signal were dominated by systematics it would likely not fit the expected signal as we can see from Figure 1, but would be more similar to Figure 2, bottom panel (solid line). Finally we have checked that cuts in magnitudes do not significantly change the results. Therefore, they are not produced by one single class of objects of the sample. We have thus some evidences that what is observed is genuinely cosmological.

Another issue, not taken into account in the error bars, is the so-called cosmic variance, that is the amount by which such a signal can vary in surveys of finite size. The complete theory of cosmic variance for such a survey is yet to be done and is likely complex (given the topology). Figure 3 shows the cosmic variance obtained from a set of 7 ray-tracing realizations of the open CDM model (Jain et al. 2000). Note that the error-bars in the different bins *are correlated*. A rough examination of the situation shows that while the cosmic variance dominates the error budget for the two-point function, observational errors and cosmic variance are of similar amplitude for the three-point function. From the examination of Fig. 1 $`\mathrm{\Lambda }`$CDM is clearly in better agreement with the data than the other two models, while $`\tau `$CDM is marginally acceptable.

## 3 Discussion

The result shown on Figure 1 is the first detection of non-Gaussian features in a cosmic shear survey. The signal is detected with a 2.4-$`\sigma `$ confidence level on 4 independent bins which gives a 4.9-$`\sigma `$ global confidence level. Such a result opens the route to break the degeneracy between $`\mathrm{\Omega }_m`$ and $`\sigma _8`$ in a way which is independent on assumptions beyond the solely hypothesis that large-scale structures grows from gravitational instability of an initial Gaussian field. We note that the amplitude of the reduced 3-point correlation function exhibits an angular dependence which is in agreement with theoretical expectations. It supports the interpretation of these results as genuine effects of the gravitational dynamics.

The signal is however still too noisy to provide reliable information on cosmological parameters. Moreover, several obstacles have yet to be overcome. It is first important to understand to which level the measurements are contaminated by systematics. This can be done through consistency checks yet to be invented (a statistic which cancels the signal in a non-trivial way, like the $`45`$ degrees rotation test for the aperture mass will have to be found for the three-point function). We have already checked with the anisotropy of the stars that our signal is unlikely to be dominated by systematics since it would otherwise exhibits a totally different angular dependence. In the next stages, the scientific interpretation of the three-point function measurements will require a significant improvement over several issues:

-The knowledge of the redshift distribution of the sources is crucial for the third order statistics (see BvWM97);

-The source clustering effect might bias the measurement is a significant way (Bernardeau 1998, Hamana et al. 2002) if the width of the source distribution is too large;

-Intrinsic alignment of galaxies have a completely unknown effect on the non-Gaussian properties of the shear field.

Resolving these issues will require progress from both the theoretical/simulation side and from the observations, which are already on their way.

###### Acknowledgements.

We thank the VIRMOS and Terapix teams who got and processed the VIRMOS-DESCART data, B. Jain and M. White for the use of their ray-tracing simulations as well as T. Hamana for useful comments. This work was supported by the TMR Network “Gravitational Lensing: New Constraints on Cosmology and the Distribution of Dark Matter” of the EC under contract No. ERBFMRX-CT97-0172. FB thanks IAP for hospitality.


text:
# Random Distance Distribution for Spherical Objects: General Theory and Applications to 𝑛-Dimensional Physics

## I Introduction

In two recent papers D. Schleef, M. Parry, S. J. Tu, B. Woodahl, and E. Fischbach (1999); Parry and Fischbach (2000), geometric probability techniques were developed to calculate the functions $`P_3\left(s\right)`$ which describe the probability density of finding a random distance $`s`$ separating two random points distributed in a uniform sphere and in a uniform ellipsoid. As discussed in Refs. D. Schleef, M. Parry, S. J. Tu, B. Woodahl, and E. Fischbach (1999); Parry and Fischbach (2000); Deltheil (1919); Hammersley (1950); Lord (1954); Kendall and Moran (1963); Santaló (1976); Overhauser ; Fischbach (1996); E. Fischbach, D. E. Krause, C. Talmadge, and D. Tadić (1995); Tu (2001); Fedotov (1998); Jansons (1991); Mecke (1998); Miles (1980), these results are of interest as tools in mathematical physics, and have numerous applications in other fields as well. Specifically, it was demonstrated in Refs. D. Schleef, M. Parry, S. J. Tu, B. Woodahl, and E. Fischbach (1999); Parry and Fischbach (2000); Fischbach (1996); E. Fischbach, D. E. Krause, C. Talmadge, and D. Tadić (1995); Tu (2001) that knowing the random distance distribution in a spherical object greatly facilitates the calculation of self-energies for spherical matter distributions arising from electromagnetic, gravitational, or weak interactions.

As an example we calculate the total electrostatic energy $`W_3`$ of a collection of $`Z`$ charges uniformly distributed within the same spherical volume of radius $`R`$. For illustrative purposes, we assume that $`Z`$ is a large number. For each pair of charges the potential energy due to the Coulomb interaction in Gaussian units is

$$V_3=\frac{e^2}{r_{12}}=\frac{e^2}{\left|\stackrel{}{r}_2\stackrel{}{r}_1\right|},$$
(1)
where $`e`$ is the elementary charge, $`\stackrel{}{r}_1`$ ($`\stackrel{}{r}_2`$) is the coordinate of the first (second) charge. The total Coulomb energy $`W_3`$ can then be expressed as

$`W_3`$ $`=`$ $`{\displaystyle \frac{1}{2}}\times \rho ^2{\displaystyle r_1^2𝑑r_1\mathrm{sin}\theta _1d\theta _1𝑑\varphi _1r_2^2𝑑r_2\mathrm{sin}\theta _2d\theta _2\frac{1}{r_{12}}𝑑\varphi _2}`$ (2)
$`=`$ $`{\displaystyle \frac{Z^2}{2}}\times {\displaystyle \frac{6}{5}}{\displaystyle \frac{e^2}{R}},`$
where $`\rho =3Ze/4\pi R^3`$ Jackson (1999). We note that Eq. (2) requires evaluating a six-dimensional integral, and using the additional theorem for spherical harmonics Arfken and Weber (1995),

$$P_l\left(\mathrm{cos}\gamma \right)=\frac{4\pi }{2l+1}\underset{m=l}{\overset{l}{}}\left(1\right)^mY_l^m(\theta _1,\varphi _1)Y_l^m(\theta _2,\varphi _2).$$
(3)
Alternatively we can use the probability density function $`P_3\left(s\right)`$ giving the random distance distribution for a sphere with a uniform density to calculate $`W_3`$. For a collection of $`Z`$ charges there are $`Z\left(Z1\right)/2`$ such pairs, and hence the total Coulomb energy $`W_3`$ is:

$$W_3=\frac{Z\left(Z1\right)}{2}\times _0^{2R}P_3\left(s\right)\frac{e^2}{s}𝑑s=\frac{Z\left(Z1\right)}{2}\times \frac{6}{5}\frac{e^2}{R}\frac{Z^2}{2}\times \frac{6}{5}\frac{e^2}{R},$$
(4)
where Deltheil (1919)

$$P_3\left(s\right)=3\frac{s^2}{R^3}\frac{9}{4}\frac{s^3}{R^4}+\frac{3}{16}\frac{s^5}{R^6}.$$
(5)
We note that by using geometric probability techniques we can simplify the expression for $`W_3`$ from a six-dimensional non-trivial integral (Eq. (2)) to a one-dimensional elementary integral (Eq. (4)). Generalizing to $`n`$-dimensions, a calculation of the electrostatic energy $`W_n`$ for a collection of $`Z`$ charges uniformly distributed within the same $`n`$-dimensional spherical volume of radius $`R`$, can be greatly simplified by utilizing the $`n`$-dimensional random distance distribution. This reduces the complexity of calculating $`W_n`$ from a $`2n`$-dimensional intractable integral involving $`n`$-dimensional spherical harmonics, to a simple $`1`$-dimensional integral.

The probability density function $`P_n(s)`$ for the distribution of the random distance $`s`$ between two random points in a uniform spherical $`n`$-ball is well known. Hence the object of the present paper is to generalize the results of Refs. Kendall and Moran (1963); Santaló (1976); Overhauser  to the case of an arbitrarily non-uniform density distribution by using a new method which we present below. Notice that the sample space $`𝐁_n`$ for the random points is a spherical $`n`$-ball of radius $`R`$ defined as

$$𝐁_n=\{(x_1,x_2,\mathrm{},x_n)𝐑^n:x_1^2+x_2^2+\mathrm{}+x_n^2R^2\},$$
(6)
where

$$𝐑^n=\{𝒙:𝒙=(x_1,x_2,\mathrm{},x_n)\}$$
(7)
represents $`n`$-dimensional Euclidean space.

To illustrate our formalism, we begin by deriving the probability density function (PDF) for a spherical $`n`$-ball with a uniform density distribution, and compare our results to those obtained earlier by other means Kendall and Moran (1963); Santaló (1976); Overhauser . We then extend this technique to a spherical $`n`$-ball with an arbitrary density distribution, and this leads to a general-purpose master formula for $`P_n\left(s\right)`$ given in Eq. (39). The outline of this paper is as follows. In Sec. II we present our formalism and illustrate it by rederiving the well-known results for a circle and for a sphere of uniform density. In Sec. III we extend this formalism to the case of non-uniform but spherically symmetric density. In Sec. IV we develop the formalism for the most general case of an arbitrary density distribution. In Sec. V we present some applications of our results. These include the $`m`$th moment $`s^m`$ for a spherical space with a uniform and Gaussian density distribution, the evaluation of the self-energy arising from $`\nu \overline{\nu }`$-exchange interactions, the probability density functions for a sphere with multiple shells that arises in neutron star model Pines (1980); Shapiro and Teukolsky (1983), and finally a new proposed computational scheme for testing random number generators in $`n`$-dimensions.

## II Uniform Density Distributions

In this section we illustrate our formalism by deriving the PDF for a circle of radius $`R`$ having a spatially uniform density characterized by a density function $`\rho `$, where $`\rho `$ is an arbitrary constant. For two points randomly sampled inside the circle located at $`\stackrel{}{r}_1`$ and $`\stackrel{}{r}_2`$ measured from the center, define a random vector $`\stackrel{}{s}=\stackrel{}{r}_2\stackrel{}{r}_1`$ and a random distance $`s=\left|\stackrel{}{s}\right|`$, where $`0s2R`$. To simplify the discussion, we translate the center of the circle to the origin so that the equation for the circle is $`x^2+y^2=R^2`$. It is sufficient to initially consider those vectors $`\stackrel{}{s}`$ which are aligned in the positive $`\widehat{x}`$ direction, since rotational symmetry can eventually be used to extend our results to those vectors $`\stackrel{}{s}`$ with arbitrary orientations. We begin by identifying those pairs of points, $`\stackrel{}{r}_1`$ and $`\stackrel{}{r}_2`$, which satisfy $`\stackrel{}{s}=s\widehat{x}`$. One set of random points for $`\stackrel{}{r}_1`$ is located in $`A_1`$ and the other set of random points for $`\stackrel{}{r}_2`$ is located in $`A_2`$ as shown in Fig. 1. We observe that $`A_2`$ is the overlap area between the original circle $`C_1`$ and an identical circle $`C_2`$ whose center is shifted from $`(0,0)`$ to $`(\left|\stackrel{}{s}\right|,0)`$ as shown in Fig. 2. Since the areas of $`A_1`$ and $`A_2`$ are equal, it follows that the probability density of finding a given $`s=\left|s\widehat{x}\right|`$ in a circle of uniform density is proportional to the area of $`A_2`$:

$$_{sR}^{\frac{s}{2}}𝑑x_{\sqrt{R^2\left(xs\right)^2}}^{\sqrt{R^2\left(xs\right)^2}}𝑑y+_{\frac{s}{2}}^R𝑑x_{\sqrt{R^2x^2}}^{\sqrt{R^2x^2}}𝑑y.$$
(8)
Using rotational symmetry, this result can apply to any orientation of $`\stackrel{}{s}`$, where $`0\varphi 2\pi `$, and hence the probability density $`P_2\left(s\right)`$ for a circle of uniform density can be factorized as

$$P_2\left(s\right)=2\pi s\times f\left(s\right),$$
(9)
where $`f\left(s\right)`$ is a function to be determined. If we impose the normalization requirement

$$_0^{2R}P_2(s)𝑑s=1,$$
(10)
we then have

$`P_2(s)`$ $`=`$ $`{\displaystyle \frac{2\pi s_{\frac{s}{2}}^R𝑑x_{\sqrt{R^2x^2}}^{\sqrt{R^2x^2}}𝑑y}{_0^{2R}\left(2\pi s_{\frac{s}{2}}^R𝑑x_{\sqrt{R^2x^2}}^{\sqrt{R^2x^2}}𝑑y\right)𝑑s}}`$ (11)
$`=`$ $`{\displaystyle \frac{2s}{R^2}}{\displaystyle \frac{s^2}{\pi R^4}}\sqrt{4R^2s^2}{\displaystyle \frac{4s}{\pi R^2}}\mathrm{sin}^1\left(s/2R\right).`$ (12)
Equation (12) is identical to the results obtained in Refs. Kendall and Moran (1963); Santaló (1976) by other means.

The conclusion that emerges from this formalism is that the probability density of finding two random points separated by a random vector $`\stackrel{}{s}`$ in a circle of uniform density can be derived by simply calculating the overlap region of that circle with an identical circle obtained by shifting the center from the origin to $`\stackrel{}{s}`$. In the following discussion, we show that this result generalizes to higher dimensions, and provides a simple way of calculating $`P_n(s)`$ for $`n3`$.

The above formalism for a circle of uniform density can be extended to a sphere of uniform density. For a given $`s`$, we select the positive $`\widehat{z}`$ direction and study the distribution of the random vectors $`\stackrel{}{s}`$ in this direction. It then follows that the probability density of finding $`s\widehat{z}`$ is proportional to

$$_{sR}^{\frac{s}{2}}𝑑z_{\sqrt{R^2\left(zs\right)^2}}^{\sqrt{R^2\left(zs\right)^2}}𝑑x_{\sqrt{R^2\left(zs\right)^2x^2}}^{\sqrt{R^2\left(zs\right)^2x^2}}𝑑y+_{\frac{s}{2}}^R𝑑z_{\sqrt{R^2z^2}}^{\sqrt{R^2z^2}}𝑑x_{\sqrt{R^2z^2x^2}}^{\sqrt{R^2z^2x^2}}𝑑y.$$
(13)
In a $`3`$-dimensional space, we note that the direction of the random vector $`\stackrel{}{s}`$ can have the following range: $`0\theta \pi `$ and $`0\varphi 2\pi `$. Following the previous discussion, we thus arrive at the following expression for a sphere with a uniform density distribution:

$`P_3(s)`$ $`=`$ $`{\displaystyle \frac{4\pi s^2\left(_{\frac{s}{2}}^R𝑑z_{\sqrt{R^2z^2}}^{\sqrt{R^2z^2}}𝑑x_{\sqrt{R^2z^2x^2}}^{\sqrt{R^2z^2x^2}}𝑑y\right)}{_0^{2R}4\pi s^2\left(_{\frac{s}{2}}^R𝑑z_{\sqrt{R^2z^2}}^{\sqrt{R^2z^2}}𝑑x_{\sqrt{R^2z^2x^2}}^{\sqrt{R^2z^2x^2}}𝑑y\right)𝑑s}}`$ (14)
$`=`$ $`3{\displaystyle \frac{s^2}{R^3}}{\displaystyle \frac{9}{4}}{\displaystyle \frac{s^3}{R^4}}+{\displaystyle \frac{3}{16}}{\displaystyle \frac{s^5}{R^6}}.`$ (15)
The result in Eq. (15) agrees exactly with the expression obtained previously in Refs. Kendall and Moran (1963); Santaló (1976); Overhauser ; Deltheil (1919); Hammersley (1950); Lord (1954).

The present formalism can be readily generalized to express $`P_n\left(s\right)`$ for a sphere of uniform density in $`n`$-dimensions as

$$P_n(s)=\frac{s^{n1}_{s/2}^R𝑑x_n_{\sqrt{R^2x_n^2}}^{\sqrt{R^2x_n^2}}𝑑x_1\mathrm{}_{\sqrt{R^2x_n^2\mathrm{}x_{n2}^2}}^{\sqrt{R^2x_n^2\mathrm{}x_{n2}^2}}𝑑x_{n1}}{_0^{2R}\left(s^{n1}_{s/2}^R𝑑x_n_{\sqrt{R^2x_n^2}}^{\sqrt{R^2x_n^2}}𝑑x_1\mathrm{}_{\sqrt{R^2x_n^2\mathrm{}x_{n2}^2}}^{\sqrt{R^2x_n^2\mathrm{}x_{n2}^2}}𝑑x_{n1}\right)𝑑s}.$$
(16)
We find that if $`n`$ is an even number,

$$P_n(s)=\frac{n\times s^{n1}}{R^n}\left[\frac{2}{\pi }\mathrm{cos}^1\left(s/2R\right)\frac{s}{\pi }\underset{k=1}{\overset{\frac{n}{2}}{}}\frac{(n2k)!!}{(n2k+1)!!}\left(R^2s^2/4\right)^{\frac{n2k+1}{2}}R^{2k2n}\right],$$
(17)
where $`0!=0!!=1`$. If $`n`$ is an odd number,

$$P_n(s)=\frac{n\times s^{n1}}{R^n}\frac{n!!}{(n1)!!}\underset{k=0}{\overset{\frac{n1}{2}}{}}\frac{(1)^k}{2k+1}\frac{\left(\frac{n1}{2}\right)!}{k!\left(\frac{n1}{2}k\right)!}\left[1\left(s/2R\right)^{2k+1}\right].$$
(18)
The functional forms of Eq. (16) for $`n=1`$, $`2`$, $`3`$, and $`4`$ are shown in Fig. 3.

The cumulative distribution function (CDF) Feller (1967); Gardiner (1985); van Kampen (1992) for $`P_n(s)`$ is given by

$`D_n(x)`$ $`=`$ $`{\displaystyle _0^x}P_n(s)𝑑s`$ (19)
$`=`$ $`{\displaystyle \frac{x^n}{R^n}}{\displaystyle \frac{B_\alpha (\frac{1}{2},\frac{n}{2}+\frac{1}{2})}{B(\frac{1}{2},\frac{n}{2}+\frac{1}{2})}}\times {\displaystyle \frac{x^n}{R^n}}+2^n\times {\displaystyle \frac{B_\alpha (\frac{n}{2}+\frac{1}{2},\frac{n}{2}+\frac{1}{2})}{B(\frac{1}{2},\frac{n}{2}+\frac{1}{2})}},`$
where $`0x2R`$, $`\alpha =x^2/4R^2`$, and $`B_\alpha `$ is the incomplete beta function.

We summarize three important representations for the probability density function $`P_n\left(s\right)`$ for a spherical $`n`$-ball of radius $`R`$ with a uniform density distribution as follows:

1. Integral representation:

$$P_n(s)=\frac{s^{n1}_{\frac{s}{2}}^R\left(R^2x^2\right)^{\frac{n1}{2}}𝑑x}{\frac{1}{2n}B(\frac{n}{2}+\frac{1}{2},\frac{1}{2})R^{2n}}.$$
(20)
2. Generating function representation:

$$P_n(s)=\frac{s^{n1}\frac{1}{n!}\left(\frac{}{h}\right)_{h=0}^n\frac{1}{\sqrt{1h^2}}\left[\mathrm{sin}^1h\mathrm{sin}^1\left(\frac{2h\sqrt{4s^2}}{2h\sqrt{4s^2}}\right)\right]}{\frac{1}{2n}B(\frac{n}{2}+\frac{1}{2},\frac{1}{2})},$$
(21)
where $`\left|h\right|<1`$ and $`R=1`$. We note that one can obtain a number of identities and recursion relations for $`P_n\left(s\right)`$ from the generating function representation in Eq. (21).
3. Hypergeometric function representation:

$$P_n(s)=\frac{2n}{B(\frac{n}{2}+\frac{1}{2},\frac{1}{2})}\times \left[{}_{2}{}^{}F_{1}^{}(a,b,c,\alpha )R\frac{s}{2}_2F_1(a,b,c,\beta )\right]\times \frac{s^{n1}}{R^{n+1}},$$
(22)
where $`a=1/2`$, $`b=1/2n/2`$, $`c=3/2`$, $`\alpha =1`$, $`\beta =s^2/4R^2`$, and $`{}_{2}{}^{}F_{1}^{}\left(\mathrm{}\right)`$ is the hypergeometric function Arfken and Weber (1995). Notice that we can obtain the orthogonality relations for $`P_n\left(s\right)`$ from the hypergeometric function representation as shown in Eq. (22).

In Ref. Tu (2001) additional identities and recursion relations for $`P_n\left(s\right)`$ are discussed in greater detail.

## III Spherically Symmetric Density Distributions

In this section we extend the previous results to the case of a circle with a variable (but spherically symmetric) density characterized by a density function $`\rho \left(r\right)`$. Following the derivation presented in the previous section, we note that for any random vector $`\stackrel{}{s}=s\widehat{x}`$, if the second random point $`\stackrel{}{r}_2`$ carries the density information $`\rho (x,y)`$, then the first random point $`\stackrel{}{r}_1`$ should have the density information $`\rho (xs,y).`$ It then follows $`P_2(s)`$ can be expressed as

$`P_2(s)`$ $`=`$ $`{\displaystyle \frac{s_{\frac{s}{2}}^R𝑑x_{\sqrt{R^2x^2}}^{\sqrt{R^2x^2}}\rho (xs,y)\times \rho (x,y)𝑑y}{_0^{2R}\left(s_{\frac{s}{2}}^R𝑑x_{\sqrt{R^2x^2}}^{\sqrt{R^2x^2}}\rho (xs,y)\times \rho (x,y)𝑑y\right)𝑑s}}.`$ (23)
A general formula for $`P_n\left(s\right)`$ for an $`n`$-dimensional spherical ball of radius $`R`$ having a spherically symmetric density can be derived, and we find

$$P_n(s)=\frac{s^{n1}_{\frac{s}{2}}^R𝑑x_n_{\sqrt{R^2x_n^2}}^{\sqrt{R^2x_n^2}}𝑑x_1\mathrm{}_{\sqrt{R^2x_n^2x_1^2\mathrm{}x_{n1}^2}}^{\sqrt{R^2x_n^2x_1^2\mathrm{}x_{n1}^2}}\rho _1\times \rho _2𝑑x_{n1}}{_0^{2R}\left(s^{n1}_{\frac{s}{2}}^R𝑑x_n_{\sqrt{R^2x_n^2}}^{\sqrt{R^2x_n^2}}𝑑x_1\mathrm{}_{\sqrt{R^2x_n^2x_1^2\mathrm{}x_{n1}^2}}^{\sqrt{R^2x_n^2x_1^2\mathrm{}x_{n1}^2}}\rho _1\times \rho _2𝑑x_{n1}\right)𝑑s},$$
(24)
where

$`\rho _1`$ $`=`$ $`\rho (x_1,x_2,\mathrm{},x_ns),`$
$`\rho _2`$ $`=`$ $`\rho (x_1,x_2,\mathrm{},x_n).`$ (25)
Some analytical results for a sphere with various spherically symmetric density distributions can be found in Ref. Tu (2001).

As an application of Eq. (24) we consider the case of an $`n`$-dimensional spherical space of radius $`R\mathrm{}`$ with a Gaussian density given by

$$\rho _n(r)=\frac{N}{(2\pi )^{\frac{n}{2}}\sigma ^n}e^{\frac{1}{2}\frac{r^2}{\sigma ^2}},$$
(26)
where

$$N=\underset{R\mathrm{}}{lim}n\frac{\pi ^{\frac{n}{2}}}{\mathrm{\Gamma }\left(\frac{n}{2}+1\right)}_0^R\rho _n(r)r^{n1}𝑑r.$$
(27)
In Eq. (27) $`r`$ is measured from the center, and the integral is over all space. The PDF for an $`n`$-dimensional spherical Gaussian space can then be obtained as

$$P_n(s)=\frac{s^{n1}e^{\frac{s^2}{4\sigma ^2}}}{2^{n1}\mathrm{\Gamma }\left(\frac{n}{2}\right)\sigma ^n}.$$
(28)
Fig. 4 displays the functions in Eq. (28) for $`n=3`$, $`4`$, $`5`$, and $`6`$, where $`\sigma =1`$. Finally, we note that the maximum probability density, denoted by $`s_{\mathrm{max}}`$, occurs at

$$s_{\mathrm{max}}=\sqrt{2(n1)}\sigma .$$
(29)
## IV Arbitrary Density Distributions

We consider in this section the probability density functions for a spherical $`n`$-ball having an arbitrary density characterized by a density function $`\rho `$. We begin with a circle and use the conventional notation for polar coordinates: $`x=r\mathrm{cos}\varphi `$ and $`y=r\mathrm{sin}\varphi `$. In a $`2`$-dimensional space, the random vector $`\stackrel{}{s}`$ can be characterized by an angle $`\varphi `$ in the range $`0\varphi 2\pi `$. Associate each random unit vector $`\widehat{s}\left(\varphi \right)`$ with a rotation operator $``$ such that

$$|\widehat{x}=|\widehat{s}\left(\varphi \right).$$
(30)
To ensure that the product of $`\rho \left(\stackrel{}{r}_1\right)`$ and $`\rho \left(\stackrel{}{r}_2\right)`$ maintains the correct density information, we use a $`2\times 2`$ matrix

$$R_{2\times 2}(\varphi )=\left[\begin{array}{cc}\mathrm{cos}\varphi & \mathrm{sin}\varphi \\ \mathrm{sin}\varphi & \mathrm{cos}\varphi \end{array}\right]$$
(31)
to characterize this particular operator $``$ such that

$$\rho (x,y)\rho (\mathrm{cos}\varphi x\mathrm{sin}\varphi y,\mathrm{sin}\varphi x+\mathrm{cos}\varphi y).$$
(32)
Notice that $`R_{2\times 2}(\varphi )`$ is an orthogonal matrix which satisfies $`R_{2\times 2}^1(\varphi )=R_{2\times 2}^T(\varphi )`$ and its determinant is $`+1`$, where $`T`$ denotes the transpose.

We can then express $`P_2\left(s\right)`$ for a circle with an arbitrary density distribution as

$$P_2(s)=\frac{s_0^{2\pi }𝑑\varphi _{s/2}^R𝑑x_{\sqrt{R^2x^2}}^{\sqrt{R^2x^2}}\rho (x^{},y^{})\times \rho (x^{\prime \prime },y^{\prime \prime })𝑑y}{_0^{2R}\left(s_0^{2\pi }𝑑\varphi _{s/2}^R𝑑x_{\sqrt{R^2x^2}}^{\sqrt{R^2x^2}}\rho (x^{},y^{})\times \rho (x^{\prime \prime },y^{\prime \prime })𝑑y\right)𝑑s}$$
(33)
where

$`\rho (x^{},y^{})`$ $`=`$ $`\rho (\mathrm{cos}\varphi (xs)\mathrm{sin}\varphi y,\mathrm{sin}\varphi (xs)+\mathrm{cos}\varphi y),`$
$`\rho (x^{\prime \prime },y^{\prime \prime })`$ $`=`$ $`\rho (\mathrm{cos}\varphi x\mathrm{sin}\varphi y,\mathrm{sin}\varphi x+\mathrm{cos}\varphi y).`$
Figure 5 exhibits $`P_2(s)`$ when $`R=1`$, and illustrates the agreement between the Monte Carlo simulation and the analytical result.

The preceding discussion can be generalized to the case of a spherical $`n`$-ball with an arbitrary density. Define the following hyperspherical coordinates Fleming (1977),

$`x_1`$ $`=`$ $`r\mathrm{sin}\theta _{n2}\mathrm{sin}\theta _{n3}\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{sin}\theta _2\mathrm{sin}\theta _1\mathrm{cos}\varphi ,`$
$`x_2`$ $`=`$ $`r\mathrm{sin}\theta _{n2}\mathrm{sin}\theta _{n3}\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{sin}\theta _2\mathrm{sin}\theta _1\mathrm{sin}\varphi ,`$
$`\mathrm{}`$
$`x_i`$ $`=`$ $`r\mathrm{sin}\theta _{n2}\mathrm{sin}\theta _{n3}\mathrm{}\mathrm{}\mathrm{sin}\theta _{i1}\mathrm{cos}\theta _{i2},`$
$`\mathrm{}`$
$`x_{n1}`$ $`=`$ $`r\mathrm{sin}\theta _{n2}\mathrm{cos}\theta _{n3},`$
$`x_n`$ $`=`$ $`r\mathrm{cos}\theta _{n2}.`$ (34)
where $`\theta _i[0,\pi ]`$ and $`\varphi [0,2\pi ]`$. Associate each random unit vector $`\widehat{s}(\theta _1,\mathrm{},\varphi )`$ with a rotation operator $``$ such that

$$|\widehat{x}_n=|\widehat{s}(\theta _1,\mathrm{},\varphi ).$$
(35)
The $`n\times n`$ matrix representation for the rotation operator $``$ in Eq. (35) can be expressed as

$$R_{n\times n}(\theta _{n2},\mathrm{},\theta _1,\varphi )=R_{n\times n}(\varphi )\times R_{n\times n}(\theta _1)\times \mathrm{}\times R_{n\times n}(\theta _{n2}),$$
(36)
where $`n3`$. The various matrices on the right-hand side of Eq. (36) are defined as follows:

1. $`R_{n\times n}(\varphi )`$: The matrix elements are $`a_{11}=\mathrm{cos}\varphi `$, $`a_{12}=\mathrm{sin}\varphi `$, $`a_{21}=\mathrm{sin}\varphi `$, $`a_{22}=\mathrm{cos}\varphi `$, for $`i[1,2]`$ and $`j[3,n]`$ $`a_{ij}=0`$, for $`i[3,n]`$ and $`i[1,2]`$ $`a_{ij}=0`$, for $`i[3,n]`$ and $`j[3,n]`$ $`a_{ij}=\delta _{ij}`$.
2. $`R_{n\times n}\left(\theta _1\right)`$: The matrix elements are $`a_{11}=\mathrm{cos}\theta _1`$, $`a_{12}=0`$, $`a_{13}=\mathrm{sin}\theta _1`$, $`a_{21}=0`$, $`a_{22}=1`$, $`a_{23}=0`$, $`a_{31}=\mathrm{sin}\theta _1`$, $`a_{32}=0`$, $`a_{33}=\mathrm{cos}\theta _1`$, for $`i[1,3]`$ and $`j[4,n]`$ $`a_{ij}=0`$, for $`i[4,n]`$ and $`j[1,3]`$ $`a_{ij}=0`$, for $`i[4,n]`$ and $`j[4,n]`$ $`a_{ij}=\delta _{ij}`$.
3. $`R_{n\times n}\left(\theta _k\right)`$ and $`k[2,n3]`$: The matrix elements are $`a_{k+1,k+1}=\mathrm{cos}\theta _k`$, $`a_{k+1,k+2}=\mathrm{sin}\theta _k`$, $`a_{k+2,k+1}=\mathrm{sin}\theta _k`$, $`a_{k+2,k+2}=\mathrm{cos}\theta _k`$, for $`i[1,k]`$ and $`j[1,k]`$ $`a_{ij}=\delta _{ij}`$, for $`i[k+3,n]`$ and $`j[k+3,n]`$ $`a_{ij}=\delta _{ij}`$, for $`i[1,k]`$ and $`j[k+1,n]`$ $`a_{ij}=0`$, for $`i[k+1,n]`$ and $`j[1,k]`$ $`a_{ij}=0`$, for $`i[k+1,k+2]`$ and $`j[k+3,n]`$ $`a_{ij}=0`$, for $`i[k+3,n]`$ and $`j[k+1,k+2]`$ $`a_{ij}=0`$.
4. $`R_{n\times n}\left(\theta _{n2}\right)`$: The matrix elements are $`a_{n1,n1}=\mathrm{cos}\theta _{n2}`$, $`a_{n1,n}=\mathrm{sin}\theta _{n2}`$, $`a_{n,n1}=\mathrm{sin}\theta _{n2}`$, $`a_{n,n}=\mathrm{cos}\theta _{n2}`$, for $`i[1,n2]`$ and $`j[1,n2]`$ $`a_{ij}=\delta _{ij}`$, for $`i[1,n2]`$ and $`j[n1,n]`$ $`a_{ij}=0`$, for $`i[n1,n]`$ and $`j[1,n2]`$ $`a_{ij}=0`$.

Notice that all the matrices in Eq. (36) are orthogonal and their determinants are $`+1`$. The matrix elements for $`R_{n\times n}(\theta _{n2},\mathrm{},\theta _1,\varphi )`$ can be summarized as follows:

1. The matrix elements for the $`1`$st column are $`a_{11}=\mathrm{cos}\theta _1\mathrm{cos}\varphi `$, $`a_{21}=\mathrm{cos}\theta _1\mathrm{sin}\varphi `$, $`a_{31}=\mathrm{sin}\theta _1`$, and for $`i[4,n]`$ $`a_{i1}=0`$.
2. The matrix elements for the $`2`$nd column are $`a_{12}=\mathrm{sin}\varphi `$, $`a_{22}=\mathrm{cos}\varphi `$, and for $`i[3,n]`$ $`a_{i2}=0`$.
3. The matrix elements for the $`j`$th column where $`j[3,n1]`$ are $`a_{j+1,j}=\mathrm{sin}\theta _{j1}`$, for $`i[j+2,n]`$ $`a_{ij}=0`$, and for $`i[1,j]`$ $`a_{ij}=\mathrm{cos}\theta _{j1}\times x_i`$, where $`x_i=\widehat{r}\times \widehat{x}_i`$ and $`x_i`$ is the $`i`$th Cartesian coordinate component for a unit vector $`\widehat{r}`$ in $`j`$-dimensions. For example:

$`a_{14}`$ $`=`$ $`\mathrm{cos}\theta _3\mathrm{sin}\theta _2\mathrm{sin}\theta _1\mathrm{cos}\varphi ,`$
$`a_{24}`$ $`=`$ $`\mathrm{cos}\theta _3\mathrm{sin}\theta _2\mathrm{sin}\theta _1\mathrm{sin}\varphi ,`$
$`a_{34}`$ $`=`$ $`\mathrm{cos}\theta _3\mathrm{sin}\theta _2\mathrm{cos}\theta _1,`$
$`a_{44}`$ $`=`$ $`\mathrm{cos}\theta _3\mathrm{cos}\theta _2.`$
4. The matrix elements for the $`n`$th column are for $`i[1,n]`$ $`a_{in}=x_i`$, where $`x_i=\widehat{r}\times \widehat{x}_i`$ and $`x_i`$ is the $`i`$th Cartesian coordinate component for a unit vector $`\widehat{r}`$ in $`n`$-dimensions.

Additionally, it is convenient to define the following transforms:

$`x_i^{}`$ $`=`$ $`{\displaystyle \underset{j=1}{\overset{n}{}}}a_{ij}\left(x_j\delta _{jn}s\right),`$ (37)
$`x_i^{\prime \prime }`$ $`=`$ $`{\displaystyle \underset{j=1}{\overset{n}{}}}a_{ij}x_j,`$ (38)
where $`a_{ij}`$ are the matrix elements for $`R_{n\times n}(\theta _{n2},\mathrm{},\theta _1,\varphi )`$ defined in Eq. (36). The master equation of $`P_n(s)`$ for a spherical $`n`$-ball with an arbitrary density characterized by a density function $`\rho \left(𝒙\right)=\rho (x_1,x_2,\mathrm{},x_n)`$ can then be formulated as

$$P_n(s)=\frac{s^{n1}\times 𝐅_n(x,\theta )}{_0^{2R}\left[s^{n1}\times 𝐅_n(x,\theta )\right]𝑑s},$$
(39)
where

$`𝐅_n(x,\theta )`$ $`=`$ $`{\displaystyle _0^\pi }\mathrm{sin}^{n2}\theta _{n2}d\theta _{n2}{\displaystyle _0^\pi }\mathrm{sin}^{n3}\theta _{n3}d\theta _{n3}\mathrm{}`$
$`\mathrm{}\mathrm{}{\displaystyle _0^\pi }\mathrm{sin}\theta _1d\theta _1{\displaystyle _0^{2\pi }}𝑑\varphi {\displaystyle _{s/2}^R}𝑑x_n{\displaystyle _{\sqrt{R^2x_n^2}}^{\sqrt{R^2x_n^2}}}𝑑x_1\mathrm{}`$
$`\mathrm{}\mathrm{}\mathrm{}{\displaystyle _{\sqrt{R^2x_n^2\mathrm{}x_{n2}^2}}^{\sqrt{R^2x_n^2\mathrm{}x_{n2}^2}}}\rho \left(𝒙^{}\right)\times \rho \left(𝒙^{\prime \prime }\right)𝑑x_{n1},`$
and

$`\rho \left(𝒙^{}\right)`$ $`=`$ $`\rho (x_1^{},x_2^{},\mathrm{},x_n^{}),`$
$`\rho \left(𝒙^{\prime \prime }\right)`$ $`=`$ $`\rho (x_1^{\prime \prime },x_2^{\prime \prime },\mathrm{},x_n^{\prime \prime }).`$
The probability density function for the random distance distribution including Euclidean distance and geodesic distance for an $`n`$-dimensional sphere (i.e. the boundary of a spherical $`n`$-ball) with an arbitrary surface density distribution will be discussed in Ref. Tu (a).

## V Applications

### V.1 $`m`$th Moment of the Distance

In some applications the $`m`$th moment of the distance, rather than the distance itself, is of interest. As an example, for a collection of nucleons interacting via simple harmonic oscillator potentials, $`s^2`$ may be of interest rather than $`s`$ itself. We then calculate the $`m`$th moment $`s^m`$ for the case of a spherical $`n`$-ball of uniform density, where

$$s^m=_0^{2R}s^mP_n(s)𝑑s.$$
(40)
Using Eq. (20), the $`m`$th moment $`s^m`$ has the general form

$$s^m=2^{n+m}\left(\frac{n}{n+m}\right)\frac{B(\frac{n}{2}+\frac{1}{2},\frac{n}{2}+\frac{1}{2}+\frac{m}{2})}{B(\frac{n}{2}+\frac{1}{2},\frac{1}{2})}R^m,$$
(41)
where $`B(p,q)`$ is the beta function, and $`m\left(n1\right)`$.

Additionally, $`s^m`$ can be evaluated for a spherical space having a Gaussian density distribution where $`\rho e^{r^2/2\sigma ^2}`$. From Eq. (28) we find,

$$s^m=\underset{R\mathrm{}}{lim}_0^{2R}s^mP_n(s)𝑑s=(2\sigma )^m\frac{\mathrm{\Gamma }\left(\frac{n+m}{2}\right)}{\mathrm{\Gamma }\left(\frac{n}{2}\right)}.$$
(42)
In some applications (such as in nuclear physics) involving low-energy interactions among nucleons the lower limit (zero) should be replaced by the hard-core radius $`r_c0.5\times 10^{13}`$ cm Bohr and Mottelson (1969). In such cases the expressions for $`P_n(s)`$ and $`s^m`$ can be expressed as follows:

$$P_n(s)=\frac{s^{n1}_{\frac{s}{2}}^R\left(R^2x^2\right)^{\frac{n1}{2}}𝑑x}{C(2R,0,n)C(r_c,0,n)}$$
(43)
and

$$s^m=_{r_c}^{2R}s^mP_n(s)𝑑s=\frac{H(R,r_c,m,n)}{H(R,r_c,0,n)},$$
(44)
where

$$C(a,m,n)=_0^as^{m+n1}𝑑s_{s/2}^R\left(R^2x^2\right)^{\left(n1\right)/2}𝑑x,$$

and

$`H(R,r_c,m,n)`$ $`=`$ $`{\displaystyle \frac{(2R)^{n+m}}{n+m}}\left[B({\displaystyle \frac{n}{2}}+{\displaystyle \frac{1}{2}},{\displaystyle \frac{n}{2}}+{\displaystyle \frac{1}{2}}+{\displaystyle \frac{m}{2}})B_x({\displaystyle \frac{n}{2}}+{\displaystyle \frac{1}{2}},{\displaystyle \frac{n}{2}}+{\displaystyle \frac{1}{2}}+{\displaystyle \frac{m}{2}})\right]`$ (45)
$`{\displaystyle \frac{r_c^{n+m}}{n+m}}\left[B({\displaystyle \frac{1}{2}},{\displaystyle \frac{n}{2}}+{\displaystyle \frac{1}{2}})B_x({\displaystyle \frac{1}{2}},{\displaystyle \frac{n}{2}}+{\displaystyle \frac{1}{2}})\right],`$
with $`x=r_c^2/4R^2`$, and $`m`$ an integer. We note that Eqs. (43) and (44) are the first known analytical results for $`P_n(s)`$ and $`s^m`$ to incorporate the hard-core radius $`r_c`$.

### V.2 Neutrino-Pair Exchange Interactions

A second example of interest is the $`\nu \overline{\nu }`$-exchange (neutrino-pair exchange) contribution to the self energy of a nucleus or a neutron star. For two point masses the $`2`$-body potential energy is given by Fischbach (1996); Feinberg and Sucher (1968); G. Feinberg, J. Sucher, and C. K. Au (1989); Hsu and Sikivie (1994); Woodahl and Fischbach (2000a, b, 2001)

$$V_{\nu \overline{\nu }}\left(\left|\stackrel{}{r}_i\stackrel{}{r}_j\right|\right)=\frac{G_F^2a_ia_j}{4\pi ^3\left|\stackrel{}{r}_i\stackrel{}{r}_j\right|^5},$$
(46)
where $`a_i`$ and $`a_j`$ are coupling constants which characterize the strength of the neutrino coupling to fermions $`i`$ and $`j`$ ($`i`$, $`j`$ = electron, proton, or neutron). In the standard model Partical Data Group (2000),

$`a_e`$ $`=`$ $`{\displaystyle \frac{1}{2}}+2\mathrm{sin}^2\theta _W=0.964`$
$`a_p`$ $`=`$ $`{\displaystyle \frac{1}{2}}2\mathrm{sin}^2\theta _W=0.036`$
$`a_n`$ $`=`$ $`{\displaystyle \frac{1}{2}}.`$
As an example, we find for the case of a uniform density distribution of radius $`R`$ containing $`N`$ neutrons,

$$W_3=\frac{N(N1)}{2}\left(\frac{3}{2r_c^2R^3}\frac{9}{4r_cR^4}+\frac{9}{8R^5}\frac{3r_c}{16R^6}\right)\frac{G_F^2}{4\pi ^3},$$
(47)
where $`r_c`$ is the hard-core radius. The analogous result for a Gaussian density distribution is

$$W_3=\left[\frac{1}{r_c^2}e^{r_c^2/4\sigma ^2}\frac{\mathrm{\Gamma }(0,r_c^2/4\sigma ^2)}{4\sigma ^2}\right]\frac{N(N1)G_F^2}{32\sigma ^3\pi ^{7/2}},$$
(48)
where $`\mathrm{\Gamma }(a,b)`$ is the incomplete gamma function.

### V.3 Neutron Star Models

Another application of current interest is the self-energy of neutron star arising from the exchange of $`\nu \overline{\nu }`$ pairs Woodahl and Fischbach (2000a, b, 2001). Here we evaluate the probability density functions in $`3`$ dimensions for neutron stars with a multiple-shell density distribution, which is what is typically assumed in neutron star models Pines (1980); Shapiro and Teukolsky (1983). For illustrative purposes, we discuss a spherically symmetric model with $`2`$ spherical shells, each of uniform density, where for simplicity we assume shells of equal thickness. Some other multiple-shell models and their $`n`$-dimensional probability density functions can be found in Tu (2001).

For a $`2`$-shell model with a uniform density in each shell define $`\rho =\rho _1`$ for $`0rR/2`$ and $`\rho =\rho _2`$ for $`R/2rR`$, where $`\rho _1`$ and $`\rho _2`$ are arbitrary constants and $`r`$ is measured from the center of the neutron star. Using the preceding formalism we can show that $`P_3\left(s\right)`$ has $`4`$ different functional forms specified by $`4`$ regions:

1. $`0s\frac{1}{2}R:`$

$$P_3(s)=\frac{24(\rho _1^2+7\rho _2^2)s^2}{\left(\rho _1+7\rho _2\right)^2R^3}\frac{36(\rho _1^22\rho _1\rho _2+5\rho _2^2)s^3}{\left(\rho _1+7\rho _2\right)^2R^4}+\frac{12(\rho _1^22\rho _1\rho _2+2\rho _2^2)s^5}{\left(\rho _1+7\rho _2\right)^2R^6},$$
(49)
2. $`\frac{1}{2}RsR:`$

$`P_3(s)`$ $`=`$ $`{\displaystyle \frac{81(\rho _1\rho _2)\rho _2s}{2\left(\rho _1+7\rho _2\right)^2R^2}}+{\displaystyle \frac{24\rho _1s^2}{\left(\rho _1+7\rho _2\right)R^3}}{\displaystyle \frac{36\rho _1(\rho _1+3\rho _2)s^3}{\left(\rho _1+7\rho _2\right)^2R^4}}`$ (50)
$`+{\displaystyle \frac{12\rho _1^2s^5}{\left(\rho _1+7\rho _2\right)^2R^6}},`$
3. $`Rs\frac{3}{2}R:`$

$`P_3(s)`$ $`=`$ $`{\displaystyle \frac{81(\rho _1\rho _2)\rho _2s}{2\left(\rho _1+7\rho _2\right)^2R^2}}+{\displaystyle \frac{24(9\rho _1\rho _2)\rho _2s^2}{\left(\rho _1+7\rho _2\right)^2R^3}}{\displaystyle \frac{36(5\rho _1\rho _2)\rho _2s^3}{\left(\rho _1+7\rho _2\right)^2R^4}}`$ (51)
$`+{\displaystyle \frac{12(2\rho _1\rho _2)\rho _2s^5}{\left(\rho _1+7\rho _2\right)^2R^6}},`$
4. $`\frac{3}{2}Rs2R:`$

$$P_3(s)=\frac{192\rho _2^2s^2}{\left(\rho _1+7\rho _2\right)^2R^3}\frac{144\rho _2^2s^3}{\left(\rho _1+7\rho _2\right)^2R^4}+\frac{12\rho _2^2s^5}{\left(\rho _1+7\rho _2\right)^2R^6}.$$
(52)

We observe that the probability density functions defined in adjacent regions are continuous across the boundaries separating the regions.

### V.4 RIPS: A New Test for Random Number Generators in $`n`$-Dimensions

Another interesting application of the present work is as a new test of random number generators in $`n`$-dimensions. Our statistical method follows by applying the probability density functions for the random distance distribution to evaluate the expectation values of $`\stackrel{}{r}_{12}\stackrel{}{r}_{23}`$ in $`n`$-dimensions, where $`\stackrel{}{r}_{12}=\stackrel{}{r}_2\stackrel{}{r}_1`$, $`\stackrel{}{r}_{23}=\stackrel{}{r}_3\stackrel{}{r}_2`$, and $`\stackrel{}{r}_1`$, $`\stackrel{}{r}_2`$, and $`\stackrel{}{r}_3`$ are three random points independently sampled from a spherical $`n`$-ball. The quantity $`\stackrel{}{r}_{12}\stackrel{}{r}_{23}_n`$ is one of the geometric probability constants Tu (2001) and has applications in many areas of physics Fischbach (1996); Tu (2001). It follows from the preceding formalism that for a spherical $`n`$-ball of radius $`R`$ with a uniform density distribution,

$$\stackrel{}{r}_{12}\stackrel{}{r}_{23}_n=\frac{n}{n+2}R^2.$$
(53)
Note that

$$\underset{n\mathrm{}}{lim}\stackrel{}{r}_{12}\stackrel{}{r}_{23}_n=R^2.$$
(54)
In a spherical Gaussian space where $`\rho e^{r^2/2\sigma ^2}`$ and $`R\mathrm{}`$,

$$\stackrel{}{r}_{12}\stackrel{}{r}_{23}_n=n\sigma ^2.$$
(55)
We refer to this test as RIPS (Random Inner Product in a Sphere).

We now apply the RIPS test to check three popular random number generators frequently used in Monte Carlo simulations Knuth (1981); W. H. Press, S. A. Teukolsky, W. T. Vetterling, and B. P. Flannery (1992); I. Vattulainen, T. Ala-Nissila, and K. Kanlaala (1994); B. L. Holian, O. E. Percus, T. T. Warnock, and P. A. Whitlock (1994); K. V. Tretiakov and K. W. Wojciechowski (1999). The random number generators tested are:

1. RAN$`0`$ W. H. Press, S. A. Teukolsky, W. T. Vetterling, and B. P. Flannery (1992) which is a linear congruential generator and uses the following algorithm

$$I_n=16807\times I_{n1}\mathrm{mod}m=2^{31}1.$$
(56)
2. R$`31`$ I. Vattulainen, T. Ala-Nissila, and K. Kanlaala (1994) which uses the generalized feedback shift register (GFSR) method

$$x_n=x_{np}x_{nq,}$$
(57)
where $`p=31`$, $`q=3`$, and $``$ is the bitwise exclusive OR operator.
3. NWS B. L. Holian, O. E. Percus, T. T. Warnock, and P. A. Whitlock (1994); K. V. Tretiakov and K. W. Wojciechowski (1999) which uses the nested Weyl sequence

$$Y_n=\left\{n\left\{n\alpha \right\}\right\},$$
(58)
where $`\left\{y\right\}`$ is the fractional part of $`y`$ and $`\alpha `$ is an irrational number.

The results are shown in Table 1. We have checked several initialization methods and the results are not affected. We note that the NWS generator using the nested Weyl sequence method fails to pass our randomness tests in all dimensions selected. Hence caution should be exercised in using this particular random number generator for Monte Carlo simulations, especially in molecular dynamics simulations. More statistical tests for a variety of random number generators including the PSLQ algorithm using the binary digits of $`\pi `$ Bailey and Ferguson , and other newly proposed computational schemes for $`n`$-dimensions, will be discussed in greater detail in Ref. Tu (b).

The geometric probability constants, $`\stackrel{}{r}_{12}\stackrel{}{r}_{23}_n`$, can be added to the family of various computational tests for random number generators. They can then serve to investigate the quality of random number generators for questions of randomness, especially in higher dimensions ($`n>3`$) where few results are currently available. The possibility that some random number generators such as NWS pass other tests but not ours, may indicate that the properties of random points in a sphere provide a more sensitive test of randomness than is otherwise available. These and other issues will be discussed in more detail in Ref. Tu (b).

## VI Conclusions

A formalism has been presented in this paper for evaluating the analytical probability density function of the random distance distribution for a spherical $`n`$-ball with an arbitrary density distribution. We show that the random distance distribution technique can reduce otherwise difficult calculations from the complexity of $`2n`$-dimensional integrals to just a $`1`$-dimensional integral, even when the nucleon-nucleon hard-core radius $`r_c`$ is included. Our formalism has applications to the currently active area of research surrounding string-inspired theories of higher dimensional physics, and has numerous potential applications to other fields as well. Specifically the results presented here are of interest in the context of recent work on the modifications to the Newtonian inverse-square law arising from the existence of extra spatial dimensions N. Arkani-Hamed, S. Dimopoulos, and G. Dvali (1999). We have also presented a new computational method to test random number generators in $`n`$-dimensions used in Monte Carlo simulations.

## VII Acknowledgments

The authors wish to thank Chinh Le, Michelle Parry, David Schleef, Dave Seaman, Christopher Tong, and A. W. Overhauser for helpful discussions and the Purdue University Computing Center for computing support. This work was supported in part by the U.S. Department of Energy under Contract No. DE-AC02-76ER01428.