text:
# Strong gravitational lensing: Why no central black holes?

## 1 Introduction

Cold Dark Matter (CDM) has become the standard theory of cosmological structure formation. The $`\mathrm{\Lambda }`$CDM variant of CDM with $`\mathrm{\Omega }_m=1\mathrm{\Omega }_\mathrm{\Lambda }0.3`$ appears to be in good agreement with the available data on large scales (Primack primack (2002)). On smaller (sub-galactic) scales, there seem to be various discrepancies, such: N-body CDM simulations which give cuspy halos with divergent profiles towards the center (Navarro, Frenk and White nfw96 (1996), nfw97 (1997), NFW hereafter); bar stability in high surface brightness spiral galaxies which also demands low-density cores; CDM models which yield an excess of small scale structures; formation of disk galaxy angular momentum, which is much too small in galaxy simulations. Issues that have arisen on smaller scales have prompted people to propose a wide variety of alternatives to CDM, such as warm dark matter (WDM) and self-interacting dark matter (SIDM). Now that problems arise from galaxy-size halos and centers of all dark matter halos, high-resolution simulations and observations are the final criterion. Recent highest-resolution simulations appear to be consistent with NFW (Klypin, klypin (2002); Power et al. power (2002)) until scales smaller than about 1 kpc. Meanwhile, a large set of high-resolution optical rotation curves has recently been analyzed for low surface brightness (LSB) galaxies. One can also conclude that the NFW profile is a good fit down to about 1 kpc. Although further simulations and observations, including measurement of CO rotation curves (Bolatto et al., bolat (2002)), may help to clarify the nature of the dark matter, it now appears that WDM and SIDM are both probably ruled out, while the small-scale predictions of $`\mathrm{\Lambda }`$CDM may be in better agreement with the latest data than appeared to be the case as recently as a year ago.

In addition to direct simulations and observations, gravitational lensing provides another powerful probe of mass distribution in the universe. Since mass within small scales only deflect light rays slightly, it is difficult to extract mass information from a single lensing event, and thus statistical gravitational lensing is needed even for “strong” gravitational lensing of small halos(Turner, Ostriker & Gott turner (1984); Narayan & White narayan (1988); Cen et al. cen (1994); Kochanek kochanek (1995); Wambsganss et al. wambs95 (1995); Wambsganss, Cen, & Ostriker wambs98 (1998); Porciani & Madau porci (2000); Keeton & Madau keeton (2001)). Li & Ostriker (li (2002)) first used the semi-analytical approach to analyze gravitational lensing of remote quasars by foreground dark halos in various cold dark matter cosmologies. The mass function of dark halos they used is alternatively given by singular isothermal sphere (SIS), the NFW profile, or the generalized NFW profile. They found that none of these models can completely explain the current observations: the SIS models predict too many large splitting lenses, while the NFW models predict too few small splitting lenses, so they proposed that there must be at least two populations of halos in the universe: small mass halos with a steep inner density slope and large mass halos with a shallow inner density slope. The author conclude that a combination of SIS and NFW halos can reasonably reproduce the current observations. Similarly, Sarbu et al. (sarbu (2001)) investigated the statistics of gravitational lenses in flat, low-density cosmological models with different cosmic equations of state $`\omega `$. It was found that COBE-normalized models with $`\omega >0.4`$ produce too few arcsecond-scale lenses in comparison with the JVAS/CLASS radio survey, a result that is consistent with other observational constraints on $`\omega `$.

When attention is attracted to alternatives of CDM dark matter density profile at small scales, another kind of dark matter — super-massive black holes in the centers of most galactic halos is forgotten or ignored in this case, although the idea of detecting supermassive compact objects by their gravitational lensing effects was proposed very early (Press & Gunn press73 (1973), Wilkinson et al. wilki (2001)) and the lensing effects of Schwarzschild black holes in the strong field regime have been discussed in detail (e.g., Virbhadra & Ellis virbh (2000); Frittelli, Kling & Newman fritt (2000); Bozza et al. bozza (2001)). On the other hand, cosmological voids can form directly after the collapse of extremely large wavelength perturbations into low-density black holes or cosmological black holes; such black holes can also be detected through their weak and strong lensing effects (Stornaiolo storn (2001)). The observational evidence presented so far suggests the ubiquity of BHs in the nuclei of all bright galaxies, regardless of their activity, and BH masses correlate with masses and luminosities of the host spheroids and, more tightly, with stellar velocity dispersions (Magorrian et al. magor (1998); Ferrarese & Merritt ferra (2000); Ravindranath et al. ravin (2001); Merritt & Ferrarese 2001a , 2001b ; Wandel wandel (2002); Sarzi et al. sarzi (2002)). Most recent high-resolution observational data gives $`M_{}/M_{\mathrm{bulge}}10^3`$(Merritt & Ferrarese 2001c ). Ferrarese (ferra02 (2002)) further gave the relation between masses $`M_{}`$ of supermassive black holes and the total gravitational mass of the dark matter halo in which they presumably formed

$$\frac{M_{}}{10^8M_{\mathrm{}}}0.046\left(\frac{M_{\mathrm{DM}}}{10^{12}M_{\mathrm{}}}\right)^{1.57}.$$
(1)
In this paper, we investigate the contributions of galactic central black holes to lensing probabilities at small image separations. Since $`\mathrm{\Lambda }`$CDM cosmology and NFW profile are in good agreement with the available data of structure formation on almost all scales as mentioned above, we only chose these two models respectively as cosmology and mass density function in our calculations. We model the lenses as a population of dark matter halos with an improved version of the Press-Schechter (press74 (1974), PS) mass distribution function, and central BHs are considered for galaxy-size halos.

The paper is organized as follows: the lensing equation is given in Sect. 2, lensing probabilities are calculated in Sect. 3, and discussion and conclusions are provided in Sect. 4.

## 2 Lensing equation

The NFW profile is

$$\rho _{\mathrm{NFW}}=\frac{\rho _\mathrm{s}r_\mathrm{s}^3}{r(r+r_\mathrm{s})^2}$$
(2)
where $`\rho _\mathrm{s}`$ and $`r_\mathrm{s}`$ are constants. We can define the mass of a halo to be the mass within $`r_{200}`$ (which is the radius of a sphere around a dark halo within which the average mass density is $`200`$ times the critical mean mass density of the universe),

$$M_{\mathrm{DM}}=4\pi _0^{r_{200}}\rho r^2𝑑r=4\pi \rho _\mathrm{s}r_\mathrm{s}^3f(c_1),$$
(3)
with $`c_1=r_{200}/r_\mathrm{s}`$ the concentration parameter, the value of which is chosen to be $`7/(1+z)`$ (Bartelmann el al. barte (1998)). And

$$f(c_1)=_0^{c_1}\frac{x^2dx}{x(1+x)^2}=\mathrm{ln}(1+c_1)\frac{c_1}{1+c_1}.$$
(4)
In flat $`\mathrm{\Lambda }`$CDM cosmology, the constants $`\rho _\mathrm{s}`$ and $`r_\mathrm{s}`$ can then be expressed as (Li & Ostriker li (2002)),

$$\rho _\mathrm{s}=\rho _{\mathrm{crit}}\left[\mathrm{\Omega }_\mathrm{m}(1+z)^3+\mathrm{\Omega }_\mathrm{\Lambda }\right]\frac{200}{3}\frac{c_1^3}{f(c_1)},$$
(5)
$$r_\mathrm{s}=\frac{1.626}{c_1}\frac{M_{15}^{1/3}}{\left[\mathrm{\Omega }_\mathrm{m}(1+z)^3+\mathrm{\Omega }_\mathrm{\Lambda }\right]^{1/3}}h^1\mathrm{Mpc}.$$
(6)
where $`\rho _{crit}`$ is the present value of the critical mass density of the universe, and $`M_{15}`$ is the reduced mass of a halo defined as $`M_{15}=M_{\mathrm{DM}}/(10^{15}\mathrm{h}^1M_{\mathrm{}})`$.

The surface mass density for NFW profile is

$`\mathrm{\Sigma }_{\mathrm{NFW}}(x)`$ $`=`$ $`2\rho _\mathrm{s}r_\mathrm{s}\times \{\begin{array}{cc}\frac{\mathrm{ln}x\sqrt{1x^2}\mathrm{ln}(1\sqrt{1x^2})}{(1x^2)^{3/2}},\hfill & (x>1),\hfill \\ \frac{1}{3},\hfill & (x=1),\hfill \\ \frac{\mathrm{arcsin}(1/x)+\sqrt{x^21}\pi /2}{(x^21)^{3/2}},\hfill & (0<x<1).\hfill \end{array}.`$ (10)
Where $`x=|x|`$ and $`x=\xi /r_\mathrm{s}`$, $`\xi `$ is the position vector in the lens plane. The galactic central black holes are assumed to be point masses, and we consider fist there is only a single black hole with mass $`M_{}`$ for each galaxy. So the surface mass density for galactic halos each with a single central black hole can be written as

$$\mathrm{\Sigma }_{\mathrm{galaxy}}(x)=M_{}\delta ^2(x)+\mathrm{\Sigma }_{\mathrm{NFW}}(x),$$
(11)
where $`\delta ^2(x)`$ is the two dimensional Dirac-delta function. The lensing equation with galactic central black holes considered then is

$$y=x\mu _\mathrm{s}\frac{f_{\mathrm{BH}}+g(x)}{x},$$
(12)
where $`y=|y|`$, $`\eta =yD_\mathrm{S}^\mathrm{A}/D_\mathrm{L}^\mathrm{A}`$ is the position vector in the source plance, in which $`D_\mathrm{S}^\mathrm{A}`$ and $`D_\mathrm{L}^\mathrm{A}`$ are angular-diameter distances from the observer to the source and to the lens respectively. And

$$\mu _\mathrm{s}=\frac{4\rho _\mathrm{s}r_\mathrm{s}}{\mathrm{\Sigma }_{\mathrm{cr}}},$$
(13)
where $`\mathrm{\Sigma }_{\mathrm{cr}}=(c^2/4\pi G)(D_\mathrm{S}^\mathrm{A}/D_\mathrm{L}^\mathrm{A}D_{\mathrm{LS}}^\mathrm{A})`$ is the so called critical surface mass density, in which $`D_{\mathrm{LS}}^\mathrm{A}`$ is the angular-diameter distance from the lens to the source. And

$`g(x)`$ $`=`$ $`\mathrm{ln}{\displaystyle \frac{x}{2}}+\{\begin{array}{cc}\frac{\mathrm{arctan}\sqrt{x^21}}{\sqrt{x^21}}\hfill & (x>1),\hfill \\ 1\hfill & (x=1),\hfill \\ \frac{\mathrm{arctanh}\sqrt{1x^2}}{\sqrt{1x^2}}\hfill & (0<x<1).\hfill \end{array}`$ (17)
In Eq.(12), the term $`f_{\mathrm{BH}}`$ stands for the contribution of a black hole, and by using Eq.(1) and Eq.(3), it has the form

$$f_{\mathrm{BH}}=2.78\times 10^4f(c_1)M_{15}^{0.57}.$$
(18)
Since there is always more than one black holes in a bulge, and thus the bulge itself can act like a black hole, we can treat a bulge as a point mass in this paper as an extreme case. The mass of a black hole correlates linearly with that of its host bulge as $`M_{}/M_{\mathrm{bulge}}10^3`$, so we can simply multiply the term $`f_{\mathrm{BH}}`$ by $`10^3`$ to stand for the contribution of a bulge. However, some light rays from the source will definitely travel across the bulge, so there must exist a kind of “effective” black hole with mass larger than a single “real” black hole but less than the bulge. In order to investigate the tendency of image separations contributed by different point mass, we can multiply $`f_{\mathrm{BH}}`$ by, for example, $`10^2`$, etc.

The lensing equations for three cases are plotted in Fig. 1 according to Eq.(12), where we have extended $`x`$ and $`y`$ to their opposite values because of symmetry. The full line, dashed line and dotted line, respectively, represent the NFW lens with $`f_{\mathrm{BH}}=0`$, ‘NFW+BH’ lens with $`f_{\mathrm{BH}}=2.78\times 10^4f(c_1)M_{15}^{0.57}`$ and ‘NFW+bulge’ with $`f_{\mathrm{BH}}=2.78\times 10^2f(c_1)M_{15}^{0.57}`$, in which $`\mu _\mathrm{s}=0.49`$ and $`f(c_1)=0.91`$.

We find that, as point masses, both a single central black hole and a bulge can more often produce small separation images than the case when no central black holes are considered. This result will be further confirmed by the lensing probability given in next section.

## 3 Lensing probability

