text:
# 1 Annotation

## 1 Annotation

The analysis results of the isotopic yields for products obtained in reactions of protons at energies 0.66, 1, 8.1 GeV with separated tin isotopes are presented. The experimental data are compared in form of reactions of element cross sections for two targets of different isotope compositions. The scaling behaviour have been shown in the dependence of this ratios from neutron and proton composition at the reaction products. The discussion of scaling parameters allows to explain the mechanism of light and medium mass fragment production.

## 2 Introduction

Our representation about of mechanism of nuclear interactions constrained with the characteristics of the emitting products in the nuclear collision. The complex fragments were generated in the nuclear reactions of intermediate and hight energies related to the various processes dynamics of those don’t described exactly.

Recently some scientific nuclear centers have been carrying out research with ion beams to study the effect of the isotopic composition of fragments on their production probability. Note that despite a lot of measurements of multiplicities of light and intermediate-mass fragments and their characteristics, little is known about the effect of target and projectile nucleon composition on the formation probability for various reaction products. Some authors - investigated ratios between cross sections for formation of light fragments in targets with different neutron-to-proton composition:

$`R_{21}(N,Z)=Y_2(N,Z)/Y_1(N,Z)`$ (1)
where $`Y_i(N,Z)`$ is the formation cross section for a fragment with Z protons and N neutrons, $`Y_2(N,Z)`$ is for yields from targets with a greater number of neutrons than $`Y_1(N,Z)`$. Systematization of experimental data  according to isotopic ratios of formation cross sections for light nuclei with Z = 1-8 made it possible to reveal clear dependence on the nucleon composition of products and to parameterize relation (1) in the form of a function of the number of neutrons and protons in the emitted fragments

$`R_{21}(N,Z)=c\mathrm{exp}(\alpha N+\beta Z)`$ (2)
where $`c`$ is the normalizing factor, $`\alpha `$ and $`\beta `$ are the empirical parameters. In , the use of three parameters to describe this dependence for a variety of reactions was considered as an isotopic scaling.

In the ideal experiment information on primary interaction can be obtained from the distributions of prefragments resulting from nuclear collision of about $``$ $`10^{21}`$ s. Yet, the real time scale of nuclear fragment measurement is $`10^{10}`$ s. The analysis , with various model representations showed that when considering isotopic relation (1), one can ignore the effect of the slow nucleon emission step for various types of reaction (multifragmentation, evaporation, deep spallation) provided that formation of the products in question is of statistical character and the temperature of excited nuclear systems is approximately identical.

In the model considered, the coefficients $`\alpha `$ and $`\beta `$, which define factorization of (2), characterize dependence of the measured distribution on the hot emitter properties, including proton and neutron density in various transition phases, on nuclear matter asymmetry allowed for by entering a corresponding term into the equation of nuclear state, on difference in nucleon evaporation conditions, etc.

Thus, the method of isotopic relations allows studying the early reaction step of the decay composite fragmenting system reaction.

Most of the data analyzed in , relate to nucleus-nucleus processes, which are characterized by considerable effects resulting from strong deformation of nuclear matter in collision, such as compression, rotation, etc. In this sense the study of reaction induces by light nuclei and particles makes it possible to isolate the thermal component in nuclear excitation.

According to , the multifragmentation process suggests decay of a hot nuclear system with formation of light and residual nuclei with $`3Z30`$ . Light nuclei can be measured and identified and their isotopic properties and formation channels can be established only with modern detecting facilities including telescopic systems of high-resolution counters in the 4$`\pi `$ geometry. Identification of heavy product nuclei is a more difficult problem. Activation analysis is one of the methods allowing the nucleon composition of heavy products to be determined exactly.

The purpose of the presented work was to analyze the ratios of product yields from reactions induced by 0.66 and 8.1-GeV protons in targets of separated tin isotopes by the induced activity method.

Our earlier investigations of the relationship between the formation cross sections for residual nucleus and the nucleon composition of the residual and target nuclei revealed correlation between the yields of these nuclei and their third projection of the isotopic spin -. Therefore, the isotopic ratio of the product formation cross sections can be considered in the form

$`R_{21}(N,Z)=Y_2(N,Z)/Y_1(N,Z)=\mathrm{exp}(C+(NZ)/2)`$ (3)
where $`C`$ and $`B`$ are the fitting parameters and $`(NZ)/2=t_3`$ is the third projection of the product isotopic spin. This expression is in general similar to formula (2) in , but the yield ratios were systematized according to the isospin of the resulting fragment.

During the element-by-element analysis of formation cross sections for various isotopic states dependence of the ratio $`R_{21}`$ on the number of protons in the product is excluded from (3) and the neutron content effect defined by the parameter $`B/2`$ corresponds to the parameter $`\alpha `$ in (2) . According to the theoretical evaluation ,, , the value of $`\alpha `$ depends on a few factors, including the neutron chemical potential difference $`\mathrm{\Delta }\mu _n`$ in both targets, and on the value of the asymmetrical energy term appearing in the equation of nuclear state .

## 3 Reactions at 8.1 GeV

The analysis of experimental data from reactions obtained by protons at 8.1 GeV energy on different tin isotopes are presented below. The yields of the reaction products in range $`19Z37`$ considerably increased at this energy -, which agrees in general with the predictions of the statistical multifragmentation model (SMM) . A distinctive feature of our measurements  was a possibility of investigating dependence of reaction product yields on the neutron content of the targets used ($`{}_{}{}^{124}Sn,^{120}Sn,^{118}Sn,^{112}Sn`$ ).

We investigated the dependence of $`R_{21}(N,Z)`$ for the measured reaction products on the isotopic spin in various target pairs. Approximating these data in the form of the exponential dependence (3), we got a set of parameters $`B`$ characterizing dependence of the ratio $`R_{21}(t3)`$ both on the isotopic spin of the product and on the difference in neutron content of the targets ($`\mathrm{\Delta }N=12,8,6,4`$). The results of calculations for some reaction products are presented in Table 1; the dependence in question for $`Br`$ isotopes is displayed in Fig. 1. (For clarity, the calculated curves are shown on the relative scale. The factors introduced for convenience are given in the figure captions.) It is evident that within the calculation error all $`R_{21}(t3)`$ are characterized by the identical slope to the $`(NZ)/2`$ axis for the targets with a certain neutron excess and that the dependence curves become steeper as the difference in the neutron content of targets increases.

