text:
# High performance Beowulf computer for lattice QCD

## 1 INTRODUCTION

The ZhongShan University Computational Physics group’s interests  cover such topics as lattice QCD quantum instanton and quantum chaos. Most of these topics can be studied through Monte Carlo simulation, but can be quite costly in terms of computing power. In order to do large scale numerical investigations of these topics, we required a corresponding development of our local computing resources.

The demarcation between super computers and personal computers has been further blurred in recent years by the high speed and low price of modern CPUs and networking technology and the availability of low cost or free software. By combining these three elements - all readily available to the consumer - one can assemble a true super computer that is within the budget of small research labs and businesses.

We document the construction and performance of a Beowulf cluster of PCs, configured to be capable of parallel processing.

## 2 SYSTEM

Our cluster consists of ten PC type computers, each with two Pentium III-500 CPUs inside. The logic behind dual CPU machines is that one can double the number of processors without the expense of additional, cases, power supplies, motherboards, network cards, et cetera. Also, the inter-node communication speed is faster for each pair of processors in the same box as compared to communication between separate computers. Each computer has an 8GB EIDE hard drive, 128 MB of memory, a 100Mbit/s Ethernet card, a simple graphics card a floppy drive and a CDROM. In practice the CDROM, the floppy drive, and even the graphics card could be considered extraneous, as all interactions with the nodes could be done through the network. One computer has a larger hard disk (20 GB), and a SCSI card for interaction with a tape drive. For the entire cluster we have only one console consisting of a keyboard, mouse and monitor.

A fast Ethernet switch handles the inter-node communication. The switch has 24 ports so there is ample room for future expansion of the cluster to up to a total of 48 processors. Of course it is possible to link multiple switches or use nodes with more that two processors, so the possibilities for a larger cluster are nearly limitless.

We have installed a Red Hat Linux 6.1 distribution. It automatically supports dual CPU computers. It is also able to support a Network File System (NFS), allowing all of the nodes in the cluster to share hard disks, and a Network Information System (NIS), which standardizes the usernames and passwords across the cluster.

We can use the the cluster for parallel processing by using the message passing interface (MPI), a library of communication functions and programs that allow for communication between processes on different CPUs. The programmer must design the parallel algorithm so that it appropriately divides the task among the individual processors. He or she must then include message passing functions in the code which allow information to be sent and received by the various processors.

## 3 PERFORMANCE

As we primarily developed the cluster for numerical simulations of lattice QCD, we have also performed a benchmark which specifically tests the performance in a parallel lattice QCD code. The algorithm can conveniently divide the lattice and assign the sections to different processors.

Hioki and Nakamura  provide comparison performance data on SX-4 (NEC), SR2201 (Hitachi), Cenju-3 (NEC) and Paragon (Intel) machines. Specifically, we compare the computing time per link update in microseconds per link and the inter-node communication speed in MB/sec. The link update is a fundamental computational task within the QCD simulation and is therefore a useful standard. The test was a simulation of improved pure gauge action ($`1\times 1`$ plaquette and $`1\times 2`$ rectangle terms) on a $`16^4`$ lattice. In each case the simulation was run on 16 processors. The results are summarized in Table 1.

A widely used QCD program is the MILC code . It has timing routines provided so that one can use the parallelized conjugate gradient routine in the simulation as a benchmark. Furthermore, as this code is very versatile and is designed to be run on a wide variety of computers and architectures. This enables quantitative comparison of our cluster to both other clusters and commercial supercomputers. In the MILC benchmark test we ran to a convergence tolerance of $`10^5`$ per site. For consistency with benchmarks performed by others, we simulated Kogut-Susskind fermions.

We have run the benchmark test for different size lattices and different numbers of processors. It is useful to look at how performance is affected by the number of CPUs, when the amount of data per CPU is held fixed, that is each CPU is responsible for a section of the lattice that has $`L^4`$ sites. For one CPU, the size of the total lattice is $`L^4`$. For two CPUs it is $`L^3\times 2L`$. For four CPUs the total lattice is $`L^2\times (2L)^2`$; for eight CPUs, $`L\times (2L)^3`$, and for 16 CPUs the total size of the lattice is $`(2L)^4`$.

Note that the falloff in performance with increased number of CPUs is dramatic. This is because inter-processor message passing is the slowest portion of this or any MPI program and from two to sixteen CPUs, the amount of communication per processor increases by a factor of four. Table 2 shows that for a lattice divided into $`2^j`$ hypercubes, each of size $`L^4`$, there will be $`j`$ directions in which the CPUs must pass data to their neighbors. The amount of communication each processor must perform is proportional to the amount of interface per processor. As this increases, per node performance decreases until $`j=4`$ and every lattice dimension has been divided (for a $`d=4`$ simulation). and the per-processor performance should remain constant as more processors are added. The shape of this decay is qualitatively consistent with $`1/j`$ falloff.

Of course there are other ways to divide a four-dimensional lattice. The goal of a particular simulation will dictate the geometry of the lattice and the therefore the most efficient way to divide it up (generally minimizing communication). A four-CPU simulation using a $`4L\times L^3`$ lattice has the four hypercubic lattice sections lined up in a row (as opposed to in a $`2\times 2`$ square for a $`L^2\times (2L)^2`$ lattice) and has the same amount of communication per CPU as does the $`L^3\times 2L`$ two-CPU simulation. In a benchmark test the per-CPU performance was comparable to the performance in the two-CPU test.

For a single processor, there is a general decrease in performance as $`L`$ increases, as shown in Table 3. This is well explained in  as due to the larger matrix size using more space outside of the cache memory, causing slower access time to the data.

For multiple CPUs there is in performance improvement as $`L`$ is increased. The explanation for this is that the communication bandwidth is not constant with respect to message size. For very small message sizes, the bandwidth is very poor. It is only with messages of around 10kB or greater that the bandwidth reaches the full potential of the fast Ethernet hardware, nearly 100Mbit/sec. With a larger $`L`$, the size of the messages is also, improving the communication efficiency. The inter-node communication latency for our system is 102$`\mu s`$. As inter-node communication is the slowest part a parallel program this far outways the effect of cache misses.

To summarize, a parallel cluster of PC type computers is an economical way to build a powerful computing resource for academic purposes. On an MPI QCD benchmark simulation it compares favorably with other MPI platforms. The price/performance ratio is $`\$7/Mflop`$. It is drastically cheaper than commercial supercomputers for the same amount of processing speed. It is particularly suitable for developing research groups in countries where funding for pure research is more scarce. We have been doing large scale calculations of hadron and glueball spectrum.

X.Q.L. is supported by National Science Fund for Distinguished Young Scholars (19825117), Guangdong Provincial Natural Science Foundation (990212), Ministry of Education, Gosun Ltd., and the Foundation of Zhongshan Univ. Advanced Research Center. We also thank S. Hioki and S. Gottlieb for their benchmark codes.


format: text

text:
# Gamma-Ray Imaging by Silicon Detectors in Space: the AREM Method

## 1 Introduction

High energy gamma-ray astrophysics is foreseen to be one of the more challenging fields of study in the coming years. Previous gamma-ray missions such as SAS-2 , COS-B , and especially the recent EGRET experiment  on the Compton Gamma Ray Observatory, left us with a large amount of exciting results and open questions. Indeed, of the nearly 300 gamma-ray sources detected so far only a small fraction ($`30\%`$) has been identified . The discovery of gamma-ray blazars, pulsars, high-energy gamma-ray bursts (GRBs) and of a large amount of unidentified sources, many of which are strong high-energy transients, has provided clear evidence for the need of a next generation of gamma-ray experiments with increased field of view and improved angular resolution. In this context, the AGILE mission , planned to be operational during the year 2003, is integrated toward the GLAST mission , planned for the year 2006. In this paper, we consider the EGRET spark chamber experiment as the direct AGILE predecessor, and we will refer to it for comparisons.

We recall here that AGILE (Astro-rivelatore Gamma a Immagini LEggero)111 More information about the project can be found at the site: http://www.ifctr.mi.cnr.it/Agile.  is a Small Scientific Mission of ASI (Agenzia Spaziale Italiana) with a tracking system based on the state-of-the-art Silicon strip technology . The AGILE GRID (Gamma-Ray Imaging Detector) consists of a Silicon-Tungsten Tracker, a Cesium-Iodide Mini-Calorimeter and a segmented Anticoincidence system of plastic scintillators and is sensitive to 30 MeV - 50 GeV photons. Thanks to the fast readout electronics and to the segmented Anticoincidence system, AGILE will have, among other features, an unprecedently large field of view (FOV $`3`$ sr), larger than previous gamma-ray experiments by a factor $`5`$. Furthermore the Silicon Tracker has a very good intrinsic spatial resolution, comparable or even less than the microstrip pitch $`<121\mu `$m, by using the analog information on the released charge distribution between strips. For comparison, we recall that the EGRET spark chamber wire spacing was equal to $`820\mu `$m.

The AGILE goal is to obtain the best sensitivity ever reached for off-axis events (up to $`60^{}`$), and an on-axis sensitivity comparable to that of EGRET, despite the smaller dimensions and effective area. Therefore, the optimization of the angular resolution algorithms becomes a crucial point to fulfil the mission scientific objectives.

