text:
# References

UCRHEP-T342

August 2002

The All-Purpose Neutrino Mass Matrix
Ernest Ma
Physics Department, University of California, Riverside, California 92521

## Abstract

A four-parameter Majorana neutrino mass matrix is proposed, which is exactly diagonalizable with two mixing angles, one of which is $`\pi /4`$ and the other large but less than $`\pi /4`$. The form of this mass matrix allows seven realizations of the present data on atmospheric and solar neutrino oscillations: three have the normal hierarchy, two have the inverse hierarchy, and two have three nearly degenerate masses. The possible inclusion of CP violation is also discussed.

With the recent addition of the SNO (Sudbury Neutrino Observatory) neutral-current data , the overall picture of solar neutrino oscillations  is becoming quite clear. Together with the well-established atmospheric neutrino data , the neutrino mixing matrix is now determined to a very good first approximation by

$$\left(\begin{array}{c}\nu _1\\ \nu _2\\ \nu _3\end{array}\right)=\left(\begin{array}{ccc}\mathrm{cos}\theta & \mathrm{sin}\theta /\sqrt{2}& \mathrm{sin}\theta /\sqrt{2}\\ \mathrm{sin}\theta & \mathrm{cos}\theta /\sqrt{2}& \mathrm{cos}\theta /\sqrt{2}\\ 0& 1/\sqrt{2}& 1/\sqrt{2}\end{array}\right)\left(\begin{array}{c}\nu _e\\ \nu _\mu \\ \nu _\tau \end{array}\right),$$
(1)
where $`\nu _{1,2,3}`$ are Majorana neutrino mass eigenstates. In the above, $`\mathrm{sin}^22\theta _{atm}=1`$ is already assumed and $`\theta `$ is the solar mixing angle which is now known to be large but not maximal , i.e. $`\mathrm{tan}^2\theta 0.4`$. The $`U_{e3}`$ entry has been assumed zero but it is only required to be small , i.e. $`|U_{e3}|<0.16`$. A nonzero $`U_{e3}`$ with a possible complex phase for CP violation will be discussed later.

Denoting the masses of $`\nu _{1,2,3}`$ as $`m_{1,2,3}`$, the solar neutrino data  require that $`m_2^2>m_1^2`$ with $`\theta <\pi /4`$, and in the case of the favored large-mixing-angle solution ,

$$\mathrm{\Delta }m_{sol}^2=m_2^2m_1^25\times 10^5\mathrm{eV}^2.$$
(2)
The atmospheric neutrino data  require

$$|m_3^2m_{1,2}^2|2.5\times 10^3\mathrm{eV}^2,$$
(3)
without deciding whether $`m_3^2>m_{1,2}^2`$ or $`m_3^2<m_{1,2}^2`$.

The big question now is what the neutrino mass matrix itself should look like. It is proposed here that it should be of the form

$$_\nu =\left(\begin{array}{ccc}a+2b+2c& d& d\\ d& b& a+b\\ d& a+b& b\end{array}\right).$$
(4)
This matrix has 4 parameters which are assumed first to be all real. In that case, it is $`exactly`$ diagonalized by Eq. (1) with its eigenvalues $`exactly`$ given by (I):

$`m_1`$ $`=`$ $`a+2b+c\sqrt{c^2+2d^2},`$ (5)
$`m_2`$ $`=`$ $`a+2b+c+\sqrt{c^2+2d^2},`$ (6)
$`m_3`$ $`=`$ $`a,`$ (7)
with

$$\mathrm{tan}\theta =\frac{\sqrt{2}d}{\sqrt{c^2+2d^2}c},$$
(8)
for $`a+2b+c>0`$, $`c<0`$, and (II):

$`m_1`$ $`=`$ $`a+2b+c+\sqrt{c^2+2d^2},`$ (9)
$`m_2`$ $`=`$ $`a+2b+c\sqrt{c^2+2d^2},`$ (10)
$`m_3`$ $`=`$ $`a,`$ (11)
with

$$\mathrm{tan}\theta =\frac{\sqrt{c^2+2d^2}c}{\sqrt{2}d},$$
(12)
for $`a+2b+c<0`$, $`c>0`$.

Note that $`\theta `$ depends only on the ratio $`d/c`$, which must be of order unity. This shows the advantage for adopting the parametrization of Eq. (4). The constraints of Eqs. (2) and (3) are then realized by the following 7 different conditions on $`a`$, $`b`$, and $`c`$.

(1) $`||a+2b+c|\sqrt{c^2+2d^2}|<<|a+2b+c|<<|a|`$, i.e. $`|m_1|<<|m_2|<<|m_3|`$.

(2) $`\sqrt{c^2+2d^2}<<|a+2b+c|<<|a|`$, i.e. $`|m_1||m_2|<<|m_3|`$.

(3) $`|a+2b+c|<<\sqrt{c^2+2d^2}<<|a|`$, i.e. $`|m_1||m_2|<<|m_3|`$.

(4) $`|a|,\sqrt{c^2+2d^2}<<|a+2b+c|`$, i.e. $`|m_3|<<|m_1||m_2|`$.

(5) $`|a|,|a+2b+c|<<\sqrt{c^2+2d^2}`$, i.e. $`|m_3|<<|m_1||m_2|`$.

(6) $`\sqrt{c^2+2d^2}<<||a+2b+c||a||<<|a|`$, i.e. $`|m_1||m_2||m_3|`$.

(7) $`|a+2b+c|<<\sqrt{c^2+2d^2}|a|`$, i.e. $`|m_1||m_2||m_3|`$.

Cases (1) to (3) have the normal hierarchy. Cases (4) and (5) have the inverse hierarchy. Cases (6) and (7) have 3 nearly degenerate masses. The versatility of Eq. (4) has clearly been demonstrated.

The above 7 cases encompass all models of the neutrino mass matrix that have ever been proposed which also satisfy Eq. (1). They are also very useful for discussing the possibility of neutrinoless double beta ($`\beta \beta _{0\nu }`$) decay in the context of neutrino oscillations . The effective mass $`m_0`$ measured in $`\beta \beta _{0\nu }`$ decay is $`|a+2b+2c|`$. However, neutrino oscillations constrain $`|a+2b+c|`$ and $`\sqrt{c^2+2d^2}`$, as well as $`|d/c|`$. Using

$$|a+2b+2c|=||a+2b+c|\pm |c||=||a+2b+c|\pm \mathrm{cos}2\theta \sqrt{c^2+2d^2}|,$$
(13)
the following conditions on $`m_0`$ are easily obtained:

$`(1)`$ $`m_0\mathrm{sin}^2\theta |m_2|\mathrm{sin}^2\theta \sqrt{\mathrm{\Delta }m_{sol}^2},`$ (14)
$`(2)`$ $`m_0|m_{1,2}|<<\sqrt{\mathrm{\Delta }m_{atm}^2},`$ (15)
$`(3)`$ $`m_0\mathrm{cos}2\theta |m_{1,2}|<<\mathrm{cos}2\theta \sqrt{\mathrm{\Delta }m_{atm}^2},`$ (16)
$`(4)`$ $`m_0\sqrt{\mathrm{\Delta }m_{atm}^2},`$ (17)
$`(5)`$ $`m_0\mathrm{cos}2\theta \sqrt{\mathrm{\Delta }m_{atm}^2},`$ (18)
$`(6)`$ $`m_0|m_{1,2,3}|,`$ (19)
$`(7)`$ $`m_0\mathrm{cos}2\theta |m_{1,2,3}|.`$ (20)
If $`m_0`$ is measured  to be significantly larger than 0.05 eV, then only Cases (6) and (7) are allowed. However, as Eqs. (19) and (20) show, the true mass of the three neutrinos is still subject to a two-fold ambiguity, which is a well-known result.

The underlying symmetry of Eq. (4) which results in $`U_{e3}=0`$ is its invariance under the interchange of $`\nu _\mu `$ and $`\nu _\tau `$. Its mass eigenstates are then separated according to whether they are even $`(\nu _{1,2})`$ or odd $`(\nu _3)`$ under this interchange, as shown by Eq. (1). To obtain $`U_{e3}0`$, this symmetry has to be broken. One interesting possibility is to rewrite Eq. (4) as

$$_\nu =\left(\begin{array}{ccc}a+2b+2c& d& d^{}\\ d& b& a+b\\ d^{}& a+b& b\end{array}\right),$$
(21)
where $`a,b,c`$ are real but $`d`$ is complex. This reduces to Eq. (4) if $`Imd=0`$, but if $`Imd0`$, then $`U_{e3}0`$.

To obtain $`U_{e3}`$ in a general way, first rotate to the basis spanned by $`\nu _e,(\nu _\mu +\nu _\tau )/\sqrt{2}`$, and $`(\nu _\tau \nu _\mu )/\sqrt{2}`$, i.e.

$$_\nu =\left(\begin{array}{ccc}a+2b+2c& \sqrt{2}Red& \sqrt{2}iImd\\ \sqrt{2}Red& a+2b& 0\\ \sqrt{2}iImd& 0& a\end{array}\right)$$
(22)
Whereas $`_\nu `$ is diagonalized by

$$U_\nu U^T=\left(\begin{array}{ccc}m_1& 0& 0\\ 0& m_2& 0\\ 0& 0& m_3\end{array}\right),$$
(23)
$`_\nu _\nu ^{}`$ is diagonalized by

$$U(_\nu _\nu ^{})U^{}=\left(\begin{array}{ccc}|m_1|^2& 0& 0\\ 0& |m_2|^2& 0\\ 0& 0& |m_3|^2\end{array}\right).$$
(24)
Here

$$_\nu _\nu ^{}=\left(\begin{array}{ccc}(a+2b+2c)^2+2|d|^2& 2\sqrt{2}(a+2b+c)Red& 2\sqrt{2}i(a+b+c)Imd\\ 2\sqrt{2}(a+2b+c)Red& (a+2b)^2+2(Red)^2& 2iRedImd\\ 2\sqrt{2}i(a+b+c)Imd& 2iRedImd& a^2+2(Imd)^2\end{array}\right).$$
(25)
To obtain $`U_{e3}`$ for small $`Imd`$, consider the matrix

$$A=_\nu _\nu ^{}[a^2+2(Imd)^2]I,$$
(26)
where $`I`$ is the identity matrix. Now $`A`$ is diagonalized by $`U`$ as well and $`U_{e3}`$ is simply given by

$$U_{e3}\frac{A_{e3}}{A_{ee}}=\frac{2\sqrt{2}i(a+b+c)Imd}{(a+2b+2c)^2a^2+2(Red)^2}$$
(27)
to a very good approximation and leads to

$`(1),(2),(3)`$ $`U_{e3}{\displaystyle \frac{\sqrt{2}iImd}{a}},`$ (28)
$`(4),(6)`$ $`U_{e3}{\displaystyle \frac{iImd}{\sqrt{2}b}},`$ (29)
$`(5)`$ $`U_{e3}{\displaystyle \frac{\sqrt{2}iImd}{c}},`$ (30)
$`(7)`$ $`U_{e3}{\displaystyle \frac{\sqrt{2}i(a+c)Imd}{c^2a^2+2(Red)^2}}.`$ (31)
In all cases, the magnitude of $`U_{e3}`$ can be as large as the present experimental limit  of 0.16 and its phase is $`\pm \pi /2`$. Thus the CP violating effect in neutrino oscillations is predicted to be maximal by Eq. (21), which is a very desirable scenario for future long-baseline neutrino experiments.

The above analysis shows that for $`U_{e3}=0`$ and $`\mathrm{sin}^22\theta _{atm}=1`$, the seven cases considered cover all possible patterns of the $`3\times 3`$ Majorana neutrino mass matrix, as indicated by present atmospheric and solar neutrino data. Any successful model should predict Eq. (4) at least as a first approximation. One such example already exists , where $`b=c=d=0`$ corresponds to the non-Abelian discrete symmetry $`A_4`$, i.e. the finite group of the rotations of a regular tetrahedron. This leads to Case (6), i.e. three nearly degenerate masses, with the common mass equal to that measured in $`\beta \beta _{0\nu }`$ decay. It has also been shown recently  that starting with this pattern, the correct mass matrix, i.e. Eq. (21) with the complex phase in the right place, is automatically obtained with the most general application of radiative corrections. In particular, if soft supersymmetry breaking is assumed to be the origin of these radiative corrections, then the neutrino mass matrix is correlated with flavor violation in the slepton sector, and may be tested in future collider experiments.

