S Touimi

J L Jauberteau

) I Jauberteau

J Aubreton

Plasma chemistry and diagnostic in an Ar-N 2 -H 2 microwave expanding plasma used for nitriding treatments

This paper reports on mass spectrometry analysis performed downstream a microwave discharge in Ar-N 2 -H 2 gas mixture in nitriding conditions. Investigations are focused on the main simple radicals NH 2 , NH and N, and on the molecular species NH 3 and N 2 H 2 produced.

Because of wall desorptions due to catalytic effects, we must develop a specific method taking into account both wall desorption and dissociative ionization effects, in order to correct the mass spectrometer signal intensity. The relative concentrations of the previous species are studied in various gas mixtures. Correlations are made between the plasma chemistry and plasma parameters (electron density and energy electron distribution function), measured by means of Langmuir probes spatially resolved within the plasma expansion. These results show the efficiency of ternary gas mixtures (Ar-N 2 -H 2 ) to produce electrons and N x H y species used in plasma nitriding process.

introduction.

Microwave discharge sustained in Ar-N 2 -H 2 gas mixtures are used in a wide range of industrial applications. Like ammonia synthesis /1/, nitriding processes for surface treatment of metallic or polymer pieces /2,3,4/ and for the elaboration of optical or electronic devices /5/. The efficiency of such plasma processes depends on the knowledge of the various elementary mechanisms involved in the gas phase and at the surface of the reactor wall or of the work pieces. In previous papers /6, 7/, we have shown that the surface oxide reduction of thin molybdenum films, in the nitriding treatment, depends on the gas composition in the plasma expansion and on the distance between the work piece surface and the discharge centre. We confirm correlations between nitrided layers characteristics, discharge expansion

conditions, gas composition and Langmuir probe measurements performed within the plasma expansion.

The purpose of this paper is to investigate the various chemical species produced within the discharge and in the expansion in nitriding treatment conditions. We make correlation between the gas compositions and the plasma parameters in the working conditions.

Investigations are performed by means of mass spectrometer downstream a microwave discharge and by means of spatially resolved Langmuir probe in the plasma expansion in nitriding conditions 2. Experimental set-up.

The experimental setup is shown on Figure 1. It consists of a microwave discharge (SAIREM generator working at 2.45 GHz) produced in a quartz tube (internal diameter 16mm and external diameter 19mm). In order to keep a constant wall temperature, the quartz tube is refreshed by means of air cooling. A Roots pumps (70-700 m 3 /hour) with an adjustable rotational speed is used and a gas flow velocity ranges from 10 to 30 m/s is maintained constant in the stainless tube (inner diameter 50mm) located above the discharge and collecting species produced within the plasma. The flow is considered as a "plug" flow and the velocity has been previously measured by means of absorption spectroscopy measurements on Ar( 3 P 2 ) metastable species produced in an Ar microwave discharge. The total pressure in the discharge tube is held constant and ranges from 10 to 100 Pa. The mass spectrometer (QMG 421 Balzers) is fixed above this stainless steel tube at 30cm from the discharge. The ionization of neutral species is performed using the electron impact technique by means of a cross beam source. Before each experiment, the reactor is heated and pumped at 10 -4 Pa and maintained at this pressure for several hours by means of a molecular pump in order to clean the reactor wall. Because the quadrupole cannot operate at pressure lower than 10 -3 Pa, a specific two-stage differential pumping device has been made and is detailed in ref /8/. It consists of two independent coaxial chambers fitted into each other and equipped with separated Turbo pumps units. The quadrupole is located at the middle of the inner chamber (see Figure1). A charged particles collector (Langmuir probe) is located in the grounded stainless steel tube. Results show that no ion or electron is crossing this area located 20cm below the mass spectrometer sample hole. Thus, investigations are performed in a region without charged particle and no reaction between charged particles and neutral species occurs in this tube between the discharge and the mass spectrometer. The present study gives information of reactivity between neutral species only. Even though, in the microwave expanding plasma used for molybdenum nitriding process, reactive processes between charged particles and neutral species are efficient probably until the substrate holder /7/.

3. Mass spectrometry analysis downstream the discharge.

In discharges sustained in the ternary gas mixture Ar-N 2 -H 2 , NH x species partly results of plasma catalysis processes on the reactor wall /9, 10/ or are directly produced in the discharge bulk. A large part of the mass spectrometer signal intensity is due to desorbed species from the reactor wall. This effect occurs when the discharge is off or on. So, it is necessary to take into account this trouble in order to measure the real effect of the discharge on the gas composition.

The main species detected

The main species detected in the Ar-N 2 -H 2 discharge afterglow are observed at m/q=1, 2, 14, For all species (i) detected, the signal intensity measured downstream the discharge, where ion concentration is negligible, is the sum of two contributions. One is due to the direct ionization of species (i) and the second is due to the dissociative ionization of larger species (j). This effect can be written /8/,

I i = I(i + /i) + Σ j I(i + /j)=T i .σ (ι+/i) (ε e ).n i + T i . Σ j (σ (ι+/j) (ε e ).n j ) (1)
Where, σ (ι+/i) (ε e ) is the ionization cross section due to the direct ionization of i species and σ (i+/j) (ε e ) is the dissociative ionization cross section of j species producing i ions. T i is the transmission factor. It depends on m/q, on the pressure and on the electron energy in the ionization chamber. When the measurements are performed within the discharge, we must add a term to the signal given by Equation 1, which is due to ions (i) produced within the discharge.

The main problem is now to separate in the total signal intensity the part due to the direct ionization of i species from the part due to the dissociative ionization processes of j species larger than i species. In the next part, we will investigate these different contributions in the case of the different species detected.

3.1.1. The case of NH 3 .

This species is detected at m/q=17 and consequently, the signal measured can be mixed with that of OH + ion, which is due to the dissociative ionization process of H 2 O (residual water) in the ionization chamber. The contribution of each dissociation process can be estimated, considering the change of the signal intensity versus the electron energy used in the ionization chamber. This method was previously reported by Where T( m/q=18 ) is the transmission factor for m/q=18, (H 2 O + /H 2 O) is the ionization cross section for the direct ionization of the water and n H2O is the residual water concentration.

Mixing these two last equations, the contribution to the signal intensity due to the direct ionization of ammonia is given by,

I(NH 3 + /NH 3 )= I (m/q=17) -[T (m/q=17) / T(m/q=18) (OH + /H 2 O)/ (H 2 O + /H 2 O) I (m/q=18) )] (2)
In this equation, the ratio T (m/q=17) /T (m/q=18) is calculated using results reported in ref (8).

Consequently, the part due to the direct ionization of NH 3 , can be calculated by means of the signal intensity measured at m/q=17 and m/q=18 and using Equation 2. Results are reported on Figure 2. It can be seen a good agreement between results of Equation 2 and the previous ones calculated considering the direct ionization cross section values given in the literature in the case of deuterated species (Tarnovsky or Rejoub /17-18/). The part due to the direct ionization of NH 3 ranges from 70 to 90% of the total signal intensity measured for m/q=17, when electron energy ranges from 25 to 70 eV. The part of the signal results mainly of the dissociative ionization of water and can be easily suppressed considering the difference between the signal measured when the discharge is on and off. This simple way for the signal correction will be used in the next part of the article in order to suppress the contribution due to the wall desorption in the afterglow.

3.1.2. The case of NH x radicals.

The same considerations have been performed in the case of the simple radicals NH 2 + , NH + and N + detected in the discharge afterglow. quite different. In one case (Tarnovsky), we can consider that there is no contribution due to the dissociative ionization of ammonia. In the second case (Märk), this contribution would be quite obvious. Assuming that this last case gives an overestimation of the dissociative ionization effect, we will consider it in the following part of this paper, to correct the signal of NH 2 . In these conditions (using Märk et al cross section values), results show that the part of the direct ionization of NH 2 is larger than 70% of the total signal, when electron energy ranges from 25 to 70eV.

As previously done for NH part remains lower than 30%, in our experimental conditions. In the case of NH and N, this correction can be neglected. In the case of NH 3 , the dissociative ionization is mainly due to desorbed water from the wall, which can be easily withdrawn considering the difference between the signal intensity measured with and without discharge. In case of NH 2 , the part due to the dissociative ionization is mainly due to NH 3 and the correction factor depends on the electron energy and is displays on Figure7. In the case of N 2 H 2 , correction is not possible because ionization cross section values are not available.

In a previous work /10/, we report a more accurate method taking into account the different dissociative processes contributing to the signal intensity measured for each species.

Polynomial equations are calculated using the relative concentrations values, NH 2 /NH 3 , NH/NH 3 and N/NH 3 . The coefficients of these polynomial equations depend on the cross sections values of the different ionization processes involved. This should be a more accurate method that the previous one. However, because of the strong uncertainty on the cross section values reported in the literature, 15% for the parent ionization cross section values and 18%

for the dissociative ionization cross section values/17/, the different coefficient of the polynomial equations are determined with a too large uncertainty (typically from 30% to 60%, according to the coefficient). The error on relative density is very large and can be over 100% according to the experimental conditions, especially in gas composition with high N 2 % and H 2 %, producing a large amount of NH 3 . For these reasons in the present works corresponding to experiments performed in nitriding conditions, this method cannot be used because of the too large error and the previous one is preferable.

Mass spectrometer diagnostics in nitriding conditions.

In previous works, it was shown that a (Ar-25%N 2 -30%H 2 ) plasma exposure at 8cm from the centre of the discharge seems the most efficient ternary gas mixture to reduce the oxide layers remaining at the surface of a molybdenum film /7/. Correlations had already been performed

with electron density and plasma propagation conditions. As previously related this gas composition lead to a better plasma expansion due to a larger electron density than in pure nitrogen or in other ternary gas mixture with lower H 2 contents /7,8/. The present work gives complementary results about the chemical gas composition in such ternary gas mixture used in expanding microwave plasma.

Two gas mixtures have been investigated to carry out this study. They are:

The gas mixture 1:

It contains (0 to 30%) H 2 injected in the binary mixture (64%Ar-36%N 2 ). It corresponds to a gas flow containing 180 sccm Ar, 100 sccm N 2 .

The gas mixture 2:

It contains (0 to 15%) H 2 injected in the binary mixture (91%Ar-9%N 2 ). It corresponds to a mixture with 1000 sccm Ar and 100 sccm N 2.

