text:
# Using the DIFF Command for Natural Language Processing

## 1 Introduction

Diff, a software program that detects differences between two sets of data, has numerous applications in natural language processing. In this paper, we first describe diff and then give several examples of how diff can be used in natural language processing. These examples include the detection of differences, the merging of two sets of different data, the extraction of rewriting rules, and the best matching of two different sets of data.111This papar is an English rough translation of our papers .

We summarize the value of this paper below.

* Since diff is one of the basic programs that comes with UNIX systems, it is readily available and very easy to use. As this paper shows, diff is also applicable to various kinds of studies in the field of natural language processing. This combination makes the program particularly valuable.
* Paraphrasing has been become an active area of study . Section 3.2 shows that we can use diff to examine the differences between a spoken language and the corresponding written language and in acquiring the rules for the transformations between the spoken and written languages. Diff is not only useful for such studies, but is also applicable to other kinds of studies in the field of paraphrasing. This paper will provide a basis for such further application of diff.
* In general, diff is used to detect differences between two data. However, it can also be used to merge data and to match data optimally. Section 4 covers several interesting studies of the use of diff in merging and best-matching tasks. These examples are the alignment of a paper and a corresponding presentation and a question-answering system. This paper thus has originality in terms of topics researched and approach to the research.

## 2 Diff and Mdiff

In this section, we describe diff. Diff is a software program of Unix systems that is used to compare files. The program displays line-by-line differences between a pair of text files while retaining the order of the data. For example, suppose we have the following two files:

```
    File 1:               File 2:
    I                     I
    go                    go
    to                    to
    school.               university.

```

When we give these data to diff, the difference is displayed in the following way.

```
    < school.
    > university.

```

Diff has a -D option, which is very useful. When we use diff with this option, common parts as well as differences are displayed. That is, files can be merged by using this option. However, the output of diff -D is in a form which is used for C preprocessor such as “Ifdef” and this is difficult for people to read. Therefore, in this paper, we give the differences in the following way:

```
    ;===== begin =====
    (The parts which only exist in the first file)
    ;-----------------
    (The parts which only exist in the second file)
    ;=====  end  =====

```

where, “`;===== begin =====`” indicates the beginning of the differences, “`;====` `= end =====`” indicates the end of the differences, and “`;-----------------`” indicates the boundary between the two sets of data. In this paper, we refer to diff in the case where the files are merged by using the `-D` option and the differences are displayed in the above form as mdiff.

When we give our earlier pair of files to mdiff, we obtain the following result.

```
    I
    go
    to
    ;===== begin =====
    school.
    ;-----------------
    university.
    ;=====  end  =====

```

“I go to” matchs while “school” and “university” are differences. The output of mdiff is easy to examine and understand because it also displays, unlike diff, the common parts.

We can reproduce the two original files from the output of mdiff. When we take the common part and the upper part of the differences, we obtain the content of the first file. When we take the common part and the lower part of the differences, we obtain the content of the second file. We can reproduce all of the original data in this way.

Since mdiff only displays the common part of the data once, it is able to reduce the amount of data. Since it is possible to fully reproduce the original data from mdiff’s output, we are able to say that mdiff compresses the data while retaining the original information.

Since the output of diff is difficult to read and the output of mdiff contains all information output by diff, we use mdiff for our explanations in the following sections. In the following sections, let’s look at some actual examples of the application of mdiff to natural language processing.

## 3 Detecting differences and acquiring transformation rules

In this section, we describe the studies where we used mdiff to detect differences and then acquired transformation rules from the differences. In more concrete terms, we give the following two examples.

* Detection of differences between the outputs of multiple systems
* Examination of differences to acquire transformation rules

### 3.1 Detection of differences between the outputs of multiple systems

We have been using multiple morphological analyzers to improve the quality of the results of morphological analysis.

Suppose that we use two morphological analyzers to analyze “We like apples.”; the results are as follows.

```
    System 1                 System 2
    We       Noun            We       Noun
    like     Verb            like     Preposition
    apples   Noun            apples   Noun

```

The word, “like” is ambiguous in terms of part of speech (POS) and can take “Verb” or “Preposition”. The POS of “like” in the above sentence is “Verb” and the analysis by System 2 was wrong. Here, when we give the two results to mdiff, we obtain the following results.

```
    We      Noun
    ;===== begin =====
    like    Verb
    ;-----------------
    like    Preposition
    ;=====  end  =====
    apples  Noun

```

Mdiff makes it easy for us to detect differences between the results produced by multiple systems. In this case, we notice that there are differences in the line for “like”. Here, we determine beforehand that a worker corrects the outputs and judges which is correct when such differences are detected. When the output of the first system is correct he does nothing and when the output of the second system is correct he marks the beginning of “;—————–” with an “x”. When we do so, we can automatically select the better result of each difference from the results marked by the worker. We simply delete the lower part of the difference where there is no “x” and delete the upper part of the difference when there is an “x”. As a result, we produce a more accurate overall result than either of the two original results. There are many cases when both parts of a difference are incorrect. In such cases, we manually rewrite the data correctly above “;—————–”.

When we use this method, the only incorrect data that is not corrected are the data on which both the systems have performed exactly the same wrong analysis. Many errors can thus be corrected. What we want to say here is that we must prepare multiple systems that have different characteristics. If we use two systems that often produce the same wrong analysis, we will fail to correct a large number of errors.

Diff3, which is another unix program, is available for analysing the outpus of three systems. Diff3 detects differences among three files.

Although we have given an example of morphological analysis, we can also use mdiff to detect differences between the results of other forms of analysis as long as the results are expressed in a line-by-line form.

Although we have described a study where we detected differences between the outputs of multiple systems, we are also able to detect and correct errors in tagged corpora by comparing the tagged corpora and the results of analysis of the corpora by a certain system.222There are studies on the correction of errors in tagged corpora .

### 3.2 Examining differences and acquiring transformation rules

In this section, we describe our study of the use of mdiff to detect differences between a spoken language and the corresponding written language. In this study, we used the alignment between data on the spoken language and on the written language to examine the differences between the spoken and written languages on the basis of the results of mdiff. We thus acquired rules for transforming the spoken language to the written language along with the inverse transformation rules. We used presentations at academic conferences as examples of the spoken language and used papers which had the same content as the presentations as examples of the writtten language.

For example, we supposed that data in the spoken and written languages are given as in Table 1.333Although we have given English sentences in the table, we actually performed these experiments on Japanese-language materials. The data in the table are transformed so that each line has one word before we apply mdiff. Applying mdiff to the data produces the results shown in Table 2. By extracting the differences from these results, we obtain the results shown in Table 3.

The results show us that “uh” is inserted in the spoken language and that “performed” can, in this case, be paraphrased as “done”. We can use mdiff to detect differences between the spoken language and the written language and then examine them linguistically. We can also consider the detected differences as rules for the transformation between written and spoken languages. For example, the occurrence of “uh” can be regarded as indicating a rule that we may insert “uh” into the written-language data as part of transforming it into the spoken language. The difference in terms of whether “performed” or “done” is used can be regarded as a rule that we transform “performed” to “done” when we transform data from the written to the spoken language. We thus find that we are able to use mdiff to detect the transformation rules.

In this section, we used data in the written and spoken language. However, we can study a variety of data in a similar way. For example, when we apply mdiff to written texts before and after spelling and grammar are checked, we can find out how we should modify the sentences to eliminate spelling and grammatical errors and thus acquire rules for checking spelling and grammar. When we use mdiff on texts and the results of their automatic summarization, we can clearly see how the texts were summarized and can thus acquire the rules used in automatic summarization. In other cases of paired textual data, too, we will be able to apply mdiff to examine the data in the various ways that such pairs can be obtained and to obtain the rules for transformation.

## 4 The merging and best matching of two different sets of data

In this section, we describe studies of the use of diff’s function of optimally merging data.444In diff, best matching is performed by maximizing the extent of the common parts. We describe the following two examples.

* The alignment of presentations and corresponding papers
* A question-answering system using the function of finding the the best match

### 4.1 Alignment of a paper with the corresponding presentation

In this section, we align a paper with the corresponding presentation . The paper and corresponding presentation is the same example as was used in describing the acquisition of the rules for transformation from written to spoken-language data (Section 3.2). The presentation was made at an academic conference and the paper corresponds to this presentation. When we are able to align each part of a presentation with a corresponding part of the paper, this form of alignment is very useful. For example, when we listen to a presentation, we are able to refer to the part of the paper that corresponds to the part of we are now listening. When we read a paper, we are able to refer to the part of the presentation that corresponds to the part of the paper we are now reading . In this section, we consider the use of mdiff to align a paper with the corresponding presentation.

Here, let’s try to determine the parts of the presentation to which each chapter of the paper corresponds by using mdiff. We suppose that the content is laid out in the same order in the paper and in the presentation. In advance, we place symbols such as `<Chapter 1>`, as shown in Figure 1, into the paper. This is so that we are easily able to recognize the chapters of the paper. By applying mdiff to the data set that consists of a presentation after both have been transformed so that each line has one word, we obtained the results shown in Figure 2. Next, we obtained the results shown in Figure 3 by eliminating the upper parts of the differences, i.e., those that correspond to the paper, and leaving symbols such as `<Chapter 1>` in place. The symbols such as `<Chapter 1>` are only inserted in the data of the presentation. We are then able to recognize the chapter of the paper that corresponds with a given part of the presentation.

