text:
# The Faint Optical Afterglow and Host Galaxy of GRB 020124: Implications for the Nature of Dark Gamma-Ray Bursts

## 1. Introduction

One of the main observational results stemming from five years of $`\gamma `$-ray burst (GRB) follow-ups at optical wavelengths is that about $`60\%`$ of well-localized GRBs lack a detected optical afterglow, (“dark bursts”; Taylor et al. 2000; Fynbo et al. 2001; Reichart & Yost 2001; Lazzati, Covino, & Ghisellini 2002). In some cases, a non-detection of the optical afterglow could simply be due to a failure to image quickly and/or deeply enough. However, there are two GRBs for which there is strong evidence that the optical emission should have been detected, based on an extrapolation of the radio and X-ray emission (Djorgovski et al. 2001a; Piro et al. 2002). One interpretation in these two cases is that the optical light was extinguished by dust, either within the immediate environment of the burst or elsewhere along the line of sight (e.g. Groot et al. 1998). An alternative explanation is a high redshift, leading to absorption of the optical light in the Ly$`\alpha `$ forest. However, the redshifts of the underlying host galaxies of these GRBs are of order unity (Djorgovski et al. 2001a; Piro et al. 2002).

Several authors have recently argued that a large fraction of the dark bursts are due to dust extinction within the local environment of the bursts (e.g. Reichart & Yost 2001; Lazzati et al. 2002; Reichart & Price 2002), but other scenarios have also been suggested (e.g. Lazzati et al. 2002). Moreover, it has been noted that regardless of the location of extinction within the host galaxy, the fraction of dark bursts is a useful upper limit on the fraction of obscured star formation (Kulkarni et al. 2000; Djorgovski et al. 2001b; Ramirez-Ruiz, Trentham, & Blain 2002; Reichart & Price 2002).

However, from an observational point of view, we must have a clear understanding of the diversity of afterglow properties before extracting astrophysically interesting inferences from dark bursts. For example, afterglows which are faint or fade rapidly (relative to the detected population) would certainly bias the determination of the fraction of truly obscured bursts. In this vein, Fynbo et al. (2001), noting the faint optical afterglow of GRB 000630, argued that some dark bursts are due to a failure to image deeply and/or quickly enough, rather than dust extinction.

Here we present optical and radio observations of GRB 020124, an afterglow that would have been classified dark had it not been for rapid and deep searches. Furthermore, GRB 020124 is an example of an afterglow, which is dim due to the combination of intrinsic faintness and a relatively fast decline, and not strong extinction.

## 2. Observations

### 2.1. Ground-Based Observations

GRB 020124, localized by the HETE-II satellite on 2002, Jan 24.44531 UT, had a duration of $`70`$ s and a fluence ($`6400`$ keV) of $`3\times 10^6`$ erg cm-2 (Ricker et al. 2002). Eight minutes after receiving the coordinates111This corresponds to 1.6 hours after the burst detection. we observed the error box with the dual-band ($`B_M`$, $`R_M`$) MACHO imager mounted on the robotic 50-in telescope at the Mount Stromlo Observatory (MSO). We also observed the error box with the Wide-Field Imager on the 40-in telescope at Siding Spring Observatory (SSO). We were unable to identify a transient source within the large error box (Price, Schmidt, & Axelrod 2002).

We subsequently observed the error box refined by the Inter-Planetary Network (Hurley et al. 2002) with the Palomar 48-in Oschin Schmidt using the unfiltered NEAT imager. PSF-matched image subtraction (Alard 2000) between the MACHO and NEAT images revealed a fading source (Price et al. 2002), which was $`R18`$ mag at the epoch of our first observations, and not present in the Digitized Sky Survey. Two nights later we observed the afterglow using the Jacobs CAMera (JCAM; Bloom et al. 2002b) mounted at the East arm focus of the Palomar 200-in telescope (Bloom et al. 2002a). The position of the fading source is $`\alpha `$(J2000)=9h32m50.78s, $`\delta `$(J2000)=$`11^{}31^{}10.6^{\prime \prime }`$, with an uncertainty of about $`0.4`$ arcsec in each coordinate (Fig. 1).

Using the Very Large Array (VLA222The VLA is operated by the National Radio Astronomy Observatory, a facility of the National Science Foundation operated under cooperative agreement by Associated Universities, Inc.) we observed the fading source at 8.46 and 22.5 GHz (see Table 3). We detect a faint source, possibly fading, at 8.46 GHz located at $`\alpha `$(J2000)=9h32m50.81s, $`\delta `$(J2000)= $`11^{}31^{}10.6^{\prime \prime }`$, with an uncertainty of about $`0.1`$ arcsec in each coordinate. Given the positional coincidence between the fading optical source and radio detection we suggest this source to be the afterglow of GRB 020124.

The optical images were bias-subtracted and flat-fielded in the standard manner. To extract the photometry we weighted the aperture with a Gaussian equivalent to the seeing disk (”weighted-aperture photometry”), using IRAF/wphot. The photometric zero-points were set through photometry of calibrated field stars (Henden 2002) with magnitudes transformed to the appropriate system (Bessell & Germany 1999; Smith et al. 2002). The photometry is summarized in Tab. 1

### 2.2. Hubble Space Telescope Observations

We observed the afterglow with the Hubble Space Telescope (HST) using the Space Telescope Imaging Spectrograph (STIS) on 2002 Feb. 11.09, 18.30, and 25.71 UT (Bloom et al. 2002a), as part of our large HST Cycle 10 program (GO-9180, PI: Kulkarni). The HST observations consisted of 750–850 sec exposures. The HST data were retrieved after “On-The-Fly” pre-processing. Using IRAF we drizzled (Fruchter & Hook 2002) each image onto a grid with pixels smaller than the original by a factor of two and using pixfrac of 0.7.

We found an astrometric tie between the HST and JCAM images using IRAF/geomap with nine suitable astrometric tie objects in common between the images. The rms of the resultant mapping is 133 mas (RA) and 124 mas (Dec). Using this mapping and IRAF/geoxytrans we transfered the afterglow position on the JCAM image to the HST images. The rms of the transformation is 604 mas (RA) and 596 mas (Dec), and is dominated by the uncertainty in the JCAM position.

The source ”S1” (Fig. 2) coincides with the afterglow position within the astrometric uncertainty. We performed differential photometry at the position of S1 by registering the images of epochs 1 and 2 using a cross-correlation of a field of size 10 arcsec centered on S1 (using IRAF/crosscor and shiftfind). We used IRAF/center and the FWHM of a relatively bright point source (”PSF star”; Fig. 1) to fix the position of S1 in each of the final images, and to determine the uncertainty in the position.

We photometered the source (and the PSF star) in epoch 1 using IRAF/phot, and determined a count-rate of $`0.0921\pm 0.013`$ e- s-1 (corrected by $`17\%`$ for the loss of flux from an infinite aperture radius). Using IRAF/synphot and assuming a source spectrum of $`f_\lambda \lambda ^{1.4}`$ (see below), we find that the source was $`R=28.55_{0.14}^{+0.16}`$ mag at the time of epoch 1. The photometry of the three epochs is summarized in Tab. 2. Please note that this more careful analysis supersedes our preliminary report (Bloom et al. 2002c).

There are no obvious persistent sources within 1.75 arcsec of the OT down to $`R29.5`$ mag. To date, all of the GRBs localized to sub-arcsecond accuracy have viable hosts brighter than this level within $`1.3`$ arcsec of the OT position (Bloom, Kulkarni, & Djorgovski 2002). The faintest host to date is that of GRB 990510, $`R28.5`$ mag (z=1.619; Vreeswijk et al. 2001). Thus, the host of GRB 020124 may be at a somewhat higher redshift; however, $`z4.5`$ since the afterglow was detected in the $`B_M`$ filter.

## 3. Modeling of the Optical Data

In Figure 3 we plot the optical lightcurves of GRB 020124, including a correction for Galactic extinction, $`E(BV)=0.052`$ mag (Schlegel et al. 1998). The optical lightcurves are usually modeled as $`F_\nu (t,\nu )=F_{\nu ,0}(t/t_0)^\alpha (\nu /\nu _0)^\beta `$. However, as can be seen in Fig. 3, the $`R`$-band lightcurve cannot be described by a single power law. Restricting the fit to $`t<2`$ days we obtain ($`\chi _{\mathrm{min}}^2=15`$ for 14 degrees of freedom) $`\alpha _1=1.60\pm 0.04`$, $`\beta =1.43\pm 0.14`$, and $`F_{\nu ,0}=2.96\pm 0.25`$ $`\mu `$Jy; here $`F_{\nu ,0}`$ is defined at the effective frequency of the $`R_M`$ filter and $`t=1`$ day. For $`t>2`$ days we get $`\alpha _21.9`$.

To account for the steepening we modify the model for the $`R`$-band lightcurve to:

$$F_\nu (t,\nu )=F_{\nu ,0}(\nu /\nu _0)^\beta [(t/t_b)^{\alpha _1n}+(t/t_b)^{\alpha _2n}]^{1/n},$$
(1)
where, $`\alpha _1`$ is the asymptotic index for $`tt_b`$, $`\alpha _2`$ is the asymptotic index for $`tt_b`$, $`n<0`$ provides a smooth joining of the two asymptotic segments, and $`t_b`$ is the time at which the asymptotic segments intersect. We retain the simple model for the $`R_M`$ and $`B_M`$ lightcurves since they are restricted to $`t0.13`$ days (i.e. well before the observed steepening).

We investigate two alternatives for the observed steepening in the framework of the afterglow synchrotron model (e.g. Sari, Piran, & Narayan 1998). In this framework, $`\alpha _1`$, $`\alpha _2`$, and $`\beta `$ are related to each other through the index ($`p`$) of the electron energy distribution, $`N(\gamma )\gamma ^p`$ (for $`\gamma >\gamma _{\mathrm{min}}`$). The relations for the models discussed below, as well as the resulting closure relations, $`\alpha _1+b\beta +c=0`$, are summarized in Tab. 4.

### 3.1. Cooling Break

The observed steepening, $`\mathrm{\Delta }\alpha \alpha _2\alpha _10.3`$, can be due to the passage of the synchrotron cooling frequency, $`\nu _c`$, through the $`R`$-band. This has been inferred, for example, in the afterglow of GRB 971214, at $`t0.2`$ days (Wijers & Galama 1999). If the steepening is due to $`\nu _c`$, this rules out models in which the ejecta expand into a circumburst medium with $`\rho r^2`$ (hereafter, Wind), because in this model $`\nu _c`$ increases with time ($`t^{1/2}`$; Chevalier & Li 1999), and one expects $`\mathrm{\Delta }\alpha =0.25`$.

There are two remaining models to consider in this case: (i) spherical expansion into a circumburst medium with constant density (hereafter, $`\mathrm{ISM}_\mathrm{B}`$; Sari, Piran, & Narayan 1998), and (ii) a jet with $`\theta _{\mathrm{jet}}<\mathrm{\Gamma }_{t0.06\mathrm{d}}^1`$ (i.e. a jet break prior to the first observation at $`t0.06`$ days; hereafter, $`\mathrm{Jet}_\mathrm{B}`$). The subscript B indicates that $`\nu _c`$ is blueward of the optical bands initially. In both models we use Eqn. 1 for the $`R`$-band lightcurve, with $`t_b`$ defined as the time at which $`\nu _c=\nu _R`$, and $`\alpha _2\alpha _11/4`$.