In the different series N 2 gas flow is equal to 100 sccm and the total pressure Ar+N 2 is kept constant and is equal to 100 Pa. The pressure is controlled by changing the rotational speed of the Roots blower pumps. The second gas mixture contains more Argon than the first one, and recombination processes of NH x radicals are expected to be less important than in the first one because of the dilution. For the two gas mixture, the microwave power used for the chemical analysis is 150 W. considering the two previous gas mixtures respectively. In these figures the main radical produced is NH 2 then NH and N. The ratio NH x /NH 3 slowly decreases with increasing % H 2 .

In the first gas mixture, the ratio NH 2 /NH 3 ranges from 0.1 to 0.5 all over the %H 2 under investigation. The ratio NH/NH 3 and N/NH 3 ranges roughly from 0.1 to 0.01 and lower than 0.01 for the first ratio and the second respectively. In the gas mixture2, the relative NH 2 concentration produced is larger than in the gas mixture 1, and is over 1. The two other ratios N/NH 3 and NH/NH 3 are almost identical to those obtained in the former gas mixture. This increase in the relative NH 2 concentration can be explained by the larger Ar flow rate used in the mixture 2. Radicals produced in such conditions are more diluted and the recombination process is decreasing. The effect of the dilution is more evident on NH 2 that on the two other radicals which concentrations are one or two magnitude orders lower.

Figure 10, shows the relative concentration of ammonia compared to the argon (NH 3 /Ar), versus the %H 2 mixed to Ar+N 2 , in the two gas mixtures (64%Ar-36%N 2 ) and (91%Ar-9%N 2 ). It can be seen that for the two binary gas mixtures under investigation, the relative concentration of NH 3 remains nearly the same and roughly linearly depends on H 2 % injected.

No saturation is observed on the figure. Comparison with results displayed on Figure 8 and9 prove that the concentration of radicals NH 2 , NH and N also increase with % H 2 increasing in both Ar-N 2 binary mixtures. On Figure 10 are also reported three ternary gas mixtures which will be discussed later in the text. The admixture of H 2 in the gas composition increases the coadsorbed hydrogen density (on a Ni(100) surface), and stabilizes the NH intermediate, increasing the surface concentration of imide in the hydrogen-deficient surface environment in the 200K-450K /24/. Consequently the diimide concentration is significantly increased in the presence of post-adsorbed hydrogen even in a 10 -7 Torr hydrogen flow /24/. For larger %H 2 injected in the gas mixture, the surface sites are saturated with hydrogen and the recombination surface mechanism of NH s producing N 2 H 2 is now decreasing. Thus, there is a maximum of efficiency to produce N 2 H 2 by means of the previous surface reaction. In our results, this maximum is observed at about 5% H 2 injected.

Correlation between plasma chemistry and plasma parameters.

In order to make correlations between the plasma parameters in the expanding plasma and the gas composition, investigations have been performed by means of Langmuir probe downstream the discharge in the reactor previously used for plasma nitriding treatment /7, 8/.

The Langmuir probe is fixed in the expanding plasma along the expansion axis and can be moved according to this axis from the tube discharge exit to the substrate holder. This experimental set-up is detailed in ref 7 and 8. The electron energy distribution function has been measured using the Druyvesteyn equation and the numerical simulation of harmonic component method previously reported /27/. Measurements have been performed at injected powers equal to 400 W and 200 W and in three gas mixtures, corresponding to Ar-30%N 2 -

12%H 2 , Ar-8%N 2 -10%H 2 , Ar-25%N 2 -30%H 2 .
As previously reported /7/, a ternary mixture containing Ar-N 2 -H 2 can be used to produce plasma expansion below the discharge tube exit in the reactor even with a large amount of H 2 and N 2 (more than 50%) mixed in Ar. Whereas in a binary mixture containing Ar-N 2 or Ar-H 2 , the plasma expansion is strongly shrinking in the discharge tube for H 2 or N 2 concentration larger than (1% for H 2 and 5% for N 2 ). Figure 12 shows the electron density measured at 400W, versus the distance in the plasma expansion.

Measurements have been performed in the three previous mixtures, Ar-30%N 2 -12%H 2 , Ar-8%N 2 -10%H 2 and Ar-25%N 2 -30%H 2 , in the binary mixtures, Ar-1%N 2 , Ar-50%N 2 and in pure Ar. The electron density is calculated using the electron current value measured at the plasma potential /27/. Binary mixtures containing Ar-H 2 are strongly disturbed even at a low %H 2 injected (<1%) and accurate results are not available. As expected the electron density decreases with increasing distance from the tube discharge exit. It can be seen that when N 2 is injected in argon, the electron density decreases with the % N 2 increasing, the same effect could be observed with H 2 . This behaviour can be ascribed to the increase of inelastic collisions in the gas mixture. Surprisingly, in the ternary gas mixtures Ar-25%N 2 -30%H 2 and Ar-8%N 2 -10%H 2 , the electron density increases compared to results obtained in pure argon discharge even for a (N 2 +H 2 ) admixture larger than 50%. Nevertheless in the case of Ar-30%N 2 -12%H 2 , the electron density remains very low, lower than in Ar-50%N 2 .

Figure 13 displays the electron distribution function measured at 200 W, 3cm below the discharge exit in the three ternary gas mixtures under investigation (Ar-30%N 2 -12%H 2 , Ar-8%N 2 -10%H 2 and Ar-25%N 2 -30%H 2 ). It can be seen that for Ar-30%N 2 -12%H 2 , the electron distribution function is smaller than those obtained for the two other gas mixtures and in Ar-8%N 2 -10%H 2 , the electron distribution function is the largest. As expected, these results agree with the previous ones displayed on Figure 12, the electron density can also be calculated using the electron distribution functions /27/. These results prove the good efficiency of the second and third gas mixtures for the plasma expansion even at 200W, and the difficulty encountered with the first one to sustain the plasma expansion at a power lower than 400W.