To put this simply, we place information to indicate chapters in the data of presentations. We then use the merging function of mdiff to match the parts of the paper and presentation. Then, by eliminating the content of the paper, we are left with information on the chapters. We can easily align a paper with a corresponding presentation at the level of chapters by using mdiff.

### 4.2 A question-answering system using best matching

In this section, we describe our question-answering system that utilizes mdiff’s best-matching function .

A question-answering system receives a question and outputs an answer. For example, when “Where is the capital of Japan?” is given, it returns “Tokyo”. When we consider now much knowledge is in written natural language form, we see that the system only has to match a given question with a sentence (a knowledge sentence) that includes the knowledge then select the word in the sentence that corresponds to the interrogative pronoun of the question sentence. For example, in the case of the previous question, the system detects this knowledge sentence in a knowledge database: “Tokyo is the capital of Japan.” “Tokyo” is output as the answer since corresponds to the interrogative pronoun in the sentence. Here, we consider the use of mdiff in answering questions.

Firstly, we transform the interrogative pronoun in the question sentence into X and change “?” into “.” at the end of the sentence. We then obtain “X is the capital of Japan.”. We thus obtain “Tokyo is the capital of Japan.” from the knowledge base. We obtain the results in Figure 4 (a) by applying mdiff to the two sentences. We consider the part that is paired with X in the differences as the answer and correctly obtain “Tokyo”.

Even if there are differences between the compared sentences, we can still use mdiff to correctly extract the answer. For example, suppose that the knowledge sentence is “Tokyo is the capital in Japan.” In this case, the result of mdiff is as in Figure 4 (b). Although the number of differences increased, the part that corresponds to X remains “Tokyo” and the correct answer has been detected.

The question-answering system we have proposed repeats to transform the question sentence and knowledge sentence so that the two sentences are as similar as possible. It then extracts the answer by matching the sentences when they are at this point. Since we use the similarity between two sentences in that time, we have to define the degree of the similarity between sentences. Since we can recognize the common parts and the different parts by using mdiff, we can calculate the degree of similarity as (the number of characters in the common part)/(the number of all characters).555Here, we use mdiff to calculate the degree of similarity of the sentences. Mdiff is also useful in this way. Here, if we suppose that we have a rule for transforming “in” to “of”, we match the data after transforming “Tokyo is the capital in Japan.” to “Tokyo is the capital of Japan.” and we are able to obtain the answer more reliably by decreasing the number of differences.

## 5 Conclusion

In this paper, we describe a lot of examples of the use of diff in various problems of natural language processing. In Section 3.1, we describe how we apply diff to a system that combines multiple systems and thus obtain higher precisions than from any of individual systems. Studies on combination of systems are very well known. We describe how such combination can easily be performed by diff. In Section 3.2, we showed that diff is easily able to examine the differences between a spoken and written language and acquire rules for the transformation between the spoken and written language. We handled the paraphrasing of written language into spoken language. However, the study on paraphrasing covers a wide range: automatic extraction of compressed sentences, sentence polishing-up to modify sentences correctly, and transformations from difficult sentences to plain sentences. In these studies, too, it will be easy to apply diff to various forms of examination and to the extraction of various transformation rules. This paper will become a basis for studies of paraphrasing, an area in which there has been increasing activity. In Section 4, we showed examples of the merging of data and the optimal matching of data. Since diff is generally used to detect differences, Section 4, which shows applications of diff to merging and matching, presents a new and original approach. In the section, we showed that it is easy to apply diff to interesting studies of two kinds, the alignment of a paper with a corresponding presentation and a question-answering system.

In this paper, we have described many examples of the application of diff to natural language processing. We hope that diff will be applied to an ever-widening range of studies.

## Acknowlegement

We were given data for the experiments on alignment of a paper with a corresponding presentation by Kiyotaka Uchimoto of the Communications Research Laboratory.


format: text

text:
# Population synthesis of wide binary millisecond pulsars

## 1 Introduction

In close binaries containing a neutron star, the evolution of the system is driven by the nuclear evolution of the companion and the loss of angular momentum through processes as magnetic braking and gravitational radiation. When the combination of these processes causes the companion to fill its Roche lobe, mass and angular momentum are transferred to the neutron star which is then spun up to form a recycled pulsar, unless transient behaviour prevents effective accretion. At the end of the mass-transfer phase, the binary emerges as a binary millisecond pulsar (BMSP) consisting of a rapidly rotating neutron star in a nearly circular orbit around a low-mass white dwarf, which used to be the core of the Roche-lobe overflowing companion (for a review see Bhattacharya & van den Heuvel 1991, Phinney & Kulkarni 1994). If the mass transfer phase takes place during the companion’s ascent of the red giant branch, the relation between the core mass and the radius of the star and the relation between the Roche-lobe radius and the orbital separation, lead to a correlation between the orbital period of the BMSP and the mass of the white dwarf. This correlation has been studied extensively by Joss, Rappaport & Lewis (1987), Savonije (1987), Rappaport et al. (1995), Ritter (1999), and Tauris & Savonije (1999).

The precise evolutionary path followed by a BMSP progenitor depends mainly on the orbital period and on the mass of the donor star at the onset of Roche-lobe overflow after the supernova explosion of the primary. Systems with orbital periods longer than $`20`$ days and donors more massive than $`2.0M_{}`$ are likely to evolve through a common-envelope phase during which the envelope of the donor is expelled and the neutron star spirals in towards the donor star’s core (Kalogera & Webbink 1998). If the orbital period is shorter, the common-envelope phase can be avoided for donor stars up to $`4.0M_{}`$ (Tauris, van den Heuvel & Savonije 2000; Kolb et al. 2000; Podsiadlowski, Rappaport & Pfahl 2002). Binaries with donor stars less massive than $`2.0M_{}`$ will become low-mass X-ray binaries, for which the evolution depends on whether the orbital period is longer or shorter than the bifurcation period separating diverging from converging systems (Pylyser & Savonije 1988, 1989). Only low-mass X-ray binaries with orbital periods longer than the bifurcation period will eventually evolve into wide BMSPs.

At present, 33 BMSPs with orbital periods longer than one day have been found in the Galactic disc. Their orbital periods, eccentricities, and white dwarf masses are listed in Table 1 (Ritter, private communication). Inspection of the observed orbital periods shows a paucity of systems with periods between 30 and 60 days, and an absence of systems with periods longer than 200 days. A theoretical interpretation of the possible period gap was given by Tauris (1996) and Taam, King & Ritter (2000), while an explanation for the upper limit was put forward by Ritter & King (2002).

With this paper, we initiate a systematic study based on population synthesis techniques to investigate the different evolutionary channels leading to the formation of wide BMSPs and to compare the theoretical with the observed orbital period distribution of wide BMSPs. We particularly investigate whether a set of standard evolutionary parameters can be found to reproduce the observed orbital period distribution without including observational selection effects or pulsar lifetime issues. Due to the large uncertainties involved in the determination of the masses of the white dwarfs in the observed BMSPs (see Table 1) we presently do not consider the white dwarf mass distribution as good a tool as the orbital period distribution to confine stellar and binary evolution and formation parameters.

Since BMSPs located in globular clusters are thought to be formed through tidal encounters of neutron stars with primordial binaries (e.g. Rappaport, Putney & Verbunt 1989), we limit ourselves to BMSPs in the Galactic disc. We also leave aside systems with orbital periods shorter than 1 day since they are likely to follow evolutionary paths different from those of the BMSPs investigated here. In particular, the evolution of the short-period systems is dominated by angular momentum losses for which quantitative descriptions are still uncertain (Ergma & Sarna 1996; Kalogera, Kolb & King 1998; Ergma, Sarna & Antipova 1998; Rasio, Pfahl & Rappaport 2000; Dewi et al. 2002; and in particular Podsiadlowski et al. 2002 and references therein).

The plan of the paper is as follows. In Sect. 2, we briefly summarise the basic concepts of the binary evolution code used in this investigation. In Sect. 3, the evolutionary channels leading to the formation of wide BMSPs are discussed. In Sect. 4, we explore the parameter space occupied by the BMSP progenitors during various phases of their evolution. The effects of the input parameters adopted in our binary evolution calculations on the orbital period distribution of wide BMSPs and the relative contributions of the different formation channels to the total population of wide BMSPs are investigated in Sect. 5. The final section is devoted to concluding remarks.

## 2 The binary evolution code

### 2.1 Single star evolution

We developed a rapid binary evolution code based on the analytical approximations for the evolution of single stars with masses ranging from $`0.1`$ to $`100M_{}`$ and metallicities between $`Z=0.0001`$ and $`Z=0.03`$ derived by Hurley, Pols & Tout (2000). The formulae express the luminosity, radius, and core mass of a star as a function of its age, mass, and metallicity for all phases from the zero-age main sequence up to and including the remnant stages. For massive stars ending their life as a neutron star, we simplified the prescription for the mass of the remnant so that all neutron stars are born with a mass of $`1.4M_{}`$. The simplification is justified by the fact that most neutron stars with accurately known masses have a mass close to this value (Thorsett & Chakrabarty 1999). We furthermore limit ourselves to Population I stellar compositions.

The single star evolution scheme takes into account mass loss from the envelope due to stellar winds. Stars beyond the main sequence are subjected to the mass-loss rates given by Kudritzki & Reimers (1978), the superwind phase on the asymptotic giant branch (AGB) is modelled using the prescription of Vassiliadis & Wood (1993), massive stars are subjected to mass loss over the entire Hertzsprung-Russell diagram as outlined by Nieuwenhuizen & De Jager (1990) and Kudritzki et al. (1989), and naked helium stars are subjected to a Wolf-Rayet-type mass-loss rate. For more details, we refer to Hurley et al. (2000).

