text:
# What causes a neuron to spike?

## 1 Introduction

There are two distinct approaches to characterizing a neuron based on spike trains recorded during stimulus presentation. One can take the view of a downstream neuron, and ask how the spike train should be decoded to reconstruct some aspect of the stimulus. Alternatively, one can take the neuron’s own view, and try to determine how the stimulus is encoded, i.e. which aspects of the stimulus trigger spikes. This paper is concerned with the latter problem: the identification of relevant stimulus features from neural data.

The earliest attempts at neural characterization, including classic electrophysiological experiments like those of Adrian , presented the neural system with simple, highly stereotyped stimuli with one or at most a few free parameters; this makes it relatively straightforward to map the input/output relation, or tuning curve, in terms of these parameters. While this approach has been invaluable, and is still often used to advantage, it has the shortcoming of requiring a strong assumption about what the neuron is ‘looking for’. As this assumption is relaxed by probing the system using stimuli with a larger number of parameters, the length of the experiment needed to probe neural response grows combinatorially.

In some cases, white noise can be used as an alternative stimulus which naturally contains ‘every’ signal . The stimulus feature best correlated with spiking is then the spike triggered average (STA), i.e. the average stimulus history preceding a spike . Such reverse correlation methods have been further developed and applied to characterize the motion-sensitive identified neuron H1 in the fly visual system . By extending reverse correlation to second order, one can identify multiple directions in stimulus space relevant to the neuron’s spiking ‘decision’ . These techniques allow the experimenter to probe neural response stochastically, using minimal assumptions about the nature of the relevant stimulus.

In this paper, we will examine in detail the relation between the stimulus features extracted through white noise analysis and the computation actually performed by the neuron. In the usual interpretation of white noise analysis, the features extracted using reverse correlation are considered to be linear filters which, when convolved with the stimulus, provide the relevant inputs to a nonlinear spiking decision function. While the features recovered by reverse correlation are optimal for stimulus reconstruction , we will show that they are distinct from the filter or filters relevant to the neuron’s spiking decision. In particular, reconstruction filters necessarily depend on the stimulus ensemble, while the filters relevant for spiking—at least, for simple neurons like leaky integrate-and-fire—are properties of the model alone.

The main difficulty in interpreting the spike triggered average is that spikes interact, in the sense that the occurrence of a spike affects the probability of the future occurrence of a spike. Both stimulus reconstruction and the identification of triggering features are affected by this issue. Some of the effects of correlations between spikes have been treated in . To separate the influence of previous spikes from the influence of the stimulus itself, we propose to analyze only ‘isolated’ spikes—spikes sufficiently distant in time from previous spikes that they can be considered to be statistically independent. This procedure allows us to recover exactly the relevant stimulus features. However, it introduces additional complications. By requiring an extended silence before a spike, we bias the prior statistical ensemble. This bias will emerge in the white noise analysis, and we will have to identify it explicitly.

We carry out this program for the simplest model neuron with temporal dynamics: the leaky integrate-and-fire neuron. We choose this model because, unlike more complex models, we know exactly which stimulus feature causes spiking, and can make a precise comparison with the output of the analysis.

The methods presented and validated in the present work using the integrate-and-fire model have been applied to the more biologically interesting Hodgkin–Huxley neuron , resulting in new insights into its computational properties. Application of the same analytic techniques to real neural data is currently underway.

## 2 Characterizing neural response

In this section we briefly review material that has appeared elsewhere .

### 2.1 The linear/nonlinear model

Consider a neuron responding to a dynamic stimulus $`I(t)`$. While here we will consider $`I`$ to be a scalar function, in general, it may depend on any number of additional variables, such as frequency, spatial position, or orientation. We can suppress these dependencies without loss of generality. Even without additional dependent variables, $`I(t)`$ is a high-dimensional stimulus, as the occurrence of a spike at time $`t_0`$ generally depends on the stimulus history, $`I(t<t_0)`$. The effective dimensionality of this input is determined by the temporal extent of the relevant stimulus, and the effective sampling rate, which is imposed by biophysical limits. For the computation of the neuron to be characterizable at all, the dimensionality of the stimulus relevant to spiking must generally be much lower than this total dimensionality. The simplest approach to dimensionality reduction is to approximate the relevant dimensions as a low-dimensional linear subspace of $`I(t)`$.

If the neuron is sensitive only to a low dimensional linear subspace, we can define a small set of signals $`s_1,s_2,\mathrm{},s_K`$ by filtering the stimulus,

$$s_\mu =_0^{\mathrm{}}𝑑\tau f_\mu (\tau )I(t_0\tau ),$$
(1)
so that the probability of spiking depends only on these few signals,

$$P[\mathrm{spike}\mathrm{at}t_0|I(t<t_0)]=P[\mathrm{spike}\mathrm{at}t_0]g(s_1,s_2,\mathrm{},s_K).$$
(2)
Classic characterizations of neurons in the retina, for example , typically assume that these neurons are sensitive to a single linear dimension of the stimulus; the filter $`f_1`$ then corresponds to the (spatiotemporal) receptive field, and $`g`$ is proportional to the one-dimensional tuning curve. A number of other neurons, particularly early in sensory pathways, have been similarly characterized, e.g. . Here we will concentrate primarily on the problem of finding the filters $`\{f_\mu \}`$; once they are found, $`g`$ can be obtained easily by direct sampling of $`P[\mathrm{spike}\mathrm{at}\mathrm{t}_0|s_1,s_2,\mathrm{},s_K]`$, assuming that the relevant dimensionality $`K`$ is small.

### 2.2 Reverse correlation methods

While one would like to determine the neural response in terms of the probability of spiking given the stimulus, $`P[\mathrm{spike}\mathrm{at}t_0|I(t<t_0)]`$, we follow  in considering instead the distribution of signals conditional on the response, $`P[I(t<t_0)|\mathrm{spike}\mathrm{at}t_0]`$. This quantity can be extracted directly from the known stimulus and the resulting spike train. The two probabilities are related by Bayes’ rule,

$$\frac{P[\mathrm{spike}\mathrm{at}t_0|I(t<t_0)]}{P[\mathrm{spike}\mathrm{at}t_0]}=\frac{P[I(t<t_0)|\mathrm{spike}\mathrm{at}t_0]}{P[I(t<t_0)]}.$$
(3)
$`P[\mathrm{spike}\mathrm{at}t_0]`$ is the mean spike rate, and $`P[I(t<t_0)]`$ is the prior distribution of stimuli. If $`I(t)`$ is Gaussian white noise, then this distribution is a multidimensional Gaussian; furthermore, the distribution of any filtered version of $`I(t)`$ will also be Gaussian.

The first moment of $`P[I(t<t_0)|\mathrm{spike}\mathrm{at}t_0]`$ is the spike-triggered average,

$$\mathrm{STA}(\tau )=[dI]P[I(t<t_0)|\mathrm{spike}\mathrm{at}t_0]I(t_0\tau ).$$
(4)
We can also compute the covariance matrix of fluctuations around this average,

$$C_{\mathrm{spike}}(\tau ,\tau ^{})=[dI]P[I(t<t_0)|\mathrm{spike}\mathrm{at}t_0]I(t_0\tau )I(t_0\tau ^{})\mathrm{STA}(\tau )\mathrm{STA}(\tau ^{}).$$
(5)
In the same way that we compare the spike-triggered average to some constant average level of the signal in the whole experiment, an important advance of  was to compare the covariance matrix $`C_{\mathrm{spike}}`$ with the covariance of the signal averaged over the whole experiment,

$$C_{\mathrm{prior}}(\tau ,\tau ^{})=[ds]P[s(t<t_0)]s(t_0\tau )s(t_0\tau ^{}).$$
(6)
If the probability of spiking depends on only a few relevant features, the change in the covariance matrix, $`\mathrm{\Delta }C=C_{\mathrm{spike}}C_{\mathrm{prior}}`$, has a correspondingly small number of outstanding eigenvalues. Precisely, if $`P[\mathrm{spike}|\mathrm{stimulus}]`$ depends on $`K`$ linear projections of the stimulus as in Eq. (2), and if the inputs $`I(t)`$ are chosen from a Gaussian distribution, then the rank of the matrix $`\mathrm{\Delta }C`$ is exactly $`K`$. Further, the eigenmodes associated with nonzero eigenvalues span the relevant subspace. A positive eigenvalue indicates a direction in stimulus space along which the variance is increased relative to the prior, while a negative eigenvalue corresponds to reduced variance.

## 3 Leaky integrate-and-fire

The leaky integrate-and-fire neuron is perhaps the most commonly used approximation to real neurons. Its dynamics can be written:

$$C\frac{dV}{dt}=I(t)\frac{V}{R}CV_c\underset{i}{}\delta (tt_i),$$
(7)
where $`C`$ is a capacitance, $`R`$ is a resistance and $`V_c`$ is the voltage threshold for spiking. The first two terms on the right hand side model an RC circuit; the capacitor integrates the input current as a potential $`V`$ while the resistor dissipates the stored charge (hence ‘leaky’). The third term implements spiking: spikes are defined as the set of times $`\{t_i\}`$ such that $`V(t_i)=V_c`$. When the potential reaches $`V_c`$, a restoring current is instantaneously injected to bring the stored charge back to zero, resetting the system. This formulation emphasizes the similarity between leaky integrate-and-fire and the more biophysically realistic Hodgkin–Huxley equation for current flow through a patch of neural membrane; in the Hodgkin–Huxley system, however, the nonlinear spike-generating term is replaced with ionic currents of a more complex (and less singular) form. Hence the simple properties of the leaky integrate-and-fire model—a single degree of freedom $`V`$, instantaneous spiking, total resetting of the system following a spike, and linearity of the equation of motion away from spikes—are only rough approximations to the properties of realistic neurons.