We choose the most generally accepted values of the parameters for flat $`\mathrm{\Lambda }`$CDM cosmology, for which, with usual symbols, the matter density parameter, vacuum energy density parameter and Hubble constant are respectively: $`\mathrm{\Omega }_\mathrm{m}=0.3`$, $`\mathrm{\Omega }_\mathrm{\Lambda }=0.7`$, $`h=0.75`$.

The quasars of redshift $`z_\mathrm{s}=1.5`$ are lensed by foreground CDM halos of galaxy clusters and galaxies, the lensing probability with image separations larger than $`\mathrm{\Delta }\theta `$ is (Schneider, Ehlers, & Falco schne (1992))

$$P(>\mathrm{\Delta }\theta )=_0^{z_\mathrm{s}}\frac{dD_\mathrm{L}(z)}{dz}𝑑z_0^{\mathrm{}}\overline{n}(M,z)\sigma (M,z)𝑑M.$$
(19)
Where $`D_\mathrm{L}(z)`$ is proper distance from the observer to the lens located at redshift $`z`$

$$D_\mathrm{L}(z)=\frac{c}{H_0}_0^z\frac{dz}{(1+z)\sqrt{\mathrm{\Omega }_\mathrm{m}(1+z)^3+\mathrm{\Omega }_\mathrm{\Lambda }}},$$
(20)
here $`c`$ is the speed of light in vacuum and $`H_0`$ is the current Hubble constant. The physical number density $`\overline{n}(M,z)`$ of virialized dark halos of masses between $`M`$ and $`M+dM`$ is related to the comoving number density $`n(M,z)`$ by $`\overline{n}(M,z)=n(M,z)(1+z)^3`$; the latter was originally given by Press & Schechter (press74 (1974)), and the improved version is

$$n(M,z)dM=\frac{\rho _0}{M}f(M,z)dM,$$
(21)
where $`\rho _0`$ is the current mean mass density of the universe, and

$$f(M,z)=\sqrt{\frac{2}{\pi }}\frac{\delta _c(z)}{M\mathrm{\Delta }}\frac{d\mathrm{ln}\mathrm{\Delta }}{d\mathrm{ln}M}\mathrm{exp}\left[\frac{\delta _c^2(z)}{2\mathrm{\Delta }^2}\right]$$
(22)
is PS mass function. In Eq.(22) above, $`\mathrm{\Delta }^2(M)`$ is the present variance of the fluctuations in a sphere containing a mean mass $`M`$,

$$\mathrm{\Delta }^2(M)=\frac{1}{2\pi ^2}_0^{\mathrm{}}P(k)W^2(kr_\mathrm{M})k^2𝑑k,$$
(23)
where $`P(k)`$ is the power spectrum of density fluctuations, $`W(kr_\mathrm{M})`$ is the Fourier transformation of a top-hat window function

$$W(kr_\mathrm{M})=3\left[\frac{\mathrm{sin}(kr_\mathrm{M})}{(kr_\mathrm{M})^3}\frac{\mathrm{cos}(kr_\mathrm{M})}{(kr_\mathrm{M})^2}\right],$$
(24)
and

$$r_\mathrm{M}=\left(\frac{3M}{4\pi \rho _0}\right)^{1/3}.$$
(25)
In Eq.(22), $`\delta _c(z)`$ is the over density threshold for spherical collapse by redshift $`z`$ (Navarro, Frenk, & White nfw97 (1997)):

$$\delta _c(z)=\frac{1.68}{D(z)},$$
(26)
where $`D(z)`$ is the linear growth function of density perturbation (Carroll, & Press carroll (1992))

$$D(z)=\frac{g(\mathrm{\Omega }(z))}{g(\mathrm{\Omega }_\mathrm{m})(1+z)},$$
(27)
in which

$$g(x)=\frac{5}{2}x\left(\frac{1}{70}+\frac{209x}{140}\frac{x^2}{140}+x^{4/7}\right)^1,$$
(28)
and

$$\mathrm{\Omega }(z)=\frac{\mathrm{\Omega }_\mathrm{m}(1+z)^3}{1\mathrm{\Omega }_\mathrm{m}+\mathrm{\Omega }_\mathrm{m}(1+z)^3}.$$
(29)
We use the fitting formulae for CDM power spectrum $`P(k)`$ given by Eisenstein & Hu (eisen (1999))

$$P(k)=AkT^2(k),$$
(30)
where $`A`$ is the amplitude normalized to $`\sigma _8=\mathrm{\Delta }(r_\mathrm{M}=8h^1\mathrm{Mpc})=0.95`$, and

$$T=\frac{L}{L+Cq_{\mathrm{eff}}^2},$$
(31)
with

$$L\mathrm{ln}(e+1.84q_{\mathrm{eff}}),$$
(32)
$$q_{\mathrm{eff}}\frac{k}{\mathrm{\Omega }_\mathrm{m}h^2\mathrm{Mpc}^1},$$
(33)
$$C14.4+\frac{325}{1+60.5q_{\mathrm{eff}}^{1.11}}.$$
(34)
We need to know the cross-sections in Eq.(19). Since we are interested in the lensing probabilities with image separations larger than a certain value $`\mathrm{\Delta }\theta `$ (ranging from $`010`$ arcseconds, for example), the cross-section is defined under the condition that the multiple images can be created. For the lenses with NFW profile, one can see from Fig.1 that multiple images can be produced only if $`|y|y_{\mathrm{cr}}`$, where $`y_{\mathrm{cr}}`$ is the maximum value of $`y`$ when $`x<0`$, which is determined by $`dy/dx=0`$ when $`f_{\mathrm{BH}}=0`$ in Eq.(12). For galaxy-size halos, the mass of which is confined to be less than $`10^{14}M_{\mathrm{}}`$ through out this paper, a central black hole or a bulge as a point mass is considered (see Eq.(12)). In this case, multiple images will always exist: when the source is close to the point caustic, i.e., when $`y`$ is small, there are three images, two of which are within the Einstein circle and the third one the outside which Einstein circle; when $`y`$ is large enough, there are two images, the weaker one is close to the center of the Einstein circle and the brighter one locates outside of the Einstein circle. So another condition is needed to define the cross-section, for which we use the brightness ratio between the brighter and weaker images just mentioned, and it is enough to set the ratio to be $`10`$.

The brightness ratio $`r`$ for the two images is just the ratio of the corresponding absolute values of magnifications (Schneider, Ehlers, & Falco schne (1992)),

$$r=\left|\frac{\mu _+}{\mu _{}}\right|,$$
(35)
where

$$\mu _+(y(x))=\left(\frac{y}{x}\frac{dy}{dx}\right)_{x>0},$$
(36)
$$\mu _{}(y(x))=\left(\frac{y}{x}\frac{dy}{dx}\right)_{x<0}.$$
(37)
Once the source position $`y_{\mathrm{cr}}`$ is determined by

$$|\mu _+(y_{\mathrm{cr}})|=10|\mu _{}(y_{\mathrm{cr}})|,$$
(38)
the cross-section can be calculated, both with and without central black holes, as

$$\sigma (M,z)=\pi y_{\mathrm{cr}}^2r_\mathrm{s}^2\vartheta (MM_{\mathrm{min}}),$$
(39)
where $`\vartheta (x)`$ is a step function, and $`M_{\mathrm{min}}`$ is determined by lower limit of image separation

$$\mathrm{\Delta }\theta =\frac{r_\mathrm{s}\mathrm{\Delta }x}{D_\mathrm{L}^\mathrm{A}}\frac{2x_0r_\mathrm{s}}{D_\mathrm{L}^\mathrm{A}}$$
(40)
and Eq.(6) as

$$M_{\mathrm{min}}=8.927\times 10^8M_{15}\left(\mathrm{\Omega }_\mathrm{m}(1+z)^3+\mathrm{\Omega }_\mathrm{\Lambda }\right)\left(\frac{c_1D_\mathrm{L}^\mathrm{A}\mathrm{\Delta }\theta }{x_0}\right)^3.$$
(41)
In Eq.(40), we have approximated the image separation $`\mathrm{\Delta }x`$ to be $`2x_0`$, where $`x_0`$ is the positive zero position of function $`y(x)`$ , both when $`f_{\mathrm{BH}}=0`$ (for NFW lens only) and $`f_{\mathrm{BH}}0`$ (for galactic, NFW+BH/bulge lenses) in Eq.(12), since image separation is insensitive to the source position $`y`$ (Li & Ostriker, li (2002)).

We plot the lensing probability with image separations larger than $`\mathrm{\Delta }\theta `$ in Fig. 2. In order to show the tendency of contributions to the lensing probability for different fraction of the bulge mass (the so called ‘effective’ black hole), we take the term $`f_{\mathrm{BH}}`$ in Eq.(12) to be $`f_{\mathrm{BH}}=2.78\times 10^1f(c_1)M_{15}^{0.57}`$, $`f_{\mathrm{BH}}=2.78\times (5.0\times 10^2)f(c_1)M_{15}^{0.57}`$ and $`f_{\mathrm{BH}}=2.78\times 10^2f(c_1)M_{15}^{0.57}`$; they are represented by the first three lines from top down, respectively. When central black holes or bulges are included and treated as point masses, the mass of their host halos is confined to be less than $`10^{14}M_{\mathrm{}}`$, because Eq.(1) strictly applies in the range $`10^6<M_{}<2\times 10^9M_{\mathrm{}}`$ and $`10^{14}M_{\mathrm{}}`$ is the upper-limit of galaxy mass (Ferrarese 2002).

## 4 Discussion and conclusions

Our numerical results for lensing probability with image separations larger than $`\mathrm{\Delta }\theta `$ in five different cases are shown in Fig. 2. In all cases, lensing probabilities keep nearly constant until $`\mathrm{\Delta }\theta 0.1`$ arc seconds, and obvious dropdown takes place at about $`1`$ arc second if central black holes are included, which, of course, does not mean that the main lensing events have image separations larger than $`0.1`$ arc seconds. As a matter of fact, in the NFW case (without galactic central black holes, the full line in Fig. 2), the lensing probability drops quite slowly in the whole range of image separations: $`\mathrm{\Delta }\theta `$ 0—10 arc seconds; such a tendency would extend even to $`30`$ arc seconds if it is plotted beyond this range, which implies a uniform distribution of lensing probability for its $`log`$ value among image separations. However, note that in the single black hole case (dotted line), the lensing probability drops to the same value of NFW at 2 arc seconds, which gives the influence range of a single black hole. In the range of $`00.1`$ arc seconds, the lensing probability for the single black hole case is about 3 times that for NFW. So, clearly, the contributions from central black holes cannot be omitted, although such contributions alone are indeed not enough to explain the observational data. As we have pointed out, there is always more than one black hole in a galactic bulge, and the collector of black holes would make a bulge itself ‘act like’ a black hole. On the other hand, not all the mass of a bulge is concentrated in black holes, so if we treat a whole bulge as an extreme black hole, such a model would produce too many lenses at image separations larger than 3 arc seconds compared with the JVAS/CLASS radio survey. As mentioned above, we have sufficient reason to tune the fraction of a bulge mass to produce a ‘right’ profile of lensing probability at larger image separations required by observational data, but this ‘sufficient reason’ seems not make us produce sufficient lensing probabilities at smaller image separations, as shown by the dash-dot line in Fig. 2.

However, we can attribute sufficient lens events at small image separations to galactic central black holes or the bulge. On the one hand, since this paper focuses on whether galactic central black holes would contribute considerably to the lensing probability, we have not considered the effect of magnification bias, which would increase the final result provided here at all image separations. On the other hand, we have used an improved version of the PS halo mass function but not the ‘best’ version. The shape of the mass function predicted by standard PS theory (the improved version) is in reasonable agreement with what is measured in numerical simulations of hierarchical clustering from Gaussian initial conditions only for massive halos,; less massive halos are more strongly clustered or less anti-biased than the standard PS predicted. Sheth & Tormen (st (1999), ST) proposed a model that provides a reasonably good fit to the bias relation of less massive haloes as well as to that of massive halos. Note that central black holes are only found in galactic bulges,; ST’s correction for mass function in the range of less massive halos would definitely change the lensing probability discussed in this paper. Also note that one of the two images produced by a galactic central black hole is close to the lens center and very faint; however, VLBI experiments can detect its existence (Hirabayashi hirab (1990); Ulvestad ulves (1999)), and further radio lensing surveys would have the ability to identify high flux density ratio of the two images. How and to what extent lensing magnification bias, flux ratio and modified PS mass function may change the final result will be discussed in another paper.

###### Acknowledgements.

I thank Professor Xiang-Ping Wu for his original idea for this paper and the anonymous referee for useful comments and suggestions. This work was supported by the National Natural Science Foundation of China.


format: text

text:
# PARALLEL DELAUNAY REFINEMENT: ALGORITHMS AND ANALYSES

