text:
# Analyze This! A cosmological constraint package for cmbeasy

## 1 Introduction

The wealth of recent precision measurements in cosmology  places stringent constraints on any model of the Universe. Typically, such a model is given in terms of a number of cosmological parameters. Numerical tools, such as Cmbfast , Camb  and Cmbeasy , permit to calculate the prediction of a given model for the observational data. While these tools are comparatively fast, scanning the parameter space for the most likely model and confidence regions can become a matter of time and computing power.

The cost of evaluating models on a n-dimensional grid in parameter space increases exponentially with the number of parameters. In contrast, the Markov Chain Monte Carlo (MCMC) method scales roughly linearly with the number of parameters . The MCMC method has already been used to constrain various models . A popular tool for setting up MCMC simulations is the Cosmo-Mc package  for the Camb code, an improved proposal distribution for the local Metropolis algorithm has been proposed in .

In this paper, we introduce the AnalyzeThis package111It is part of the cmbeasy v2.0 release. for Cmbeasy. It includes a parallel MCMC driver, as well as routines to calculate the likelihood of a model with respect to various data sets. We took special care in designing a step-proposal strategy that leads to fast convergence and mixing of the chains. This strategy is applied during the early stages of the simulation. As soon as the likelihood contour has been roughly explored, the adaptive step proposal freezes in. This ensures that the MCMC results are not contaminated by the adaptive steps. At the same time, the automated step optimization considerably improves performance and is rather convenient. The raw data files can be processed from within a graphical user interface (gui). Using the gui, one can marginalize, visualize and print one and two dimensional likelihood surfaces (see figure 1).

The plan of this paper is as follows: we describe the MCMC method and our implementation in section 2. A brief introduction to the software is given in section 3. We analyze the constraints on cosmological parameters from observational data sets in section 4. In section 5 we present our conclusions, while the format of the MCMC data files is defined in A.

## 2 Markov Chain Monte Carlo simulation

In the following, we will review the basic ideas of Markov Chain Monte Carlo simulation  and describe our implementation.

Suppose that $`\theta `$ denotes a vector of all model parameters, $`X`$ represents some observed data and $`L(X|\theta )`$ is the likelihood of observing $`X`$ given parameters $`\theta `$. Specifying a prior distribution $`P(\theta )`$ for the parameters, Bayes’ theorem yields the posterior distribution $`\pi (\theta |X)`$ of $`\theta `$ given the observed data $`X`$:

$$\pi (\theta |X)=\frac{P(\theta )L(X|\theta )}{P(\theta )L(X|\theta )𝑑\theta }.$$
(1)
Using $`\pi (\theta |X)`$, one can compute expectation values

$$E[f(\theta )|X]=\frac{f(\theta )\pi (\theta |X)𝑑\theta }{\pi (\theta |X)𝑑\theta },$$
(2)
as well as confidence levels.

The idea of the MCMC method is to directly draw samples from the posterior $`\pi (\theta |X)`$. The statistical properties of $`\pi (\theta |X)`$ may then be estimated using this sample.222For instance, one can infer the mean of parameter $`\theta ^{(i)}`$ from the sample of parameter values $`\theta _1,\mathrm{},\theta _n`$ (the ergodic average):
$$E[\theta ^{(i)}]=\frac{1}{n}\underset{j=1}{\overset{n}{}}\theta _j^{(i)}.$$
(3)   To accomplish the sampling of $`\pi (\theta |X)`$ one uses a Markov Chain, which is a stochastic process $`\{\theta _0,\theta _1,\mathrm{},\theta _n\}`$ where $`\theta _n`$ only depends on $`\theta _{n1}`$. The idea is to choose the next point in the chain based on the previous point such that $`\pi (\theta |X)`$ becomes the stationary distribution of the chain

$$\text{Dist}\{\theta _0,\mathrm{},\theta _n\}\pi (\theta |X)\text{ as }n\mathrm{}.$$
(4)
There are several methods to accomplish this. We will concentrate on the Metropolis algorithm  and its implementation in Cmbeasy.

### 2.1 The Metropolis Algorithm

The algorithm is defined as follows (for an illustration see Fig. 2):

1. Choose starting parameter vector $`\theta _0`$.
2. Compute the likelihood $`L_0(X|\theta _0)`$ of observing the experimental data given the parameters $`\theta _0`$.
3. Obtain a new parameter vector by sampling from a “proposal distribution” $`q(\theta _{i1},\theta _i)`$ (see section 2.3).
4. Compute the likelihood $`L_i(X|\theta _i)`$.
5. If $`L_i>L_{i1}`$ then save $`\theta _i`$ as new point in the chain (“take the step”) and go to (3).
6. If $`L_i<L_{i1}`$ then generate a random variable $`u`$ from $`[0,1]`$. If $`u<L_i/L_{i1}`$ take the step as in (5). If $`u>L_i/L_{i1}`$ then reject $`\theta _i`$, save $`\theta _{i1}`$ as new point in the chain and go to (3).

This algorithm assumes flat priors $`P(\theta )`$ and a symmetric proposal distribution $`q(\theta _{i1},\theta _i)=q(\theta _i,\theta _{i1})`$. Note that we assign likelihood zero to any parameter set that has at least one point outside its prior. In this version of the Metropolis algorithm, all parameters change with every step, a strategy called global Metropolis.111 One may also change one parameter (or a subset of parameters) at a time. Such a local Metropolis algorithm is implemented in Cosmo-Mc.

### 2.2 Convergence Testing

At the beginning, the chain migrates from its random starting point to regions of higher likelihood. Points during this “burn-in” do not constitute a sample from $`\pi (\theta |X)`$ and should be eliminated. In principle, it may be difficult to tell from a single chain if it has converged towards the underlying $`\pi (\theta |X)`$. In MCMC, one therefore uses several chains with random starting points and monitors mixing and convergence.

The convergence test of Gelman and Rubin  monitors the variance of a parameter between the chains. To be precise, consider using the last $`n`$ points of each of $`m`$ chains for the test. Let $`\psi _{ij}`$ label one entry of the parameter vector $`\theta `$ at point $`j=1,\mathrm{},n`$ in chain $`i`$ with $`\overline{\psi }_i`$ the mean for chain $`i`$ and $`\overline{\psi }`$ the mean of all chains. The variance between chains $`B`$ and the within-chain variance $`W`$ are then given by

$`B`$ $`=`$ $`{\displaystyle \frac{n}{m1}}{\displaystyle \underset{i=1}{\overset{m}{}}}(\overline{\psi }_i\overline{\psi })^2,`$ (5)
$`W`$ $`=`$ $`{\displaystyle \frac{1}{m}}{\displaystyle \underset{i=1}{\overset{m}{}}}s_i^2,\text{ where }s_i^2={\displaystyle \frac{1}{n1}}{\displaystyle \underset{j=1}{\overset{n}{}}}(\psi _{ij}\overline{\psi }_i)^2,`$ (6)
and the quantity

$$R=\frac{\frac{n1}{n}W+\frac{1}{n}B}{W}$$
(7)
should converge to one.333In a realistic situation, the numerator is an overestimate whereas the denominator is an underestimate of the variance of the stationary distribution of $`\psi `$. A value of $`R<1.2`$ for all parameters indicates that the chain is sampling from $`\pi (\theta |X)`$.444Indicating that the step sizes and directions will freeze in in our implementation. From this point onwards one may use the chain points for parameter estimates.

When do we have enough points for parameter estimation? This question is not easy to answer, since it is depends on the model, the used data sets and the desired accuracy how many chain points are needed for a robust estimate of parameters. Therefore, in our implementation the MCMC simulation runs indefinetly. However, any “breaking-criterion” may be implemented easily, and the chains may be monitored with external programs during the run.

### 2.3 Adaptive Step Size Gaussian Sampler

The number of steps needed for good convergence and mixing depends strongly on the step proposal distribution. If the proposed steps are too large, the algorithm will frequently reject steps, giving slow convergence of the chain. If, on the other hand, the proposed steps are too small, it will take a long time for the chain to explore the likelihood surface, resulting in slow mixing. In the ideal case the proposal distribution should be as close to the posterior distribution $`\pi (\theta |X)`$ as possible – which unfortunately is not known a priori. While a simple Gaussian proposal distribution with step sizes $`\sigma _k`$ is sufficient, it is not optimal in terms of computing costs if cosmological parameters are degenerate. A naive Gaussian sampler would move slowly along the degeneracy direction, unaware of any degeneracy (see figure 3).

Instead of using a naive Gaussian proposal distribution, we sample from a multivariate Gaussian distribution with covariance matrix estimated from the previous points in the chains. By taking into account the covariances among the parameters, we effectively approximate the likelihood contour in extent and orientation – the Gaussian samples are taken along the principal axis of the likelihood contour. Denoting $`\theta _{i1}\theta _i=\mathrm{u}_\mathrm{i}`$ for notational convenience, the proposal distribution we use for the steps is

$$q(\theta _{i1},\theta _i)N\mathrm{exp}\left[\frac{1}{2}\mathrm{u}_\mathrm{i}^\mathrm{T}\mathrm{S}^\mathbf{}1\mathrm{u}_\mathrm{i}\right],$$
(8)
where $`N=\left(2\pi \right)^{K/2}(\text{det}\mathrm{S}\mathbf{)}^{\mathbf{}1\mathbf{/}2}`$ and $`\mathrm{S}`$ is the covariance matrix

$$\mathrm{S}\mathbf{=}\mathbf{\left(}\begin{array}{cccc}𝝈_1^2& 𝝆_{12}& \mathbf{}& 𝝆_{1\mathrm{K}}\\ 𝝆_{21}& 𝝈_2^2& \mathbf{}& 𝝆_{2\mathrm{K}}\\ \mathbf{}& & \mathbf{}& \mathbf{}\\ 𝝆_{\mathrm{K}1}& \mathbf{}& 𝝆_{\mathrm{K}\mathbf{}1\mathrm{K}}& 𝝈_\mathrm{K}^2\end{array}\mathbf{\right)}\mathbf{.}$$
(9)
The sampling is most easily performed by diagonalizing the covariance matrix

$$\mathrm{T}^\mathrm{T}\mathrm{S}\mathrm{T}\mathbf{=}\mathrm{D}\mathbf{}\mathrm{T}^\mathrm{T}\mathrm{S}^\mathbf{}1\mathrm{T}\mathbf{=}\mathrm{D}^\mathbf{}1\mathbf{,}$$
(10)
where $`\mathrm{T}`$ is an orthogonal matrix. Using this, equation (8) becomes

