text:
# Anomalous relaxation and self-organization in non-equilibrium processes

\[

## Abstract

We study thermal relaxation in ordered arrays of coupled nonlinear elements with external driving. We find, that our model exhibits dynamic self-organization manifested in a universal stretched-exponential form of relaxation. We identify two types of self-organization, cooperative and anti-cooperative, which lead to fast and slow relaxation, respectively. We give a qualitative explanation for the behavior of the stretched exponent in different parameter ranges. We emphasize that this is a system exhibiting stretched-exponential relaxation without explicit disorder or frustration.

\]

The world around us is full of non-equilibrium and non-stationary processes, many of which are robust and easy to measure. However, there is no a priori reason to believe, that these non-equilibrium processes should fulfill universal laws, in a way equilibrium systems fulfill laws of thermodynamics and statistical mechanics. It is, therefore, important to recognize examples of complex non-equilibrium processes organizing themselves in a simple and universal way, a phenomenon termed “self-organization” or “emergent behavior” . Kolmogorov turbulence  is an example of such non-equilibrium process: Liquid, forced at macroscopic length scales, produces a flow of energy from large to small length scales, and this flow organizes itself into a stationary universal distribution. Turbulence is a steady, though non-equilibrium, process; the flow of energy is constant and the resulting self-organized probability distribution does not change in time. It is important to establish, whether boundaries of self-organization can be expanded to include non-equilibrium and non-stationary processes. In this paper we show that dynamic self-organization can indeed be found in relaxational dynamics of extended systems and is manifested in stretched-exponential dependence on time of physical quantities.

Stretched-exponential ($`\mathrm{exp}\left[(t/t_{rel})^\alpha \right]`$) relaxation laws have been observed in a large variety of physical and biological processes, such as recombination of carriers in semiconductors and polymers , protein relaxation  and folding , ligand binding to myoglobin , relaxation in magnetic clusters , superconducting vortices  and charge density waves , dynamics of alloys  and glasses . In the case of glasses, the stretching exponent $`\alpha `$ defines a glass transition temperature, $`T_g`$, i.e. $`\alpha =1`$ for $`T>T_g`$, and $`\alpha <1`$ for $`T<T_g`$. The non-exponential dynamics of glasses has long been related to the high degree of disorder, which leads to the existence of a large number of metastable states .

Observations of non-exponential behavior in simpler systems, such as magnetic clusters  and proteins , suggests that a high degree of disorder is not a necessary requirement for a system to display anomalous relaxation. Especially notable are recent observations of folding dynamics in proteins, yeast phosphoglycerate kinase (PGK) and ubiquitin mutant . These proteins fold according to a “downhill folding” scenario , meaning that the folding path between the unfolded and folded states is free of deep metastable minima, so the process of folding is “downhill” relaxation along this path. Nevertheless, the number of folded proteins, as a function of time, displays stretched-exponential behavior over a large time interval .

In this Letter we present a simple minimal model which is translationally invariant and without disorder or explicit frustration , which displays perfect stretched-exponential relaxation over wide time intervals. The stretched exponent $`\alpha `$ changes continuously from slow relaxation $`\alpha <1`$ to fast relaxation $`1\alpha 2`$, as parameters of the system vary. We identify dynamic self-organization as the origin of the stretched-exponential relaxation, and show, that slow and fast relaxation are caused by anti-cooperative and cooperative behavior, respectively. We provide a theoretical explanation of the discrete and continuous limiting cases. We also present a qualitative theory, which accounts for behavior of the stretched exponent in the intermediate range of parameters.

Real life systems \[4-9\] are much more complex than our simple model. However, we conjecture that at least one underlying reason for stretched-exponential relaxation is universal. Namely, self-organization, with fast and slow relaxation corresponding to cooperative and anti-cooperative behavior. We suggest that proteins belong to the anti-cooperative universality class. Our argument is, that initially the protein is loose and deforms easily. When folding into a particular local pattern occurs, this part becomes stiff, making folding of neighboring parts more difficult, and therefore behaving in an anti-cooperative way. Similar anti-cooperative behavior occurs in our model system for negative values of the coupling, as we discuss below.

Let us now specify the model. We consider a chain of nonlinear bistable elements. Each element is described by the order parameter $`u_n`$. The local energy $`E(u)`$ has two minima, one of which is a metastable state (local minimum of energy), the other one is the absolutely stable state (absolute minimum of energy). We assume overdamped dynamics of $`u_n`$ at the presence of delta-correlated thermal noise. Equations of motion of the system are:

$$\frac{u_n}{t}=\beta (u_{n+1}+u_{n1}2u_n)+F(u_n)+f_n^{\mathrm{stc}}.$$
(1)
Here $`\beta `$ represents linear coupling of neighboring sites, $`F(u)=dE/du=u(uu_0)(u1)`$ is a forcing term, corresponding to the bistable polynomial potential $`E(u)=(1/4)u^4(1/3)(u_0+1)u^3+(1/2)u_0u^2`$. The potential has two minima, $`u=0`$ and $`u=1`$, separated by a barrier, $`u=u_0`$. At $`u_0=1/2`$ the minima have equal energy. For $`u_0<1/2`$ the minimum $`u=1`$ becomes absolutely stable and $`u=0`$ metastable, with energy difference given by $`\mathrm{\Delta }E=1/121/6u_0`$. The stochastic term $`f^{\mathrm{stc}}`$ is given by a delta-correlated Langevin force

$$<f_n^{\mathrm{stc}}(t_1)f_m^{\mathrm{stc}}(t_2)>=T\delta _{mn}\delta (t_1t_2),$$
(2)
where $`T`$ defines the temperature in our system . Eq. (1) can be viewed as discrete one-dimensional Ginzburg-Landau equation with noise.

Let us first review equilibrium properties of our system. At $`T=0`$ the system is in the absolutely stable phase. When temperature is introduced, some (small) number of particles overcome the barrier because of equilibrium thermal fluctuations. The relative number of particles in the metastable phase is determined by Boltzmann statistics, and is exponentially small for temperatures much less than the energy difference between potential wells.

Our goal is to study non-equilibrium and non-stationary process of relaxation from the metastable phase to the absolutely stable phase. We assumes that all sites are initially in the metastable phase, and then introduces the thermal noise. Due to fluctuations, particles start to overcome the barrier, and the number of particles in the absolutely stable phase increases. We monitor relaxation by introducing the function $`n(t)=N(t)/N_0(T)`$, which describes the ratio between the concentration of particles in the absolutely stable phase at time $`t`$ and their equilibrium number $`N_0(T)`$. In the thermodynamic limit of infinite chain length, $`n(t)`$ is a well defined smooth function, which satisfies the conditions $`n(0)=0`$ and $`n(\mathrm{})=1`$.

We study relaxation numerically by integrating the Langevin equations (1) defined on long chain segments ($`N=200`$), and averaging the resulting functions $`n(t)`$ over many noise realizations. To improve numerical convergence, we use an implicit integration scheme. Details of the numerical method will be reported elsewhere .

A typical function $`n(t)`$ is given in Fig. 1. Important global features of the relaxation, valid in all ranges of parameters, are:

i) After a short transient time $`t_0`$ (Fig. 1), the system self-organizes and $`n(t)`$ starts to obey the stretched-exponential form $`n(t)=1\mathrm{exp}[(t/t_{rel})^\alpha ]`$. Note, that there are only two parameters in our fit, which are the stretched exponent $`\alpha `$ and the relaxation time $`t_{rel}`$. The precision of the stretched-exponential fit is extremely high, and is about $`0.1\%`$ for $`t>t_0`$;

ii) The stretched exponent $`\alpha `$ changes continuously, as temperature $`T`$ and coupling $`\beta `$ vary, see (Figs. 2 and 3). The observed relaxation is fast ($`\alpha >1`$) for positive values of $`\beta `$ and slow ($`\alpha <1`$) for negative values of $`\beta `$. As temperature $`T`$ is increased, the stretched exponent $`\alpha `$ approaches the Arrhenius law ($`\alpha =1`$), corresponding to the intrinsic frustration being overcome thermally;

iii) The stretched exponent $`\alpha `$ approaches the value $`\alpha =2`$ in the continuous limit of large $`\beta `$, representing a binary relaxation channel.

The dependence of $`\alpha `$ on $`\beta `$ is given in Fig. 2. At zero $`\beta `$ the system represents a set of uncoupled nonlinear sites. Each site is described by a one-dimensional Fokker-Planck equation. In this case relaxation is known to be described by an exponential (Arrhenius) law, with decay time determined by the lowest excited state of the Fokker-Planck operator . Therefore, for $`\beta =0`$ one has $`\alpha =1`$.

Let us now explain, why introducing positive/negative $`\beta `$ corresponds to cooperative/anti-cooperative behavior, respectively. If $`\beta `$ is positive, one may think of the inter-site coupling $`\beta `$ as an elastic chain connecting neighboring sites $`n`$ and $`n+1`$. If a particle overcomes the barrier, it attracts its neighbors and makes it, therefore, easier for them to jump over the barrier. Contrary to this, negative coupling corresponds to a repulsive force, which makes relaxation of nearest neighbors more difficult. (Indeed in the limit of large negative $`\beta `$, our model produces a strong staggered dimerization of the lattice). We conclude from Fig. 2 that the cooperative behavior leads to fast relaxation, $`\alpha >1`$, and anti-cooperative behavior leads to slow relaxation, $`\alpha <1`$.

We will now discuss the behavior of the stretched exponent $`\alpha `$ at large positive $`\beta `$. In this limit the system becomes continuous and is described by the continuous Ginzburg-Landau equation.

$$\frac{u}{t}=\beta \frac{^2u}{x^2}+F(u)+f(x,t)^{\mathrm{stc}}.$$
(3)
In the absence of thermal noise, topological excitations of this equation are fronts (kinks), which separate the absolutely stable and metastable phases. Due to the energy difference between the phases, the fronts propagate at finite velocity $`v`$, increasing the size of the absolutely stable phase. When temperature is introduced, local fluctuations of the order parameter give birth to kink-antikink pairs (Fig. 4). These pairs counter-propagate, replacing the metastable phase by the absolutely stable one. Let us now estimate the probability for the order parameter $`u(0,t)`$ at $`x=0`$ to stay in the metastable phase after time $`t`$. This probability is approximately equal to the probability $`P(t,l)`$, that no kink-antikink pair will be created by fluctuations during time $`t`$ at a distance $`lvt`$ from the origin. If such a pair is created, then the newly-born kink has enough time to reach $`x=0`$, before the time interval $`t`$ elapses, and annihilate the metastable phase; see Fig. 4. Since fluctuations are local, one can estimate $`P(l,t)`$ as $`\mathrm{exp}[lt/\eta ]`$, where $`\eta `$ is a constant. Therefore, the probability to stay in the metastable phase after time $`t`$ is approximately $`\mathrm{exp}[vt^2/\eta ]`$. Since the number of particles in the metastable phase after time $`t`$ is proportional to this probability, we conclude that in the continuous limit the stretched exponent $`\alpha =2`$. Our numerical data (see Fig. 2) are in a good agreement with this prediction.

