Search for a small egg by spermatozoa in restricted geometries

The search by swimmers for a small target in a bounded domain is ubiquitous in cellular biology, where a prominent case is that of the search by spermatozoa for an egg in the uterus. This is one of the severest selection processes in animal reproduction. We present here a mathematical model of the search, its analysis, and numerical simulations. In the proposed model the swimmers’ trajectories are rectilinear and the speed is constant. When a trajectory hits an obstacle or the boundary, it is reflected at a random angle and continues the search with the same speed. Because hitting a small target by a trajectory is a rare event, asymptotic approximations and stochastic simulations are needed to estimate the mean search time in various geometries. We consider searches in a disk, in convex planar domains, and in domains with cusps. The exploration of the parameter space for spermatozoa motion in different uterus geometries leads to scaling laws for the search process.


Introduction
Searches for a small target by random particle or by random switching dynamics depend critically on the local geometry of the target and the global geometry of the domain, as revealed recently by analytical, numerical methods and by stochastic simulations of molecular motion, viral trafficking, or prey search (Holcman and Schuss 2004;Redner 2001;Benichou et al. 2007;Schuss 2013, 2014;Ward et al. 2010;Kurella et al. 2015;Cheviakov et al. 2012;Schuss 2010). Asymptotic formulas for the expected search time were found for diffusive motion , but much less is known for other search processes. The problem under consideration here is that of estimating the mean search time of spermatozoa for an egg lodged in the uterus or in the fallopian tubes. This search is one of the most severe selection processes of one among hundreds of millions spermatozoa. How the winner is selected remains unclear and may depend on many parameters of the spermatozoa and the uterus. It is also unclear how to quantify geometrical parameters of the uterus. For example, how does a drop in the initial number of sperms affect fertility?
Between 1989 and 2005 the concentration of sperm cells in human semen has dropped significantly and continuously at an average rate of 1.9 %/year, leading to a reduction of 32.2 % in sperm count over 16 years (Rolland et al. 2013). Do these changes in concentration really matter? It is well documented that reducing the concentration of spermatozoa by four leads to infertility, but how such a result can be explained? because it still remains millions Reynaud et al. (2015). We still do not understand the role of the huge redundancy in sperm population. Motivated by these questions, we develop here a mathematical model of spermatozoa motion based on a simplified dynamics and on some elements of the uterus geometry. Our goal is to estimate the expected arrival time of a spermatozoon (a sperm cell) to a given neighborhood of a small egg. In that neighborhood, the motion of a spermatozoon is modified, because it responds to chemotaxis gradient generated by the egg (Armon et al. 2014;Armon and Eisenbach 2011;Teves et al. 2009). This terminal phase of the search is not studied here.
Classical models of spermatozoa motion include the beating of flagella in viscoelastic fluids (Fu et al. 2008) and attraction to a flat wall due to hydrodynamic interactions of the swimmer with the surface (Elgeti et al. 2010;Berke et al. 2008;Smith et al. 2011;Gaffney et al. 2011;Kantsler et al. 2014Kantsler et al. , 2013. For high spermatozoa concentration, collective modes of locomotion, different from those displayed by isolated cells, have been described by long-time kinematics of their relative locomotion (Michelin and Lauga 2010). The study of asymmetric flagellar bending was based on cytoplasmic calcium dynamics in the flagellum (Olson et al. 2011). The trajectories of spermatozoa can also be influenced by fluid motion (Marcos et al. 2012). However none of the present studies have addressed the generic question of the search of a small egg in the context of the uterus. Indeed, we shall explore here the role of the uterus geometry, quantified by various parameters such as the height, width, the radius of the cervix δ, local curvature near the fallopian tubes or the aspect ρ = L/W (as shown in Fig. 4) that measures the non-convexity, all should be key parameters in directing the spermatozoa toward the egg, yet this possibility has not been explored so far. However, contrary to the mentioned references, in the present work, we will use a very crude model of the spermatozoa motion, approximated as a ballistic directed motion, to derive asymptotic formula for the search process in convex geometries and use numerical simulations for non-convex ones. The novelty and difficulty here are in geometry of the uterus-like domain, where the motion occurs and the search for the small egg target.
We recall that the adult female uterus (Gray et al. 1974) is pear-shaped and is about 7.5 cm long, has a maximum diameter of 5 cm and a height of 3.4 cm, with a mean volume of tens of cm 3 . It is a hollow thick-walled non-convex muscular organ. On its upper part, the uterine tubes open, one on either side, while below, its cavity connects to that of the vagina. After an egg is released from the ovaries and moves inside the uterine cavity through the uterine tubes, it waits for fertilization (see Fig. 1). In summary, fertilization occurs most likely between the junction at the end of the uterus and somewhere in the fallopian tube, but not inside the uterus. We thus define the position of the target for this search process as the entrance of the fallopian tube, modeled here as a small gap between straight line and quarter-ellipses (see Fig. 1).
In estimating the expected search time, we model the spermatozoa motion as rectilinear far away from the egg, with random reflection when a trajectory hits the boundary of the uterus. Indeed, we adopt this model here following recent in vivo movies that we observed (courtesy of Reynaud K, private communication), where spermatozoa motions were mostly ballistic for distances of millimeters, as long as they do not encounter any obstacle. After hitting a small obstacle, the new direction was much correlated to the previous one, which is a significant deviation from a classical random walk. However, after hitting a wall and for trajectories that immediately returned to the bulk, the reflected angle was not correlated with the initial one. Other reflection behaviors include motion along the boundary, as observed in microfluidic chambers (Denissenko et al. 2012). It is however unclear at this time what exactly is the spermatozoa motion in vivo at long distances of centimeter sizes. We further construct numerical simulations of the search to explore a part of the phase-space, especially for non-convex shape geometry, where we have not yet found a fruitful analytical approach. Mathematically, finding a small target by this search process with deter-ministic reflection is analogous to the escape from a billiard table through one and two small holes (Bunimovich and Dettmann 2005).
This manuscript is organized as follows. We start with describing a crude rectilinear model for the swimmer motion and use it to estimate the time to find a small target (one of the two small fallopian entrances, which are two well-separated small targets, see Fig. 1). The swimmer motion is modeled as that of a point moving at constant speed in a fixed direction. It is reflected at the boundary in a random direction and absorbed when it hits the small target (of size ε). We estimate analytically the arrival time to the small target for different situations: starting with the unit disk, followed by the three-dimensional ball and then for more general two-dimensional convex domains (that are small deviation of a disk). We explore non-convex geometries by numerical simulations in dimension two. For non-convex domains, there are regions, including the one where the target is located, which are not accessible to direct rays, leading to a difficulty in computing the probability to find the target in finite time. The present approach provides a framework for the study of various parameters involved in the search process and reveals the role of the geometry in defining the search time.

Model simplifications
Spermatozoa motion results from the flagella beating attached to the round head. The ensemble head plus flagella allows the spermatozoa to swim in a complex medium, inside or on the surface of a domain. At a micrometer scale, the details of the spermatozoa geometry are relevant to the prediction of its motion, to account for the variety of displacements and of its orientation. For example, spermatozoa are reflected along a corner or at a surface with different properties, with random angles that are different from the mechanical and optical rules of the Snell-Descartes law of reflection (Kantsler et al. 2013;Ishimoto and Gaffney 2014).
In the present study, we focus on long-scale motion on the order of centimeters. We neglect the effect of any chemotaxis gradient that could direct the cell toward a specific target. Thus the motion of a spermatozoon is simplified to that of a point, which away from any boundary, is ballistic or rectilinear with a constant velocity (Prez-Cerezales et al. 2015). We postulate that when a spermatozoon hits a surface it is reflected at a random angle, uniformly distributed in the tangent half-space. Thus the simplified model of spermatozoon motion in a bounded domain is assumed rectilinear with constant velocityẊ where u is a fixed vector chosen on the unit sphere from a uniform distribution. Upon hitting an obstacle at a boundary point X 0 , the velocity changes tȯ where v is chosen on the unit sphere in the supporting half space at X 0 from a uniform distribution, independently of u. This model captures the crude sperm dynamics, neglecting the additional motion induced by the flagella and will be implemented in the simulations performed here.

Expected search time in a disk
When is a planar circular disk of radius R, the velocity is v 0 , and the target is an arc of length ε on the boundary, the probability p(ε) to hit the target in one step, starting the search at any point on the circle, is The probability to hit the target in exactly k steps is a geometric distribution q k = (1 − p(ε)) k−1 p(ε). Thus the expected number of steps to the target is To determine the expected search time, we assume that the target is centered at x = 0, y = R and that k reflections occur at the uniformly distributed angles (θ 1 , . . . θ k ) ( Fig. 2), where the angle of reflection is θ = P QS. Each ray is represented by its projection Re iθ k on the circle of radius R. The distance between the reflection points is The expectation of the right hand side is given by The time spent on a single ray is τ k = d k /v 0 , therefore the mean time to exit in N steps is Search for a small target in a disk and a ball. a Schematic representation of the dynamics in a disk: the motion is reflected on the boundary at a random angle. The target size ε is on the boundary. b A path with many random reflection escapes ultimately at the target (green dot). c Probability density function of the arrival time obtained from stochastic simulations. d The expected search time is the MFPT of the trajectory to the target, obtained from stochastic simulations (blue circles) and compared to the analytical (9) (red). e Search in a ball. f Comparison of the analytical curve (red) with the stochastic simulations (blue circles). Non-dimensional parameters: R = 1, v 0 = 1 (color figure online) Thus the mean first passage time to the target is where S = π R 2 . Formula (9) shows that the present search process is asymptotically of the order 1 ε , as ε goes to zero, much longer than the one for a free Brownian particle with diffusion coefficient D, searching for a target of size ε, in a two-dimensional domain of surface S, for which the Narrow escape theory (Cheviakov et al. 2012;Schuss et al. 2007) gives The results of numerical simulations of a search for a small arc (Fig. 2a) in a circular disk are shown in Fig. 2b, where the search starts at a point uniformly distributed on the boundary (red point). The tail of the exit time density decays exponentially (Fig. 2c). The increase at short times corresponds to a fraction of trajectories that shoot-up directly to the target. Finally, the numerical approximation for the expected search time is well matched by the asymptotic formula (9), where The distribution of the number of search steps prior to hitting the target decays exponentially, because Thus, the probability density of the search time τ e is The exit time event {τ e = t} can always be decomposed as τ e = kτ m + t * , where the residual time t * is distributed as where τ m = 2R/v 0 (the maximum search time) and the normalization constant is α = 1 τm 0 Pr{t * =s}ds . Indeed the velocity is constant and the initial point is uniformly distributed on the circle. Thus, as the number of steps before exit decays exponentially, we have that Pr{τ e = kτ m + t * = t, t * = s} ∼ Pr{τ e = kτ m = t − s} ∼ Pr{τ e = kτ m = t} for t 1. Thus, where

Expected search time in the three-dimensional ball
To compute the expected search time of a swimmer in a ball, we assume that the target is a small disk centered at the north pole of the sphere and that the dynamics are (1) and (2). The apex of the target is an angle β = 2ε on the sphere, the solid angle at any point on the sphere at position φ is in general a function of both variables φ and ε, denoted (φ, ε) (see Fig. 2), θ, φ are spherical coordinates on the sphere. At each step of hitting the boundary, the probability for the trajectory to find the target is p(θ ) = (φ, ε)/2π . The solid angle (φ, ε) depends on the position and the size of the hole (Paxton 1959) and the exact expression for formula (62) is given in Appendix A.2. However, in contrast to the case of a disk, the expected number of steps to hit the target, or the search time, depends on the history, due to the dependence of the probability on the variable φ. Indeed, the probability to hit the target in k steps, with the history angles φ 1 , .., φ k , is Thus the expected number of hits prior to escape is where the expectation is conditioned on all trajectories that are reflected with angles Thus, because all angles variables are i.i.d. random variables where the mean solid angle is defined by The solid angle at which the target is seen from the sphere center is 0 = 2π(1−cos ε).
To calculate the expected search time, we first compute the mean Euclidean distance along a ray between two consecutive reflecting points P k−1 = (φ k−1 , θ k−1 ) and P k = (φ k , θ k ) on the sphere, where (φ, θ ) are the spherical angles uniformly distributed in Thus, The integral (22) is evaluated numerically (Mathematica) to give an approximation for the mean length on a sphere It follows that the expected hitting time after exactly N hits is where v 0 is the moving speed of the swimmer, which is set to a constant value here.
Using (19), we finally get the general formula for the expected search time Hence using expression (21), it follows the estimate ). An asymptotic formula for an elliptic target can be derived from formula 62 (Appendix A.2) involving elliptic integrals.