We find that in the $`\mathrm{ISM}_\mathrm{B}`$ model $`t_c0.4`$ days, while in the $`\mathrm{Jet}_\mathrm{B}`$ model $`t_c0.65`$ days. Moreover, in both models the closure relations can only be satisfied by including a contribution from dust extinction within the host galaxy, $`A_V^{\mathrm{host}}`$. We estimate the required extinction using the parametric extinction curves of Cardelli, Clayton, & Mathis (1989) and Fitzpatrick & Massa (1988), along with the interpolation calculated by Reichart (2001). Since the redshift of GRB 020124 is not known we assume $`z=0.3,\mathrm{\hspace{0.17em}1},\mathrm{\hspace{0.17em}3}`$, which spans the range of typical redshifts for the long-duration GRBs. The inferred values of $`A_V^{\mathrm{host}}`$ are summarized in Tab. 4, and range from 0.2 to 0.9 mag.

### 3.2. Jet Break

An alternative explanation for the steepening is a jet expanding into: (i) an ISM medium with $`\nu _c`$ blueward of the optical bands (J-$`\mathrm{ISM}_\mathrm{B}`$), (ii) a Wind medium with $`\nu _c`$ blueward of the optical bands (J-$`\mathrm{Wind}_\mathrm{B}`$), and (iii) an ISM or Wind medium with $`\nu _c`$ redward of the optical bands (J-$`\mathrm{ISM}/\mathrm{Wind}_\mathrm{R}`$). We note that the J-$`\mathrm{ISM}_\mathrm{B}`$ model is different than the $`\mathrm{ISM}_\mathrm{B}`$ model (§3.1) since previously it was defined such that the jet break is later than the last observation. In these models, $`t_bt_{\mathrm{jet}}`$ is the time at which $`\mathrm{\Gamma }(t_{\mathrm{jet}})\theta _{\mathrm{jet}}^1`$.

From the closure values we note that the J-$`\mathrm{ISM}/\mathrm{Wind}_\mathrm{R}`$ requires no extinction within the host galaxy, while the J-$`\mathrm{ISM}_\mathrm{B}`$ and J-$`\mathrm{Wind}_\mathrm{B}`$ models require values of about $`0.05`$ to $`0.3`$ mag.

We find $`t_{\mathrm{jet}}1020`$ days, corresponding to $`\theta _{\mathrm{jet}}10^{}`$. Using the measured fluence (§2.1) we estimate the beaming-corrected $`\gamma `$-ray energy, $`E_\gamma 5\times 10^{50}`$ erg, assuming a circumburst density of 1 cm-3 and $`z=1`$ ($`E_\gamma `$ is a weak function of $`z`$). This value is in good agreement with the distribution of $`E_\gamma `$ for long-duration GRBs (Frail et al. 2001).

## 4. Discussion and Conclusions

Regardless of the specific model for the afterglow emission, the main conclusion of §3 is that the optical afterglow of GRB 020124 suffered little or no dust extinction. Still, this afterglow would have been missed by typical searches undertaken even as early as 12 hours after the GRB event. As shown in Fig. 4, about $`70\%`$ of the searches conducted to date would have failed to detect an optical afterglow like that of GRB 020124.

This is simply because the afterglow of GRB 020124 was faint and exhibited relatively rapid decay. From Fig. 5 we note that GRB 020124 is one of the faintest afterglows detected to date (normalized to $`t=1`$ day), and while it is not an excessively rapid fader, it is in the top $`30\%`$ in this category.

Thus, the afterglow of GRB 020124, along with that of GRB 000630 (Fynbo et al. 2001; Fig. 5), indicates that there is a wide diversity in the brightness and decay rates of optical afterglows. In fact, the brightness distribution spans a factor of about 400, while the decay index varies by more than a factor of three. Coupled with the low dust extinction in the afterglow of GRB 020124, this indicates that some dark bursts may simply be dim, and not dust obscured.

Given this wide diversity in the brightness of optical afterglows, it is important to establish directly that an afterglow is dust obscured. This has only been done in a few cases (§1). Therefore, while statistical analyses (e.g. Reichart & Yost 2001) point to extinction as the underlying reason for some fraction of dark bursts, and may even account for an afterglow like that of GRB 020124, it is clear that observationally the issue of dark bursts is not settled, and the observational biases have not been traced fully.

Since progress in our understanding of dark bursts will benefit from observations, we need consistent, rapid follow-up of a large number of bursts to constrain the underlying distribution, as well as complementary techniques which can directly measure material along the line of sight. This includes X-ray observations which allow us to measure the column density to the burst (Galama & Wijers 2001), and thus infer the type of environment, and potential extinction level. Along the same line, radio observations allow us to infer the synchrotron self-absorption frequency, which is sensitive to the ambient density (e.g. Sari & Esin 2001); the detection of radio emission, as in the case of GRB 020124, implies a density $`n10^2`$ cm-3. Finally, prompt optical observations, as we have carried out in this case, may uncover a larger fraction of the dim optical afterglows, and provide a better constraint on the fraction of truly obscured bursts.

J. S. B. is a Fannie and John Hertz Foundation Fellow. F. A. H. acknowledges support from a Presidential Early Career award. S. R. K. and S. G. D. thank the NSF for support. R. S. is grateful for support from a NASA ATP grant. R. S. and T. J. G. acknowledge support from the Sherman Fairchild Foundation. J. C. W. acknowledges support from NASA grant NAG59302. K. H. is grateful for Ulysses support under JPL contract 958056 and for IPN support under NASA grants FDNAG 5-11451 and NAG 5-17100. Support for Proposal HST-GO-09180.01-A was provided by NASA through a grant from the Space Telescope Science Institute, which is operated by the Association of Universities for Research in Astronomy, Inc., under NASA contract NAS5-26555.


format: text

text:
# System Size Coherence Resonance

## Abstract

We show the existence of a system size coherence resonance effect for an ensemble of globally coupled excitable systems. Namely, we demonstrate numerically that the regularity in the signal emitted by an ensemble of globally coupled Fitzhugh-Nagumo systems, under excitation by independent noise sources, is optimal for a particular value of the number of coupled systems. This resonance is shown through several different dynamical measures: the time correlation function, correlation time and jitter.

Noise induced resonance is a topic that has attracted a lot of attention in the last years. In particular, it has been unambiguously shown that the response of some systems to an external perturbation can be enhanced by the presence of noise (stochastic resonance). A different effect is that of coherence resonance in which an excitable system shows a maximum degree of regularity in the emitted signal in the presence of the right amount of fluctuations (or the related one of stochastic resonance without the need of an external forcing) Coherence resonance has also been studied in dynamical systems close to the onset of a bifurcation, as well as in other bistable and oscillatory systems. It has also been analyzed in different neuronal models such as the FitzHugh–Nagumo and Hodgkin-Huxley models. It has been observed experimentally in electronic circuits, either excitable or chaotic, and in lasers operating in an excitable regime.

In an important recent paper , Pikovsky et al. have shown that when one considers an ensemble of coupled bistable systems subjected to an external periodic forcing (and in the presence of a constant amount of noise), it turns out that an optimal response is obtained for an appropriate value of the number $`N`$ of coupled systems. In other words, that there is a resonance with respect to the number of coupled elements, rather to the usual one that involves the noise level. The authors speculate that this system size resonance might be relevant to neuronal dynamics, in which the neuronal connections or the coupling strengths between neurons can be tuned in order to achieve maximum sensitivity to external signals.

In this paper, we extend the previous result by considering an ensemble of globally coupled excitable systems, each one under the influence of its own noise with a fixed intensity, but without an external forcing. We show that there is a coherence resonance effect as a function of the number $`N`$ of coupled systems. More specifically, we show that the excitable systems pulse on average with a regularity which is optimal for a specific value of $`N`$.

Motivated by the biological applications suggested in  we consider an ensemble of $`i=1,\mathrm{},N`$ coupled FitzHugh–Nagumo systems, each one described by the activator, $`x_i`$, and inhibitor, $`y_i`$, variables. The FitzHugh–Nagumo model provides the simplest representation of firing dynamics and has been widely used as a prototypic model for spiking neurons as well as for cardiac cells. The dynamical FitzHugh–Nagumo equations modified to account for the global coupling are as follows:

$`ϵ\dot{x_i}`$ $`=`$ $`x_i{\displaystyle \frac{1}{3}}x_i^3y_i+{\displaystyle \frac{K}{N}}{\displaystyle \underset{j=1}{\overset{N}{}}}(x_jx_i)`$ (1)
$`\dot{y_i}`$ $`=`$ $`x_i+a+D\xi _i(t)`$ (2)
where independent noises of intensity $`D`$ have been added to the slow variables $`y_i`$ as in ref. . The $`\xi _i(t)`$ are white noises with Gaussian distribution of zero mean and correlations $`\xi _i(t)\xi _j(t^{})=\delta _{ij}\delta (tt^{})`$. The difference in the time scales of $`x_i`$ and $`y_i`$ is measured by $`ϵ`$, a small number. The systems are globally coupled, as indicated by the last term of Eq.(1), where $`K`$ is the coupling strength.

We consider the excitable regime, $`a>1`$. In the absence of coupling, $`K=0`$, each FitzHugh–Nagumo system emits pulses. The pulses are trajectories that exit the basin of attraction of the stable fixed point and are triggered by the influence of the noise term. For zero coupling, coherence resonance exists when the regularity of the time between pulses is optimal for a certain value of the noise intensity $`D`$. This behavior is the same for all the systems although the response of each one is uncorrelated with any other.

Let us now study the collective response of the coupled system and compare it with the individual responses. For the collective response, we introduce the average values of the activator and inhibitor variables as

$$X(t)=\frac{1}{N}\underset{i=1}{\overset{N}{}}x_i(t)Y(t)=\frac{1}{N}\underset{i=1}{\overset{N}{}}y_i(t)$$
(3)
By following the approach by Desai and Zwanzig (see also reference ) it is possible to reach an approximate effective equation for these average values of the form:

$`ϵ\dot{X}=F(X)Y`$ (4)
$`\dot{Y}=X+a+{\displaystyle \frac{D}{\sqrt{N}}}\xi (t)`$ (5)
where $`\xi (t)`$ is a white noise source. Although the exact form of the function $`F(X)`$ and the analysis of the approximations assumed in the derivation will be presented elsewhere, we need only to remark here that in the (exact) equation for $`Y(t)`$ the noise intensity appears rescaled as $`D/\sqrt{N}`$. Therefore, this approximation suggests that the optimal effective noise intensity for the appearance of coherence resonance can be achieved by varying the number of coupled elements $`N`$, as in the case of stochastic resonance for the bistable system considered in . To go beyond this approximation, we numerically integrate the equations of motion (1) and (2).

Figure 1 (left panel) shows the time trace for the variable $`X(t)`$ while 1 (right panel) shows the time trace for the variable $`x_i(t)`$ of one of the elements chosen randomly, for three different values of the number of coupled elements (see the caption of the figure for details of the parameters). Notice that for $`N=160`$ the regularity of the emitted pulses is better than that corresponding to larger or smaller values of $`N`$. This is a clear signature of coherence resonance. Moreover, it can be seen that the regularity in the averaged variable $`X(t)`$ is better than in one of the individual elements, showing that the coupling allows for a smoothness of the trace. It is worth noting that the peaks in the collective variable $`X(t)`$ and $`x_i(t)`$ are very well synchronized in time indicating that the individual systems are pulsing synchronously in time. In Figure 2 (left panel) we plot the time trace for the slow variable $`Y(t)`$, as well as a time trace for a single $`y_i(t)`$ (right panel). At variance with the fast variable $`X`$, it turns out that the averaged $`Y(t)`$ shows a very nice regular behavior for an intermediate number of elements while the individual traces $`y_i(t)`$ do not.

We have computed two indicators commonly used to quantify this effect. First, we have computed the time correlation function $`C_X(t)`$ of the averaged $`X`$ variable, defined as