The fact, that the stretched exponent $`\alpha =2`$ in the continuous case, is striking and has important consequences. Namely, space discreteness is necessary to observe continuous dependence of $`\alpha `$ on parameters. We can rephrase this statement as follows. It is known, that critical properties of equilibrium systems are described by continuous models, and these models can be organized into universality classes, such that each universality class has a set of universal critical exponents . We conjecture, that effective models for non-equilibrium, non-stationary self-organization, in particular relaxation, must include discreteness, in order to incorporate the continuous dependence of stretched exponents on parameters.

Let us now consider the nontrivial dependence of the stretched exponent $`\alpha `$ on the coupling $`\beta `$ and temperature $`T`$, for the intermediate values of $`\beta `$, given in Figs. 2 and 3 (a more detailed discussion of these effects will be given in ). In the absence of noise, the topological excitations of system (1) undergo a sequence of bifurcations, as described in . The most important of these bifurcations is the propagation failure bifurcation, at which fronts cease to propagate, being pinned by the lattice. We denote the corresponding value of the coupling $`\beta `$ as $`\beta _c`$. An analytical expression for $`\beta _c`$ was derived by us in . The second important bifurcation point is $`\beta _1`$. At this point, a one-site nucleus of the globally stable phase, representing a bound state of a kink and antikink pair, becomes unstable, and ”bursts” into counter-propagating kink and antikink. In  we have shown, that there exists an infinite number of bifurcation points $`\beta _n`$, such that $`\beta _c\mathrm{}\beta _n\beta _{n1}\mathrm{}\beta _1`$, corresponding to the instabilities of a one-site, two-site, three-site, etc. nuclei. In Fig 2 we see the corresponding change of shape in the coupling-dependence of the stretched exponent $`\alpha `$, in the interval $`\beta _c<\beta <\beta _1`$. Accordingly, Fig. 3 demonstrates a transition in the temperature-dependence of $`\alpha `$, with $`\beta `$ crossing $`\beta _c`$. We conclude from Figs. 2 and 3, that this region of bifurcations is manifested as a characteristic feature of the non-equilibrium relaxation of the system.

In conclusion, we have introduced a model of dynamic self-organization, and conjectured that it should be applicable to stretched-exponential relaxation in many biological and physical systems. We studied properties of self-organization in a simple model, and showed that fast/slow relaxation is related to cooperative/anti-cooperative behavior. We showed, that neither explicit disorder nor frustration are necessary requirements for anomalous relaxation laws. In view of this, an experimental evidence of anomalous relaxation in a physical system should not be considered as automatic proof of disorder or frustration in the system. Appropriate nonlinearity (arising, for example, from coupled degrees-of-freedom with feedback) can alone produce both the ingredients of local structures (”intrinsic disorder”) and competition (of length-scales), leading to a distribution of metastable states and anomalous dynamics . The detailed analytical description of dynamic self-organization is a promising direction for future study. This will include extensions to higher dimensions as well as specific physical models. We also note, that stretched exponents were also observed in equilibrium correlation functions , and it is important to understand, whether they can be related to their non-equilibrium counterparts.

We thank R. B. Laughlin for fruitful discussions and for sharing his views on self-organization in non-equilibrium processes. We are also grateful to J. Pearson, for his valuable comments. This work was supported in part (K.K.) by the Otto Hahn Fellowship of the Max Planck Society, Germany. Work at Los Alamos was supported by US DOE under contract No. W-7405-ENG-36.


format: text

text:
# Comparison of the Two Followup Observation Strategies for Gravitational Microlensing Planet Searches

## 1 Introduction

Experiments to detect microlensing events by monitoring millions of stars located in the Magellanic Clouds and the Galactic bulge have been and are being carried out by several groups (Alcock et al., 1993; Aubourg et al., 1993; Udalski et al., 1993; Alard & Guibert, 1995; Abe et al., 1997). With their efforts, the total number of candidate events now reaches up to a thousand (P. Popowski 1999, private communication).

Although the primary goal for these experiments was to investigate the nature of Galactic dark matter, it turns out to be that microlensing can be very useful in various other fields of astronomy. One of the important applications of microlensing is the detection of extra-solar planets. Planet detection by using microlensing is possible because the event caused by a lens system with a planet can exhibit detectable anomalies in the light curve when the source passes close to the lens caustics (Mao & Paczyński, 1991; Gould & Loeb, 1992; Bolatto & Falco, 1994). For this lens system, there are 2 or 3 disconnected sets of caustics depending on the projected separation between the planet and the primary lens (planet-lens separation). Among them, one is located very close to the primary lens, ‘central caustic’, and the other(s) is (are) located relatively away from the primary lens, ‘planetary caustics’ (Griest & Safizadeh, 1998). Accordingly, there exist two different types of anomalies in the light curves; one affected by the planetary caustic(s), ‘type I anomaly’, and the other affected by the central caustic, ‘type II anomaly’ (Covone et al., 2000). Due to the characteristics of the central caustic, type II anomalies occur near the peak of the light curves of high amplification events.

Compared to the frequency of type I anomalies, type II anomalies occur with a relatively low frequency due to the smaller size of the central caustic than the corresponding planetary caustic. However, the efficiency of detecting type II anomalies can be high because intensive monitoring of events is possible due to the predictable time of anomalies (Griest & Safizadeh, 1998). We call the planet search strategy of intensive monitoring near the peak of high amplification events as the ‘type II’ strategy, while the strategy of continuous monitoring for all events detected from lensing survey programs as the ‘type I’ strategy. Then a naturally arising questions is whether the planet detection rate (not the detection efficiency per event) under the type II strategy can make up a substantial fraction of the total rate. In this paper, we answer to this question by estimating the relative detection rates under the two strategies of planet search followup observations.

## 2 Types of Planet-induced Anomalies

The lens system with a planet is described by the formalism of the binary lens system with a very low-mass companion. When lengths are normalized to the combined angular Einstein ring radius, which is equivalent to the angular Einstein ring radius of a single lens with a mass equal to the total mass of the binary, the lens equation in complex notations for a binary lens system is represented by

$$\zeta =z+\frac{m_1}{\overline{z}_1\overline{z}}+\frac{m_2}{\overline{z}_2\overline{z}},$$
(1)
where $`m_1`$ and $`m_2`$ are the mass fractions of individual lenses (and thus $`m_1+m_2=1`$), $`z_1`$ and $`z_2`$ are the positions of the lenses, $`\zeta =\xi +i\eta `$ and $`z=x+iy`$ are the positions of the source and images, and $`\overline{z}`$ denotes the complex conjugate of $`z`$ (Witt, 1990). The combined angular Einstein ring radius is related to the physical parameters of the lens system by

$$\theta _\mathrm{E}=\left(\frac{4GM}{c^2}\frac{D_{\mathrm{ls}}}{D_{\mathrm{os}}D_{\mathrm{ol}}}\right)^{1/2},$$
(2)
where $`M`$ is the total mass of the binary lens system and $`D_{\mathrm{ol}}`$, $`D_{\mathrm{ls}}`$, and $`D_{\mathrm{os}}`$ represent the separations between the observer, lens, and source star, respectively. The amplification of each image is given by the Jacobian of the transformation (1) evaluated at the image position, i.e.

$$A_i=\left(\frac{1}{|\mathrm{det}J|}\right)_{z=z_i};\mathrm{det}J=1\frac{\zeta }{\overline{z}}\frac{\overline{\zeta }}{\overline{z}}.$$
(3)
Then the total amplification of a binary lens event is given by the sum of the amplifications of the individual images, i.e. $`A=_iA_i`$. The set of source positions with infinite amplifications, i.e. $`\mathrm{det}J=0`$, form closed curves called caustics. As a result, whenever a source passes very close to the caustic, the resulting light curve deviates significantly from the standard Paczyński (1986) curve.

The caustic structure of the lens system with a planet varies depending on the planet-lens separation $`d`$ normalized by $`\theta _\mathrm{E}`$. A lens system with a wide planet ($`d>1.0`$) forms 2 disconnected sets of caustics. One caustic is located very close to the center of mass (and thus referred as the central caustic), and the other planetary caustic is located away from the center of mass on the planet-side (with respect to the center of mass) axis connecting the primary lens and the planet. A system with a close planet ($`d<1.0`$) also has a single central caustic, but has two planetary caustics, which are located on the opposite side of the planet with respect to the center of mass and not on the planet-lens axis. The sizes of both the central and planetary caustics are maximized when the planet is located in the lensing zone of $`0.618d1.618`$ (Griest & Safizadeh, 1998), and decrease as the planet-lens separation becomes either smaller or larger. However, regardless of the separation the central caustic of a lens system is always smaller than the corresponding planetary caustic(s).

In Figure 1, we illustrate the central and planetary caustics of an example lens system with a planet and the corresponding type I and II anomalies in the resulting light curves. Presented in upper panel are the contour map of the fractional amplification excess, which is computed by

$$ϵ=\frac{\left|AA_\mathrm{s}\right|}{A_\mathrm{s}},$$
(4)
where $`A`$ and $`A_\mathrm{s}`$ represent the amplifications with and without the existence of the planet, respectively. The single lens event amplification is related to the separation (normalized by the angular Einstein ring radius of the primary lens) between the source and the primary lens, $`u`$, by