Integrating Eq. (7) away from spikes, we obtain

$$CV(t)=_{\mathrm{}}^t𝑑\tau \mathrm{exp}\left(\frac{\tau t}{RC}\right)I(\tau ),$$
(8)
assuming initialization of the system in the distant past at zero potential. The right hand side can be rewritten as the convolution of $`I(t)`$ with the causal exponential kernel

$$f_{\mathrm{}}(\tau )=\{\begin{array}{cc}0\hfill & \text{if }\tau 0\hfill \\ \mathrm{exp}(\tau /RC)\hfill & \text{if }0<\tau \text{,}\hfill \end{array}$$
(9)
so the condition for a spike at time $`t_1`$ is that this convolution reach threshold:

$$CV_c=_{\mathrm{}}^{\mathrm{}}𝑑\tau f_{\mathrm{}}(\tau )I(t_1\tau ).$$
(10)
The interaction between spikes complicates this picture somewhat. If the system spiked last at time $`t_0<t_1`$, we must replace the lower limit of integration in Eq. (8) with $`t_0`$, as no current injected before $`t_0`$ can contribute to the accumulated potential. However, if we wish always to evaluate stimulus projections by integrating over the entire current history $`I(t)`$ as in Eq. (10), then we must replace $`f_{\mathrm{}}`$ with an entire family of kernels dependent on the time to the previous spike, $`\mathrm{\Delta }t=t_1t_0`$:

$$f_{\mathrm{\Delta }t}(\tau )=\{\begin{array}{ccc}0\hfill & \text{if }\tau 0\hfill & \\ \mathrm{exp}(\tau /RC)\hfill & \text{if }0<\tau <\mathrm{\Delta }t\hfill & \\ 0\hfill & \text{if }\mathrm{\Delta }t\tau \text{.}\hfill & \end{array}$$
(11)
This equation completes an exact model of the leaky integrate-and-fire neuron.

## 4 Response to white noise

We now consider the statistical behavior of the integrate-and-fire neuron when driven with Gaussian white noise current of standard deviation $`\sigma `$,

$`I(t)`$ $`=`$ $`0,`$
$`I(t)I(t^{})`$ $`=`$ $`\sigma ^2\delta (tt^{}).`$ (12)
For the following analysis, we take $`R=10`$ k$`\mathrm{\Omega }`$, $`C=1`$ $`\mu `$F, $`V_c=10`$ mV and $`\sigma =\sqrt{200}`$ $`\mu `$A, and use Euler integration with a time step of 0.05 msec. These values were chosen to mimic biophysically plausible parameters for a 1 $`\mathrm{cm}^2`$ patch of membrane, but none of our results depend qualitatively on this choice.

### 4.1 Unqualified reverse correlation

At first glance, it would appear that the neuron is perfectly described by a reduced model in one dimension: the current filtered by the exponential kernel of Eq. (9). If so, the spike-triggered average under white noise stimulation should be proportional to this filter.111In possible linear combination with a derivative-like filter, arising from the additional constraint that the threshold must be crossed from below. This point will be addressed in more detail later. However, as is evident in Fig. 1, the STA is not exponential, nor does it asymptote to the appropriate decay rate $`1/RC`$. One reason is that, in the presence of previous spikes, the integrate-and-fire neuron is not low-dimensional. Recall that the relative timing $`\mathrm{\Delta }t`$ of the previous spike selects a filter $`f_{\mathrm{\Delta }t}`$ (Eq. (11)); these filters are truncated versions of the original exponential $`f_{\mathrm{}}`$. The STA will therefore include contributions from each of these filters, weighted by their probability of occurrence (as determined by the interspike interval distribution).222This is not a prescription for recovering $`f_{\mathrm{}}`$ from the STA: the STA is further influenced by other effects, including the average stimulus histories leading up to previous spikes. It is readily shown that the orthogonal component of each filter $`f_{\mathrm{\Delta }t}`$ is a $`\delta `$-function at $`\mathrm{\Delta }t`$; so the set of relevant filters spans the entire stimulus space.

Covariance analysis reveals this high dimensionality very clearly. In Fig. 2 we show the spectrum of eigenvalue magnitudes of $`\mathrm{\Delta }C`$ as a function of the number of spikes summed to accumulate the matrix; with decreasing noise, an arbitrary number of significant eigenmodes emerge. As shown in Fig. 3, none of these eigenmodes correspond to an exponential filter. The first mode most closely resembles an exponential, but as shown in the right half of Fig. 3, it is more peaked near $`t=0`$, and its time constant is very different from $`RC`$.

### 4.2 Isolated spikes

It is clear that making progress requires first considering the causes of a single spike independently of the spiking history. Therefore we will include in our analysis only spikes which follow previous spikes after a sufficiently long interval. The appropriate interval can be estimated from the interspike interval distribution, Fig. 4.

As in many real neurons, the interval distribution has three main features: a ‘refractory period’ during which spikes are strongly suppressed, a peak at a preferred interval, and an exponential tail. Intervals in the exponential tail follow Poisson statistics; for intervals in this range, we can take spikes to be statistically independent. With reference to Fig. 4, we set our isolated spike criterion—very conservatively—as 75 msec of previous ‘silence’.

As shown in Fig. 5, even with this constraint, we do not recover a spike-triggered average matching $`f_{\mathrm{}}`$. We should note that this is no longer, strictly speaking, a spike-triggered average, but an event-triggered average, where the event includes the preceding silence. Thus we can expect the silence itself to contribute to any average we construct. To first order, during silence there is a constant negative bias, since certain positive currents which would have caused a spike are excluded from the average. Similarly, we expect that the covariance of the current will be altered from that of the Gaussian prior during extended silences. Because extended silence is not locked to any particular moment in time, the covariance matrix is stationary during silence.

We construct the covariance matrix as before: $`\mathrm{\Delta }C_{\mathrm{iso}}=C_{\mathrm{isolated}\mathrm{spike}}C_{\mathrm{prior}}`$. An entire series of modes once again appears, Fig. 6. This time, they arise because silence is, by definition, translationally invariant, resulting in a Fourier-like spectrum. Modes associated with the spike, on the other hand, are locked to a particular time, and should therefore not have Fourier-like modes. In particular, these spiking modes should have temporal support only in the time immediately before the spike. This implies that two types of modes—silence- and spike-associated—should emerge independently in the covariance analysis, and should exhibit a clear difference in their domain of support.

In Fig. 7, we consider the fraction of the energy of each mode over the interval $`t[65,45]`$ msec. At $`45`$ msec, the exponential filter $`f_{\mathrm{}}`$ has decayed to 1% of its peak value, so we are well away from the spike at $`t=0`$. The 75 msec silence criterion also puts the lower boundary of this interval well inside the stationary silence (at least 10 msec after the previous spike). As we see, exactly two modes emerge as local, with their energy fraction decaying with the noise in the matrix until reaching a stable value well below the rest of the modes.

The first localized mode, shown in Fig. 8, is the long-awaited exponential filter $`f_{\mathrm{}}`$. Using our modified second-order reverse correlation analysis on ‘experimental data’ we have thus finally been able to recover the stimulus feature to which the leaky integrate-and-fire neuron is sensitive.

### 4.3 Constraints

The second mode is a consequence of an implicit constraint in the spiking model. Naïvely, we know that at the moment when $`V=V_c`$, $`\dot{V}>0`$; that is, the threshold must be crossed from below. Hence

$$_{\mathrm{}}^{\mathrm{}}𝑑\tau I(t\tau )\frac{d}{d\tau }f_{\mathrm{}}(\tau )<0,$$
(13)
making the neuron also sensitive to the time derivative of $`f_{\mathrm{}}`$. The integrate-and-fire neuron is peculiar in that $`f_{\mathrm{}}`$ and $`df_{\mathrm{}}/d\tau `$ are almost linearly dependent:

$$\frac{d}{d\tau }f_{\mathrm{}}(\tau )=\{\begin{array}{cc}\delta (\tau )\hfill & \text{if }\tau 0\text{.}\hfill \\ (RC)^1\mathrm{exp}(\tau /RC)\hfill & \text{if }0<\tau \text{.}\hfill \end{array}$$
(14)
The orthogonal part of this filter is therefore only $`\delta (\tau )`$, and we can rewrite this extra condition more simply as $`I(0)>V_c/R`$.

However, we see that although the second spiking mode is indeed very sharply peaked at $`\tau =0`$, it is not a $`\delta `$-function. While the above argument is correct in constraining the possible values of $`I(0)`$, we have also required no threshold crossings for an extended time prior to the spike. Hence certain trajectories $`I(t)`$ which satisfy the two basic conditions ($`V(0)=V_c`$ and $`\dot{V}(0)>0`$) are still disallowed, because they would have led to a spike earlier. The second mode is therefore an extended version of the derivative condition. Although it appears to have very extended support, this is due to orthogonalization with respect to the first mode, $`f_{\mathrm{}}`$. The true constraint is positive definite and highly peaked, with virtually all of its energy concentrated in the $`5`$ msec interval prior to the spike.

We can more easily understand this constraint by considering the time-dependent distribution of $`V(t)`$ leading up to an isolated spike, shown in the left panel of Fig. 9. The process $`V(t)`$ leading up to a spike is a centrally biased random walk with an absorbing barrier at $`V_c`$ and a capacitatively determined correlation time. If we trace the distribution of $`V`$ backward in time from the spike at $`t=0`$, we see a constrained diffusive process from a $`\delta `$-function at the spike time to the steady-state silence distribution. The constraint on $`V`$ enforces a corresponding constraint on $`I`$, whose time-dependent distribution is shown in the right panel of Fig. 9. Notice that the $`I`$ distribution is only noticeably distorted from its silence steady-state at times very near the spike.

### 4.4 Silence modes

Finally, let us return to the spectrum of silence-associated modes in Fig. 7. Several of these are shown overlaid in Fig. 10. It is immediately clear that they have the expected Fourier-like structure over the period of silence. Their behavior near the spike is less obvious, however: each mode appears to execute an FM chirp, culminating in a sharp feature at the spike time itself. Clearly, while it is true by construction that spike-associated covariance modes have no support in periods of extended silence, the converse is not true. This structure is due to the essential non-stationarity of silence near spikes: at some point, the silence must come to an abrupt end.

It is clear that in order to implement our procedure for separating silence from spikes, we could not have taken the approach of looking for modes with significant energy near the spike. All of the silence modes have interesting structure near the spike, and could easily be mistaken for additional spike-generating features. Like the spike-associated modes, these silence modes have negative eigenvalues, indicating reduced stimulus variance in certain dimensions. Yet they do not cause spiking, nor do they inhibit it; they are in no meaningful sense ‘suppressive’ modes. They simply show us the second-order statistics of silences, which emerge from the enforced absence of spikes. This underscores the necessity of considering a sufficiently extended silence before isolated spikes for the translationally invariant part of the silence modes to emerge, distinguishing them from spiking modes.

This structure also highlights the reason why we could not, as one might expect, simply subtract from $`C_{\mathrm{isolated}\mathrm{spike}}`$ a prior obtained from extended silences instead of $`C_{\mathrm{prior}}`$ to obtain a simple, low-dimensional $`\mathrm{\Delta }C_{\mathrm{iso}}`$. While such a subtraction results in a local matrix around the spike, it cannot take into account the non-stationary part of the silence, and leaves the spike-related and non-stationary silence-related aspects of the covariance near the spike mixed. This again results in a high-dimensional matrix with incorrect (mixed) modes.

## 5 Adaptation

It has previously been observed that the STAs of some neurons depend on the stimuli used to probe their response . White noise stimulus was originally introduced to circumvent this difficulty by sampling stimulus space in an unbiased way; however, one must still choose a noise variance. While one would hope that this does not affect the STA, it has been shown that, at least in some sense, even the integrate-and-fire neuron ‘adapts’ to the stimulus variance , showing both a change in the overall STA and in its measured nonlinear decision function. The reason for this form of adaptation is simple: the statistics of spike intervals depend on the stimulus variance, as the time to reach threshold depends on the variance. By construction, however, neither the spike-triggering feature $`f_{\mathrm{}}`$ nor the decision function $`g`$ (here simply the threshold) of the integrate-and-fire neuron actually depends on the stimulus statistics. As we have discussed, the unqualified STA is a linear combination of the filters $`f_{\mathrm{\Delta }t}`$ conditioned on the time to the last spike; the contribution of each filter to the overall STA thus depends on the probability $`P(\mathrm{\Delta }t)`$ of the interspike interval $`\mathrm{\Delta }t`$. As different stimulus ensembles produce different interspike interval distributions, these coefficients are functions of the stimulus ensemble. Therefore even in the absence of adaptation in the linear filters $`f`$, the overall STA will be stimulus-dependent. Because sampling the input/output relation normally involves convolving the stimulus with the STA, the measured nonlinearity will also appear to be stimulus-dependent, despite the fact that the neuron’s nonlinear decision function $`g`$ (based on convolution of the stimulus with the fixed filters $`f`$) is also fixed. Hence this form of ‘adaptation’ does not reflect any changing property or internal state of the neuron in response to the stimulus statistics. Rather, it arises from the nonlinearity inherent in spiking itself.

This emphasizes the significance of our methods: we have introduced a way to extract from data the intrinsic, unique spike-triggering feature, removing the effects of the stimulus-dependent interspike interval statistics. These statistics are produced by the model when driven by a particular stimulation regime; they are not part of the model itself. If a neuron does in fact exhibit ‘true’ adaptation in the sense that the functions $`f`$ and $`g`$ are variance-dependent, then the isolated spike analysis presented here will reveal this. Such variance dependence can come about as a result of explicit gain control mechanisms, e.g. , or because nonlinearities in the subthreshold neural response select non-constant stimulus features, an effect which can play an important role in more complex models such as the Hodgkin–Huxley neuron .

## 6 Discussion

While the interspike interval statistics are not an intrinsic model property, they do depend on an aspect of the complete characterization of the neuron that we have explicitly neglected here: the interspike interaction. There are two senses in which addressing only isolated spikes appears to be a limitation of our method. One is practical: in most situations, this will substantially reduce the number of spikes considered in the analysis. Further, much of the work on white noise analysis to date has focussed on the problem of stimulus reconstruction, and our method does not provide a recipe for decoding non-isolated spikes. With respect to the second point, there are good reasons for considering isolated spikes first. Evidence from several systems suggests that isolated spikes contain more information about the stimulus than the spikes that follow them . This is clear from first principles, as non-isolated spikes are triggered only partly by the stimulus, and partly by the timing of the previous spike. Understanding the complex interactions between spikes and the stimulus to create spike patterns is obviously an important next step in understanding the neural code, but the logical starting point is to establish the causes of isolated spikes. If we attempt stimulus reconstruction using reverse correlations averaged over all spikes, then the reconstruction will be incorrect; similarly, predicting spike timing without knowledge of the recent spiking history is generally impossible. A step along the path of treating spike interactions explicitly has recently been made in the modelling of retinal ganglion cell spike trains .

The separation of spikes into isolated and non-isolated also allows a more efficient use of the data—ensemble averages will not be complicated by the inclusion of effects from interspike interaction. Thus although we are using less data, it is probably the case that we are using the data more effectively, by summing stronger covariance signals from a simpler, lower-dimensional space. While in a simulation the amount of data can be unbounded, and we have used up to millions of spikes to demonstrate our points here, we note that the significant effects all emerge very clearly after $`10^4`$ spikes. This can be seen from the evolution of the eigenvalues and silence energy fractions with increasing spike count (Figs. 2, 6 and 7). This number is in a range often accessible to the experimenter.

It should be noted that the need to consider the effect of silence arises even in the standard analysis of non-isolated spikes. Because all spiking neurons have a refractory period, there is always an implicit silence or isolation constraint, albeit sometimes of only a few milliseconds. The second spike-associated mode identified here (Fig. 8) is localized precisely in the disallowed interval, and in fact emerges almost identically in the unqualified covariance analysis. Its contribution to the highly peaked STA is very significant, since both the spike-conditional variance and mean along this stimulus dimension are substantially altered from the prior. Recall that this mode is related to a derivative condition expressing upward threshold crossing, but it is extended, as there is always a finite period of silence prior to a spike. An analogous mode should arise for all spiking neurons, though its shape will depend strongly on the shape of the spike-generating filter or filters. Reduced stimulus variance in the direction of this derivative-like feature should not be thought of as part of the cause of a spike; it is a necessary concomitant to spike generation conditioned on the first filter alone. This inevitable derivative-like condition is relevant to the inverse problem of stimulus reconstruction from a spike train, but not to the forward problem of causally predicting spike times from the stimulus. Notice that this implies that, in general, convolving the stimulus with the STA—either unqualified or for isolated spikes only—is not an ideal way to predict spike times.

In summary, we have shown that standard white noise analysis is unable to recover the known filter properties of the leaky integrate-and-fire neuron. We show that this is due to the influence of interactions between spikes, and we demonstrate that by removing this influence, requiring spikes to be ‘isolated’, we are able to recover the computation performed by the neuron. However, in so doing, we introduce a complication: we must address not only the isolated spikes, but the silences that precede them. Silences produce complex structure in the covariance analysis, which is nonetheless statistically independent and clearly distinguishable from structure tied to the spikes themselves. In a sense, the distinction between spike-related and silence-related features is artificial, for silence is simply the complement to spiking: silence is structured by the absence of features which would produce a spike. The distinction is useful, however, in allowing us to use the isolated spikes to identify the spike-generating feature or features unambiguously.

We have shown that failing to consider spikes in isolation leads to distortion in the filters obtained through standard white noise reverse correlation. Our conclusions are relevant not just for the simple case of the leaky integrate-and-fire neuron, but for any system in which the presence of a spike has an influence on the generation of subsequent spikes. The new isolated spike method presented should improve the accuracy of white noise analysis of neural spike trains with interspike interactions. In some cases, it may significantly change our picture of what stimulus features a neuron is sensitive to.

> We would like to thank William Bialek for his advice, discussion, and encouragement to study this problem.


format: text

text:
# On the Selection of Photometric Planetary Transits

## 1. Introduction

In recent years the discoveries of large numbers of planets via high precession radial velocity studies (Marcy & Bulter 2000) has fostered a surge of activity aimed at the discovery of low mass companions to nearby stars. However, with the radial velocity technique alone it is not possible to uniquely determine the mass of a planetary candidate since the orbital inclination ($`i`$) is poorly determined. The uncertainty in $`i`$ leads to a $`sin(i)`$ degeneracy in the mass of the orbiting object. In contrast, planetary transit searches measure light curves which put strong constraints on the inclinations of the planetary orbits. This information can be combined with radial velocity measurements to determine a transiting object’s mass and mean density. Furthermore, once a transiting planet or brown dwarf has been confirmed by radial velocity measurements, additional high accuracy spectral and photometric observations make it possible to detect the presence of specific chemical elements in a transiting planet’s atmosphere (Seager & Sasselov 2000, Brown 2001, Charbonneau et al. 2002).

The radial velocity survey based discovery of a planetary companion to HD 209458 and the subsequent discovery of a planetary transit (Henry et al. 2000, Charbonneau et al. 2000, Mazeh et al. 2000) has led to a great deal of interest in the detection of extra-solar planets via photometric transits. The planetary object, HD 209458b, has a mass of 0.69MJ, a radius of 1.35RJ, and an orbital distance of 0.047AU (Cody & Sasselov 2002). A number of other extra-solar planets have been discovered with small separations from their parent stars in radial velocity surveys (51-Peg-b, $`\tau `$-Boo-b, HD 187123b, etc.), but no photometric transits have been observed for these.

Within the last year planetary transit candidates have been discovered in the data from the OLGE-III (Udalski et al. 2002a,b) and the EXPLORE (Mallén-Ornelas et al. 2002) projects as well as the Vulcan campaign (Jenkins et al. 2002). However, when considering whether a photometric transit is due to a planet or a small star, such as a late M-dwarf, it is necessary to know the mass of the transiting object. This mass can be determined with radial velocity measurements by using large telescopes with high resolution spectrographs (such as VLT or Keck). However, the determination of the mass of a single object may require a number of observations taken on multiple nights. As there are currently more than 60 planetary transit candidates and many hundreds more expected from space missions and other searches, it will soon be impractical to measure the radial velocity profiles of all candidates.

A method of selecting planetary transit candidates has been put forward by Seager & Mallén-Ornelas (2002). Among other things, this selection favors planets that have circular orbits, unblended parent stars and produce flat bottomed eclipse shapes. However, it is clear that any of these criteria may be restrictive against bona-fide planet transits. Circular orbits seem probable for transiting planets since they are strongly influenced by gravitational tidal forces due to their close proximity to their parent stars. However, this may not always be true since significant eccentricity has been observed among the planets found in radial velocity surveys (Udry et al. 2002). The second selection criteria can also fail to produce all good candidates because in some cases (such as OGLE-III) dense stellar fields are monitored to maximize the planet detection efficiency. Furthermore, many blends are likely to occur because of multiple star systems (Mallén-Ornelas et al. 2002). The selection of flat bottomed light curves excludes planets transiting near stellar limbs. These may not exhibit any flat regions whatsoever. Although the number of transits expected at these inclinations is generally small, the exact fraction of limb transits is strongly dependent on the as-yet-unknown sizes of tightly orbiting planets. Furthermore, the effect of limb darkening is strong in photometric bands blueward of I-band. Thus, a large fraction of bona-fide planets may not exhibit flat bottomed eclipses. The effect of limb darkening on the shape of planetary transits is clear from the HST data of planetary transit HD 209458 (Brown et al. 2001).

The selection of good planetary candidates can be improved if one knows the size of the parent star $`R_{}`$. This can be estimated with some accuracy by determining the spectral type of the parent star. Since the fractional drop in flux during the eclipse, $`\mathrm{\Delta }F`$, is related to the stellar and planetary radii by,

$$\mathrm{\Delta }F=\left(\frac{R_p}{R_{}}\right)^2,$$
(1)
a determination of the candidate planet’s radius $`R_p`$ can be obtained. However, limb darkening changes the measured $`\mathrm{\Delta }F`$ by a small amount which depends on the stellar type and the passband used.

The evolutionary theory of isolated planets by Baraffe et al. (1998, 2002) predicts that evolved extra-solar planets should have radii similar to that of Jupiter. Dreizler et al. (2002) used these models to select the best transit candidates from Udalski et al. (2002a,b), for follow-up with radial velocity measurements based on their sizes. However, currently there are little empirical data about the actual radii of gaseous planets at small separations from stars. For example, presently the only extra-solar planet with a measured radius is HD 209458b and as Dreizler et al. (2002) conceded, their adopted model is inconsistent with HD 209458b. It would seem that evolutionary models of isolated planets are poor representations of real planets. The giant extra-solar planet models of Burrows et al. (2000), Bodenheimer et al. (2001) and Guillot & Showman (2002) suggest that large extra-solar planets can exist at small distances from main-sequence stars. Guillot & Showman (2002) point out that, if only 1% of the flux incident on HD 209458b from its parent star was transformed to kinetic energy in the planetary atmosphere, a planet can maintain its size. In such models the incident radiation suppresses the typical energy loss by radiative diffusion and hence the contraction of such planets. As many planets have been found at small distances from their host stars in radial velocity surveys, it seems that gaseous planets at small orbits could have larger radii than HD 209458b. Furthermore, it seems possible that transiting brown dwarfs may also have large radii. In light of the current level of uncertainty about the radii of extra-solar planets and brown dwarfs in tight orbits, it would be useful to select candidates based on tracers for the transiting objects’ masses as well as their sizes. One tracer for the mass of the transiting companion in tight binary systems is the presence of gravity darkening effects in eclipsing binary light curves.

In the next section we outline the nature of the ellipsoidal variations observed in binary light curves due to gravity darkening. In §3, we show how the effect can be used to select bona-fide planet/brown dwarf candidates. In §4, we test some public photometric transit data for the presence of the gravity darkening effect. Following this in §5, we discuss the problems with planetary transit selection due to photometric blending. A partial solution to the blending problem will be discussed in §6, followed by concluding remarks.

## 2. Gravity darkening

When two stars exist in a tight orbit, their gravitational potential stretches the envelopes of the stars into ellipsoids (von Zeipel 1924). These distortions leads to sinusoidal modulations of the kind observed in the light curves of stars such as W UMa binary systems (Kitamura & Nakamura 1988). When these binary systems rotate their brightnesses vary depending on the observed luminosity and cross section of the stars. The temperature of the surface elements is seen to act in proportion to the effective gravity, such that

$$T_{eff}^4g^{\beta _1}.$$
(2)
It is clear that the extent of the gravity darkening effect observed depends on the exponent $`\beta _1`$. This parameter varies strongly between radiative stars (where it is 1) and convective stars (where it is typically around 0.25). The exact value of $`\beta _1`$ varies with the stellar mass and has been modeled for a range of stellar parameters by Alencar & Vaz (1997) and Claret (1998, 2000). However, there is only a small amount of empirical data available, such as Rafert & Twigg (1980).

The sinusoidal modulation due to ellipsoidal geometry of the stars in a binary system varies at twice the orbital frequency. Since the distortion of the primary is independent of whether or not the secondary is luminous, it is an ideal tracer for an unseen massive companion to a luminous star (Beech 1985). Furthermore, when both stars have similar luminosities the reflection effect can also be observed in the light curves. The ellipsoidal effect can be observed in binary systems even if an eclipse is not observed, because the ellipsoidal effect can be observed over a larger range of orbital inclinations than a transit. In the next section, we will show that, although accurate values of $`\beta _1`$ are yet to be determined empirically for many stellar types, even the slightest signal of a modulation can point to a secondary object with a significant mass.

## 3. Selecting planetary candidates

Typical planetary transits are expected to last a few hours for orbital periods of a few days. These dips in flux during the transit will be $`1\%`$. Therefore, the detection of a planet transit requires a well sampled light curves with high accuracy data. Transit detections can be made by either making frequent observations of candidate stars (dedicated transit surveys), or by folding light curves and searching for significant periods in data taken less frequently over a longer time (such as microlensing surveys). Because planet transits are short, roughly 95% of the photometry points taken in a survey will occur outside the transit region for any given light curve. These measurements are of little use to the determination of a planet’s size or shape. However, these data make it possible to measure variations in the baseline of the light curve to much greater accuracy than the transit depth.

The OGLE-III transit search data set consists of 800 photometry points (Udalski et al. 2002a,b) with an accuracy of better than 1.5%. A simple Poisson argument would suggest that at a 3 sigma level, a modulation amplitude of $`1.6\times 10^3`$ magnitudes should be measurable. For the Explore-I project the best candidates ($`15I17`$) have 1600 photometry points with better than 0.5% photometry. In this case, we expect that a variation should be detectable at the level of $`3.75\times 10^4`$ magnitudes. However, in reality it is likely that there are systematic contributions to the uncertainty in the photometry points. Systematic errors can occur due to observations of a field being taken nightly at a similar airmass range. The observed scatter will not behave in a Poissonian way, so the actual measured baseline uncertainties will always exceed this limit. Indeed, many of the transit light curves presented by the OGLE-III project appear to exhibit time dependent fluctuations which are either due to systematic noise contributions or real variability of the parent stars.

In Figure 1 we present the modulation amplitudes expected due to the gravity darkening effect. These values have been determined using the Nightfall111See http://www.lsw.uni-heidelberg.de/$``$rwichman/Nightfall.html for details eclipsing binary analysis program. This program can be used to model and fit the shapes of eclipsing binary light curves, line profiles and radial velocities, and includes effects such as limb darkening, star spots, reflections and third light . The values for the gravity darkening exponent $`\beta _1`$ are taken from the non-grey models of Alencar & Vaz (1997). For the curves in this particular example we have taken an orbital period of 2 days. This is consistent with a large number of the OGLE-III candidates. We expect that the spectral types of primary stars will generally be within the range of types shown on Figure 1. For example, the OGLE-III candidates measured by Dreizler et al. (2002) have spectral types ranging from A3 V to K4 V. Stars with spectral types later than K5 V are very faint and are, therefore, difficult to photometer accurately, while types earlier than A5 V (and giant types) will only have small planetary transit signals due to their sizes.

Current planetary transit surveys are attempting to find planets and brown dwarfs with masses in the range $`\mathrm{M}_2<0.08\mathrm{M}_{\mathrm{}}`$. From Figure 1 we see that (for most types of primary stars), the detection of modulations greater than $`>0.5`$ mmag in amplitude is sufficient to reject candidates for orbital periods longer than 2 days. For larger orbital periods, the separations of the stars and candidate planets are greater and the gravity induced modulations are smaller. Thus, the presence of measurable modulations at longer periods provides even stronger evidence that transiting objects have significant masses. For the most accurately measured stars it may be possible to limit the secondary’s mass to within the brown dwarf mass range $`<0.08\mathrm{M}_{\mathrm{}}`$ using photometry alone. This will be very useful for reducing the number of possible planetary candidates in space missions such as Kepler.

The transits of giant-type stars by main-sequence stars can be distinguished from those of planets eclipsing main-sequence stars, by using the parent star’s stellar density derived from the transit parameters (Seager & Mallén-Ornelas 2002). These stars can also be distinguished with spectra or multi-band photometry (Bessell & Brett 1988). Furthermore, the ellipsoidal effect is larger for giant stars than main-sequence stars since they are more readily distorted due to their low surface gravities.

If a candidate planetary transit is due to the partial eclipse of two normal stars, the observed gravity darkening effect may be large even though the observed eclipse may be very small. In some cases an observed ellipsoidal effect could also be due to the presence of additional massive object which is not the transiting object. Such a third object would have an orbital period different from that of the transiting object and thus would appear as a source of correlated noise in the phased transit light curve. In other cases a transit light curve may exhibit sinusoidal modulations due to effects other than the ellipsoidal effect. For instance, one can imagine cases where an observed modulation was due to the concentrations of star spots on one hemisphere of a star. However, this is unlikely to properly mimic the gravity darkening effect since this would require that the rotational period of the star was exactly twice the orbital period of a planet. Furthermore, the number and distribution of spots is likely to change over time.

Because there is some uncertainty in the value of $`\beta _1`$, the exact amplitudes presented in Figure 1 are also quite uncertain. In any case, it seems prudent to reject candidate planetary transits exhibiting sinusoidal modulations whether they are due to a massive transiting object or some other phenomena.

## 4. Examination of the OGLE-III data

To test whether the ellipsoidal effect can really be used to select planetary transit candidates we retrieved the freely distributed OGLE-III transit data from the OGLE-III archive222ftp://bulge.princeton.edu/ogle/ogle3/transits. We phased these light curves with the most recently derived periods which are given on the OGLE-III web site333http://bulge.princeton.edu/$``$ogle. Of the OGLE-III transit candidates, numbers 43 to 46 do not have known periods so could not be phased. On examination of the transit candidate OGLE-TR-39 we found it to have, not only a strong sinusoidal modulation but additional features visible at a period of 2.44565 days (three times the OGLE-III period). The presence of this additional signal would invalidate any sinusoidal fit, so we did not analyze this candidate. For each of the remaining set of light curves we subtracted the data during the period of the transit. Simply fitting the data with the inclusion of the transit dip would clearly bias the results toward higher values. In the absence of any ellipsoidal or other systematic effects, the data outside the transit region should be constant.

Many of the OGLE-III transit light curves show the clear signs of systematic trends with time. The exact origins of these features are unknown. However, such features are common in the photometry of crowded fields with large cameras. These occur because of blending between neighbouring stars, flat fielding errors, inaccurate airmass corrections, changes in transmission (due to atmospheric dust), etc. To obtain a reliable estimate of whether the effect of gravity is significant in any light curve, it is necessary to have a good estimate of the photometric scatter. In this analysis we attempted to determine the significance of the ellipsoidal variations in the presence of the real systematic noise. To do this, firstly, we fitted the phased light curves with the ellipsoidal modulation approximated by $`\mathrm{a}_1\mathrm{cos}(2\varphi )`$ plus a constant (where $`a_1`$ is the gravity darkening amplitude). We note that Beech (1985) and others have shown that, although additional sinusoidal terms exist, they are small in comparison to the $`\mathrm{cos}(2\varphi )`$ term. The initial fits resulted in large $`\chi ^2`$ values and very significant sinusoidal amplitudes because of the underestimated uncertainties in the data points. Next, we subtracted the fitted values from the data to remove the gravity darkening trend. Based on a Gaussian uncertainty model, we determined the time averaged uncertainty in the residuals for each light curve. We scaled the error bars to match the observed scatter and determined the parameter uncertainties in the fit terms. In many cases the fit $`\chi _r^2`$ value remained high because of real variability structure within the light curves. As an additional test of the significance of the gravity darkening term we perform an F-distribution test for the significance of this extra parameter. That is, we compared the fit $`\chi _1^2`$ value for a constant baseline with the $`\chi _2^2`$ value for a constant term plus sinusoid. The F-test statistic is given by

$$F=\frac{\chi _1^2\chi _2^2}{\chi _{r2}^2},$$
(3)
where $`\chi _{r2}^2`$ is the reduced chi-square value of the fit with the extra term. This statistic obeys the F-distribution for $`\nu _1=1`$ and $`\nu _2=Nm`$ ($`N`$ points, $`m`$ parameters) degrees of freedom. An F-distribution value larger than 6.668 (for 750 data points) is expected in only 1% of experiments and larger than 10.91 for 0.1%. The larger this value the more significant the decease in the fit $`\chi ^2`$ is with the additional term. In Table 1 we present the fit values for the OGLE-III transit candidates with amplitudes, $`a_1`$, greater than zero at a 3.5 sigma level and F-statistic values greater than 6.668. We believe that these transit candidates are the most likely to have stellar mass secondaries.

We note that the planetary transit candidate OGLE-TR-3 was selected by Dreizler et al. (2002) from the OGLE-III candidates as a likely planetary transit. This object exhibits a sinusoidal modulation at the $`4\sigma `$ level. The light curve also appears to exhibit the presence of a second dip in the phased light curve, although this dip is not very significant ($`2\sigma `$). The presence of a second dip in a binary light curve is the clear sign that transiting object is luminous. In Figure 2 we present the fits to two of the OGLE-III planetary transit candidates, OGLE-TR-18 and OGLE-TR-30. The presence of a sinusoidal modulation is clearly seen in these light curves.

## 5. Blending

A common type of contaminant in searches for planetary transits are cases where an ordinary eclipsing binary star is blended with one or more additional stars. The light curves of these candidates exhibit a dip in flux where one of the stars is eclipsed. This may mimic a planetary transit as only a small dip in the group flux is observed. The additional flux from the blended stars causes the size of eclipsed object to be underestimated from the light curve. As an example, take an M5 V spectral type star eclipsing a K2 V type star. Without any blended flux the complete stellar eclipse would be too long and too deep to be due to a planet. However, if these stars were blended with the light from an F5 V type star, the transit time and the dip in flux would be consistent with a $`1.3R_J`$ planet transiting the F5 type star near its limb. In this example we have assumed that the F5, K2 and M5 stars have the following typical parameters, $`\mathrm{M}_\mathrm{V}=`$ ($`3.5`$, $`6.4`$, $`12.3`$) and $`\mathrm{R}_{\mathrm{}}=`$ ($`1.3`$, $`1.1`$, $`0.798`$), respectively.

Seager & Mallén-Ornelas (2002) have shown that the percentage of blended events can be derived statistically by using the ratio of the flat part of the eclipse to the total transit time. However, this statistical information does not determine which of the light curves are blended in a given survey. Alternatively, if a spectrum has been taken of the parent star, it is possible to assume an approximate mass and radius based on the spectral type. With the additional information it is possible to compare the density of the parent star with that derived from transit parameters, timescales (total and flat), depth, and orbital period (Seager & Mallén-Ornelas 2002). However, this process is complicated when the blended star is much brighter than the stars involved in the eclipse. In such cases, the measured spectra and density will be for the blended star rather than the eclipsing stars.

In principal it is also possible to determine the presence of blending astrometrically since the location of the photo-center of the blended system will vary from the location of the eclipsing stars. The location of the eclipsing system can be found by using difference image analysis (Alcock et al. 1999). This technique allows the determination of the photo-center of a variable object or binary system by subtracting the constant flux component (the blended baseline flux). An accurate location can be found by combining difference images taken at times during the eclipse. The position of the eclipse is the location of the residual flux. A significant offset between the position of the residual and the photo-center of the blended system can prove that candidate is blended (Alcock et al. 2000, 2001). However, in some cases the separation between the blended centroid and the eclipse centroid will be indistinguishable. Furthermore, many surveys use very large pixel sizes to survey large areas of sky. In such cases astrometry is probably impractical.

The gravity darkening effect occurs whether or not a planetary transit candidate is blended. However, the amplitude of the observed ellipsoidal light variations depends on the brightness of the eclipsing stars relative to the additional flux component. In our example above, the ellipsoidal modulation is very small for periods longer than a day and is washed out by the flux from the much brighter blended star. If the blending is caused by an association of stars, the number of planets is limited to the few stable orbital configurations in such systems. Alternatively, if the additional blended flux comes from the transiting object, such as when the planet candidate is due to a partial stellar eclipse, ellipsoidal variations will be detectable with good photometry.

## 6. Using colour information

In section 5 we demonstrated how a blended stellar system can mimic a planetary transit. We will now consider how colour information can be used to select bona-fide transits. Our blending example consists of F5, K2 and M5 type stars. An F5 V star has a $`(VI)_{KC}`$ colour of $`0.53`$ while a K2 V star has a colour of $`0.98`$. The difference between the colour of these two stars will led to a variation in the transit depth when observed in multiple passbands. For our test case the maximum observed transit depths will be $`0.011`$ magnitudes in $`V`$-band and $`0.007`$ magnitudes in $`I`$-band. Although this difference is very small, it is a large fraction of the eclipse depth and is measurable at the accuracy levels achieved in current transit surveys.

The colours of stars on the main-sequence vary in relation to their sizes. The difference in the transit-depths observed in two passbands will increase as the difference in the sizes of the blended stars increases. Therefore, as the amount of blending increases the fractional difference in observed transit depth will generally increase. The variations in transit depth are most easily observed when observations are taken at widely separated wavelengths. With this in mind, it may well be worth pursuing an observational search strategy where filters are swapped between each observation. Such a strategy can achieve a sampling rate similar to a single colour survey (depending on the observation times and passbands chosen), and gain additional information about colour evolution.

Photometric planetary transit searches require light curves exhibiting a single eclipse since the planet does not contribute to the luminosity of the system. When the orbit of a binary is non-circular there are inclinations where a partial eclipse can appear like a planetary transit. In such cases the orbital eccentricity can cause only a single eclipse to be observed. The presence of a second eclipse is also significantly reduced when there are large differences in the luminosities of the stars in a binary. Therefore, candidate planetary transits are also biased towards binary stars with large luminosity differences. If both stars are on the mainsequence they are likely to have different colours. This colour information can be used to remove candidate planetary transits due to partial stellar eclipses. The detection of a significant colour change during an eclipse is evidence of blending. Colour changes have been observed in searches for planetary transits candidates in MACHO project data (Drake et al. 2003).

Testing planetary transit candidate light curves for colour changes during eclipses may be the easiest and most robust way of finding which events are due to blending. In many cases, a few high signal-to-noise photometric follow up observations, inside and outside the eclipse, will detect the presence of blending. However, in some situations, radial velocity measurements will still be necessary to determine whether or not blending is responsible. It is clear that none of the techniques listed above can give a clear indication of whether blending is present in all cases. Still, the combination of these techniques should reduce the number of blends in planetary transit searches.

## 7. Conclusions and Discussion

We have shown that the presence of modulations in the light curves of a planetary transit candidates is likely to be due to the gravity darkening effect. This photometric information can be used as an effective way of reducing the number of spurious candidates in current and future planet searches. Unlike the transit timescale, the gravity darkening modulations are directly related to the secondary’s mass. Such a size-independent relation is important because of current uncertainty about the sizes of planets and brown dwarfs in close stellar orbits.

In many cases when very good photometry ($`<1\%`$) is obtained, it is possible to rule out transit candidates due to M-dwarfs to white dwarfs. If greater than 1000 highly accurate ($`<0.5\%`$) photometry points are obtained, it may be possible to select objects in the planet/brown dwarf regime with photometry alone. The Kepler mission is expected to obtain a photometric precision of 90 micromags (Sahu & Gilliland 2002). However, radial velocity follow-up will still be necessary to determine accurate masses of individual objects and thus to separate brown dwarfs from planets.

White dwarf secondaries are expected to mimic small extra-solar planets in transit surveys due to the transit depths they will cause. It has been shown by Marsh (2001) that microlensing of a transiting white dwarf can cause a magnification during the transit. It has also been reported by Sahu & Gilliland (2002) that the presence of this lensing can be used to break the similarity between planetary and white dwarf transits. However, as white dwarfs are expect to have masses around 0.6 $`\mathrm{M}_{\mathrm{}}`$, the gravity darkening effect will be strong. Therefore, at small orbital distances ($`<0.03AU`$) white dwarf transits will not appear similar to planets. At larger orbital distances microlensing becomes important. However, few transits are expected with long periods since the transit probability decreases as $`R_{}/a`$, where $`R_{}`$ is the primary star’s radius and $`a`$ is the orbital separation.

In this work we have shown that many of the Udalski et al. (2002a,b) OGLE-III transit candidates exhibit the presence of ellipsoidal modulations in their light curves with $`>3.5\sigma `$ significance. In the coming years many hundreds of planetary candidates are expected from on-going transit surveys. Such large numbers may necessitate the use of a lower selection threshold. However, any selection is tentative without a clear understanding of the systematic uncertainties in the light curves.

Udalski et al. (2002a,b) selected main-sequence stars as targets under the assumption that the transiting objects passed before the center part of their parent stars. Under this assumption it is possible to derive the parent star’s radius and mass from the transit light curve parameters using the equations presented by Seager & Mallén-Ornelas (2002). However, Udalski et al. (2002a,b) instead assumed that all candidates had a mass of $`1\mathrm{M}_{\mathrm{}}`$ and simply noted that there was a small scaling for parent star and planet with the true mass. Although this assertion is correct, the list of OGLE-III candidates contains many stars with radii which are inconsistent with a standard main-sequence mass-radius relationship. Since the radii of transiting objects are derived from the sizes primary stars, these values are biased by the $`1\mathrm{M}_{\mathrm{}}`$ assumption.

When a bright star is blended with an eclipsing binary system the resultant light curve can mimic the signal of a planet transiting the bright star. In such cases it may be difficult to separate bona-fide candidates from blends with spectra or based on the presence of sinusoidal baseline variations. However, it may be possible to detect the blend by measuring the depth of the transit in more than one passband.

I wish to thank the members of OGLE-III project for making their data publicly available. I also wish to thank Gabriela Mallen-Ornelas and Sara Seager for their comments which improved this paper.


text:
# Charge ordering in oxides : a conundrum solved by resonant diffraction

Present address: \]Dept. of Physics, Brookhaven National Laboratory, Upton, NY 11973, USA. grenier@bnl.gov

## Abstract

We show that in mixed-valence $`3d`$ transition metal oxides undergoing a structural transition, the low temperature phase results from an effective ordering of the charge. This arrangement and the quantitative evaluation of the atomic charges are determined by using resonant x-ray scattering experiments further analyzed with the help of ab initio calculations of the corresponding scattering factors. We have found that this reorganization concerns only a small fraction of electron and is necessary to reconcile all the experimental data.

There is an issue in understanding the charge state in the high and low temperature phase of mixed-valence transition metal compounds. Few experimental techniques can access the actual charge state of a given atom and most of the accumulated knowledge results from theoretical models 1  and from indirect results based on transport property measurements 2  or from heuristic bond-valence calculations on the crystal structure 3 . Nevertheless, the answer to very simple questions would turn out to be of large relevance in the context of metal-insulator or insulator-insulator phase transitions issued from a charge ordering (CO), chiefly if one could deduce a general pattern. For instance, given an average valence state as deduced from charge neutrality considerations of the chemical composition, would the ground state be an effective CO ? and if so how are the corresponding electrons arranged ? Are they located in the very localized $`3d`$ orbital of the transition metal, or in the extended $`4s`$ or even on the oxygen ligands ? In this paper we claim that resonant X-ray diffraction is able to answer quantitatively these questions. In our study, we have found a small but effective ordering occurring mainly in the $`3d`$ transition metal orbitals. We think that the same trend occurs in many compounds exhibiting CO which range in diversity and physical properties from colossal magnetoresistance manganites, to the astonishingly complex Verwey 4  transition in magnetite, to the CO inducing Spin-Peierls phase transition in NaV2O5 5 .

X-ray absorption spectroscopies (XAS) and resonant X-ray scattering are extensively used techniques to explore and elucidate the electronic properties of materials 6 ; 7 . In these spectroscopies the incoming photons are absorbed by promoting a core electron to an empty intermediate excited state $`f`$. When the electron subsequently decays to the same core-hole, it emits a second photon with the same energy than the incoming one. This process is virtual and preserves the coherence of the electromagnetic wave. Set in diffraction conditions the associated experimental technique is known as diffraction anomalous fine structure spectroscopy or resonant diffraction. When the absorption process is considered alone, and therefore loosing the coherence properties of the X-ray field, the experimental technique is called XAS or X ray absorption near edge spectroscopy (XANES) when considering the structure close to the edge. In this low energy range a higher sensitivity to the electronic configuration around the photoabsorber is expected because the photoelectron probes allowed electronic states very close to the Fermi level as required by the selection rules resulting from the Fermi golden rule. Indeed by playing with the polarization of the incoming (and eventually outgoing) photon beam, with the different channels of transition (dipole, quadrupole) and with the orientation of the crystal with respect to the photon polarization direction, the probability of transition becomes highly angle dependent. This phenomena is illustrated by a tensor form of the scattering structure factors, and the analysis of each term of the tensors gives access to specific projections of the electronic states. Although complete calculations are performed we will restrict to the dipole approximation and the atomic scattering factor (ASF) reads,

$`f^{}+if^{\prime \prime }`$ $`=`$ $`{\displaystyle \underset{\alpha \beta }{}}ϵ_\alpha ^oϵ_\beta ^iD_{\alpha \beta }`$ (1)
$`\alpha `$ and $`\beta `$ are the index of the coordinates in the orthogonal basis. $`ϵ^i`$ and $`ϵ^o`$ stand for the incoming and outgoing polarizations and $`D_{\alpha \beta }`$ is given by:

$`D_{\alpha \beta }`$ $`=`$ $`m_e\omega ^2{\displaystyle \underset{fg}{}}{\displaystyle \frac{\psi _g|r_\alpha |\psi _f\psi _f|r_\beta |\psi _g}{E_fE_g\mathrm{}\omega +i\frac{\mathrm{\Gamma }(E_f)}{2}}}`$ (2)
where $`E_f`$, $`E_g`$ and $`\mathrm{}\omega `$ are respectively the energies of the intermediate state, the core state and the photon; $`m_e`$ is the electron mass and $`\mathrm{\Gamma }(E_f)`$ includes the contribution of the inverse lifetime of all the states involved in the resonant process. The summation over $`f`$ is limited to the unoccupied states. In diffraction, the intensity of the reflections is proportional to the modulus squared of a weighted sum of the ASF:

$$I(\stackrel{}{Q};\mathrm{}\omega )\left|\underset{a}{}e^{i\stackrel{}{Q}\stackrel{}{R_a}}\left(f_{0a}(\stackrel{}{Q})+f_a^{}(\omega )+if_a^{\prime \prime }(\omega )\right)\right|^2$$
(3)
where $`\stackrel{}{R_a}`$ stands for the atom $`a`$ position and $`\stackrel{}{Q}`$ is the diffraction vector. The atomic scattering amplitude is the sum of the Thomson factor $`f_{0a}`$ and the ASF.

Of particular importance are the *forbidden* Bragg peaks that are extinct off resonance conditions because of the presence of glide planes or screw axis, if allowed by the space group of the given material. For these reflections, the diffracted intensity is related to the difference between ASF of the resonant (or anomalous) atoms. Therefore, these reflections are supposed to be much more sensitive to slight modifications of the crystallographic as well as the local electronic structure. In the absence of new superlattice peaks originating from a phase transition (typically ferrodistortive phase transitions), these *forbidden* Bragg peaks might yield relevant informations on the corresponding CO.

This phenomenology has been extensively used in the study of the interplay between orbital, spin and charge degrees of freedom in transition metal oxides 8 ; 9 . In particular it is specially suited to study charge 10  and orbital orderings 11 ; 12 ; 13 . Nevertheless apparently contradictory results have led to opposite conclusions on the presence or not of real CO. Recent resonant diffraction experiments at *forbidden* reflections in Fe3O4 have yielded no change in the resonant signal below and above the ordering temperature and it has been concluded that no CO occurs at the Verwey transition 14 . Conversely, resonant scattering experiments in NaV2O5 15  and La1/2Sr1/2MnO3 16  have shown that the energy line shape of some reflections retains close similarities with the derivative of the resonant X-ray scattering factor, a fact which has been interpreted as a signature for the presence of two different V (or Mn) charge states ordered in the lattice. In this letter we shall show that all these experimental evidences can actually be interpreted within a unique framework. In order to accomplish this task, we shall proceed in a quantitative analysis of the resonant diffraction experiments with the help of *ab initio* simulations. It will then be possible to evaluate the influence of the different parameters and estimate the effective charge redistribution resulting from the CO transition.

To illustrate our analysis we shall discuss very recent resonant X-ray diffraction results in NaV2O5 and all the experimental details can be found in ref. 15 . Among the numerous diffraction peaks measured as a function of the incident photon energy and under different photon polarization conditions, we select three of them to illustrate our results. The first one is the (1, 0, 0)c *forbidden* reflection which does not exhibit any change through the transition. The incoming photon polarization is along $`c`$, the outgoing one is in the scattering plane ($`\sigma \pi `$ polarization condition). The two other reflections under consideration are superlattice peaks, (15/2, 1/2, 1/4)b,c, for which the incoming and outgoing photon polarizations were kept parallel ($`\sigma \sigma `$ condition). According to the index, the polarization directions are along either the $`b`$ or $`c`$axis. The superstructure peaks appear due to the phase transition and the energy dependence of the scattered intensity exhibits peculiar features that we examine in detail below. *Ab initio* monoelectronic calculations are performed using the FDMNES code. 17 ; 18  Starting from the atomic positions and an associated electronic density, this code solves the Schrödinger or Dyson equations to determine the intermediate states probed by the photoelectron. The resulting structure factors are evaluated by taking into account the different excitation-deexcitation channels (dipole-dipole, dipole-quadrupole and quadrupole-quadrupole). Diffracted intensities are calculated with the crystallographic structures corresponding to the high ($`Pmmn`$ 19 ) and low ($`A112`$ 20 ) temperature phases with and without CO. Because the purpose of the present paper is to fit the electronic structure, it is important to define easy to handle parameters. We choose the valence orbital occupancy rates (VOOR). The radial shape of the atomic orbitals remains fixed. For any VOOR set, the total electronic density is obtained by simple superposition of the electronic orbital densities weighted by the VOOR. First-guess VOOR values corresponding to the high temperature (HT) phase are obtained by a fit on the electron density calculated by a LAPW 21  band structure procedure. These orbital occupancies, which reflect the location of the charge in space, are further refined in order to accurately reproduce the experimental curves in the low temperature phase. The vanadium and oxygen VOOR are then the only adjustable parameters of our calculations. The HT starting values are respectively 2.3 and 1.7 electrons for the $`3d`$ and $`4s`$ V atoms and 4.6 electrons for the $`2p`$ O atoms. If one uses the ionic radii as a definition of the atomic spheres, the $`4s`$ charge (average radius is 1.5Å, to be compared with the $`3d`$ average radius, 0.6Å) lies inside the oxygen atoms. Thus, the V formal charge is 2.7+ in this context. The sodium atom has no electron on its $`3s`$ valence orbital and remains 1+ in all the study. With this VOOR set, a good agreement is obtained at HT between the calculated (100) reflection, the polarized XANES line shape and the corresponding experimental spectra. To get such agreement, the finite difference procedure is used to solve the Schrödinger equation. This technique permits to avoid the muffin-tin approximation on the potential shape. That point is absolutely necessary in such distorted systems.

We now discuss the low temperature (LT) phase and the starting point is the crystallographic structure as determined in ref. 19 . The most important result is that an explicit CO is absolutely necessary to get a satisfactory agreement with the experiment (Fig. 1). The fit of the experimental data yields two equivalent solutions : *(a)* The charge exchange is located at the $`4s`$ vanadium and amounts $`\delta =\pm `$ 0.2 electron, and *(b)* it is located at the $`3d`$ vanadium and amounts only $`\delta =\pm `$ 0.04 electron. Thus the CO takes place either in the intermediate region between the oxygen and the vanadium where an atomic assignment of the electron is merely formal (see Fig. 2) or closer to the atom centers with a smaller magnitude. In any case, the bridging oxygen located in the rung of the ladder, being the atom that experiences the largest displacement at the phase transition (0.05Å towards the V5+, that is 5 to 10 times more than the other atoms), is suspected to participate to the electronic rearrangement as well. Thus, the electronic reorganization we propose does not concern isolated atoms but all the rung of the ladder, with the center of gravity of the electronic density displaced in the opposite direction than that of the bridging oxygen, i.e., toward the V4+ atom. Taking into account the experimental resolution, the choice among the two proposed solutions cannot be safely done. Nevertheless the second one with a CO on the $`3d`$ orbital is in agreement with that obtained by Suaud and Lepetit 22  on the basis of ab initio calculations using large configurations interaction method on embedded fragments. These authors found a $`\pm `$ 0.04 CO on the $`3d`$, in excellent agreement with the present study. They have also shown that the modification of the magnetic properties at the phase transition is due to the redistribution between the different $`d`$orbitals, keeping the total $`d`$ charge only slightly modified. We want to end this paragraph by underlying the extreme sensitivity of resonant diffraction to observe such a small amount of charge disproportionation.

A second important result, with straightforward implications in other CO compounds, is that the CO has very little influence on the (1, 0, 0) *forbidden* reflection whereas it is at the origin of the derivative line shape for the non integer reflections. A charge disporportionation in the vicinity of the absorbing atom greatly affects localized orbitals, and thus both the $`1s`$ and the $`3d`$ energy levels will undergo an analogous shift. However the $`4p`$ states, less localized, are hardly affected joly . Therefore the dipole-dipole contribution to the $`K`$edge of an atom experiencing a CO will be shifted towards higher or lower energies as compared to that found in the HT phase whereas the quadrupole pre-edges are little shifted during the oxidation processes. In a more general way, one has to proceed very cautiously on applying the derivative effect concept to $`L`$ and $`M`$edges if the CO is expected to take place in $`d`$ and $`f`$ final states, respectively. We have computed the effect of the associated geometrical distortion onto the ASF and we have concluded that it remains negligible ($`\pm `$ 0.1 eV) as compared to the observed energy shift and to the effect of the charge disproportionation ($`\delta E=\pm 0.8`$ eV). A first order Taylor expansion of the scattering factor is justified by the small value of $`\delta E`$. It reads as $`f^{+\delta }f+f/E\delta E`$ and $`f^\delta ff/E\delta E`$ where the $`f`$, $`f^{+\delta }`$ and $`f^\delta `$ are the ASF tensors for the vanadium atoms with the HT charge, and the $`+\delta `$ and $`\delta `$ charge disproportionation, respectively. The derivative effect condition is best achieved for diffraction wave vectors with direction and amplitude such that the putative disproportionated vanadium ASF subtract. Conversely, for the *forbidden* reflections -as the (1, 0, 0)- the derivative term remains negligible in front of the combination of the ASF tensors weighted by its phase factors, and the sum of phase factors of V-atoms of a given charge state almost cancels out.

By using the multiple scattering theory (MST) we have checked that, except for the energy shift, the ASF are very little modified through the transition. Therefore our assumption above is plainly justified and we use energy shifted HT ASF in the LT structure. It can also be noticed that the XANES spectra show little sensitivity to the charge occupancy. The reason is exactly the same as the non-sensitivity of the (1, 0, 0) reflection: XANES adds the contributions of the different atoms, the small shifts due to the CO remain smeared within the experimental resolution.

We have shown that a charge-ordering transition model quantitatively explains the resonant scattering experiments in NaV2O5. The method presented here has the advantage to allow for a direct evaluation of the actual ordered charge disproportionation in the $`3d`$-transition metal oxides. The charge ordering concerns a small fraction of an electron mainly the $`d`$orbitals even if a contribution of the $`s`$orbitals cannot be excluded from the present study. This conclusion, although heavily suspected from the results of other works, points to reconsider electronic interactions in strongly correlated electron systems and to the role of small lattice distortions in accommodating extra charges. Finally with respect to compounds such as manganites and magnetite, where resonant diffraction experiments have yielded no sign of CO, we suspect that the studied diffraction peaks are hardly sensitive to an actual CO. Our results in NaV2O5 give evidence that spectroscopic effects associated to a very small CO are not at all negligible and that they should be observable in resonant diffraction experiments if convenient meaningful peaks are measured and the results contrasted with the help of appropriate simulations codes.


text:
# The Jet Angular Profile and the Afterglow Light Curves

## 1. Introduction

Most GRB jet models consider a uniform (or ‘top hat’) jet, where the Lorentz factor, $`\mathrm{\Gamma }`$, and the energy per unit solid angle, $`ϵ`$, do not depend on the angle $`\theta `$ from the jet axis, within some finite well defined opening angle, $`\theta _j`$, and drop sharply outside of $`\theta _j`$. The possibility that GRB jets can display an angular structure, where $`ϵ`$ and $`\mathrm{\Gamma }`$ are power laws in $`\theta `$, was proposed by Mészáros, Rees & Wijers (1998). Recently, several different groups have analyzed afterglow observations within the frame work of the ‘top hat’ jet model, and have inferred a relatively narrow distribution both for the total energy output in gamma-rays (Frail et al. 2001) and in the initial kinetic energy of the relativistic outflow (Panaitescu & Kumar 2001; Piran et al. 2001). These results may alternatively be interpreted as GRB jets having a universal structure, which is intrinsically the same for all GRBs, and the observed differences between different GRBs are a result of different viewing angles, $`\theta _v`$, w.r.t the jet symmetry axis (Lipunov, Postnov & Prokhorov 2001; Rossi, Lazzati & Rees 2002; Zhang & Mészáros 2002). Whereas in the ‘top hat’ jet interpretation, the jet break time, $`t_j`$, depends mainly on the initial opening angle of the jet, $`\theta _j`$, in the universal ‘structured’ jet interpretation, $`t_j`$ depends mainly on the viewing angle, $`\theta _v`$, and the light curve is roughly similar to that for a ‘top hat’ jet with $`\theta _j=\theta _v`$ and $`ϵ=ϵ(\theta _v)`$.

While the evolution of ‘top hat’ jets and their light curves has been widely investigated, much less work has been done on ‘structured’ jets. Here we describe the main results of Kumar & Granot (2003, KG hereafter) and Granot & Kumar (2003, GK hereafter), and refer the reader to these works for more details.

## 2. The Jet Dynamics & Afterglow Light Curves

We begin with a numerical hydrodynamic calculation of the evolution of a relativistic jet (KG). The hydrodynamic equations are reduced to 1-D by assuming axial symmetry and integrating over the radial profile of the flow, thus considerably reducing the computation time. We examined initial conditions where $`ϵ`$ and $`\mathrm{\Gamma }`$ are power laws in $`\theta `$, outside of a core angle, $`\theta _c`$: $`ϵ(\theta ,t_0)=ϵ_0\mathrm{\Theta }^a`$ and $`\mathrm{\Gamma }(\theta ,t_0)=1+(\mathrm{\Gamma }_01)\mathrm{\Theta }^b`$, where $`\mathrm{\Theta }[1+(\theta /\theta _c)^2]^{1/2}`$, as well as a Gaussian profile: $`[\mathrm{\Gamma }(\theta ,t_0)1]/(\mathrm{\Gamma }_01)=ϵ(\theta ,t_0)/ϵ_0=\mathrm{exp}(\theta ^2/2\theta _c^2)`$. For simplicity, we show results only for a uniform ambient medium. The hydrodynamic evolution is shown in Fig. 1. For power law jet profiles the lateral energy transfer is small, and as long as $`\mathrm{\Gamma }(\theta =0)>`$a few, $`ϵ(\theta ,t)ϵ(\theta ,t_0)`$. For the Gaussian profile, which is the steepest, a shock forms in the lateral direction, but most of the energy still remains at $`\theta <\theta _c`$. The lateral velocity in the comoving frame is found to be $`v_\theta ^{}c/\mathrm{\Gamma }c_sc/\sqrt{3}`$ for $`\mathrm{\Gamma }1`$ (in KG we show analytically that $`v_\theta ^{}c/(\mathrm{\Gamma }\delta \theta )`$ where $`\delta \theta `$ is the angle on which $`\mathrm{\Gamma }`$ or $`ϵ`$ change significantly).

Now we examine two simple models, where either: (1) $`ϵ(\theta ,t)=ϵ(\theta ,t_0)`$, or (2) $`ϵ`$ is averaged over the region to which a sound wave can propagate (this simulates the maximal lateral energy transfer that is consistent with causality). We assume initial power law profiles, as described above, and calculate the afterglow light curves for observers at different viewing angles, $`\theta _v`$, assuming synchrotron emission (see Fig. 2). For $`(a,b)=(0,2)`$ the light curve initially rises \[before the deceleration time $`t_{\mathrm{dec}}(\theta _v)`$\], and there is no jet break, which is quite different from observations for most afterglows. For $`(a,b)=(2,2),(2,0)`$ we find a jet break at $`t_j`$ when $`\mathrm{\Gamma }(\theta _v)\theta _v^1`$. For $`(a,b)=(2,2)`$ the value, $`\alpha _1`$, of the temporal decay slope $`\alpha d\mathrm{log}F_\nu /d\mathrm{log}t`$ at $`t<t_j`$, increases with $`\theta _v`$, while $`\alpha _2=\alpha (t>t_j)`$ decreases with $`\theta _v`$. This effect is more prominent in model 1, and appears to a lesser extent in model 2, but is prominent in the light curves from the numerical calculation of the jet dynamics (KG). This suggests that $`\delta \alpha =\alpha _1\alpha _2`$ should increase with $`t_j`$, which is not supported by observations. For $`(a,b)=(2,0)`$, there is a flattening of the light curve just before the jet break, for $`\theta _v>3\theta _c`$. Again, this effect is larger in model 1, compared to model 2. The light curves for the simulated jet dynamics (KG) show a reasonably sharp jet break at $`\theta _v>(23)\theta _c`$ and flattening of the light curve just before the jet break becomes strong at $`\theta _v>(57)\theta _c`$. This leaves a factor of $`3`$ in $`\theta _v/\theta _c`$ for which there is a sharp jet break not preceded by a flattening, as is typically observed in afterglow light curves. Since a larger range of $`\theta _v`$ ($`\theta _j`$ for ‘top hat’ jets) is inferred ($`220^{}`$), this poses a problem for a universal jet structure.

## 3. Conclusions

We have described a numerical scheme for calculating the dynamics of structured relativistic jets. The lateral velocity in the local frame is found to be smaller than the sound speed, except for very sharp jet angular profiles. This causes the energy per unit solid angle, $`ϵ`$, to remain close to its initial value in the first few days, suggesting that a simple model (model 1) where $`ϵ`$ is constant in time, and each segment of the jet evolves as if it were part of a spherical flow, provides a good approximation for the jet dynamics and light curves for $`t<`$a few days.

A universal jet profile, where all GRB jets are intrinsically identical, and only our viewing angle changes between different GRBs, has difficulty explaining all the afterglow light curves observed to date. This was demonstrated for $`ϵ(t_0)\theta ^2`$, which is needed to reproduce the recent results of a roughly constant energy in GRB jets. Furthermore, we find that the jet break for structured jets in an external density $`r^2`$ are much smoother compared to a uniform density, in agreement with the result of Kumar & Panaitescu (2000) for a ‘top hat’ jet. Finally, the light curves for a Gaussian initial jet profile are found to be similar to those for a ‘top hat’ jet, and (for $`\theta _v<\theta _c`$) produce sharp jet breaks, that are compatible with afterglow observations.

## References

Frail, D. et al., 2001, ApJ, 562, L55

Granot, J., & Kumar, P. 2002, submitted to ApJ (astro-ph/0212540) (GK)

Kumar, P., & Granot, J. 2002, submitted to ApJ (KG)

Kumar, P., & Panaitescu, A. 2000, ApJ, 541, L9

Lipunov, V.M., Postnov, K.A., & Prokhorov, M.E. 2001, Astron. Rep., 45, 236

Mészáros, P., Rees, M.J, & Wijers, R.A.M.J. 1998, ApJ, 499, 301

Panaitescu, A., & Kumar, P. 2001, ApJ, 560, L49

Piran, T., Panaitescu, A., Kumar, P., & Piro, L. 2001, ApJ, 560, L167

Rossi, E., Lazzati, D., & Rees, M.J. 2002, MNRAS, 332, 945

Zhang, B., & Mészáros, P. 2002, ApJ, 571, 876