If $`m_0`$ turns out to be of order $`\sqrt{\mathrm{\Delta }m_{atm}^2}`$ or less, then Cases (6) and (7) would be ruled out, and further progress in understanding the pattern of neutrino oscillations will depend on having a rationale for choosing one of the remaining 5 cases. Consider for example Case (1) with $`a+2b+c=\sqrt{c^2+2d^2}`$ and $`c=d/2<0`$, then Eq. (4) becomes

$$_\nu =\left(\begin{array}{ccc}d& d& d\\ d& b& b+2d\\ d& b+2d& b\end{array}\right),$$
(32)
which results in $`m_1=0`$, $`m_2=3d`$, $`m_3=2b2d`$, and $`\mathrm{tan}^2\theta =0.5`$.

As another example, consider Case (4) with $`a=0`$, $`c=d/2<0`$, then Eq. (4) becomes

$$_\nu =\left(\begin{array}{ccc}2bd& d& d\\ d& b& b\\ d& b& b\end{array}\right),$$
(33)
which results in $`m_1=2b2d`$, $`m_2=2b+d`$, $`m_3=0`$, and $`\mathrm{tan}^2\theta =0.5`$. Thus $`m_2^2m_1^212db`$ and $`m_{1,2}^2m_3^24b^2`$, where $`d<<b`$ has been used.

Both examples are simple acceptable two-parameter solutions of the present atmospheric and solar neutrino data. This shows the utility of Eq. (4) as a way to find possible patterns which may even turn out to be supported by some symmetry. As further examples, the case $`c=0`$ was discussed in Ref.  and the case $`a+2b+c=0`$ was discussed in Ref. .

In conclusion, recent experimental progress on neutrino oscillations points to a neutrino mixing matrix which can be understood in a systematic way in terms of an all-purpose neutrino mass matrix, i.e. Eq. (4), and its simple extension, i.e. Eq. (21), to allow for a nonzero and $`imaginary`$ $`U_{e3}`$, i.e. Eq. (27). Seven possible cases have been identified, each with a different prediction for $`\beta \beta _{0\nu }`$ decay, i.e. Eqs. (14) to (20). The present analysis should be useful as both a checklist and a roadmap for the further theoretical understanding of the origin of neutrino masses and mixing.

This work was supported in part by the U. S. Department of Energy under Grant No. DE-FG03-94ER40837.


format: text

text:
# The Probability Distribution of Binary Pulsar Coalescence Rates. I. Double Neutron Star Systems in the Galactic Field

## 1 INTRODUCTION

The detection of the double neutron star (NS–NS) prototype PSR B1913+16 as a binary pulsar Hulse & Taylor (1975) and its orbital decay due to emission of gravitational waves Taylor et al. (1979) has inspired a number of quantitative estimates of the coalescence rate of NS–NS binaries Clark et al. (1979); Narayan et al. (1991); Phinney (1991); Curran & Lorimer (1995). In general, the coalescence rate of NS–NS binaries can be calculated based on: (a) our theoretical understanding of their formation (see Belczynski & Kalogera 2001 for a review and application of this approach); (b) the observational properties of the pulsars in the binary systems and the modeling of pulsar survey selection effects (see e.g. Narayan 1987). Interest in these coalescence derives from an intrinsic motivation of understanding their origin and evolution and their connections to other NS binaries. However, significant interest derives from their importance as gravitational-wave sources for the upcoming ground-based laser interferometers (such as LIGO) and their possible association with $`\gamma `$-ray burst events (Popham et al. 1998 and references therein).

The traditional way of calculating the coalescence rate based on observations involves an estimate of the scale factor, an indicator for the number of pulsars in our Galaxy with the same spin period and luminosity Narayan (1987). Corrections must then be applied to these scale factors to account for the faint end of the pulsar luminosity function, the beamed nature of pulsar emission, and uncertainties in the assumed spatial distribution. The estimated total number in the Galaxy can then be combined with estimates of their lifetimes to obtain a coalescence rate, $``$. This method was first applied by Narayan et al. (1991) and Phinney (1991) and other investigators who followed Curran & Lorimer (1995); van den Heuvel & Lorimer (1996). Various correction factors were (or were not) included at various levels of completeness. Summaries of these earlier studies can be found in Arzoumanian et al. (1999) and Kalogera et al. (2001; hereafter KNST). The latter authors examined all possible uncertainties in the estimates of the coalescence rate of NS–NS binaries in detail, and pointed a small-number bias that introduces a large uncertainty (more than two orders of magnitude) in the correction factor for the faint-pulsar population that must be applied to the rate estimate. They obtained a total NS–NS rate estimate in the range $`=10^65\times 10^4`$ yr-1, with the uncertainty dominated by the small-number bias. Earlier studies, which made different assumptions about the pulsar properties (e.g. luminosity and spatial distributions and lifetimes), are roughly consistent with each other (given the large uncertainties). Estimated ranges of values until now were not associated with statistical significance statements and an “all-inclusive” estimated Galactic coalescence rate lies in the range $`10^710^5`$ yr-1.

The motivation for this paper is to update the scale factor calculations using the most recent pulsar surveys, and present a statistical analysis that allows the calculation of statistical confidence levels associated with rate estimates. We consider the two binaries found in the Galactic disk: PSR B1913+16 Hulse & Taylor (1975) and PSR B1534+12 Wolszczan (1991). Following the arguments made by Phinney (1991) and KNST, we do not include the globular cluster system PSR B2127+11C Prince et al. (1991); this system will be the subject of a later paper. Radio-pulsar-survey selection effects are taken into account in the modeling of pulsar population. As described in what follows, the small-number bias and the effect of a luminosity function are implicitly included in our analysis, and therefore a separate correction factor is not needed. For each population model of pulsars, we derive the probability distribution function of the total Galactic coalescence rate weighted by the two observed binary systems. In our results we note a number of important correlations between $`_{\mathrm{peak}}`$ and model parameters that are useful in generalizing the method. We extrapolate the Galactic rate to cover the detection volume of LIGO and estimate the detection rates of NS–NS inspiral events for the initial and advanced LIGO.

The plan for the rest of this paper is as follows. In §2, we describe our analysis method in a qualitative way. Full details of the various pulsar population models and survey selection effects are then given in §3 and §4 respectively. In §5, we derive the probability distribution function for the total Galactic coalescence rate and calculate the detection rate of LIGO. In §6, we summarize our results and discuss a number of intriguing correlations between various physical quantities. Finally, in §7, we discuss the results and compare them with previous studies.

## 2 BASIC ANALYSIS METHOD

Our basic method is one of “forward” analysis. By this we mean that we do not attempt to invert the observations to obtain the total number of NS–NS binaries in the Galaxy. Instead, using Monte Carlo methods, we populate a model galaxy with NS–NS binaries (that match the spin properties of PSR B1913+16 and PSR B1534+12) with pre-set properties in terms of their spatial distribution and radio pulsar luminosity function. Details about these “physical models” are given in §3.

For a given physical model, we produce synthetic populations of different total numbers of objects ($`N_{\mathrm{tot}}`$). We then produce a very large number of Monte Carlo realizations of such pulsar populations and determine the number of objects ($`N_{\mathrm{obs}}`$) that are observable by all large-scale pulsar surveys carried out to date by detailed modeling of the detection thresholds of these surveys. This analysis utilizes code to take account of observational selection effects in a self-consistent manner, developed and described in detail by Lorimer et al. (1993; hereafter LBDH) which we summarize in §4. Performing this analysis for many different Monte Carlo realizations of the physical model allows us to examine the distribution of $`N_{\mathrm{obs}}`$. We find, as expected and assumed by other studies, that this distribution closely follows Poisson statistics, and we determine the best-fit value of the mean of the Poisson distribution $`\lambda `$ for each population model and value of $`N_{\mathrm{tot}}`$.

The calculations described so far are performed separately for each of PSR B1913+16 and PSR B1534+12 so that we obtain separate best-fit $`\lambda `$ values for the Poisson distributions. Doing the analysis in this way allows us to calculate the likelihood of observing just one example of each pulsar in the real-world sample. Given the Poissonian nature of the distributions this likelihood is simply: $`P(1;\lambda )=\lambda \mathrm{exp}(\lambda )`$. We then calculate this likelihood for a variety of assumed $`N_{\mathrm{tot}}`$ values for each physical model.

The probability distribution of the total coalescence rate $`_{\mathrm{tot}}`$ is derived using the Bayesian analysis and the calculated likelihood for each pulsar (described in detail in §5). The derivation of this probability distribution allows us to calculate the most probable rate as well as determine its ranges of values at various statistical confidence levels. Finally, we extrapolate the Galactic rate to the volume expected to be reached by LIGO and calculate the detection rate, $`_{\mathrm{det}}`$ (see §6).

## 3 MODELS FOR THE GALACTIC PULSAR POPULATION

Our model pulsar populations are characterized by a Galactocentric radius ($`R`$), vertical distance ($`Z`$) from the Galactic plane and radio luminosity ($`L`$). Assuming that the distributions of each of these parameters are independent, the combined probability density function (PDF) of the model pulsar population can be written as:

$$f(R,Z,L)dRdZdL=\psi _\mathrm{R}(R)\mathrm{\hspace{0.17em}2}\pi RdR\psi _\mathrm{Z}(Z)dZ\varphi (L)dL,$$
(1)
where $`\psi _R(R)`$, $`\psi _Z(Z)`$ and $`\varphi (L)`$, are the individual PDFs of $`R`$, $`Z`$ and $`L`$, respectively. In all models considered, we assume azimuthal symmetry about the Galactic center.

The spatial distribution of pulsars is rather loosely constrained, but we find that it does not affect the results significantly for a wide range of models. For the radial and the vertical PDFs, we consider Gaussian and exponential forms with different values of the radial $`R_0`$ and the vertical $`Z_0`$ scale. In our reference model, we assume a Gaussian PDF for the radial component and an exponential PDF for the vertical component. Hence, the combined spatial PDF is given by:

$$f(R,Z)\mathrm{exp}\left(\frac{R^2}{2R_0^2}\frac{|Z|}{Z_0}\right),$$
(2)
We set $`R_0=4.0`$ kpc and $`Z_0=1.5`$ kpc as standard model parameters. Following Narayan et al. (1991), these and other values considered reflect the present-day spatial distribution of the NS–NS binary population after kinematic evolution in the Galactic gravitational potential.

Having assigned a position of each pulsar in our model galaxy, for later computational convenience, we store the positions as Cartesian $`x,y,z`$ coordinates, where the Galactic center is defined as (0.0,0.0,0.0) kpc and the position of the Earth is assumed to lie 8.5 kpc from the center along the $`xy`$ plane, i.e. (8.5,0.0,0.0). From these definitions, the distance $`d`$ to each pulsar from the Earth can be readily calculated, as well as the apparent Galactic coordinates $`l`$ and $`b`$.

For the luminosity PDF, we follow the results of Cordes & Chernoff (1997) and adopt a power-law function of the form

$$\varphi (L)=(p1)L_{\mathrm{min}}^{\mathrm{p}1}L^\mathrm{p},$$
(3)
where $`LL_{\mathrm{min}}`$ and $`p>1`$. The cut-off luminosity, $`L_{\mathrm{min}}`$, and the exponent $`p`$ are the model parameters. Cordes & Chernoff (1997) found $`L_{\mathrm{min}}=1.1_{0.5}^{+0.4}`$ mJy kpc2 and $`p=2.0\pm 0.2`$ at 68% confidence. We set $`p=2.0`$ and $`L_{\mathrm{min}}=1.0`$ mJy kpc2 for our reference model. Throughout this paper, luminosities are defined to be at the observing frequency $`\nu =400`$ MHz.