$`q(\theta _{i1},\theta _i)N\mathrm{exp}\left[\mathrm{u}_\mathrm{i}^\mathrm{T}\mathrm{T}\mathrm{T}^\mathrm{T}\mathrm{S}^\mathbf{}1\mathrm{T}\mathrm{T}^\mathrm{T}\mathrm{u}_\mathrm{i}\right]`$ (11)
$`=N\mathrm{exp}\left[{\displaystyle \frac{1}{2}}\mathrm{v}_\mathrm{i}^\mathrm{T}\mathrm{D}^\mathbf{}1\mathrm{v}_\mathrm{i}\right],`$ (12)
where $`\mathrm{v}_\mathrm{i}\mathbf{}\mathrm{T}^\mathrm{T}\mathrm{u}_\mathrm{i}`$. Thus, the procedure for obtaining a sample $`\mathrm{u}_\mathrm{i}`$ is as follows:

1. Find the eigenvalues $`\stackrel{~}{\sigma }_j^2`$ and eigenvectors of $`\mathrm{S}`$. Construct the transformation matrix $`\mathrm{T}`$ from the eigenvectors.
2. Draw Gaussian samples with variances $`\stackrel{~}{\sigma }_j^2`$, thereby obtaining the vector $`\mathrm{v}_\mathrm{i}`$.
3. Then $`\mathrm{u}_\mathrm{i}\mathbf{=}\mathrm{T}\mathrm{v}_\mathrm{i}`$ is the desired sample from the multivariate Gaussian with covariance matrix $`\mathrm{S}`$.

The convergence can be further improved by scaling the covariance matrix $`\mathrm{S}`$ with a variable factor $`\alpha `$. Using $`\alpha `$, we can cope better in situations where the projected likelihood takes on banana shapes such as in . It also improves the convergence during the early stages when the low number of points available limits the estimate of the covariance matrix $`\mathrm{S}`$. We dynamically increase $`\alpha `$ if a chain takes steps too often,222Frequent acceptance means that the likelihood at the next step is roughly comparable to the current one. This happens when the chain rarely takes steps larger than $`1\sigma `$. while we decrease $`\alpha `$ if the acceptance rate is too low.555Rare acceptance means that the chain rarely explores points with the same likelihood, i.e. neighboring points within $`1\sigma `$. The positive effect of our scheme on the convergence is illustrated in figure 4.

One can show that modifying the proposal distribution based on previous chain data during the run may lead to a wrong stationary distribution $`\pi ^{}(\theta |X)`$. Therefore, we only apply the dynamical strategy of finding an optimal step proposal during the early stages of the simulation. When the convergence is better than $`R=1.2`$ and the chain has calculated a certain number of points, we freeze in the step proposal distribution. All points in the chain before this freeze in should be discarded (see also A).

## 3 The Software

The package is part of Cmbeasy and consists of two main components. The first one is a MCMC driver using LAM/MPI  for parallel execution of each chain. The second one is the AnalyzeThis class which is designed to evaluate the likelihood of a given model with respect to various data sets. These sets include the latest data of WMAP TT and TE , ACBAR , CBI , VSA  , 2dFGRS , SDSS , the SNe Ia compilations of Riess et al. , Tonry et al.  and Knop et. al.  as well as the IfA Deep Survey SNe Ia data . Data files for all experiments are included for convenience. New data sets are added continuously to the code.

### 3.1 The MCMC Driver

The example MCMC driver consists of two routines: master() and slave(). Using LAM/MPI for parallel computing, one master and up to ten slaves may be started. The master() will determine the initial random starting position for each chain. In a never ending loop, it then sends the parameters to the slave()’s and collects the results when the computation is finished. Whenever a step has been successful, it stores the parameters and likelihoods of the last step together with the number of times the chain remained at the same point in parameter space in a new line of one file per chain. The format of these raw MCMC data files is defined in A. The master() monitors convergence and mixing and determines the next step for the slave(). Before freeze in, the covariance is estimated and the step proposal accordingly modified. After freeze in, the proposal distribution remains unchanged.

### 3.2 The AnalyzeThis Class

The AnalyzeThis class provides several routines concerning CMB, SNe Ia and Large Scale Structure measurements. It also contains routines for marginalizing and plotting the Monte Carlo data.

#### 3.2.1 WMAP

The WMAP routines are (slightly modified) C++ ports of the likelihood code  available at the LAMBDA  web site. When the WMAP routines are called for the first time, the covariance matrices provided by WMAP will be converted to a binary format to speed up future use. A routine for WMAP normalization of the $`C_l`$ spectrum using the binned TT may be used instead of the old COBE normalization. (After the quick normalization, a best fit normalization may be called which uses the full likelihood routine provided by .)

#### 3.2.2 ACBAR, CBI and VSA

We use the procedures for likelihood computation described by the ACBAR, CBI and VSA collaborations , using window functions and calibration uncertainty. One can de-select data bins from each of these datasets. It is, for instance, possible to calculate the likelihood for the low $`l`$ data only, thus speeding up the computation of a model, because the window functions of the high $`l`$ data require multipoles up to $`l=4000`$. One may also wish to exclude $`l<800`$ data if one is using the WMAP observations in order to keep the sets independent.

#### 3.2.3 2dFGRS

For this Large Scale Structure dataset one compares the data with the theoretical power spectrum at $`z=0.17`$, the effective redshift of the survey, multiplied with the window function. We only include the region with $`k/h<0.15\text{ Mpc}^1`$ since at smaller scales nonlinear effects need to be taken into account . For these values, the bias is nearly constant . One may either specify the bias or marginalize over it.

#### 3.2.4 SDSS

The theoretical power spectrum to be compared with the data should be evaluated at the effective redshift of this survey, at $`z=0.1`$. We use the data given in Table 3 of  and the appropriate window functions. One may select the maximum $`k`$-value to be included in the likelihood estimate, but including data beyond $`k/h>0.15\text{ Mpc}^1`$ requires nonlinear corrections. Again, one can specify the bias or marginalize over it.

#### 3.2.5 Supernovae Ia

We include four routines for calculating the likelihood with respect to SNe Ia data. Please note that the sets of Riess et al., Tonry et al. and Knop et al. are not independent.

##### Riess et. al.

One can use the full dataset, subset selection of the “gold” set as described in  is possible. Likelihood computation as given in this paper.

##### Tonry et. al.:

From the supernovae compilation of Tonry et al.  one can use the full data set of 230 Supernovae or, alternatively, one may use a restricted set of 172 supernovae, where supernovae with $`z<0.01`$ and with excess reddening have been omitted as suggested in . In any case, we have taken the particular velocity uncertainty to be $`v=500\text{ km/s}`$ corresponding to $`\mathrm{\Delta }z=0.00167`$ and computed the corresponding uncertainty in the luminosity distance to obtain the likelihood:

$`\chi ^2={\displaystyle \underset{i=1}{\overset{N}{}}}{\displaystyle \frac{(\mathrm{log}d_l^{exp}\mathrm{log}d_l^{th})^2}{\sigma _{tot}^2}},`$ (13)
$`\text{with }\sigma _{tot}^2=\sigma _{exp}^2+\left({\displaystyle \frac{d\mathrm{log}d_l^{th}}{dz}}\right)^2\mathrm{\Delta }z^2.`$ (14)
Here, $`d_l^{exp}`$ and $`d_l^{th}`$ is the experimental and theoretical luminosity distance, respectively.

##### Knop et. al.:

The 54 SNe Ia presented in Tables 3-5  have been included in the code. We use the data with strech and extinction correction applied, subsample selection as discussed in  is also possible.

##### IfADS survey:

One may choose to include all 23 SNe Ia given in , or one may exclude the supernovae with excess reddening and those not unambigously identified as SNe Ia as described in .

### 3.3 The Graphical User Interface

The gui may be used to process the raw output files of the MCMC chains. After starting cmbeasy, the first step is to “distill” the chain data files. By distilling we mean the merging of the raw chain output files into one file. In addition, distilling removes all burn-in data. To get started immediately, we include raw data from an example MCMC run in the resources directory of cmbeasy. The four chains are called montecarlo\_chain$`\alpha `$.dat. They are runs for a $`\mathrm{\Lambda }`$CDM model and need to be distilled first. Two and one dimensional marginalized likelihoods may then be plotted and printed from within the gui (see figure 1). Please see the “howto-montecarlo” document shipped with cmbeasy for an introduction.

## 4 Cosmological constraints from the data sets

Having described the method and software, we can now proceed to investigate the impact of the different data sets on the distribution of parameters for a given cosmological model. For illustrative purposes we limit ourselves to a flat $`\mathrm{\Lambda }`$CDM cosmology with five parameters. We take the reduced baryon and matter densities $`\mathrm{\Omega }_bh^2`$ and $`\mathrm{\Omega }_mh^2`$, the Hubble parameter $`h=H_0/(100\text{km s}^1\text{Mpc}^1)`$, the optical depth to the last scattering surface $`\tau `$ and the spectral index of the initial power spectrum $`n_s`$ as parameters. We neglect any tensor contributions and we marginalize over the amplitude of the initial power spectrum. Thus, the amplitude is treated as a “nuisance” parameter that is integrated out.

In table 1 we display the cosmological parameters and the flat priors used.

The limits of the parameter ranges are chosen in accordance with previous results .

First, we will investigate constraints from CMB-only data, subsequently adding Supernovae Ia and large scale structure data sets. For a detailed description of how we implemented the datasets we refer the reader to section 3.2 of this paper. We have performed three MCMC simulations with CMB only, CMB+ SNe Ia and CMB + SNe Ia + LSS data sets. Convergence was reached in each run after about $`2000`$ iterations, each simulation was run until 55 000 models had been computed.

For the CMB-only run we use the measurements of WMAP (TT and TE spectra), CBI, VSA and ACBAR up to $`l=2000`$, removing data points from CBI, VSA and ACBAR where we have WMAP measurements with comparable error bars to keep the data sets independent. The entire set we used is displayed in figure 5.

For CMB + SNe Ia we add the data set of Riess et al.  (“gold sample”), and finally for CMB + SNe Ia + LSS we add the large scale structure measurements of the SDSS collaboration using all points with $`k/h<0.15\text{Mpc}^1`$ where the perturbations are still linear . The results are presented in table 2, marginalized likelihoods for the parameters are displayed in table 3 and figure 6.