## 1 Introduction

Delaunay refinement is a popular and practical technique for generating well-shaped unstructured meshes . The first step of a Delaunay refinement algorithm is the construction of a constrained or conforming Delaunay triangulation of the input domain. This initial Delaunay triangulation need not be well-shaped. Delaunay refinement then iteratively adds new points to the domain to improve the quality of the mesh and to make the mesh respect the boundary of the input domain. A sequential Delaunay refinement algorithm typically adds one new vertex per iteration, although sometimes one may prefer to insert more than one new vertex at each iteration. Each new point or set of points is chosen from a set of potential candidates — the circumcenters of poorly conditioned simplices (to improve mesh quality) and the diameter-centers of boundary simplices (to conform to the domain boundary). Ruppert  was the first to show that the proper application of Delaunay refinement produces well-shaped meshes in two dimensions whose size is within a small constant factor of the best possible. Ruppert’s result was then extended to three dimensions by Shewchuk  and Li and Teng . Efficient sequential Delaunay refinement software has been developed both in two  and three dimensions . Chrisochiedes and Nave  and Okusanya and Peraire  developed parallel software using Delaunay refinement, for which they have reported good performance. Recently, Nave et al.  presented a parallel Delaunay refinement algorithm and proved that it produces well-shaped meshes. The complexity of their algorithm as well as the size of the mesh it outputs remains unanalyzed.

In this paper, we study the parallel complexity of a natural parallelization of Delaunay refinement. One of the main ingredients of our parallel method is a notion of independence among potential candidates for Delaunay insertion at each iteration. Our parallel Delaunay method performs the following steps during each iteration.

1. Generate an independent set of points for parallel insertion;
2. Update the Delaunay triangulation in parallel.

Our independent sets have the following properties:

* Their insertion can be realized sequentially by Ruppert’s method in 2D and Shewchuk’s in 3D. Hence, an algorithm that inserts all their points in parallel will inherit the guarantees of Ruppert’s and Shewchuk’s methods that the output mesh be well-shaped and have size optimal up to a constant.
* The independent sets can be generated efficiently in parallel. In addition, they are “large enough” so that the number of parallel iterations needed is $`O(\mathrm{log}^2(L/s))`$, where $`L`$ and $`s`$ are the diameter of the domain and the smallest edge in the output mesh, respectively.
* When a quasi-uniform mesh is desired as in Chew’s method, the number of iterations can be reduced to $`O(\mathrm{log}(L/s))`$.

We should emphasize here that our analysis focuses on the number of parallel iterations of Delaunay refinement. The independence of the new points do not necessarily imply a straightforward parallel insertion scheme at each iteration. There are several existing parallel Delaunay triangulation algorithms that we can employ at each iteration. For example, in 2D we can use the divide-and-conquer parallel algorithm developed by Blelloch et al.  for Delaunay triangulation. Their algorithm uses $`O(n\mathrm{log}n)`$ work and $`O(\mathrm{log}^3n)`$ parallel time. We can alternatively use the randomized parallel algorithms of Reif and Sen , or by Amato et al. , in both two and three dimensions. Both of these randomized parallel Delaunay triangulation algorithms have expected parallel running time $`O(\mathrm{log}n)`$. Using one of these adds a logarithmic factor to our worst-case total parallel time complexity analysis. To the best of our knowledge, these are the first provably polylog$`(L/s)`$ parallel time Delaunay meshing algorithms that generate well-shaped meshes of size optimal to within a constant.

### 1.1 Motivation and Related Work

This work is motivated by the observation that both sequential and parallel implementations of Delaunay refinement algorithms seem to produce the best meshes in practice. However, improvements in the speed of parallel numerical solvers are creating the need for comparable speedups in meshing software: Löhner and Cebral  have reported that improvements in parallel numerical solvers  have resulted in the simulation time of numerous practical systems being dominated by the meshing process.

Quadtree-based methods are an alternative to Delaunay refinement. They also generate well-shaped meshes whose size is within a constant factor of the best possible . In practice, however, they often generate meshes larger than Delaunay refinement on the same input. The parallel complexity of the quadtree-based methods is nevertheless better understood.

Several parallel mesh generation algorithms have been developed. On the theoretical extreme, Bern, Eppstein and Teng  gave a parallel $`O(\mathrm{log}n)`$ time algorithm using $`K/\mathrm{log}n`$ processors to compute a well-shaped quadtree mesh, where $`K`$ is the final mesh size. There is also a simple level-by-level quadtree-based method that is used in practice . One can easily show that this level-by-level based method takes $`O(\mathrm{log}(L/s)+K/p)`$ parallel time, using $`p`$ processors .

Building upon , Miller et al.  developed a parallel sphere-packing based Delaunay meshing algorithm that generates well-shaped Delaunay meshes of optimal size in $`O(\mathrm{log}n)`$ parallel time using $`K/\mathrm{log}n`$ processors. Their method uses a parallel maximal independent set algorithm  to directly generate the set of final mesh points, and then constructs the Delaunay mesh using parallel Delaunay triangulation. As this algorithm has not been implemented, we do not know how the meshes it produces will compare.

Various parallel Delaunay refinement methods have been implemented and been seen to have good performance . These methods address some important issues such as how to partition the domain so as to minimize the communication cost among the processors. Our new analysis on the number of parallel iterations of Delaunay refinement could potentially provide provable bounds on their parallel complexity.

Our work also helps explain the performance of some sequential implementations of Delaunay refinement, especially those which use a Delaunay triangulator as a black-box. In such situations, it is often desirable to minimize the number of calls to the black-box Delaunay triangulator by inserting multiple points at each iteration. Our bounds on the number of iterations provide a bound on the number of calls to the Delaunay triangulator.

We omit the proofs of Lemmas 8, 9, 16, 10, and 18 and Theorems 13, 14, 21, and 22 in this version due to page limitation. A full version of the paper is available at http://www.cs.duke.edu/~ungor/abstracts/parallelDelRef.html.

## 2 Preliminaries

### 2.1 Input Domain

In 2D, the input domain $`\mathrm{\Omega }`$ is represented as a planar straight line graph (PSLG)  — a proper planar drawing in which each edge is mapped to a straight line segment between its two endpoints. The segments express the boundaries of $`\mathrm{\Omega }`$ and the endpoints are the vertices of $`\mathrm{\Omega }`$. The vertices and boundary segments of $`\mathrm{\Omega }`$ will be referred to as input features of $`\mathrm{\Omega }`$. A vertex is incident to a segment if it is one of the endpoints of the segment. Two segments are incident if they share a common vertex. In general, if the domain is given as a collection of vertices only, then the boundary of its convex hull is taken to be the boundary of the input.

Miller et al.  presented a natural extension of PSLGs, called piecewise linear complexes (PLCs), to describe domains in three and higher dimensions. In three dimensions, the domain $`\mathrm{\Omega }`$ is a collection of vertices, segments, and facets where (i) all lower dimensional elements on the boundary of an element in $`\mathrm{\Omega }`$ also belongs to $`\mathrm{\Omega }`$, and (ii) if any two elements intersect, then their intersection is a lower dimensional element in $`\mathrm{\Omega }`$. In other words, a PLC in $`d`$ dimensions is a cell complex with polyhedral cells from $`0`$ to $`d`$ dimensions.

### 2.2 Delaunay Triangulation

Let $`P`$ be a point set in $`^d`$. A simplex $`\tau `$ formed by a subset of $`P`$ points is a Delaunay simplex if there exists a circumsphere of $`\tau `$ whose interior does not contain any points in $`P`$. The Delaunay triangulation of $`P`$, denoted $`Del(P)`$, is a PLC that contains all Delaunay simplices. If the points are in general position, that is, if no $`d+2`$ points in $`P`$ are co-spherical, then $`Del(P)`$ is a simplicial complex.

The Delaunay triangulation of a point set can be constructed in $`O(n\mathrm{log}n)`$ time in 2D  and in $`O(n^{d/2})`$ time in $`d`$ dimensions . A nice survey of these algorithms can be found in .

One way to obtain a triangulation that conforms to the boundary of a PSLG domain is to use a Constrained Delaunay triangulation. Let $`P`$ be the set of vertices of a PSLG $`\mathrm{\Omega }`$. Two points $`p`$ and $`q`$ in $`P`$ are said to be visible from each other if the line segment $`pq`$ does not intersect the interior of any segment in $`\mathrm{\Omega }`$. Three points form a constrained Delaunay triangle if the interior of their circumcircle contains no point from $`P`$ that is visible from all three points. The union of all constrained Delaunay triangles forms a constrained Delaunay triangulation $`CDT(\mathrm{\Omega })`$. Chew developed an algorithm for computing constrained Delaunay triangulations .

A Delaunay triangulation $`T`$ of input and Steiner points is a conforming Delaunay triangulation of a PLC $`\mathrm{\Omega }`$ if every face of $`\mathrm{\Omega }`$ is a union of faces of $`T`$. In 2D, Edelsbrunner and Tan proved that $`O(n^3)`$ additional points are sufficient to generate a conforming triangulation of a PSLG of complexity $`n`$ . A 2D solution proposed by Saalfeld  is extended to 3D by Murphy et al.  and Cohen-Steiner et al. . However, it remains open whether the size of their output is polynomial in the input size or local feature size. The definition of local feature size will be given in Section 3. When the angle between the faces of a PLC is bounded from below, say for example by $`\pi /2`$, then one can apply Delaunay refinement to generate well-shaped conforming triangulations whose size is close to optimal both in two  and three dimensions .

## 3 2D Sequential Delaunay Refinement

In this section, we recall Ruppert’s and Chew’s algorithms for constructing Delaunay meshes of PLSGs in 2D. Following Ruppert , we assume that the angle between two adjacent input segments is at least $`\pi /2`$. Boundary treatments that relax this assumption are discussed in .

In the process of Delaunay refinement, one could either maintain a constrained Delaunay triangulation, or one just keeps track of the set of input segments that are not respected. The first approach does not extend to three dimensions because, in 3D, some PLCs do not have a constrained Delaunay triangulation. We therefore use the second approach.

At each iteration, we choose a new point for insertion from a set of candidate points. There are two kinds of candidate points: (1) the circumcenters of existing triangles, and (2) the midpoints of existing boundary segments.

Let the diametral circle of a segment be the circle whose diameter is the segment. A point is said to encroach a segment if it is inside the segment’s diametral circle. (See Figure 1.)

At iteration $`i`$, the circumcenter of a triangle is a potential candidate for insertion if the triangle is poorly shaped. For example, in Ruppert’s algorithm, a triangle is considered poorly shaped if the ratio of its circumradius to the length of its shortest side is larger than a pre-specified constant $`\beta _R\sqrt{2}`$. Let $`\dot{𝒞}^{(i)}`$ denote the set of all potential candidate circumcenters that do not encroach any segment. Let $`𝒞^{(i)}`$ denote their corresponding circumcircles. Similarly, let $`\dot{}^{(i)}`$ denote the set of all potential candidate circumcenters that do encroach some segment. Let $`^{(i)}`$ denote their corresponding circumcircles.

The midpoint of a boundary segment is a candidate for insertion if (1) the segment is not part of the current Delaunay triangulation, that is, its diametral circle is encroached by some existing mesh points, or (2) its segment is encroached by a circumcenter in $`\dot{}`$. In the latter case, this potential circumcenter candidate is rejected from insertion. Let $`\dot{𝒟}_T^{(i)}`$ be all midpoint candidates of type (1) and let $`\dot{𝒟}_{}^{(i)}`$ be all midpoint candidates of type (2).

The points inserted by the Delaunay refinement are often called Steiner points.

If a quasi-uniform mesh, such as that produced by Chew’s method, is desired , then we use the following notion of poorly shaped triangle: A triangle is poorly-shaped if the ratio of its circumradius to the length of the shortest edge in the current Delaunay triangulation $`T`$ is more than a pre-specified constant $`\beta _C\sqrt{2}`$.

Figure 2 shows the output of the Delaunay refinement illustrating the difference between Chew’s and Ruppert’s refinement. We call these two variants of the Delaunay refinement algorithm Chew’s algorithm and Ruppert’s algorithm.

In their original papers , Chew and Ruppert presented their Delaunay refinement algorithms as particular variations of Algorithm 1 —they specified how to choose the next point at each iteration from the set of candidates. In this paper, we will consider the following variation of Algorithm 1 which is more aggressive in adding boundary points — we choose this variation to parallelize because its analysis is relatively simpler to present.