Expected search time in a two-dimensional convex domain
Next, we consider the search for a small target in a planar starlike domain . We calculate the probability of a swimmer's trajectory to hit a small target on the boundary ∂ of the domain, which is parameterized in polar coordinates by r (θ ) so that in the complex plane, Q(θ ) = r (θ )e iθ .

The probability to hit a small window
When a trajectory starts at Q 0 ∈ ∂ it ends at another point Q 1 ∈ ∂ (see Fig. 3), the probability that the angular deviation of the particle moving path Q 0 Q 1 (assuming the path is a straight line) will lie between angle ψ and ψ + dψ is P(d Q 1 |Q 0 ) = dψ/π, in which ψ is the angle between the moving path and the tangent line at point Q 0 . Based on triangle geometry in the plane, we obtain π/2 + β = ψ + α, in which β is the angle between the line Q 0 Q 1 and the vertical direction (see Fig. 3) and α is the angle of the tangent line at Q 0 , then it turns out that dψ = dβ after applying differentiation to the equation (α is constant). The line Q 0 Q 1 , can be written in polar coordinate.
r cos(θ − β) = constant, in which (r, θ) are polar coordinates in the plane. Therefore the coordinates of the points Q 0 = (r 0 , θ 0 ) and Q 1 = (r 1 , θ 1 ) satisfy the equation where r 0 = r (θ 0 ) and r 1 = r (θ 1 ). Expanding (26) and rearranging, we obtain Differentiating (27), we get Therefore, the probability density of hitting the target in one step, P(d Q 1 |Q 0 ), can be written as In contrast to the case of a circle, here the density depends on the initial position.