We will now discuss the impact of the different sets on the parameters. Consider figure 6. The CMB-only data set does not well constrain $`\mathrm{\Omega }_mh^2`$. Adding supernovae improves the bounds on the total matter content considerably, because the lumiosity distance depends sensitively on $`\mathrm{\Omega }_m`$. Adding large scale structure data breaks the degeneracy between matter contribution and initial power spectrum amplitude, hence we obtain still tighter bounds. Note that the distribution of $`\mathrm{\Omega }_mh^2`$ obtained from CMB-only and CMB + SNe Ia + LSS data sets differ somewhat. This has already been noted in . Even from CMB measurements alone we can infer the presence of dark matter at high significance. In contrast, the baryon contribution $`\mathrm{\Omega }_bh^2`$ is already well constrained by the CMB alone, adding SNe Ia and LSS data only improves the bounds slightly. The CMB is most sensitive to this parameter since what one observes are essentially oscillations in a photon-baryon plasma, the density of baryons is a critical parameter for the shape of the CMB spectrum. SNe Ia and LSS measurements are more or less insensitive to this parameter, though hints at oscillations in the power spectrum have been detected .

The Hubble parameter $`h`$, is only slightly constrained by CMB measurements alone. Adding SNe Ia data considerably improves the bounds, adding large scale structure data does not improve the bounds very much. The bounds obtained for the Hubble parameter for the full data set are consistent with the HST key project value $`h=0.72\pm 0.08`$  derived from measuring the local Hubble flow. Even though the optical depth does not directly influence SNe Ia predictions, adding SNe Ia data tightens the bound on the optical depth $`\tau `$ as seen in figure 6. The mechanism is somewhat indirect: SNe Ia data tightens the bound on $`\mathrm{\Omega }_m`$ and thus limits the range of allowed values for $`\tau `$, $`h`$ and $`n_s`$ from CMB measurements. As we marginalize over the bias of the SDSS data, large scale structure does not add further to the bounds on $`\tau `$.

## 5 Conclusions

We have introduced the AnalyzeThis package, which can be used to constrain cosmological models using observational data sets. The AnalyzeThis class provides many functions to compute the likelihood of a model with respect to measurements of the CMB, SNe Ia and Large Scale Structure.888We would be happy to include contributions of routines for recent and future measurements. So if the readers favorite measurement is not included yet, please send in a few lines of code.

In order to constrain models of the Universe with a substantial number of parameters, we include a Markov Chain Monte Carlo driver. As the MCMC step strategy determines the convergence speed of the chains, we implemented a multivariate Gaussian sampler with an additional dynamical scaling. We stop the adaptive improvement of the step proposal density as soon as a good level of convergence is reached and discard all data calculated before. Using this approach, we combine the advantage of a single-run automatic optimization of the step strategy with the demand of a static step proposal density.

The output of the Monte Carlo chains is in a human readable format and may be processed by any software, even during the run. For convenience however, one may use Cmbeasy’s gui to marginalize, plot and print two and one dimensional likelihood surfaces.

Finally, we discussed the impact of the different data sets on the parameters for the case of a $`\mathrm{\Lambda }`$CDM cosmology with five parameters.

We would like to thank Robert R. Caldwell, Robert A. Knop, Havard B. Sandvik, Max Tegmark, Adam G. Riess and Timothy J. Pearson for helpful discussions. M. Doran is supported by NSF grant PHY-0099543. C.M. Müller is supported by GRK grant 216/3-02.

## Appendix A Files and Formats

### A.1 MCMC chain output data

The parameters and likelihoods of each model are stored in files with the naming convention “montecarlo\_chain$`𝜶`$.dat”, where “$`𝜶`$” is an integer. Each file corresponds to a slave() process (and hence to an independent chain). Each line in such a file represents a new point in parameter space. The format of one such line is:

$$p_0|p_1|\mathrm{}|p_n|L|l_0|\mathrm{}|l_m|a_1|\mathrm{}|a_k|M,$$

where $`p_i`$ are the parameters of the model, $`L`$ is the overall likelihood, the $`l_i`$ are likelihoods of different experiments, $`a_i`$ are auxiliary fields (to store $`\sigma _8`$, for instance). Finally, $`M`$ determines the weight of this parameter set, i.e. the “time” spend before leaving this point in parameter space. When the step sizes and the covariance estimates are frozen in, one line in each file gets a zero $`M`$. All data before this line should be regarded as “burn-in”. The graphical user interface (using the distillChains() routine of AnalyzeThis) will automatically discard all models before freeze-in.

### A.2 Monitoring Progress and Covariance

The master() routine outputs the convergence and some more information into a file called “progress.txt”. Until the step proposal distribution is frozen in, the covariance matrix is in addition output to “covMatrix.txt”. Messages of any errors occuring in slave() will be appended to the file “errorlog.txt”.


format: text

text:
# X-ray observations of Be/X-ray binaries in the SMC Based on observations with XMM-Newton, an ESA Science Mission with instruments and contributions directly funded by ESA Member states and the USA (NASA)

## 1 Introduction

X-ray observations of the Small Magellanic Cloud (SMC) with recent observatories like ASCA, BeppoSAX, ROSAT and RXTE revealed a large number of Be/X-ray binary systems. Most of them were discovered through the detection of X-ray pulsations (e.g. see summary of ASCA X-ray sources in the SMC Yokogawa et al. 2003), indicating the spin period of the neutron star which orbits an early type star. Optical follow-up observations are then required to determine the spectral type of the high mass star and in the SMC only Be stars were identified (with the only secure exception of the supergiant system SMC X-1), making the SMC very different to the Milky Way and the Large Magellanic Cloud in absolute numbers of high mass X-ray binaries with respect to the galaxy mass and also in their ratio of Be to supergiant OB systems.

Haberl & Sasaki (2000, hereafter HS00) searched for Be/X-ray binaries in the ROSAT catalogues of Haberl et al. (2000) and Sasaki et al. (2000) by correlating them with catalogues of optical H$`\alpha `$ emission line objects (Meyssonnier & Azzopardi 1993; Murphy & Bessell 2000). HS00 presented a list of 47 Be/X-ray binaries and candidates with X-ray positions typically better than 10″. Since the launch of Chandra and XMM-Newton new X-ray pulsars and Be/X-ray binary candidates were discovered (Sasaki et al. 2001, 2003; Macomb et al. 2003). Together with eleven pulsars detected by RXTE with very uncertain X-ray positions we know about sixty high mass X-ray binaries (HMXBs) and candidates in the SMC, from which more than thirty are known as X-ray pulsars (Yokogawa et al. 2003; Mereghetti 2001).

The position uncertainties of X-ray sources in the SMC obtained by previous missions of up to tens of arc-seconds (including the most precise ROSAT positions) made it in many cases impossible to uniquely correlate sources from the different X-ray missions, nor allowed exact optical identifications. We therefore started to systematically analyse archival XMM-Newton data of HMXBs and candidates in the SMC in order to confirm their proposed nature. First results were presented in Sasaki et al. (2003). Here we continue this work analyzing the XMM-Newton EPIC data of two fields in the SMC which contain in total fifteen known HMXBs and candidates.

## 2 XMM-Newton observations and analysis

XMM-Newton observed the fields around the SMC supernova remnants SNR 0047$-$73.5 and SNR 0103$-$72.6 on October 15, 2000 and on October 17, 2000 with EPIC-MOS (MOS1 and MOS2, Turner et al. 2001) for 26.9 ks and 34.8 ks and EPIC-pn (Strüder et al. 2001) for 23.0 ks and 30.9 ks. For both observations the three EPIC cameras were operated in imaging mode covering the full field of view of $``$14′ radius, MOS in full-frame (2.6 s time resolution) and pn in extended full-frame (0.2 s) CCD readout mode. For optical light blocking the medium filter was used in all cameras.

The data were processed using the XMM-Newton analysis package SAS version 5.4.1 to produce the photon event files and binned data products like images, spectra and light curves. To avoid artificial photon arrival time jumps in the EPIC pn data an updated SAS development version from 2003, May 16 was used to create pn event files.

The observation of SNR 0103$-$72.6 suffered from strong background flaring activity (Fig. 1) which hampered the analysis of the data. We applied two different levels of background screening depending on the goals of the analysis. To obtain clean images which also show faint sources we used only data from times of low background ($`<`$10 cts ks-1 arcmin-2 in EPIC-pn). To search for periodicities in the brighter sources we used two different background cuts of 10 and 50 cts ks-1 arcmin-2. While the strong screening removes the flares which have typical time scales of a few 1000 s, the latter method avoids too many data gaps. For the spectral analysis of point sources the cut level of 50 cts ks-1 arcmin-2 was applied. For the MOS data corresponding levels of 2 and 5 cts ks-1 arcmin-2 were defined, which result in similar good time intervals in the contemporaneous parts of MOS and pn data.

Fig. 2 shows the EPIC-pn images of the two observed SMC fields. The investigated Be/X-ray binaries are marked with source numbers given in Table 3 which we use throughout the paper. The fields of view of the MOS cameras are slightly shifted with respect to pn and do not cover the 345 s pulsar SAX J0103.2$-$7209 (source 44).

### 2.1 Source positions

Source detection based on sliding window and maximum likelihood methods available in the SAS package was applied to the MOS and pn images, simultaneously for four energy bands (0.3-1.0 keV, 1.0-2.0 keV, 2.0-4.5 keV and 4.5-7.5 keV), but separately for the three instruments. Dividing the data into several energy bands enhances the signal to noise ratio for sources which have their major part of emission in a relatively narrow band, like e.g. super-soft X-ray sources or highly absorbed sources. In Table 1 we list the detected Be/X-ray binaries with XMM-Newton names, X-ray positions, count rates as obtained from the pn images (except for source 5 which is located near a CCD border in the pn camera). Coordinate uncertainties include a 4″ systematic error (Barcons et al. 2002). We correlated our source list with the compilation of Be/X-ray binaries from HS00 and with the H$`\alpha `$ emission line object catalogue of Meyssonnier & Azzopardi (1993). Source 47 (Fig. 2) detected in ROSAT HRI data (Sasaki et al. 2000) was not detected by the EPIC instruments. Also source 4 (Yokogawa et al. 2003) was not detected. The association of the 9.1 s pulsar AX J0049$``$732 (source 8) with either of source 7 or 9 (Yokogawa et al. 2003; Filipović et al. 2000b) could not be confirmed by a detection of the pulsations from either of the two sources in the EPIC data (see below). AX J0049$``$732 may be an uncorrelated source detected during a brighter state during the ASCA observations. HS00 classified RX J0105.7$``$7226 (\[HS00\]54) as candidate Be/X-ray binary due to the positional coincidence with the emission line object 1544 in Meyssonnier & Azzopardi (1993). Our newly derived X-ray position fully agrees with the ROSAT position, but the smaller uncertainty now excludes an association with the emission line object. Since no other evidence argues for RX 0105.7-7226 as Be/X-ray binary, we remove it from the list of candidates. Similarly, source 48 is now very unlikely to be associated with \[MA93\]1517. This makes the situation concerning the identification of the possible pulsar AX J0105$``$722 with a ROSAT source more unclear. However, the spectral analysis of RX J0105.1$``$7211 (see below) supports the HMXB nature of the source and we therefore keep it as a candidate.