$$A_\mathrm{s}=\frac{u^2+2}{u(u^2+4)^{1/2}}.$$
(5)
The contours (thin solid curves) around the individual types of caustic (closed figures drawn by the thick solid line) are drawn at the levels of $`ϵ=1\%`$, 5%, and 10% with increasing distance from the caustics. The example lens system has a mass ratio and a projected planet-lens separation of $`q=10^3`$ (corresponding to a Jupiter-mass planet around a $`1M_{}`$ star) and $`d=1.2`$. The positions of the primary lens and the planet are chosen such that the center of mass is at the origin and all lengths are scaled by $`\theta _\mathrm{E}`$. The lower panels show the light curves resulting from the two source trajectories marked by straight lines in the upper panel. The two curves in each panel represent those expected in the presence (solid curve) and absence (dotted curve) of the planet. Times are in units of the Einstein ring radius crossing time (Einstein time-scale). From the figure, one finds that the region of large type II deviations is smaller than the corresponding region of type I deviations. One also finds that the type II anomaly in the light curve occurs near the peak amplification.

## 3 Type I versus Type II Strategies

Due to the difference in the characteristics of the central and planetary caustics and the resulting anomalies, detecting planets by using the type II strategy has both advantages and disadvantages. The greatest advantage of the type II strategy is that the time of anomalies, which occurs near the peak amplification, can be predicted in advance, and thus high time-resolution monitoring of the event is possible (Griest & Safizadeh, 1998). In addition, more accurate photometry can be performed because more photons will be available for high amplification events. On the other hand, the disadvantage is that type II anomalies occur with a relatively low frequency due to the smaller size of the central caustic. Then, the question is whether the planet detection rate under the type II strategy can be high enough to be comparable to the rate under the type I strategy. In this section, we answer to this question by estimating the relative planet detection rates expected under the two different types of planet search strategies.

To estimate the detection rates, we assume the following observational conditions and detection criteria. For the type I strategy, all events with $`A>A_{\mathrm{th}}=1.34`$ (i.e. $`uu_{\mathrm{th}}=1.0`$) are assumed to be monitored in a round-the-clock manner during 8 hours per night on average. Here $`u_{\mathrm{th}}`$ and $`A_{\mathrm{th}}`$ represent the threshold lens-source separation and the corresponding threshold amplification, which are required for the event to become the target of followup observations. On the other hand, if the events are suspected to have high amplifications with $`AA_{\mathrm{th}}=10`$ (i.e. $`u<u_{\mathrm{th}}=0.1`$), they become targets for intensive montoring by the type II strategy. We assume that the type I strategy is employed during $`0.4t_\mathrm{E}t_{\mathrm{obs}}2.0t_\mathrm{E}`$, while the type II strategy is employed only during the time of high amplifications with $`AA_{\mathrm{th}}`$. For planet detection, it is assumed that one should detect light curve deviations greater than a threshold value, $`ϵ_{\mathrm{th}}`$, at least more than 5 times. The adopted threshold deviations are $`ϵ_{\mathrm{th}}=5\%`$ for the type I and $`1\%`$ for the type II strategies, respectively. We adopt the lower value of $`ϵ_{\mathrm{th}}`$ for the type II strategy because the events monitored by this strategy will have higher amplifications. For more sophisticated computations of the planet detection probability including not only more refined models of photometric precision but also the effects of finite source size and blending, see Gaudi & Sackett (2000). According to the current microlensing detection rate, there are on average several dozens of on-going events. Then, considering the average total observation time of $`8`$ hours per night and $`5`$ – 10 minutes of required time per event, one can observe on average $`2`$ – 3 times for each event under the type I strategy. However, since the frequency can be increased by employing multiple telescopes like the current followup observation teams (Albrow et al., 1998; Rhie et al., 1999), we test various monitoring frequencies per night of $`f_{\mathrm{obs}}=2`$, 5, and 10 for the type I strategy. For the type II strategy, we test $`f_{\mathrm{obs}}=10`$, 20, and 50.

In Figure 2, we present the distributions of the planet detection probability expected under the individual strategies with various monitoring frequencies as a function of the planet-lens separation for events caused by a planetary lens system with $`q=10^3`$. The events are assumed to have a common Einstein time scale of $`t_\mathrm{E}=20\mathrm{days}`$, which corresponds to that of the bulge self-lensing event with a lens mass $`1M_{}`$. The region enclosed by the dot-dashed lines represents the lensing zone. Note that our probability under the type II strategy are normalized by the total number of events with $`A1.34`$ unlike those of Griest & Safizadeh (1998), who normalized by the number of only the intensively monitored high amplification events. Since we applied $`u_{\mathrm{th}}=0.1`$, our probability is equivalent to approximately one tenth of their value. In Table 1, we list the fractions of planet detection rate by the type II strategy out of the total rate (type I + type II) for various combinations of $`f_{\mathrm{obs}}`$ of the individual strategies. We estimate two sets of fractions. The first set is for all planetary systems with separations $`d5.0`$ and the second set is for those with planets located within the lensing zone only.

The findings from Figure 2 and Table 1 are summarized as follows:

1. As pointed out by Griest & Safizadeh (1998), the planet detection efficiency per event, which is ten times of the probability presented in Figure 2, expected under the type II strategy is very high. We note, however, that our estimate is somewhat lower than their estimate. We suspect the difference is caused because they assume that the type II monitoring strategy is applied during the entire duration of the event, while we assume that intensive monitoring is performed only during the time of high amplifications. We actually obtain similar efficiency distributions to theirs under the observational conditions they assumed.
2. In addition to the higher efficiency, the planet detection rate under the type II strategy is expected to comprise a substantial fraction of the total rate. We find that the fraction will be $`1/4`$ – 1/2 depending on the combinations of the individual strategies’ monitoring frequencies.
3. Compared to the type I strategy, the type II strategy is more efficient in detecting planets located outside of the lensing zone. This is because while the planetary caustic lies at the position with increasing distance from the primary lens as the planet-lens separation becomes smaller or large than $`d1.0`$ (Wambsganss, 1997), the central caustic is located always near the primary lens regardless of the planet-lens separation (see Figure 3). As a result, while the type II efficiency decreases only by the decrease of the (central) caustic size, the type I efficiency is decreased additionally by the increasing distance of the (planetary) caustic location from the primary lens.
4. The optimal monitoring frequency of the type II strategy will be $`f_{\mathrm{obs}}20`$. This is because observations with frequencies higher than this value will result in very slight increase in the planet detection rate (see Table 1). This optimal frequency can be easily achieved by the current followup programs even with a single telescope.

## 4 Conclusion

We compare the planet detection rates expected by the type I and II followup observation strategies. From this comparison, we find that the rate under the type II strategy will be substantial, comprising $`1/4`$ – 1/2 of the total rate. It is also found that compared to the type I strategy, the type II strategy is more efficient in detecting planets located outside of the lensing zone. Therefore, in addition to yielding high planet detection efficiency for each high amplification event, employment of the type II strategy will allow one not only to increase the total planet detection rate but also to probe the existence of planets with a wider range of separations. Despite these benefits, implementation of the strategy requires only a modest improvement in the monitoring frequencies.

This work was supported by the grant from Korea Science & Engineering Foundation (KOSEF).


text:
# Theory of Ideal Four-Wave Mixing in Bose-Einstein Condensates

## 1  Introduction

The experimental realizations of Bose-Einstein condensation (BEC) in trapped dilute alkali atomic vapours  and the atom “laser” have made possible the experimental study of four-wave mixing (4WM) of coherent matter-waves. As a typical effect of nonlinear atom optics, the possibility of 4WM with matter waves of Bose-Einstein condensate (BEC) can be deduced from the similarity between the Maxwell equations for lasing action in nonlinear optical media and the Gross-Pitaevskii(GP) equation widely used in studying the properties of macroscopic quantum systems such as BECs.

In the past few years, 4WM in BECs and related problems have been studied intensively . Early in 1995, Goldstein et al. had started to investigate the problem of phase conjugation in BECs, including the case of multi-component BECs. A 4WM experiment using three colliding BEC wavepackets with different central momenta was proposed by Trippenbach et al.. Recently, Deng et al. successfully carried out 4WM experiment in sodium condensates by means of Bragg diffraction technology. In addition, stimulated by the work of Ho and coworkers and Ohmi and Machida on the Bose condensates with internal degrees of freedom, Law et al. studied quantum spin-mixing in spinor BECs using an algebraic method developed in problems related to cavity QED, and Goldstein and Meystre developed a quantum theory of atomic 4WM involving both the internal and spatial degrees of freedom. Rza̧żeewski et al. investigated the fluctuations in the populations of atoms in 4WM wavepackets due to the quantum mechanical nature of the mean-field BEC wavefunctions. Although the experimental results turn out to be consistent with the numerical simulations by Trippenbach et al. based on the GP equation, further theoretical studies are still needed in order to understand 4WM in BECs more clearly. Without invoking the undepleted pump approximation, Wu et al. calculated analytically the 4WM in the experimental condition of Ref.. Based on the slowly-varying-envelope approximation, Trippenbach et al. developed a three-dimensional quantum mechanical description for 4WM of wavepackets created from BEC.

The main purpose of the present work is to study the ideal 4WM with coherent matter waves from BEC in which only four plane waves with definite wavevectors and frequencies are involved. We organize the paper as follows. A set of time evolution equations for 4WM are derived in Sec.2 starting from a Hamiltonian for many-particle systems in the second quantized form rather than the original GP equation. Section 3 gives a short discussion on the phase-matching condition for ideal 4WM. In Sec.4, we present results demonstrating the effects of relative phase difference and the collapse and revival (C&R) behaviour of ideal 4WM. The problem of conversion efficiency and the main characters of recent 4WM experiment are investigated in Sec.5. Sec.6 gives a brief summary.

## 2 Time evolution equations for 4WM

To study 4WM with matter waves, we consider the dynamics of $`N`$ interacting identical bosons. For simplicity, we omit the internal degrees of freedom. In standard experimental setups for studying 4WM, the confining potential is removed, and the Hamiltonian can be generally written in terms of bosonic creation and annihilation operators as