The probability of hitting the target after n steps
The probability for a swimmer starting at Q 0 ∈ ∂ to hit a small arc ∂ a of length |∂ a | = ε, centered at Q 1 ∈ ∂ a , is The probability to hit the target after bouncing n times in the boundary There is no general simplifying form for this multiple integral. We consider, therefore, a perturbation of a circle.

Perturbation of a circle
A perturbation of the unit circle can be written as which for η 1 is a small amplitude deviation, as long as ρ is a function that keeps the domain convex. Substituting (32) in (29) and expanding to leading order in powers of η we obtain which allows to approximate the exit probability after n th iterations (relation (31)) as For η 1 and n ≥ 2 and whereP This can be further simplified tõ where The operators K 1 ( f ρ ) and K 2 ( f ρ , ∂ a ) are well defined, because

Number of hits before hitting the target
The expected number of steps E[N hit ] before hitting the target is expressed in terms of the probability P n (Q 0 ) to hit the target in n-steps, given the initial point Q 0 , The expansion (34) for small ε and η gives to leading order We decompose with to obtain for η ε, It follows that the expected number of steps prior to hitting the target depends on the integral K 1 ( f ρ ), which accounts for certain global properties of the boundary.
To compute the mean time before escape, we estimate the distance between two points P(θ 1 ) and P(θ 2 ) using the perturbation condition (32) Thus, Because the time spent on a single ray is τ k = d k /v 0 , the expected search time in N steps is Finally, (38) gives the expected search time as This expansion is valid for η ε. It shows that the first perturbation term depends on two integrals involving a generalized curvature function over the boundary. Although we shall not use directly these estimations in the numerical analysis below, these formulas reveal the geometrical features contributing to the search time. It would certainly be interesting to obtain a full expansion to order ε n .

