text:
# An accelerator mode based technique for studying quantum chaos

## Abstract

We experimentally demonstrate a method for selecting small regions of phase space for kicked rotor quantum chaos experiments with cold atoms. Our technique uses quantum accelerator modes to selectively accelerate atomic wavepackets with localized spatial and momentum distributions. The potential used to create the accelerator mode and subsequently realize the kicked rotor system is formed by a set of off-resonant standing wave light pulses. We also propose a method for testing whether a selected region of phase space exhibits chaotic or regular behavior using a Ramsey type separated field experiment.

Chaos in quantum mechanics is still a relatively poorly defined concept. Typically it is taken to refer to the behavior of a quantum system which in the classical limit exhibits an exponential sensitivity to initial conditions. A much studied example of such a system is the delta-kicked rotor Lichtenberg and Lieberman ; Casati . A realization of its quantum version, in the guise of cold atoms exposed to the ac-Stark shift potential of an off-resonant standing wave of light Raizen 1st QDKR , has elucidated many of the concepts associated with quantum chaos ResonancePRL ; Tunnelling2 ; Quantum-Classical ; Multi-dimensions . However, one substantial problem remains: there is no clear way of distinguishing regular from chaotic dynamics. In the quantum regime it is not possible to define a chaotic region of phase space as being one in which two initially similar states have an overlap which decreases exponentially with time. The unitary nature of any interaction necessarily implies that the overlap between such states remains unchanged. It has been suggested by Peres Peres  that an alternative definition is needed for quantum systems. Peres’ proposal involves examining the evolution of a state which can interact with two potentials which differ very slightly in form. If the overlap between states produced by an evolution under each potential decreases exponentially as a function of time, then the region of phase space occupied by the initial state can be said to be chaotic in the quantum sense. The justification for this definition is that perturbation theory only converges when the region under examination has regular properties Feingold . More recently Gardiner et al. Gardiner97  have proposed an experiment with a trapped ion which would implement the essential features of Peres’ idea.

In this paper we address two problems related to the experimental determination of a phase space stability map for the atom optical version of the QDKR. Firstly, we demonstrate a method for preparing cesium atoms in a restricted region of phase space using the same standing wave light pulses which create the rotor. Secondly, we show how by examining the overlap between the wavefunctions of the two hyperfine levels of the $`6^2S_{1/2}`$ cesium ground state it should be possible to test the type of dynamics exhibited by the prepared atoms.

The basis of our atom optical version of the QDKR is exposure of laser cooled cesium atoms to pulses (duration $`t_\mathrm{p}`$ and separation time $`T`$) of a standing wave of off-resonant light. The pulses are short enough to allow the effect of the atoms’s kinetic energy to be neglected during a pulse. This places our experiment in the Raman-Nath regime in which the spatially periodic ac-Stark shift potential created by the light acts as a thin phase grating for the atoms AccModePRL . Thus an incident plane de Broglie wave is diffracted into a series of “orders” separated in momentum by $`\mathrm{}G`$, where $`G=2k_\mathrm{L}`$ ($`k_\mathrm{L}`$ is the light wavevector). One of the striking features of this system is the existence of quantum resonances ResonancePRL ; newfishman . These resonances occur when the pulse interval is a multiple of the half-Talbot time, $`T_{\frac{1}{2}}=2\pi m/\mathrm{}G^2`$. During these special times all of the diffraction orders formed from incident plane waves with certain momenta will freely evolve phases which are multiples of $`2\pi `$. For cesium, the first quantum resonance occurs at $`T=67`$ $`\mu `$s.

In addition to being a paradigm of experimental quantum chaos, this system can also produce quantum accelerator modes AccModePRL . This is achieved by adding a potential of the form $`U_\mathrm{A}(x)=max`$, where $`m`$ is the atomic mass, $`a`$ is an applied acceleration and $`x`$ is the position along the standing wave. In our experiment the standing wave is oriented in the vertical direction, so $`a`$ is the acceleration due to gravity. Quantum accelerator modes are characterized by a fixed momentum transfer of (on average) $`\xi \mathrm{}G`$ during each standing wave pulse (see Fig. 1 of Ref. AccModePRL ). Typically the accelerator modes are formed when the value of $`T`$ is near to a quantum resonance. For simplicity, we will henceforth confine our attention to pulse repetition times near to the first quantum resonance. One way of modelling the accelerator modes is to make the approximation that they consist of just a few diffraction orders (say $`q1`$, $`q`$ and $`q+1`$, where $`q`$ is an integer) which in the time between two light pulses accumulate a phase difference which is very close to an integer multiple of $`2\pi `$. At the next pulse this makes it possible for interference to occur between the diffraction orders in such a way that the population of the three orders centered on $`q+\xi `$ non-integer am  is enhanced. By using this rephasing condition, it can be shown AccModePRA  that after $`N_\mathrm{p}`$ pulses the central momentum of the accelerator mode is given by $`q=\frac{N_\mathrm{p}}{\gamma }\frac{\alpha ^2}{1\alpha }=N_\mathrm{p}\xi `$, where $`\gamma =\mathrm{}^2G^3/2\pi m^2a`$ and $`\alpha =T/T_{\frac{1}{2}}`$. The same rephasing condition also leads to the conclusion that only incident plane waves with momenta

$$p_{\mathrm{init}}=\left(\alpha ^1l+\mathrm{const}\right)\mathrm{}G$$
(1)
(where $`l`$ is any integer) can ever participate in an accelerator mode. Thus in momentum space an accelerator mode resembles a comb with a tooth spacing of $`\alpha ^1\mathrm{}G`$. Furthermore, since plane waves separated by this momentum behave identically under the action of the kicks, the accelerator mode effectively contains only a single momentum. Figure 1(a) shows the theoretical momentum distribution of the accelerator mode after one set of pulses, as calculated using a model based on diffraction AccModePRL . To reflect the experimental situation the starting distribution was gaussian with a width of 12 $`\mathrm{}k_\mathrm{L}`$ at FWHM. Although there is good qualitative agreement between Eq. (1) and the numerically derived momentum distribution, an explanation for the finite widths of the comb elements is needed. To understand this effect recall that each comb element must give rise to a set of diffraction orders which are always spaced by exact multiples of $`\mathrm{}G`$. To determine what the real momentum distribution looks like we must add together the diffraction from all the different comb elements. These have a spacing of $`\alpha ^1\mathrm{}G`$, which for pulse interval times just less than the Talbot time is slightly greater than $`\mathrm{}G`$. The most obvious effect of including the diffraction orders is that the width of each comb element of the resultant distribution becomes non-zero and increases with $`|\alpha 1|`$. Additionally, when all the diffraction orders are weighted by the intensity of the comb element from which they originated, the nearest-neighbor spacing of the peaks in the momentum distribution is reduced as one moves away from the center of the distribution. The opposite effect is produced when the pulse interval is slightly greater than the Talbot time. Although the change in spacing is not immediately obvious in Fig. 1(a), it has been confirmed in a detailed analysis of this data.

In Fig. 1(b) we show the spatial form of the accelerator mode wavefunction as calculated with our numerical model. Such a distribution is similar to that deduced from the assumption that the mode consists of the sum of three plane waves. Importantly, the wavefunction is periodically localized in position with maxima occurring every $`\lambda _{\mathrm{spat}}=2\pi /G`$. Since this is the wavelength of the standing wave potential, points having this separation behave equivalently and the wavefunction effectively has a spread in position which is less than $`\lambda _{\mathrm{spat}}`$. From Fig. 1 we expect the width of one element of the momentum comb to be approximately 0.4 $`\mathrm{}G`$ and the extent of each region of strong spatial localization to be $`\lambda _{\mathrm{spat}}/3`$. For comparison, the extent in momentum of a unit cell of the classical phase space is $`(T_{1/2}/T)\mathrm{}G`$ while that in position is $`\lambda _{\mathrm{spat}}`$. Hence the accelerator mode isolates a restricted region in phase space which is about 10% of the overall area. Additionally, since the effect of the accelerator mode is to produce a large momentum offset, it should be straightforward to isolate these atoms.

We now discuss the experimental observation of the distributions described above. Our atomic source was a cesium magneto-optic trap containing approximately 107 atoms. After the trap was switched off the atoms were cooled in an optical molasses to a temperature of 5 $`\mu `$K, corresponding to a momentum width of 12 $`\mathrm{}k_\mathrm{L}`$ at FWHM. Following release, the atoms were exposed to a series of pulses from a vertically oriented standing light wave, detuned $`30`$ GHz below the $`F=4F^{}=3`$ line in the D1 transition at 894 nm. The pulses had a duration of $`t_\mathrm{p}=0.5\mu `$s and a peak intensity of approximately 20 W/cm2. The atoms then fell freely until they passed through a probe laser beam resonant with the D2 cycling transition. By measuring the absorption of the probe light as a function of time we were able to determine the momentum distribution of the atoms in the $`F=4`$ level. Additional details of our apparatus can be found in Ref. AccModePRA .

The momentum resolution of our experiment was approximately 2 photon recoils, determined by the initial spatial extent of the optical molasses and the thickness of the probe laser beam. Hence direct observation of the accelerator mode momentum comb was not possible. Instead we used two sets of temporally separated light pulses to infer its existence. The first set contained 20 pulses and was used to create an accelerator mode. The resultant atomic distribution was then translated in momentum by allowing the atoms to fall for a variable amount of time $`t_{\mathrm{wait}}`$. Finally, a second set of pulses, identical to the first, was applied. Figure 2 shows the resultant momentum distributions as a function of $`t_{\mathrm{wait}}`$. Both panels contain a large fraction of unaccelerated atoms near $`p=0`$ (dotted line), a group of atoms that have been accelerated by one set of pulses (dot-dashed line), and atoms that have been accelerated by both pulse sets (dashed line). In Fig. 2(a) the pulses had an interval of 60 $`\mu `$s and the accelerator mode imparted momentum in the negative direction (with the convention gravity is negative), while in Fig. 2(b) the interval was 74 $`\mu `$s and the accelerator mode imparted positive momentum. The effect of $`t_{\mathrm{wait}}`$ is to allow gravity to translate the distribution in momentum space. If any of this distribution overlaps with the momentum comb of Eq. (1) when the second set of pulses occurs then a further acceleration takes place. The periodic variation in the doubly accelerated population with $`t_{\mathrm{wait}}`$ is just the length of time required for gravity to accelerate the atoms by the momentum separation of adjacent comb elements, $`\alpha ^1\mathrm{}G`$. In the T=60 $`\mu `$s case this gives a theoretical value of 753 $`\mu `$s between accelerator mode revivals. Along the dashed line of Fig. 2(a) we observe 763 $`\pm `$ 4 $`\mu `$s, the discrepancy with the calculated value most likely due to the slight variation in comb spacing across the accelerator mode discussed previously. A similar level of agreement is found for T=74 $`\mu `$s. As regards the width of each momentum comb element, we note that numerical simulations have provided almost identical results to those shown in Fig. 2 so it is reasonable to assume that the actual extent of each momentum element is as shown in Fig. 1(a).