In this variation, $`^{(i)}`$, $`𝒞^{(i)}`$, and $`𝒟_T^{(i)}`$ are the same as in Algorithm 1. The set $`𝒟_{}^{(i)}`$ is built incrementally. At iteration $`i`$, we compute $`^{(i)}`$ first and let $`𝒟^{(i)}`$ be the set of diametral circles that are encroached by some circumcenters of $`^{(i)}`$. We then set $`𝒟_{}^{(i)}=𝒟_{}^{(i1)}𝒟^{(i)}`$.

In other words, if a segment is encroached by a circumcenter of a poorly-shaped Delaunay triangle, its midpoint will be added to the set of candidate midpoints and remains candidate thereafter. This is in contrast with Algorithm 1, in which an encroached midpoint is added to the set of candidate midpoints only for the next iteration. If another candidate is chosen that is in a circumcircle whose center encroaches the segment, the circumcircle will no longer be in the Delaunay triangulation at the end of the iteration, and hence the segment might not be encroached in the triangulation at the end of the iteration. So, its midpoint might not be a candidate in future iterations.

Assuming that the angle between two adjacent input segments is at least $`\pi /2`$, Chew’s algorithm terminates with well-shaped quasi-uniform meshes, while Ruppert’s algorithm  terminates with a well-shaped Delaunay mesh of the input domain whose elements adapt to the local geometry of the domain. The number of triangles in the mesh generated by Ruppert’s algorithm is asymptotically optimal up to a constant. The proofs of Ruppert’s and Chew’s  that their algorithms terminate with a well-shaped mesh of size within a constant factor of optimal can be easily extended to our variation of Algorithm 1 discussed above. We refer interested readers to  and . Here we give a high level argument and introduce an important concept that will be used in Section 4.2.3 for preprocessing an input domain in parallel.

Given a domain $`\mathrm{\Omega }`$, the local feature size of each point $`x`$ in $`\mathrm{\Omega }`$, denoted by $`\mathrm{lfs}_\mathrm{\Omega }(x)`$, is the radius of the smallest disk centered at $`x`$ that touches two non-incident input features. Ruppert showed that every Delaunay triangle in the final mesh is well-shaped and that the length of the longest edge in each Delaunay triangle is within a constant factor of $`\mathrm{lfs}_\mathrm{\Omega }(x)`$ for each $`x`$ in the interior of the triangle.

Suppose $`M`$ is a mesh generated by our variation of Algorithm 1. Let $`\mathrm{\Omega }^{}`$ be the domain obtained from $`\mathrm{\Omega }`$ by adding to $`\mathrm{\Omega }`$ all mesh points in $`M`$ that are on the boundary segments of $`\mathrm{\Omega }`$. Then we can show (i) for all $`x`$ in $`\mathrm{\Omega }`$, $`\mathrm{lfs}_\mathrm{\Omega }(x)`$ and $`\mathrm{lfs}_\mathrm{\Omega }^{}(x)`$ are within a small constant factor of each other; and (ii) $`M`$ can be obtained by applying Ruppert’s (or Chew’s) variations of Algorithm 1 to $`\mathrm{\Omega }^{}`$. Therefore, the mesh produced by our variation of Algorithm 1 has size within a small constant factor of the one generated by Ruppert’s (or Chew’s) refinement method.

## 4 Parallel 2D Delaunay Refinement

To better illustrate our analysis of parallel Delaunay refinement, we first focus on the case in which the input is a periodic point set (PPS) as introduced by Cheng et al. . See also . We will then extend our results to produce boundary conforming meshes when the input domain is a PSLG.

### 4.1 Input Domain: Periodic Point Sets

If $`P`$ is a finite set of points in the half open unit square $`[0,1)^2`$ and $`^2`$ is the two dimensional integer grid, then $`S=P+^2`$ is a periodic point set . The periodic set $`S`$ contains all points $`p+v`$, where $`pP`$ and $`v`$ is an integer vector. The Delaunay triangulation of a periodic point set is also periodic.

As $`P`$ is contained in the unit square, the diameter of $`P`$ is $`L\sqrt{2}`$. When we refer to the diameter of a periodic point set, we will mean the diameter of $`P`$.

#### 4.1.1 A generic parallel algorithm (PPS)

For a periodic point set, the only candidates for insertion are the circumcenters of poorly shaped triangles. We need a rule for choosing a large subset of the candidates with the property that a sequential Delaunay refinement algorithm would insert each of the points in the subset. Our rule is derived from the following notion of independence among candidates.

###### Definition 1 (Independence)

Two circumcenters $`\dot{c}`$ and $`\dot{c}^{}`$ (and also the corresponding circles $`c`$ and $`c^{}`$) are conflicting if both $`c`$ and $`c^{}`$ contain each other’s center. Otherwise, $`\dot{c}`$ and $`\dot{c}^{}`$ (respectively $`c`$ and $`c^{}`$) are said to be independent.

If two candidates conflict, at most one of them can be inserted. Our rule is to insert a maximal independent set (MIS) of candidates at each iteration. We will show that if an algorithm follows this rule, then it will terminate after a polylogarithmic number of rounds.

In the next few subsections, we will discuss how to generate the independent sets used by the algorithm. But first, we prove that regardless of how one chooses the independent set, our parallel algorithm can be sequentialized. This implies that the algorithm inherits the guarantee of its sequential counterpart that it generates a well-shaped mesh of size that is within a constant factor of optimal.

###### Theorem 2

Suppose $`M`$ is a mesh produced by an execution of the Generic Parallel Delaunay Refinement algorithm. Then $`M`$ can be obtained by some execution of one of the sequential Delaunay refinement algorithms discussed in Section 3.

* Let $`_1,_2,\mathrm{},_k`$ be the sets of vertices inserted by the parallel algorithm above at iterations $`1,\mathrm{},k`$, respectively. We describe a sequential execution that inserts all the points in $`_i`$ before any point of $`_j`$ for $`i<j`$. For each independent set $`_i`$, we insert the candidates according to their circumradius in the order from largest to smallest. For any two circumcenters $`\dot{a},\dot{b}_i`$, assume that the radius of $`a`$ is larger than the radius of $`b`$. This implies that $`\dot{a}`$ can not be in the circumcircle of $`\dot{b}`$, because $`\dot{a}`$ and $`\dot{b}`$ are independent. Therefore, the insertion of $`\dot{a}`$ will not eliminate the triangle of $`\dot{b}`$.

Furthermore, observe that in any sequential execution, the insertion of point $`\dot{p}_i`$ can not eliminate the triangle corresponding to $`\dot{q}_j`$ for any $`i<j`$, for otherwise, $`\dot{q}`$ would not exist in the $`j^{th}`$ iteration of the parallel execution.

Therefore, the parallel and sequential executions terminate with the same Delaunay mesh. $`\mathrm{}`$

To minimize the number of iterations, intuitively, we should choose a maximal independent set of candidates at each iteration. In Section 4.1.3, we will give a geometric algorithm that computes a maximal independent set of candidates efficiently in parallel. Our algorithm makes use of the following observation.

###### Lemma 3

Suppose $`c_a`$ and $`c_b`$ are two conflicting circumcircles at iteration $`i`$, and let $`r_a`$ and $`r_b`$ be their circumradii. Then $`r_b/2<r_a<2r_b`$.

#### 4.1.2 Parallelizing Chew’s Refinement (PPS)

In this section, we show that our parallel implementation of Chew’s refinement only needs $`O(\mathrm{log}(L/s))`$ iterations. The basic argument is very simple — we will show that the radius of the largest Delaunay circle reduces by a factor of 3/4 after some constant number (e.g., 98) of iterations. Because the largest circumradius initially is $`O(L)`$ and the largest circumradius in the final mesh is $`\mathrm{\Omega }(s)`$, the iteration bound of $`O(\mathrm{log}(L/s))`$ follows immediately.

###### Lemma 4

For all $`i`$, let $`r_i`$ be the largest circumradius of a triangle in the Delaunay triangulation at the end of iteration $`i`$. For all $`k98`$, $`r_k3r_{k98}/4`$.

* We assume by way of contradiction that $`r_k>3r_{k98}/4`$. Let $`i=k98`$. Let $`c_k`$ be a circumcircle with radius $`r_k`$ after iteration $`k`$. Let $`\dot{c}_k`$ be the center of $`c_k`$.

For $`jk`$, it is clear that $`c_k`$ is also an empty circle in iteration $`j`$, because the refinement process only adds new points. But, $`c_k`$ might not be a circumcircle at iteration $`j`$. We now show that for each iteration $`j`$, where $`ijk`$, there exists a circumcircle $`c_j^{}`$ with center $`\dot{c}_j^{}`$, and radius $`r_j^{}`$ such that (1) $`\dot{c}_j^{}\dot{c}_k3r_i/4`$ and (2) $`r_j^{}3r_i/4`$.

Let $`p_k`$, $`q_k`$ and $`t_k`$ be the vertices of the Delaunay triangle at iteration $`k`$ that defines $`c_k`$. We will alter $`c_k`$ in three stages to produce a suitable $`c_j^{}`$ that is the circumcircle of three points that exist at stage $`j`$: $`p_j`$, $`q_j`$ and $`t_j`$.

1. Dilate $`c_k`$ until it touches a mesh point $`p_j`$. Note that $`p_j`$ might well be $`p_k`$, $`q_k`$, or $`t_k`$, so, $`c_k`$ might not actually expand at all during this step.
2. Grow the circle by moving its center away from $`p_j`$ along the ray $`\stackrel{}{p_j\dot{c}_k}`$, and maintaining the property that $`p_j`$ lies on the boundary of the circle, until it touches a mesh point $`q_j`$.
3. Continue to grow the circle, maintaining its contact with $`p_j`$ and $`q_j`$, moving its center away from the chord $`\overline{p_jq_j}`$, until it touches a vertex $`t_j`$.

The resulting circle $`c_j^{}`$ is a circumcircle of a Delaunay triangle $`p_jq_jt_j`$ at iteration $`j`$. Moreover, $`p_jq_jt_j`$ is a poorly-shaped triangle because its circumradius $`r_j^{}`$ is at least $`r_k`$. Thus, its center $`\dot{c}_j^{}`$ is a candidate at iteration $`j`$. Note also that $`r_j^{}r_k3r_i/4`$.

Consider the triangle $`p_j\dot{c}_k\dot{c}_j^{}`$, which is non-acute at vertex $`\dot{c}_k`$. Let $`x=|\dot{c}_kp_j|`$ and $`y=|\dot{c}_k\dot{c}_j^{}|`$. Since $`|\mathrm{}\dot{c}_j^{}\dot{c}_kp_j|`$ is non-acute $`(r_j^{})^2x^2+y^2`$. As $`r_j^{}`$ is the radius of a Delaunay triangle and $`ji`$, $`r_j^{}r_i`$. Combining this fact with $`xr_k>3r_i/4`$, we find $`r_j^{}<x+r_i/4`$. So we can write, $`(x+r_i/4)^2x^2+y^2`$. By simplifying this inequality to $`xr_i/2+r_i^2/16y^2`$ and substituting $`xr_i`$, we derive $`9r_i^2/16y^2`$. Hence, $`y=\dot{c}_j^{}\dot{c}_k3r_i/4`$.

Because $`c_j^{}`$ is empty at the end of iteration $`j`$, we know $`\dot{c_j^{}}`$ was not chosen during iteration $`j`$. Because the independent set of candidates that we select is maximal, there must be another circumcircle $`c_j^{\prime \prime }`$ chosen in iteration $`j`$ that conflicts with $`c_j^{}`$. By Lemma 3, the radius of $`c_j^{\prime \prime }`$ is at least one half of the radius of $`c_j^{}`$, and so is at least $`|r_j^{}|/23r_i/8`$. Moreover, the radius of $`c_j^{\prime \prime }`$ is at most $`r_j^{}`$ and hence at most $`r_i`$. So $`\dot{c}_j^{\prime \prime }\dot{c}_j^{}r_i`$. Hence,

$$\dot{c}_j^{\prime \prime }\dot{c}_k\dot{c}_j^{\prime \prime }\dot{c}_j^{}+\dot{c}_j^{}\dot{c}_kr_i+3r_i/47r_i/4.$$

Let $`\dot{C}^{\prime \prime }=\{\dot{c}_{i+1}^{\prime \prime },\dot{c}_{i+2}^{\prime \prime },\mathrm{},\dot{c}_j^{\prime \prime },\mathrm{},\dot{c}_k^{\prime \prime }\}`$, and let $`C^{\prime \prime }`$ be the corresponding set of circumcircles. As $`\dot{c}_l^{\prime \prime }`$ is inserted during round $`l`$ for each $`l`$, each circle $`c_j^{\prime \prime }C^{\prime \prime }`$ is empty of all the centers $`\dot{c}_l^{\prime \prime }`$ for $`l<j`$. So the centers in $`\dot{C}^{\prime \prime }`$ are pairwise at least $`3r_i/8`$ away from each other. Thus, one can draw disjoint circles of radius $`3r_i/8`$ around each of these points. The annulus containing these disjoint circles has area at most $`\pi (7/4+3/16)^2r_i^2\pi (3/43/16)^2r_i^2`$. So, one can pack at most