Having defined the position and luminosity of each pulsar in our model Galaxy, the final step in defining the model population is to calculate a number of derived parameters required to characterize the detection of the model pulsars: dispersion measure (DM), scatter-broadening time ($`\tau `$) and sky background temperature ($`T_{\mathrm{sky}}`$). To calculate DM and $`\tau `$, we use the software developed by Taylor & Cordes (1993) to integrate their model of the free-electron column density along the line of sight to each pulsar defined by its model Galactic coordinates $`l`$ and $`b`$ out to its distance $`d`$. Frequency scaling of $`\tau `$ to different survey frequencies is done assuming a Kolmogorov turbulence spectrum with a spectral index of –4.4111Although recent studies suggest a variety of spectral indices for $`\tau `$ Löhmer et al. (2001), the effects of scattering turn out to be negligible in this study since the detections of NS–NS binaries are limited by luminosity to nearby systems.. Finally, given the model Galactic coordinates of each pulsar, the sky background noise temperature at 408 MHz ($`T_{\mathrm{sky}}`$) is taken from the all-sky catalog of Haslam et al. (1981). Scaling $`T_{\mathrm{sky}}`$ to other survey frequencies assumes a spectral index of –2.8 Lawson et al. (1987).

## 4 PULSAR SURVEY SELECTION EFFECTS

Having created a model pulsar population with a given spatial and luminosity distribution, we are now in a position to determine the fraction of the total population which are actually detectable by current large-scale pulsar surveys. To do this, we need to calculate, for each model pulsar, the effective signal-to-noise ratio it would have in each survey and compare this with the corresponding detection threshold. Only those pulsars which are nominally above the threshold count as being detectable. After performing this process on the entire model pulsar population of size $`N_{\mathrm{tot}}`$, we are left with a sample of $`N_{\mathrm{obs}}`$ pulsars that are nominally detectable by the surveys. Repeating this process many times, we can determine the probability distribution of $`N_{\mathrm{obs}}`$ which we then use to constrain the population and coalescence rate of NS–NS binaries. In this section we discuss our modeling of the various selection effects which limit pulsar detection.

### 4.1 Survey Parameters

The main factors affecting the signal-to-noise ratio ($`\sigma `$) of a pulsar search can be summarized by the following expression

$$\sigma \frac{S_\nu G}{T}\sqrt{\frac{P\mathrm{\Delta }\nu t}{w_e}},$$
(4)
where $`S_\nu `$ is the apparent flux density at the survey frequency $`\nu `$, $`G`$ is the gain of the telescope, $`T`$ is the effective system noise temperature (which includes a contribution $`T_{\mathrm{sky}}`$ from the sky background described in the previous section), $`P`$ is the pulse period, $`\mathrm{\Delta }\nu `$ is the observing bandwidth, $`t`$ is the integration time and $`w_e`$ is the pulse width. More exact expressions are given in the detailed description of the survey selection effects in §2 of LBDH (in particular see their eqns. 14–18) which we adopt in this work. In what follows, we describe the salient points relevant to this study.

For each model pulsar, with known 400-MHz luminosity $`L`$ and distance $`d`$, we calculate the apparent 400-MHz flux density $`S_{400}=L/d^2`$. Since not all pulsar surveys are carried out at 400 MHz, we need to scale $`S_{400}`$ to take account of the steep radio flux density spectra of pulsars. Using a simple power law of the form $`S_\nu \nu ^\alpha `$, where $`\alpha `$ is the spectral index, we can calculate the flux $`S_\nu `$ at any frequency $`\nu `$ as:

$$S_\nu =S_{400}\left(\frac{\nu }{400\mathrm{MHz}}\right)^\alpha ,$$
(5)
Following the results of Lorimer et al. (1995), in all simulations, spectral indices were drawn from a Gaussian PDF with a mean of –1.6 and standard deviation 0.4.

The telescope gain, system noise, bandwidth and integration time are well-known parameters for any given survey and the detailed models we use take account of these. In addition to the surveys considered by LBDH, we also model surveys listed by Curran & Lorimer (1995), and more recent surveys at Green Bank (Sayer, Nice & Taylor 1997) and Parkes Lyne et al. (2000); Manchester et al. (2001); Edwards et al. (2001). A complete list of the surveys considered, and the references to the relevant publications is given in Table 1.

Up to this point in the simulations, the model parameters are identical for both PSRs B1913+16 and B1534+12. Since we are interested in the individual contributions each of these systems make to the total Galactic merger rate of NS–NS binaries similar to these systems, we fix the assumed spin periods $`P`$ and intrinsic pulse widths $`w`$ to the values of each pulsar and perform separate simulations over all physical models considered. The assumed pulse widths are 10 ms and 1.5 ms respectively for PSRs B1913+16 and B1534+12. The effective pulse width $`w_e`$ required for the signal-to-noise calculation must take into account pulse broadening effects due to the interstellar medium and the response of the observing system. The various contributions are summarized by the quadrature sum:

$$w_e^2=w^2+\tau ^2+t_{\mathrm{samp}}^2+t_{\mathrm{DM}}^2+t_{\mathrm{\Delta }\mathrm{DM}}^2,$$
(6)
where $`\tau `$ is the scatter-broadening timescale calculated from the Taylor & Cordes (1993) model, $`t_{\mathrm{samp}}`$ is the data sampling interval in the observing system, $`t_{\mathrm{DM}}`$ is the dispersive broadening across an individual frequency channel and $`t_{\mathrm{\Delta }\mathrm{DM}}`$ is the pulse broadening due to dedispersion at a slightly incorrect dispersion measure. All of these factors are accounted for in our model described in detail in LBDH.

### 4.2 Doppler Smearing

For binary pulsars, we need to take account of the reduction in signal-to-noise ratio due to the Doppler shift in period during an observation. This was not considered in LBDH since their analysis was concerned only with isolated pulsars. For observations of NS–NS binaries, however, where the orbital periods are of the order of 10 hours or less, the apparent pulse period can change significantly during a search observation causing the received power to be spread over a number of frequency bins in the Fourier domain. As all the surveys considered in this analysis search for periodicities in the amplitude spectrum of the Fourier transform of the time series, a signal spread over several bins can result in a loss of signal-to-noise ratio. To take account of this effect in our survey simulations, we need to multiply the apparent flux density of each model pulsar by a “degradation factor”, $`F`$.

To calculate the appropriate $`F`$ values to use, we generate synthetic pulsar search data containing signals with periods and duty cycles similar to PSR B1913+16 and PSR B1534+12. These data are then passed through a real pulsar search code which is similar to those in use in the large-scale surveys (see Lorimer et al. 2000 for details). For each of the two pulsars, we first generate a control time series in which the signal has a constant period and find the resulting signal-to-noise ratio, $`\sigma _{\mathrm{control}}`$, reported by the search code. We then generate a time series in which the pulses have identical intensity but are modulated in period according to the appropriate orbital parameters of each binary pulsar. From the resulting search signal-to-noise ratio, $`\sigma _{\mathrm{binary}}`$, the degradation factor $`F=\sigma _{\mathrm{binary}}/\sigma _{\mathrm{control}}`$. Significant degradation occurs, therefore, when $`F1`$. Since accumulated Doppler shift, and therefore $`F`$, is a strong function of the orbital phase at the start of a given observation, for both binary systems, we calculate the mean value of $`F`$ for a variety of starting orbital phases appropriately weighted by the time spent in that particular part of the orbit.

A similar analysis was made by Camilo et al. (2000) for the millisecond pulsars in 47 Tucanae. In this paper, where we are interested in the degradation as a function of integration time, we generate time series with a variety of lengths between 1 minute and 1 hour using sampling intervals similar to those of the actual surveys listed in Table 1. The results are summarized in Figure The Probability Distribution of Binary Pulsar Coalescence Rates. I. Double Neutron Star Systems in the Galactic Field, where we plot average $`F`$ versus integration time for both sets of orbital parameters. As expected, surveys with the longest integration times are most affected by Doppler smearing. For the Parkes Multibeam survey Lyne et al. (2000); Manchester et al. (2001), which has an integration time of 35 min, mean values of $`F`$ are 0.7 and 0.3 for PSR B1913+16 and PSR B1534+12 respectively222In order to improve on the sensitivity to binary pulsars, the Parkes Multibeam survey data are now being reprocessed using various algorithms designed to account for binary motion during the integration time (Faulkner et al. 2002). The greater degradation for PSR B1534+12 is due to its mildly eccentric orbit ($`e0.3`$ versus 0.6 for PSR B1913+16) which results in a much more persistent change in apparent pulse period when averaged over the entire orbit. For the Jodrell Bank and Swinburne surveys Nicastro et al. (1995); Edwards et al. (2001), which both have integration times of order 5 min, we find $`F0.9`$ for both systems. For all other surveys, which have significantly shorter integration times, no significant degradation is seen, and we take $`F=1`$.

## 5 STATISTICAL ANALYSIS

In this section we describe in detail the derivation of the probability distribution of the Galactic coalescence rate $``$. The analysis method makes use of Bayesian statistics and takes into account the rate contributions of both observed NS–NS binaries. At the end of the section we derive the associated detection rates for LIGO.

### 5.1 The Rate Probability Distribution for Each Observed NS–NS Binary

As already mentioned in § 2, for each of the two observed NS–NS binaries (PSR B1913+16 or PSR B1534+12) we generate pulsar populations in physical and radio luminosity space with pulse periods and widths fixed to the observed ones and with different absolute normalizations, i.e., total number $`N_{\mathrm{tot}}`$ of pulsars in the Galaxy. We generate large numbers of “observed” pulsar samples by modeling the pulsar survey selection effects (see § 4) and applying them on these model populations of PSR B1913+16–like and PSR B1534+12–like pulsars separately (see § 3). For a fixed value of $`N_{\mathrm{tot}}`$, we use these “observed” samples to calculate the distribution of the number of objects in the samples. One might expect that the number of observed pulsars in a sample $`N_{\mathrm{obs}}`$ follows very closely a Poisson distribution:

$$P(N_{\mathrm{obs}};\lambda )=\frac{\lambda ^{N_{\mathrm{obs}}}e^\lambda }{N_{\mathrm{obs}}!},$$
(7)
where, by definition, $`\lambda <N_{\mathrm{obs}}>`$. We confirm our expectation by obtaining excellent formal fits to the Monte Carlo data using such a distribution and calculating the best-fitting value of $`\lambda `$. We vary $`N_{\mathrm{tot}}`$ in the range $`1010^4`$ and find that $`\lambda `$ is linearly correlated with $`N_{\mathrm{tot}}`$:

$$\lambda =\alpha N_{\mathrm{tot}},$$
(8)
where $`\alpha `$ is a constant that depends on the properties (space and luminosity distributions and pulse period and width) of the Galactic pulsar population. Examples of the Poisson fits and the linear correlations are shown in Figures The Probability Distribution of Binary Pulsar Coalescence Rates. I. Double Neutron Star Systems in the Galactic Field and The Probability Distribution of Binary Pulsar Coalescence Rates. I. Double Neutron Star Systems in the Galactic Field, respectively, for our reference model 1 (see Table 2 and § 6).

The main step in deriving a rate probability distribution for each of the observed systems is to first derive the probability distribution of the total number of pulsars like the observed ones in the Galaxy. We obtain the latter by applying Bayes’ theorem:

$$P(H|DX)=P(H|X)\frac{P(D|HX)}{P(D|X)},$$
(9)
where $`P(H|DX)`$ is the probability of a model hypothesis $`H`$ given data $`D`$ and model priors $`X`$, $`P(D|HX)`$ is the likelihood of the data given a model hypothesis and priors, $`P(H|X)`$ is the probability of a model hypothesis in the absence of any data information, and $`P(D|X)`$ is the model prior probability, which acts as a normalization constant.

In the present work, we make following identifications: $`D:\mathrm{is}\mathrm{the}\mathrm{real}\mathrm{observed}\mathrm{sample}H:\mathrm{is}\lambda proportionaltoN_{\mathrm{tot}}X:\mathrm{is}\mathrm{the}\mathrm{population}\mathrm{model}(\mathrm{space}\mathrm{and}\mathrm{luminosity}\mathrm{distributions}\mathrm{and}\mathrm{pulse}\mathrm{period}\mathrm{and}\mathrm{width})`$ With these identifications, $`P(D|HX)`$ is the likelihood of the real observed sample (one “PSR B1913+16–like” and one “PSR B1534+12–like” pulsar) and is obtained by the best-fitting Poisson distribution:

$$P(D|HX)=P(1;\lambda (N_{\mathrm{tot}}),X)=\lambda (N_{\mathrm{tot}})e^{\lambda (N_{\mathrm{tot}})}.$$
(10)
In the absence of any data information, the absolute normalization of the model population, i.e., the total pulsar number $`N_{\mathrm{tot}}`$ and hence $`\lambda `$ is expected to be independent of the shape of the population distributions and properties represented by $`X`$. Therefore, the probability of $`\lambda `$ given a set of assumptions $`X`$ for the model Galactic population is expected to be flat:

$$P(H|X)=P(\lambda (N_{\mathrm{tot}})|X)=\mathrm{constant},$$
(11)
and is essentially absorbed by the model prior probability $`P(D|X)`$ as a normalization constant. The probability distribution of $`\lambda `$ then is given by:

$$P(\lambda |DX)=\frac{P(\lambda |X)}{P(D|X)}P(1;\lambda ,X)=\mathrm{constant}\times P(1;\lambda ,X).$$
(12)
We impose the normalization constraint: $`_0^{\mathrm{}}P(\lambda |DX)𝑑\lambda =1`$ and find that $`P(H|X)/P(D|X)=1`$ and

$$P(\lambda |DX)=P(1;\lambda ,X)=\lambda e^\lambda .$$
(13)
Note that based on the above expression, the maximum value of $`P(\lambda )`$ equal to $`e^1`$ always occurs at $`\lambda =1`$ or at $`N_{\mathrm{tot}}=\alpha ^1`$ (see eq.(8)). It is straightforward to calculate $`P(N_{\mathrm{tot}})`$:

$`P(N_{\mathrm{tot}})`$ $`=`$ $`P(\lambda )\left|{\displaystyle \frac{d\lambda }{dN_{\mathrm{tot}}}}\right|`$ (14)
$`=`$ $`\alpha ^2N_{\mathrm{tot}}e^{\alpha N_{\mathrm{tot}}}`$
For a given total number of pulsars in the Galaxy, we can calculate their rate using estimates of the associated pulsar beaming correction factor $`f_\mathrm{b}`$ and lifetime $`\tau _{\mathrm{life}}`$:

$$=\frac{N_{\mathrm{tot}}}{\tau _{\mathrm{life}}}f_\mathrm{b}.$$
(15)
Note that this estimate is equivalent to previous studies (KNST and references therein) where the concept of the scale factor is used instead of our calculated $`N_{\mathrm{tot}}`$. We can write equivalently for our calculation:

$$\frac{N_{\mathrm{tot}}}{N_{\mathrm{obs}}}=\frac{_{\mathrm{V}_\mathrm{G}}f(R,Z,L)𝑑V𝑑L}{_{\mathrm{V}_\mathrm{D}}f(R,Z,L)𝑑V𝑑L},$$
(16)
where $`f(R,Z,L)`$ is the probability distribution of the pulsar population (§ 3), $`V_\mathrm{G}`$ and $`V_\mathrm{D}`$ are the Galactic volume and the detection volume (for pulsars with pulse period and width similar to each of the two observed pulsars), respectively. Note that we do not fix the luminosity of the pulsar population to the observed values, but instead, we estimate $`N_{\mathrm{tot}}`$ for a distribution of radio luminosities. Since we consider separately PSR B1913+16–like and PSR B1534+12–like populations, $`N_{\mathrm{obs}}`$= 1.

For pulsar beaming fractions, we adopt the estimates obtained by KNST: 5.72 for PSR B1913+16 and 6.45 for PSR B1534+12. For the lifetime estimates we also follow KNST. Our adopted values for the pulsar lifetimes are $`3.65\times 10^8`$ yr for PSR B1913+16 and $`2.9\times 10^9`$ yr for PSR B1534+12.

Using eq. (15) we calculate $`P()`$ for each of the two observed pulsars:

$`P()`$ $`=`$ $`P(N_{\mathrm{tot}})\left|{\displaystyle \frac{dN_{\mathrm{tot}}}{d}}\right|`$ (17)
$`=`$ $`\left({\displaystyle \frac{\alpha \tau _{\mathrm{life}}}{f_\mathrm{b}}}\right)^2e^{\left(\frac{\alpha \tau _{\mathrm{life}}}{f_\mathrm{b}}\right)}.`$
### 5.2 The Total Galactic Coalescence Rate

Once the probability distributions of the rate contributions of the two observed pulsars are calculated, we can obtain the distribution functions of the total coalescence rate $`_{\mathrm{tot}}`$. We define the following two coefficients for each observed system:

$$A\left(\frac{\alpha \tau _{\mathrm{life}}}{f_\mathrm{b}}\right)_{1913}\mathrm{and}B\left(\frac{\alpha \tau _{\mathrm{life}}}{f_\mathrm{b}}\right)_{1534}$$
(18)
and rewrite the coalescence rate for each binary system:

$$P_{1913}()=A^2e^A\mathrm{and}P_{1534}()=B^2e^B.$$
(19)
One can confirm that each distribution function satisfies the normalization condition

$$_0^{\mathrm{}}P()𝑑=1.$$
(20)
We then define two new variables: $`_+_{1913}+_{1534}`$ and $`_{}_{1913}_{1534}`$. Since $`_{1913}0`$ and $`_{1534}0`$, we have $`_+_{}_+`$. For convenience, we rename $`_{1913}=_1`$ and $`_{1534}=_2`$ and perform a two-dimensional probability distribution transformation:

$$P(_+,_{})=P(_1,_2)\left|\begin{array}{cc}\frac{d_1}{d_+}& \frac{d_2}{d_{}}\\ \frac{d_1}{d_{}}& \frac{d_2}{d_+}\end{array}\right|=\frac{1}{2}P(_1,_2).$$
(21)
The probability distribution of the total rate $`_{\mathrm{tot}}_+`$ is obtained after integrating $`P(_+,_{})`$ over $`_{}`$:

$$P(_+)=_{_{}}P(_+,_{})𝑑_{}=\frac{1}{2}_{_{}}P(_1,_2)𝑑_{}.$$
(22)
Since the probability distributions for the rate contributions of each of the two observed pulsars are independent of each other, their two-dimensional distribution is simply the product of the two.

$$P(_1,_2)=P(_1)P(_2).$$
(23)
After we replace all variables to $`_+`$ and $`_{}`$, we have

$`P(_+_{\mathrm{tot}})`$ $`=`$ $`{\displaystyle \frac{A^2B^2}{2^3}}e^{\left(\frac{A+B}{2}\right)_+}\left[_+^2{\displaystyle __+^{+_+}}𝑑_{}e^{\left(\frac{BA}{2}\right)_{}}{\displaystyle __+^{+_+}}𝑑_{}_{}^2e^{\left(\frac{BA}{2}\right)_{}}\right]`$ (24)
$`=`$ $`\left({\displaystyle \frac{AB}{BA}}\right)^2\left[_+\left(e^{A_+}+e^{B_+}\right)\left({\displaystyle \frac{2}{BA}}\right)\left(e^{A_+}e^{B_+}\right)\right].`$
We have confirmed that the above function satisfies the normalization $`_0^{\mathrm{}}P(_{\mathrm{tot}})𝑑_{\mathrm{tot}}=1`$.

Having calculated the probability distribution of the Galactic coalescence rate, we can take one step further and also calculate ranges of values for the rate $`_{\mathrm{tot}}`$ at various confidence levels (CL). The lower ($`_\mathrm{a}`$) and upper ($`_\mathrm{b}`$) limits to these ranges are calculated using:

$$__\mathrm{a}^_\mathrm{b}P(_{\mathrm{tot}})𝑑_{\mathrm{tot}}=\mathrm{CL}$$
(25)
and

$$P(_\mathrm{a})=P(_\mathrm{b}).$$
(26)
In all our results, we quote coalescence rates for 68%, 95% and 99% CL.

### 5.3 The Detection Rate for LIGO

NS–NS binary systems are expected to emit strong gravitational waves during their inspiral phase, the late stages of which may be detected by the ground-based gravitational wave detectors. In this paper we estimate expected detection rates for LIGO (initial and advanced) using the derived probability distribution of the Galactic rate and its values at the maximum of the distribution, for different physical models of pulsar properties. To calculate the detection rate, we need to extrapolate the Galactic rate to the volume detectable by LIGO. We use the ratio between the B-band luminosity density of the Universe and the B-band luminosity of our Galaxy as the scaling factor (Phinney 1991; KNST). This is based on the assumptions that (i) the B-band luminosity (corrected for dust absorption) correlates with the star–formation rate in the nearby universe and hence the coalescence rate, (ii) the B-band luminosity density is constant in the nearby universe. The detection rate, $`_{\mathrm{det}}`$, is calculated by the following equation:

$$_{\mathrm{det}}=ϵ_{\mathrm{tot}}V_{\mathrm{det}},$$
(27)
where $`ϵ`$ is the scaling factor assumed to be $`10^2`$ Mpc-3 (for details see KNST). $`V_{\mathrm{det}}`$ is the detection volume defined as a sphere with a radius equals to the maximum detection distance Dmax for the initial ($`20`$ Mpc) and advanced LIGO ($`350`$ Mpc) Finn (2001).

## 6 RESULTS

We have calculated the probability distribution of the Galactic coalescence rate $`_{\mathrm{tot}}`$, its most likely value $`_{\mathrm{peak}}`$ and ranges at different statistical confidence levels, and the most likely expected detection rates for the initial and advance LIGO for a large number of model pulsar population properties (see Table 2). We have chosen one of them to be our reference model based on the statistical analyses and results presented by Cordes & Chernoff (1997).

For our reference model ($`L_{\mathrm{min}}=1.0`$ mJy kpc2, $`p=2.0`$, $`R_\mathrm{o}=4.0`$ kpc, and $`Z_0=1.5`$ kpc), we find the most likely value of $`N_{\mathrm{tot}}`$, to be $`390`$ pulsars for the “PSR B1913+16-like” population, and $`350`$ pulsars for the “PSR B1534+12-like” population. Using eq. (24), we evaluate the total Galactic coalescence rate of NS–NS binaries for this reference case (model 1 in Table 2). The most likely value of the coalescence rate is $`8`$ Myr-1 and the ranges at different statistical confidence levels are: $`320`$ Myr-1 at 68%, $`130`$ Myr-1 at 95%, and $`0.740`$ Myr-1 at 99%.

In Figure The Probability Distribution of Binary Pulsar Coalescence Rates. I. Double Neutron Star Systems in the Galactic Field, $`P(_{\mathrm{tot}})`$ along with $`P(_{1913})`$ and $`P(_{1534})`$ are shown for the reference model. It is evident that the total rate distribution is dominated by that of PSR B1913+16. At first this appears to be in contradiction to most other studies of the NS–NS coalescence rate (Narayan et al. 1991; Phinney 1991; Curran & Lorimer 1995; van den Heuvel & Lorimer 1996; Stairs et al. 1998; Arzoumanian et al. 1999; KNST). However, it turns out that this difference is due to the fact that earlier studies restricted the scale factor calculation to the actual observed pulsar luminosity and the rate estimates were dominated by the low luminosity of, and hence large scale factor for, PSR B1534+12. Any corrections to the rate estimate which take into account the range of pulsar luminosities were applied as subsequent upward corrections. However, this effect is eliminated in our case because we calculate the two rate contributions having relaxed the luminosity constraint and instead having allowed for a range in luminosity for the pulsars. In this case any differences in the two separate rate contributions depend only on differences in pulse periods, and widths. Given that the latter are rather small, it makes sense that, for example, the most likely values of $`N_{\mathrm{tot}}`$ for the two pulsars come out to be very similar (e.g. $`390`$ and $`350`$, for PSR B1913+16 and PSR B1534+12, respectively, in model 1). Consequently any difference in the rate contributions from the two populations is due to the difference in lifetimes (about a factor of 10) for the two observed pulsars (note that the two do not only have similar $`N_{\mathrm{tot}}`$ estimates, but also similar beaming correction factors). Since the lifetime estimate for PSR B1913+16 is much smaller, the total rate distribution is dominated by its contribution.