$$=\underset{k}{}\epsilon _k\widehat{a}_k^{}\widehat{a}_k+\underset{𝐤_4+𝐤_3=𝐤_1+𝐤_2}{}\frac{U_q}{2𝒱}\widehat{a}_{k_4}^{}\widehat{a}_{k_3}^{}\widehat{a}_{k_2}\widehat{a}_{k_1},$$
(1)
where $`\epsilon _k=\mathrm{}^2k^2/2M`$ is the kinetic energy of a particle in a single-particle state of momentum $`𝐩=\mathrm{}𝐤`$, $`𝐪𝐤_4𝐤_1=𝐤_2𝐤_3`$ denotes the momentum transfer in two-body scattering, and $`M`$ denotes the atomic mass. Momentum conservation in the second term results from the calculation of the interaction matrix elements $`U_{𝐤_4𝐤_3𝐤_2𝐤_1}=U_q\delta _{𝐤_4+𝐤_3,𝐤_1+𝐤_2}`$, where $`U_qd^3𝐫U(r)\mathrm{exp}(i𝐪𝐫)`$ is the Fourier transform of the interaction potential $`U(r)`$. The field operator can be expanded as $`\widehat{\mathrm{\Psi }}(𝐫,t)=_k\widehat{a}_k\mathrm{exp}(i𝐤𝐫)/\sqrt{𝒱}`$ with $`𝒱`$ denoting the volume of the system.

Within the Heisenberg picture, the wavefunction is independent of time and the dynamical evolution of the system is given by that of the field operators. The annihilation operators satisfy equations of the form

$$i\mathrm{}\frac{\widehat{a}_{k_4}}{t}=(\epsilon _{k_4}+U_0\widehat{n})\widehat{a}_{k_4}+\underset{\stackrel{𝐤_1+𝐤_2𝐤_3=𝐤_4}{q0}}{}\frac{U_q}{𝒱}\widehat{a}_{k_3}^{}\widehat{a}_{k_2}\widehat{a}_{k_1},$$
(2)
where the subscripts $`k_1,k_2,k_3`$ and $`k_4`$ can be any allowed values of wavevectors (or particle momenta), and the operator of the averaged number density of atomic gas is a fixed constant defined by $`\widehat{n}(1/𝒱)_k\widehat{a}_k^{}\widehat{a}_k`$.

Guided by results obtained through the experimental realization of BEC and the technique of Bragg diffraction, we assume that the creation and annihilation operators $`\widehat{a}_k^{}`$ and $`\widehat{a}_k`$ of a certain mode with macroscopic population of atoms can be treated as a pair of conjugated complex numbers $`a_k^{}`$ and $`a_k`$, and those for the other modes vanished. Therefore the field operator becomes a macroscopic wavefunction satisfying the GP equation, except that $`U_q`$ may depend on $`q`$ for scatterings with large momentum transfer. Introducing a new variable $`A_ka_k/\sqrt{𝒱}`$ and assuming for simplicity a $`q`$-independent $`U_q4\pi a_0\mathrm{}^2/M`$ with $`a_0`$ denoting the s-wave scatteing length, the slowly varying envelope function of the field, $`\stackrel{~}{A}_kA_k\mathrm{exp}[i(\epsilon _knU_0)t/\mathrm{}]`$ with $`n=_ka_k^{}a_k/𝒱`$, satisfies a set of coupled nonlinear equations

$$i\frac{\stackrel{~}{A}_{k_4}}{t}=\frac{2a_0h}{M}\underset{\stackrel{𝐤_1+𝐤_2𝐤_3=𝐤_4}{q0}}{}\stackrel{~}{A}_{k_3}^{}\stackrel{~}{A}_{k_2}\stackrel{~}{A}_{k_1}e^{i\mathrm{\Delta }\omega t}$$
(3)
with $`\mathrm{\Delta }\omega =(\epsilon _{k_4}+\epsilon _{k_3}\epsilon _{k_2}\epsilon _{k_1})/\mathrm{}`$. The momentum conservation indicated in the summation is a consequence of the translational invariance of the system. In contrast, collision processes with $`\mathrm{\Delta }\omega 0`$ are not forbidden in principle. As the system has a fixed number of atoms, the averaged number density is a constant, i.e., $`_k|\stackrel{~}{A}_k|^2=n`$. We can then define the normalized field amplitudes as $`\xi _k=\stackrel{~}{A}_k/\sqrt{n}`$, which satisfy the normalization condition $`_k|\xi _k|^2=1`$, to describe the relative distribution of atoms more conveniently. The time evolution equations for these normalized amplitues take on the more compact form of

$$i\frac{\xi _{k_4}}{\tau }=\frac{1}{2}\underset{\stackrel{𝐤_1+𝐤_2𝐤_3=𝐤_4}{q0}}{}\xi _{k_3}^{}\xi _{k_2}\xi _{k_1}e^{i\mathrm{\Omega }\tau },$$
(4)
where the re-scaled time $`\tau t/T`$ is dimensionless with the characteristic time $`TM/4na_0h`$ being inversely proportional to the number density, and the phase mismatch or ‘energy-loss’ parameter is defined as $`\mathrm{\Omega }T\mathrm{\Delta }\omega `$. Before going into details, some important conclusions can be drawn on the common properties of 4WM with matter waves. The number density of atoms can affect the time evolution of 4WM only by modifying the characteristic time $`T`$, implying that 4WM evolves faster in situations of higher number density of atoms. In actual experimental situations, the characteristic time is estimated to be about $`3ms`$ (see Sec.5).

The general form of time evolution equations for multi-wave mixing (MWM) is a rather complicated set of coupled nonlinear equations because the subscripts $`k_1,k_2,k_3`$ and $`k_4`$ in Eq.(4) can take on any allowed values of wavevectors (or particle momenta). However it can be simplified for an ideal 4WM process to

$$\frac{\xi _1}{\tau }=i\xi _2^{}\xi _3\xi _4e^{i\mathrm{\Omega }\tau },\frac{\xi _2}{\tau }=i\xi _1^{}\xi _3\xi _4e^{i\mathrm{\Omega }\tau },$$
(5)
$$\frac{\xi _3}{\tau }=i\xi _4^{}\xi _2\xi _1e^{i\mathrm{\Omega }\tau },\frac{\xi _4}{\tau }=i\xi _3^{}\xi _2\xi _1e^{i\mathrm{\Omega }\tau }.$$
(6)
This set of equations is analogous to those in optical 4WM except that the spatial variable is replaced by the re-scaled time. In what follows, we will study 4WM of matter waves based on the simplified model described by Eqs.(5) and (6). All numerical calculations are carried out using the 4th order Runge-Kutta method.

## 3 Phase-matching condition for 4WM

The recent 4WM experiment shows that three colliding wavepackets with arbitrary central momenta cannot always create a fourth wavepacket, and a possible 4WM process must satisfy the phase-matching condition, i.e. momentum-energy conservation. The Hamiltonian (Eq.(1)) describes a system satisfying the conservation of total energy, momentum and particle-number. Therefore the phase-matching condition is always understood as $`𝐤_4+𝐤_3=𝐤_1+𝐤_2`$ and $`\epsilon _{k_4}+\epsilon _{k_3}=\epsilon _{k_2}+\epsilon _{k_1}`$. Note that from the dispersion relation of massive particles $`\epsilon _k=(\mathrm{}k)^2/2M`$, the phase mismatch parameter can be written as $`\mathrm{\Omega }=(𝐤_3𝐤_2)(𝐤_3𝐤_1)/8\pi na_0`$. The phase-matching condition, i.e. $`\mathrm{\Omega }=0`$, for matter waves can be shown graphically as in Fig.1.

It has also been emphasized, however, in the discussion after Eq.(1) that the scattering process of bosons does not require the conservation of energy. So the constraint of energy conservation deserves some more explanations. In fact, numerical calculations show that the parameter $`\mathrm{\Omega }T\mathrm{\Delta }\omega `$ affects 4WM processes, but only large enough phase mismatch can suppress the creation of fourth wave. Figure 2 shows the population $`N_4`$ of the created wave as a function of $`\tau `$ for different values of $`\mathrm{\Omega }`$. For small $`\mathrm{\Omega }`$, $`N_4(\tau )`$ is similar to that of an ideal 4WM process satisfying the phase-matching condition exactly. The conversion efficiency is lower for higher values of $`\mathrm{\Omega }`$, and no 4WM signal can be observed when $`|\mathrm{\Omega }|\pi `$.

## 4  Dynamics of Ideal 4WM

The time evolution of an ideal 4WM system can be described by Eqs.(5) and (6). This kind of process satisfies not only particle-number conservation but also the conditions

$$\frac{|\xi _4|^2}{\tau }=\frac{|\xi _3|^2}{\tau }=\frac{|\xi _2|^2}{\tau }=\frac{|\xi _1|^2}{\tau },$$
(7)
implying that the two pumping waves must provide atoms to amplify the probing and created waves at the same rate. When one of the two pumping waves is exhausted, the another one will take on the role of the probing wave. Figure 3 shows the calculated results for typical collapse and revival (C&R) behaviour of 4WM. Initially, the $`N`$ atoms are distributed in three states corresponding to the pumping wave 1, pumping wave 2 and probing wave 3. Due to two-body interaction and the bosonic enhancement effect, the atoms in the pumping waves begin to transfer into the states 3 and 4. After a short period of time, all the atoms are transferred out of state 2 and no more atoms are transferred into state 4. The second phrase of 4WM begins at around $`\tau 5`$. Now waves 3 and 4 play the roles of the pumping waves, and wave 1 acts as a probing wave, wave 2 will become appreciable again until the pumping wave 4 vanishes at about $`\tau 10.5`$. This completes one period of ideal 4WM.

Another special dynamical character in 4WM processes is the effect of relative phase difference. It is not the phase of each individual wave that matters. Instead if we write $`\xi _i=\sqrt{I_i}e^{i\phi _i}`$, it is the combination of the phases $`\alpha =\phi _4+\phi _3\phi _2\phi _1\mathrm{\Omega }\tau `$ that affects 4WM in the way that

$$\frac{I_1}{\tau }=\frac{I_2}{\tau }=\frac{I_3}{\tau }=\frac{I_4}{\tau }=2𝒜\mathrm{sin}\alpha ;$$
(8)
$$\frac{\alpha }{\tau }=\mathrm{\Omega }+𝒜(\frac{1}{I_1}+\frac{1}{I_2}\frac{1}{I_3}\frac{1}{I_4})\mathrm{cos}\alpha ,$$
(9)
where the amplification factor is defined as $`𝒜\sqrt{I_1I_2I_3I_4}`$. When one of the four waves vanishes, we can set the combined phase $`\alpha `$ to be an arbitrary value. Therefore the initial values of the pumping and probing waves do not affect the time evolution of atomic distributions.

## 5 Conversion efficiency of 4WM