There remain eleven Be/X-ray binaries and candidates in the two SMC fields for which we present a more detailed temporal and spectral analysis in the following.

### 2.2 X-ray pulsations

We searched for pulsations in the X-ray light curves (from different energy bands and background screening levels) of all eleven Be/X-ray binaries and candidates using Fast Fourier Transform (FFT) and Rayleigh Z$`{}_{}{}^{2}{}_{n}{}^{}`$ methods. The power spectra obtained from broad band data (0.3$``$7.5 keV) of the three known pulsars in the field of SNR 0047$-$73.5 (sources 11, 12 and 17) clearly reveal their pulse periods. This - together with the accurate source positions \- now allows the unique identification of the ASCA sources (which showed pulsations) with ROSAT sources (good X-ray positions), which could previously be done on positional coincidence only. In particular RX J0051.9$``$7311 (source 17), which was suggested to be the counterpart of the 16.6 s pulsar XTE J0050$``$732#1 (Lamb et al. 2002b) or the 172 s pulsar AX J0051.6$``$7311 (Yokogawa et al. 2003), is consistent in position with XMMU J005152.2$``$731033, which shows 172.2 s pulsations. In addition we find 263 s pulsations from XMMU J004723.7$``$731226=RX J0047.3$``$7312 which is likely associated with AX J0047.3$``$7312 (see also Yokogawa et al. 2003). In Fig. 3 we present FFT power spectra and light curves folded at the pulse period of the four pulsars. As indicated by the presence of the first harmonic in the power spectrum of RX J0049.7$``$7323, the pulse profile of this long period pulsar is double peaked. The determined pulse periods with 1$`\sigma `$ error are listed in Table 1.

From SAX J0103.2$-$7209 (source 44, Israel et al. 2000) which is located at the rim of the field of SNR 0103$-$72.6, we detect pulsations with a period of 341.2 s. This indicates that the period of large spin-up reported by (Israel et al. 2000) is now lasting more than 4.5 years (Figs. 4 and 5). In Fig. 5 we also plot period measurements of the two long period pulsars RX J0050.8$``$7316 and RX J0049.7$``$7323. The two pulsars with periods measured at more than two epochs, SAX J0103.2$-$7209 and RX J0050.8$``$7316=AX J0051$``$733, show very similar pulse periods and also (linear) period changes $`\dot{P}`$ of $``$1.6 s yr-1 and $``$1.0 s yr-1, respectively. A $`\dot{P}`$ of -14 s yr-1 is derived for the 750 s pulsar XMMU 004942.3$``$732313 but is based on only two period measurements, one of them with large error.

For six of the eleven sources we do not detect significant pulsations. From a Z$`{}_{}{}^{2}{}_{1}{}^{}`$ analysis no period with more than 2.5 $`\sigma `$ significance was found. In particular in the data of observation 0110000101 with low background, modulations with a similar pulsed fraction as seen from the detected pulsars (Fig. 3) should be easily picked up. On the other hand, the sensitivity for period detection is lower in the high background observation 0110000301. The pulse period of SAX J0103.2$-$7209 was detected with a formal 3.2 $`\sigma `$ significance (for an unknown period) and a slightly weaker modulation would not be detected for sources 43 and 48.

### 2.3 X-ray spectra

EPIC-pn and -MOS spectra were extracted for the eleven Be/X-ray binaries and candidates and simultaneously fit. An absorbed power-law model describes the data sufficiently well (with typical reduced $`\chi ^2`$ values around 1.0) for the given statistical quality of the spectra. Normalization factors were allowed to vary between the individual instruments to account for cross-calibration effects. Fig. 6 shows an example fit to the spectra of XMMU J004942.3$``$732313 (source 11), the 750 s pulsar located in the field of SNR 0047$-$73.5. Table 2 summarizes the spectral fit results for all eleven sources. The power-law photon indices $`\gamma `$ show a narrow distribution with an average value of 1.0 and a standard deviation of 0.16. The total absorption column densities are above the Galactic foreground value of 5.74$`\times 10^{20}`$ cm-2 as observed towards the SMC (Dickey & Lockman 1990) and indicates absorption by interstellar matter in the SMC and/or local matter around the X-ray sources. For the absorbing material solar element abundances were used in the spectral fits. If an SMC abundance of 0.2 solar is assumed (Russell & Dopita 1992), the contribution to the column density by SMC matter (total NH reduced by the Galactic foreground) needs to be multiplied by a factor of 5.0. In Table 2 we list observed fluxes and source intrinsic luminosities in the 0.2-10.0 keV energy band.

### 2.4 Long term variability

To investigate intensity variations on time scales of years we list in Table 2 the range of X-ray luminosities in which the investigated sources were detected in the past. We do this separately for instruments covering an energy band similar to the EPIC band (ASCA, BeppoSAX) and for ROSAT with the soft band of 0.1-2.4 keV. Most broad band luminosities were reported from ASCA observations by Yokogawa et al. (2003, 0.7-10.0 keV) which were corrected by a factor of 1.05 to account for the slightly larger (0.2-10.0 keV) energy band. Only source 44 (SAX J0103.2$-$7209) was observed with maximum luminosity by BeppoSAX (Israel et al. 2000). In this case the correction factor is 1.23. The correction factors were derived for an unabsorbed power-law spectrum with photon index of 1.0. Source 5 was not detected by ASCA and for sources 7 and 9 the association with the nearby ASCA source 8 is unclear. We also do not give an ASCA flux for source 48 because of contamination of the spectrum by nearby sources due to the large ASCA point spread function. To obtain 0.2-10 keV luminosities based on ROSAT measurements we used the spectral parameters inferred from the EPIC spectra, computed expected PSPC count rates by folding the spectral model through the PSPC detector response and scaled the observed count rates to derive the conversion factor. Measured PSPC count rates were taken from the Results Archive (RRA) PSPC source catalogue from pointed ROSAT observations (ROSAT Consortium 2000) which contains multiple detection entries for a source if observed more than once during the ROSAT mission. The last column in Table 2 shows the ratio between maximum and minimum luminosity which is at most a factor of 13. The intensity variations are relatively moderate for this kind of objects and can probably be explained by variations of the mass accretion rate along a neutron star orbit with moderate eccentricity within the circum-stellar matter of the Be star. While variations up to a factor of $``$20 are observed from many Be/X-ray binaries (e.g. Haberl & Pietsch 1999), ‘giant’ outbursts (factor 100-1000) are rare and are thought to be explained by increased mass outflow through the Be stellar wind disc (Apparao 1994). No giant outbursts were so far detected from any of the eleven systems.

## 3 The HMXB population of the SMC

The detection of a periodic 263.6 s X-ray flux modulation from XMMU J004723.7$``$731226 increases the number of likely HMXB pulsars in the SMC to at least thirty-seven. In Table 3 we give an updated list of HMXBs and candidates in the SMC which is based on the work of HS00 and add transient pulsars detected by RXTE. We do not include two additional known pulsars which may be of different type: AX J0043$``$737 for which confirmation of a 87 ms period is needed and which could be a Crab-like pulsar (Yokogawa & Koyama 2000; Yokogawa et al. 2003), and CXOU J0110043.1$``$721134 a possible anomalous X-ray pulsar in the SMC (Lamb et al. 2002a).

The number of Be/X-ray binaries and candidates is the highest known from any galaxy, including the Milky Way where about forty such systems (Liu et al. 2000) are known today. The SMC Be/X-ray binaries, all at a similar distance, constitute an ideal sample for statistical studies. Luminosity distributions were derived by HS00 and Sasaki et al. (2003). Compared to the Milky Way the SMC harbours a higher fraction of low-luminosity systems ($``$$`10^{35}`$ erg s-1), which may be a selection effect due to a more difficult identification of such objects in the Galactic Plane (HS00). The XMM-Newton results from SMC HMXBs allow to compare the X-ray spectra from already 16 sources \[11 from this work and sources 33, 34, 35, 41 and 45 from Sasaki et al. (2003)\] obtained with the same instruments and processed in a coherent way. In Fig. 7 (upper panel) the distribution of the power-law photon index is presented. The distribution is strongly peaked at 1.0, but includes a weak wide contribution with indices between 0.65 and 1.45.

Although the column densities derived from the spectra have relatively large errors, the NH distribution (Fig. 7 bottom panel) exhibits two maxima. An excess of objects with little or no absorption may either indicate spectra with a soft component which is not recognized because of insufficient statistics or a special spatial distribution of the sources on the near side of the SMC (star formation arm). The second peak may be explained by interstellar absorption in the SMC with the majority of objects located in the inner parts of the SMC. A comparison with similar data obtained from background AGN behind the SMC should allow to derive constraints on the total absorbing column density through the SMC on a statistical basis.

## 4 Summary

We performed temporal and spectral analyses of eleven Be/X-ray binaries and candidates observed serendipitously by the EPIC instruments on board XMM-Newton in two observations of the supernova remnants SNR 0047$-$73.5 and SNR 0103$-$72.6. From five of the objects pulsations in their X-ray flux are detected, four of them were previously known pulsars and in one case, XMMU J004723.7$``$731227, the XMM-Newton data revealed 263 s pulsations for the first time. XMMU J004723.7$``$731227 coincides in position with RX J0047.3$``$7312 which was suggested to be a Be/X-ray binary due to its likely identification with an H$`\alpha `$ emission line object (HS00). The detection of the pulse period strongly supports this identification.

The detection of pulsations, together with the good X-ray positions allowed to securely identify sources detected by various X-ray space missions, like e.g. the 172 s pulsar AX J0051.6$``$7311 with XMMU J005152.2$``$731033 (via pulsations) and RX J0051.9$``$7311 (via accurate position). The latter was previously also suggested to be the counterpart of XTE J0050$``$732#1, a 16.6 s pulsar. The 345 s pulsar SAX J0103.2$-$7209 was found to continue its large spin down with a rate of $``$1.6 s yr-1 extending now over 4.5 years. The 323 s pulsar AX J0051$``$733 was found to exhibit very similar properties as SAX J0103.2$-$7209 with respect to spin period, period derivative and X-ray luminosity.

The large number of Be/X-ray binaries in the SMC allows first statistical studies of their spectral and temporal properties. XMM-Newton and Chandra are ideally suited to detect pulsars among the Be-/X-ray binary candidates at low source fluxes and future observations are required to further increase the sample and allow to analyse it in a homogeneous way.

###### Acknowledgements.