In Table 2 we list the population parameters and results of Galactic coalescence and LIGO detection rates for a large number of models (but still a subset of the models we have investigated in detail). Model 1 is our reference model defined above and the rest of the models are used to explore the sensitivity of our results on the pulsar population properties in luminosity and space distributions. The variation of the luminosity function parameters is within ranges that correspond to 68% confidence levels from the statistical analysis of Cordes & Chernoff (1997). The variation of the space scale lengths $`R_0`$ and $`Z_0`$ extends to values that are not favored by our current understanding, only to allow us to examine the presence of any correlations. Scrutiny of the results presented in Table 2 reveals that rate estimates are modestly sensitive (and in some cases, e.g., $`Z_0`$, insensitive) to most of these variations, except in cases of very low values of $`L_{\mathrm{min}}`$, down to 0.3 mJy kpc2, and unphysically low values of $`R_0`$, down to $`23`$ kpc. Assumptions about the shape of the space distribution (exponential or Gaussian) are not important. The most important model parameter seems to be the slope of the luminosity function. The most likely values for the rates are found in the range $`322`$ Myr-1 for all models with luminosity-function parameters consistent with a 68% confidence level estimated by Cordes & Chernoff (1997) and $`R_0`$ in the range $`48`$ kpc. Within a given model the range of estimated values at 68% confidence level is typically broad by a factor of 5, while at the 95% level the ranges broaden by another factor of about 5, reaching an uncertainty of $`25`$ for the rate estimates.

In terms of qualitative variations, it is clear that models with increasing fraction of low-luminosity or very distant pulsars lead to increasing rate estimates, as expected. A demonstration of this kind of dependence is shown in Figures The Probability Distribution of Binary Pulsar Coalescence Rates. I. Double Neutron Star Systems in the Galactic Field and The Probability Distribution of Binary Pulsar Coalescence Rates. I. Double Neutron Star Systems in the Galactic Field. We have found that there is a strong correlation between the peak value of the total Galactic coalescence rate $`_{\mathrm{peak}}`$ estimated by eq. (24) and the cut-off luminosity $`L_{\mathrm{min}}`$ and the power index p. As seen in these Figures $`_{\mathrm{peak}}`$ increases rapidly with decreasing $`L_{\mathrm{min}}`$ (Fig. The Probability Distribution of Binary Pulsar Coalescence Rates. I. Double Neutron Star Systems in the Galactic Field) or with increasing p (Fig. The Probability Distribution of Binary Pulsar Coalescence Rates. I. Double Neutron Star Systems in the Galactic Field). Scale lengths of the spatial distribution (either $`R_0`$ or $`Z_0`$) do not show a correlation with $`_{\mathrm{peak}}`$ as strong as that of the luminosity-function parameters. In Figure The Probability Distribution of Binary Pulsar Coalescence Rates. I. Double Neutron Star Systems in the Galactic Field the most likely rate estimate is plotted as a function of the radial scale length $`R_0`$. It is evident that the rate is relatively insensitive to $`R_0`$ variations unless $`R_0`$ becomes very small ($`3`$ kpc). Small values of the radial scale for the Galactic distribution imply a large fraction of very distant (and hence undetectable) pulsars in the Galaxy and lead to a rapid increase of the estimated coalescence rate. However, such small $`R_0`$ values are not consistent with our current understanding of stellar populations in the Galaxy Lyen et al. (1985).

## 7 DISCUSSION

In this paper we present a new method of estimating the total number of pulsars in our Galaxy and we apply it to the calculation of the coalescence rate of double neutron star systems in the Galactic field. The method implicitly takes into account the small number of pulsars in the observed double-neutron-star sample as well as their distribution in luminosity and space in the Galaxy. The modeling of pulsar survey selection effects is formulated in a “forward” way, by populating the Galaxy with model pulsar populations and calculating the likelihood of the real observed sample. This is in contrast to the “inverse” way of the calculation of scale factors used in previous studies. The formulation presented here allows us to: (a) calculate the probability distribution of coalescence rates; (b) assign statistical significance to these estimates; and (c) quantify the uncertainties associated with them.

As originally shown by KNST, the most important uncertainties originate from the combination of a small-number observed pulsar and a pulsar population dominated by faint objects. The probability distribution covers more than 2 orders of magnitude in agreement with the uncertainties in excess of two orders of magnitudes asserted by KNST. However, for the reference model, even at high statistical confidence level (99%) the uncertainty is reduced to a factor of $`60`$. At confidence levels of 95% and 68%, the uncertainty is further reduced to just $`25`$ and $`5`$, respectively. We use our results to estimate the expected detection rates for ground-based interferometers, such as LIGO. The most likely values are found in the range $`(130)\times 10^3`$ yr-1 and $`4140`$ yr-1, for the initial and advanced LIGO.

The statistical method developed here can be further extended to account for distributions of pulsar populations in pulse periods, widths, and orbital periods. Most importantly the method can be applied to any type of pulsar population with appropriate modifications of the modeling of survey selection effects. Currently we are working on assessing the contribution of double neutron stars formed in globular clusters as well as the formation rate of binary pulsars with white dwarf companions that are important for gravitational-wave detection by LISA, the space-based interferometer planned by NASA and ESA for the end of this decade.

This work is partially supported by NSF grant PHY-0121420 to VK. DRL is a University Research Fellow funded by the Royal Society. DRL is also grateful for the hospitality and support of the Theoretical Astrophysics Group at Northwestern U. VK would also like to acknowledge partial support by the Clay Fellowship at the Smithsonian Astrophysical Observatory where this work was initiated.


text:
# On non-abelian homomorphic public-key cryptosystems

## 1 Introduction

In what follows all the groups are presented in some natural way depending on the problem. For example, the special constructions of Section 2 are based on the groups $`_n^+`$ and $`_n^{}`$ just given via $`n`$, whereas the general construction of Section 3 requires only that elements of a group in question can be generated and moreover, the multiplication and taking the inverse in the group can be performed efficiently. In the latter case the groups can be presented by generators and relations or even by generic algorithms (see e.g. ).

There is a lot of public-key cryptosystems using groups (see e.g. ) but only a few of them have a homomorphic property in the sense of the following definition (cf. also ).

###### Definition 1.1

Let $`H`$ be a finite non-identity group, $`G`$ a finitely generated group and $`f:GH`$ an epimorphism. Suppose that $`R`$ is a set of distinct representatives of the right cosets of $`G`$ with respect to $`ker(f)`$, $`A`$ is a set of words in some alphabet and a mapping $`P:AG`$ such that $`im(P)=ker(f)`$. A triple $`𝒮=(R,A,P)`$ is called a homomorphic cryptosystem over $`H`$ with respect to $`f`$, if the following conditions are satisfied:

1. one can get random elements (of the sets $`A,G,H`$), compute the inverse of an element and the product of two elements (in the group $`G`$ or $`H`$) in polynomial in $`N`$ probabilistic time where $`N`$ is the size of presentations of $`G,H`$ and $`A`$;
2. $`|R|=|H|`$ and for any element $`gR`$ its image $`f(g)`$ as well as for any element $`hH`$ its unique preimage $`gR`$ such that $`f(g)=h`$ can be computed in polynomial in $`N`$ probabilistic time;
3. the mapping $`P`$ is a trapdoor function.

###### Remark 1.2

We require that the set $`R`$ is given explicitly by a list of elements of $`G`$. So, condition (H2) implies that without loss of generality one can assume that the group $`H`$ is represented by its multiplication table.

Condition (H3) means (see ) that the values of $`P`$ can be computed in polynomial in $`N`$ probabilistic time, whereas finding of the inverse mapping $`P^1`$ is a hard computational problem which can be solved with the help of some additional secret information (for instance, knowing some invariant of the group $`G`$). In a homomorphic cryptosystem $`𝒮`$ the elements of $`H`$ are (publically) encrypted in a probabilistic manner by the elements of $`G`$, all the computations are performed in $`G`$ and the result is decrypted to $`H`$. More precisely:

Public Key: $`G,H,R,A,P,f|_R`$.

Secret Key: finding $`P^1`$.

Encryption: given a plaintext $`hH`$ take $`rR`$ such that $`f(r)=h`$ (invoking (H2)) and a random element $`aA`$; the ciphertext of $`h`$ is the element $`P(a)r`$ of $`G`$ (the element $`a`$ as well as the product $`P(a)r`$ is computed by means of (H1)).

Decryption: given $`gG`$ find the elements $`rR`$ and $`aA`$ such that $`rg^1=P(a)`$ (for computing $`P(a)`$ see (H3)); set the plaintext of $`g`$ to be $`f(g)=f(r)`$ (the element $`f(r)`$ is computed by means of (H2)).

One can see that the encryption procedure can be performed by means of public keys efficiently. However, the decryption procedure is a secret one in the following sense. To find the element $`r`$ one has to solve in fact, the membership problem for the subgroup $`ker(f)`$ of the group $`G`$. We assume that a solution for each instance $`g^{}ker(f)`$ of this problem must have a “proof”, which is actually an element $`aP^1(g^{})`$. Thus, the secrecy of the system is based on the assumption that finding an element in the set $`P^1(g^{})`$ is an intractable computation problem. On the other hand, our ability to compute $`P^1`$ enables us to efficiently implement the decryption algorithm. One can treat $`P`$ as a proof system for $`ker(f)`$ in the sense of . Moreover, in case when $`A`$ is a certain group and $`P`$ is a homomorphism we have the following exact sequence of group homomorphisms

$$A\stackrel{P}{}G\stackrel{f}{}H\{1\}$$

(recall that the exact sequence means that the image of each homomorphism in it coincides with the kernel of the next one).

In the present paper the group $`H`$ being an alphabet of plaintext messages is always finite (and rather small) and given by its multiplication table, while the group $`G`$ of ciphertext messages could be infinite but being always finitely generated. However, the infinitness of $`G`$ is not an obstacle for encrypting (and decrypting) since an element from $`H`$ is encrypted by a finite word in generators of $`G`$. For example, in  for an (infinite, non-abelian in general) group $`H`$ given by $`m`$ generators and relations a natural epimorphism $`f:F_mH`$ from a free group $`F_m`$ is considered. Thus, for any element of $`H`$ one can produce its preimages (encryptions) by inserting in a word (being already a produced preimage of $`f`$) from $`F_m`$ any relation defining $`H`$. In other terms, decrypting of $`f`$ reduces to the word problem in $`H`$. The main difference with our approach is that we consider free products over groups of number-theoretic nature like $`Z_n^{}`$ (rather than given by generators and relations). This allows one to provide evidence for difficulty of decryption.

###### Definition 1.3

$`𝒢_{crypt}`$ is the class of all finite groups $`H`$ for which there exists a homomorphic cryptosystem over $`H`$.

In the context of our definition of a homomorphic cryptosystem the main problem we study in this paper is to prove that the class $`𝒢_{crypt}`$ contains all solvable nonidentity groups (see Theorem 3.6).

To our knowledge all known at present homomorphic cryptosystems are more or less modifications of the following one. Let $`n`$ be the product of two distinct large primes of size $`O(\mathrm{log}n)`$. Set $`G=\{g_n^{}:𝐉_n(g)=1\}`$ where $`𝐉_n`$ is the Jacobi symbol, and $`H=_2^+`$. Then given a non-square $`g_0G`$ the triple $`(R,A,P)`$ where

$$R=\{1,g_0\},A=_n^{},P(g):gg^2,$$

is a homomorphic cryptosystem over $`H`$ with respect to the natural epimorphism $`f:GH`$ with $`ker(f)=\{g^2:g_n^{}\}`$ (see ). We call it the quadratic residue cryptosystem. It can be proved (see ) that in this case finding $`P^1`$ is not easier than factoring $`n`$, whereas given a prime divisor of $`n`$ the computation of $`P^1`$ can be performed in polynomial time in $`\mathrm{log}n`$.

It is an essential assumption (being a shortcoming) in the quadratic residue cryptosystem as well as other cryptosystems cited below that its security relies on a fixed a priori (proof system) $`P`$. Indeed, it is not excluded that adversary could verify whether an element of $`G`$ belongs to $`ker(f)`$ avoiding making use of $`P`$, for example, in case of the quadratic residue cryptosystem that would mean verifying that $`gG`$ is a square without providing a square root of $`g`$. Although, there is a common conjecture that verifying for an element to be a square (as well as some power) is also difficult.