The $`\gamma `$-ray photon direction reconstruction is based on the physical process of pair production, and is obtained from the identification and the detailed analysis of the electron/positron tracks stemming from a common vertex. The direction reconstruction should take into account the effect of multiple Coulomb scattering and the distribution of the total energy of the incident photon between the $`e^+/e^{}`$ particles. Until now, this was done by analyzing separately the two tracks projections in the ZX and ZY views. We will show in section 2 that the “2-D projection method” is a good approximation only for nearly on-axis events, but it induces two kinds of systematic error in the $`\gamma `$-ray direction reconstruction for off-axis events because of the intrinsic ambiguity in the proper identification of 3-D tracks and of the problem of the correct reconstruction of the true 3-D photon direction. Finally, we point out the importance of the choice of a proper track weighting scheme.

In this paper we emphasize that, contrary to previous $`\gamma `$-ray experiments, the simplified 2-D projection method would not be a good approximation for AGILE because of its large field of view and very good intrinsic spatial resolution. As we discuss in the next section, the first source of error follows from the ambiguity in the association of coordinate pairs $`(X,Y)`$ corresponding to the hits in the Tracker in the two orthogonal ZX and ZY views. This ambiguity can be resolved in all those cases in which one of the two tracks becomes distinguishable from the other in both views of the Tracker, as for example when it stops or exits the Tracker, or when it suffers a more significant multiple scattering. Regarding the second point, we notice that in general, and expecially for off-axis events, the 3-D reconstruction is not equivalent to the composition of reconstructed directions in each projected view. So we are forced to consider an intrinsically three-dimensional strategy in order to reconstruct the true incident photon direction. The third point keeps into account the fact that, in general, the photon energy is not evenly divided between the two pair particles. It turns out that the incident photon direction is closer to that of the most energetic particle, and an “energy-weighted” direction reconstruction is necessary.

The AGILE REconstruction Method (AREM) presented in this paper, provides a general baseline, to be optimized for each particular $`\gamma `$-ray instrument, to keep into account these three points of the analysis.

## 2 The AREM method

### 2.1 The conversion plane problem

Let us first illustrate and discuss the conversion plane problem in the simplified case of an even energy sharing between the electrons and positrons i.e., when the photon direction coincides with the bisector. We restrict the discussion to the idealized case of absence of $`\delta `$-rays or other secondaries in the Tracker. In general, each Si-Tracker plane of a $`\gamma `$-ray detector is composed of a converter layer and by two orthogonal strips layers. Charged particles entering the Tracker will produce an electric signal in the front-end electronics which would correspond to the hit readout strips in the two active X and Y layers. As long as a single charged particle crosses the Tracker, this signal will identify a unique point in 3-D space for each hit plane. However, when two separated particles hit simultaneously the active layers (as in the case with $`e^+/e^{}`$ pairs) the signal could correspond to two possible couples of points in 3-D space, as shown in Fig. 1.

The two possible couples of points are: pair LL/RR = $`[(X_L,Y_L),(X_R,Y_R)]`$, open symbols in the figure, or pair LR/RL = $`[(X_L,Y_R),(X_R,Y_L)]`$, filled symbols. This gives rise to an intrinsic ambiguity in the conversion plane identification. In terms of projected views this “conversion plane problem” can be phrased as: “Does the track to the left in the ZX view correspond to the left track in the ZY view - pair LL/RR? Or does it correspond to the track to the right - pair LR/RL?”.

Note that for on-axis events the “pyramid” in Fig. 1 becomes isosceles and the conversion plane problem apparently disappears, since both alternatives - pair LL/RR or LR/RL - would lead to the same bisector. However, we recall that we are considering here the particular and unrealistic case of an even energy sharing. In general the photon energy is not evenly divided between the two particles (see Sect. 2.3), and the photon direction does not coincide with the bisector. Hence the ambiguity for the identification of the “energy-weighted” reconstructed direction applies also to on-axis events.

AREM algorithms solve the intrinsic ambiguity in the identification of the 3-D $`e^+/e^{}`$ tracks by identifying a primary and a secondary track in both ZX and ZY views. The primary track is defined as that one carrying most of the incident energy, and subject to less multiple scattering. In practice, we identify the primary track as the one for which (starting from the third plane) the quantity $`|X_{CLUSTER}X_{EXTRAPOLATED}|`$ is on average minimized. This means that the primary track should have the smallest mean distance between the selected clusters and the position obtained by extrapolation of the trajectory from the previous planes. The 3-D primary track is then obtained by associating the primary track in the ZX view with the primary track in the ZY view (analogously for the secondary one).

### 2.2 The projection problem

According to the simplified 2-D projection method, the first step is to identify the two projected tracks in each view. The following step is to take their (in general weighted) bisectors and compose them to obtain the reconstructed gamma direction. This procedure is not correct: the true 3-D bisector has ZY and ZY projections which do not correspond to the two bisecting lines in each view. This is illustrated in Fig. 2. The only case in which this equality holds true is for on-axis events. As shown in ref. , this systematic effect increases for increasing off-axis angles and large opening angles up to values of $`0.5^{}`$, hence the 2-D projection method is not acceptable for AGILE and similar detectors222We warmly thank the EGRET collaboration, and in particular D.L. Bertsch and D. Thompson, for allowing us to perform a test of our reconstruction algorithms on their calibration data. In the case of EGRET data this systematic error can be considered almost negligible, since it is hidden by the relatively low spark chamber intrinsic resolution..

The proper identification with AREM of the 3-D pair production plane and photon direction reconstruction can be obtained by several completely equivalent geometrical methods, such as: spherical triangle formulae, normalized 3-D vector components, or double rotation of reference system. In the following we use the latter method.

### 2.3 Track weighting scheme

In the pair production process the total energy of the incident photon is split between the $`e^+/e^{}`$ particles. We recall that the probability distribution of having an electron, or a positron, with energy $`E_{e^+(e^{})}`$ when the primary photon has energy $`E_\gamma `$, is almost flat in the “low” energy range $`E_\gamma =10÷40`$ MeV. However, for increasing incident energy it is more and more probable that one of the two secondary particles carries most of the total energy (see for example Fig. 2.19.2 of ref. ). The RMS emission angle $`\alpha _{e^+(e^{})}`$ between the original photon direction and the direction of the electron, or positron, roughly depends on the inverse of the photon energy : $`\alpha _{e^+(e^{})}1/E_\gamma `$. This implies that for most events in the AGILE energy range the main information on the original photon direction will be carried by the most energetic $`e^+(e^{})`$ particle. It follows that the choice of a proper track weighting scheme to keep into account the energy distribution between the $`e^+/e^{}`$ particles plays a fundamental rôle in the angular reconstruction algorithms.

Can we obtain a satisfactory track energy estimate relying mainly on Tracker data? Multiple scattering spreads the information on the original particle directions, but on the other side it can be exploited to obtain a track energy estimate. According to the Molière theory , the small planar deflection angle of multiple scattering $`\beta `$ has an approximately Gaussian distribution with high tails (for a detailed analysis see also ref. ). The Molière formula for the RMS scattering angle can be written as:

$`\beta ={\displaystyle \frac{13.6\mathrm{MeV}}{pv}}\sqrt{{\displaystyle \frac{X}{X_0}}}\left[1+0.038\mathrm{ln}({\displaystyle \frac{X}{X_0}})\right],`$ (1)
where $`X/X_0`$ is the thickness of the scattering medium in radiation lengths. In our case, the relation between particle momentum and energy is: $`pcE_{e^+(e^{})}`$, and the velocity is such that: $`v/c1`$, since we consider $`E_{e^+(e^{})}>>m_ec^2`$. There is then an (approximate) inverse relation between particle energy and RMS deflection angle:

$`E_{e^+(e^{})}{\displaystyle \frac{\sqrt{cos\theta }}{\beta }},`$ (2)
where here $`\theta `$ is the angle between the considered track and the vertical axis. This means that the most energetic particle suffers less multiple Coulomb scattering.

The AREM choice of weighting scheme is based on Molière multiple scattering theory. In principle, each 3-D track should be weighted by its energy to some power. According to eq. (2), for each projected track we define a weight as:

$`w(E_{e^+(e^{})})^r(\sqrt{cos(\theta )}/\beta )^r,`$ (3)
where $`r`$ is a parameter to be optimized for each particular $`\gamma `$-ray instrument. In practice, the quantity $`|(X_{CLUSTER}X_{EXTRAPOLATED})|`$ is used to estimate the RMS multiple scattering angles $`<\beta >`$, and hence the particle energies $`E_{e^+(e^{})}`$.

### 2.4 The AREM flow

Summing up, we show in Fig. 3 the general flow chart for the software corresponding to the AREM method described in the previous subsections.

## 3 Applications to the AGILE Tracker

### 3.1 Instrumental and geometrical set-up

The AGILE baseline configuration  is the following:

* Silicon Tracker: the AGILE Tracker consists of 14 detection planes $`42\times 42`$ cm2 total area with $`1.6`$ cm interplane distance. Each plane is composed by two layers of microstrip Silicon detectors with orthogonal strips (X and Y view). The first 12 planes, starting from the top, have also a Tungsten conversion layer each, with thickness equal to $`0.07X_0`$ radiation length. The last two planes have no Tungsten layers, since the readout trigger requires at least three Si-layers providing signals. The active area of the Si-detectors is made of $`4\times 4`$ tiles where $`4\times 768`$ microstrips of pitch of $`121\mu \mathrm{m}`$ are implanted. Only half of the strips in every tile are readout strips. However, it is possible to detect a signal even if the particle hits a “floating” strip, by taking into account the effect of the capacitive coupling between contiguous strips333The “floating strip” configuration has been chosen to achieve an excellent spatial resolution while keeping under control the number of readout channels and hence the detector power consumption..
* Mini Calorimeter: it consists of two orthogonal planes, each containing 16 CsI bars, for a total radiation length of $`1.5X_0`$.
* Anticoincidence System: it is made of lateral segmented Plastic Scintillator planes composed of 3 panels for each Tracker side and of a top plane of thickness of $`0.5`$ cm.

