text:
# Bayesian detection of planetary transits

## 1 Introduction

The search for rocky, terrestrial planets around other stars is a key research topic in astrophysics for the next decade. Following the first exo-planet detection around a sun-like star (Mayor & Queloz, 1995), gaseous giants around other solar-type stars have been shown to be relatively common (Butler et al., 2001). The mass function of the current crop of extra-solar planets grows rapidly toward the lower masses (Butler et al., 2001), showing that low-mass planets must be common. However, the radial velocity technique, which has resulted in the detection of the exo-planets detected so far, is limited to planetary masses somewhat smaller than Saturn, and cannot reach the domain of terrestrial planets. This is due to astrophysical effects, such as microturbulence in the star’s atmosphere, rather than instrumental limitations.

The most promising approach for the detection of (significant numbers) of terrestrial planets around stars other than the Sun appears to be the search for planetary transits, i.e. dips in the light curve of the parent stars caused by the planet transiting in front of the stellar disk. The flux dip caused by the transit is also small, $`\mathrm{\Delta }F/F=(R_\mathrm{p}/R_{})^2`$, which for the transit of an Earth-Sun system gives $`\mathrm{\Delta }F/F=10^4`$. This is well below the scintillation noise caused by the Earth’s atmosphere (Favata & the Eddington Science Team, 2000, see), so that high-accuracy space-based photometry will be needed for the detection of such events. The probability of occurrence of a transit depends on the inclination of the planetary orbit relative to the line of sight (which must be close to $`i=90`$ degrees), and is relatively small (for a set of randomly oriented Sun-Earth systems $`p0.5\%`$), so that searches for planetary transits must be based on observation of large samples of target stars. A typical transit duration will be of order $`\mathrm{\Delta }t10`$ hours, and the transit periodicity will be the same as the orbital period of the planet, typically several months.

A number of space missions wholly or partially dedicated to the search for planetary transits are either in development or in the planning stage. The CNES/European satellite COROT is planned for launch in 2005, while the ESA mission *Eddington* and the NASA mission *Kepler* are planned for launch in 2007. Given the intrinsically statistical nature of planetary transit searches, these missions will acquire large number of stellar light curves, ranging from thousands for COROT to hundreds of thousands for *Eddington* and *Kepler*. Also, some smaller searches are being conducted for limited time periods (and concentrating on larger planets) using e.g. HST (Gilliland et al., 2000) or ground-based telescopes (e.g. Doyle et al., 2000).

The analysis of data from such searches, and in particular the detection of transits with a high degree of certainty and a low false alarm rate, is a challenging task. The transit signal is weak ($`\mathrm{\Delta }F/F=10^4`$), and concentrated in a small fraction of the total signal: for a habitable planet orbiting a K5V star the orbital period will be roughly 4 months, so that for a 1 year light curve three events will be present. As each transit lasts $``$ 10 hours, the transit signal is present in only $`0.3\%`$ of the total light curve. In the Euclidean regime, the number of stars in a given field increases toward fainter magnitudes by a factor of $`4`$ per magnitude. This is the case for the range of magnitudes and the low Galactic target latitudes of interest for currently planned missions . Therefore, most of the detected planets will be in the light curves of the fainter (and thus statistically noisier) stars, impying the need for effective robust data analysis algorithms able to reliably detect transits “hidden in the noise”. At the same time, the large number of light curves which will need to be analyzed, each with a large number of points (of order $`\mathrm{10\hspace{0.17em}000}`$ points for a year of data) makes the use of efficient algorithms necessary, and rules out brute force approaches.

Some ground- (Doyle et al., 2000) and HST-based (Gilliland et al., 2000) transit searches, which deal with relatively small numbers of light curves, use a detection approach based on comparing large numbers of model transits to the light curves and minimising a $`\chi ^2`$ statistic (or a linear statistic in the case of Doyle). These approaches are computationally very intensive, and thus may be unsuitable for the routine processing of the large number of light curves which will be produced by upcoming space missions.

As an alternative, transit detection algorithms based on Bayesian methods have recently been the subject of some attention. They have the advantage of combining computational efficiency with flexibility. While a global statistic can be used for the detection, information is directly available to reconstruct the detected signal if wanted, therefore providing a tool to discriminate between planetary transits and other types of periodic signals (Defaÿ et al., 2001a), as well as directly measuring additional planetary characteristics such as the planet’s radius.

In the present paper we present a novel algorithm for the detection of planetary transits based on the method developed by Gregory & Loredo (1992) (hereafter referred to as GL method) for the search of pulsed emission from pulsars in X-ray data. While the algorithm was developed to be “general purpose”, we have tuned it with the parameters of the upcoming *Eddington* planet finding mission in mind. The present paper discusses the characteristics of the algorithm on the basis of extensive simulations for the case in which the light curve is dominated by photon noise. Its performance in the case in which stellar activity is the dominating noise source will be the subject of a future paper.

Bayesian algorithms for the detection of planetary transits are also being developed in the context of the COROT mission. In particular, an approach based on expansion of the light curve into a truncated Fourier series is being investigated (Defaÿ et al., 2001b). Perfoming the detection in the Fourier domain can make the algorithm computationally sensitive to data gaps and sampling rates. Here we explore a more robust *direct space* approach.

The GL (Gregory & Loredo, 1992) method, was initially developed for the detection of X-ray pulsars (where Poisson statistics dominate) and later extended to the Gaussian noise case (Gregory, 1999). At the flux levels of interest for the transit searches for *Eddington*, the photon shot noise per detection element (which is Poissonian) can be very well represented by Gaussian noise (see Sect. 2.1). The original formulation of the GL algorithm is well-suited to the detection of periodic signals of unknown shape. However, in the planetary transit problem we have strong prior information about the transit shape. In this paper we modify the GL algorithm to perform more optimally for planetary transit detection. We do this by allowing one of the bins to have a variable width, to represent the out of transit constant signal level. This formulation also permits the phase of the transits to be identified, a task the original GL method is not suited for (see Sect. 3.1). The fitted parameters are the period, duration and phase of the transit. The shape of the transit can then be reconstructed from the phase-folded light curve.

The simulated light curves are described in Sect. 2. The algorithm is outlined in Sect. 3 and compared with the original GL algorithm in Sect. 4. Sect. 5 describes the evaluation of the algorithm’s performance by determining the number of false alarms and missed detections in a large sample of simulated light curves with and without transits. Conclusions and options for future work are presented in Sect. 7.

## 2 The light curves

### 2.1 Transits

Given the presence of limb darkening in stellar photospheres, planetary transits are not perfectly “flat bottomed” (nor are they, strictly speaking, truly grey). To simulate transits in a realistic way, the Universal Transit Modeler (UTM) software written by H.J. Deeg (Deeg, 1999) was used. Limb darkening coefficients were taken from Van Hamme (1993). Two types of transits were simulated, one representing a Jupiter-type planet in a short orbit around a Sun-like star and another representing a Earth-like planet in an habitable orbit around a K5V star.

The input characteristics of the system for the Jovian transit were:

* Time step $`t_{\mathrm{unit}}=15`$ minutes.
* Radius of star $`R_{}=R_{}`$.
* Luminosity of star $`L_{}=L_{}`$.
* Radius of planet $`R_\mathrm{p}=R_{\mathrm{}}`$. The ratio of the duration of the ingress / egress, to that of the ‘flat bottom’ of the transit (affected only by limb darkening) is roughly $`2\mathrm{R}_\mathrm{p}/(\mathrm{R}_{}\mathrm{R}_\mathrm{p})`$, in this case $`0.22`$.
* Period of transit $`P=2880\times t_{\mathrm{unit}}`$ ($`1`$ month). This is at the short end of the range of periods of interest for *Eddington*, but extrapolation to longer periods is to some extent possible (when activity is not included) by acting on the number of transits in the light curve.
* The distance star-planet, which is used to determine the transit duration, was varied between $`d=15.3\times R_{}`$ (resulting in a transit duration of 15 hours) and $`d=45.9\times R_{}`$ (5 hours). (N.B.: for the same period different distances would correspond to different planetary to stellar mass ratios).
* The duration of light curve was varied between $`D=3\times P`$ and $`D=5\times P`$.
* The phase of the transit was randomly varied in the different simulations. The posital phase (between 0 and 1) is used in the course of the present paper.

Light curves were normalized to the photon count level expected for a star of a given $`V`$ magnitude (between 7 and 17 typically), based on the throughput expected for the baseline *Eddington* mission design (Favata & the Eddington Science Team, 2000), i.e. a collecting area of $`0.6`$ m2 and a total system throughput of $`70\%`$. With these instrument parameters a $`V=21.5`$ G2V star will yield $`50`$ detected photons/sec. Gaussian noise was then added with variance defined by the number of detected photons per pixel.

## 3 The algorithm

### 3.1 A Bayesian method

The method employed consists in calculating the likelihood of the data given a certain number of parameters, varying the parameters over a given range and identifying the value of each parameter whose probability is maximized according to Bayes’ theorem:

$$p(\theta |data,I)=p(\theta |I)\times \frac{p(data|\theta ,I)}{p(data|I)}$$
(1)
where :

* $`\theta `$ is a set of parameter values (i.e. a hypothesis).
* $`data`$ is the dataset.
* $`I`$ represents information about the ensemble of hypotheses considered i.e. the type of model used and knowledge about the other models. For the remainder of this section $`I`$ will be implicit in likelihood expressions.
* $`p(data|I)`$ is a prior for the type of model used.
* $`p(\theta |I)`$ is the combined prior for the parameters.

An excellent description of the theory on which the present algorithm is based is given in Gregory & Loredo (1992). In the present paper we will give a brief outline of the calculations, detailing only those aspects in which our work differs from the discussion of Gregory & Loredo (1992). As a starting point we constructed an algorithm following exactly the GL prescription, and we tested it on sets of 10 simulated light curves containing transits with varying characteristics. This benchmark was later used to ensure that the modifications in our algorithm were indeed improvements.

The GL algorithm employs a family of stepwise models to describe the periodic signal plus background. Each member of the family resembles a histogram, with $`m`$ equal sized bins per period $`P`$. The family members are distinguished by the value of $`m`$ which is varied in the range from 2 to some upper limit (typically 15 for X-ray pulsar detection work). Such a model is capable of approximating a light curve of essentially arbitrary shape, which is desirable for detecting periodic signals of unknown shape, in contrast to the current planetary transit problem, for which the shape is known a priori. GL also employs a phase parameter $`\varphi `$. If the time offset $`o`$ is defined as the time elapsed between the start of the first bin and the start of the data, $`\varphi =2\pi (o/P)`$. The parameters fitted by the GL method are then $`P`$, $`m`$, $`\varphi `$ (the flux level in each bin of the model is marginalised over). In the case of planetary transits, it is not desirable to let the model vary outside the transit. We therefore have a slightly different type of model. The number of steps in the step function is $`n+1`$. Bin $`0`$ is the ‘out of transit’ bin and lasts for a large fraction of any given period, and bins $`1`$ to $`n`$ are ‘in transit’, each lasting $`d/n`$ where $`d`$ is the duration of the model transit. We have also adopted a different definition for the time offset, as we have a significant event – the transit – which we can use to determine the start of a new period. Defining the time offset $`o`$ as the time from the start of the data to the start of the next transit, the phase is then related to the offset in the same way as before. The parameters required are now $`P`$, $`\varphi `$ and $`d`$. These parameter definitions are illustrated for both methods in Fig. 1