Although we have accounted for the accelerator mode revivals, we have said nothing about their sloping orientation (as seen in Fig. 2), nor about this slope’s dependence on pulse interval. An explanation for this effect can be found by returning to the earlier observation that the momentum intervals between the peaks of Fig. 1(a) change as one moves away from the center of the accelerator mode. This implies that not all parts of the distribution have to be translated by the same amount to regain the accelerator mode condition of Eq. (1). At 60 $`\mu `$s, the outer peaks bulge towards the center of the distribution. Assuming that the whole distribution is accelerating under gravity, the first part to regain the accelerator mode condition will be the component with the greatest positive momentum. For 74 $`\mu `$s the opposite will be the case and the direction of the revival’s slope will flip.

When $`t_{\mathrm{wait}}`$ is scanned in much smaller time steps the previous experiment can also be used to demonstrate localization in position. These experiments are only meaningful if they are performed when the wavefunction has some degree of spatial localization, that is when $`t_{\mathrm{wait}}`$ is close or equal to an integer multiple of the pulse separation time. Furthermore, for these wait times all the elements of the momentum comb described by Eq. (1) translate by an integer number of standing wave periods (or equivalently acquire a phase of $`2\pi \times integer`$). As we have seen, the momentum distribution produced by the accelerator mode is not a comb of delta functions. Thus we do not expect an exact rephasing of the distribution even at the special wait times. Experimentally we have found that for wait times greater than three pulse separation times very little rephasing is observed. Figure 3 shows the momentum distributions observed when the pulse separation time is 60 $`\mu `$s and the wait time is scanned near to 120 $`\mu `$s. The accelerator mode reappears at the dashed line whenever the atoms prepared by the first set of pulses have moved a whole number of standing-wave wavelengths. We can quantify what the period of the reappearance or revival should be by determining the distance which the atoms move during $`t_{\mathrm{wait}}`$. This distance is given by $`\mathrm{\Delta }x=v_0t_{\mathrm{wait}}+\frac{1}{2}at_{\mathrm{wait}}^2`$, where $`v_0`$ is the velocity of the atoms after the end of the first set of pulses. We have calculated this velocity to be 177 mm/s revival period , giving a revival period of 2.53 $`\mu `$s. This is in good agreement with the value of $`2.56\pm 0.03`$ $`\mu `$s that is observed along the dashed line in Fig. 3. Another important feature of this figure is that except for the exact multiple of the pulse separation time, not all parts of the accelerator mode produced by the first set of pulses are able to simultaneously undergo acceleration by the second set of pulses. This can be explained by noting that the distance moved by an atom during the wait time depends upon its initial velocity. Thus the revival period for atoms that have been accelerated the most is the smallest. This picture also provides an alternative way of understanding why we only see clear signatures of spatial localization near the multiples of $`T`$; the revivals from the different comb elements in the accelerator mode all come into phase at these times.

The preceding results demonstrate that by using a pulsed periodic potential it is possible to select restricted regions of phase space. We now suggest a technique for investigating the dynamical properties of these regions under the action of a chaotic potential. Our proposal is to take two atomic states which have been prepared by an accelerator mode to have identical external components. We would expose each state to chaotic potentials which differ slightly in strength and then measure the overlap between the states in order to determine the dynamics. The $`F=3`$ and $`F=4`$ ground state hyperfine levels of cesium offer a convenient means of realizing such a procedure through a technique analogous to that of a Ramsey separated field experiment Ramsey . After selecting the desired region of phase space and preparing atoms in the $`F=4`$, $`m_\mathrm{F}=0`$ state, a 50:50 superposition of the two hyperfine levels can be created by employing a $`\frac{\pi }{2}`$ pulse of microwave radiation resonant with the 9.19 GHz $`F=3`$, $`m_\mathrm{F}=0F=4`$, $`m_\mathrm{F}=0`$ magnetic dipole transition. The QDKR light pulses can then be applied using the same laser beam and pulse envelope as the accelerator mode, but with the key change of moving the standing wave to make $`a=0`$ ResonancePRL . Importantly, each of the states in the superposition would experience the standing wave potential with a different strength. For example, light detuned 30 GHz below resonance for atoms in the $`F=4`$ level is detuned approximately 40 GHz for $`F=3`$ atoms. Since the light shift scales inversely with detuning, the component of the superposition in the $`F=3`$ state sees a kicking potential which is $`0.75`$ the size experienced by the component in $`F=4`$. The final step of our proposal would be to expose the superposition to a second $`\frac{\pi }{2}`$ microwave pulse which differs in phase by an amount $`\mathrm{\Delta }\varphi `$ from the first pulse. If the overlap between the two final states is good then a sinusoidal variation in the population of either state would be observed as $`\mathrm{\Delta }\varphi `$ were systematically changed. The exact value of the overlap could be determined by measuring the population at $`\mathrm{\Delta }\varphi =0`$ and $`\mathrm{\Delta }\varphi =\frac{\pi }{2}`$ and then using a method similar to that outlined in Ref. Gardiner97 .

To summarize, we have demonstrated a new method of preparing atoms in a narrow region of phase space from an initially broad distribution. This was achieved using quantum accelerator modes, a process in which atoms are exposed to pulses of an off-resonant standing light wave which is accelerating relative to the atoms’ frame of reference. Since the experimental configuration is very similar to that used for experiments which have studied the QDKR, it should be relatively straightforward to apply this technique to studies of quantum chaos. Previous experiments have concentrated on the ensemble behavior of the atoms. The work presented here paves the way for a more detailed examination of the dynamics. We have also suggested an experiment which would allow the production of a quantum phase space stability map. By using a method based on Ramsey separated fields we can find the amount of overlap between identical initial states which have interacted with chaotic potentials of slightly different strengths. We hope that the combination of these two ideas will eventually lead to a fuller understanding of the concepts of quantum chaos.

###### Acknowledgements.

We thank A. Buchleitner, K. Burnett, S. Fishman, S. A. Gardiner, I. Guarneri, M. K. Oberthaler and S. Wimberger for stimulating discussions. This work was supported by the UK EPSRC, the Paul Instrument Fund of The Royal Society, the EU as part of the TMR ‘Cold Quantum Gases’ network, contract no. HPRN-CT-2000-00125 and the ESF BEC2000+ program.


format: text

text:
# Effect of initial correlations on short–time decoherence

## Appendix A

In this appendix we collect for convenience the exact expressions of the quantities appearing in the propagator (12). The derivation of these quantities can be found in Ref. . The functions $`K(t)`$ and $`N(t)`$ are given by

$`K(t)`$ $`=`$ $`\nu \text{ cotan }\nu t,`$ (34)
$`N(t)`$ $`=`$ $`{\displaystyle \frac{\nu \mathrm{exp}[\gamma t]}{\mathrm{sin}\nu t}},`$ (35)
where $`\nu =(\omega _0^2\gamma ^2)^{1/2}`$. We work in the limit of an underdamped oscillator where the friction coefficient $`\gamma `$ is smaller than the frequency $`\omega _0`$ of the harmonic oscillator. For a free particle, $`\omega _0=0`$, one has to replace $`\nu i\gamma `$. The variance of the position in equilibrium, $`\kappa =x^2_{eq}`$, is given by

$$\kappa =\frac{2\gamma }{\pi }_0^{\mathrm{}}𝑑\omega \frac{\omega \mathrm{coth}(\beta \omega /2)}{(\omega ^2\omega _0^2)^2+4\gamma ^2\omega ^2}.$$
(36)
Its asymptotic behavior at high ($`T\omega _0`$) and low temperature ($`T\omega _0`$) is respectively given by $`\kappa 1/\beta \omega _0^2`$ and $`\kappa 1/2\omega _0`$. The functions $`C(t)`$ and $`E(t)`$ are, on the other hand, of the form

$$f(t)=\frac{\gamma }{\pi }_0^{\omega _c}𝑑\omega \omega \mathrm{coth}(\beta \omega /2)f_\omega (t),$$
(37)
where $`f_\omega (t)`$ has respectively to be replaced by

$`C_\omega (t)`$ $`=`$ $`{\displaystyle \frac{1}{\mathrm{sin}^2\nu t}}{\displaystyle _0^t}𝑑t^{}{\displaystyle _0^t}𝑑t^{\prime \prime }\mathrm{sin}\nu (tt^{})\mathrm{cos}\omega (t^{}t^{\prime \prime })\mathrm{sin}\nu (tt^{\prime \prime })\mathrm{exp}[\gamma (t^{}+t^{\prime \prime })]`$ (38)
$`+`$ $`{\displaystyle \frac{4\gamma \omega ^2}{(\omega _0^2\omega ^2)^2+4\gamma ^2\omega ^2}}{\displaystyle \frac{1}{\mathrm{sin}\nu t}}{\displaystyle _0^t}𝑑t^{}\mathrm{cos}\omega t^{}\mathrm{sin}\nu (tt^{})\mathrm{exp}[\gamma t^{}]`$ (39)
$``$ $`{\displaystyle \frac{2\omega (\omega ^2\omega _0^2)}{(\omega _0^2\omega ^2)^2+4\gamma ^2\omega ^2}}{\displaystyle \frac{1}{\mathrm{sin}\nu t}}{\displaystyle _0^t}𝑑t^{}\mathrm{sin}\omega t^{}\mathrm{sin}\nu (tt^{})\mathrm{exp}[\gamma t^{}]+{\displaystyle \frac{\omega ^2}{(\omega _0^2\omega ^2)^2+4\gamma ^2\omega ^2}},`$ (40)
$`E_\omega (t)`$ $`=`$ $`{\displaystyle \frac{2}{\mathrm{sin}\nu t}}{\displaystyle \frac{1}{(\omega _0^2\omega ^2)^2+4\gamma ^2\omega ^2}}{\displaystyle _0^t}𝑑t^{}\mathrm{sin}\nu (tt^{})\mathrm{exp}[\gamma t^{}]\left((\omega _0^2\omega ^2)\mathrm{cos}\omega t^{}2\gamma \omega \mathrm{sin}\omega t^{}\right).`$ (41)
## Appendix B

