text:
# New Statistical Methods for Analysis of Large Surveys: Distributions and Correlations

## 1 INTRODUCTION

One of the important ways of testing the models of AGN, or any other astrophysical source, is through the investigations of the distributions, ranges, and more importantly the correlations among, the relevant physical characteristics, such as luminosity, spectrum, redshifts or distances. A reliable determination of these features requires large samples. As evident from papers presented in this proceedings the samples are becoming larger and larger. Combining the samples, however, is a very challenging task, because different samples are obtained by different instruments and techniques, so that they suffer from different and varied selection biases and data truncations. Overcoming these biases requires care.

The primary goal of this paper is to describe some of the relatively new methods we have developed at Stanford over the past decade (Efron & Petrosian 1992, 1999). These methods are very general and are applicable to any data with well defined but arbitrary truncations. We have applied these to various astrophysical data such as solar flares (Lee, Petrosian, & McTiernan 1993, 1995), gamma-ray burst (see e.g. Lloyd, Petrosian, & Mallozzi 2000) and quasars (Meloney & Petrosian 1999). Instead of using abstract mathematical symbols, the method will be demonstrated using the luminosity function of AGNs and its cosmological evolution, i.e. its variation with redshift $`z;\mathrm{\Psi }(L,z)`$. Without loss of generality, we can write the luminosity function as

$$\mathrm{\Psi }(L,z)=\rho (z)\psi (L/g(z),\alpha _i)/g(z),$$
(1)
where $`\rho (z)`$ describes the co-moving density evolution and $`g(z)`$ (with $`g(0)=1`$) describes the luminosity evolution of the population with $`L_o=L/g(z)`$ as the luminosity adjusted to its present epoch value; $`\psi (L_o,\alpha _i)`$ gives the local luminosity function. Here I explicitly include the shape parameters $`\alpha _i`$, which could also vary with redshift. A surprising result has been the absence of evidence for a strong shape evolution. In this paper I ignore such effects and concentrate on the determination of the the density and luminosity evolution functions $`\rho (z)`$ and $`g(z)`$.

In the next section we describe various kind of truncations and in §3 give a brief summary of some of the past methods used for this kind of analysis. The new methods and their accuracy are described in §4 and a sample of results from application to AGNs are summarized in §5.

## 2 TYPES OF TRUNCATIONS

The left panel of Figure 1 shows a set of arbitrary data points labeled as luminosity $`L`$ and redshift $`z`$ and several generic truncations. The distribution may be truncated parallel to the axis as shown by the dotted lines. This only limits the observed ranges of the variables but does not introduce any bias. This kind of data will be referred to as untruncated data. However, this is rarely the case for astronomical data, and in general, one is dealing with cases where the truncation is not parallel to the axis. The simplest and most common case is when the data suffers a one-sided truncation from below as shown by the solid curve. This is the case for magnitude or flux limited samples; $`L4\pi d_L^2(\mathrm{\Omega }_i,z)f_{min}`$, where $`f_{min}`$ is the limiting flux and $`d_L(\mathrm{\Omega }_i,z)`$ is the luminosity distance at $`z`$ for an assumed cosmological model represented by $`\mathrm{\Omega }_i`$. In some cases the data may be truncated from above as shown by the dashed curve. The statistical methods do not distinguish between truncation from above or below. However, the data analysis is affected when there are truncations both from above and below. This is the case for some AGN samples. The situation becomes even more complex when the truncation boundaries are not monotonic, or when one tries to combine samples with different truncations, say different upper and lower flux limits. Such variation may be present even within a given catalog where data taken at different times and directions in the sky may have different limits. The most general truncation is when each data point, say \[$`L_i,z_i`$\], has its individual upper and lower limits, $`L_i^{}<L_i<L_i^+`$ and $`z_i^{}<z_i<z_i^+`$, as shown by the large cross for one point on Figure 1. The methods we have developed can treat this most general truncation.

## 3 A BRIEF HISTORICAL REVIEW

The process of determination of the distribution of physical characteristics from truncated data has a long history starting with first and rudimentary observations of stars in the disk of our galaxy. Here I will touch upon some of the relevant highlights. A more detailed discussion can be found in the references cited below and in a review article (Petrosian 1992). I will limit the discussion to works aimed at determination of the luminosity function and spatial distribution of sources; $`\mathrm{\Psi }(L,r)`$ from a magnitude limited data (one sided truncation) as shown by the right panel in Figure 1. Almost all the methods I will describe ignore possible correlations between the variables and assume that these are independent. In the case of the luminosity function this means that $`g(z)=1`$, so that the bivariate distribution is separable; $`\mathrm{\Psi }(L,r)=\psi (L)\rho (r)`$. Unfortunately this often unjustified simplification is prevalent even today. This unnecessary assumption often can lead to erroneous conclusions. As stressed below the first task must be testing the bivariate distribution for correlations between the variables.

In general, most of the methods can be divided into two categories, parametric or non parametric. In the former one assumes a parametric form for the two functions $`\psi (L)`$ and $`\rho (r)`$. In the non parametric methods one often ends up with a description for the cumulative functions

$$\varphi (L)=_L^{\mathrm{}}\psi (L^{})𝑑L^{},\sigma (r)=_0^r(dV/dr)\rho (r^{})𝑑r^{},$$
(2)
where $`V(r)`$ is the volume of space (included in the observation) from the origin to distance $`r`$. In these relation $`r`$ could be any measure of distance including redshift, look back time etc.

### 3.1 Parametric Methods

The first methods were developed for investigation of the luminosity function and spatial distribution of stars perpendicular to the galactic disk. These and subsequent applications to other sources have come to be known as correction for Malmquist Bias after Malmquist (1922); Eddington (1915, 1940) also describes this method. A more recent description can be found in Trumpler & Weaver (1953). These early works dealing with the distribution of stars assume a Gaussian distribution of absolute magnitudes (i.e. a log-normal luminosity function) and a Gaussian spatial distribution perpendicular to the disk. Because of the truncation there is absence of low luminosity stars at large distance so that the raw observed distributions of $`L`$ and $`r`$ are biased. The essence of the method was to correct for this bias. It turns out that the method used was sound but the final results were incorrect because of the erroneous assumptions about the forms of the distributions. We know today that the luminosity function of stars is best represented by a broken power law and that the fall of the stars perpendicular to the disk is exponential and not Gaussian. This demonstrates the shortcoming of this and other parametric methods. Similar parametric methods have been used for extragalactic sources (galaxies and quasars) first by Neymann & Scott (1959) and in numerous works ever since.

### 3.2 Non parametric Methods

One of the most commonly non parametric methods used is the so-called $`V/V_{max}`$ method first described by Kafka (1966) soon after the discovery of the quasars but used most successfully by Schmidt (1968). Independence is again assumed and the presence or absence of density evolution is tested by a single moment, namely the average value of $`V/V_{max}`$, where $`V_{max}`$ is the volume up to the maximum redshift (or distance) that a source of luminosity $`L`$ can be visible given the limiting flux $`f_{min}`$; $`L=4\pi d_L^2(\mathrm{\Omega }_i,z_{max})`$ (see Figure 1). In the absence of evolution one expects a value of 0.5 for this average. A more general method was described later by Avni & Bahcall (1980). Of course, one need not be limited only to one moment of the distribution. More information can be obtained by examining the distribution of $`(V/V_{max})_i`$ of the whole data set.

Schmidt (1968) also described a method for determination of $`\varphi (L)`$, which was later dubbed as the ”Schmidt Estimator” by Felten (1976). It is straightforward to show that in absence of evolution, i.e. for a uniform spatial distribution, the contribution to this cumulative luminosity function of each source is proportional to $`V_{max,i}^1`$. It is also easy to show that if there is evolution this contribution is proportional to $`\sigma _{max,i}^1`$ defined in equation (2); $`\sigma `$ was denoted as $`V^{}`$ by Schmidt.

Most other methods employ binning, which simplifies the problem conceptually but has several shortcomings, the primary being loss of data points in the incomplete bins at the truncation boundaries. Examples of these are anlysis by Nicole & Segal (1978, 1983), Turner (1979) and Choloniewski (1986, 1987). As shown by Petrosian (1986), it turns out that all these procedures, in the limit of one source per bin reduce to Lynden-Bell’s (1979) $`C^{}`$ method. For a detailed comparison of these methods see Petrosian (1992).

### 3.3 The General Method

The right panel of Figure 1 depicts two boxes. Let us assume that number of data points in the narrow box is $`n(L)dL`$ and the number in the big box, excluding the narrow region is $`N(L)`$. It is easy to show that if the variables are independent, i.e. the luminosity function is separable in $`L`$ and $`z`$, then

$$\frac{n(L)dL}{N(L)}=\frac{\sigma (z_{max})\psi (L)dL}{\sigma (z_{max})\varphi (L)}=d\mathrm{ln}\varphi (L).$$
(3)
In the limiting case of one object per (narrow) bin, $`n(L)=\delta (LL_i)`$ and $`N(L)=N_i+\mathrm{\Theta }(LL_i)`$, where $`\mathrm{\Theta }`$ is the Heviside step function. Integration of the above equation then shows that the cumulative luminosity function $`\varphi (L)`$ increases by $`\delta \mathrm{ln}\varphi (L)=\mathrm{ln}(1+N_i^1)`$ going across the source at $`L_i`$. Thus, the cumulative luminosity function can be build in increments starting, say from the highest observed luminosity $`L_1`$, as

$$\varphi (L)=\varphi (L_1)\underset{L_i>L}{}\mathrm{ln}(1+N_i^1).$$
(4)
The set of $`N_i`$ sources in the big box are referred to as the comparable or the associated set of the source with luminosity $`L_i`$. It is clear that because of the complete mathematical symmetry between the two variables we can get an identical expression for the cumulative distribution $`\sigma (z)`$. The (different) associated sets are defined in a similar manner.

In the next section I will describe the extension of this method to more complex truncations. I will consider only bivariate distributions. The generalizations to multivariate distributions is straightforward.

## 4 THE NEW METHODS

The complete description of the distributions is a two step process. The first step is to determine whether the two variables are correlated or they are independent, and if correlated, then find a way to account for this correlation. The latter step can be done only parametrically. For the luminosity function it entails the determination of the form of the luminosity evolution function $`g(z)`$ such that a redefined luminosity $`L_0=L/g(z)`$ is uncorrelated with, or is independent of, $`z`$. The second step is to determine the univariate distributions of the independent variables $`z`$ and $`L_0`$.

### 4.1 Untruncated Data

If $`z`$ and $`L`$ are independent then the rank $`R_i`$ of $`z_i`$ (or $`L_i`$) in an untruncated sample (i.e. a sample truncated parallel to the axes; $`z>z_{min}`$ and $`L>L_{min}`$) will be distributed uniformly between $`1`$ and $`N`$ with an expected mean $`E=\frac{1}{2}(N+1)`$ and variance $`V=\frac{1}{12}(N^21)`$. We may then normalize $`R_i`$ to have a mean of $`0`$ and a variance of $`1`$ by defining the statistic $`T_i=(R_iE)/V`$. The hypothesis of independence is then rejected or accepted using a statistics based on the distribution of the $`T_i`$. The quantity

$$\tau =\frac{\underset{i}{}(R_iE)}{\sqrt{_iV}}$$
(5)
is one choice of such a test statistic with a mean of $`0`$ and a variance of $`1`$. The hypothesis of independence is rejected if $`|\tau _{data}|`$ is too large (e.g. $`|\tau _{data}|1`$ for rejection of independence at the $`1\sigma `$ level). This $`\tau `$ is equivalent to Kendell’s $`\tau `$ statistic (see, e.g. Press et al. 1990)