Let us mention that a cryptosystem from  over $`H=_n^+`$ (for the same assumptions on $`n`$ as in the quadratic residue cryptosystem) with respect to the homomorphism $`f:GH`$ where $`G=_{n^2}^{}`$ and $`ker(f)=\{g^n:gG\}`$, in which $`A=G`$ and $`P:gg^n`$, is not homomorphic in the sense of Definition 1.1 because condition (H3) of it does not hold. (Since $`|G||H|^2`$, one can inverse $`P`$ in a polynomial time in $`|H|`$.) By the same reason the cryptosystem from  over $`H=_p^+`$ with respect to the homomorphism $`f:GH`$ where $`G=_{p^2q}^{}`$ and $`ker(f)=\{g^{pq}:gG\}`$ (here the integers $`p,q`$ are distinct large primes of the same size) is also not homomorphic (besides, in this system only a part of the group $`H`$ is encrypted). Some cryptosystems over certain dihedral groups were studied in .

We note in addition that an alternative setting of a homomorphic (in fact, isomorphic) encryption $`E`$ (and a decryption $`D=E^1`$) was proposed in . Unlike Definition 1.1 the encryption $`E:GG`$ is executed in the same set $`G`$ (being an elliptic curve over the ring $`_n`$) treated as the set of plaintext messages. If $`n`$ is composite, then $`G`$ is not a group while being endowed with a partially defined binary operation which converts $`G`$ in a group when $`n`$ is prime. The problem of decrypting this cryptosystem is close to the factoring of $`n`$. In this aspect  is similar to the well-known RSA scheme (see e.g. ) if to interprete RSA as a homomorphism (in fact, isomorphism) $`E:Z_n^{}Z_n^{}`$, for which the security relies on the difficulty of finding the order of the group $`Z_n^{}`$.

We complete the introduction by mentioning some cryptosystems using groups but not being homomorphic in the sense of Definition 1.1. The well-known example is a cryptosystem which relies on the Diffie-Hellman key agreement protocol (see e.g. ). It involves cyclic groups and relates to the discrete logarithm problem ; the complexity of this system was studied in . Some generalizations of this system to non-abelian groups (in particular, the matrix groups over some rings) were suggested in  where secrecy was based on an analog of the discrete logarithm problems in groups of inner automorphisms. Certain variations of the Diffie-Hellman systems over the braid groups were described in ; here several trapdoor one-way functions connected with the conjugacy and the taking root problems in the braid groups were proposed. Finally it should be noted that a cryptosystem from  is based on a monomorphism $`_m^+_n^{}`$ by means of which $`x`$ is encrypted by $`g^x(modn)`$ where $`n,g`$ constitute a public key; its decrypting relates to the discrete logarithm problem and is feasible in this situation due to a special choice of $`n`$ and $`m`$ (cf. also ).

## 2 Homomorphic cryptosystems over cyclic groups

In this section we present an explicit homomorphic cryptosystem over a cyclic group of a prime order $`m`$ whose decription is based on taking $`m`$-roots in the group $`_n^{}`$ for a suitable $`n`$. It can be considered in a sense as a generalization of the quadratic residue cryptosystem over $`_2^+`$. Throughout this section given $`n`$ we denote by $`|n|`$ the size of the number $`n`$.

Given $`m,N`$ set $`𝒯_N=\{(p,q):p,q\mathrm{are}\mathrm{primes},|p|=|q|=N,p<q\}`$ and

$$D_{N,m}=\{n:n=pq,(p,q)𝒯_N,m|p1,GCD(m,q1)=1\}.$$

From the Dirichlet’s theorem on primes in arithmetic progressions  it follows that given an odd prime $`m`$, the set $`D_{N,m}`$ is not empty for sufficiently large numbers $`N`$.

Let $`nD_{N,m}`$ for some natural number $`N`$ and an odd prime $`m`$. Then the group $`G=_n^{}`$ has a (normal) subgroup $`G_0=\{g^m:gG\}`$ the factor by which is isomorphic to the group $`H=_m^+`$. Denote by $`f`$ the corresponding epimorphism from $`G`$ to $`H`$. The mapping

$$P:GG,gg^m$$
(1)
is obviously a polynomial time computable homomorphism such that $`im(P)=ker(f)`$. Next, any element of the set

$$R_{m,n}=\{RG:|f(R)|=|R|=m\}$$

is a system of distinct representatives of the cosets of $`G`$ by $`G_0`$. We observe that given the decomposition $`n=pq`$ one can find an element $`RR_{m,n}`$ in probabilistic time $`|n|^{O(1)}`$. Indeed, since $`m`$ is a prime, it suffices to compute a random element $`s_p_p^{}`$ such that $`s_p^{(p1)/m}1`$, and an element $`s_q_q^{}`$, then find by the Chineese reminder theorem the unique element $`s_q^{}`$ such that $`s=s_p(modp)`$, $`s=s_q(modq)`$, and set $`R=\{s^it_i^m:i=0,\mathrm{},m1\}`$ for arbitrary elements $`t_i_n^{}`$.

We claim that the triple $`𝒮_{N,m,n}=(R,A,P)`$ with arbitrary chosen set $`RR_{m,n}`$, $`A=G`$ and $`P`$ defined by (1) is a homomorphic cryptosystem over the group $`H`$ with respect to the epimorphism $`f`$ whenever the following statement is true:

Assumption (\*). For an odd prime $`m`$ the problem $`𝒫(m)`$, of finding the $`m`$-root in $`_n^{}`$ with $`nD_{N,m}`$ given an element $`RR_{m,n}`$ is not easier than the same problem without any such $`R`$.

Let us present the group $`G`$ by the number $`n`$ and the group $`H`$ by the set of its elements. Then for the triple $`S_{N,m,n}`$ conditions (H1) and (H2) of Definition 1.1 are trivially satisfied (the image of the above element $`s^it_i^m`$ with respect to the homomorphism $`f`$ equals $`i_m^+`$). In fact, condition (H3) would follow from the next lemma.

###### Lemma 2.1

Let $`N`$, $`m`$ be an odd prime and $`nD_{N,m}`$. Then

1. given primes $`p`$ and $`q`$ such that $`n=pq`$ and an element $`g_n^{}`$ one can verify whether $`g`$ is an $`m`$-power and if it is the case one can find an $`m`$-root of $`g`$ in probabilistic polynomial time in $`N`$;
2. the factoring problem for $`n`$ is probabilistic polynomial time reducible to the problem of finding an $`m`$-root in $`_n^{}`$.

Proof. Throughout the proof we will use the canonical decomposition $`_n^{}=_p^{}\times _q^{}`$. To prove statement (1) we make use of Rabin’s probabilistic polynomial-time algorithm for finding roots of polynomials over finite prime fields (see ). Namely, given the primes $`p,q`$ and $`g_n^{}`$ we proceed as follows:

* Step 1. Find the elements $`g_p_p^{}`$ and $`g_q_q^{}`$ such that $`g=g_p\times g_q`$, i.e. $`g_p=g(modp)`$, $`g_q=g(modq)`$.
* Step 2. By Rabin’s algorithm (for a prime field) find some roots $`h_p_p^{}`$ and $`h_q_q^{}`$ of the polynomials $`x^mg_p`$ and $`x^mg_q`$, respectively.
* Step 3. Output $`h=h_p\times h_q`$.

Observe that the described algorithm fails (at Step 2) if and only if $`g`$ is not an $`m`$-power. Since, obviously, $`h^m=h_p^m\times h_q^m=g_p\times g_q=g`$, statement (1) of the lemma is proved.

To prove statement (2) suppose that we are supplied with a probabilistic polynomial-time algorithm $`Q_n`$ that given $`g_n^{}`$ computes an $`m`$-root $`Q_n(g)`$ of $`g`$. The following procedure using well-known observations  shows how $`Q_n`$ helps to find the numbers $`p`$ and $`q`$.

* Step 1. Randomly choose $`x_n^{}`$.
* Step 2. Set $`y=Q_n(x^m)`$. If $`x=y`$, then go to Step 1.
* Step 3. Output $`q=GCD(xy,n)`$ and $`p=n/q`$.

Let $`x=x_p\times x_q`$ and $`y=y_p\times y_q`$ where $`x_p,y_p_p^{}`$ and $`x_q,y_q_q^{}`$. From Step 2 it follows that $`x_q^m=y_q^m`$. On the other hand, since $`nD_{N,m}`$, we have $`GCD(q1,m)=1`$. Thus $`x_q=y_q(modq)`$ and hence

$$x=x_q=y_q=y(modq).$$

So, $`xy0(modn)`$ is a multiple of $`q`$. To complete the proof we note that since $`m=O(1)`$, the loop of Steps 1,2 terminates with a large probability after a polynomial number of iterations.

Unfortunately, we don’t know how to apply this lemma without assumption (\*) because in our case the system $`𝒮_{N,m,n}`$ includes the set $`RR_{m,n}`$. However, from it we obtain the following statement.

###### Theorem 2.2

Under assumption (\*) the triple $`𝒮_{N,m,n}`$ for an odd prime $`m`$ is a homomorphic cryptosystem over $`_m^+`$; in particular, the class $`𝒢_{crypt}`$ contains each cyclic group of a prime order.

We complete the section by mentioning that $`𝒮_{N,m,n}`$ can be slightly modified to avoid the applying of Rabin’s algorithm for finding roots of polynomials over finite fields. In principle, to implement the decryption algorithm it suffices to determine whether a given number $`gG`$ belongs to the group $`G_0`$ or not. However, this can be done by observing that $`gG_0`$ iff $`g_p^{(p1)/m}=1(modp)`$ where $`g_p`$ is the component of $`g`$ in the factor $`_p^{}`$ of $`G=_p^{}\times _q^{}`$.

## 3 Homomorphic cryptosystems using free products

3.1. Throughout the section for a set $`X`$ we denote by $`W(X)`$ the set of all words in the alphabet $`X`$. For an element $`wW(X)`$ we denote by $`|w|`$ the length of $`w`$.

Let $`G_1,\mathrm{},G_m`$ be a set of $`m1`$ pairwise disjoint finite groups. For $`i=1,\mathrm{},m`$ set

$$X_i=G_i\{1_{G_i}\},_i=\{xyzW(X_i):x,y,zX_i,z^1=xy\}.$$

Then $`G_i=X_i;_i`$, i.e. $`G_i`$ is the group given by the set $`X_i`$ of generators and the set $`_i`$ of relations. Set $`X_G=_{i=1}^mX_i`$ and $`_G=_{i=1}^m_i`$. The group

$$G=G_1\mathrm{}G_m=X_G;_G$$

is called the free product of the groups $`G_1,\mathrm{},G_m`$ (see ). From the definition it follows that each element of $`G`$ can be represented by the uniquely determined (canonical) word of $`W_G=W(X_G)`$ such that no two adjacent letters of it belong to the same set among the sets $`X_i`$. This enables us to identify $`G`$ with the subset of $`W_G`$ consisting of all such words. Thus $`G=\{\overline{w}W_G:wW_G\}`$ where $`\overline{w}`$ is the canonical word corresponding to a word $`w`$. In particular, $`\overline{g}=g`$ for all $`gG`$. Due to identifying the groups $`G_1,\mathrm{},G_m`$ and $`G`$ with the corresponding subsets of the set $`W_G`$, we will assume below that the identities of these groups are equal to the empty word of $`W_G`$.

Suppose we are given epimorphisms $`f_i:G_iK_i`$, $`i=1,\mathrm{},m`$. Assuming the groups $`K_1,\mathrm{},K_m`$ being pairwise disjoint we set $`K=K_1\mathrm{}K_m`$ and $`W_K=W(X_K)`$ where $`X_K=_{i=1}^m(K_i\{1_{K_i}\})`$. Then the natural surjection $`W_GW_K`$ replacing the elements of $`X_G`$ by their images in $`X_K`$ with respect to corresponding $`f_i`$, induces an epimorphism