In this appendix we evaluate the function $`\stackrel{~}{C}(t)=\mathrm{sin}^2\nu tC(t)`$ \[the contributions coming from the function $`E(t)`$ are of higher order and can therefore be neglected in the limit considered in the paper\]. We begin by computing the first double integral in Eq. (25). Introducing the new variables $`u=(t^{}+t^{\prime \prime })/2`$ and $`v=t^{}t^{\prime \prime }`$, we have

$`a_\omega (t)`$ $`=`$ $`{\displaystyle _0^t}𝑑t^{}{\displaystyle _0^t}𝑑t^{\prime \prime }\mathrm{sin}\nu (tt^{})\mathrm{sin}\nu (tt^{\prime \prime })\mathrm{cos}\omega (t^{}t^{\prime \prime })\mathrm{exp}[\gamma (t^{}+t^{\prime \prime })]`$ (42)
$`=`$ $`{\displaystyle _0^t}𝑑u{\displaystyle _0^t}𝑑v\mathrm{sin}\nu (tuv/2)\mathrm{sin}\nu (tu+v/2)\mathrm{exp}[2\gamma u]\mathrm{cos}\omega v.`$ (43)
The integral over $`u`$ is readily obtained as,

$$_0^t𝑑u\mathrm{sin}\nu (tuv/2)\mathrm{sin}\nu (tu+v/2)e^{2\gamma u}=\frac{1}{4}\left[\frac{(e^{2\gamma t}1)\mathrm{cos}\nu v}{\gamma }+\frac{\gamma \mathrm{cos}2\nu t\nu \mathrm{sin}2\nu t\gamma e^{2\gamma t}}{\gamma ^2+\nu ^2}\right],$$
(44)
while the integral over $`v`$ can be calculated using

$$_0^{\mathrm{}}𝑑\omega \omega \mathrm{coth}(\beta \omega /2)\mathrm{cos}\omega v2T_0^{\mathrm{}}𝑑\omega \mathrm{cos}\omega v=2\pi T\delta (v),$$
(45)
at high temperature, and

$$_0^{\mathrm{}}𝑑\omega \omega \mathrm{coth}(\beta \omega /2)\mathrm{cos}\omega v_0^{\mathrm{}}𝑑\omega \omega \mathrm{cos}\omega v=\frac{1}{v^2},$$
(46)
at low temperature (the latter should be interpreted as a principal value ). Furthermore, by using the following two expressions,

$`{\displaystyle _0^t}𝑑t^{}\mathrm{cos}\omega t^{}\mathrm{sin}\nu (tt^{})e^{\gamma t^{}}={\displaystyle \frac{\nu (\omega _0^2\omega ^2)\mathrm{cos}\nu t\gamma (\omega ^2+\omega ^2)\mathrm{sin}\nu t+\nu e^{\gamma t}[(\omega _0^2\omega ^2)\mathrm{cos}\omega t+2\gamma \omega \mathrm{sin}\omega t]}{(\omega _0^2\omega ^2)^2+4\gamma ^2\omega ^2}},`$ (47)
$`{\displaystyle _0^t}𝑑t^{}\mathrm{sin}\omega t^{}\mathrm{sin}\nu (tt^{})e^{\gamma t^{}}={\displaystyle \frac{2\gamma \nu \mathrm{cos}\nu t+\omega (\gamma ^2\nu ^2+\omega ^2)\mathrm{sin}\nu t+\gamma e^{\gamma t}[2\gamma \omega \mathrm{cos}\omega t+(\omega _0^2\omega ^2)\mathrm{sin}\omega t]}{(\omega _0^2\omega ^2)^2+4\gamma ^2\omega ^2}},`$ (48)
we can rewrite the remaining three terms in Eq. (25) in the compact form

$$b_\omega (t)=\frac{2\omega \mathrm{sin}\omega t\nu \mathrm{sin}\nu t\mathrm{exp}[\gamma t]\omega ^2\mathrm{sin}^2\nu t}{(\omega _0^2\omega ^2)^2+4\gamma ^2\omega ^2}.$$
(49)
The function $`\stackrel{~}{C}(t)`$ is eventually given by

$$\stackrel{~}{C}(t)=\frac{\gamma }{\pi }_0^{\omega _c}𝑑\omega \omega \mathrm{coth}(\beta \omega /2)[a_\omega (t)+b_\omega (t)].$$
(50)


text:
# Extrasolar Terrestrial Planets: Can We Detect Them Already?

## 1. Introduction

All extrasolar planets orbiting main-sequence stars known to date are giant planets ranging in mass from $`m\mathrm{sin}i=0.12\mathrm{M}_{\mathrm{Jupiter}}`$ to $`17\mathrm{M}_{\mathrm{Jup}}`$ (with $`i`$ the unknown angle between the orbital plane and the sky). It is virtually impossible to detect terrestrial planets (rocky objects of 1 to a few $`\mathrm{M}_{}`$) around F, G and K-type stars using the high-precision $`RV`$ technique, since Earth-mass planets induce only negligible reflex-motions on their host stars. This scenario changes in the faint and low-mass regime of the Hertzsprung-Russell diagram. In the case of M-dwarf stars the low stellar primary mass leads to detectable $`RV`$ amplitudes even for planets of a few $`\mathrm{M}_{}`$ and less in short-period orbits. Fig.1 shows the $`RV`$ signatures of planets with $`m\mathrm{sin}i=1`$ to $`2.5\mathrm{M}_{}`$ orbiting an M5V star with a mass of $`0.1\mathrm{M}_{}`$. Due to its intrinsic faintness the habitable zone is located very close to the star (Kasting, Whitmire, & Reynolds 1993). This, too, favors the detection by the $`RV`$ technique, since shorter periods also mean higher $`RV`$ amplitudes. Although these nearby planets are probably tidally locked into synchronous rotation, a 3-D climate study by Joshi, Haberle, & Reynolds (1997) demonstrated that these planets are still likely to be habitable.

## 2. Searching for terrestrial planets inside the habitable zone of <br>Proxima Centauri

For the past two years we have used the ESO VLT and the UVES spectrograph to monitor a sample of M-dwarf stars in the southern hemisphere to search for planetary companions. Our $`RV`$ results for Proxima Cen (M5V, $`V=11.05`$), obtained with an enhanced version of our Austral radial velocity code (Endl, Kürster, & Els 2000), are displayed in Fig.2. Over the time span of more than one year the star experiences a linear acceleration of unknown origin. The rms-scatter around this trend is $`2.5\mathrm{m}\mathrm{s}^1`$.

### 2.1. Detection sensitivity

We then take the 22 RV-measurements (after subtraction of the linear trend) and determine our planet detection threshold in the period range for orbits inside the habitable zone (2 to 16 days), using the technique from Endl at al. (2001). For the primary mass we adopt a value of $`0.11\mathrm{M}_{}`$ after Henry et al. (1999). Fig.3 shows our momentary detection capability based on the $`RV`$ data collected so far. With a confidence of $`>99\%`$ we could have already detected all planets with $`m\mathrm{sin}i`$ values higher than $`4`$ to $`6\mathrm{M}_{}`$ inside the habitable zone. With a chance probability of $`50\%`$ (i.e. only half of the test signals were recovered successfully with a confidence of $`>99\%`$) we would have even found signals of planets down to $`m\mathrm{sin}i2.5\mathrm{M}_{}`$. These new constraints on planets inside the habitable zone of Proxima Cen supplement the existing limits for giant planets around this star by our ESO CES planet search program (Kürster et al. 1999) and the HST Fine Guidance Sensor astrometric results of Benedict et al. (1999).

## 3. Conclusions

Measuring radial velocities of M-dwarfs with the appropriate high level of precision allows us to detect extremely low-mass planets in short period orbits. With the 22 $`RV`$ measurments we obtained so far for Proxima Cen, we could have already detected all planets with $`m\mathrm{sin}i=4`$ to $`6\mathrm{M}_{}`$ inside the habitable zone. So the answer to the title question is: not yet, but we get pretty close, and our sensitivity will further improve by extending the monitoring time span.

#### Acknowledgments.

This work is supported by NASA Grant NAG5-9227 and by NSF Grant AST-9808980. We thank the ESO OPC and DDTC for generous allocation of observing time to our VLT M-dwarf planet search program (Proposals 69.C-0722, 68.C-0415, 67.C-5700).

## References

Benedict G.F., McArthur B., Chappell D.W., et al. (13 authors) 1999, AJ, 118, 1086

Endl M., Kürster M., & Els S. 2000, A&A, 362, 585

Endl M., Kürster M., Els S., Hatzes A.P., & Cochran W.D. 2001, A&A, 374, 675

Henry T.J., Franz O.G., Wasserman L.H., Benedict G.F., Shelus P.J., Ianna P.A., Kirkpatrick J.D., & McCarthy D.W. 1999, ApJ, 512, 864

Joshi M.M., Haberle R.M., & Reynolds R.T. 1997, Icarus, 129, 450

Kasting J.F., Whitmire D.P., & Reynolds R.T. 1993, Icarus, 101, 108

Kürster M., Hatzes A.P., Cochran W.D., Döbereiner S., Dennerl K., & Endl M. 1999, A&A, 344, L5


text:
# 1 Introduction

## 1 Introduction

Uncovering the origin, composition and energy spectrum of the highest energy cosmic rays is one of the biggest challenges in astroparticle physics. The Pierre Auger project is the next step in the search for answers to intriguing questions about the origin of these particles . The Auger Southern Observatory will consist of 1600 water Cerenkov detector stations (each 10 m2 x 1.2 m deep) on a hexagonal grid of 1500 m spacing overlooked by four detectors capable of detecting the fluorescence light emitted by the nitrogen molecules excited by the shower. The array covers a ground area of approximately 3000 km2 at a mean altitude of 875 g cm-2 ($`1400`$ m), near Malargüe in Mendoza State, Argentina (lat = -35.2, long = -69.2). The Auger Southern Observatory will be able to measure the energy of the incoming cosmic ray using fluorescence detectors (FD) and thus calibrate the energy inferred from the surface detectors.