If it turns out that $`|\tau _{data}|1`$ and that the variables are independent, then the determination of the univariate distributions in each variable is obtained ignoring the value of the other. For example, the cumulative luminosity function will be described by the histogram $`\varphi (L_{i+1}>L>L_i)=i`$. However, if the variables are correlated, one must then carry out a transformation to remove the correlation by some parametric function, say $`L_0=L/g(z)`$. It is important to note that the transformed data will now appear truncated in the $`L_0z`$ plane with $`L_{0,min}=L_{min}/g(z)`$, which is no longer parallel to the $`z`$ axis. The use of the method described by equation (3) is then required to obtain the univariate distributions.

### 4.2 Data with One-sided Truncation

A straightforward application of the above method to a truncated data will clearly give a false correlation signal. Efron & Petrosian (1992), and independently Tsai (1990), describe how this method can be applied to data with one-sided truncation. The above procedure is modified as follows. For each object define a new comparable or associated set

$$J_i=\{j:L_j>L_i,L_j^{}<L_i\},$$
(6)
where $`L_j^{}=4\pi d_L^2(\mathrm{\Omega }_i,z_j)f_{min}`$. It is easy to see that this is the same set defined above (Figure 1) as the big box containing $`N_i`$ sources. (Note that the set defined in Efron & Petrosian (1992) includes the object $`i`$ in question.) This is the largest subset of luminosity and volume limited data that can be constructed for each point $`(L_i,z_i)`$. If $`z`$ and $`L`$ are independent then we expect the rank $`R_i`$ of $`z_i`$ (or $`L_i`$) in this limited set, not in the whole sample to be uniformly distributed between $`1`$ and $`N_i`$. The rest of the procedure follows the steps described above.

Similarly the determination of the univariate distributions will require the use of the method described in the previous section for the transformed and uncorrelated variables. Note again that the truncation boundary also gets transformed in case there is a correlation between $`L`$ and $`z`$.

### 4.3 Complex Truncations

A generalization of the above method to doubly (or multiply) truncated data was developed by Efron & Petrosian (1999), which is valid for the most general truncation $`L_i^{}<L_i<L_i^+`$ and $`z_i^{}<z_i<z_i^+`$. The method is equivalent to the previous method, with the associated set defined as

$$J_i=\{j:L_j>L_i,L_i(L_j^{},L_j^+)\}.$$
(7)
In this case, however, the distribution of the rankings (or of $`\tau `$) is unknown. If the data are uncorrelated then $`\tau `$ must still have a mean of zero. But a bootstrap method using simulations based on the $`\psi (L)`$ obtained from the data, as described below, is required for the purpose of the estimation of the variance.

For doubly (or more complexly) truncated data the comparable set is not completely observed, thus a simple analytic method such as the one described in equation (4) is not possible. However, it turns out that a simple iterative procedure can lead to a maximum likelihood estimate of the distributions. Efron and Petrosian (1999) give a thorough description of this method; for a more brief and transparent description see Maloney & Petrosian (1999).

### 4.4 Tests of the Correlation Algorithm

Lloyd et al. (2000) describe two simulations which test how well the above procedures reproduce known distributions. In one simulation the method was applied to a single sided truncation of an uncorrelated bivariate parent distribution. The rank test applied to the untruncated simulated points gave a value of $`\tau =0.9`$ But when applied to the truncated data without considering the effects of the truncation resulted in $`\tau >5.0`$, indicating the presence of a strong (of course false) correlation that is introduced by the truncation. However, when the test was carried out correctly by accounting for the effects of the truncation (as described above) it was found that $`\tau =0.6`$, recovering the fact that the variables were independent. In a second set the simulated parent sample had a strong anti correlation, with variables $`y`$ and $`x`$ correlated with average value of $`\overline{y}x^{1.5}`$. Application of the correlation to the untruncated randomly selected sample, as expected, gave a negative value for the statistic; $`\tau <5`$. The data was truncated and the method applied blindly. This resulted in the value of $`\tau =+2.5`$, or a (false) positive correlation. However, when applied properly it gave a value of $`\tau =4.0`$ indicating the presence of the anti correlation at 4 sigma level. Assuming a correlation of $`\overline{y}x^\delta `$, the data was transformed accordingly and the method applied to the new data. Figure 2 shows the variation of $`\tau `$ with the exponent $`\delta `$, from which we deduce a mean value and one sigma range of $`\delta =1.51\pm 0.11`$. This is a strong support for the accuracy of the method.

Further simulations are required to test the methods for determination of the distribution and correlations for more complex truncations.

## 5 SOME RESULTS FROM AGN DATA

Our application of these methods to a combine set of several optically selected samples of quasars are described in Maloney & Petrosian (1999). Here is a brief summary.

$``$ We found a strong correlation between luminosity and redshift, indicating the presence of a rapid luminosity evolution.

$``$ The parametric model of luminosity evolution $`(1+z)^k^{}`$ provides a better description of the data than the model $`e^{kt(z)}`$, where $`t(z)`$ is the look back time. Neither parameterization perfectly removes the correlation in all areas of the $`Lz`$ plane. In order to better model this evolution future analyses of quasar evolution could consider parametric forms, with more than one free parameter. Some of the more complex forms of evolution suggested in the literature, e.g. the used by La Franca & Cristiani (1996) are equivalent with a luminosity evolution form that contains two independent free parameters.

Given the form of the luminosity evolution we make the simple transformation of all luminosities to their hypothetical present epoch values, $`L_0=L/g(z)`$, so that $`L_0`$ and $`z`$ are uncorrelated. This allows us to use our methods to determine the univariate distributions of $`z`$ (the density evolution) and $`L_0`$ (the local luminosity function).

$``$ We find that the co-moving density of quasars also evolves, but its extent depends on the cosmological model. For example, for the Einstein$``$de Sitter model $`\rho (1+z)^{2.5}`$ for low redshifts and rapidly declines as $`\rho (1+z)^5`$ for $`z>2`$. This is much slower evolution than is obtained when one (incorrectly) assumes a pure density evolution model; $`g(z)=1`$ (Schmidt 1968; Miyaji, Hasionger & Schmidt 1998).

$``$ The cumulative local luminosity function $`\varphi (L_o)`$ has the double power law form found previously (Caditz & Petrosian 1990, La Franca & Cristiani 1996), with a break luminosity of $`L_{}=6\times 10^{29}`$ erg / (sec Hz), in the Einstein$``$de Sitter model. The power-law indices at the low and high luminosity ends are -1.5 and -2.3 in rough agreement with previous estimates. There is however, some evidence for evolution of the shape of the luminosity function. More data is required for a quantitative description of this evolution.


format: text

text:
# 2MASSI J1315309-264951: An L Dwarf with Strong and Variable H𝛼 Emission

## 1 Introduction

For many years the coolest, lowest-mass stars known were M dwarfs, but in recent years the L and T dwarfs have extended the stellar sequence to even lower temperatures and masses (Kirkpatrick et al. 1999). These three classes of low-mass stars and brown dwarfs exhibit distinctly different spectral features due to the decrease in photospheric temperature from M through L to T. In addition, M dwarfs often show H$`\alpha `$ in emission. Down to the early M dwarfs, H$`\alpha `$ activity correlates with rotation and thus decreases with age as stars lose angular momentum over time via stellar winds. Beyond objects of spectral type M8, however, it appears that H$`\alpha `$ activity is stronger in more massive objects, even if they are older (Gizis et al. 2000). The frequency of H$`\alpha `$ emission peaks around type M7 and declines for later-type L and T dwarfs, reaching zero at L5 (Gizis et al. 2000). However, contrary to this trend, Burgasser et al. (2000) reported the discovery of a T dwarf with strong H$`\alpha `$ emission. Here I report the discovery of an L3 dwarf with similarly strong H$`\alpha `$ emission.

## 2 Spectroscopy

2MASSI J1315309$``$264951 (hereafter 2MASS 1315$``$2649) was targeted for spectroscopy as part of a color-selected survey of AGN candidates from the Two-Micron All-Sky Survey Second Incremental Data Relase (2MASSI). That survey will appear elsewhere (Hall 2002, in preparation).

Potential AGN candidates were selected from the 2MASSI Point Source Catalog by requiring $`\alpha <5^h`$ or $`\alpha >15^h`$, $`55\mathrm{°}<\delta <+5\mathrm{°}`$, $`|b|>25\mathrm{°}`$, and that they not be affected by blending, saturation or confusion with other objects, known artifacts or minor planets. Actual AGN candidates were then selected from the list of potential candidates using color criteria based on those of Barkhouse & Hall (2001). 2MASS 1315$``$2649 was included in a subsample of objects with $`JK_s`$$`>`$1.5 and $`K_s`$$`<`$13.5.

Spectra of 2MASS 1315$``$2649 were secured during two runs with CSPEC at the CTIO 1.5m telescope. Both runs used a 300 lines/mm grating blazed at 6750 Å, for a resolution of 8.6 Å given the 1$`\stackrel{}{\mathrm{.}}`$5 slit. Wavelength coverage was 5780$``$9270 Å at 2.89 Å/pixel. Two 12-minute exposures were begun at UT 06:09 on 2001 March 30. A 12-minute and a 24-minute exposure were obtained through cirrus beginning at UT 23:39 on 2001 August 15. The two-dimensional spectra were debiased and flatfielded in the usual manner using IRAF.111The Image Reduction and Analysis Facility (IRAF) is distributed by NOAO, which is operated by AURA, Inc., under contract to the NSF. Each pair of reduced two-dimensional spectra was then coadded with exposure time weighting and 5$`\sigma `$ cosmic ray rejection.

The two-dimensional spectra show a close pair of objects. 2MASS 1315$``$2649 at $`\alpha `$=13:15:30.94, $`\delta `$=$``$26:49:51.3 (J2000) is an optically faint L dwarf with H$`\alpha `$ emission, and USNO J131531.23$``$264953.0 an optically brighter star slightly south of east. Both objects were visible on the guide camera, and the east-west slit was positioned as best as possible to obtain spectroscopy of both. However, unlike USNO J131531.23$``$264953.0, 2MASS 1315$``$2649 is faint enough that it may have drifted off the slit without being immediately noticed. The east-west separation between the objects in the two-dimensional spectra is 5$`\stackrel{}{\mathrm{.}}`$5$`\pm `$0$`\stackrel{}{\mathrm{.}}`$3, as determined from the centroids of the spatial profiles of the two objects in the spectra from the two independent epochs. This is significantly larger than the 3$`\stackrel{}{\mathrm{.}}`$9 east-west separation expected for the 2MASS-USNO separation of 4$`\stackrel{}{\mathrm{.}}`$23 at PA 114.5° E of N. This suggests that one or both of the objects exhibited proper motion between the 2MASS epoch of 1998.411 and the spectroscopic observations. However, better data are needed to verify this since CSPEC delivers a 2-3 pixel FWHM spectrum at a spatial scale of 1$`\stackrel{}{\mathrm{.}}`$3/pixel.

The spectra were optimally extracted using IRAF, with the USNO star’s trace used for both objects. Flux calibration was provided by GD108 (Oke 1990) in March and LTT1788 (Hamuy et al. 1994) in August. The spectral type of USNO J131531.23$``$264953.0 must be no later than K5, since its spectrum lacks strong features. The spectra of 2MASS 1315$``$2649 are discussed next.

## 3 Analysis

Figure 1 shows the March and August spectra of 2MASS 1315$``$2649. By comparison to the spectral atlas of Kirkpatrick et al. (2000), I estimate a spectral class between L2 and L4 on the Kirkpatrick et al. (1999) system, and adopt L3. Due to the poor quality of the spectra, this classification is not based on single features, but on the agreement between numerous features. The absence of strong TiO and VO features at 7000-8400 Å puts the class no earlier than L2. The flatness of the spectrum from 8400-9000 Å and the continuum level there compared to 7500 Å are inconsistent with spectral types L5 and higher.