Comparison with previous results obtained by means of mass spectrometer in the same gas mixtures (see the Figure 8 drastically increases the electron density, since the total ionization cross section of N 2 is larger than that of H 2 and the dissociation of H 2 is enhanced to generate more H radicals, which contribute to the etching of organic low k film. In our experiments, the electron density increases with NH x concentration increasing. This could also contributes to increase NH x dissociation processes and to enhance oxide reduction at the surface layer, as it was previously related /7-8/.

Conclusion.

Investigations have been performed downstream a microwave discharge sustained in Ar-N 2 -H 2 gas mixtures, by means of mass spectrometry. In a first time we develop a method in order to correct the measured signal disturbed by wall desorption and dissociative ionization processes. In a second part we study the relative concentration of the main radicals and molecules produced in various gas compositions previously used for thin molybdenum layers nitriding treatment. Then we make correlations between the plasma chemistry and the plasma parameters measured by means of Langmuir probe spatially resolved in the plasma expansion.

Results show that with ternary mixtures (Ar-N 2 -H 2 ) a large amount of N 2 +H 2 can be mixed in the discharge (up to 50%), even at low power (200W) without plasma shrinking. In binary gas mixture plasma expansion cannot be obtained at such power for concentration larger than 5% for N 2 or 3% for H 2 mixed to Ar. This behaviour can be explained by the large amount of NH x radicals, and NH 3 produced. Due to the low ionization threshold values of these species compared to those of N 2 , H 2 and Ar, the production of these species gives rise to an increase of ionization process number and of the electron density. Consequently, the plasma expansion in the reactor is improved. Moreover, the electron density increase gives rise to an increase of dissociation processes, producing NH x radicals. As it has been already reported, in the case of H atoms, these radicals could act on the substrate layers, reducing oxide /7, 8/.

Figure 3

 3 Figure 3 displays the signal intensity of NH 2 taking into account the change of the emission current intensity versus the electron energy. The results are compared to the values calculated assuming only the direct ionization of NH 2 and using the ionization cross-section values given by Tarnovsky et al /17/. They are also compared to the calculated values obtained assuming the dissociative ionization process of ammonia and using the ionization cross section values given by Tarnovsky /17/ and Märk et al /16/. Results obtained using these two references are

3 and NH 2 ,

 2 Figures 4 and 5 give results obtained in the case of NH and N. Experimental values are compared versus the electron energy to the signal intensity calculated assuming only the direct ionization of the two radicals. We use the cross section values measured by Tarnovsky et al /17/ for NH and Brook et al /20/ for N. It can be seen a good agreement between experimental and calculated values for both species. This shows that the signal measured for m/q=14 and 15 result mainly of the direct ionization of NH and N.

3. 3 . 1 .

 31 The case of NH x species.Analysis, corresponding to different nitriding conditions, are performed downstream the discharge. Figure8and 9 displays the relative signal intensity measured by means of mass spectrometer, N/NH 3 , NH/NH 3 and NH 2 /NH 3 , versus %H 2 injected in Ar+N 2 binary mixture,

Figure 1 .

 1 Figure 1. Experimental set-up.

Figure 2 .

 2 Figure 2. Mass spectrometer signal intensity measured for NH 3 + (m/q=17) versus electron energy in the ionization chamber. Comparison with the signal intensity calculated using cross section values given in the literature in the cases of the direct ionization of NH 3 (Märk, Tarnovsky, Rejoub) and of the dissociative ionization of water producing OH + (m/q=17) (Straub).

Figure 3 .

 3 Figure 3. Mass spectrometer signal intensity measured for NH 2 + (m/q=16) versus electron energy in the ionization chamber. Comparison with the signal intensity calculated using cross section values given in the literature in the cases of the direct ionization of NH 2 (Tarnovsky)and of the dissociative ionization of NH 3(Märk, Tarnovsky).

Figure 4 .

 4 Figure 4. Mass spectrometer signal intensity measured for NH + (m/q=15) versus electron energy in the ionization chamber. Comparison with the signal intensity calculated assuming the direct ionization of NH using cross section values given in the literature (Tarnovsky).

Figure 5 .

 5 Figure 5. Mass spectrometer signal intensity measured for N + (m/q=14) versus electron energy in the ionization chamber. Comparison with the signal intensity calculated assuming the direct ionization of N using cross section values given in the literature (Brook).

Figure 6 .

 6 Figure 6. Mass spectrometer signal intensity measured for N 2 H 2 + (m/q=30) versus electron energy in the ionization chamber. Comparison with the signal intensity calculated assuming the dissociative ionization of hydrazine, using cross section values given in the literature (Syage).

Figure 7 .

 7 Figure 7. Contribution to the total signal intensity measured at m/q=17 and 16 of the direct ionization of NH 3 and NH 2 , versus electron energy.

Figure 8 .

 8 Figure 8. The relative signal intensity measured by means of mass spectrometer N/NH 3 , NH/NH 3 and NH 2 /NH 3 , versus %H 2 injected in Ar+N 2 binary mixture, in the case of the first gas mixture.

Figure 9 .

 9 Figure 9. The relative signal intensity measured by means of mass spectrometer N/NH 3 , NH/NH 3 and NH 2 /NH 3 , versus %H 2 injected in Ar+N 2 binary mixture, in the case of the second gas mixture.

Figure 10 .

 10 Figure 10. The relative concentration of ammonia compared to the argon (NH 3 /Ar), versus %H 2 mixed to Ar+N 2 , in the two gas mixtures (64%Ar-36%N 2 ) (circle) and (91%Ar-9%N 2 ) (square). The ternary mixtures pointed on the graph bring back to comments displayed in section 3.4.

Figure 11 .

 11 Figure 11. The relative concentration ratio (N 2 H 2 /NH 3 ) versus % H 2 injected in the binary mixtures (64%Ar-36%N 2 ) (full circles) and (91%Ar-9%N 2 ) (empty squares).

Figure 12 .

 12 Figure 12. The electron density measured by means of a langmuir probe versus distance from the tube discharge exit. Measurements are performed at 400W, in different gas mixtures.

Figure 13 .FIG. 1 .

 131 Figure 13. The electron energy distribution function measured by means of a langmuir probe and numerical simulation of harmonic components method, in the case of the three ternary gas mixtures under investigation.

  FIG.3.

  FIG.4.

  FIG.5.

2 FIG. 6 .

 26 FIG.6.

  FIG.7.

2 FIG. 8 .

 28 FIG.8.

  FIG.9.

  FIG.11.

  FIG.12.

  FIG.13.

  Figure 6 displays the signal intensity measured for m/q=30 versus the electron energy. As previously, intensities values are corrected taking into account the change due to the emission current in the ionization chamber.Comparison is performed with calculated values assuming only the dissociative ionization of Hydrazine (N 2 H 4 ), and using ionization cross section values given by Syage et al /21/. A good agreement is observed between the two curves. This show that the signal intensity measured for N 2 H 2 + could be ascribed to the dissociative ionization of hydrazine. However, N 2 H 4 + ion is not detected in these experiments and is expected to be easily detected at these electron energies when N 2 H 4 is produced /22/. Although, the signal intensity measured for N 2 H 2 + is the same that the value calculated assuming only the dissociative ionization of hydrazine, the previous remarks show that N 2 H 4 is not significantly produced in the present experiments and the signal measured for m/q=30 is probably mainly due to the direct ionization of N 2 H 2 . Comparison of experimental values to calculations assuming only the direct ionization of N 2 H 2 , is necessary in order to confirm the previous conclusion. However, no data about this ionization process is available in the literature. to the direct ionization of the neutral radical. However, in the case of NH 2 and NH 3 , it results of two contributions: The direct ionization and the dissociative ionization of larger species. These contributions depend on the electron energy in the ionization chamber. Before any analysis by means of mass spectrometer in nitriding conditions, it is necessary to correct the signal intensity measured for each species, removing the part due to the dissociative ionization of larger species. As previously shown on Figures 2, 3, 4 and 5, this

	The figure 7 displays the part due to the direct ionization of NH 3 and NH 2 compared to the
	total signal intensity measured for m/q=16 and 17, versus electron energy, and considering the
	cross section values given by Märk et al /16/ for NH 2 (case of an overestimation). For both
	species, the dissociative ionization contributes to less than 30% to the total signal for electron
	energy lower than 70 eV.

3.1.3. The case of N 2 H 2 (diimide)

According to Stothard et al /23/, the main channel reaction of NH 2 +NH 2 at low pressure and ambient temperature is producing N 2 H 2 and H 2 . The reaction rate constant is rather significant k=1.3.10 -12 cm 3 s -1 . But, the energy difference between NH 2 +NH 2 and N 2 H 4, is too large and a three body channel reaction must be considered in order to produce N 2 H 4 and to dissipate the energy excess. The third body could be no reactive species in the gas volume (k=7.0.10 -14 cm 3 s -1 in He at 1mbar) or at low pressure it is more probably the wall reactor. Considering a Ni(100) surface at temperature ranges from 200 to 450K, Huang et al /24/ have shown that the primary mechanism for diimide formation is the recombination of an adsorbed NH surface intermediate.

This first study shows that in the case of NH and N the signal measured by mass spectrometry is mainly due 3.2. The correction of the mass spectrometry signal.

  It is worth noting that in this reactive scheme, NH 3 is supposed produced by catalytic plasma reaction /9/ and reacts with the species produced in the plasma expansion because of the recirculation flow produced within the reactor.In the case where the surface reaction dominates, NH produced within the discharge, reacts on Figure11shows the relative concentration ratio (N 2 H 2 /NH 3 ) versus % H 2 injected in the two previous binary mixtures Ar+N 2 . It can be seen that the relative density of N 2 H 2 produced remains the same in both mixtures and that a maximum is observed at about 5% H 2 .According to Huang et al /24/, the primary mechanism for diimide formation is the recombination of an adsorbed NH surface intermediate. This mechanism can be written, NH s +NH s N 2 H 2.

	the reactor wall (stainless steel) according to /9, 10/,
	Two different reaction channels are involved to produce N x H y species in an Ar-N 2 -H 2 NH s + H NH 2s
	discharge. They are plasma catalysis processes (heterogeneous reactions) /9/ or volume
	reactions (homogeneous reactions) /25/. In the present experiments, the concentration of As previously, NH 3 results of the surface reactions,
	NH 2 >(NH or N). In a previous work /10/, we have shown that NH>NH 2 at low % N 2 injected
	(98.7%Ar-1.3%N 2 ). For larger N 2 %, NH concentration is similar to NH 2 concentration or NH 2s +H NH 3
	And even lower in N reactions, NH s +H 2 NH 3 .
	N+H 2 NH + H (k=4.10 -16 m 3 s -1 )
	NH 3 +H NH 2 + H 2 ((k=4.10 -18 m 3 s -1 at 1750K).
	Because of the large activation energy needs for the last reaction, radical concentrations are
	NH>NH 2 in the bulk plasma, where homogenous reaction dominate. However, the reaction
	rate can be enhanced by means of fast H atom injected within the plasma /25/. Similar

2 -H 2 gas mixture. This agrees with the present experiments. However, Helden et al /25/ obtained in N 2 -H 2 a NH concentration larger than NH 2 and they show that the NH concentration decreases and NH 2 concentration increases with %N 2 decreasing. This is in contradiction with our results. Such different behaviour can be understood considering a change in the dominant reactive mechanism involved. In the present experiments, the inner reactor diameter is 50mm and the pressure is 100Pa, whereas in the case of Helden et al /25/, the plasma is expanding supersonically in the reactor, which inner diameter is 400mm and the pressure is 20KPa in the cascaded arc plasma source and 20 Pa to 100 Pa in the reactor. Because of the low pressure and of the low reactor diameter, the dominant reaction channels are probably heterogeneous reaction in our experiments, whereas homogeneous mechanisms are dominant in the case of Helden et al /25/, where NH and NH 2 are mainly produced by the processes with fast ions have been related in the case of ion-molecule reaction /26/.

Where, NH s is adsorbed NH species.

As previously related, at low %N 2 injected, NH>NH 2 . This behaviour as been observed on Ni(100) surface by Huang et al /24/ and it can be explained by the large adsorption of H 2 on the wall, which stabilise NH and stop the surface reaction of NH s with H s , producing NH 2.

3.3.2. The case of N 2 H 2:

  , 9, 10), show that the two first gas mixtures correspond to plasma containing low density of NH x radicals and NH 3 molecules. Whereas the third gas mixture corresponds to plasma containing a large amount of NH x radicals and NH 3 molecules. In the particular case of the second gas mixture, N 2 and H 2 and NH x species are largely diluted in Ar, compared to the two other ones. These remarks show that the plasma expansion in the reactor is efficient in plasma containing N 2 or H 2 strongly diluted in Ar or in plasma containing a large concentration of N 2 and H 2 (more than 50%) and producing a large amount of ammonia and radicals NH x . This behaviour can be explained looking at the ionization threshold values of the different species observed in these plasmas. The ionization threshold values of Ar, N 2 and H 2 are 15.75, 15.58 and 15.426 eV respectively /28/. The ionization threshold of NH 3 is 10.166 eV, that of NH 2 is 11.22 eV and NH is 12.8 eV /28/. Moreover the ionization cross section values of NH x species are larger than ionization cross section values of N 2 and H 2 /16, 17, 18, 28/. Consequently, NH x radicals or NH 3 , increases the efficiency of ionization processes within the plasma, and therefore increase the electron density.Consequently, the plasma expansion is easier in gas mixture containing large amounts of NH x species. H. Nigai et al /29/, have observed that the admixture of a little amount of N 2 in H 2 ,

Acknowledgements: This work has been supported by La Région Limousin.


Michael Gb Blum
email: michael.blum@imag.fr

Mattias Jakobsson

Michael G B Blum

Deep divergences of human gene trees and models of human origins

Keywords: human origins, time to the most recent common ancestor, TMRCA, archaic admixture, African bottleneck, coalescent Running head: Deep divergences of human gene trees

Two competing hypotheses are at the forefront of the debate on modern human origins. In the first scenario, known as the recent Out-of-Africa hypothesis, modern humans arose in Africa about 100,000-200,000 years ago, and spread throughout the world by replacing the local archaic human populations. By contrast, the second hypothesis posits substantial gene flow between archaic and emerging modern humans. In the last two decades, the young time estimates -between 100,000 and 200,000 years -of the most recent common ancestors for the mitochondrion and the Y-chromosome provided evidence in favor of a recent African origin of modern humans. However, the presence of very old lineages for autosomal and X-linked genes has often been claimed to be incompatible with a simple, single origin of modern humans. Through the analysis of a public DNA sequence database, we find, similar to previous estimates, that the common ancestors of autosomal and X-linked genes are indeed very old, living, on average, respectively 1,500,000 and 1,000,000 years ago. However, contrary to previous conclusions, we find that these deep gene genealogies are consistent with the Out-of-Africa scenario provided that the ancestral effective population size was approximately 14,000 individuals. We show that an ancient bottleneck in the Middle Pleistocene, possibly arising from an ancestral structured population, can reconcile the contradictory findings from the mitochondrion on the one hand, with the autosomes and the X-chromosome on the other hand.

Introduction

The process by which modern humans arose has been the subject of much debate in paleoanthropology [START_REF] Stringer | Modern human origins -progress and prospects[END_REF]. Especially the extent of admixture between anatomically modern humans and archaic populations of Homo has been vigorously debated [START_REF] Wolpoff | Multiregional, not multiple origins[END_REF][START_REF] Templeton | Out of Africa again and again[END_REF][START_REF] Plagnol | Possible ancestral structure in human populations[END_REF][START_REF] Garrigan | Reconstructing human origins in the genomic era[END_REF][START_REF] Fagundes | Statistical evaluation of alternative models of human evolution[END_REF]. At one end of the spectrum, the recent Out-of-Africa hypothesis posits that a group of modern humans, arising in Africa about 100,000-200,000 years ago, spread throughout the world by replacing, without admixture, the local archaic human populations [START_REF] Mellars | Why did modern humans populations disperse from Africa ca. 60,000 years ago[END_REF]. At the other end, the Multi-Regional model posits substantial gene flow between local archaic humans and the emerging modern humans, even though it does not exclude Africa as the cradle of modern humans [START_REF] Wolpoff | Multiregional, not multiple origins[END_REF]. In between these two hypotheses, alternative scenarios assume archaic admixture restricted to Africa before the emerging modern humans eventually colonized the globe (Harding and McVean 2004;[START_REF] Gunz | Early modern human diversity suggests subdivided population structure and a complex out-of-Africa scenario[END_REF]).

Substantial genetic evidence has been put forward in support of the recent Outof-Africa hypothesis. For example, a continuous decrease of genetic diversity with increasing distance from Africa has been observed for autosomal microsatellites [START_REF] Prugnolle | Geography predicts neutral genetic diversity of human populations[END_REF][START_REF] Ramachandran | Support from the relationship of genetic and geographic distance in human populations for a serial founder effect originating in Africa[END_REF] and SNPs [START_REF] Li | Worldwide human relationships inferred from genome-wide patterns of variation[END_REF], as well as a continuous increase of linkage disequilibrium with distance from Africa [START_REF] Jakobsson | Genotype, haplotype and copy-number variation in worldwide human populations[END_REF]. These world-wide patterns are consistent with a serial founder model in which migrant populations are formed from a subset of the previous population in the migration wave outward from Africa [START_REF] Degiorgio | Explaining worldwide patterns of human genetic variation using a coalescent-based serial founder model of migration outward from Africa[END_REF]). Sex-linked markers have also provided evidence for a recent African origin since the common ancestor of all contemporary mitochondrial haplotypes existed as recently as ∼200,000 years ago [START_REF] Cann | Mitochondrial DNA and human evolution[END_REF], and the ancestor of all Y-chromosomes lived ∼100,000 years ago [START_REF] Thomson | Recent common ancestry of human Y chromosomes: evidence from DNA sequence data[END_REF][START_REF] Wilder | Genetic evidence for unequal effective population sizes of human females and males[END_REF].