In Fig 2 dependence of the parameter $`B`$ on the relative neutron excess of targets is plotted for the target pairs $`{}_{}{}^{118}Sn/^{112}Sn,^{120}Sn/^{118}Sn,^{124}Sn/^{118}Sn,^{124}Sn/^{112}Sn`$ (solid curve). The relative neutron excess in targets ($`\mathrm{\Delta }N`$) was used as variable. The parameter is seen to vary practically linearly with increasing neutron content difference, the curve slope being 0.11 $`\pm `$ 0.02.

It is evident from Table 1 that the average value of the parameter $`B`$ is 1.1$`\pm `$0.22 ($`\mathrm{\Delta }N=12`$), which corresponds to $`\alpha `$=0.58$`\pm `$0.06 . This value ($`\alpha =0.6`$)was used in  to approximate data for evaporation products.

The calculations by the SMM  used by us to describe this type of reactions reflect to an extent the behavior of isospin characteristics but no full agreement is found on closer examination (dashed line in Fig. 1).

Cross sections for formation of light fragments in the mass number region $`7A39`$ obtained by us at the proton energy 8.1 GeV  were also analyzed by formula (3). As there were not enough data for the element-by-element analysis, we used all yields of light products in the above mass region for approximation. The results of the calculations are presented in Table 1 and displayed in Fig. 3. The parameter B is seen to vary with the isotopic spin of both the product and the target. When isoscaling in reactions of this type is discussed , , the parameter value $`\alpha `$=0.37 ($`\mathrm{\Delta }N=12`$) is used to describe multifragmentation processes. We obtained the value $`B/2=\alpha 0.28\pm 0.02`$ for $`{}_{}{}^{124}Sn/^{112}Sn`$ (Fig. 3), which points to multifragmentation rather than evaporation mechanism for formation of these nuclei at the proton energy 8.1 GeV. Calculations by the SMM  showed that the isotopic dependence in general manifests itself in relative cross sections for formation of light fragments but particular comparison does not show satisfactory agreement.

Dependence of the parameter B on the neutron excess in the targets is shown by a curve with a slope of 0.049$`\pm `$0.01 (dashed curve in Fig. 3).

## 4 Reactions at 0.66 GeV

The experimental data in this energy region are not enough to carry out a thorough element-by-element analysis, as in the case of 8.1 GeV. Figure 4 (Table 2) present some cross section ratios from which it is evident that the isospin dependence takes place at this energy also, including the effect of target neutron excess. If we take into account closeness of the charge states of the detected products ($`Z=4045`$) and build an averaged dependence curve (Fig. 5), all values for $`\mathrm{\Delta }N=12`$ and $`\mathrm{\Delta }N=6`$ can be placed within the regions characterized by the slope parameters $`\alpha =0.48\pm 0.04`$ and $`\alpha =0.295\pm 0.015`$ respectively. In general, the parameter values yielded by our analysis agree within the determination error with the data used in  to explain formation of reaction products through the evaporation. It is seen that for $`\mathrm{\Delta }N=12`$ the effect of the neutron excess in targets also manifests itself. The calculations by the cascade model  also yield the isospin dependence of reaction cross sections, but the values of the coefficients in the approximating functions turned out to be larger than in the above analysis (dashed line in Fig. 4).

To get more information on generalities of the isospin-cross section relationship in nuclear reactions, we carried out a similar analysis using the experimental data on formation of light fragments from targets of different isotopic composition in reactions induced by 0.66 and 1-GeV protons , . The cited papers point to distinct dependence of the formation probability for various isotopes of light nuclei on the isotopic state of target nuclei. Analysis by the scheme proposed in  yielded that $`R_{21}=Y_2(t_3)/Y_1(t_3)`$ for the $`H,He,Li,Be,B`$, and $`C`$ produced in $`{}_{}{}^{124}Sn`$ and $`{}_{}{}^{112}Sn`$ targets ($`\mathrm{\Delta }N=12`$) was satisfactorily defined by the expression $`\mathrm{exp}(0.6N0.82Z)`$ within the 15-20% accuracy. The same character of dependence and coefficient values were taken for description of evaporation processes in . Thus, we may assume that fragmentation in proton-nucleus interactions at the energy of 1 GeV and below has predominantly of the evaporative character.

## 5 Conclusion

Analysis of the experimental data on nuclear reactions in tin targets of various isotopic composition ($`1.24N/Z1.48`$ ) in the region of intermediate and high energies showed distinct isotopic dependence of product formation cross sections. The systematics proposed in - allows the assumption that this dependence is determined by the initial step of the interaction and the observed isoscaling is due to the statistical nature of fragment formation. This dependence is reflected in general in the statistical models used for description of fragments. Analysis of the coefficients in the expression for the scaling function on the basis of the data -, ,, showed that yields of light fragments at the proton energy 8.1 GeV corresponded to the multifragmentation process ($`\alpha =0.28`$ ) and at 1 GeV and below ($`\alpha =0.6`$ ) to the evaporation process. Yields of medium-heavy and heavy products are due to evaporation and multifragmentation. Measurements of energy distributions and other experimental data corroborate the presence of two sources of their formation. Observation of the dependence in question in a wide range of energies and reaction types indicates that a universal property of nuclear interactions manifests itself.