All of the previous works on the 4WM of matter waves have focussed on the initial rate of growth of the fourth wave because the time of wavepackets interaction in experiment is very short (determined by the size and individual momentum of wavepackets). Initially the rate of growth of the fourth wave is proportional to the product of the intensities of pumping and probing waves. If the experimental setup can be changed in such a way that 4WM can proceed for a much longer period, we need more precise knowledge about the conversion efficiency of 4WM.

In the experiment  of 4WM with matter waves in BEC with $`N1.7\times 10^6`$ sodium atoms ($`M=3.848\times 10^{23}g`$ and $`a_0=28\AA `$), the pumping (1,2) and probing (3) waves are prepared by Bragg scattering. It was estimated that the population numbers in the three states took on the initial values of $`N_1^06.8\times 10^5`$, $`N_2^06.9\times 10^5`$ and $`N_3^03.3\times 10^5`$, respectively. Assuming that 4WM took place in a region of volume $`𝒱=L^3`$ with $`L0.1mm`$, the characteristic time will be $`TML^3/4Na_0h3ms`$.

The four curves in Fig. 4 show typical results of the conversion efficiency $`N_4/N`$ as a function of the re-scaled time $`\tau `$ obtained by Eqs.(5) and (6), hence assuming an ideal 4WM process with $`\mathrm{\Omega }=0`$, for different sets of initial conditions. Curve (i) shows the situation in which $`N_1^0=N_2^0N_3^0`$, which achieves a maximum efficiency of $`50\%`$ if 4WM can proceed for a time long enough, e.g., $`t>10T`$. Curve (ii) is calculated using a set of initial conditions similar to the experimental situation. Comparing with curve (i), it demonstrates a much greater initial rate of growth of the fourth wave and a, however, smaller saturation efficiency. Note that a conversion efficiency of $`10\%`$ is achieved at $`t1.9T5.8ms`$, which is consistent with the results reported in Ref.. The initial conditions for curves (iii) and (iv) are chosen so that both of them have the same initial rate of growth of the fourth wave. Even so, the two curves demonstrate very different evolution processes as a function of $`\tau `$. Therefore, carefully prepared initial atomic distributions are crucial in achieving a higher initial rate of increase in the efficiency and at the same time a high enough saturation efficiency.

While the initial conditions in Ref. are close to those values in achieving a high efficiency, the observed value of $`610\%`$ is still rather low. A higher conversion efficiency can then be achieved by increasing the interaction time and/or by shortening the characteristic time by increasing the number density of atoms.

## 6 Summary

In conclusion, a set of nonlinear equations for the time evolution of multi-wave mixing in BECs is dervied. The formalism leads to a physically transparent picture of 4WM with matter waves. The ideal 4WM model is successfully applied to analyse results of 4WM experiments. Factors affecting the value of the conversion efficiency are discussed. The GP equation in momentum space provides us with an efficient way to study the wavepacket effects of 4WM because the four wavepackets are limited to four spots of the momentum space even though they may spread widely in real space.

## 7 Acknowledgments

This research is supported in part by the NSF of China for “Climbing Project in Theoretical Physics”, the Guangdong Provincial Natural Science Foundation of China, and a Direct Grant for Research at the Chinese University of Hong Kong.

FIGURE CAPTIONS

* Fig.1 The momentum sphere of phase-matching condition for 4WM. The center is located at $`(𝐤_1+𝐤_2)/2`$, and the radius of the sphere is defined as the magnitude of the vector $`(𝐤_1𝐤_2)/2`$. If the momentum of the probing wave points exactly at a point on the surface of the sphere, the largest conversion efficiency of 4WM is achieved.
* Fig.2 The conversion coefficiency as a function of interaction time for five different values of the phase-mismatch paramter $`\mathrm{\Omega }=T\mathrm{\Delta }\omega `$. The time is given in units of the characteristic time as defined in the text. The five curves from top to bottom correspond to $`\mathrm{\Omega }/\pi =0,\mathrm{\hspace{0.33em}0.1},\mathrm{\hspace{0.33em}0.2},\mathrm{\hspace{0.33em}0.4},\mathrm{\hspace{0.33em}1}`$, respectively. The initial conditions for the pumping and probing waves are chosen to be $`\xi _1^0,\xi _2^0,\xi _3^0=0.68,\mathrm{\hspace{0.33em}0.58},\mathrm{\hspace{0.33em}0.45}`$.
* Fig.3 Typical time evolution curves of the populations of the matter waves involved in a 4WM process. The initial conditions are chosen to be $`\xi _1^0,\xi _2^0,\xi _3^0=0.68,\mathrm{\hspace{0.33em}0.58},\mathrm{\hspace{0.33em}0.45}`$.
* Fig.4 Typical results of the conversion efficiency as a function of the interaction time for different initial conditions. The initial conditions for curve (i) are $`\xi _1^0,\xi _2^0,\xi _3^0=0.707,\mathrm{\hspace{0.33em}0.703},\mathrm{\hspace{0.33em}0.243}`$, which lead to a maximum efficiency of $`50\%`$ but with a low initial increasing rate. Curve (ii) has initial conditions similar to those in experimental situation: $`\xi _1^0,\xi _2^0,\xi _3^0=0.637,\mathrm{\hspace{0.33em}0.632},\mathrm{\hspace{0.33em}0.441}`$. Curves (iii) and (iv) show that in order to achieve an appreciable efficiency, the probing wave can be weak but the pumping waves must be strong. The initial conditions are: (iii) $`\xi _1^0,\xi _2^0,\xi _3^0=0.804,\mathrm{\hspace{0.33em}0.553},\mathrm{\hspace{0.33em}0.217}`$ and (iv) $`\xi _1^0,\xi _2^0,\xi _3^0=0.804,\mathrm{\hspace{0.33em}0.217},\mathrm{\hspace{0.33em}0.553}`$.


text:
# Exchange Frequencies in the 2d Wigner crystal

## Abstract

Using Path Integral Monte Carlo we have calculated exchange frequencies as electrons undergo ring exchanges in a “clean” 2d Wigner crystal as a function of density. The results show agreement with WKB calculations at very low density, but show a more rapid increase with density near melting. Remarkably, the exchange Hamiltonian closely resembles the measured exchanges in 2d 3He. Using the resulting multi-spin exchange model we find the spin Hamiltonian for $`r_s175\pm 10`$ is a frustrated antiferromagnetic; its likely ground state is a spin liquid. For lower density the ground state will be ferromagnetic.

PACS Numbers: 73.20Dx, 75.10-b, 67.8Jd

The uniform system of electrons is one of the basic models of condensed matter physics. In this paper, we report on the first exact calculations of the spin Hamiltonian in the low density 2 dimensional Wigner crystal (2dWC) near melting. This system is realized experimentally with electrons confined at an semiconductor MOSFET’s and heterostructures, and for electrons on the surface of liquid helium.

A homogeneous charged system is characterized by two parameters: the density given in terms of $`r_s=a/a_0=(m^{}/mϵ)(\pi a_0^2\rho )^{1/2}`$ and the energy in effective Rydbergs $`Ry^{}=(m^{}/m_eϵ^2)Ry`$ where $`m^{}`$ is the effective mass and $`ϵ`$ the dielectric constant. Figure 1 summarizes the 2d phase diagram. At low density (large $`r_s`$) the potential energy dominates over the kinetic energy and the system forms a perfect triangular lattice, the Wigner crystal. Tanatar and Ceperley determined that melting at zero temperature occurs at $`r_s37\pm 5`$. Recent calculations have shown that the low temperature phase is free of point defects for densities with $`r_s50`$ but defects may be present very near melting. At densities for $`r_s100`$ the melting is classical, and occurs for temperatures $`T_{melt}=2Ry/(\mathrm{\Gamma }_cr_s)`$ where $`\mathrm{\Gamma }_c137`$.

We determine the spin-spin interaction in the Wigner crystal, using Thouless’ theory of exchange. According to this theory, in the absence of point defects, at low temperatures the spins will be governed by a Hamiltonian of the form:

$$_{spin}=\underset{P}{}(1)^PJ_P\widehat{𝒫}_{spin}$$
(1)
where the sum is over all cyclic (ring) exchanges described by a cyclic permutation $`P`$, $`J_P`$ is its exchange frequency and $`\widehat{𝒫}_{spin}`$ is the corresponding spin exchange operator. Path Integral Monte Carlo (PIMC) as suggested by Thouless and Rogerhas proved to be the only reliable way to calculate these parameters. The theory and computational method have been tested thoroughly on the magnetic properties of bulk helium obtaining excellent agreement with measured properties. Rather surprisingly, it has been found that in both 2d and 3d solid 3He, exchanges of 2, 3 and 4 particles have roughly the same order of magnitude and must all be taken into account. This is known as the multiple spin exchange model(MSE).

A WKB calculation of the exchange frequencies in the 2dWC by Roger predicted that the three electron $`J_3`$ nearest neighbor exchange would dominate, leading to a ferromagnetic(F) ground state. It has been recently pointed out that if the 2dWC can be stabilized by disorder to higher densities, a transition to an antiferromagnetic(AF) phase will occur. Although the 1/r interaction is characterized as ’soft’, in the 2dWC the two particle pair correlation function for the electron system is quite similar to that of solid 3He supporting the idea that multiple exchanges could be important in the 2dWC near melting.

The method of determining the exchange frequency in quantum crystals has been given earlier. One computes the free energy necessary to make an exchange beginning with one arrangement of particles to lattice sites $`Z`$ and ending on a permuted arrangement $`PZ`$:

$$F_P(\beta )=Q(P,\beta )/Q(I,\beta )=\mathrm{tanh}(J_P(\beta \beta _0)).$$
(2)
Here $`Q(P,\beta )`$ is the partition function corresponding to an exchange $`P`$ at a temperature $`1/\beta `$. $`I`$ is the identity permutation. Note that these paths are of “distinguishable” particles since Fermi statistics are implemented through the spin Hamilitonian in Eq. (1). We determine the function $`F_P(\beta )`$ using a method which directly calculates free energy differences and thereby determines $`J_P`$ and $`\beta _0`$. The only new feature with respect to the 3He calculations is the method of treating the long-range potential. We use the standard Ewald breakup and treat the short-range part using the exact pair action and the long-range k-space term using the primitive approximation. We used a hexagonal unit cell with periodic boundary conditions, most calculations containing 36 electrons. Checks with up to 144 electrons did not change the results within the statistical errors. The number of particles is not as important in the 2dWC as in solid 3He because the $`1/r`$ interaction suppresses the long wavelength charge fluctuations.