Search for a small target in a uterus-like geometry: a simulation approach
Consider the search problem for a small target in a two-dimensional non-convex domain with cusps ( Fig. 4) and assume the dynamics (1)-(2). The two-dimensional domain in Fig. 4a captures some of the main features of the planar cross section of the uterine cavity (see discussion of the three-dimensional approximation below). The planar domain consists of two quarter-ellipses (black thick lines) underneath the upper part which is made of a straight line (black line on the top). The left and right sides are symmetric. The gaps between the straight line and the quarter-ellipses are the narrow openings.

Estimate of the expected the search time from numerical simulations
The simulation of equations (1)-(2) with random reflections at the boundary, except when a trajectory hits one of the two ending cusps, where it is terminated, gives the results shown in Fig. 4b-d (see Appendix A.1 for the details of the simulations and Table 1: we use the following non-dimensional parameters L = 5, W = 2.5, = 0.1, δ = 0.3 and velocity v = 1). A typical trajectory before reaching a small target is shown in Fig. 4a. The probability density function (pdf) Pr{τ ε = t} of the arrival time to either one of the cusp opening (size ε = 0.1) is shown in Fig. 4b. The tail of the distribution is well approximated by an exponential where the rate is λ = E[τ ε ] −1 . The results of numerical simulations give the following value for the mean time E[τ ε ] = 84.9. The exponential approximation for the tail seems valid for t > 100, but deviates slightly when 20 < t < 100. The short-time  distribution is not studied here for t < 20. We show the different pdfs as ε varies in Fig. 4c. After rescaling the pdfs by their mean values, they all align into one universal curve (Fig. 4c, inset), suggesting that a possible power scaling law for the pdf. To explore this possibility, we estimated the search time for various values of ε (Fig. 4d). Using a log-log plot, we compare the arrival time computed from stochastic simulations (blue dots) and the approximation by a power law (read line) E[τ ] ∼ Cε −α . An optimal fitting procedure shows that α = 0.658 ≈ 2/3. At this stage, we obtain the following approximation The power law is an approximation for small ε, but deviates as ε becomes of order O(1) (ε > 0.2 in Fig. 4e).
To convert in dimensional units, we use that a spermatozoa mean velocity is v = 75 μms −1 (Reynaud K, private communication) and that for the length in Fig. 4 are L = 7.5 cm, W = 2.5 cm, ε = 0.03 cm, δ = 0.2 cm (where the non-dimensional velocity was v = 1) and the non-dimensional arrival time was E[τ ] ≈ 238 (Fig. 4). By scaling with the velocity, we obtain the dimensional search time E[τ DU ] = 238·130 = 30,940 s, which is about 8.6 h. This two-dimensional estimate provides a time scale for the arrival to one of the two targets. Indeed, an egg is positioned at one of these places, from which we obtain the expected search time as E[τ target ] ≈ 17 h. This time estimation is compatible to previous experimental report (Chang 1951). Finally, the shortest time, which can be obtained by the concatenation the two straight lines (see discussion below) gives an estimate of 1333 s, which is about 22 min.