H$`\alpha `$ is visible in both the two-dimensional and one-dimensional spectra from both observing runs, but it is clear that at least its equivalent width (EW) varied between the two runs. The total H$`\alpha `$ flux may not have been measured accurately in either run because the east-west slit was comparable in size to the seeing and was not centered precisely on 2MASS 1315$``$2649. Also, the August spectrum was obtained through cirrus and so its absolute fluxing is untrustworthy. On the other hand, EW measurements from both runs will be accurate unless the H$`\alpha `$ emission is extended, for which there is no evidence in the two-dimensional spectra, or the contamination of the L dwarf spectrum by the USNO star spectrum was different.

To estimate the contamination of the L dwarf spectrum by the USNO star spectrum, in each two-dimensional spectrum the L dwarf aperture was reflected in the spatial direction around the centroid of the USNO star aperture and the flux in the mirrored aperture was extracted. This contaminating flux is negligible at $``$8000 Å in both epochs, but is larger at $``$8000 Å in the August spectrum, yielding the stronger continuum seen there in Figure 1b. This makes the August EW measurement an underestimate, but does not affect the spectral typing of the object.

The H$`\alpha `$ EW was 121$`\pm `$31 Å in March but only 25$`\pm `$10 Å in August. The H$`\alpha `$ flux was measured to be 3.23$`\pm `$0.82 $`\times `$ 10-16 ergs s-1 cm-2 in March. Even including a correction factor of 1.42 for the nonphotometric conditions of the August measurement estimated by scaling the flux from 8000-9260 Å to match the March data, the August H$`\alpha `$ flux is only 1.49$`\pm `$0.58 $`\times `$ 10-16 ergs s-1 cm-2, about a factor of two below the March measurement. Thus it does appear that the absolute flux as well as the equivalent width of the H$`\alpha `$ emission varied.

The ratio of the H$`\alpha `$ luminosity (or flux) to the bolometric luminosity (or flux) is a useful quantity since it can constrain the emission mechanism. The bolometric flux of 2MASS 1315$``$2649 was estimated using the bolometric correction from the $`K_s`$ band of BCK=3.33 found by Tinney et al. (1993) for the L4 dwarf GD165B.222The bolometric correction is not always negative, as erroneously stated on p. 381 of Cox (2000), because the zero point of the magnitude system in which it is calculated can be different from that of the bolometric magnitude system. This yields an apparent bolometric magnitude of $`m_{\mathrm{bol}}=16.79`$. Since $`f_{\mathrm{bol}}=2.48\times 10^5`$ ergs s-1 cm-2 for $`m_{\mathrm{bol}}=0`$ (Cox 2000), $`f_{\mathrm{bol}}=4.77\times 10^{12}`$ ergs s-1 cm-2 for 2MASS 1315$``$2649. Thus $`\mathrm{log}(L_{\mathrm{H}\alpha }/L_{\mathrm{bol}})4.17`$ in March, and $`4.51`$ in August (after accounting for nonphotometric conditions). Both values are lower limits due to the use of a narrow slit with uncertain placement relative to the object.

## 4 Discussion

How does the H$`\alpha `$ emission in 2MASS 1315$``$2649 compare to that in other late-type dwarf stars? Gizis et al. (2000) studied the H$`\alpha `$ emission properties of a sample of nearby M and L dwarfs and found that $``$20$`\pm `$10% of L2-L4 dwarfs show H$`\alpha `$ emission. However, the strongest H$`\alpha `$ emission object in their sample3332MASSW J0326137+295015, whose spectrum can be seen in Figure 16 of Kirkpatrick et al. (1999). has only $`\mathrm{log}(L_{\mathrm{H}\alpha }/L_{\mathrm{bol}})5`$. Subsequently, Burgasser et al. (2000) reported H$`\alpha `$ emission with $`\mathrm{log}(L_{\mathrm{H}\alpha }/L_{\mathrm{bol}})4.3`$ in the T dwarf 2MASSI J1237392+652615. Possible variability in this object’s H$`\alpha `$ flux was reported at the 2.8$`\sigma `$ level in one spectrum, but more detailed analysis shows that the spectra are consistent with no variability (Burgasser et al. 2002).

Thus even the August spectrum of 2MASS 1315$``$2649 has stronger H$`\alpha `$ emission than reported in any other L dwarf, and the variable H$`\alpha `$ emission in this object is the first reported for any L dwarf. The study of Gizis et al. (2000) indicates that as an L dwarf with H$`\alpha `$ emission, 2MASS 1315$``$2649 is likely massive enough (and old enough) to have burned lithium, but the physical explanation for this correlation with mass is unknown. The various explanations considered by Burgasser et al. (2000) for H$`\alpha `$ emission in the T dwarf 2MASS 1237+6526 can also be considered for 2MASS 1315$``$2649. Acoustic heating is ruled out as a possible dominant energy source for strong H$`\alpha `$ emission in L or T dwarfs (§3.3 of Burgasser et al. 2000), but several other possibilities remain viable.

### 4.1 Accretion In An Interacting Binary System

Burgasser et al. (2000) discuss the possibility of sustained Roche lobe overflow in close brown dwarf binaries. Such accretion might explain steady H$`\alpha `$ emission without a strong accompanying thermal spectrum if the accretion produces H$`\alpha `$ emission from ionization at a shock front or from magnetic field lines streaming onto the pole of the primary, rather than from an accretion disk. Burgasser et al. (2000) suggest that the T dwarf 2MASS 1237+6526 and the M9.5e dwarf PC 0025+0447 could be examples of such systems. Such an explanation cannot be ruled out for 2MASS 1315$``$2649, but its confirmed variable H$`\alpha `$ would require variability in some part of the accretion process. Such variability has not been convincingly observed in the other two objects and might not be expected in the case of accretion from sustained Roche lobe overflow.

### 4.2 A Strong Magnetic Field

Magnetic fields are believed to drive H$`\alpha `$ emission from F to early M stars via an internal dynamo which is stronger for faster-rotating stars, such as the $`\alpha `$-$`\mathrm{\Omega }`$ dynamo (Parker 1955). Such dynamos require a radiative/convective boundary to anchor flux lines and therefore break down as stars become fully convective ($``$0.3 $`M_{}`$, spectral type M4). However, the observed H$`\alpha `$ activity level in M dwarfs shows no sign of this transition, remaining constant at $`\mathrm{log}(L_{\mathrm{H}\alpha }/L_{\mathrm{bol}})3.8`$ (Hawley, Gizis, & Reid 1996). This suggests that some other mechanism, probably a turbulent dynamo (Durney, De Young, & Roxburgh 1993), contributes substantially to the magnetic field in very late type dwarfs.

The value of $`\mathrm{log}(L_{\mathrm{H}\alpha }/L_{\mathrm{bol}})`$ in 2MASS 1315$``$2649 is an order of magnitude larger than in most L dwarfs, but is still lower than the average value of $``$3.8 measured for M dwarfs (see Figure 7 of Gizis et al. 2000). Thus it is at least plausible that the same mechanism at work in late M dwarfs (whether a turbulent dynamo or something else) could provide the magnetic field required to explain the H$`\alpha `$ flux in 2MASS 1315$``$2649 and possibly the T dwarf 2MASS 1237+6526. This would mean that they just have stronger than average magnetic fields for L and T dwarfs.

If this is the case, then the variable H$`\alpha `$ emission in 2MASS 1315$``$2649 requires at least a slow variation in its magnetic field. However, no late type dwarf is known to show long-term variability of more than a factor of two at X-ray or extreme ultraviolet wavelengths (§ 4 of Drake et al. 1996), and variability in H$`\alpha `$ is not likely to be greater than variability at those wavelengths. On the other hand, flares are known to occur in M dwarfs and have been detected at radio wavelengths in the L3.5 dwarf 2MASSW J0036159+182110 (Berger 2002). Thus while the spectra presented here do not rule out a slow variation in its magnetic field as the cause of the H$`\alpha `$ variability in this L dwarf, a more rapid variation — a flare — seems much more likely, as discussed next.

### 4.3 Flaring

Flares (seen in H$`\alpha `$, the radio, the UV, and X-rays) on late type stars are powered by the energy released during a sudden reconfiguration of the magnetic field structure. Detection of a flare therefore confirms the existence of a stellar magnetic field. However, strong flares in M dwarfs typically last only a few hours (Hawley & Pettersen 1991) and occur unpredictably, making them difficult to observe.

It may very well be that 2MASS 1315$``$2649 was flaring when the discovery spectra were obtained. The H$`\alpha `$ emission did not vary more than 1$`\sigma `$ (25%) on a timescale of $``$15 minutes between successive two-dimensional spectra in March, but this is consistent with the rate of change of the H$`\alpha `$ flux in the M9.5 dwarf 2MASS 0149+2956 during the flare observed by Liebert et al. (1999). Given that H$`\alpha `$ flaring is known to occur in M dwarfs and radio flaring in at least one L dwarf (Berger 2002), flaring seems the most likely explanation for the H$`\alpha `$ variability in 2MASS 1315$``$2649. However, the other possibilities discussed earlier cannot yet be ruled out. The current data does not prove that the variation in H$`\alpha `$ occurred on a timescale short enough to be called a flare, and radio, H$`\alpha `$ and X-ray and activity may be decoupled for spectral types later than M8 (Berger 2002), Further monitoring of 2MASS 1315$``$2649 is needed to firmly determine the energy source for its H$`\alpha `$ emission.

If its variable H$`\alpha `$ emission is due to flaring, 2MASS 1315$``$2649 is probably similar to the M9.5 dwarf BRI 0021$``$0214 (Reid et al. 1999). Both objects show only H$`\alpha `$ in emission and have peak $`\mathrm{log}(L_{\mathrm{H}\alpha }/L_{\mathrm{bol}})`$ values lower than the average for M dwarfs. However, the spectra do not rule out the possibility that 2MASS 1315$``$2649 is a flare star with emission from lines other than H$`\alpha `$, similar to 2MASS 0149+2956 (Liebert et al. 1999). There is no evidence for such lines in March to limits of 0.5$`f_{\mathrm{H}\alpha }`$, but in 2MASS 0149+2956 the strongest other lines had 0.1$`f_{\mathrm{H}\alpha }`$. There is also no evidence in the March spectra for the stronger continuum and veiling of molecular bands seen during the flare in 2MASS 0149+2956. Only changes of a factor of two or more in the continuum level can be ruled out, but this is sufficient to place 2MASS 1315$``$2649 at either an earlier or later stage in any such flare, as follows. The H$`\alpha `$ EW in 2MASS 0149+2956 increased (from 200 Å to 330 Å) during the observed portion of the flare because the continuum flux at H$`\alpha `$ decreased more quickly (from ten to three times the quiescent flux) than the H$`\alpha `$ line flux did. In 2MASS 1315$``$2649 the H$`\alpha `$ EW was 121 Å in March, but the continuum flux at H$`\alpha `$ was at most twice the quiescent flux (estimated from the August spectrum). This H$`\alpha `$ EW and continuum level cannot be simultaneously reproduced by simply scaling any of the 2MASS 0149+2956 observations, even accounting for the weaker continuum at H$`\alpha `$ in an L3 dwarf. Thus, if they were the same type of flare, and if the relationship between the continuum and H$`\alpha `$ line flux in this type of flare is independent of its luminosity, then 2MASS 1315$``$2649 must have been observed at a later (or perhaps earlier) stage of its flare. Any model for such flares must then incorporate the observed limit of $``$25% variation in H$`\alpha `$ on $``$15 minute timescales at such a stage in the flare.

If the H$`\alpha `$ emisssion in 2MASS 1315$``$2649 in March was from a flare, regardless of what type, then the large range in $`\mathrm{log}(L_{\mathrm{H}\alpha }/L_{\mathrm{bol}})`$ for flaring M dwarfs extends to L dwarfs as well. As pointed out by Reid et al. (1999), this suggests a wide variation in intrinsic magnetic field strength, efficiency of energy transport to the stellar chromosphere during a flare, or both.

## 5 Conclusion