Direct evidence against the recent Out-of-Africa hypothesis can potentially come from comparisons between ancient DNA of archaic humans, such as Neanderthals, and DNA of present-day modern humans [START_REF] Noonan | Neanderthal genomics and the evolution of modern humans[END_REF]. Recently, when releasing the draft sequence of a Neanderthal genome, [START_REF] Green | A Draft Sequence of the Neandertal Genome[END_REF] found 1 -4% of Neanderthal introgression in the gene pool of non-Africans; however, previous comparisons involving ancient Neanderthal DNA did not provide evidence in favor of admixture between humans and Neanderthals [START_REF] Krings | Neandertal DNA sequences and the origin of modern humans[END_REF][START_REF] Noonan | Sequencing and analysis of Neanderthal genomic DNA[END_REF][START_REF] Wall | Inconsistencies in neanderthal genomic dna sequences[END_REF]. Genetic evidence that does not involve ancient DNA, in particular elevated levels of Linkage Disequilibrium (LD) in modern humans, was also found to be indicative of ancient admixture [START_REF] Plagnol | Possible ancestral structure in human populations[END_REF][START_REF] Wall | Detecting ancient admixture and estimating demographic parameters in multiple human populations[END_REF]). It has additionally been argued that the presence of very old lineages, or deep divergences, for autosomal genes and genes on the X-chromosome is incompatible with a simple, single origin of modern humans, and that these deep divergences are instead evidence in favor of archaic admixture [START_REF] Harris | X chromosome evidence for ancient human histories[END_REF]Harding and McVean 2004;[START_REF] Garrigan | Deep haplotype divergence and long-range linkage disequilibrium at Xp21.1 provide evidence that humans descend from a structured ancestral population[END_REF][START_REF] Evans | Evidence that the adaptive allele of the brain size gene microcephalin introgressed into homo sapiens from an archaic homo lineage[END_REF][START_REF] Garrigan | Reconstructing human origins in the genomic era[END_REF][START_REF] Hayakawa | Fixation of the humanspecific CMP-N-acetylneuraminic acid hydroxylase pseudogene and implications of haplotype diversity for human evolution[END_REF][START_REF] Patin | Deciphering the ancient and complex evolutionary history of human arylamine N-acetyltransferase genes[END_REF][START_REF] Cox | Testing for archaic hominin admixture on the x chromosome: Model likelihoods for the modern human RRM2P4 region from summaries of genealogical topology under the structured coalescent[END_REF][START_REF] Kim | Population genetic analysis of the n-acylsphingosine amidohydrolase gene associated with mental activity in humans[END_REF].

A measure of the (deepest) divergence of a gene tree is the Time to the Most Recent Common Ancestor (TMRCA). The TMRCA is the time at which the Most Recent Common Ancestor (MRCA) of all existing copies of a given gene lived. A genome-wide frequency distribution of the TMRCAs has been reported by curating the literature [START_REF] Garrigan | Reconstructing human origins in the genomic era[END_REF] but no systematic and consistent analysis has been performed in a single genome-wide data set. We report the first genomewide estimation of the TMRCAs of anatomically modern humans, and we investigate if different scenarios of human evolutionary history are supported by this estimate.

In particular, we investigate to what extent the ages of the autosomal and X-linked lineages are compatible with the recent Out-of-Africa hypothesis.

Materials and Methods

Sequence data.