Effect of changing some geometrical parameters on the search time
We now explore the consequences of changing various parameters on the arrival time, this includes looking at the geometry of the uterus-like domain. We vary here several parameters such as the target size ε, the length L and the width W and the entrance size δ, which represents the cervix outside radius. First, by varying the size δ (in range [0.1 − 0.8]), while keeping the other parameters constant (Fig. 5a), we obtain that the pdfs of arrival time show only a slight difference (less than 20 %). Thus, the numerical simulations show that the cervix size δ has little influence on the search time.
We next varied (Fig. 5b), the aspect ratio as indicated in Fig. 5 (the length of the quarter of an ellipse L and W are defined in Fig. 4). For large aspect ratio ρ, the domain has a long narrow shape, while it is flat and wide for small aspect ratio. The pdf of the search time for different aspect ratios, while keeping the circumference of the quarter-ellipses fixed, is shown in Fig. 5b. Note that curves (normalized by the maximum) obtained for different parameters, align to a single curve (see inset of Fig. 5b), suggesting that the pdf has a scaling law structure. Actually, all normalized pdfs for the search time with different opening size δ and aspect ratios ρ follow a universal curve (Fig. 5c), indicating that a general scaling law is expected as a function of the two parameters. Finally, we estimated numerically, the mean search time as a function of the aspect ratio (Fig. 5d): long narrow or flat geometries are associated with short mean search time, while a unit aspect ratio (a quarter-circle) has the maximum search time, four times larger. There is no theory today to explain this maximum.