| Residual | Values of the fitting parameter B($`E_p=8.1`$) | | | | |
| --- | --- | --- | --- | --- | --- |
| nuclei | $`\mathrm{\Delta }N=4`$ | $`\mathrm{\Delta }N=6`$ | $`\mathrm{\Delta }N=8`$ | $`\mathrm{\Delta }N=12`$ | $`\alpha _{\mathrm{\Delta }N=12}`$ |
| $`{}_{}{}^{4347}Sc_{21}`$ | 0.33$`\pm `$0.13 | 0.53$`\pm `$0.08 | 0.81$`\pm `$0.16 | 1.23$`\pm `$0.20 | 0.61$`\pm `$0.10 |
| $`{}_{}{}^{7074}As_{33}`$ | 0.36$`\pm `$0.01 | 0.52$`\pm `$0.01 | - | 1.10$`\pm `$0.22 | 0.55$`\pm `$0.11 |
| $`{}_{}{}^{7477}Br_{35}`$ | 0.34$`\pm `$0.08 | 0.63$`\pm `$0.10 | 0.81$`\pm `$0.04 | 1.16$`\pm `$0.05 | 0.58$`\pm `$0.05 |
| $`{}_{}{}^{7679}Kr_{36}`$ | - | - | - | 1.10$`\pm `$0.03 | 0.55$`\pm `$0.17 |
| $`{}_{}{}^{8184}Rb_{37}`$ | 0.27$`\pm `$0.10 | 0.61$`\pm `$0.06 | 0.94$`\pm `$0.10 | 1.23$`\pm `$0.13 | 0.61$`\pm `$0.06 |
| $`{}_{}{}^{9396}Tc_{43}`$ | 0.28$`\pm `$0.12 | - | - | 1.22$`\pm `$0.12 | 0.61$`\pm `$0.06 |
| $`{}_{}{}^{94103}Ru_{44}`$ | - | 0.56$`\pm `$0.07 | 0.69$`\pm `$0.04 | 1.18$`\pm `$0.08 | 0.59$`\pm `$0.04 |
| $`{}_{}{}^{103113}Ag_{47}`$ | 0.37$`\pm `$0.03 | - | - | - | - |
| the mean | - | - | - | B=1.17$`\pm `$0.36 | $`\alpha `$=0.58$`\pm `$0.06 |
| Light fragments ($`4Z17`$) | | | | | |
| $`\mathrm{\Delta }N=2`$ | $`\mathrm{\Delta }N=4`$ | $`\mathrm{\Delta }N=6`$ | $`\mathrm{\Delta }N=8`$ | $`\mathrm{\Delta }N=12`$ | $`\alpha _{\mathrm{\Delta }N=12}`$ |
| 0.11$`\pm `$0.02 | 0.13$`\pm `$0.05 | 0.27$`\pm `$0.10 | 0.44$`\pm `$0.05 | 0.56$`\pm `$0.02 | 0.28$`\pm `$0.02 |

| Residual | Values of the fitting parameter B ($`E_p=0.66`$) | | |
| --- | --- | --- | --- |
| nuclei | $`\mathrm{\Delta }N=4`$ | $`\mathrm{\Delta }N=6`$ | $`\mathrm{\Delta }N=12`$ |
| $`{}_{}{}^{9096}Nb_{41}`$ | 0.60$`\pm `$0.10 | - | 1.10$`\pm `$0.14 |
| $`{}_{}{}^{9196}Tc_{43}`$ | - | 0.55$`\pm `$0.11 | - |
| $`{}_{}{}^{94105}Ru_{44}`$ | - | $`0.67\pm `$0.05 | - |
| $`{}_{}{}^{99105}Ru_{45}`$ | - | $`0.67\pm `$0.06 | 1.30$`\pm `$0.07 |
| $`{}_{}{}^{103110}Ag_{47}`$ | - | $`0.60\pm `$0.12 | - |


format: text

text:
# Component separation for Cosmic Microwave Background data: A blind approach based on spectral diversity

## 1 Introduction

The precise measurement of Cosmic Microwave Background (CMB) anisotropies is one of the main objectives of modern cosmology. The small temperature fluctuations of the CMB with respect to the pointing direction on the sky reflect the primordial density perturbations in the young Universe. The spatial power spectrum of those fluctuations depend on a set of important parameters, known as cosmological parameters. The accurate estimation of CMB power spectrum is thus of prime importance in cosmology. The Planck satellite (to be launched in 2007) will map whole sky emission with unprecedented signal to noise ratio.
The accuracy of CMB power spectrum measurements are limited by the other astrophysical emissions present in the sub-millimeter range of the spectrum. Those emissions, called foregrounds (as they are emitted in front of the CMB), are of different origin. Some of them originate from within our own Galaxy, as the dust and synchrotron emissions, others are of extra-galactic origin as the Sunyaev-Zel’dovich effects. All depend on the wavelength.

It is thus important for a precise measurement of the CMB power spectrum to deal with the various foregrounds present with the CMB anisotropies. The availability of several detectors operating in several bands (10 for Planck ranging from 30 to 850 GHz) allows to distinguish the various contributions. Component separation methods has been adressed by a number of authors . The standard approach consists in producing the cleanest maps of the various components, followed by estimation of the power spectrum of the CMB from the separated CMB map. This approach is not fully satisfactory for two reasons. First, component separation requires prior knowledge of the electromagnetic spectra of the components which are not all very well known. Second, it would be preferable to estimate the CMB power spectrum in one step by jointly analysing the observation maps.
A new approach to process multi-detector multi-component data has been developed in papers  and . The method is based on likelihood maximization in the Whittle approximation and takes additive noise into account (in the CMB observations, a significant amount of noise is expected). The spectral diversity of the various components is used. The method has been implemented on the domain of small sky maps, in the flat sky approximation.
In this paper, we have refined the multi-detector multi-component (MDMC) spectral matching method in several aspects to account for distinctive CMB observation features. First, as current observations cover as large a portion of the sky as possible, we have adapted the method for a spherical harmonic expansion of all sky maps in order to exploit all the information contained in the maps. Second, we account for the finite spatial resolution of the detectors, to take into account the fact that the sky is not seen with the same resolution by all detectors. Finally, our method allows for the inclusion of some physical knowledge about the components. Fixing some parameters, the likelihood is maximized over the other parameters. This operation allows, for example, to break degeneracies due to components having similar spatial power spectra.
We have implemented the method on full sky simulated Planck observations. We compare the estimated CMB power spectra in the blind approach and in the ideal case where one has a perfect knowledge of the component emission laws.