$$f^{}:GK,f^{}|_{X_i}=f_i,i=1,\mathrm{},m.$$
(2)
Moreover, since the conditions $`f^{}|_{X_i}=f_i`$ define the images of the generators of $`G`$, the epimorphism $`f^{}`$ is the unique epimorphism from $`G`$ onto $`K`$ satisfying these conditions.

3.2. Let us study the kernel of the epimorphism $`f^{}`$. To do this suppose that $`K_i`$ is a cyclic group of a prime order and $`(R_i,A_i,P_i)`$ is a homomorphic cryptosystem over the group $`K_i`$ with respect to the epimorphism $`f_i:G_iK_i`$, described in Section 2, $`i=1,\mathrm{},m`$. Let $`w_1,w_2W_G`$ and $`xG_i`$ for some $`i`$. By an elementary transformation of the word $`w_1xw_2W_G`$ we mean replacing of $`x`$ by a certain word $`x_1P_j(a_j)x_2`$ where $`x_1,x_2G_i`$ with $`x=\overline{x_1x_2}`$ and $`a_jA_j`$ for some $`j`$:

$$w_1xw_2w_1(x_1P_j(a_j)x_2)w_2.$$
(3)
Denote by $`W_0`$ the set of all words of $`W_G`$ which can be obtained from the empty word by a sequence of elementary transformations and set $`G_0=\{\overline{w}:wW_0\}`$.

###### Lemma 3.1

In the above notations, $`ker(f^{})=G_0`$.

Proof. First, let $`w=w_1xw_2`$ where $`w_1,w_2W_G`$ and $`xG_i`$ for some $`i`$, to be a word of $`W_G`$. Suppose that a word $`w^{}W_G`$ is obtained from $`w`$ by the elementary transformation (3). Then the words $`\overline{w},\overline{w^{}}`$ belong or not to the group $`ker(f^{})`$ simultaneously. Indeed,

$$f^{}(\overline{w}^{})=f^{}(\overline{w}_1\overline{x}_1\overline{P_j(a_j)}\overline{x}_2\overline{w}_2)=f^{}(\overline{w}_1)f^{}(\overline{x}_1)f^{}(\overline{x}_2)f^{}(\overline{w}_2)=$$

$$f^{}(\overline{w}_1)f^{}(\overline{x})f^{}(\overline{w}_2)=f^{}(\overline{w_1xw_2})=f^{}(\overline{w}).$$
(4)
So, the inclusion $`ker(f^{})G_0`$ follows by the induction on the number of elementary transformation used for constructing an element of $`G_0`$.

Conversely, let $`w^{}ker(f^{})`$. Let us prove that $`w^{}G_0`$ by the induction on $`|w^{}|`$. If $`|w^{}|=0`$, then the statement is obvious. Suppose $`|w^{}|>0`$. Then $`w^{}=w_1xw_2`$ for some $`w_1,w_2W_G`$ and $`xker(f_i)`$ for some $`i`$. (Indeed, otherwise since $`\overline{w^{}}=w^{}`$, we conclude that $`|w^{}|=|f^{}(w^{})|`$. So, $`|f^{}(w^{})|>0`$ which contradicts the fact that $`w^{}ker(f^{})`$.) So, $`w^{}`$ is obtained from the word $`w=w_1w_2`$ of $`W_G`$ by the elementary transformation (3). Since $`w^{}ker(f^{})`$ from (4) it follows that $`\overline{w}ker(f^{})`$. On the other hand, it is easy to see that $`|\overline{w}||w|<|w^{}|`$. So, by the induction hypothesis we conclude that $`\overline{w}G_0`$. By the definition of $`\overline{w}`$ this implies that $`w^{}G_0`$. Thus $`ker(f^{})G_0`$ and we are done.

Let $`gker(f^{})`$. Then from Lemma 3.1 it follows that $`g`$ can be obtained from the empty word by a sequence of elementary transformations. Moreover, the proof of this lemma implies that there exists such a sequence consisting of at most $`|g|`$ elementary transformations. Any such sequence is called a proof for $`g`$ (more precisely, a proof of the membership of $`gker(f^{})`$, cf. (H3) in Definition 1.1). It is easy to see that any elementary transformation (3) is uniquely determined by the following data: the position of the letter $`xG_i`$, the word $`x_1x_2x^1_j`$ and the element $`a_jA_j`$. Thus any proof for the element $`g`$ can be represented by a word $`p`$ in the alphabet $`\times _G\times (_iA_i)`$. One can see that in this case $`|p|`$ is bounded by a polynomial in $`|g|`$.

We define $`A^{}`$ to be the set of all proofs for the elements of $`ker(f^{})`$. It should be stressed that $`A^{}`$ includes only “short” (consisting of at most $`|g|`$ elementary transformations) proofs for an element $`gker(f^{})`$ and does not contain “proofs” for all words of $`W_0`$. For a given $`s`$ one can generate a random element of $`A^{}`$ of the length $`s`$ in time $`s^{O(1)}`$. Indeed, due to the definition of the elementary transformation (3) it suffices to choose randomly positions in a current word and elements of $`A_i`$ for all $`i=1,\mathrm{},m`$. However, this can be done with the help of the algorithms of the homomorphic cryptosystem $`(R_i,A_i,P_i)`$ over $`K_i`$ (see condition (H1) of Definition 1.1).

###### Lemma 3.2

The image of the mapping $`P^{}:A^{}G`$ defined by $`P^{}(a)=g`$ iff $`a`$ is a proof for $`g`$, equals $`ker(f^{})`$. Moreover, the following statements hold:

1. given $`aA^{}`$ the element $`P^{}(a)`$ can be found in polynomial time in $`|a|`$,
2. if for each $`i\{1,\mathrm{},m\}`$ there is an oracle $`Q_i`$ which for any element $`g_iker(f_i)`$ produces a certain $`a_iP_i^1(g_i)`$, then given $`gker(f^{})`$ a proof $`aA^{}`$ for $`g`$ can be found by means of at most $`|g|^2`$ calls of oracles $`Q_i`$ for $`P_i^1(g_i)`$, $`i=1,\mathrm{},m`$, $`g_iker(f_i)`$,
3. for each $`i\{1,\mathrm{},m\}`$ and $`gker(f_i)`$ the problem of finding an element in $`P_i^1(g)`$ is polynomial time reducible to the problem of finding an element in $`(P^{})^1(g)`$.

Proof. The equality $`im(P^{})=ker(f^{})`$ follows from Lemma 3.1. Statement (i1) follows from the fact that any elementary transformation (3) of $`gW_G`$ is reduced to finding $`P_j(a_j)`$ which can be done in polynomial in $`|a_j||a|`$ time. To prove statement (i2) one can apply the following obvious procedure testing membership of a word $`wW_G`$ to the set $`W_0`$.

* Step 1. Using multiplications in the groups $`G_i`$, $`i=1,\mathrm{},m`$, find the canonical word $`\overline{w}`$ of the word $`w`$.
* Step 2. Using the oracles $`Q_i`$, $`i=1,\mathrm{},m`$, delete any letter $`xker(f_i)`$ from the (current) word $`\overline{w}`$. If there was at least one deletion, then go to Step 1.
* Step 3. If the resulting word is empty, then $`wW_0`$.

In fact, this procedure is the algorithmic version of the proof of Lemma 3.1. To find a proof for arbitrary $`gker(f^{})`$, it suffices to apply the above procedure to the word $`gW_G`$ and to collect all results at Steps 1 and 2. Since the number of them is at most $`|g|`$, and the number of calls the oracles $`Q_i`$ at Step 2 is also at most $`|g|`$, statement (i2) follows.

To prove statement (i3) let $`i\{1,\mathrm{},m\}`$ and $`gG_i`$. Then since obviously $`gker(f_i)`$ iff $`gker(f^{})`$, one can test whether $`gker(f_i)`$ by means of an algorithm finding $`(P^{})^1`$. Moreover, if $`gker(f_i)`$, then this algorithm yields a proof from $`A^{}`$ for the element $`g`$. Set $`T`$ to be the set of all elements $`a_jA_i`$ of elementary transformations (3) belonging to this proof. Then $`g=_{a_jT}P(a_j)`$. Since the set $`A_i`$ is an Abelian group and the mapping $`P_i:A_iG_i`$ is a homomorphism, this implies that $`a=_{a_jT}a_j`$ is a proof for $`g`$ and we are done.

Let us describe a special system of distinct representatives of the right cosets of $`G`$ with respect to $`ker(f^{})`$. Set $`W_R=W(_iR_i)`$. Then $`W_RW_G`$ and the set

$$R^{}=GW_R$$
(5)
is a system of distinct representatives of the right cosets of $`G`$ with respect to $`ker(f^{})`$. Indeed, $`R^{}`$ consists of all the words of $`W_R`$ which are canonical words of $`W_G`$. So, no two elements of $`R^{}`$ belong to the same right coset of $`G`$ with respect to $`ker(f^{})`$. Thus our claim follows from the fact that the restriction of the mapping $`f^{}`$ on $`R^{}`$ is the bijection from $`R^{}`$ onto $`K`$ coinciding with $`f_i`$ on $`R_i`$, $`i=1,\mathrm{},m`$.

3.3. We need one more special homomorphism of a free product. To do this we recall some facts on semidirect products of groups (see e.g. ). Suppose that $`K_1,K_2`$ are groups and $`\phi :K_2Aut(K_1)`$ is a homomorphism. Then the set $`K_1\times K_2`$ forms a group with the multiplication given by

$$(k_1,k_2)(l_1,l_2)=(k_1(l_1)^{\phi (k_2^1)},k_2l_2),k_1,l_1K_1,k_2,l_2K_2.$$

This group is called a semidirect product of $`K_1`$ and $`K_2`$ (with respect to the homomorphism $`\phi `$) and is denoted by $`H=\mathrm{\Pi }(K_1,K_2)`$. One can see that it contains the subgroup $`K_i^{}`$ (isomorphic to $`K_i`$) consisting of all pairs with $`1_{K_{3i}}`$ as the $`(3i)`$-th coordinate, $`i=1,2`$. Moreover, $`K_1^{}`$ is a normal subgroup of $`H`$, $`K_1^{}K_2^{}=\{1_H\}`$ and $`H=K_1^{}K_2^{}`$. In general, if an arbitrary group $`H`$ have such two subgroups $`K_1^{},K_2^{}`$, then it is isomorphic to the semidirect product of them with respect to the homomorphism $`\phi `$ induced by the action of $`K_2^{}`$ on $`K_1^{}`$ by the conjugation. In what follows we shall identify $`K_i`$ with $`K_i^{}`$, $`i=1,2`$. We also extend the definition of the semidirect product to arbitrary number of factors by means of setting for $`m3`$

$$\mathrm{\Pi }(K_1,K_2,\mathrm{},K_m)=\mathrm{\Pi }(K_1,\mathrm{\Pi }(K_2,\mathrm{},K_m))$$

with respect to the suitable homomorphisms $`\phi `$. Thus $`\mathrm{\Pi }(K_i,\mathrm{},K_m)=\mathrm{\Pi }(K_i,K^{(i)})`$ where $`K^{(i)}=\mathrm{\Pi }(K_{i+1},\mathrm{},K_m)`$, for all $`i=1,\mathrm{},m1`$ (for $`i=m1`$ we adopt that $`\mathrm{\Pi }(K_m)=K_m`$). In what follows the group $`\mathrm{\Pi }(K_1,K_2,\mathrm{},K_m)`$ will be “small” and presented by its multiplication table. Thus, its subgroups $`K_1,\mathrm{},K_m`$ are also small and for a given $`i`$ the homomorphism $`\phi _i:K^{(i)}Aut(K_i)`$ can be presented by indicating the permutations $`\phi _i(k)`$ of the set $`K_i`$ for all $`kK^{(i)}`$.

###### Lemma 3.3

Let $`H=\mathrm{\Pi }(K_1,\mathrm{},K_m)`$ and $`K=K_1\mathrm{}K_m`$ for a set of pairwise disjoint finite groups $`K_1,\mathrm{},K_m`$. Then there exists an epimorphism $`Q:KH`$ such that given $`kK`$ one can find the element $`Q(k)`$ in time polynomial in $`|k|`$ and $`|H|`$.

Proof. Due to our assumptions we see that the set $`H`$ as well as the set