$$C_X(t)=\frac{\delta X(t^{})\delta X(t+t^{})}{\delta X(t^{})^2}\delta X(t)=X(t)X(t^{})$$
(6)
and similarly for the correlation function $`C_Y(t)`$ for the averaged $`Y`$ variable. Here the averages $``$ are with respect to the time $`t^{}`$, after a small transient has been neglected. Figure 3 shows this correlation function for both the $`X`$ and $`Y`$ variables. It can be seen that the correlations extend further in time for an intermediate value, neither very large nor very small, of the number of coupled systems $`N`$. To obtain a quantitative indicator of this effect, we define the characteristic correlation times $`\tau _X`$ and $`\tau _Y`$ for each variable as

$$\tau _{X,Y}=_0^{\mathrm{}}|C_{X,Y}(t)|𝑑t$$
(7)
In practice, the upper limit of the integral is replaced by a value $`t_{max}`$ such that the correlation function can be considered as decayed to its asymptotic value $`C_{X,Y}=0`$ ($`t_{max}=50`$ for the data shown in figure 3). We have plotted these two correlation times in the left panel of figure 4. Both times reach a maximum at approximately the same value $`N160`$, indicating that, for the set of parameters chosen, the maximum extent of the time correlation occurs for this number of coupled excitable systems.

Another common indicator for the regularity of the emitted pulses can be obtained by the jitter of the time between pulses. A pulse in the $`X(t)`$ variable is defined when $`X(t)`$ exceeds a certain threshold value $`X_0`$ (taken arbitrarily as $`X_0=0.3`$, although other values yield similar results). The jitter $`R_X`$ is defined as the root mean square of the time $`T_X`$ between two consecutive pulses normalized to its mean value:

$$R_X=\frac{\sigma [T_X]}{T_X}$$
(8)
and an equivalent definition for the jitter $`R_Y`$ of the $`Y`$ variable. The smaller the value of $`R_{X,Y}`$, the larger the regularity of the pulses (a value of $`R_{X,Y}=0`$ indicates a perfectly periodic signal). It is shown in the right panel of figure 4 that indeed the jitter in both variables have a well defined minimum at a value of $`N80`$, again showing the existence of the system size resonance. When comparing with the results of the correlation time, it is not uncommon that the two indicators (the correlation time $`\tau `$ and the jitter $`R`$) have their optimal values at different values of the system parameters.

In summary, we have shown that an ensemble of globally coupled FitzHugh–Nagumo excitable systems subjected to independent noises pulse on average with a regularity that is maximum for a given value of the number $`N`$ of coupled systems. An approximate calculation indicates that the collective variable $`Y(t)`$ is subjected to a noise of effective intensity $`D/\sqrt{N}`$. Therefore, even in the presence of a large amount of noise ($`D`$ large), it is possible to couple the right number of systems in order to optimize the periodicity of the emitted pulses. Since the FitzHugh–Nagumo system has been used previously to model some biological systems, we believe that our results can be relevant when analyzing the collective response of such systems in a noisy environment and can help to explain the observed size of some ensembles of excitable cells in living organisms.

Acknowledgments: This work is supported by the Ministerio de Ciencia y Tecnología (Spain) and FEDER, projects BFM2001-0341-C02-01 and BFM2000-1108, and NSF grant DMR9813409. J.D.G. acknowledges IBERDROLA for financial support.


text:
# Spin decay and quantum parallelism

## I Introduction

The interest in electron spin dynamics in semiconductor structures has remarkably increased in the recent years generating the emerging field of spintronics . This key word summarizes the efforts to use the spins of quantum objects rather than or in combination with their charge for information processing, or, even more ambitious, for quantum information processing. Meanwhile several proposals for quantum information processing using (electron or nuclear) spins have been put forward .

In quantum information processing the coherence of quantum bits is crucial. This issue becomes particularly important in solid state systems where the qubits are ususally assumed to be affected by much more and much stronger perturbating influences than in other experimental setups related to quantum information processing in various other fields of physics such as atomic physics, quantum optics, or NMR experiments. An important advantage however solid state systems have is that they offer the possibility scalability once individual qubits and elementary gate operations between them are estabilshed. Such a perspective is usually not given in other quantum computation scenarios.

Motivated by these developments we investigate in this work the time evolution of single spin $`\stackrel{}{S}`$ that is coupled inhomogeneously to a non-interacting environment of other spins. A natural realization of such a system is given by the spin of a single electron in a semiconductor quantum dot interacting with surrounding nuclear spins via hyperfine coupling . Alternatively one can think of an electron bound to a phosphoros atom implantated into a silicon matrix , or of other combinations of impurities and host materials such as Si in a GaAs or Ge matrix. The time evolution of an electron spin under such circumstances is of fundamental interest in its own right and of particular relevance to the quantum computing proposal of Ref. . In fact, very recently a series of studies of electron spin dynamics related to the present one has appeared . Here we build on recent work in Ref. , where the dynamics of an electron spin due to the hyperfine interaction with nuclear spins in a semiconductor quantum dot was investigated. This scource for spin decay can be assumed to be the dominant one in a quantum dot geometry where other mechanisms induced by spin-orbit interaction are believed to be suppressed, although this issue has not been clarified entirely yet .

Our approach here is based on numerical simulations of the full quantum mechanical spin dynamics in sufficiently small systems. These studies complement the earlier work in Ref. , where approximation-free analytical results where obtained for the case of a fully polarized system of nuclei. For a more general initial condition such as an unpolarized nuclear spin system low-order time-dependent perturbation theory was employed which unfortunately suffers from divergent terms in higher order.

In our numerical simulations we observe a decay of the electron spin as measured in terms of the expectation values of its components. We study this phenomenon as a function of the initial nuclear polarization, and the type of the initial state of the nuclear spin system. We compare the time evolutions of initial states where the nuclear system is in a simple tensor product state with the situations where the initial nuclear state is randomly correlated. A major results is that the time evolution depends very significantly on the type of initial state of the spin environment. The time evolution of simple tensor product states can be quite individual, while randomly correlated (and therefore highly entangled) states show very reproducible dynamics that mimics the average over the time evolutions of all possible tensor product states. This observation is an example of quantum parallelism in a decoherence phenomenon.

The spin decay is adccompanied by the generation of quantum correlations between the electron spin and the nuclear spins, illustrating a general concept of quantum information theory where the decoherece of a quantum bit (here the decay of the electron spin) is viewed as the result of the generation of entanglement (i.e. quantum correlations) between the qubit and its environment. We quantify this entanglement using well-established methods and concepts of quantum information theory. By this we also hope that studies of this kind will faciliate fruitful interactions between the communities of solid state physics and quantum infromation.

Finally we compare the results for the full quantum mechanical dynamics with simulations of a classical spin model that arises as the classical limit of the underlying Hamiltonian. This comparision shows that the spin decay observed in the quantum system depends crucially on two properties of the system: (i) the inhomogeneity of the hyperfine contact interaction induced by the spatial variation of the electron wave function, and (ii) the quantum mechanical nature of the dynamics allowing for nontrivial correlations (entanglement) between the electron spin and the nuclear spins.

This paper is organized as follows: In section II we briefly describe the details of our modeling and technical aspects of our numerical simulations. In section III we report on our numerical results. We discuss the role of different initial conditions for the nuclear spin system, and the connection between decoherence and the generation of entanglement observed in our simulations. We close with conclusions in section IV.

## II The model

We consider single spin $`\stackrel{}{S}`$ which is coupled inhomogeneously to a non-interacting environment of other spins,

$$=\stackrel{}{S}\underset{i}{}A_i\stackrel{}{I}_i.$$
(1)
The coupling is inhomogeneous since the constants $`A_i`$ vary among the environment spins $`\stackrel{}{I}_i`$. For an electron spin residing in semiconductor quantum dot the coupling constants $`A_i`$ are given by $`A_i=Av_0|\mathrm{\Psi }(\stackrel{}{r}_i)|^2`$ where $`A`$ is an overall coupling parameter and $`v_0`$ the inverse density of nuclei in the material. $`\mathrm{\Psi }(\stackrel{}{r}_i)`$ is the electron envelope wave function at location $`\stackrel{}{r}_i`$. This factor induces a spatial dependence of the cupliong constants $`A_i`$ which is crucial for the spin dynamics. For simplicity we will consider in the following nuclear spin of length $`I=\frac{1}{2}`$ in a spherical quantum dot In our simulations a given number $`N`$ of nuclear spins is contained in a sphere of radius $`R=(3N/4\pi n_0)^{1/3}`$ where $`n_0=1/v_0`$ is the density of nuclei. The electron wave function is given by

$$|\mathrm{\Psi }(\stackrel{}{r})|^2=\left(\frac{1}{\pi (R/a)^2}\right)^{\frac{3}{2}}e^{\frac{r^2}{(R/a)^2}}$$
(2)
where the parameter $`a`$ describes the confinement of the electron in the dot. In the following we shall use $`a=2`$ and the material parameters of gallium arsenide with $`n_0=45.55\mathrm{nm}^3`$. Therefore a typical quantum dot contains about $`N=10^5`$ nuclei. For the alternative scenario of an electron bound to a $`{}_{}{}^{31}\mathrm{P}`$ in silicon the number $`N`$ of nuclear spins effectively coupled to the electron spin is smaller. The Bohr radius of the hydrogen-like electron state is about 3nm, and with the lattice constant of silicon and the natural abundance of $`{}_{}{}^{29}\mathrm{Si}`$ this leads to values for $`N`$ of about a few hundred.

To mimic their sperical distribution also in systems of smaller size used in our simulations we choose the radial coordinate $`r_i`$ of the $`i`$-th nucleus according to $`r_i=[3(i1/2)/(4\pi n_0)]^{1/3}`$ with $`i`$ ranging from $`1`$ to $`N`$.

The Hamiltonian (1) does not include the direct dipolar interaction between nuclear spins. This interaction is by orders of magnitude weaker than the scale $`A`$ of the hyperfine coupling which is of order $`10^5\mathrm{eV}`$ in GaAs . In this material the characteristic time $`T_{2N}`$ for the nuclear spin decay due to dipolar interaction is of order $`10^4\mathrm{s}`$ while the time scales considered in this work will be at least two orders of magnitude smaller. We also mention a recent interesting numerical study by Dobrovitski et al.  on spin dynamics stressing the role of entropy. There a central spin is coupled inhomogeneously to an essentially non-interacting spin environment, where, differently from the present study, an Ising-like coupling was used. To allow for nontrivial dynamics the authors of Ref.  introduced a magnetic field perpendicular to the $`z`$-direction of the Ising coupling.

The model (1) was specifically studied recently in Ref. , where approximation-free analytical results were obtained for the case of a fully polarized system of nuclei. For a more general initial conditions such as an unpolarized nuclear spin system low-order time-dependent perturbation theory was employed which unfortunately suffers from divergent terms in higher order. In this work we choose a different route and perform finite-size exact diagonalizations from which we obtain the full time evolution.

Since our Hamiltonian conserves the total spin $`\stackrel{}{J}=\stackrel{}{S}+_i\stackrel{}{I}_i`$ it is convenient to work in a subspace of given $`J^z`$ having a dimension of