The proposal to use the Auger observatories to search for very inclined showers induced by ultrahigh energy neutrinos  has led to an investigation of the characteristics of high energy showers at large zenith angles, i.e. showers arriving at zenith angles larger than 600 . These showers would not be very different from vertical showers except for the fact that they develop in the upper part of the atmosphere. As a result the electromagnetic part of the shower, produced mainly from $`\pi ^0`$ decay, is mostly absorbed well before the shower front reaches ground level. The muon front propagates through the atmosphere and it gets attenuated more slowly through pair production, bremsstrahlung, and hadronic interactions that reduce the muon energy and increase their probability to decay in flight. Therefore the muon energy spectrum at ground will have a low energy cutoff which increases as the zenith angle rises. As a result the average muon energy at ground level also increases. Although the bulk of the overall increase in the average muon energy with zenith is due to the rise of the cutoff, there is a smaller contribution due to the rise of the pion critical energy (the energy at which the pions are more likely to decay than to interact). These energetic muons travel making a small angle to the incoming cosmic ray direction but their trajectories are deflected by the Earth’s magnetic field. As a result the muon density patterns at ground are different from typical densities measured in vertical showers, reflecting the structure of the shower core and the geomagnetic field which acts as a kind of “natural magnetic spectrometer”.

This work outlines the sensitivity of the Auger Observatory to inclined air showers induced by ordinary cosmic rays, that is cosmic rays produced by protons, heavy nuclei, and photons. They constitute the background to neutrino detection but their observation also provides a significant increase in the aperture of the array and may improve mass composition studies . Below, we predict the expected rate for the Auger Surface Array for zenith angles above $`60^{}`$ assuming different primary mass compositions: proton, iron and photons. Primary photons can interact with the geomagnetic field before reaching the top of the atmosphere . This will have a major impact on the expected rate of detection for those primaries, as we will show. We can take advantage of this effect to outline a way to identify photon primaries at energies as high as 1020 eV.

The article is organized as follows: in section 2 we describe the procedure to calculate the expected rate for the Auger array above 60 for primary protons, giving the expected values for the energy resolution, core error reconstruction, and multiplicity of stations as a function of the energy: a detailed explanation of the procedure can also be found in ; in section 3 we describe the effect of the geomagnetic interactions in inclined air showers initiated by photon primaries; in section 4 we present the expected rate for three primary mass compositions (protons, iron and photons). Finally in section 5 we end up with some conclusions.

## 2 A simulation of inclined events and their analysis

The effect of the magnetic field on inclined showers makes the density patterns at ground level dependent on zenith and azimuth angles, in addition to their dependence on assumptions that must be made about energy, mass composition and interaction model as is usual in the the case for vertical showers. The work presented here relies on the model of the muon density patterns produced by inclined showers under the influence of the Earth’s magnetic field described in . This model approximates the magnetic distortions of circularly symmetric muon density patterns generated by simulation in the absence of magnetic effects. Moreover the model has revealed relatively simple energy scaling functions for the density patterns once the zenith and azimuthal angles are determined . The model provides a simple and continuous average muon density function which simplifies fitting procedures and allows us to simulate the expected sensitivity of the Auger Surface Array for showers with zenith angle above 60. The magnetic field at Malargüe has been taken to have an intensity $`B=26\mu `$T, an inclination of $`\beta _i=33.18^{}`$ (pointing in the upward direction) and a declination $`\beta _d=6.88^{}`$ (eastward). A northern site with a larger intensity and higher inclination would have a different distortion pattern and results are expected to be slightly different particularly in relation to azimuthal angle distributions.

We have simulated a subarray containing 91 detectors arranged in five concentric hexagons and showers are assumed to fall in the central region to save computing time. Results obtained in this subarray are then scaled up to the total array surface of 3000 km2. This approximate approach has two important implications. The subarray “radius” ranges between $`r_{sub}6.5`$ and 7.5 km. This implies that in the transverse plane the maximum distance to shower axis is $`r_{sub}\mathrm{cos}\theta `$ which can be rather small for high zenith angles (for example for a zenith of 75 the subarray stops sampling at about $`1.8`$ km). As a result the number of tanks that have a signal in the total array will be severely underestimated for high zenith and high energy showers. Another implication is that no edge effects are taken into consideration but this is not expected to affect much neither the shower rate nor most of our results because a large fraction of the showers will fall well within the array.

With the help of AIRES  (version 2.2.1) and this model, we have generated libraries with average muon density maps for inclined showers between 60 and 88 in steps of 2, and for all azimuth angles in steps of 5, at a constant energy of 1019 eV. Muon density profiles for different primary energies are obtained using the appropriate scaling factor which depends on primary composition and hadronic model considered . In this work we consider the QGSJET98  as the high energy hadronic interaction model. To account for the detector response to muons of different energy and impact angle we have used the WTANK  program based on the GEANT CERN package . As was previously shown in  for Haverah Park detectors, the signal of very inclined muons in water Čerenkov tanks can be greatly enhanced by direct light and muon interactions, which are included in these simulations. Between 60 and 70 the electromagnetic component, from $`\pi ^0`$ decay, was directly calculated using Monte Carlo simulation with AIRES and treated as a small correction to the tank signals, which is ignored for $`\theta >70^{}`$.

Using these muon density libraries, and assuming an isotropic cosmic ray distribution consisting of protons, and a given parameterization of the measured spectrum, we have simulated artificial events in the energy range 1018.5-1021 eV and in the zenith angle range 60-89. We have assumed a simple trigger condition demanding 5 stations with a density greater than 0.4 vertical equivalent muons (VEM) m-2. For each artificial event the densities at each detector are fluctuated according to the Poisson fluctuations in the number of muons, and also according to the distributions of light fluctuations in the tank. These zenith angle dependent distributions include the response to the geometry of the muon tracks, the direct light contribution and the catastrophic muon energy losses, namely hard bremsstrahlung, pair production or hadronic interactions which become more important as the average muon energy rises, as obtained in detailed simulations using WTANK.

Following the steps of  the shower parameters of these artificial events ($`\theta `$, $`\varphi `$, impact point position, and energy) are reconstructed as if they were real data. The procedure is based on a two-step fitting process using a maximum likelihood method which is iterated. Silent stations are taken into account in the fitting procedure. The arrival direction angles are first fitted to the signal times in the “triggered” detectors and then the energy and impact point are varied to make the average muon density map fit the actual “measured” densities in the tank. Curvature corrections to the arrival times are implemented once there is an estimate of the core location: the process is iterated three times for convergence.

An example of a reconstructed event in the plane perpendicular to the shower direction is shown in Fig. 1, together with the contours of densities that best fit the data. In the figure the array is rotated in the shower plane such that the $`y`$-axis is aligned with the component of the magnetic field perpendicular to the shower axis. The asymmetry in the density pattern due to the geomagnetic field is apparent. As anticipated , in spite of a relatively moderate energy, common events can be rather spectacular having over twenty tanks triggered.

We have estimated the detection rate as a function of primary energy integrating over all zenith and azimuth angles, using the parameterization of the energy spectrum of . This parameterization lies between the energy spectra recently reported by the AGASA group and the HiRes group ,. Assuming these parameterizations of the energy spectrum our predicted rate will increase (decrease) up to 20 %. For self-consistency, the energy spectrum to be used in a practical application of the ideas outlined here should be the one derived from the Auger experiment using the fluorescence technique, as the spectrum derived by this technique has the advantage of being independent of assumptions of primary mass composition.

The results are shown in Fig. 2. The dashed line is the integral rate that could be expected in a 3000 km2 array from the cosmic ray spectrum that has been assumed in this work. The continuous line represents the reconstructed spectrum. We have made quality cuts to the reconstructed events, namely the downward error in the energy was required to be less than 60% and the $`\chi ^2`$ probability for the energy and time fits have been forced to be greater than 1% for the accepted events. The difference in normalization, which amounts to $``$ $`20\%`$, is due to the combined effect of the quality cuts, systematic effects and the effect of the finite energy resolution of the detector.

The number events expected per year, assuming a pure proton composition, is $`1000`$ above $`10^{19}`$ eV and $`18`$ above $`10^{20}`$ eV after making these quality cuts. It should be noted that these numbers are obtained using the parameterization of the energy spectrum of , where no Greisen-Zatsepin-Kuzmin (GZK) cutoff is indicated.

The multiplicity of detectors with a signal above 0.4 VEM m-2 as a function of energy for the range of angles simulated is also shown in Fig. 2. For a given energy the multiplicity is much larger than for vertical showers. Two factors contribute to this effect: the reduction of the effective spacing between stations in the shower plane, and the flatter density profile of horizontal air showers compared to vertical.

As a shower becomes inclined the earliest and latest tanks that are above threshold can be quite far apart corresponding to different depths in shower development. Due to the additional atmosphere traversed an attenuation can be expected for the late arriving particles. This was studied for muon distributions which is the dominant component in inclined showers in  and for electron distributions in . If the number density of muons changes significantly over such a distance it can represent a significant correction to the densities derived in our work. It should be remarked that no attenuation of the shower particles across the array, has been taken into account in this work. As it is shown in , such attenuation effects become important at extreme zenith angles $`>87^{}`$, but for most of the zenith angle interval discussed in this article the attenuation, which is less important than in the vertical case, can be ignored in a first approximation. As an example, simulations show that for a 30 shower, water Cerenkov density at 1000 m from the shower core in the shower plane is attenuated by a $``$ 30 % when reaching ground level. At 700 this factor is less than 15 %. The main reason for this is that in the vertical direction the signal is dominated by electromagnetic particles from $`\pi ^0`$ decay, while in the horizontal direction the signal at ground level is dominated by high energy muons and secondary electromagnetic particles from decaying muons. The scale of the attenuation is roughly determined by the interaction length, which for photons and electrons is orders of magnitude smaller than that for muons.

We should also point out that in spite of using average density maps, statistical fluctuations of the shower densities are considered in this work. These are complicated effects that include fluctuations in the ratios of energy transported by charged and neutral pions in the hadronic interactions as well as fluctuations in depth of maximum. As obtained in simulations it typically represents a $`20\%`$ effect on the overall muon number density normalization for proton primaries and it is rather independent of primary energy . Much of the fluctuation corresponds to the first interaction. In the left panel of Fig. 3 we show the energy resolution for two energy ranges. Two factors contribute to the energy resolution calculated in this work: the variation of the number of muons at ground due to shower to shower fluctuations ($`20\%`$ for all energies) and the muon size reconstruction error (which evolves from 20 % at 1019 eV to 10 % at 1020 eV). Hence, at 1020 eV the energy resolution is dominated by shower to shower fluctuations. The errors in the reconstruction of the impact points of the analysed events are also shown in Fig. 3. As an example it is $`70`$ m at 60 and it has a small variation with zenith angle due to the different effective spacing of the detectors when projected on to the transverse shower plane.

## 3 Inclined showers induced by primary photons

