text:
# First Passage Time Statistics For Systems Driven by Long Range Gaussian Noises

## I Introduction

The arrival of a random process at a particular state often triggers some important behavior. Among the large number of instances that one can name are the firing of a neuron, the nucleation of a phase associated with a phase transition, the triggering of an alarm, the occurrence of a major earthquake, and the crossing of an activation barrier by a reaction coordinate that converts reactants to products in a chemical reaction. These examples range from the macroscopic to the microscopic, indicating that the arrival problem is of interest on all scales. In many instances it is not only important to know the statistics of occurrence of such events, but more specifically the statistics of first occurrence. This then leads to the study of the statistical properties of the time that it takes a random process to reach a specified state for the first time, that is, the mean first passage time. The first passage time literature is enormous and extends over many decades stratonovich ; wax ; gardiner ; risken ; redner .

For simplicity we restrict our attention in this work to one-dimensional random processes $`x(t)`$, although the generalization to vector processes or even to fields is conceptually straightforward. The statistical properties of the temporal evolution of the random process are typically described in one of two ways. One is by way of an evolution equation for the probability distribution $`P(x,t)`$ for the random variable to take on the value $`x(t)=x`$ at time $`t`$ (the random variable and the values that it can take on are often denoted by the same symbol, $`x`$ in this case, although the distinction should be kept in mind). The other is by way of an evolution equation (“Langevin equation”) for the random variable itself. From the solution of the appropriate evolution equation (via steps that are in general non-trivial) one can then obtain the first passage time statistics for $`x(t)`$ to reach a prescribed value for the first time. In either case, one often thinks of the evolution of $`x(t)`$ as being driven by fluctuations $`\eta (t)`$ of prescribed statistical properties from which the statistics of $`x(t)`$ then follow.

Let us then turn to the “driving fluctuations” $`\eta (t)`$ and their typical statistical properties. By far the most common assumption is that fluctuations have a Gaussian distribution, although there has always been a great deal of interest in fluctuations that are distributed in some other fashion. The most common examples include Lévy distributions ditlevsen ; grigolini1  (distributions with very long tails), and dichotomous processes ourbook ; grigolini2  (where the random variable can take on only two values). The former typically arise when the fluctuations are a result of many multiplicative inputs; the latter often serve as an analytically tractable prototype.

In order to specify the fluctuations fully, one must also explicitly state their correlation properties. If they are Gaussian,

$$P(\eta )=\frac{1}{\sqrt{2\pi }\sigma }e^{\eta ^2/2\sigma ^2},$$
(1)
then only the average $`\eta (t)`$ (usually and here as well taken to be zero) and two-time correlation function (assumed stationary) $`\gamma (tt^{})\eta (t)\eta (t^{})`$ need to be specified, since all other correlation properties then follow. In Eq. (1) $`\sigma ^2=\gamma (0)`$. For other distributions it is in general necessary to also specify higher order correlations functions, even though these are often difficult to determine. The early literature dealt primarily with $`\delta `$\- correlated or “white noise” processes,

$$\gamma (t)=2D\delta (t)$$
(2)
and $`D`$ (the integral of the correlation function) measures the intensity of the noise. In this case the second moment $`\sigma ^2`$ diverges as $`D/\mathrm{\Delta }t`$, where $`\mathrm{\Delta }t0`$ is the “width” of the $`\delta `$-function and therefore the shortest time scale in the problem. The evolution equation $`P(x,t)`$ for a random process driven by Gaussian $`\delta `$-correlated fluctuations $`\eta (t)`$ is the familiar Fokker-Planck equation stratonovich ; risken , and the associated first passage time properties are well understood stratonovich ; risken ; FPfpt .

About two decades ago a great deal of attention began to be directed toward understanding the effects of “colored noise”, that is, of driving fluctuations that are not $`\delta `$-correlated. Most of the attention was focused on exponentially correlated noise,

$$\gamma (t)=\frac{D}{\tau }e^{|t|/\tau },$$
(3)
where $`\tau `$ is the correlation time of the noise. The associated evolution equation for $`P(x,t)`$ and related first passage time properties as a function of the correlation time were studied in detail and are also well understood reviews . For instance, it is firmly established that the mean first passage time from one potential well to another of a process $`x(t)`$ driven by exponentially correlated noise increases with increasing correlation time. It is also understood that the particular exponential form of the correlation function is not crucial in the qualitative features of the first passage time statistics. The most important feature determining the qualitative features of the first passage time statistics is the fact that there is a finite correlation time associated with the fluctuations.

More recently there has been considerable interest in (typically Gaussian) fluctuations that display long-range power-law correlations with an infinite correlation time. Such highly correlated fluctuations have been considered in a broad array of circumstances ranging from the biological to the physical, from the economical to the atmospheric, and encompassing theoretical and experimental studies makse96 ; Eva1 . There are two separate issues that need to be addressed with such highly correlated noise. One is the numerical generation of the correlated noise itself. Only recently have the traditional Fourier filtering methods been superseded by far more efficient procedures makse96 ; nises ; pang ; Barabasi ; Romero99  that will briefly be reviewed below in the context of our systems. The second issue is that of the effect of such a driving noise on the first passage time properties of the system that is driven by these highly correlated fluctuations, and this is our principal interest in this paper. Herein we examine the mean first passage time of a particle, driven by highly correlated fluctuations, in three different physical situations. In the first case, we discuss the arrival of a free particle (moving superdiffusively) at an absorbing boundary. As a second case, we study the decay of a particle from an initial unstable state. Finally, we discuss the time for a particle to cross a barrier separating one well from the other in a double well potential.

In Sec. II we summarize the main steps in the generation of long ranged correlated Gaussian noises. In Secs. III, IV, and V we present respectively the arrival of the superdiffusive particle at an absorbing boundary, the decay of an unstable state, and the barrier crossing problem. Finally, we conclude with a summary in Sec. VI.

## II Generation of long range correlated Gaussian noise

The main steps to generate noises with a Gaussian distribution and arbitrary correlation properties can be found in Refs. makse96 ; nises ; pang ; Barabasi ; Romero99 . Here we present a summary of these “spectral methods” including some details relevant to our applications. The goal is to generate a Gaussian noise $`\eta (t)`$, with correlation function $`\gamma (t)`$ defined by:

$$\eta (t)\eta (t^{})=\gamma (tt^{}),$$
(4)
and with a Fourier transform

$$\stackrel{~}{\gamma }(\omega )=_{\mathrm{}}^{\mathrm{}}𝑑te^{i\omega t}\gamma (t).$$
(5)
This correlation function may be specified analytically or numerically. In the $`\omega `$ Fourier-space, the transformed noise $`\stackrel{~}{\eta }(\omega )`$ has a correlation function

$$\stackrel{~}{\eta }(\omega )\stackrel{~}{\eta }(\omega ^{})=2\pi \gamma (\omega )\delta (\omega +\omega ^{}).$$
(6)
The algorithm for the noise generation can be summarized as follows. First, the time interval $`(0,t)`$ is discretized into $`N=2^n`$ intervals of mesh size $`\mathrm{\Delta }t`$. This time interval has to be much smaller than any other characteristic time of the system (and hence $`N`$ must be sufficiently large) because it is used as the time integration step in the numerical simulation of the stochastic evolution equations. The intervals in time will be denoted by a roman index and the resulting frequency intervals in Fourier space are denoted by a greek index.

In discrete Fourier space the noise has a correlation function given by

$$\stackrel{~}{\eta }(\omega _\mu )\stackrel{~}{\eta }(\omega _\mu ^{})=N\mathrm{\Delta }t\stackrel{~}{\gamma }(\omega _\mu )\delta _{\mu +\mu ^{},0},$$
(7)
where $`\stackrel{~}{\eta }(\omega _\mu )`$ can be constructed as follows:

$$\begin{array}{ccc}\stackrel{~}{\eta }(\omega _\mu )=\sqrt{N\mathrm{\Delta }t\stackrel{~}{\gamma }(\omega _{\mu )})}\alpha _\mu ,\hfill & & \mu =0,\mathrm{},N1\hfill \\ & & \\ \stackrel{~}{\eta }(\omega _N)=\stackrel{~}{\eta }(\omega _0),\hfill & & \omega _\mu =\frac{2\pi \mu }{N\mathrm{\Delta }t}.\hfill \end{array}$$

Here the $`\alpha _\mu `$ are Gaussian random numbers with zero mean and with correlation

$$\alpha _\mu \alpha _\nu =\delta _{\mu +\nu ,0}.$$
(8)
This type of $`\delta `$-anticorrelated noise can be generated rather easily if the symmetry properties of real periodic series in the Fourier space ($`\alpha _i`$) are used nises ; Romero99 . Owing to the periodicity of $`\stackrel{~}{\gamma }(\omega _\mu )`$, the index $`\mu `$ can be made to run from $`N/2`$ to $`N/2`$, and the maximum frequency is $`\omega _{max}=\pi /\mathrm{\Delta }t`$ Romero99 . The Fourier components are related by $`\alpha _{\mu +pN}=\alpha _\mu `$ for any integer $`p`$, and $`\alpha _\mu =\alpha _\mu ^{}`$ because the original random numbers in time are real. The anticorrelated random numbers can then be constructed as $`\alpha _\mu =a_\mu +ib_\mu `$ with $`b_0=0`$, $`a_0^2=1`$, and $`a_\mu `$ and $`b_\mu `$ are each Gaussian random numbers with zero mean and a variance of $`1/2`$ for $`\mu 0`$.

The discrete inverse transform of any sequence $`\stackrel{~}{\eta }(\omega _\mu )`$ is then numerically calculated by a standard Fast Fourier Transform algorithm. The result is a string of $`N`$ numbers, $`\eta (t_i)`$ which, by construction, have the proposed time correlation (4). However, due to the symmetries of the Fourier transform only $`N/2`$ of these values are actually independent and the remaining numbers are periodically correlated to them. Thus we have constructed a true random process from $`t=0`$ to a maximum time $`T_{cutoff}=\mathrm{\Delta }tN/2`$.

This algorithm is sufficiently general to allow generation of noise of any given correlation function. For instance, in Romero99  we considered Gaussian noise with a Gaussian correlation function,

$$\gamma (t)=\frac{2D}{\tau \sqrt{2\pi }}e^{t^2/2\tau ^2}$$
(9)
whose Fourier transform can be obtained analytically,

$$\stackrel{~}{\gamma }(\omega )=2De^{\tau ^2\omega ^2/4}.$$
(10)
According to our prescription, if we generate a discrete field of random numbers according to

$$\stackrel{~}{\eta }(\omega _\mu )=\left[2DN\mathrm{\Delta }t\mathrm{exp}\left(\frac{\tau ^2}{(\mathrm{\Delta }t)^2}(\mathrm{cos}(2\pi \mu /N)1)\right)\right]^{1/2}\alpha _\mu $$
(11)
for sufficiently large $`N`$ the resulting correlation function upon Fourier inversion should match Eq. (9) (it does). Similarly, for exponentially correlated Gaussian noise, cf. Eq.(3), whose Fourier transform is

$$\stackrel{~}{\gamma }(\omega )=\frac{2D}{1+\tau ^2\omega ^2},$$
(12)
we generate the discrete random numbers according to