2MASSI J1315309$``$264951 is an L3 dwarf with strong H$`\alpha `$ emission which decreased in strength by about a factor of two between two epochs separated by 137 days, the first reported variable H$`\alpha `$ emission in an L dwarf. The H$`\alpha `$ emission in 2MASSI J1315309$``$264951 must be powered either by magnetic activity or by accretion in a binary system. Accreting binaries are rare, so that hypothesis is unlikely. The spectra presented here do not rule out a slow variation in H$`\alpha `$ strength, but slow variations of the observed amplitude are rare among M dwarfs. Since flaring powered by reconnection of magnetic fields is common in M dwarfs and a radio flare has been detected in the L3.5 dwarf 2MASSW J0036159+182110 (Berger 2002), a flare is the logical explanation for the H$`\alpha `$ variability in 2MASSI J1315309$``$264951,

The value of $`\mathrm{log}(L_{\mathrm{H}\alpha }/L_{\mathrm{bol}})>4.17`$ observed at the discovery epoch is larger than that of any other L dwarf but comparable to the value of $`4.3`$ observed for 2MASSI J1237392+652615, the only reported T dwarf with H$`\alpha `$ emission. However, both these values lie well below the average $`\mathrm{log}(L_{\mathrm{H}\alpha }/L_{\mathrm{bol}})=3.8`$ observed in M dwarfs.

Only two L dwarfs and one T dwarf are known to exhibit H$`\alpha `$ emission of strength $`\mathrm{log}(L_{\mathrm{H}\alpha }/L_{\mathrm{bol}})>5`$. Thus perhaps two percent of L or T dwarfs exhibit H$`\alpha `$ emission this strong at any given time. Given the small number statistics, this is consistent with the duty cycles observed for H$`\alpha `$ flares ($``$7%; Gizis et al. 2000) and radio flares (2-10%; Berger 2002) among late-type dwarfs. Spectroscopic or narrow-band H$`\alpha `$ monitoring of L and T dwarfs on timescales of hours to days is needed to determine if the frequency of strong H$`\alpha `$ emission is governed by flaring, by the upper envelope of the magnetic field strength distribution, or by accretion in binary systems.

I thank an anonymous referee and A. Bugasser for helpful comments, CNTAC for observing time, A. Alvarez and A. Guerra for operating the CTIO 1.5m and putting up with my musical tastes, CTIO for excellent observing support, and Fundación Andes and Chilean FONDECYT grant #1010981 for financial support. The Two Micron All Sky Survey (2MASS) is a joint project of the University of Massachusetts and the Infrared Processing and Analysis Center/California Institute of Technology, funded by the National Aeronautics and Space Administration and the National Science Foundation.


text:
# 1 Introduction

## 1 Introduction

This communication has been inspired by a report on the theory of a perfect lens by Pendry . This lens is supposedly constructed of a material whose permittivity and permeability, respectively, are exactly negative of the permittivity and the permeability of free space. In a predecessor paper , that material was postulated as the anti–vacuum.

## 2 Canonical Formulation for a Lens

As geometrical optics is applicable for lenses, the canonical formulation for a lens merely involves a linear homogeneous material confined to the region between two parallel planes. Let us, however, generalize the situation in order to understand the issue at greater depth by involving two linear homogeneous materials and four interfaces, as shown in Figure 1. The regions $`0zd_1`$ and $`d_1+d_2zd_1+d_2+d_3`$ are occupied by a material labelled $`a`$, and the region $`d_1zd_1+d_2`$ by a material labelled $`b`$. Both materials are linear, homogeneous, bianisotropic and necessarily dispersive, their frequency–domain constitutive relations being as follows:

$$\begin{array}{c}𝐃(x,y,z,\omega )=ϵ_0\left[\underset{¯}{\underset{¯}{ϵ}}^{a,b}(\omega )\text{ }\text{}\text{ }𝐄(x,y,z,\omega )+\underset{¯}{\underset{¯}{\alpha }}^{a,b}(\omega )\text{ }\text{}\text{ }𝐇(x,y,z,\omega )\right]\hfill \\ 𝐁(x,y,z,\omega )=\mu _0\left[\underset{¯}{\underset{¯}{\beta }}^{a,b}(\omega )\text{ }\text{}\text{ }𝐄(x,y,z,\omega )+\underset{¯}{\underset{¯}{\mu }}^{a,b}(\omega )\text{ }\text{}\text{ }𝐇(x,y,z,\omega )\right]\hfill \end{array}\}.$$
(1)
In these relations, the dielectric properties are delineated by $`\underset{¯}{\underset{¯}{ϵ}}^{a,b}(\omega )`$, the magnetic properties by $`\underset{¯}{\underset{¯}{\mu }}^{a,b}(\omega )`$, and the magnetoelectric properties by $`\underset{¯}{\underset{¯}{\alpha }}^{a,b}(\omega )`$ as well as $`\underset{¯}{\underset{¯}{\beta }}^{a,b}(\omega )`$, all of these dyadics being functions of the angular frequency $`\omega `$. The permittivity and permeability of free space (i.e., vacuum) are denoted by $`ϵ_0`$ and $`\mu _0`$, respectively.

Figure 1: Schematic for the canonical lens formulation.

Without loss of essential generality, we can take the spatial Fourier transformations of all electromagnetic phasors with respect to $`x`$ and $`y`$; thus,

$$𝐄(x,y,z,\omega )=𝐞(z,\kappa ,\psi ,\omega )\mathrm{exp}\left[i\kappa (x\mathrm{cos}\psi +y\mathrm{sin}\psi )\right],$$
(2)
etc., suffice for the present purposes. Wave propagation in the two materials can then be cast in terms of 4$`\times `$4 matrix ordinary differential equations as follows :

$$\frac{d}{dz}[\underset{¯}{f}(z,\kappa ,\psi ,\omega )]=i[\underset{¯}{\underset{¯}{P}}_{a,b}(\kappa ,\psi ,\omega )]\text{ }\text{}\text{ }[\underset{¯}{f}(z,\kappa ,\psi ,\omega )].$$
(3)
Here, $`[\underset{¯}{f}][e_x;e_y;h_x;h_y]^T`$ is a column 4–vector with the superscript T denoting the transpose, while $`[\underset{¯}{\underset{¯}{P}}_{a,b}]`$ are 4$`\times `$4 matrixes. Accordingly, we obtain the basic relation

$$[\underset{¯}{f}(d_1+d_2+d_3,\kappa ,\psi ,\omega )]=[\underset{¯}{\underset{¯}{M}}(d_1+d_2+d_3,\kappa ,\psi ,\omega )]\text{ }\text{}\text{ }[\underset{¯}{f}(0,\kappa ,\psi ,\omega )],$$
(4)
where

$`[\underset{¯}{\underset{¯}{M}}(d_1+d_2+d_3,\kappa ,\psi ,\omega )]`$ $`=`$ $`\mathrm{exp}\left\{id_3[\underset{¯}{\underset{¯}{P}}_a(\kappa ,\psi \omega )]\right\}\text{ }\text{}\text{ }\mathrm{exp}\left\{id_2[\underset{¯}{\underset{¯}{P}}_b(\kappa ,\psi \omega )]\right\}`$ (5)
$`\text{ }\text{}\text{ }\mathrm{exp}\left\{id_1[\underset{¯}{\underset{¯}{P}}_a(\kappa ,\psi \omega )]\right\}.`$
Within the confines of continuum electromagnetics, the canonical lens problem involves finding the material $`b`$ and the thickness $`d_2`$ for a specified material $`a`$ and thicknesses $`d_1`$ and $`d_3`$, such that

$$[\underset{¯}{\underset{¯}{M}}(d_1+d_2+d_3,\kappa ,\psi ,\omega )]=\left[\begin{array}{cccc}1& 0& 0& 0\\ 0& 1& 0& 0\\ 0& 0& 1& 0\\ 0& 0& 0& 1\end{array}\right]$$
(6)
for all $`\kappa `$, $`\psi `$ and $`\omega `$.

## 3 Analysis and Discussion

Obviously, that would be a fruitless endeavor in practice. Hence, lens designers effectively settle for some acceptable ranges of $`\kappa `$, $`\psi `$ and $`\omega `$ in which (6) holds, when material $`a`$ is air. Deviations from an ideal match introduce aberrations .

Mathematically, and at first glance, an excellent candidate for ideal match is the following:

$$\begin{array}{cc}\underset{¯}{\underset{¯}{ϵ}}^b(\omega )=\underset{¯}{\underset{¯}{ϵ}}^a(\omega ),\hfill & \underset{¯}{\underset{¯}{\alpha }}^b(\omega )=\underset{¯}{\underset{¯}{\alpha }}^a(\omega )\hfill \\ \underset{¯}{\underset{¯}{\beta }}^b(\omega )=\underset{¯}{\underset{¯}{\beta }}^a(\omega ),\hfill & \underset{¯}{\underset{¯}{\mu }}^b(\omega )=\underset{¯}{\underset{¯}{\mu }}^a(\omega )\hfill \\ d_2=d_1+d_3\hfill & \end{array}\}.$$
(7)
Close inspection of (6), however, shows that (7) is suitable for all $`\kappa `$, only if material $`a`$ has orthorhombic symmetry, i.e.,

$$\underset{¯}{\underset{¯}{\xi }}^a=\left[\begin{array}{ccc}\xi _1^a& 0& 0\\ 0& \xi _2^a& 0\\ 0& 0& \xi _3^a\end{array}\right],\xi =ϵ,\alpha ,\beta ,\mu .$$
(8)
A sandwich of equal thicknesses of materials $`a`$ and $`b`$ then constitutes a planar realization of a medium named nihility earlier in this journal . Nihility is the postulated electromagnetic nilpotent, with the following constitutive relations:

$$\begin{array}{c}𝐃(x,y,z,\omega )=\mathrm{𝟎}\hfill \\ 𝐁(x,y,z,\omega )=\mathrm{𝟎}\hfill \end{array}\}.$$
(9)
Wave propagation cannot occur in nihility, because $`\times 𝐄(x,y,z,\omega )=\mathrm{𝟎}`$ and $`\times 𝐇(x,y,z,\omega )=\mathrm{𝟎}`$ in the absence of sources therein. Whereas the phase velocity and the wavevector of a plane wave in vacuum/anti–vacuum are co–parallel/anti–parallel, the directionality of the phase velocity relative to the wavevector in nihility is a non–issue.

Physically, (7) and (8) are still deficient because the principle of energy conservation has not been considered. With the normal constraint that material $`a`$ be passive, it follows from (7) that material $`b`$ has to be active. Although their electromagnetic response can be simulated via composites containing active circuit elements , I do not think that active materials will provide a realistic option in the near future. Therefore, the passivity constraint on material $`b`$ would lead to aberrations due to absorption. To those must be added chromatic aberrations, which would arise from the non–fulfilment of (7) outside some limited range of $`\omega `$.

Suppose next that material $`a`$ is isotropic and non–magnetoelectric, i.e., $`\underset{¯}{\underset{¯}{ϵ}}^a=ϵ^a\underset{¯}{\underset{¯}{I}}`$, $`\underset{¯}{\underset{¯}{\mu }}^a=\mu ^a\underset{¯}{\underset{¯}{I}}`$, and $`\underset{¯}{\underset{¯}{\alpha }}^a=\underset{¯}{\underset{¯}{\beta }}^a=\underset{¯}{\underset{¯}{0}}`$, where $`\underset{¯}{\underset{¯}{I}}`$ is the identity dyadic and $`\underset{¯}{\underset{¯}{0}}`$ is the null dyadic. Then, a reasonable match in some frequency range could conceivably be provided by materials that supposedly possess a negative index of refraction , so long as absorption is acceptably low.

