text:
# The intrinsic colour dispersion in Type Ia supernovae

## 1 Introduction

In the last few years Type Ia supernovae (SNe Ia) have proved to be excellent distance estimators and have been successfully used to investigate the fate of the universe (Perlmutter et al., 1999; Riess et al., 1998). Despite the broad use of these objects by cosmologists, the current knowledge of the nature of SNe Ia is rather limited. Thus, there is a strong demand for further understanding to assess important issues about these explosions. For cosmological implications, the main concerns are related to the possible evolution of the SN properties with redshift. Moreover, critical tests for extinction along the line of sight based on supernova colours require good knowledge of the intrinsic properties of these objects.

In this paper, a statistical study on 48 well observed nearby SNe Ia is carried out. In particular, the intrinsic dispersion in SN colours is investigated using published $`BVRI`$ data. We also focus on the time correlation of intrinsic optical colours. This information is needed to address the possible host galaxy or intergalactic extinction by dust of supernovae used for cosmological tests (see e.g. Riess et al. (2000); Nobili et al. (2003)), and also to probe for other exotic sources of dimming at high-z with differential extinction (Mörtsell et al., 2002).

## 2 The data set

Published $`BVRI`$ lightcurves of well observed nearby SNe Ia were analysed. The considered sample consists of 48 SNe Ia from 2 different sets, the Calan/Tololo data published by Hamuy et al. (1996), and the set in Riess et al. (1999), usually referred to as the CfA data. The list of SNe is given in Table 1, along with the observed filter data available for each of them, their redshift and the $`B`$-band lightcurve “stretch”, s, as defined in Perlmutter et al. (1997) and Goldhaber et al. (2001).
The selected samples include a broad variety of SNe, well distributed in stretch factor, s. This parameter, related to $`\mathrm{\Delta }m_{15}`$ (Phillips, 1993), has been found to correlate with the supernova luminosity. Thus, a sample well distributed in stretch should imply a broad distribution in luminosity. Fig. 1 shows the distribution of the $`B`$-band light curve stretch factor and SN redshifts ($`0.003z0.12`$), for both samples. The timescale stretch parameter was determined from the lightcurve fits as in Goldhaber et al. (2001) using their $`B`$-band lightcurve template with a parabolic behavior for the earliest epoch after explosion.

$`K`$-corrections were applied to account for the small cosmological redshift as described in e.g. Goobar et al. (2002), using the ($`s=1`$) spectroscopic template of Nugent et al. (2002) as a starting point. The results of this analysis were used to improve the spectral template of SN Ia’s (Sect. 6) and we iterated the analysis once re-calculating the $`K`$-corrections with the improved template. Note that even though the $`K`$-corrections for the used data set are small, typically of the order of a few hundreds of a magnitude, for some of the more distant objects in the sample considered, they reach up to $`0.5`$ mag.

The SNe light curves were corrected for both Milky Way and host galaxy extinction as in Phillips et al. (1999) using the method first proposed by Lira (1995) for estimating host galaxy extinction using late epoch light curves. There is empirical evidence that the $`BV`$ colors of SNe Ia show extremely small scatter for the period between 30 and 90 days post $`B`$-band maximum, despite any difference in the light curve shapes at earlier epochs. As the set of SNe used in this article is a subset of the one analysed in Phillips et al. (1999), the host galaxy extinctions listed in Table 2 of their paper were used. These were derived combining the late time $`BV`$ colour with information on the $`BV`$ and $`VI`$ at maximum light.

The Cardelli et al. (1989) relation, modified by O’Donnell (1994), was used to compute the extinction in other colours given $`E(BV)`$. Spectral templates of Type Ia SNe were used to compute the evolution of the extinction with the supernova epoch.

The extinction corrected lightcurves were further screened to exclude the most peculiar SNe, as the main emphasis of this work is to establish the properties of “normal” supernovae. Fig. 2 shows the difference of the $`B`$-band lightcurve maximum, $`B_{\mathrm{max}}`$ and the $`V`$-band light curve maximum $`V_{\mathrm{max}}`$ plotted against the decline rate parameter, $`\mathrm{\Delta }m_{15}`$, as reported in Phillips et al. (1999). A $`3\sigma `$ clipping rejection criteria was applied, iterating until no data points were further rejected. At least 5 SNe in the 2 data sets deviate significantly from the expected linear relation derived in Phillips et al. (1999) from an independent set of 20 non-reddened SNe (solid line in Fig. 2):

$$B_{\mathrm{max}}V_{\mathrm{max}}=0.07(\pm 0.012)+0.114(\pm 0.037)(\mathrm{\Delta }m_{15}1.1)$$

The outliers are SN 1993ae, SN 1995bd, SN 1996ai, and SN 1996bk for the CfA data, and SN 1992K for the Calan/Tololo set. Moreover SN 1995E and SN 1995ac have been excluded since they are similar to the peculiar SN 1991T (Branch et al., 1993).111Including these 2 SNe in the analysis does not change the results significantely.

## 3 The average optical colours

Figs. 38 show the time evolution of the extinction and $`K`$-corrected $`BV`$, $`VR`$, $`RI`$, $`BI`$ and $`VI`$, colours. The plotted errors include the uncertainty on the host galaxy exinction correction. The time axis ($`t^{}`$) has been corrected by the SN redshift and the $`B`$-band stretch factor $`s`$, to account for the dependence of colours on the stretch. This rescaling of the time axis was found very effective in reducing the measured intrinsic dispersion (Sect. 4)

$$t^{}=\frac{tt_{B_{\mathrm{max}}}}{s(1+z)}$$
(1)
For each colour lightcurve we have spline-interpolated the weighted mean values computed in four days wide, non-overlapping bins. The aim of this procedure was to find a curve that describes the average time evolution of the colours. This parameterization will be referred as a “model” in the following discussion. Other methods were used to fit the data, as for example a least squares cubic spline fit. These result in curves that differ from the model typically by 0.01 magnitudes. The models show some systematic discrepancies between the two data sets, especially in $`RI`$ for which the C-T model is always redder than the CfA model. In Fig. 9 we investigate the differences of each of the models from the one built on both sets together. The differences are usually of the order of a few hundreds of a magnitude. The largest deviation was found for $`RI`$ in the Calan/Tololo set, resulting in a difference of about 0.2 mag at maximum with respect to the model built up using data from both sets. Note, however, that the statistics in the Calan/Tololo set for $`RI`$ colour, all along the evolution and in particular around the time of $`B_{\mathrm{max}}`$, is extremely poor, as shown in Fig. 5. Due to the smaller quoted observational error bars, the CfA set dominates the weighted average used to build the model out of both data sets.

A source of uncertainty in this analysis is the ability of the observers to convert the instrumental magnitudes from the used filter+CCD system transmission into the standard $`BVRI`$ system. The apparent systematic effects for $`RI`$ in Fig. 9 may be indicative of this.

In order to assess systematic effects in builing the colour models, we checked whether the host galaxy extinction was over-corrected. Thus, we compared the models with the data of those SNe that, according to Phillips et al. (1999), suffered no extinction from the host galaxy. The comparison, shown in Fig. 7, exhibits no obvious deviations and we may conclude that the extinction corrected colours are consistent with the uncorrected sub-sample.

A possible remaining dependence of colours on the stretch factor $`s`$ was investigated. For data points in a 5 days broad bin around time of B-maximum and around day 15, a linear weighted fit was done on the residuals from the models not showing any further significant dependence on $`\mathrm{\Delta }m_{15}`$ or equivalently on the stretch parameter $`s`$. Note that this is not in contradiction with Fig. 2, as we are considering $`(BV)_{t_{B_{\mathrm{max}}}}`$ instead of $`B_{\mathrm{max}}V_{\mathrm{max}}`$. Moreover, our analysis includes data points in a broad bin around day zero, and not a fitted lightcurve maximum as in Fig. 2.

## 4 Intrinsic colour dispersion

The observed colour dispersion around the models is overall inconsistent with a statistical random distribution of the data with the reported measurement errors. The $`\chi ^2`$ value per degree of freedom computed for each model is typically between 3 and 6, giving indication of a very poor fit of the data. Thus, we conclude that the analysed data supports the existence of an intrinsic colour dispersion of Type Ia supernovae 222This has been further investigated as described at the end of this section..

The intrinsic colour dispersion was computed on the residuals of each of the data sets from the corresponding model, built as explained above. Labeling $`XY_i`$ any of the measured colours, $`BV`$, $`VR`$, $`RI`$, $`BI`$ and $`VI`$, the residual with respect to the model expectation is referred as $`R_{XY_i}`$, i.e.,

$$R_{XY_i}=XY_iXY_i^{\mathrm{model}}$$
(2)
Next, the time axis was binned into 5 days wide bins to compensate for the modest statistics. The bins are centered at days $`d=5n`$, with $`n`$ ranging from 0 to 8. For each bin $`k`$, we compute the weighted average for the data $`XY`$, (see Appendix A)

$$XY^k=\frac{_{i=1}^{N_k}w_iXY_i}{_{i=1}^{N_k}w_i}$$
(3)
where $`N_k`$ is the number of points in bin $`k`$ and $`w_i`$ is the inverse of the uncertainty on the $`i`$-th measurement squared, $`1/\sigma _i^2`$, where $`\sigma _i`$ include both measurement errors and uncertainty on the host galaxy extinction corrections. The weighted sample standard deviation was computed as the square root of the weighted second moment:

$$s_{XY}^k=\sqrt{m_{w_2}}=\sqrt{\frac{_{i=1}^Nw_iR_{XY_i}^2}{_{i=1}^Nw_i}}$$
(4)
The uncertainty on the expression in Eq. (4) is given by the square root of the variance, $`V[s_{XY}^k]`$, computed as:

$$V[s_{XY}^k]=\frac{V[m_{w_2}]}{4m_{w_2}}$$
(5)
where $`V[m_{w_2}]`$ is the variance of the weighted second moment. We can consider the result of Eq. (4) as an estimate of the intrinsic dispersion in each bin, for each $`XY`$ colour: $`\mathrm{\Delta }_{XY}=s_{XY}^k`$. This has been applied at each one of the sets separately and at the whole set of data. The results of the method for the whole set of data are given in Table 2.

Fig. 10 shows a comparison of the intrinsic dispersion computed for each of the sets separately and their combination. In most of the cases the differences exceed the statistical uncertainties, seemingly arising from systematic differences between the two data sets. The computed intrinsic dispersion for the CfA sample was found to be smaller than for the CT data set and the whole (combined) sample in most cases, again pointing at systematic differences in the reported magnitudes.

Note that the adopted method leads to an overestimation of the intrinsic colour dispersion due to the contribution from the measurement errors. However, the weighting procedure ensures that the most accurate measurements dominate the result. To asses the impact of the measurement accuracy, we run a Monte-Carlo simulation to generate a synthetic colour data set with a dispersion given by the measurement uncertainties alone, i.e. no intrinsic dispersion. Three hundred data sets, with the same distributions in epochs and formal error bars as the CfA and Calan/Tololo were simulated and the weighted standard deviation (and its error) were computed, according to Eqs. (4) and (5). Note that the simulation, for simplicity, generates gaussian distributed and completely uncorrelated data. The averages, $`\delta `$, were used to disentangle the contribution of the intrinsic dispersion from the measurement errors. First, an hypothesis test was run to verify whether the simulated data and the measured data had the same dispersion; e.g. implying null intrinsic dispersion:

| $`H_0:\mathrm{\Delta }=\delta `$ |
| --- |
| $`H_1:\mathrm{\Delta }\delta `$ |

A level of significance $`\alpha =0.01`$ was set for rejecting the null hypothesis, (Cowan, 1998). Only 10 cases were not rejected, indicated by a (a) in Table 2. For all the other cases, for which the $`H_0`$ hypothesis was rejected, the intrinsic dispersion was computed as:

$$\mathrm{\Delta }^{\mathrm{corr}}=\sqrt{\mathrm{\Delta }^2\delta ^2}$$
(6)
and a lower limit on its value was set at a 99% confidence limit. The cases for which the null hypothesis was not rejected, were considered as compatible with no intrisic dispersion, and an upper limit on its value was set at a 99% confidence level. The corrected intrinsic dispersions are listed in the 5:th column of Table 2, together with upper and lower limits. We notice that the narrowest colour dispersion happens for $`VR`$, especially at late times. At 25 days after $`B_{\mathrm{max}}`$ and later, this colour is compatible with no intrinsic spread at all. Further, it should be noted that at day 35, all the colours but $`VI`$ are consistent with vanishing intrinsic dispersion.

## 5 Correlation

The correlation between optical colours at different epochs was also estimated. The property that was tested is whether a supernova that is blue at a certain epoch for example, say at maximum, stays blue at all epochs. In Riess et al. (2000), the authors argue that data measurements more than 3 days apart may be considered as uncorrelated estimators of colour. Our analysis does not support that assumption.333For high-z objects this is even more questionable when some of the main sources of uncertainty is the subtraction of a common image of the host galaxy and the $`K`$-corrections. We find significant correlations for data points up to a month apart, as shown below. The method followed is essentially the one used to compute the intrinsic dispersion. One can summarize the following steps:

* Bin the data in time
* Select only the SNe present in all the time bins
* For each bin compute the weighted average of the measurements belonging to the same SN.
* Compute the linear correlation coefficient between bins as in equation A-8.
* Test the correlation coefficient significance (Appendix A.2)

We refer to appendix A.2 for what follows. The correlation coefficient between different epochs h and k is:

$$r_{hk}=\frac{_{i=1}^n(R_{XY_i}^h\overline{R}_{XY}^h)(R_{XY_i}^k\overline{R}_{XY}^k)}{\sqrt{_{i=1}^n(R_{XY_i}^h\overline{R}_{XY}^h)^2_{j=1}^n(R_{XY_j}^k\overline{R}_{XY}^k)^2}}$$
(7)
where the summation is on the i-th SN, which by construction is present in all the bins. The uncertainty on $`r_{hk}`$ has been computed converting it into the normally distributed variable z, as described in (Appendix A.2). The results, given in Table 3, show that the correlation is important and non-zero all along the time evolution. Fig. 11 shows the weighted mean colours for the selected SNe in the bins centered at day 0, 15 and 30. Note that the supernovae selected are different in different colours, but, by construction, are the same in all the bins for each of the colour. It appears that supernovae that deviate from the average colour at a certain epoch are likely to keep their colour excess all along the 30 days evolution considered here.

## 6 Supernovae template

Nugent et al. (2002) analysed the relation between colours of SNe Ia and $`K`$-corrections. In particular they showed how $`K`$-corrections are mainly driven by the overall colour of the SN rather than by peculiarities of single features. Using spectra and colour light curves they give a recipe to build a SN Ia template, to be used for computing $`K`$-corrections of $`s1.0`$ SNe. As the sample used in our work is larger than the one considered in their paper, we have used our results to improve the spectral template. Note that this affects only the $`BVRI`$ magnitudes, which are the only bands treated here. Referring to Nugent et al. (2002), we proceed as follows:

* correct the $`BVRI`$ of their list of $`UBVRIJHK`$ magnitudes using our $`BV`$, $`VR`$ and $`RI`$ models
* correct the corresponding spectral templates with this colours444The same code used by Nugent et al. (2002) was used for this.
* iterate the whole analysis described in this paper, using the newly created templates to compute the $`K`$-corrections for both Calan/Tololo and CfA data set
* correct again the templates using the most recent version of the models

Even though the template itself is modified considerably from the original one, iterating once results in small corrections ($``$ 0.001-0.005 mag) in the intrinsic dispersion in all colours at any epoch, and a few percent in the colour models. As the correction obtained are small and well within the given uncertainties, we do not iterate further. The results of this analysis should not be considered definitive and may change as SNe Ia and their colour evolutions are studied in even more detail. However, it can be considered a good estimate of the average SN Ia template (in $`BVRI`$), as it is constructed from a quite broad sample of SNe. The final corrected template is available upon request.

## 7 Discussion

Throughout our analysis all data points were treated as independent measurements. This is particularly important for the cases with significant host galaxy light underneath the supernova where a single reference image was used, introducing a correlation between the data points not considered here.

A first estimate of the intrinsic dispersion was calculated for each data set separately and the two data sets together. A comparison of the results shows that the colours extracted from the CfA sample have smaller scatter indicating that the contribution from measurement errors is not negligible. Further, in the analysis the measurements are assumed to be Gaussian distributed and the weighted standard deviation has been taken as an estimate of the intrinsic dispersion. However, we noticed that 2 SNe, SN 1993H and SN 1992bc, seem to be rather deviant in $`BI`$ and $`VI`$ for the C-T set, as shown in Figs. 6 and 8. The effect of the “outliers” is particularly important aroud 20 days after B-maximum. For comparison, we recalculated the intrinsic spread in $`BI`$ and $`VI`$ excluding these 2 SNe. The results are shown on Table 4. The agreement between the intrinsic dispersion between the data sets improves when these 2 SNe are excluded. We emphasize that there are systematic differences between the Calan/Tololo and CfA data sets, and they might have been introduced while converting from the instrumental system used to the standard $`BVRI`$ system.

An attempt of disentangling the intrinsic dispersion from the contribution of the uncertainties has been done resulting in upper limit values for the intrinsic dispersion in some cases. The most relevant is the case of $`VR`$, that seems compatible with null intrinsic dispersion for most of the epochs at 99% C.L. Moreover this analysis brings out an important feature at day 35, when all the colours but $`VI`$ are consistent with zero intrinsic dispersion. This indicates that further studies will be needed to investigate this intriguing finding.

When computing the correlation coefficients we considered data up to 30 days after $`B`$-band maximum, even though later epoch data are available for several supernovae. However, the necessary condition of each SN being observed in all the time bins reduces the statistics if later epochs are introduced. This limitation is specific to the sample used and can be overcome with the use of a more extensively observed sample, such as what will be provide by SNfactory (Aldering et al., 2002).

The intrinsic dispersion sets constrains on the ability to determine the host galaxy extinction,$`A_V`$. This will depend on the colour, as the intrinsic dispersion is different for different colours. As an example we used the Cardelli et al. (1989) relation at the effective wavelength for each bandpass to compute the expected uncertainty in the extinction, neglecting any dependence on the supernova phase.