$$\frac{\pi [(7/4+3/16)^2(3/43/16)^2]r_i^2}{\pi (3/16)^2r_i^2}=97$$

disjoint circles of radius $`3r_i/16`$ in this region. This implies $`|C^{\prime \prime }|=ki97`$, a contradiction. $`\mathrm{}`$

###### Theorem 5

Our parallel implementation of Chew’s refinement algorithm takes at most $`98\mathrm{log}_{4/3}(L/s)`$ iterations.

#### 4.1.3 Parallel Computation of MIS

One can use the parallel maximal independent set algorithm of Luby  to compute a parallel independent set of candidates for each iteration in $`O(\mathrm{log}^2n)`$ parallel time. In this section, we will explain how we can exploit the geometric structure of the independence relation to compute a maximal independent set in a constant parallel time.

We will make extensive use of the result of Lemma 3 that two circumcircles are conflicting at iteration $`j`$ only if their radii are within a factor of 2 of each other.

###### Lemma 6

At iteration $`j`$, if there are $`n_j`$ circumcircles, then a maximal independent set of candidates for Delaunay refinement can be computed in constant parallel time using $`n_j`$ processors.

* Let $`C_h^j`$ be the set of circumcircles of radius more than $`L/2^{h+1}`$ and less than or equal to $`L/2^h`$, where $`h`$ ranges from $`0`$ to $`\mathrm{log}(L/s_j)`$ and $`s_j`$ is the smallest circumradius at iteration $`j`$. Note that a circumcircle in $`C_h^j`$ does not conflict with any circumcircle in $`C_l^j`$ if $`l>h+1`$.

To compute a maximal set of non-conflicting candidates, we first in parallel find a maximal independent sets of circumcircles in $`C_h^j`$, independently for all even $`h`$. We will show below that a maximal independent set of circumcircles in $`C_h^j`$ can be computed in constant time in parallel. Let $`I_{even}^j`$ be the set of independent circumcircles computed. Then in one parallel step, we can eliminate all conflicting circumcircles in $`_{h:odd}C_h^j`$. We then compute a maximal independent set for remaining circumcircles in $`C_h^j`$ for all odd $`h`$. Let this set be $`I_{odd}^j`$. Then $`I_{even}^jI_{odd}^j`$ is a maximal independent set of circumcircles for iteration $`j`$.

Note that all circumcircles in $`C_h^j`$ have radius between $`L/2^{h+1}`$ and $`L/2^h`$. If we divide the square containing all circumcenters into a $`2^h`$-by-$`2^h`$ grid, then any circumcenter that is conflict with a circumcenter in the grid box $`(x,y)`$ must lie either in grid box $`(x,y)`$ or one of its eight grid neighbors.

We color grid boxes $`(x,y)`$ with color $`(xmod3,ymod3)`$. We then cycle through the 9 color classes and, by a method we will explain momentarily, find a maximal independent set of the candidates in each grid-box of the current color in parallel. We then eliminate in parallel the conflicting circumcenters that are in the color classes that have not yet been processed.

Finally, we explain how to compute a maximal independent set among the candidates that lie in a given grid-box. First notice that any maximal independent set of candidates in a grid-box can have at most a constant number of members, and hence a maximal independent set can be found by a constant number of parallel selection-elimination operations: choose a center that has not been eliminated, and in parallel eliminate any centers with which it conflicts.

In a parallel system that supports primitives such fetch\_and\_add, test\_and\_set, or parallel scan, we can use such a primitive to select in constant time a candidate in a grid box. The processor that holds this candidate becomes a “leader” in that round and broadcasts its candidate so that the conflicting candidates can be eliminated. With these primitives, our algorithm can be implemented in parallel constant time. However, if the parallel system does not support these primitives, then for each grid cell we can emulate parallel scan to select a leader in $`O(\mathrm{log}n)`$ time, where $`n`$ is the number of candidate centers in the cell.

In general, many grid cells are empty and there is no need to generate them at all. We can use hashing to select grid cells that are not empty. The idea is very simple, each candidate center can compute the coordinates of its grid cell from the coordinates of its center and its radius. We can hash grid cells using their coordinates and therefore, all candidate centers belonging to a grid cell can independently generate the hash identity of the cell. We can then use parallel primitives discussed in the paragraph above to support the computation of a maximal independent set of candidates for all non-empty grid cells. $`\mathrm{}`$

#### 4.1.4 Parallelizing Ruppert’s Refinement (PPS)

In this section, we show that our parallelization of Ruppert’s method for periodic point sets in 2D takes $`O(\mathrm{log}^2(L/s))`$ iterations. For simplicity, we give an analysis for the case $`\beta _R=\sqrt{2}`$, although our analysis can be easily extended to the case when $`\beta _R=1+ϵ`$, for any $`ϵ>0`$. We recall that $`\beta _R`$ is the threshold of the ratio of the radius to shortest edge-length defining a poorly shaped triangle. Thus, for $`\beta _R=\sqrt{2}`$, inserting the circumcenter of a poorly shaped triangle whose shortest edge is $`h`$ introduces new Delaunay edges of length at least $`\sqrt{2}h`$.

Let $`s`$ be the length of the shortest edge in the initial Delaunay triangulation. At each iteration, we assign an edge to class $`_i`$ if its length is in $`[\sqrt{2}^{i1}s,\sqrt{2}^is)`$. Similarly, we assign a Delaunay triangle to $`_i`$ if its shortest edge has length in $`[\sqrt{2}^{i1}s,\sqrt{2}^is)`$. There are at most $`\mathrm{log}_\sqrt{2}(L/s)`$ of such classes. Using this definition, we can state and analyze the Parallel Ruppert’s Refinement Algorithm.

###### Theorem 7

Given a periodic point set in 2D of diameter $`L`$, the Parallel Ruppert’s Refinement Algorithm takes $`O(\mathrm{log}^2(L/s))`$ iterations.

* Lemmas 8 and 9 prove that after the $`i`$th iteration of the outer loop, each Delaunay triangle touching an edge in class $`_i`$ will be well-shaped, and successive iterations cannot degrade the shape of the Delaunay triangles touching that edge. Lemma 10 implies that during each iteration of the outer loop, the inner loop of the algorithm will execute at most $`O(\mathrm{log}(L/s))`$ times. As the outer loop is executed $`O(\mathrm{log}(L/s))`$ times, the whole algorithm takes at most $`O(\mathrm{log}^2(L/s))`$ iterations. $`\mathrm{}`$

###### Lemma 8

During the $`i`$th iteration of the outer loop of the Parallel Ruppert’s Refinement Algorithm, no Delaunay edges are added to or removed from class $`_i`$.

###### Lemma 9

Suppose $`e`$ is an edge in $`_j`$ where $`ji`$. Then during the $`i`$th outer loop, the radius-edge ratio of triangles containing $`e`$ does not increase.

###### Lemma 10

Let $`e_i`$, and let $`r_l`$ be the radius of the larger of the two circumcircles containing $`e`$ at the end of the $`l^{th}`$ iteration of the inner loop during the $`i`$th iteration of the outer loop. Then, at the end of iteration $`k=l+81`$ of the inner loop, either (1) both Delaunay triangles containing $`e`$ are well-shaped, or (2) $`r_k3r_l/4`$ where $`r_k`$ is the radius of the larger of the two circumcircles containing $`e`$ after iteration $`k`$.

### 4.2 Input Domain: PSLG

In this subsection, we extend our parallel algorithm for generating a Delaunay mesh from a domain given by a periodic point set to a domain defined by a planar straight-line graph. Following Ruppert, we assume that the angle between two adjacent input segments is at least $`\pi /2`$. A key step in Delaunay refinement for a domain specified by a PSLG is to properly add points to the boundary segments so that the Delaunay mesh is conforming to the boundary. In our parallel algorithm, we make our mesh conform to the boundary in two steps: First, we give, in Section 4.2.3, an $`O(\mathrm{log}L/s)`$ time parallel preprocessing algorithm to insert points to input segments so that the initial Delaunay mesh is conforming to the boundary and no diametral circle intersects any other non-incident input features. Second when a segment is encroached during parallel Delaunay refinement, we include its midpoint as candidate for insertion.

The preprocessing step might not be needed to implement our parallel algorithm: one could probably add points to the boundary as needed. However, the preprocessing step simplifies our analysis in this Section by greatly reducing the number of cases in the analysis.

#### 4.2.1 A generic parallel algorithm (PSLG)

After applying the preprocessing step, the initial Delaunay triangulation is conforming to the boundary and no diameter circle contains any point of the triangulation. We will maintain this invariant in our algorithm.

In order to perform parallel refinement, as in Section 4.1.1, we need a rule of independence among candidates for refining boundary segments and poorly shaped triangles. We first recall the set of candidates for insertion defined in Section 3.

Let $``$ be the set of circumcircles of poorly shaped triangles whose centers $`\dot{}`$ encroach some boundary segments. Let $`𝒞`$ be the set of circumcircles of poorly shaped triangles whose centers $`\dot{𝒞}`$ don’t encroach any boundary segments. Let $`𝒟`$ is the set of diametral circles that are encroached by some centers in $`\dot{}`$. So, $`\dot{𝒞}\dot{𝒟}`$ are candidate points for insertion.

We will still apply Definition 1 to determine whether two circumcenters from $`\dot{𝒞}`$ are independent. Because the angle between two adjacent input segments is at least $`\pi /2`$, after preprocessing, any two diametral circles from $`\dot{𝒟}`$ are not overlapping. Every two diametral centers from $`\dot{}`$ are independent.

We will use the following definition of independence between a diametral center in $`\dot{𝒟}`$ and a circumcenter in $`\dot{𝒞}`$. Note that because a circumcenter in $`\dot{𝒞}`$ does not encroach any boundary segment, a diametral circle of $`𝒟`$ does not contain any center in $`\dot{𝒞}`$.

###### Definition 11

A circumcenter $`\dot{c}\dot{𝒞}`$ and a diametral center $`\dot{d}\dot{𝒟}`$ are conflicting if (i) $`\dot{d}`$ is inside $`c`$; and (ii) the radius of $`c`$ is smaller than $`\sqrt{2}`$ times the radius of $`d`$. Otherwise, $`\dot{c}`$ and $`\dot{d}`$ (also $`c`$ and $`d`$) are independent.

This definition of independence is motivated by the following lemma first proved by Ruppert .

###### Lemma 12

If a circumcircle $`c`$ of radius $`r_c`$ encroaches a diametral circle $`d`$ of radius $`r_d`$, then $`r_dr_c/\sqrt{2}`$.

The following theorem extends Theorem 2 for domains given by PSLGs.

###### Theorem 13

For a domain $`\mathrm{\Omega }`$ specified by a PSLG, suppose $`M`$ is a mesh produced by an execution of the parallel algorithm above. Then $`M`$ can be obtained by some execution of one of the sequential Delaunay refinement algorithms discussed in Section 3.

#### 4.2.2 Parallelizing Chew’s Refinement (PSLG)

To parallelize Chew’s algorithm for domain defined by a PSLG, we apply Algorithm 4 and use a maximal independent set of the candidates at each iteration. In addition, because each pair of diametral centers in $`\dot{𝒟}`$ is independent, we include all centers $`\dot{𝒟}`$ in the independent set. The parallel algorithm of Section 4.1.3 can be used to construct the maximal independent set.

###### Theorem 14

Our parallel implementation of Chew’s refinement algorithm takes $`O(\mathrm{log}(L/s))`$ iterations for a domain given by a PSLG, where $`L`$ is the diameter of the domain and $`s`$ is smallest local feature size.

#### 4.2.3 Parallel Preprocessing

In the algorithm and proof presented in the last subsection, we assume that the boundary of the domain has been preprocessed to satisfy the following property.

###### Definition 15 (Strongly Conforming)

A domain $`\mathrm{\Omega }`$ specified by a PSLG is strongly conforming if no diametral circle contains any vertex or intersects any other non-incident input features.

Clearly, if $`\mathrm{\Omega }`$ is strongly conforming, then the Delaunay triangulation of the vertices of $`\mathrm{\Omega }`$ is conforming to $`\mathrm{\Omega }`$.

We will use the following parallel method to preprocess a domain $`\mathrm{\Omega }`$ to make it strongly conforming. This method repeatedly adds midpoints to boundary segments whose diametral circles intersect non-incident input features.

###### Lemma 16

Parallel Boundary Preprocessing terminates in $`O(\mathrm{log}(L/𝐬))`$ iterations.