## 2 Model of sky emission

The key assumption is that the sky emission at a given frequency is the linear superposition of astrophysical components. In addition, we assume that the emission laws of the various components do not depend on the position on the sky. The signal measured by a detector is the sky emission convolved by a beam shape (depending on the detector), plus an additive noise. Assuming a symmetric beam shape, the observation by a detector $`d`$ is given by :

$$x_d(\stackrel{}{r})=\underset{c=1}{\overset{N_c}{}}A_{dc}.b_d(|\stackrel{}{r}\stackrel{}{r^{}}|).s_c(\stackrel{}{r^{}})d\stackrel{}{r^{}}+n_d(\stackrel{}{r})$$
(1)
$`\stackrel{}{r}`$ is the direction of observation on the sky, $`s_c`$ is the emission template for source $`c`$, $`n_d`$ represents the noise of the detector $`d`$, $`b_d(|\stackrel{}{r}\stackrel{}{r^{}}|)`$ represents the beam, and depends only on $`|\stackrel{}{r}\stackrel{}{r^{}}|`$ (or on $`\stackrel{}{r}.\stackrel{}{r^{}}`$) for symmetric beams and $`A`$ is the mixing matrix. Each element of the mixing matrix results from the integration of the emission law of one component over one detector frequency band.

A natural basis for the application of the MDMC spectral matching method, exploiting the components spectral diversity is the spherical harmonics basis. We write the decomposition of the signal in this basis :

$$x_d(\stackrel{}{r})=\underset{l=0}{\overset{\mathrm{}}{}}\underset{m=l}{\overset{l}{}}x_d(l,m)Y_l^m(\stackrel{}{r})$$
(2)
where $`Y_l^m(\stackrel{}{r})`$ are the spherical harmonic functions. Coefficients $`x_d(l,m)`$ can be easily computed using the orthogonality of spherical harmonic functions :

$$x_d(l,m)=_{4\pi }x_d(\stackrel{}{r})Y_l^m(\stackrel{}{r})^{}𝑑\mathrm{\Omega }$$
(3)
The convolution between a symmetric beam and the signal in real space becomes a product in the spherical harmonic basis. Then, we obtain the observation coefficients combining equations 1 and 3

$$x_d(l,m)=b_d(l)\underset{c=1}{\overset{N_c}{}}A_{dc}s_c(l,m)+n_d(l,m)$$
(4)
where $`b_d(l)`$ is the Legendre Polynomial expansion of
$`\mathrm{b}_\mathrm{d}(\mathrm{cos}\mathrm{\Theta })`$, $`\mathrm{\Theta }`$ being the angular distance from the center of the beam so that $`\mathrm{cos}\mathrm{\Theta }=\stackrel{}{r}.\stackrel{}{r^{}}`$. For a Gaussian beam, $`b(l)\mathrm{exp}(\sigma _b^2l(l+1)/2)`$ and $`\sigma _b=\mathrm{\Theta }_{\mathrm{beam}}/\sqrt{8\mathrm{l}\mathrm{n}2}`$, where $`\mathrm{\Theta }_{\mathrm{beam}}`$ is the full width at half maximum of the beam.
Let us consider the diagonal matrix $`B(l)`$ such that the diagonal element $`B_{dd}(l)=b_d(l)`$. We define the new coefficients $`x^{}(l,m)=B(l)^1x(l,m)`$. It is useful to write these “deconvolved” observations coefficients in a matrix form :

$$x^{}(l,m)=As(l,m)+B(l)^1n(l,m)$$
(5)
The introduction of $`x^{}`$ will be justified later.

Spectral statistics

We now need to compute the spectral statistics of the observations $`x^{}(l,m)`$, given by $`R_x^{}(l,m)=x^{}(l,m)x^{}(l,m){}_{}{}^{}`$, where $`^{}`$ denotes transpose-conjugation.

$$R_x^{}(l,m)=AC(l,m)A^{}+B(l)^1N(l,m)B(l)^{}$$
(6)
$`C(l,m)`$ and $`N(l,m)`$ are respectively the component and the noise covariance matrices. Statistical independence between components implies that $`C_{lm}`$ is a diagonal matrix. We also assume that the noise is white and independent across detectors, so that : $`N(l,m)=diag(\sigma _1^2,\mathrm{},\sigma _d^2)`$. The main reason to work with variable $`x^{}`$ is that we can average the covariance matrices $`R_x^{}(l,m)`$ over bins and preserve the simple structure of equation (6) for the signal part. Let us define the following particular averaging over bins:

$$R_x^{}(q)=\frac{1}{n_q}\underset{l=l_{\mathrm{min}}(q)}{\overset{l_{\mathrm{max}}(q)}{}}\underset{m=l}{\overset{l}{}}R_x^{}(l,m)$$
(7)
Here $`q=1,Q`$ is the spectral bin index, the modes $`(l,m)`$ contributing to $`q`$ are such that $`l_{\mathrm{min}}(q)<l<l_{\mathrm{max}}(q)`$. Since $`m`$ can vary between $`l`$ and $`l`$, the number of modes in each bin $`q`$ is $`n_q=_{l=l_{\mathrm{min}}(q)}^{l_{\mathrm{max}}(q)}(2l+1)`$. The reason for the choice of such averaging is that the spectral covariance matrices of isotropic components on the sky (we have strong reasons to think that cosmological components are isotropic) do not depend on the parameter $`m`$, so $`C(l,m)=C(l)`$. The average covariance matrices are :

$$R_x^{}(q)=AC(q)A^{}+M(q)$$
(8)
where $`M(q)=1/n_q_{l=l_{\mathrm{min}}(q)}^{l_{\mathrm{max}}(q)}(2l+1)B(l)^2N`$ is a diagonal matrix.