In addition to the stellar evolution of the binary components, we also take into account the evolution of the binary configuration due to gravitational radiation, magnetic braking, mass loss and accretion from stellar winds, and conservative or non-conservative Roche-lobe overflow. For the implementation of magnetic braking, we adopt a Skumanich-type parametrisation modified to take into account the dependency on the mass of the convective envelope as described by Hurley et al. (2000). We furthermore assume the orbit of the binary to be circular and keep the rotation of the stars synchronised with their orbital motion at all times. The recipe followed for the implementation of the binary evolution is essentially along the lines laid down by Hurley, Tout & Pols (2002). In the following paragraphs we therefore restrict ourselves to summarising the points where our treatment deviates from the outline given by these authors.

### 2.2 Roche-lobe overflow

When the combined effects of stellar and orbital evolution cause a component of a close binary to become larger than its Roche lobe, mass is transferred from the Roche-lobe overflowing star to the companion. The stability of the mass transfer process depends on the adiabatic and thermal radius-mass exponents, $`\zeta _{\mathrm{ad}}`$ and $`\zeta _{\mathrm{th}}`$, and the Roche-lobe index $`\zeta _\mathrm{L}`$ (for a definition see, e.g., Webbink 1985). For stars on the main sequence or in the Hertzsprung gap, we use the tabulated values for $`\zeta _{\mathrm{ad}}`$ and $`\zeta _{\mathrm{th}}`$ given by Hjellming (1989). In all other cases we adopt the prescription given by Hurley et al. (2002). For the determination of the Roche-lobe index $`\zeta _\mathrm{L}`$, we follow the procedure outlined by Kolb et al. (2001).

If mass transfer proceeds on either the nuclear ($`\zeta _\mathrm{L}<\zeta _{\mathrm{ad}}`$ and $`\zeta _\mathrm{L}<\zeta _{\mathrm{eq}}`$) or the thermal ($`\zeta _{\mathrm{eq}}<\zeta _\mathrm{L}<\zeta _{\mathrm{ad}}`$) timescale of the Roche-lobe overflowing star, the companion is assumed to accrete a fraction $`1\gamma `$ of the transferred mass determined by the ratio of the mass-transfer timescale to the accretor’s thermal timescale (for details see Hurley et al. 2002). The remaining fraction is assumed to be expelled from the system, carrying away the specific orbital angular momentum of the companion. In the particular case of a neutron star accretor, an average accretion rate is determined by imposing an upper limit of $`\left(\mathrm{\Delta }M_{\mathrm{NS}}\right)_{\mathrm{max}}`$ on the amount of mass that can be accreted (for details see Appendix A). Once this limit has been reached, mass accretion onto the neutron star is switched off completely to ensure that the neutron star accretes no more than $`\left(\mathrm{\Delta }M_{\mathrm{NS}}\right)_{\mathrm{max}}`$. We impose this upper limit to mimic the low accretion efficiency neutron stars are thought to have (e.g. King & Ritter 1999, Ritter & King 2002). In addition, we also limit the mass-accretion rate to the Eddington rate at all times.

In the case of dynamically unstable mass transfer ($`\zeta _{\mathrm{ad}}<\zeta _\mathrm{L}`$) from a giant-like star, the binary is assumed to enter a common-envelope phase during which the orbital separation is reduced and the envelope is ejected from the system. The phase is modelled by equating the binding energy of the giant’s envelope to the change in orbital energy of the two stars:

$$\frac{G\left(M_\mathrm{c}+M_\mathrm{e}\right)M_\mathrm{e}}{\lambda R_\mathrm{L}}=\alpha [\frac{GM_\mathrm{c}M_2}{2a_\mathrm{f}}\frac{G\left(M_\mathrm{c}+M_\mathrm{e}\right)M_2}{2a_\mathrm{i}}].$$
(1)
In this equation, $`G`$ is the gravitational constant, $`M_\mathrm{c}`$ and $`M_\mathrm{e}`$ are the core and envelope mass of the Roche-lobe overflowing star, $`R_\mathrm{L}`$ is the radius of the donor star’s Roche lobe, $`M_2`$ is the mass of the companion, $`a_\mathrm{i}`$ and $`a_\mathrm{f}`$ are the orbital separations of the binary at the start and at the end of the common-envelope phase, $`\lambda `$ is a dimensionless structure parameter determining the binding energy of the giant’s envelope, and $`\alpha `$ is the fraction of the orbital energy that is transferred to the envelope. The parameters $`\lambda `$ and $`\alpha `$ are usually assumed to be constant and equal to 0.5 and 1.0, respectively (e.g. de Kool 1992, Politano 1996). In fact, the value of $`\lambda `$ depends on the evolutionary stage of the donor star at the onset of the mass-transfer phase (Dewi & Tauris 2000, 2001). However, since no elaborate set of $`\lambda `$-values has been published for the full range of stellar masses considered in our investigation, we here treat $`\lambda `$ as a constant.

### 2.3 Supernova explosions

Whenever a component of a close binary undergoes a supernova explosion leading to the birth of a neutron star, we use the expressions established by Kalogera (1996) to derive the post-supernova orbital parameters. If the birth of the neutron star is accompanied by a kick, we assume its magnitude $`v_{\mathrm{kick}}`$ to be distributed according to the Maxwellian distribution

$$P\left(v_{\mathrm{kick}}\right)=\sqrt{\frac{2}{\pi }}\frac{v_{\mathrm{kick}}^2}{\sigma ^3}\mathrm{exp}\left(\frac{v_{\mathrm{kick}}^2}{2\sigma ^2}\right),$$
(2)
where $`\sigma `$ is the dispersion of the kick velocities. If the supernova explosion causes the binary to become eccentric, we instantaneously circularise the orbit while conserving the total orbital angular momentum of the binary.

## 3 Evolutionary channels

Before proceeding to the statistics of the BMSP population and the effects of the various assumptions made in our treatment of stellar and binary evolution and formation, we explore the different evolutionary channels leading to the formation of wide BMSPs.

We followed the evolution of a large number of binaries starting with zero-age main sequence components up to a maximum evolutionary age of 15 Gyr. The initial masses of the binary components are taken in the interval between $`0.1M_{}`$ and $`60M_{}`$, and the initial orbital periods in the range from 10 to 10 000 days. We used equidistantly spaced grids consisting of 40 grid points for the logarithm of the initial masses and 200 grid points for the logarithm of the initial periods. For symmetry reasons only binaries with $`M_1>M_2`$ are evolved. Each time a supernova explosion leads to the birth of a neutron star, 200 random kick velocities are generated and the evolution of the surviving binaries is followed until the imposed age limit of 15 Gyr. In order to identify as many formation channels as possible, we do not apply any weighting to the initial masses and the initial orbital periods at this stage. The weighted contributions of the different formation channels to the population of wide BMSPs will be discussed in Sect. 5.

For the results presented in this section, a standard set of input parameters was adopted for the binary evolution calculations: common-envelope phases are modelled using the values $`\lambda =0.5`$ and $`\alpha =1.0`$ for the binding-energy parameter and the envelope-ejection efficiency, the occurrence of a delayed dynamical instability is determined using the values tabulated by Hjellming (1989), supernova kicks are generated with a Maxwellian velocity dispersion $`\sigma =190\mathrm{km}/\mathrm{s}`$, and an upper limit $`\left(\mathrm{\Delta }M_{\mathrm{NS}}\right)_{\mathrm{max}}=0.2M_{}`$ is imposed on the mass that can be accreted by a neutron star. These assumptions will serve as a reference for the population synthesis calculations presented in Sect. 5, where we will refer to them as model A.

We find four major evolutionary channels leading to the formation of wide BMSPs. In three of these channels, the primary undergoes a common-envelope phase (CE) which ejects the primary’s envelope from the system and reduces the orbital separation. For primaries more massive than $`20M_{}`$ the common-envelope phase takes place when the primary crosses the Hertzsprung gap, while for primaries less massive than $`20M_{}`$ the common-envelope phase takes place when the primary ascends the AGB. Both cases lead to similar evolutionary paths. For the continuation of the evolution, it is necessary to distinguish between binaries with low- and intermediate mass secondaries.

The main evolutionary phases of a BMSP progenitor with a low-mass secondary are illustrated in Fig. 1, where the initial orbital period is 1400 days and the binary components have initial masses of $`15M_{}`$ and $`1M_{}`$. The primary fills its Roche lobe after 14 Myr, at which time a stellar wind has reduced the mass of the primary to $`10.3M_{}`$ and increased the orbital period to $`2690`$ days. The primary exits the common-envelope phase as a massive naked helium star which evolves further until it explodes in a type Ib/c supernova and gives birth to a $`1.4M_{}`$ neutron star. After circularisation of the post-supernova orbit, the binary consists of a $`1.0M_{}`$ main-sequence star which orbits the neutron star in a period of 5.7 days. The binary subsequently evolves as a detached system until the secondary expands and fills its Roche lobe during its ascent of the giant branch. The resulting mass-transfer phase is both dynamically and thermally stable and takes place on the nuclear timescale of the giant star. During this phase, mass accretion spins up the neutron star to high rotation rates until it becomes a recycled pulsar. At the end of the mass-transfer phase, the binary emerges as a wide BMSP consisting of a rapidly rotating neutron star in a 64 day orbit around a $`0.3M_{}`$ helium white dwarf, the former core of the giant star. This formation channel has been discussed previously by Bhattacharya & van den Heuvel (1991) and Phinney & Kulkarni (1994), among others. We will refer to it with the label He N.