The recent analysis of inclined data obtained at Haverah Park has opened a new window to establishing bounds on photon composition at ultrahigh energies . As photon showers have far fewer muons than hadron showers, the expected rate of inclined showers created by photons detected with a surface array is reduced relative to the rate from protons or nuclei primaries. This is because the detection of horizontal showers induced by photons and hadrons by surface detectors is practically only sensitive to the muon component. The determination of the spectrum of cosmic rays through vertical shower measurements using the fluorescence technique, which is fairly independent of the primary composition together with the inclined shower rate thus allow measurements of composition .

Following a similar procedure to that described in the previous section we obtain in this section the predicted integral spectrum for photon primaries for the Auger Observatories. There are a number of potential problems that could compromise the capabilities of the Auger Observatories to establish bounds or measure photon composition through inclined event rates, namely the Landau-Pomerančuck-Migdal  effect, the rise in the photoproduction cross section  and the interactions of photons with the geomagnetic field of the Earth . We discuss these effects in detail to show to what extent it can be expected that the Auger detector will be sensitive to photon primaries.

### 3.1 General features of inclined photon showers

The dominant process for muon production in an electromagnetic cascade is photoproduction:

$$\gamma +nucleushadrons$$
(1)
When photoproduction occurs the reaction products are essentially like those of a pion-nucleus interaction. Muons originate from the decay of the produced pions and kaons and their progeny in the resulting hadronic subshowers. The cross section for photoproduction has been measured up to $`10^4`$ GeV  for the incident photons in the laboratory frame. Above the resonance region (see Fig. 4) the cross section is about 100 $`\mu `$barns per nucleon, and rises slowly for photon energies above $`10`$ GeV. The corresponding cross section on an air nucleus is $``$ 1.1 mb, obtained by scaling $`\sigma _{\gamma A}A^\alpha \sigma _{\gamma p}`$ where $`\alpha 0.9`$. The ratio of the pair production to the photohadronic cross section gives the probability of a hadronic interaction to occur. This ratio is at 10 GeV:

$$R=\frac{\sigma _{\gamma hadrons}}{\sigma _{\gamma e+e}}2.810^3$$
(2)
The ratio grows with the incident photon energy, because of the rise in the photoproduction cross section, and it is expected to be $`10^2`$ at $`10^{19}`$ eV. Fig. 4 shows the measured photoproduction cross sections at different energies, together with the standard parameterization used by the AIRES code  and a recent parameterization given by Block et al.  from a combined analysis of $`pp`$, $`p\gamma `$ and $`\gamma \gamma `$ interactions.

As pointed out in ref., the photoproduction cross section is close to the saturation limit at the HERA energies and cannot grow with energy faster than the currently used parameterizations. Therefore, the Block et al. parameterization may be considered as an upper limit to the photoproduction cross section. Other potential sources of muons are smaller than this. Muon pair production, in which a pair of electrons is replaced by a pair of muons is suppressed by a factor $`m_e^2/m_\mu ^2210^5`$ . Hadron production by electrons also contributes less than photoproduction because the process has to occur through exchange of a virtual photon, and the energies of the produced hadrons are small.

The probability of photoproduction versus pair production is further enhanced at high energies when the Landau-Pomeranchuk-Migdal (LPM) effect takes place . The LPM effect is a collective effect of the electric potential of several atoms and it tends to suppress the pair-production and bremsstrahlung cross sections for energies above a given value $`E_{LPM}`$, which is inversely proportional to the medium density. It starts to become relevant for the development of photon induced air showers at primary energies above $`10^{19}`$ eV for vertical showers , a value somewhat above $`E_{LPM}`$. This energy rises for large zenith angles corresponding to shower development higher up in the atmosphere where the air is thinner. We have found from simulations using AIRES code that for the zenith angle range being discussed, the LPM effect has a mild impact on the total number of muons produced by photon showers ($`<`$ 25 % at a zenith angle of 60).

It has been shown before that to a very good approximation the number of muons produced by a proton and iron primary rises with energy according to a simple scaling law :

$$N_\mu =N_0(E/10^{19})^\beta ,$$
(3)
where $`N_0`$ (the normalization, here taken at $`10^{19}`$ eV) and $`\beta `$ (a constant parameter of order 0.9) are fixed numbers for a given zenith angle, hadronic interaction model and mass composition. The approximation is also valid for photon primaries although a number of differences are worth a brief discussion.

We have simulated batches of 100 photon induced air showers using AIRES code at different zenith angles and energies, in the absence of magnetic field, using the two different extrapolations of the photoproduction cross sections which are shown in Fig. 4, namely the standard parameterization of AIRES and the parameterization of ref . In Fig. 5 we plot the average number of muons generated by a photon as a function of energy for different zenith angles as obtained in the simulation. We note that the absolute number of muons of a photon shower is less than that of a hadron shower as expected from the ratio of probabilities in Eq. 2. However the number of muons in a photon-induced shower rises with energy faster than in a proton shower so that differences between the two become smaller as the energy rises. This is also qualitatively expected from the behaviour of the photoproduction cross section with energy. We also show for comparison an approximate scaling law using $`\beta =1.2`$ which we have used before as a conservative estimate .

We can observe in Fig. 5 a steepening kink of the curve at the highest energies. This is due to the LPM effect which makes the showers develop deeper and further enhances the ratio $`R`$. The rise in $`R`$ is a potential problem for the method. If the average number of muons in a photon shower becomes sufficiently close to that of a proton shower, no conclusion could be extracted regarding photon composition since a photon shower that photoproduces in the first interaction is practically undistinguishable from a proton shower of the same energy.

On the other hand, the difference in the number of muons produced for the two parameterizations used for the photoproduction cross section is small which gives an idea of the uncertainty involved in this calculation (see Fig. 5). This is not surprising since it is clear that the number of muons in a photon shower is dominated by processes at lower energies closer to the region where data constrain the cross section.

Fluctuations in the number of muons in photon-induced showers are rather different from ordinary cosmic ray showers. If the first interaction of the incident photon happens to be hadronic (probability R$``$ 0.01 at $`10^{19}`$ eV) then the shower is indistinguishable from a hadronic shower. We can thus expect a distribution for the total number of muons in a photon shower with a long tail in the region of showers with large number of muons, close to those of a proton shower. These non-gaussian tails can affect the detection rates and for that reason we have implemented fluctuations in the number of muons in photon shower using the distributions obtained in the simulations.

### 3.2 Interactions with the geomagnetic field

It was pointed out in  that photons above $`10^{19}`$ eV have a large probability to convert into an electron-positron pair in the presence of the Earth magnetic field before entering the atmosphere. The strength of this interaction depends on $`EB_{}`$ where $`E`$ is the photon energy and $`B_{}`$ is the transverse magnetic field. Therefore the probability of interaction will depend on the direction of the incoming photon with respect to the Earth frame. For a given zenith angle it is a maximum when the photon trajectory intercepts the Earth’s magnetic axis, which roughly corresponds to the geomagnetic south at the Malargüe site. This is because at a given distance to the center of the Earth the dipole field is strongest along the magnetic axis. For the same reason it is clear that the field will be largest for higher zenith angles. This azimuthal asymmetry was also pointed out by Stanev and Vankov  and by Bertou et al.  to be a potential way to establish a photon contribution in the highest energy cosmic rays.

The point we want to make is that this interaction reduces the potential problems that could arise because of the increase in the photoproduction cross section and therefore “safeguards” the method of using the inclined shower rate for composition studies also for the highest energies. If such an interaction is deemed to happen an electron-positron pair is created that will radiate strongly by synchrotron radiation. A large number of photons is thus produced and some of them may give secondary pairs provided they still have sufficient energy. Instead of a single particle entering the upper atmosphere we end up with a spectrum of gammas and one or a few electron positron pairs adding up to the primary photon energy.

To study this effect we have devised a Monte Carlo to simulate the cascading of photons in the geomagnetic field. The photon is injected at 20000 km from the top of the atmosphere along the incoming direction and tracked in steps of 100 km. The values of the interaction length for pair production are taken from . Fig. 6 shows the probability of conversion for photons at different energies, zenith and azimuth angles. When a photon is converted, the resulting electron and positron are tracked in steps of 2 km, and the number and energy of the photons emitted via synchrotron radiation is calculated using the spectral distribution in . The photons produced by magnetic bremsstrahlung are tracked and the probability of conversion calculated in each step. In this way, we obtain the spectrum of electrons, positrons and gammas at the top of the atmosphere.

The relevant parameter for the inclined shower rate is the total number of muons in the shower. When the photon converts in the magnetic field the number of muons at ground will be given by a sum over the contribution of the corresponding particle distribution:

$$N_\mu =\underset{i}{}N_0(E_i/10^{19})^\beta ,$$
(4)
where $`E_i`$ is the energy of particle $`i`$. Meanwhile for an unconverted photon it is given by Eq. 3. Since the slope parameter in this equation, $`\beta `$, is greater than one the number of muons is reduced when the photon interacts with the Earth’s magnetic field relative to when it does not.

It is of some interest to comment on a few facts about this conversion procedure. The strength of the Earth’s magnetic field is such that efficient conversion starts at energies of order $`10^{20}`$eV. When the photon energy is well above this energy, the photon is always converted, giving rise to a sort of universal photon spectrum. The shape of the resulting spectrum is governed by magnetic bremsstrahlung and it is independent of the primary energy and only mildly dependent on the arrival direction; only the total number of particles changes as can be seen in Fig. 7 for two different primary energies. The spectrum displays a low energy tail and a maximum at around $`10^{19}`$eV and can be parameterized by:

$$\frac{dN_{em}}{d\mathrm{log}E}=k\left[\frac{E}{E_c}\right]^\gamma \left[1+\frac{E}{E_c}\right]^\delta ,$$
(5)
where $`\gamma `$=0.24, $`\delta `$=-1.7, $`E_c`$=1019.76 eV, and $`k`$ is a normalization constant which due to energy conservation scales linearly with the primary photon energy. At very high energies there is a sharp cutoff in the spectrum due to the fact that all photons of such energies convert in the magnetic field. Thus at very high energies, when the photons convert, we can equivalently talk about an effective slope parameter of (Eq. 3) $`\beta _{\mathrm{eff}}=1`$.

At this point we can use the method described in the previous section to generate artificial events in the Auger array. We use the conversion probability for the photon showers to predict their behavior. These events are then reconstructed as if they were protons, and an equivalent energy ($`E_p`$) is calculated.