Pendry  actually took material $`a`$ to be air, which is optically indistinguishable from vacuum ($`\underset{¯}{\underset{¯}{ϵ}}^a=\underset{¯}{\underset{¯}{I}}`$, $`\underset{¯}{\underset{¯}{\mu }}^a=\underset{¯}{\underset{¯}{I}}`$, and $`\underset{¯}{\underset{¯}{\alpha }}^a=\underset{¯}{\underset{¯}{\beta }}^a=\underset{¯}{\underset{¯}{0}}`$) for most purposes. Thus, the requirements on material $`b`$ became very simple, viz., $`\underset{¯}{\underset{¯}{ϵ}}^b=\underset{¯}{\underset{¯}{I}}`$, $`\underset{¯}{\underset{¯}{\mu }}^b=\underset{¯}{\underset{¯}{I}}`$, and $`\underset{¯}{\underset{¯}{\alpha }}^b=\underset{¯}{\underset{¯}{\beta }}^b=\underset{¯}{\underset{¯}{0}}`$. (In other words, material $`b`$ has to be anti–vacuum.) But Ziolkowski  has recently concluded from two–dimensional computer simulations that even these simple requirements (in some narrow frequency range) cannot be met by realistic meta–materials. Recent correspondence between Pendry and others adds to the debate .

A perfect lens remains unrealizable in my opinion, as does nihility.


text:
# Crystal experiments on efficient beam extraction

(Presented at Int. Conf. on Atom. Coll. in Solids
(ICACS-19: Paris, July 29 - August 3, 2001))

## Abstract

Silicon crystal was channeling and extracting 70-GeV protons from the U-70 accelerator with efficiency of 85.3$`\pm `$2.8% as measured for a beam of $``$1012 protons directed towards crystals of $``$2 mm length in spills of $``$2 s duration. The experimental data follow very well the prediction of Monte Carlo simulations. This success is important to devise a more efficient use of the U-70 accelerator in Protvino and provides a crucial support for implementation of crystal-assisted collimation of gold ion beam in RHIC and slow extraction from AGS onto E952, now in preparation at Brookhaven Nat’l Lab. Future applications, spanning in the energy from sub-GeV (medical) to order of 1 GeV (scraping in the SNS, extraction from COSY) to order of 1 TeV and beyond (scraping in the Tevatron, LHC, VLHC), can benefit from these studies.

The technique of bent crystal channeling to steer particle beams, with applications in extraction from accelerator and in beam collimation, has progressed rapidly thanks to efforts at IHEP Protvino, CERN, and FNAL. Beams of up to 106 proton/s were extracted from CERN SPS and Tevatron by Si crystals of just 4 cm in length, with typical efficiencies on the order of 10-20%. It was predicted that efficiency of crystal channeling extraction can be boosted to much higher values by multiple particle encounters with a shorter crystal installed in a circulating beam. To clarify this mechanism a new experiment was started at IHEP at the end of 1997, with intention to test very short crystals and achieve very high efficiencies of extraction.

The benefits of a crystal-assisted extraction are fourfold. In hadron colliders this mode of extraction can in general be made compatible with the colliding mode of operation. The time structure of the extracted beam is practically flat, since the extraction mechanism is resonance-free. The size of the extracted beam is smaller and more round than in a resonant extraction. Finally polarized beams can be extracted without detrimental effects on the polarization.

The classic two-stage collimation system for loss localisation in accelerators typically uses a small scattering target as a primary element and a bulk absorber as secondary element. The role of the primary element is to give a substantial angular kick to the incoming particles in order to increase the impact parameter on the secondary element, which is generally placed in the optimum position to intercept transverse or longitudinal beam halo. An amorphous primary target scatters the impinging particles in all possible directions. Ideally, one would prefer to use a ”smart target” which kicks all particles in only one direction: for instance, only in radial direction, only outward, and only into the preferred angular range corresponding to the centre of the absorber (to exclude escapes). A bent crystal is the practical implementation for such a smart target: it traps particles and conveys them into the desired direction .

These two possible applications of crystal channeling in modern hadron accelerators, extraction and halo collimation, can be exploited in a broad range of energies, from sub-GeV cases (i.e. for medical accelerators) to multi-TeV machines (for high-energy research). Indeed, several projects are in progress to investigate them. Crystal collimation is studied at RHIC (100-250 GeV), considered for the Tevatron (1000 GeV) and the LHC, for the Spallation Neutron Source (1 GeV), whilst crystal-assisted slow extraction is considered for AGS (25 GeV) and COSY (1-2 GeV). In all cases, the critical issue is the channeling efficiency.

In the last two years, we demonstrated crystal channeling with 50% efficiency. We also showed that, when properly aligned, these crystals could be efficiently used as primary collimators, thereby reducing by a factor two the radiation level measured downstream of the collimation region of U-70. To continue our investigations, we installed and tested several new crystals in different straight sections of the U-70 ring. Three of them were produced by different manufacturers with a new shape. They were made with narrow strips of Si(111), about 40 mm long vertically and a fraction of mm long in the radial direction. Their azimuthal length was only a few mm. They were bent by a metallic holder providing deflections of 0.8 to 1.5 mrad.

To make a crystal deflector this short (factor of $``$20 shorter than the crystals used at CERN and FNAL), a new design has been developed, Fig.1. The crystal strip has the shape of a saddle, being bent in both vertical and radial directions. The advantages of ”new-generation” crystals are threefold: (a) they can be made shorter along the beam than in previous designs, (b) they have no straight ends as the bending mechanism is continuous, and (c) they have no amorphous material close to the beam (like the ”legs” of U- and O-shaped deflectors used earlier in CERN and IHEP).

Two crystals were assembled in Protvino: one $`2`$ mm long was bent by 0.9 mrad, the other $`4`$ mm long was bent by 1.5 mrad. The third crystal, $`1.8`$ mm long, 0.8 mrad bent, was prepared at the University of Ferrara and chemically polished at optical level to remove the defects induced during diamond slicing. The two Russian crystals were used in extraction mode, whilst the Italian one was tested as a primary collimator. The three crystals were exposed to 70 GeV proton beams and used to channel and extract halo particles.

Figure 2 illustrates the beneficial effect of crystals when used as primary collimators. It shows beam profiles in the radial direction downstream of the crystal as measured on the entry face of the collimator. Four cases are reported. First, an amorphous collimator is used as primary target whilst the crystal is kept outside of the beam envelope. As expected, the beam profile is peaked at the collimator edge, Fig.2(a). In the second case, Fig.2(b), the crystal is used as the primary scraper but is not aligned to the beam. When properly aligned, Fig.2(c), the crystal channels most of the incoming beam into the depth of the collimator. In the last case, Fig.2(d), the beam is simply kicked by a magnet towards the secondary collimator, whilst the crystal is retracted.

The extraction efficiency is given by the ratio of the extracted beam intensity as measured in the external beam line to all the beam loss as measured in the entire ring. We obtained very high efficiencies in each of the three new crystals: namely, both the 1.8 and 2 mm long crystals reached 85% efficiency, whilst the 4mm long crystal reached 68% efficiency. These striking results were obtained in a steady manner over many runs.

In Fig.3 we plot the predicted (by Monte Carlo code) and the measured extraction efficiencies together with the data obtained earlier. The agreement between measurements and simulations is excellent.

Beside the channeling efficiency, the ability to withstand a high beam intensity and the crystal lifetime also important for a practical application. Crystals located in the region upstream of the U-70 cleaning area were irradiated with the entire circulating beam, spilt out in a rather short time durations to simulate very dense halo collimation. Analysis has shown that our crystals were irradiated up to $`2\times 10^{14}`$ particles per spill of $``$1 s duration. When averaged over machine cycles, the irradiation rate was as high as $`2\times 10^{13}`$ proton hits/s. Notice that this irradiation rate compares with the expected beam loss rate at the Spallation Neutron Source. Indeed, the SNS Accumulator Ring should generate a 1 GeV proton flux of 60$`\times `$2$`\times `$1014 per second. At the expected rate of beam loss of 0.1–1% the halo flux will be (1.2–12)$`\times `$1013 protons/s. Crystals may well tolerate high flux of protons as large as those expected to hit the SNS beam collimation system.

Another crystal, $`5`$ mm long, was exposed for several minutes to even higher radiation flux. Beams of up to 1013 70-GeV protons (resulting in $``$1014 proton hits) were directed towards the crystal in spills of 50 ms every 9.6 s. After this extreme exposure, the channeling properties of the crystal were tested in an external beam. The deflected beam observed with photo-emulsion  was perfectly normal, without breaks, nor significant tails eventually produced by dechanneled particles. Several crystals in use in U-70 have been exposed to high intensity beams for months, thereby accumulating very large integrated doses. After the estimated irradiation of $``$1020p/cm2 the initial channeling efficiency of about 43% was practically unaffected. This irradiation is still below the world highest results. The CERN experiment showed that at 5$`\times `$1020 p/cm2 the crystal lost only 30% of its deflection efficiency, which means $``$100 years lifetime in the intense beam of NA48 experiment.

On the same location in U-70 ring with the same 1.8-mm crystal of Si(111) positioned $``$20 m upstream of the ring collimator, we have repeated the crystal collimation experiment at the injection flattop of U-70, proton kinetic energy of 1.3 GeV. With the crystal aligned to the incoming halo particles, the radial beam profile at the collimator entry face showed a significant channeled peak far from the edge, Fig.6. About half of the protons intercepted by the collimator jaw, have been channeled there by a crystal; i.e., the crystal has doubled the amount of particles intercepted by the jaw. As only part (about 34%) of all particles scattered off the crystal have reached the jaw, we estimate the crystal deflection efficiency as 15-20%. The observed figure of efficiency could be well reproduced in computer simulations. This figure is orders of magnitude higher than previous world data for low-GeV energy range. It is remarkable that the same crystal was efficiently channeling and really helping in cleaning halo particles both at 70 GeV and at 1.3 GeV, thus demonstrating to be operational in a very wide energy range.

In summary, the crystal channeling efficiency has reached unprecedented high values both at top energy and at injection energy. The same 2 mm long crystal was used to channel 70 GeV protons with an efficiency of 85.3$`\pm `$2.8% during several weeks of operation and 1.32 GeV protons with an efficiency of 15-20% during some test runs. Crystals with a similar design were able to stand radiation doses over $`10^{20}`$ proton/cm2 and irradiation rates of $`2\times 10^{14}`$ particles incident on crystal in spills of $``$2 s duration without deterioration of their performances.

The efficiency results well match the figures theoretically expected for ideal crystals. As simulations show, extraction and collimation with channeling efficiencies over 90-95% is feasible. The obtained high figures provide a crucial support for the ideas to apply this technique in beam cleaning systems, for instance in RHIC and Tevatron. Earlier Tevatron scraping simulations have shown that a crystal scraper can reduce accelerator-related backgrounds in CDF and D0 experiments by a factor of $``$10. This year, first experimental data is expected from RHIC where a crystal collimator is installed. The technique presented here is potentially applicable also in LHC for instance to improve the efficiency of the LHC cleaning system by embedding bent crystals in the primary collimators. This work was supported by INTAS-CERN grant 132-2000, RFBR grant 01-02-16229, and by the ”Young researcher Project” of the University of Ferrara.


text:
# Dynamical ultrametricity in the critical trap model

## I Introduction

Notable theoretical progress in our understanding of the ubiquitous aging phenomena in glassy systems was made possible by the recognition that the discussion of correlation and response functions requires two times: the waiting time $`t_w`$ and the total elapsed time $`t_w+t`$. This appears very clearly in the framework of mean-field spin glass models, domain growth (coarsening models), or more phenomenological trap models . Although the basic phenomenology of all these models are rather similar, the underlying physical picture is completely different. For example, aging in domain growth models is associated to the growth of a coherence length. In mean field or trap models, space is absent and aging is related to the structure of the energy landscape, but here again the intuition is completely different. In mean field models, the system never reaches the bottom of an energy valley and there are no activated processes involved in the dynamics. Rather, the dynamics slows down because saddles with less and less ‘descending’ directions are visited as the system ages. Conversely, in the trap model, activation is the basic ingredient of the model, and aging is associated to the fact that deeper and deeper valleys are reached as the system evolves. Dynamics in the latter case in fundamentally intermittent: either nothing moves, or there is a jump between two traps. This must be contrasted with mean field dynamics which is continuous in time. However, mean field and trap dynamics can be shown to correspond to two successive time regimes in a particular class of models .