For BMSP progenitors with intermediate-mass secondaries, the formation channel becomes a bit more complex. Two possible evolutionary paths are displayed in Figs. 2 and 3, in the case of a binary with an initial orbital period of 1400 days and initial component masses of $`30M_{}`$ and $`3.5M_{}`$. Again, the primary emerges from the common-envelope phase as a massive naked helium star which gives birth to a $`1.4M_{}`$ neutron star in a type Ib/c supernova explosion. The further evolution depends on the circularised post-supernova orbital separation.

If the post-supernova orbital period is smaller than $`2.5`$ days (Fig. 2), the binary enters a case A thermal-timescale mass-transfer phase after 214 Myr and detaches again after 217 Myr when the donor mass has decreased to $`1.2M_{}`$. This short thermal-timescale mass-transfer phase is followed by a stable nuclear timescale mass-transfer phase after 939 Myr when the secondary approaches the end of the main sequence. This stable mass transfer phase continues all the way up to the giant branch and leads to the formation of a binary consisting of a rapidly rotating neutron star in a 10 day orbit around a $`0.2M_{}`$ helium white dwarf. We label this formation channel He T.

If the post-supernova orbital period is larger than $`2.5`$ days (Fig. 3), the secondary enters a case B thermal-timescale mass-transfer phase when it crosses the Hertzsprung gap. Due to the high mass-transfer rate, the neutron star rejects most of the mass transferred by the companion so that the spin-up process is probably less effective and thus leads to more slowly rotating pulsars than the spin-up process during a slow stable mass-transfer phase (e.g. Li 2002, and references therein). The semi-detached state ends when the hydrogen envelope of the secondary is completely stripped away and its helium core is exposed as a low-mass naked helium star. During the following 135 Myr, the secondary burns helium in the core while a Wolf-Rayet-type stellar wind slowly removes the outer helium layers. The system ends up as a $`0.5M_{}`$ carbon/oxygen white dwarf in a 10.4 day orbit around a neutron star which may or may not be spun up to millisecond periods. This channel was first found independently by King & Ritter (1999) and by Podsiadlowski & Rappaport (2000) \[see also Tauris et al. (2000)\]. A possible BMSP progenitor that is currently at the end of the early case B mass transfer phase in this channel is the persistent X-ray binary Cygnus X-2 (King & Ritter 1999, Kolb et al. 2000, Podsiadlowski & Rappaport 2000). We will refer to this formation channel with the label CO T.

The fourth and final evolutionary channel for the formation of wide BMSPs corresponds to the direct supernova mechanism discussed by Kalogera (1998). It applies to binaries with low-mass secondaries and initial orbital separations that are wide enough to avoid any kind of Roche-lobe overflow before the supernova explosion of the primary. An example of such a system is displayed in Fig. 4, where the initial orbital period of the binary is 1600 days and the initial masses of the primary and the secondary are $`9M_{}`$ and $`1M_{}`$, respectively. The wide orbit allows the primary to evolve all the way up to the AGB where it explodes in a type II supernova. After the circularisation of the post-supernova orbit, the binary consists of a $`1.4M_{}`$ neutron star orbiting a fairly unevolved main-sequence star with a long orbital period of 1100 days. The secondary then in turn evolves up to the AGB where it fills its Roche lobe after 12.4 Gyr. At this point, its mass is reduced to $`0.7M_{}`$ due to strong stellar winds on the giant branch and the AGB. The following mass-transfer phase is dynamically and thermally stable, and lasts for about 0.3 Myr. At the end of the mass-transfer phase, the binary consists of a rapidly rotating neutron star in a 1440 day orbit around a $`0.5M_{}`$ carbon/oxygen white dwarf. We label this formation channel CO N.

## 4 Parameter space

Next, we examine the parameter space occupied by systems forming wide BMSPs at the start of their evolution, at the instant just before the supernova explosion of the primary, at the onset of Roche-lobe overflow after the supernova explosion, and at the birth of the BMSP.

The initial orbital periods $`P_{\mathrm{orb}}`$ and secondary masses $`M_2`$ of binaries that will evolve into BMSPs are shown in the left-hand panels of Fig. 5. Dark and light dots indicate systems evolving into BMSPs containing a helium white dwarf and systems evolving into BMSPs containing a carbon/oxygen white dwarf, respectively, while big and small dots distinguish between systems undergoing a nuclear or thermal-timescale mass-transfer phase after the supernova explosion of the primary. The bulk of the BMSP progenitors occupy a rather narrow band of orbital periods around 1000 days. The upper limit of the band corresponds to the longest orbital periods for which the primary fills its Roche lobe and undergoes a common-envelope phase prior to its supernova explosion, while the lower limit is formed by the shortest orbital periods for which mergers or contact configurations can be avoided. Besides the band around 1000 days, there is also a small group of systems with orbital periods in the range from 10 to 30 days. These systems evolve along the same lines as illustrated in Figs. 2 (the He T channel) and 3 (the CO T channel), except that the primary undergoes a rapid thermal-timescale mass-transfer phase instead of a common-envelope phase prior to its supernova explosion. Since these systems do not contribute significantly to the total population of BMSPs we do not consider them separately in our present investigation.

The right-hand panels of Fig. 5 show the orbital periods and secondary masses of the BMSP progenitors just before the supernova explosion of the primary. With the exception of a few rare cases, the orbital separations of the initial binaries are substantially reduced due to the common-envelope phase of the primary. The few systems managing to avoid the common-envelope phase are those evolving through the direct supernova mechanism (Fig. 4, the CO N channel). Their orbital periods have increased slightly with respect to their initial orbital periods due to mass loss from the system caused by the stellar wind of the primary. The systems evolving through the direct supernova mechanism into BMSPs containing a helium white dwarf are included in the He N channel, since their post-supernova evolutionary path is similar to that shown in Fig. 1.

In Fig. 6, we show the parameter space occupied by the progenitors of the wide BMSPs after the supernova explosion of the primary, at the onset of Roche-lobe overflow from the secondary. The lines labelled ZAMS, TMS, and BGB represent the orbital periods at which the secondary would fill its Roche lobe if it was located at the zero-age main sequence, the terminal main sequence, or the base of the giant branch, respectively. The thick lines delimit the different stability regions for mass transfer from the secondary. The dark grey region indicates dynamically unstable or delayed dynamically unstable mass transfer, the light grey region indicates thermal-timescale mass transfer, and the intermediate grey region indicates thermally and dynamically stable mass transfer. For the construction of both thin and thick lines, the mass of the neutron star was assumed to be $`1.4M_{}`$.

The post-supernova orbital periods of the BMSP progenitors depend on the magnitude of the kick velocity imparted to the neutron star at birth. For secondaries with a mass larger than $`1.3M_{}`$ the dynamical instability of mass transfer sets an upper limit for the post-supernova orbital period. The lower limit on the orbital periods observed in Fig. 6 is a consequence of the lower limit on the post-common-envelope orbital separations resulting from the requirement that the systems exit the common-envelope phase as detached binaries. Lower and upper limits on the mass of the secondary are caused by the imposed age limit of 15 Gyr and by the delayed dynamical instability. Secondaries with a mass smaller than $`0.9M_{}`$ evolve too slowly to fill their Roche lobe within 15 Gyr, while for secondaries with a mass larger than $`4.04.5M_{}`$ the initial thermal-timescale mass-transfer phase will evolve into a dynamically unstable phase shortly after the onset of Roche-lobe overflow.

Furthermore, four groups of systems can be distinguished in Fig. 6, in accordance with the four identified formation channels. Systems evolving into BMSPs containing a helium white dwarf either have low-mass secondaries which enter a stable mass-transfer phase when they cross the Hertzsprung gap or ascend the giant branch (Fig. 1, the He N channel), or intermediate-mass secondaries for which mass transfer initially takes place on a thermal timescale (Fig. 2, the He T channel). Systems evolving into BMSPs containing a carbon/oxygen white dwarf undergo a thermal-timescale mass-transfer phase when the secondary crosses the Hertzsprung gap (Fig. 3, the CO T channel) or a slow stable mass-transfer phase on the AGB if they manage to avoid the common-envelope phase of the primary prior to the supernova explosion (Fig. 4, the CO N channel).

The parameter space occupied by the BMSPs in the $`\mathrm{log}P_{\mathrm{orb}}`$ \- $`M_{\mathrm{WD}}`$ plane, where $`M_{\mathrm{WD}}`$ is the mass of the white dwarf, at the time of their formation is displayed in Fig. 7. The positions of the observed BMSPs listed in Table 1 are indicated by diamonds. The arrows indicate the uncertainties on the observed white dwarf masses. The orbital periods of the BMSPs containing helium white dwarfs are correlated with the masses of the white dwarfs. The correlation occurs for both the evolutionary channels He N and He T, and arises due to the relation between the core mass and the radius of the giant which is the progenitor of the white dwarf, and due to the relation between the Roche-lobe radius of the giant and the semi-major axis of the orbit (e.g., Rappaport et al. 1995, Tauris & Savonije 1999). The BMSPs forming through the evolutionary channel CO T are located below the systems containing helium white dwarfs, while the ones forming through the direct supernova mechanism CO N form an extension of the $`P_{\mathrm{orb}}`$ \- $`M_{\mathrm{WD}}`$ relation for systems containing a helium white dwarf.