### 3.2 Monte Carlo simulations

The Monte Carlo simulations of AGILE performances are done with the GEANT 3.21 code . The GEANT package reproduces all possible interactions of radiation and particles through matter, such as pair production, multiple scattering, Compton scattering, etc.

In the case of the AGILE Tracker, CERN testbeams indicate that the experimental uncertainty in the identification of the position of a cluster of adjacent strips providing signal on a track can be reduced at the level of $`40\mu `$m . This result is achieved by using the analog readout which gives information on the charge distribution released in the Si-microstrips, and the floating strip configuration. However, since the experimental validation of the parametrization of the capacitive coupling effect at large off-axis angles is still in progress, in order to implement and test the AREM method we use here a simplified version of the AGILE simulation code, originally developed by the AGILE simulation software group . In this simplified model, which does not yet include capacitive coupling and floating strips, all the strips are readout strips and the digitization uncertainty is simply determined by the pitch. This detector model is very accurate and detailed enough to be fully suitable for our analysis.

The simulated events are characterized by the incident photon energy, $`E`$, and by a fixed incident direction ($`\theta `$ , $`\varphi `$) with respect to the detector sistem of reference. The algorithms which generate flux sources originating from a fixed direction are optimized to correctly simulate a plane wave front from an infinite distance .

### 3.3 Analysis of simulated events

From our Monte Carlo simulations, we obtain the complete detector response to the incident $`\gamma `$-ray flux, at several energies and incident angles. We recall that the standard data analysis based on Kalman filter algorithms , which make use of the information from all hit planes, is optimized for charged particle track identification, but not specifically for $`\gamma `$-ray direction reconstruction. Hence, as a first step to test our reconstruction method, we perform an analysis of the first n-Planes Resolution, i.e. using only information from the first hit planes, for which the information on the original photon direction is less disturbed by multiple scattering. In this preliminary study we present in particular our results regarding AREM 2PR and 3PR (2PR = two planes resolution, 3PR = three planes resolution).

In view of our introductory discussion, multiple scattering has a two-fold effect: it disperses the information on the original direction of the particle going through subsequent converter layers; on the other hand, it can be exploited to obtain a particle energy estimate. Integrating AREM with the high efficiency algorithms based on the Kalman filters will be the next step to improve the energy determination and weighting scheme, and to further optimize the AGILE angular resolution .

### 3.4 Event selection

According to AGILE Level-1 trigger conditions, there must be at least three consecutive Tracker planes which give a signal in both views444 The actual on-board first level trigger is given by the coincidence of 3 out of 4 consecutive Tracker planes, to take into account the possibility of one-plane failure. Simulations show that the final difference on the total number of rejected events at subsequent trigger level is of order of $`1\%`$ .. With respect to trigger Level-1, we perform a further event selection by asking that accepted events show only one cluster of strips on the first hit Tracker plane in both views. This amounts to require that it be possible to identify a “clean” conversion vertex in both Tracker views. We also require no more than two clusters on the second hit Tracker plane, always in both ZX and ZY views. In such a way we reject all events with spurious hits on the second plane, i.e., we require a “clean” second plane. We do not impose any further conditions on the third and subsequent planes.

### 3.5 Two plane resolution

We have analized first the so-called two plane resolution (2PR) of the Tracker, which takes into account only the information coming from the first two hit planes. The 2PR gives no further information on the deviation from the initial track direction obtained from the firts two planes and we are forced to neglect multiple scattering and energy distribution effects. We can only take a random choice regarding the conversion plane problem described in Sect. 2.1, and we get a coordinate ambiguity error for about $`50\%`$ of the photons. Regarding the projection problem described in Sect. 2.2, we can still utilize the 3-D method to guarantee the proper reconstruction of the bisector.

### 3.6 Three plane resolution

We then include in our analysis the information coming from the third hit plane. For the three plane resolution (3PR), the event selection criteria remain the same as in the 2PR case, and are the ones described above. This implies that the efficiency of the two reconstruction algorithms are equal. In this case we are in the condition to follow the complete AREM flow and to correctly address the intrinsic track ambiguity and the projection problem. Having identified a primary and a secondary track, we can assign a weight to each track as described in Sect. 2. In the following, we chose as a particular weighting scheme: $`wE_{e^+(e^{})}`$, i.e. we set $`r=1`$ in eq. (3), and the track energy is estimated from the inverse of square root of the variance of the track linear fit.

### 3.7 Results

In Fig. 4 we show the 3-D integral Point Spread Function (PSF) for events at $`\theta =10^{}`$, for the energies $`E=1\mathrm{GeV}`$, $`E=400\mathrm{MeV}`$, $`E=200\mathrm{MeV}`$ and $`E=100\mathrm{MeV}`$. In that figure we compare the PSF profiles obtained by 2PR and 3PR. The 3PR $`68\%`$ containment radius (c.r.) is systematically better than the 2PR one for $`E\stackrel{>}{}100\mathrm{MeV}`$. Since the efficiency of the two reconstruction algorithms is the same, the improvement of the angular resolution obtained with AREM 3PR is in part due to the solution of the coordinate track ambiguity and of the adopted weighting scheme. For the event set considered in Fig. 4 this effect ranges from $`10\%`$ up to $`30\%`$ at higher energy. This is also shown in Fig. 5, were we compare the 3-D PSF at $`68\%`$ c.r. as a function of energy in the case of EGRET on-axis and of AGILE 2PR and 3PR at incidence angle $`\theta =10^{}`$.

The comparison between the results obtained with the 2-D and 3-D reconstruction methods, both applied to the 3PR algorithm, is shown if Fig. 6 which summarizes the effect of AREM on the AGILE resolution at incidence angle $`\theta =30^{}`$.

In Fig. 7 we report the final results obtained with AREM 3PR for AGILE angular resolution near on-axis and at $`\theta =30^{}`$ off-axis. The efficiency of the AREM 3PR reconstruction (same as 2PR) is shown in Fig. 8 and it is compatible with the one obtained by EGRET . Summing up, we find that the AGILE 3PR resolution is better than that of EGRET by a factor of $`2`$ above 400 MeV. These results are remarkable considering the fact that in this preliminary study we are disregarding all the information on the following hit planes, and that we make no use yet of the Mini-Calorimeter information.

## 4 Conclusions

We presented here the Agile REcontruction Method and the results of a preliminary study of its application using only information from the first three hit planes, for which the information on the original photon direction is less dispersed by multiple scattering. The preliminary results of AREM 3PR provide a satisfactory PSF, that turns out to be better than that of EGRET by a factor of $`2`$ at energies above or equal $`400`$ MeV, up to incidence angles of $`\theta =30^{}`$ off-axis. We expect to improve this result by integrating the AREM method with the standard Kalman filter algorithms, to take properly into account the information from all hit planes, and with the Mini Calorimeter data, in order to optimize the track energy determination and the AREM weighting scheme . Moreover, a further improvement of the AGILE PSF should be obtained by applying the AREM method to a model of the AGILE instrument which includes the forthcoming results of a more accurate experimental study of the charge deposition in the Si-microstrips. This analysis will allow an optimal use of the analog readout, which is one of the main characteristics of the AGILE Tracker.

## ACKNOWLEDGMENTS

We thank D.L. Bertsch and D. Thompson for many discussions and collaborative work on EGRET calibration data. We also thank V. Cocco, A. Giuliani, F. Longo, S. Mereghetti, A. Pellizzoni, S. Vercellone and D. Zanello for many interesting discussions and scientific exchanges. We wish to thank the Astroparticle Physics WIZARD Group and INFN Section of Roma II for kind hospitality. C.P. thanks in particular M.P. De Pascale for useful comments after a careful reading of the manuscript, A. Morselli for collaboration at the initial stage of this work and P. Picozza for having set up the Cosmic-Ray research area at the “ Tor Vergata” University of Rome.

Work carried out under the auspices of the Agenzia Spaziale Italiana.


text:
# Spectroscopic Evidence for Starspots on the Secondary Star of SS Cygni

## 1 Introduction

Cool stars with a spectral type later than F5 (cooler than about 6500K) have convective outer layers, which when coupled with rotation, create a dynamo that amplifies magnetic fields and brings them to the stellar surface. Radial magnetic field lines trap the stellar wind, forcing it to co-rotate with the star out to large radii and removing angular momentum from the star. Weber & Davies (1967) first applied such a model to the solar angular momentum loss, and variants of their model can successfully explain the observed rotational spin-down of cool stars in open clusters of increasing age (Krishnamurthi et al., 1997; Bouvier et al., 1997; Tinker et al., 2000; Sills et al., 2000).