The data comprises of 40 re-sequenced independent intergenic regions from the autosomes and the X-chromosome [START_REF] Wall | A novel DNA sequence database for analyzing human demographic history[END_REF]. Each region encompasses approximately 20kb and consists of three 2kb sequence fragments, separated by 7kb of unsequenced DNA ("locus trio design", see [START_REF] Garrigan | Deep haplotype divergence and long-range linkage disequilibrium at Xp21.1 provide evidence that humans descend from a structured ancestral population[END_REF]. We present the results for 78 individuals from three African populations (Mandenka, Biaka, and San from Namibia), one European population (Basque), one Asiatic population (Han) and one population from Oceania (Melanesia). Two common chimpanzee sequences, available in the database, were used as outgroups.

TMRCA estimation.

We use the method of [START_REF] Thomson | Recent common ancestry of human Y chromosomes: evidence from DNA sequence data[END_REF] for estimating TMRCAs and we reiterate the basic motivation for the estimator (for an extensive description, see [START_REF] Thomson | Recent common ancestry of human Y chromosomes: evidence from DNA sequence data[END_REF]Hudson 2007). Thomson's estimator requires an outgroup that provides information on the ancestral state at every polymorphic position and assumes an infinite sites model.

Let T be the time to the MRCA for a sample of n lineages, and let x i be the number of mutations that have occurred between the MRCA and lineage i. We assume that the number of mutations along a branch, x i , is Poisson distributed with mean equal to the product of the mutation rate (u) and the branch length (T ). Then, Thomson's estimator of TMRCA is

T = n i=1 x i /(nu).
To estimate the mutation rate u, we assume a molecular divergence of 6 millions years between human and chimpanzee [START_REF] Glazko | Estimation of divergence times for major lineages of primate species[END_REF], and a generation time of 25 years. Computing the mean number of nucleotide differences between two chimp sequences and the human sequences, we find a mean mutation rate of 9.90×10 -10 /bp/year, on the same order as [START_REF] Fagundes | Statistical evaluation of alternative models of human evolution[END_REF].

To compute Thomson's estimator T , we reduce the data set to human polymorphic sites. Since we only consider sites that are polymorphic within humans, the mutations unique to the outgroup sequences are excluded from the computation of the TMRCA estimator (and have no effect on the TMRCA estimates). We reconstruct the ancestral state of each polymorphic site and, under the infinitely-many-sites mutation model, the state of the most recent common ancestor (MRCA) of the within-species sample is assumed to be equal to the outgroup allelic state. For 20 out of 1588 sites, the ancestral state could not be deduced, and, for these sites, the most frequent allele was assumed to be the ancestral variant. The choice of the infinitely-many-sites model is motivated by the small estimated mutation rate (u = 9.90 × 10 -10 /bp/year), which makes the probability of two mutations hitting the same site negligible.

We investigate the effect of recombination on Thomson's estimator using simulations, and we consider an effective recombination rate ranging from 0 to 10. This range encompasses previous estimates of the recombination rate in humans (The International HapMap Consortium 2005;[START_REF] Voight | Interrogating multiple aspects of variation in a full resequencing data set to infer human population size changes[END_REF]Coop et al. 2008b). Using, for instance, a standard estimate of 1cM/Mb (i.e. 10 -8 /bp) for the recombination rate (The International HapMap Consortium 2005) gives an effective recombination rate of ρ = 8 for a 20kb region, assuming an effective population size of 10, 000 individuals. To quantify the effect of recombination on T , we compute the relative bias and mean square-error of T as a function of the recombination rate (Figure 1).

We also investigate if Thomson's estimator T remains accurate when the population experienced demographic changes, and we consider three different demographic models: a) a population with a constant size of N = 10, 000 individuals; b) an expanding population where the population experienced a 10-fold expansion starting at a time distributed uniformly between 0 and 200, 000 years ago from a population-size of N = 10, 000 individuals; and c) a population splitting model. For the last model c, we assume that an ancestral population (of size N = 10, 000 individuals) split at a time distributed uniformly between 0 and 200, 000 years ago into two subpopulations. One of the two subpopulations contains N = 10, 000 individuals and the other population contains N = 10, 000 × p individuals with p chosen from a uniform distribution between 0 and .5.

Approximate Bayesian Computation (ABC).

We use ABC to find the range of demographic parameters that yield TMRCAs similar to the empirical estimates of TMRCAs. For each scenario of human evolution (described in the Results and Discussion Section), the ABC statistical procedure can be described as follows:

• Generate the demographic parameters according to the prior distributions given in Supplementary Table 1.

• Simulate sequence data with the software ms [START_REF] Hudson | Generating samples under a Wright-Fisher neutral model of genetic variation[END_REF]. One simulation comprises of generating 20 autosomal and 20 X-linked sequence-regions with the same number of samples and the same sequence lengths as in the empirical data.

• Compute the summary statistics (see below) for the simulated sequences and compute the Euclidean distance between observed and simulated summary statistics.

After performing a total of 100,000 simulations, we retain the 500 simulations that provide the best match to the data. We use an Epanechnikov kernel to assign larger weights to simulations that provide the best match [START_REF] Beaumont | Approximate Bayesian computation in population genetics[END_REF]. To account for the imperfect match between simulated and observed summary statistics, we then use regression adjustment as described by [START_REF] Blum | Non-linear regression models for approximate Bayesian computation[END_REF]. After completion of the algorithm, the posterior distribution of the parameters consists of the set of accepted parameters after adjustment.

To compute the summary statistics, we consider, for each of the 40 sequence-regions, the mean number of mutations, n i=1 x i /n, between the set of sequences and the ancestral sequence. To reduce the number of summary statistics, we compute, separately for the X-linked and the autosomal markers, the three quartiles of the 20 mean numbers of mutations. This procedure results in a total of 6 summary statistics.

Choice of priors for the mutation and recombination rates.

For the mutation rate, we choose an empirical Bayes approach in which the prior depends on the data [START_REF] Casella | An introduction to empirical Bayes data analysis[END_REF]. We choose a Gamma distribution for the mutation rate and we estimate the parameters so that the Gamma distribution fits the empirical distribution of the 40 estimated mutation rates. We obtain a Gamma distribution with a shape parameter of 15.18 and a scale parameter of 6.50 × 10 -11 . This results in a median mutation rate of 9.65 × 10 -10 mutations/bp/year and 95% of the simulated mutation rates are between 5.54 × 10 -10 and 1.54 × 10 -9 mutations/bp/year. For the crossing-over rate, we consider a log-normal distribution with a mean and standard deviation (on a log scale) of -18.148 and 0.5802 [START_REF] Voight | Interrogating multiple aspects of variation in a full resequencing data set to infer human population size changes[END_REF]. We assume homogeneous cross-over rate along each sequence-region.

Computation of the relative model probabilities.

We introduce an indicator variable Y = {1, 2, 3, 4} that indicates, for each simulation, which one of the four different demographic scenarios (described in the Results and Discussion Section) generated the data. We perform the same number of simulations for each scenario. We then regress the indicator variable Y by the 6 summary statistics Posterior predictive simulations.

To perform goodness-of-fit of the scenarios of human evolution (Figure 2), we perform posterior predictive checks [START_REF] Gelman | Bayesian Data Analysis, Second Edition[END_REF]. We sample, with replacement, 10,000 multivariate demographic parameters at random from the posterior distribution obtained with the ABC algorithm. Using ms [START_REF] Hudson | Generating samples under a Wright-Fisher neutral model of genetic variation[END_REF], we simulate, for each multivariate demographic parameter, gene trees along a 20kb sequence, and compute the median of the (potentially different) simulated TMRCAs found along the sequence.

For each scenario, this results in a total of 10,000 median TMRCAs that are displayed in Figure 2.

Results and Discussion

Genome-wide estimation of the TMRCAs.

The data comprise 40 re-sequenced independent intergenic regions from the autosomes and the X-chromosome provided by a public DNA sequence database that has been designed for the purpose of analyzing human prehistory [START_REF] Wall | A novel DNA sequence database for analyzing human demographic history[END_REF]). To compute genome-wide estimation of the TMRCA, we consider the statistic of [START_REF] Thomson | Recent common ancestry of human Y chromosomes: evidence from DNA sequence data[END_REF], which has been used for instance to date the ancestor of the human Ychromosome [START_REF] Thomson | Recent common ancestry of human Y chromosomes: evidence from DNA sequence data[END_REF] and the human FOXP2 gene (Coop et al. 2008a).

For each DNA fragment, the TMRCA estimate is obtained by computing the number of mutations between each sample (gene-copy) and a reconstructed ancestral sequence,

and averaging across the gene-copies of the sample. Generally, the bias of this estimator has been shown to be small (Hudson 2007), and we demonstrate that it is also robust to recombination (see Figure 1). To affix a time scale in years to the TMRCA estimates, we assume a molecular divergence of 6 millions years between human and chimpanzee [START_REF] Glazko | Estimation of divergence times for major lineages of primate species[END_REF], and a generation time of 25 years. For each region of 20kb, we compute one TMRCA estimate so that this estimate captures an average TMRCA for the possibly different TMRCAs found along the 20kb region.

Figure 2 displays the distribution of the TMRCAs for the 20 autosomal loci (Figure 2A) and the 20 loci located on the X-chromosome (Figure 2B). The median of the TMRCA is approximately 1,500,000 years for the autosomes (first quartile = 950,000, third quartile = 1,700,000) and approximately 1,000,000 years for the X-chromosome (first quartile = 700,000, third quartile = 1,350,000). These numbers are close to what [START_REF] Garrigan | Reconstructing human origins in the genomic era[END_REF] found by collecting TMRCA estimates from the literature. However, they are at odds with numerical predictions obtained by [START_REF] Fagundes | Statistical evaluation of alternative models of human evolution[END_REF] for the recent Out-of-Africa model. They found that the distribution of autosomal TMRCAs should peak around the time when modern humans emerged (100,000 to 200,000 years ago) and that 50% of the TMRCAs should be more recent than 650,000 years. In contrast, we find that only two out of 20 autosomal TMRCAs are more recent than 650,000 years.

A number of authors have argued that deep genealogical histories are incompatible with the recent Out-of-Africa hypothesis, and instead claimed that these deep genealogies are evidence in favor of archaic admixture [START_REF] Harris | X chromosome evidence for ancient human histories[END_REF]Harding and McVean 2004;[START_REF] Garrigan | Deep haplotype divergence and long-range linkage disequilibrium at Xp21.1 provide evidence that humans descend from a structured ancestral population[END_REF][START_REF] Evans | Evidence that the adaptive allele of the brain size gene microcephalin introgressed into homo sapiens from an archaic homo lineage[END_REF][START_REF] Garrigan | Reconstructing human origins in the genomic era[END_REF][START_REF] Hayakawa | Fixation of the humanspecific CMP-N-acetylneuraminic acid hydroxylase pseudogene and implications of haplotype diversity for human evolution[END_REF][START_REF] Patin | Deciphering the ancient and complex evolutionary history of human arylamine N-acetyltransferase genes[END_REF][START_REF] Cox | Testing for archaic hominin admixture on the x chromosome: Model likelihoods for the modern human RRM2P4 region from summaries of genealogical topology under the structured coalescent[END_REF][START_REF] Kim | Population genetic analysis of the n-acylsphingosine amidohydrolase gene associated with mental activity in humans[END_REF]. In the following sections, we investigate to what extent the genome-wide distribution of the TMRCAs is compatible with the recent Out-of-Africa model. We also consider alternative models that assume archaic admixture and check if they provide a better fit to the TMRCA distribution.

TMRCA distributions predicted by different scenarios of human evolution.

We compare observed TMRCAs to simulated TMRCAs for four different scenarios of modern human origin (Figure 3 and Supplementary Table 1). The two first scenarios are versions of the Out-of-Africa model, and the last two scenarios are versions of the Multi-Regional model: 3. 'Recent Admixture' A Multi-Regional scenario in which archaic and modern humans were isolated during 300-600,000 years and admixed recently in Eurasia, 30-70,000 years ago [START_REF] Plagnol | Possible ancestral structure in human populations[END_REF], and 4. 'Long-Standing Admixture' A Multi-Regional scenario with continuous and long-standing admixture between archaic and modern humans.

To find the range of demographic parameter-values that yield TMRCAs similar to the 40 estimates from the empirical data, we use Approximate Bayesian Computation (ABC) [START_REF] Beaumont | Approximate Bayesian computation in population genetics[END_REF][START_REF] Csilléry | Approximate Bayesian Computation in practice[END_REF]) and coalescent simulations [START_REF] Hudson | Generating samples under a Wright-Fisher neutral model of genetic variation[END_REF]. The empirical genetic data is summarized by six summary statistics that capture the divergence of gene trees. For each DNA fragment, we compute, for the human polymorphic sites, the mean number of mutations between the gene-copies of the sample and the reconstructed ancestral sequence (see Materials and Methods). We compute -separately for the X-chromosome and the autosomes -the three quartiles of these averaged number of mutations.

We estimate the ancestral effective population size of archaic humans (Figure 4 and Supplementary Table 2) and find, for the recent Out-of-Africa scenario, that the most likely value is 14,500 (95% CI = 12,000-17,000) similar to previous estimates [START_REF] Takahata | Allelic genealogy and human evolution[END_REF][START_REF] Harding | Archaic African and Asian lineages in the genetic ancestry of modern humans[END_REF][START_REF] Wall | Estimating ancestral population sizes and divergence times[END_REF][START_REF] Voight | Interrogating multiple aspects of variation in a full resequencing data set to infer human population size changes[END_REF]). In scenario 2, which assumes a structured population in Africa before the emergence of modern humans, we find a slightly lower estimate on the order of 8,000 individuals (95% CI = 5,000-15,000)

reflecting that several archaic African populations contribute to the modern gene pool.