In the scheme above, we can grow a quadtree level by level to support the query of whether the diametral circle of a segment intersects another non-incident feature. The number of levels of the quadtree that we need to grow is at most $`\mathrm{log}(L/s)`$. As shown in , one can use balanced quadtree to approximate local feature size function of $`\mathrm{\Omega }`$ to within a constant factor. Therefore, using a balanced quadtree as , we can preprocess the domain in $`\mathrm{log}(L/s)`$ parallel time so that the preprocessed domain is strongly feature conforming as defined below.

###### Definition 17 (Strongly Feature Conforming)

Let $`\alpha >2`$ be a constant. A domain $`\mathrm{\Omega }`$ specified by a PSLG is strongly feature conforming with parameter $`\alpha `$ if it is strongly conforming, and in addition, the length of each segment is no more than $`\alpha `$ times the local feature size of its midpoint.

In the next subsection, we will present a parallel implementation of Ruppert’s algorithm for domains that are strongly feature conforming and show that it terminates in $`O(\mathrm{log}^2(L/s))`$ iterations.

We use the following lemma to show that the size optimality of our results are not affected much by the preprocessing.

###### Lemma 18

Let $`\mathrm{\Omega }`$ and $`\mathrm{\Omega }^{}`$ denote the input before and after preprocessing, respectively. Then, for any point $`x`$ in these domains, $`\mathrm{lfs}_\mathrm{\Omega }(x)/3\mathrm{lfs}_\mathrm{\Omega }^{}(x)\mathrm{lfs}_\mathrm{\Omega }(x)`$.

#### 4.2.4 Parallelizing Ruppert’s Refinement (PSLG)

In this section, we show that our parallelization of Ruppert’s method for a domain given by a PSLG takes $`O(\mathrm{log}^2(L/s))`$ iterations. Again, for simplicity, we will only give an analysis for the case when $`\beta _R=\sqrt{2}`$.

The parallel algorithm follows basic steps of the parallel Ruppert’s Refinement presented earlier in Section 4.1.4. But first, we apply the parallel preprocessing algorithm of Section 4.2.3 so that the preprocessed domain is strongly feature conforming. So below we can assume that $`\mathrm{\Omega }`$ is strongly conforming.

Let $`s`$ be smallest local feature of $`\mathrm{\Omega }`$. At each iteration, we assign an edge to class $`_i`$ if its length is in $`[\sqrt{2}^{i1}s,\sqrt{2}^is)`$. Similarly, we assign a Delaunay triangle to $`_i`$ if its shortest edge has length in $`[\sqrt{2}^{i1}s,\sqrt{2}^is)`$. There are at most $`\mathrm{log}_\sqrt{2}(L/s)`$ of such classes.

###### Theorem 19

Given a domain specified by a PSLG, the Parallel Ruppert’s Refinement Algorithm takes $`O(\mathrm{log}^2(L/s))`$ iterations.

The proof of Theorem 19 is essentially the same as the proof of Theorem 7 where we need to address the following two issues.

1. The center of a circumcircle could potentially encroach a boundary segment whose length is much larger than that the circumradius.
2. The insertion of a midpoint on the boundary could potentially introduce smaller edges.

To address the first issue, we apply parallel processing algorithm of Section 4.2.3 and hence assume $`\mathrm{\Omega }`$ is strongly feature conforming. Hence if a circumcenter encroaches a boundary segment, the circumradius and the length of the segment are with a constant factor of each other. In addition, because each boundary segment can only be split at most a constant times in the refinement, it can not introduce smaller edges too many times.

## 5 3D Delaunay Refinement

A 3D domain is specified by a PLC (see Section 2.1). In this section, we assume that the angle between any two intersecting elements, when one is not contained in the other, is at least $`90^{}`$. There are three kinds of spheres associated with a 3D Delaunay mesh that we are interested: the circumspheres, the diametral spheres, and the equatorial sphere given below.

###### Definition 20

The equatorial sphere of a triangle in 3D is the smallest sphere that passes through its vertices. A triangular subfacet of a PLC is encroached if the equatorial sphere is not empty.

Chew’s algorithm extends naturally to 3D. In , Shewchuk developed a 3D extension of Ruppert’s algorithm. In Shewchuk’s refinement, given below, a tetrahedron is bad if the ratio of its circumradius to its shortest edge, referred as the radius-edge ratio, is more than a pre-specified constant $`\beta _S2`$.

### 5.1 Parallel 3D Delaunay Refinement

In this subsection, we show that our results for a domain given by a periodic point set can be extended from two dimensions to three dimensions to parallelize both Chew’s and Shewchuk’s algorithm. So far, we have not completed the analysis for domains specified by PLCs, although we think similar results can be obtained.

The following is a parallel Delaunay refinement algorithm for domains specified by 3D periodic point sets.

To parallelize Chew’s 3D refinement, we use a maximal independent set of candidate centers in Algorithm 8. With almost the same proof as we have presented in Section 4.1.2, we can show that the number of iterations needed is $`1076\mathrm{log}(L/s)`$.

###### Theorem 21

Suppose $`M`$ is a mesh produced by an execution of the 3D Generic Parallel Delaunay Refinement algorithm. Then $`M`$ can be obtained by some execution of the sequential Delaunay refinement algorithm.

#### 5.1.1 Parallelizing Shewchuk’s Refinement

We will present our analysis for the case when $`\beta _R=\sqrt{2}`$, although our analysis can be easily extended to the case when $`\beta _R=1+ϵ`$, for any $`ϵ>0`$. Thus, for $`\beta _R=\sqrt{2}`$, inserting the circumcenter of a poorly shaped triangle whose shortest edge is $`h`$ introduces new Delaunay edges of length at least $`\sqrt{2}h`$.

Let $`s`$ be the length of the shortest edge in the initial Delaunay triangulation. At each iteration, we assign an edge to class $`_i`$ if its length is in $`[\sqrt{2}^{i1}s,\sqrt{2}^is)`$. Similarly, we assign a Delaunay tetrahedra to $`_i`$ if its shortest edge has length in $`[\sqrt{2}^{i1}s,\sqrt{2}^is)`$. There are at most $`\mathrm{log}_\sqrt{2}(L/s)`$ of such classes.

Our parallel implementation of Shewchuk’s algorithm is analogous to our parallel implementation of Ruppert’s algorithm. In addition, our proof in 3D is also analogous to the proof in 2D.

###### Theorem 22

For a given periodic point set $`P`$ in $`^3`$ of diameter at most $`L`$, if the length of the shortest edge in the mesh generated by Shewchuk’s refinement is $`s`$, then parallel Shewchuk refinement takes $`O(\mathrm{log}^2(L/s))`$ iterations to generate a bounded radius-edge ratio mesh.

## 6 Discussion

Polylogarithmic upper bounds on the number of parallel iterations presented in Sections 4 and 5.1 constitutes the main component of the analyses of our parallel algorithms. At each iteration, our algorithms perform two main operations: i) compute a maximal independent set of points for parallel insertion; ii) update the Delaunay triangulation inserting all these points. For the first one, we proposed a new constant time parallel algorithm. For the second, we suggested to use an existing logarithmic time parallel Delaunay triangulation algorithm. These immediately imply polylogarithmic total time complexity for our parallel Delaunay refinement algorithms.

We opted for simplicity in our analyses. So, the constants in lemmas 4 and 10 are probably not optimal and likely to be much smaller in practice than 98 and 81.

The 3D extension of Chew’s and Shewchuk’s algorithms do not always guarantee that the resulting mesh has an aspect-ratio bounded by a constant. However, they both guarantee a constant bound on the ratio of the circumradius to the length of the shortest edge (the radius-edge ratio) of any tetrahedra in the final mesh. So, the meshes these two algorithms generate might potentially contains slivers, which are elements with close to zero aspect-ratio but with a constant radius-edge ratio. Several quality enhancing and guaranteeing meshing algorithms  have been developed recently. Cheng et al.  and Edelsbrunner et al.  have already given parallel complexity of their sliver removal algorithms. Our framework can be used to analyze parallel complexity of the other two algorithms, by Chew  and Li and Teng .

We conclude the paper with two conjectures.

* There is a parallel implementation of Ruppert’s  and Shewchuk’s  algorithm that runs in $`O(\mathrm{log}(L/s))`$ iterations.
* There is a parallel Ruppert’s  and Shewchuk’s  algorithm that runs in $`O(\mathrm{log}n)`$ time where $`n`$ is the input complexity. Notice that Bern et al.  showed that the quadtree algorithm can be implemented in $`O(\mathrm{log}n)`$ time with $`K`$ processors.

We would also like to see results that establish the parallel complexity of other mesh generation algorithms such as sink insertion .

## 7 Acknowledgments

We thank Jeff Erickson, Sariel Har-Peled and Dafna Talmor for helpful conversations and comments on the paper. We also thank Jonathan Shewchuk for helpful emails about boundary assumptions and boundary preprocessing for Ruppert’s algorithm.


text:
# 1 Introduction

## 1 Introduction

The $`bs\gamma `$ transition proceeds by an effective flavor-changing neutral current. The Standard Model prediction, in which the $`t`$ penguin loop gives the largest contribution, has been calculated to next–to–leading order . Measurements of the inclusive branching fraction have been used to constrain new physics contributions to the decay amplitude .

Since $`bs\gamma `$ is a two–body decay process, the photon energy, $`E_\gamma `$, in the $`B`$ rest frame, is related to the recoil hadronic mass, $`M_{Had}`$, by:

$$E_\gamma =\frac{M_B^2M_{Had}^{}{}_{}{}^{2}}{2M_B}.$$
(1)
In our analysis we fully reconstruct twelve exclusive $`bs\gamma `$ decays, and use the hadronic mass spectrum to measure the photon energy spectrum in the $`B`$ rest frame. Note that in a fully inclusive analysis , where only the photon is measured, the gamma energy is smeared by the resolution of the calorimeter and by the motion of the $`B`$ meson in the $`\mathrm{{\rm Y}}(4S)`$ rest frame. With our semi-exclusive approach, the gamma energy resolution depends only on the $`M_{Had}`$ resolution, which is a few $`\mathrm{MeV}/c^2`$.

We fit the $`M_{Had}`$ spectrum with the model proposed by Kagan and Neubert , which predicts the shape of the gamma energy spectrum using heavy quark effective theory (HQET). There are two main parameters, an effective $`b`$ quark mass, $`m_b`$, and $`\lambda _1`$, which is related to the kinetic energy of the $`b`$ quark in the $`B`$ meson. The model introduces the known $`K^{}\gamma `$ contribution  using local quark-hadron duality to convert the portion of the spectrum below a cutoff mass into a $`K^{}(892)`$ Breit-Wigner shape. The cutoff mass is a free parameter in the model, but is expected to be between 1.0 and 1.2 $`\mathrm{GeV}/c^2`$. Above this cutoff the inclusive model does not explicitly include any higher resonances.

We find a correlated band of allowed values for the parameters $`m_b`$ and $`\lambda _1`$, and a cutoff mass of about 1.11 $`\mathrm{GeV}/c^2`$. By measuring the first moment of the photon energy spectrum we constrain the HQET parameter $`\overline{\mathrm{\Lambda }}`$ and hence $`m_b`$. A fit to the $`M_{Had}`$ spectrum with this constraint applied to $`m_b`$ gives the inclusive branching fraction and a range of values for $`\lambda _1`$. There is currently much interest in the HQET parameters, since they are needed to extract $`V_{ub}`$ and $`V_{cb}`$ from semileptonic $`B`$ decays .

## 2 The BABAR detector and dataset

This paper describes preliminary results obtained from 20.7 $`\text{ fb}^1`$ of data recorded with the BABAR detector at the PEP-II asymmetric $`e^+e^{}`$ storage ring during 1999-2000. The data correspond to a total of $`(22.74\pm 0.36)10^6`$ $`B\overline{B}`$ pairs collected on the $`\mathrm{{\rm Y}}(4S)`$ resonance.

The BABAR detector is described in detail elsewhere . We briefly summarize the detector systems most relevant to the current paper. The BABAR detector contains a five–layer silicon vertex tracker (SVT) and a forty–layer drift chamber (DCH) situated in a 1.5 T solenoidal magnetic field. These devices detect charged particles and measure their momentum and ionization energy loss ($`dE/dx`$). The transverse momentum resolution is $`\sigma _{p_t}/p_t=(0.13\pm 0.01)\%p_t+(0.45\pm 0.03)\%`$, with $`p_t`$ in $`\mathrm{GeV}/c`$. Photons are detected in a CsI(Tl) crystal electromagnetic calorimeter (EMC). The EMC detects photons with energies as low as 20 $`\mathrm{MeV}`$. The nominal EMC resolution for photons and electrons is $`\sigma _E/E=2.3\%/E^{1/4}1.9\%`$, with $`E`$ in $`\mathrm{GeV}`$. The charged particle identification (PID) combines $`dE/dx`$ measurements in the SVT and DCH with particle velocity measurements, obtained with an internally reflecting ring-imaging Cherenkov detector (DIRC) of quartz bars surrounding the DCH.