$$\left(\begin{array}{c}N+1\\ \frac{N+1}{2}J^z\end{array}\right).$$
(3)
To obtain the time evolution of the total spin system with the initial state lying in a subspace with a given value $`J^z`$ we diagonalize the Hamitonian within this subspace and compute the time evolution of certain expectation values from the eigensystem data. For initial states having non-zero overlap in several subspaces with different $`J^z`$ the time evolutions obtained in the different spaces have to be superimposed. The fact that the full eigensystem of the Hamiltonian is required in this procedure is different from most other numerical investigations of spin systems where, for reason of the physical questions being investigated, it is sufficient to concentrate on the ground state and some low-lying excitations. In our case we need the full eigensystem and are therefore restricted to system sizes $`N`$ and values of $`J^z`$ yielding dimensions of not more than a few thousands. However, as we will explain in more detail below, our findings neither depend on the restriction to a certain value of $`J^z`$ nor on the specific coupling constants induced by the quantum dot geometry. In fact, qualitatively the same results are obtained when working in the full Hilbert space, or if the coupling parameters are chosen at random from a uniform distribution. On the other hand it is essential that the coupling is inhomogeneous leading to a time evolution which is for all practical time scales aperiodic and in this sense irreversible. The case of homogeneous coupling with all $`A_i`$ being the same is readily solved analytically and generates periodic dynamics with recurrence time $`T=4\pi N\mathrm{}/A`$.

## III Results for the time evolution of spins

In the time evolutions to be discussed below the initial state will always be a simple direct product of the state of the nuclear or environmental spin system and the electron spin with the latter pointing downwards along the $`z`$-direction. Therefore the nuclear spins and the electron spin are initially uncorrelated. For the nuclear spin system itself we consider two types of initial conditions which give rise to significantly different time evolutions.

### A Product states versus randomly correlated states

We investigate two types of initial states for the nuclear spin system which differ crucially in their corrrelation properties and, as a we shall see below, as a consequence of this also in their time evolution: (i) The nuclear spins are initially in a simple tensor product state. If $`J^z`$ is fixed to a certain value such a state consists of eigenstates of $`I_i^z`$ for each nuclear spin $`i`$. If this restriction is not applied tensor product states consisting of more general spin-coherent states are possible. (ii) The nuclear spin state $`|\chi _N`$ is initially a linear superposition, $`|\chi _N=_T\alpha _T|T`$, where the sum goes over all tensor product states $`|T`$ consisting of eigenstates $`I_i^z`$, $`i\{1,\mathrm{},N\}`$, and is, for fixed $`J^z`$, restricted to the appropriate subspace. The coefficients $`\alpha _T`$ in this entangled pure state are subject to a normalization condition and chosen either at random or coherently (for example, can have the same phases).

As we shall see shortly, a single tensor product state on the one hand and a randomly correlated nuclear state on the other hand generate strikingly different time evolutions for the electron spin. Fig. 1 shows numerical data for time-evolved expectation value $`S^z(t)`$ for an initially randomly correlated system and different degrees of its polarization (characterized by $`J^z`$). In all cases, $`S^z(t)`$ decreases in magnitude. With decreasing polarization the decay becomes more pronounced, and the oscillations accompanying this process get suppressed. Note that it is the decay of the envelope in these graphs but not the fast oscillations itself that signals the decay of the spin. The distance between two neighboring maxima of the oscillations can depend slightly on the initial state and the coupling constants in the Hamiltonian. However, a good estimate for this effective period is usually given by $`T=4\pi \mathrm{}/A`$ since $`A/2`$ is an estimate (neglecting quantum fluctuations) for the width of the spectrum, i.e. the difference between the largest and the smallest eigenvalue of the Hamiltonian.

When the nuclear spin system is initially in a randomly correlated state the time evolution of $`S^z(t)`$ is very reproducible in the sense that it depends only very weakly on the particular representation of the initial random state. This is illustrated in Fig. 2 where the results of different initial random configurations are compared for two different system sizes and degrees of polarization.

This behavior of randomly correlated initial states sharply contrasts with the time evolution of simple tensor product nuclear spin state. The upper two panels of Fig. 3 show the time evolution of the electron spin for two initial tensor product states. In the lower right panel we compare this data with the time evolution of a representative of the randomly correlated initial condition. In the former case the time evolution depends significantly on the concrete initial tensor product state, and the decay of the electron spin occurs typically clearly more slowly than in the case of an initially randomly correlated nuclear spin system.

In the lower left panel of Fig. 3 we show the time evolution of the electron spin averaged over all nuclear tensor product states . Comparing the two lower panels one sees that this data is very close to the time evolution of a randomly correlated state. This observation is also made for other system sizes and degrees of polarization and constitutes an example of quantum parallelism : The time evolution of each initially uncorrelated (and therefore classical-like) nuclear state is present in the evolution of a linear superposition of all such states. In other words, the time evolutions of all uncorrelated classical-like states are perfomed in parallel in the time evolution of the randomly correlated state. An experimental consequence of this observation is that if the electron spin dynamics would be detected on an array of independent quantum dots one could not distinguish whether the nuclear spin system in each dot was initially randomly correlated or in an uncorrelated tensor product state. In other words, the spin dynamics of a randomly correlated pure state of the nuclear system in a single dot cannot be distinguished from a mixed state of an ensemble of dots.

The observation that the time evolution of a randomly correlated state quite closely mimics the average over all tensor product initial conditions relies on the cancellation of off-diagonal terms $`\alpha _T^{}\alpha _T^{},T|\stackrel{}{S}(t)|,T^{}`$, $`TT^{}`$, due to the randomness in the phases of the coefficients $`\alpha _T`$. In this sense our system has a self-averaging property. This can be checked explicitly by reducing this randomness. The left panel of Fig. 4 shows the time evolution of a randomly correlated state where the amplitudes $`\alpha _T`$ are restricted to have a non-negative real and imaginary part. This time evolution turns out to be similarly reproducible as before, i.e. it does not depend on the concrete realization of the initial random state, but it is clearly different from the former case since the cancellation of off-diagonal contributions is inhibited . For comparison we show in the right panel data where the amplitudes in the initial nuclear spin state have a random phase but are restricted to have the same modulus. Here the proper averaging process takes place again. The results described so far were obtained in certain subspaces of $`J^z`$ and for the form of coupling constants $`A_i`$ as induced by the quantum dot geometry. However, our findings do not depend on these choices. We have also performed simulations were the initial state has overlap in the full Hilbert space. For a randomly correlated initial nuclear spin state the only difference is that now also transverse components $`S^x(t)`$, $`S^y(t)`$ of the electron spin evolve. However, these are tiny in magnitude and oscillate around zero. For an initial tensor product states these transverse components can become sizable, and the time evolution again strongly depends on the concrete initial tensor product state. Moreover, as mentioned earlier, the exact form of the coupling constants is also not crucial as long as they are sufficiently inhomogeneous. For instance, we obtain qualitatively the same results if we choose the coupling parameters randomly from a uniform distribution.

We also note that coupling a magnetic field to the electron spin has only a quantitative influence on our results. Here again the time dependence of tensor product initial nuclear state is very individual, while a randomly correlated states gives very reproducible results that mimic closely the average over tensor product states.

### B Decoherence and the generation of entanglement

In circumstances of quantum information processing the decay of a qubit is ususally viewed as some ’decoherence’ process due to the environment attacking the quantum information. As seen above, the spin decay is generically slower if the spin environment is initially in a uncorrelated state. This finding suggests that it is advantageous for protecting quantum information to disentangle the environment that unavoidably interacts with the qubit system.

A ’decoherence’ process of the above kind can be viewed as the generation of entanglement between a qubit and its environment. The system investigated here provides an illustrative example for this statement. The entanglement in the total state $`|\mathrm{\Psi }(t)`$ between the central electron spin and its environment can be measured by the von-Neumann entropy of the partial density matrix where either the electron or the environment has been traced out from the pure-state density matrix $`|\mathrm{\Psi }(t)\mathrm{\Psi }(t)|`$ . Tracing out the nuclear system we have

$$\rho _{el}(t)=\left(\begin{array}{cc}\frac{1}{2}+S^z(t)& S^+(t)\\ S^{}(t)& \frac{1}{2}S^z(t)\end{array}\right).$$
(4)
This matrix has eigenvalues $`\lambda _\pm =1/2\pm |\stackrel{}{S}(t)|`$, and the measure of entanglement reads $`E(|\mathrm{\Psi }(t))=\lambda _+\mathrm{log}\lambda _+\lambda _{}\mathrm{log}\lambda _{}`$. Thus, the formation of expectation values $`|\stackrel{}{S}(t)|1/2`$ (or, in the case of fixed $`J^z`$, just $`|S^z(t)|1/2`$), is a manifestation of the entanglement between the electron spin and the nuclear spin system. The maximum entanglement, $`E=\mathrm{log}2`$, is achieved if the electron spin has decayed completely as measured by the expectation values of its components, $`\stackrel{}{S}(t)=0`$. The generation of quantum entanglement between the electron spin and the nuclear spin system signalled by a reduced value of $`\stackrel{}{S}(t)`$ is a main and crucial difference between the quantum system studied here and its classical ’counterpart’ described by a system of Landau-Lifshitz equations. These equations can be obtained from the Heisenberg equations of motion for the quantum system, $`\stackrel{}{S}/t=i[,\stackrel{}{S}]/\mathrm{}`$, $`\stackrel{}{I}_i/t=i[,\stackrel{}{I}_i]/\mathrm{}`$, by performing expectation values of both hand sides within spin-coherent states and assuming that the expectation values of all operator products factorizes to products of expectation values. This procedure becomes exact in the classical limit . The resulting equations do not contain operators any more but just describe the dynamics of three-component vectors (classical spins) of fixed length. We have performed simulations of such a classical spin system by solving the Landau-Lifshitz equation via the fourth-order Runge-Kutta scheme. As a result, the central classical spin performes an irregular chaotic motion which does not show any similarity to the results for the quantum spin-$`\frac{1}{2}`$ case. In particular all qualitative features of quantum effects such as the generation of entanglement (signalled by a decay of spins as measured by their expectation values) are not present in such a time evolution. Therefore the Landau-Lifshitz equation provides only a rather poor description of the underlying quantum system.

Let us now briefly discuss how the different initial conditions can be prepared experimentally. A tensor product initial state can be produced by applying a magnetic field and having the underlying crystal lattice at a temperature high enough such that spin-lattice relaxation processes to the nuclear spins are efficient. These interaction with the phonon environment will effectively perform projection-type measurements on each spin and force the system to be in a state close to a tensor product of nuclear states pointing in each of the two direction along the field axis. Another possibilty is the use of all-optical NMR techniques as described in Ref. . A randomly correlated nuclear state on the other hand can be achieved by cooling down the lattice to temperatures where phonon processes are suppressed. Then the highly anisotropic and long-ranged dipolar interaction will produce a sufficiently ’disordered’ state with a highly irregular pattern of amplitudes when expressed in the tensor product basis, as we have confirmed by explicit simulations of a system of 8 nuclear spin placed on the edges of a cube. The highly correlated (or entangled) character of these states can be detected by following the individual nuclear spins in terms of their expectation values $`|\stackrel{}{I}_i(t)|`$. This quantity decays from its initial value of $`1/2`$ (in a tensor product initial state) on a time scale determined by the dipolar interaction  to values typically close to zero. According to the entanglement measure $`E`$ discussed above this indicates strong entanglement between each nuclear spin $`\stackrel{}{I}_i`$ and its environment of all other nuclear spins. In both cases the initial state of the full system can be prepared by injecting the electron to the quantum dot from an external lead, or the electron state can be prepared by cooling in a magnetic field and ESR techniques.

We finally consider the nuclear spin correlator $`C(t)=I^z(t)I^z(0)`$, $`\stackrel{}{I}=_i\stackrel{}{I}_i`$, which can be measured directly by local NMR like measurements such as magnetic resonance force microscopy . In a subspace of given $`J^z`$ and the electron spin pointing downwards initially this quantity reads $`C(t)=(J^zS^z(t))(J^z+1/2)`$. A realistic initial state will have its dominant weight in a series of subspaces with neighboring $`J^z`$ centered around some value. Then the time evolution of $`S^z(t)`$ is very similar in these subspaces, and the dynamics of the total nuclear spin can be mapped out by measuring the electron spin and vice versa.