The band-averaged spectral covariance matrices are estimated by :

$$\stackrel{~}{R}_x^{}(q)=\frac{1}{n_q}\underset{l=l_{\mathrm{min}}(q)}{\overset{l_{\mathrm{max}}(q)}{}}\underset{m=l}{\overset{l}{}}x^{}(l,m)x^{}(l,m)^{}$$
(9)
which is real valued since $`x^{}(l,m)=x^{}(l,m)^{}`$ for real data.

## 3 The MDMC method

The aim of the MDMC spectral matching method is to obtain an estimate of different parameters of the model, of relevance in astrophysics and cosmology, without the need of prior information. Those parameters are the mixing matrix $`A`$, the band averaged component power spectra $`C_q`$ and the noise covariance $`N`$. They are collectively referred to as
$`\theta =\{A,C(q),N\}`$.
The method is based on minimazing the mismatch between the empirical spectral covariance matrices of the observations (equation 9) and their expected values which depends on the mixing parameters (equation 8). The mismatch is quantified by the average divergence measure between two matrices :

$$\mathrm{\Phi }(\theta )=\underset{q=1}{\overset{Q}{}}n_qD(\stackrel{~}{R}_x^{}(q),R_x^{}(q))$$
(10)
where $`D(R_1,R_2)`$ is a mesure of the divergence (the Kullback divergence) between two positive matrices and $`n_q`$ is the number of modes in each band $`q`$. Assuming that the spherical harmonics coefficients of the components are random realizations of Gaussian field of variance $`R_x^{}(q)`$ and are uncorrelated, the log-likelihood (up to irrelevant factor) take the same form as in equation 10 (in the frame of the Whittle approximation). The divergence is given in this case by :

$$D(R_1,R_2)=\mathrm{tr}(R_1R_2^1)\mathrm{log}\mathrm{det}(R_1R_2^1)m_d$$
(11)
The parameter estimate is given by $`\widehat{\theta }=\mathrm{argmin}_{/\theta }\mathrm{\Phi }(\theta )`$. The connection with the log-likelihood guarantees good statistical properties of the estimates (at least asymptotically).

Optimization algorithm:

The optimization is made using an Expectation-Maximization (EM) algorithm, and completed by quasi-Newton algorithm (BFGS) in order to accelerate the convergence. The formal algorithm described in  is used, with minor changes (the equations change slighly when we account for the beam smoothing effect).

Degeneracies

As seen from equation 8, one can exchange a factor between the mixing matrix $`A`$ and the component power spectra $`C(q)`$ without changing the result of the equation. We fix this degeneracy in the EM step by fixing the norm of each column of $`A`$ to unity; the power spectra are adjusted accordingly. In the quasi-Newton step, a penalty term is added to $`\mathrm{\Phi }`$ which penalizes the deviation of each column of $`A`$ with a minimum penalization at unity.

## 4 Application

We now turn to the application of the MDMC spectral matching method on synthetic data.

### 4.1 Simulated Planck data

We use full sky Planck simulations provided by the Planck consortium. The maps are generated at the ten frequencies of Planck instruments (30, 44, 70, 100 GHZ for the low frequency instrument and 100, 143, 217, 353, 545, 857 GHz for the high frequency instrument). Five components and white noise at the nominal level of Planck instruments are mixed according to equation 4 (the beam sizes are by increasing order in frequency : 33, 23, 14, 10, 10.6, 7.4, 4.9, 4.5, 4.5, 4.5 arcminutes). The components are the CMB, the kinetic and thermal SZ effect from galaxy clusters, and the galactic dust and synchrotron. They are obtained as follows:
The CMB component is randomly generated using a power spectrum $`C_l`$ predicted by CMBFast  with standard cosmological parameters. The galactic components are generated using observation templates at different frequencies, the galactic dust is modeled using the DIRBE-IRAS 100$`\mu \mathrm{m}`$ maps and the galactic synchrotron is a destriped version of the 408 MHz Haslam survey with additional small scale structures (see ). The Thermal and the Kinetic templates are simulated with a gas dynamics code . Note that the kinetic SZ effect (always sub-dominant) and the CMB component have similar emission laws. The simulations are performed up to the resolution of 3.5 arcminutes. Figure 3 shows the simulated Planck observation maps at all the frequencies.

### 4.2 Estimated parameters

We have applied our method on the above synthetic data. It is necessary to fix the number of components assumed to be present in the data. We choose to characterize 4 components. One of the components in the simulations, the kinetic SZ, is negligible at all frequencies. Moreover, it can not be separated from the CMB directly with our approach as the consequence of their proportional emission law (CMB and kinetic SZ form one component).
Spherical harmonics are computed up to multipole $`l`$=3000 (corresponding to 3.5 arcminutes), and we choose bins of width $`\mathrm{\Delta }l=10`$. We investigate two different approaches :
$``$ First we adopt a quasi-blind approach. We estimate all the elements of the mixing matrix except three elements, corresponding to three of the four elements at 857GHz, which are fixed to zero (we expect the presence of only one component at this frequency). We will discuss the reason for this choice in subsection 4.4. All the other parameters are estimated, including component power spectra and the noise variances. The total number of parameters is $`10\times 43+4\times 300+10=1247`$, compared to $`300\times 10\times (10+1)/2=16500`$ data elements
$``$ In a second approach, we assume that we know the electromagnetic spectra of the components perfectly, so we fix all the elements of the mixing matrix to their true values. We estimate the component spatial power spectra and the noise variances.

### 4.3 Results

In the blind case, all the parameters of the mixing matrix corresponding to the CMB and dust components are estimated with a very good accuracy. Table 1 gives the ratio between the recovered and the true mixing elements of the CMB. The mixing matrix elements corresponding to the thermal SZ are estimated with a good precision. The galactic synchrotron emission law is very well constrained at lower frequencies. Therefore, we show that strong constraints can be put with our method on the component emission law, in particular for the galactic components at frequencies far from their maximum emission.