In spite of these important differences, many predictions are common to the latter two pictures, such as:

* A short time singularity of the response function in the aging regime, which leads to an (aging) low frequency noise.
* Non trivial violations of the Fluctuation Dissipation Theorem (fdt), first pointed out within mean field models, but that also exists in trap models.
* The possibility of rejuvenation and memory, which involves the existence of different degrees of freedom with different time scales.

Furthermore, a certain class of mean field models (that includes the Sherrington-Kirkpatrick model) has been conjectured to possess ultrametric dynamical properties, that very precisely reflect and encode the ultrametric nature of the static solution. This dynamical ultrametricity is associated to an infinite number of time scales (which diverges with the age of the system), in the following sense: if $`C(t_2,t_1)`$ is the correlation function between times $`t_1<t_2`$, then in the limit of large times:

$$C(t_3,t_1)=\mathrm{min}(C(t_2,t_1),C(t_3,t_2)),t_2[t_1,t_3].$$
(1)
This means that either $`t_2`$ is close enough to $`t_3`$, and then no further dynamics takes place between $`t_2`$ and $`t_3`$, or $`t_2`$ is close enough to $`t_1`$ but then the age of the system hardly changes between $`t_1`$ and $`t_2`$. This property of the correlation function has been shown to hold for the aging solution of the dynamical (Mode Coupling) equations describing the dynamics of ‘continuous’ spin-glasses. This property is furthermore invariant under reparametrizations of time, where $`th(t)`$ with an arbitrary monotonous function $`h(t)`$. Testing whether or not dynamical ultrametricity also holds in realistic disordered systems is made difficult because this property is only expected in the limit of asymptotically large times, and corrections are expected on finite times. How large are these corrections?

In this paper, we show that exact ultrametricity holds at the critical point of the single-level trap model (or random energy model). We give the explicit form of the correlation function and discuss finite time (or finite forcing) corrections. Note that in this single-level trap model, the dynamics is ultrametric although the statics is not. The issue of finite time fdt plot is also adressed. We discuss multi-level extensions of the trap model and argue that dynamical ultrametricity should be generic at ‘short times’, i.e. at the beginning of the aging region. We show thermoremanent magnetization data that support this idea. The relation with $`1/f`$ noise, already discussed in this context , is recalled.

## II The model

The trap model, introduced in the context of aging in  and further developed in , is one of the simplest soluble model exhibiting a dynamical glass transition. In this model, one considers a particle which is trapped in low energy states $`i`$ of depth $`E_i>0`$, where the $`E_i`$ are random variables distributed according to $`\rho (E)=\frac{1}{T_c}e^{E/T_c}`$. The dynamics is chosen to be activated: each particle stays in trap $`i`$ an exponential random time, equal on average to $`\tau _0e^{E_i/T}`$. The quantity $`\tau _0`$ is a microscopic time scale which we shall take as the time unit in the following. When the particle leaves the trap, it chooses at random a new one among all the others. As a consequence, at high temperature ($`T>T_c`$), the particle spends most of its time in the small traps, because the number of these traps (the entropic factor) dominates the Boltzmann factor and the system equilibrates. On the contrary, for $`T<T_c`$, the Boltzmann factor is dominant and the particle explores deeper and deeper traps, so that the system never equilibrates (in the limit of an infinite number of traps). In this regime, the dynamics ages: correlation and response functions are no longer time translation invariant, but depend both on the waiting time $`t_w`$ and the total time $`t_w+t`$.

## III Correlation function

Correlation functions are useful tools to characterize the dynamics and to compare several models. The simplest (but nevertheless non trivial) correlation in the trap model is $`\mathrm{\Pi }(t_w+t,t_w)`$ defined as the probability to remain in the same trap during the time interval $`[t_w,t_w+t]`$. As we are considering the infinite dimensional (or fully connected) model, this probability is also equal to the probability $`P(t_w+t,t_w)`$ to be in the same trap at $`t_w+t`$ as at $`t_w`$, since the probability to go back to the same trap vanishes in the limit of an infinite number of traps. (This is not true in finite dimension. For instance, in one dimension, $`\mathrm{\Pi }(t_w+t,t_w)`$ and $`P(t_w+t,t_w)`$ scale in a different way with $`t_w`$).

The correlation function $`\mathrm{\Pi }(t_w+t,t_w)`$ has been calculated in the general case in . We shall now focus on the critical case $`T=T_c`$, for which we have:

$$\mathrm{\Pi }(t_w+t,t_w)=_1^{t_w}𝑑u\frac{1\mathrm{exp}[(t_w+tu)]}{(\mathrm{log}u)(t_w+tu)}$$
(2)
for large $`t_w`$. We note that this expression simplifies further in the limit $`t,t_w\mathrm{}`$. In particular, the exponential term vanishes since $`t_w+tu>t`$. Writing $`u=t_w(1v)`$, we get

$`\mathrm{\Pi }(t_w+t,t_w)`$ $``$ $`{\displaystyle _0^{1\frac{1}{t_w}}}{\displaystyle \frac{dv}{[\mathrm{log}t_w+\mathrm{log}(1v)][v+\frac{t}{t_w}]}}`$ (3)
$``$ $`{\displaystyle \frac{1}{\mathrm{log}t_w}}{\displaystyle _0^{\frac{t_w}{t}}}{\displaystyle \frac{dz}{z+1}}+𝒪\left({\displaystyle \frac{1}{(\mathrm{log}t_w)^2}}\right)`$ (4)
with the new integration variable $`z=\frac{t_w}{t}v`$. This expression is readily integrated to give:

$$\mathrm{\Pi }(t_w+t,t_w)\frac{\mathrm{log}(1+\frac{t_w}{t})}{\mathrm{log}t_w}$$
(5)
This relation shows that $`\mathrm{\Pi }(t_w+t,t_w)`$ is not a function of $`\frac{t}{t_w}`$, at variance with the results that hold the whole low temperature phase $`T<T_c`$ . On the contrary, taking $`tat_w^\alpha `$, we obtain, in the limit $`t_w`$ going to $`+\mathrm{}`$:

$$\mathrm{\Pi }(t_w+t,t_w)=𝒞(\alpha )$$
(6)
with $`𝒞(\alpha )`$ given by:

$`𝒞(\alpha )`$ $`=`$ $`1\alpha (\alpha <1)`$ (7)
$`=`$ $`0(\alpha 1).`$ (8)
Note that $`𝒞(\alpha )`$ is a monotonous decreasing function of $`\alpha `$. The above result is only true in the limit $`t\tau _0`$, as well as $`\mathrm{log}(t_w/\tau _0)1`$. Indeed, all the finite time results reported in this paper should be understood as first order terms in a $`1/\mathrm{log}t_w`$ expansion.

An important remark is that the correlation function $`\mathrm{\Pi }(t_w+t,t_w)`$ is a function of $`\alpha =\mathrm{log}t/\mathrm{log}t_w`$ that cannot be written as $`h(t_w+t)/h(t_w)`$. The latter ratio naturally appears (with an unknown function $`h`$) in the aging part of the solution of the dynamical equation corresponding to one step Replica Symmetry Breaking (rsb) mean field spin glass models. However, for full rsb models where dynamical ultrametricity is indeed expected, the correlation function is given by an infinite sum of contributions coming from different times sectors :

$$C(t_w+t,t_w)=\underset{i}{}𝒞_i\left(\frac{h_i(t_w+t)}{h_i(t_w)}\right),$$
(9)
where $`h_i`$ are unknown (monotonous) functions defining the $`i`$th time scale, and the $`𝒞_i`$ are monotonously decaying to zero for large arguments.

A useful (but up to now unjustified theoretically) form for $`h_i(t)`$, that allows one to give some flesh to the above formula, is :

$$h_i(u)=\mathrm{exp}\left[\frac{u^{1\mu _i}}{1\mu _i}\right],0\mu _i1.$$
(10)
It is easy to see that for this choice of $`h_i`$, the time scale on which the ratio $`h_i(t_w+t)/h_i(t_w)`$ varies significantly is precisely $`t_w^{\mu _i}`$. The choice $`\mu =0`$ therefore corresponds to stationnary dynamics, whereas $`\mu =1`$ gives full aging. Now, take $`t=t_w^\alpha `$ (with $`0<\alpha <1`$) in Eq. (9) and take the limit $`t_w\mathrm{}`$. All the sectors such that $`\mu _i<\alpha `$ have relaxed to zero, whereas the sectors corresponding to $`\mu _i>\alpha `$ have not decayed at all. Introducing a continuum of different values of $`\mu `$, we find that the correlation function is given by:

$$C(t_w+t_w^\alpha ,t_w)=_\alpha ^1𝑑\mu \rho (\mu )𝒞_\mu (1),$$
(11)
where $`\rho (\mu )`$ is the ‘density’ of time sectors of order $`t_w^\mu `$ and $`𝒞_\mu (1)`$ is the initial value of the correlation function in this sector. From this result, one sees that $`C(t_w+t,t_w)`$ indeed becomes a function of $`\alpha =\mathrm{log}t/\mathrm{log}t_w`$ in the long time limit. Therefore, interestingly, the superposition of an infinite number of subaging contributions defined by (10) naturally leads to a correlation function that depends on $`\mathrm{log}t/\mathrm{log}t_w`$, for which the dynamical ultrametricity property is explicit. The critical trap behaviour corresponds to a uniform contribution of all time sectors, i.e. $`\rho (\mu )𝒞_\mu (1)=1`$, $`\mu `$.

## IV Dynamical ultrametricity

As recalled in the Introduction, Cugliandolo and Kurchan have defined dynamical ultrametricity for the correlation function $`C`$ if the following property is true: take three times $`t_1<t_2<t_3`$, then:

$$C(t_3,t_1)=\mathrm{min}(C(t_2,t_1),C(t_3,t_2)),t_2[t_1,t_3].$$
(12)
Let us show now that $`\mathrm{\Pi }(t_w+t,t_w)`$ at the critical temperature is ultrametric in the sense defined hereabove. It will be useful to introduce the following notations:

$`\mathrm{\Pi }(t_2,t_1)=C_1`$ (13)
$`\mathrm{\Pi }(t_3,t_2)=C_2`$ (14)
$`\mathrm{\Pi }(t_3,t_1)=C_3`$ (15)
From the monoticity of correlation functions, the inequality $`C_3\mathrm{min}(C_1,C_2)`$ holds in general. We simply have to check that (at least) one of the two correlations $`C_1`$ and $`C_2`$ is equal to $`C_3`$. In order to take the infinite time limit, we need to specify how $`t_2`$ and $`t_3`$ scale with $`t_1`$. A natural parametrization is the following:

$`t_3t_1`$ $``$ $`at_1^\alpha `$ (16)
$`t_2t_1`$ $``$ $`bt_1^\beta `$ (17)
We can now look at various cases:

* If $`\beta <\alpha `$, then one has, for large $`t_1`$:

$$t_3t_2at_1^\alpha bt_1^\beta at_1^\alpha $$
(18)
so that $`C_1=𝒞(\beta )`$ and $`C_2=C_3=𝒞(\alpha )<𝒞(\beta )`$.
* Assuming now $`\beta =\alpha `$, we get $`C_1=C_3=𝒞(\alpha )`$, as well as $`C_2C_3`$.

As a result, we have shown that the relation $`C_3=\mathrm{min}(C_1,C_2)`$ always holds, which implies that dynamical ultrametricity is satisfied in this model.

The appearance of dynamical ultrametricity can be considered as a signature of the existence of many time scales involved in the dynamics. This is indeed the case in the present model, even though there is no static ultrametricity to account for this hierarchy of time scales. This property in fact arises naturally from the exact balance, at the critical point, of the Boltzmann weight $`e^{E/T}`$ and of the entropic factor $`\rho (E)`$: the probability for the particle to have a given energy (or equivalently, a given trapping time, in logarithmic scale) is essentially uniform on the interval $`[0,T_c\mathrm{log}t_w]`$. In other words, dynamical ultrametricity here is a consequence of the critical scale invariance.