$$\begin{array}{ccc}\sigma _{A_V}^{\mathrm{B}\mathrm{V}}\hfill & =\hfill & 3.1\mathrm{\Delta }_{\mathrm{BV}}^{\mathrm{corr}}\hfill \\ \sigma _{A_V}^{\mathrm{V}\mathrm{R}}\hfill & =\hfill & 6.2\mathrm{\Delta }_{\mathrm{VR}}^{\mathrm{corr}}\hfill \\ \sigma _{A_V}^{\mathrm{R}\mathrm{I}}\hfill & =\hfill & 4.1\mathrm{\Delta }_{\mathrm{RI}}^{\mathrm{corr}}\hfill \\ \sigma _{A_V}^{\mathrm{B}\mathrm{I}}\hfill & =\hfill & 1.4\mathrm{\Delta }_{\mathrm{BI}}^{\mathrm{corr}}\hfill \\ \sigma _{A_V}^{\mathrm{V}\mathrm{I}}\hfill & =\hfill & 2.5\mathrm{\Delta }_{\mathrm{VI}}^{\mathrm{corr}}\hfill \end{array}$$
(8)
where $`R_V`$ was assumed equal to 3.1. Table 5 shows the results of the equation 8 for the epochs for which the intrinsic dispersion was calculated. The $`\mathrm{\Delta }_{XY}^{\mathrm{corr}}`$ was used for this. The results indicate that, with the present knowledge, extinction by dust with $`R_V=3.1`$ may only be determined to $`\sigma _{A_V}>0.10`$ with Type Ia restframe optical data within the first 40 days after $`B`$-band lightcurve maximum, for the colours and epochs with non-zero intrinsic dispersion. However observations in $`VR`$ are preferable to other colours to set limit on the extinction of an observed supernova. To account for possible different extinction parameters $`R_V`$ in the different supernova host galaxies (as noticed in e.g. (Riess et al., 1996) and (Krisciunas et al., 2000)), we repeated the analysis considering a Gaussian uncertainty $`\sigma _{R_V}=1`$ on $`R_V`$, propagating this scatter on the host galaxy extinction corrections. This changes the values of the intrinsic dispersion, but typically within the quoted errors in Table 2.

## 8 Conclusion

A statistical analysis of colours of SNe Ia using a sample of 48 nearby SNe was performed. Of this sample 7 SNe have been excluded based on their large extinction, or peculiar behavior of their colours at maximum. With the present knowledge and data quality we computed the average colour evolution for $`BV`$, $`VR`$, $`RI`$, $`BI`$ and $`VI`$ with time and the derived intrinsic scatter. We find that the correlation of colours during the first 30 days after restframe $`B`$-band maximum is not negligible, i.e. arguing against the assumptions made in the analysis of SN 1999Q (z=0.46; Riess et al., 2000) where five measurements of restframe $`BI`$ along its lightcurve were treated as independent estimates of extinction. According to our findings, their limit on the presence of intergalactic grey dust must be revised. A reanalysis of the data from SN 1999Q is in preparation (Nobili et al., 2003). With the data at hand, host galaxy extinction corrections from restframe optical colours within the first 30 days after maximum light are generaly limited to $`\sigma _{A_V}>0.1`$ due to the intrinsic variation of Type Ia colours at those epochs, with the possible exception of extinction corrections derived from the rest-frame $`VR`$ colour. The results of this analysis have been used to correct spectroscopic templates which are available upon request.

###### Acknowledgements.

S. Nobili wishes to thank Christian Walck for his valuable help and constructive discussions on the statistical analysis presented in this paper and Don Groom, Chris Lidman and Vallery Stanishev for their careful reading and useful suggestions. SN is supported by a graduate student grant from the Swedish Research Council. AG is a Royal Swedish Academy Research Fellow supported by a grant from the Knut and Alice Wallenberg Foundation. PEN acknowledges support from a NASA LTSA grant and by the Director, Office of Science under U.S. Department of Energy Contract No. DE-AC03-76SF00098.

## Appendix A Appendix

### A.1 Algebraic moments

In what follows we refer to Cowan (1998) and Kendall & Stuart (1958) (see Volume 1, Chapter 10), our aim here is to give uniformity to the notations used. Given n independent observations of a variable, $`x_1,x_2,\mathrm{},x_n`$, the r-th moment or algebraic moment is given by:

$$m_r=\frac{1}{n}\underset{i=1}{\overset{n}{}}x_i^r$$
(A-1)
The expectation value, or mean value, and the variance of $`m_r`$ are:

| $`E[m_r]=\mu _r`$ |
| --- |
| $`V[m_r]=\frac{1}{n}(\mu _{2r}\mu _r^2)`$ |

(A-2)
where $`\mu _{2r}`$ is the $`2r`$-th moment. The A-2 are exact formulae as long as one knows $`\mu _{2r}`$ and $`\mu _r`$. However this is not always the case, and one has to use their estimators $`m_{2r}`$ and $`m_r`$ from the sample itself. The variance of the standard deviation can be computed using standard error propagation:

$$V[\sigma ]=V[\sqrt{m_2}]=\frac{V[m_2]}{4m_2}=\frac{m_4m_2^2}{4nm_2}$$
(A-3)
In experimental situation each observation $`x_i`$ is often attached to a certain weight $`w_i`$. Supposing that the weight themselves are known without errors, one can define the following formula for the r-th weighted moments:

$$m_{w_r}=\frac{_{i=1}^nw_ix_i^r}{_{i=1}^nw_i}$$
(A-4)
In order to evaluate the variance on the weighted standard deviation we simply extended the A-3 for the case of weighted moments, obtaining:

$$V[\sigma _w]=V[\sqrt{m_{w_2}}]=\frac{V[m_{w_2}]}{4m_{w_2}}=\frac{m_{w_4}m_{w_2}^2}{4n_{eff}m_{w_2}}$$
(A-5)
where the $`m_{w_4}`$ and $`m_{w_2}`$ are the 4th and 2nd weighted moments respectively and $`n_{eff}=(_{i=1}^nw_i)^2/_{i=1}^nw_i^2`$.

A simple Monte Carlo simulation has been run to verify the accuracy of the approximated formula A-5.

### A.2 Correlation coefficient

Given 2 random variables x and y the correlation coefficient is defined as

$$\rho _{xy}=\frac{V_{xy}}{\sigma _x\sigma _y}$$
(A-6)
The unbiased estimator of the covariance $`V_{xy}`$ is:

$$\widehat{V}_{xy}=\frac{1}{n1}\underset{i=1}{\overset{n}{}}(x_i\overline{x})(y_i\overline{y})$$
(A-7)
so that the estimator of the correlation coefficient will be:

$$r_{xy}=\frac{\widehat{V}_{xy}}{s_xs_y}=\frac{_{i=1}^n(x_i\overline{x})(y_i\overline{y})}{\sqrt{_{i=1}^n(x_i\overline{x})^2_{j=1}^n(y_j\overline{y})^2}}$$
(A-8)
For the A-8 to be a good estimator of the correlation coefficient there are a few caveats to check. First condition is that the samples are randomly defined from the population, that is to make sure that the samples x and y are not selected in some ways that would operate as to increase or decrease the value of r. Even though one has random samples it is possible to compute the errors due to sampling. Commonly this is computed as $`(1r^2)/\sqrt{n}`$. Unfortunately this is just an approximation. Moreover r’s for successive samples are not distributed normally unless n is large and the true value $`\rho `$ is near zero. This yields to a distinction:

* if $`n>30`$
in order to know whether the value calculated for r is significantly different from zero, one can compute its standard error as:

$$\sigma _r=1/\sqrt{n1}$$
(A-9)
If $`r/\sigma `$ is greater than 2.58, one can conclude that the universe value of r is likely to be greater than zero.
* if $`n<30`$
the variable:

$$t=r\frac{\sqrt{n2}}{\sqrt{1r^2}}$$
(A-10)
follows the t-distribution with dof=n-2. This can be used only for testing the hypothesis of zero correlation.

R. A. Fisher developed a technique to overcome these difficulties. The variable r is transformed into another variable that is normally distributed. This is especially useful for high value of r, when none of the above test can be safely applied. The transformation to the variable z:

$$z=\frac{1}{2}\mathrm{ln}(1+r)\frac{1}{2}\mathrm{ln}(1r)$$
(A-11)
allows some important semplifications. The distribution of z’s for successive samples does not depend on the universe value $`\rho `$ and the distribution of z for successive samples is so near to normal that it can be treated as such without any loss of accuracy, see for example (Volume 1, Chapter 16 and Volume 2,Chapter 26) of Kendall & Stuart (1958). Moreover the standard error for z is independent on its $`\sigma `$:

$$\sigma _z=1/\sqrt{n3}$$
(A-12)
The way to proceed is then very simple:

* Compute r according to ( A-8)
* Transform r into z according to ( A-11)
* Compute $`\sigma _z`$ according to ( A-12)
* In order to state a 1 sigma confidence limit for r, transform the 2 values $`z\pm \sigma _z`$ back to r, using the inverse of ( A-11)


format: text

text:
# On signals of new physics in global event properties in pp collisions in the TeV energy domain

## 1 Introduction