Low-mass secondary stars in cataclysmic variables (CVs) should exhibit strong (0.1-0.3T) magnetic fields since they have deep convection zones and are rapidly rotating as a result of being tidally locked with their companion. (Tidal interaction ensures that the secondary star is brought into synchronous rotation with the orbit.) The currently accepted hypothesis is that magnetic braking by the secondary star is responsible for removing angular momentum from the system, perpetuating the mass transfer (Verbunt & Zwaan, 1981; Rappaport et al., 1983). Magnetic braking is also crucial in the standard explanation for the dearth of CVs with orbital periods between 2 and 3 hours (the “period gap”). As CVs evolve towards shorter orbital periods the secondary star becomes less massive, until finally, at masses of 0.25M and periods around 3 hours, it should be fully convective (Robinson et al., 1981). At this point it is said that either magnetic activity may cease (Spruit & Ritter, 1983; Rappaport et al., 1983), or perhaps the magnetic field changes topology (Taam & Spruit, 1989), and thus magnetic braking becomes ineffectual. In turn, mass transfer will cease, and the secondary star, which is slightly bloated due to the effect of mass loss, shrinks to within its Roche-lobe. Gravitational radiation is now the only angular momentum loss mechanism, and only when the binary reaches a period of two hours does mass transfer begin again.

There are two worrying problems with this scenario. First that it relies on applying single star models to semi-detached binaries at somewhat shorter rotation periods. The danger of this assumption is underlined by the case of the W UMa binaries, whose magnetic fields are a factor of order four or five below those for analogous single stars (Stȩpień et al., 2001). The second problem is the lack of good observational evidence for magnetic fields in CV secondaries. Strong (0.1-0.3T) magnetic fields covering a large fraction of the stellar surface are well established in rapidly rotating single and binary G to M-type stars (e.g. Saar & Linsky, 1985). These fields can be detected directly by Zeeman splitting (Johns-Krull & Valenti, 1996) and inferred indirectly from the presence of cool starspots on the stellar surface, or from the non-radiative magnetic heating in outer atmospheres diagnosed by coronal X-ray emission (O’Neal et al., 1996; Haisch & Schmitt, 1996; Strassmeier & Rice, 1998). Applying these techniques to CVs is difficult, due in part to the variability of the accretion disc masking that of the secondary star, and in part due to their intrinsic faintness. There are, however, other indications that spots and/or magnetic activity may be present in CVs. Magnetic activity cycles are a plausible explanation for changes seen on timescales of decades in the overall brightness, inter-outburst periods and orbital periods of CVs. A critical account of these is given in Richman et al. (1994), and given that other explanations are possible, they do not amount to a strong case for magnetic activity. More convincing is to use starspots as an explanation for subtle differences in secondary star line strengths, either in comparison with single stars of similar spectral type, or as a function of orbital phase (e.g. Howell et al., 2000; Dhillon et al., 2002). However, there are always other possible explanations for these effects, including light from the accretion disc, or the peculiar evolutionary history of the secondary star. In summary, the evidence to date for magnetic activity in CV secondary stars amounts to using them as possible post-facto explanations of extant observational data. What is required is an observational test of a prediction based on the idea that the secondary stars have strong magnetic fields.

## 2 The technique

Starspots are the key to the technique we present here. Their presence can be deduced from rotational modulation of stellar brightness or more directly from Doppler imaging of the surfaces of rapidly rotating stars (e.g. Barnes et al., 2001). Donati & Collier Cameron (1997) have used Zeeman Doppler imaging to show that starspots do indeed coincide with regions of strong magnetic field, and starspots are now known to exist in both single cool stars and non-interacting binaries which have rotation periods from about 12 hours to a few weeks.

O’Neal et al. (1996) and O’Neal et al. (1998) have developed a new way to measure the spot area and temperature, which relies on a typical starspot being approximately 1000K cooler than the unspotted photosphere. The total spectrum from the star is then the flux-weighted mean from the spotted and unspotted areas. In particular, when one observes magnetically active K stars, molecular absorption features due to TiO are seen, which can only come from much cooler M-type atmospheres. By modeling the total spectrum as the sum of template K-type and M-type stars, an estimate of the spot area and temperature is recovered, assuming that the star is uniformly covered by spots. In this letter we present the first results using a similar method for a CV secondary star.

Although the application of this technique to CVs is straightforward in principle, there are observational complications which make it difficult in practice. First, unlike O’Neal et al. (1996), we cannot normalise the spectra such that the continuum flux is unity, as the presence of non-stellar components (e.g. the accretion disc) whose flux varies with wavelength, would make the relative fluxes in the lines incorrect. Thus we must flux our spectra, ensuring that the relative levels across the spectra are accurate to approximately one percent. It is well known that normal techniques can introduce problems at this level, and worse still the TiO bands are contaminated by absorption from the Earth’s atmosphere. However, preserving this flux information means that our fitting technique can utilize the overall continuum shape of the template spectra, in addition to the strengths and widths of the absorption lines. For the latter information to be useful, very high signal-to-noise spectra are required.

## 3 Observations and Results

We chose to investigate the 6.6 hour orbital period CV, SS Cygni. Since the secondary star in SS Cyg is an early K star (Harrison et al. (2000) and references therein), the mere detection of M-type TiO features would be a positive indication of spot activity on its surface. Four SS Cyg spectra totaling 24 mins exposure time were taken within 0.3 in phase of superior conjunction of the white dwarf (i.e. when the secondary star is at its closest to the Earth), using the Intermediate Dispersion Spectrograph on the 2.5m Isaac Newton Telescope. With the 1.6 arcsec slit used, the spectrograph yielded a resolution of 2.5Å. The spectra covered the wavelength range 6820-8190Å at a dispersion of 1.22Å per pixel. However, we only used the region 7000-7590Å, since this was well corrected for telluric absorption, and contains the TiO band of interest. The reduced spectra were shifted so that the spectral lines lay at their rest wavelength and added, using a weighted mean, to give a spectrum with a signal-to-noise of about 200. We also observed 35 single K and M stars to use as template spectra, which were also shifted so that the spectral lines lay at their rest wavelength positions. We ensured good relative flux calibration by observing an O or B-star close to each target, after the target had been observed. This hot star was itself fluxed (using data from flux standards), and the resulting spectrum (excluding telluric absorption bands) was fitted with a quadratic (a good approximation to the flux distribution of hot stars in this region of the spectrum), which removed any residual fluxing errors on scales of less than 300Å. We next divided the unfluxed spectrum of the target by the unfluxed spectrum of the O-star, thus removing the telluric absorption lines. Finally we multiplied this spectrum by the quadratic fit to the hot star spectrum to obtain a fluxed spectrum of the target.

Our mean SS Cyg spectrum is shown in Figure 1. There are two strong emission lines at 7065Å and 7280Å (both HeI), whose broadness implies they originate in the accretion disc. (They are broadened both by the rotation of the disc, and our shifting of the spectra into the rest frame of the secondary star.) The rest of the spectrum is a host of weak, rotationally broadened absorption lines from the secondary star. It is immediately evident, even without modeling, that there is a broad absorption trough between 7100Å and 7200Å, which is normally associated with TiO (although the bandhead itself, at 7055Å is obscured by an emission line). However, the secondary star in this CV is of spectral type early K (4800-5300K), and is thus too hot for TiO to exist in significant quantity. Indeed TiO only appears in stars of spectral type K7 and later, where the temperature is below 4300K. It would appear, therefore, that we are observing a K-type star, with cooler, M-type regions, which are presumably starspots.

We modeled the mean SS Cyg spectrum using varying proportions of our template K and M, giant and dwarf stars, with the addition of a shallow quadratic to represent the non-stellar components. The spectral lines of the template stars were broadened to match the lines of SS Cyg which are rotationally broadened to about 90km s<sup>-1</sup> (Martinez-Pias et al., 1994). We restricted the fit to the region 7000-7590Å, since it contained one TiO band and was well corrected for telluric absorption. To begin with we tried a fit which contained just one template star, and the quadratic (Figure 1). The best fitting spectral type was K4V. When we added a second star (Figure 2) $`\chi _\nu ^2`$ fell by a factor of 2.6. (There were 126 degrees of freedom for the two star fit, which returned a $`\chi _\nu ^2`$ of 3.9.) The two star model not only fits the area around TiO well, it also fits the region around 7400Å better than the single star model. This second fit gave spectral type of K4V (4800K) for the unspotted parts of the photosphere, and a spot spectral type M2V (3900K). The hemisphere we observed had approximately 22% of its area covered in spots. To assess the errors in these figures we fitted a two star model to the individual spectra which had been combined to create our mean SS Cyg spectrum. The resulting scatter implied our error is about a sub-type for both the “normal” and spotted spectral types. About 50% of the light was calculated to come from the secondary star. Surprisingly there was little variation in the implied spot coverage with changes in either the normal or spotted spectral type. This is because the fall in flux (as opposed to equivalent width) across the TiO bandhead changes very little with spectral type, when calculated per unit area of atmosphere. The fall in continuum flux being matched almost exactly by the increasing depth of the absorption. Using the data of Jones et al (1994) and Reid et al. (1995) we find a 10 percent change in the flux drop across the TiO bandhead between M2 and M4.

To prove the method does not have any propensity to introduce TiO features where there should not be any, we performed the same experiment on the known spotted star HD 82443. In this star the temperature of the unspotted areas is known to be 4925K, with a spot temperature of 4300K (Messina et al., 1999); thus this object should not show detectable TiO. We fitted the spectrum of this star in an identical fashion to SS Cyg, allowing both spectral types to run free. We found no evidence for TiO, and the fit returned a figure of less than 5% for the non-stellar contribution.

## 4 Caveats

The simplest explanation of our observations is that we have observed starspots on the surface of the secondary star in SS Cyg. However, there are three other obvious possibilities which should be considered: gravity darkening of the secondary star; a peculiar abundance pattern, or an effect caused by the non-stellar component.