## V Finite time analysis

An interesting analysis was introduced by Cugliandolo and Kurchan in the context of the dynamical analysis of the Sherrington-Kirkpatrick (SK) model . These authors made the assumption that, in the limit of large times, there exists a certain function $`f(x,y)`$, not necessarily smooth, such that $`C_3=f(C_1,C_2)`$. We have shown in the previous section that such a function indeed exists in the present case and is given by $`f(x,y)=\mathrm{min}(x,y)`$.

A useful representation, proposed in , is to plot in the $`(C_1,C_2)`$-plane the curves of constant $`C_3=f(C_1,C_2)`$, which reduces for our case to two straight lines ($`C_1=C_3`$ and $`C_2=C_3`$) at right angle. For finite times, the function $`f`$ has to include a time scale as third argument, so that $`C_3=F(C_1,C_2;t_1)`$, where $`f(x,y)=lim_{t_1\mathrm{}}F(x,y;t_1)`$. The $`(C_1,C_2)`$-plane representation is then a good way to visualize the convergence towards the asymptotic function $`f(C_1,C_2)`$. It has been used with numerical data to test if dynamic ultrametricity holds in realistic systems, with rather inconclusive results , except in a few cases, like in a recent study of the 4-dimensional Edwards-Anderson model .

Let us apply this procedure to the critical trap model. In order to deal with finite time expressions, we shall come back to Eqn. 5, restated as follows:

$$C(t_j,t_i)=\frac{\mathrm{log}(1+\frac{t_i}{t_jt_i})}{\mathrm{log}t_i}$$
(19)
Inverting these relations so as to express $`t_2`$ and $`t_3`$ as functions of $`t_1`$, $`C_1`$ and $`C_3`$, we can write an explicit expression for $`F(C_1,C_2;t_1)`$:

$`C_3`$ $`=`$ $`F(C_1,C_2;t_1)`$ (20)
$`=`$ $`{\displaystyle \frac{1}{\mathrm{log}t_1}}\mathrm{log}\left(t_1^{C_1}+t_1^{C_2}(1t_1^{C_1})^{1+C_2}\right)`$ (21)
In order to plot the constant $`C_3`$ curves, it will be useful to express also $`C_2`$ as a certain function $`G(C_1,C_3;t_1)`$:

$$C_2=G(C_1,C_3;t_1)=\frac{\mathrm{log}(1t_1^{C_1})\mathrm{log}(t_1^{C_3}t_1^{C_1})}{\mathrm{log}t_1\mathrm{log}(1t_1^{C_1})}$$
(22)
The resulting plots in the $`(C_1,C_2)`$-plane are displayed in Fig. 1, for $`C_3=0.1,0.3,0.5,0.7`$ and for three different times ($`t_1=10^5`$, $`10^{10}`$ and $`10^{15}`$). Note that the convergence is very slow close to the infinite time singularity.

Cugliandolo and Kurchan also introduced in  the notion of correlation time scales through the representation of $`f(C,C)`$ versus $`C`$. As already mentionned before, $`f(C,C)C`$ in general. There may exist some special ‘fixed’ points $`C^{}`$ such that $`f(C^{},C^{})=C^{}`$. Each of these fixed points has been shown to be associated with a correlation time scale. If dynamical ultrametricity holds in a particular time sector, all $`C`$’s belonging to a certain interval $`[C^{},C^{\prime \prime }]`$ are fixed points. In our case, ultrametricity holds over the full correlation interval $`[0,1]`$. But in our model, we can go beyond the infinite time analysis and quantify the convergence of $`F(C,C;t_1)`$ with $`t_1`$ towards the asymptotic function $`f(C,C)=C`$. We find:

$$F(C,C;t_1)=C\frac{\mathrm{log}(1+(1t_1^C)^{1+C})}{\mathrm{log}t_1}$$
(23)
For $`C>0`$, this expression simplifies further at large times:

$$F(C,C;t_1)C\frac{\mathrm{log}2}{\mathrm{log}t_1}$$
(24)
Interestingly, the leading correction does not depend on $`C`$, and it is valid only if $`C`$ is not too close to $`0`$. Fig. 2 displays the plots of $`F(C,C;t_1)C`$, for the same values of $`t_1`$ as in Fig. 1.

This last result may also be interpreted in the framework of Fig. 1. For a given value of $`C_3`$, the point $`C(t_1)`$ defined by the relation $`F(C(t_1),C(t_1);t_1)=C_3`$ converges to the infinite time right angle singularity $`C_1=C_2=C_3`$ as (see Eqn. 24):

$$C(t_1)C_3+\frac{\mathrm{log}2}{\mathrm{log}t_1}$$
(25)
This logarithmic correction may explain why it seems so difficult to observe the convergence towards the ultrametric relation in experimental or numerical data, where only a few decades (usually between four and six) are available.

## VI The effect of ‘shear’

The effect of an external ‘shear’ (or power injection) on aging was investigated in the context of mean field models in , and in the context of the trap model in . In both models, aging is interrupted by the shear beyond a time scale $`\tau _r`$ which diverges as the shear rate $`\dot{\gamma }`$ tends to zero. In the model considered in , this time scale is given by:

$$\tau _r\frac{1}{\dot{\gamma }}\left(\mathrm{log}\frac{1}{\dot{\gamma }}\right)^{\frac{1}{2}}.$$
(26)
For $`\tau \tau _r`$, the (power-law) distribution of trapping times is in a first approximation unaffected by the shear, whereas for longer times, the distribution decays exponentially. In the limit where the waiting time $`t_w`$ is much larger than $`\tau _r`$, the dynamics of the model becomes stationary and one finds, for $`t\tau _r`$, the same result as above with $`t_w`$ replaced by $`\tau _r`$:

$$C(t+t_w,t_w)=𝒞\left(\frac{\mathrm{log}t}{\mathrm{log}\tau _r}\right)1\frac{\mathrm{log}t}{\mathrm{log}\tau _r}.$$
(27)
As discussed in , dynamical ultrametricity in this context manifests itself by the appearance of an infinity of time scales in the limit $`\tau _r\mathrm{}`$ (i.e. $`\dot{\gamma }0`$): the time needed for the correlation to decay to a certain value $`c`$ diverges as $`\tau _r^{1c}`$ (see  for a further discussion).

## VII The fluctuation-dissipation ratio

It is interesting to study the fluctuation dissipation ratio $`X`$ in the trap model at the critical point. This ratio is defined as:

$$X(t_2,t_1)=\frac{TR(t_2,t_1)}{C(t_2,t_1)/t_1},$$
(28)
where $`R(t_2,t_1)`$ is the response of the system at time $`t_2>t_1`$ to a small bias field applied at time $`t_1`$ (see  for details). In the trap model, where the field only changes the trapping time of the starting site, one finds that the following relation holds in general:

$$TR(t_2,t_1)=\frac{C(t_2,t_1)}{t_2}.$$
(29)
Using this result, one sees that in the ‘liquid’ phase $`T>T_c`$ where all two time functions only depend on time differences, $`X1`$: the usual Fluctuation-Dissipation Theorem (fdt) holds. In the glass phase, on the other hand, one finds that $`X(t_2,t_1)=t_1/t_2`$: the value of $`X`$ is non trivial in the whole scaling regime $`t_2t_1t_1`$. Right at the critical point $`T=T_c`$, one can express $`X`$ as:

$$X(t_2,t_1)=\frac{(t_1^C1)^2}{t_1^C(t_1^C+1C)}CC(t_2,t_1)>0.$$
(30)
So for any fixed $`C>0`$, $`X`$ tends to $`1`$ in the asymptotic limit $`t_1\mathrm{}`$. We show in Fig. 3 the now famous plot of the integrated response versus $`C`$, which should yield a straight line of slope $`1/T_c`$ when the fdt holds. Again, one sees a very slow convergence towards the asymptotic result for small values of $`C`$.

## VIII The multilevel trap model and discussion

The simple trap model described above can be considered as a one-step rsb model. In order to generalize the model to a full rsb one, it has been proposed in  (and further studied in ) to follow Parisi’s procedure for the static solution of the SK model. Roughly speaking, it means that each trap is recursively subdivided into a new series of traps, in a hierarchical manner. Each level $`k`$ of traps is characterized a certain overlap between states $`q_k`$ and by an exponential probability distribution of the energy barriers, with a critical temperature $`T_c^k`$ depending on the level index $`k`$. The critical temperatures are related to Parisi’s function $`x_k=x(q_k)`$ as $`T_c^k=T/x_k`$, and satisfy the relations $`T_c^k<T_c^{k1}`$. At any temperatures, the levels of the tree corresponding to $`q>q_{EA}(T)`$, where $`q_{EA}`$ is the Edwards-Anderson parameter, are such that $`x_k>1`$, so that these levels are equilibrated (i.e. $`T_c^k<T`$).

In the single level trap model, the correlation function $`\mathrm{\Pi }(t_w+t,t_w)`$ in the aging phase $`T<T_c`$ behaves at short times as:

$$\mathrm{\Pi }(t_w+t,t_w)1\frac{\mathrm{sin}\pi x}{\pi (1x)}\left(\frac{t}{t_w}\right)^{1x}1tt_w$$
(31)
with $`x=\frac{T}{T_c}`$. In the multi-level model with a finite number $`M`$ of levels, the total correlation function is defined as:

$`C(t_w+t,t_w)`$ $`=`$ $`{\displaystyle \underset{k=0}{\overset{M}{}}}q_k[\mathrm{\Pi }_k(t_w+t,t_w)\mathrm{\Pi }_{k+1}(t_w+t,t_w)]`$ (32)
$`=`$ $`q_0+{\displaystyle \underset{k=1}{\overset{M}{}}}(q_kq_{k1})\mathrm{\Pi }_k(t_w+t,t_w)`$ (33)
$`q_k`$ being the $`k^{th}`$ level overlap, and $`\mathrm{\Pi }_k(t_w+t,t_w)`$ is the probability that the process has never jumped beyond the $`k^{\mathrm{th}}`$ layer of the tree between $`t_w`$ and $`t_w+t`$ with the convention that $`\mathrm{\Pi }_0(t_w+t,t_w)=1`$ and $`\mathrm{\Pi }_{M+1}(t_w+t,t_w)=0`$ (see  for details). From this definition, we see that $`C(t_w+t,t_w)`$ is dominated at short times by the levels $`k`$ with $`x_k`$ close to $`1`$, for which the short time singularity is strongest. (We assume that $`T<T_c^0`$, i.e. that at least one level is aging). Therefore, we expect to observe the dynamical ultrametricity associated to the level $`k^{}`$ for which $`x_k^{}=1`$ in the ‘short’ time regime $`\mathrm{log}t/\mathrm{log}t_w<1`$, before the $`t/t_w`$ regime associated to the levels $`k<k^{}`$ sets in. (Note that if $`C(t_w+t,t_w)`$ is a function of $`t/t_w`$, then the function $`f(x,y)`$ defined above cannot be equal to $`\mathrm{min}(x,y)`$).

Interestingly, ‘short time’ dynamical ultrametricity exists for the hierarchical tree model in the whole low temperature phase, but has no relation with the static ultrametricity built in the tree structure which encodes Parisi’s rsb solution. Thus the origin of dynamical ultrametricity in the generalized trap model is again very different from the dynamical ultrametricity found in mean field spin glass models, which in the latter case is deeply related to the Parisi function $`x(q)`$ which encodes the structure of the tree. As mentioned in the introduction, the physical interpretation of aging in the two pictures are radically different, although some of the phenomenology is very similar.