We conclude this section by noting that the ranges of orbital periods and secondary masses occupied by the BMSP progenitors at the onset of Roche-lobe overflow after the supernova explosion of the primary and at the formation of the BMSP found from our population synthesis code are in good agreement with those marked out by Tauris et al. (2000) and Podsiadlowski et al. (2002) on the basis of detailed numerical calculations.

## 5 Population synthesis

### 5.1 Model assumptions

Now that we have established the different formation channels leading to the birth of wide BMSPs and the parameter space occupied by their progenitors, we turn our attention to the distribution of the final BMSP orbital periods. In particular, we investigate how the distribution changes with the assumptions made in the binary evolution calculations. We therefore repeated the calculations discussed in the previous sections for a variety of models consisting of different input parameters. Our reference model, model A, corresponds to the input parameters adopted in Sects. 3 and 4. In models CE1 and CE2, the common-envelope efficiency parameter $`\alpha `$ is varied. The role of supernova kicks is explored in models K0, KM50, KM100, and KM400: in model K0 no kicks are imparted to the neutron star at birth, while in models KM50, KM100, and KM400, the Maxwellian velocity dispersion takes the values 50, 100, and 400 km/s, respectively. In model NS2, the maximum amount of mass that can be accreted by a neutron star is increased from $`0.2M_{}`$ to $`0.6M_{}`$. The influence of the critical mass ratio $`q_{\mathrm{c},\mathrm{ddi}}`$ for the development of a delayed dynamical instability is investigated in models QCDD and QCDD2. In model C, finally, the limiting case of conservative mass transfer is considered. Details of the differences between the various models are summarised in Table 2.

Each time the evolution of an initial binary leads to the formation of a wide BMSP, its contribution to the orbital period distribution is weighted according to the probability density distributions of its initial parameters. The initial mass of the primary (the initially most massive star) is assumed to be distributed according to the normalised initial mass function