Possible scenarios for pp collisions in the TeV energy domain, based on extrapolations from data and fits in the GeV region to charged particle multiplicity distributions, $`P_n`$, and on the weighted superposition mechanism of two classes of events, i.e., soft (without minijets) and semihard (with minijets), each described by a negative binomial multiplicity distribution (NBMD) with characteristic parameters ($`\overline{n}_{\text{soft}},k_{\text{soft}}`$ and $`\overline{n}_{\text{sh}},k_{\text{sh}}`$, respectively) have been discussed both in full phase space  and rapidity intervals . (For an experimental investigation of the substructures see Ref. .)

It has been pointed out that in the semihard component either with a strongly KNO scaling violating mechanism (parameter $`k_{\text{sh}}`$ decreases with c.m. energy $`\sqrt{s}`$ as $`1/\mathrm{log}\sqrt{s}`$) or with a QCD inspired KNO scaling violating behaviour, (which is milder than the previous one as parameter $`k_{\text{sh}}`$ decreases as $`1/\sqrt{\mathrm{log}s}`$), the average number of clans, $`\overline{N}_{\text{sh}}`$, becomes smaller as the c.m. energy increases (from $`23`$ and 22, resp., at 900 GeV, to $`11`$ and 18, resp., at 14 TeV) and the corresponding average number of particles per clan, $`\overline{n}_{c,\text{sh}}`$, much larger (from $`2.5`$ and 2.6, resp., at 900 GeV, to $`7`$ and 5, resp., at 14 TeV), favouring clan aggregation and higher average particle population within clans. This behaviour was considered suggestive and indicative of a class of events with high particle density.

Motivated by this remark, we decided to explore in the same framework the consequences of an eventual more dramatic decrease of the average number of clans towards 1, i.e., towards a class of events with the maximum clan and particle aggregation, a situation which implies a $`k_{\text{sh}}`$ value less than one but which is unfortunately met in the semihard component only asymptotically at an extraordinary high c.m. energy in the KNO scaling violating extrapolations of $`k_{\text{sh}}`$. This result contrasts sharply with the main properties of the class of events belonging to the semihard component at non-asymptotic energies, where $`k_{\text{sh}}`$ is always larger than or equal to one. It raises the intriguing question whether a very high clan and particle density is only an asymptotic behaviour of the semihard component, or the property of an effective new class of events, different from (maybe harder than) the semihard one, and whose onset might happen already at 14 TeV c.m. energy: if this is the case, its contribution to the total charged particle multiplicity distribution, $`P_n^{(\text{total})}`$, in terms of a new weighted NB(Pascal)MD, should be added to those of the soft and semihard ones of Ref.  and . The claim is that the benchmark of this third new class of events (called ‘th’ from now on) is $`k_{\text{th}}<1`$.

It should be recalled that the average number of clans is in general a non trivial function of the average number of particles of each component, $`\overline{n}_i`$, and of the parameter $`k_i`$,

$$\overline{N}_i=k_i\mathrm{ln}(1+\overline{n}_i/k_i)$$
(1)
with $`i`$=soft, sh and th, enumerating the three classes. Whereas $`k_{\text{soft}}`$ is taken constant in our scenarios in the new energy domain and therefore $`\overline{N}_{\text{soft}}`$ depends only on the average multiplicity $`\overline{n}_{\text{soft}}`$ as the c.m. energy increases, the decrease of $`k_{\text{sh}}`$ with c.m. energy suggested by our assumptions on KNO scaling violations remains at non-asymptotic energies always larger than or equal to one in the semihard component and is only in part contrasted by the increase of $`\overline{n}_{\text{sh}}`$, a fact which outlines $`k_{\text{sh}}`$ dominance in the general $`\overline{N}_{\text{sh}}`$ behaviour. As already mentioned the request of an average number of clans, $`\overline{N}_{\text{th}}`$, of few units in the third component leads to values of $`k_{\text{th}}`$ less than one whereas that of $`\overline{N}_{\text{th}}1`$ implies $`k_{\text{th}}0`$ for large $`\overline{n}_{\text{th}}`$. Notice that in this last extreme situation the average number of particle per clan, $`\overline{n}_{c,\text{th}}`$ coincides of course with the average number of particles of the new component $`\overline{n}_{\text{th}}`$. In addition when $`\overline{n}_{\text{sh}}k_{\text{sh}}`$ —a quite normal situation at large c.m. energies— the NB(Pascal)MD describing the final charged particle MD of the component is well approximated by a log-concave gamma MD for $`k_{\text{sh}}`$ larger than one which becomes an exponential when $`k_{\text{sh}}`$ is equal to one: the maximum of the gamma distribution occurs at $`n/\overline{n}_{\text{th}}=1k_{\text{th}}^1`$. (See Fig. 1). This behaviour should be compared with what happens in the third component under the condition $`k_{\text{th}}\overline{n}_{\text{th}}`$: in this case the NBMD is in fact well approximated by a log-convex gamma distribution which for a certain value of the above parameter close to zero leads to the average number of clans $`1`$ and to the total MD well described by a logarithmic one. This result is consistent with the standard interpretation of the occurrence of the negative binomial MD for final charged multiplicity at hadron level as a two step process in which the independently emitted clans contain logarithmically distributed charged particles .

Finally in consequence of the assumed quite extreme clan and particle aggregation into few clans, events with quite large forward-backward multiplicity correlations close to the maximum allowed leakage from one hemisphere to the opposite one  are also expected. At partonic level the mentioned remarks would suggest for the third component high parton density clan production with huge colour exchange processes originated from a relatively small number of high virtuality ancestors, which would indicate probably an emission mechanism harder than that seen in the semihard component.

In conclusion it might well be that one could observe already at 14 TeV c.m. energy in pp collisions three classes of events or components instead of two, the first and the second class being those examined in Refs  and  and the third one fully characterised by the reduction of the the average number of clans to the minimum allowed by the condition $`k_{\text{th}}<1`$. The classes of events or components contributing to the total $`n`$ charged particle multiplicity distributions would therefore be the following:

I) the class of soft events (events with no minijets), which in clan structure analysis would be characterised by quite large values of $`\overline{N}_{\text{soft}}`$ and quite small $`\overline{n}_{c,\text{soft}}`$; $`P_n^{(\text{soft})}`$ obeys KNO scaling and $`k_{\text{soft}}`$ is assumed to be constant from the GeV region to the new energy domain;

II) the class of semihard events (events with minijets), for which $`\overline{N}_{\text{sh}}`$ decreases quickly with c.m. energy and $`\overline{n}_{c,\text{sh}}>\overline{n}_{c,\text{soft}}`$; $`k_{\text{sh}}`$ decreases also with c.m. energy but its value is larger than or equal to one and KNO scaling is violated;

III) the new class of events (events generated probably by quite hard partons) with $`k_{\text{th}}<1`$, $`\overline{n}_{\text{th}}k_{\text{th}}`$ (and quite small $`\overline{N}_{\text{th}}`$), with large forward-backward multiplicity correlations.

Accordingly, the total charged particle MD will have the following expression in terms of the weighted composition of the three NBMD’s, one for each component:

$$P_n^{(\text{total})}=\alpha _{\text{soft}}P_n^{(\text{soft})}+\alpha _{\text{sh}}P_n^{(\text{sh})}+(1\alpha _{\text{soft}}\alpha _{\text{sh}})P_n^{(\text{th})}$$
(2)
where $`\alpha _{\text{soft}}`$ and $`\alpha _{\text{sh}}`$ are the fraction of soft and semihard events, respectively.

At 14 TeV c.m energy the weight of the third component is expected to be of course quite small, in our calculations we assumed it to be between 1 and 3 percent of the total number of events.

Since the quantitative properties of the classes of soft and semihard events have been extensively discussed in Refs.  and  we focus our attention in this paper on the global properties of the new class of events, which we would be tempted to call hard events.

## 2 A new class of events, a third component in pp collisions?

Assuming for simplicity that the new class of events is described —as are the soft and semihard ones— by a NBMD with parameters $`\overline{n}_{\text{th}}`$ and $`k_{\text{th}}`$, we decided to examine under such conditions possible signals of new physics which should be easily detectable already at 14 TeV c.m. energy in pp collisions by studying the general behaviour of the various sets of parameters characterising the final charged particles MD of the component $`P_n(\overline{n}_{\text{th}},k_{\text{th}})`$. In order to stress the properties of the new class of events, the extreme case $`\overline{N}_{\text{th}}=1`$ will be discussed in the following.

### 2.1 Consequences in terms of the $`\overline{n}_{\text{th}}`$, $`k_{\text{th}}`$ parametrisation of the NBMD.

The request that $`\overline{N}_{\text{th}}=1`$, through its definition Eq. (1), leads to the straightforward relation between $`\overline{n}_{\text{th}}`$ and $`k_{\text{th}}`$:

$$\overline{n}_{\text{th}}=k_{\text{th}}(e^{1/k_{\text{th}}}1)$$
(3)
with (being $`\overline{n}_{\text{th}}k_{\text{th}}`$)