It is interesting to note that, as a coincidence, the energy at which photon conversion becomes efficient is also that at which the LPM effects start to show up in shower development, $`10^{19}`$ eV  for a vertical air shower. When a photon is converted a quick look at the spectrum shown in Fig. 7 shows that most of the particles that will reach the top of the atmosphere and initiate the shower will be below the LPM threshold. The magnetic field of the Earth shields the atmosphere and the probability that showers develop in the atmosphere with strong LPM suppression becomes very small. The effective shielding is strongest at high zenith angles because the energy at which the LPM effect starts to take up is higher.

This is fortunate because, as Fig. 5 shows, if a high energy photon reaches the top of the atmosphere, the number of muons produced in the corresponding shower could become rather large and close to that of proton induced showers.

## 4 Composition analysis with HAS

In Fig. 8 we show the integral energy spectra obtained from the artificial events under three different assumptions for the primary composition (protons, iron, and photons) using the cosmic ray parameterization given in . In all the cases the artificial events are generated with muon density maps obtained for each specific primary, and then reconstructed as if they were protons, obtaining an equivalent energy ($`E_p`$). For primary photons we have calculated the spectra switching on (off) the interaction with the geomagnetic field. We show the prediction assuming the Block et al. photoproduction cross section and a parameterization used in . This parameterization is based on AIRES simulations but the authors set the value of $`\beta `$ to be 1.2 to be conservative.

Three important features are clear in this picture:

* The expected rate for primary photons above 1020 eV is reduced 30 % when geomagnetic interactions are taken into account.
* The two photoproduction cross sections used predict a rate for photon primaries $``$ 10 times smaller compared to the predictions for other hadronic primaries.
* Geomagnetic interactions modify the slope of the spectrum for primary photons at energies above 1019.6 eV.

According to current wisdom and partly because of the geomagnetic effect the uncertainty on the photoproduction cross sections at high energies has a small effect on the muon number at ground level. The predicted distortion of the longitudinal shower development for photon primaries because of the LPM effect is reduced because of the geomagnetic interactions. Assuming that the vertical flux is well measured independently of mass composition, as expected to be achieved with a fluorescence detector, the rate of inclined showers could be used to impose strong constraints on primary photon composition assuming a given hadronic model.

The hadronic model uncertainty can change the normalization of the spectra shown in Fig. 8. Therefore this method is not adequate to impose constraints on both mass composition and hadronic models at the same time. It should be used in combination with vertical measurements on mass composition using the same hadronic model. However, the differences in muon densities between the different hadronic models available are $`20\%`$, which is comparable to the difference between densities produced by proton and iron primaries, and therefore, these primaries cannot be resolved by using the overall rate measurements. Nevertheless the difference between photon and hadronic primaries are much larger than 20 % so that this method can, on its own, impose severe constraints on the photon content of the highest energy cosmic rays.

Currently the ultra high energy data suggest a mostly hadronic composition , and the inclined showers rate is well described by proton primaries at energies above $`10^{19}`$ eV. If there is a change of composition to a photon dominance at higher energies, the rate of inclined showers will be much lower. After a few years of operation of the Auger Observatory these features will help to establish bounds on the flux of photons at energies as high as 1020 eV.

## 5 Conclusions

Inclined showers will be seen in the Auger Observatory as spectacular events with as many as 30 or 40 hit detectors. We have calculated the approximate rate of inclined showers with zenith angle exceeding $`60^{}`$ expected to be observed at the Pierre Auger Observatory. This rate increases the aperture of the observatory by almost a factor 2. Assuming a pure proton composition, there will be over 1000 well reconstructed events above $`10^{19}`$ eV, with a mean error energy $`25`$ %. The rate is sensitive to composition. If photons were dominant at high energy, the rate would be an order of magnitude smaller than if they were protons or nuclei allowing for a clear discrimination of the two cases. Uncertainties in the physics at very high energies have implications for our results on the detailed quantitative predictions but these uncertainties have little impact on the previous conclusions.

There are other signatures of the presence of photons (see for instance ref. ). Surely the combination of the inclined shower rate measurements together with vertical flux determinations and detailed analysis of the expected photon signatures will be a great step in the establishment of the overall photon rate at very high energies with the forthcoming data from the Auger experiment.

## Acknowledgments

We thank F. Halzen, M. Teshima, and K. Shinozaki for discussions and P. Billoir and S. Sciutto for very useful comments after careful reading of the manuscript. We also thank F. Halzen for bringing up the issue of photoproduction uncertainty. This work was partly supported by a joint grant from the British Council and the Spanish Ministry of Education (HB1997-0175), by Xunta de Galicia (PGIDT00PXI20615PR), by CICYT (AEN99-0589-C02-02 and) by MCYT (FPA 2001-3837) and by PPARC(PPA/9/5/1998/00453). The work of R.A.V. is supported by the ”Ramón y Cajal” program. We thank the ’Centro de Supercomputación de Galicia” (CESGA) for computer resources.


text:
# References

Title: Mechanisms of noise-resistance in genetic oscillators
Authors: José M. G. Vilar1,2, Hao Yuan Kueh1, Naama Barkai3, and Stanislas Leibler1,2
1Howard Hughes Medical Institute, Departments of Molecular Biology and Physics, Princeton University, Princeton, NJ 08544
2The Rockefeller University, 1230 York Avenue, New York, NY 10021
3Departments of Molecular Genetics and Physics of Complex Systems, Weizmann Institute of Science, Rehovot, Israel

Abstract

A wide range of organisms use circadian clocks to keep internal sense of daily time and regulate their behavior accordingly. Most of these clocks use intracellular genetic networks based on positive and negative regulatory elements. The integration of these “circuits” at the cellular level imposes strong constraints on their functioning and design. Here we study a recently proposed model \[N. Barkai and S. Leibler, *Nature*, 403:267–268, 2000\] that incorporates just the essential elements found experimentally. We show that this type of oscillator is driven mainly by two elements: the concentration of a repressor protein and the dynamics of an activator protein forming an inactive complex with the repressor. Thus the clock does not need to rely on mRNA dynamics to oscillate, which makes it especially resistant to fluctuations. Oscillations can be present even when the time average of the number of mRNA molecules goes below one. Under some conditions, this oscillator is not only resistant to but paradoxically also enhanced by the intrinsic biochemical noise.

The environment changes in a highly periodic manner. There are, among others, daily cycles of light and dark as well as annual cycles of changing climates and physical conditions. Such environmental periodicity may create the necessity for organisms to develop internal time-keeping mechanisms to accurately anticipate these external changes and modify their state accordingly . In particular, a wide range of organisms, as diverse as cyanobacteria and mammals, have evolved circadian rhythms —biological clocks with a period of about twenty four hours that evoke and regulate physiological and biochemical changes to best suit different times of the day.

Recent findings show that the molecular mechanisms upon which these clocks rely share many common features among species . The main characteristic is the presence of intracellular transcription regulation networks with a set of clock elements that give rise to stable oscillations in gene expression. A positive element activates genes coupled to the circadian clock. It simultaneously promotes the expression of a negative element, which in turn represses the positive element. The cycle completes itself upon degradation of the negative element and re-expression of the positive element.

A crucial feature of circadian clocks is the ability to maintain a constant period over a wide range of internal and external fluctuations . Such robustness ensures that the clock runs accurately and triggers the expression of clock-dependent genes at the appropriate time of the day. For instance, fluctuations in temperature affect chemical reaction rates and may perturb oscillatory behavior. Another source of fluctuations may be the presence of internal noise due to the stochastic nature of chemical reactions . Low numbers of molecules may be responsible for random fluctuations that can destabilize the oscillatory behavior of the biochemical network . Yet, circadian clocks maintain a fairly constant period amidst such fluctuations.

Description of the model

To study possible strategies, or principles, that biological systems may use to minimize the effect of stochastic noise on circadian clocks, we examine a minimal model based on the common positive and negative control elements found experimentally . This model is described in Figure 1. It involves two genes, an activator $`A`$ and a repressor $`R`$, which are transcribed into mRNA and subsequently translated into protein. The activator $`A`$ binds to the $`A`$ and $`R`$ promoters, which increases their transcription rate. Thus, $`A`$ acts as the positive element in transcription, whereas $`R`$ acts as the negative element by sequestering the activator.

The deterministic dynamics of the model is given by the set of reaction rate equations

$$\begin{array}{ccc}dD_A/dt\hfill & =& \theta _AD_A^{}\gamma _AD_AA\hfill \\ dD_R/dt\hfill & =& \theta _RD_R^{}\gamma _RD_RA\hfill \\ dD_A^{}/dt\hfill & =& \gamma _AD_AA\theta _AD_A^{}\hfill \\ dD_R^{}/dt\hfill & =& \gamma _RD_RA\theta _RD_R^{}\hfill \\ dM_A/dt\hfill & =& \alpha _A^{}D_A^{}+\alpha _AD_A\delta _{M_A}M_A\hfill \\ dA/dt\hfill & =& \beta _AM_A+\theta _AD_A^{}+\theta _RD_R^{}\hfill \\ & & A(\gamma _AD_A+\gamma _RD_R+\gamma _CR+\delta _A)\hfill \\ dM_R/dt\hfill & =& \alpha _R^{}D_R^{}+\alpha _RD_R\delta _{M_R}M_R\hfill \\ dR/dt\hfill & =& \beta _RM_R\gamma _CAR+\delta _AC\delta _RR\hfill \\ dC/dt\hfill & =& \gamma _CAR\delta _AC,\hfill \end{array}$$
(1)
where the variables and constants are as described in the caption for Figure 1. This simple model is not intended to reproduce the particular details of each organism but to grasp the properties that the core design confers. As in any general model, the parameters of the values we use are typical ones. For instance, the rates for bimolecular reactions are all in the range of diffusion limited reactions.

The preceding equations would be strictly valid in a well-stirred macroscopic reactor. At the cellular level, a more realistic approach has to consider the intrinsic stochasticity of chemical reactions . This can be done by transforming the reaction rates into probability transition rates and concentrations into number of molecules. One then obtains the so-called master equation which gives the time evolution of the probability of having a given number of molecules. There is no general procedure to solve this type of equation analytically, but it is the starting point to simulate the stochastic behavior of the system. The basic idea behind such simulations is to perform a random walk through the possible states of the system, which are defined by the numbers of molecules of the different reacting species. Starting from a state with given numbers of molecules, the probability of jumping to other state with different numbers of molecules (i.e. the probability for an elementary reaction to happen) can be computed from the master equation. One can pick up a state and the jumping time according to that probability distribution, then consider the resulting state as a new initial state, and repeat this procedure until some final state or time is reached. In this way, the numbers of molecules change in time with the statistical properties given by the master equation. There are several algorithms to implement this. The main difference among them is the specific way in which they compute the probabilities and select the states. For chemical reactions with few components, it is customary to use the so-called Gillespie algorithm .