To examine the possibility that gravity darkening may be responsible for the cooler region(s), we modeled a star with a pole temperature of 4800K, observed at an inclination of 60 degrees (probable for SS Cyg), using the code of Shahbaz et al (1993). Gravity darkening will have its greatest effect at phase 0.5, when the inner Lagrangian point is facing the observer. This was our reason for preferring to observe phases close to zero, but to ensure our modeling is conservative, we modeled phase 0.3, the closest that we observed to 0.5. We found that only four percent of the flux weighted area lay below the temperature of 4500K where TiO first becomes visible. Only if the TiO band depth at this and lower temperatures were 100 percent, could our observed depth of 4 percent be explained. In fact the depth at these temperatures is only a few percent.

Peculiar abundances are also unlikely to cause the observed effect. The problem is not merely changing the relative strengths of certain lines, but producing simultaneously the features of an early K star and TiO, which only occurs in cooler atmospheres. Modeling a range of abundance possibilities to prove this lies outside the scope of this letter, but to be credible, such modeling would have to be backed with an explanation as to why the star has the implied abundances.

Finally, could the subtraction of the non-stellar component cause the TiO features? Although the exact nature of the non-stellar component is controversial, it does seem to be smooth on scales of 1000Å and more (Harrison et al., 2000). Given that our residuals are also smooth, it seems that only if the non-stellar component itself has TiO could it explain our results.

## 5 Conclusions

Our conclusion is therefore that the most convincing explanation of our data is that starspots do exist on the surfaces of CV secondary stars. These are likely to be associated with equipartition magnetic fields of strength 0.1-0.3T. This puts the idea that magnetically controlled wind braking is responsible for CV evolution on a much firmer footing and opens the way for further studies of these magnetic fields in ways similar to those being pursued in rapidly rotating single stars. It also shows that magnetic fields and spot activity do not disappear at extremely high rotation rates, that are an order of magnitude faster than the claimed dynamo or angular momentum saturation threshold seen in the K stars of young open clusters (Stauffer et al., 1994). However, much crucial work still needs to be done, including establishing whether spots disappear at rotation periods of 2-3 hours, thus providing a natural explanation for the CV period-gap. CV secondaries at these periods are themselves M-type stars and so the mere detection of TiO bands will not be sufficient to prove the presence of starspots.