## IV Conclusions

In summary we have studied the dynamics of a single spin coupled inhomogeneously to a spin environment. As a main result the time evolution depends sensitively on the type of initial state of the spin environment. While the time evolution of simple tensor product states can be quite individual, randomly correlated states show very reproducible dynamics that mimics the average over the time evolutions of all possible tensor product states. This observation constitutes an example of quantum parallelism in a decoherence phenomenon. This effects is clearly seen for all finite system sizes studied here and can therefore be expected to be also present in realistic quantum dot systems containing about $`N=10^5`$ nuclei, and also in the thermodynamic limit $`N\mathrm{}`$.

The decay of the single spin in terms of its expectation values is due to the formation of entanglement between this spin and its environment. Since this decay is generally slower if the spin environment is initially in a simple tensor product state (i.e. no entanglement among the environmental spins) our results suggest that it is advantagous for protecting quantum information to disentangle the environment. We expect this result to be of quite general nature, i.e. it should also be valid for other systems consisting of some central quantum object coupled to a bath of other quantum degrees of freedom.

###### Acknowledgements.

This work has been supported by the Swiss NSF, NCCR Nanoscience, DARPA, and ARO.


text:
# Ballistic Magnetoresitance over 4000% in Ni-Ni Electrodeposited Nanocontacts

## Abstract

This paper reports ballistic magnetoresistance values over 4000% measured in electrodeposited Ni-Ni nanocontacts with T geometry previously developed . Over the time, after several magnetic field cycles, the ballistic magnetoresistance relaxed to a 400%. While that the magnetoreistance of a contact could rise indefinitely; relaxtion and reproducibility are, however, the main issue. We find that the tip ending radius conforming the contacts appears not to play the main role.

PACS:72.15.Gd, 75.70.-i

* corresponding author

email: nicolas.garcia@fsp.csic.es

Tel: 34-91-561 88 06

Fax: 34-91-563 15 60

Ballistic magnetoresitance (BMR) of 300% was first reported by Garcia, Muñoz and Zhao, in atomic contacts of resistances larger than 1000 $`\Omega `$  and was explained as resulting from scattering with very sharp magnetic domain walls . Later on it was possible to obtain, by electrodeposition techniques, contacts of a few Ohms resistance with a $`T`$ shape configuration (see Fig.1) having BMR of 700% , that relax to a final ´equilibrium´value of 400% after field cycling the sample during a few days. Contrary to the case of atomic size contacts, the magnetoresistance of these few Ohms contacts (10-30 $`nm`$ size) cannot be attributed  to domain wall scattering . This has been reported by us in Ref. 4. We postulated that the physical cause of BMR in the large nanocontacts was the presence of a dead magnetic layer at the contact point. Thus if the polarization at the Fermi level could be modified to a value near unity, then the resulting BMR could rise up to infinity.

Recently BMR values over 3000% in a Ni-Ni nanocontact with $`T`$ geometry have been reported . The large increase of BMR over the values previously obtained by Garcia $`etal.`$  in similar contacts (3150% against 700%) was attributed to the fact that the tip ending had been sharpened to 40 $`nm`$ by electrochemical etching in KCl  (see Fig.1).

In this communication, we report Ni-Ni contacts made by electrodeposition with the same $`T`$ shape configuration that have 4000% BMR. Both the tip I and the wire II (Fig. 1) are treated electrochemically with KCl, but the tip end was 1000 $`nm`$ in size, and not 40 $`nm`$ size as in the experiment above , which implies that the increase of BMR is not due to the tip size but depends, instead, on the electrochemical processing.

We arranged the Ni wires in a $`T`$ configuration as shown in Fig.1. The applied field during magnetoresistance measurements is in the direction of the Ni wire labeled ´I´in Fig.1. This arrangement, originally conceived by us, is well suited for magnetoresistance across the Ni nanocontact. The tip of the wire ´I´in Fig.1 was positioned to within a few microns to few tens of microns to the Ni wire labeled ´II´prior to electro-deposition of the Ni nanocontact. The Ni wires (except for the region in the immediate vicinity of the contact) were insulated by a fast curing resin epoxy in order to limit the deposition to the gap between the wires. The resin epoxy also firmly holds the Ni wires to the underlying dielectric substrate. The Ni nanocontacts were electrochemically deposited at room temperature. Electro-deposition was performed from a saturated nickel sulphate (Ni2SO4) electrolyte (PH=3.5). We used a cathode potential of -1.2 $`V`$ against a saturated calomel electrode. Deposition times were typically less than one minute. The magnetoresistance measurements were performed at room temperature for magnetic fields up to H=3000 Oe in the current-in-plane / field-in-plane(CIP/FIP) confiuration. The process of fabrication of the tip ´I´was as follows: we first made a tip by breaking mechanically a Ni wire of diameter (0.125 $`mm`$). Then we used an electrochemical etching technique to get the tip radius down to 600-1000 $`nm`$. The Ni wire was deeped into a cell filled with 2M KCl and a 2$`V`$ constant voltage was applied. The tip is shown in Fig.1a. Etching takes place according to the anodic reaction $`Ni(s)+2Cl^{}=NiCl_2+2e^{}`$; the reaction occurring at the Platium cathode is $`2H_2O+2e^{}=H_2+2OH^{}`$. It should be stressed that electrochemical sharpening of the tip was not performed, contrarily to the case in Ref. 6. It would appear that sharpening the tip to 40 $`nm`$ is not the determining step for obtaining contacts with high BMR values, as stated in Ref. 6. We show here that one can get 4000% BMR with a tip sharpened to 1000 nm diameter. Instead, we believe that the additional electrochemical treatment in the KCl cell is the cause for the high BMR.

Fig.2 shows consecutive magnetoresistance loops in a sample whose initial zero-field contact resistance was 15 $`\Omega `$ after electro-deposition. This contact resistance R determines the diameter d=$`\sqrt{1000/R(\Omega )}`$ (in $`nm`$)  of the nanocontact, being equal to 8 $`nm`$ for this sample. As seen in Fig. 2a, the resistance increases rapidly with increasing the field strength. At saturation the resistance rises to 634 $`\Omega `$, after which it remains essentially unchanged with further increases in field strength. This represents a 4000% BMR at room temperature in a field about 800 Oe. The peak of the BMR curve corresponds to about 280 Oe. Fig. 2b shows the BMR curve of the same sample measured after several hours. It can be seen that, although the low field resistance increases somewhat and the BMR value decreases significantly, the high field saturation resistance remains virtually unchanged. This relaxation process needs to be understood and is very important, cannot be neglected, as was already stressed in our previous work  .

In conclusion we have reported BMR over 4000% and have shown that the size of the tip conforming the contact does not play a relevant role because tips 1 micrometer size, as reported here, give also BMR values much larger than the ones previously reported in  . In fact, there is not a single insight in the physics of the BMR process that would point to the importance of sharpening the tip to achieve very high BMR values. The whole thing is regulated by the nanocontacts and on the products O, S, Cl, etc, or density of states at Fermi, that can be segregated at the contact. This has been reported in a very recent paper (cond-mat/0207323)  indicating that the localization, through bonds, of the Ni sp-electrons can lead to full polarization of d-electrons contributing to current and therefore to an indefinitely grow of the BMR.Notice that :

BMR(%)= D$`\times `$2P2/(1-P2)$`\times `$100 (1)

Where D is the dynamical factor approximately equal to unity for very thin domain wall  fixed by the dead magnetic layer. This is the factor that gives the spin conservation in the current or nonadiabaticity. P is the polarization of the electrons contributing to the current if the sp electrons in the case at hand of Ni are bonded by Cl this will not contribute to the current and only the fully polarized d-electrons will participate to the current and thus P~1. Inserting these values in equ. 1 it can be seen that BMR can grow indefinitely (see Ref. 9).

We thank Prof. J. A Rausell for reviewing the manuscript and M. Muñoz for discussions. This work has been supported by the Spanish DGICyT.

$`References:`$

$`FigureCaptions:`$

Fig.1 (a) Optical microscope view of the tip conforming the T contact in 1(b). The tip and the wire have been given an electrochemical treatment in KCl (see text). This seems to enhance BMR to very large values, irrespective of tip ending dimensions.

Fig.2 (a) BMR cycles showing values up to 4014%. However, a new cycle was done after a time and the BMR was reduced to 354% (see 2(b)). This happens always to our samples and is the important point to understand .


text:
# Exact Results for the Reactivity of a Single-File System

## Abstract

We derive analytical expressions for the reactivity of a Single-File System with fast diffusion and adsorption and desorption at one end. If the conversion reaction is fast, then the reactivity depends only very weakly on the system size, and the conversion is about 100%. If the reaction is slow, then the reactivity becomes proportional to the system size, the loading, and the reaction rate constant. If the system size increases the reactivity goes to the geometric mean of the reaction rate constant and the rate of adsorption and desorption. For large systems the number of nonconverted particles decreases exponentially with distance from the adsorption/desorption end.

Zeolitic structures like mordenite, ZSM-22, and $`\mathrm{AlPO}_4`$-o5 are industrially immensely important. They are examples of single-file systems (SFS).Meier et al. (1996) They have one-dimensional pores through which molecules can diffuse, but with cross sections that are too small to allow passing. As for many other 1D systems, this leads to very interesting kinetic effects.Privman (1997) Research on SFS zeolites has focussed on the diffusion, because the mean-square displacement in a SFS is proportional to the square root of time, and not proportional to time.Fedders (1978); van Beijeren et al. (1983); Hahn and Kärger (1996) Molecular Dynamics,Keffer et al. (1995); Hahn and Kärger (1996); Schuring et al. (2000) Dynamic Monte Carlo (DMC),Rödenbeck et al. (1995); Nelson and Auerbach (1999); Nedea et al. (2002) and reaction-diffusion equationsKrishna (2000) have mainly been used. Few studies have included reactions.Tsikoyiannis and Wei (1991); Kärger et al. (1992); Nedea et al. (2002) We have argued that this is unfortunate, because if diffusion is fast the temporal dependence of the mean-square displacement is not very relevant, but there are still effects from the non-passing of the particles in the pores.

In a previous paper we have investigated steady-state properties of a SFS for different assumptions of the reactive site distribution using DMC as well as analytical techniques.Nedea et al. (2002) In the current paper we derive exact results for the case that the entire pore is reactive. These results give detailed insight in the relationships of the various system parameters and the reactivity.

Our model consists of $`S`$ sites forming a one-dimensional finite chain. Each site is numbered consecutively from 1 for the site on one end, to $`S`$ for the site on the other end. Each site is vacant or is occupied by a particle. We have two types of particles: A and B. An A can be converted into a B on any site. Adsorption and desorption of particles can only occur on site 1. Only A’s are adsorbed, but both A’s and B’s desorb. Both types of particle diffuse by making random hops to neighboring sites if vacant.

The evolution of the system is described by a master equationGelten et al. (1999); van Kampen (1981)

$$\frac{dP_\alpha }{dt}=\underset{\beta }{}\left[W_{\alpha \beta }P_\beta W_{\beta \alpha }P_\alpha \right],$$
(1)
where $`\alpha `$ and $`\beta `$ refer to the configuration of the adlayer (a particular distribution of particle over the sites), the $`P`$’s are the probabilities of the configurations, $`t`$ is real time, and the $`W`$’s are constants that give the rates with which reactions change the occupations of the sites. $`W_{\alpha \beta }`$ corresponds to the reaction that changes $`\beta `$ into $`\alpha `$. The rate constants in our model are $`W_{\mathrm{ads}}`$ for adsorption of an A onto site 1 if vacant, $`W_{\mathrm{des}}`$ for desorption of a particle from site 1, $`W_{\mathrm{rx}}`$ for the conversion of an A into a B on any site, and $`W_{\mathrm{diff}}`$ for a hop of a particle to a vacant neighboring site.