$$\xi \left(M_1\right)=\{\begin{array}{cc}0\hfill & M_1/M_{}<0.1,\hfill \\ 0.38415M_1^1\hfill & 0.1M_1/M_{}<0.75,\hfill \\ 0.23556M_1^{2.7}\hfill & 0.75M_1/M_{}<\mathrm{},\hfill \end{array}$$
(3)
so that the probability that a star has a mass between $`M_1`$ and $`M_1+dM_1`$ is given by $`\xi \left(M_1\right)dM_1`$. This distribution function is a simplified version of the initial mass function used by Hurley et al. (2002) and is based on the initial mass function derived by Kroupa, Tout & Gilmore (1993). We recall that for our purpose, we limit ourselves to masses between $`0.1M_{}`$ and $`60M_{}`$.

For a given mass of the primary, the initial mass of the secondary (the initially least massive star) is determined by the mass ratio $`q=M_2/M_1`$. We assume the latter to be distributed according to

$$n(q)=\{\begin{array}{cc}1\hfill & 0<q1,\hfill \\ 0\hfill & q>1,\hfill \end{array}$$
(4)
so that the masses of the component stars are correlated with each other. The existence of such a correlation was suggested, among others, by Eggleton, Fitchett & Tout (1989).

The distribution of the initial orbital separations is assumed to be uniform in the logarithm of the semi-major axis $`a`$. Following Hurley et al. (2002), we adopt the normalised distribution

$$\chi (a)=\{\begin{array}{cc}0\hfill & a/R_{}<3\text{ or }a/R_{}>10^4,\hfill \\ 0.12328a^1\hfill & 3a/R_{}10^4.\hfill \end{array}$$
(5)
After weighting the contribution of each newly formed BMSP to the total population, the resulting distribution function is convolved with a constant star-formation rate, on the assumption that the BMSPs under consideration do not evolve any further once they are formed. In particular, we presently do not take into account the spin-down time and thus the finite lifetime of the recycled pulsars. After the convolution, the distribution function is normalised so that the integral over all systems found is equal to one.

### 5.2 The orbital period distribution

We considered the distribution of the orbital periods separately for each of the evolutionary channels discussed in Sect. 3. The most illustrative results are shown in Fig. 8, where the contribution of each channel to the orbital period distribution is displayed using a different shade of grey. For display purposes, the distribution functions of the BMSPs forming through the evolutionary channel He N are multiplied by a factor of 5, and those of the BMSPs forming through the evolutionary channel CO N by a factor of 100. The orbital period distribution of the observed BMSPs listed in Table 1 is represented by the black solid line.

The orbital period distribution of the BMSPs forming through the He N channel (black) shows a peak between $`50`$ and $`100`$ days, depending on the envelope-ejection efficiency adopted for the common-envelope phase of the primary and on the velocity dispersion adopted for the magnitude of the kick velocity imparted to the neutron star at birth.

A larger envelope-ejection efficiency implies less spiral-in and thus a larger orbital separation at the end of the common-envelope phase. If the binary furthermore survives the subsequent supernova explosion of the primary, the circularised post-supernova orbital separation must be smaller than twice the orbital separation before the supernova explosion (e.g. Kalogera 1996). Increasing the common-envelope efficiency therefore shifts the peak in the orbital period distribution to longer orbital periods. If we identify the peak at long orbital periods in the observed orbital period distribution with the distribution resulting from the He N channel, the best agreement with the observed distribution is found for moderate to large values of the envelope-ejection efficiency ($`\alpha 1.05.0`$). An efficiency parameter larger than unity was also suggested by Tauris & Bailes (1996) in an investigation on the origin of the millisecond pulsar velocities and by Tauris & Sennels (2000) in a study on the formation of the binary pulsars PSR B2303+46 and PSR J1141-6545. Such a large value of $`\alpha `$ may be in conflict with the principle of energy conservation.

The velocity dispersion adopted for the magnitude of the kick velocity imparted to the neutron star at birth determines the spread of the orbital periods after the supernova explosion of the primary. Small velocity dispersions lead to a narrow distribution of orbital periods close to the upper limit beyond which the binary is disrupted by the supernova explosion, while larger velocity dispersions produce a larger spread with more systems at shorter orbital periods \[see Fig. 7 in Kalogera (1996)\]. This tendency becomes particularly clear if one compares the orbital period distribution in the absence of kicks (model K0) with the orbital period distribution in the case of large velocity dispersions (model KM400). The comparison with the observed long orbital period peak is most satisfactory in the absence of kicks or when the kick velocity dispersion is not excessively high ($`\sigma 100\mathrm{km}/\mathrm{s}`$).

Increasing the amount of mass that can be accreted by the neutron star from $`0.2M_{}`$ to $`0.6M_{}`$ does not significantly change the orbital period distribution of the BMSPs. In the limiting case of conservative mass transfer (not shown), the peak in the distribution is located around the orbital period of 15 days, which is far too short in comparison to the position of the long-period peak in the observed distribution. Systems with longer orbital periods do not appear in the conservative case because their progenitors have donor stars with masses below $`0.9M_{}`$, which evolve too slowly to fill their Roche lobe within the imposed age limit of 15 Gyr. Changes in the critical mass ratio for the delayed dynamical instability, finally, do not influence the orbital period distribution for systems evolving through this formation channel because the mass ratio of the progenitor systems at the onset of Roche-lobe overflow is always well below the critical mass ratio (see Fig. 6).

In all models considered, few systems are found with orbital periods longer than $`200`$ days. The decrease of the simulated distribution functions at longer orbital periods is, however, a lot less sharp than the decrease in the observed orbital period distribution, except for model K0. The paucity of systems beyond $`200`$ days also occurred in a test simulation where we allowed the neutron star to accrete at the Eddington rate without imposing an upper limit on the amount of mass that can be accreted. The lack of long period systems in the simulated orbital period distributions is therefore not due to the inefficiency of the accretion process resulting from, e.g., disc instabilities as proposed by Ritter & King (2002), but due to the upper limit on the post-supernova orbital periods which in turn is caused by the upper limit on the initial orbital periods (see Figs. 5 and 6 and the associated discussion in Sect. 4). The upper limit is preserved by the mass-transfer process after the supernova explosion due to the limited range of donor star masses ($`0.9M_{}M21.3M_{}`$) available for stable nuclear timescale mass transfer. This is consistent with the results presented by Podsiadlowski, Rappaport & Pfahl (2001).

The orbital period distribution of the BMSPs descending from the evolutionary channel He T (dark grey) typically has a peak at orbital periods smaller than 10 days, followed by a long-period tail which can extend up to $`60`$ days. The position of the peak is very robust with regard to changes in the common-envelope efficiency and the kick velocity dispersion, except in the case of the low envelope ejection efficiency $`\alpha =0.2`$ when the simulated distribution has a significant secondary peak located in the period gap of the observed orbital period distribution.

Variations in the assumptions regarding the mass-accretion rate onto the neutron star again do not significantly influence the distribution, except in the limiting case of conservative mass transfer (not shown). In the latter case, mass accretion from donor stars more massive than $`1.8M_{}`$ pushes the mass of the neutron star over the imposed upper limit of $`3.0M_{}`$, causing the neutron star to collapse into a black hole. Since mass transfer for systems evolving through this formation channel initially takes place on a thermal timescale, the orbital period of the system at the end of the Roche-lobe overflow phase increases with decreasing mass of the donor star at the onset of Roche-lobe overflow (e.g. King & Ritter 1999). Hence, the lack of systems with initial donor star masses above $`1.8M_{}`$ in model C shifts the peak in the orbital period distribution towards orbital periods that are too long in comparison to the position of the short orbital period peak in the observed distribution. A similar shift takes place when the critical mass ratio for the delayed dynamical instability is decreased from the values tabulated by Hjellming (1989) to $`q_{\mathrm{c},\mathrm{ddi}}=2.5`$ (model QCDD, not shown).

The BMSPs forming through the evolutionary channel CO T (light grey) generally dominate the simulated orbital period distribution at short orbital periods. The distribution function shows a peak around $`10`$ days, and a cut-off at $`100`$ days. In the conservative mass-transfer model (not shown), neutron stars with Roche-lobe filling companions more massive than $`1.8M_{}`$ again collapse into a black hole, so that no BMSPs are formed through this channel (see also Fig. 6). A decrease of the critical mass ratio for the delayed dynamical instability to $`q_{\mathrm{c},\mathrm{ddi}}=2.5`$ shifts the peak of the distribution function towards the period gap of the observed orbital period distribution. All other model parameters do not significantly affect the orbital period distribution of the BMSPs forming through this channel.

Comparison of the simulated orbital period distributions resulting from the formation channels He N and CO T shows that for comparable post-supernova orbital periods, case B Roche-lobe overflow from a low-mass donor star can lead to very different final orbital periods than case B Roche-lobe overflow from an intermediate-mass donor star. This behaviour was put forward by Taam et al. (2000) as a possible explanation for the period gap in the orbital period distribution of the observed Galactic BMSPs. We note that we here also find a contribution to the distribution function below the period gap from systems with intermediate-mass donors evolving through an initial thermal-timescale case A mass-transfer phase followed by a stable case A/case B Roche-lobe overflow phase.

Systems evolving into BMSPs through the formation channel CO N (medium grey) typically have orbital periods around $`1000`$ days and provide the smallest contribution to the total number of BMSPs in our simulated samples. The contribution of these systems to the total population of wide BMSPs is largest for kick velocities with a dispersion of $`\sigma 50`$ km/s. It is however clear that the orbital periods of these BMSPs are too long to contribute to the orbital period distribution of the observed Galactic BMSPs listed in Table 1. We therefore did not consider these systems in more detail for the present investigation.

### 5.3 Relative contributions

The relative contributions of the four formation channels to the total population of wide BMSPs are given in Table 3 for each of the models considered in our investigation. For comparison, the relative contributions to the observed orbital period distribution of systems below and above the period gap are about 75% and 25%, respectively.

The simulated populations of wide BMSPs are dominated by the systems forming through the evolutionary channel CO T, except when mass transfer is treated conservatively. However, due to the high mass-transfer rates inherent to donor stars located in the Hertzsprung gap, the question may be raised whether neutron stars in binaries evolving through this evolutionary channel can accrete enough matter to be spun up to a millisecond pulsar or not. If we impose a lower limit of $`0.05M_{}`$ on the mass that needs to be accreted by the neutron star (see, e.g., Burderi et al. 1999), only binaries with donor stars initially less massive than $`2.8M_{}`$ end up as BMSPs. The contributions of the four formation channels to the population of wide BMSPs then change as indicated in Table 4, so that the dominant contribution now stems from the systems evolving through the evolutionary channel He T. The contribution of the systems forming through the evolutionary channel CO T decreases even further if a lower limit of $`0.1M_{}`$ is imposed on the mass that needs to be accreted by the neutron star. In this case, the mass of the donor star at the onset of Roche-lobe overflow must be smaller than about $`2.4M_{}`$ for the binary to end up as a BMSP.

The relative contributions of the four formation channels also depend on the distribution adopted for the mass ratio of the binaries at the start of their evolution. Increasing the weight of binaries with smaller initial mass ratios results in a larger fraction of systems evolving through the evolutionary channel He N, so that more BMSPs are formed with longer orbital periods. Increasing the weight of binaries with larger initial mass ratios has the opposite effect. For illustration, the relative contributions in the case of the mass ratio distribution $`n(q)=1/q`$, for $`0<q1`$, are given in Table 5. In this case, the relative contributions of the systems evolving through the evolutionary channel He N approximately increase by a factor of 3. The overall behaviour of the orbital period distributions associated with the four formation channels is similar to that described in Sect. 5.2.

As a test, we also considered the case where the mass of the secondary is weighted independently from the mass of the primary according to the initial mass function given by Eq. (3). For the models considered, the corresponding orbital period distributions only show a peak at short orbital periods if no kicks are imparted to the neutron star at birth or if the velocity dispersion of the Maxwellian kicks is small ($`\sigma 50`$ km/s). However, as can be seen from Table 6, the contribution of the systems forming through the evolutionary channel CO N is then significantly too large in comparison to the observed orbital period distribution, so that an independent choice of the mass of both the primary and the secondary from the adopted initial mass function seems unlikely in this context.

### 5.4 A ’best-fit’ solution?

From the results presented in the previous subsections, it is clear that there are still too many uncertainties in some phases of stellar and binary evolution to come up with a unique set of input parameters giving the best possible representation of the observed orbital period distribution. An example of a possible combination of input parameters reproducing the main features of the observed distribution is shown in Fig. 9. The simulated distribution is obtained from combining model KM100 with a lower limit $`\mathrm{\Delta }M_{\mathrm{NS}}=0.1M_{}`$ for the mass that needs to be accreted to spin the neutron star up to millisecond periods and a mass ratio distribution $`n(q)=1/q`$, for $`0<q1`$. Since systems forming through the CO N channel only accrete about $`0.05M_{}`$, they here do not meet the criteria to be classified as BMSPs.

The combination of parameters reproduces the observed short- and long-period peaks as well as the observed period gap, but it fails to reproduce the BMSPs with orbital periods between 1 and 3 days. A possible reason is that the progenitors of these systems are likely to converge during the post-supernova mass-transfer phase, so that they evolve through a different evolutionary channel than those discussed here. The model furthermore only predicts BMSPs containing carbon/oxygen white dwarfs at orbital periods above the period gap, which may pose a problem to explain some of the observed systems with more massive white dwarfs and orbital periods shorter than 20 days (see Fig. 7). This potential inconsistency could be resolved if the uncertainties in the accretion process onto neutron stars during the rapid thermal-timescale mass-transfer phase of the systems evolving through the CO T channel is better understood.

The relative contributions of the systems forming through the He T, He N, and CO T channels are 72.51%, 21.65%, and 5.84%, respectively, which agrees fairly well with the relative number of systems below and above the period gap in the observed orbital period distribution. The agreement is even better if the observed distribution is restricted to orbital periods longer than 3 days, in which case 70% of the observed systems have orbital periods below the gap and 30% have orbital periods above the gap.

## 6 Concluding remarks

We have investigated the formation of wide BMSPs through four evolutionary channels using a rapid binary evolution code based on the analytical approximation formulae for the evolution of single stars derived by Hurley et al. (2000).

In three of the channels, dynamically unstable mass transfer from the primary prior to the supernova explosion results in a common-envelope phase during which the orbital separation is reduced and the helium core of the primary is exposed as a naked helium star. The primary continues its evolution until it becomes a neutron star in a type Ib/c supernova explosion. The further evolution depends on the mass of the secondary and on the post-supernova orbital period after circularisation:

* If the mass of the secondary is smaller than $`1.6M_{}`$, the binary undergoes a stable case B mass-transfer phase when the secondary reaches the giant branch. The binary becomes a BMSP consisting of a rapidly rotating neutron star orbiting a helium white dwarf with an orbital period determined by the mass of the white dwarf.
* If the mass of the secondary is larger than $`1.6M_{}`$ and the orbital period shorter than $`2.5`$ days, the binary undergoes a thermal-timescale case A mass-transfer phase, followed by a stable phase of case A and/or case B Roche-lobe overflow. The system again evolves into a BMSP containing a helium white dwarf whose mass determines the orbital period.
* If the mass of the secondary is larger than $`1.6M_{}`$ and the orbital period longer than $`2.5`$ days, the binary undergoes a thermal-timescale case B mass-transfer phase when the secondary crosses the Hertzsprung gap. The system evolves into a binary containing a neutron star and a carbon/oxygen white dwarf without developing a tight correlation between the orbital period and the mass of the white dwarf. The accretion onto the neutron star in these binaries may not be efficient enough to spin them up to millisecond pulsars.

The fourth evolutionary channel corresponds to the direct supernova mechanism where no interaction between the binary components takes place until the primary gives birth to a neutron star in a type II supernova explosion on the AGB. If the resulting post-supernova orbit is wide enough to avoid any interaction until the secondary becomes an AGB star, stable case C mass transfer leads to the formation of a BMSP containing a carbon/oxygen white dwarf whose mass is correlated with the orbital period.

Our main goal in this investigation was to compare the orbital period distribution of simulated samples of BMSPs resulting from different sets of input parameters for the binary evolution calculations with the orbital period distribution of wide BMSPs observed in the Galactic disc. The observed long orbital period peak is identified with the BMSPs containing a helium white dwarf whose progenitor underwent a mass-transfer phase that was stable at all times. The observed short orbital period peak on the other hand is identified with the BMSPs containing a helium white dwarf whose progenitor evolved through a thermal-timescale mass-transfer phase on the main sequence or in the Hertzsprung gap.

The simulated distribution functions all show a rapid decrease for orbital periods longer than 200 days, irrespective of the accretion efficiency of neutron stars. The lack of longer period systems is a consequence of the upper limit on the initial orbital periods beyond which the binary remains detached instead of going through a common-envelope phase prior to the supernova explosion of the primary. The agreement between the simulated and the observed orbital period distributions is best for models with highly non-conservative mass transfer, common-envelope efficiencies equal to or larger than unity, a critical mass ratio for the delayed dynamical instability larger than 3, and no or moderate supernova kicks at the birth of the neutron star.

The first results of the population synthesis study presented here are encouraging, but many uncertainties remain to be solved. The lack of BMSPs containing carbon/oxygen white dwarfs at orbital periods shorter than 20 days in our ’best-fit’ solution, for instance, surely poses a problem that needs to be resolved. A more detailed treatment of the accretion process onto neutron stars and the potential effects of pulsar turn-on, evaporation, and spin-down are of particular importance to improve the population synthesis study of wide BMSPs. Preliminary calculations also reveal a dependency of the orbital period distribution on the mass-loss rates from stellar winds, but the results are still inconclusive.

Furthermore, our population synthesis study implicitly incorporates a model for the Galactic population of neutron star X-ray binaries. The long-standing apparent conflict between the observationally implied birth rate of BMSPs and their LMXB progenitors (e.g. Ruderman, Shaham & Tavani 1989; Podsiadlowski et al. 2002 and references therein) may point towards non-standard evolutionary effects in the X-ray binary phase which are not considered in our study. It is therefore desirable to consider the populations of LMXBs and BMSPs simultaneously.

We will address these matters in more detail in future investigations.

## Acknowledgements

We are grateful to Hans Ritter for providing the data of the observed BMSPs and to Jarrod Hurley, Onno Pols, and Chris Tout for sharing their SSE software package and a copy of the paper describing their binary evolution algorithm prior to publication. We also thank Jarrod Hurley, Chris Tout, Simon Portegies-Zwart, and Firoza Sutaria for useful discussions. We thank the referee, Philipp Podsiadlowski, for useful comments. This research was supported by the British Particle Physics and Astronomy Research Council (PPARC).

## Appendix A Accretion onto neutron stars

In this appendix, we briefly outline the treatment of mass accretion by neutron stars in semi-detached binaries adopted in our binary evolution code. We assume mass transfer to be non-conservative and denote by $`\gamma `$ the fraction of the transferred mass that is lost from the system. The mass-accretion rate onto the neutron star is then related to the mass-transfer rate from the donor star by

$$\dot{M}_{\mathrm{NS}}=(1\gamma )|\dot{M}_2|,$$
(6)
where $`M_{\mathrm{NS}}`$ is the mass of the neutron star, and $`M_2`$ the mass of the Roche-lobe overflowing companion.

If we furthermore assume that the neutron star accretes a maximum mass $`\left(\mathrm{\Delta }M_{\mathrm{NS}}\right)_{\mathrm{max}}`$ during the mass transfer process, the average mass-accretion rate onto the neutron star is given by

$$\dot{M}_{\mathrm{NS}}=\left(\mathrm{\Delta }M_{\mathrm{NS}}\right)_{\mathrm{max}}\frac{|\dot{M}_2|}{M_2},$$
(7)
so that

$$\gamma =1\frac{\left(\mathrm{\Delta }M_{\mathrm{NS}}\right)_{\mathrm{max}}}{M_2}.$$
(8)
In order to ensure that the neutron star accretes not more than $`\left(\mathrm{\Delta }M_{\mathrm{NS}}\right)_{\mathrm{max}}`$, we then adopt the artificial rate given by (7) as the mass-accretion rate onto the neutron star.


text:
# Particle density fluctuations

## 1 Introduction

At high temperatures or high baryon number density, Quantum Chromo-Dynamics (QCD) describes a world of weakly interacting quarks and gluons very different from the hadronic world in which we live. This suggests the possibility of a phase transition as the temperature or density is increased in which there is a transition from a state of matter where quarks are confined inside hadrons to one where quarks are free to move around (deconfined ) within a large volume - the quark gluon plasma (QGP). This can be addressed through experimental studies involving relativistic heavy-ion collisions. Experimental searches have focused on isolating signatures of two types of phase transitions which might occur in extremely hot and/or dense nuclear matter. One is related to the deconfinement of quarks while the other is related to chiral symmetry restoration. There is also an interesting possibility of the existence of a tri-critical point in the phase diagram , where the transition changes from first to second order. At the tri-critical point, one would observe singularities in several thermodynamical variables, such as the specific heat and matter compressibility. The thermodynamical variables are related to event-by-event fluctuations in experimental observables like particle multiplicity, transverse energy, and mean transverse momentum. For example, the total heat capacity is related to transverse momentum fluctuations and matter compressibility to multiplicity fluctuations . Hence the existence of the tri-critical point can be probed in an experiment by varying the control parameters of the reaction (impact parameter, rapidity acceptance, beam energy, and system size) and studying the event-by-event fluctuations in the global observables. It is also believed that in high energy heavy-ion collisions there is a possibility of creating a chiral symmetry restored phase. One of the possible interesting consequences of chiral symmetry restoration, is the formation of disoriented chiral condensates (DCC) . The detection and study of a DCC state is expected to provide valuable information about the chiral phase transition and vacuum structure of strong interactions. The probability distribution of the neutral pion fraction in a DCC domain follows the relation :

$$P(f)=1/2\sqrt{f}\mathrm{where}f=N_{\pi ^0}/N_\pi .$$
(1)
Since the majority of charged particles consist of charged pions and majority of photons originate from $`\pi ^0`$ decays, DCC formation in a given domain would be associated with large correlated event-by-event fluctuations in the multiplicities of charged particles and photons. Here we present the experimental results of multiplicity fluctuations and correlations in relativistic heavy-ion (Pb+Pb) collisions as measured by the WA98 experiment  at the CERN SPS.

## 2 Multiplicity fluctuations

Recently, much theoretical interest has been generated on the subject of event-by-event fluctuations, primarily motivated by the nearly perfect Gaussian distributions of several observables ($`N_\gamma `$, $`N_{ch}`$, $`E_T`$, and $`<p_T>`$) for a fixed centrality bin . One may define the relative fluctuation as,

$$\omega _X=\frac{\sigma _X^2}{<X>}$$
(2)
where $`\sigma _X^2`$ is the variance of the distribution. The fluctuation is then studied by varying the control parameters of the reaction, such as the centrality and rapidity acceptance.

The WA98 experiment at the CERN SPS has carried out a detailed analysis of photon and charged particle multiplicity  fluctuations. The photon multiplicity is measured using the photon multiplicity detector (PMD)  which has the pseudo-rapidity coverage from 2.9 to 4.2. The charged particle multiplicity is obtained from the silicon pad multiplicity detector (SPMD)  with pseudo-rapidity coverage from 2.35 to 3.75. The centrality of the reaction is defined through the measurement of transverse energy ($`E_T`$) with the mid-rapidity calorimeter (MIRAC) having a pseudo-rapidity coverage from 3.5 to 5.5.

### 2.1 Centrality dependence of multiplicity fluctuations

It is very important to control the centrality selection carefully in fluctuation studies so that the impact parameter fluctuations are kept to a minimum and the distributions are good Gaussians. With this in mind we have used narrow centrality selections with bin widths of 2% in cross section, as determined from the measured total transverse energy, such as, 0-2%, 2-4%, 4-6%, 6-8%, etc. The fluctuations calculated for these centrality bins using Eqn. 1 are plotted in Fig. 1 as a function of the number of participants. The results from data for both charged particles and photons are compared to those obtained from the VENUS event generator and a simple participant model . In the participant model, where the total particle multiplicity in an event is the sum of the number of particles produced by each participant, the multiplicity fluctuations, $`\omega _N`$, can be expressed as

$$\omega _N=\omega _n+n\omega _{N_{\mathrm{p}art}}$$
(3)
The value of $`\omega _n`$, the fluctuation in the number of particles falling within the detector acceptance ($`n`$) produced per participant can be obtained from nucleon-nucleon data. The value of $`\omega _{N_{\mathrm{p}art}}`$, the fluctuation in the number of participants, can be obtained from simulations. From Fig. 1 we observe that the measured fluctuations are reasonably well reproduced by statistical models.

### 2.2 Acceptance dependence of multiplicity fluctuations

Fluctuations in the particle multiplicity have also been studied by varying the rapidity acceptance. The results are shown in Fig. 2. One observes, contrary to naive expectation that the fluctuations should increase as the multiplicity decreases with decrease in acceptance, that instead the fluctuations decrease. This can be explained through a simple statistical model, based on the assumption that particles are accepted following a binomial distribution. The details of this model are given in Ref . As per this model, the fluctuations ($`\omega _n`$ ) in a small acceptance region are related to the fluctuations ($`\omega _m`$) in a larger acceptance as

$$\omega _n=1f+f\omega _m$$
(4)
where $`f`$ is the ratio of the average number of particles in the smaller acceptance to the average number of particles in the larger acceptance. The results from the model shown in Fig. 2 are found to be in reasonably good agreement with the data.

## 3 $`N_\gamma `$ vs. $`N_{ch}`$ fluctuations

The main motivation to study event-by-event fluctuations in the photon vs. charged particle multiplicity is to search for possible formation of DCCs. For this one would like a photon multiplicity detector and a charged particle multiplicity detector with as large a common coverage in $`\eta \varphi `$ phase space as possible. The main DCC search result from the SPS has come from the WA98 experiment . First we discuss the results for top $`5\%`$ most central events. For the pseudo-rapidity and azimuthal region common to both PMD and SPMD used for this analysis, the average photon and charged particle multiplicity is 335 and 323, respectively. The experimental results are compared to simulated events and various types of mixed events to investigate possible DCC formation. The simulated events were generated by passing the VENUS output through a detector simulation package (GEANT), which incorporates the WA98 experimental setup. Four different types of mixed events were generated from the data with equivalent event sample sizes as for real events. The details of the construction of the mixed events can be found in Ref. . In Table 1 we summarize the construction and the physics issue probed by each type of mixed event.

A common analysis was carried out on the data, simulated, and mixed events to investigate the source of any observed fluctuations.

### 3.1 $`N_\gamma `$ vs. $`N_{ch}`$ correlation analysis

Analysis of the correlation between $`N_\gamma `$ and $`N_{\mathrm{ch}}`$ is useful to search for DCC-type fluctuations. The details of this can be found in Ref. . From the event-by-event correlation between $`N_{\gamma \mathrm{like}}`$ and $`N_{\mathrm{ch}}`$ in various $`\varphi `$-segments (obtained by dividing the $`\varphi `$-space into 2, 4, 8, and 16 bins) and a common correlation axis ($`Z`$), one can obtain the closest distance ($`D_Z`$) of the data points to the correlation axis. The correlation plot is shown in Fig. 3.

In order to compare the fluctuations for different $`\varphi `$ bins on a similar footing, a scaled variable, $`S_Z=D_Z/s(D_Z)`$, is used, where $`s(D_Z)`$ represents the RMS deviation of the $`D_Z`$ distribution for VENUS events analyzed in the same manner. The width of the distribution of $`S_Z`$ represents the relative fluctuations of $`N_{\gamma \mathrm{like}}`$ and $`N_{\mathrm{ch}}`$ from the correlation axis at any given $`\varphi `$ bin. Since the width of the $`S_Z`$ distribution quantifies the amount of fluctuation, the RMS deviations of these distributions for data are compared with those from mixed events and simulations in Fig. 4.

The RMS deviations of the $`S_Z`$ distributions for mixed events and data agree for 1 bin in $`\varphi `$ (by construction) and for 16 bins in $`\varphi `$ within the quoted errors. However, for the other bins in $`\varphi `$ one observes that results from M1 mixed events are lower than those from data. This indicates the presence of localised non-statistical fluctuations in the data. The source of the additional fluctuations is understood by the comparison to the M3-type of mixed events. Comparison shows that the excess fluctuations have contributions both from $`N_\gamma `$ and $`N_{\mathrm{c}h}`$. However the RMS deviations of the M2-type of mixed events closely follow those from data, suggesting the absence of event-by-event correlated DCC-like $`N_\gamma `$ vs. $`N_{\mathrm{c}h}`$ fluctuations.

### 3.2 Multi-resolution analysis based on discrete wavelet transformation

A multi-resolution analysis using discrete wavelet transformations (DWT) has been shown to be quite powerful in the search for localized domains of DCC . For the present DWT analysis the full azimuthal region is divided into smaller bins in $`\varphi `$, the number of bins at a given scale $`j`$ being $`2^j`$. The input to the analysis is a spectrum of the sample function at the smallest bin in $`\varphi `$ corresponding to the highest resolution scale, $`j_{max}`$(= 5 here). In the present case, the sample function is chosen to be the photon fraction, given as,

$$f^{}(\varphi )=N_{\gamma \mathrm{like}}(\varphi )/(N_{\gamma \mathrm{like}}(\varphi )+N_{\mathrm{ch}}(\varphi ))$$
(5)
The output of the DWT consists of a set of wavelet or father function coefficients (FFC) at each scale, from $`j=1`$,…,($`j_{max}1`$). The coefficients obtained at a given scale, $`j`$, are derived from the distribution of the sample function at one higher scale, $`j+1`$. The FFCs quantify the deviation of the bin-to-bin fluctuations in the sample function at that higher scale relative to the average behavior. The presence of localized non-statistical fluctuations will increase the RMS deviation of the distribution of FFCs and may result in non-Gaussian tails . Once again, comparing the RMS deviations of the FFC distributions of data, mixed events, and VENUS events may allow to draw inference about the presence of localized fluctuations.

The RMS deviations of these FFC distributions are summarized in Fig. 5. The RMS deviations of the FFC distributions for the data, VENUS, and mixed events are found to be close to each other (within quoted errors) for the case of $`32`$ bins in $`\varphi `$. While the values for M2 mixed events are found to closely follow those of the data for all bins in $`\varphi `$, the RMS deviations for the M3 mixed events lie between those of the data and M1 mixed events. These results are consistent with those obtained from the analysis of the $`S_Z`$ distributions. These observations indicate the absence of event-by-event localized correlated fluctuations (DCC-like) between $`N_{\gamma \mathrm{like}}`$ and $`N_{\mathrm{ch}}`$. However, they do indicate the presence of localized independent fluctuations in both photon and charged particle multiplicities for intermediate bin sizes in azimuth.

### 3.3 Centrality dependence of $`N_\gamma `$ vs. $`N_{ch}`$ fluctuations

It is interesting to study the centrality dependence of $`N_\gamma `$ vs. $`N_{ch}`$ fluctuations . This study has been carried out for four centrality classes corresponding to the top $`5\%`$, $`5\%`$ \- $`10\%`$, $`15\%`$ \- $`30\%`$, and $`45\%55\%`$ of the minimum bias cross section. The correlation and DWT analysis was carried out on events from each of the centrality classes and the results in terms of RMS deviations of $`S_Z`$ and FFC distributions for data, various sets of mixed events, and simulations were compared. For all four classes, the RMS deviations of the $`S_Z`$ and FFC distributions of data and M2 mixed events agreed reasonably well with each other. However, comparison to M1 and M3 mixed events showed the presence of localised uncorrelated event-by-event non-statistical fluctuations in both photon and charged particle multiplicities. In order to quantify the strength of the total localised $`N_{\gamma \mathrm{like}}`$ and $`N_{\mathrm{ch}}`$ fluctuations for various bins in $`\varphi `$ and for different centrality classes, we define a sensitivity parameter $`\chi `$ as :

$$\chi =\frac{\sqrt{(s^2s_1^2)}}{s_1}$$
(6)
where $`s_1`$ and $`s`$ correspond to the RMS deviations of the FFC distributions of the M1 mixed events and real data, respectively. The results are shown in Fig. 6 as a function of the number of bins in $`\varphi `$ for the four different centrality classes. Qualitatively similar results are obtained when $`\chi `$ is calculated using the RMS deviations of the $`S_Z`$ distributions. The shaded portion indicates the region of $`\chi `$ where $`s`$ is one $`\sigma `$ greater than the RMS deviation FFC distributions for M1 events, where $`\sigma `$ is the total error on the M1 event RMS deviation. It represents the limit above which a signal is detectable. The result shows that the strength of the fluctuations decreases as the number of bins in $`\varphi `$ increases, with a strength which decreases to below detectable level (within the quoted errors) for 16 and 32 bins. It is also observed that the strength of the signal decreases with decreasing centrality for 4 and 8 bins in azimuthal angle, although the tendency is not very strong.

### 3.4 Upper Limit on DCC production at SPS

If the amount of DCC-like fluctuations in the experimental data were large, then the RMS deviations of $`S_Z`$ and FFC distributions for data would have been larger than those of M2 events. Since this is not the case, we have extracted upper limits on the probability of DCC-like fluctuations at the 90% confidence level following the standard procedure as discussed in Ref. . This is done within the context of the simple DCC model described in Ref. . We give the upper limits for the two most central event classes ($`05\%`$ and $`510\%`$). The $`90\%`$ CL upper limit contour has been calculated as $`\chi `$ \+ 1.28$`e_\chi `$, where $`\chi `$ is calculated using Eqn.(6). Here $`s_1`$ and $`s`$ correspond to the RMS deviations of the FFC (or $`S_Z`$) distributions for M2 mixed events and real data, respectively, and $`e_\chi `$ is the error in $`\chi `$ from the FFC (or $`S_Z`$) analysis. To relate the measured upper limit on the size of the fluctuations in terms of DCC domain size and frequency of occurrence we take use the simple simulated DCC model described in Ref. . The upper limit is set at that value of frequency of occurrence for a fixed DCC domain size at which the $`\chi `$ value from the DCC model matches with that of the $`\chi `$ \+ 1.28$`e_\chi `$ upper limit from the experimental data. The results for centrality classes 1 and 2 are shown in Fig. 7.

## 4 Summary

A detailed study of centrality and rapidity acceptance dependence of multiplicity fluctuations carried out by the WA98 experiment at the CERN SPS shows an absence of significant non-statistical fluctuations in the photon and charged particle multiplicity. A model independent study of event-by-event correlated fluctuations (DCC-type) in photon and charged particle multiplicity using a robust mixed event technique reveals an absence of significant DCC-like correlated fluctuations at the SPS. However, the analysis indicates the presence of event-by-event uncorrelated fluctuations in both $`N_\gamma `$ and $`N_{\mathrm{ch}}`$ beyond those observed in simulated and mixed events for limited region of azimuthal angle, with the strength of fluctuation increasing with increase in centrality. Using the results from the data, mixed events, and a simple model of DCC formation, an upper limit on DCC production in Pb+Pb collisions at SPS energies has been set.

#### Acknowledgments

One of us (B.M) is grateful to the Board of Research on Nuclear Science and Department of Atomic Energy, Government of India for financial support in form of the Dr. K.S. Krishnan fellowship. This work was supported jointly by the German BMBF and DFG, the U.S. DOE, the Swedish NFR and FRN, the Dutch Stichting FOM, the Polish KBN the Grant Agency of the Czech Republic the DAE, DST, CSIR, UGC of the Government of India, the PPE division of CERN, the Swiss National Fund, The Ministry of Education, Science and Culture, the University of Tsukuba Special Research Projects, and the JSPS Research Fellowships for Young Scientists.