We performed calculations for the densities $`r_s=\{45,50,60,75,100,140,200\}`$. We found accurate results using a “time-step” for the descretized imaginary time path integrals of $`\tau 0.3r_s^{3/2}`$ and extrapolated to the $`\tau =0`$ limit using $`J_p(\tau )=J_p(0)+J_p^{}\tau ^3`$. The inverse temperature $`\beta `$ in Eq. (2 )must be larger than the values determined by melting and twice the exchange “time” $`\beta _05r_s^{3/2}`$. We have determined the exchange frequencies with an accuracy between 1.5% and 6% independent of their magnitude or the number of exchanging electrons though the computer time increases with the number of exchanging electrons. Breakdown of Thouless’ theory caused by many states contributing to the ratio of partition functions would be signaled by $`F_P(\beta )`$ not described by Eq. (2). Except for $`r_s<50`$ where our calculations are too unstable to make definite predictions, we observed no problems of convergence.

Figure 2 shows the ring exchanges considered here. Except near melting these exchanges give rise to most of the thermodynamic properties. Note that we consider the 6-particle parallelogram (6p) exchange, which is not taken into account in solid 3He. We have also calculated several 2-5 particle exchanges having next-nearest neighbor exchanges and all possible 6 particle nearest-neighbor exchanges for $`50r_s75`$, but because their magnitudes are much smaller, we do not report those results.

Our calculated exchange frequencies vary rapidly with density as shown in fig. (3). One can see that they are much less than the zero point energy of the electrons, thus justifying the use of Thouless’ theory. The WKB method where one approximates the Path Integral in Eq. (2) by the single most probable path, explains most of this density dependence. In the 2dWC, the WKB expression for the exchange frequency is:

$$J_P=A_P(r_s)b_P^{1/2}r_s^{5/4}e^{b_Pr_s^{1/2}}.$$
(3)
Here $`b_Pr_s^{1/2}`$ is the minimum value of the action integral along the path connecting $`PZ`$ with $`Z`$. The 3 particle exchange exponent is the smallest indicating that as $`r_s\mathrm{}`$, $`J_3`$ will dominate and the system will have a ferromagnetic ground state. However, note that in fig. (3) $`J_2>J_3`$ for $`r_s90`$. We observe a more rapid increase in all the exchange frequencies than predicted by Eq. (3) for $`r_s<100`$ caused by non-linear fluctuations about the classical path. Near melting the spread in exchange frequencies is much less than in the low density limit.

Figure (4) shows density dependence of the prefactor, $`A_P`$. The prefactors are larger than unity as expected (our definition of $`A_p`$ differs from ref.  by a factor of $`1.63/\sqrt{2\pi }`$). We observe that the prefactor is roughly constant for $`r_s>100`$ but rises rapidly as the system approaches the melting density. We are unaware of any work that definitively establishes the form of the correction to the WKB formula. A good fit was obtained with the function $`A_P(r_s)=A_P(0)(1+(r_P/r_s)^3)`$ to determine $`b_P`$, $`A_P(0)`$ and $`r_P`$. The exponents are close to the WKB values, only differing significantly for the 2-particle exchange.

The prefactor for $`J_2`$ is twice as large as the other exchanges. One might think that since one cannot observe the difference in pair exchanges in the initial and final spins, the sense of the exchange would not matter. Our results show that the 2 orientations of 2 particle exchange, clockwise and counter-clockwise, should be counted separately. A recent calculation directly evaluates $`J_2`$ in the high density regime, $`5r_s45`$, by numerically solving for the difference between the even and odd parity 2-electron energies. In those calculations, the spectator electrons were fixed at their lattice sites and the 2 exchanging electrons had vanishing wavefunctions outside a rectangle centered around the exchange. At $`r_s=45`$ a direct comparison shows their result is 2.3 times larger than that obtained with PIMC. This result is expected because of the additional localization caused by the spectator electrons fluctuating into the exchanging region.

Having determined the exchange frequency, one is left with the spin Hamiltonian of Eq. (1). This is a non-trivial many-body problem which we will not discuss in detail here. For spin 1/2 systems, $`J_2`$ and $`J_3`$ contribute only with a nearest neighbor Heisenberg term: $`J_2^{\mathrm{eff}}=J_22J_3`$. This term is negative (ferromagnetic) but approaches zero near melting. For convenience we use $`J_4`$ as a reference to fix the overall scale of the magnetic energy. Neglecting $`J_{6p}`$, the Hamiltonian has 3 remaining parameters $`J_2^{\mathrm{eff}}/J_4`$, $`J_5/J_4`$ and $`J_{6h}/J_4`$. The dependence of these ratios on density is shown in fig. (5).

High temperature series expansions  determine the the specific heat $`C_V`$ and magnetic susceptibility $`\chi _0/\chi `$ for temperatures $`k_BTJ_P`$. The susceptibility is given by: $`\chi _0/\chi =T\theta +B/T\mathrm{}`$ and the specific heat $`C_V/Nk_B=(3J_c/2T)^2+\mathrm{}`$ where the Curie-Weiss constant is given by $`\theta =3(J_2^{\mathrm{eff}}+3J_45J_5+5/8J_{6h}+15/4J_{6p})`$ with a quadratic expression of the $`J`$’s for $`J_c`$. These two constants which set the scale of the temperature where exchange is important are shown as a dotted and dashed lines on fig. (1). Both $`\theta `$ and $`J_c`$ decrease very rapidly at low density showing that experiments must be done at $`r_s100`$ if spin effects are to be at a reasonable temperature, e.g. $`T_c>0.1mK`$ (assuming $`ϵ=m^{}=1`$). Note that $`\theta `$ changes from positive to negative at $`r_s130`$.

The zero temperature state can be studied by exact diagonalization (ED) of an $`N`$ site system. (The present limitation is $`N36`$.) Here we focus on the zero temperature transition from a ferromagnetic to an antiferromagnetic state. An estimate  of the critical transition parameters is shown on fig.(5). The ferromagnetic phase is obtained only at very low density: we estimate the the F-AF transition for the 2dWC will occur at $`r_s=175\pm 10`$. (Note that this estimate does not include the effect of $`J_{6p}`$; this will increase the stability of the antiferromagnetic region to roughly $`r_s=200`$.)

At higher density, the frustration between large cycle exchanges (4-6 body loops) results in a disordered spin state. For example, the point ($`J_2^{\mathrm{eff}}/J_4=2`$, $`J_5=0`$, $`J_{6h}=0`$), close to the parameters at $`r_s=100`$, is a spin liquid with a gap to all excitations: The spin liquid properties can be understood from a resonance valance bond model with ordering into spin-1 diamond plaquettes. The breakup of the plaquettes is responsible for a low temperature peak on the specific heat. A second high temperature peak develops at a temperature $`TJ_c`$, shown as a dotted line on fig. (1).

We find that at higher densities ($`r_s<100`$), the trajectory of the MSE models parallels the F-AF phase line, with the possibility of a re-entrant ferromagnetic phase for $`r_s<40`$. Note that QMC calculations of the normal fermi liquid at $`r_s=30`$ show that the ferromagnetic phase has a slightly lower energy that the unpolarized phase. Hence both the high density 2dWC and the low density electron fluid are characterized by a spin Hamiltonian which is nearly ferromagnetic.

We note a remarkable similarity between the exchange parameters of the 2dWC to those extrapolated from measurements of the second layer of 3He absorbed on grafoil as shown in fig. (5). The existence of these two related spin systems should allow a fuller understanding of this remarkable spin liquid. Such an equivalence could arise from an underlying virtual vacancy-interstitial (VI) mechanism giving rise to ring exchanges. In this model the prefactor of the exchange is controlled by the rate of VI formation (non-universal) but the ratios of the various ring exchanges arises from geometry of the triangular lattice and from the attraction of point defects (universal). We have recently determined that the VI formation energy vanishes at melting in the 2dWC. This is consistent with the fact that the various $`J_P`$ increase rapidly near melting.

In summary, we find that the magnetic interactions are characterized by a frustrated spin order. Application of a magnetic field transforms the exchange frequencies to $`J_pe^{2\pi ieB_{}a_p/h}`$ where $`a_p`$ is the area of the exchange (see table I) and $`B_{}`$ the magnetic field. Experiments with magnetic fields will allow exploration of this Aharonov-Bohm effect and thereby provide direct information on ring exchanges.

The semiconductor realizations of the 2dWC have significant disorder which can stabilize localized electronic states at higher densities than in the clean system. One can also stabilize the Wigner crystal at higher densities using bilayers . Exchange frequencies in those systems including effects of layer thickness and the exchange properties of point defects could be calculated with the present method.

This research was funded by NSF DMR 98-02373, the CNRS-University of Illinois exchange agreement, support by Fundação de Amparo à Pesquisa do Estado de São Paulo (FAPESP) and the Dept. of Physics at the University of Illinois. We used the computational facilities at the NCSA.


text:
# Snapshot Identification of Gamma Ray Burst Optical Afterglows1footnote 11footnote 1Dedicated to the memory of Jan van Paradijs.

## 1 Introduction

The search for counterparts to gamma ray bursts (GRBs) at longer wavelengths has been recognized from the beginning as a key to understanding the bursts’ nature. This motivated models that predicted gamma ray burst afterglows from hot material in relativistically expanding GRB remnants (Paczyński & Rhoads 1993; Katz 1994; Mészáros & Rees 1997). The afterglow was predicted to arise from synchrotron radiation, and its spectrum is therefore well described by a broken power law. Both the prediction of fading counterparts with power law spectra and the general expectation that GRB counterparts would tell us a great deal about the bursts have been amply borne out by observation since the first successful observations of GRB afterglows in 1997 (Costa et al 1997; van Paradijs et al 1997; Frail et al 1997).

The greatest rewards so far have come from detailed studies of individual afterglows with the best time and spectral coverage. Rapid afterglow identifications are a key to achieving such data sets, because some essential types of data cannot be obtained before the afterglow position is known to $`1`$ arcminute (near-infrared and submillimeter photometry) or even $`1`$ arcsecond (spectroscopy). Traditionally, afterglows are found through a variability search within GRB error boxes. This requires multiple epoch observations whose separation in time $`\mathrm{\Delta }t`$ must be a reasonable fraction ($`0.3`$) of the time elapsed from the GRB to the first epoch. This lag is usually further increased by telescope scheduling issues to $`1`$ day. Moreover, afterglow light curves occasionally include “plateau” phases (e.g. GRB 970508 at $`t<2`$ days \[Pedersen et al 1998\]; GRB 971214 at $`t<0.3`$ days \[Gorosabel et al 1998\]; and GRB 000301C at $`t6`$ days \[Bernabei et al 2000; Masetti et al 2000\]). Two-epoch observations of such afterglows could miss them entirely.