$$k_{\text{th}}<1.$$
(4)
In Fig. 2 are plotted $`\overline{n}_{\text{th}}P_n`$ vs $`n/\overline{n}_{\text{th}}`$ in KNO form for $`\overline{N}_{\text{th}}=1`$ and different average multiplicities, with the respective values of $`k_{\text{th}}`$ obtained via Eq. (3), which highlights the general trend of the exponential behaviour of the MD and $`k_{\text{th}}<1`$ dominance. In addition we notice that for $`k_{\text{th}}<1`$, as for the gamma distribution corresponding to the limit $`\overline{n}k`$ of the NBMD, we can distinguish two regions: $`n>\overline{n}_{\text{th}}`$ and $`n<\overline{n}_{\text{th}}`$. In the former a smooth behaviour of $`\overline{n}_{\text{th}}P_n`$ is seen as $`n`$ increases (events with large multiplicities) whereas in the latter a dramatic decrease of $`\overline{n}_{\text{th}}P_n`$ is visible (events with low multiplicities). The fact that the distribution decreases so slowly for $`n>\overline{n}_{\text{th}}`$ and so fast for $`n<\overline{n}_{\text{th}}`$ produces more easily events with very large or very small charged multiplicity, which reminds us of what occurs in cosmic rays for centauro and anti-centauro events. A provocative result to be tested in experiments.

### 2.2 Consequences in terms of clan structure analysis

In Fig. 3 are plotted $`\overline{N}`$ and $`\overline{n}_c`$ for the three components as a function of c.m. energy in full phase space: the first two components are those discussed in , notice that $`k_{\text{sh}}`$ behaviour corresponds to that suggested by a strong KNO scaling violation. The lack of experimental data on the third component at lower c.m. energies and the consequent impossibility to guess its behaviour in terms of extrapolations as it was done for the soft and semihard components in Ref. , led us to show a band of possible values. A conservative guess that the total multiplicity variation due to the third component over the extrapolation of Ref.  is limited to 10% leads to $`\overline{n}_{\text{th}}`$ values from 3 to 10 times larger than the total multiplicity, i.e., to a third component weight ranging from 1 to 3%. Notice that to a small variation of $`k_{\text{th}}`$ according to Eq. (3) corresponds a quite large variation of $`\overline{n}_{\text{th}}`$. This is clearly shown in the figure, where $`\overline{N}_{\text{th}}=1`$ whereas $`\overline{n}_{c,\text{th}}`$ is 3 times larger in one case with respect to the other. Last but not least, $`\overline{n}_{c,\text{th}}`$ is, as assumed, much larger than in the two other components.

### 2.3 Consequences in terms of parameters $`a_{\text{th}},b_{\text{th}}`$ of the NBMD.

Since

$$\frac{(n+1)P_{n+1}^{(\text{th})}}{P_n^{(\text{th})}}=a_{\text{th}}+b_{\text{th}}n,$$
(5)
with

$$a_{\text{th}}=\frac{k_{\text{th}}\overline{n}_{\text{th}}}{k_{\text{th}}+\overline{n}_{\text{th}}},$$
(6)
and

$$b_{\text{th}}=\frac{\overline{n}_{\text{th}}}{k_{\text{th}}+\overline{n}_{\text{th}}},$$
(7)
it follows that

$$a_{\text{th}}=\overline{N}_{\text{th}}P_{\text{log}}(1)=\overline{N}_{\text{th}}\frac{b_{\text{th}}}{\mathrm{ln}(1b_{\text{th}})},$$
(8)
(where $`P_{\text{log}}`$ is the logarithmic distribution) which for $`\overline{N}_{\text{th}}=1`$ leads to

$$a_{\text{th}}=b_{\text{th}}/\mathrm{ln}(1b_{\text{th}})=P_{\text{log}}(1).$$
(9)
It should be noticed that in the limit $`k_{\text{th}}0`$ (it corresponds to $`a_{\text{th}}0`$ and $`b_{\text{th}}1`$) a truncated NBMD (i.e., with zero multiplicity missing) and with constant $`\overline{n}_{\text{th}}/k_{\text{th}}`$ (i.e., constant $`b_{\text{th}}`$) leads to a logarithmic MD, i.e., to the MD of a single average clan:

$$P_n^{(\text{th})}\frac{(b_{\text{th}})^{n1}}{n}P_{\text{log}}(1)\stackrel{a_{\text{th}}0}{=}\frac{(b_{\text{th}})^n}{n\mathrm{ln}(1b_{\text{th}})}=P_{\text{log}}(n).$$
(10)
This result is of course consistent with the standard interpretation of the occurrence of the NBMD as a two step process (independently produced clans with at least one particle ancestor decay according to a logarithmic MD) when the average number of clans is reduced to one.

In Fig. 4 are plotted the c.m. energy dependence of $`a`$ and $`b`$ for the three components which highlight the general trend of $`a_{\text{th}}0`$ and $`b_{\text{th}}1`$ for $`\overline{N}_{\text{th}}1`$. Notice that in general, for $`\overline{n}_ik_i`$, one has $`a_ik_i`$.

### 2.4 Clan aggregation and correlations.

Being the probability of two particles of the hard component to join the same clan $`1/k_{\text{th}}`$ times larger than the probability to join two different clans  it is clear that $`k_{\text{th}}<1`$ is enhancing particle aggregation properties within clans.

Being $`\overline{n}_{\text{th}}^2/k_{\text{th}}`$ much larger than in the semihard component, from the relation

$$\overline{n}_{\text{th}}^2/k_{\text{th}}=C_2(\eta _1^{},\eta _2^{\prime \prime })𝑑\eta _1^{}𝑑\eta _2^{\prime \prime }$$
(11)
two particle correlations are also expected to be larger. Since cumulants are dependent on $`1/k_{\text{th}}`$, which is again much larger than $`1/k_{\text{sh}}`$, higher order cumulants are also expected to be quite large in the third component.

Concerning forward-backward multiplicity correlations and leakage parameter behaviour, it should be pointed out, as discussed in Ref. , that the forward-backward multiplicity correlation strength is given for each component by

$$b_{\text{FB},i}=\frac{2b_ip_i(1p_i)}{12b_ip_i(1p_i)},$$
(12)
where $`p_i`$ is the average fraction of particles, within each clan, which remain in the same hemisphere where the clan was emitted and do not ‘leak’ to the opposite hemisphere. In the third component, in presence of only one clan, the leakage is expected to be maximum ($`p_{\text{th}}1/2`$); being at the same time, as we have seen, $`b_{\text{th}}1`$, one obtains $`b_{\text{FB},\text{th}}1`$.

Accordingly forward-backward multiplicity correlations in the third component are expected to be quite stronger than in the semihard component.

### 2.5 What Monte Carlo event generators say on the third component at 14 TeV c.m. energy in pp collisions

It is interesting to point out that the total charged particle MD, $`P_n`$, of the events generated with the last version of the Pythia Monte Carlo model  both in full phase space and in a restricted rapidity interval ($`|\eta |<0.9`$) cannot be fitted in terms of the weighted superposition of two NBMD’s and that the plot of $`P_n`$ vs $`n`$ shows two shoulders (see Fig. 5). These two shoulders appearing in the total multiplicity distribution come, in our framework, from the weighted superpositions of the first with the second and of the second with the third component. This is true also in the pseudo-rapidity interval. It would be tempting to identify the third component with the hard one discussed above, but striking differences between the third class of events generated with Pythia Monte Carlo and that one discussed above prevent us from doing this identification.

An “elbow structure” is expected in $`P_n^{(\text{total})}`$ as the result of the weighted superposition of the NB MD of the semihard component with the gamma log-convex MD of the hard one, as shown in Fig. 6.

In the figure are shown in addition to the weighted superposition of the soft with the semihard components already seen in Ref  the results of the weighted superposition of the semihard and hard components at 14 TeV c.m. energy.

### 2.6 The case of few clans

It should be pointed out that similar results to those discussed so far for the extreme case $`\overline{N}_{\text{th}}=1`$ can be obtained with a small average number of clans, provided that $`k_{\text{th}}<1`$, as illustrated in Table 1. The closer $`k_{\text{th}}`$ is to 1, the less evident becomes the elbow structure discussed above.

## 3 Conclusions

Signals of new physics could be visible in the tail of the $`n`$ charged particle multiplicity distribution $`P_n`$ in complete disagreement with Pythia Monte Carlo predictions under the assumption that the aggregation of the average number of clans in the semihard component reaches its minimum not asymptotically but at lower c.m. energy. The occurrence of this minimum would suggest the onset of a third (maybe hard) component to be added to the soft and semihard ones discussed in previous work.

In the paper we examined the main consequences of the occurrence of the mentioned situation already at 14 TeV c.m energy in pp collision, assuming that the new class of events is between 1 and 3 percent of the total.

## Acknowledgements

We would like to thank M. Monteno for providing the Pythia Monte Carlo simulation results.

## References


text:
# Particle Transport in Tangled Magnetic Fields and Fermi Acceleration at Relativistic Shocks

## 1 Introduction

Fermi acceleration at relativistic shocks is an important topic for understanding the formation of spectra of ultrarelativistic particles and radiation in relativistic flows such as those observed in active nuclei, microquasars, $`\gamma `$ray bursts and pulsar wind nebulae (see Kirk & Duffy 2001 and references therein). Of particular interest is the acceleration timescale that can be as short as a Larmor time for relativistic Fermi acceleration; the smaller return probability to the shock for downstream particles, as compared to the non-relativistic regime, is compensated by the much larger energy gain at each cycle. However the study of Fermi acceleration in the relativistic limit is more involved than in the non-relativistic regime due to the increased importance of the anisotropy of the distribution function (Gallant & Achterberg 1999).