If the diffusion is infinitely fast, we can derive a simpler master equation for the number of particles in the system.

$`{\displaystyle \frac{dP_N}{dt}}`$ $`=`$ $`W_{\mathrm{ads}}\left[1{\displaystyle \frac{N1}{S}}\right]P_{N1}W_{\mathrm{ads}}\left[1{\displaystyle \frac{N}{S}}\right]P_N`$ (2)
$`+`$ $`W_{\mathrm{des}}{\displaystyle \frac{N+1}{S}}P_{N+1}W_{\mathrm{des}}{\displaystyle \frac{N}{S}}P_N`$
where $`P_N`$ is the probability that there are $`N`$ particles in the system. This is a master equation of a one-step Markov process.van Kampen (1981)

We are interested in the probability distribution $`f_{MK}(t)`$ that if at time $`t=0`$ the number of particles is $`K`$, this number becomes $`M`$ for the first time at time $`t`$ with $`M<K`$. If there are $`N`$ particles in the system and at time $`t=0`$ a particle adsorbs, then $`f_{N,N+1}`$ is the probability distribution for the time that this particle desorbs. This is based on the property of a SFS that particles can not pass each other. So $`f_{N,N+1}`$ is the probability distribution for the residence time of a particle that adsorbs in a system with $`N`$ particles.

Let $`P_{NK}`$ be the solution of Eq. (2) with $`N`$ the number of particles and the initial condition $`P_{NK}(0)=\delta _{NK}`$. Let $`Q_{NK}`$ also be a solution with $`Q_{NK}(0)=\delta _{NK}`$, but now for the master equation with an adsorbing boundary at $`M`$: i.e., we remove the term in Eq. (2) that corresponds to an adsorption then there are $`M`$ particles in the system. With $`NM`$ we have

$$P_{NK}(t)=Q_{NK}(t)+_0^t𝑑t^{}P_{NM}(tt^{})f_{MK}(t^{}).$$
(3)
This equation is called the renewal equation.van Kampen (1981) If we take $`N=M`$ in the renewal equation, then we have $`Q_{MK}(t)=0`$ by definition. So we get the integral equation

$$P_{MK}(t)=_0^t𝑑t^{}P_{MM}(tt^{})f_{MK}(t^{})$$
(4)
for $`f_{MK}`$.

A particle that adsorbs, stays in the system for a period $`t`$, and then desorbs, has a probability $`\mathrm{exp}(W_{\mathrm{rx}}t)`$ that it is not converted during the time it was in the system. Because the probability distribution for the time that it stayed in the system is given by $`f_{N,N+1}(t)`$, the probability that the particle desorbs without being converted is

$$_0^{\mathrm{}}𝑑t\mathrm{exp}(W_{\mathrm{rx}}t)f_{N,N+1}(t).$$
(5)
This is equal to $`\widehat{f}_{N,N+1}(W_{\mathrm{rx}})`$ with $`\widehat{f}_{N,N+1}`$ the Laplace transform of $`f_{N,N+1}`$. The reason for using this Laplace transform is that is that it is related to the Laplace transform of the solution of the master equation through the renewal equation. Laplace transforming Eq. (4) yields

$$\widehat{P}_{MK}(s)=\widehat{P}_{MM}(s)\widehat{f}_{MK}(s).$$
(6)
If we write the master equation (2) in matrix-vector notation as $`\dot{P}=𝐖P`$, and Laplace transform it we get

$$\underset{M}{}(s\delta _{NM}𝐖_{NM})\widehat{P}_M(s)=P_N(0)$$
(7)
With the initial condition for $`P_{MK}`$ this yields

$$\widehat{P}_{MK}(s)=[(s𝐖)^1]_{MK},$$
(8)
so that

$$\widehat{f}_{MK}(s)=\frac{[(s𝐖)^1]_{MK}}{[(s𝐖)^1]_{MM}}.$$
(9)
We define the reactivity $`B_{\mathrm{prod}}`$ as the number of particles that is being converted per unit time. To get $`B_{\mathrm{prod}}`$ we have to multiply the probability that a particle is converted by the rate of adsorption and make a weighted average over the number of particles in the system. The probability that a particle is converted equals $`1\widehat{f}_{N,N+1}(W_{\mathrm{rx}})`$, and the rate of adsorption equals $`W_{\mathrm{ads}}(1N/S)`$. At steady state the number of particles in the system is given byNedea et al. (2002)

$$P_N^{(ss)}=\left(\begin{array}{c}S\\ N\end{array}\right)\left[\frac{W_{\mathrm{des}}}{W_{\mathrm{ads}}+W_{\mathrm{des}}}\right]^S\left[\frac{W_{\mathrm{ads}}}{W_{\mathrm{des}}}\right]^N.$$
(10)
Combining this results in

$$B_{\mathrm{prod}}=W_{\mathrm{ads}}\underset{N=0}{\overset{S1}{}}P_N^{(ss)}\left[1\frac{N}{S}\right]\left[1\widehat{f}_{N,N+1}(W_{\mathrm{rx}})\right].$$
(11)
Figure 1 shows this reactivity for two finite system sizes and three loadings. The loading $`\theta `$ is defined as the probability that a site is occupied at steady state. It’s equal toNedea et al. (2002)

$$\theta =\frac{W_{\mathrm{ads}}}{W_{\mathrm{ads}}+W_{\mathrm{des}}}.$$
(12)
The reactivity is compared to the (average) rate of adsorption and desorption $`\lambda `$, which at steady state is given by

$$\lambda =W_{\mathrm{ads}}(1\theta )=W_{\mathrm{des}}\theta .$$
(13)
If the rates of adsorption and desorption are small compared to the reaction rate constant $`W_{\mathrm{rx}}`$, then $`B_{\mathrm{prod}}\lambda `$. This means that almost every particle that enters the system is converted. If the rates of adsorption and desorption are much larger, then

$$B_{\mathrm{prod}}S\theta W_{\mathrm{rx}}$$
(14)
for $`\lambda /W_{\mathrm{rx}}\mathrm{}`$ and $`S`$ not too large. Figure 1 shows the $`\theta `$ dependence for small $`S`$ and the $`S`$ dependence for small $`\theta `$. It does not show this $`S`$ and $`\theta `$ dependence for $`S`$ and $`\theta `$ both large. $`\lambda `$ needs to be much larger than the values in the figure for that. For $`S=100`$ we see that the reactivity may be proportional to $`\theta `$ for $`\theta 0.5`$, but not around $`\theta =0.9`$. For the latter value the reactivity is even smaller than for $`\theta =0.5`$.

In practice only $`W_{\mathrm{ads}}`$ can usually be varied independently, either through the pressure or the concentration of A’s outside the system. The limit $`W_{\mathrm{ads}}\mathrm{}`$ corresponds to $`\theta =1`$ and $`\lambda =W_{\mathrm{des}}`$. So we see that there generally is an upper bound on $`\lambda `$.

Comparing $`S=5`$ and $`S=100`$ in Fig. 1 for some finite $`\lambda `$ shows that $`B_{\mathrm{prod}}`$ varies less with $`\theta `$ for larger $`S`$. We can show that in the limit $`S\mathrm{}`$ the reactivity depends only on $`\lambda `$ and $`W_{\mathrm{rx}}`$. If the system becomes infinitely large, then the fluctuations become small with respect to the number of sites. If the system is then at steady state we can write the master equation as

$$\frac{dP_N}{dt}=\lambda (P_{N+1}+P_{N1}2P_N).$$
(15)
Formally we can let the number of particles $`N`$ run from $`\mathrm{}`$ to $`+\mathrm{}`$. The matrix $`s𝐖`$ can easily be diagonalized. The eigenvalues are $`s+2\lambda (1\mathrm{cos}k)`$ with $`\pi <k\pi `$, and the corresponding eigenvector has components $`\mathrm{exp}[ikN]`$. We can use this to calculate elements of the inverse of the matrix $`s𝐖`$.

$`[(s𝐖)^1]_{NM}`$ $`=`$ $`{\displaystyle \frac{1}{2\pi }}{\displaystyle _\pi ^\pi }𝑑k{\displaystyle \frac{\mathrm{cos}(|NM|k)}{s+2\lambda (1\mathrm{cos}k)}}`$ (16)
$`=`$ $`{\displaystyle \frac{1}{ys}}\left[{\displaystyle \frac{y1}{y+1}}\right]^{|NM|}.`$
with

$$y=\sqrt{1+4\lambda /s}.$$
(17)
This then yields for the probability that a particle adsorbs and desorbs at a later time without being converted the expression

$$\widehat{f}_{N,N+1}(W_{\mathrm{rx}})=\frac{y1}{y+1}.$$
(18)
(Here and in the rest of the paper we set $`s=W_{\mathrm{rx}}`$ in the definition of $`y`$.) For the reactivity $`B_{\mathrm{prod}}`$ we have

$$B_{\mathrm{prod}}=\lambda [1\widehat{f}_{N,N+1}(W_{\mathrm{rx}})]=\frac{1}{2}W_{\mathrm{rx}}(y1).$$
(19)
We see that the reactivity no longer depends on the loading $`\theta `$ if the system becomes large, but only on the reaction rate and on $`\lambda `$. For small rates of adsorption and desorption we again find $`B_{\mathrm{prod}}\lambda `$. If the rates of adsorption and desorption are large, then

$$B_{\mathrm{prod}}\left(\lambda W_{\mathrm{rx}}\right)^{1/2}$$
(20)
for $`\lambda /W_{\mathrm{rx}}\mathrm{}`$. The approach to the limit $`S\mathrm{}`$ becomes very slow when $`\lambda `$ is large as can be seen in Fig. 1.

The procedure for the first-passage problem above can also be used to derive the A and B profiles: i.e., the distribution of the A’s and B’s in the chain. We will first deal with the question what the probability is that the $`n`$th particle, counting from site 1, is an A. Then we will answer the question that site $`m`$ is occupied by an A. The answers to the same questions about B’s follow trivially from the ones for the A’s.

We are only interested in the steady state. In this case we have detailed balance: i.e., the number of transitions per unit time from $`N`$ to $`N+1`$ particles is equal to those of $`N+1`$ to $`N`$. This means that for each sequence in the number of particles

$$N_0\stackrel{0}{}N_1\stackrel{\mathrm{\Delta }t_1}{}N_2\stackrel{\mathrm{\Delta }t_2}{}\mathrm{}\stackrel{\mathrm{\Delta }t_{T1}}{}N_T=N_0,$$
(21)
where after a time lapse $`\mathrm{\Delta }t_i`$ the number of particles changes from $`N_i`$ to $`N_{i+1}`$, there is another sequence

$$N_0=N_T\stackrel{0}{}N_{T1}\stackrel{\mathrm{\Delta }t_{T1}}{}N_{T2}\stackrel{\mathrm{\Delta }t_{T2}}{}\mathrm{}\stackrel{\mathrm{\Delta }t_1}{}N_0.$$
(22)
The second sequence is the time reversed of the first one. Moreover, both sequences are equally likely because of the detailed balance. Consequently, the probability distribution that the number of particles in the system is $`K`$ at time $`t=0`$ and $`M`$ with $`M<K`$ at time $`t`$ for the last time equals $`f_{MK}(t)`$. This means that if there are $`N`$ particles in the system the probability distribution that particle $`n`$ is in the system for a time $`t`$ equals $`f_{Nn,N}(t)`$. The probability that that particle has not been converted is then $`\widehat{f}_{Nn,N}(W_{\mathrm{rx}})`$ following the reasoning after Eq. (5). This probability can be calculated using Eq. (9).