As there is no feature in a step function of unconstrained shape with equal duration steps which can mark the beginning of a period in the data, the concept of phase is not well defined. Any bin in the step function could be the first. Thus we do not expect the GL method to enable phase determination directly. Only after reconstruction of the entire light curve could the position of the transit be pin-pointed relative to the start of the data. By introducing a transit feature in the model, the phase is built into the model function and we expect it to be detected effectively by the modified algorithm.

Models with a lower number $`n`$ of ‘in transit’ bins will incur a lower Occam penalty factor, as emphasised in Gregory & Loredo (1992). In general, $`n`$ should be chosen to be the lowest value possible. For pure detection purposes, given that transits are relatively simple events, $`n=1`$ should suffice. For transit reconstruction purposes, a higher value can be used.

Despite the modifications we made to the GL models, we followed the method outlined in Gregory (1999) to calculate the likelihoods.

### 3.2 Likelihood calculation

The likelihood is initially calculated for a given set of parameters $`P`$ (period), $`d`$ (duration), $`o`$ (offset). For convenience the results were sometimes expressed in terms of posital phase: $`\mathrm{ph}=\varphi /2\pi =o/P`$.

Due to the different type of model function used, Eq. (6) in Gregory (1999), which describes the assigmnent of a bin number $`j`$ to each data point $`y_i`$ taken at time $`t_i`$, was replaced by the following:

$$j(t_i)=\{\begin{array}{ccc}t_{\mathrm{mod}}\hfill & :& \text{if }0<t_{\mathrm{mod}}n\hfill \\ 0\hfill & :& \text{otherwise}\hfill \end{array}$$
(2)
where:

$$t_{\mathrm{mod}}=\mathrm{int}\left(\frac{(t_i+Po)\mathrm{mod}(P)}{d/n}+1\right)$$
(3)
$`n`$ is the number of bins per transit, $`\mathrm{int}(x)`$ is the nearest integer lower than or equal to $`x`$ and $`(a)\mathrm{mod}(b)`$ is the remainder of $`a`$ divided by $`b`$.

At time $`t_i`$, the observed flux count $`y_i`$ can be written as $`y_i=y(t_i)+e_i`$ where $`y(t_i)`$ is the value predicted by the model for time $`t_i`$ and $`e_i`$ is a noise component. The noise is assumed to have a Gaussian distribution (see Gregory, 1999 and references therein) with variance $`\sigma _i^2`$. In the present case it is appropriate and clearer to use the same value of $`\sigma `$ for all data points111In Gregory (1999) a noise parameter $`b`$ is introduced to account for incomplete knowledge of $`\sigma `$, and is then marginalised over. We have not made use of this parameter in this work.. Strictly speaking the noise in the *Eddington* case is Poisson distributed (being photon shot noise), however given the large number of photons in each time bin used for the transit search this is indistinguishable from a Gaussian noise distribution. The likelihood is therefore given by:

$$p(data|P,d,o)=\underset{i=1}{\overset{N}{}}\left[\frac{\sigma ^1}{\sqrt{2\pi }}\times \mathrm{exp}\left[\frac{\left(y_iy(t_i)\right)^2}{2\sigma ^2}\right]\right]$$
(4)
where $`N`$ is the total number of data points.

Re-expressed in terms of the $`n+1`$ bins of the model:

$$\begin{array}{ccc}p(data|P,d,o)\hfill & =& \sigma ^N(2\pi )^{N/2}\hfill \\ & \times & \underset{j=0}{\overset{n}{}}\mathrm{exp}\left[\left(\underset{i=1}{\overset{n_j}{}}y_ir_j\right)^2/2\sigma ^2\right]\hfill \end{array}$$
(5)
where $`n_j`$ is the number of data points in bin $`j`$ and $`r_j`$ is the model value in bin $`j`$.

As shown in Gregory (1999) the argument of the exponential can be reduced to:

$$\frac{\left({\displaystyle \underset{i=1}{\overset{n_j}{}}}y_ir_j\right)^2}{2\sigma ^2}=\frac{W_j\left(r_j\overline{d_{W_j}}\right)^2+\chi _{W_j}^2}{2}$$
(6)
This allows the marginalization over the $`r_j`$’s to be performed, which we do identically to Gregory, to obtain:

$$\begin{array}{ccc}p(data|P,d,o)\hfill & =& \sigma ^N(2\pi )^{N/2}\left(\mathrm{\Delta }_r\right)^{(n+1)}\left(\frac{\pi }{2}\right)^{(n+1)/2}\hfill \\ & & \\ & \times & \mathrm{exp}\left[\underset{j=0}{\overset{n}{}}\chi _{W_j}^2/2\right]\hfill \\ & & \\ & \times & \underset{j=0}{\overset{n}{}}W_j^{1/2}\left[\text{erfc}(y_{j,\mathrm{min}})\text{erfc}(y_{j,\mathrm{max}})\right]\hfill \end{array}$$
(7)
where:

* $`\mathrm{\Delta }_r=r_{\mathrm{max}}r_{\mathrm{min}}`$ is the range of values the model step function is allowed to take;
* the quantities $`W_j`$, $`\chi _{W_j}^2`$, $`y_{j,\mathrm{min}}`$ and $`y_{j,\mathrm{max}}`$ are taken directly from Eqs. (11) to (16) in Gregory (1999) ;
* $`\text{erfc}(y)`$ is the complementary error function.

### 3.3 Odds ratio calculation

In order to use the likelihood to determine a given parameter all the other parameters must be marginalized, by multiplying by the corresponding prior and integrating over the parameter’s range of values. When marginalizing the phase, in order to minimise the computing time, we incremented the phase by steps of $`\pi d/P`$ where $`P`$ is the period and $`d`$ the duration of the transit, corresponding to time offset increments of $`d/2`$. We used the same priors for each parameter as Gregory, with a flat prior for the new parameter $`d`$. Although we worked in terms of period rather than frequency this does not change the calculations.

Odds ratios $`O_{P,c}`$ were then computed by comparing the probabilities as defined in Eq. (7) for a range of period values, integrating out all other parameters, to the probability obtained with a constant model denoted by the subscript $`c`$. These odds ratios can then used to check for evidence of a periodic signal over the entire frequency range before proceeding to determine individual parameters, as described in Gregory (1999). The a posteriori probabilities needed for parameter estimation can be directly evaluated from the odds ratios by multiplying by the relevant prior and normalizing.

In Gregory & Loredo (1992), a global odds ratio is calculated for each light curve by marginalising over all the parameters, in order to determine whether there is evidence for a periodic signal. If the global odds ratio is larger than $`1`$, the answer is yes. In that case, posterior probability distributions for individual parameters are used to determine the optimal parameter values.

### 3.4 Weighting factor to compensate for uneven distribution into the bins

When the number of periods is low such that one bin might be represented four times while another only three times, or if there are gaps in the data which may not be evenly distributed over the bins, Gregory & Loredo (1992) noted that some of their initial assumptions may fail, leading to the appearance of an erroneous trend in the posterior probability for the period.

In an appendix to Gregory & Loredo (1992), a solution to this problem was proposed. A weighting factor $`s_j`$ is applied to each bin:

$$s_j=\left(\frac{n_jm}{N}\right)^{n_j}$$
(8)
This factor is derived in the context of Poisson statistics and does not apply to the present, Gaussian noise case.

Despite the low number of periods in our light curves it was found that no weighting factor was required in the benchmark algorithm that reproduced the GL identically. However it is clear that the problem is more acute in the modified algorithm. The ‘out of transit’ bin contains many more data points than the others, and therefore has a much larger effective weight. A weighting factor is required to compensate for this problem. The expression given above for $`s_j`$ is only appropriate in the photon count context in which it was derived, not in the Gaussian noise case adopted here. A different weighting factor can be heuristically derived by considering Eq. (6). The contribution of each model level to the likelihood is a $`\chi ^2`$ sum. The variance of a $`\chi ^2`$ distribution is given by the number of degrees of freedom $`\nu `$. In each bin there are $`n_j`$ data points and $`\mathrm{n}_{\mathrm{param}}`$ parameters to adjust. As $`n_j\mathrm{n}_{\mathrm{param}}`$, $`\nu =n_j\mathrm{n}_{\mathrm{param}}n_j`$. Weighting each bin by a factor $`1/n_j`$ is therefore equivalent to weighting by the variance. In practice this is achieved by maintaining the expressions for $`\overline{d_{W_j}}`$ and $`\overline{d_{W_j}^2}`$, given in Gregory (1999) in terms of $`d_i`$ and $`\sigma `$, but replacing $`W_j`$ by $`W_j/n_j`$.

This modification was implemented in our algorithm and found to give more robust results.

### 3.5 Minimizing the computing time

For a given set of parameters, the calculation of the likelihood involves summing over each element in each bin. The time required to compute the likelihood for a given set of $`P`$, $`d`$, $`o`$ therefore scales linearly with the number of points in the light curve. It also increases with the number of bins, but this is a slow increase. It does not depend on the individual parameter values.

The overall computing time also depends, of course, on how tightly the parameter space is sampled. It is necessary to minimise the number of trial values for each parameter without missing potentially localised likelihood maxima. Because of the relative sharpness of the peak in the posterior probability for the period, the period increment needs to be kept fairly small (typically once or twice the time step between data points). Attention was therefore concentrated on what increment was suitable in terms of phase. The results are not significantly worsened by increasing the *posital* phase increment from $`1/P`$ (i.e. shifting the model by 1 sampling time at each increment) to $`d/2nP`$ (i.e. shifting the model by half the duration of an in-transit bin at each increment). Further increase leads to sharp steps in the posterior probability distribution (analogous to Shannon’s sampling theorem).

However, the computing time is inversely proportional to the increment, and the steps in the distribution are effectively removed by dividing it by the equivalent distribution for an entirely flat light curve with the same duration, sampling and data gaps as the light curve. We call this dividing function the “window function” 222This also has the advantage of ironing out any residual effects of the uneven bin duration not removed by the weighting factor.. We therefore used a phase increment of $`d/2P`$ and performed the division before analyzing the results. As the window function only needs to be calculated once per set of parameters, this is much faster than using a smaller increment (see Sect. 5).

Note that due to the use of this window function one should not strictly speaking use the word ’posterior probability’ when talking about the output of the algorithm. In the rest of this paper we will refer to ‘modified posterior probability’ to mean ‘posterior probability distribution divided by the window function’. This also implies that the global odds ratios mentioned in Sect. 3.3 cannot be used to directly measure the ratio of the probabilities for a periodic model compared to a constant model. Instead, we use bootstrap simulations (see Sect. 5.1) to set a threshold value of the detection statistic above which a detection is accepted.

## 4 Comparing the modified algorithm with the original version

In order to establish a reference point and to have a preliminary estimate of the modified algorithm’s performance, some tests were run on both the original and the modified version. Table 1 summarizes the names and meanings of the various parameters in each method.

### 4.1 A few qualitative tests

From a typical light curve described below a number of parameters were varied one by one and the odds ratios were plotted as a function of period and as a function of phase. The base light curve lasted $`\mathrm{11\hspace{0.17em}500}\times 15`$ minutes (119.8 days), contained a transiting giant planet with a period of 30 days, a duration of 15 hours and a posital phase (i.e. phase in radians divided by $`2\pi `$) of 0.25. The magnitude of the parent star was $`10.0`$, which for *Eddington* corresponds to a signal to noise ratio of roughly $`1400`$ over 15 minutes, so that the depth of the transit for the Jupiter-sized planet is 14 times the noise standard deviation. It was analyzed with $`m=10`$ in the case of the GL method, and $`n=4`$ in the case of our method333The possibility of using $`n=1`$ for detection only purposes, then a larger value of $`n`$ for transit reconstruction, will be the subject of investigations in a further paper.. In order to sample the transit as well with the GL method as with the modified method, a much higher value of $`m`$ would need to be used, but this would be too computationally expensive. Instead the values of $`m`$ and $`n`$ we chosen such that the computing times were similar. The results obtained for this benchmark case are shown in Fig. 2.

Each of the parameters (be they associated with the light curve or with the model) was varied over a small range of representative values. These one-off tests on a small parameter space confirmed some expected trends:

* for a given light curve duration the detection is less precise for longer periods as the light curve contains less transits;
* as expected, the unmodified GL method is not well suited to detecting the phase as there is no way of labeling one bin the first one. A detection is still possible by folding the posterior probability for the phase over the number of bins used. On the other hand the phase is very successfully recovered with the modified version, and the precision does not vary with the phase itself;
* the larger the value of $`m`$ (GL method), the sharper the detection. However $`m=10`$ appeared sufficient for our purposes;
* increasing the value of $`n`$ (modified method) does not necessarily improve the detection ability since one starts to fit the noise inside the transits, which is not periodic. When fitting Gaussian profiles it is standard to require a minimum of $`2`$ bins per FWHM. The shape of the transit is not Gaussian but it is relatively simple, hence we multiplied by a safety factor of $`2`$, leading to $`n=4`$ in further calculations. However when dealing with a particular value of $`d`$ it is advantageous to choose $`n`$ so $`d`$ is a multiple of it to avoid introducing extra noise by splitting individual data points across bin boundaries(see footnote 3);
* although the modified method should in principle allow us to determine the duration of the transit, in practice this is not successful. The program may be fitting a much wider region than the transit itself. In the GL method, as there are only 10 to 20 bins per period, with $`P`$ of order several hundred sampling times or more, the bin in which the transit falls is much larger than the transit itself. We have seen that the loss of information this implies does not prevent the detection of the period by the GL method. The modified algorithm is likely to overestimate the transit duration because fitting a region larger than the transit does not significantly reduce the likelihood. For now the duration of the transit was simply marginalized over; once the presence of a transit is asserted and its period known, phase folding should allow a fairly quick determination of the shape and duration;
* for a given set of parameters, with $`m=10`$ and $`n=4`$, such that both algorithms have similar computing times, the detection peaks are much sharper with the modified version.

## 5 Performance of the algorithm

### 5.1 Method

To evaluate the performance of the algorithm, we used the same method as Doyle et al. (2000). For each set of trial parameters the algorithm was run first on a set of one hundred simulated light curves containing only Gaussian noise and no transits. Subsequently it was run on another set of one hundred simulated light curves containing Jovian-type planetary transits with the characteristics described in Sect. 2.1, with the same level but different realizations of the photon noise, and with uniformly distributed random phases

For each simulation, the modified posterior probabilities were plotted versus period and the value of the maximum was noted. This maximum is our ‘detection statistic’, on the basis of which we determine whether there is a transit or not. We then plot a histogram of the detection statistics measured from running the algorithm over all the light curves with transits and one histogram for all the light curves with noise only. In other words, one histogram corresponds to the cases where the transit hypothesis is correct and one to the cases where the null hypothesis is correct. Ideally, the two distributions would be completely separated, with no overlap, and choosing a detection threshold located between the two histograms would guarantee a 100% detection rate and a 0% false alarm rate. In practice, for the cases of real interest, close to the noise level, the two histograms will show an overlap. A compromise has to be found by choosing a threshold which minimises a penalty factor designed to take into account both false alarm and missed detection rate. This is illustrated in Fig. 3.

Depending on the circumstances, it may be more important to minimise the false alarm rate than the missed detection rate. This is the approach followed by Jenkins et al. (2002), on the basis that detections from space experiments are hard to follow-up from the ground. An alternative view is any real transit that is rejected is a loss of valuable scientific information. As long as the false alarm rate is kept to a manageable level, further analysis of the light curves will prune out the false events. We have opted here for an intermediate position, and our penalty factor is simply the sum of the missed detection rate $`\mathrm{N}_{\mathrm{MD}}`$ and the false alarm rate $`\mathrm{N}_{\mathrm{FA}}`$:

$$\mathrm{F}_{\mathrm{penalty}}=\mathrm{N}_{\mathrm{FA}}+\mathrm{N}_{\mathrm{MD}}$$
(9)
However, the marginalised detection algorithm yields modified posterior probabilities as a function of period, and also as a function of phase. The simultaneous use of the two detection statistics $`\mathrm{S}_{\mathrm{per}}`$ and $`\mathrm{S}_{\mathrm{ph}}`$ (plotting 2-D rather than 1-D distributions) increases the discriminating power of the algorithm, (as long as the two distributions do not have secondary maxima in 2-D space). This is shown when comparing the false alarm and missed detection rates obtained from period and phase information separately and together. The threshold in the 2-D case takes the form of a line: $`\mathrm{S}_{\mathrm{ph}}=a+b\times \mathrm{S}_{\mathrm{per}}`$. Here the optimal values of $`a`$ and $`b`$ were found by trial and error, although standard discriminant analysis techniques can be used to determine them automatically.

### 5.2 Results

#### 5.2.1 An ideal case

In Defaÿ et al. (2001a), analysis performed on the basis of 200 bootstrap samples for the COROT observations of a star with magnitude 13 and an Earth-sized planet showed, with 6 transits lasting 5 hr each, a probability of true detection of around 0.3. We performed the simulations described in Sect. 5.1 for a similar case: Earth-sized planet orbiting a K5V type star with $`V=13`$ with a period of $`932\times 15`$ minutes and a transit duration of 5 hr. The light curve is sampled with 15 minute bins. The noise is different from the COROT case, as we concentrate uniquely on the photon noise expected for *Eddington*.

The results are shown in Fig. 4 for period and phase separately. As the distributions for the noise only and transit light curves are completely separated, each parameter alone is sufficient to determine a threshold ensuring null false alarm and missed detection rates.

#### 5.2.2 Performance of the algorithm at the noise limit

Given that the key scientific goal of *Eddington* in the field of planet-finding is the detection of habitable planets, the performance of the algorithm was extensively tested for habitable planets at (or close to) the noise limit of *Eddington*. The case of a Earth-like planet orbiting a K dwarf in a habitable orbit was used as benchmark. The light curve was simulated for a system with the following parameters:

* the star is a K5 dwarf ($`R_{}=0.8`$ $`R_{}`$) with a range of apparent V-band magnitudes $`V=14.0,14.5,15.0`$;
* planet with radius $`R_\mathrm{p}=R_{}`$ and a period of 4 months, orbiting the star at a distance of $`0.64`$ A.U. (leading to a transit duration of $`10.5`$ hours);
* light curve duration of 16 months, containing 4 transits. The light curves were sampled every hour.

An example of a light curve is shown in Fig. 5. The resulting transit event has a depth $`\mathrm{\Delta }F/F=1.4\times 10^4`$. For the *Eddington* baseline collecting area a star at $`V=14`$ will result in a photon count of $`1.8\times 10^8`$ per hour, so that the Poisson noise standard deviation will be $`1.34\times 10^4`$. The $`S/N`$ of the transit event in each 1 hour bin will thus be $`1.88`$. Following the same reasoning for the $`V=15`$ case, the $`S/N`$ of the of transit event in a single one hour bin is $`1.19`$. As there are 4 transits lasting 10 hours each in the light curves considered, the overall transit signal has a $`S/N`$ of $`\sqrt{40}\times 1.197.5`$.

With the results of the simulations, an example of which is shown in Fig. 6, the analysis described in Sect. 5.1 was performed for all three magnitudes, confirming that the combined use of the two statistics improves the results. This is illustrated for the $`V=14.5`$ case in Figs. 7 & 8 (for this particular case 1000 rather than 100 runs were computed to improve precision).

As illustrated in Fig. 9, a mean error rate444i.e. the mean of the false alarm and missed detection rates. of $`<`$ 3% can be achieved up to magnitude 14.5. This magnitude is therefore taken as the performance limit for the algorithm for an Earth-sized planet around a K5V-type star. However this analysis is not complete enough to allow a precise determination of the limit, as the noise treatment is incomplete (photon noise only being considered) and one would need more runs per simulations to compute meaningful errors on the false alarm and missed detection rates (sets of 1000 runs, as was done for the limiting $`V=14.5`$ case, should be computed for all cases).

The asymetric shape of the distributions shown in Figs. 4, 7 & 8 implies that, even though the thresholds are chosen to minimise false alarms and missed detections equally, the optimal threshold results in more false alarms than missed detections. This could easily be avoided, if needed, by replacing Eq. 9 by:

$$\mathrm{F}_{\mathrm{penalty}}=A\times \mathrm{N}_{\mathrm{FA}}+\mathrm{N}_{\mathrm{MD}}$$
(10)
where $`A`$ is a factor greater than $`1`$. Alternatively one could keep the penalty factor unchanged but set a strict requirement on the maximum acceptable false alarm rate.

As in any unbiased search for periodicity in a time-series, the inclusion of a larger range of periods in the search will lead to a higher chance of finding a spurious (noise-induced) period signal in the data. The simulations used here to assess the algorithm’s performance are based on a search through a relatively small range of periods. In practice, lacking any a priori knowledge of the possible periodicity of planetary orbits around the star being observed, one will have to test a large range of periods, ranging from few days (the physical limit of the period of planetary orbits) all the way to the duration of the data set (searching for individual transit events).

#### 5.2.3 Data gaps

Any realistic data set will suffer from gaps in the data. While the orbits of both *Eddington* and *Kepler* have been chosen to minimize gaps, 100% availability is not realistic, and gaps will be present due to e.g. telemetry dropouts, spacecraft momentum dumping maneuvers, showers of solar protons during large solar flares, etc. For this reason any realistic algorithm must be robust against the presence gaps in the data, showing graceful degradation as a function of the fraction of data missing from the time series.

We have therefore tested the algorithm discussed here using simulated light curves with 5%, 10% and 20% data gaps, randomly distributed in the data, i.e. 5% of the points in the time series are selected randomly with a uniform distribution and removed from the light curve. The gaps will probably not be randomly distributed in reality, but as the typical gap duration is expected to be of order 1 or 2 hours, simulated random gaps can already be used to test the algorithm’s robustness. For reasons of computing time, to avoid having to recalculate the “window function” at each run, the distribution of the data gaps is the same for all runs of a simulation. As the gaps are chosen one by one there are rarely gaps of more than two consecutive time steps, i.e. 2 hours. Note that e.g. the *Eddington* mission is designed to produce light curves with a duty cycle $`90\%`$, so that the case with 20% data gaps represents a worst case analysis.

The results are shown in Fig. 10. There is visibly very little degradation up to 20% data gaps. When using $`\mathrm{S}_{\mathrm{per}}`$ alone or the two statistics combined there is no perceptible difference. We can therefore say this this algorithm is robust at least for data gaps of the type likely to occur due to e.g. telemetry dropouts, which last only a few hours. One would also expect the algorithm to perform well in the presence of longer gaps: the effect of gaps is to render the number of samples per bin uneven, and this is already the case for this particular method with no gaps at all.

#### 5.2.4 Number of transits in the data

The planetary transits detection phase of the *Eddington* mission is planned to last 3 years with a single pointing for the entire duration of that phase. There will therefore be three or four transits in the light curve for a typical habitable planet. However, other missions such as COROT are planned with shorter (5 months) pointings and its is of interest for this type of mission to study the degradation of the algorithm’s performance as the number of transits in the light curve reduces. If the algorithm performs well with 2 or less transits, in the context of *Eddington* it may also allow the detection of “cool Jupiters”, i.e. Jupiter-sized planets with orbits more similar to those of the gaseous giants in our solar system. This would be of relevance to the question of how typical our solar system is.

Sets of 100 runs with the characteristics specified in Sect. 5.2.2 for a star of magnitude 14.5 were computed for light curve durations of 4, 8, 12, 16 and 20 months, containing between 1 and 5 transits. The results are shown in Fig. 11. The degradation only becomes significant when less than three transits are present. However, even mono-transits could be detectable for larger planets at that magnitude.

Defaÿ et al. (2001a) compared a matched filter approach with a Bayesian method based on the decomposition of the light curve into its Fourier coefficients. Their results suggest that the performance degradation in the low number of transits case is faster for the Bayesian method than for the matched filter. This is because the matched filter makes use of assumptions about the transit shape. It is also shown that when the Bayesian method fails to detect a transit, it can still reconstruct it if the detection is performed using a matched filter. Our algorithm has not been directly compared to a matched filter. Its very design is based on the search for a short periodic signal in an otherwise flat lightcurve, which is itself an assumption about the shape of the signal. The matched filter makes use of more detailed knowledge of the transit shape and is therefore likely to perform better in the low transit number limit. However our algorithm with $`n=1`$ may provide already a very good approximation to the relatively simple shape that is a transit, and therefore perform nearly as well.

#### 5.2.5 Differences in the two statistics

The two a posteriori probabilities show a different behavior. In general the phase statistic is far more discriminatory than the period statistic. The period statistic’s lesser effectiveness may be explained in the following way. If the phase is wrong, even if the period is right, it is likely none of the transits will be matched. If the phase is right, whatever the period, at least the first transit will be matched by the model. First we consider the likelihood distribution a function of phase, normalised over all periods. For an incorrect phase the contribution from the correct period is nil as all transits are missed, but for the correct phase all trial periods produce a non-negligible contribution (the correct period of course contributing most). The likelihood distribution as a function of phase is therefore sharply peaked. Then we consider the likelihood distribution as a function of period, normalised over all phases. The contribution from the correct phase is non-negligible whatever the period. When the period is correct, the contribution from the correct phase is washed out by the contributions from all the incorrect phases. The likelihood distribution as a function of period is therefore less sharply peaked.

However the combined use of the two parameters is more successful than the phase statistic alone. The reason for this is illustrated in Fig. 8: in 2-D space the two distributions are aligned on a diagonal, such that no single value cutoff is optimal in either direction, compared to the line shown. In an upcoming paper, the direct use of a combined statistic shall be investigated. The global odds ratio described in Sect. 3.3 could be used for such a purpose. We have noted in Sect. 3.5 that the global odds ratio for a given lightcurve cannot be used as an absolute statitstic in the context of the present method. It can however be used as relative detection statistic, like $`\mathrm{S}_{\mathrm{per}}`$ & $`\mathrm{S}_{\mathrm{ph}}`$, combined with bootstrap simulations.

## 6 Discussion

Efficient data processing is one of the challenges for the upcoming generation of large scale searches for exo-planets through photometric transits. While radial velocity searches concentrate on limited number of stars, transit searches will investigate simultaneously large numbers of stars, and produce large amounts of data (photometric light curves) for each of them. A computationally efficient and robust algorithm for the processing of these data sets is necessary to make transit searches feasible. It is likely that the photometric time series which represent the observational product of the transit searches will be analyzed in different stages, using more than a single approach. In particular, a first level of processing (after instrumental effects have been removed) should concentrate on singling out high-probability transit candidates, while efficiently pruning out the large number (more than 90%, even if all stars have planets, due to the low probability of transit events) of light curves in which no transits are present. In this first stage of analysis the ability to efficiently screen real transits in the data – even at the price of a moderate number of false alarms – is a key requirement for the algorithm. The candidate light curves in which a transit is suspected will then later be subject to a more detailed processing, which can then afford to be computationally less efficient (given it has to operate on a much smaller amount of data).

The algorithm we have developed and discussed here is able to detect transit events at the limit of the photon noise present in the light curve. It shows a graceful degradation of its performance as function of different parameters of interest, e.g. the noise level in the data, as well as the presence of data gaps and the number of transits actually observed. Its strong sensitivity to the phase of periodic transits supplies significant additional information to be then used by further steps of processing for e.g. the reconstruction of the transit parameters. Thus, while little used in astronomy, Bayesian algorithms appear to be a powerful tool in the processing of transit data.

## 7 Conclusions and future work

A novel algorithm to detect transits due to extra-solar planets in stellar light curves has been developed and tested. The algorithm, based on a Bayesian approach, has proved successful in the tests performed so far, which include the effects of photon noise and data gaps. Using the photometric accuracy and throughput expected for the *Eddington* mission, we are able to detect an Earth-sized planet orbiting a K5V-type star with a period of 4 months down to an apparent stellar magnitude of $`V14.5`$. Randomly distributed data gaps lasting up to two hours each and covering up to 20% of the light curve do not significantly affect the performance of the algorithm. The minimum number of transits in one light curve required for high confidence detections is three, however the algorithm’s performance degrades gracefully for small number of transits, so that detections are possible for individual transits, albeit at a lower confidence level. This will allow for the detection of larger planets in long-period orbits (analogous to the gaseous giants of our solar system), likely to transit only once in the three year planet detection phase planned for the *Eddington* mission.

The most serious additional noise source to perturb planetary transit detections from space, is likely to be intrinsic stellar micro-variability (mostly activity-induced). At the moment it is also the least well investigated. The consequences of activity on the detection efficiency (using simulated light curves based on the solar light curves recorded by the VIRGO instrument on board SOHO, which spans all solar activity levels, from solar minimum to solar maximum) will be the subject of a future paper, in which the feasibility and effectiveness of using color information, as well as a number of pre-processing techniques such as whitening, will also be investigated.

The algorithm we have developed and discussed here has the potential to form part of a powerfull, multi-stage approach to analysing transit lightcurves. A more optimised processing method will be discussed in a separate paper. It will include a variability filtering stage, followed by distinct detection and parameter estimation stages, using a combination of a matched filter approach and of the present algorithm.

The performance of the algorithm presented here shows that the search of planetary transits with amplitudes comparable to the intrinsic noise level of the data set is fully feasible, and thus represents an important element in the development of the future generation of transit-based planet finding missions.

###### Acknowledgements.

This work was supported by the European Space Agency’s Young Graduate Trainee program. The authors are grateful to Hans Deeg for making his transit modelling package UTM available. We wish to thank Mike Irwin, David Katz & Rolf Jansen for advice, and Celine Defaÿ for advice and access to her algorithm. We are also grateful for the thorough and helpful comments of the referee.


format: text

text:
# Renormalization group approach to energy level statistics at the integer quantum Hall transition

## I Introduction

It has been realized long ago that, alongside with the change in the behavior of the eigenfunctions, a localization-delocalization transition manifests itself in the statistics of the energy levels. In particular, as the energy is swept across the mobility edge, the shape of the level spacing distribution (LSD) crosses over from the Wigner-Dyson distribution, corresponding to the appropriate universality class, to the Poisson distribution.

Moreover, finite-size corrections to the critical LSD exactly at the mobility edge allow to determine the value of the correlation length exponent,ShkSSL93  thus avoiding an actual analysis of the spatial extent of the wave functions. For this reason, the energy level statistics constitutes an alternative to the MacKinnon-KramerPicS81a ; PicS81b ; MacK81 ; MacK83  and to the transmission-matrixLan70 ; FisL81  approaches to the numerical study of localization.

Another reason why a large number of numerical simulations of the LSD at the transition HofS94b ; Eva94 ; VarHSP95 ; KawOSO96 ; BatSZK96 ; ZhaK97 ; Met98b ; Eva95 ; SchZ95 ; FeiAB95 ; OhtO95 ; BatS96 ; MetV98 ; BatSK98 ; Met99  were carried out during the past decade is the controversy that existed over the large-spacing tail of the critical LSD. Conclusive demonstrationBatSZK96 ; ZhaK97 ; BatSK98  that this tail is Poissonian, i.e., that there is no repulsion between the levels with spacings much larger than the mean value,ShkSSL93  rather than super-Poissonian,KraLAA94  implying that repulsion is partially preserved, required a very high accuracy of the simulations.ZhaK97 ; BatSK98  The bulk of numerical work on the level statistics at the transition was carried out for three-dimensional systems HofS94b ; Eva94 ; VarHSP95 ; KawOSO96 ; BatSZK96 ; ZhaK97 ; Met98b  for which there exists a mobility edge separating localized and extended states. In two dimensions all the states are localized in the absence of a magnetic field. In the presence of a magnetic field, localization-delocalization transitions in two dimensions (quantum Hall transitions) are infinitely sharp. Still the reasoning of Ref. ShkSSL93,  applies. Numerical studies have established a Poissonian tail of the LSD.BatS96  It was also demonstrated BatS96  that the procedure of extracting the localization length exponent from the finite-size corrections yields a value close to $`\nu =2.35`$ found from large-size simulations of the wave functions.HucK90 ; Huc92 ; Huc95 ; LeeWK93

Recently, a semianalytical description of the integer quantum Hall transition, based on the extension of the scaling ideas for the classical percolation StaA92  to the Chalker-Coddington (CC) model of the quantum percolation,ChaC88  has been developed.GalR97 ; AroJS97  The key idea of this description, a real-space-renormalization group approach (RG), is the following. Each RG step corresponds to a doubling of the system size. The RG transformation relates the conductance distribution of the sample at the next step to the conductance distribution at the previous step. The fixed point of this transformation, yields the distribution of the conductance, $`P_\text{c}(G)`$ of a macroscopic sample at the quantum Hall transition. This universal distribution describes the mesoscopic properties of a fully coherent quantum Hall sample. Analogously to the classical percolation,StaA92  the correlation length exponent, $`\nu `$, was extracted from the RG procedure CaiRSR01  using the fact that a slight shift of the initial distribution with respect to the fixed-point, $`P_\text{c}(G)`$, drives the system to the insulator upon renormalization. Then the rate of the shift of the distribution maximum determines the value of $`\nu `$. Remarkably, both $`P_\text{c}(G)`$ and the critical exponent obtained within the RG approachCaiRSR01 ; WeyJ98 ; JanMMW98  agree very well with the “exact” results of the large-scale simulations.Huc92 ; LeeWK93 ; WanJL96 ; WanLS98 ; AviBB99

The goal of the present paper is two-fold. Firstly, we extend the RG approach to the level statistics at the transition in order to subject its validity to yet another test. Secondly, we apply the method analogous to the finite-size-corrections analysis to extract $`\nu `$ from the LSD obtained within the RG approach. This method yields $`\nu =2.37\pm 0.02`$, which is even closer to the most precise large-scale simulations result $`\nu =2.35\pm 0.03`$ Huc92  than the value $`\nu =2.39\pm 0.01`$ inferred from the conductance distribution.CaiRSR01  The latter result is by no means trivial. Indeed, the original RG transformation CaiRSR01  related the conductances, i.e., the absolute values of the transmission coefficients of the original and the doubled samples, while the phases of the transmission coefficients were assumed random and uncorrelated. In contrast, the level statistics at the transition corresponds to the fixed point in the distribution of these phases. Therefore, the success of the RG approach for conductances does not guarantee that it will be equally accurate quantitatively for the level statistics.

Within both RG transformations, for the magnitudes and for the phases of the transmission coefficients, an initial deviation from the critical distribution drives the system towards an insulator with zero transmission and Poissonian LSD. Thus the procedures of the extraction of $`\nu `$ from both transformations are technically different, but conceptually similar. In fact, the shape of the critical LSD, obtained from the RG approach, shows systematic deviations from the large-scale simulation resultsFeiAB95 ; OhtO95 ; BatS96 ; MetV98 ; BatSK98 ; Met99 ; OnoOK96 ; KleM97  which yield the body of LSD very close to the Gaussian unitary random matrix ensemble (GUE).Meh90  However, the RG flow of the LSD towards the insulator appears to be robust.

The paper is organized as follows. First, in Sec. II we review the real-space RG approach GalR97 ; AroJS97 ; CaiRSR01  and adjust it to the computation of the energy levels and the LSD. In Sec. III we present our numerical results for the LSD. The finite-size scaling (FSS) analysis of the obtained LSD at the QH transition is reported in Sec. IV. Concluding remarks are presented in Sec. V.

## II Model and RG method for the LSD

### II.1 RG approach to the conductance distribution

A detailed description of the RG approach to the conductance distribution can be found in Refs. GalR97, ; AroJS97, ; CaiRSR01, . It is based on the RG unit shown in Fig. 1. The unit is a fragment of the CC network consisting of five nodes. Each node, $`i`$, is characterized by the transmission coefficient $`t_i`$, which is an amplitude to deflect an incoming electron along the link to the left. Analogously, the reflection coefficient $`r_i=(1t_i^2)^{1/2}`$ is the amplitude to deflect the incoming electron to the right. Doubling of the sample size corresponds to the replacement of the RG unit by a single node. The RG transformation expresses the transmission coefficient of this effective node, $`t^{}`$, through the transmission coefficients of the five constituting nodesGalR97

$$t^{}=\left|\frac{t_1t_5(r_2r_3r_4e^{i\mathrm{\Phi }_2}1)+t_2t_4e^{i(\mathrm{\Phi }_3+\mathrm{\Phi }_4)}(r_1r_3r_5e^{i\mathrm{\Phi }_1}1)+t_3(t_2t_5e^{i\mathrm{\Phi }_3}+t_1t_4e^{i\mathrm{\Phi }_4})}{(r_3r_2r_4e^{i\mathrm{\Phi }_2})(r_3r_1r_5e^{i\mathrm{\Phi }_1})+(t_3t_4t_5e^{i\mathrm{\Phi }_4})(t_3t_1t_2e^{i\mathrm{\Phi }_3})}\right|.$$
(1)
Here $`\mathrm{\Phi }_j`$ are the phases accumulated along the closed loops (see Fig. 1). Within the RG approach to the conductance distribution, information about electron energy is incorporated only into the values of $`t_i`$. The energy dependence of phases, $`\mathrm{\Phi }_j`$, is irrelevant; they are assumed completely random. Due to this randomness, the transmission coefficients, $`t_i`$, for a given energy, are also randomly distributed with a distribution function $`P(t)`$. Then the transformation (1) allows, upon averaging over $`\mathrm{\Phi }_j`$, to generate the next-step distribution $`P(t^{})`$. Therefore, within the RG scheme, a delocalized state corresponds to the fixed point, $`P_\text{c}(t)`$, of the RG transformation. Due to the symmetry of the RG unit, it is obvious that the critical distribution, $`P_\mathrm{c}(t^2)`$, of the power transmission coefficient, $`t^2=G`$, which has the meaning of the two-terminal conductance, is symmetric with respect to $`t^2=\frac{1}{2}`$. In other words, the RG transformation respects the duality between transmission and reflection. The critical distribution $`P_\mathrm{c}(G)`$ found in Refs. GalR97,  and CaiRSR01,  agrees very well with the results of direct large-scale simulations.

### II.2 RG approach to the LSD

Universal features of the energy level statistics in a macroscopic fully coherent sample at the quantum Hall transition complement the universality in the conductance distribution. The prime characteristics of the level statistics is the LSD – the distribution of the spacings between neighboring energy levels. In order to adjust the RG approach to the calculation of the LSD, it is necessary to “close” the sample at each RG step in order to discretize the energy levels. One of the possible variants of such a closing is shown in Fig. 1 with dashed lines.

For a given closed RG unit with a fixed set of $`t_i`$-values at the nodes, the positions of the energy levels are determined by the energy dependences, $`\mathrm{\Phi }_j(E)`$, of the four phases along the loops. These phases change by $`\pi `$ within a very narrow energy interval, inversely proportional to the sample size. Within this interval the change of the transmission coefficients is negligibly small. A closed RG unit in Fig. 1 contains $`10`$ links, and, thus, it is described by $`10`$ amplitudes. These amplitudes are related by $`10`$ equations ($`2`$ at each node). Each link is characterized by an individual phase. On the other hand, it is obvious that the energy levels are determined only by the phases along the loops. One possible way to derive the system, in which individual phases combine into $`\mathrm{\Phi }_j`$ is to exclude from the original system of $`10`$ equations all amplitudes except the “boundary” amplitudes $`\mathrm{\Psi }_j`$ (see Fig. 1). This procedure is similar to the derivation of Eq. (1). The system of equations for the remaining four amplitudes takes the form

$$\left(\begin{array}{cccc}(r_1r_2t_1t_2t_3)e^{i\mathrm{\Phi }_1}& (t_1r_2+t_2t_3r_1)e^{i\mathrm{\Phi }_1}& t_2t_5r_3e^{i\mathrm{\Phi }_1}& t_2r_3r_5e^{i\mathrm{\Phi }_1}\\ t_1r_3r_4e^{i\mathrm{\Phi }_2}& r_1r_3r_4e^{i\mathrm{\Phi }_2}& (t_4r_5+t_3t_5r_4)e^{i\mathrm{\Phi }_2}& (t_4t_5t_3r_4r_5)e^{i\mathrm{\Phi }_2}\\ t_1t_4r_3e^{i\mathrm{\Phi }_4}& t_4r_1r_3e^{i\mathrm{\Phi }_4}& (r_4r_5t_3t_4t_5)e^{i\mathrm{\Phi }_4}& (t_5r_4+t_3t_4r_5)e^{i\mathrm{\Phi }_4}\\ (t_2r_1+t_1t_3r_2)e^{i\mathrm{\Phi }_3}& (t_1t_2t_3r_1r_2)e^{i\mathrm{\Phi }_3}& t_5r_2r_3e^{i\mathrm{\Phi }_3}& r_2r_3r_5e^{i\mathrm{\Phi }_3}\end{array}\right)\left(\begin{array}{c}\mathrm{\Psi }_1\\ \mathrm{\Psi }_2\\ \mathrm{\Psi }_3\\ \mathrm{\Psi }_4\end{array}\right)=e^{i\omega }\left(\begin{array}{c}\mathrm{\Psi }_1\\ \mathrm{\Psi }_2\\ \mathrm{\Psi }_3\\ \mathrm{\Psi }_4\end{array}\right),$$
(2)
where the parameter $`\omega `$ should be set zero. Then the energy levels, $`E_k`$, of the closed RG unit are the energies for which, with phases $`\mathrm{\Phi }_j(E)=\mathrm{\Phi }_j(E_k)`$, one of the four eigenvalues of the matrix in the left-hand side of Eq. (2) is equal to one. If we keep $`\omega `$ in the right-hand side of Eq. (2), then the above condition can be reformulated as $`\omega (E_k)=0`$. Thus, the calculation of the energy levels reduces to a diagonalization of the $`4\times 4`$ matrix.

The crucial step now is the choice of the dependence $`\mathrm{\Phi }_j(E)`$. If each loop in Fig. 1 is viewed as a closed equipotential as it is the case for the first step of the RG procedure,ChaC88  then $`\mathrm{\Phi }_j(E)`$ is a true magnetic phase, which changes linearly with energy with a slope governed by the actual potential profile, which, in turn, determines the drift velocity. Thus we have

$$\mathrm{\Phi }_j(E)=\mathrm{\Phi }_{0,j}+2\pi \frac{E}{s_j},$$
(3)
where a random part, $`\mathrm{\Phi }_{0,j}`$, is uniformly distributed within $`[0,2\pi ]`$, and $`2\pi /s_j`$ is a random slope. Strictly speaking, the dependence (3) applies only for the first RG step. At each following step, $`n>1`$, $`\mathrm{\Phi }_j(E)`$ is a complicated function of $`E`$ which carries information about all energy scales at previous steps. However, in the spirit of the RG approach, we assume that $`\mathrm{\Phi }_j(E)`$ can still be linearized within a relevant energy interval. The conventional RG approach suggests that different scales in the real space can be decoupled. Linearization of Eq. (3) implies a similar decoupling in the energy space. In the case of phases, a “justification” of such a decoupling is that at each following RG step, the relevant energy scale, that is the mean level spacing, reduces by a factor of $`4`$.

With $`\mathrm{\Phi }_j(E)`$ given by Eq. (3) and fixed values of $`t_i`$, the statistics of energy levels determined by the matrix equation (2) is obtained by averaging over the random initial phases $`\mathrm{\Phi }_{0,j}`$. In particular, each realization of $`\mathrm{\Phi }_{0,j}`$ yields 3 level spacings which are then used to construct a smooth LSD. We now outline the RG procedure for the LSD. The slopes $`s_j`$ in Eq. (3) determine the level spacings at the first step. They are randomly distributed with a distribution function $`P_0(s)`$. Diagonalization of the matrix in Eq. (2) with subsequent averaging over realizations yields the LSD, $`P_1(s)`$, at the second step. Then the key element of the RG procedure, as applied to the level statistics, is using $`P_1(s)`$ as a distribution of slopes in Eq. (3). This leads to the next-step LSD and so on.

It is instructive to compare our procedure of calculating the energy levels with an approach adopted in large-scale simulations within the CC model.Met98b ; KleM97  This approach is based on the unitary network operator $`U`$.KleM97  For a single RG unit this operator acts analogously to the matrix in the left-hand side of Eq. (2). However, within the approach of Refs. Met98b,  and KleM97, , the energy dependence of phases $`\mathrm{\Phi }_j`$ in the elements of the matrix was neglected (only the random contributions, $`\mathrm{\Phi }_{0,j}`$, were kept). Then, instead of the energy levels, $`E_k`$, diagonalization of the matrix (2) yielded a set of eigenvalues, $`\mathrm{exp}(i\omega _k)`$. The numbers $`\omega _k`$ were named quasienergies, and it is the statistics of these quasienergies that was studied in Ref. KleM97, . Comparison of the two procedures for a single RG unit is illustrated in Fig. 2. Fig. 2 shows the dependence of the $`4`$ quasienergies $`\omega _k`$ on the energy $`E`$ calculated for two single sample RG units, with $`t_i`$ chosen from the critical distribution $`P_\text{c}(t)`$. The energy dependence of the phases $`\mathrm{\Phi }_j`$ was chosen from LSD of the GUE according to Eq. (3). It is seen that the dependences $`\omega (E)`$ range from remarkably linear and almost parallel (Fig. 2a) to strongly nonlinear (Fig. 2b).

## III Numerical results

### III.1 The LSD at the QH transition

As a first step of the RG procedure for the calculation of the critical LSD we chose for $`P_0(s)`$ the distribution corresponding to the GUE random matrix ensemble, since previous simulationsBatS96 ; KleM97  indicated that the LSD at the transition is close to GUE. According to $`P_0(s)`$, we pick $`s_j`$ and set $`\mathrm{\Phi }_j`$, $`j=1,\mathrm{},4`$ as in Eq. (3). For the transmission coefficients $`t_i`$, $`i=1,\mathrm{},5`$ we use the fixed-point distribution $`P_\text{c}(t)`$,comment-1  obtained previously.CaiRSR01

From the solutions of Eq. (2) corresponding to $`\omega _j(E_k)=0`$ the new LSD $`P_1(s^{})`$ is constructed using the “unfolded” energy level spacings $`s_m^{}=(E_{m+1}E_m)/\mathrm{\Delta }`$, where $`m=1,2,3`$, $`E_{k+1}>E_k`$ and the mean spacing $`\mathrm{\Delta }=(E_4E_1)/3`$. Due to the “unfolding”Haa92  with $`\mathrm{\Delta }`$, the average spacing is set to one for each sample and in each RG-iteration step we superimpose spacing data of $`2\times 10^6`$ RG units. The resulting LSD is discretized in bins with largest width $`0.01`$. In the following iteration step we repeat the procedure using $`P_1`$ as initial distribution. We assume that the iteration process has converged when the mean-square deviation of distribution $`P_n(s)`$ deviates by less than $`10^4`$ from its predecessor $`P_{n1}(s)`$. The RG iteration process converges rather quickly after only $`23`$ RG steps. The resulting LSD, $`P_\text{c}(s)`$, is shown in Fig. 3 together with an LSD for the unitary random matrix ensemble.

Although $`P_\text{c}(s)`$ exhibits the expected features, namely, level repulsion for small $`s`$ and a long tail at large $`s`$, the overall shape of $`P_\text{c}(s)`$ differs noticeably from GUE. In the previous large-size lattice simulations BatS96 ; KleM97  the obtained critical LSD was much closer to GUE than $`P_\text{c}(s)`$ in Fig. 3. This fact, however, does not reflect on the accuracy of the RG approach. Indeed, as it was demonstrated recently, the critical LSD – although being system size independent — nevertheless depends on the geometry of the samplesPotS98  and on the specific choice of boundary conditions.BraMP98 ; SchP98  Sensitivity to the boundary conditions does not affect the asymptotics of the critical distribution, but rather manifests itself in the shape of the “body” of the LSD. Recall now that the boundary conditions which we have imposed to calculate the energy levels (dashed lines in Fig. 1) are non-periodic.

There is another possibility to assess the critical LSD, namely by iterating the distribution of quasienergies. In Fig. 3 we show the result of this procedure. It appears that the resulting distribution is almost identical to $`P_\text{c}(s)`$. This observation is highly non-trivial, since, as follows from Fig. 2, there is no simple relation between the energies and quasienergies. Moreover, if instead of the linear $`E`$-dependence of $`\mathrm{\Phi }_j`$, we choose another functional form, say,

$$\mathrm{\Phi }_j(E)=\mathrm{\Phi }_{0,j}+2\mathrm{arcsin}\left(\frac{E}{s_j}2p\right),$$
(4)
where the integer $`p`$ insures that $`\left|\frac{E}{s_j}2p\right|1`$, then, the RG procedure would yield an LSD which is markedly different (within the “body”) from $`P_\text{c}(s)`$. This is illustrated in Fig. 4.

Both procedures, using quasienergies instead of real energies \[as in Ref. KleM97, \], and linearization of the energy dependence of phases \[as in Eq. (3)\] are not rigorous. Linearization is dictated by the RG concept. The coincidence of the results of the two procedures indicates that the concept of quasienergies, namely, that they obey the same statistics as real energies, is equivalent to the RG.

### III.2 Small and large $`s`$ behavior

As it was mentioned above, the general shape of the critical LSD is not universal. However, the small $`s`$ behavior of $`P_\text{c}(s)`$ must be the same as for the unitary random matrix ensemble, namely $`P_\text{c}(s)s^2`$. This is because delocalization at the quantum Hall transition implies the level repulsion.ShkSSL93 ; FyoM97  Earlier large-scale simulations of the critical LSD KawOSO96 ; BatSZK96 ; Met98b ; FeiAB95 ; OhtO95 ; BatS96 ; MetV98 ; BatSK98 ; Met99 ; OnoOK96 ; KleM97  satisfy this general requirement. The same holds also for our result, as can be seen in Fig. 5. The given error bars of our numerical data are standard deviations computed from a statistical average of $`100`$ FP distributions each obtained for different random sets of $`t_i`$’s and $`\mathrm{\Phi }_j`$’s within the RG unit. In general, within the RG approach, the $`s^2`$-asymptotics of $`P(s)`$ is most natural. This is because the levels are found from diagonalization of the $`4\times 4`$ unitary matrix with absolute values of elements widely distributed between $`0`$ and $`1`$.

The right form of the large-$`s`$ tail of $`P(s)`$ is Poissonian, $`P_\text{c}(s)\mathrm{exp}(bs)`$.ShkSSL93  For the Anderson model in three dimensions, unambiguous confirmation of this prediction in numerical simulations became possible only when very high numerical accuracy had been achieved.BatSZK96 ; ZhaK97  This is because $`P_\text{c}(s)`$ assumes the Poissonian asymtotics only at large enough $`s3\mathrm{\Delta }`$. For the quantum Hall transition, a linear behavior of $`\mathrm{ln}P_\text{c}(s)`$ with a slope corresponding to the value $`b4.1`$ has been found in Ref. BatS96,  from the analysis of the interval $`2<s/\mathrm{\Delta }<4`$. Our data, as shown in Fig. 6, has a high accuracy only for $`s/\mathrm{\Delta }2.5`$. For such $`s`$, the distribution $`P_\text{c}(s)`$ does not yet reach its large-$`s`$ tail. Thus, the value of parameter $`b`$ extracted from this limited interval is somewhat ambiguous. Namely, we obtain $`b=5.442`$ for $`s/\mathrm{\Delta }[1.5,2.0]`$ and $`b=6.803`$ for $`s/\mathrm{\Delta }[2.0,2.5]`$.

Summarizing, the accuracy of the RG approach, applied to the level statistics, is insufficient to discern the only non-trivial feature of the critical LSD, i.e., the universal Poissonian asymtotics. However, the scaling analysis of LSD clearly reveals the universal features of the quantum Hall transition as we demonstrate in the next Section.

## IV Scaling results for the LSD

### IV.1 Finite-size scaling at the QH transition

The critical exponent, $`\nu `$, of the quantum Hall transition governs the divergence of the correlation length $`\xi _{\mathrm{}}`$ as a function of the arbitrary control parameter $`z_0`$, i.e.,

$$\xi _{\mathrm{}}(z_0)|z_0z_\text{c}|^\nu ,$$
(5)
where $`z_\text{c}`$ is the critical value. The values of $`\nu `$ calculated using different numerical methods, e.g., $`\nu =2.35\pm 0.03`$,Huc92  $`2.4\pm 0.2`$,LeeWK93  $`2.5\pm 0.5`$ChaC88  agree with each other. The RG approach for the conductance distribution also yields a rather accurate value $`\nu =2.39\pm 0.01`$.CaiRSR01  In Sec. II we have introduced a complimentary RG approach to the distribution of the energy levels at the transition. It can be expected on general grounds, that the LSD obtained from the RG approach obeys scaling at small enough $`z_0z_\text{c}`$. However, it is by now means obvious, whether the value of $`\nu `$ extracted from different variants of the RG approach are consistent.

In order to extract $`\nu `$ from the LSD we employ the one-parameter-scaling analysis. This analysis is based on the rescaling of a quantity $`\alpha (N;\{z_i\})`$ — depending on (external) system parameters $`\{z_i\}`$ and the system size $`N`$ — onto a single curve by using a scaling function $`f`$

$$\alpha (N;\{z_i\})=f\left(\frac{N}{\xi _{\mathrm{}}(\{z_i\})}\right).$$
(6)
Since Eq. (5), as indicated by “$`\mathrm{}`$”, holds only in the limit of infinite system size, we now use the scaling assumption to extrapolate $`f`$ to $`N\mathrm{}`$ from the finite-size results of the computations. Once $`f`$ and $`\xi _{\mathrm{}}`$ are known, the value of $`\nu `$ can be then inferred.

In the original formulation of the RG approachGalR97  it was demonstrated that there is a natural parametrization of the transmission coefficients $`t`$, i.e., $`t=(e^z+1)^{1/2}`$. For such a parametrization, $`z`$ can be identified with a dimensionless electron energy. The quantum Hall transition occurs at $`z=0`$, which corresponds to the center of the Landau band. The universal conductance distribution at the transition, $`P_\text{c}(G)`$, corresponds to the distribution $`Q_\text{c}(z)=P_\text{c}\left[(e^z+1)^1\right]/4\mathrm{cosh}^2(z/2)`$ of parameter $`z`$, which is symmetric with respect to $`z=0`$ and has a shape close to a gaussian.CaiRSR01  The RG procedure for the conductance distribution converges and yields $`Q_\text{c}(z)`$ only if the initial distribution is an even function of $`z`$. This suggests to choose as a control parameter in Eq. (6), $`z_0`$, the position of the maximum of the function $`Q(z)`$. Then the meaning of $`z_0`$ is the electron energy measured from the center of Landau band. The fact that the quantum Hall transition is infinitely sharp implies that for any $`z_00`$, the RG procedure drives the initial distribution $`Q(zz_0)`$ towards an insulator, either with complete transmission of the network nodes (for $`z_0>0`$) or with complete reflection of the nodes (for $`z_0<0`$).

### IV.2 Scaling for $`\alpha _\text{P}`$ and $`\alpha _\text{I}`$

In principle, we are free to choose for the finite-size analysis any characteristic quantity $`\alpha (N;z_0)`$ constructed from the LSD which has a systematic dependence on system size $`N`$ for $`z_00`$ while being constant at the transition $`z_0=0`$. Because of the large number of possible choicesShkSSL93 ; HofS94b ; ZhaK97 ; BatS96 ; ZhaK95b ; ZhaK95c  we restrict ourselves to two quantities which are obtained by integration of the LSD and have already been successfully used in Refs. HofS94b,  and HofS93, , namely

$$\alpha _\text{P}=_0^{s_0}P(s)𝑑s$$
(7)
and second

$$\alpha _\text{I}=\frac{1}{s_0}_0^{s_0}I(s)𝑑s,$$
(8)
with $`I(s)=_0^sP(s^{})𝑑s^{}`$. The integration limit is chosen as $`s_0=1.4`$ which approximates the common crossing pointHofS94b  of all LSD curves as can be seen in Fig. 7.

Thus $`P(s_0)`$ is independent of the distance $`|z_0z_\text{c}|`$ to the critical point and the system size $`N`$. We note that $`N`$ is directly related to the RG step $`n`$ by $`N=2^n`$. The double integration in $`\alpha _\text{I}`$ is numerically advantageous since fluctuations in $`P(s)`$ are smoothed. We now apply the finite-size-scaling approach from Eq. (6)

$$\alpha _{\text{I,P}}(N,z_0)=f\left(\frac{N}{\xi _{\mathrm{}}(z_0)}\right).$$
(9)
Since $`\alpha _{\text{I,P}}(N,z_0)`$ is analytical for finite $`N`$, one can expand the scaling function $`f`$ at the critical point. The first order approximation yields

$$\alpha (N,z_0)\alpha (N,z_\text{c})+a|z_0z_\text{c}|N^{1/\nu }$$
(10)
where $`a`$ is a dimensionless coefficient. For our calculation we use a higher order expansion proposed by Slevin and Ohtsuki.SleO99a  In Ref. SleO99a,  the function $`f`$ is expanded twice, first, in terms of the Chebyshev polynomials of order $`𝒪_\nu `$ and, second, in Taylor series with terms $`|z_0z_c|`$ in the power $`𝒪_z`$. This procedure allows to describe the deviations from linearity in $`|z_0z_c|`$ at the transition. In addition, in Ref. SleO99a,  the contributions from an irrelevant scaling variable which leads to a shift of the transition for small system sizes was taken into account. In our case, in contrast to the Anderson model of localization, the transition point $`z_0=0`$ is known.

Therefore, we can neglect the influence of irrelevant variables. In order to obtain the functional form of $`f`$, the fitting parameters, including $`\nu `$, are evaluated by a nonlinear least-square ($`\chi ^2`$) minimization. In Fig. 8 we show the resulting fit for $`\alpha _\text{P}`$ and $`\alpha _\text{I}`$ at the transition.

The fits are chosen in such a way that the total number of parameters is kept at a minimal value, while the fit agrees well with the numerical data.NdaRS02  The corresponding scaling curves are displayed in Fig. 9. In the plots the two branches corresponding to complete reflection ($`z_0<0`$) and complete transmission ($`z_0>0`$) can be clearly distinguished. In order to estimate the error of fitting procedure we compare the results for $`\nu `$ obtained by different orders $`𝒪_\nu `$ and $`𝒪_z`$ of the expansion, system sizes $`N`$, and regions around the transition. A part of our over $`100`$ fit results together with the standard deviation of the fit are given in Table 1.

The value of $`\nu `$ is calculated as the average of all individual fits where the resulting error of $`\nu `$ was smaller than $`0.02`$. The error is then determined as the standard deviation of the contributing values. By this method we assure that our result is not influenced by local minima of the nonlinear fit. So we consider $`\nu =2.37\pm 0.02`$ as a reliable value for the exponent of the localization length at the QH transition obtained from the RG approach to LSD. This is in excellent agreement with $`\nu =2.35\pm 0.03`$ (Ref. Huc92, ), $`2.4\pm 0.2`$ (Ref. LeeWK93, ), $`2.5\pm 0.5`$ (Ref. ChaC88, ), and $`2.39\pm 0.01`$ (Ref. CaiRSR01, ) calculated previously. In addition to $`\alpha _\text{P}`$ and $`\alpha _\text{I}`$, we tested also a parameter-free scaling quantity $`_0^{\mathrm{}}s^2P(s)𝑑s`$,ZhaK95b  where the whole distribution $`P(s)`$ is taken into account. Here, due to the influence of the large $`s`$-tail a less reliable value $`\nu =2.33\pm 0.05`$ was obtained.

### IV.3 Test of consistency

Finally we address the question, how the actual form of the distribution $`Q(z)`$ affects the results for LSD and the scaling analysis. Recall that in the above calculations we have used at each step of the RG procedure the distribution $`Q(z)`$ derived from the critical conductance distribution, $`P_\text{c}(G)`$. The function $`P_\text{c}(G)`$ is shown in Fig. 10 (inset) with a full line. In order to understand the importance of the fact that $`P_\text{c}(G)`$ is almost flat, we have repeated our calculations choosing for $`P(G)`$ a relatively narrow gaussian distribution $`P(G)P_{\text{Gauß}}(G)`$ at each RG step. This distribution is shown with a dashed line in Fig. 10 (inset). The obtained LSD is presented in Fig. 10. Obviously, it agrees much worse with GUE, which can be considered as a reference point, than the LSD computed using the true $`P_\text{c}(G)`$. Our data for $`\alpha _\text{I}`$ calculated for $`P(G)=P_{\text{Gauß}}(G)`$ is plotted in Fig. 11.

The curves for small system sizes $`N`$ exhibit strong deviations, i.e., there is initially no common crossing point, while for large $`N`$ a behavior similar to Fig. 8 is observed. Therefore, small $`N`$ data are neglected in the scaling analysis. The $`\chi ^2`$ fits for $`\alpha _\text{I}`$ and $`\alpha _\text{P}`$ are carried out using $`z_0[0.05,0.05]`$ and $`N=16512`$. They yield the values $`\nu _\text{I}=2.43\pm 0.02`$ and $`\nu _\text{P}=2.46\pm 0.03`$, which are also less accurate than $`\nu `$ calculated with the critical $`P_\text{c}(G)`$. Overall, Figs. 10 and 11 illustrate the consistency of the RG approaches for the conduction distribution and for the level statistics, in the sense, that the best fixed point distribution of the level spacings corresponds to the fixed point of the conductance distribution.

## V Conclusion

Network models introduced in Ref. Sha82,  turned out to be a powerful tool to study the Anderson localization. Without magnetic field, the propagation of electron waves along each link of the network is allowed in both directions. In two dimensions the transmission coefficient of the network is zero for all parameters of the scattering matrix at the nodes,FreJM99  illustrating complete localization of electronic states. On the other hand, the two-channel network model with inter-channel mixing, that models spin-orbit interaction, exhibits a localization-delocalization transitionMerJH98  that is also in accord with the scaling theory of localization.HikLN80  However, the version of the network model that has been most widely studied, is the chiral version, i.e., the CC model,ChaC88  describing the electron motion in a disordered system in a strong magnetic field limit. Within the CC model, the scattering matrix at the node is parametrized by a single number, e.g., the transmission coefficient $`t`$. On the qualitative level, the CC model yields a transparent explanation why delocalization occurs only at a single energy, for which $`t^2=1/2`$. On the quantitative level, in addition to the exponent, $`\nu `$, more delicate characteristics of the critical wave functions were extracted from the numerical analysis of the CC model.HucK99 ; KleZ01

The fact that the RG approach, within which the correlations between different scales are neglected, describes the results of the large-scale simulations of the CC model so accurately, indicates that only a few spatial correlations within each scale are responsible for the critical characteristics of the quantum Hall transition. More precisely, the structure of the eigenstates of a macroscopic sample at the transition can be predicted from the analysis of a single RG unit consisting of only five nodes. Earlier we have demonstrated this fact for the conductance distribution.CaiRSR01  In the present paper this statement is reinforced by the study of the level statistics at the transition, which is a complimentary (to the conductance distribution) characteristics of the localization.

###### Acknowledgements.

We thank R. Klesse, L. Schweitzer and I. Zharekeshev for stimulating discussions. This work was supported by the NSF-DAAD collaborative research Grant No. INT-0003710. P.C. and R.A.R. also gratefully acknowledge the support of DFG within the Schwerpunktprogramm “Quanten-Hall-Systeme” and the SFB 393.


text:
# HU-EP-02/39LU-ITP 2002/020Edinburgh 2002/10 LTH 559 Finite Size Effects in Nucleon Masses in Dynamical QCD presented by A. Ali Khan

## 1 INTRODUCTION

For a lattice calculation of physical quantities such as hadron masses it is desirable to have a theoretical understanding of finite size effects and to be able to extrapolate to infinite volume. On lattices used in typical quenched calculations, finite size effects of hadron masses seem to be small, but they may be considerable in the dynamical case . A source of finite size effects is a virtual pion going around the boundary of the lattice before it is absorbed again. At sufficiently small pion masses and large volumes this effect can be described by chiral perturbation theory, the low-energy effective theory of nucleons and pions, in a finite box. Here, we discuss finite size effects on nucleon masses generated by the QCDSF and UKQCD collaborations.

## 2 THE SIMULATION

The simulations were done using the plaquette gauge action and $`N_f=2`$ dynamical non-perturbatively $`O(a)`$-improved Wilson fermions. We consider lattices of size $`1`$ fm, $`1.6`$ fm, and $`2`$ fm. Pion masses are in the range of about $`0.61`$ GeV. The lattice spacing is $`a0.1`$ fm. We use $`r_0=0.5`$ fm to set the scale. Valence and sea quark masses are taken to be equal.

The lattice data for the nucleon mass are plotted in Fig. 1. We see significant finite size effects. We expect $`O(a^2)`$ effects to be small. Indeed, our two mass points at $`r_0m_{PS}2`$ on the 1.6 fm lattice differ by less than a percent, whereas $`a^2`$ differs by $`O(20)\%`$.

Chiral perturbation theory in the infinite volume predicts

$`m_Nr_0`$ $`=`$ $`a+b(m_{PS}r_0)^2+c(m_{PS}r_0)^3,`$ (1)
where the cubic term reflects the leading non-analytic (LNA) behaviour of the nucleon self-energy. The chiral extrapolation is performed separately for different lattice sizes using Eq. (1). At 1.6 fm, the chiral extrapolation gives a value that is 8% or $`2\sigma `$ higher than experiment, at 2 fm no discrepancy is found. The coefficient $`c`$ disagrees significantly with the infinite volume LNA prediction.

## 3 THEORETICAL PREDICTIONS

The finite size effects are calculated from the difference of the nucleon self-energy $`\mathrm{\Sigma }`$ in finite and infinite spatial volume,

$$\delta m_N=\mathrm{\Sigma }(L_s)\mathrm{\Sigma }(\mathrm{}).$$
(2)
In the perturbative calculation, the $`x_4`$ direction is taken to be of infinite extent. We calculate the self-energy at one loop in Heavy Baryon $`\chi PT`$ according to the lattice prescription given in Ref.  adapted to the case of 2-flavours. Additionally, we include $`\mathrm{\Delta }`$ intermediate states. For this we discretise the $`N\mathrm{\Delta }\pi `$ interaction Lagrangian given in Ref. . The numerical values for $`f_\pi `$, $`g_A=𝒟+`$, and the $`N\mathrm{\Delta }\pi `$ coupling $`c_A`$ are also taken from this reference. The relevant Feynman diagrams are given in Fig. 2.

To regularize the lattice integrals we vary the cut-off $`\mathrm{\Lambda }`$ between $`\pi /a=6`$ GeV and $`\mathrm{}`$. The difference between using $`\mathrm{\Lambda }=6`$ GeV and the continuum limit is around $`10\%`$ for $`m_{PS}600`$ MeV. The values presented here are for $`\mathrm{\Lambda }30`$ GeV, which practically amounts to infinite cut-off. Results for the smallest pion mass used in the simulation at $`1.6`$ fm are given in Fig. 3. For $`L_s1.6`$ fm, the $`L_s`$ dependence can be approximated by an exponential decay,

$$\delta m_N=(c_0/L_s)\mathrm{exp}(c_1L_s),$$
(3)
where the fit with Eq. (3) describes the $`\chi PT`$ result with an accuracy of several percent.

The finite size effect of the nucleon mass has been estimated previously, relating the mass shift to the $`\pi N`$ scattering amplitude . A comparison with our one-loop result is given in Fig. 3.

In Fig. 4, we show $`m_N\delta m_N`$. We see that the masses on the smaller volumes are brought down towards a universal curve. In particular, 50% of the finite size corrections at $`L_s=1`$ fm, and 60% at $`L_s=1.6`$ fm, are accounted for by chiral perturbation theory. We consider this a remarkable result, considering that the pion mass is relatively heavy.

## 4 CONCLUSIONS AND DISCUSSION

We find indications for finite size effects in the nucleon mass on lattices of $`1.62`$ fm size. We calculated the finite volume corrections in lattice-regularised $`\chi PT`$ at one loop. At $`L_s1.6`$ fm, the correction is of the order of a few percent for the pion masses considered. $`\chi PT`$ describes about $`5060\%`$ of the finite size effects of the data. The mass dependence is exponential rather than power-like even for intermediate lattices and pion masses. Higher order corrections to the one-loop calculation may be important at present pion masses.

## ACKNOWLEDGEMENTS

We acknowledge useful discussions with Th. Hemmert. The numerical calculations were performed on the Hitachi SR8000 at LRZ (Munich), the APEmille at NIC (Zeuthen), the Cray T3E at EPCC (Edinburgh) and NIC (Jülich). We wish to thank all institutions for their support. This work has been supported in part by the European Community’s Human Potential Program under contract HPRN-CT-2000-00145, Hadrons/Lattice QCD.


text:
# Disks on a Sphere and two-dimensional Glasses

## Abstract

Talk given at the conference on ‘Inhomogeneous Random Systems’ at the University of Cergy-Pontoise, France, (23 January 2001).
I describe the classic circle-packing problem on a sphere, and the analytic and numerical approaches that have been used to study it. I then present a very simple Markov-chain Monte Carlo algorithm, which succeeds in finding the best solutions known today. The behavior of the algorithm is put into the context of the statistical physics of glasses.
keywords: Monte Carlo Methods; Packing of Circles; Glass Transition
AMS numbers: 00A72; 52C15; 82D30

What is the minimal radius $`R_N`$ of a sphere allowing $`N`$ non-overlapping circles of radius $`r=\frac{1}{2}`$ to be drawn on its surface? In a classic paper, Schütte and van der Waerden  reduced this problem to a study of planar graphs: As the sphere has minimal radius, many circles are blocked, i.e. they cannot move. These circles make up the $`MN`$ vertices of the graph. If two circles are in (blocking) contact, they are connected by an edge.

Schütte and van der Waerden considered all graphs with $`M12`$ vertices. Some of these graphs can correspond to a blocked configuration of circles on a sphere of radius $`\stackrel{~}{R}`$, with $`R=\mathrm{min}_{\text{graphs}}\stackrel{~}{R}`$. They were able to determine $`R_N`$ for $`N12`$. For this latter case $`N=12`$, e.g., the centers of the circles are located at the vertices of an icosahedron of side length $`1`$, and lie on a sphere of radius

$$R_{12}=\frac{\sqrt{10+2\sqrt{5}}}{4}=0.951\mathrm{}$$
(1)
Schütte and van der Waerden  were later able to prove that $`R>1`$ for $`N>12`$. This result (cf figure 1) implies that $`13`$ spheres of radius $`1`$ cannot be placed on top of another sphere of the same size, and solves a packing problem discussed at least since the time of Newton. The best currently known packing of spheres for $`N=13`$ concerns spheres of radius $`1`$ on top of another sphere of radius $`1.0911\mathrm{}`$. The corresponding circle arrangement is shown in figure 2, and is almost certainly optimal. Today, it should be easy to enumerate planar graphs on a computer, and to prove this conjecture.

The problem of packing disks on the surface of a sphere has appeared repeatedly in the literature, as there are important applications in mathematics and the natural sciences. Hundreds of papers have been written on this subject (for a partial list of references cf e.g. ), most of them improving upper bounds on $`R_N`$ in the range of $`13<N100`$. Only for $`N=24`$ has a new optimal solution been proven .

The problem of $`N`$ non-overlapping disks on a minimal sphere can be reformulated as a minimization problem of the potential energy of $`N`$ particles $`i=1,\mathrm{},N`$ on the surface of the sphere of radius $`1`$ with a two-body interaction

$$V=\underset{i=1}{\overset{N1}{}}\underset{j=i+1}{\overset{N}{}}\left(\frac{c}{r_{ij}}\right)^p$$
(2)
where $`r_{ij}`$ is the Euclidean distance between circles $`i`$ and $`j`$. Therefore, customary minimization programs (such as the Newton-Raphson method) can be applied to find local - and hopefully global - minima for ever increasing values of $`p`$. In our first units (circles of radius $`\frac{1}{2}`$ on a sphere of radius $`R`$), the minimal radius $`R_N`$ is given by

$$R_N=\frac{1}{\mathrm{min}_{ij}r_{ij}}$$
(3)
in the limit $`p\mathrm{}`$.

Current records were set by Kottwitz , who used values of $`p`$ in the astonishing range of $`p=80`$ to $`p=\mathrm{1\hspace{0.17em}310\hspace{0.17em}720}`$. Running such a sophisticated Newton-Raphson program is a delicate task as one has to avoid getting trapped in local minima of the potential energy. Nevertheless, the strategy proved to be successful, as a large number of new values for $`R_N`$ were found, and as these values have so far stood the test of time.

Recently, I have reconsidered the non-overlapping circle problem from the view point of Markov-chain Monte Carlo algorithms. Instead of searching for the minimum of a regularized energy, I move circles during a large number of time steps $`t=1,\mathrm{},T_{\mathrm{max}}`$. At time $`t`$, a single, randomly chosen circle $`i`$ is moved a tiny bit into a random direction. Specifically, let us denote the position of $`i`$ at time $`t`$ by $`\stackrel{}{r}_i`$, and let the displacement vector be $`\stackrel{}{\delta r}=(\delta x,\delta y,\delta z)`$ with Gaussian random numbers $`\delta x`$ etc. that have zero mean and variance $`\sigma ^21`$. Then the following move is considered

$$\stackrel{}{r}_i\stackrel{}{r}_i^{}=\frac{\stackrel{}{r}_i+\stackrel{}{\delta r}}{|\stackrel{}{r}_i+\stackrel{}{\delta r}|}R$$
(4)
If the move can be performed without violating the constraints, it is accepted, and rejected otherwise. This procedure implements a standard local Metropolis algorithm for a constant stationary probability distribution: Any non-overlapping configuration of circles should be visited with the same probability.

In addition to performing the Metropolis algorithm at fixed $`R`$, I also implemented the “simulated annealing” procedure : Every so often, the radius of the sphere is reduced by a very small amount, if this is possible without introducing overlaps between particles .

I have found that this simple algorithm reproduces the best results known up to date (cf ) for $`N=5\mathrm{}90`$, typically within a few minutes of computer time .

This observation may be of practical interest, as the Monte Carlo algorithm is much simpler than the previous approaches. However, it also underlines the quality of the bounds obtained by Kottwitz, which I have been able to reach, but not to improve on.

In physical terms, the graph giving $`R_N`$ is the ground state of the system (which does not have to be unique: for $`N=15`$, e.g., two best solutions are known ). Following Stillinger and Weber , all the other graphs that can be put on the surface of a sphere (with larger radii $`\stackrel{~}{R}>R_N`$) may be called ‘inherent structures’. These are the configurations that the simulated annealing algorithm may get trapped in. If the Monte Carlo algorithm is trapped in an inherent structure at radius $`\stackrel{~}{R}`$, it could by continuity also get trapped by the same inherent structure at a slightly larger radius $`R^{}\stackrel{>}{}\stackrel{~}{R}`$. This means that, in the strict sense, the Monte Carlo algorithm is non-ergodic at least for all $`R^{}`$ satisfying

$$\underset{\text{graphs}}{\mathrm{max}}\stackrel{~}{R}>R^{}>R_N.$$
(5)
In this and similar systems, very little is known rigorously about the number of inherent structures as a function of radius (their density of states), and enumeration methods at small $`N`$ should be very helpful. For physical applications, the thermodynamic limit $`N\mathrm{}`$ is of prime interest. In this limit, $`(R^{}\stackrel{~}{R})/\stackrel{~}{R}`$ certainly goes to zero. This means that strict ergodicity (on timescales that diverge with $`N`$) is reinstalled in the thermodynamic limit. Precise understanding of this complicated thermodynamic limit is lacking.

If, in contrast, the simulated annealing algorithm got trapped in a non-optimal state (with $`R`$ clearly larger than $`R_N`$) with probability close to one, we would say that the ‘physical’ system of disks on the surface of a sphere had a glassy phase. Such a phase seems to be absent for monodisperse circles. Our findings agree very well with what has been found in monodisperse hard disk and hard sphere systems on two- and three-dimensional tori .

In contrast, there is strong numerical evidence that polydisperse hard disks can be glassy: in a system of hard disks of varying radii on a torus , a comparable local Monte Carlo algorithm always gets trapped at finite density. Again, it is fundamentally unclear how this observation carries over into the thermodynamic limit. Empirically, diffusion times of particles seem to grow without bounds at a well-defined, finite density.

In these two-dimensional systems, it has been found  that other Monte Carlo dynamics do not suffer the same slow-down as the local algorithm. These methods also provide strong indications that the glass transition in these two-dimensional models is a purely dynamic phenomenon, which does not manifest itself in equilibrium averages over the stationary probability distribution.

Acknowledgements: I would like to thank L. Santen and S. Tanase-Nicola for helpful discussions.