Aging experiments is spin-glasses have been interpreted within the framework of the multi-level trap model in . The need for several levels comes from rejuvenation and memory in temperature shift experiments, but also from the detailed shape of the thermo-remanent magnetization (trm) relaxation at a given temperature, which shows that the short time and long time singularities are described by different exponents $`x_k`$ . It is natural to interpret the different levels of the hierarchy in terms of length scales, as is actually suggested by the rsb solution of the pinned manifold problem . Therefore, one can expect that for any temperature within the spin-glass phase, there will be a particular ‘critical’ length scale $`\mathrm{}`$ for which $`x=x(\mathrm{})=1`$, that will contribute to the correlation and response functions as a function of $`\mathrm{log}t/\mathrm{log}t_w`$. We have reanalyzed some trm data in the spirit of the present discussion. We show in Fig. (4) the decay of the trm in AgMn at $`T=0.75T_g`$ , plotted as a function of $`s(t,t_w)=\mathrm{log}(1+t_w/t)/\mathrm{log}t_w`$, as suggested by Eq. (5). This figure shows that the rescaling is very good at short times, but is violated for large $`t/t_w`$ (corresponding to small $`s(t,t_w)`$). From Eq. (32), we indeed expect to observe the sum of a contribution $`(k<k^{})`$ from the levels $`k<k^{}`$ (that only depends on $`t/t_w`$) and a contribution from $`kk^{}`$, proportional to $`s(t,t_w)`$. When $`\alpha =\mathrm{log}t/\mathrm{log}t_w<1`$, the first contribution does not vary much, since $`(t/t_w)^{1x_k}=t_w^{(\alpha 1)(1x_k)}`$ is very small for large $`t_w`$, whereas the second contribution is a function of $`\alpha `$. This suggests that in the short time regime one should observe:

$$M(t_w+t,t_w)=(k<k^{})+m_1s(t,t_w),$$
(34)
where $`m_1`$ is a certain prefactor. The quantity $`\phi =m_1/((k<k^{})+m_1)`$ measures the relative contribution of the levels $`<k^{}`$ and around $`k^{}`$ in the total decay of the signal. The dashed line shown in Fig. 4 is a linear fit of the initial decay, as a function of $`s`$, from which one extracts (in this particular example) $`\phi 0.2`$.

In all the above formulas, time is implicitely measured in units of a microscopic time $`\tau _0`$. The value of $`\tau _0`$ is not necessarily an individual flip time $`10^{12}`$ sec., since collective dynamics may exist in the vicinity of the transition point (see  for a detailed discussion). In Fig. 4, we have chosen $`\tau _0=10^5`$ sec. to achieve the best rescaling. This value is very close to the one extracted from the analysis of .

The conclusion is that the trm data is indeed compatible with a $`\mathrm{log}t/\mathrm{log}t_w`$ behaviour for short times. Note that this behaviour is tantamount to a logarithmic dependence of the a.c. susceptibility, or else to $`1/f`$ noise. We have insisted in previous papers  on the fact that the existence of levels in the vicinity of $`x=1`$ generically leads to $`1/f`$ noise for long times and low frequencies, because the contribution of other levels (both above and below $`k^{}`$) become negligible in that regime. We show here that the degrees of freedom contributing to $`1/f`$ noise also give rise to dynamical ultrametricity.

## Aknowledgments

We thank L. Berthier and E. Vincent for fruitful discussions.


text:
# Optimising the Linear Collider Luminosity: Feedback on Nanosecond Timescales

## I Introduction

The luminosity achievable in a colliding-beam accelerator may be written:

$$L=\frac{fNN_1N_2H}{4\pi \sigma _x\sigma _y}$$

where $`f`$ is the number of machine cycles per second, $`N`$ is the number of bunches per machine cycle, $`N_1`$ and $`N_2`$ are the number of particles per bunch, and $`\sigma _x`$, $`\sigma _y`$ are the transverse overlaps of the colliding bunches. For an $`e^+e^{}`$ collider $`H`$ is the beam-beam self-focussing parameter or ‘pinch enhancement’. The table shows the values of several of these parameters for the various designs of a next-generation $`e^+e^{}`$ linear collider currently under consideration. $`\mathrm{\Delta }t`$ is the time separation between the $`N`$ bunches in the train.

The nanometre-level vertical beam overlap is a particularly challenging goal for all these designs, most notably CLIC. Any source of beam motion which results in relative vertical offsets of the two beams at the interaction point (IP) at the nm level will clearly reduce the luminosity from the nominal value. In all of the collider designs stabilisation below the $`1\sigma `$ level is required to keep the luminosity loss below 10%.

The many kinds of potential beam motion may be characterised in two classes: (i) slow drifts resulting from eg. thermal excursions or component settling, with characteristic timescales varying from seconds to months; (ii) jitter on a timescale comparable with the machine repetition time. Both kinds of motion were experienced in the decade-long experience at SLC, and were dealt with by employing slow- and fast-feedback systems, respectively.

We are addressing the design of an intra-bunch-train fast-feedback (FB) system for the next-generation linear collider (LC). The system comprises a fast beam position monitor (BPM) to detect the relative misalignment of the leading electron and positron bunches at the IP, a feedback loop, and a fast kicker for kicking the trailing bunches back into collision.

The system time-response requirements for J/NLC and CLIC are clearly very different from those for TESLA. We have therefore chosen the more challenging case of J/NLC and CLIC as the primary focus of our efforts to develop a working prototype hardware system. However, from the timing point-of-view, a system which works on the 10ns scale could clearly be applied to the less demanding TESLA timescale of 300ns. With current technology it is clear that a system for J/NLC and CLIC must be based on fast analogue signal processing, rather than digital processing ingrid . However, given the rapid advances in signal processing speeds over the past decade it is quite possible that digital technology will be fast enough by the time the real LC system is deployed c. 2010; if so, the FB electronics could look almost identical for any of J/NLC, CLIC or TESLA. We aim to monitor this situation and to take advantage of technological developments as they arise.

## II Simulation Studies

During the past 18 months we have created powerful software tools for the simulation of a nanosecond timescale feedback system for correcting the relative displacement of nanometre-sized beams.

* We have imported and installed the code GUINEAPIG schulte . This allows us to simulate the beam-beam interaction between colliding electron and positron beams of arbitrary size and bunch charge. When they do not collide head-on each bunch gives the other a strong transverse EM kick which can be detected in a downstream BPM. This forms the physical input to the FB system (Figure 1).
* We have installed MATLAB and SIMULINK to create a modular FB system simulation. The response times of the BPM, feedback loop and kicker magnet, as well as cable delays, can all be chosen arbitrarily. The feedback algorithm can similarly be programmed at will.
* We have set up a GEANT model of the NLC interaction region to allow us to incorporate the FB system components. The location of the BPM and kicker directly affect the system latency (due to signal propagation times). In addition, the material of these components can contribute to knock-on backgrounds in the detector resulting from the showering of beam-produced photons and $`e^+e^{}`$. A corresponding model of the CLIC interaction region is currently being implemented.
* We have set up an equivalent FLUKA model to allow us to calculate background neutron production.

We have used these tools to simulate the luminosity loss due to beam misalignment for J/NLC, TESLA and CLIC (Figure 2). For the NLC and TESLA cases we have set up a feedback model which simulates the performance of the real hardware components: the BPM, the kicker, the FB logic/loop, with signal delay times corresponding to the actual locations of the components within the respective IR. In addition, for NLC we have evaluated the neutron and $`e^+e^{}`$ pair backgrounds in the IR that result from interactions with the material of the added components; the equivalent simulations will be done for CLIC.

We have used this software infrastructure to optimise the design of the FB system for the NLC case in terms of the minimisation of the response time and the knock-on backgrounds, subject to the constraints imposed by the locations of the other components in the crowded IR, principally the final-focus quadrupole magnets and the vertex detector. Figure 3 glen  shows the performance for the optimal hardware layout in terms of the luminosity loss vs. beam offset and kicker gain. For offsets below 10$`\sigma `$ there is a comfortable ‘valley’ in which the specific gain choice is not critical. However, for offsets significantly larger than 10$`\sigma `$ the choice of gain is delicate; an unfortunate choice could lead to catastrophic luminosity loss. We are now in a position to perform similar design optimisation studies for CLIC and TESLA.

## III FONT1: First Beam Test at the NLC Test Accelerator

The importance of FB systems in accelerator stabilisation/control and luminosity optimisation was learned painfully at SLC. No future linear collider can be physically aligned for nm-size beam-beam collisions. The Feedback on Nanosecond Timescales (FONT) system is key to successful operation of the future linear collider; without it the luminosity may be 1-2 orders of magnitude below the nominal design. It is therefore highly desirable to test experimentally a prototype FB system under conditions as close as possible to those at the LC phil .

For J/NLC or CLIC one requires an experiment with colliding trains of O(100) bunches with an inter-bunch separation of O(1 ns); no such facility exists, or is likely to exist, until the collider turns on. However, non-colliding bunch trains with similar properties can be produced. We have decided to base our first-round test, FONT1, at the NLC Test Accelerator (NLCTA), located in the Research Yard at SLAC. NLCTA is currently used for the evaluation of RF components such as Cu structures, but the downstream section is relatively ‘open’ beamline. NLCTA can be operated in a ‘long-pulse’ mode, which provides a train of bunches filled at X-band frequency, 11.4 GHz, typically of 120 ns duration. This train length is within a factor of 2 of J/NLC and CLIC requirements, and the total train charge is comparable, although the bunch spacing, 0.1 ns, is about an order of magnitude shorter; the prototype experiment is in some sense therefore ‘harder’ than the real case. We intend to use a dipole magnet to perturb the whole train, a fast BPM to monitor the displacement at a given point, an analogue FB circuit, a fast kicker to attempt to correct the displacement as fast as possible, and two downstream BPMs to measure and monitor the correction.

The key components of the system, and their status, are listed below:

* Dipole magnet: A suitable magnet has been found, a ‘type 4 linac corrector’, amongst spares at SLAC. It will produce a dipole offset (i.e. angular deviation) that corresponds to up to 5 mm position offset at the measurement BPM. The peak current is 6A and it is planned to use an existing power supply. The dipole can be pre-triggered so as to be fully ‘on’ by the time the bunch train arrives.
* Fast BPM: A prototype X-band ‘button’ BPM has been fabricated. The response of the BPM at X-band has been measured using a network analyser, as well as using the BPM test facility at Lawrence Berkeley Laboratory. Its resolution has been characterised and appears to be of the order of 10 $`\mu `$m, which is more than adequate for our test.
* Fast kicker: Three suitable fast kickers have been found amongst spares at SLAC; only one is needed for a 1-dimensional test.
* Kicker amplifier: A preliminary design study for a multi-stage tube amplifier capable of delivering up to 5 kW peak power has been made. The design employs Y690 planar triode tubes.
* Diagnostic BPMs: Two downstream X-band BPMs (see above) will be used to monitor the beam position with the FB in off/on modes.

## IV Schedule

We expect to take data on the FB system performance in spring and summer 2002. Initially a fixed offset will be applied, and we will attempt to ‘zero’ it as fast as possible within the bunch train passage time. The fixed offset will then be varied to simulate different beam offsets in the NLC, and we will characterise the performance vs. offset, investigating the effect of the choice of gain. The ultimate test will be to wire a random-number generator to the dipole and have the FB system operate ‘blind’.

Without pre-empting the results of these tests, it is possible to envision an extended R&D programme for investigation of a number of issues germane to a real FB system:

* Testing of improved analogue processing algorithms.
* Handling of correlated ‘jitter’ within the bunch train, eg. sinusoidal or other regular transverse position oscillations.
* Dynamic gain choice, either from adaptive learning or using a priori information from upstream systems.
* Interplay between position-correction and angle-correction; here we have addressed only the former.
* Correction in both planes: dealing with transverse beam coupling.
* Dealing with bunch tails in highly non-Gaussian (i.e. realistic) bunches.
* Upgrade to digital technology as signal processing speed improves.

Some of these could be carried out at the NLCTA, and/or at other facilities such as the ATF at KEK, TTF at DESY, or CTF3 which is under construction at CERN.