We now concentrate on the spatial power spectrum of the CMB component. Figure 1 shows the estimated CMB power spectrum in the blind approach. Figure 2 shows the relative errors made on the CMB power spectrum estimation given by ($`|\stackrel{~}{C}(q)C(q)|/C(q)`$) for the two approaches.

The remarkable point is that the estimation errors in both approaches are equivalent. Therefore, it seems that the emission intensity of the components at the different observation frequencies can be estimated in Planck data without loss of accuracy in the linear CMB power spectrum estimation.
In both cases, the method allows to estimate accurately the CMB power spectrum up to $`l2000`$. At smaller scales the dispersion begins to be significant. This result is not surprising since the noise and the beam smoothing effect are important at these scales for all detectors. Also, the estimated power spectrum does not seem to be contaminated by the kinetic SZ.

### 4.4 Discussion

Semi-blind approach?

As we have seen before, in the case where we estimate the mixing matrix as well as the power spectra of the components, we have fixed to zero the contribution of three components at 857GHz, the only component we assume to be present at this frequency is the galactic dust. This hypothesis, being very realistic, appears to be necessary because we expect that the galactic dust and the galactic synchrotron have quasi proportional power spectra. Without any prior, the method can not separate these two components, but the previous simple procedure of fixing some parameters allows to break the degeneracies and to estimate accurately all the components.

Partial coverage

The galactic components are very inhomogeneous on the sky. In the region of the galactic plane, their intensities are several hundred times stronger than in the outer regions. Since we model the components as homogeneous, the method, when applied on those data, is sub-optimal. However, more accurate CMB power spectrum estimation may be obtained by making a galactic cut, involving a partially covered sky processing. The coefficients obtained from the spherical harmonic decomposition of partially covered sky are correlated. Our spectral matching method can be used on such data and a good choice would be to take bin size $`\mathrm{\Delta }l`$ larger than the correlation length in the power spectrum measurement. The method, in particular, is applicable to data such that from the Archeops balloon-borne experiment in which the sky coverage is about $`30\%`$ of the whole sky.

## 5 Conclusion

We have adapted our blind MDMC spectral matching method for the processing of all sky CMB maps. The observations are modeled as a noisy linear mixture of beam-convolved components. By maximizing the likelihood of the system, we estimate the power spectra of the components, their contribution levels at each frequency as well as the noise levels. The method has been applied on full sky Planck simulations containing five components. The spherical harmonics basis was used.
The power spectrum of the CMB is accurately estimated up to $`l2000`$ in bins of size $`\mathrm{\Delta }l=10`$. The comparison between the results obtained in blind and semi-blind approaches shows that the mixing elements can be estimated without loss of accuracy in the CMB power spectrum estimation.

## 6 Acknowledgments