Various methods have been used to study the relativistic regime of Fermi acceleration (see Kirk & Duffy 2001 and references therein), starting with analytical estimates by Peacock (1981), followed by semi-analytical methods (Kirk & Schneider 1987; Gallant & Achterberg 1999; Kirk et al. 2000; Achterberg et al. 2001) and numerical Monte-Carlo techniques (Ballard & Heavens 1992; Ostrowski 1993; Bednarz & Ostrowski 1998), which in spite of their differences converge to an asymptotic spectral index $`s2.22.3`$ in the ultra-relativistic limit.

In the present Letter, we propose a new numerical Monte-Carlo method of simulation of the acceleration of test particles at relativistic shocks. The trajectories of particles in the upstream and downstream inhomogeneous magnetic fields are integrated out exactly from the entry of each particle through the shock until its first return to the shock. The law of probability of return to the shock as a function of ingress and egress pitch angles is then constructed by Monte-Carlo iteration. Finally we combine these probability laws, one defined for upstream and the other for downstream, with the energy gain formula at shock crossing to simulate the acceleration process. This latter use of the angular probability laws is similar to the method of Gallant et al. (2000) with some differences to be discussed below.

The present method appears efficient and potentially more powerful when compared to direct Monte-Carlo simulations which follow each particle through its repeated shock crossings (e.g. Ballard & Heavens 1992; Ostrowski 1993). It allows one to simulate relativistic Fermi acceleration in any magnetic configuration, albeit for test particles only. We describe the method and numerical simulations in Section 2 and then present the results for a planar shock with fully turbulent magnetic field both upstream and downstream in Section 3.

## 2 Numerical simulations

The hydrodynamic jump conditions at an adiabatic strong shock, neglecting magnetic fields, are given in Blandford & Mc Kee (1977), Kirk & Duffy (2001) and Gallant (2002). The shock Lorentz factor is $`\mathrm{\Gamma }_\mathrm{s}`$ (upstream or lab frame111Unless otherwise noted, all quantities are calculated in the upstream (lab) frame; wherever needed, quantities relative to a given frame but calculated in an other will be marked with the subscript $`|`$, e.g., $`\beta _{\mathrm{s}|\mathrm{d}}`$ refers to the shock velocity measured in the downstream rest frame.), and the relative Lorentz factor $`\mathrm{\Gamma }_\mathrm{r}`$ between upstream and downstream $`\mathrm{\Gamma }_\mathrm{r}\mathrm{\Gamma }_\mathrm{s}\mathrm{\Gamma }_{\mathrm{s}|\mathrm{d}}(1\beta _\mathrm{s}\beta _{\mathrm{s}|\mathrm{d}})`$. The downstream Lorentz factor $`\mathrm{\Gamma }_{\mathrm{s}|\mathrm{d}}`$ (and $`\mathrm{\Gamma }_\mathrm{r}`$) can be obtained as a function of $`\mathrm{\Gamma }_\mathrm{s}`$ from the relations derived from the shock jump conditions given in Gallant (2002). In particular, in the ultra-relativistic limit $`\mathrm{\Gamma }_\mathrm{s}+\mathrm{}`$, one finds the well-known results $`\beta _{\mathrm{s}|\mathrm{d}}1/3`$ and $`\mathrm{\Gamma }_\mathrm{r}\mathrm{\Gamma }_\mathrm{s}/\sqrt{2}`$.

We conduct our simulations in two steps. We first perform Monte-Carlo simulations of particle propagation in a given magnetic field structure following Casse, Lemoine & Pelletier (2001) (wherein one may find more details on the numerical procedure). These simulations are carried out separately in the downstream and in the upstream rest frames. It is possible to set up any magnetic field structure including regular and tangled components, but in the following, for the sake of simplicity, we restrict ourselves to the case of pure Kolmogorov turbulence in both downstream and upstream rest frames. The equations of motion of each particle are integrated out exactly, the magnetic field being calculated at each point of the trajectory as a sum of plane waves, using 200 modes spaced logarithmically on three decades of wavelength and whose wavevector directions are drawn at random222Using a higher number of modes does not modify the results shown here as we have checked.. Ostrowski (1993) used a similar technique to construct the magnetic field albeit with 3 modes only, while Ballard & Heavens (1992) used three-dimensional FFT methods (see Casse, Lemoine & Pelletier 2001 for a comparison of these methods).

The trajectories are integrated over 100 scattering times upstream and 1000 scattering times downstream in order not to miss possible late returns to the shocks. The laws of return probability as a function of ingress and egress pitch angles are then constructed in the following way. We draw at random a point along a simulated trajectory which defines the point of entry through the shock. The ingress pitch angle cosine to shock normal $`\mu ^\mathrm{i}`$ is recorded, the trajectory scanned to find the point of exit through the shock, and the egress pitch angle cosine $`\mu ^\mathrm{e}`$ is then recorded. Iteration of the above then yields the desired law of conditional return probability $`𝒫(\mu ^\mathrm{i};\mu ^\mathrm{e})`$ which gives the probability for a particle entering with a pitch angle cosine $`\mu ^\mathrm{i}`$ to return to the shock with a pitch angle cosine $`\mu ^\mathrm{e}`$. The simulations also give a direct measurement of the return time to the shock as a function of pitch angles.

Once the upstream and downstream laws of return probability, respectively $`𝒫_\mathrm{u}(\mu _\mathrm{u}^\mathrm{i};\mu _\mathrm{u}^\mathrm{e})`$ and $`𝒫_\mathrm{d}(\mu _\mathrm{d}^\mathrm{i};\mu _\mathrm{d}^\mathrm{e})`$ are known, the simulation of the acceleration process itself can be performed as follows. We denote by $`f_\mathrm{d}^{2n+1}(\mu _\mathrm{d},ϵ_\mathrm{d})`$ the distribution function of particles that enter the shock to downstream and that have experienced $`2n+1`$ shock crossings, and $`f_\mathrm{u}^{2n}(\mu _\mathrm{u},ϵ_\mathrm{u})`$ similarly upstream: particles are upstream for an even number of shock crossings and $`f_\mathrm{u}^0`$ represents the injection population. These distribution functions are normalized to the total number of particles $`N`$ injected such that, in the absence of escape from the acceleration site, after $`2n`$ shock crossings $`N=d\mu _\mathrm{u}dϵ_\mathrm{u}f_\mathrm{u}^{2n}(\mu _\mathrm{u},ϵ_\mathrm{u})`$, and after $`2n+1`$ shock crossings $`N=d\mu _\mathrm{d}dϵ_\mathrm{d}f_\mathrm{d}^{2n+1}(\mu _\mathrm{d},ϵ_\mathrm{d})`$. The conservation of particle number at shock crossing $`ud`$ and Lorentz transforms of pitch angles and energies imply:

$$f_\mathrm{d}^{2n+1}(ϵ_\mathrm{d},\mu _\mathrm{d}^\mathrm{i})\mathrm{d}\mu _\mathrm{d}^\mathrm{i}\mathrm{d}ϵ_\mathrm{d}=\left[_{\beta _\mathrm{s}}^1d\mu _\mathrm{u}^\mathrm{i}𝒫_\mathrm{u}(\mu _\mathrm{u}^\mathrm{i};\mu _\mathrm{u}^\mathrm{e})f_\mathrm{u}^{2n}(ϵ_\mathrm{u};\mu _\mathrm{u}^\mathrm{i})\right]\mathrm{d}\mu _\mathrm{u}^\mathrm{e}\mathrm{d}ϵ_\mathrm{u}$$
(1)
$$\mu _\mathrm{d}^\mathrm{i}=\frac{\mu _\mathrm{u}^\mathrm{e}\beta _\mathrm{r}}{1\beta _\mathrm{r}\mu _\mathrm{u}^\mathrm{e}},ϵ_\mathrm{d}=\mathrm{\Gamma }_\mathrm{r}(1\beta _\mathrm{r}\mu _\mathrm{u}^\mathrm{e})ϵ_\mathrm{u},$$
(2)
and one obtains a similar system for shock crossing $`du`$:

$$f_\mathrm{u}^{2n}(ϵ_\mathrm{u},\mu _\mathrm{u}^i)\mathrm{d}\mu _\mathrm{u}^\mathrm{i}\mathrm{d}ϵ_\mathrm{u}=\left[_1^{\beta _{\mathrm{s}|\mathrm{d}}}d\stackrel{~}{\mu }_\mathrm{d}^\mathrm{i}𝒫_\mathrm{d}(\stackrel{~}{\mu }_\mathrm{d}^\mathrm{i};\mu _\mathrm{d}^\mathrm{e})f_\mathrm{d}^{2n1}(\stackrel{~}{ϵ}_\mathrm{d};\stackrel{~}{\mu }_\mathrm{d}^\mathrm{i})\right]\mathrm{d}\mu _\mathrm{d}^\mathrm{e}\mathrm{d}\stackrel{~}{ϵ}_\mathrm{d}$$
(3)
$$\mu _\mathrm{u}^\mathrm{i}=\frac{\mu _\mathrm{d}^\mathrm{e}+\beta _\mathrm{r}}{1+\beta _\mathrm{r}\mu _\mathrm{d}^\mathrm{e}},ϵ_\mathrm{u}=\mathrm{\Gamma }_\mathrm{r}(1+\beta _\mathrm{r}\mu _\mathrm{d}^\mathrm{e})\stackrel{~}{ϵ}_\mathrm{d},$$
(4)
The terms within brackets in Eqs. (1) and (3) correspond to the distribution function upon exit from upstream and downstream respectively. The return probability to the shock $`P_{\mathrm{ret}}(\mu _\mathrm{d}^\mathrm{i})`$ as a function of ingress pitch angle can be obtained as: $`P_{\mathrm{ret}}(\mu _\mathrm{d}^\mathrm{i})d\mu _\mathrm{d}^\mathrm{e}𝒫_\mathrm{d}(\mu _\mathrm{d}^\mathrm{i};\mu _\mathrm{d}^\mathrm{e})`$. The corresponding return probability for upstream is obviously unity. After each cycle $`udu`$, a fraction $`f_{\mathrm{out}}^{2n+1}(ϵ)=d\mu _\mathrm{d}^\mathrm{i}[1P_{\mathrm{ret}}(\mu _\mathrm{d}^\mathrm{i})]f_\mathrm{d}^{2n+1}(\mu _\mathrm{d}^\mathrm{i};ϵ)`$ of the particle population has escaped downstream and accumulates to form the outgoing particle population $`f_{\mathrm{out}}(ϵ)=_{n=0}^{n=+\mathrm{}}f_{\mathrm{out}}^{2n+1}(ϵ)`$. By following each shock crossing, and using Eqs. (1),(2),(3), (4) one can follow the evolution of $`f_\mathrm{d}`$, $`f_\mathrm{u}`$ and $`f_{\mathrm{out}}`$, starting from a mono-energetic and isotropic initial injection distribution upstream. The distribution $`f_{\mathrm{out}}(ϵ)`$ eventually provides the accelerated particle population. A similar formal development of the acceleration process by repeated shock crossings has been proposed independently by Vietri (2002): the flux of particles crossing the shock in the stationary regime, noted $`J_{\mathrm{in}}`$ in Vietri (2002) is related to the above as $`J_{\mathrm{in}}=C_{n=0}^{n=+\mathrm{}}f_\mathrm{d}^{2n+1}`$ with $`C`$ a normalization constant.

This method assumes that the angular probability laws do not depend on rigidity. This is true in the diffusive limit but one might expect some weak dependence to appear in the relativistic limit $`\mathrm{\Gamma }_\mathrm{s}1`$. Indeed we have found numerically such a weak dependence of $`𝒫_\mathrm{d}`$ and $`𝒫_\mathrm{u}`$ on the particle rigidity. However it remains weak, and the change in $`𝒫`$ averages to a few percent when the rigidity changes by two orders of magnitude. In terms of energy spectral index $`s`$, this dependence introduces an error of $`\delta s=\pm 0.02`$ for $`\mathrm{\Gamma }_\mathrm{s}=2`$ up to $`\delta s=\pm 0.04`$ for $`\mathrm{\Gamma }_\mathrm{s}=100`$. Therefore in the following we neglect the dependence on rigidity but keep the above errors as uncertainties on our results. Note that one can in principle incorporate this dependence on rigidity in our method at the expense of having to calculate the probability laws $`𝒫`$ for a wide range of values of the rigidity.

The present technique has significant advantages when compared to standard Monte-Carlo techniques which follow the particle trajectories on both sides of the shock through the whole acceleration process. Indeed, provided one neglects the dependence on rigidity of $`𝒫`$, one can simulate the trajectories of particles of high rigidity only (near the end of the resonance range) which are must faster to integrate than the trajectories of particles of small rigidity since the ratio of scattering time to Larmor time decreases with increasing rigidity. The direct Monte-Carlo methods also suffer from the problem of a small dynamic range of the magnetic fields as compared to the wide dynamic range of particles momenta. The present method also offers a significant gain in signal as will be obvious in the following. Finally our method differs from Gallant et al. (2000) as they use Itô differential techniques to simulate scattering downstream and analytical methods for scattering in a regular magnetic field upstream assuming $`\mathrm{\Gamma }_\mathrm{s}1`$. In contrast, the present simulations can be applied to any shock Lorentz factor and magnetic field configuration. Furthermore they use Monte-Carlo methods to simulate the acceleration process after constructing the laws of return probability while we directly fold over repeatedly the probability distributions in conjunction with the energy gain formula.

## 3 Results

The downstream return probability to the shock as a function of ingress pitch angle cosine is shown in Fig. 3 for $`\mathrm{\Gamma }_\mathrm{s}=2,100`$. The return probability for $`\mathrm{\Gamma }_\mathrm{s}=100`$ is an asymptote which is reached to within a percent as early as $`\mathrm{\Gamma }_\mathrm{s}=5`$. In Fig. 3, we show the average energy gain $`gϵ_f/ϵ_i`$ per cycle $`udu`$ (diamonds) and $`dud`$ (triangles) for $`\mathrm{\Gamma }_\mathrm{s}=100`$. This energy gain is defined as the ratio of the average energies at the end $`ϵ_f`$ and at the beginning of the cycle $`ϵ_i`$, with $`ϵd\mu dϵϵf(\mu ,ϵ)/d\mu dϵf(\mu ,ϵ)`$. The average energy gain in each cycle subsequent to first shock crossing is $`1.93`$ per cycle u–d–u for $`\mathrm{\Gamma }_\mathrm{s}=100`$, and this asymptotic value is reached immediately after the first cycle. This is a rather dramatic confirmation of the analytical expectations of Gallant & Achterberg (1999) and Achterberg et al. (2001) which had argued that only the first cycle should yield a gain $`\mathrm{\Gamma }_\mathrm{s}^2`$ since the anisotropy of the distribution function upstream is so pronounced that the gain in subsequent cycles is reduced to $`2`$.

An example of the spectrum of accelerated particles for $`\mathrm{\Gamma }_\mathrm{s}=100`$ after 20 cycles $`udu`$ is shown in Fig. 3; the thin lines in this figure show the fractions of particles $`f_{\mathrm{out}}^{2n+1}`$ that escape after $`2n+1`$ shock crossings. One clearly sees in this figure the piling up of populations of particles of ever decreasing size (due to finite escape probability) and ever increasing energy which gives rise to the accelerated population $`f_{\mathrm{out}}=_nf_{\mathrm{out}}^{2n+1}`$. The spectral index of the escaping population for $`\mathrm{\Gamma }_\mathrm{s}=100`$ is here $`s=2.26\pm 0.04`$ (incorporating the error due to dependence of $`𝒫`$ on rigidity), in excellent agreement with previous results by Bednarz & Ostrowksi (1998), Kirk et al. (2000) and Achterberg et al. (2001).

Finally, in Fig. 3, we give the average return probabilities (open squares), average asymptotic energy gains (diamonds) and spectral indices (filled circles) for values of $`\mathrm{\Gamma }_\mathrm{s}`$ comprised between 2 and 100. The average return probabilities shown in this figure are the return probabilities shown in Fig. 3 weighted by the corresponding asymptotic downstream ingress pitch angle distribution, i.e. $`P_{\mathrm{ret}}=lim_{n+\mathrm{}}𝒫_{\mathrm{ret}}\left(\mu _\mathrm{d}^\mathrm{i}\right)f_\mathrm{d}^n(\mu _\mathrm{d}^\mathrm{i})d\mu _\mathrm{d}^\mathrm{i}/f_\mathrm{d}^n(\mu _\mathrm{d}^\mathrm{i})d\mu _\mathrm{d}^\mathrm{i}`$. A naive unweighted average of the return probability shown in Fig. 3 for $`\mathrm{\Gamma }_\mathrm{s}=100`$ would give 0.33, whereas the weighted average gives 0.40: the difference is directly related to the strong anisotropy at shock crossing.

The standard non-relativistic formula for the (energy) spectral index $`s`$ (Bell 1978), $`s=1\mathrm{log}\left(P_{\mathrm{ret}}\right)/\mathrm{log}\left(g\right)`$, with $`g`$ the average energy gain, is in relatively good agreement with the slopes obtained, provided one uses the weighted average for the return probability as above, see Fig. 3. A more general formula which includes relativistic effects has been proposed by Vietri (2002): $`P_{\mathrm{ret}}g^{s1}=1`$. One can derive this formula and variants of it by using our Eqs. (1),(2),(3), (4). For instance, insert Eq. (3) into Eq. (1) then sum over $`n`$ (shock crossing number) to go to the stationary regime and consider an energy range where $`ϵϵ_0`$, $`ϵ_0`$ being the maximal injection energy. There one expects that $`_nf_\mathrm{d}^{2n+1}ϵ_\mathrm{d}^s\varphi (\mu _\mathrm{d})`$, i.e. the distribution factorizes out in an energy power law times a function of pitch angle, and indeed this property is verified numerically to high accuracy. Then one introduces the energy gain per cycle $`g(\mu _\mathrm{d}^\mathrm{e},\mu _\mathrm{d}^\mathrm{i})=ϵ_\mathrm{d}/\stackrel{~}{ϵ}_\mathrm{d}`$ and integrates over $`\mu _\mathrm{d}^i`$ both sides of the equation. Finally, dividing one side by the other yields the following relation, which is a variant of the formula of Vietri (2002):