For the profile we need the probability that site $`m`$ is occupied by particle $`n`$. This probability is given by

$$P_{\mathrm{occ}}(n,N;m,S)=\left(\begin{array}{c}m1\\ n1\end{array}\right)\left(\begin{array}{c}Sm\\ Nn\end{array}\right)/\left(\begin{array}{c}S\\ N\end{array}\right).$$
(23)
The probability $`\mathrm{A}_m`$ that site $`m`$ is occupied by a particle that has not been converted is then given by

$`\mathrm{A}_m`$ $`=`$ $`{\displaystyle \underset{N=0}{\overset{S}{}}}P_N^{(ss)}\left(\begin{array}{c}S\\ N\end{array}\right)^1`$
$`\times `$ $`{\displaystyle \underset{n=1}{\overset{m}{}}}\left(\begin{array}{c}m1\\ n1\end{array}\right)\left(\begin{array}{c}Sm\\ Nn\end{array}\right)\widehat{f}_{Nn,N}(W_{\mathrm{rx}})`$
where the first summation averages over the number of particles in the system and $`P_N^{(ss)}`$ is given by Eq. (10)

The expression above can be simplified and interpreted more readily for an infinite system. With Eq. (16) we have

$$\widehat{f}_{Nn,N}=\left(\frac{y1}{y+1}\right)^n$$
(25)
for the probability that particle $`n`$ has not been converted. We see that this probability decreases exponentially.

The probability $`P_{\mathrm{occ}}(n,N;m,S)`$ becomes $`P_{\mathrm{occ}}(n,m;\theta )`$ with $`\theta =N/S`$ for $`S\mathrm{}`$. This limit of the combinatorial factors yields

$$P_{\mathrm{occ}}(n,m;\theta )=\left(\begin{array}{c}m1\\ n1\end{array}\right)\theta ^n(1\theta )^{mn}.$$
(26)
Substituting this expression and (25) in (Exact Results for the Reactivity of a Single-File System) yields

$$\mathrm{A}_m=\frac{\theta (xs)(x+s2\theta s)^{m1}}{(x+s)^m}.$$
(27)
Again we find an exponential decrease but now with a characteristic length of $`\mathrm{\Delta }=1/[\mathrm{ln}(y+1)\mathrm{ln}(y+12\theta )]`$. For high loadings ($`\theta 1`$) we find $`\mathrm{\Delta }y/2`$. For low loading ($`\theta 0`$) we get $`\mathrm{\Delta }(y+1)/(2\theta )`$. So the characteristic length is larger for low loadings and higher rate of adsorption and desorption.

Figure 2 shows some typical profiles. The straight line corresponds to the exponential decrease of Eq. (27). The result for finite system sizes can be understood from the fact that smaller systems are less reactive because there are fewer sites at which conversion can take place. As the reactivity must be equal to the number of desorbing converted particles, $`W_{\mathrm{des}}\mathrm{B}_1`$ must be smaller for smaller $`S`$, and because $`\mathrm{B}_1=\theta \mathrm{A}_1`$ the curves in Fig. 2 must start out at larger values of $`\mathrm{A}_1`$. The reactivity also equals $`W_{\mathrm{rx}}_{n=1}^S\mathrm{A}_n`$. This means that $`_{n=1}^S\mathrm{A}_n`$ must be smaller for smaller $`S`$. Because $`\mathrm{A}_1`$ is larger for smaller $`S`$, $`\mathrm{A}_n`$ must decrease faster.

Figure 2 confirms the slow convergence of the system to the limit $`S\mathrm{}`$. At $`S=100`$ the values of $`\mathrm{A}_n`$ drop down to $`210^5`$. At such low value one would not expect an influence of these sites on the kinetics, but one clearly sees differences between $`S=100`$ and $`S\mathrm{}`$ for all sites.

We would finally like to comment on three extensions of our model: adsorption and desorption at both ends, finite diffusion, and other reactions. A system with $`2S`$ sites and adsorption and desorption at both ends has a reactivity that is twice the one of a system with $`S`$ sites and adsorption and desorption at just one end provided that $`S\mathrm{\Delta }`$. DMC simulation show that if $`S`$ is near or smaller than $`\mathrm{\Delta }`$, then the system open at both ends is relatively more reactive.

Finite diffusion is expected to be equivalent to infinite diffusion as long as it is much faster than the other reactions. If this is not the case, then it may become a rate determining step. The system then becomes transport limited. As a consequence the reactivity can only go down. We have seen that in our simulations.Nedea et al. (2002)

We expect that our results will change only little when we change the details of the reaction or when interaction between the particles are included.Sholl and Fichthorn (1997) Our results are a consequence of the non-passing of the particles and of the fact that particles that stay longer in the system have a higher probability of being converted. We therefore also do not expect our results to change when we change our model of discrete sites to some continuous model.


text:
# Frequency modulation and pulse compression by coherent multimode molecular motion

## Abstract

We study beating of a probe field with a time-varying susceptibility in a coherently prepared Raman medium. We consider the general case of an arbitrary variation of susceptibility, which corresponds to a superposition of an arbitrary number of excited Raman transitions. We derive a general analytical solution and conservation relations for this process. We show that the interference between Raman polarizations may substantially affect frequency modulation and pulse compression for the probe field.

We have recently presented a detailed analysis of ultrashort laser pulse compression by parametric beating with a sinusoidal molecular oscillation Kien02 . In this Brief Report, we extend the previous work to the general case of an arbitrarily complex molecular motion. We assume a wave-like molecular excitation, produced by external fields in a dispersionless medium. Molecular motion results in a time-varying susceptibility, which (in the absence of dispersion) leads to time variations in phase and group velocities for a probe field Kien02 . We analyze the effect of these variations on pulse propagation, derive an analytical solution in terms of an integral for a time transformation, and use it to obtain conservation relations.

It was shown recently that coherently driven molecular oscillations can produce frequency modulation with ultrabroad bandwidth Modulation ; all experiments , and result in subfemtosecond some theory ; Sokolov01 ; korn02  and subcycle HarrisSeries  pulse compression. A molecular oscillation can be either sinusoidal (a single coherently excited Raman transition) Modulation ; all experiments ; some theory ; Sokolov01 ; korn02  or more complex (multimode vibrational or rotational molecular wave packets, which correspond to a coherent phased excitation of many Raman transitions) Ivanov02 ; Kapteyn01 . Molecular wave packets can be excited either impulsively korn02 ; Ivanov02 ; Kapteyn01 , or by applying several quasi-monochromatic laser fields HarrisSeries . Generalization of our previous results Kien02  to the important case of a complex molecular motion will provide insight into the Raman pulse compression techniques.

The essence of our result is simple: A pulse is compressed or stretched by time-varying phase and group velocities, such that envelope compression or stretching is accompanied by frequency increase or decrease, respectively. Pulse amplitude and oscillation frequency are changed by reciprocal factors, such that the photon number, the pulse area, and the number of optical oscillations are conserved, while the pulse energy is not. When the resultant change in the pulse bandwidth is small compared to the optical carrier frequency, pulse deformation by the complex non-sinusoidal molecular motion can be obtained as a multiplication of the effects of individual harmonic components of the motion. However, for stronger modulation this intuitive thinking can not be applied. In order to emphasize this point, in the second part of our Report we consider an example of a biharmonic molecular excitation and show that frequency modulation and pulse compression are significantly affected by the interference between two Raman polarizations.

We study the propagation of a probe field $`E`$ in a Raman medium, characterized by a polarization $`P`$. In the local time $`\tau =tz/c`$, the reduced wave equation reads

$$\frac{E}{z}=\frac{1}{2ϵ_0c}\frac{P}{\tau }.$$
(1)
We assume that dispersion is negligible. The instantaneous susceptibility is defined as $`\chi =P/ϵ_0E`$. We assume that $`\chi =\chi _0+\chi _m(\tau z/v)`$, where $`\chi _0`$ is the linear time-independent susceptibility, $`\chi _m`$ is the Raman time-varying susceptibility, and $`v=2c/\chi _0`$ is the average phase and group velocity in the local time coordinates.

We use the reduced local time $`\eta =\tau z/v`$. In the coordinates $`z`$ and $`\eta `$, Eq. (1) becomes

$$\frac{}{z}E(z,\eta )=\frac{1}{2c}\frac{}{\eta }\chi _m(\eta )E(z,\eta ).$$
(2)
The solution to Eq. (2) is found to be

$$E(z,\eta )=E_{\mathrm{in}}(s)G(\eta ),$$
(3)
where $`E_{\mathrm{in}}(s)=E(z=0,s)`$ is the input field, $`s`$ is the input time determined from the output time $`\eta `$ by

$$_s^\eta \frac{d\theta }{\chi _m(\theta )}=\frac{z}{2c},$$
(4)
and $`G(\eta )`$ is the compression factor given by

$$G(\eta )=\frac{\chi _m(s)}{\chi _m(\eta )}.$$
(5)
Differentiating Eq. (4) with respect to $`\eta `$, we obtain

$$\frac{ds}{d\eta }=G(\eta ).$$
(6)
According to Eqs. (3) and (6), the height and duration of optical oscillations in the probe field are changed by the reciprocal factors $`E(z,\eta )/E_{\mathrm{in}}(s)=G(\eta )`$ and $`d\eta /ds=1/G(\eta )`$, respectively. The value $`G(\eta )>1`$ ($`G(\eta )<1`$) indicates pulse compression (stretching) in the vicinity of $`\eta `$. From Eqs. (3) and (6), we find the relation

$$_{\eta _1}^{\eta _2}E(z,\eta )𝑑\eta =_{s_1}^{s_2}E_{\mathrm{in}}(s)𝑑s,$$
(7)
which describes the conservation of the pulse area.

We find from Eq. (6) that the instantaneous oscillation frequency of the output field is

$$\omega _{\mathrm{osc}}(\eta )=G(\eta )\omega _0,$$
(8)
where $`\omega _0`$ is the input frequency. As seen, $`\omega _{\mathrm{osc}}(\eta )`$ is modulated in time by the compression factor $`G(\eta )`$. Using Eqs. (3), (6), and (8), we find

$$\frac{cϵ_0}{2}_{\eta _1}^{\eta _2}\frac{E^2(z,\eta )}{\mathrm{}\omega _{\mathrm{osc}}(\eta )}𝑑\eta =\frac{cϵ_0}{2}_{s_1}^{s_2}\frac{E_{\mathrm{in}}^2(s)}{\mathrm{}\omega _0}𝑑s.$$
(9)
This relation describes the conservation of the photon number, which is always satisfied for Raman processes.

We introduce the mean frequency $`\overline{\omega }(\eta _1,\eta _2)=(\eta _2\eta _1)^1_{\eta _1}^{\eta _2}\omega _{\mathrm{osc}}(\eta )𝑑\eta `$ for the oscillations in the time interval $`(\eta _1,\eta _2)`$. Then, we find from Eqs. (6) and (8) that $`\overline{\omega }(\eta _1,\eta _2)(\eta _2\eta _1)=\omega _0(s_2s_1)`$, that is, the product of the pulse length and the mean frequency is constant during the propagation process. Furthermore, Eq. (3) says that $`E(z,\eta )=0`$ if $`E_{\mathrm{in}}(s)=0`$ and vice versa, that is, a zero of the input field at an input time $`s`$ leads to a zero of the output field at the corresponding output time $`\eta `$. Hence, the number of optical oscillations is conserved.