The Isaac Newton Telescope is operated on the island of La Palma by the Isaac Newton Group in the Spanish Observatorio del Roque de los Muchachos of the Instituto de Astrofisica de Canarias. This work made use of Peter van Hoof’s Atomic Line List (http://www.pa.uky.edu/ peter/atomic). The majority of this work was carried out when NAW was supported by a PPARC studentship at Keele University, and TN was in receipt of a PPARC Advanced Fellowship at the same institute. We are grateful to Steve Saar for comments on an earlier version of this paper.


text:
# Some Observational Consequences of Brane World Cosmologies

## I INTRODUCTION

Recent results based on the magnitude-redshift relation for extragalactic type Ia supernovae (SNe Ia) suggest an accelerating universe dominated by some kind of negative-pressure dark component, the so-called *quintessence* perlmutter ; riess . The existence of this dark energy has also been confirmed, independently of the SNe Ia analyses, by combining the latest galaxy clustering data with cosmic microwave background (CMB) measurements efs . Together, these results seem to provide the remaining piece of information connecting the inflationary flatness prediction ($`\mathrm{\Omega }_\mathrm{T}=1`$) with astronomical observations. Such a state of affairs has stimulated the interest for more general models containing an extra component describing this dark energy and, simultaneously, accounting for the present accelerated stage of the Universe. The simplest and most theoretically appealing possibility is the vacuum energy (cosmological constant). Because of their observational successes, flat models with a relic cosmological constant are considered nowadays our best description of the observed Universe. However, we face at least a serious problem when one considers a nonzero vacuum energy: in order to dominate the dynamics of the Universe only at recent times, a very small value for the cosmological constant ($`\mathrm{\Lambda }_o10^{56}\mathrm{cm}^2`$) is required from observations, while naive estimates based on quantum field theories are 50-120 orders of magnitude larger, thereby originating an extreme fine tunning problem weinberg ; sahni  or making a complete cancellation (from an unknown physical mechanism) seem more plausible.

If the cosmological constant is null, something else must be causing the Universe to speed up. Several possible dark energy candidates have been discussed in the literature, e.g., a vacuum decaying energy density, or a time varying $`\mathrm{\Lambda }`$-term ozer , a relic scalar field peebles , an extra component, the so-called “X-matter”, which is simply characterized by an equation of state $`p_x=\omega \rho _x`$, where $`1\omega <0`$ turner  and that includes, as a particular case, models with a cosmological constant ($`\mathrm{\Lambda }`$CDM) or still models based on the framework of brane-induced gravity dvali ; dvali1 ; deff ; deff1 ; deffZ . In this last case, the basic idea is that our 4-dimensional Universe would be a surface or a brane embedded into a higher dimensional bulk space-time to which gravity can spread. Despite the fact that there is no cosmological constant on the brane, such scenarios explain the observed acceleration of the Universe because the bulk gravity sees its own curvature term on the brane as a negative-pressure dark component and accelerates the Universe deff1 .

Brane world cosmologies have been discussed in different contexts. For example, the issue related to the cosmological constant problem has been addressed ccp as well as the evolution of cosmological perturbations in the gauge-invariant formalism brand , cosmological phase transitions cpt , inflationary solutions cpt1 , baryogenesis dvali99 , stochastic background of gravitational waves hogan1 , singularity, homogeneity, flatness and entropy problems aaa , among others (see hogan  for a discussion on the different perspectives of brane world models). In the observational front, some analyses deffZ  have shown that such models are in agreement with the most recent cosmological observations (see, however, avelino ; dnew ). In this case, constraints from SNe Ia + CMB data require a flat universe with $`\mathrm{\Omega }_\mathrm{m}=0.3`$ and $`\mathrm{\Omega }_{\mathrm{r}_\mathrm{c}}=0.12`$, where $`\mathrm{\Omega }_{\mathrm{r}_\mathrm{c}}`$ is the density parameter associated with the crossover distance between the 4-dimensional and 5-dimensional gravities (see deff1  for details).

In the present work we focus our attention on these kinds of cosmologies. Following deffZ , we study models based on the framework of the brane-induced gravity of Dvali et al. dvali ; dvali1  that have been recently proposed in Refs. deff ; deff1 . We will also consider only the case in which the bulk space-time is 5-dimensional. Our main purpose is to investigate some observational consequences of such scenarios with emphasis on the constraints provided by observations of the angular size of high-$`z`$ milliarcsecond radio-sources.

We structured this paper as follows. In Section II the basic field equations and distance formulas are presented. We also briefly discuss the predicted age of the Universe. In Section III we use measurements of the angular size of high-$`z`$ milliarcsecond radio sources gurv1  to constrain the free parameters of the model. We show that a good agreement between theory and observation is possible if $`\mathrm{\Omega }_\mathrm{m}0.38`$, $`\mathrm{\Omega }_{\mathrm{r}_\mathrm{c}}0.29`$ and $`\mathrm{\Omega }_\mathrm{m}0.09`$, $`\mathrm{\Omega }_{\mathrm{r}_\mathrm{c}}0.29`$ ($`68\%`$ c.l.) for values of the characteristic length of the sources between $`l20h^130h^1`$ pc, respectively. In particular we find that a slightly closed, accelerated model with $`\mathrm{\Omega }_\mathrm{m}=0.06`$, $`\mathrm{\Omega }_{\mathrm{r}_\mathrm{c}}=0.28`$, and $`l=27.06h^1`$ pc is the best fit for these data.

## II Field equations, distance formulas and the age of the universe

The geometry of our 4-dimensional Universe is described by the Friedmann-Robertson-Walker (FRW) line element ($`c=1`$)

$$ds^2=dt^2R^2(t)\left[\frac{dr^2}{1kr^2}+r^2(d\theta ^2+\mathrm{sin}^2\theta \mathrm{d}\varphi ^2)\right],$$
(1)
where $`k=0`$, $`\pm 1`$ is the curvature parameter of the spatial section, $`r`$, $`\theta `$, and $`\varphi `$ are dimensionless comoving coordinates, and $`R(t)`$ is the scale factor. The Friedmann’s equation for the kind of models we are considering is deff1 ; deffZ

$$\left[\sqrt{\frac{\rho }{3M_{pl}^2}+\frac{1}{4r_c^2}}+\frac{1}{2r_c}\right]^2=H^2+\frac{k}{R(t)^2},$$
(2)
where $`\rho `$ is the energy density of the cosmic fluid, $`H`$ is the Hubble parameter, $`M_{pl}`$ is the Planck mass, and $`r_c=M_{pl}^2/2M_5^3`$ is the crossover scale defining the gravitational interaction among particles located on the brane ($`M_5`$ is the 5-dimensional reduced Planck mass). For distances smaller than $`r_c`$ the force experienced by two punctual sources is the usual 4-dimensional gravitational $`1/r^2`$ force whereas for distances larger than $`r_c`$ the gravitational force follows the 5-dimensional $`1/r^3`$ behavior comm .

Equation (2) implies that the normalization condition is given by

$$\mathrm{\Omega }_k+\left[\sqrt{\mathrm{\Omega }_{\mathrm{r}_\mathrm{c}}}+\sqrt{\mathrm{\Omega }_{\mathrm{r}_\mathrm{c}}+\mathrm{\Omega }_\mathrm{m}}\right]^2=1$$
(3)
where $`\mathrm{\Omega }_\mathrm{m}`$ and $`\mathrm{\Omega }_k`$ are the matter and curvature density parameters, respectively and

$$\mathrm{\Omega }_{\mathrm{r}_\mathrm{c}}=1/4r_c^2H^2,$$
(4)
is the density parameter associated to the crossover radius $`r_c`$.

The deceleration parameter, usually defined as $`q_o=R\ddot{R}/\dot{R}^2|_{t_o}`$, now takes the following form

$`q_o`$ $`=`$ $`{\displaystyle \frac{3}{2}}\mathrm{\Omega }_\mathrm{m}\left[1+\sqrt{{\displaystyle \frac{\mathrm{\Omega }_{\mathrm{r}_\mathrm{c}}}{\mathrm{\Omega }_{\mathrm{r}_\mathrm{c}}+\mathrm{\Omega }_\mathrm{m}}}}\right]`$
$`\left[\sqrt{\mathrm{\Omega }_{\mathrm{r}_\mathrm{c}}}+\sqrt{\mathrm{\Omega }_{\mathrm{r}_\mathrm{c}}+\mathrm{\Omega }_\mathrm{m}}\right]^2.`$
For $`\mathrm{\Omega }_k=0`$ (flat case), the above expression reduces to

$$q_o=\frac{3}{2}\mathrm{\Omega }_\mathrm{m}\left[1+\sqrt{\frac{\mathrm{\Omega }_{\mathrm{r}_\mathrm{c}}}{\mathrm{\Omega }_{\mathrm{r}_\mathrm{c}}+\mathrm{\Omega }_\mathrm{m}}}\right]1.$$
(6)
Figure 1 shows the behavior of the deceleration parameter as a function of redshift for selected values of $`\mathrm{\Omega }_\mathrm{m}`$ and $`\mathrm{\Omega }_{\mathrm{r}_\mathrm{c}}`$. As explained above, although there is no cosmological constant on the brane, brane-world models allows periods of accelerated expansion because the bulk gravity sees its own curvature term on the brane as a negative-pressure dark component deff1 . The best fit $`\mathrm{\Lambda }`$CDM case is also showed for the sake of comparison. Note that at late times brane-cosmologies with $`\mathrm{\Omega }_{\mathrm{r}_\mathrm{c}}=0.3`$ accelerates slower than $`\mathrm{\Lambda }`$CDM models with $`\mathrm{\Omega }_\mathrm{\Lambda }=0.7`$ and the same value of $`\mathrm{\Omega }_\mathrm{m}`$. For the best fit model found in Ref. deffZ , i.e., $`\mathrm{\Omega }_\mathrm{m}=0.3`$ and $`\mathrm{\Omega }_{\mathrm{r}_\mathrm{c}}=0.12`$, the accelerated expansion begins at $`z0.5`$ whereas for $`\mathrm{\Lambda }`$CDM we find $`z0.7`$. For our best fit, found in Section 3, we see that the Universe always accelerates at a faster rate than the best fit $`\mathrm{\Lambda }`$CDM model. In this case, the Universe begins to accelerate at $`z2.3`$ (see also sahini  for a more detailed discussion on this topic).

From Eqs. (1) and (2), it is straightforward to show that the comoving distance $`r(z)`$ is given by

$$r(z)=\frac{1}{R_oH_o|\mathrm{\Omega }_k|^{1/2}}\left[|\mathrm{\Omega }_k|^{1/2}_x^{}^1\frac{dx}{x^2f(\mathrm{\Omega }_j,x)}\right],$$
(7)
where the subscript $`o`$ denotes present day quantities, $`x^{}=\frac{R(t)}{R_o}=(1+z)^1`$ is a convenient integration variable and the function $`(r)`$ is defined by one of the following forms: $`(r)=\text{sinh}(r)`$, $`r`$, and $`\text{sin}(r)`$, respectively, for open, flat and closed geometries. The dimensionless function $`f(\mathrm{\Omega }_j,x)`$ takes the following form:

$$f(\mathrm{\Omega }_j,x)=\left[\mathrm{\Omega }_kx^2+\left(\sqrt{\mathrm{\Omega }_{\mathrm{r}_\mathrm{c}}}+\sqrt{\mathrm{\Omega }_{\mathrm{r}_\mathrm{c}}+\mathrm{\Omega }_\mathrm{m}x^3}\right)^2\right]^{1/2},$$
(8)
where $`j`$ stands for $`m`$, $`r_c`$ and $`k`$.

Similarly, the predicted age of the Universe as a function of the redshift can be written as

$$t_z=H_o^1_0^x^{}\frac{dx}{xf(\mathrm{\Omega }_j,x)}.$$
(9)
As one may check from Eqs. (2), (5), (7) and (9), for $`\mathrm{\Omega }_{\mathrm{r}_\mathrm{c}}=0`$, the standard relations are recovered. In Fig. 2 we show the dimensionless age parameter $`H_ot_o`$ as a function of $`\mathrm{\Omega }_\mathrm{m}`$ for several values of $`\mathrm{\Omega }_{\mathrm{r}_\mathrm{c}}`$. Note that for a fixed value of $`\mathrm{\Omega }_\mathrm{m}`$ the predicted age of the Universe is larger for larger values of $`\mathrm{\Omega }_{\mathrm{r}_\mathrm{c}}`$, thereby showing, similarly to what happens in the $`\mathrm{\Lambda }`$CDM context, that the class of models studied here is efficient to solve the “already” classical age of the Universe problem. For example, if $`\mathrm{\Omega }_\mathrm{m}=0.3`$, as sugested by dynamical estimates on scales up to about $`2h^1`$ Mpc calb , and $`\mathrm{\Omega }_{\mathrm{r}_\mathrm{c}}=0.15`$ we find $`t_o13`$ Gyr ($`H_o=70\mathrm{K}\mathrm{m}.\mathrm{s}^1.\mathrm{Mpc}^1`$) or, in terms of the age parameter, $`H_ot_o0.93`$, a value that is compatible with the most recent age estimates of globular clusters carreta ; chab  as well as very close to some determinations based on SNe Ia data riess ; tonry .

## III Constraints from high-$`z`$ angular size measurements

In this section we study the constraints from angular size measurements of high-$`z`$ radio sources on the free parameters of the model. In the following we briefly outline our main assumptions for this analysis. Our approach is based on Ref. alcaniz .

The angular size-redshift relation for a rod of characteristic length $`l`$ can be written as sand

$$\theta (z)=\frac{D(1+z)}{r(z)}.$$
(10)
In the above expression $`D=100lh`$ is the angular size scale expressed in milliarcsecond (mas) for $`l`$ measured in parsecs (compact sources).

In order to constrain the parameters $`\mathrm{\Omega }_\mathrm{m}`$ and $`\mathrm{\Omega }_{\mathrm{r}_\mathrm{c}}`$ we use the angular size data for milliarcsecond radio sources recently compiled by Gurvits et al. gurv1 . This data set, originally composed by 330 sources distributed over a wide range of redshifts ($`0.011z4.72`$), was reduced to 145 sources with spectral index $`0.38\alpha 0.18`$ and total luminosity $`Lh^210^{26}`$ W/Hz in order to minimize any possible dependence of angular size on spectral index and/or linear size on luminosity. This new subsample was distributed into 12 bins with 12-13 sources per bin. Two points, however, should be stressed before discussing the resulting diagrams. First of all, the determination of cosmological parameters is strongly dependent on the characteristic length $`l`$ (see, e.g., alcaniz ). In the absence of a statistical study describing the intrinsic length distribution of the sources, we follow alcaniz  and, instead of assuming a specific value for the mean projected linear size, we have worked on the interval $`l20h^130h^1`$ pc, i.e., $`lO(40)`$ pc for $`h=0.65`$, or equivalently, $`D=1.42.0`$ mas (see also gurv  for a detailed discussion on this topic). Second, following Kellermann kelle , we assume that compact radio sources are free of the evolutionary and selection effects that have bedevilled attempts to use extended double radio source in this context (see, for exemple, buc ), as they are deeply embedded in active galact nuclei, and, therefore, their morphology and kinematics do not depend considerably on the changes of the intergalactic medium. Moreover, these sources have typical ages of some tens of years, i.e., it is reasonable to suppose that a stable population is estabilished, characterized by parameters that do not change with the cosmic epoch jack .

Following a procedure similar to that described in alcaniz , we determine the cosmological parameters $`\mathrm{\Omega }_\mathrm{m}`$ and $`\mathrm{\Omega }_{\mathrm{r}_\mathrm{c}}`$ through a $`\chi ^2`$ minimization for a range of $`\mathrm{\Omega }_\mathrm{m}`$ and $`\mathrm{\Omega }_{\mathrm{r}_\mathrm{c}}`$ spanning the interval  in steps of 0.02,

$$\chi ^2(l,\mathrm{\Omega }_\mathrm{m},\mathrm{\Omega }_{\mathrm{r}_\mathrm{c}})=\underset{i=1}{\overset{12}{}}\frac{\left[\theta (z_i,l,\mathrm{\Omega }_\mathrm{m},\mathrm{\Omega }_{\mathrm{r}_\mathrm{c}})\theta _{oi}\right]^2}{\sigma _i^2},$$
(11)
where $`\theta (z_i,l,\mathrm{\Omega }_\mathrm{m},\mathrm{\Omega }_{\mathrm{r}_\mathrm{c}})`$ is given by Eqs. (7) and (10) and $`\theta _{oi}`$ is the observed values of the angular size with errors $`\sigma _i`$ of the $`i`$th bin in the sample.

In Fig. 3 we show the binned data of the median angular size plotted as a function of redshift for several values of $`\mathrm{\Omega }_\mathrm{m}`$ and $`\mathrm{\Omega }_{\mathrm{r}_\mathrm{c}}`$. For comparison we also show the standard prediction (thick line). Fig. 4 displays the $`95\%`$ and $`68\%`$ c.l. limits from angular size data on the $`\mathrm{\Omega }_\mathrm{m}\mathrm{\Omega }_{\mathrm{r}_\mathrm{c}}`$ plane for the interval $`l20h^130h^1`$ pc. Note that the limits on the plane are more restrictive for increasing values of the characteristic length $`l`$. It happens because for $`z2`$ (where most of our data points are concentrated) the parameter $`\mathrm{\Omega }_{\mathrm{r}_\mathrm{c}}`$ has a behavior similar to a cosmological constant or quintessence, i.e., it increases the distance between two different redshifts. In this way, according to Eq. (10), for the same $`\theta _{oi}`$’s the larger the value of $`l`$ the larger the value of $`r(z)`$ that is required or, equivalently, the smaller the value of $`\mathrm{\Omega }_\mathrm{m}`$. For $`l=20.58h^1`$ pc (D = 1.4 mas) the peak of likelihood is located at $`\mathrm{\Omega }_\mathrm{m}=0.22`$ and $`\mathrm{\Omega }_{\mathrm{r}_\mathrm{c}}=0.18`$. This assumption provides $`\mathrm{\Omega }_\mathrm{m}0.38`$ and $`\mathrm{\Omega }_{\mathrm{r}_\mathrm{c}}0.29`$ at 1$`\sigma `$. In the subsequent panels of the same figure similar analyses are displayed for $`l=23.53h^1`$ pc (D = 1.6 mas), $`l=26.47h^1`$ pc (D = 1.8 mas) and $`l=29.41h^1`$ pc (D = 2.0 mas). For an analysis independent of the choice of the characteristic length $`l`$, i.e., minimizing Eq. (11) for $`\mathrm{\Omega }_\mathrm{m}`$, $`\mathrm{\Omega }_{\mathrm{r}_\mathrm{c}}`$ and $`l`$, we obtain $`\mathrm{\Omega }_\mathrm{m}=0.06`$, $`\mathrm{\Omega }_{\mathrm{r}_\mathrm{c}}=0.28`$ and $`l=27.06h^1`$ pc (D = 1.84 mas) as the best fit for these data with $`\chi ^2=4.25`$ and 9 degrees of freedom. We also remark that although not discussed here it is possible to determine the influence of $`\mathrm{\Omega }_{\mathrm{r}_\mathrm{c}}`$ on the critical redshift $`z_m`$ at which the angular size takes its minimal value. However as shown elsewhere jailS , this test cannot discriminate among world models since different scenarios provide similar values of $`z_m`$.

An elementary combination of our best fit with Eq. (4) enables us to estimate $`r_c`$ (the crossover distance between 4-dimensional and 5-dimensional gravity) in terms of the Hubble radius $`H_o^1`$. One obtains,

$$r_c0.94H_o^1.$$
(12)
Such a value is slightly different from that one found by Deffayett et al. deff  using SNe Ia and CMB data. In their analysis it was found as a concordance model for these two tests a flat model with $`\mathrm{\Omega }_\mathrm{m}=0.3`$ and $`\mathrm{\Omega }_{\mathrm{r}_\mathrm{c}}=0.122`$, leading to an estimate of the crossover radius of $`r_c1.4H_o^1`$. Naturally, with new results from different cosmological tests it will be possible to delimit the $`\mathrm{\Omega }_\mathrm{m}\mathrm{\Omega }_{\mathrm{r}_\mathrm{c}}`$ plane more precisely. An analysis on the observational constraints from statistics of gravitational lenses will appear in a forthcoming communication alc .

## IV conclusion

Based on a large body of observational evidence, a consensus is beginning to emerge that we live in a flat, accelerated universe composed of $``$ 1/3 of matter (barionic + dark) and $``$ 2/3 of a negative-pressure dark component. However, since the nature of this dark energy is still not well understood, an important task nowadays in cosmology is to investigate the existing possibilities in the light of the current observational data. In this paper we have focused our attention on some observational aspects of brane world cosmologies. These models, inspired on superstring-M theory, explain the observed acceleration of the Universe through a large scale modification of gravity arising from a gravitational leakage into an extra dimension deffZ . We showed that their predicted age of the Universe is compatible with the most recent age estimates of globular clusters and, therefore, that there is no age crisis in the context of these models. By using a large sample of milliarcsecond radio sources recently updated and extended by Gurvits et al. gurv1  we obtained, as the best fit for these data, a slightly closed, accelerated universe with $`\mathrm{\Omega }_\mathrm{m}=0.06`$ and $`\mathrm{\Omega }_{\mathrm{r}_\mathrm{c}}=0.28`$. Such values lead to an estimate of the crossover radius of $`r_c0.94H_o^1`$.

###### Acknowledgements.

It is my pleasure to acknowledge helpful discussions with Professor J. A. S. Lima and Professor C. J. Hogan. I also thank Professor L. I. Gurvits for sending his compilation of the data and J. V. Cunha for producing Figure 4. This research is supported by the Conselho Nacional de Desenvolvimento Científico e Tecnológico (CNPq - Brazil) and CNPq(620053/01-1-PADCT III/Milenio)


text:
# Brillouin propagation modes in optical lattices: interpretation in terms of nonconventional stochastic resonance

## Abstract

We report the first direct observation of Brillouin-like propagation modes in a dissipative periodic optical lattice. This has been done by observing, in both theoretical and experimental work, a resonant behaviour of the spatial diffusion coefficient in the direction corresponding to the propagation mode with the phase velocity of the moving intensity modulation used to excite these propagation modes. Furthermore, we show theoretically that the amplitude of the Brillouin mode is a nonmonotonic function of the strength of the noise corresponding to the optical pumping, and discuss this behaviour in terms of nonconventional stochastic resonance.

The last decade has witnessed dramatic progress in laser cooling techniques and nowadays in several laboratories around the world atoms are routinely trapped and cooled at very low temperatures and high densities . Most of the current efforts within the cold atoms community are directed to reaching the regime of quantum degeneracy in both bosonic and fermionic samples, in order to investigate the properties of the quantum gases thus obtained, and realizing an atom laser, the matter wave analog of the laser. Cold atomic samples also constitute an ideal system for the study of complex nonlinear phenomena. This turns out to be especially true if the cold atoms are ordered by the light fields in periodic structures, so-called optical lattices . These are obtained by the interference of two or more laser fields: the light imposes its order on the matter via the dipole force , creating a periodic structure of atoms.

Among the most significant studies of nonlinear dynamics in optical lattices, we recall here the observation of mechanical bistability in a strongly driven dissipative optical lattice  and the realization of the kicked rotor and corresponding detection of chaotic motion in a far detuned lattice . Furthermore the macroscopic transport of atoms in an asymmetric optical lattice without the application of external forces has been observed . This corresponds to the realization of an optical motor, i.e. a ratchet for cold atoms, a well-controlable model system for the molecular combustion motor . Brillouin-like propagation modes of atoms in a dissipative optical lattice have also been theoretically studied to explain the nonlinear optical properties of optical lattices . In this Letter, we report on the first direct observation of these modes. Furthermore we will discuss the propagation mechanism associated with these modes, completely different from the one encountered in dense fluids or solid media. Indeed in dilute optical lattices the interaction between the different atoms is completely negligible, therefore the mechanism for the propagation of atoms cannot be ascribed to any sound-wave-like mechanism. While a sound wave corresponds to a propagating density wave without a net transport of atoms, the Brillouin-like resonances analyzed in this work consist of a net motion of the atoms. In fact, the propagation of atoms through the lattice are determined here by the interaction with the light. The light fields determine both the potential wells where the atoms can oscillate, and whose vibrational frequency determines the velocity of the propagation modes, and the escape from the potential wells, which allows the propagation of atoms. We also show theoretically that the amplitude of the Brillouin mode is a nonmonotonic function of the strength of the noise corresponding to the optical pumping, and discuss this behaviour in terms of nonconventional stochastic resonance .

We consider a three dimensional (3D) lin$``$lin near resonant optical lattice, as in previous work . The periodic structure is determined by the interference of four linearly polarized laser beams, arranged as in Fig. 1. This arrangement results in a periodic modulation of the light polarization and light intensity, which produces a periodic modulation of the light shifts of the different ground states of the atoms (optical potentials) . The optical pumping between the different atomic ground states combined with the spatial modulation of the light shifts leads then to the cooling of the atoms  and to their localization  at the minima of the optical potentials, thus producing a periodic array of atoms.

After the cooling phase, characterized by a significant reduction of the atomic kinetic temperature and by the creation of a periodic spatial order, the atoms keep interacting with the light undergoing optical pumping cycles. The optical pumping may transfer an atom from a potential well to a neighbouring one, giving rise to a variety of transport phenomena . Among these, there are modes which correspond to the propagation of atoms through the optical lattice in a given direction. They consist of a sequence in which one half oscillation in a potential well is followed by an optical pumping process to the neighbouring well, and so on. One can estimate their velocity by $`v_i=\lambda _i\mathrm{\Omega }_i/(2\pi )`$ where $`\lambda _i`$ is the lattice constant and $`\mathrm{\Omega }_i/(2\pi )`$ the vibrational frequency in the $`i`$direction . These modes were first identified through Monte-Carlo simulations in Ref.  and shown to produce resonance lines in the nonlinear optical response of optical lattices. However up to now no direct observation of these modes has been reported. This is achieved in the present work by observing a resonant behaviour of the spatial diffusion coefficient in the direction corresponding to the propagation mode with the phase velocity of the moving intensity modulation used to excite these propagation modes.

The modulation scheme for the excitation of the propagation modes is completely analogous to the one used in previous investigations of the nonlinear optical response of optical lattices . An additional laser field linearly polarized along the $`y`$-axis is introduced with the $`z`$-axis as propagation direction. This probe field interferes with the copropagating lattice beams, creating an intensity modulation. The interference pattern consists of two propagating intensity waves moving with phase velocities $`\stackrel{}{v}_j=\stackrel{}{n}_j\delta /|\mathrm{\Delta }\stackrel{}{k}_j|`$ ($`j=1,2`$) with $`\stackrel{}{n}_j=\mathrm{\Delta }\stackrel{}{k}_j/|\mathrm{\Delta }\stackrel{}{k}_j|`$, and $`\mathrm{\Delta }\stackrel{}{k}_j=\stackrel{}{k}_j\stackrel{}{k}_p`$ the difference between the $`j`$-th lattice beam and the probe ($`p`$) wavevectors . Here $`\delta =\omega _p\omega _L`$ is the detuning between the probe ($`\omega _p`$) and the lattice ($`\omega _L`$) frequencies. According to the numerical simulations for the atomic trajectories presented in Ref. , for $`\delta =\pm \mathrm{\Omega }_x`$, the propagation modes along $`x`$ are excited by the driving field, with the atoms effectively dragged by the moving intensity modulation . Intuitively, the dragging of atoms by the two propagating intensity modulations should result in an increase of the spatial diffusion coefficient $`D_x`$ in the $`x`$-direction. Therefore it should be possible to detect these Brillouin propagation modes by monitoring $`D_x`$ as a function of the detuning $`\delta `$. The propagation modes are then revealed by a resonance in $`D_x`$ around $`\delta =\pm \mathrm{\Omega }_x`$. We tested the validity of this reasoning with the help of semiclassical Monte-Carlo simulations . Taking advantage of the symmetry between the $`x`$ and $`y`$ directions (see Fig. 1), we restricted the atomic dynamics in the $`xOz`$ plane. Our calculations are for a $`J_g=1/2J_e=3/2`$ transition, as customary in numerical analysis of Sisyphus cooling, of an atom of mass $`M`$. We expect our 2D calculations to reproduce the dependencies of the different quantities associated with the real 3D atomic dynamics to within a scaling factor corresponding to the difference in dimensionality . In the numerical simulations, we monitored the variance of the atomic position distribution at a given value of the probe field detuning. We verified that the spatial diffusion is normal, i.e. the atomic square displacements $`\mathrm{\Delta }x^2`$ and $`\mathrm{\Delta }z^2`$ increase linearly with time. Accordingly, we derived the spatial diffusion coefficients $`D_x`$ and $`D_z`$ by fitting the numerical data with $`\mathrm{\Delta }x_i^2=2D_{x_i}t`$ $`(x_i=x,z).`$ Results for the spatial diffusion coefficients as functions of the probe detuning $`\delta `$ are shown in Fig. 2. Two narrow resonances, centered approximately at $`\delta =\pm \mathrm{\Omega }_x`$, appear clearly in the spectrum of the diffusion coefficient along the $`x`$-axis. In contrast, $`D_z`$ does not show any resonant behaviour with the driving field detuning. This demonstrates the validity of the detection scheme based on the measurement of the diffusion coefficients.

In the experiment 85Rb atoms are cooled and trapped in a magneto-optical trap (MOT). The MOT laser beams and magnetic field are then suddenly turned off. Simultaneously the four lattice beams are turned on and after $`10`$ ms of thermalization of the atoms in the lattice the probe laser field is introduced along the $`z`$-axis. We studied the transport of atoms in the optical lattice by observing the atomic cloud expansion with a Charge Coupled Device (CCD) camera. The procedure to derive the diffusion coefficients has been described in detail in previous work , and we recall here only the basic idea. For a given detuning of the probe field we took images of the expanding cloud at different instants after the atoms have been loaded into the optical lattice. From the images of the atomic cloud we derived the atomic mean square displacement along the $`x`$\- and $`z`$-axes.

We verified that the cloud expansion corresponds to normal diffusion and derived the diffusion coefficients $`D_x`$ and $`D_z`$. The procedure has been repeated for several different values of the detuning $`\delta `$ of the probe field. Results for $`D_x`$ and $`D_z`$ as functions of $`\delta `$ are shown in Fig. 3. The probe transmission spectrum is also reported to allow the comparison of the position of the resonances in the spectrum of the diffusion coefficients and in that of the probe transmission. We observe two narrow resonances in the diffusion coefficient along the $`x`$-axis centered at $`\delta =\pm \mathrm{\Omega }_x`$. In contrast, no resonant behaviour of $`D_z`$ with $`\delta `$ is observed. This is in agreement with the numerical simulations and constitutes the first direct observation of Brillouin-like propagation modes in an optical lattice.

We turn now to the analysis of the mechanism behind these propagation modes. Brillouin-like propagation modes have been widely studied in condensed matter and dense fluids . However in the present case the mechanism associated with these modes is clearly of a different nature, as in dilute optical lattices the interaction between atoms is negligible and therefore sound-wave-like propagation modes cannot be supported. On the contrary, in a dilute optical lattice the propagation of the atoms is determined by the synchronization of the oscillation within a potential well with the optical pumping from a well to a neighbouring one, as first identified in the numerical analysis of Ref. . This dynamics can be interpreted in terms of noise-induced resonances: the probe field induces a large scale moving modulation of the periodic potential of the four-beam optical lattice, with the optical pumping constituting the noise source which allows transfer from a well to a neighbouring one. It is then natural to investigate the dependence of the amplitude of the Brillouin mode on the strength of the noise, i.e. on the optical pumping rate. We studied, via semiclassical Monte-Carlo simulations, the atomic cloud expansion for a given depth of the potential well at different values of the optical pumping rate $`\mathrm{\Gamma }_0^{}`$, proportional to the rate $`\mathrm{\Gamma }_{esc}^{^{}}`$ of escape from the well . This has been done by varying the lattice intensity $`I`$ and detuning $`\mathrm{\Delta }`$ so as to keep the depth of the potential wells $`U_0I/\mathrm{\Delta }`$ constant while varying $`\mathrm{\Gamma }_0^{}I/\mathrm{\Delta }^2`$. The diffusion coefficient in the $`x`$-direction has been calculated both for a probe field at resonance ($`|\delta |=\mathrm{\Omega }_x`$) and for a probe field far off-resonance ($`|\delta |\mathrm{\Omega }_x`$). The two diffusion coefficients will be indicated by $`D_x`$ and $`D_x^0`$ respectively. To characterize quantitatively the response of the atomic system to a noise strength variation, we introduce the enhancement factor $`\xi `$ defined as

$$\xi =\frac{D_xD_x^0}{D_x^0}.$$
(1)
Numerical results for the enhancement factor $`\xi `$ as a function of the optical pumping rate at a given value of the potential well depth (i.e. for fixed light shift per beam $`\mathrm{\Delta }_0^{}`$) are shown in Fig. 4. At small pumping rates, $`\xi `$ increases abruptly with $`\mathrm{\Gamma }_0^{}`$; then a maximum is reached, corresponding to the synchronization of the oscillation of the atoms within a well with the escape from a well to the neighbouring one; finally at larger pumping rates this synchronization is lost and $`\xi `$ decreases. This dependence recalls the typical behaviour of stochastic resonance , with the noise enhancing the response of the atomic system to the weak moving modulation. It should be noted that the system analyzed here has one important peculiarity with respect to the model usually considered in the analysis of stochastic resonance. Stochastic resonance is in general understood as the noise-induced enhancement of a weak periodic signal with a frequency much smaller than the intrawell relaxation frequency within a single metastable state. In contrast, in the present case, the noise synchronizes precisely with the intrawell motion of the atoms. This corresponds to a nonconventional stochastic resonance scenario .

In summary, in this Letter we introduced a scheme for the detection of Brillouin propagation modes in optical lattices and we reported on their direct observation. Furthermore, we studied via Monte-Carlo simulations the amplitude of the Brillouin mode, as characterized by an increase of the diffusion coefficients due to the presence of the probe field, as a function of the rate of escape from the potential wells. The Brillouin modes examined in this work differ from their counterparts in solid state or dense fluids as they are sustained by a medium of noninteracting particles. From our analysis it turns out that in the presence of noise the Brownian motion of a system of particles in a periodic potential can be turned in a motion at a well-defined velocity by the application of a weak moving modulation. This represents a quite unusual situation in statistical physics and may constitute a model for many biological phenomena, like the transmission of weak signals in neuronal systems .

We thank Franck Laloë for his continued interest in our work and David Lucas for comments on the manuscript. This work was supported by CNRS, the European Commission (TMR network ”Quantum Structures”, contract FMRX-CT96-0077) and by Région Ile de France under contract E.1220. Laboratoire Kastler Brossel is an ”unité mixte de recherche de l’Ecole Normale Supérieure et de l’Université Pierre et Marie Curie associée au Centre National de la Recherche Scientifique (CNRS)”.