In addition to parameter inference, the ABC approach also offers a convenient way to assign a probability to each of the scenarios [START_REF] Fagundes | Statistical evaluation of alternative models of human evolution[END_REF][START_REF] François | Demographic History of European Populations of Arabidopsis thaliana[END_REF][START_REF] Verdu | Origins and genetic diversity of pygmy huntergatherers from Western Central Africa[END_REF]. We find that the four different models are almost equally supported by the divergence of gene trees since the four posterior probabilities range between 20% and 30% (Supplementary Figure 1). These even probabilities reflect that the relatively ancient lineages found in the autosomes and X-linked genes neither favor nor disfavor the models with archaic admixture (models 2-4).

To check if the different scenarios of human evolution provide a good fit to the data, we compare the empirical TMRCA estimates to the TMRCAs predicted by the different scenarios. All four models predict, on average, lineages as old as 1,500,000 years for autosomal fragments (Figure 2A) and as old as 1,000,000 years for X-linked fragments (Figure 2B). In short, we find that both the simple replacement model and the models with archaic admixture are perfectly compatible with the deep divergences found in the empirical data.

However, not all aspects of the empirical TMRCAs are well captured by the modeled scenarios of human evolution. The variance of the empirical TMRCAs is larger than the variance predicted by three of the four different models of human evolution (see Figure 2 and Supplementary Table 3), and this large variance has been interpreted as the result of archaic sub-structure in Africa (Harding and McVean 2004). Indeed, the Multiple Archaic Populations' (scenario 2) shows similar variance of TMRCAs as the empirical data, but the inflated variance of the empirical TMRCA estimates can also be due to variation in mutation or recombination rate across the 40 sequence-regions [START_REF] Mcvean | The fine-scale structure of recombination rate variation in the human genome[END_REF]).

The mitochondrion and the Y-chromosome.

We also investigate the distribution of TMRCAs that is expected for the Y-chromosome and the mitochondrial chromosome (Figure 2C). The models of human evolution typically predict older TMRCAs compared to the estimated 170,000 years for mtDNA [START_REF] Ingman | Mitochondrial genome variation and the origin of modern humans[END_REF] and the upper estimate of 100,000 years for the Y-chromosome [START_REF] Tang | Frequentist estimation of coalescence times from nucleotide sequence data using a tree-based partition[END_REF][START_REF] Wilder | Genetic evidence for unequal effective population sizes of human females and males[END_REF][START_REF] Shi | A Worldwide Survey of Human Male Demographic History Based on Y-SNP and Y-STR Data from the HGDP-CEPH Populations[END_REF]. For mtDNA, a TMRCA of 170,000 years is within the range of values predicted by the 'Multiple Archaic Populations' scenario (P(TMRCA < 170, 000) = 0.21), but the mitochondrial TMRCA estimate is difficult to reconcile with the remaining three scenarios (P < 4 × 10 -2 ). For the Y-chromosome, a TMRCA of 100,000 years is clearly at odds with three of the models (P < 6 × 10 -4 ), but for the 'Multiple Archaic Populations' scenario with archaic African admixture, the proportion of simulated gene trees with TMRCAs younger than 100,000 years is larger than for the other three models, albeit quite small (P = 1.5×10 -2 ). Although assuming an archaic structured population in Africa, as in scenario 2, reduces the Y-chromosome and mtDNA TMRCAs, the model cannot fully explain a Y-chromosome ancestor living as recently as 100,000 years ago. In scenario 2, we consider three archaic populations, and increasing the number of archaic populations will further decrease the effective population size of each sub-population, which will decrease the predicted haploid TM-RCAs, bringing them in line with the empirical estimates. However, since the mtDNA and the Y-chromosome are non-recombining units, their young TMRCAs can also be explained by recent selective sweeps, caused by directional selection at any gene within the non-recombining units [START_REF] Kreitman | Methods to detect selection in populations with applications to the human[END_REF][START_REF] Jobling | The human Y chromosome: an evolutionary marker comes of age[END_REF].

A bottleneck when anatomically modern humans emerged

An alternative explanation for young haploid TMRCA involves a demographic bottleneck concomitant with the emergence of anatomically modern humans. The earliest known suite of derived morphological traits associated with modern humans appears in fossil remains from Ethiopia dated to 150-190 kya [START_REF] White | Pleistocene homo sapiens from middle awash, Ethiopia[END_REF][START_REF] Mcdougall | Stratigraphic placement and age of modern humans from kibish, Ethiopia[END_REF]. A model of the origin and spread of modern humans proposed by [START_REF] Lahr | Multiple dispersals and modern human origins[END_REF] assumes that the emerging modern humans experienced bottlenecks within Africa during the penultimate glacial age (130-190 kya) when cold, dry climates caused isolation of populations within Africa (see also [START_REF] Ambrose | Late Pleistocene human population bottlenecks, volcanic winter, and differentiation of modern humans[END_REF]). Here we consider a bottleneck that occurred 150,000 years ago, and we measure the bottleneck intensity by the inbreeding coefficient during the bottleneck

F = 1 -(1 - 1 2bN A ) D ,
where bN A is the diploid population size during the bottleneck and D corresponds to the duration (in generations) of the bottleneck. To give a reference value, the inbreeding coefficient F corresponding to the Out-of-Africa bottleneck was inferred at F = 0.175 [START_REF] Akey | Population history and natural selection shape patterns of genetic variation in 132 genes[END_REF][START_REF] Voight | Interrogating multiple aspects of variation in a full resequencing data set to infer human population size changes[END_REF]. Considering an uniform prior between 0 and 1 for F , we find that there is a large range of parameter values compatible with the autosomal and X-linked TMRCAs (Figure 5). For instance, both pairs of parameter values (F = 0, N A = 14, 000) and (F = 0.4, N A = 18, 000) are clearly within the range of the bivariate posterior distribution found with our ABC approach. The relatively large bivariate posterior range of the inbreeding coefficient F and of the ancestral size N A shows that a strong bottleneck can still produce an average autosomal TMRCA of 1.5 million year provided that the ancestral size was large enough before the bottleneck.

To investigate if a bottleneck 150,000 years ago in Africa can account for both recent haploid and old autosomal ancestors, we simulate TMRCAs for the different chromosomes by sampling the demographic parameters from the posterior distribution that was obtained using ABC. We find that the recent mtDNA TMRCA is clearly within the range of predicted values (P=0.37, Figure 2C) as well as the old X-linked and autosomal TMRCAs (Figure 2A-B). The fact that a bottleneck can accommodate an 8-fold discrepancy between autosomal and mtDNA TMRCA can be seen by plotting TMRCAs as a function of F and N A (Figure 6). Setting the inbreeding coefficient at F = 0.4 for instance, Figure 6 shows that the mean mtDNA TMRCA is smaller than 200,000 years, for a large range of values of N A , whereas the 20kb-averages of TMRCAs, for autosomal and X-linked markers, are older than 1.5 and 1 million years when N A > 16, 000 individuals. Finally, as for the other models of human evolution, the TMRCA of 100,000 years for the Y-chromosome remains unexpectedly young (P = 4 × 10 -4 ).

We find that both the 'multiple archaic population' model and a sudden bottleneck, 150,000 years ago in Africa, can account for the 8-fold discrepancy between autosomal and mtDNA TMRCA. Although modeling different patterns of human evolution, these two scenarios are different versions of a bottleneck in the human lineage that arose before the succeeding migration out of Africa. Previous attempts to detect an ancestral African bottleneck have often been inconclusive and the genetic evidence in favor of the Out-of-Africa bottleneck are more salient [START_REF] Marth | The allele frequency spectrum in genome-wide human variation data reveals signals of differential demographic history in three large world populations[END_REF][START_REF] Voight | Interrogating multiple aspects of variation in a full resequencing data set to infer human population size changes[END_REF]).

However, it is more difficult to detect this potential middle-Pleistocene bottleneck compared to the Out-of-Africa bottleneck since it is more ancient [START_REF] Depaulis | Power of neutrality tests to detect bottlenecks and hitchhiking[END_REF].

Additionally, the unexpected large levels of African LD that has been interpreted as evidence of archaic admixture [START_REF] Plagnol | Possible ancestral structure in human populations[END_REF][START_REF] Wall | Detecting ancient admixture and estimating demographic parameters in multiple human populations[END_REF]) may potentially be explained by an ancient bottleneck [START_REF] Schmegner | Genetic variability in a genomic region with long-range linkage disequilibrium reveals traces of a bottleneck in the history of the european population[END_REF]. We therefore suggest that an ancestral bottleneck, possibly following ancestral sub-structure and admixture, is important to consider when investigating demographic models of human evolution [START_REF] Schaffner | Calibrating a coalescent simulation of human genome sequence variation[END_REF][START_REF] Gutenkunst | Inferring the joint demographic history of multiple populations from multidimensional SNP frequency data[END_REF][START_REF] Laval | Formulating a historical and demographic model of recent human evolution based on resequencing data from noncoding regions[END_REF].

A simple population genetic prediction