This delay in afterglow identification is not necessary, because afterglows are characterized not only by fading behavior but also by power law spectra. Such spectra are relatively rare in the optical sky, which is dominated by Galactic stars at bright flux levels and by galaxies at faint flux levels. Relative to a stellar spectrum, a GRB afterglow is blue at blue wavelengths ($`0.4\mu m`$), and red at red wavelengths ($`0.6\mu m`$). A color-color plot should therefore isolate a GRB afterglow from nearby stars efficiently, provided it spans a sufficient wavelength range to detect spectral curvature (or the Balmer jump) in stars of any surface temperature. Under many circumstances it will also be possible to identify afterglows using simpler though less robust single color criteria (either ultraviolet excess \[“UVX”\] or infrared excess \[“IRX”\]). Moreover, absolute photometric calibration is not required so long as a sufficient solid angle is observed to empirically determine the stellar locus in color space.

Section 2 describes model color-color plots that show the expected loci of GRB afterglows, stars, galaxies, and quasars. Section 3 demonstrates empirically that the color-color method could have identified the afterglow of GRB 000301C in a single epoch. Finally, section 4 considers advantages and limitations of the method, including confusion limits as a function of flux level, likely selection effects in the sample of afterglows found through color-based searches, and applicability of the method to “orphan afterglow” searches.

## 2 Color-Color Space Models

To explore this method in detail, I have calculated synthetic color space locations for GRB afterglows, stars, galaxies, and quasars. Model calculations were performed for the widely available Johnson U, B, and V filters and Cousins R and I filters, though color-color searches will work for any set of filters spanning a sufficient wavelength range. Filter transmission curves were taken from Bessell (1990). The photometric zero point for all Johnson-Cousins system synthetic colors was set to the Vega model atmosphere of Kurucz (1979).

Model spectra of objects at high redshift were attenuated using the mean intergalactic absorption given by Madau (1995)222The photoelectric absorption was treated using the approximation to Madau’s equation 16 from his footnote 3, which is within 5% of the full calculation.. For redshifts $`z2`$, the flux in the bluest filter can be greatly reduced by neutral hydrogen absorption in the intergalactic medium. This absorption consists of a superposition of Lyman $`\alpha `$ forest lines, higher order Lyman lines, and bound-free Lyman continuum absorption (Madau 1995). The net effect is a shift in the bluest color measured. The relevant redshift depends on the choice of bluest filter. The Lyman $`\alpha `$ line is redshifted to the peak transmission wavelength of the Johnson U band at $`z=2.04`$, and the Johnson B filter around $`z=2.5`$.

In addition, I calculated reddening vectors for dust both in the Milky Way and in afterglow host galaxies at a range of redshifts, using the empirical fits of Pei (1992) to the extinction laws of the Milky Way and the Magellanic Clouds.

In the following subsections, I describe the color calculations for each class of object considered, and then summarize the results.

#### Stars:

Stellar colors were derived from the spectrophotometric catalog of Gunn and Stryker (1983), which spans spectral types from O to M and luminosity classes from the main sequence to supergiants.

#### GRB Afterglows:

Afterglows were taken as pure power law spectra, $`f_\nu \nu ^\alpha `$, with indices $`0.5\alpha 1.5`$. Such spectra describe synchrotron emission from a power law distribution of electrons away from break frequencies. Deviations from power law spectra occur near spectral breaks (see Sari, Piran, & Narayan 1998) and in the case of interstellar absorption (section 4). The range of spectral indices considered corresponds to the commonly measured range of electron energy indices, $`2<p<3`$. Here $`p`$ describes the number-energy distribution of electrons accelerated at the GRB remnant’s external shock, $`N()^p`$, and I have considered spectra both above and below the cooling frequency.

#### Galaxies:

Galaxy spectra were taken from the GISSEL stellar population synthesis code (Bruzual & Charlot 1993) for a few representative star formation histories. No attempt was made to include galaxy evolution effects, which may result in elliptical galaxy colors that are slightly too red at high redshift (see, e.g., Bruzual ). The resulting galaxy colors agree with those computed by Fukugita, Shimasaku, & Ichikawa (1995) to within $`0.1\text{mag}`$ for most filters, galaxy spectral types, and redshifts; this is acceptable given the differences in the input spectra in the two works.

Galaxy spectra at optical wavelengths are a superposition of diverse stellar spectra, modified slightly by nebular emission and dust effects. This superposition implies that galaxies have less sharply peaked spectra than individual stars and so lie on the same side of the color-color space stellar locus as power law spectra. Indeed, the spectrum of a star forming galaxy can be a good power law over most of the optical window. Fortunately, the $`4000`$Å break due to stellar atmospheric features interrupts this approximate power law and can be used to distinguish almost all galaxies at $`z1`$ from afterglow candidates given sufficient color information. With observations in three filters, there is some redshift where the $`4000`$Å break lies in the middle filter, and galaxy colors can match those of a pure power law. The addition of a fourth filter breaks this degeneracy. For example, in the $`UB`$, $`RI`$ plane, an irregular galaxy at $`z0.4`$ lies on the power law sequence; however, its $`BR`$ color is much redder than the corresponding power law. In addition to color criteria, it is of course possible to distinguish galaxies from GRB afterglows morphologically. The fraction of galaxies that appear extended depends on the available spatial resolution, with the a few galaxies appearing pointlike even in Hubble Space Telescope images, but subarcsecond resolution will in practice greatly reduce confusion due to galaxies.

#### Quasars:

Quasar colors were calculated using the composite quasar spectrum by Francis et al (1991), extrapolated to $`\lambda >0.6\mu m`$ with $`f_\lambda \lambda ^{2.18}`$. This extrapolation will result in modest color errors for low redshifts where the H$`\alpha `$ line falls in a measured filter, but should not affect our conclusions greatly.

Quasar spectra are reasonably well approximated by power laws (plus broad emission lines) across a wide wavelength range. Color-based criteria have been used to identify quasars for nearly 40 years (see the review by Burbidge 1967), and the use of far-red filters to supplement bluer UBV passbands was suggested as early as 1967 (Braccesi 1967; Braccesi, Lynds, & Sandage 1968). Quasars are therefore extremely difficult to distinguish from GRB afterglows using only broad band colors. The main discriminant available is that quasars tend to be somewhat bluer than typical afterglows. In practice, we expect confusion with quasars to set the maximum area over which an unambiguous afterglow identification can be expected. This confusion can be resolved with low resolution spectra, which identify most quasars unambiguously by their strong, broad emission lines.

### 2.1 Results

Sample color-color diagrams are shown for the $`UV`$$`VI`$ plane in figure 1, and the $`BV`$$`RI`$ plane in figure 2. These have been chosen as representative cases where the two-color method works. Additional combinations of filters are also effective provided that at least three filters are used, the bluest filter is at least as blue as the Johnson B band, and the reddest is at least as red as Cousins I. Unfortunately, the observationally easiest filters (BVR) do not yield an adequate separation between afterglow and stellar colors to expect reliable identification of afterglows, barring unusually high photometric precision.

Given a sufficient filter set, the afterglows are well separated from the stellar locus. Moreover, the reddening vector for Milky Way dust runs essentially parallel to both the stellar sequence and the power law spectral sequence. Thus, even though reddened afterglow spectra are no longer strict power laws, they remain distinct in color-color space.

Of greater concern is dust in the GRB host galaxy. Because this dust is at high redshift, its $`2175`$Å feature can enter the observed wavelength range, changing the direction of the reddening vector in color-color space. For redshifts $`2.2z3.3`$, this feature falls in the reddest filter of the color-color plot, thus making the long-wavelength color (e.g., V-I) slightly bluer. A moderate amount of dust ($`A_V=0.5`$) can then move the afterglow locus onto the stellar locus in color-color space. However, this effect is only important for Milky-Way type dust. The most important consequence of dust for color-based afterglow searches is therefore not reddening but extinction, which can reduce the received flux below the detection limit of the search whatever the reddening law. However, GRBs can destroy dust at distances up to $`10pc`$ (Waxman & Draine 2000) or beyond (Fruchter, Krolik, & Rhoads 2001), reducing these concerns substantially for bursts at high Galactic latitude.

## 3 Empirical Tests

Archival multicolor data is available for a few GRB fields from the U.S. Naval Observatory’s GRB followup team (see Henden et al 2000). I have used this data for the field of GRB 000301C, together with published photometry of the afterglow (Jensen et al 2000), for an empirical demonstration of color-based optical afterglow searches.

Two color-color plots for this field are shown in figure 3. In both cases, we see that the afterglow is an outlier in color-color space. In the $`UV`$$`VI`$ figure, there are a few other outliers, all of which are systematically redder than the afterglow in both colors. In the $`BV`$$`RI`$ plot, the afterglow is the single most dramatic outlier and would be the first target chosen for prompt spectroscopy or other rapid narrow-field followup observations.

## 4 Discussion

Color-based searches have the potential to identify GRB afterglows faster than a two-epoch variability search, and hence to allow more thorough followup of many afterglows. There are three basic color selection criteria available for afterglows: IR excess (“IRX”), UV excess (“UVX”), and Lyman break. The appropriate combination of these depends on both the redshift of the burst and the available instrumentation. The primary focus of this paper is on color-color plane methods, which combine the IR excess and UV excess criteria and are therefore more robust than either criterion alone.

Still, as the observations of GRB 000301C demonstrate, potentially occupied regions of color-color space remain vacant in practice for high-latitude fields observed to a limiting magnitude $`V21`$. In particular, there are no very hot stars, few or no high redshift galaxies, and few or no quasars. Under these circumstances, a simple UV Excess (“UVX”) test could identify the afterglow on the basis of a single color measurement (preferably including one filter at $`\lambda <4000`$Å).

Likewise, at near-infrared wavelengths, most stars and galaxies will have spectral slopes near the Rayleigh-Jeans value ($`2`$) and will therefore be bluer than afterglows. This is again demonstrated by GRB 000301C, for which we found $`JK^{}1.45`$ (Rhoads & Fruchter 2001), about $`0.3`$ magnitudes redder than even the coolest stars tabulated by Johnson (1966). However, UVX and IRX methods are vulnerable to the presence of very hot stars and highly reddened sources respectively, while the two color method is not.