This intrinsic probabilistic behavior in the evolution of the number of molecules gives rise to fluctuations that are usually referred to as noise. In general, the term noise is used for any disturbance interfering with a signal or with the operation of system. In the case of chemical reactions, the signal would be the average production of the reacting species whereas the disturbance would arise as a consequence of the fluctuations around that average value. We use term noise rather than fluctuations to emphasize the disturbing effect that these fluctuations can have. Thus, although related, both terms do not mean the same. For instance, there can be large fluctuations in some molecular species but, if their characteristic time is very short compared to those of other processes that take place, they would introduce little noise.

In Figure 2 we compare the results of the stochastic and deterministic approaches. We show the levels of $`A`$ protein and $`R`$ protein over time for the set of parameter values and initial conditions given in the caption of Figure 1. The deterministic results were obtained from numerical integration of Eqs. (1), whereas the stochastic results were obtained by computer simulation using the Gillespie algorithm. The main difference between the deterministic and stochastic time courses is the presence of random fluctuations in the latter. In the deterministic model every circadian cycle is identical to the previous one. The stochastic model shows some variability in the numbers of molecules and the period length, corresponding to the intrinsic fluctuations of the biochemical network. For these values of the parameters, both stochastic and deterministic approaches give similar results. We have also used different initial conditions and in all the cases we have observed that the behavior of the long term solution is the same.

Model simplification

To gain further insight into the essential elements that are responsible for the oscillations, we will simplify as much as possible the deterministic rate equations. By making various quasi-steady state assumptions , it is possible to simplify the set of Eqs. (1) into a two variable system with the repressor $`R`$ and the complex $`C`$ as the two slow variables:

$`{\displaystyle \frac{dR}{dt}}`$ $`=`$ $`{\displaystyle \frac{\beta _R}{\delta _{M_R}}}{\displaystyle \frac{\alpha _R\theta _R+\alpha _R^{}\gamma _R\stackrel{~}{A}(R)}{\theta _R+\gamma _R\stackrel{~}{A}(R)}}\gamma _C\stackrel{~}{A}(R)R+\delta _AC\delta _RR`$
$`{\displaystyle \frac{dC}{dt}}`$ $`=`$ $`\gamma _C\stackrel{~}{A}(R)R\delta _AC`$ (2)
where

$$\stackrel{~}{A}(R)=\frac{1}{2}(\alpha _A^{}\rho (R)K_d)+\frac{1}{2}\sqrt{(\alpha _A^{}\rho (R)K_d)^2+4\alpha _A\rho (R)K_d}$$

with $`\rho (R)=\beta _A/\delta _{M_A}(\gamma _CR+\delta _A)`$ and $`K_d=\theta _A/\gamma _A`$. Notice that the nonlinearity in the equations enters trough the quasi-equilibrium value of $`A`$, $`\stackrel{~}{A}(R)`$, which is a function of $`R`$. The main idea behind these approximations is that there are fast and slow variables. Fast variables are assumed to be at an effective equilibrium whereas slow variables are responsible for the dynamics of the system. Thus, given the set of Eqs. (1), we assume that all the derivatives except $`dR/dt`$ and $`dC/dt`$ are zero.

In Figure 3 we show that for the values of the parameters we use the numerical solutions for the trajectories of the two-variable system \[Eqs. (2)\] agree closely with the solutions of the full system \[Eqs. (1)\], except for quantitative differences in the peak levels and times at the beginning of each cycle. These differences arise because the time scale separation between fast and slow variables is not sufficiently large for quasi-steady state assumptions to be exact. These results indicate, nevertheless, that the dynamics of the system is mainly determined by two component concentrations: those of the complex and the repressor. The other components are driven mainly by these two elements and their effects enter the system through effective parameters. It is worth emphasizing that the reduced two-variable model is aimed just to offer insights into the qualitative behavior of the system and to show how the properties that one observes in the full system are already present in a simple two variable model. Thus, whenever we present simulation results for the deterministic system, except if otherwise stated, we are referring to the full system. Regarding the validity of the two-variable model, it is a good approximation when the dynamics of mRNA and the activator is faster than that of the complex and repressor. For instance, it will remain valid if $`\delta _R`$ is decreased or some of $`\delta _A`$, $`\delta _{M_A}`$, $`\delta _{M_R}`$, $`\theta _A`$, and $`\theta _R`$ are increased with respect to the parameters of the caption for Figure 1.

Limit cycle oscillations and stability analysis

The existence of oscillations in the two-variable system can be inferred from application of the Poincaré-Bendixson theorem . This theorem states that a two-dimensional system of the type we are considering exhibits limit cycles if it is confined in a closed bounded region that does not contain any stable fixed points. The trajectory of the system is confined since the number of molecules cannot reach infinite values. The fixed points and their stability can be determined by following a standard linear stability analysis. There is a single fixed point for positive concentrations. In our case, the signs of the real parts of the eigenvalues of the matrix describing the linearized dynamics around this point are given by

$$\tau =\left|\frac{\stackrel{~}{A}(R)}{R}\right|\left[\gamma _CR\frac{\beta _R}{\delta _{M_R}}\frac{(\alpha _R^{}\alpha _R)\theta _R\gamma _R}{[\theta _R+\gamma _R\stackrel{~}{A}(R)]^2}\right][\gamma _C\stackrel{~}{A}(R)+\delta _R+\delta _A],$$
(3)
where $`\tau `$ is the trace of the matrix. All the quantities are evaluated at the fixed point $`(R_0,C_0)`$. When $`\tau `$ is positive, the real part of the eigenvalues is positive, the fixed point is unstable, and there is a limit cycle in the system, which gives sufficient conditions for the existence of oscillations. Evaluation of Eq. (3) shows that $`\tau `$ is indeed positive for the set of parameters we are using. The domain in which $`\tau `$ is positive is rather broad. For instance, the function $`\tau `$ remains positive when $`\theta _A`$ and $`\theta _R`$ are multiplied by a factor $`K`$ with $`0.024<K<10.7`$; when all transcription rates ($`\alpha `$ and $`\alpha ^{}`$) are multiplied by $`K`$ with $`K>0.08`$; or when protein and mRNA degradation rates are multiplied by $`K`$ with $`0.0009<K<3.5`$. When $`\tau `$ is negative, the fixed point is stable and, in principle, no conclusion about the existence of limit cycles can be drawn. In this case, the presence of oscillations could also depend on the initial conditions. For the full model, the ranges of parameters that give rise to oscillations are not exactly the same but remain very close to the previous ones.

The mechanism responsible for oscillations is illustrated in Figure 4 through the phase portrait of the two-variable model. Starting with low numbers of initial molecules near the origin of the phase plane, the trajectory rapidly shoots upwards along the $`\dot{R}=0`$ nullcline (the dot over a variable means its time derivative). Here, the high levels of $`A`$, present due to auto-activated transcription, rapidly induce the formation of the complex $`C`$. Reaching the maximum of the nullcline, the trajectory ‘falls off’ the edge and moves rapidly diagonally right and downwards, corresponding to a drop in $`C`$ and a rise in $`R`$. The trajectory curves around the $`\dot{R}=0`$ nullcline and hits the $`\dot{C}=0`$ nullcline, where it slowly returns to the left and approaches the fixed point $`(R_0,C_0)`$. When approaching the fixed point, $`\dot{R}`$ decreases sharply, taking the trajectory past the fixed point and sending it back upwards to initiate a new cycle.

The trajectory in Figure 4 comprises a fast phase corresponding to the rapid production of $`C`$ and $`R`$, and a slow phase corresponding to the slow degradation of $`R`$. These two distinct phases are characteristic of excitable systems, the classic example of which is the Fitz Hugh-Nagumo model for action potential transmission in neurons . The fast and slow leg correspond to the excitable and refractory phase of the system, respectively. Thus, the system oscillates as it avoids the fixed point $`(R_0,C_0)`$ and hits the $`\dot{R}=0`$ nullcline on the left to begin the excitable phase of a new cycle.

As we have already pointed out, the deterministic analysis can be useful to grasp the main properties of the system under certain conditions. Unfortunately such conditions are not known a priori without a stochastic analysis. Surprisingly enough, we have found that parameter values that give rise to a stable steady state in the deterministic limit continue to produce reliable oscillations in the stochastic case, as shown in Figure 5. Therefore, the presence of noise not only changes the behavior of the system by adding more disorder, but can also lead to marked qualitative differences.

How can the system continue to produce oscillations even when deterministic rate equations predict a stable steady state? The system always evolves towards a stable fixed point, as sketched in Figure 6. However, a perturbation of sufficient magnitude near the fixed point, e.g. as illustrated by the dotted arrow in Figure 6 (notice that the figure has not been drawn to scale and that the size of the arrow is not representative of the actual size of the perturbation), may initiate a new cycle. For low numbers of molecules, the intrinsic fluctuations of chemical reactions can be large enough to continually send the system into the fast phase after each cycle and thus produce sustained oscillations. In the deterministic limit (or close to it), there are no perturbations (or the perturbations are too small) to initiate a new cycle and the trajectory stays close to the fixed point. In this case, the system performs better if enough noise is present in the system. This situation is analogous to that observed in the Fitz Hugh–Nagumo model, where an optimal amount of noise maximizes the reliability of the oscillations . It is important to realize that the effects that noise may have on non-linear systems can be difficult to predict and rather paradoxical . Therefore, the smaller number of molecules does not necessarily imply more irregular behavior of the system, as one might intuitively assume .

Significant parameters and noise resistance

The mechanism responsible for oscillations involves only two variables. This has some important consequences for the functioning of the clock. If we consider the deterministic limit, a two-dimensional dynamical system of this type either oscillates regularly or does not oscillate at all. In two dimensions, since trajectories cannot cross, fixed points and periodic orbits are the only possible attractors. Other more complicated behaviors such as chaos or quasi periodicity are not allowed . On the other hand, the intrinsic stochastic fluctuations of the remaining variables are effectively averaged and do not significantly affect the performance of the system.