The XMM-Newton project is supported by the Bundesministerium für Bildung und Forschung / Deutsches Zentrum für Luft- und Raumfahrt (BMBF / DLR), the Max-Planck-Gesellschaft and the Heidenhain-Stiftung.


text:
# On the definition of temperature in FPU systems.

## 1 Introduction.

The Fermi–Pasta–Ulam system consists of a chain of $`N`$ nonlinear oscillators with certain given boundary conditions, tipically fixed ends. It is well known (see  and ) that, for energies below a certain threshold $`E_\mathrm{c}`$, if the energy is initially given to a few low frequency modes, equipartition of energy among the modes is eventually attained only after an extremely long time, while at intermediate times a kind of metaequilibrium state is attained, in which the energy is shared essentially within a packet of low frequency modes.

An interesting and much discussed problem, is whether the specific energy threshold $`E_\mathrm{c}/N`$ vanishes or not in the limit of an infinite number of particles. Here we leave this problem aside: we will suppose for example that the number $`N`$ of particles be fixed, so that the threshold certainly exists. We address instead the problem raised by the fact that below the threshold one meets, as in the theory of glasses, with time scales to thermal equilibrium which are very long, even longer than any fixed observational time scale. Does this lack of thermalization have any consequences on the relevant thermodynamical quantities? or even, is it possible to correctly (i.e. uniquely) define the quantities of interest? In other words, is it still possible to have a thermodynamics below the threshold?

In the literature, the discussion is usually focused on the specific heat, because from heuristic arguments it is suggested that to less chaotic motions there correspond smaller specific heats, with eventually zero specific heat for totally ordered motions (i.e. for integrable systems). Thus it is expected that by lowering the energy below threshold the specific heat should diminish; such a property, in turn, should be considered as a good indicator of the weakening of chaos.

The papers  and , which aim at evaluating the specific heat of Fermi–Pasta–Ulam systems below threshold by numerical computations, reach two opposite conclusions: in  the value of the specific heat remains constant (as would follow from the equipartition principle) even below the threshold, while in  the specific heat indeed begins to fall down, below the threshold, and finally vanishes as temperature goes to zero. This striking difference is apparently due to the different methods used in the two papers in estimating the specific heat. Actually, in both papers the Fermi–Pasta–Ulam system is kept isolated (fixed ends), so that a direct measurement of the specific heat is precluded (because a direct measurement requires at least one heat bath). The specific heat is then estimated from the fluctuations of energy of a subsystem through the well known relation between specific heat and mean square deviation of energy, which holds in the canonical ensemble. The two papers differ in the choice of the subsytem: in the paper , one considers the energy fluctuations of a small piece of the total chain, while in  one considers the energy fluctuations of a small packet of nearby modes. As the energy of each mode remains nearly constant below the threshold while the energy of a piece of chains still present large fluctuations, this indeed explains why the two papers reach opposite conclusions. Now, at most only one of the conclusions can be correct, if a right conclusion does exist at all; indeed it is not clear whether the specific heat can be defined in an unambiguous way below the threshold.

A different approach was proposed in the paper . In short, the basic remark is the following one. As the above mentioned relation between specific heat and energy fluctuations is obtained from the equilibrium Gibbs ensemble, then its validity below the threshold is in doubt just because, up to the considered times, the system has not yet thermalized. So, it is argued that in measuring the specific heat one should revert to the conventional method which makes use of a heat bath at a given temperature $`T`$ coupled to the Fermi–Pasta–Ulam system, with the corresponding familiar calorimetric expression for the specific heat. Namely, the energy exchange $`\mathrm{\Delta }Q`$ between the bath and the Fermi–Pasta–Ulam system is measured when the temperature is varied by $`\mathrm{\Delta }T`$, then the ratio $`\mathrm{\Delta }Q/\mathrm{\Delta }T`$ is computed, and (in principle) the limit is taken for vanishing $`\mathrm{\Delta }T`$. However, even with such a method, one is still confronted with a delicate problem, because the amount of exchanged energy $`\mathrm{\Delta }Q`$ does depend on how much time one has waited in making the measurement (this is the so–called waiting time of the theory of glasses). The curve predicted by equipartition is recovered for infinitely long waiting times, while for finite waiting times the specific heat is expected to vanish at sufficiently low temperatures.

The question is thus: does there exist a natural way to choose a definite waiting time? Another problem then arises, due to the fact that, in the ratio $`\mathrm{\Delta }Q/\mathrm{\Delta }T`$ defining the specific heat, one should insert in the denominator the variation of temperature of the Fermi–Pasta–Ulam system and not that of the bath. The question is then whether the temperature of the Fermi–Pasta–Ulam system is the same as that of the heat bath. The very fact that the quantity $`\mathrm{\Delta }Q`$ depends on the waiting time actually shows that this is not the case, just because the equality of the two temperatures would imply $`\mathrm{\Delta }Q=0`$. On the other hand, if one were able to identify the temperature of the Fermi–Pasta–Ulam system, then the question of the waiting time would have a quick answer: one should wait until the heat bath temperature and that of the Fermi–Pasta–Ulam system have become equal, and only at that time should one measure the corresponding energy exchange. So the possibility of having available well defined thermodynamic quantities on short time scales is based on the possibility of providing a good notion of temperature for the FPU system before complete equipartition be achieved. From this point of view, the identification usually made of the temperature of the Fermi–Pasta–Ulam system with the mean kinetic energy of its particles is not the correct one, because in such a case the temperatures of the two systems (Fermi–Pasta–Ulam system and heat bath) would remain different for extremely long times.

The identification of the mean kinetic energy with temperature is so deeply rooted in our minds, that the existence of another quantity playing that role seems hardly conceivable. The aim of this paper is to show instead that this is possible.

In Section 2 we give a preliminary discussion of the zeroth law for states of metaequilibrium such as those of the Fermi–Pasta–Ulam system below threshold, in Section 3 we describe the model we employ for measuring of the temperature of the Fermi–Pasta–Ulam system through heat baths by numerical computations, and the numerical results are given in Section 4.

## 2 Zeroth law and temperature in states of metaequilibrium.

One of the basic features of thermal equilibrium is the so–called zeroth law, which essentially amounts to the transitivity of the equilibrium. From this follows (see ) that for any system there exists a function of its macroscopic (the so called *empirical temperature*) which has the same value for bodies in equilibrium.

However, it is not granted that, for a given macroscopic state, the equilibrium is unique if some of the internal degrees of freedom are dynamically frozen. We are thinking typically of the case of polyatomic molecules (see ) for which it is known that the exchanges of energy between the center of mass and the internal degrees of freedom are so slow that the number of effective degrees of freedom depends on the time of observation. This is actually the general situation that occurs in states of metaequilibrium.

We thus address the problem whether it is possible to have a zeroth principle, and so also an empirical temperature, in situations of metaequilibrium, in which the physical quantities are changing only on a very long time scale. So, if we put our Fermi–Pasta–Ulam system in heat contact with another body, and observe that at first there is a rapid relaxation to a certain state, while a later evolution to a final equilibrium would take place on a time scale much longer than our observational scale, we can think of our system as if it were equilibrated. Obviously, one is not granted that in such a situation the zeroth law, i.e. the transitivity of this metaequilibrium state, holds. But, if this is the case, an empirical temperature can be defined. In other terms, if we put the Fermi–Pasta–Ulam system in contact with a thermometric body which, after a short transient, appears to have reached a temperature $`T`$ (not evolving on our time scale), and if later, after having put the Fermi–Pasta–Ulam system in contact with another body at the same temperature $`T`$, nothing seems to occur (i.e. there is no exchange of energy in mean between the bodies), then the metastable state reached does have the transitive property, and we are authorized to assign to the Fermi–Pasta–Ulam system the temperature $`T`$ reached by the thermometric body.

An equivalent arrangement, which we have actually implemented in our numerical simulations to be described below, is the following one: the Fermi–Pasta–Ulam system is put at the same time in contact both with a heat bath and with a thermometric body (obviously, with no direct connection between the two external bodies). In such a situation, the transitive property reduces to the property that in a short time the thermometer attains the same temperature of the heath bath; the subsequent evolution to a final equilibrium should take place later, at a much smaller rate. Then, *by definition*, the temperature of the Fermi–Pasta–Ulam system in the metaequilibrium state is the one reached by the thermometer after the short–time relaxation.

This definition can appear satisfactory from an operative point of view. However, as it stands, it still lacks a clear connection with the properties of the Fermi–Pasta–Ulam system itself. Indeed there remains the problem of understanding, how the zeroth law can hold even if the Fermi–Pasta–Ulam system did not yet thermalize. In this connection, we make reference to a known phenomenon  concerning the isolated Fermi–Pasta–Ulam system, namely the fact that, below threshold, for initial excitations of the low frequency modes the energy turns out to remain confined to modes below a certain critical mode $`k_{\mathrm{cr}}`$, while the higher modes are not significantly involved in the dynamics. Notice, moreover, that an analogous phenomenon, i.e. a dynamical involvement restricted to the modes of sufficiently low frequency, is know to occur also when a system is coupled to an external body, for the case of polyatomic molecules (see ). So, it is known that, on a short time scale, the high frequencies modes (above $`k_{\mathrm{cr}}`$) do not get dynamically involved, neither by the internal nonlinarities nor by an interaction with external bodies. In both cases a packet of low frequency modes is formed which are in mutual equilibrium, and moreover are active in the process of thermalization with external bodies. Only after a much larger time scale there follows a relaxation of the Fermi–Pasta–Ulam system to the true equilibrium state. Before, it appears as if there existed an adiabatic partition (of a dynamical nature) between low and high frequency modes.

If this is the correct picture, it is clear to what property of the Fermi–Pasta–Ulam system should our definition of temperature correspond: namely, to the mean energy of each of the low frequency modes (those below $`k_{\mathrm{cr}}`$). In the rest of the paper we will illustrate the results of some numerical computations, which, in our opinion, strongly support the fact that such a metaequilibrium is transitive, and that the empirical temperature thus defined coincides with the energy of the low frequency modes.

## 3 The model.

As mentioned above, the model is constituted by a Fermi–Pasta–Ulam system in contact with two bodies, acting, the one as a thermometer and the other one as a heat bath. We make the simplest choice, in which the two bodies are perfect gases. Each gas is modeled as a system of point particles having no interactions among them, while interacting with the Fermi–Pasta–Ulam system through some smooth force between each of the gas particles and one of the edge Fermi–Pasta–Ulam particles. In more detail, concerning the Fermi–Pasta–Ulam system we denote as usual by $`x_i`$, $`i=1,\mathrm{},N`$, the distance of the $`i`$–th particle from its equilibrium position, by $`p_i`$ its conjugate momentum, and consider the familiar “$`\beta `$–model” Hamiltonian