$$_H=\{x^1yxy^{}W_K:xK^{(i)},yK_i,y^{}=(x^1yx)^1,i=1,\mathrm{},m1\}$$

are the subsets of the group $`K`$. From the definition of the free product it follows that the elements of $`H`$ are distinct elements of the quotient group $`\overline{K}`$ obtained from $`K`$ by imposing the set $`_H`$ of relations. On the other hand, from the definition of $`_H`$ it follows that any element $`\overline{k}\overline{K}`$ can be represented by an element of $`K`$ of the form $`k_1\mathrm{}k_m`$ for some $`k_iK_i`$, $`i=1,\mathrm{},m`$, moreover this representation is unique, since otherwise the equality of two such representations one could deduce from the relations $`_H`$ which hold in the group $`H`$ as well, but in the group $`H`$ any two such representations differ. After identifying $`\overline{K}`$ with the set of all such elements, we see that the mapping

$$\overline{K}H,\overline{k}k_1\mathrm{}k_m$$

is an isomorphism. Denote by $`Q`$ the composition of the natural epimorpism $`K\overline{K}`$ with this isomorphism. Then the mapping $`Q:KH`$ is an epimorphism and given $`kK`$ the computation of $`Q(k)`$ consists in a reduction of $`k`$ modulo the relations of $`_H`$ to the form $`k_1\mathrm{}k_m`$. This can be done by means of the following procedure

* Step 1. If $`k`$ is the empty word, then output $`Q(k)=k`$.
* Step 2. Using the relations $`x^1yxy^{}_H`$ with arbitrary $`xK^{(1)}`$ and $`yK_1`$ reduce $`k`$ to the form $`k_1\stackrel{~}{k}`$ where $`k_1K_1`$ and $`\stackrel{~}{k}\stackrel{~}{K}`$ with $`\stackrel{~}{K}=K_2\mathrm{}K_m`$.
* Step 3. Applying the procedure recursively to $`\stackrel{~}{k}\stackrel{~}{K}`$ and $`\stackrel{~}{H}=K^{(2)}`$, output $`Q(k)=k_1Q(\stackrel{~}{k})`$.

First, we observe that the length of any intermediate word in the above procedure is at most $`|k|`$. Next, the number of recursive calls (at Step 3) is at most $`|m|`$. Thus the procedure can be done in time polynomial in $`|k|`$ and $`|H|`$. Lemma is proved.

3.4. We are ready to describe the main construction of this section. Let $`K_1,\mathrm{},K_m`$ be a set of pairwise disjoint cyclic groups of prime orders and $`H=\mathrm{\Pi }(K_1,\mathrm{},K_m)`$. Suppose that we are given a homorphic cryptosystem $`(R_i,A_i,P_i)`$ over the group $`K_i`$ with respect to the homomorphism $`f_i:G_iK_i`$ from Section 2, $`i=1,\mathrm{},m`$. Without loss of generality we assume that the groups $`G_i`$ are pairwise disjoint. Set $`G=G_1\mathrm{}G_m`$. Then from the definition of $`f^{}`$ (see formula (2)) and Lemma 3.3 it follows that the mapping $`f=Qf^{}`$ from $`G`$ to $`H`$ is an epimorphism.

###### Theorem 3.4

The triple $`(R,A,P)`$ where $`R`$ is an arbitrary set of distinct representatives of the right cosets of $`G`$ with respect to $`ker(f)`$, provided that $`R`$ fulfils the condition (H2) of Definition 1.1,

$$A=\{(a,r)A^{}\times R^{}:f(r)=1_H\},P:AG,(a,r)\overline{P^{}(a)r}$$

with $`A^{}`$, $`P^{}`$ defined in Subsection 3 and $`R^{}`$ defined in (5), is a homomorphic cryptosystem over the group $`H`$ with respect to the homomorphism $`f:GH`$.

Proof. From the definition of $`R^{}`$ it follows that given $`gG`$ there exist uniquely determined $`g_0ker(f^{})`$ and $`rR^{}`$ such that $`g=g_0r`$. So, $`f(g)=Q(f^{}(g_0)f^{}(r))=Q(f^{}(r))=f(r)`$. Thus $`gker(f)`$ iff $`f(r)=1`$. By Lemma 3.2 this implies that $`im(P)=ker(f)`$. To check the condition (H1) of Definition 1.1 we have to show how to get random elements of $`A`$. This follows from the remarks before Lemma 3.2 for random generating elements of $`A^{}`$, whereas the sets $`R_i`$, $`i=1,\mathrm{},m`$, are given explicitly. Thus it remains to verify that $`P`$ is a trapdoor function (i.e. the condition (H3)).

First, we observe that by statement (i1) of Lemma 3.2 and by Lemma 3.3 the mapping $`P`$ is polynomial time computable. Second, by condition (H3) for homomorphic cryptosystems $`(R_i,A_i,P_i)`$ there exists an algorithm that given $`i\{1,\mathrm{},m\}`$ and $`g_iG_i`$ efficiently finds an element of the set $`P_i^1(g_i)`$. Let us show that it suffices to invert $`P`$. Indeed, in this case given $`gG`$ the element $`f^{}(g)K`$ can be found efficiently. Since $`g=g_0r`$ for uniquely determined $`g_0ker(f^{})`$ and $`rR^{}`$, and $`f^{}(r)=f^{}(g)`$, one can compute the element $`r`$ and hence the element $`g_0=gr^1`$ within the same time. By statement (i2) of Lemma 3.2 we can also find an element $`aA^{}`$ such that $`P^{}(a)=g_0`$. Thus to invert $`P`$ it suffices to test whether $`f(r)=1_H`$ holds (if $`f(r)=1_H`$, then $`(a,r)A`$ and $`P(a,r)=g`$). We have $`f(r)=Q(f^{}(r))`$. Next, by condition (H2) for homomorphic cryptosystems $`(R_i,A_i,P_i)`$ we can find the element $`f^{}(r)`$ and so by Lemma 3.3 the element $`Q(f^{}(r))`$, and finally test the equality $`f(r)=1_H`$.

Suppose that one can invert $`P`$ efficiently. Let $`gG`$. If $`gker(f)`$, then obviously $`gker(f^{})`$. Let now $`gker(f)`$ and $`(a,r)A`$ be a proof for $`g`$, i.e. $`\overline{P^{}(a)r}=g`$. Since $`r`$ belongs to the right transversal $`R^{}`$ of $`ker(f^{})`$ in $`G`$, it follows that $`gker(f^{})`$ iff $`r=1_G`$. Moreover, if $`r=1_G`$, then obviously $`P^{}(a)=g`$. Thus the problem of finding $`(P^{})^1`$ is polynomial time reducible to the problem of finding $`P^1`$. So by statement (i3) of Lemma 3.2 the problem of finding $`P_i^1`$, $`i=1,\mathrm{},m`$, is polynomial time reducible to the problem of finding $`P^1`$. Thus $`P`$ is a trapdoor function which completes the proof.

Observe that one can explicitly produce a set $`R`$ satisfying the condition of Theorem 3.4 (i.e. the condition (H2)). Namely, for each element $`hH`$ find a representation $`h=k_1\mathrm{}k_m`$ where $`k_iK_i`$ (cf. Subsection 3.2), and take $`r_iR_i`$ such that $`f_i(r_i)=k_i`$, $`i=1,\mathrm{},m`$. Then the set of all elements $`r_1\mathrm{}r_m`$ for all $`hH`$ can be chosen as the set $`R`$.

From Theorems 2.2 and 3.4 we immediately obtain the following statement.

###### Corollary 3.5

Let $`K_1,\mathrm{},K_m`$ be a set of pairwise disjoint cyclic groups of prime orders, $`m1`$. Then $`\mathrm{\Pi }(K_1,\mathrm{},K_m)𝒢_{crypt}`$.

3.5. The special cases of a semidirect product are the direct and wreath products. Indeed, in the latter case the resulting group is a semidirect product of the direct power of the first group (with the number of the factors being equal to the order of the second group) by the second group which acts on the product by permutations of direct factors, see e.g. . Thus as an immediate consequence of Theorem 3.4 we conclude that the class $`𝒢_{crypt}`$ contains direct and wreath products of cyclic groups of prime orders (cf. Corollary 3.5). Using this fact we can prove the main result of the paper.

###### Theorem 3.6

Any solvable nonidentity group belongs to the class $`𝒢_{crypt}`$.

Proof. It is a well-known fact that any solvable group can be constructed from a cyclic group of prime order by a sequence of cyclic extensions. On the other hand, from \[10, Theorem 6.2.8\] it follows that any extension of one group by another one is isomorphic to a subgroup of the wreath product of them. So it suffices to verify (cf. Corollary 3.5) that any nonidentity subgroup of a semidirect product of cyclic groups of prime orders belongs to the class $`𝒢_{crypt}`$.

To do this let $`H𝒢_{crypt}`$ be such a group. Then there exists a homomorphic cryptosystem $`𝒮=(R,A,P)`$ over $`H`$ with respect to some epimorphism $`f:GH`$. Without loss of generality we assume that $`𝒮`$ is the homomorphic cryptosystem from Theorem 3.4. Given explicitely a non-identity subgroup $`H^{}`$ of $`H`$ set $`G^{}=f^1(H^{})`$, $`f^{}=f|_G^{}`$ and $`𝒮^{}=(R^{},A^{},P^{})`$ where $`R^{}=Rf^1(H^{})`$,

$$A^{}=\{(a,r)A^{}\times (R^{})^{}:f(r)=1_H\},P^{}:A^{}G^{},(a,r)\overline{P^{}(a)r}$$

and $`(R^{})^{}=\{rR^{}:f(r)H^{}\}`$. Then $`𝒮^{}`$ is a homomorphic cryptosystem over $`H^{}`$ with respect to the homomorphism $`f^{}:G^{}H^{}`$. Indeed, we present the group $`G^{}`$ as a subgroup of $`G`$ generated by the sets $`im(P^{})`$ and $`R^{}`$. (In this presentation of $`G^{}`$ we would be unable to recognize its elements in $`G`$, but we do not need this.) Now the first two conditions of the Definition 1.1 are satisfied for $`𝒮^{}`$ because they are satisfied for $`𝒮`$ (to generate a random element of $`A^{}`$, it suffices to generate a random element $`r^{}`$ of $`(R^{})^{}`$ and for this purpose one can generate a random $`rR^{}`$ and set $`r^{}=\overline{r\stackrel{~}{r}}`$ where $`\stackrel{~}{r}`$ is the element of $`R^{}`$ such that $`f(r)f(\stackrel{~}{r})^1H^{}`$). Since $`ker(f^{})=ker(f)`$, we have $`im(P^{})=ker(f^{})`$ and condition (H3) is also satisfied for $`𝒮^{}`$ because it is satisfied for $`𝒮`$.

It should be remarked that the construction of a homomorphic cryptosystem over a solvable group $`H`$ described in this section is rather theoretical. The computational complexity of the underlying algorithms is bounded by a polynomial the degree of which is a function of $`|H|`$. Besides, the size of representing $`G`$ could be exponential in $`|H|`$ due to involving wreath products. However, it seems that more careful implementation can be developed for small groups.

From Theorem 3.6 it follows that there exists a homomorphic cryptosystem over the group $`Sym(n)`$ for $`n4`$. It would be interesting to construct a homomorphic cryptosystem over an arbitrary symmetric group, because such a system would provide secret computations with any permutation group and moreover an implementation of any boolean circuit in the sense of . In this connection we remark that every boolean circuit of logarithmic depth can be implemented by a polynomial-time computation in an arbitrary nonsolvable group (see ). On the other hand, it was proved in  that over an arbitrary nilpotent group not any boolean circuit can be implemented. If the group is not nilpotent but is solvable, then only an exponential size implementation is known  and it is conjectured that one is unable to do better. Thus, if the latter conjecture was wrong, then combining with Theorem 3.6 would enable us to encrypt any boolean circuit.

Acknowledgements. The authors would like to thank the Max-Planck Institut fuer Mathematik (Bonn) during the stay in which this paper was done; also Eberhard Becker who has drawn the authors attention to homomorphic cryptosystems and Igor Shparlinski for useful discussions.