The interactions of particles traversing the detector are simulated using GEANT4 . Variations in detector conditions and beam-induced backgrounds are taken into account in these simulations. Signal and generic background samples are used to study the effect of the event selection criteria and to estimate the backgrounds. The generic background simulations consist of $`e^+e^{}q\overline{q}`$ $`(q=u,d,s,c)`$ and $`B\overline{B}`$ events, where the $`bs\gamma `$ signal events have been removed.

## 3 Event Selection

We reconstruct 12 final states111Charge conjugate modes are assumed throughout this paper. formed from a high–energy photon, a charged or neutral kaon, and 1-3 pions, one of which can be neutral:

$`K^+\pi ^{}`$, $`K_S\pi ^0`$, $`K^+\pi ^{}\pi ^0`$, $`K_S\pi ^+\pi ^{}`$, $`K^+\pi ^{}\pi ^+\pi ^{}`$, $`K_S\pi ^+\pi ^{}\pi ^0`$

$`K^+\pi ^0`$, $`K_S\pi ^+`$, $`K^+\pi ^{}\pi ^+`$, $`K_S\pi ^+\pi ^0`$, $`K^+\pi ^{}\pi ^+\pi ^0`$, $`K_S\pi ^+\pi ^{}\pi ^+`$

These make up about 50% of the total $`bs\gamma `$ rate for $`M_{Had}`$ $`<`$ 2.4 $`\mathrm{GeV}/c^2`$.

We select events with at least one photon with an energy $`1.5<E_\gamma ^{}<3.5\mathrm{GeV}`$ in the $`e^+e^{}`$ centre-of-mass (CM) frame. If the photon can be combined with another photon with energy $`>`$ 50 $`\mathrm{MeV}`$ (250 $`\mathrm{MeV}`$) to make a $`\pi ^0`$ ($`\eta `$) with a mass between 115 $`\mathrm{MeV}/c^2`$ (508 $`\mathrm{MeV}/c^2`$) and 155$`\mathrm{MeV}/c^2`$ (588 $`\mathrm{MeV}/c^2`$), it is vetoed.

In the hadronic system, a $`\pi ^0`$ is made from two photons, with a minimum photon energy of 30 $`\mathrm{MeV}`$, and a minimum $`\pi ^0`$ energy of 200 $`\mathrm{MeV}`$. Well reconstructed tracks are used for charged pions, but charged kaons must satisfy additional kaon identification criteria. Candidates for $`K_S\pi ^+\pi ^{}`$ decays are selected by requiring a mass between $`489`$ and $`507\mathrm{MeV}/c^2`$, and a decay length $`>2`$ mm. These kaon selections reduce any $`bd\gamma `$ background contributions to a negligible level.

The mass of the hadronic system is restricted to:

$$0.6<M_{Had}<2.4\mathrm{GeV}/c^2$$

where the upper limit corresponds to a requirement on $`E_\gamma >2.094`$ $`\mathrm{GeV}`$. Most of the background is accounted for by events in which the photon comes from initial state radiation (ISR) or from $`\pi ^0`$/$`\eta `$ decays that survive the veto. We suppress these by requiring $`\mathrm{cos}\theta _T^{}<0.7`$, where $`\mathrm{cos}\theta _T^{}`$ is the angle in the CM system between the high–energy photon candidate and the thrust axis of the rest of the event calculated after removing the $`B`$ candidate. Further suppression is achieved by combining several event shape variables in a Fisher discriminant. The variables include the ratio of the second and zeroth Fox Wolfram moments  computed in the CM and in the photon recoil system, the CM direction of the $`B`$ candidate, and energy flow in a 20 cone along the photon direction and in a 40 cone opposite the photon.

The hadronic system is combined with a high–energy photon, and accepted as a $`B`$ candidate if the energy in the CM system, $`E_B^{}`$, is within 150 $`\mathrm{MeV}`$ of the CM beam energy, $`E_{beam}^{}`$, :

$$|\mathrm{\Delta }E|=|E_B^{}E_{beam}^{}|<150\mathrm{MeV}$$
(2)
and the beam-energy substituted mass $`m_{\mathrm{ES}}`$:

$$m_{\mathrm{ES}}=\sqrt{E_{beam}^2p_B^2}>5.21\mathrm{GeV}/c^2$$
(3)
where $`p_B^{}`$ is the $`B`$ momentum in the CM frame.

If there is more than one $`B`$ candidate left in an event, we choose the candidate with the minimum $`|\mathrm{\Delta }E|`$. Then, we rescale the measured photon energy to give $`\mathrm{\Delta }E=0\mathrm{GeV}`$, and recalculate $`m_{\mathrm{ES}}`$. This corrects for energy leakage and improves the $`m_{\mathrm{ES}}`$ resolution, particularly for candidates without a $`\pi ^0`$.

## 4 Background Subtraction

The background is mainly from continuum, although $`B\overline{B}`$ backgrounds become significant at high $`M_{Had}`$ values. The continuum background $`m_{\mathrm{ES}}`$ distribution can be fit by an ARGUS  function. When divided into bins in $`M_{Had}`$, the generic continuum Monte Carlo samples show a slight dependence of the shape parameter of the ARGUS function on $`M_{Had}`$, which we take into account in our fits.

The $`B\overline{B}`$ background is mainly composed of final states with high multiplicity and high hadronic mass, where the photon is the daughter of a $`\pi ^0`$ or $`\eta `$ that decayed asymmetrically. We fit the generic $`B\overline{B}`$ Monte Carlo samples with an ARGUS function for the non-peaking background and a Crystal Ball function  for the peaking background, such as $`BD^{()}\rho ^{}`$ where only a low–energy photon is missing. The shape parameters of the ARGUS functions for the continuum and non-peaking $`B\overline{B}`$ backgrounds are sufficiently similar that it is possible to use a combined ARGUS function to describe the sum of the two backgrounds.

Figure 1 shows the continuum and non-peaking $`B\overline{B}`$ backgrounds for $`m_{\mathrm{ES}}>5.21\mathrm{GeV}/c^2`$, in bins of hadronic mass. The level of $`B\overline{B}`$ background is only significant for $`M_{Had}>1.8`$ $`\mathrm{GeV}/c^2`$.

We treat $`bs\gamma `$ events as “cross-feed” background when they are reconstructed with the correct photon, but with the wrong hadronic system, since this gives an incorrect measurement of $`M_{Had}`$. Cross-feed events come from three sources, which are, in decreasing order of importance:

* Final states that are not considered in the analysis, e.g. $`K_L`$, $`2\pi ^0`$, $``$ 5–body hadronic final states, and events with $`M_{Had}>2.4\mathrm{GeV}/c^2`$;
* Final states which are considered, but in which one of the final state particles is not detected; and
* Multiple candidate events where the wrong candidate is chosen.

The cross-feed $`m_{\mathrm{ES}}`$ distribution is fit with an ARGUS function plus the Crystal Ball function to allow for a possible peaking contribution. The cross-feed ARGUS function rises towards the signal region. Note that the level of cross-feed background is proportional to the inclusive branching fraction. The peaking component of the cross-feed background contributes just a few events.

Fits to the $`m_{\mathrm{ES}}`$ distributions of data events are performed in each of the nine $`M_{Had}`$ bins using a Crystal Ball function for the contributions that peak in the signal region and two ARGUS functions for the non-peaking backgrounds. The plots of the fits to the eight $`M_{Had}`$ bins from 0.8 to 2.4 $`\mathrm{GeV}/c^2`$ can be seen in Figure 2.

Table 1 lists the peaking yields for the sum of all final states in data, $`B\overline{B}`$ Monte Carlo and cross–feed, normalized to our measured inclusive branching fraction.

## 5 Signal Efficiency

The signal efficiency is determined from Monte Carlo samples using the yield for true signal events after cross-feed has been removed. In the first two $`M_{Had}`$ bins, 0.6-1.0 $`\mathrm{GeV}/c^2`$, the signal is modeled by the four exclusive $`K^{}\gamma `$ modes. In the range 1.0–2.4 $`\mathrm{GeV}/c^2`$, we use the inclusive $`X_s\gamma `$ model of Kagan and Neubert with $`m_b=4.65`$ $`\mathrm{GeV}/c^2`$. This inclusive model uses non-resonant JETSET  fragmentation of the $`X_S`$ to produce the hadronic final states. Changing the modeling of the 1.0-1.2 $`\mathrm{GeV}/c^2`$ bin to $`K^{}\gamma `$, and varying the parameters of the inclusive model form part of our systematic studies.

The overall efficiency in each $`M_{Had}`$ bin is an average over the individual final states, weighted by the fractional contribution of that final state to the total in that $`M_{Had}`$ bin. These efficiencies have to be corrected by $`10`$% for the differences between data and Monte Carlo detection efficiencies for $`\pi ^+`$, $`K^+`$, $`\pi ^0`$, $`\gamma `$ and $`K_S`$. There are significant differences between the efficiencies for the individual final states and a strong dependence of the overall efficiency on $`M_{Had}`$. To illustrate this, the efficiencies can be understood in terms of a 1/$`M_{Had}`$ dependence of the efficiency for each final state, a 60% efficiency ratio between $`\pi ^0`$ and $`\pi ^+`$ final states, a 55% efficiency ratio between $`K_S`$ and $`K^+`$ final states, and a 60% efficiency factor between 3 and 2–body, and between 4 and 3–body hadronic systems.

To check the fractions of each final state in the non-resonant Monte Carlo, we divide the sample above $`M_{Had}=1.0\mathrm{GeV}/c^2`$ into three sets of final state categories: with and without a $`\pi ^0`$, with and without a $`K_S`$, and 2, 3, or 4–body hadronic systems. A check of the $`K_S/K^+`$ samples shows that the data and Monte Carlo values agree, and the ratio $`K_S/K^+`$ is consistent with 0.5 as expected by isospin arguments. However, we find that the fraction of $`\pi ^0`$ final states needs to be increased by a factor 1.5 and the fraction of 2–body final states decreased by a factor 0.4 in order to obtain a reasonable agreement between non-resonant Monte Carlo simulation and data. When these adjustments are applied to the fractions of final states generated by JETSET, they lead to a correction of $`15`$% in the overall efficiency. Table 2 presents the signal efficiencies after all corrections have been applied.

## 6 Branching Fractions

Partial branching fractions ($`\mathrm{\Delta }`$) are calculated for each $`M_{Had}`$ bin for the sum of the twelve final states reconstructed in this analysis. The calculation uses the data and $`B\overline{B}`$ yields from Table 1 and the signal efficiencies from Table 2. The peaking cross–feed background has to be normalized to our measured inclusive branching fraction ($``$). We take care of this dependence of the peaking cross–feed background by including effective terms that modify the signal efficiencies in each bin, and report the equivalent peaking cross–feed yields in Table 1. Figure 3 shows the data yields after subtraction of the peaking $`B\overline{B}`$ and cross–feed backgrounds (signal yields).

The inclusive branching fraction for $`bs\gamma `$ in each $`M_{Had}`$ bin can be obtained from the partial branching fractions by taking into account the fraction, $`f`$, of final states that are considered in the analysis. The $`f`$s are taken from the generator level JETSET fragmentation of the inclusive $`bs\gamma `$ signal. Since there are discrepancies in the modeling of our twelve reconstructed final states as already noted above, we know that the values of $`f`$ are not completely reliable. There are also uncertainties coming from other final states, e.g. with $`\eta `$ and $`\varphi `$ mesons, which have not been looked for and therefore cannot be compared with data.

We divide up the missing final states from our generator–level Monte Carlo into categories with a $`K_L`$ (25%), more than one $`\pi ^0`$ (4-18%), more than 4–body hadronic systems (0-14%), an $`\eta /\eta ^{}`$ (5-11%), a $`s\overline{s}s`$ system (0-7%), and “other”(0-4%), where the missing fractions increase from $`M_{Had}=1.0`$ to $`M_{Had}=2.4`$ $`\mathrm{GeV}/c^2`$. The $`K_L`$ fraction is known to be 25% from isospin and the measured $`K_S/K^+`$ ratio in data. For all the other fractions we allow a variation between 0.5 and 1.5 times the generator fraction, and add these variations in quadrature to give an overall systematic uncertainty on the missing fractions. Table 3 summarizes the partial branching fractions with their statistical errors, the fraction of studied final states with the corresponding errors, the resulting branching fractions with their statistical and systematic errors and the cumulative $`bs\gamma `$ branching fraction, with statistical and systematic errors, in bins of $`M_{Had}`$.