Although we perform coalescent simulations, standard population genetic theory can predict the TMRCAs found for the autosomes and the X-chromosome in the Out-of-Africa model when there is no ancestral bottleneck or admixture. For a population of constant diploid size N in which the effective number of males and females is the same, the expected waiting time (in generations) before the coalescence of a (reasonably large) sample of genes is approximately 4N for the autosomes and 3N for the X-linked genes (see e.g. [START_REF] Hein | Gene Genealogies, Variation and Evolution[END_REF]. Using a generation time of 25 years and an effective population size of 14,000 individuals, the computation leads to an average TMRCA of 1,400,000 years for autosomal genes and 1,050,000 years for X-linked genes, which are both very close to our estimates from the sequence data. This theoretical argument

shows that the difference between the TMRCAs of the autosomes and the X-linked genes is easily explained by the difference in effective population size. The same argument yields an average TMRCA of 350,000 years for both the mtDNA and Y-chromosome, clearly deviating from the estimates based on empirical data [START_REF] Ingman | Mitochondrial genome variation and the origin of modern humans[END_REF][START_REF] Wilder | Genetic evidence for unequal effective population sizes of human females and males[END_REF][START_REF] Shi | A Worldwide Survey of Human Male Demographic History Based on Y-SNP and Y-STR Data from the HGDP-CEPH Populations[END_REF].

Conclusion

We provide 20 autosomal and 20 X-linked estimates of TMRCAs of a sample of contemporary humans and find that the autosomal ancestors of modern humans lived ∼1,500,000 years ago and that the X-linked ancestors lived ∼1,000,000 years ago. The ranges of values for the TMRCAs are quite large: 450,000-2,400,000 years for the autosomes, and 380,000-2,000,000 for the X-chromosome. These values are in the same range as previous estimates for autosomal and X-linked genes (see e.

  using local multinomial logistic regression to obtain the model probabilities P (Y |s) as a function of the summary statistics s (Fagundes et al. 2007; Beaumont 2008). Local logistic regression differs from standard logistic regression because larger weights are given to the simulations for which the summary statistics are close to the observed summary statistics. By computing the logistic regression equation for the observed summary statistics P (Y |s = s obs ), we obtain the relative model probabilities. To perform local multinomial logistic regression, we use the R package VGAM (Yee 2010).

  1. 'Single Origin Population' A recent Out-of-Africa scenario in which modern humans descended from one subpopulation of archaic humans that was a separate population for a long time in Africa, 2. 'Multiple Archaic Populations' A recent Out-of-Africa scenario in which different archaic African populations were connected by gene flow, even though only one archaic population eventually colonized the globe(Harding and McVean 2004;[START_REF] Garrigan | Reconstructing human origins in the genomic era[END_REF][START_REF] Campbell | African genetic diversity: Implications for human demographic history, modern human origins, and complex disease mapping[END_REF]

  g.[START_REF] Templeton | Out of Africa again and again[END_REF][START_REF] Garrigan | Reconstructing human origins in the genomic era[END_REF][START_REF] Tishkoff | Human origins within and out of Africa[END_REF]. We investigate to what extent the recent Out-of-Africa model reproduces the pattern of estimated TMRCAs, and when setting the ancestral effective size of humans to ∼14,000, this model reproduces the old TMRCAs of the empirical data. Deep divergences in human gene trees are therefore not incompatible with the recent Out-of-Africa hypothesis and the observation of deep gene genealogies should not be taken as evidence for the Multi-Regional hypothesis. Finally, we show that an ancestral bottleneck in Africa, possibly arising in a structured population, can account for the unexpectedly large discrepancy between young mtDNA and Y-chromosome ancestors and old autosomal and X-linked ancestors.

Figure 1 :Figure 2 :Figure 4 :Figure 5 :Figure 6 :

 12456 Figure1: Average errors of the TMRCA estimate for n = 100 individuals. We compare estimated TMRCAs T to true averaged TMRCAs T av . The true TMRCAs are averaged using the median along the 20kb sequences. The relative bias is defined as the average over simulations of ( T -T av )/T av and the relative quadratic error is the average of ( T -T av ) 2 /T 2 av . The effective recombination rate is equal to 4 times the ancestral population size times the recombination rate. In all simulations, we use the mutation rate estimated from the data.

Acknowledgments

This work was supported by a Grant from the Swedish Foundation for International Cooperation in Research and Higher Education (STINT) to MJ and MGBB. MJ would like to thank the Swedish Research Council and the Erik Philip-Sörensen's foundation for financial support.

recent Out-of-Africa scenario in which modern humans descended from one subpopulation of archaic humans that was a separate population for a long time in Africa. The recent Out-of-Africa scenario potentially includes a bottleneck before the exodus from Africa, 150,000 years ago. 2) 'Multiple Archaic Populations' A recent Out-of-Africa scenario in which different archaic African populations were connected by gene flow, even though only one archaic population eventually colonized the globe (Harding and McVean 2004;[START_REF] Garrigan | Reconstructing human origins in the genomic era[END_REF]. 3) 'Recent Admixture' A Multi-Regional scenario in which archaic and modern humans were isolated during 300-600,000 years and admixed recently in Eurasia, 30-70,000 years ago [START_REF] Plagnol | Possible ancestral structure in human populations[END_REF]. 4) 'Long-Standing Admixture' A Multi-Regional scenario with continuous and long-standing admixture between archaic and modern humans. Ancestral population size: N A , time of the migration out of Africa: T , inbreeding coefficient during the bottleneck: F , time of structuring of archaic African population: T + T S + T M , ending of structured archaic African population: T + T S , time of archaic humans exiting Africa: T 0 , time of admixture: T admix , and migration rate: M .


Éric Chassande-Mottin

Miriam Miele

Satya Mohapatra

Laura Cadonati

Detection of GW bursts with chirplet-like template families

Gravitational Wave (GW) burst detection algorithms typically rely on the hypothesis that the burst signal is "locally stationary", that is it changes slowly with frequency. Under this assumption, the signal can be decomposed into a small number of wavelets with constant frequency. This justifies the use of a family of sine-Gaussian wavelets in the Omega pipeline, one of the algorithms used in LIGO-Virgo burst searches. However there are plausible scenarios where the burst frequency evolves rapidly, such as in the merger phase of a binary black hole and/or neutron star coalescence. In those cases, the local stationarity of sine-Gaussians induces performance losses, due to the mismatch between the template and the actual signal. We propose an extension of the Omega pipeline based on chirplet-like templates. Chirplets incorporate an additional parameter, the chirp rate, to control the frequency variation. In this paper, we show that the Omega pipeline can easily be extended to include a chirplet template bank. We illustrate the method on a simulated data set, with a family of phenomenological binary black-hole coalescence waveforms embedded into Gaussian LIGO/Virgo-like noise. Chirplet-like templates result in an enhancement of the measured signal-to-noise ratio.

Motivations

Current searches for gravitational wave transients in LIGO-Virgo data focus on two signal classes: short unmodelled bursts and longer quasi-periodic signals from inspiralling black hole and/or neutron star binaries as predicted by post-Newtonian approximations. To account for intermediate scenarios, we consider "chirping burst" GW target signals that exhibit characteristics from both the above categories: a short duration and a "sweeping" frequency.

We propose here an extension of the Omega pipeline [START_REF] Chatterji | The search for gravitational-wave bursts in data from the second LIGO science run[END_REF] (originally known as Q-pipeline) that searches for chirping bursts. The Omega pipeline projects the data over a family of sine-Gaussian wavelets with fixed frequency. The idea is to replace these templates by frequency varying waveforms, referred to as chirplets.

In this paper, we first define chirplets and the related chirplet transform. We discuss the implementation of the chirplet transform and its insertion into the Omega pipeline, with attention to how the chirplet template bank is built. Finally, we present a few examples using simulated data.

ψ(τ ) ≡ A exp - (2πf ) 2 Q 2 (τ -t) 2 exp 2πi f (τ -t) + d/2 (τ -t) 2 , (1)
with A = (8πf 2 /Q 2 ) 1/4 , a normalization factor ensuring that |ψ| 2 = 1. t and f are the center time and frequency, respectively and Q is the dimensionless quality factor. See in Fig. 1 for an example of a chirplet. The main difference from a sine-Gaussian wavelets is the chirp rate, an additional term in the phase denoted d that changes the chirplet frequency linearly in time as

f (τ ) = f + d(τ -t).
The chirp rate controls the slope of the frequency evolution. When d = 0, we retrieve standard sine-Gaussians. Chirplets are thus associated with a four-dimensional parameter space instead of three for sine-Gaussians. In the sequel, we will concatenate all those parameters into a descriptor θ ≡ {t, f, Q, d}.

Chirplet transform

The chirplet transform T is obtained by correlating the data with the chirplets defined in the previous section. In the frequency domain, it reads:

T [x; θ] = X(ξ)Ψ * (ξ; θ)dξ 2 , (2)
where X(•) and Ψ(•; θ) denotes the Fourier transform of the (whitened) data stream x(•) and chirplet ψ(•) with descriptor θ resp. The chirplet Fourier transform can be expressed as

Ψ(ξ; θ) = A exp - Q2 4 (ξ -f ) 2 f 2 , (3)
where A = [( Q4 /Q2 )/(2πf 2 )] 1/4 is written in terms of a "complex-valued" quality factor Q = Q √ z/|z| where z = 1 + id∆ 2 t with the chirplet duration

1 ∆ t = Q/ (2 √ πf ).

Building template banks with chirplets

By varying the chirplet parameters, we obtain a continuous signal space. In this space, we need to select a finite-size family of representatives which will be used to analyze the data. The coverage of the chirplet space has to meet two conflicting goals i.e., satisfy a worst case mismatch with a minimum number of templates. We adopt the method proposed in [START_REF] Owen | Matched filtering of gravitational waves from inspiraling compact binaries: Computational cost and template placement[END_REF][START_REF] Chatterji | The search for gravitational-wave bursts in data from the second LIGO science run[END_REF] which consists of sampling the space with equi-spaced templates using the intrinsic metric deduced from the mismatch T ψ(θ ); θ between two neighboring templates with a small discrepancy δθ ≡ θ -θ in their parameters. The metric results from the second-order expansion of the mismatch µ θ (δθ) ≡ 1 -T [ψ(θ + δθ); θ] ≈ δs 2 when δθ → 0 and leads to 2 :

δs 2 = Q 4 d 2 + 16π 2 f 4 4Q 2 f 2 δt 2 + 2 + Q 2 4f 2 δf 2 + δQ 2 2Q 2 + Q 4 128π 2 f 4 δd 2 - Q 2 d 2f 2 δtδf - δf δQ Qf . (4)
There are several differences and additional terms from the sine-Gaussian metric, due to the non-zero chirping rate. Along the time axis and for small f Q √ d, the sampling step δt ∝ f/(Qd) is finer than that of the sine-Gaussian case δt ∝ Q/f . We note also that the sampling step along the chirping rate axis scales with δd ∝ (f/Q) 2 . We thus expect to get many chirplets in the low-frequency band and for large values of Q.

The chirplet space equipped with the above metric (off-diagonal terms being neglected for simplicity [START_REF] Chatterji | The search for gravitational-wave bursts in data from the second LIGO science run[END_REF]) can be discretized by a cubic lattice with templates placed at the vertices. The worst case occurs when the real signal is farther apart from all vertices, at the center of the cube. Let us denote δs = µ 1/2 , this worst-case distance, which corresponds to the half-length of the cube diagonal and assign a maximum value µ 1/2 max that we can tolerate. Since we are in a four-dimensional space, the length of the cube edge is equal to that of its half-diagonal. Therefore, we must have ≤ √ µ max . The discretization along each axis of the parameter space which results directly from this condition ensures that µ θ (θ -θ n ) ≤ µ max for any θ with θ n the closest vertex of the lattice. In the following, we set the maximum mismatch to the value µ max = 0.2 typically used when applying standard Omega. Fig. 2 shows an example of a chirplet template bank resulting from this template placement scheme.

In Fig. 3, we apply the same scheme in two different settings. In both cases we computed the number of templates necessary to cover the signal space in the sine-Gaussian (standard Omega) and chirplet (chirpletized Omega) cases. This computation is done at a fixed time t. We compare the result to the estimate given by the ratio of the whole space volume V = |µ | 1/2 d3 θ (where θ = {f, Q, d} and δs 2 = |µ | denotes the metric in Eq. ( 4) without the components associated to the time axis) to the size of a cubic element of the lattice. We find 3 :

N ≡ V/ 3 ∝ f -2 min Q 3 max d max , (5)
where we assume that for each coordinate the lower boundaries (min) are much smaller than the higher boundaries (max).

It is important to note that both the count and estimate are obtained assuming an infinite bandwidth. Since the data are sampled, we are restricted to a limited Nyquist frequency. Chirplets with frequency exceeding this limits are aliased and have to be discarded. Fig. 3 also show the number of non-aliased chirplets. This number is about a factor of 10 larger than the number of sine-Gaussians required to cover the entire parameter space. Two comments can be made at this point: the larger number of templates indicate that with chirplets, we define and explore a much larger signal space than with sine-Gaussians (we investigate this question further in Sec. 5). As the computing cost scales approximately linearly with the number of templates, analyzing the data with chirplets requires with a ten-fold increase in computing resources (a factor that will be rapidly absorbed by the exponential growth of computing power).

1 By definition, ∆t ≡ 2 √ π R (τ -t) 2 ψ 2 (τ )dτ .

From "standard" to "chirpletized" Omega

In this session we discuss other aspects of the analysis pipeline, in addition to the implementation of the chirplet template bank.

Filtering

The modulus of Ψ(•) in Eq. ( 3) is a Gaussian function as in the sine-Gaussian case. This allows the use of the same filtering scheme as in the standard Omega pipeline to generate the chirplet transform. The Omega scheme [START_REF] Chatterji | The search for gravitational-wave bursts in data from the second LIGO science run[END_REF] operates in the frequency domain following Eq. ( 2). It consists in multiplying the Fourier transform of the data, computed with the FFT algorithm, with that of the templates and take the inverse Fourier transform of the product. Omega uses a bi-square frequency window that approximates the Gaussian shape. The compact support of the bi-square window prevents aliasing.

This scheme can be applied to the chirplet case with two differences. First, the template Ψ is now complex, thus we need to multiply the data spectrum both in modulus and phase. Second, the template bandwidth now results from the quadratic sum ∆ 2 f = (∆ finite size f . The results are compared to the size of the sine-Gaussian template bank used by standard Omega. See text for a detailed discussion. All template banks are generated using the same maximum mismatch µ max = 0.2. components, one due to the finite size of the chirplet ∆ finite size f = 1/∆ t and the other due to its sweeping frequency ∆ chirp f = d∆ t , where ∆ t is the chirplet duration (defined above). According to [START_REF] Chatterji | The search for gravitational-wave bursts in data from the second LIGO science run[END_REF], the width of the bi-square window should be set to the chirplet frequency bandwidth ∆ f rescaled by a factor of √ 11.

) 2 +(∆ chirp f ) 2 of two

Pre-and post-processing

In this paper we focus on a single-detector network, where most of the pre-and post-processing can be adopted from the standard Omega pipeline. Pre-processing consists of whitening the input data stream. Post-processing consists of selecting among the chirplets with partial time and frequency overlap to the one with maximum correlation with the data. Each chirplet is associated with a time-frequency tile, defined by [t ± ∆ t /2, f ± ∆ f /2] where ∆ t and ∆ f are the chirplet duration and bandwidth, respectively. Two chirplets overlap if their time-frequency tiles overlap.

Performances of Chirpletized Omega

In this section we present a comparison between the standard version of the Omega pipeline, which uses sine-Gaussian wavelets, and its chirpletized version. We configure the pipelines with identical values for the parameters they have in common (frequency and Q range and maximum mismatch). We identify cases where we can expect advantages from analyzing the data with chirplets.

Analyzing chirplets with sine-Gaussians

The signal space associated with sine-Gaussians is contained in the larger space associated with chirplets. We estimate the signal-to-noise ratio (SNR) loss occurring when analyzing a chirplet by correlating this signal against a sine-Gaussian template bank. The chirplet parameters have been set to Q = 50, f = 256 Hz and d = 2048 Hz/s. Those parameters correspond to observable physical signals in the LIGO/Virgo frequency band (for instance, the selected chirping rate is approximately that of an inspiralling binary chirp with total mass M ∼ 3M -assuming equal masses -at f = 256 Hz according to the Newtonian model). Fig. 4 presents the result of this analysis. Consistently to the metric estimate, the loss is √ 128πf 2 /(dQ 2 ) ∼ 50% in the present case. Note also that the maximum correlation is shifted to lower Q which may lead to a possible bias in the estimation of this parameter.

Analyzing inspiralling binary chirps with chirplets

One case study -As an illustration, we show here results from chirpletized Omega on simulated Gaussian noise, colored with the spectral characteristics of LIGO/Virgo noise with simulated gravitational wave signals. The signal we consider here results from the phenomenological approximation introduced in of the coalescing binary black-hole chirp signals and includes the inspiral, merger and ringdown parts of the coalescence.

In Fig. 5, we compare the results from the standard and chirpletized Omega pipelines obtained for a black-hole binary chirp embedded in simulated Gaussian noise at large SNR. Chirplets with a positive slope are preferred to sine-Gaussian with constant frequency: the correlation of the most significative chirplet is, in this example, ∼ 45% larger than the most significative sine-Gaussian. Work is currently in progress to understand how the background in chirpletized Omega is different from standard Omega. Preliminary studies in Gaussian noise suggest that the background rates are comparable, so that we can expect an increase of ∼ 30-40% in distance reach by using chirplets.

Note that the chirplet slope provides indication of the frequency evolution of the observed signal and thus may be very useful in the a posteriori interpretation of an event.

Systematic study -We also performed a more systematic comparison over a population of inspiralling binaries. We considered a total of 5500 binaries with equal mass components. The total mass M is extracted from a flat distribution in the 4-100 M range. The signal amplitudes are scaled so that the SNR is distributed over an interval ranging from ∼ 10 to ∼ 10 3 . We obtain an estimate of the injected SNR from the amplitude of the most significant template. In the ideal case where signal and template are identical, the estimate equals the injected value. In Fig. 6, we show the relative difference of the SNR estimated by chirpletized Omega and standard Omega. Chirpletized Omega estimates a higher SNR across the mass range. However, there are two regimes: for the high-mass range M 60M , the SNR improvement is small (∼ 5%) while it is more pronounced (∼ 20%) for the low-mass range M 60M . The improvement may go upto ∼ 40% for M 20M .

Generally speaking, the spectrum of the GW chirp is shifted toward low frequencies when the binary mass increases. The frequency associated with the innermost stable circular orbit (ISCO), which corresponds to the transition between inspiral and merger phase of the coalescence, is below 70 Hz for masses M 60M . In this condition, the chirp phase of the waveform with a spectral content at frequencies below ISCO is outside the detector sensitive band4 and thus does not contribute significantly to the SNR. Sine-Gaussian waveform provides a good enough fit of the remaining few waveform cycles associated with the merger and ringdown parts of the coalescence. This explains the two regimes in Fig. 6.

Status and future plans

We introduced a new chirplet-based extension of the Omega pipeline. We show preliminary results using coalescing binary mergers waveforms. Versatility (robustness to signal model uncertainty) and algorithmic simplicity are two advantages of this methodology when compared to more standard approaches for the detection of such waveforms.

The single-detector network search code is ready and it can be downloaded [1] and used to produce chirpletized Omega scans similar to the one we show in Fig. 5.

We continue to study the response of the code to real noise and we aim at a complete, operating pipeline using chirplets as templates, and new clustering strategies tailored to these templates, as well as a multi-detector network strategy.

Figure 1 .

 1 Figure 1. Example of a chirplet

Figure 2 .

 2 Figure 2. Example of a chirplet family resulting from the template placement procedure presented in Sec. 3. In this graph, each box represents the time-frequency tile associated with a chirplet. The oblique boxes are associated to non-zero values of the chirping rate d. The slope of the upper and lower edges equals d.

Figure 3 .

 3 Figure 3. Size of the chirplet template bank in two cases: (left) assuming chirp rate limits between ±d max for any frequencies; (right) assuming frequency dependent limits consistent with the Newtonian model of the inspiralling binary chirp: CM 5/3 min f 11/3 d CM 5/3 max f 11/3 . The results are compared to the size of the sine-Gaussian template bank used by standard Omega. See text for a detailed discussion. All template banks are generated using the same maximum mismatch µ max = 0.2.

Figure 4 .

 4 Figure 4. Correlation measurement between a single chirplet and a sine-Gaussian template bank. The parameters of the analyzed chirplet is indicated with a black cross in this diagram.

Figure 5 .

 5 Figure 5. (top/left) Inspiralling black-hole binary (with masses m 1 = 14M , m 2 = 16M and non-precessing spin parameters χ 1 = -0.68 and χ 2 = -0.48) signal in simulated Gaussian LIGO/Virgo-like noise. (top/right) Spectrogram (bottom/left) Significant time-frequency tiles for standard Omega (using sine-Gaussian wavelets only) (bottom/right ) Significant timefrequency tiles for chirpletized Omega (using chirplets).

Figure 6 .

 6 Figure 6. SNR enhancement from standard Omega to chirpletized Omega. The largest improvements about ∼ 30 to 40% are mainly due to binaries in the low mass range M < 20M .

This calculation assumes that the detector noise has a flat spectrum. Contrarily to the sine-Gaussian case, this approximation has significant effect since the chirplet frequency varies across the detector bandwidth.

This result is valid both when neglecting or retaining the off-diagonal terms of the metric. The scaling is actually exact when the off-diagonal terms are included and valid to a good approximation when Qmax 2 in the other case.

This statement is valid for the LIGO detector noise curve which was used for the systematic study. The mass cut-off is higher for the Virgo detector as its sensitive band extends to lower frequencies.

Acknowlegments

We are grateful to the LIGO Scientific Collaboration for the use of its algorithm library available at https://www.lsc-group.phys.uwm.edu/daswg/projects/lalapps.html. This work is supported in part by NSF grant PHY-0653550. M. Miele is supported by an ERASMUS fellowhip from the University of Sannio at Benevento (Italy).