Discussion and conclusion
In this study, we presented the search process in a reduced domain presenting some aspects of uterus geometry, such as opening sizes of oviduct and cervix. We varied some parameters such as the aspect ratio ρ = L W defining the shape of the uterus and estimated the search time of a spermatozoon to the egg target site, in a cusp that represents the entrance of the oviduct.
In the absence of guide mechanisms, random reflection is the key element that determines the search time. We computed the expected search time in two-and threedimensional balls and domains that resemble the uterus geometry. We found that the search time depends on the geometry of the domain. In two dimensions, we found for ε 1 that where K is a constant, which we determined for a sphere and some convex domains. When the domain contains ends of cusp geometry, where K C is now a dimensional constant. We could only obtain the exponent α = 0.66 numerically. The exact relation with the geometry remains an open question. The formula for a domain of arbitrary shape is unknown. Although the arrival time decreases with the oviduct opening size, the opening size of cervix (parameter δ) had little effect. According to our numerical simulations, we found here that long narrow shape (high aspect ratio ρ) leads to arrival time much shorter than for a round uterus shape. These results should be considered as predictions because there are little available data on this subject. The search time for the rectilinear motion studied here, describing spermatozoa motion, is much longer than for a diffusion particle, which is O(log 1 ε ) in dimension 2 and O( 1 ε ) for a cusp . We conclude that spermatozoa stay much longer in the domain before they can find the target, and when there are many swimmers, this long period can be seen as a selection process based on intrinsic spermatozoa properties, to select the fittest one. In three dimensions, we found that E[τ ε ] = K V v 0 ε 2 in a ball of volume V . We conjecture that this formula is valid for convex domains. However, for a domain with cusps, the scaling law may depend on the local cusp property. Another surprising result is that the power law for the search time in a cusp domain is smaller than 1, thus the search time is faster in a domain with a cusp compared to the disk. This situation is exactly the opposite for Brownian motion, where the search time becomes exponentially longer in a cusp (Schuss et al. 2007).
The present modeling approach allows for the first time to obtain laws for the expected search time. Most of the previously published numbers were obtained from empirical data, where spermatozoa were counted after opening the uterus at different stages (Chang 1951). Modeling is a different approach that accounts for the local motion and for the entire global geometrical structure of the uterus. Our approach should certainly be complemented by in vivo endoscopy studies. But there are still some difficulties in reconstructing the search process from the local data generated by the endoscopy probe that can only access local domains of few microns.
Another consequence of the present analysis is that the number of swimmers matters when we consider the first one that finds a neighborhood of the egg. By increasing the number of trajectories (Fig. 6), when computing the arrival time of the first one, the trajectories that are almost optimizing the search time are located in a neighborhood around the optimal ones (white dashed lines in Fig. 6). The optimal trajectory is indeed the geodesic that joins the initial point to the final egg location. The shortest search time is achieved by the geodesic that minimizes the functional where d(s) is the Euclidean distance at time s and as indicated by the results of our simulations (Fig. 6). Possible future directions include the exploration of other parameters of the search process, such as the fitness of spermatozoa, the time window when they can sense the presence of the egg and chemotaxis processes, and so on. These elements should certainly be included in future models, as well as a possible killing field, which should account for sperm degradation during their sojourn time.
In summary, isometric shapes are not optimal for the search process and thus for fertilization. The uterus of most animals including human and rabbit, have long narrow shape, optimal for the search time and possibly fertilization, however, the optimization search process may also depend on other factors, such as body size, dynamics, surface of the uterus, uterus contraction and many more. The aspect ratio of the human uterus is around 2 to 3, which does not give the optimal ratio according to our simulations (Fig. 5) and we speculate that it might be associated with a lower fertility rate compared to other species.
where trajectories are reflected with a random angle with respect to the tangent. In case of the uterus-like domain, there are two small targets located at the end of the two symmetrical cusps. The target can only be found by trajectories moving inside de the domain.
The escape time of a trajectory from the domain is computed from the total length before exiting, as the velocity is constant. For each case, the statistic is obtained for 30,000 runs. We did not see any differences by using additional simulations such as 50,000.

A.2 Probability that a trajectory enter the solid angle of the target in the three dimensional ball
In this section, we compute the solid angle (φ, ε) starting from any point (φ, θ ) located on a sphere. The solid angle subtended by a circular disk (Paxton 1959) is summarized here.
By definition, the solid angle is given by the integral = n P · ds P R 2 , where n P · ds P is the area of the projection of the surface element ds P onto the plane perpendicular to the radius R = |N P|, from the north pole to point P (see Fig. 7). Specifically, Integral 55 can be written as: = ds cos θ R 2 = ρdβdρ cos θ R 2 = sin θ dθ dβ (56) Fig. 7 Solid angle subtended at points a over the interior or over the periphery of disk (r 0 ≤ r m ); b outside disk boundary (r 0 > r m ) (these figures are screenshots from Poxtan's paper; Paxton 1959) in which ρ = L tan θ and L/R = cos θ . When r 0 ≤ r m (see Fig. 7a), Eq. 56 can be written as 2 = π 0 θ s 0 sin θ dθ dβ = π − π 0 cos θ s dβ (57) After a few more steps and re-arrangements, Eq. 57 turns into where K (k) and 0 (ξ, k) are the complete elliptic integrals of the first and third kind, respectively with k 2 = 1 − (R 1 /R max ) 2 R 2 1 = L 2 + (r 0 − r m ) 2 (59) ξ = sin −1 (L/R 1 ).
When r 0 > r m (see Fig. 7b Similarly, after certain steps and re-arrangements, Eq. 60 can be written as follows: In summary, the expression for the solid angle is (ξ, k) is a function of ε and φ as we shall see now. We apply formula 62 to the case of a small hole of radius ε located at the north pole of the three-dimensional sphere. We obtain for a point of coordinates (θ, φ) We then used Mathematica to compute the φ>ε/R (φ) sin(φ)dφ. We obtain for small ε that φ>ε/R (φ) sin(φ)dφ = π (1 − cos ε) .