To get insight into the behavior of the compression factor $`G`$, we derive an explicit, approximate expression for this factor. For this purpose, we consider a particular case where the deviation of an output time $`\eta _i`$ from its corresponding input time $`s_i`$ is small compared to the characteristic variation time $`T_m`$ of $`\chi _m`$, that is, $`|s_i\eta _i|T_m`$. We take $`\eta `$ close to $`\eta _i`$, and use the approximation $`\chi _m(\theta )=\chi _m(\eta _i)+\chi _m^{}(\eta _i)(\theta \eta _i)`$ to calculate the integral in Eq. (4). Then, we find

$$s\eta _i+\frac{\chi _m(\eta _i)}{\chi _m^{}(\eta _i)}=\left[\eta \eta _i+\frac{\chi _m(\eta _i)}{\chi _m^{}(\eta _i)}\right]\mathrm{exp}\left[\frac{z}{2c}\chi _m^{}(\eta _i)\right].$$
(10)
Hence, Eq. (6) yields

$$G(\eta _i)=\frac{ds}{d\eta }|_{\eta =\eta _i}=\mathrm{exp}\left[\frac{z}{2c}\chi _m^{}(\eta _i)\right].$$
(11)
Equation (11) shows that the factor $`G`$ is approximately determined by the time derivative of the Raman susceptibility $`\chi _m`$, multiplied with the propagation length $`z`$. Note that, around the chosen time $`\eta _i`$, the factor $`G`$ is multiplicative, that is, pulse deformation by the complex molecular motion can be obtained as a multiplication of the effects of individual components of the motion. The assumption $`|s_i\eta _i|T_m`$ requires

$$\left|\frac{\mathrm{exp}\left[(z/2c)\chi _m^{}(\eta _i)\right]1}{\chi _m^{}(\eta _i)}\right|\frac{T_m}{|\chi _m(\eta _i)|}.$$
(12)
Condition (12) is satisfied when $`(z/c)|\chi _m^{}(\eta _i)|1`$ and $`(z/c)|\chi _m(\eta _i)|T_m`$, i.e., when the Raman susceptibility and its modulation are small. In this case, we have $`G(\eta _i)1(z/2c)\chi _m^{}(\eta _i)1`$, indicating that the pulse compression is weak. Condition (12) is also satisfied when $`\eta _i`$ is a zero of $`\chi _m`$.

In general, $`G`$ is not multiplicative with respect to the individual susceptibilities. To show this, we consider the case where the Raman susceptibility $`\chi _m`$ is a sum of two sinusoids corresponding to two individual Raman transitions, namely, $`\chi _m=\chi _a+\chi _b`$, where $`\chi _j=\chi _j^{(0)}\mathrm{sin}(\omega _j\eta +\phi _j)`$ for $`j=a,b`$. The compression factor corresponding to the individual component $`\chi _j`$ is given by the periodic comb function Kalosha and Kien ; Kien02

$$G_j=\frac{1}{e^{\alpha _jz}\mathrm{cos}^2\frac{\omega _j\eta +\phi _j}{2}+e^{\alpha _jz}\mathrm{sin}^2\frac{\omega _j\eta +\phi _j}{2}},$$
(13)
where $`\alpha _j=(\omega _j/2c)\chi _j^{(0)}`$. In general, the total compression factor $`G`$ is different from the product of the individual compression factors $`G_a`$ and $`G_b`$. The reason is that the compression factor $`G`$ is produced by the time-varying Raman susceptibility $`\chi _m`$ via a nonlinear mechanism, see Eqs. (4) and (5). Due to the nonlinearity of this mechanism, the individual susceptibility components may interfere with each other in beating with the field.

To see the interference between the individual susceptibility components in beating with the field, we illustrate the functions $`G`$ and $`G_aG_b`$ in Figs. 1(a) and 1(b), respectively. The plots are calculated for the parameters $`\alpha _az=0.8`$, $`\alpha _bz=0.6`$, $`\omega _b=0.07\omega _a`$, and $`\phi _a=\phi _b=0`$. The time is normalized to the Raman period $`T_a=2\pi /\omega _a`$. Both $`G`$ and $`G_aG_b`$ reveal a sequence of teeth corresponding to the oscillation of the susceptibility component $`\chi _a`$, which has a higher frequency. The peaks of these teeth are modulated in accordance with the oscillation of the component $`\chi _b`$, which has a lower frequency. The differences between $`G`$ and $`G_aG_b`$ are clearly observed.

To see in detail the difference between $`G`$ and $`G_aG_b`$, we magnify in Fig. 2 the curves of Fig. 1 for the time windows $`1\eta /T_a2`$ (a) and $`7.5\eta /T_a8.5`$ (b). We also plot in Fig. 2 the output field for the cases where the input peak time is $`\eta _p=0.83T_a`$ (a) and $`8.8T_a`$ (b). The carrier frequency and pulse length of the input probe field are $`\omega _0=15.2\omega _a`$ and $`T=0.08T_a`$, respectively. We observe $`G<1<G_aG_b`$ and $`G>1>G_aG_b`$ around the times $`\eta /T_a=1.5`$ and 8.1, respectively. These relations imply that, due to the interference between $`\chi _a`$ and $`\chi _b`$, the total compression factor $`G`$ and the product of individual factors $`G_aG_b`$ may give opposite indications on pulse stretching and compression. Thus, the interference between the susceptibility components can substantially affect frequency modulation and pulse compression.

In the above, the total compression factor $`G`$ has been compared with the product $`G_aG_b`$ of the individual compression factors. However, $`G_aG_b`$ cannot be rigorously considered as a compression factor for the case where there is no interference between $`\chi _a`$ and $`\chi _b`$. Below we will compare the compression factor of a combined system with the compression factor of a sequence of two individual systems.

We call $`M_aM_b`$ a cascade in which the probe field is sent through, first, a cell $`M_a`$ (with susceptibility $`\chi _a`$ and length $`z_a`$) and, then, a cell $`M_b`$ (with susceptibility $`\chi _b`$ and length $`z_b`$). For this configuration, we have $`E_{\mathrm{out}}(\eta )=E_{\mathrm{in}}(s)G_{ab}(\eta )`$ and $`ds/d\eta =G_{ab}(\eta )`$, where

$$G_{ab}(\eta )=G_a(\eta ^{},z_a)G_b(\eta ,z_b).$$
(14)
Here $`\eta ^{}`$ and $`s`$ are determined by the equations Kalosha and Kien ; Kien02

$`\mathrm{tan}[(\omega _b\eta ^{}+\phi _b)/2]=e^{\alpha _bz_b}\mathrm{tan}[(\omega _b\eta +\phi _b)/2],`$
$`\mathrm{tan}[(\omega _as+\phi _a)/2]=e^{\alpha _az_a}\mathrm{tan}[(\omega _a\eta ^{}+\phi _a)/2].`$ (15)
Note that $`G_{ab}`$ is not symmetric with respect to the indices $`a`$ and $`b`$, i.e., $`G_{ab}G_{ba}`$. In other words, the cascades $`M_aM_b`$ and $`M_bM_a`$ produce different compression factors. This feature is different from the results of Ref. HarrisSeries , which are valid for the case of limited bandwidths.

Comparison between Eqs. (5) and (14) shows that the compression factor $`G`$ of the combined system is different from the factors $`G_{ab}`$ and $`G_{ba}`$ of the cascades. To see the difference between the mixture and cascades of the individual systems $`M_a`$ and $`M_b`$, we plot the factors $`G`$, $`G_{ab}`$, and $`G_{ba}`$ in Figs. 3(a), 3(b), 3(c), respectively, for the time window $`4\eta /T_a5`$. We also plot the output field and the output spectrum for the input peak time $`\eta _p=3.1T_a`$. The other parameters are the same as for Fig. 2. We observe $`G_{ab}<G1<G_{ba}`$ around the time $`\eta /T_a=4.5`$. These relations as well as the temporal and spectral profiles of the output field show clearly the differences between the combined and cascade systems in frequency modulation and pulse compression.

In the case where $`\alpha _jz_j1`$ for $`j=a,b`$, we find to lowest order in $`\alpha _jz_j`$ the compression factor $`G_{ab}=1_j\alpha _jz_j\mathrm{cos}(\omega _j\eta +\phi _j)`$ and the input-output time relation $`s=\eta _j(\alpha _jz_j/\omega _j)\mathrm{sin}(\omega _j\eta +\phi _j)`$. Hence, with the use of the generating function $`e^{i\xi \mathrm{sin}\theta }=_{k=\mathrm{}}^{\mathrm{}}J_k(\xi )e^{ik\theta }`$ of the Bessel functions $`J_k`$, we can expand an input oscillation $`e^{i\omega _0s}`$ into a series of output harmonics $`\mathrm{exp}[i(\omega _0+_jq_j\omega _j)\eta ]`$. Then, we can approximate the spectrum of the field at the output of a series of $`M_a`$ and $`M_b`$ as the product of the Bessel-function spectra of the individual cells HarrisSeries . The condition $`\alpha _jz_j1`$ means that the pulse bandwidths $`\gamma _jz_j\omega _j`$ Kien02 ; some theory , with $`\gamma _j=\alpha _j\omega _0/\omega _j`$, produced by the individual cells, are small compared to the optical carrier frequency $`\omega _0`$. The same results are also obtained for the case of mixed cells (except that individual cell lengths $`z_j`$ should be replaced by a common length $`z`$). Under the condition of small pulse bandwidths, there is no interference between $`\chi _a`$ and $`\chi _b`$, and therefore, no difference between the mixture and cascades of $`M_a`$ and $`M_b`$ in frequency modulation and pulse compression.

Finally, we demonstrate a numerical example for a realistic system, namely, for a cell containing mixed H2 and D2 molecules. We take $`\omega _a=587`$ cm-1 and $`\omega _b=179`$ cm-1 so as to correspond to the rotational transitions of H2 and D2. We use a probe with a carrier frequency $`\omega _0=\mathrm{20\hspace{0.17em}000}`$ cm-1, a pulse length $`T=4.5`$ ps, and a peak time $`\eta _p=0`$. We assume the comb depths $`\alpha _az=0.587`$ and $`\alpha _bz=0.179`$, which correspond to the modulation depths $`\gamma _az=\gamma _bz=20`$. Similar to the results of Ref. HarrisSeries  for the case of a cascade, a broad spectrum with a large number of sidebands, ranging from 8000 cm-1 to 46 000 cm-1, is generated. Unlike the results of Ref. HarrisSeries , the spectrum in Fig. 4(a) is asymmetric. Such asymmetry is due to the deviation of the spectrum from the Bessel-function spectrum, and occurs when the comb depths $`\alpha _jz`$ are not too small compared to unity. The output field in Fig. 4(b) is a long train of slightly compressed sections. As has been shown in Ref. HarrisSeries  for the case of a cascade, single-cycle pulses can be synthesized when the sidebands are phase-corrected by a phase compensator. We illustrate in Fig. 4(c) a single-cycle pulse synthesized from the phase-corrected spectrum of Fig. 4(a).

In summary, we have studied beating of a probe field with a time-varying susceptibility of a Raman medium. We have derived a general analytical solution and conservation relations for this process. We have shown that the interference between Raman polarizations may substantially affect frequency modulation and pulse compression. We emphasize that the analysis of this work can be directly applied only when dispersion is negligible. When dispersion becomes substantial, our analytical solution and conservation relations still provide a useful insight, but a numerical simulation of pulse propagation is necessary in order to obtain an exact solution.

We thank A. K. Patnaik for helpful discussions. A. V. S. acknowledges support from the Welch Foundation.