The systematic errors will be described in Section 7. It can be seen that the branching fraction for $`M_{Had}<1.0`$ $`\mathrm{GeV}/c^2`$ agrees with the exclusive $`K^{}\gamma `$ branching fraction of about $`410^5`$ , and that the integral for $`M_{Had}>2.0`$ $`\mathrm{GeV}/c^2`$ is higher than previous measurements of the inclusive $`bs\gamma `$ branching fraction , although with large uncertainties. Note that these results are essentially independent of the modeling of the $`bs\gamma `$ spectrum.

Figure 4 shows the branching fractions as a function of $`M_{Had}`$. The same results are shown as a function of $`E_\gamma `$ in Figure 5.

## 7 Systematic Errors

The first category of systematic errors (general systematic) are those which do not depend on the hadronic final state. The systematic error on the detection of the high–energy photon is 3.7%, and the systematic error on the total number of $`B`$ mesons in our data sample is 1.6%. There is also a systematic uncertainty of 1% coming from the use of the $`\mathrm{cos}\theta _T^{}`$ and Fisher Discriminant requirements to suppress continuum backgrounds. These errors give a total contribution of 4.2% independent of the hadronic mass bin.

A second category of errors are due to the fitting procedure and signal definition (fit systematic). These errors are evaluated by varying the fixed parameters used in the fits, and by varying the fit procedures. The width of the signal Crystal Ball function is varied by $`\pm `$0.2 $`\mathrm{MeV}/c^2`$ for both signal and peaking backgrounds. This introduces a systematic uncertainty of 5% on the signal yield. Varying the Crystal Ball shape to reflect our understanding of the radiative tail introduces a further systematic uncertainty of 5%, and varying the peak position introduces a systematic uncertainty of 3%. The ARGUS shape parameters are varied within the ranges allowed by the Monte Carlo samples. They give a contribution to the systematic uncertainty of 10%, which is limited by the Monte Carlo statistics. We also include a systematic error of 20% on the peaking $`B\overline{B}`$ yield in the region $`M_{Had}`$= 1.8-2.4 $`\mathrm{GeV}/c^2`$ due to the method of fitting the $`B\overline{B}`$ yield, and to uncertainties in the content of the generic $`B\overline{B}`$ Monte Carlo simulation. This translates into a 10% systematic error in the last three $`M_{Had}`$ bins.

Next we consider errors in the detection efficiencies for the hadronic final states (efficiency systematic). We assign systematic errors based on the uncertainties in the corrections to the efficiencies that we have applied. The uncertainties coming from the single particle detection differences between data and Monte Carlo are between 4.0 and 4.6% as a function of $`M_{Had}`$, and the fragmentation differences give uncertainties between 2.1 and 6.3%.

A final category of errors is due to the generator and modeling we use (generator systematic). The errors due to the missing final states not considered in the analysis were discussed in the previous section and shown in Table 3. Also, we consider contributions due to changing the modeling of the bin 1.0-1.2 $`\mathrm{GeV}/c^2`$ to $`K^{}\gamma `$ and varying the parameters of the inclusive model in our signal Monte Carlo.

The contributions to the systematic errors coming from all sources are summarized in Table 4.

## 8 Fits to the Spectrum

So far our results have had no significant dependence on the $`bs\gamma `$ model. We now fit the hadronic mass spectrum with the shape predicted by the Kagan and Neubert model . This extrapolates from the measured range of $`M_{Had}`$ to give the inclusive branching fraction for $`bs\gamma `$.

In the fits we vary the parameters $`m_b`$ and $`\lambda _1`$, the transition point between the $`K^{}`$ and non-resonant contributions, and the total normalization, which corresponds to our integrated branching fraction. In the plane $`\lambda _1`$ and $`m_b`$, Figure 6 shows the value of the $`\chi ^2`$, having fit the transition point (resulting around 1.1 $`\mathrm{GeV}/c^2`$) and the normalization for each point. There is a strong correlation between these parameters, which is expected in this model, and a shallow minimum at $`m_b=4.65`$ $`\mathrm{GeV}/c^2`$ and $`\lambda _1=0.48[\mathrm{GeV}/\mathrm{c}^2]^2`$. Figure 7(a) shows the correlation curve between $`\lambda _1`$ and $`m_b`$ along which we have the minimum value of $`\chi ^2`$, and Figure 7(b) shows the values of $`\chi ^2`$ as a function of $`m_b`$. Figure 7(c) shows the inclusive branching fraction and the branching fraction for $`M_{Had}`$ $`<`$ 2.4 $`\mathrm{GeV}/c^2`$ as a function of $`m_b`$, where the result for the minimum value of $`m_b`$ found in the free fit is marked with a line. Since our spectrum contains a $`K^{}`$ peak with a known branching fraction of about $`410^5`$, a reduction in $`m_b`$ leads to an increase in the non-resonant component of the spectrum, and hence to a higher inclusive branching fraction. This increase can be seen both in the branching fraction restricted to the measured range, and in the extrapolation over all $`M_{Had}`$.

The free fit to the model parameters only allows us to set an upper bound on $`m_b`$, since the $`\chi ^2`$ rises rather slowly on the other side of the minimum. We obtain a lower bound on $`m_b`$ by measuring the first moment of the $`E_\gamma `$ distribution, $`<E_\gamma >`$, and then extracting $`\overline{\mathrm{\Lambda }}(\alpha _s^2,1/M_B^3)`$ using HQET calculations  with the coefficients relevant for our minimum value of $`E_\gamma =2.094`$ $`\mathrm{GeV}`$ . Then, we truncate the HQET expressions to order $`O(\alpha _s,1/M_B^2)`$ and recompute $`m_b=M_B\overline{\mathrm{\Lambda }}(\alpha _s,1/M_B^2)`$ to obtain the parameter required in the Kagan and Neubert model. There are significant theoretical uncertainties associated with this method coming from varying the scale at which $`\alpha _s`$ is determined between $`m_b/2`$ and $`2m_b`$, from the truncation of higher–order terms, and from the dependence of $`\overline{\mathrm{\Lambda }}`$ on the choice of the minimum $`E_\gamma `$ .

The value of the first moment, $`<E_\gamma >|_{E_\gamma >2.094\mathrm{GeV}}`$, is:

$`<E_\gamma >|_{E_\gamma >2.094\mathrm{GeV}}=2.35\pm 0.04(stat)\pm 0.04(syst)`$ $`\mathrm{GeV}`$,

which gives:

$`\overline{\mathrm{\Lambda }}=0.37\pm 0.09(stat)\pm 0.07(syst)\pm 0.10(model)`$ $`\mathrm{GeV}/c^2`$.

where the theoretical uncertainties have been included $`(model)`$. With the truncation to first order in $`\alpha _s`$, $`\overline{\mathrm{\Lambda }}=0.49`$ $`\mathrm{GeV}/c^2`$, which gives:

$`m_b=4.79\pm 0.08(stat)\pm 0.10(syst)\pm 0.10(model)`$ $`\mathrm{GeV}/c^2`$.

## 9 Final Results

We perform final fits to the hadronic mass spectrum with $`m_b=4.79\pm 0.16\mathrm{GeV}/c^2`$ taken from the first moment analysis. These fits give the inclusive branching fraction:

$`(bs\gamma )=4.3\pm 0.5(stat)\pm 0.8(syst)\pm 1.3(model)10^4`$

and a range for the parameter $`\lambda _1`$:

$`\lambda _1=0.24_{0.04}^{+0.03}(stat)\pm 0.02(syst)_{0.21}^{+0.15}(model)[\mathrm{GeV}/c^2]^2`$.

Figure 8 shows the superposition of the theoretical spectrum on the data distribution for the best fits with $`m_b=4.79`$ $`\mathrm{GeV}/c^2`$ and $`m_b=4.65`$ $`\mathrm{GeV}/c^2`$.

Our values for $`\overline{\mathrm{\Lambda }}`$, at $`O(\alpha _s^2,1/M_B^3)`$, and $`\lambda _1`$, at $`O(\alpha _s,1/M_B^2)`$, are consistent with other results , but our inclusive branching fraction is somewhat higher . Our error bars are still larger than the previous measurements, but they can be reduced significantly with the addition of more data. It is possible to reduce the systematic errors by improvements in the treatment of backgrounds, and by studying a larger number of final states to reduce the missing fraction errors. The model dependence will also be reduced by fitting the spectrum to determine the model parameters.

## 10 Acknowledgments

We would like to thank Alexander L. Kagan, Zoltan Ligeti and Mark B. Wise, for useful discussions and for having provided us with numbers and results from their works. We are grateful for the extraordinary contributions of our PEP-II colleagues in achieving the excellent luminosity and machine conditions that have made this work possible. The success of this project also relies critically on the expertise and dedication of the computing organizations that support BABAR. The collaborating institutions wish to thank SLAC for its support and the kind hospitality extended to them. This work is supported by the US Department of Energy and National Science Foundation, the Natural Sciences and Engineering Research Council (Canada), Institute of High Energy Physics (China), the Commissariat à l’Energie Atomique and Institut National de Physique Nucléaire et de Physique des Particules (France), the Bundesministerium für Bildung und Forschung and Deutsche Forschungsgemeinschaft (Germany), the Istituto Nazionale di Fisica Nucleare (Italy), the Research Council of Norway, the Ministry of Science and Technology of the Russian Federation, and the Particle Physics and Astronomy Research Council (United Kingdom). Individuals have received support from the A. P. Sloan Foundation, the Research Corporation, and the Alexander von Humboldt Foundation.


text:
# To the problem of Poincaré recurrences in generic Hamiltonian systems

\[

## Abstract

We discuss the problem of Poincaré recurrences in area-preserving maps and the universality of their decay at long times. The work is related to to the results presented in Refs. .

\]

The authors of Comment  question the existence of universal power-law decay of Poincaré recurrences discussed in our Letter  for generic Hamiltonian systems. There are two aspects in their quarry:

The first one states that the asymptotic decay $`P(\tau )1/\tau ^3`$, which is expected from the scaling theory of universal phase-space structure in the vicinity of critical golden curve (see Refs. 1-2, 12-18 in ), is not valid. Such a conclusion is based on Fig. 1 in  obtained from numerical simulations of the standard map with the total number of map iterations $`N_{tot}`$ being only by a factor 10 larger than for the similar case of Fig. 2 in  where $`N_{tot}=10^{12}`$. New data for $`8<\mathrm{log}\tau 9`$ indeed demonstrate noticeable deviations from the theoretical estimates for $`P(\tau )`$ obtained from the exit times $`\tau _n`$ from golden resonance scales $`r_n`$ (Fig. 1 in ). These deviations are rather intriguing. Indeed, the local properties of critical golden curve are known to be self-similar and universal (e.g. phase-space structure, size of resonances and local diffusion rate, see e.g. Refs. 2, 15 in ). Moreover, the exit times $`\tau _n`$, found in  for a first time, give an example of a non-local characteristic which follows the universal scaling law. The arguments based on this scaling lead to the asymptotic decay $`P(\tau )1/\tau ^3`$ which however can start after extremely long times $`\tau >\tau _a`$. In  we have shown that at least $`\tau _a>\tau _g2\times 10^5`$ and it is not excluded that this time scale is still much longer , and is not yet reached even in the simulations presented in .

If to assume this then there is an interesting possibility to see if $`P(\tau )`$ would have some universal properties on the presently available (intermediate) time scales ($`\tau 10^9`$). This is the second aspect of the quarry on which the authors of the Comment tend to give a negative answer. To this end we show in Fig. 1 an example of another map with critical golden curve which power-law decay $`P(\tau )1/\tau ^p`$ has the exponent $`p1.5`$. This average value of exponent $`p`$ in the range $`10^2<\tau <10^8`$ agrees, indeed, with that of the decay in the standard map in the range $`10^5<\tau <10^9`$ shown in . There is also a clear similarity for the decay of $`P(\tau )`$ in two maps (see Fig. 1). For other values of parameter $`\lambda `$ the decay $`P(\tau )`$ in the separatrix map is no longer a simple power law but rather some irregular oscillations around that which represent a peculiar phenomenon of the so-called renormalization chaos (Ref. 12 in ). However, averaging over several $`\lambda `$ values smoothes away the signs of renormalization chaos and reveals the underlying picture of the power-law distribution with the same exponent $`p1.5`$ (see Ref. 9 in  and ). Certainly, further studies are required for the understanding of the origin of this intermediate asymptotics and for the estimates of the time scale $`\tau _a`$ after which a transition to $`p=3`$ is expected.