For instance, one variable that usually plays a prominent role in many circadian rhythm models is the number of mRNA molecules . In our case, however, mRNA does not enter directly into the dynamics. It is just an intermediate step in the production of the proteins. Thus, if protein production remains unaltered the system will oscillate regardless of the number of mRNA molecules involved. This can be achieved, for instance, by increasing simultaneously the degradation rates of mRNA and the translation rate of the proteins. In the deterministic limit of the two-variable model this has no effect on protein dynamics at all. In the stochastic case the effects are negligible. Figure 7 shows the time evolution of repressor mRNA and protein levels in the system for $`\beta _A`$, $`\beta _R`$, $`\delta _{M_A}`$ and $`\delta _{M_R}`$ multiplied by a factor of hundred. The system essentially alternates between having zero and one mRNA molecule in the cycle, and the proteins continue to exhibit remarkably good oscillations.

There are also parameters that do affect the properties of the oscillations. In the deterministic limit, oscillations are always regular (provided that the two-variable model is a good approximation). When fluctuations are taken into account, the reliability of the oscillations may depend on those parameters. One such parameter, as we have seen in the previous section, is the repressor degradation rate $`\delta _R`$. This parameter affects the period of the oscillations (compare e.g. Figures 2 and 5) and also the stability of the fixed point. One can infer from Eq. (3) that for high or low values of $`\delta _R`$ the fixed point becomes stable. In such cases, the deterministic system stops oscillating but this does not need to be so for the stochastic one, which may continue producing reliable oscillations.

Notice that the positive feedback is a key element in the clock dynamics. Its most obvious use is the generation of the instability that gives rise to oscillations. But it has other not so obvious role that is closely related to the resistance to noise. In general, gene regulation is a slow process (with typical characteristic times of about one hour) and, as a such, is prone to be affected by fluctuations. This problem gets even worse if the dynamics relies in several coupled transcriptional feedbacks since the effects of the fluctuations are then amplified. The positive feedback gives a fast transcriptional switch, allowing to move fast from low to high transcription rates. In this way, the time in which the system is prone to fluctuations is greatly reduced.

Conclusions

The presence of noise in transcriptional and enzymatic networks is a fundamental consequence of the stochastic nature of biochemical reactions. The ability to function effectively and consistently amidst such random fluctuations is a major issue in gene expression and network behavior. In this paper we have studied how different factors affect a simple model for circadian rhythms that exhibits noise resistance. We found that the oscillations in this model are driven mainly by two components: a repressor protein and an activator-repressor complex. This fact is responsible for the reliability of the oscillations. First, a two-dimensional dynamical system of this kind has a very simple behavior: in the deterministic limit, it either oscillates or goes to a steady state. Second, noise and perturbations in the other variables affect the system only slightly since they do not enter directly the dynamics. Finally, resistance to noise is achieved as the number of molecules of any of the two key components reaches small values only for short periods of time or when they are not driving the dynamics of the system. In this way, even though some molecular species may be present in very low numbers, the intrinsic stochasticity of biochemical reactions can be bypassed.

It is important to emphasize that organisms have evolved networks to function in the extremely noisy cellular environment. Suitable network designs, as those that are now emerging from the experimental data , can confer resistance against such noise. In addition, some of these networks may not only be resistant to but could also be taking advantage of the cellular noise to perform their functions under conditions in which it would not be possible by deterministic means.

Acknowledgements

This work was supported by the U.S. National Institutes of Health.

Figure Captions

* Biochemical network of the circadian oscillator model. $`D_A^{}`$ and $`D_A`$ denote the number of activator genes with and without $`A`$ bound to its promoter respectively; similarly, $`D_R^{}`$ and $`D_R`$ refer to the repressor promoter; $`M_A`$ and $`M_R`$ denote mRNA of $`A`$ and $`R`$; $`A`$ and $`R`$ correspond to the activator and repressor proteins; and $`C`$ to the the inactivated complex formed by $`A`$ and $`R`$. The constants $`\alpha `$ and $`\alpha ^{}`$ denote the basal and activated rates of transcription, $`\beta `$ the rates of translation, $`\delta `$ the rates of spontaneous degradation, $`\gamma `$ the rates of binding of $`A`$ to other components, and $`\theta `$ the rates of unbinding of $`A`$ from those components. Except if otherwise stated, along this paper we have assumed the following values for the reaction rates: $`\alpha _A=50\text{ hr}^1`$, $`\alpha _A^{}=500\text{ hr}^1`$, $`\alpha _R=0.01\text{ hr}^1`$, $`\alpha _R^{}=50\text{ hr}^1`$, $`\beta _A=50\text{ hr}^1`$, $`\beta _R=5\text{ hr}^1`$, $`\delta _{MA}=10\text{ hr}^1`$, $`\delta _{MR}=0.5\text{ hr}^1`$, $`\delta _A=1\text{ hr}^1`$, $`\delta _R=0.2\text{ hr}^1`$, $`\gamma _A=1\text{ mol}^1\text{ hr}^1`$, $`\gamma _R=1\text{ mol}^1\text{ hr}^1`$, $`\gamma _C=2\text{ mol}^1\text{ hr}^1`$, $`\theta _A=50\text{ hr}^1`$, $`\theta _R=100\text{ hr}^1`$, where mol means number of molecules. The initial conditions are $`D_A=D_R=1\text{ mol}`$, $`D_A^{}=D_R^{}=M_A=M_R=A=R=C=0`$, which require that the cell has a single copy of the activator and repressor genes: $`D_A+D_A^{}=1\text{ mol}`$ and $`D_R+D_R^{}=1\text{ mol}`$. The cellular volume is assumed to be the unity so that concentrations and number of molecules are equivalent. Notice that we assume that the complex breaks into $`R`$ due to the degradation of $`A`$ and therefore the parameter $`\delta _A`$ appears twice in the model.
* Oscillations in repressor and activator protein numbers obtained from numerical simulations of the deterministic ($`a`$,$`b`$) and stochastic ($`c`$,$`d`$) descriptions of the model.
* Time evolution of the quantities $`R`$ (continuous line line) and $`C`$ (dashed line) for the system reduced to two variables ($`a`$) by various quasi-steady state assumptions and for the complete system ($`b`$).
* Phase portrait of the two variable oscillator \[Eqs. (2)\] for the parameter values given in the caption for Figure 1 (the drawing is not to scale). The thick line illustrates the trajectory of system. $`(R_0,C_0)`$ is the fixed point of the system, and $`\dot{R}dR/dt=0`$ and $`\dot{C}dC/dt=0`$ are the $`R`$ and $`C`$ nullclines respectively. The solid arrows give the orientation of the direction field on the nullclines.
* Time evolution of $`R`$ for the deterministic \[Eqs. (1)\] ($`a`$) and stochastic ($`b`$) versions of the model. The values of the parameters are as in the caption of Figure 1 except that now we set $`\delta _R=0.05\text{ hr}^1`$. For these parameter values $`\tau <0`$ so that the fixed point is stable.
* Phase portrait as in Figure 4 but for a situation in which the system falls into the stable fixed point $`(R_0,C_0)`$. The dotted arrow to the left of the fixed point illustrates a perturbation that would initiate a single sweep of the (former) oscillatory trajectory.
* Stochastic time evolution of the number of activator ($`a`$) and repressor ($`c`$) molecules; and activator ($`b`$) and repressor ($`d`$) mRNA molecules. The values of the parameters are as in the caption for Figure 1 but now with $`\beta _A=5000\text{ hr}^1`$, $`\beta _R=500\text{ hr}^1`$, $`\delta _{MA}=1000\text{ hr}^1`$, and $`\delta _{MR}=50\text{ hr}^1`$.

FIGURE 1

FIGURE 2

FIGURE 3

FIGURE 4

FIGURE 5

FIGURE 6

FIGURE 7


text:
# 1 Introduction

## 1 Introduction

The microquasar XTE J1550-564 was discovered with the All Sky Monitor on board RossiXTE  in September 1998  and to date it showed four outbursts. The compared analysis of all these outbursts allows to relate different behaviours of the source with different values of flux reached during the outburst. The outburst here analyzed is quite different from the first two and similar to the 2001 one. Radio observations were reported by . A more complete analysis is described in .

## 2 Energy spectra

We analyzed spectra in a range from 3 to 150 keV from PCA  and HEXTE  instruments corresponding to channels 4-52 for the PCA and 10-50 for HEXTE. We used a model consisting of a power law with a high-energy cutoff, a smeared Fe edge and a Fe emission line. The hydrogen column density was fixed at 8.5 $`\times `$ 1021 cm-2 and a systematic error of 0.75% was quadratically added to the data. For the first seven observations the line was not necessary while for the remaining the edge was not required. The absence of a soft thermal component may be due to the lack of data at energies lower than 3 keV. The unabsorbed 3-150 keV flux decreased from 1.2 $`\times `$ 10-8 erg cm-2 s-1 to 7.1 $`\times `$ 10-10 erg cm-2 s-1, corresponding to a luminosity at 6 kpc of 5.0-0.2 $`\times `$ 1037 erg s-1.

The values found for the power law photon index are tipical of a low/hard state. The cutoff energy decreases substantially during the observations (from $``$340 keV to $``$64 keV) while the photon index does not increase much (from $``$1.4 to $``$1.5).

## 3 Power density spectra

The variability observed in the PCA light curve is tipical of the low/hard state. For each observation we created a Power Density Spectrum (PDS) in the following way: we produced individual Leahy-normalized  PDS from 512-second stretches with time resolution of 512 s-1 from PCA data ranging from channel 0 to 35 (roughly corresponding to energies from 2 to 15 keV). These PDS were averaged and the contribution due to Poisson statistic was subtracted . The averaged PDS were normalized to squared fractional rms .

Each PDS was fitted with a model consisting of a sum of Lorentzian components : we always needed only two broad Lorentzian to reach good values of the chi-squared. In Figure 1 all the PDS extracted are shown, together with their best fit Lorentzian components.

## 4 Discussion

From the analysis of energy spectra we can see how the presence of a high energy cutoff in the power law component is consistent with the presence of a thermal distribution of electrons in the inner region. In such a model the decrease of the cutoff energy implies a decrease of electron temperature: since the photon index does not change correspondingly the optical depth of the Comptonizing cloud must have increased at the same time.

Interpreting the two Lorentzian components present in the PDS as Lb and L (see ), we can verify how the values found are consistent with existing correlations found in other sources, like the one concerning $`\nu _b`$ vs. $`\nu _{\mathrm{}}`$  and $`\nu _b`$ vs. rms2 . The timing analysis results confirm that XTE J1550-564 is in a low/hard state.

The general behaviour showed by XTE J1550-564 during this outburst is compatible with the one observed in other BHC in their low/hard state, indicating that at lower fluxes (and therefore possibly lower accretion rates) the source behaves like “normal” BHC.