$$_1^{\beta _{\mathrm{s}|\mathrm{d}}}d\mu _\mathrm{d}^\mathrm{i}_{\beta _{\mathrm{s}|\mathrm{d}}}^1d\mu _\mathrm{d}^\mathrm{e}\stackrel{~}{P}_{\mathrm{ret}}(\mu _\mathrm{d}^\mathrm{e})\overline{𝒫}_\mathrm{u}(\mu _\mathrm{d}^\mathrm{e},\mu _\mathrm{d}^\mathrm{i})g^{s1}(\mu _\mathrm{d}^\mathrm{e},\mu _\mathrm{d}^\mathrm{i})=1.$$
(5)
In this equation, $`\overline{𝒫}_\mathrm{u}(\mu _\mathrm{d}^\mathrm{e},\mu _\mathrm{d}^\mathrm{i})`$ simply corresponds to $`𝒫_\mathrm{u}`$ expressed in terms of downstream pitch angle cosines, and $`\stackrel{~}{P}_{\mathrm{ret}}(\mu _\mathrm{d}^\mathrm{e})d\mu _\mathrm{d}^\mathrm{i}𝒫_\mathrm{d}(\mu _\mathrm{d}^\mathrm{i};\mu _\mathrm{d}^\mathrm{e})`$. Equation (5) is indeed verified to within the numerical noise ($`<1`$%).

Our simulations also provide a direct measurement of the acceleration timescale, which can be defined as the cycle time in the upstream rest frame when $`\mathrm{\Gamma }_\mathrm{s}1`$: $`t_{\mathrm{acc}}(ϵ)t_{\mathrm{u}|\mathrm{u}}(ϵ)+\mathrm{\Gamma }_\mathrm{r}t_{\mathrm{d}|\mathrm{d}}(ϵ/\mathrm{\Gamma }_\mathrm{r})`$, where $`t_{\mathrm{u}|\mathrm{u}}(ϵ)`$ and $`t_{\mathrm{d}|\mathrm{d}}(ϵ)`$ denote the upstream and downstream return times measured in their respective rest frames for a particle of energy $`ϵ`$. For $`t_{\mathrm{d}|\mathrm{d}}`$ we find to within the noise of the simulations: $`t_{\mathrm{d}|\mathrm{d}}1.5\beta _{\mathrm{s}|\mathrm{d}}^1t_{\mathrm{scatt}|\mathrm{d}}`$, with $`t_{\mathrm{scatt}|d}`$ the scattering time downstream. The scattering time is given as a function of Larmor time $`t_\mathrm{L}`$ in Casse, Lemoine & Pelletier (2001), and for pure turbulence one finds: $`t_{\mathrm{scatt}}/t_\mathrm{L}0.4\rho ^{2/3}`$ for $`\rho 0.1`$ and $`t_{\mathrm{scatt}}/t_\mathrm{L}2`$ for $`0.1\rho 1`$. Here $`\rho k_{\mathrm{min}}r_\mathrm{L}`$ denotes the rigidity, with $`k_{\mathrm{min}}`$ the smallest wavenumber of the magnetic field modes. The above result for $`t_{\mathrm{d}|\mathrm{d}}`$ can be understood as follows. In the non-relativistic limit one derives $`t_{\mathrm{d}|\mathrm{d}}(2/3\beta _{\mathrm{s}|\mathrm{d}}^2)(1P_{\mathrm{ret}})/P_{\mathrm{ret}}`$ and we find $`P_{\mathrm{ret}}(1\beta _{\mathrm{s}|\mathrm{d}})/2`$ in the relativistic limit, which gives $`t_{\mathrm{d}|\mathrm{d}}t_{\mathrm{scatt}|\mathrm{d}}/\beta _{\mathrm{s}|\mathrm{d}}`$.

Gallant & Achterberg (1999) have conjectured $`t_{\mathrm{u}|\mathrm{u}}t_{\mathrm{L}|\mathrm{u}}/\mathrm{\Gamma }_\mathrm{s}`$ for $`\mathrm{\Gamma }_\mathrm{s}1`$, corresponding to deflection by an angle of order $`1/\mathrm{\Gamma }_\mathrm{s}`$ in a regular magnetic field. We confirm this result up to a small residual dependence on scattering time/rigidity: $`t_{\mathrm{u}|\mathrm{u}}5\mathrm{\Gamma }_\mathrm{s}^1\rho _\mathrm{u}^{0.15}t_{\mathrm{L}|\mathrm{u}}`$ for $`\mathrm{\Gamma }_\mathrm{s}5`$. In the mildly relativistic case, for $`\mathrm{\Gamma }_\mathrm{s}=2`$, we find $`t_\mathrm{u}t_{\mathrm{scatt}|\mathrm{u}}`$, which indicates that the particles have time to wander before being caught by the shock in this case.

The final acceleration time depends on how the magnetic field downstream $`B_\mathrm{d}`$ is related to that upstream $`B_\mathrm{u}`$ (Gallant & Achterberg 1999). If one assumes that $`B_\mathrm{d}=\mathrm{\Gamma }_BB_\mathrm{u}`$ with $`\mathrm{\Gamma }_B\mathrm{\Gamma }_\mathrm{s}`$, and the turbulence is isotropic downstream but the length scale has been contracted by a factor $`\mathrm{\Gamma }_B`$, i.e. $`k_{\mathrm{min}|\mathrm{d}}\mathrm{\Gamma }_Bk_{\mathrm{min}|\mathrm{u}}`$, one finds, for $`\rho _{\mathrm{d}|\mathrm{d}}1`$, $`t_{\mathrm{acc}}(5\mathrm{\Gamma }_\mathrm{s}^{0.15}+9)t_{\mathrm{L}|\mathrm{u}}/\mathrm{\Gamma }_\mathrm{s}`$. The acceleration timescale can thus be as short as a fraction of Larmor time in the ultra-relativistic limit.

This is a most interesting aspect of relativistic Fermi acceleration, as it implies that the particles can reach the energy confinement limit $`ϵ_{\mathrm{cl}}=Ze\mathrm{\Gamma }Br`$ when the acceleration is limited by expansion losses or by the age of the shock wave. Here, $`r=\mathrm{min}(l,t/\kappa c)`$, $`l`$ is the size of the accelerating region, $`t`$ the age or losses timescale and $`\kappa =t_{\mathrm{acc}}/t_\mathrm{L}<1`$, all quantities being expressed in the comoving frame. This is particularly relevant for the generation of ultra-high energy cosmic rays in relativistic winds such as $`\gamma `$ray bursts (Waxman 1995, Vietri 1995, Gallant & Achterberg 1999, Gialis & Pelletier 2003). In particular our results for $`\mathrm{\Gamma }_\mathrm{s}=2`$ are of direct relevance to the acceleration of protons and electrons in internal shocks of $`\gamma `$ray bursts, while the results for $`\mathrm{\Gamma }_\mathrm{s}1`$ can be applied directly to the external shock model of $`\gamma `$ray bursts.

To sum up, our simulations confirm the results of Bednarz & Ostrowski (1998) concerning the spectral index and those of Achterberg et al. (2001) concerning the spectral index, the acceleration time scale and energy gains. Very recently, Ellison & Double (2002) found a similar index by taking into account the backreaction of cosmic rays on the shock structure. In the near future, the present method will be applied to more general magnetic field configurations including parallel, transverse, subluminal and superluminal shocks.

We would like to thank Y. Gallant for fruitful discussions.


text:
# Reply to the comment by I.Ispolatov and M.Karttunen, cond-mat/0303564

The description of self-gravitating stellar systems as equilibrized objects in the microcanonical thermodynamics applied also to ”Small” systems gross174  with sharp total energy and angular momentum demands a spherical container centered in the center of mass of the total system. Of course in reality there is no such container. The container is necessary to avoid particle loss by evaporation. The statistical equilibrium thus obtained is a metastable equilibrium which lives over times such that the mass loss due to evaporation is negligible, and makes sense as long as nuclear hydrogen burning stabilizes the stars against collapse.

In order for the calculated equilibrium configuration to have anything in common with the real free stellar system the effect of the container on the obtained mass distribution must be minimized. I.e., only configurations with the center of mass fixed at the center of the container are admitted. Moreover, the density near the container walls should be small. This was achieved in our calculation gross187 ; gross190 ; gross191 ; gross195  by setting the contributions of the spherical harmonics $`b_{l=1,m}(x)x^2𝑑x=0`$.

In the molecular-dynamics simulation by Ispolatov and Karttunen this necessary constraint was forgotten. Therefore the spatial distribution they show in their Figure 1 is irrelevant for reality because the container has a non-negligible momentum and the center of mass of the gas is not at rest. Certainly after removing this crucial restriction it is impossible to find any realistic configuration like binary stars rotating around their common center of mass and with a negligible density of matter at the artificial container walls, as the one we describe e.g. in gross187 . On the other hand, if in our setting we release the constraint on the functions $`b_{1,m}`$ we get configurations similar to those they find.

A possible way to enforce preservation of center of mass within molecular dynamics may be to take only those initial conditions with positive parity, where for each particle at $`(r_i,p_i)`$ there is a partner starting at opposite initial point in phase space $`(r_i,p_i)`$. Such pairs of particles would hit the container simultaneously at opposite sides and the center of mass of the whole system without the container would not move. This should lead to symmetric binary configurations as the ones found by us in gross187 .