$$H_{\mathrm{FPU}}=\underset{i=1}{\overset{N}{}}\frac{p_i^2}{2m}+\underset{i=0}{\overset{N}{}}\frac{\mathrm{\Omega }^2}{2}(x_{i+1}x_i)^2+\frac{\beta }{4}(x_{i+1}x_i)^4,$$
(1)
involving two positive parameters $`\beta `$ and $`\mathrm{\Omega }`$, with fixed boundary conditions $`x_0=0`$, $`x_{N+1}=0`$. Now, to the ‘’left” of the Fermi–Pasta–Ulam chain we place a perfect gas which acts as a heat bath: denoting by $`y_i`$, $`\pi _i`$ the $`i`$–th gas particle’s position and momentum respectively, we have $`L<y_i<x_1`$, with $`L>0`$ playing the role of the volume of the gas. The motion of each particle is thus free apart from the fact that it suffers an elastic reflection as $`y_i=L`$, and that it moreover interacts with the first Fermi–Pasta–Ulam particle $`x_1`$ through a short range potential, which we choose as

$$V=V_0\frac{e^{(y_ix_1)/l_0}}{(y_ix_1)/l_0},$$

$`l_0`$ and $`V_0`$ denoting its range and strength respectively. In agreement with the bound given above, due to the singularity of the potential at $`y_i=x_1`$, the solutions $`y_i(t)`$ of the equations of motion cannot cross the point $`x_1(t)`$, i.e. for all times $`t`$ one has $`L<y_i(t)<x_1(t)`$. The Hamiltonian of the heat bath is thus

$$H_\mathrm{B}=\underset{i=1}{\overset{N}{}}\frac{\pi _i^2}{2m}+V_0l_0\frac{e^{(y_ix_1)/l_0}}{(y_ix_1)},$$
(2)
supplemented by the boundary condition that the particles are reflected at $`y_i=L`$.

To the ‘’right” of the Fermi–Pasta–Ulam chain we place the thermometer, which is taken again as a perfect gas, with Hamiltonian

$$H_\mathrm{T}=\underset{i=1}{\overset{N}{}}\frac{\stackrel{~}{\pi }_i^2}{2m}+V_0l_0\frac{e^{(\stackrel{~}{y}_ix_N)/l_0}}{(\stackrel{~}{y}_ix_N)},$$
(3)
(plus a reflection condition at $`\stackrel{~}{y}_i=L`$) where $`\stackrel{~}{y}_i`$, $`\stackrel{~}{\pi }_i`$ are the positions and momenta of the gas particles respectively; each of the particles interacts only with the last particle $`x_N`$ of the Fermi–Pasta–Ulam system via the same potential as for the heath bath.

In our simulations we chose an equal number of particles for the three systems, while in principle the number of bath particles should be larger than that of the Fermi–Pasta–Ulam system, and this in turn should be larger than that of the thermometer. Our choice is dictated only by the computational power available: we cannot take the total number of particles too large, but at the same time the number of particles in each system cannot be too small if a good statistics has to be insured. Taking the same number $`N=100`$ of particles for each of the three subsystems, seemed to us a good compromise.

We took as units of mass, length and energy the values $`m`$, $`l_0`$ and $`V_0`$, which were thus put equal to one in our computations. The values of the parameters $`\mathrm{\Omega }`$ and $`\beta `$ were set equal to $`\mathrm{\Omega }=400`$ and $`\beta =3742`$ respectively. Such strange values come from the following consideration: the intermolecular interaction in a crystal is well represented by a Lennard–Jones potential, whose relevant parameters (the range and the strength) are of order one with our choice of units. On the other hand the Fermi–Pasta–Ulam potential should just be a Taylor expansion of the Lennard–Jones potential around the equilibrium position. Performing such a Taylor expansion and putting the parameters equal to one, the indicated values for $`\mathrm{\Omega }`$ and $`\beta `$ are found.

Finally in our numerical simulations we took $`L=25`$; this in order to ensure a sufficient total number of collisions (of the order $`10^6`$ in our actual integrations), while at the same time letting the gas particles be free for a large part of their paths.

## 4 Numerical results.

The integration step was taken equal to a twentieth of the shortest period $`\tau _\mathrm{f}=\pi /\mathrm{\Omega }`$ of the Fermi–Pasta–Ulam chain, and the numerical solutions were computed up to times of order $`210^7\tau _\mathrm{f}`$.

The numerical experiments were performed in the following way. For the bath we chose a temperature $`T_1`$ and took random initial conditions extracted from a Maxwellian at the chosen temperature $`T_1`$ (we also checked that the value of the mean kinetic energy should not deviate too much from the expected one, in order to avoid too large fluctuations); for the Fermi–Pasta–Ulam system we chose initial data at equipartition with a temperature $`T_1/10`$ and random phases; and finally for the thermometer we chose initial data in the same way as for the heath bath, but at a temperature $`T_1/10`$. We let the system evolve for a time $`10^4\tau _\mathrm{f}`$, and then began to compute the time averages of the kinetic energies of the gases and of the harmonic energy of the Fermi–Pasta–Ulam system. The results of the computations for two representative cases are shown in Figures 1 and 2, where we report, versus time, the temperatures (i.e. twice the kinetic energies per particle) of the gases and the harmonic energy per particle of the Fermi–Pasta–Ulam system (actually, time averages of such quantities are reported). Let us recall that such three quantities should be equal according to the equipartition principle, i.e. for sufficiently long times.

In Figure 1, we started with a temperature $`T_1=1`$. One sees that, after a time of order $`10^5\tau _\mathrm{f}`$, the temperatures of two gases and the specific harmonic energy of the Fermi–Pasta–Ulam system have become essentially equal, although still presenting significant fluctuations. This case should correspond to a situation of thermal equilibrium.

Figure 2 refers instead to the choice of $`T_1=0.4`$. One sees that the heat bath and the thermometer still reach the same temperature (although after a time almost one order of magnitude larger than before), while the specific harmonic energy per particle of the Fermi–Pasta–Ulam system remains well below the common temperature of the two gases, up to the observation time. Actually the curve is so flat that the Fermi–Pasta–Ulam system can be expected to possibly reach the equilibrium only on a totally different time scale. The global system seems indeed to be in a situation of metaequilibrium.

It appears, however, that the zeroth law can still be valid, and that the temperature “measured” by the thermometer is a good empirical one, because the temperatures of the two gases have become equal. To understand to which quantity of the Fermi–Pasta–Ulam system does this measured temperature correspond, in Figures 3 and 4 we report the spectra (time–averaged energies of the modes versus mode number) of the Fermi–Pasta–Ulam system at the end of the two runs. Figure 3 refers to the case of complete thermalization, and correspondingly a complete equipartition among the modes is found, as expected. More interesting is Figure 4: here equipartition obtains only among modes of sufficiently low frequency, say below $`k_{\mathrm{cr}}=10`$, while the energy starts decreasing for larger values of $`k`$, going down, say for $`k>25`$, to the initial equipartition value $`0.04`$. It does not appear as a surprise to observe that the mean energy of the low frequency modes essentially agrees with the common temperature of the two gases. This seems to indicate that, for our metastable state, the “good” definition of temperature of the Fermi–Pasta–Ulam system is the mean energy of the (sufficiently) low frequency modes.

## 5 Conclusions.

In conclusion, we hope to have shown, through our numerical study of a Fermi–Pasta–Ulam system in contact with two gases, that there are cases of metastable equilibrium for which a notion of temperature can be defined. However, at variance with the familiar case of equilibrium, such a temperature does not coincide with the “canonical” one, namely twice the mean kinetic energy per particle.

As a further comment, we would like to add that the metaequilibrium states met in Fermi–Pasta–Ulam systems present characteristics which are somehow opposite to those of glasses. Indeed, in the latter case the lack of thermalization is ascribed to the low frequency modes and correspondingly the thermometer measures the mean energy of the high frequency modes, which are the one being in mutual equipartition.

Finally, we would like to mention that the possibility of having a thermodynamics for situations of metaequilibrium, typically involving the presence of adiabatic invariants, was amply discussed in the second part of a very interesting paper of Poincaré , which appears to have been almost completely forgotten. In fact we became aware of such work only after completing the present work, through a conference of V. Kozlov . In fact, V. Kozlov was addressing only the problem dealt with in the first part of the paper of Poincaré, namely how it occurs that the fast variables of an integrable hamiltonian system approach equilibrium, notwithstanding the reversibility and the return property of the system. In the second part of his paper, Poincaré was instead considering a situation in which one has at first a quick relaxation to a “provisional equilibrium” while a “definitive equilibrium” would be attained after a much larger time, i.e. one is concerned, in his very words, with “very long times of first order” and “very long times of second order”, which is a situation analogous to the one discussed in the present paper.


text:
# Brain neurons as quantum computers: in vivo support of background physics

## 1 Introduction

Many authors have argued that consciousness can be understood as a macroquantum effect. In particular, Penrose  proposed that this takes place in microtubules, the ubiquitous hollow cylinders that among other things help cells maintain their shapes. It has been argued that microtubules can process information like a cellular automaton , and Penrose suggests that they operate as a quantum computer (see for a review ). On the other hand, a number of other authors have conjectured that environment-induced decoherence will rapidly destroy the quantum macrosuperpositions in the brain (see for instance , and references therein). Indeed, the decoherence (as well as dissipation and heating) is a serious obstacle to all applications exploiting quantum coherence. Recently, considerable effort has been devoted to designing strategies able to counteract the undesired effects of the coupling with an external environment. Notable examples of these strategies in the field of quantum information are quantum error correction codes, error avoiding codes , and ”parity kicks” method . It is only beginning of such activity. It is quite possible that their exist still unknown strategies (and processes ,) which can eliminate in principle any undesired effect of the environment. Nature can utilize these strategies (processes) in neurons. Therefore, the simple calculations of the decoherence characteristic times (without taking into account such corrective possibilities) cannot be considered as a final verdict.

Methods of statistical physics are actively used for investigation of probabilistic properties of different neurons -. In the present paper we will study the data obtained in vivo for so-called Ventral Tegmental Area (VTA).

The Ventral Tegmental Area (VTA) is a midbrain nucleus consisting of the dopaminergic cells. VTA is known as a part of the limbic brain that corresponds to such high brain functions as cognition, learning, rewarding and emotional behavior. VTA plays a key role in the generation of pleasure and in development of psychological drug addiction. This area is also involved in control of the gonadal hormones .