$$\stackrel{~}{\eta }(\omega _\mu )=\left(\frac{2DN\mathrm{\Delta }t}{1+\left(\frac{2\tau }{\mathrm{\Delta }t}\mathrm{sin}(\pi \mu /N)\right)^2}\right)^{1/2}\alpha _\mu $$
(13)
and again reproduce the correlation function (3) upon Fourier inversion. In both cases $`2D`$ is the intensity of the noise and $`\tau `$ its correlation time. We find that it does not matter whether we discretize the $`\omega `$ using the function $`\mathrm{cos}(2\pi \mu /N)`$ as in (11) or $`\mathrm{sin}(\pi \mu /N)`$ as in (13).

A more difficult noise to generate numerically is one with a memory of inverse power law form,

$$\gamma (t)=\frac{\epsilon }{(1+|t|/t_0)^\beta }.$$
(14)
Here $`t_0`$ is an adjustable (small) parameter. For reasons described below and related to our particular algorithm, the value of $`t_0`$ turns out to depend on $`\beta `$. When the power law decay exponent $`\beta >2`$ then not only is the intensity of the fluctuations finite (and given by $`2\epsilon t_0/(\beta 1)`$, which increases with decreasing $`\beta `$), but they also have a finite correlation time (which also grows as $`\beta `$ decreases. The effects of this sort of noise are essentially the same as those of any Gaussian noise with a finite correlation time, although direct comparison with most of the literature would have to be done carefully reviews  because correlation time effects have usually been studied in the context of noise of a fixed intensity. When $`1<\beta <2`$ then the fluctuations have a finite intensity but an infinite correlation time. We have not analyzed this case. Even more persistent correlations, and the ones of interest to us here, occur when $`0<\beta <1`$.

In order to implement the spectral method, we would first have to know the Fourier transform of the correlation function. The form (14) does not have an closed-form Fourier transform note , and so instead we postulate or guess the form

$$\stackrel{~}{\gamma }(\omega _\mu )=\frac{\epsilon \beta \pi \omega _{max}^\beta }{\left[\frac{2}{\mathrm{\Delta }t}\mathrm{sin}(\pi \mu /N)+\omega _0\right]^{(1\beta )}}$$
(15)
where $`\omega _{max}`$ was defined earlier as $`\pi /\mathrm{\Delta }t`$ and where $`\omega _0`$ is a low frequency cutoff that is chosen to control the low-frequency behavior and avoids a zero-frequency divergence. To check whether Eq. (15) is indeed an appropriate choice one has to perform the numerical inverse transform and compare with (14). For a given exponent $`\beta `$ this comparison involves the three parameters $`\mathrm{\Delta }t`$ (the discretization time step), $`\omega _0`$, and $`t_0`$. The value of $`\omega _0`$ that leads to the best agreement turns out to depend on the time step. In our calculations we have mostly used $`\mathrm{\Delta }t=0.01`$ but in our barrier-crossing calculations in Sec. V we take $`\mathrm{\Delta }t=0.02`$. For this range of $`\mathrm{\Delta }t`$ we find that the best choice (giving the best concordance between the desired correlation function and the assumed Fourier transform) is $`\omega _0=0.0001`$. This is the value used throughout this analysis.

From Eq.( 15) we obtain an important relation,

$$\gamma (0)\frac{1}{\pi }_0^{\omega _{max}}\stackrel{~}{\gamma }(\omega )𝑑\omega =ϵ+𝒪(\omega _0/\omega _{max}),$$
(16)
consistent with Eq. (14). Therefore $`\gamma (0)`$ is independent of $`\beta `$ and is a direct measure of the noise control parameter $`\epsilon `$. This provides a first check for our algorithm.

A second test of our algorithm is presented in Figs. 1 and 2, where we exhibit the numerically generated correlation functions as well as the analytic counterparts with $`\beta =0.3`$ and $`0.6`$. Within the resolution of these figures, the simulated results agree well with the expected power laws given by Eq.(14) for the particular parameter choices shown in the caption. The value of $`\epsilon `$ has been assumed to be selectable independently of $`\beta `$ and according to Eq. (16). Figure 3 shows that this is indeed the case, i.e., that $`\gamma (0)`$ as generated from our simulations is essentially independent of $`\beta `$. The small variations are due to the discreteness of the time variable and the stochasticity of the problem. More complicated is the choice of the parameter $`t_0`$: for a fixed $`\epsilon `$, the choice of $`t_0`$ that leads to a best match between the numerically generated correlation function and the analytic form (14) depends on $`\beta `$, a dependence that is exhibited in Fig. 4. The parameter $`t_0`$ also varies with the time step $`\mathrm{\Delta }t`$, a dependence we have not exhibited because we hold $`\mathrm{\Delta }t`$ fixed in our analysis. In any case, implementation of our algorithm requires attention to these dependences. The values of $`t_0`$ indicated in Figs. 1 and 1 result from this best fit.

## III Superdiffusion of a free particle

The first problem we study is of interest in the analysis of the mean exit time of a free Brownian particle from a domain bounded by absorbing walls. We consider a free particle that moves in one dimension under the influence of a long range correlated noise until it covers, on average, the distance to one of the absorbing barriers.

A dynamical equation that describes this type of dynamics is the very simple stochastic differential equation (Langevin equation)

$$\frac{dx}{dt}=\eta (t),$$
(17)
where $`\eta (t)`$ is a Gaussian power law-correlated random noise as described in section II. Note that here and in all subsequent examples we “start the clock” at $`t=0`$, that is, in addition to specifying an initial value or distribution of $`x`$ at $`t=0`$ we also assume that the noise is “turned on” at that time. Otherwise one would need to be concerned about the role of the time $`t=0`$ in the evolution of the noise, which in turn would affect the description of its correlation function.

We are particularly interested in the dependence of the mean arrival (absorption) time on the parameters $`\beta `$ and $`ϵ`$ and on the length $`2L`$ of the interval. The essential scaling features (exponents) of the dependences on these parameters can be obtained from a formal integration of Eq.(17), which immediately leads to the mean squared displacement relation

$$x^2(t)\epsilon t^{2\beta },$$
(18)
characteristic of superdiffusive behavior Romero99 . Superdiffusion here arises because the arrival at the boundaries is more likely to occur via essentially ballistic motion than in the white noise case, where frequent reversals of $`dx/dt`$ cause the net displacement to be slow. The time at which the mean squared displacement is proportional to $`L^2`$ can be obtained by inverting this expression. The difference between this time and the mean time $`T`$ for first arrival at an absorbing boundary is due to recrossings of the boundary and only affects the prefactors, not the exponents. That is, we can use Eq. (18) to write

$$T\epsilon ^{1/(2\beta )}L^{2/(2\beta )}.$$
(19)
Note that the (Gaussian) statistics of the noise play no role in this result, which is determined entirely by the correlation function of the noise.

In order to test this theoretical prediction we have performed numerical simulations for different values of the interval length $`2L`$, the exponent $`\beta `$, and the coefficient $`\epsilon `$. Throughout this paper all Langevin equations are numerically integrated using the Heun method, which is an extension of a second order Runge-Kutta algorithm for stochastic differential equations Toral95 . The results presented below are averages over $`2000`$ realizations.

A comparison of the theoretical exponent ($`2/(2\beta )`$) in Eq. (19) and the simulation results is summarized in Table I. The agreement is clearly excellent. In Fig. 5 the average exit time is presented as function of the interval length on a log-log plot. The exit time $`T`$ follows the predicted power law behavior with the slopes of Table I. In the figure we also show the result for white noise with its characteristic exponent of $`2`$.

Fig. 6 shows the dependence of the mean first passage time on the noise intensity parameter $`\epsilon `$. The simulation results are in excellent agreement with the predictions of Eq. (19).

## IV Decay of an unstable state

The decay of unstable states is of course easily triggered by any disturbance, and in many physical instances the decay is caused by fluctuations. In its simplest rendition this process can be modeled by a particle initially placed at the top ($`x=0`$) of an inverted parabolic potential and subject to fluctuations,

$$\frac{dx}{dt}=xx^3+\eta (t),x(0)=0.$$
(20)
The decay at the early stages is dominated by the linear term and the fluctuations, so that in this regime the Langevin equation can be further simplified to sancho89

$$\frac{dx}{dt}=x+\eta (t).$$
(21)
Formal integration of this linear equation immediately yields

$$x(t)=h(t)e^t,h(t)_0^t𝑑t^{}e^t^{}\eta (t^{}).$$
(22)
The mean first passage time to, say, $`x=\pm 1`$ is identified as the time at which $`x^2(t)=1`$. Recrossings are even less relevant here than in the previous section; larger values of $`x`$ require consideration of the nonlinear term in the evolution equation. If the time of interest is $`t1`$ (a condition that will be seen to be satisfied by the mean first passage time) then the upper limit in the integral in Eq. (22) can be extended to infinity, that is, we set

$$x^2(t)=h^2e^2t,h=_0^{\mathrm{}}𝑑t^{}e^t^{}\eta (t^{}).$$
(23)
Thus the fluctuations $`h^2`$ simply play the role of a random initial condition on $`x^2(t)`$. The random variable $`h`$ has zero mean and variance $`\sigma ^2`$ given by caceres ,

$`\sigma ^2=h^2={\displaystyle _0^{\mathrm{}}}𝑑t^{}e^t^{}\eta (t^{}){\displaystyle _0^{\mathrm{}}}𝑑t^{\prime \prime }e^{t^{\prime \prime }}\eta (t^{\prime \prime })=`$
$`{\displaystyle _0^{\mathrm{}}}e^{2t^{}}𝑑t^{}{\displaystyle _0^{\mathrm{}}}𝑑t^{\prime \prime }e^{t^{\prime \prime }}\gamma (s)={\displaystyle \frac{ϵ}{2}}t_0^\beta e^{t_0}\mathrm{\Gamma }(1\beta ,t_0).`$ (24)
The stationarity of the noise has been used to calculate the variance; $`\mathrm{\Gamma }(1\beta ,t_0)`$ is the incomplete Gamma function.

The mean first passage time to $`x^2=1`$ for this process for a generic colored noise was calculated in Ref. sancho89 :

$$T=\frac{1}{2}\mathrm{ln}\left(\sigma ^2\right)+\frac{\gamma }{2},$$
(25)
where $`\gamma =0.57721\mathrm{}`$ is the Euler constant. This result together with (24) gives an explicit expression for the mean first passage time. Figures 7 and 8 show excellent agreement between the theory and numerical simulations for the dependence of $`T`$ on $`\epsilon `$ and on $`\beta `$.

## V The Barrier Crossing Problem

A problem that is more complex than the essentially linear ones posed so far is the “barrier crossing problem.” Here we consider the diffusion of an overdamped particle in a double-well potential. The mean first passage time of interest is the time it takes the particle to go from one of the potential minima to the other when the transition is driven by power-law-correlated Gaussian fluctuations. The particle dynamics is modeled by the following Langevin equation:

$$\frac{dx}{dt}=\frac{dV(x)}{dt}+\eta (t),$$
(26)
where the double-well potential is,

$$V(x)=\frac{1}{2}x^2+\frac{1}{4}x^4.$$
(27)
Our simulations follow the dynamics of each particle starting in the left well at $`x(t=0)=1.0`$ until it arrives at the right well at $`x(t=T)=1.0`$ for the first time. We then calculate the mean first passage time $`T`$, that is, the average of $`T`$ over many ($`10,000`$) realizations.

Two questions are of interest: (1) If two noises $`\eta (t)`$, one $`\delta `$-correlated and the other power-law correlated, lead to the same average transition time from one well to the another, are there other properties that allow a clear distinction between them? (2) For a power-law correlated noise $`\eta (t)`$, how does the transition time depend on the correlation function parameters? We address both of these questions below.

The numerical barrier crossing problem is fundamentally different from the problems considered in the previous section. In those, no matter how small or persistent a realization of the noise $`\eta (t)`$, the process $`x(t)`$ will eventually reach one of the boundaries of interest. How long a simulation must run in order to make sure that all processes have reached the boundaries is essentially a matter of insuring that those with the lowest value of $`\eta `$ as given in the discretization scheme arrive there. In the barrier crossing problem, on the other hand, small-$`\eta `$ realizations will remain in one well until the value of $`\eta `$ changes to a sufficiently high value (easy to estimate reviews ) so as to effectively eliminate the barrier. The time it takes to effect this change is of course large (infinite on average) for the fluctuations considered herein. In our simulations we have carefully chosen run times that insure that all $`10,000`$ particles in our ensemble have made the passage from one well to the other. As we decrease $`\epsilon `$ (thus weakening the noise) and/or increase $`\beta `$ (thus in effect also weakening the noise), the simulation time has to be increased accordingly. Note that these remarks address the simulation problem. A complete theory of this process would have to deal with the relative magnitudes of a number of infinities (noise intensity, noise correlation, and passage over the barrier of the slowest processes).

Fig. 9 shows two typical trajectories of the barrier crossing process for two different Gaussian noises. One is white noise, the other is power-law correlated, and the parameters have been chosen so that the mean first passage times for both are essentially the same. The effect of the long range correlations is certainly not evident from these trajectories. This observation suggests a more detailed study of the dynamical properties. The first passage time distribution and, in particular, the mean first passage time for the Gaussian white noise case, are well known analytically (see reviews  and the many original references therein). Approximate theories for exponentially correlated noise also abound reviews , but to our knowledge there are no results for the infinite correlation time case. To deal with this problem we note that the trajectories (in either the white or the correlated noise cases) consist of two distinct components characterized by two different time scales. One component describes the rapid fluctuations in each well, and the other describes the much slower switching events between the two wells. We implement this observation by writing the trajectory $`x(t)`$ as a sum of these two contributions:

$$x(t)=x_D(t)+x_B(t),$$
(28)
where $`x_B(t)`$ represents the rapid random motion inside each well, and $`X_D(t)`$ is a dichotomous random process between the values $`x=1`$ and $`x=1`$, with a characteristic time controlled by the mean transition time $`T`$. The dynamics of $`x_B(t)`$ can be approximated by expanding Eq. (26) around the minimum of one of the wells, e.g. as $`x(t)=1+x_B(t)`$ in the well. Then $`x_B(t)`$ evolves according to the Langevin equation

$$\frac{dx_B}{dt}=2x_B+\eta (t),$$
(29)
which evolves on a characteristic time scale $`\tau =0.5`$ that is much smaller than $`T`$.

With this decomposition, the power spectrum of $`x(t)`$, denoted by $`S(\omega )`$, is just the sum of the spectra associated with $`x_D(t)`$ and $`x_B(t)`$. These can be calculated explicitly:

$$S(\omega )=\frac{4T}{4+\omega ^2T^2}+\frac{\tau ^2\gamma (\omega )}{1+\omega ^2\tau ^2}.$$
(30)
The first contribution is the spectrum of a dichotomous process governed by the time scale $`T`$, and the second contribution comes from the exact solution of Eq. (29) with the noise spectrum $`\gamma (\omega )`$. In the white noise case $`\gamma (\omega )=2D`$, and for the long ranged noise the spectrum is obtained from Eq. (15) in the continuum limit,

$$\gamma (\omega )=\frac{ϵ\beta \pi \omega _{max}^\beta }{(\omega +\omega _0)^{1\beta }}$$
(31)
It is interesting to note that the main difference between the white and power law spectra appear in the high-frequency (short-time) regime. In the white noise case, $`S(\omega )\omega ^2`$, whereas in the correlated case we observe the more rapid decay $`S(\omega )\omega ^{(3\beta )}`$.

Figure 10 shows the spectra for the two trajectories shown in Fig. 9 as well as the prediction (30) for each. As noted earlier, the parameters have been chosen so as to lead to the same mean first passage time $`T`$, but other than that there are no adjustable parameters. The agreement between our predictions and the numerical results is quite good in both cases. The exponent obtained from the numerical data in the high frequency regime for the correlated noise case is 2.49, which agrees with the theoretical prediction of 2.5. It is also clear that a distinction between the two types of fluctuations is difficult even at this level of detail for a given mean transition time.

In the above analysis $`T`$ is the adjustable parameter, that is, we have fixed $`\beta `$ and $`\epsilon `$ so as to yield a particular value of $`T`$. The dependence of this quantity on the model parameters $`\epsilon `$ and $`\beta `$ is yet to be explored.

The first panel in Fig. 11 is a plot of $`\mathrm{log}T`$ as a function of $`1/\epsilon `$. The straight line behavior shows that the transition time is described by a law of the Kramers form,

$$T\mathrm{exp}\left(\frac{A}{\epsilon }\right).$$
(32)
Fitting the simulation data in the figure with a straight line we obtain $`A=1.933`$, very far from the value $`A=2/27`$ obtained for exponentially correlated noise reviews . We stress that the “activation parameter” $`A`$ here depends on $`\beta `$, that is, a plot exactly like that of the first panel Fig. 11 but for a different value of $`\beta `$ leads to a different value of $`A`$. The values of $`A`$ vs $`\beta `$ collected in this way are shown in the second panel in Fig. 11. In contrast, for noises with a finite (albeit large) correlation time $`A`$ is essentially a constant independent of the noise intensity and of the correlation time. Comparison with those results is therefore not appropriate.

The dependence of the mean transition time on $`\beta `$ qualitatively parallels that of the activation parameter but is even more complex because it is non-monotonic. The typical behavior is shown in Fig. 12. Although we have no detailed theory to account for this non-monotonicity, a qualitative explanation is possible in terms of two competing effects, one that causes a decrease in $`T`$ with increasing $`\beta `$ and the other that causes $`T`$ to increase with increasing $`\beta `$. At very small $`\beta `$ the mean first passage time is dominated by those realizations of very slow passage for which $`\eta (0)`$ is not large enough to cause a transition and a very long wait is involved before $`\eta (t)`$ changes to a value that is sufficiently large. With increasing $`\beta `$ the wait for change decreases, thus leading to faster passage. On the other hand, increasing $`\beta `$ in effect decreases the intensity of the noise. To see this suppose that we calculate an effective intensity by integrating the correlation function only up to a maximum time $`t_{max}`$ ($`t_0`$), the running time of our simulations:

$`{\displaystyle _0^{t_{max}}}𝑑t\gamma (t)`$ $`=`$ $`{\displaystyle \frac{\epsilon t_0}{1\beta }}\left[\left(1+{\displaystyle \frac{t_{max}}{t_0}}\right)^{1\beta }1\right]`$ (33)
$`=`$ $`{\displaystyle \frac{\epsilon t_0}{1\beta }}\left({\displaystyle \frac{t_{max}}{t_0}}\right)^{1\beta }\left[1+𝒪\left({\displaystyle \frac{(1\beta )t_0}{t_{max}}}\right)\right].`$
The $`1\beta `$ in the exponent of the large ratio $`t_{max}/t_0`$ dominates the $`1\beta `$-dependence in the denominator, so that for a fixed $`t_{max}`$ this decreases rapidly with increasing $`\beta `$. This lower effective intensity of the noise leads to slower passage. We conjecture that the interplay of the two effects leads to the non-monotonic behavior observed in Fig. 12, although the specific value and specific location of the minimum depend on $`\epsilon `$.

In order to verify this minimum and to exhibit more information on the statistics of $`T`$, we have also obtained numerical results for the probability distribution $`P(T)`$ of the exit times (of which $`T`$ is the average). In Figs. 13 and 14 we have plotted the numerical histogram of $`P(T)`$ for white noise and for a power-law correlated noise, again chosen so as to correspond to the same means. The scatter of numerical points around the analytic forms discussed below (solid lines) is quite large and could be decreased by doing many more realizations for much longer times.

In the case of white noise the first passage time distribution is known to be exponential risken ; FPfpt ,

$$P(T)=\frac{1}{T}\mathrm{exp}\left(\frac{T}{<T>}\right).$$
(34)
Our numerical results agree with this result. For the power law case we conjecture a stretched exponential distribution for the transition time,

$$P(T)=\frac{\theta }{T_s\mathrm{\Gamma }(1/\theta )}\mathrm{exp}\left(\left[\frac{T}{T_s}\right]^\theta \right),$$
(35)
where the prefactor has been chosen to ensure normalization ($`\mathrm{\Gamma }(z)`$ is the Gamma function) and where the parameters $`\theta `$ and $`T_s`$ are obtained from the simulation data. The mean first passage time in terms of these parameters is

$$T=T_s\frac{\mathrm{\Gamma }(2/\theta )}{\mathrm{\Gamma }(1/\theta )}=\frac{T_s}{\sqrt{2\pi }}2^{\frac{2}{\theta }\frac{1}{2}}\mathrm{\Gamma }(\frac{1}{\theta }+\frac{1}{2}).$$
(36)
We do not display the $`\epsilon `$-dependence of the parameters because it is the $`\beta `$ dependence of the first passage time that requires further understanding. Our fitting procedure for the case of $`ϵ=1.15`$ and $`\mathrm{\Delta }t=0.02`$ gives the $`\beta `$-dependences observed in Fig. 15. Although the dependences of the parameters on $`\beta `$ are monotonic, the dependence of $`T`$ on $`\theta `$ is not monotonic (note that $`T_s`$ also varies with $`\theta `$) and gives rise to the observed minimum. Using the data shown in Fig. 15 to construct $`T`$ according to Eq. (36) one again obtains the minimum as a function of $`\beta `$ observed in Fig. 12. Another point worth noting is the approach of the exponent $`\theta `$ to unity with increasing $`\beta `$ that is seen in Fig. 15. The exponent is exactly unity for Gaussian white noise and would also be unity for $`\beta >2`$.

From the proposed distribution (35) one can evaluate the second moment of the first passage time distribution. We obtain

$$\frac{T^2}{T^2}=\frac{\mathrm{\Gamma }(1/\theta )\mathrm{\Gamma }(3/\theta )}{\mathrm{\Gamma }(2/\theta )^2},$$
(37)
For white Gaussian noise this ratio is equal to $`2`$. We see in Fig. 16 that for correlated noise the ratio is greater than $`2`$ reflecting the greater width of the stretched exponential distribution, but with increasing $`\beta `$ the approach toward $`2`$ is evident. The figure shows the moment ratio $`T^2/T^2`$ obtained from numerical simulations (full circles) and from the analytic expression Eq. (37).

## VI Summary and Conclusions

We have presented a numerical study and derived (in some cases) or conjectured (in others) analytic results for the first passage time problem in systems driven by long-range correlated noises. We specifically considered Gaussian noise, as opposed to noises with Lévy or other long-tailed distributions, or noises defined in terms of long waiting times between events, that may also lead to long-range correlations. There is an extensive literature on these non-Gaussian problems klafter , but very little on the Gaussian counterpart makse96 ; nises ; Romero99 . We have considered highly correlated noise with an inverse power law form for the correlation function, Eq. (14), with $`0<\beta <1`$. Not only does this noise have no finite correlation time, but the correlation function is not integrable. The parameters of the problem are the exponent $`\beta `$ and the noise control parameter $`\epsilon `$.

First we discussed a numerical algorithm to generate this type of noise. This is in itself not a trivial problem. Then we considered several systems driven by Gaussian power-law correlated noises with a special focus on the problem of first passage to a prescribed boundary or set of boundaries. Two of these problems, the arrival of a free particle in one dimension at either end of a finite interval, and the decay of an unstable state, admit straightforward analytic solution. The dependence of the respective mean first passage times on the noise parameters were calculated explicitly and compared favorably with numerical simulations.

Our final application is more complex: Here we considered the first passage of a particle that evolves in a double well potential from one well, over the barrier, to the other well. We specifically considered two issues. One is a comparison of the trajectories of such a particle subject to Gaussian white noise and to Gaussian power-law correlated noise, with parameters chosen in such a way that the mean first passage times for crossing from one well to the other in both cases is the same. It is difficult to find a clear difference that would reveal in an experiment which type of noise the particle was subjected to. Only the high frequency decay of the respective power spectra shows some differences. The other issue we explored is the dependence of the first passage time distribution on the parameters $`\beta `$ and $`\epsilon `$. We found that the dependence of the mean first passage time on the noise control parameter is of the usual “activated” form, but with an effective activation parameter that depends on $`\beta `$. We also found that the mean first passage time exhibits a minimum as a function of $`\beta `$, which we explain on the basis of competing mechanisms, one of which leads to an increase of the mean first passage time with increasing $`\beta `$ and the other to a decrease. This non-monotonicity is a signature characteristic of these fluctuations. We also found that a stretched exponential form for the first passage time distribution describes the numerical results, with a stretched exponent $`\theta `$ and a characteristic time $`T_s`$ that both depend on $`\beta `$. This stretched exponential form is consistent with the non-monotonic behavior of the mean first passage time, and also gives results in agreement with the simulation outcomes for the second moment of the first passage time distribution. The stretched exponential first passage time distribution, and the dependence of the parameters of this distribution on the noise parameters $`\epsilon `$ and $`\beta `$, are features that have no analog in Gaussian noises with a finite correlation time.

## Acknowledgements

The work was supported in part by the U.S. Department of Energy under Grant No. DE-FG03-86ER13606, by FONDECYT, Chile under Grant 1010988, by the Millenium Project under Grant P99-135F, and by the Comisión Interministerial de Ciencia y Tecnología (Spain) Project No. BFM2000-0624.


format: text

text:
# The Y-Band at 1.035 𝜇m: Photometric Calibration and the Dwarf Stellar / Sub-stellar Color Sequence

## 1 Introduction

The traditional near-infrared photometric bandpasses that are observable from the ground through windows of moderate atmospheric transparency are JHKLMNQ (and their slight variants K’, Ks, L’, M’, N1, N2, N3, etc). This set of filters defined broadly as in Johnson (1966), with the H-band being devised soon thereafter by Becklin and similarly the Q-band by Low, neglected to include one rather clean window through the atmosphere near 1.035 $`\mu `$m. This region of the electromagnetic spectrum is in between the traditional “optical” and “infrared” wavelength regimes, which are usually divided from one another based on detector technology (nowadays CCDs vs InSb or HgCdTe arrays) around 1 $`\mu `$m. We discuss in this paper a filter called “Y” that takes advantage of this window and has particular interest for studies of low-mass and extremely cool stars and substellar objects, i.e. brown dwarfs.

The Y filter bandpass was designed by S.E.P. and K.M. after visually inspecting the USAF’s LOWTRAN (Selby & McClatchey, 1975) atmospheric transmission curve and also referring to Manducca and Bell (1979). The half-power points of the filter were aimed at the 75% transmission points of the broad atmospheric absorption features near 0.955 and 1.112 $`\mu `$m. The filters were manufactured by Barr Associates (Worchester, MA). We name the filter “Y” in publication as a compromise between its current discrepant designations (X and Z as described below). Further, we choose Y to distinguish this filter from z and Z filters in use which are quantitatively different.

Nearly identical Y filters reside in the Keck Near-InfraRed Camera (NIRC, where it is called “Z”) and the Palomar 60” cassegrain InfraRed Camera (P60/IRC, where it is called “X”), as well as in the Palomar 200” D-78 and D-80 cameras. At 77 K, the Keck/NIRC filter has $`\lambda _c=1.032\mu `$m and $`\delta \lambda _c=0.157\mu `$m with $`>`$50% transmission from 0.965 - 1.115 $`\mu `$m while the P60/IRC filter has $`\lambda _c=1.04\mu `$m and $`\delta \lambda _c=0.15\mu `$m with $`>`$50% transmission from 0.9540-1.1105 $`\mu `$m. The two filter profiles are nearly identical in shape with the Keck/NIRC filter 3-5% less transmissive. It was found to be difficult to produce the specified filter with peak transmission much above 70%; transmission profiles are given in Table 1. In Figure 1, filter transmission curves for the P60/IRC filter set are shown relative to a model of the atmosphere above Mauna Kea, produced using the program IRTRANS4 as presented on the UKIRT worldwide web pages 111http://www.ukirt.jach.hawaii.edu/JACpublic/UKIRT/astronomy/calib/atmos-index.html. for an altitude of 4200 km assuming 1.2 mm precipitable water vapor and zenith pointing. The Y-band is a relatively clean atmospheric window, both in terms of molecular absorption (Figure 1) and in terms of O2 and OH night sky emission lines which are insignificant in the Y-band compared to J-band and H-band. Thermal background at Y-band is also insignificant.

The only other filter located in this atmospheric window is, ironically, the F1042M medium/narrowband filter in WFPC-2 aboard HST which has $`\lambda _c=1.015\mu `$m and $`\delta \lambda _c=0.03\mu `$m. Since the original preparation of this manuscript we have also become aware of the UFTI “Z” filter as described in Leggett et al. (2002). This filter has half-power points at 0.851 and 1.056 $`\mu `$m (their Table 2), slightly blueward of and broader than our Y filter, and inclusive of the strong and broad atmospheric water absorption around 0.93 $`\mu `$m (see Figure 1). The optical “Gunn z” and its successor the “Sloan z” filters (Fukugita et al. 1996) as well as the “Bessel z”, by contrast, are defined as long-pass filters and hence have their effective profiles determined by a rapidly falling CCD response curve towards 1 $`\mu `$m (Figure 1). The effective bandpass of these optical z filters is centered near 9250 $`\AA `$ but can vary substantially from instrument to instrument depending on the CCD (e.g. SITe detectors vs the new red-sensitive LBNL devices as compared in NOAO Newsletter 67, p. 3).

The Y filter introduced here has a well-defined bandpass, as opposed to the long-pass optical z filters, and takes advantage of the relatively clean atmospheric window centered near 1.035 $`\mu `$m, in moderate contrast to the UFTI Z filter. Y-band photometry therefore should be relatively independent of the specific detector and relatively insensitive to site-to-site differences in the transparency vs wavelength and the night-to-night variability of water vapor at a given site.

The motivating forces behind the present study are the instrumental/atmospheric advantages of the 1.035 $`\mu `$m window just mentioned, combined with the peak brightness of cool and ultra-cool dwarfs slightly longward of 1 $`\mu `$m. As shown in Figure 2, wavelengths $``$1$`\mu `$m represent the absolute peak in the spectral energy distributions of low-mass stars with spectral types $``$M5, and a strong local peak for stars and young brown dwarfs with spectral types $`>`$L6 (Burrows et al. 2001). In these very late type dwarfs the atmospheric opacity redistributes flux far from expectations based on blackbody radiation (Marley et al., 1996). Bolometric corrections become positive in the near-infrared bands, and J-H and H-K colors become bluer with decreasing effective temperature (Burrows et al., 1997).

Our goal here is to investigate the behavior of near-infrared colors involving the Y-band as a function of spectral type, and to determine the slope of the interstellar reddening vector. Our specific interest is in assessing the utility of the Y-filter in studies of low-mass objects in star-forming regions.

## 2 Observations

Imaging observations were obtained by L.A.H. on the nights of 1999, July 30 and August 1 (UT) using the cassegrain InfraRed Camera (P60/IRC; Murphy et al. 1995) at the Palomar 60” telescope under photometric sky conditions. The instrument contains a 256$`\times `$256 HgCdTe NICMOS-3 array with 40 $`\mu `$m pixels and produces a final platescale of 0.62”/pixel. Observations at Y (called “X” in P60/IRC), H, and K were obtained using a 9-point dither pattern centered on target with 10” separation between adjacent grid points.

Dwarf optical standard stars spanning a range of B-V color from -0.32 (spectral class O) to 2.2 (spectral type M6.5) were selected from Landolt (1992) and supplemented at the red end with late-type objects from Kirkpatrick, Henry, & McCarthy (1991), Tinney (1993), Kirkpatrick, Henry, & Simons (1995), Reid, Hawley, & Gizis (1995), Hawley, Gizis, & Reid (1996), Kirkpatrick et al. (1999, 2000), and Burgasser et al. (2002). The latest spectral type observed was L8V. Infrared standard stars used for calibration of the H- and K-bands were selected from Persson et al. (1998) and include both unreddened G-type stars and reddened objects of unknown spectral type. A set of 8 A0 stars was observed in order to define the instrumental offsets between the Y filter and the H and K filters. The star GSPC S813-D (Lasker et al., 1988) was used to map the atmospheric extinction curve at Y-band. Finally, a moderately reddened field ($`A_V`$ = 10-20 mag) in the Ophiuchus dark cloud was selected from the large area survey of Barsony et al. (1997) in order to study the interstellar extinction curve at Y-band. On 2001, September 9, some additional data for this program was obtained using the same equipment, by A. Burgasser and M. Brown, focussing on T dwarfs.

Also, observations that pre-date those above were made of 6 late-M- and L-type dwarfs using the Keck I telescope and NIRC (Matthews & Soifer 1994) by L.A.H. and J. Carpenter on 1999 February 10 (UT). Details are given in Hillenbrand & Carpenter (2000). Data of relevence to the present discussion were taken using a 5-point dither pattern in the Y (called “Z” in Keck/NIRC) and K filters with appropriate flat-field and dark current calibration frames also obtained.

## 3 Data Reduction

Our data reduction steps were performed within IRAF and include the standard procedures of applying a linearity correction, masking bad pixels, flat-fielding for correction of sensitivity variations across the array, sky-subtracting, and, finally, mosaicking the dither pattern to produce a final image for photometry. We describe these steps in detail for the P60/IRC observations though the same basic procedures were followed for the Keck/NIRC data.

Non-linearity of the detector array was assessed from a series of 2 to 20 second exposures of the back of the telescope’s primary mirror cover taken at the end of the observing run. A central section of the array (1/4 of the total area) that is relatively defect-free was used to measure ADU vs exposure time. The count rate became significantly non-linear above 36,000 ADU, but to be safe images with a maximum pixel value above 30,000 ADU were discarded.

Below 30,000 counts the response was almost, but not quite, linear. The advertised correction (on the Palomar 60” web-page) linearizes data according to:

$$xl=x+x^21.85\times 10^6(1+\delta /EXPTIME)$$

where $`xl`$ represents linear ADU and $`x`$ raw ADU. For our data, $`\delta `$ = 3.385 for an array speed of 1000 msec and 2.6 for a speed of 743 msec; hence integration times shorter than 5 seconds were avoided during the observations (sometimes by slightly defocussing the telescope when the observing plan called for observations of bright ($`<`$7.5) stars and integrations longer than 10 seconds preferred. Different linearization factors were examined and the correction with the lowest chi-squared relative to a linear curve was

$$xl=x+x^22.1\times 10^6$$

without a dependence on integration time for EXPTIME $`>`$5 seconds. Both the magnitude of this correction and its difference from the standard correction given above are quite small.

Raw bias and flatfield frames were used to verify that the gain and readnoise of the detector were similar to the values in Murphy et al. (gain = 8e-; readnoise = 30e-). Bias frames were acquired as 1.365 second exposures with the filter wheel in the closed position. Ten frames were averaged after rejecting among the frames the highest and lowest value for each pixel. On each night a sequence of twilight exposures through each filter was used to construct flat field images. The minimum and maximum values for each pixel were rejected and the individual frames averaged after subtracting the bias image. The sky flat was then normalized to have an overall average value equal to unity. A bad pixel list was created based on a histogram of the normalized flatfield images by flagging all pixels below 0.65 and above 1.2. This list had considerable overlap with the advertised bad pixel list and the two lists were merged into a single bad pixel mask.

The images of astronomical targets were proccessed first by constructing sky images as averages of the 9 data frames in each dither sequence and subtracting this sky image from each individual data frame, then by determining the relative offsets in fractional pixels between each component of the dither pattern for the point source, and finally by applying the offsets and combining the individual frames with bad pixel masking into a mosaic. These steps were performed in the NOAO/IRAF environment.

## 4 Photometry and Calibration

We performed aperture photometry with the IRAF/PHOT task. An aperture radius of 12 pixels (15” diameter on the sky) was chosen uniformly for all objects after looking at FWHM values (2-4) pixels and nearest neighbor distances for the entire dataset. The bright A0 stars which tended to be in crowded fields due to proximity to the galactic plane, set the maximum aperture size. A growth curve of magnitude vs aperture size peaked in the range 11-14 pixels, beyond which the photometric errors increased substantially as the magnitudes decreased for a sky annulus defined from 20-25 pixels. Only for a few (1-2) stars with FWHM just slightly less than 4 pixels (2.5”) were the fluxes possibly under-estimated using an aperture only $``$3 times their FWHM. The empirical aperture correction for these objects was difficult to find in many cases, but it was examined for several of the largest objects and found to be $``$0.007 mag. Due to the difficulty in determining this number and its relatively small magnitude this correction was not applied. Uncertainty in the derived instrumental magnitudes was underestimated by IRAF/PHOT. A better estimate was obtained by examining the photometry of sources in the individual unmosaicked frames where the statistical scatter was observed to be $``$0.01 mag. We therefore assumed a minimum photometric error of 0.01 mag and characterized our instrumental error as errorIRAF/PHOT+0.01 mag.

The relationship between airmass and (standard - instrumental) magnitude was used to derive the zeropoints and atmospheric extinction curves for the K and H filters (Figure 3). Persson et al. (1998) standard stars that were observed over the airmass range 1.0-2.5, with airmass taken from the middle frame in the 9-position dither pattern, defined these relationships. The star GSPC S813-D was observed at 5-6 different airmasses each night for the explicit purpose of deriving the atmospheric extinction curve for the Y filter without having to know the standard star magnitudes at Y (Figure 4). Because conditions were photometric, data from both nights were combined in order to obtain a more tightly constrained fit. We now describe the fitting procedure.

While the slopes from each of the two nights were approximately the same at Y, H, and K, the zero point for the second night was systematically higher by a few hundredth’s of a magnitude at all wavelengths. In combining data from the two nights, various vertical offsets of night 2 data to night 1 data were applied and the value which minimized the $`\chi ^2`$ in the linear fit of the extinction curve was retained. The linear fits assumed errors in airmass given by the 1/2 width of the airmass range of the observations, and errors in the ordinate given by a quadratic sum of our photometry errors (errorIRAF/PHOT+0.01 mag) plus an arbitrarily assumed 0.01 mag error in the standard star magnitudes.

We thus calibrated the K and H photometry according to:

$$K_{true}=K_{instr.}(0.065\pm 0.010)\times AM_K+\left\{\begin{array}{cc}20.354\pm 0.015(n1)\hfill & \\ 20.389\pm 0.015(n2)\hfill & \end{array}\right\}$$

$$H_{true}=H_{instr.}(0.029\pm 0.010)\times AM_H+\left\{\begin{array}{cc}20.925\pm 0.015(n1)\hfill & \\ 20.948\pm 0.015(n2)\hfill & \end{array}\right\}$$

as shown in Figure 3. The second term in each equation contains a numerical value and error for the atmospheric extinction, $`\kappa _{K,H}`$, and the third term values and errors for the zero point, $`ZP_{K,H}`$. The extinction coefficients we derive for Palomar, $`\kappa _K`$ = 0.065 and $`\kappa _H`$ = 0.029 agree well with those from the 2MASS survey 222http://www.ipac.caltech.edu/2mass/releases/first/doc/sec4\_8.html at Mt. Hopkins a few hundred miles away and a few thousand feet higher, $`\kappa _{Ks}`$ = 0.061 and $`\kappa _H`$ = 0.031. Note that our quoted zero points are for A.M.=0, i.e. outside the atmosphere, assuming linear extrapolation of the fit derived from A.M.$`>`$1; see Figure 3 of Tokunaga, Simons, & Vacca (2002) for the (in)appropriateness of this exptrapolation in various near-infrared bandpasses.

To calibrate the Y filter we assume by definition that unreddened A0 stars have zero color. The A0 stars observed for this program are lightly reddened due to the requirement that they be faint enough to observe with P60/IRC (implying K$`>`$7.5 mag) which is essentially a minimum distance requirement given the absolute magnitude of an A0 star, and hence a minimum reddening requirement given the properties of the local interstellar medium. We therefore include an interstellar extinction term in our calibration equations. We used both the K-band and the H-band data independently to derive the zeropoint at Y, from:

$$K_{instr.}+Y_{instr.}=ZP_KZP_Y(\kappa _K\kappa _Y)\times AM(A_K/A_VA_Y/A_V)\times A_V$$

and a similar equation extrapolating to the Y-band from H-band instead of K-band, with the difference between the two approaches being less than 0.015 mag and within the $`1\sigma `$ A0 star photometry errors. Yinstr. is the instrumental magnitude, ZPY is the zero point, $`\kappa _Y`$ is the slope of the atmospheric extinction curve, and $`A_Y`$ is the interstellar extinction. To find $`\kappa _Y`$ we combined the data for GSPC S813-D from both nights, shifting the night 2 data by 0.035 mag to obtain the best fit to a straight line, giving $`\kappa _Y=0.047`$. Figure 4 shows this slope determination. For the interstellar extinction term we use literature or Tycho B-V photometry and assume AV = 3.1 $`\times `$ E(B-V), AY = 0.38 $`\times `$ AV (anticipating our derivation of this result below in section 5.2), and A$`{}_{H}{}^{}=0.155\times A_V`$, and A$`{}_{K}{}^{}=0.090\times A_V`$ (Cohen et al., 1981). The two A0 stars with less well-determined B-V photometry (see Table 2) were left out of the final analysis since they introduced larger scatter and a systematic shift in the results between the H and K derivations.

With $`\kappa _Y`$ and $`A_Y`$ in hand, the above equation and an analogous one for H are solved for ZPY; errors are calculated by summing in quadrature the errors from each term of the equations. From this solution for ZPY we then correct the 10-12 data points for GSPC S813-D and combine these measurements with the A0 points to derive the final Y-band calibration equation:

$$Y_{true}=Y_{instr.}(0.047\pm 0.018)\times AM_Y+\left\{\begin{array}{cc}20.835\pm 0.030(n1)\hfill & \\ 20.875\pm 0.030(n2)\hfill & \end{array}\right\}.$$

Although the errors in the final K-band transformations are lower than those in the final H-band transformations, the zero point at Y derived from K and from H agree to well within the estimated errors, differing only by 0.015 mag. The $`\chi ^2`$ values for the fits are 0.8 and 0.7 for the extrapolation from K and H, respectively, a little less than unity which indicates a good assessment of errors. The difference in zero point between the two nights is 0.040 mag whereas we had shifted the data for S813-D between the two nights by 0.035 to derive the slope of atmospheric extinction, all consistent, and suggestive that only a single iteration of the procedure is needed.

The final Y, Y-H, and Y-K photometry for all objects is presented in Table 2 along with errors and ancillary information. These measurements are in the natural system of the camera, which should be that defined by the Persson et al. (1998) standards set up in part using this same camera. Magnitude uncertainties are given by:

$$\sigma _{magnitude}=\sqrt{\sigma _{ZP}^2+airmass^2\sigma _\kappa ^2+\sigma _{instrumental}^2}$$

which contains 1-$`\sigma `$ error terms for the zero point, atmospheric extinction assuming no uncertainty in the recorded airmass, and instrumental errors, $`\sigma _{instrumental}`$ = ($`\sigma _{IRAF/PHOT}`$+0.01 mag). Color errors are assumed to be the root-sum-square of the relevant magnitude errors which means that we are ignoring the effects of any correlated errors in our photometry. Multiple measurements of the same star were averaged, mostly occurring in the case of GSPC standards. Also given in Table 2 are optical photometry and spectral types from the literature, and 2MASS photometry from the Second Incremental Release for comparison with our H and K measurements.

## 5 Analysis

### 5.1 Colors

Color-color diagrams were made for different color combinations using the H and K results presented in this paper and B, V, and I photometry obtained from the literature where available. These plots were superposed on standard color-color lines (e.g. Bessell & Brett 1988 transformed to the CIT photometric system) and good agreement was seen between our H and K photometry and expected I-K and H-K colors. Moving now to the Y filter, Figure 5 shows the infrared Y-K color compared to optical B-V and V-I colors. These colors redden together with increasing spectral type to mid-M in the case of B-V and mid-L in the case of V-I, at which point the optical colors become saturated and/or turn bluer while the infrared color continues to grow redder. Note that optical colors are not available for all of the objects in Table 2, particularly the later type, redder stars and brown dwarfs.

Plots of Y-H versus H-K, Y-K versus H-K, and Y-H versus Y-K were then investigated, as seen in Figures 6 and 7. The range of spectral types extends from sdO through T5V with photometry for the T6.5V brown dwarf Gl229b (Matthews et al., 1996) shown for comparison. Gaps in the data around Y-K = 0.3 and 0.8 mag are due to the lack of F and K stars in our sample. These plots reveal a relatively tight relationship between the YHK colors as they grow redder together through spectral type $``$L6V. Beyond L6V the colors turn blueward and from L7V to at least T6.5V the colors are nearly degenerate with those of warmer L and M stars.

There are six objects in common between this study and that of Leggett et al. (2002): LHS 3003 (M7V), LHS 2924 (M9V), 2MASS J1439284+192915 (L1V), DENIS-P J0205.4-1159 (L7V), 2MASS 1632291+190441 (L8V), and SDSSp 015141.69+124429.6 (T1) Our Y-band magnitudes are 0.3-0.6 mag brighter than their UFTI Z magnitudes and our colors accordingly bluer for these objects, consistent with the nature of the filter center and width differences described in our introduction.

### 5.2 Extinction

In order to study the effects of interstellar extinction on colors involving the Y-band filter, a moderately reddened field towards the edge of the Ophiuchus dark cloud was observed. Eleven stars were identified in all three filter images, but only 7 were bright enough in Y to produce acceptable photometry. L 134 and L 547 from Persson et al. (1998) were also included in this analysis. The locations of these 9 stars in YHK color-color diagrams did not enable us to deduce the magnitude of the reddening vector, as they lay basically along an extrapolation of the color-color line defined by the dwarf stellar/sub-stellar sequence (see Figures 6 and 7).

In order to quantify how objects redden in Y-band relative to V-band we retrieved J, H, and Ks photometry for the 7 Ophiuchus and the 2 Persson “red objects” from the 2MASS survey’s Second Incremental Release. We then dereddened these 9 stars using the 2MASS J-H and H-Ks colors and the Cohen et al. (1981) interstellar extinction curve until they intersected the dwarf color sequence in a J-H vs H-Ks diagram over the color range corresponding to spectral types K7-M6. Fortunately these stars all appear to be reddened background stars and not young stellar objects associated with the Ophiuchus (and other, in the case of the Persson et al. “red stars”) clouds, although we can not prove this. If they were indeed located within clouds they would likely have circumstellar dust/gas and hence colors that reflect substantial near-infrared excess in addition to the effects of reddening, invalidating our approach. The estimated H-K and J-H color excess due to interstellar reddening was converted to $`A_V`$.

We then employed our Y-K vs H-K diagram to experiment with different values of $`A_Y`$ and find the value that best put the 9 reddened stars on the empirically defined color-color line. 2MASS gives only upper limits to the J magnitudes for 2 of the 9 stars, implying lower limit $`A_V`$’s derived from the J-H vs H-Ks diagram; these stars were given little weight in the $`A_Y`$ determination. The relation $`A_Y`$ = (0.38$`\pm `$0.03)$`\times A_V`$ was found to best bring the ensemble of stars onto the color sequence defined by unreddened stars. The derived value of total extinction at Y-band relative to V-band is within the rough expectations from interpolation between adjacent atmospheric windows, given $`A_J`$=0.265$`\times A_V`$ and $`A_I`$ = 0.617$`\times A_V`$. For comparison, $`A_{z(gunn)}`$=0.472$`\times A_V`$ and $`A_{z(sloan)}`$=0.453$`\times A_V`$ have been quoted in the literature (Schlegel et al. 1998). No information is available to the authors regarding AZ(UFTI).

## 6 Discussion

Soon after incorporation of the H-band into standard filter sets it was shown by Bahng (1969) that the colors of GKM stars differ from blackbody expectations in the near-infrared. The excursions from monotonic reddening with advancing spectral type are largely due to the influence of H- opacity in the H-band. Deviation from blackbody colors is apparent in the now-familiar J-H vs H-K diagram in the growth of J-H relative to H-K in a manner that is “faster” than blackbody from spectral types G through K, and then in the peak and blueward turnover of J-H near spectral type K7 as H-K continues slowly to get redder. After spectral type M4, J-H again grows redder with increasing H-K through $``$L6.

In our data we see that colors involving Y also steadily redden out to Y-J $``$ 1 mag, Y-H $``$ 2.5 mag, and Y-K $``$ 3 mag at spectral type L6V. Then, Y-H and Y-K colors both turn sharply blueward in the L6-L8V temperature range with later T-type objects having colors similar to those of warmer L and M stars. Y-J colors turn only softly blueward, remaining relatively constant from L6-T4V and then overlapping dwarf stars only as warm as M8-M9 by spectral type T5V.

Cool M- and L-type objects initially become redder with advancing spectral type because of increased opacity sources at the shorter Y-band, mostly H20 absorption at first, then H2-H2 collision induced absorption. When CH4 opacity begins to contribute substantially in the T dwarfs, the H and K fluxes decrease, sending the colors bluer. This bluing behavior is consistent with that exhibited in J-H vs H-K diagrams of late-type objects (e.g. Burgasser et al. 2002) where the transition between L-dwarfs and T-dwarfs takes colors from the reddest occupied by stars, all the way though the heart of the color-color diagram past every known astronomical object, to the bluest colors occupied by both hot stars and cool brown dwarfs. In the case of YHK colors, however, instead of crossing the entire color sequence the T dwarf sequence overlaps the only slightly warmer L and M dwarfs (see Figures 6 and 7). The latest T dwarfs ($`>`$T5) appear to fall more blueward in Y-K than Y-H and hence still may be uniquely distinguishable in this color-color plane. Overall, Y-J appears to be the most robust diagnostic of late-type cool stars and brown dwarfs.

The sharp turnaround in Y-K and Y-H occurring in the T dwarfs stands in contrast to T dwarf colors at longer and shorter wavelengths. K-L’ colors appear to redden out to at least L7, become flat with advancing spectral type from L7 to T1, and then redden again from T1 to at least T6 (Stephens et al., 2001; Leggett et al. 2002). However, hampering the utility of this color dignostic is that L-band data on stellar/sub-stellar photospheres is hard to come by due to high thermal backgrounds and resulting low sensitivity. I-z colors have been studied through only the L spectral types, with Steele & Howells (2000) claiming that I-z colors redden through the M and early L range, flatten from L1-L5, and then redden again towards even later spectral types. However, I-band photometry is not readily obtainable for many of the nearby L and T objects due to source faintness (MI = 15 mag at M9 but $`>`$18 mag by L7). We thus advocate the use of Y-K, Y-H, and Y-J colors for identifying very cool objects.

The downside of the Y-K vs Y-H diagram is the effect of extinction, which drives stars nearly directly along the dwarf locus. In other color-color planes such as I-Y or Y-J vs Y-K and J-K, however, there is more leverage of extinction vs the stellar/substellar color squence.

## 7 Conclusions

We present the first exploration of the colors of normal dwarf stars (and some brown dwarfs) spanning spectral types sdO through T5V using the Y-band at 1.035 $`\mu `$m. The behavior of the earth’s atmosphere is relatively clean in this window of near unit transparency and few emission lines, especially compared to the other near-infrared J, H, and K windows. Atmospheric extinction at Y-band is $`<`$0.05 mag per airmass. The relationship between interstellar extinction in the Y-band and that at V-band is $`A_Y`$ = (0.38$`\pm `$0.03) $`\times A_V`$.

Our photometric results, including errors, appear in Table 2. Although for most objects only one measurement was taken, photometric variability should not be an issue since our stars are mostly accepted photometric standards in either the optical or the near-infrared. The exception to this are the late type dwarfs (late M, L, and T spectral types) for which small-amplitude near-infrared variability has been reported in some objects (e.g. Bailer-Jones & Mundt, 2001).

The large range in color exhibited by low-mass stars and brown dwarfs makes colors involving the Y-band particularly diagnostic of these extremely cool objects. The I-Y, Y-J, Y-H, and Y-K colors all grow redder with increasing spectral type through at least L6V. At later spectral types YJHK colors all saturate and turn blueward to somewhat differing degrees, due to the strong peak in the SEDs of extremely cool objects in the 1.1-1.3 $`\mu `$m region. This bluing is caused by H2O, CH4, and collision-induced H2 opacity supressing first the K-band and then the H-band flux as temperatures grow cooler (e.g. Marley et al. 2002). Infrared spectroscopy of Gl229b and Gl570d, the coolest known T dwarfs, appear to confirm these statements for most of the T spectral type range. Colors involving the Y, J, H, and M bands can uniquely identify extremely cool brown dwarfs and so-called “planetary mass objects” apart from other astronomical targets, and should form the basis for deep surveys of field and cluster samples in the future.

We acknowledge with appreciation Adam Burgasser for obtaining P60/IRC (“X”) Y-band data for us of several T-dwarfs. LAH thanks John Carpenter for allowing her to talk him into observing with the Keck/NIRC (“Z”) Y-band filter during data acquisition for our H/K study of Orion Nebula Cluster faint stars and brown dwarfs. This publication makes use of data products from the Two Micron All Sky Survey (2MASS) and has benefited from information provided by the SIMBAD database.


text:
# 1 Introduction

## 1 Introduction

One of the common goals in the forthcoming experiments is the search for new phenomena. In the forthcoming high energy physics experiments (LHC, TEV22, NLC, …) the main goal is the search for physics beyond the Standard Model (supersymmetry, $`Z^{}`$-, $`W^{}`$-bosons, …) and the Higgs boson discovery as a final confirmation of the Standard Model. In estimation of the discovery potential of the planned experiments (to be specific in this paper we shall use as an example CMS experiment at LHC ) the background cross section (the Standard Model cross section) is calculated and for the given integrated luminosity $`L`$ the average number of background events is $`n_b=\sigma _bL`$. Suppose the existence of a new physics leads to additional nonzero signal cross section $`\sigma _s`$ with the same signature as for the background cross section that results in the prediction of the additional average number of signal events 111It should be noted that the existence of new physics can also lead to the decrease of the cross section due to destructive interference or some nonlocal formfactors. In this paper we consider the case when the new physics existence leads to additional positive contribution to the background cross section. The consideration of the opposite case is straightforward. $`n_s=\sigma _sL`$ for the integrated luminosity $`L`$.

The total average number of the events is $`<n>=n_s+n_b=(\sigma _s+\sigma _b)L`$. So, as a result of new physics existence, we expect an excess of the average number of events. In real experiments the probability of the realization of $`n`$ events is described by Poisson distribution

$$f(n;\lambda )=\frac{\lambda ^n}{n!}e^\lambda .$$
(1)
Here $`\lambda =<n>`$ is the average number of events. Remember that the Poisson distribution $`f(n;\lambda )`$ gives  the probability of finding exactly $`n`$ events in the given interval of (e.g. space and time) when the events occur independently of one another at an average rate of $`\lambda `$ per the given interval. For the Poisson distribution the variance $`\sigma ^2`$ equals to $`\lambda `$. So, to estimate the probability of the new physics discovery we have to compare the Poisson statistics with $`\lambda =n_b`$ and $`\lambda =n_b+n_s`$. Usually, high energy physicists use the following “significances” for testing the possibility to discover new physics in an experiment:

* “significance” $`S_1={\displaystyle \frac{n_s}{\sqrt{n_b}}}`$ ,
* “significance” $`S_2={\displaystyle \frac{n_s}{\sqrt{n_s+n_b}}}`$ ,
* “significance” $`2S_{12}=2(\sqrt{n_s+n_b}\sqrt{n_b})`$ .

A conventional claim is that for $`S_1(S_2)5`$ we shall discover new physics (here, of course, the systematic uncertainties are ignored). For $`n_bn_s`$ the significances $`S_1`$ and $`S_2`$ coincide (the search for Higgs boson through the $`h\gamma \gamma `$ signature). For the case when $`n_sn_b`$, $`S_1`$ and $`S_2`$ differ. Therefore, a natural question arises: what is the correct definition for the significance $`S_1`$, $`S_2`$ or anything else ?

It should be noted that there is a crucial difference between the planned experiment and the real experiment. In the real experiment the total number of events $`n_{obs}`$ is a given number (already has been measured) and we compare it with $`n_b`$ when we test the validity of the standard physics. So, the number of possible signal events is determined as $`n_s=n_{obs}n_b`$ and it is compared with the average number of background events $`n_b`$. The fluctuation of the background is $`\sigma _{fb}=\sqrt{n_b}`$, therefore, we come to the $`S_1`$ significance as the measure of the distinction from the standard physics. In the conditions of the planned experiment when we want to search for new physics, we know only the average number of the background events and the average number of the signal events, so we have to compare the Poisson distributions $`f(n;n_b)`$ and $`f(n;n_s+n_b)`$ to determine the probability to find new physics in planned experiment.

In this paper we describe a method for estimation of the discovery potential and exclusion limits on new physics in planned experiments. The effective significance of signal for given probability of observation is proposed for planned experiments instead of the usual significances $`S_1={\displaystyle \frac{n_s}{\sqrt{n_b}}}`$ and $`S_2={\displaystyle \frac{n_s}{\sqrt{n_s+n_b}}}`$, where $`n_s`$ and $`n_b`$ are the average numbers of signal and background events. We also estimate the influence of systematic uncertainties related to nonexact knowledge of signal and background cross sections on the probability to discover new physics in planned experiments. An account of such systematics is very essential in the search for supersymmetry at LHC.

The organization of the paper is the following. In the next Section mainly due to completeness we discuss the case of real experiment. In Section 3 we describe a method for the estimation of new physics discovery potential in planned experiment. Section 4 deals with estimation of exclusion limit. In Section 5 we estimate the influence of the systematics related to nonexact knowledge of the signal and background cross sections on the probability to discover new physics and set up exclusion limits on new physics in planned experiments. Section 6 contains concluding remarks.

## 2 New physics discovery in real experiment

In this section well known situation with real experiment is reminded to pedagogical reasons. Consider the case when the average number $`\lambda `$ of the events in the Poisson distribution (1) is big $`\lambda 1`$. In this case the Poisson distribution (1) approaches the Gaussian distribution

$$f_G(n;\mu ,\sigma )=_{n0.5}^{n+0.5}P_G(x;\mu ,\sigma ^2)𝑑x,$$
(2)
with $`P_G(x;\mu ,\sigma ^2)=\frac{1}{\sigma \sqrt{2\pi }}e^{\frac{(x\mu )^2}{2\sigma ^2}}`$, $`\mu =\sigma ^2`$, $`\mu =\lambda `$ and $`n0`$. Note that for the Poisson distribution the mean equals to the variance. According to common definition  new physics discovery corresponds to the case when the probability that background can imitate signal is less than $`5\sigma `$ (here of course we neglect any possible systematic uncertainties). Suppose we have observed some excess of events $`n_{obs}>n_b`$. The probability that for the background $`n_b`$ we shall observe events with $`nn_{obs}`$ is determined by standard formula

$$P(nn_{obs})=_{n_{obs}}^{\mathrm{}}P_G(x;n_b,n_b)𝑑x=\frac{1}{\sqrt{2\pi }}_{S_1}^{\mathrm{}}exp(x^2/2)𝑑x,$$
(3)
where

$$S_1=\frac{n_{obs}n_b}{\sqrt{n_b}}\frac{n_s}{\sqrt{n_b}}$$
(4)
According to common definition for $`S_15`$ the Standard Model is excluded (the probability that background can imitate signal is less than $`2.8510^7`$) and we have new physics discovery 222Here we neglect any possible systematic uncertainties..

Suppose some model with new physics predicts $`\lambda =n_s+n_b`$ and

$$|n_{obs}n_sn_b|k_1\sigma _{s+b},$$
(5)
then for $`k_1=1.64`$ ($`k_1=1.96`$) the model with “new physics” agrees with experimental data at 90% C.L. ( 95% C.L.). Here $`\sigma _{s+b}^2=n_s+n_b`$.

For $`S_12`$ ($`S_13`$ ) in formula (3) the probability that background can imitate signal is less than 2.28% ( 0.14%) and according to our definition we have weak (strong) evidence in favor of new physics.

Suppose that the measured number of events is such that

$$|n_{obs}n_b|(k_1=1.96)(\sigma _b=\sqrt{n_b})$$
(6)
It means that the Standard Model agrees at 95 % C.L. with experimental measurement. In this case we can also obtain exclusion limit on new physics (limit on the average number of signal events $`n_s`$) . Namely, for $`\lambda =n_s+n_b`$ we require that

$$|n_{obs}n_sn_b|(k_1=1.96)(\sigma _{s+b}=\sqrt{n_s+n_b})$$
(7)
From the equation (7) we obtain 95% C.L. upper bound on the average number $`n_s`$ of signal events.

Consider now the case of the Poisson distribution (1). Suppose we have measured the number of events $`n_{obs}>n_b`$ (an excess of events). We define the statistical significance $`s`$ of a signal  in the Standard Model by

$$\frac{1}{\sqrt{2\pi }}_s^{\mathrm{}}exp(x^2/2)𝑑x=\underset{k=n_{obs}}{\overset{\mathrm{}}{}}f(k;n_b).$$
(8)
The formula (8) is nothing but the probability to observe $`nn_{obs}`$ of background events in an identical independent experiments. Note that $`s`$ is a function on $`n_{obs}`$ and $`n_b`$$`s=s(n_{obs},n_b)`$. If $`s5`$ then by common definition we have new physics discovery. For $`s2`$ ($`s3`$) according to our definition we have weak(strong) evidence in favor of new physics. If the model with additional $`n_s`$ signal events obeys the inequality

$$|s(n_{obs},n_s+n_b)|(k_1=1.96)$$
(9)
then the model with new physics agrees at 95% C.L. with an experiment.

Suppose that we have observed the number of events compatible at 95% C.L. with the Standard Model, i.e.

$$s(n_{obs},n_b)(k_1=1.96).$$
(10)
In this case one can obtain at 95% C.L. exclusion limit on the average number of signal events $`n_s`$ from the inequality

$$s(n_{obs},n_s+n_b)(k_1=1.96).$$
(11)
## 3 Planned experiments.

As it has been mentioned in the introduction the crucial difference between planned experiment and real experiment is that in real experiment we know the number of observed events, therefore we can compare the Standard Model with experimental data directly, whereas in the case of planned experiment we know only the average number of background events $`n_b`$ and the average number of signal events (for the case when we have new physics in addition to the Standard Model). Therefore in the case of planned experiment an additional “input” parameter is the probability of the discovery. Suppose we test two models: the Standard Model with the average number of events $`\lambda =n_b`$ and the model with new physics and the average number of events $`\lambda =n_s+n_b`$.

To discover new physics we have to require that the probability $`\beta (\mathrm{\Delta })`$ of the background fluctuations for $`n>n_0(\mathrm{\Delta })`$ is less than $`\mathrm{\Delta }`$, namely

$$\beta (\mathrm{\Delta })=\underset{n=n_0(\mathrm{\Delta })+1}{\overset{\mathrm{}}{}}f(n;n_b)\mathrm{\Delta }$$
(12)
The probability $`1\alpha (\mathrm{\Delta })`$ that the number of events in a model with new physics will be bigger than $`n_0(\mathrm{\Delta })`$ is equal to

$$1\alpha (\mathrm{\Delta })=\underset{n=n_0(\mathrm{\Delta })+1}{\overset{\mathrm{}}{}}f(n;n_s+n_b)$$
(13)
It should be stressed that if $`\mathrm{\Delta }`$ is a given number then $`\alpha (\mathrm{\Delta })`$ is a function of $`\mathrm{\Delta }`$ or vice versa we can fix the value of $`\alpha `$ in formula (13) then $`\mathrm{\Delta }`$ is a function of $`\alpha `$. The meaning of the probability of the discovery $`1\alpha `$ is the probability that in the case of new physics an experiment will measure the number of events bigger than $`n_0`$ such that the probability that the Standard Model can reproduce such number of events is rather small ($`\beta `$).

In other words we choose the critical value $`n_0`$ for hypotheses testing333A simple statistical hypothesis $`H_0`$ (new physics is present, i.e. $`\lambda =n_s+n_b`$) against a simple alternative hypothesis $`H_1`$ (new physics is absent, i.e. $`\lambda =n_b`$. about observability of new physics requiring that Type II error $`\beta \mathrm{\Delta }`$. Then we calculate the Type I error $`\alpha `$ and the probability of discovery (or evidence) $`1\alpha `$.

For fixed value of $`\alpha `$ and known values of $`n_s`$, $`n_b`$ we can calculate $`\beta `$ using formulae (12,13). In our numerical calculations we take $`\alpha =0.5;0.25;0.1;0.05`$. Consider now the limiting case $`n_b1`$ when Poisson distribution approaches Gaussian distribution. The equations (12,13) take the form

$$\beta _{n_0}^{\mathrm{}}P_G(x;n_b,n_b)𝑑x$$
(14)
$$1\alpha _{n_0}^{\mathrm{}}P_G(x;n_s+n_b,n_s+n_b)𝑑x$$
(15)
Consider at first the most simple case when $`\alpha =0.5`$ (see Figs.1-2 for an illustration). For $`\alpha =0.5`$ parameter $`n_0`$ in formula (15) is equal to $`n_0=n_s+n_b`$. The equation (14) takes the form

$$\beta _{S_1}^{\mathrm{}}P_G(x;0,1)𝑑x,$$
(16)
where

$$S_1=\frac{n_s}{\sqrt{n_b}}$$
(17)
The significance $`S_1`$ is determined by the formula (17) and it is often used in experiment proposals .

For $`1\alpha >0.5`$ (see Fig.3 for an illustration) the parameter $`n_0`$ in formula (14) is equal to

$$n_0=n_s+n_bk(\alpha )\sqrt{n_s+n_b},$$
(18)
where $`k(\alpha )`$: $`k(0.5)=0`$; $`k(0.25)=0.66`$; $`k(0.1)=1.28`$; $`k(0.05)=1.64`$ (as an example, Tab.28.1 ). The effective significance $`s`$ in the equation (8) (i.e. corrected significance $`S_1`$, corresponding the discovery probability $`1\alpha `$) has the form

$$s=\frac{n_s}{\sqrt{n_b}}k(\alpha )\sqrt{1+\frac{n_s}{n_b}}$$
(19)
So, we see that the asymptotic formula (17) for the significance $`s`$ is valid only for $`1\alpha =0.5`$.

As it has been shown in refs.  the more proper of the significance in planned experiments is $`2S_{12}`$. The generalization of this significance to the case of $`1\alpha >0.5`$ looks very attractive for approximate estimation of discovery potential

$$s=2(\sqrt{n_s+n_b}\sqrt{n_b})k(\alpha ).$$
(20)
The comparison of formulae (19,20) is shown in Fig.4.

It should be stressed that very often in the conditions of planned experiment the average numbers of background and real events are not very big and we have to solve the equations (12, 13) directly to construct $`5\sigma `$ discovery, $`3\sigma `$ strong evidence and $`2\sigma `$ weak evidence curves. Our numerical results are presented in Figs. (5 - 10).

As an example consider the search for standard Higgs boson with a mass $`m_h=110GeV`$ using the $`h\gamma \gamma `$ decay mode at the CMS detector. For total luminosity $`L=310^4pb^1(210^4pb^1)`$ one can find  that $`n_b=2893(1929),n_s=357(238)`$, $`S_1={\displaystyle \frac{n_s}{\sqrt{n_b}}}=6.6(5.4)`$. Using the formula (19) and Table of the standard normal probability density function  we find that $`1\alpha (\mathrm{\Delta }_{dis})=0.93(0.60)`$. It means that for total luminosity $`L=310^4pb^1(210^4pb^1)`$ the CMS experiment will discover at $`5\sigma `$ level standard Higgs boson with a mass $`m_h=110GeV`$ with a probability 93(60) percents 444In other words let us suppose that we have constructed 100 identical CMS detectors. At $`5\sigma `$ level the Higgs boson will be discovered at 93(60) CMS detectors.

For the case when we are interested in estimation of the lower bound on number $`n_s`$ of signal events (bound on new physics) we can use the equations

$$1\alpha (\mathrm{\Delta })=\underset{n=n_0(\mathrm{\Delta })+1}{\overset{\mathrm{}}{}}f(n;n_b+n_s)$$
(21)
$$\beta (\mathrm{\Delta })=\underset{n=n_0(\mathrm{\Delta })+1}{\overset{\mathrm{}}{}}f(n;n_b)\mathrm{\Delta }$$
(22)
## 4 Exclusion limits.

It is important to know the range in which the planned experiment can exclude the presence of signal at the given confidence level ($`1ϵ`$). It means that we will have uncertainty in future hypotheses testing about non observation of signal equals to or less than $`ϵ`$. In refs. different methods to derive exclusion limits in prospective studies have been suggested. As is seen from Fig.11 the essential differences in values of the exclusion limits take place. Let us compare these methods by the use of the equal probability test . In order to estimate the various approaches of the exclusion limit determination we suppose that new physics exists, i.e. the value $`n_s`$ equals to one of the exclusion limits from Fig.11 and the value $`n_b`$ equals to the corresponding value of expected background. Then we apply the equal probability test ($`f(n_0;n_s+n_b)=f(n_0;n_b)`$) to find critical value $`n_0`$ for hypotheses testing in planned measurements (Fig.12). Here a zero hypothesis $`H_0`$ is the statement that new physics exists and an alternative hypothesis $`H_1`$ is the statement that new physics is absent. After the calculation of the Type I error $`\alpha `$ (the probability that the number of the observed events will be equal or less than the critical value $`n_0`$) and the Type II error $`\beta `$ (the probability that the number of the observed events will bigger than the critical value $`n_0`$ in the case of the absence of new physics) we can compare the methods. For this purpose the relative uncertainty  $`\stackrel{~}{\kappa }={\displaystyle \frac{\alpha +\beta }{2(\alpha +\beta )}}`$ which will take place under hypotheses testing $`H_0`$ versus $`H_1`$ is calculated. This relative uncertainty $`\stackrel{~}{\kappa }`$ in case of applying the equal-probability test is a minimal relative value of the number of wrong decisions in the future hypotheses testing for Poisson distributions. It is the uncertainty in the observability of the new phenomenon. Note the $`1\stackrel{~}{\kappa }`$ (the relative number of correct decisions) is a distance between two distributions (the measure of distinguishability of two Poisson processes) in frequentist sense.

In Table 1 the result of the comparison is shown. As is seen from this Table the ”Typical experiment” approach  gives too small values of the exclusion limit. The difference in the 90% CL definition is the main reason of the difference between our results and the results from ref. . We require that $`\stackrel{~}{\kappa }`$ equals to $`ϵ`$, i.e. we use only one parameter ($`ϵ`$) as measure of uncertainty in the hypotheses testing. In ref  the criterion $`\beta <\mathrm{\Delta }`$ and $`{\displaystyle \frac{\alpha }{1\beta }}<ϵ`$ for determination of the exclusion limits has been applied. It means the experiment will observe with probability at least $`1\mathrm{\Delta }`$ at most a number of events such that the limit obtained at the $`1ϵ`$ confidence level excludes the corresponding signal. In this case two parameters ($`ϵ,\mathrm{\Delta }`$) are required to construct exclusion limits. Nevertheless, we have got results close to  555The using $`\kappa =\alpha +\beta `$ as measure of uncertainty  gives a somewhat different results..

## 5 An account of systematic uncertainties related to nonexact knowledge of background and signal cross sections

In the previous sections we took into account only statistical fluctuation in the number of signal and background events and did not take into account other uncertainties. There is considered the systematic uncertainties 666In ref.  the systematic uncertainty is the uncertainty in the sensitivity factor. This uncertainty has statistical properties which can be measured or estimated. The systematic effects in ref.  as supposed has stochastic behaviour too. The account for statistical uncertainties due to statistical errors in determination of values $`n_b`$ and $`n_s`$  implies the existence of conditional probability for parameter of Poisson distribution. We consider here forthcoming experiments to search for new physics. In this case the systematic uncertainties has theoretical origin without any statistical properties.  due to imperfect knowledge of the background and signal cross sections.

In this section we investigate the influence of the systematic uncertainties related to nonexact knowledge of the background and signal cross sections on the discovery potential in planned experiments. Denote the Born background and signal cross sections as $`\sigma _b^0`$ and $`\sigma _s^0`$. An account of one loop corrections leads to $`\sigma _b^0\sigma _b^0(1+\delta _{1b})`$ and $`\sigma _s^0\sigma _s^0(1+\delta _{1s})`$, where typically $`\delta _{1b}`$ and $`\delta _{1s}`$ are $`O(0.5)`$ for the LHC. Two loop corrections for most reactions at present are not known. So, we can assume that the uncertainty related with nonexact knowledge of cross sections is around $`\delta _{1s}`$ and $`\delta _{1b}`$ correspondingly. In other words we assume that exact cross sections lie in the intervals $`(\sigma _b,\sigma _b(1+\delta _b))`$ and $`(\sigma _s,\sigma _s(1+\delta _s))`$, where $`\sigma _b`$ and $`\sigma _s`$ are calculated at Born or one loop level of the accuracy. The average number of background and signal events lie in the intervals $`(n_b,n_b(1+\delta _b))`$, $`(n_s,n_s(1+\delta _s))`$, where $`n_b=\sigma _bL`$, $`n_s=\sigma _sL`$.

To determine the new physics discovery potential we again have to compare two Poisson distributions with and without new physics. Contrary to the Section 3 we have to compare the Poisson distributions in which the average numbers lie in some intervals. It means that we have to find the critical value $`n_0`$ and to estimate the influence of systematic uncertainty on the discovery probability. A priori the only thing we know is that the average number of background and signal events lie in some intervals but we do not know the exact values of the average background and signal events. Moreover we can not say anything about probability distributions of possible values $`n_b`$ and $`n_s`$ in this interval. Such distribution is absent.

An account of uncertainties related to nonexact knowledge of background and signal cross sections is straightforward and it is based on the results of Section 3. Suppose uncertainty in the calculation of exact background cross section is determined by parameter $`\delta `$, i.e. the exact cross section lies in the interval $`(\sigma _b,\sigma _b(1+\delta ))`$ and the exact value of average number of background events lies in the interval $`(n_b,n_b(1+\delta ))`$. Let us suppose $`n_bn_s`$. In this instance the discovery potential most sensitive to the systematic uncertainties. Because we know nothing about possible values of average number of background events, we consider the worst case. Taking into account formulae (12) and (13) we have the formulae 777Formulae (23,24) realize the worst case when the background cross section $`\sigma _b(1+\delta )`$ is the maximal one, but we think that both the signal and the background cross sections are minimal

$$\beta (\mathrm{\Delta })=\underset{n=n_0(\mathrm{\Delta })+1}{\overset{\mathrm{}}{}}f(n;n_b(1+\delta ))\mathrm{\Delta }$$
(23)
$$1\alpha (\mathrm{\Delta })=\underset{n=n_0(\mathrm{\Delta })+1}{\overset{\mathrm{}}{}}f(n;n_b+n_s)$$
(24)
This approach allows estimate the scale of influence of background uncertainty to observability of signal. As an application of formulae (23,24) consider the case $`n_b=n_s=100`$ (typical case for the search for supersymmetry at LHC). For such values of $`n_s`$ and $`n_b`$ and for $`\delta `$= 0., 0.1, 0.25, 0.5 we find that $`1\alpha (\mathrm{\Delta }_{dis})`$ = 0.9996, 0.9924, 0.8476, 0.137, correspondingly. So, we see that the uncertainty in the calculations of background cross section is extremely essential for the determination of the LHC discovery potential. Some other examples are presented in Tables 2-7 and in Figs.13-14.

## 6 Conclusions

In this paper we have described a method to estimate the discovery potential and exclusion limits on new physics in planned experiments where only the average number of background $`n_b`$ and signal events $`n_s`$ is known. We have found that in this case the more proper definition of the significance (for $`\alpha =0.5`$) is $`2S_{12}=2(\sqrt{n_s+n_b}\sqrt{n_b)}`$ in comparison with often used expressions for the significances $`S_1={\displaystyle \frac{n_s}{\sqrt{n_b}}}`$ and $`S_2={\displaystyle \frac{n_s}{\sqrt{n_s+n_b}}}`$. For $`1\alpha >0.5`$ we have additional additive contribution to the significance (see approximate formulae 19-20 of estimation of discovery potential for given discovery probability). As a result, the effective significance $`s`$ of signal for given probability of observation is proposed. The results of direct calculations of dependences $`n_s`$ versus $`n_b`$ for different discovery probabilities and significances are presented. We estimate the influence of systematic uncertainty related to nonexact knowledge of signal and background cross sections on the probability to discover new physics in planned experiments. An account of such kind of systematics is very essential in the search for supersymmetry and leads to an essential decrease in the probability to discover new physics in future experiments. The texts of programs can be found in http://home.cern.ch/bityukov.

We are grateful to Fred James, V.A.Matveev and V.F.Obraztsov for the interest and useful comments. S.B. thanks Bob Cousins, George Kahrimanis, James Linnemann, Louis Lyons, Tony Vaiciulis, Alex Read, Byron Roy and Pekka Sinervo for very useful discussions. S.B. would like to thank James Stirling, Mike Whalley and Linda Wilkinson for having organized this interesting Conference which is the wonderful opportunity for exchange of ideas. This work has been supported by CERN-INTAS 99-0377. The work of S.B. is also supported by CERN-INTAS 00-0440.