We thank the Planck collaboration and in particular V. Stolyarov and R. Kneissl for the full sky simulated maps. Guillaume Patanchon would like to thanks M.A.J. Ashdown for useful discussions while he was visiting Cavendish Laboratory. The HEALPix package  (see http://www.eso.org/ science/healpix/) was used for the spherical harmonics decomposition of the input maps.


text:
# Evidence for anisotropy in the distribution of short-lived gamma-ray bursts

## 1 Introduction

The Burst And Transient Source Experiment (BATSE) on board the Compton Gamma Ray Observatory collected more than 2000 Gamma-Ray Bursts (GRBs) in its 9-year activity. The analysis of such a large sample has deeply characterized the temporal and spectral properties of GRBs (Fishman & Meegan 1995). Despite the multitude and complexity of the light curves observed, the burst duration, typically defined by $`T_{90}`$ as the time interval during which 90% of the fluence is accumulated (Kouvelitou et al. 1993), indicates the existence of two main classes: short ($`T_{90}<2`$ s) and long ($`T_{90}<2`$ s) bursts. This distinction is also reinforced by the associated spectral shape since short bursts are typically harder than long ones (Fishman & Meegan 1995; Mitrofanov et al. 1998). At present, it is not known if this dichotomy is produced by different progenitors or by different physical mechanisms originating from the same object or, again, from different observing conditions. Nonetheless, it has been suggested (MacFadyen & Woosley 1999) that long bursts might be produced by the core collapse of a single massive stars (the Collapsar model – Woosley 1993; Paczynski 1998) whereas short bursts might be more efficiently generated by double compact object mergers (Goodman 1997; Narayan, Piran & Kumar 2001) or double neutron star heating (Salmonson & Wilson 2002). A distinctive feature of these two scenarios is the expected GRB distance: compact object mergers should have a lower median redshift than Collapsars (Fryer, Woosley & Hartmann 1999) which would be tracers of the high redshift universe (Lamb 2001).

A cosmological origin of long bursts was first suggested by the isotropic sky distribution observed by BATSE (Briggs 1993; Meegan et al. 1992; Paczynski 1995) but the definitive proof came with the first redshift measurement for GRB 970508 (Metzger et al. 1997) at $`z=0.835`$. At present, the distribution of known redshifts extends up to z=4.5 with a mean $`z1.27`$ (Andersen et al. 2000; Djorgovski et al. 2001). However very little is presently known about short GRBs because no radio, optical or X-ray counterpart has been associated with any burst lasting less than 2 s (Hurley et al. 2002), although the possible existence of short burst afterglows has been suggested (Lazzati, Ramirez-Ruiz & Ghisellini 2001). Thus, the distance scale of short GRBs is still an open issue and its solution would have important implications for the understanding of their progenitors.

The aim of this Letter is to investigate the nature of short GRBs by means of statistical techniques such as the angular two-point correlation function, since different hypotheses are expected to leave different imprints on their distribution on the sky. In particular we will concentrate on the two distinct possibilities of burst repetition and large-scale structure.
In order to do so, we have considered sources from the public BATSE on line catalogue111http://gammaray.msfc.nasa.gov/batse/grb/catalog/current/ which contains 2702 bursts from GRB 910421 to GRB 000526. For our analysis, a key ingredient is the estimate of the duration $`T_{90}`$, which has been derived only for 2041 events of this catalogue. Also, the flux and fluence are only available for 2135 events. The maximum number of bursts that we could extract with measured durations and fluences is 2035, of which 497 are short events with $`T_{90}<2`$ s.

Our work is presented as follows: in Section 2 we briefly recall the results of previous works examining evidence for GRB anisotropies, while in Section 3 we introduce the method and show the results on the angular correlation function. Section 4 discusses some possible models for the observed distribution of short GRBs and summarizes our conclusions.

## 2 Evidence of GRB anisotropies

First studies on the angular distribution of GRBs by Hartmann & Blumenthal (1989) including 80 GRBs detected by the IPN and 160 GRBs collected by KONUS, showed a possible deviation from isotropy on scales $`5^{}`$. A more refined analysis (Quashnock & Lamb 1993) based on the first BATSE (1B) catalogue (260 bursts) confirmed this evidence but with a cut on the burst positional errors $`\mathrm{\Delta }<10^{}`$ (but see cautions in Narayan & Piran 1993; Brainerd 1996).

The analysis of joint spatial and temporal clustering of the bursts in the 1B catalog (Wang & Lingenfelter 1995) showed an excess of bursts with angular separation lower than $`4^{}`$ within 4 days with a small probability ($`1\times 10^5`$) of random occurrence. However, less significant results are claimed by Petrosian & Efron (1995). In particular, the need of larger samples for investigation of the possible burst repetition by means of the two–point correlation function is discussed by Brainerd (1996).

Evidence of possible burst clustering was recently claimed for both the populations of short and intermediate GRBs from the analysis of larger samples (Balazs et al. 1998; Meszaros et al. 2000). The possible different nature of the short and long events, in fact suggests the separate analysis and comparison of their clustering properties.

## 3 The Method

One of the possible statistical tests for the detection of anisotropy is provided by the two-point angular correlation function $`w(\theta )`$, which gives the excess probability, with respect to a random Poisson distribution, of finding two sources in the solid angles $`\delta \mathrm{\Omega }_1`$ $`\delta \mathrm{\Omega }_2`$ separated by an angle $`\theta .`$ We applied this technique to the distinct samples of long and short GRBs in order to investigate for possible deviations from isotropy in their spatial distribution.

Given the large uncertainties affecting the determination of the positions of GRBs, especially in the case of short-lived ones (see Figure 1), we have decided to only consider those objects which have angular coordinates with associated errors $`\mathrm{\Delta }10^{}`$. The choice for this value results as the best compromise between the requirements of precision-position measurements and large-enough datasets to allow for sensible statistical analysis. By applying the above cut, we end up respectively with $`N_{\mathrm{long}}=1486`$ (out of 1538) and $`N_{\mathrm{short}}=407`$ (out of 497) sources in the case of long-lived and short-lived GRBs.

Random catalogues have subsequently been generated for each of the two subsets, with positions of the objects uniformly distributed on the celestial sphere. In order to minimize Poisson fluctuations due to the limited number of sources in the original datasets, random catalogues were required to include $`10\times N_{\mathrm{long}}`$ objects in the case of long-lived GRBs and $`20\times N_{\mathrm{short}}`$ objects in the case of short-lived GRBs. We then counted the number of distinct data-data pairs ($`DD`$), data-random pairs ($`DR)`$, and distinct random-random pairs ($`RR`$) as a function of angular separation, and calculated the two-point angular correlation function $`w(\theta )`$ using the estimator (Hamilton 1993)

$`w(\theta )=`$ $`{\displaystyle \frac{4DDRR}{(DR)^2}}`$ $`1`$ (1)
on the angular scales $`1\theta 20`$ degrees.

The open dots in Figs. 2 and 3 illustrate the results respectively for the population of long-lived and short-lived GRBs; the error bars indicate Poisson estimates for the points. We note that, while in the first case values of $`w(\theta )`$ at all angular scales seem to be consistent with an isotropic distribution, the angular correlation function of short-lived GRBs shows a significant excess of power on scales $`24`$ degrees.

In order to test the significance of this result, we have furtherly examined two reduced datasets, obtained from the original one by only considering those sources with positions known to better than $`\mathrm{\Delta }5^{}`$ (1239 long and 196 short GRBs). Solid squares in both Figs. 2 and 3 illustrate the values of $`w(\theta )`$ obtained in this case. Once again, measurements of the angular correlation function for long-lived GRBs are entirely consistent with absence of intrinsic non-randomness, in agreement with earlier results (e.g. Meegan et al. 1995, 1996; Bennett & Rhie 1996; Tegmark et al. 1996) obtained for samples which mainly included this population of gamma-ray emitters, and with more recent findings (e.g. Balazs et al. 1998; Meszaros et al. 2000) concerning the population of long-lived GRBs alone.

On the other hand, the distribution of short-lived, $`\mathrm{\Delta }5^{}`$ GRBs – even though with larger errors due to the small number of objects considered in the analysis – shows a departure from isotropy on small angular scales even stronger than it was for the previously-analysed $`\mathrm{\Delta }10^{}`$ subsample. We can then conclude that the distribution of short-lived gamma-ray bursts exhibits a $``$2$`\sigma `$ deviation from isotropy on scales $`24`$ degrees.

## 4 DISCUSSION

Given the above result on the anisotropic distribution of short-lived GRBs, it is worthwhile investigating (and possibly excluding) all the systematics which could produce a non-zero signal in the angular correlation function.
One of the possible causes for systematic deviations from isotropy is related to the non-uniform sky-exposure of BATSE (Fishman et al. 1994; Band 1996), which in principle could lead to measurements of anisotropy even for randomly distributed sources. We note however that such an effect would affect the observed distribution of both long and short GRBs in a similar fashion, since there is no a priori reason why different sub-classes should be observed by the instrument in different ways. The fact that long-lived sources exhibit a $`w(\theta )`$ which is highly consistent with the hypothesis of isotropy at all the angular scales probed by our analysis then tends to exclude responsibility from instrumental effects on the observed anisotropy of short-lived GRBs (see also Balazs et al. 1998; Meszaros et al. 2000).

Another possible cause for the observed anisotropic distribution can be attributed to intrinsic Poissonian fluctuations of short-lived GRBs due to their small number. Such an effect would however produce random fluctuations of the measured correlation function at all angular scales, while the data in our case show that an appreciable deviation from zero values is only obtained on small distances; measurements of $`w(\theta )`$ at all $`\theta >5^{}`$ are remarkably consistent with a null clustering signal. As a further test, we have re-binned sources belonging to the sub-class of short-lived GRBs and evaluated the two-point correlation function by means of equation (1) in the angular range $`1\theta 10`$ degrees.
Results from this analysis are presented in Figure 4: once again a marked deviation from isotropy is only detected on scales $`\theta 24`$ degrees, in full agreement with the findings of the previous Section.

Having excluded the effects of systematics on the detected anisotropy, we are then left with intrinsic non-randomness in the distribution of short-lived GRBs. This intrinsic anisotropic distribution would arise if e.g. some of the sources were repeating bursts, as indeed expected in several GRB models (Quashnock & Lamb 1993; Wang & Lingenfelter 1995; Meegan et al. 1995; Bennett & Rie 1996; Tegmark et al. 1996).

If we define $`f`$ as the fraction of observed bursts that can be labeled as repeaters and $`\nu `$ as the average number of observed events per repeating source ($`\nu 2`$) then, for positional errors drawn from a Gaussian distribution with standard deviation $`\mathrm{\Sigma }=\mathrm{\Delta }/(2\sqrt{2\mathrm{l}\mathrm{n}2})<<1`$ ($`\mathrm{\Delta }5^{}`$, see Figure 1), the correlation signal arising from the repeating-burst assumption can be written as (see Meegan et al. 1995)

$$w(\theta )=\frac{f(\nu 1)}{(N_{\mathrm{short}}1)}\left[\frac{2}{\mathrm{\Sigma }^2}\mathrm{exp}\left(\frac{\theta ^2}{2\mathrm{\Sigma }^2}\right)1\right].$$
(2)
Equation (2) was used to fit the data presented in Figure 4 for different values of the product $`f(\nu 1)`$. It turns out that the best-fit value is given by $`f(\nu 1)=0.07_{0.04}^{+0.06}`$, result which rejects the hypothesis of non-repetitivity at the $`2\sigma `$ confidence level.
By then fixing the average number of repeating events, we obtain the following constraints on the fraction of sources expected to repeat:

$`\nu =2f=0.07_{0.04}^{+0.06}`$
$`\nu =3f=0.035_{0.02}^{+0.035}`$
$`\nu =4f=0.023_{0.013}^{+0.02}`$
In other words, the data indicate that up to 13 per cent of short GRBs may have burst on average twice during the years corresponding to the lifetime of BATSE. This fraction obviously decreases for an increasing number of repetitions, going down to a more negligible $`2`$ per cent in the case of $`\nu =4`$.

The above explanation is however not the only one capable to account for the observed anisotropy on a few degrees scale. Another alternative which needs to be taken into account is the presence of large-scale structure, arising under the hypothesis that short GRBs trace the distribution of the galaxies in which they reside (note that here we are explicitly assuming GRBs to be extragalactic sources). In this framework, the angular correlation function can be modeled (to a first approximation) as the power-law $`w(\theta )=A\theta ^{(1\gamma )}`$ (see e.g. Peebles 1980).

Fits to the data presented in Figure 4 give values for the amplitude of the correlation function $`A0.3\pm 0.2`$, independent (within the errors) of the different values for the slope $`\gamma =1.6,2.0,1.8`$, chosen so to mimic the clustering properties respectively of star-forming and early-type galaxies, and of a mixture of both. The case for $`\gamma =1.8`$ is illustrated in Figure 4 by the dotted line. Despite the large errors associated to its determination, the allowed range of values for $`A`$ is in good agreement with results e.g. from the APM Survey (Loveday et al. 1995) based on samples of $`z<0.2`$ galaxies divided by morphological type. On the other hand, clustering analyses performed at higher redshifts find amplitudes which are much smaller than the values indicated by the GRB data ($`A<10^2`$ for $`z11.5`$ EROs, Roche et al. 2002; $`A10^3`$ for $`z3`$ Lyman Break galaxies, Porciani & Giavalisco 2002; $`A510^3`$ for $`z2.53.5`$ galaxies found in the Hubble Deep Field, Magliocchetti & Maddox 1999). If we then use large-scale structure to explain the clustering measurements presented in this work, we can conclude that – if confirmed – they imply short-lived GRBs to be hosted by relatively (possibly $`z<0.5`$) local galaxies, in agreement with the median redshift 0.5–0.8 predicted by neutron star/neutron star or black hole/neutron star merger models (Fryer et al. 1999). In fact (see e.g. MacFadyen & Woosley 1999), the alternative to produce short bursts in the context of the Collapsar model would meet some difficulties at $`z<0.8`$ since the intense mass loss from the progenitor - due to a higher metallicity with respect to that of stars at higher redshifts - could prevent the collapse of the helium core.

Obviously, one would like to discern between an anisotropic signal caused by burst repetition and one due to large-scale structure. Unfortunately, the still small number of detected sources does not allow to probe angular scales $`\theta <1`$ degrees where these two effects would greatly differentiate from one another. The solution to this problem is however within reach, since the goal of collecting hundreds of GRBs per year will be soon achieved by SWIFT (Burrows et al. 2002), which will also allow for precise position measurements in order to perform optical to X-ray spectrophotometry, with the aim to determine the origin of GRBs and their afterglows and to use bursts to investigate the early Universe. Moreover, the likely detection of several short bursts (which was severely limited in the BATSE experiment due to its lower sensitivity to extremely fast transients) by INTEGRAL and their location with sufficient accuracy through the IPN system will lead to a better understanding of the still mysterious population of short bursts.