The VTA dopaminergic cells fire in irregular manner with a mean rare rate between 0.5 to 10 Hz. Firing patterns of the healthy VTA dopaminergic cells include single spikes and short bursts containing 3-7 spikes. It is believed that the bursting manner of firing is probably pulled on by rewarding stimuli accepted from the glutamate neurons originating from prefrontal cortex and hippocampus. Bursting activity of the normal VTA dopaminergic cells results in dopamine release in limbic brain. Dopamine release from VTA dopaminergic cells accompanies rewarding behavior. VTA dopaminergic cells are autoregulated via dopamine autoreceptors expressed on their somas. Additionally, activity of the VTA dopaminergic cells is regulated by serotonin, noradrenalin and acetylcholine.

Specific nonlinear properties of the VTA dopaminergic neurons are actively studied in recent years (see, for instance , and references therein). Relationship between these properties and rewarding function of VTA is one of the topical problems. Therefore, comparative analysis of the nonlinear properties of the VTA signals generated by normal brains and ones generated by brains with genetically suppressed rewarding function of VTA seems to be of significant interest. For this purpose we use genetically defined rat model of depression (Flinders Sensitive Rat Line - FSL). Moreover, the genetically defined depression of the FSL rats (which suppresses the bursting activity of the dopaminergic neurons) and anaesthesia used in our experiment (see below) allows us extract a background signal generated by the dopaminergic neurons. Using universality of the background signal we will combine the signals taken from 5 single dopaminergic neurons (3 FSL rats) in order to obtain statistically representative data set (about 10000 spikes).

Usually, the time intervals, $`T`$, between neighboring spikes are used in order to study statistical properties of neuron signals. In the present paper we will use inverse interspike intervals (interspike frequencies) $`\omega =T^1`$. While for probability density analysis there is no advantage in such choice, for a moment (multifractal) analysis the interspike frequencies turned out to be more informative (see below). In particular, we will show that the FSL VTA neuron signals exhibit multifractal properties for the interspike frequencies on the scales where the healthy VTA dopaminergic neurons exhibit bursting activity. The obtained generalized dimensions spectrum corresponds to the multifractal Bernoulli distribution and for the high moments this spectrum coincides with generalized dimensions spectrum calculated for a spectral measure of a quantum system (so-called kicked Harper model, actively used as a model of quantum chaos).

## 2 Experimental methods and materials

Male Sprague-Dawley tars were used in all experiments. Animals were anaesthetized with chloral hydrate (400 mg/kg., i.p.) and mounted in stereotaxic apparatus. The hole was drilled 4.2 mm. anterior from the interaural line and 1.0 mm lateral from the medial line. Extracellular recordings were processed by an electrode from VTA (8.0-8.6 mm dorsal from the lambda). The constant level of the anaesthesia was checked by EKG and chloral hydrate was added as necessary. Single unit recording was carried out by amplitude discrimination. Each recording from the single cell included at least 2000 spike events. After each experiment, the recording site was marked by a lesion caused by 15 mA DC for 10 sec. Brains were removed and stained with formalin before histological examination. Frozen sections were cut at 50 mm intervals. Microscopic examination of the sections was carried out aiming to verify that the electrode tip was placed in VTA.

## 3 Multifractal analysis

The interspike frequencies: $`\omega (n)=T(n)^1`$, have obvious meaning in the case of simple periodic signals with constant period $`T`$. For a varying interspike interval $`T(n)`$ (where $`n`$ is number of a spike in the spike series) meaning of $`\omega (n)`$ is not so clear, especially for a random-like $`T(n)`$. Therefore, it is useful to introduce a measure based on the interspike frequencies through a moving average

$$\omega _r=\frac{1}{r}\underset{i=n}{\overset{i=n+r}{}}\omega (i)$$
$`(1)`$
This measure has a simple meaning of an average frequency for $`r`$ spikes. Then, one can try to analyze this measure on existence of scaling properties using its moments:

$$\omega _r^pr^{\mu _p}$$
$`(2)`$
Usually, normalized (dimensionless) moments $`\omega _r^p/\omega _r^p`$ are used for this purpose. In our case, however, we have $`\mu _1=0`$ directly from the definition. Therefore, the scaling exponents of the dimensionless moments (if existst) coincide with those of the moments (2). If the scaling exists one can use the exponents $`\mu _p`$ in order to calculate the generalized dimensions $`D_p`$ :

$$D_p=1\frac{\mu _p}{p1}$$
$`(3)`$
Figure 1 shows composite signal, $`\omega (n)`$, obtained from 5 singular dopaminergic (VTA) neurons of 3 FSL rats. Because the genetically defined depression suppresses the bursting activity of the dopaminergic neurons and because of the anaesthesia used in our experiment we can expect a background signal from the FSL dopaminergic neurons. In particular, we can expect an universality of the individual signals (i.e. a statistic similarity even between the signals obtained from the dopaminergic neurons which belong to different FSL rats). The combined signal contains of about 10000 data points. A few data points with relatively short (less than 0.01 s) interspike intervals were excluded from the combined data set (that is about 1$`\%`$ of the entire number of the data points).

Figure 2 shows the moments (2) calculated for the data set shown in figure 1, for $`r=313`$. Log-log scales are chosen for comparison with the scaling equation (2). The straight lines (the best fit) are drawn to indicate scaling (2) in these scales.

Figure 3 shows the generalized dimensions $`D_p`$ (3) (calculated using the exponents $`\mu _p`$ extracted from figure 2). The scales: $`D_p`$ versus $`\mathrm{ln}p/(p1)`$ are chosen in figure 3 for comparison with the theoretical prediction for the multifractal Bernoulli distribution :

$$D_p=D_{\mathrm{}}+c\frac{\mathrm{ln}p}{(p1)}$$
$`(4)`$
The straight line is drawn in figure 3 to indicate the multifractal Bernoulli distribution asymptotic (4) with multifractal specific heat $`c0.41`$.

## 4 Discussion

The observed scaling interval: $`r=313`$, covers the scales where bursting activity of the healthy dopaminergic neurons takes place (see Introduction). Analogous analysis of the control (healthy) rats did not reveal the multifractal behavior in this interval of scales. This observation can be interpreted as multifractality of the background signal of the dopaminergic (VTA) neurons. Moreover, clear tendency to the multifractal Bernoulli distribution, seen in figure 3 for the high moments (the straight line), indicates a fundamental underlying physics (kinetics) . In order to proceed further in this direction we also show in figure 3 values of generalized dimensions calculated for a spectral measure of the kicked Harper model actively used as a model of quantum chaos . The kicked Harper model is obtained upon quantization of the following map:

$$p_{n+1}=p_n+K\mathrm{sin}(x_n),x_{n+1}=x_nL\mathrm{sin}(p_{n+1})$$

where $`K`$ and $`L`$ are some parameters. Canonical quantization thus leads to the one-period evolution operator

$$U_{L,K}=\mathrm{exp}[i\frac{L}{h}\mathrm{cos}(h\nu )]\mathrm{exp}[i\frac{K}{h}\mathrm{cos}(x)]$$

where operator $`\nu =i/x`$, and $`h/2\pi `$ has to be considered as an effective Planck constant, playing a role similar to an incommensurability parameter in a quasi-periodic system. Quasi-energy eigenvalues and eigenvectors are determined by

$$U_{L,K}\psi _\omega =\mathrm{exp}[2\pi i\omega ]\psi _\omega $$

and, if we denote by $`\psi _0`$ the state in which the system is prepared at $`t=0`$, the corresponding spectral measure will be indicated by $`d\eta _{\psi _0}(\omega )`$: its support is contained in the unit interval. The autocorrelation function is easily expressed in terms of the spectral measure as follows:

$$C(t)=\psi _0|U_{L,K}|\psi _0=_0^1𝑑\eta _{\psi _0}(\omega )\mathrm{exp}[2\pi i\omega t]$$
$`(5)`$
Information on the spectral measure is thus obtained by inversion of (5).

Multifractal analysis requires a sequence of approximations to the asymptotic spectral measure $`\rho _{i(l)}`$, where $`\rho _{i(l)}`$ is the probability attached to the $`i`$-th ball (of size $`l`$) covering the support of the spectral measure. The generalized dimensions is given by

$$D_p=\underset{k\mathrm{}}{lim}\frac{1}{(p1)}\frac{\mathrm{log}\chi _k(p)}{\mathrm{log}l_k}$$

where

$$\chi _k(p)=\underset{j=1}{\overset{N_k}{}}\rho _j(l_k)^p,l_k=1/(2N_k+1)$$

Phase diagram of the kicked Harper model is roughly divided into three regions , characterized, respectively , by a pure point spectrum (region A), a purely continuous spectrum (region B), and a mixed spectrum (region C). Within region A scaling features were consistently not observed. Region B is, on the other side, characterized by good scaling properties, and a converging, non-trivial spectrum of generalized dimensions were obtained in . The generalized dimensions for the parameters $`K=2`$, $`L=5`$ (with $`N=6400`$) are shown as crosses in figure 3.

Reported above results can be considered as an indication of existence of a background dopaminergic neuron signal with universal multifractal properties for the scales where healthy dopaminergic neurons exhibit bursting behavior. This signal has multifractal Bernoulli distribution as high moments asymptotic that may be related to an underlying physics (kinetics) of mesoscopic systems . In particular, the multifractal properties of a spectral measure of the kicked Harper model are in good quantitative correspondence with the observed multifractal properties of the dopaminergic neuron signals. The kicked Harper model is a prominent example of quasi-periodically driven quantum systems with a chaotic classical analogue, that also can be of significant interest for the discussion mentioned in Introduction.

The authors are grateful to K.R. Sreenivasan for cooperation and to the Machanaim Center (Jerusalem) for support.


text:
# Experimental investigation of coherent Smith-Purcell radiation from a ”flat” grating

## 1 Introduction

The diffraction radiation (DR) is stipulated by means of dynamic perturbation of electronic shells of atoms and free electrons by an electromagnetic field of a charged particle moving near a conductive target. An effective transversal size of a relativistic electron field for its component with wavelength $`\lambda `$ is about $`\gamma \lambda `$ where $`\gamma `$ \- is the Lorentz factor of an electron. This value can reach macroscopic sizes in ultrarelativistic case, allowing the distance between an electron beam and the target to be increased. A transition radiation, in contrast, is generated during passing of an initial particle through material of a target. In case of the periodic target to be used, there maybe a constructive interference of radiation from different elements of a target and the radiation becomes resonant in this case. The special case of resonant diffraction radiation, when the electron moves in parallel to a periodic target plane is referred as Smith-Purcell radiation (SPR) . There are two types of targets for this geometry of experiment: a volume grating obtained by a periodic deformation of a continuous surface, and a flat periodic target obtained by alternation of metal strips and vacuum (dielectric) gaps. Authors of the work  suggested to use CSPR for bunch length measurements. Changing the observation angle (in other words, the CSPR wavelength) one may detect a sharp increase in CSPR intensity when the wavelength becomes more than a bunch length. But the usage of a traditional lamellar grating for this aim has some disadvantages: for a set of observation angles there may exist deep minima in the CSPR angular distribution . It leads to significant troubles in interpretation of experimental data. Another method of the beam diagnostics suggested in  is based on the rotation of the grating in the grating plane to change an effective grating period and needs the developed theory of SPR effect for a rotated grating that is absent now. As with of grating rotation the spectra of CSPR became wider and this technique should be investigated both theoretically and experimentally. The use of a flat grating for the beam diagnostic technique allows one to avoid the angular distribution unevens and to increase the radiation yield .