The Lyman break (or “dropout”) method uses intergalactic absorption to find high redshift afterglows whose strongly attenuated blue flux places them in a region of color space occupied by essentially all classes of high redshift object, on the opposite side of the stellar locus from low redshift afterglows. The observational requirements to find these bursts from Lyman break colors are more exacting. Long integrations in the bluest filter are required to conclusively identify “dropout” afterglows, while lower redshift afterglows are unusually bright in the bluest filter and require correspondingly shorter integrations. Even so, most afterglows are brighter than their host galaxies at discovery, which means that a high redshift afterglow may be an easier dropout to identify than a high redshift galaxy.

Color-based searches are subject to certain selection effects, which must be considered in interpreting their results. These effects include: (1) Redshift: The onset of intergalactic hydrogen absorption in the bluest filter observed determines the highest redshift where UV excess methods can work and the lowest redshift where Lyman break methods can work. The UVX method will fail at $`z2`$ for U band and $`z2.5`$ for B band, while the Lyman break method first becomes practical in the U band at $`z3.3`$. Thus, two-color optical searches will not identify afterglows between about $`z=2.5`$ and $`z=3.3`$. The IR excess will continue to be practical to quite high redshift, e.g. to $`z15`$ for a $`JK`$ color, since the onset of neutral hydrogen absorption simply accentuates an IR excess until the absorption reaches the redder filter. Because there is considerable variance in the hydrogen column density along different lines of sight (Madau 1995), the redshift boundaries discussed here are not sharp demarcations. Rather, the detection efficiency declines smoothly, and the characteristic redshifts given above correspond to about $`50\%`$ of peak efficiency.

(2) Heavily obscured afterglows can drop below the search’s detection threshold. Extinction affects any search method, though it is exacerbated for UVX and two-color methods by the increase of extinction towards bluer optical-UV wavelengths and the necessity of observations in observer frame blue light. Additionally, there will be a mild selection against finding afterglows in galaxies with Milky Way type dust extinction at redshifts $`2.2z3.3`$ using two-color optical searches. However, this overlaps the stronger selection due to intergalactic Hydrogen absorption to such an extent as to be almost unimportant.

(3) For gamma ray bursts occurring near dense molecular gas, absorption by excited molecular hydrogen can occur for rest wavelengths as long as $`1650`$Å (Draine 2000). Absorption between $`1650`$ and $`1300`$Å can reach 50% for H2 column densities around $`10^{19}`$ cm-2. This will act to exclude afterglows in molecular clouds from UVX or two-color samples at redshifts beyond $`1.3`$ (using U band) or $`1.65`$ (using B).

(4) Because afterglow spectra are broken power laws rather than pure power laws, there will be times in the evolution of each afterglow when its colors deviate from the models in section 2. Breaks in synchrotron spectra are not absolutely sharp, so this is unlikely to be a major selection effect for purely optical searches, where the wavelength range spanned is only a factor of $`2.5`$. Searches employing larger wavelength coverage (UV or IR) will be more sensitive to this effect.

(5) Afterglows located atop brighter host galaxies may be hard to find because the dominant color is that of the galaxy. This effect will depend on instrumental spatial resolution as well as the intrinsic properties of the afterglow and host galaxy. (Variability searches can be similarly compromised unless image subtraction methods are used.)

(6) An advantage of this method is that short-term “plateaus” in afterglow light curves do not affect it. It may therefore find some afterglows that traditional variability searches miss.

Even afterglows that lie in the predicted region of color space may be difficult to identify by these methods in cases where the error box is so large that a substantial number of other sources occupy the same region of color space. The critical error box size where one “confusing” source is expected depends on the magnitude limit of the search. Ultimately, the quasar number-flux relation gives a lower limit to the confusion level. Additionally, for three filter searches, star forming galaxies at $`z0.5`$ can add to the confusion. To achieve a unique afterglow identification, it is necessary to match the magnitude limit of the search to the size of the error box, and to observe the field while the afterglow remains brighter than this magnitude limit. A unique identification may not always be needed, however. In particular, a list of several viable candidates may be enough to obtain a same-night afterglow spectrum using a multiobject spectrograph with a sufficient field of view and rapid setup procedure. The correct candidate could then be identified later based on either spectra or variability.

The quasar number-magnitude relation for $`z<2.2`$ and magnitude $`B<24`$ is predicted by Malhotra & Turner (1995) for both ($`\mathrm{\Omega }=1`$, $`\mathrm{\Lambda }=0`$) and ($`\mathrm{\Omega }=0.1`$, $`\mathrm{\Lambda }=0.9`$) cosmologies, based on observations to $`B<22`$ by Boyle et al (1987) and Boyle, Shanks, & Peterson (1988).

The approximate range of redshift where star forming galaxies overlap the GRB afterglow locus in color-color space is $`0.35z0.6`$ for the $`UV`$, $`VI`$ diagram, or $`0.45z0.6`$ for the $`BV`$, $`RI`$ diagram. I estimate a number-magnitude relation for such galaxies using both the local luminosity function from Loveday et al (1992) and direct counts of galaxies with $`0.55<VI<1.28`$ and $`0.4<z<0.6`$ from the Canada-France Redshift Survey (Lilly et al 1995a; Le Fevre et al 1995; Hammer et al 1995). The predicted number counts amount to approximately 10% of all galaxy number counts (see, e.g., Lilly et al 1995b, Williams et al 1996) around $`R=22`$. The quasar and galaxy count estimates (converted to R band magnitudes using the colors of an $`f_\nu \nu ^1`$ power law) are plotted together in figure 4. We see that quasars are the more important foreground for $`R<20`$. At fainter magnitudes, galaxies appear more important, but this may be misleading since they can often be identified using a third independent color and/or morphological information.

Presently, most GRB error boxes have characteristic sizes of a few tens of square arcminutes. Thus, the brightest quasar expected in the error box has a magnitude $`R22`$, while the brightest galaxy having colors similar to an afterglow might be a magnitude brighter. Both this depth and the error box area are reasonably well matched to the capabilities of moderate sized (1 to 2.5 meter) ground-based telescopes.

### 4.1 Optimizing the Observations

To realize the promise of same-night afterglow identifications, a rapid and standard data analysis procedure is desirable, and could be aided by an observing procedure optimized for quick analysis. Archival bias frames and flatfield frames for all filters could be useful, and precisely cospatial exposures in multiple filters (i.e., without offsets and with telescope guiding) could ease measurement of object colors.

An additional constraint on observing strategy arises because the target is likely to fade during the period of observations. If the multiple wavelengths cannot be observed simultaneously, it is then better to begin with the bluest filter, followed by the reddest, and only then move on to the middle filters. This way, fading behavior will accentuate the expected color offset rather than reducing it. (The choice of bluest filter before reddest is driven by the possibility that a spectral break may move through the observed wavelength range during the course of observations.)

The optimal choice of filter set will depend on local conditions. The relative difficulty of U band observations is approximately balanced by the large color-plane offsets obtained with U data ($`0.7\text{mag}`$ from the stellar locus), which relax the signal to noise requirements considerably.

### 4.2 Orphan Afterglow Searches

If gamma ray bursts are highly collimated, we expect to observe many more afterglows than GRBs (Rhoads 1997). The search for such “orphan afterglows” is an important and largely model-independent way of testing GRB collimation.

Color-based techniques will be a valuable addition to orphan afterglow searches. Afterglows remain rare phenomena even in extreme collimation scenarios, and in order to find one it will be necessary to survey a large area of sky that will contain many quasars. Thus, orphan afterglow searches must combine variability with color tests to be reliable. When enough filters are observed simultaneously in a field where an archival first epoch is available, it will be possible to identify new sources, test their colors for consistency with a power law spectrum, and start intensive followup of any viable afterglow candidates. This will allow much greater confidence in the identification of orphan afterglows than is presently possible, since good spectra and light curves can then be used to distinguish real afterglows from variable stars, high redshift supernovae, active galactic nuclei, microlensing events, and other optically variable sources.

As an illustrative example, consider Sloan Digital Sky Survey (SDSS; see Knapp et al 1999). The SDSS southern strip will cover an area $`200`$ degree2 to about 22nd magnitude in 5 broadband filters many ($`45`$) times. The known GRB rate is $`3`$ per sky per day, and bright afterglows remain above 22nd magnitude for a period $`4`$ days. Thus the southern strip has an effective coverage of $`4\times 10^4\text{degree}^2\text{day}`$. The northern SDSS survey will cover $`10^4`$ degree2 in a single imaging epoch. Fortunately, GRB afterglows would likely be targeted for spectra based on their quasar-like colors, yielding a second epoch. The northern survey thus provides another $`4\times 10^4\text{degree}^2\text{day}`$. Both parts together would be expected to detect $`6`$ GRB afterglows. If GRBs are tightly collimated (opening angles $`10^{}`$), the expected number of orphan afterglows rises to $`25`$. By using the color space properties to identify these afterglows, such a survey can place strong limits on the collimation angle and hence energy requirements of gamma ray bursts.

## 5 Conclusions

Color-based searches will allow faster identification of many gamma ray burst afterglows. Intensive followup of these afterglows can then begin earlier and at brighter flux levels, yielding a larger sample of bursts for which spectroscopic redshifts and well measured spectral energy distributions (including roughly synoptic data at X-ray, optical, near-infrared, submillimeter, and radio wavelengths) can be obtained. Such data sets are crucial to detailed studies of afterglow physics. Additionally, single-epoch color-based searches may yield a more easily characterized statistical sample of afterglow detections and nondetections than is currently available. Such samples will allow reliable studies of the population properties of gamma ray burst afterglows, and thereby open the way for new insights into the nature of gamma ray bursts.

I thank Sangeeta Malhotra for extensive discussions and comments on an early draft of this paper; Andy Fruchter for suggesting a harder look at IR excess methods; and Mauro Giavalisco, Arjun Dey, Norm Grogin, and Richard Green for additional discussions. Finally, I thank Arne Henden and the USNO GRB team for providing freely the photometric calibration data for GRB fields that allowed the first empirical test of this method. This work was supported by an Institute Fellowship at The Space Telescope Science Institute (STScI), which is operated by AURA under NASA contract NAS 5-26555.