## 2 Theoretical model

The basis of theoretical model for calculation of radiation characteristics from a flat target is the exact solution of Maxwell equations for diffraction radiation caused by relativistic electron moving close to a thicknessless ideal conductive half-plane .

As shown in  the spectral - angular distribution of SPR from a flat grating consisting of  N strips with width a and period d, separated by vacuum gaps with width  (d-a) can be written as:

$`{\displaystyle \frac{d^2W_{SP}}{d\lambda d\mathrm{\Omega }}}=4{\displaystyle \frac{\alpha }{\pi }}{\displaystyle \frac{\mathrm{}c}{\lambda ^2}}{\displaystyle \frac{exp\left[\frac{4\pi h}{\gamma \lambda }\right]}{1/\beta cos\theta }}sin\left[{\displaystyle \frac{\pi a}{\lambda }}\left(cos\theta {\displaystyle \frac{1}{\beta }}\right)\right]F_N`$ (1)
Here $`\alpha =\frac{1}{137}`$ -is a fine structure constant, $`\mathrm{}c=0.2eV\mu m`$\- conversion constant, $`h`$\- distance between the particle trajectory and the target surface, $`\theta `$ \- outgoing photon angle, $`\beta =\sqrt{1\gamma ^2}`$\- velocity of particle in units of light speed,$`\gamma `$ is Lorentz-factor, $`\lambda `$ \- is SPR wavelength, $`F_N=\frac{sin^2\left(N\varphi _0/2\right)}{sin^2\left(\varphi _0/2\right)}`$ \- coefficient, taking into account the interference of radiation from $`N`$ periodically spaced strips,

$$\varphi _0=\frac{2\pi d}{\lambda }\left(cos\theta \frac{1}{\beta }\right)$$

Fig. 1 shows the calculated dependence of radiation wavelength on the observation angle which follows from the equation for $`\varphi _0=2\pi k`$ (Smith-Purcell dispersion law). It is evidently that, with increasing the observation angle the spectral lines are shifted to large wavelength region.

With the radiation wavelength exceeding the longitudinal size of an electron bunch, all the electrons from a bunch emit coherently. In this case radiation yield, normalized to one electron, increases proportionally to the number of electrons in the bunch. The yield of coherent radiation for a bunch with number of electrons $`N_e`$ in the aperture of the detector $`\mathrm{\Delta }\mathrm{\Omega }`$ registered by the detector with bandwidth $`\mathrm{\Delta }\lambda `$ can be expressed as follows:

$`\mathrm{\Delta }W_{exp}(h,\theta ,L_{bunch})=N_e{\displaystyle \underset{\mathrm{\Delta }\lambda }{}}\left|f\left({\displaystyle \frac{\lambda }{L_{bunch}}}\right)\right|^2{\displaystyle \underset{\mathrm{\Delta }\mathrm{\Omega }}{}}{\displaystyle \frac{d^2W_{SP}}{d\lambda d\mathrm{\Omega }}}𝑑\mathrm{\Omega }\epsilon \left(\lambda \right)d\lambda `$ (2)
Here $`f\left(\frac{\lambda }{L_{bunch}}\right)`$\- is the longitudinal form-factor of a bunch, $`\epsilon \left(\lambda \right)`$ \- spectral efficiency of the detector, $`L_{bunch}`$ \- length of the bunch. In Fig. 2 is given the relation between the squared form-factor modulus and radiation wavelengths for gauss longitudinal distribution of electrons in a bunch. In figure we consider a bunchlength $`L_{bunch}`$ equal to $`3\sigma `$ in longitudinal distribution of electrons. It is clear from the figure that, for example, for $`50\%`$ level of a maximum intensity corresponds wavelength $`\lambda =0.7L_{bunch}`$. That is, the measurement of angular distribution of CSPR (in other words, measurement of spectral dependence of a CSPR yield) allows a length of an electron bunch to be determined.

## 3 Experimental setup

The experiment was carried out in NPI MSU on a linac with the following parameters of electron beam (see Table 1).

The experimental setup is shown in Fig. 3. The beam with energy 5MeV is extracted through a 20 micrometer titanium foil and moves at a distance $`h=2`$ mm from a periodic target.

Table1

| Electron energy | $`E_e=5`$ MeV |
| --- | --- |
| Bunch population | $`10^9`$ |
| Bunch length | $`L_{bunch}12`$ mm |
| Bunches in macropulse | 22 |
| Macropulse duration | $`\tau 7`$ ns |
| Distance between the target center and the detector | $`R=1`$ m |

The target, size $`50\times 50\times 1.5`$ mm consisting of 8 copper strips with 3 mm width and $`30`$ micrometer thickness, located 3 mm apart on an insulating substrate, was placed in parallel to beam in the goniometer, by which it was possible to tune the parallelism between the beam trajectory and the target and change the impact parameter $`h`$.

The detector may be rotated on the table with radius $`R=1`$ m that allows one to change angle of photon detection. The center of rotation coincides with a target center. In the experiment we registered CSPR together with transition radiation (TR) background from the exit window foil.

The detector was made on the basis of a room-temperature resonance funnel-shaped antenna with the following characteristics: resonant wavelength $`\lambda _{res}=7.2`$ mm, a band width $`\mathrm{\Delta }\lambda 2`$ mm, sensitivity: $`0.18`$ V/W, response time $`\tau _D=10`$ ns. The polar angle acceptance $`\mathrm{\Delta }\theta 12mrad0.06\gamma ^1`$. Assuming that the dielectric constant of the target insulator substrate is the same as for vacuum, we have calculated the radiation power registered by the detector for our condition and $`L_{bunch}=1.5`$mm (see formula(2)) and obtained value $`\mathrm{\Delta }W_{exp}/\tau 0.9W`$.

## 4 Experimental results

Figure 4 presents the calculated dependence of CSPR intensity on the angle of detection in the real aperture of the detector with regard to the spectral efficiency. The experimental relations between CSPR intensity and angle of observation, the impact parameter and a current of an electron beam are shown in Figures 5,6,7.

Comparing Figures 4 and 5 it is possible to notice the reasonable agreement between the theoretical model and experiment. The monochromaticity of the main maximum is defined by the number of the lattice periods and band pass of the detector $`\mathrm{\Delta }\lambda /\lambda >1/N12\%`$.

The narrow peaks in Fig. 5 for angles $`\theta =56^0`$ and $`73^0`$ can be explained by resonances of a narrow-band detector. To get dependence of CSPR intensity on the impact parameter, the detector was installed at angle $`\theta =95^0`$ that corresponded to the maximal value of CSPR.

In Fig. 6 the value $`h=0`$ corresponds to a beam passing exactly through a target and the dip at this value $`h`$ is explained by scattering the electrons on target material.

The approximation of dependence of CSPR intensity on the impact parameter by exponential function $`f=F+aexp\left[bh\right]`$, where $`b=\frac{4\pi }{\gamma <\lambda >}`$, gives the value of mean detected wavelength of CSPR $`<\lambda >=7.9\pm 1.3mm`$, that agrees well with resonant wavelength of the detector.

For confirming a coherent character of radiation we measured a dependence of its intensity on pulse current of a macrobunch, $`I_{acc}`$, proportional to the number of particles in the bunch. It canbe seen from Fig. 7 that for $`I_{acc}<I_012`$ mA the dependence of CSPR yield is $`YI_{acc}^c`$, where $`c=2.02\pm 0.07`$. With experimental error accuracy, the dependence of CSPR intensity on the current of electron beam is square-law, that confirms the coherent character of radiation. For values $`I_{acc}`$, exceeding $`12`$mA detector starts to work in a nonlinear regime. It is necessary to notice that the experiment was performed for $`I_{acc}<I_0`$ condition.

## 5 Conclusion

In conclusion we may note that using $`5`$ MeV beam of a linac NPI MSU we investigated characteristics of CSPR and showed, that the CSPR mechanism provided the high monochromaticity of radiation even for a target with $`8`$ periods only and the aperture of the detector $`\gamma ^1`$, whereas the continuous component in a spectrum connected with background processes (TR from exit window, for instance) is small enough. Also it is necessary to pay attention to the fact, that for $`\theta >>1/\gamma `$ the CSPR yield practically does not depend on $`\gamma `$ (for $`\gamma =10`$ and for $`\gamma =1000`$ difference is a few percent). The macropulse power of CSPR is proportional to $`nN_e^2`$, where $`n`$ \- the number of bunches in macropulse, (for a linac NPI MSU estimation from (2) gives $`PnN_e^2\mathrm{\Delta }W_{exp}/\tau 0.9W`$). For such accelerators as KEK-ATF and SLAC FFTB with number of electrons in a bunch $`N_e10^{10}`$, this power level is reached for a single bunch. Thus, the received power level demonstrates a possibility of measurement of a single bunch length using a flat grating and a set of detectors with responsivity: $`0.1V/W`$, operating at room temperature and with different resonance wavelengths located at different angles $`\theta `$. A resonance condition and requirement of coherency of a bunch radiation define the wavelength range accessible to the bunch length measurement by this method: $`\lambda dL_{bunch}`$. For example, for $`L_{bunch}=0.1`$ mm,$`d0.1`$ mm,$`\lambda 0.1`$ mm one can use the detector for a far infrared range.

## Figure captions

Figure 1: Calculated dependence of radiation wavelength on observation angle (d=6 mm, $`\gamma =10`$, $`k=1,3\mathrm{}`$\- diffraction orders).

Figure 2: Dependence of square of a form-factor on radiation wavelength for gauss longitudinal distribution of electrons in a bunch.

Figure 3: The experiment scheme

Figure 4: The calculated relation between CSPR intensity and angle of detection $`\theta `$ in the real aperture of the detector with account of spectral effectiveness.

Figure 5: Experimental dependence of CSPR intensity from observation angle. Dotted curve is a continuous background “substrate”.

Figure 6: CSPR intensity as a function of impact parameter

Figure 7: Experimental relation between CSPR intensity and an electron beam current.