text:
# PROPER MOTION MEASUREMENTS OF PULSAR B1951+32 IN THE SUPERNOVA REMNANT CTB 80

## 1. Introduction

Associations between pulsars and supernova remnants (SNRs) allow measurements which would not be possible on either population of object alone. Since the center of an SNR marks the presumed pulsar birth-site, the pulsar’s characteristic age, combined with the offset of the pulsar from the center of the SNR, lets us estimate the pulsar’s transverse velocity (Frail, Goss, & Whiteoak (1994)). A more difficult measurement is to actually measure such a pulsar’s proper motion. The direction of motion can confirm (or refute) the association with the SNR, while the magnitude of the motion gives an independent estimate of the pulsar’s age (Gaensler & Frail (2000)).

PSR B1951+32 is a rapidly spinning ($`P=39.5`$ ms) radio, X-ray, and $`\gamma `$-ray pulsar, located on the edge of the unusual SNR CTB 80 (G69.0+2.7) (Strom (1987); Kulkarni et al. (1988)). While B1951+32 is spinning only slightly slower than the Crab pulsar, its low period derivative implies a much larger characteristic age, $`\tau _cP/2\dot{P}=107`$ kyr, and a much lower surface magnetic field, $`B5\times 10^{11}`$ G (Fruchter et al. (1988)).

The distance to the pulsar of $`2.4\pm 0.2`$ kpc, as estimated from its dispersion measure, is consistent with the SNR’s distance of 2 kpc as measured with H i absorption (Strom & Stappers (2000)), suggesting that the two sources are physically associated. However, this system appears quite different from more typical pulsar/SNR associations, in which a young ($`\tau _c20`$ kyr) pulsar sits near the center of an approximately circular shell. This difference can be understood if one invokes an evolutionary scenario in which the pulsar was originally interior to the SNR, but as a result of its high space velocity has caught up with and begun to penetrate its SNR (Fesen, Shull, & Saken (1988); Hester & Kulkarni (1988); Koo et al. (1990)). The pulsar’s relativistic wind then interacts with the SNR shell, re-energizing and distorting it (Shull, Fesen, & Saken (1989)). This sequence can explain both the strange appearance of CTB 80 and the large age for B1951+32 compared to most pulsar/SNR associations.

While this evolutionary picture is aesthetically pleasing, it is imperative that it be verified by observations. Specifically, the model proposed for CTB 80 can be tested by measurement of the proper motion of PSR B1951+32. If the two sources are associated, then the pulsar should be moving away from the center of the SNR, as defined by the shell seen in infra-red and H i (Fesen, Shull, & Saken (1988); Koo et al. (1990)). The separation between the pulsar and the shell’s center, when combined with the pulsar’s characteristic age, lets us predict a proper motion for the pulsar of $`\mu 15`$ mas yr-1 at a position angle (PA) $`250^{}`$ (north through east), independent of the distance to the system.

PSR B1951+32 shows significant timing noise, which prevents measurement of its proper motion through time-of-arrival analysis (Foster et al. (1994)). Thus, only via interferometric measurements can this prediction be tested. In this Letter, we report on measurements of the motion of PSR B1951+32 over an 11-year baseline. In §2 we describe our observations, while in §3 we explain the method used to make our measurements and present our results. In §4 we discuss the implications of these results for the pulsar and SNR.

## 2. Observations

Observations of PSR B1951+32 were made with the Very Large Array (VLA) at epochs 1989.04, 1991.55, 1993.02 and 2000.90. Each observation was eight hours in duration and used the VLA’s A configuration. In the 1989, 1991, and 1993 observations, data were recorded at two simultaneous frequencies, 1385 & 1652 MHz, within the 20-cm band. Each frequency consisted of 15 channels spread across a 23.5 MHz bandwidth with a sampling time of 10 seconds, resulting in a useable field-of-view of radius $`10^{}`$ centered on the pulsar. For the 2000-epoch observations, we also incorporated data from the Pie Town antenna of the VLBA, which doubled the resolution of the array primarily in one dimension. The observations for that epoch were conducted at 1385 & 1516 MHz, and consisted of 15 channels across a 12.5 MHz bandwidth with 5-sec sampling, again resulting in a 10-arcmin field-of-view. The flux density scale of our observations was determined using observations of 3C 286, while the time-varying gains for each antenna were measured using regular observations of PKS B1923+210 or TXS 2013+370.

After standard editing and calibration, we produced images of the pulsar field for each epoch and frequency, using multi-frequency synthesis to mitigate bandwidth smearing, and discarding all baselines shorter than 10 km (corresponding to spatial scales larger than $`4^{\prime \prime }`$). By removing these shorter baselines, we ensured that emission from SNR CTB 80 and from the compact wind-driven nebula surrounding the pulsar were not detected. The only emission seen in our images was the pulsar itself, the “hot-spot” seen immediately adjacent (Strom (1987)), and various other point sources spread throughout the field.

We have identified seven such sources within $`10^{}`$ of the pulsar, as listed in Table 1. After deconvolving each image using the CLEAN algorithm, we applied a gaussian fit to the pulsar and to each of these seven sources to measure their position and extent at each epoch and frequency.

## 3. Measurement of the Pulsar Proper Motion

In order to accurately measure the pulsar proper motion, we adopted six of the field sources as reference sources, reserving the nearest source to the pulsar (source 1 in Table 1) as a check on our measurements. The reference sources are approximately evenly distributed throughout the field. The quality of each source as an astrometric reference was determined by measuring the vector separation between all possible pairs of sources, and then verifying that no source showed trends or significant changes in the magnitude and direction of these vectors. We also verified that in all cases the dimensions of the fitted gaussian matched those of the synthesized beam, confirming that scatter broadening and bandwidth smearing (in the case of background sources) or scintillation (in the case of the pulsar) were not affecting our positional determinations.

The ionosphere can potentially distort the measured positions, producing a frequency-dependent angular displacement $`\mathrm{\Delta }\theta =k\theta \lambda ^2`$, where $`\theta `$ is the angular separation between a source and the field center, $`\lambda `$ is the observing wavelength, and $`k`$ is a constant. By measuring $`\theta `$ for each source and each frequency, we found that $`k`$ had a value consistent with zero, demonstrating that any ionospheric effects were dominated by other uncertanties in our measurements.

McGary et al. (2001) used the VLA to measure pulsar proper motions an order of magnitude smaller than that expected here, and found a variety of other ways in which the positions of reference sources can be distorted, including relativistic effects due to the Earth’s motion, errors introduced by the VLA correlator and additional empirical corrections which had no simple explanation. At the lower precision required here, we have accounted for all these effects simply by measuring the scatter in the position of each reference source between epochs.

We began by determining the vector distances between all possible pairings of the six reference sources. The standard deviations of the components of each vector in Right Ascension and Declination were computed for the four epochs; these were taken as the uncertainty for a particular pairing. This uncertainty was then decomposed between the two sources in that pair, the relative contributions of the two sources to the joint error being weighted inversely by each source’s signal-to-noise ratio. We were thus able to derive an uncertainty in the position of each source for each pairing; these uncertainties were averaged over the five possible pairings to determine the best measurement of the true uncertainty in each coordinate for each source. This analysis was performed separately for observations at 1385 MHz and at 1516/1652 MHz. A mean reference position for each epoch and frequency was then determined by averaging together the positions of the six reference sources, weighting the contribution from each source inversely by its distance from the field center.

The separation between the pulsar’s position and this mean reference position was then measured for each epoch and frequency, combining in quadrature the calculated uncertainties in the pulsar position with those determined for the reference position. The results of these measurements are plotted in Figure 1, and show a clear motion of the pulsar to the south-west. The flux density of the pulsar was higher during the 1993 observations than at other epochs, resulting in smaller uncertainties in the pulsar position for this measurement. Note that the position of the pulsar listed in Table 1 is consistent with previous determinations (Foster et al. (1994)), but with larger errors due to the systematic errors discussed above. Such uncertainties in the absolute astrometry do not affect our proper motion measurement, which has been determined based on astrometry relative to nearby background sources.

As an independent test of our approach, we have similarly measured the proper motion for source 1. As is shown in Figure 1, no change in position is seen for this source, demonstrating that the motion measured for the pulsar is real and that the uncertainties have been realistically assessed.

By applying a weighted least squares fit to the pulsar’s position at each epoch, we find a proper motion at 1385 MHz of $`\mu =31\pm 5`$ mas yr-1, at PA $`=250^{}\pm 7^{}`$, and at 1516/1652 MHz of $`\mu =26\pm 6`$ mas yr-1 at PA $`=240^{}\pm 10^{}`$. These two measurements are consistent with each other within their uncertainties. After taking into account apparent motion of the pulsar at the level of 5 mas yr-1 due to differential Galactic rotation,111To make this calculation, we have used the rotation curve of Fich, Blitz & Stark (1989) and have assumed a 10% uncertainty in the systemic velocity of the pulsar. the combination of our two measurements yields a motion $`\mu =25\pm 4`$ mas yr-1 at PA $`=252^{}\pm 7^{}`$. The corresponding projected velocity is $`V_t=(240\pm 40)d_2`$ km s-1 for a distance $`2d_2`$ kpc.

## 4. Discussion

The transverse velocity we have inferred for PSR B1951+32 is in agreement with the value $`V_t300`$ km s-1 implied by scintillation in its dynamic spectrum (Fruchter et al. (1988)). Furthermore, as shown in Figure 2 the measured position angle of the proper motion agrees closely with the direction predicted if the adjacent “hot-spot” is interpreted as a bow-shock driven ahead of the pulsar (Strom (1987)). It can also be seen in Figure 2 that the pulsar is clearly traveling directly away from the center of its SNR. While this result is not unexpected given the already strong evidence for an association between the pulsar and SNR, it is surprising to realize that this is the first time such an effect has actually been observed. Such measurements for other young pulsars in SNRs have failed to demonstrate such motion, and have either required specific models of SNR and pulsar evolution to maintain the association (Bailes et al. (1989); Gaensler & Frail (2000)), or have caused the association to be abandoned (Thompson & Córdova (1994)).

While the direction in which the pulsar is moving is as predicted, the magnitude of the pulsar’s motion is not in agreement with the simplest expectations. By projecting the pulsar’s position back by $`\tau _c=107`$ kyr to its inferred birth-site, it can be seen from Figure 2 that if the pulsar was born near the center of its SNR, then the pulsar’s true age, $`t_p`$, must be less than $`\tau _c`$.

We can make an initial estimate of $`t_p`$ by simply identifying the offset, $`\mathrm{\Theta }`$, of the pulsar from the SNR’s center. We estimate from the morphology of the infra-red shell shown in Figure 2 that $`\mathrm{\Theta }27^{}`$, so that $`t_p=\mathrm{\Theta }/\mu 64`$ kyr. However, a high velocity progenitor and/or asymmetric SNR expansion can produce a significant offset between the pulsar birth-site and the geometric center of the SNR (Różyczka et al. (1993); Dohm-Palmer & Jones (1996); Hnatyk & Petruk (1999)). Given the complicated evolutionary history of CTB 80, it is perhaps overly optimistic to define a precise geometrical center and simply assume the pulsar was born there.

If the birth-site is not at the SNR’s geometric center, there is no reason to expect the pulsar to be traveling away from this position (Gvaramadze (2000)). However, the farther the birth-site from the center of the SNR, the less likely the pulsar should be moving in this direction by chance. We can therefore quantify the uncertainty in the pulsar’s birthplace as follows. We assume that for a given angular separation $`\varphi `$ from the SNR’s geometric center, it is equally likely that the pulsar was born at any position on a circle concentric with the SNR and of radius $`\varphi `$. Thus, for a particular value of $`\varphi `$, we can compute the fraction of possible birth-sites which are consistent with the pulsar’s direction of proper motion, PA $`=252^{}\pm 7^{}`$. By considering the full range of values $`0\varphi \mathrm{\Theta }`$, we can build up a probability distribution of possible pulsar ages, from the integral of which we can determine the $`1\sigma `$ confidence limits on $`\varphi `$. Using this approach, we find that the pulsar was born within $`\pm 8^{}`$ of the SNR’s center, and that the age implied by our proper motion measurement is $`t_p=64\pm 18`$ kyr. These estimates are independent of the distance to the system and on the site of the supernova explosion; they depend only on the assumptions that an approximate geometric center can be defined for the SNR, and that the pulsar is moving away from its birthplace. The age we have derived is consisent with the estimate of 77 kyr, derived from the expansion velocity of the H i shell coincident with the SNR (Koo et al. (1990)).

It is not surprising that this system’s true age is less than the pulsar’s characteristic age. For a pulsar of constant moment of inertia and magnetic dipole moment, it is straightforward to show that (Manchester & Taylor (1977)):

$$t_p=\frac{2\tau _c}{n1}\left[1\left(\frac{P_0}{P}\right)^{n1}\right],$$
(1)
where $`n`$ is the pulsar’s “braking index” and $`P_0`$ is its initial spin-period. For $`n=3`$ and $`P_0P`$, Equation (1) reduces to $`t_p=\tau _c`$, but in general neither such condition will be satisfied. While $`n`$ has not been measured for PSR B1951+32, we find $`1.5n3`$ for the five pulsars which have had their braking indices measured (Zhang et al. (2001) and references therein). For $`n=3`$, the range in $`t_p`$ determined above implies $`P_0=25\pm 5`$ ms for PSR B1951+32, while for $`n=1.5`$ we find $`P_0=29\pm 3`$ ms. Combining these estimates, we find the initial period most likely falls in the range $`P_0=27\pm 6`$ ms. Again, this estimate does not depend on the assumed distance to the source.

As listed in Table 2, PSR B1951+32 joins a rapidly growing list of pulsars whose initial periods have been determined from their associations with either SNRs or historical supernovae. While these measurements span a reasonable range of initial periods, the data provide good evidence that radio pulsars are born as rapid rotators, with no evidence for a population of longer initial periods as has been sometimes proposed (Vivekanand & Narayan (1981); Chevalier & Emmering (1986)). On the other hand, it is becoming abundantly clear that $`\tau _c`$ can be an unreliable estimator of a young pulsar’s age, with corresponding implications for associations with SNRs, pulsar velocities, and models for neutron star cooling. With many new pulsar/SNR associations now being identified at both radio and X-ray wavelengths, a statistically useful sample of pulsar initial spin periods should soon emerge.

We are grateful to Roger Foster for initiating this project, and to Dale Frail for making us aware of these archival data. We thank Robin McGary and Glen Monnelly for useful advice. The National Radio Astronomy Observatory is a facility of the National Science Foundation, operated under cooperative agreement by Associated Universities, Inc. J. M. M. acknowledges the support of MIT through a Presidential Fellowship for graduate studies in Physics. B. M. G. acknowledges the support of a Clay Fellowship awarded by the Harvard-Smithsonian Center for Astrophysics. E. vdS. is currently supported by the European Commission under the TMR programme, contract number ERB-FMRX-CT98-0168.


format: text

text:
# References

Cherenkov radiation and pair production
by particles traversing laser beams

I.M. Dremin

Lebedev Physical Institute, Moscow 119991, Russia

## Abstract

It is shown that Cherenkov radiation can be observed at TESLA in electron collisions with optical laser pulses. The prospects for it to be observed at SLC, LEP, LHC and RHIC are discussed. The conclusions are compared with results for pair production.

The problem of very high energy charged particles collisions with laser beams is widely discussed now, mostly in connection with the $`e^+e^{}`$-pair production (see latest references ) at SLC and TESLA X-ray laser facilities  and with some other issues of fundamental physics.

Spontaneous particle creation from vacuum induced by a strong external field was theoretically considered in many papers beginning with Refs . However, very powerful high-frequency lasers are needed for this process to be observed. First observations of this effect were done at SLC .

Here, I would like to note that optical lasers can be used for studies of Cherenkov radiation. It is crucial for its observation that the main background process of Compton scattering does not contribute to the kinematical region of Cherenkov radiation. The principal possibility of such a process was first mentioned by V. Ritus in Ref. . The results can be, in general, applied to verify our ideas about the properties of the ”photon medium” in the region where new physics concepts can become essential and to measurements of beam energy and laser bunch parameters.

X-ray lasers have been proposed for use in $`e^+e^{}`$-pair production studies because the quanta energies are high enough to reach the threshold energy which in c.m.s. is equal to $`m+2m_e`$, where $`m_e`$ is the electron mass and $`m`$ is the accelerated particle mass (equal to $`m_e`$ at SLC, LEP, TESLA, to the proton mass at LHC and to the nucleus mass at RHIC). At the same time, for studies of Cherenkov radiation other characteristics of a laser, namely, the ratio $`F/F_0`$ of the electric field $`F`$ to its ”critical” value $`F_0=m_e^2/e`$ or, equivalently, its peak power density $`S`$, are important. They determine the index of refraction $`n`$ of the ”photon medium” in laser bunches. The difference of $`n`$ from 1 is proportional to the density of photons in a laser pulse, i.e., to $`S`$. Just this difference defines the threshold of Cherenkov radiation, its emission angle and intensity111The analogous problem was considered in Ref.  for particles traversing the cosmic microwave background radiation. The density of relic photons is, however, extremely low and, therefore, the index of refraction is so close to 1 that the threshold energy is too high for this effect to be observable.. The parameters $`F/F_0`$ and $`S`$ are higher for optical lasers than for presently available X-ray lasers. That is why they would be preferred for Cherenkov radiation studies nowadays. Besides, the energy limitations also favour optical lasers for this purpose.

The necessary conditions for Cherenkov radiation to be observed are the excess of the index of refraction $`n`$ over 1, i.e.

$$\mathrm{\Delta }n=n1>0$$
(1)
and the real emission angle, given by the formula

$$\mathrm{cos}\theta =\frac{1}{\beta n},$$
(2)
where $`\beta =v/c=\sqrt{1\frac{m^2}{E^2}},m,E`$ are the particle mass and energy. For small values of $`m/E`$ and $`\mathrm{\Delta }n`$ one gets

$$\theta \sqrt{2\mathrm{\Delta }n\frac{m^2}{E^2}}=\sqrt{2\mathrm{\Delta }n\gamma ^2}.$$
(3)
Hence, the condition for the energy to exceed the threshold for Cherenkov radiation $`E_{Ct}`$ is written as

$$\gamma m=EE_{Ct}=\frac{m}{\sqrt{2\mathrm{\Delta }n}}=\gamma _{Ct}m.$$
(4)
It is easily seen that the threshold can become very high for small $`\mathrm{\Delta }n`$.

The formula (3) can be rewritten as

$$0\theta =\frac{\sqrt{\gamma ^2\gamma _{Ct}^2}}{\gamma \gamma _{Ct}}\frac{1}{\gamma _{Ct}}=\theta _{max}(\gamma \mathrm{}).$$
(5)
It is seen that the emission angles of Cherenkov radiation increase from 0 at the threshold to $`\theta _{max}=\gamma _{Ct}^1`$ for $`\gamma \mathrm{}`$. However, already at $`E=2E_{Ct}`$ this angle is very close to $`\theta _{max}`$ ($`\theta (2E_{Ct})0.866\theta _{max}`$).

The number of Cherenkov photons emitted by a single particle with the electric charge $`e`$ in the interval of frequencies $`d\omega `$ from the path length $`dl`$ is given by the common expression

$$\frac{dN_1}{d\omega dl}=2\alpha \mathrm{\Delta }n,$$
(6)
where the fine structure constant $`\alpha 1/137`$. Thus all physical characteristics of the process are determined by the value $`\mathrm{\Delta }n`$. The intensity of the radiation (6) decreases with the threshold energy (4) increase:

$$\frac{dN_1}{d\omega dl}=\frac{\alpha m^2}{E_{Ct}^2}=\frac{\alpha }{\gamma _{Ct}^2}.$$
(7)
The value of $`\mathrm{\Delta }n`$ is uniquely related to the polarization operator of $`\gamma \gamma `$ scattering. For high energy electrons (protons), the laser field can be considered as the constant crossed (or null) field. The refractivity index can be expressed in terms of the photon mass acquired in such a field. It has been calculated in Refs  and its graphical representation can be found in Refs  (for other approaches see also Refs ). According to the results of Refs , the value of $`\mathrm{\Delta }n`$ depends on the photon mass $`\mu `$ and its energy $`\omega `$ in a following way

$$\mathrm{\Delta }n=\mathrm{Re}\mu ^2/2\omega ^2.$$
(8)
In head-on collisions, the photon mass depends only on the invariant variable

$$\kappa =\frac{2\omega }{m_e}\frac{F}{F_0}.$$
(9)
The value of Re$`\mu ^2`$ is negative222We consider the value for transverse polarized photons. For other polarizations it differs by a factor less than 2, and this does not change general conclusions. in the region about

$$0\kappa 15$$
(10)
and has a minimum at $`\kappa 5`$ with $`\mathrm{Re}\mu ^20.2\alpha m_e^2`$. According to eq. (8), the refractivity index exceeds 1 in this region, and, consequently, the Cherenkov radiation is possible at these values of $`\kappa `$. The perturbation theory is still applicable  because $`\alpha \kappa ^{2/3}1`$.

At a fixed laser intensity, i.e., a fixed ratio $`F/F_0`$, the index of refraction does not depend on $`\omega `$ at low energies

$$\mathrm{\Delta }n\frac{14\alpha F^2}{45\pi F_0^2},$$
(11)
because Re$`\mu ^2`$ is proportional to $`\kappa ^2`$ at small $`\kappa ^2`$. Thus, the ratio $`F/F_0`$ defines there main features of Cherenkov radiation.

The Cherenkov threshold $`\gamma _{Ct}`$ is also completely determined by this ratio as seen from formulas (4), (11). It is the same for electrons and protons. Therefore, the threshold energies $`E_{Ct}`$ are approximately 2000 times higher for protons than for electrons. The formulas (3) and (11) show that, in principle, by measuring the angle $`\theta `$ one can get the energy of the particle beam and the strength of the laser field or its peak power density.

The magnitude of $`\mathrm{\Delta }n`$ decreases at higher values of $`\kappa `$ and becomes negative at $`\kappa >15`$ so that Cherenkov radiation is impossible there.

Even though it can again become positive at extreme energies where the hadronic channels are important, this region is completely inaccessible in collisions with laser beams.

Let us remind that the energy threshold for the $`e^+e^{}`$-pair production processes in high energy head-on collisions of a particle of mass $`m`$ with laser quanta is given by

$$E_{th}\frac{m_e(m+m_e)}{\omega _L}.$$
(12)
It depends on the energy of laser quanta $`\omega _L`$ and is much lower for X-ray lasers than for optical lasers. It is approximately 1000 times higher for protons than for electrons. In particular, the threshold for $`\gamma \gamma `$-collisions follows from (12) at $`m=0`$.

The condition for Cherenkov radiation threshold to be below the pair production threshold imposes the restriction on the laser quanta energies:

$$\omega _L<\sqrt{2\mathrm{\Delta }n}m_e(1+\frac{m_e}{m})=\gamma _{Ct}^1m_e(1+\frac{m_e}{m}).$$
(13)
The condition (13) differs for electron and proton beams only by a factor about 2 in the right hand side.

For optical and X-ray lasers, according to , we accept, correspondingly, $`\omega _L=1.2`$ eV (actually, it varies from 0.12 eV for CO2-laser to 2.35 eV for Nd:glass laser) and 3.1 keV, the ratios $`F/F_0=310^4`$ and $`10^5`$ or ,equivalently, the peak power densities $`S=310^{22}`$W/cm2=$`510^{16}`$ eV4 and $`810^{19}`$ W/cm2 (with a possible goal $`710^{29}`$ W/cm2).

At these parameters, the laser field can be treated as a constant crossed (or null) field because its invariants (see )

$$x=\frac{m_e}{\omega _L}\frac{F}{F_0};\chi =2\gamma \frac{F}{F_0}$$
(14)
are large. For optical lasers with $`F/F_0=310^4`$, one gets $`x=125`$, and at TESLA energies $`\gamma =10^6`$, $`\chi =600`$. Thus, the wavelength $`1/\omega _L`$ is much larger than the formation length $`m_e/eF`$.

Electromagnetic processes at these conditions are extremely interesting. New physics concepts can be necessary here because the effective expansion parameter  $`\alpha \chi ^{2/3}`$ exceeds 0.5. In particular, this would indicate that the constant field allows the interaction with field quanta of the arbitrarily low energies. Radiation effects should be reconsidered.

At $`F/F_0=310^4`$, the relations (9), (10) impose the upper limit $`\omega <12`$ GeV. Only this region of comparatively low energies is admissible for Cherenkov quanta in such strong fields.

Using these characteristics, one also concludes that the condition (13) is satisfied for optical lasers with $`F/F_0>310^5`$ ($`S>310^{21}`$ W/cm2) and not valid for the presently available X-ray lasers. To satisfy it for X-ray lasers, one must achieve the peak power density as high as 1027 W/cm2 which is, nevertheless, within the proclaimed goals .

It follows from eqs (6) and (11) that one should deal with most intensive laser fields to get higher intensity of Cherenkov radiation. Thus, in what follows, we discuss only optical lasers briefly referring to X-ray lasers for some estimates.

The numerical value of $`\mathrm{\Delta }n`$ for $`\gamma `$-quanta in the optical laser field with $`F/F_0=310^4`$ is given by

$$\mathrm{\Delta }n=0.6510^{10}.$$
(15)
Therefore, the typical angles and threshold $`\gamma `$-factors for Cherenkov radiation are

$$\theta _{max}=1.1410^5;\gamma _{Ct}=8.810^4.$$
(16)
This implies that the energy threshold for Cherenkov radiation is exceeded at LEP2 and TESLA since it is $`E_{Ct}^{(e)}=45`$ GeV and is close to the upper energy of SLC. Only with further increase of the laser power, it would be possible to study this process at SLC.

The pair production threshold for optical lasers is about 430 GeV. Thus SLC and LEP energies are well below it while TESLA is just close333However, notice the rather wide spread of available wavelengths for optical lasers mentioned above. to the threshold value.

For proton beams the Cherenkov radiation threshold $`E_{Ct}^{(p)}=83`$ TeV is too high even for LHC. If the optical lasers with $`F/F_0>410^3`$ (the peak power density $`S>510^{24}`$ W/cm2) will become accessible, one can hope to observe this effect there as well. This energy is much lower than the threshold for pair production at proton accelerators which is about 400 TeV.

What concerns the X-ray laser facilities, the threshold for pair production (12) is well below the energies accessible at all high energy accelerators.

Now, let us calculate the intensity of the Cherenkov radiation for electron beams and compare it with the main background process of Compton scattering444The final results are valid for any charged particles.. The total number of Cherenkov quanta emitted in the energy interval $`d\omega `$ by a particle which collides with the laser bunch of the coherent spike length $`L`$ is

$$\frac{dN_{Ch}}{d\omega }=2\alpha \mathrm{\Delta }nL=1.110^5L\frac{F^2}{F_0^2}.$$
(17)
This is the energy distribution at low energies as given by eqs (6), (8). It is almost constant at low energies as demonstrated by eq. (17) but should decrease towards the cut-off at $`\kappa 15`$. For a fixed value of the ratio $`F/F_0`$, the energies of emitted Cherenkov quanta are proportional to $`\kappa `$ and limited according to eqs (9), (10). However, already at $`\kappa 4`$ the magnitude of $`\mathrm{\Delta }n`$ and, consequently, the intensity (6) are about twice lower than at $`\kappa =0`$. Therefore, the effective values of $`\omega `$ important in the distribution can be approximately estimated according to eq. (9) as

$$\omega _{eff}<2m_e\frac{F_0}{F}.$$
(18)
For the values of the ratio $`F/F_0`$ adopted above, one gets $`\omega _{eff}^{(o)}<3.5`$ GeV for optical lasers and $`\omega _{eff}^{(X)}<100`$ GeV for X-ray lasers. One can use eq. (9) for an estimate of $`\mathrm{\Delta }n`$ in these energy regions. The threshold and effective energy of Cherenkov quanta decrease while the emission angle and the intensity of radiation increase with increase of laser fields $`F`$.

The absolute intensity can be evaluated according to the formulas (6), (17). For the coherent spike length $`L1`$ mm, the number of quanta per 1 GeV is estimated as

$$\frac{dN_{Ch}}{d\omega }5\mathrm{GeV}^1.$$
(19)
Thus the emitted energy within the effective interval should be of the order of 30 GeV per 1 mm.

To proceed with similar estimates for Compton scattering, we consider first its kinematics. This leads to the following relation between the emission angle $`\theta `$ and energy $`\omega `$ of the scattered quantum in the laboratory system:

$$\mathrm{cos}\theta =\left(1+\frac{\omega _L}{E}\frac{2\omega _L}{\omega }\right)\left(1\frac{\omega _L}{E}\right)^11\frac{2\omega _L}{\omega }.$$
(20)
The precise limits imposed by this relation on the energy of emitted quanta are given by $`\omega _L\omega E`$. In the right hand side of (20) we have taken into account that the particle beam energy is much higher than the photons energies $`E\omega >\omega _L`$. At the angles typical for Cherenkov radiation (16) the energy of the backscattered quantum obtained from eq. (20) is equal to $`\omega 37`$ GeV while Cherenkov radiation is much softer ($`\omega _{Ch}<3.5`$ GeV) due to the cut-off imposed by the behaviour of $`\mathrm{\Delta }n`$. Such ”soft” photons are emitted at larger angles at Compton scattering. Therefore, there is no overlap of the kinematical regions available for Compton and Cherenkov processes. By separating the relatively soft quanta at the angle (16) one would be able to get rid of the background due to Compton processes. These processes contribute to completely different energy ranges and, therefore, can be easily disentangled.

The multiphoton processes would lead to even harder quanta in Compton scattering. In other words, they are effective at larger angles $`\theta >2\sqrt{\omega _L/\omega }`$. Thus, even at high values of $`x`$ usual Compton effect formulas are needed to estimate its intensity at small angles. The non-linear effects of the laser field are not important at Cherenkov angles (16). Let us also note that these energies are much lower than the threshold for pair production (12) equal to 220 GeV for $`\omega _L=1.2`$ eV.

The total energy loss due to Compton scattering is much larger than for Cherenkov radiation (19). It can be as large as 16/63 of the initial electron energy . This is determined by the extremely hard backscattered quanta.

The use of different lasers inevitably leads to the change of corresponding values of $`F/F_0`$ and $`S`$, and, consequently, of the threshold energy $`E_{Ct}`$. Therefore, one should be careful in estimates. The practical feasibility of observation of such an effect at high energies should be considered in close relation to the definite conditions of a particular experiment. For example, one can not use the Nd:glass laser of SLC experiments  with $`F/F_0=2.310^6`$ at TESLA because the Cherenkov threshold energy becomes too high. The optimum choice would be the laser with highest power density (see (17)).

At the end, let us note that, in principle, the heavy ion accelerator RHIC can be also used for pair production studies with X-ray lasers because the threshold energy according to (12) is equal to 165 GeV per nucleon that is available there. The Cherenkov radiation threshold is the same as for proton accelerators if estimated per nucleon. It has been discussed above in connection with LHC and is not reachable at RHIC.

I am grateful to E.L. Feinberg, V.A Maisheev, A.I. Nikishov and V.I. Ritus for discussions and comments.

This work is supported by the RFBR grant 00-02-16101.


text:
# Spin and interaction effects in quantum dots: a Hartree-Fock-Koopmans approach

## Abstract

We use a Hartree-Fock-Koopmans approach to study spin and interaction effects in a diffusive or chaotic quantum dot. In particular, we derive the statistics of the spacings between successive Coulomb-blockade peaks. We include fluctuations of the matrix elements of the two-body screened interaction, surface-charge potential, and confining potential to leading order in the inverse Thouless conductance. The calculated peak-spacing distribution is compared with experimental results.

The traditional description of Coulomb blockade in quantum dots has been the constant-interaction (CI) model, in which the electrons occupy single-particle levels in a confining potential, and the interaction is taken as a constant charging energy. In dots with chaotic dynamics, the fluctuations of the single-particle levels and wave functions can be described by random-matrix theory (RMT). The CI plus RMT model was successful in describing some of the observed statistical properties of such dots, e.g., the conductance peak-height distribution . However, several experiments have demonstrated that other statistics, such as the distribution of the spacing between Coulomb-blockade peaks, are affected by electron-electron interactions .

At low temperatures the peak spacing is given by the second-order difference of the ground-state energy versus the number of electrons. In the CI model and for spin-degenerate levels, the peak-spacing distribution is bimodal, i.e., a superposition of a $`\delta `$ function and a (shifted) Wigner-Dyson distribution. However, none of the measured distributions are bimodal and they all deviate from Wigner-Dyson statistics .

In the spinless case, exact diagonalization for a small number of electrons, Hartree-Fock (HF)  and density functional  calculations, as well as a random interaction matrix model  explained the deviation from Wigner-Dyson statistics as an interaction effect. Spin degrees of freedom were included using exact diagonalization for dots with a small number of electrons and small values of the Thouless conductance $`g`$ . Spin effects were studied in the limit $`g\mathrm{}`$  using the so-called universal Hamiltonian . Interaction effects (in the presence of a magnetic field) were included using a Strutinsky approach for quantum dots . The resulting distributions showed qualitative differences when compared with experiments. Temperature effects were shown to be important at temperatures as low as $`T0.10.2\mathrm{\Delta }`$  ($`\mathrm{\Delta }`$ is the mean spacing between spin-degenerate levels). The spacing statistics were also studied in a spin density functional theory for dots with $`10`$ electrons .

Here we develop a theory that includes spin and interaction effects and is based on a HF-Koopmans approach. The theory is generic and does not require the actual solution of the HF equations. Rather than using the non-interacting basis, we choose as a reference state the $`n`$-electron dot with spin $`S=0`$ ($`n`$ even) and work in its HF basis. We then consider the addition energies and the energy differences between various spin configurations in Koopmans’ limit where the single-particle wave functions do not change . Previously, Koopmans’ approach was discussed for dots with spinless electrons  and was used to study spectral scrambling . Here we derive the generic statistics of the spin and peak spacings, assuming that the HF levels of the reference state satisfy statistics typical of a diffusive or chaotic dot. The theory is valid for large $`g`$, and for $`g\mathrm{}`$ it reduces to the universal Hamiltonian . For finite $`g`$, we include fluctuations of the diagonal matrix elements to leading order in $`1/g`$. Compared with the approach of Ref. , our theory requires the statistics of only a few levels around the Fermi energy, and some of its results are qualitatively different.

The Hamiltonian of the quantum dot in the “disorder” basis $`|i\sigma =a_{i\sigma }^{}|0`$ ($`i`$ denotes a spatial orbit and $`\sigma =\pm 1`$ is the spin variable) is given by

$`H={\displaystyle \underset{i\sigma }{}}ϵ_i^{(0)}a_{i\sigma }^{}a_{i\sigma }+{\displaystyle \frac{1}{2}}{\displaystyle \underset{\genfrac{}{}{0pt}{}{ijkl}{\sigma \sigma ^{}}}{}}v_{ij;kl}a_{i\sigma }^{}a_{j\sigma ^{}}^{}a_{l\sigma ^{}}a_{k\sigma },`$ (1)
where $`ϵ_i^{(0)}`$ are the single-particle energies and $`v_{ij;kl}`$ are the (spin-independent) matrix elements of the Coulomb interaction. The Hamiltonian (1) can be solved in the HF approximation. For each value of the spin projection $`S_z`$, we use Slater determinants with $`n_+`$ ($`n_{}`$) spin up (down) orbitals. Usually, the HF single-particle energies $`ϵ_{\alpha \sigma }^{(n)}`$ and orbitals $`\varphi _{\alpha \sigma }`$ depend on the spin $`\sigma `$. However for even $`n`$, the HF equations have a solution where $`ϵ_\alpha ^{(n)}`$ and $`\varphi _\alpha `$ are independent of $`\sigma `$, and the lowest $`n/2`$ levels are doubly occupied. Such a Slater determinant has good $`S=0`$. We choose this solution as our reference state, and work in its HF basis $`|\alpha `$. This $`S=0`$ state is an eigenstate of the following diagonal many-particle Hamiltonian

$$H_\mathrm{d}=\underset{\alpha \sigma }{}ϵ_{\alpha \alpha }^{(0)}\widehat{n}_{\alpha \sigma }+\frac{1}{2}\underset{\genfrac{}{}{0pt}{}{\alpha \beta }{\sigma }}{}v_{\alpha \beta }^A\widehat{n}_{\alpha \sigma }\widehat{n}_{\beta \sigma }+\underset{\alpha \beta }{}v_{\alpha \beta }\widehat{n}_{\alpha +}\widehat{n}_\beta ,$$
(2)
where $`\widehat{n}_{\alpha \sigma }`$ is the number operator of the state $`\alpha \sigma `$, and $`v_{\alpha \beta }v_{\alpha \beta ;\alpha \beta }`$, $`v_{\alpha \beta }^{\mathrm{ex}}v_{\alpha \beta ;\beta \alpha }`$, and $`v_{\alpha \beta }^Av_{\alpha \beta }v_{\alpha \beta }^{\mathrm{ex}}`$ are diagonal, exchange and antisymmetrized matrix elements, respectively.

The Hamiltonian (2) also has eigenstates with $`S0`$. We label by $`\alpha =0`$ the last occupied level of the $`S=0`$ state. The lowest energy state for each spin $`S`$ is obtained by promoting $`S`$ spin down electrons from $`\alpha =0,\mathrm{},(S1)`$ to spin up electrons in $`\alpha =1,\mathrm{},S`$. The resulting Slater determinant describes a maximal spin projection state $`S_z=S`$ and has good spin $`S`$. Using (2), the energy difference $`\delta E_n(S)E_n(S)E_n(S=0)`$ between the lowest states with spin $`S`$ and spin $`S=0`$ can be written in terms of $`ϵ_\alpha ^{(n)}`$ and a few matrix elements. For example

$`\delta E_n(S=1)`$ $`=`$ $`(ϵ_1^{(n)}ϵ_0^{(n)})v_{10}.`$ (3)
The spin $`S_n`$ of the ground state of the $`n`$-electron dot is determined by minimizing $`\delta E_n(S)`$.

Similarly, the ground-state spin $`S_{n+1}`$ of the dot with $`n+1`$ electrons can be determined from the energy differences $`\delta E_{n+1}(S)E_{n+1}(S)E_{n+1}(S=1/2)`$ for half-integer $`S`$. Assuming Koopmans’ limit, we can express $`\delta E_{n+1}(S)`$ in terms of the HF levels and matrix elements of the reference state $`(n,S=0)`$. For example

$`\delta E_{n+1}(S=3/2)=(ϵ_2^{(n)}ϵ_0^{(n)})v_{10}v_{20}+v_{21}^A.`$ (4)
The addition energy $`\mu _{n+1}E_{\mathrm{gs}}(n+1)E_{\mathrm{gs}}(n)`$ is

$$\mu _{n+1}=\mu _{n+1}(01/2)+\delta E_{n+1}(S_{n+1})\delta E_n(S_n),$$
(5)
where in Koopmans’ limit $`\mu _{n+1}(01/2)=ϵ_1^{(n)}`$.

The spacing $`\mathrm{\Delta }_2`$ between successive peaks is given by the difference in addition energies. We have to distinguish between even-odd-even (“odd”) and odd-even-odd (“even”) transitions (in particle number). We consider here the odd transition $`nn+1n+2`$ ($`n`$ even), for which $`\mathrm{\Delta }_2=\mu _{n+2}\mu _{n+1}`$ (similar results can be derived for the even transition ). $`\mu _{n+1}`$ is given by (5) and $`\mu _{n+2}`$ is calculated from $`\mu _{n+2}=\mu _{n+2}(1/20)+\delta E_{n+2}(S_{n+2})\delta E_{n+1}(S_{n+1})`$ with $`\mu _{n+2}(1/20)=ϵ_1^{(n)}+v_{11}`$. $`\delta E_{n+1}(S)`$ is given by, e.g., Eq. (4), while $`\delta E_{n+2}(S)`$ is calculated from, e.g.,

$`\delta E_{n+2}(S=1)=(ϵ_2^{(n)}ϵ_1^{(n)})v_{11}+v_{21}^A.`$ (6)
To describe the statistics of $`\mathrm{\Delta }_2`$ it is necessary to model the fluctuations of the HF levels and wave functions of the $`(n,S=0)`$ dot. The spectrum is assumed to follow RMT within $`g`$ levels around the Fermi energy, and the matrix elements are uncorrelated from the single-particle spectrum. An exception is the gap $`ϵ_1^{(n)}ϵ_0^{(n)}`$. We find its statistics by comparing the single-particle spectrum of the $`(n+2)`$-electron dot with the spectrum of the $`n`$-electron dot (both in their $`S=0`$ configuration). We have the relation $`ϵ_1^{(n)}ϵ_0^{(n)}=ϵ_1^{(n+2)}ϵ_0^{(n+2)}+v_{01}^A+v_{01}v_{11}`$, in which the levels $`ϵ_0^{(n+2)}`$ and $`ϵ_1^{(n+2)}`$ are both doubly occupied, and thus their spacing should follow Wigner-Dyson statistics. The gap distribution is then a convolution of a Wigner-Dyson distribution with a Gaussian describing the distribution of $`v_{01}^A+v_{01}v_{11}`$ (see below).

We apply our HF-Koopmans approach in a restricted single-particle space of $`g`$ levels around the Fermi energy. The long-range bare Coulomb interaction should then be replaced by an effective interaction. In the limit $`r_s<<1`$ we employ an effective RPA interaction calculated by excluding particle-hole transitions within the above strip of $`g`$ levels . This effective interaction is $`v(𝒓_1,𝒓_2)=e^2/C+v_\kappa (𝒓_1,𝒓_2)+V(𝒓_1)+V(𝒓_2)`$, where $`v_\kappa (𝒓_1,𝒓_2)`$ is a two-body screened interaction in the dot, and $`V(𝒓)`$ is a one-body potential generated by the accumulation of surface charge in the finite dot .

We decompose the interaction in (2) into its average and fluctuating parts, and denote by $`U_d=\overline{v}_{\alpha \beta }`$ and $`J_s=\overline{v}_{\alpha \beta }^{\mathrm{ex}}`$ ($`\alpha \beta `$) the average values of the direct and exchange interaction, respectively. The average interaction is invariant under a change of the single-particle basis when $`\overline{v}_{\alpha \alpha }=U_d+J_s`$ and can be written in terms of the number operator $`\widehat{n}`$ and total spin $`𝑺`$

$`\overline{V}={\displaystyle \frac{1}{2}}\left(U_d{\displaystyle \frac{J_s}{2}}\right)\widehat{n}^2\left({\displaystyle \frac{U_d}{2}}J_s\right)\widehat{n}J_s𝑺^2.`$ (7)
The average interaction is determined by the spatial correlations of the single-particle wave functions . To leading order in $`1/g`$, $`J_s=(\mathrm{\Delta }/2)(1+2\beta ^1b_1/g)`$ , where $`b_1`$ is a geometry-dependent coefficient determined from $`b_1/g=𝑑𝒓\mathrm{\Pi }(𝒓,𝒓)/𝒜`$. Here $`\mathrm{\Pi }(𝒓_1,𝒓_2)`$ is the diffuson propagator and $`𝒜`$ is the area of the dot.

Koopmans’ limit is exact when the fluctuations of the matrix elements are ignored, and our calculations of $`\mathrm{\Delta }_2`$ reduce to those in Ref. , i.e., they can be derived directly from an Hamiltonian that consists of a random one-body part plus an average interaction (7). In the limit $`g\mathrm{}`$ this Hamiltonian is just the universal Hamiltonian with $`J_s=\mathrm{\Delta }/2`$ . A better estimate of $`J_s`$ can be obtained in RPA ; it increases monotonically with $`r_s`$ but remains below $`\mathrm{\Delta }/2`$.

In the simple limit (7), the spin and spacing distributions are determined by a single parameter $`J_s/\mathrm{\Delta }`$ . This limit is demonstrated in Fig. 1 where we show the probabilities of various spin values versus $`J_s`$. The inset shows the standard deviation of $`\mathrm{\Delta }_2`$ versus $`J_s`$ for the even and odd transitions as well as the combined one (“total”).

Next we discuss the fluctuations of the interaction matrix elements , which are approximately Gaussian variables. We discuss separately the bulk screened interaction and surface-charge potential. Using the diagrammatic approach for the two-body screened interaction, we have, for $`r_s1`$ and to leading order in $`1/g`$

$`\mathrm{GOE}:\sigma (v_{\alpha \beta })`$ $`=`$ $`2\sigma _2;\sigma (v_{\alpha \beta }^{\mathrm{ex}})=\sqrt{2}\sigma _2;\sigma (v_{\alpha \alpha })=2\sqrt{2}\sigma _2`$ (8)
$`\mathrm{GUE}:\sigma (v_{\alpha \beta })`$ $`=`$ $`\sigma _2;\sigma (v_{\alpha \beta }^{\mathrm{ex}})=\sigma _2;\sigma (v_{\alpha \alpha })=\sqrt{2}\sigma _2.`$ (9)
Different matrix elements (including the direct $`v_{\alpha \beta }`$ and exchange $`v_{\alpha \beta }^{\mathrm{ex}}`$) are uncorrelated. We note that the coefficients of $`\sigma _1`$ and $`\sigma _2`$ in Eqs. (8) are different from those obtained for a zero-range interaction . The parameter $`\sigma _2`$ is given by

$$\sigma _2=\left[𝒜^2𝑑𝒓_1𝑑𝒓_2\mathrm{\Pi }^2(𝒓_1,𝒓_2)\right]^{1/2}=c_2\mathrm{\Delta }/g,$$
(10)
where $`c_2`$ is a geometry-dependent coefficient. For a disk of radius $`R`$ and boundary conditions of vanishing normal derivative, we find  $`c_2=2[_{l,m}x_{l,m}^4]^{1/2}0.67`$ where $`x_{l,m}`$ are the zeros of $`J_l^{}(x)`$ ($`J_l`$ is the Bessel function of order $`l`$) and $`g\mathrm{\Delta }=2\pi \mathrm{}D/R^2`$ ($`D`$ being the diffusion constant). Since only a few matrix elements contribute to the peak spacing, the contribution of the two-body screened interaction is parametrically of the order $`\mathrm{\Delta }/g`$, unlike the $`\mathrm{\Delta }/\sqrt{g}`$ dependence found in Ref. .

The surface-charge contribution to an interaction matrix element is $`v_{\alpha \beta }=V_\alpha +V_\beta `$, where $`V_\alpha |\psi _\alpha (𝒓)|^2V(𝒓)`$ is a diagonal matrix element of the surface-charge potential. We have

$`\sigma (V_\alpha )=(2/\beta )^{1/2}\sigma _1;\overline{V_\alpha V_\beta }\overline{V_\alpha }\overline{V_\beta }0,`$ (11)
where

$`\sigma _1`$ $`=`$ $`\left[𝒜^2{\displaystyle 𝑑𝒓_1𝑑𝒓_2V(𝒓_1)\mathrm{\Pi }(𝒓_1,𝒓_2)V(𝒓_2)}\right]^{1/2}`$ (12)
$`=`$ $`c_1\mathrm{\Delta }/\sqrt{g}.`$ (13)
For an isolated two-dimensional (2D) circular disk of radius $`R`$, the surface-charge potential can be approximated by $`V(𝒓)=(e^2/2\kappa ϵR)(R^2r^2)^{1/2}`$, where $`\kappa =2\pi e^2\nu /ϵ`$ is the inverse screening length in 2D and $`ϵ`$ is the dielectric constant . We then find $`c_1=2^{1/2}[_{m0}\mathrm{sin}^2x_{0,m}/(x_{0,m}^4J_0^2(x_{0,m}))]^{1/2}0.087`$ .

The above results can be generalized to a ballistic dot, using the ballistic supersymmetric $`\sigma `$ model obtained when a weak disorder with finite correlation length is added . In particular, if the Lyapunov length of this smooth disorder is smaller than the dot’s size, relations similar to Eqs. (8) and (11) can be derived but with the ballistic propagator $`\mathrm{\Pi }_B`$ replacing the diffusive propagator $`\mathrm{\Pi }`$ in Eqs. (10) and (12). For a circular dot we define the ballistic Thouless conductance from the inverse time it takes to cross the diameter $`2R`$ of the dot, leading to $`g=\pi k_FR/2=\pi (n/2)^{1/2}`$ ($`n`$ is the number of electrons in the dot). In the ballistic case

$$\sigma _2=c_2\mathrm{\Delta }[\mathrm{ln}(c_2^{}g)]^{1/2}/g,$$
(14)
where $`c_2=\sqrt{3}/2`$ is a geometry-independent constant. The direct propagation between $`𝒓_1`$ and $`𝒓_2`$ contributes to $`\mathrm{\Pi }_B(𝒓_1,𝒓_2)`$ a term $`1/(\pi k_F|𝒓_1𝒓_2|)`$ which at shorter distances should be replaced by its quantum counterpart $`J_0^2(k_F|𝒓_1𝒓_2|)`$ . The corresponding contribution needs to be renormalized such that its integral over $`𝒓_1`$ (or $`𝒓_2`$) vanishes . The remaining part of $`\mathrm{\Pi }_B`$ involves single or multiple scattering from the boundaries and calculating it requires knowledge of the semiclassical dynamics. It can be calculated analytically for a circular billiard with diffusive boundary scattering . Using this model, we estimate $`c_2^{}=0.81`$ . A similar estimate of (12) gives $`c_1=0.123`$.

We have studied the effects of fluctuations of the interaction matrix elements on the peak-spacing distribution. In Fig. 2 we show the standard deviation $`\sigma (\mathrm{\Delta }_2)`$ versus $`\sigma _1`$ ($`J_s=0.3\mathrm{\Delta }`$ and $`\sigma _2=0.05\mathrm{\Delta }`$) for the orthogonal (solid lines) and unitary (dashed lines) symmetries. $`\sigma (\mathrm{\Delta }_2)`$ increases with $`\sigma _1`$ and it does so faster in the odd case.

Fig. 3 shows typical peak-spacing distributions for both the orthogonal (left) and unitary (right) symmetries for $`\sigma _2=0.025\mathrm{\Delta }`$ and $`\sigma _1=0`$, $`0.03\mathrm{\Delta }`$ and $`0.06\mathrm{\Delta }`$. Signatures of the bimodality can still be observed for $`\sigma _1=0`$ but they disappear at $`\sigma _1=0.03\mathrm{\Delta }`$. Nevertheless, the distributions remain asymmetric (more so in the unitary case).

An additional contribution to the fluctuations of $`\mathrm{\Delta }_2`$ arises from the variation of the gate voltage between peaks. In general, the change of the gate voltage between two peaks leads to a spatially non-uniform change $`U(𝒓)=V(𝒓)+\stackrel{~}{V}(𝒓)`$ in the confining potential, where $`\stackrel{~}{V}`$ originates in the mutual dot-gate capacitance . This leads to scrambling of the HF levels between peaks. For example, let us consider the odd transition. As the gate voltage changes between $`V_g^{n+1}`$ and $`V_g^{n+2}`$, the reference HF levels $`ϵ_\alpha ^{(n)}`$ change by $`\delta ϵ_\alpha ^{(n)}U_\alpha =𝑑𝒓|\psi _\alpha (𝒓)|^2U(𝒓)`$. We therefore substitute $`ϵ_\alpha ^{(n)}ϵ_\alpha ^{(n)}+U_\alpha `$ in the calculation of $`\mu _{n+2}`$ and in Eq. (6) ($`\mu _{n+1}`$ is unchanged and calculated from Eq. (5)). This level scrambling can also lead to spin rearrangement in the dot. The ground-state spin of the $`n+1`$-electron dot at gate voltage $`V_g^{n+2}`$ (just below the transition) is found from (6) (and its generalization to higher values of $`S`$) after the substitution $`ϵ_\alpha ^{(n)}ϵ_\alpha ^{(n)}+U_\alpha `$.

The fluctuation properties of $`U_\alpha `$ are similar to those of $`V_\alpha `$ in Eqs. (11), except that $`V(𝒓)`$ is replaced by $`U(𝒓)`$ in Eq. (12).

We now compare our theory with the experimental results of Ref.  at the lowest measured temperature of $`T=0.22\mathrm{\Delta }`$. At this temperature it is necessary to include the effect of excited states and in particular the contribution from both lowest $`S=0`$ and $`S=1`$ states . We model the gate-voltage scrambling by a potential $`\stackrel{~}{V}`$ whose matrix elements are uncorrelated from the matrix elements of $`V`$ and have the same variance. It is difficult to estimate $`\sigma _1`$ and $`\sigma _2`$ for the ballistic dot used in the experiment. The simple estimates based on a billiard with diffusive surface scattering (see the second paragraph after Eq. (12)) for $`n340`$ electrons (i.e., $`g41`$) give $`\sigma _1=0.02\mathrm{\Delta }`$ and $`\sigma _2=0.04\mathrm{\Delta }`$. Fig. 4 compares the experimental distribution of Ref.  in the presence of a magnetic field (solid histograms) with the corresponding theoretical distribution (dashed histograms) that includes an experimental noise of $`0.1\mathrm{\Delta }`$. The theoretical distribution describes rather well the asymmetry of the experimental distribution and its width $`0.27\mathrm{\Delta }`$ is slightly below the experimental width of $`(0.29\pm 0.03)\mathrm{\Delta }`$. In chaotic billiards our estimates of $`\sigma _1`$ and $`\sigma _2`$ can be enhanced by up to a factor of $`2`$ and lead to better agreement with the data.

In conclusion, we have developed a HF-Koopmans approach to study spin and interaction effects in diffusive or chaotic quantum dots. In particular we have studied the dependence of the peak-spacing distribution on the fluctuations of the interaction matrix elements to leading order in the inverse Thouless conductance. We find good agreement with the lowest temperature data of Ref..

This work was supported in part by the U.S. DOE grant No. DE-FG-0291-ER-40608. We acknowledge useful discussions with H.U. Baranger, Y. Gefen, A. Polkovnikov and G. Usaj, and in particular with A.D. Mirlin.


text:
# Pattern Stability and Trijunction Motion in Eutectic Solidification

\[

## Abstract

We demonstrate by both experiments and phase-field simulations that lamellar eutectic growth can be stable for a wide range of spacings below the point of minimum undercooling at low velocity, contrary to what is predicted by existing stability analyses. This overstabilization can be explained by relaxing Cahn’s assumption that lamellae grow locally normal to the eutectic interface.

\]

The solidification of eutectic alloys is both a striking example of spontaneous pattern formation in nature and a metallurgical problem of widely recognized practical importance . This growth process has been traditionally studied by directional solidification experiments where a sample containing a binary alloy of near-eutectic composition is pulled with a fixed speed $`V`$ in an externally imposed temperature gradient $`G`$. This setup produces a wide range of microstructures of which the simplest is an array of lamellae of two coexisting $`\alpha `$ and $`\beta `$ solid phases growing into the metastable liquid, as shown in Figure 1a. The steady-state growth of a perfectly periodic lamellar array is described by the classic Jackson-Hunt (JH) theory  that predicts the relationship

$$\mathrm{\Delta }T(\lambda )T_ET_{av}(\lambda )=K_1V\lambda +K_2/\lambda ,$$
(1)
between the lamellar spacing $`\lambda `$ (width of one lamella pair) and the undercooling, $`\mathrm{\Delta }T(\lambda )`$, which is the difference between the eutectic temperature, $`T_E`$, at which the three phases ($`\alpha `$, $`\beta `$ and liquid) coexist in equilibrium, and the average temperature of the nonequilibrium eutectic interface spatially averaged over one spacing, $`T_{av}(\lambda )`$. The first and second terms on the right-hand-side of Eq. 1 represent the undercooling necessary to drive the diffusive transport of the two chemical components of the alloy in the liquid (coupled growth) and the capillary undercooling associated with the curvature of the solid-liquid interface, respectively; $`K_1`$ and $`K_2`$ are constants that only depend on the alloy system and the overall composition of the sample.

The JH result implies that the growth undercooling has a minimum value $`\mathrm{\Delta }T_m=(4K_1K_2)^{1/2}V^{1/2}`$ for a spacing

$$\lambda _m=(K_2/K_1)^{1/2}V^{1/2},$$
(2)
which is easily found by setting $`d\mathrm{\Delta }T/d\lambda =0`$.

Lamellar growth is well known to be unstable for spacings smaller than a critical value $`\lambda _c`$. This instability leads to the local elimination of lamellae and is the mechanism by which the array increases its average spacing during the dynamical transient that produces the final pattern. Hence, it is crucially important for understanding pattern selection in this system. Oscillatory instabilities are also known to limit the array stability at large spacing . JH have credited Cahn in Ref.  for pointing out that $`\lambda _c`$ should be equal to $`\lambda _m`$ if one assumes that lamellae grow locally normal to the envelope of the eutectic interface. Langer later formalized this result by showing that a large-scale and small-amplitude spacing modulation of a steady-state array obeys the diffusion equation

$$_t\lambda (x,t)=D_x^2\lambda (x,t),$$
(3)
where $`x`$ is the coordinate perpendicular to the growth axis $`z`$, and $`D=D_{}`$ with

$$D_{}=\frac{V\lambda _0}{G}\frac{d\mathrm{\Delta }T(\lambda )}{d\lambda }|_{\lambda =\lambda _0}=\frac{K_1\lambda _0V^2}{G}\left(1\frac{1}{\mathrm{\Lambda }_0^2}\right).$$
(4)
We have defined $`\mathrm{\Lambda }_0=\lambda _0/\lambda _m`$ where $`\lambda _0`$ is the spacing of the steady-state array being perturbed, and we have used the subscript “$``$” to stress that this expression for $`D`$ is obtained under Cahn’s assumption that lamellae grow normal to the interface. Langer’s analysis reproduces the JH-Cahn result that growth is unstable below $`\lambda _m`$ since perturbations are amplified (decay) when $`D_{}<0`$ ($`>0`$). In addition, it shows that lamella elimination is initiated by a long-wavelength diffusive instability that is generically present in one-dimensional pattern forming systems with translation symmetry along $`x`$.

In this letter, we study the steady-state and stability properties of lamellar eutectic growth by thin-sample directional solidification experiments in the transparent organic system CBr4-C2Cl6 and by two-dimensional simulations of a phase-field model, and we extract from both approaches independent accurate determinations of $`\lambda _m`$ and $`\lambda _c`$. An important and novel component of our experiments is the direct measurement of $`\mathrm{\Delta }T(\lambda )`$, which allows us to obtain $`\lambda _m`$ from the minimum of this curve rather than computing its value from Eq. 2, thereby circumventing uncertainties in materials parameters. We find that, in both experiments and simulations, $`\lambda _c`$ is substantially smaller than $`\lambda _m`$, even for typical directional solidification growth conditions where the two spacings have previously been assumed equal. Furthermore, by analyzing the decay of long-wavelength perturbations of the array in both experiments and simulations, we obtain a direct measurement of $`D`$, which allows us to shed light on the origin of the discrepancy between $`\lambda _c`$ and $`\lambda _m`$.

The experiments were made with a nearly eutectic CBr4-C2Cl6 alloy prepared with zone refined materials in thin ($`12\mu m`$ thick) glass wall samples ($`8mm`$ wide and $`60mm`$ long). The values of $`G`$ used ranged from $`40Kcm^1`$ to $`110Kcm^1`$ ($`\pm 10\%`$), and those of $`V`$ from 0.125 to 0.75 $`\mu ms^1`$($`\pm 4\%`$). Details concerning the preparation of the samples, the solidification setup, and the visualization of the front shape can be found in Refs. .

The steady-state $`\mathrm{\Delta }T(\lambda )`$ curve has never been measured directly due to the fact that $`\mathrm{\Delta }T_m`$ is usually of the order $`0.01K`$, whereas the absolute temperature is not known with a precision better than about $`0.1K`$. To overcome this difficulty, we exploit two key ingredients. Firstly, as will be described elsewhere, we are able to create a large-scale modulation of spacing where $`\lambda (x)`$ varies between two extreme values that comprise $`\lambda _m`$, as shown in Fig. 1(b). Secondly, we measure the $`z`$ coordinate of the solid-liquid interface averaged over one $`\lambda `$, which we denote by $`\overline{z}(x)`$, and compute the local front undercooling using $`T_{av}(x)=G\overline{z}(x)+T_0`$, where $`T_0`$ is an unknown constant. By eliminating $`x`$ between $`T_0T_{av}(x)`$ and $`\lambda (x)`$, we obtain $`T_0T_{av}(\lambda )`$ which we then fit to Eq. 1 expressed in the form $`T_0T_{av}(\lambda )=\mathrm{\Delta }T_m(\lambda /\lambda _m+\lambda _m/\lambda )/2\mathrm{\Delta }T_0`$, using $`\lambda _m`$, $`\mathrm{\Delta }T_m`$ and $`\mathrm{\Delta }T_0=T_ET_0`$ as adjustable parameters. A plot of $`T_0T_{av}(\lambda )`$ and its fit is shown in Fig. 2. The fit is very good for $`\lambda `$ smaller than about $`1.25\lambda _m`$. The departure observed beyond this limit is compatible with the one which exists between the numerically calculated curves $`\mathrm{\Delta }T(\lambda )`$ and the JH approximation . We performed such measurements for $`V`$ ranging from 0.125 to 0.5 $`\mu ms^1`$. We found $`\lambda _m^2V=K_2/K_1=193\pm 16\mu m^3s^1`$ and $`\mathrm{\Delta }T_m^2/V=4K_1K_2=(2.7\pm 1.3)\times 10^3K^2s\mu m^1`$. These values compare well to those calculated from the material constants of CBr4-C2Cl6 given in ref. , namely, $`\lambda _m^2V=185\pm 20\mu m^3s^1`$ and $`\mathrm{\Delta }T_m^2/V=(1.2\pm 0.2)\times 10^3K^2s\mu m^1`$.

To study experimentally the small spacing stability limit, we exploit the fact that $`\lambda _mV^{1/2}`$. Therefore, we can effectively vary $`\lambda _0/\lambda _m`$ by performing downward velocity jumps of relatively large amplitude. Namely, we start from a stable quasi-stationary periodic array of spacing $`\lambda _0`$ at a higher velocity, and then observe whether the same array at the lower velocity, and hence smaller $`\lambda _0/\lambda _m`$, remains stable or becomes unstable.

To measure experimentally the array diffusion constant $`D`$, we use the fact that the amplitude of a long-wavelength modulation of spacing of the form

$$\lambda (x,t)\lambda _0+\delta \lambda _0\mathrm{exp}(ikx+\omega _kt)$$
(5)
decays exponentially in time when the array is stable ($`\lambda _0>\lambda _c`$). Substituting this form in Eq. 3, we obtain the simple dispersion relation $`\omega _k=Dk^2`$ which is valid if the wavelength of the perturbation $`2\pi /k\lambda _0`$. Knowing $`k`$, and calculating $`\omega _k`$ by a fit of the measured amplitude of the modulation to a decaying exponential, we obtain $`D`$. This is illustrated in Fig. 3 for a case where $`2\pi /k7\lambda _0`$. We deduce from this measurement that

$$D=D_{}+D_{}$$
(6)
where $`D_{}`$ is a positive contribution responsible for the overstabilization of the array. The latter is calculated by taking the difference between $`D`$, and $`D_{}`$ evaluated via Eq. 4 using the value $`K_1=1.9\times 10^3Ks\mu m^2`$ obtained from our present experiments.

Next, we simulate a phase-field model of a two-component (AB) eutectic alloy which is completely symmetric under the exchange of $`\alpha `$ and $`\beta `$. Our goal here is not to model quantitatively the experiments but to demonstrate the generality of our results in different alloy systems. An order parameter $`\varphi `$, which distinguishes between solid ($`\varphi =+1`$) and liquid ($`\varphi =1`$), is coupled to the dimensionless concentration field, $`u(cc_E)/(\mathrm{\Delta }c/2)`$, where $`c`$ is the mole fraction of B, and $`\mathrm{\Delta }c`$ is the difference between the solid composition in $`\beta `$ and $`\alpha `$ at $`T_E`$. The non-conserved and conserved dynamics for $`\varphi `$ and $`u`$, $`\tau _t\varphi =\delta F/\delta \varphi `$, and $`_tu=\stackrel{}{}(M(\varphi )\stackrel{}{}\delta F/\delta u)`$, respectively, are derived from the functional

$$F=_V𝑑V\left[\frac{W_u^2}{2}|u|^2+\frac{W_\varphi ^2}{2}|\varphi |^2+f(\varphi ,u,T)\right]$$
(7)
where $`V`$ is the volume of the system and

$$f(\varphi ,u,T)=f_{DW}(\varphi )+\frac{1+h(\varphi )}{2}f_s+\frac{1h(\varphi )}{2}f_l,$$
(8)
is the bulk free-energy density that is the sum of the standard double-well $`f_{DW}=\varphi ^2/2+\varphi ^4/4`$, and a concentration-dependent part that interpolates between the bulk free-energy of the liquid, $`f_l=u^2/2`$, and the solid $`f_s=(u^21)^2/8(T_ET)/T_E`$, with $`h(\varphi )=3(\varphi \varphi ^3/3)/2`$. Furthermore, we used the mobility $`M(\varphi )=D_l[1(1+\varphi )^4/16]`$, where $`D_l`$ is the solute diffusivity in the liquid, that makes the diffusivity in the solid vanish and yields efficient computations . Directional growth is implemented using the frozen-temperature approximation $`T(z,t)=T_E+GzVt`$, and periodic boundary conditions in $`x`$ are used in all simulations.

Steady-state $`\mathrm{\Delta }T(\lambda )`$ curves were obtained from short simulations with two lamellae. The stability was studied with long simulations where steady-state arrays of up to 20 lamellae, constructed from the two-lamella solutions, are slightly perturbed by a long-wavelength modulation of spacing of the form of Eq. 5. By exponential fits of the amplitude of modulation vs time for different $`k`$, we obtain the stability spectrum $`\omega _k`$ as illustrated in Fig. 4, and hence $`D`$ by a quadratic fit of $`\omega `$ vs $`k`$ at small $`k`$. The stability limit $`\lambda _c`$ is then obtained by determining where $`D`$ changes sign and values of $`D_{}`$ are obtained by subtracting $`D_{}`$ from $`D`$. We find that the dimensionless ratio $`D_{}/(V\lambda _0)`$ varies with $`\mathrm{\Lambda }_0`$, but negligibly with $`G`$ and $`V`$. Moreover, we find that the simple form

$$D_{}/(V\lambda _0)A\mathrm{\Lambda }_0,$$
(9)
with $`A0.15`$ gives a reasonable fit to our phase-field simulation results as shown in Fig. 5. Remarkably, our experimentally determined values of $`D_{}`$ are close to those obtained in the phase-field simulations, even though the two alloy systems are different.

To interpret our findings, let us briefly review Langer’s analysis that yields $`D=D_{}`$. Its first ingredient is the assumption that the interface adjusts adiabatically its average temperature to the local spacing, or

$$\mathrm{\Delta }T[\lambda (x,t)]G\zeta (x,t),$$
(10)
where $`\mathrm{\Delta }T(\lambda )`$ is the same as in steady-state and $`\zeta (x)`$ is the $`z`$ coordinate of the envelope of the eutectic interface, defined as a smooth curve interpolating the positions of the three-phase junctions (trijunctions), with the origin at $`T_E`$. The second is Cahn’s assumption that lamellae grow normal to this envelope, which for a small perturbation is equivalent to

$$_ty(x,t)V_x\zeta (x,t),$$
(11)
where $`y(x,t)`$ is the lateral displacement of trijunctions from their steady-state positions. Finally, it follows from the definition of $`y`$ that $`\lambda (x,t)\lambda _0(1+_xy)`$. Differentiating both sides of this identity with respect to time and using Eqs. 1, 10, and 11, one obtains the diffusion equation (3) with $`D=D_{}`$.

We checked that Eq. 10 is indeed faithfully obeyed in the range of wavelengths that we consider here. Consequently, the discrepancy between $`D`$ and $`D_{}`$ must originate from a correction to Cahn’s normal growth assumption. It is simple to show that the modified phenomenological form

$$_ty(x,t)V_x\zeta (x,t)+D_{}_x\lambda (x,t)/\lambda _0,$$
(12)
yields the diffusion equation (3) with $`D`$ given by Eq. 6. Equation 12 implies that trijunctions also move locally “parallel” to the envelope of the eutectic interface in response to a gradient of spacing. To see physically why this lateral motion overstabilizes the pattern, consider a local depression in an array of initial spacing $`\lambda _0`$. Cahn’s normal growth assumption implies that such a depression will produce a local decrease of spacing, and hence a local increase in undercooling that will amplify this depression if $`\lambda _0<\lambda _m`$ (because $`d\mathrm{\Delta }T/d\lambda <0`$ in this case). This well-known argument yields $`\lambda _c=\lambda _m`$. In contrast, the second term on the right-hand-side of Eq. 12 implies that the lateral motion of trijunctions opposes the local decrease in $`\lambda `$, and hence helps flatten the interface.

A prediction for $`\mathrm{\Lambda }_c=\lambda _c/\lambda _m`$ can be obtained by setting $`D=D_{}+D_{}=0`$, which, using Eqs. 2, 4, and 9, yields the cubic equation

$$1\frac{1}{\mathrm{\Lambda }_c^2}+\frac{AG}{K_1V}\mathrm{\Lambda }_c=0.$$
(13)
With $`K_1=1.9\times 10^3Ks\mu m^2`$ and $`A=0.15`$, we obtain $`\mathrm{\Lambda }_c=0.70`$ for the experiment of Fig. 2; the lowest observed spacings are just above the predicted stability threshold, as it should be. For the phase-field simulations ($`K_1=0.1120`$), we find $`\mathrm{\Lambda }_c=0.942`$ for $`l_T/l_D=20`$ and $`\mathrm{\Lambda }_c=0.715`$ for $`l_T/l_D=2`$.

Two previous stability analyses have predicted that $`\lambda _c`$ should be smaller than $`\lambda _m`$. The one by Caroli et al. , however, is restricted to a large $`G`$ limit that cannot be compared with our results. The other by Chen and Davis  does not have this restriction, but predicts a departure of $`\lambda _c`$ from $`\lambda _m`$ that is about one order of magnitude smaller than found here and predicted by Eq. 13. We believe that the lateral motion of the trijunctions is due to a coupling between the diffusion field and the non-planar front geometry on the scale of the individual lamellae; such effects would appear only at higher orders in the analyses cited above. Therefore, an analytical understanding of eutectic stability at small spacings from sharp-interface models largely remains to be developed.

The present results shed light on a number of previous observations. In metallic eutectics, $`K_1`$ is generally close to $`10^2Ks\mu m^2`$ , so that the value of $`V`$ below which the departure of $`\mathrm{\Lambda }_c`$ from unity becomes significant is of about $`1\mu ms^1`$ for $`G`$ in the $`100Kcm^1`$ range. This may explain the deviation from the law $`\overline{\lambda }V^{0.5}`$, where $`\overline{\lambda }`$ is some empirically defined average eutectic spacing, that has been observed at $`V`$ lower than about $`1\mu ms^1`$ in a number of metallic eutectics . Similarly, the overstability due to the lateral motion of the trijunctions may explain why coupled growth in a peritectic system has recently been found to be stable  in a situation, analogous to that of eutectics at $`\lambda <\lambda _m`$, where the interface should be unstable according to the JH-Cahn stability arguments . Finally, our results also improve our understanding of the morphological instability that leads to the formation of eutectic colonies in the presence of a dilute ternary impurity .

This work was supported by Centre National d’Etudes Spatiales, France, and by the U.S. DOE under grant No. DE-FG02-92ER45471.


text:
# Nonlinearity and Multifractality of Climate Change in the Past 420,000 Years

\[

## Abstract

Evidence of past climate variations are stored in ice and indicate glacial-interglacial cycles characterized by three dominant time periods of 20kyr, 40kyr, and 100kyr. We study the scaling properties of temperature proxy records of four ice cores from Antarctica and Greenland. These series are long-range correlated in the time scales of 1-100kyr. We show that these series are nonlinear as expressed by volatility correlations and a broad multifractal spectrum. We present a stochastic model that captures the scaling and the nonlinear properties observed in the data.

\]

Abundant geological evidence indicates that temperatures varied from the cold of ice ages to the warmth of interglacial periods. In the last 800,000 years (800kyr) there is strong evidence for a dominant glacial-interglacial cycle of 100kyr, with weaker secondary cycles of 40kyr and 20kyr . Each 100kyr cycle consists of gradual cooling for $``$ 90kyr followed by rapid warming during $``$ 10kyr. “Milankovitch forcing”, which refers to changes in insolation (solar radiation) due to variations in the precession, obliquity (tilting), and eccentricity of Earth’s orbit  are thought to play an important role in glacial dynamics. These orbital variations are characterized by periods of 20kyr, 40kyr, and 100kyr, respectively. The 20kyr and 40kyr periods in the climate records are generally believed to be a linear response of the climate system to insolation variations. In contrast, the weakness of the variations in solar radiation at the 100kyr timescale has lead to the generally accepted conclusion that the glacial-interglacial oscillations at this timescale are most likely not a direct linear response of the climate system to external solar variations .

Many deterministic theories have been developed to explain the glacial-interglacial 100kyr variability; the majority suggest that the 100kyr period is a result of self-sustained nonlinear mechanisms (see, e.g., ). Other studies proposed that climate variations are stochastic and follow scaling laws — the Milankovitch periods are second-order perturbations (e.g. ). Importantly, both the deterministic and stochastic mechanisms still assume that the variations on time scales below 100kyr down to 10kyr are linear. The objectives of the present study are to quantify the degree of nonlinearity of climate dynamics within the time scales of 1-100kyr and to provide statistical characteristics of the proxy records which can serve as a test for distinguishing between existing climate models .

We study the correlation (scaling) properties of climate records of the past 420kyr. We show that temperature variations are long-range correlated suggesting that the Milankovitch periods are indeed secondary and (contrary to common belief ) that climate dynamics of all time scales below 100kyr down to 1kyr are highly nonlinear. In addition, we quantify the degree of nonlinearity in the climate records and suggest a possible stochastic nonlinear mechanism for our findings.

Our analysis is based on isotope records obtained from four ice cores, Vostok and Taylor Dome from Antarctica, and GISP (Greenland-Ice-Sheet-Project) and GRIP (Greenland-Ice-Project) from Greenland . Measurements of oxygen and hydrogen isotope ratios ($`\delta ^{18}O`$ and $`\delta D`$) of the ice at different depths in the core provide a proxy record of temperature  when the ice was formed (Fig. 1a). These records extend back to 100-420kyr.

Fourier analysis is the standard method for studying long-range correlations in time series. When the power spectrum follows scaling laws, $`S(f)1/f^\beta `$ ($`\beta >0`$), the series is long-range correlated . However, the power spectrum might yield an inaccurate estimation of the scaling exponent due to constant or polynomial trends that are not necessarily related to the intrinsic dynamics . We therefore use the detrended fluctuation analysis (DFA) ; the $`m\mathrm{th}`$ order DFA eliminates polynomial trends of order $`m1`$ from the data and provides a more accurate estimation of the scaling exponent . If the root mean square fluctuation function, $`F(n)`$, is proportional to $`n^\alpha `$, where $`n`$ is the window scale, the series is long range correlated ($`\beta =2\alpha 1`$). For a random series $`\alpha =0.5`$ while for correlated (or anticorrelated) series $`\alpha >0.5`$ (or $`\alpha <0.5`$). We begin our analysis with the Vostok ice core and find that temperature changes are highly correlated in the time range 1-100kyr with a scaling exponent $`\alpha 1.5`$ (Fig. 2a), consistent with the previously reported power spectrum exponent $`\beta =2`$ .

Next, we analyze the nonlinear properties of the ice core record. We define a process to be linear if it is possible to reproduce its statistical properties (such as the third moment) from the power spectrum and the probability distribution alone, regardless of the Fourier phases . This definition includes autoregression processes ($`x_n=_{i=1}^Ma_ix_{ni}+_{i=0}^Lb_i\eta _{ni}`$ where $`\eta `$ is Gaussian white noise) and fractional Brownian motion; the output, $`x_n`$, of these processes may undergo monotonic nonlinear transformations $`s_n=s(x_n)`$ and still be linear. Processes which are not linear are defined as nonlinear .

Long-range correlations in the temperature time series, $`T_i`$, reflect linear aspects of $`T_i`$. Long-range correlations in the magnitudes of temperature increments, $`|\mathrm{\Delta }T_i||T_{i+1}T_i|`$ (Fig. 1b), which we define as volatility, indicate nonlinearity of the underlying process . Linear series have uncorrelated $`|\mathrm{\Delta }T_i|`$ series while nonlinear processes that follow a scaling law exhibit long-range correlations in the magnitude series $`|\mathrm{\Delta }T_i|`$. We find that the magnitude series $`|\mathrm{\Delta }T_i|`$ is highly long-range correlated within the time range 1-100kyr (Fig. 2b). Thus, the underlying process is nonlinear . The value of the correlation exponent quantifies the degree of nonlinearity in the ice core record. Correlations in the magnitude series indicates that the magnitude series is “clustered”, i.e., large magnitude is more likely to be followed by a large magnitude, as can be seen in Fig. 1b. These clusters may be associated with abrupt warming events known as Dansgaard-Oeschger events .

To demonstrate that the correlations in the magnitude series are related to the nonlinearity of the underlying process we apply a surrogate data test for nonlinearity that preserves both the power spectrum and the histogram of the temperature increment series $`\mathrm{\Delta }T_i`$ . The surrogate series has random Fourier phases; the nonlinearities that are stored in the phases are destroyed. We find that the magnitude series obtained from the surrogate series is indeed uncorrelated (Fig. 2b) confirming that the original series is nonlinear within 1-100kyr.

Correlations in the magnitude series $`|\mathrm{\Delta }T_i|`$ can be related to the width of the multifractal spectrum . We calculate the exponents $`\tau (q)`$ of different moments $`q`$ for the ice core data and find that $`\tau (q)`$ is a nonlinear function of $`q`$ (Fig. 2c), indicating that the temperature series is multifractal. We also perform multifractal analysis on the surrogate data and find that its $`\tau (q)`$ is almost linear. The multifractal spectrum, $`D(h)`$ (Fig. 2d), is broad for the original data and narrower for the surrogate data. The broadness of the multifractal spectrum may also be used to quantify the degree of multifractality, and thus the degree of nonlinearity, in the data .

We repeat the above analysis for the other three ice cores and obtain similar results (Table I). Although the DFA exponents of the original series are smaller for the Greenland cores, the magnitude series exponents are almost the same for all cores. We thus conclude that climate dynamics is nonlinear for time scales of few thousands of years up to 100kyr.

To understand the mechanism that may contribute to the nonlinearities observed in the data we modified a model for ice-volume evolution recently suggested by Wunsch . Ice volume $`V`$ was observed to be negatively correlated with temperature $`TV`$  and thus a model for ice accumulation may serve, indirectly, as a model for temperature dynamics. The Wunsch model can be summarized as follows: The ice-sheet builds randomly up to a critical volume where it breaks up rapidly. Then, growth begins again. By construction, this model is a random walk up to a time scale corresponding to the critical ice volume, followed by a crossover to random behavior for larger time scales. However, this model does not reproduce the nonlinearity in the data as defined and found above.

The assumptions of our model (Fig. 3) are:

(i) The ice volume $`V`$ changes with steps $`\delta +b/V`$; i.e., $`V(t+dt)=V(t)+\delta (t)+b/V(t)`$.

(ii) When ice volume $`V`$ “crosses” a critical volume $`\xi `$, $`b`$ is set to be negative, $`b=b_2<0`$. Ice volume is considered to lie between $`0.01\xi `$ and $`\xi `$. When $`V=0.01\xi `$, $`b`$ becomes $`b=b_1>0`$ till $`V`$ exceeds the threshold $`\xi `$.

(iii) The ice accumulation increments $`\delta `$ are the product of two stochastic inputs, $`\delta (t)=\zeta (t)\eta _{i(t)}`$. $`\zeta `$ and $`\eta `$ are Gaussian distributed random variables with zero mean and unit variance.

(iv) Random switching between the states $`\eta _i`$’s is controlled by $`i(t)`$ which is equal to $`[l(t)]`$ where $`[]`$ stands for the closest integer value. $`l`$ is a random walk described by, $`l(t+dt)=l(t)+C\omega (t)`$, where $`C`$ is the switching range and $`\omega `$ is another Gaussian random variable with zero mean and unit variance (see Fig. 3a).

Assumptions (i) and (ii) describe the random growth (with $`b=b_1>0`$) of the ice-sheet and its rapid breakup (with $`b=b_2<0`$) after crossing the critical volume $`\xi `$. The $`1/V`$ term in assumption (i) mimics the reduced ice accumulation for large ice-volume . Ice volume changes $`\delta `$ result from two interacting random inputs \[assumption (iii)\] where one, $`\zeta `$, may represent the atmosphere, the other, $`\eta `$, the ocean, and the product, $`\zeta \eta `$, the atmosphere-ocean interaction.

In our simulations of the model we use the following values, $`dt=0.1kyr`$, $`\xi =1.5\sqrt{90kyr/dt}`$, $`b_1=1`$, $`b_2=3`$, and $`C=0.27`$. We choose the value 90kyr in $`\xi `$ so that on average $`V`$ will grow from zero to $`\xi `$ after 90kyr/dt steps. We choose the values of $`b_1`$ and $`b_2`$ such that the ice-sheet grows slowly and breakup rapidly. The values of $`\xi `$, $`b_1`$, and $`b_2`$ are constrained by the natural record. The switching parameter $`C`$ determines the number of states $`\eta `$ for a given number of steps; the number of different states is proportional to the square root of the number of steps (e.g., $`3`$ states for 4kyr and $`5`$ for 10kyr).

An example of an arbitrary 400kyr time series obtained by the model is shown in Fig. 3b. The scaling of the model’s $`V`$ series (Fig. 4a) indicates random walk behavior with exponent $`\alpha =1.5`$ (as for the Vostok core). The magnitude series $`|\mathrm{\Delta }V|`$ is highly correlated with exponent $`0.8`$ (Fig. 4b) as for the ice core data (Fig. 2b and Table I). The surrogate data test applied to the $`\mathrm{\Delta }V`$ series changes the magnitude series into an uncorrelated one, indicating the nonlinearity of the model. This nonlinearity is mainly due to the product of the inputs $`\eta `$, $`\zeta `$ in assumption (iii). The multifractal spectrum is broad where, as with the ice core data (Fig. 2c,d), the exponents for negative moments, $`\tau (q<0)`$, mainly contribute to its broadness (Fig. 4c,d). After the surrogate data test the series becomes linear and statistically different from the original data.

This simple model reproduces the statistical characteristics of the ice core data under consideration. Although the natural system is undoubtedly more complex, we conjecture that the model variables may be associated with specific aspects of Earth’s climate system although our model cannot uniquely identify them. One of the random inputs, $`\eta `$, thus may represent the influence of the deep ocean on ice accumulation since the state of the deep ocean is known to have impact on glaciation (e.g., ). The other random input, $`\zeta `$, may represent the net atmospheric influence affecting ice accumulation (resulting from, e.g., variations in eddy transport, cloudiness, incoming solar radiation, ablation). We assume that the ocean has several states with a tendency to return to previous states, as does $`\eta `$ in the model; a possible example for such “switching” mechanism is the deep ocean circulation which has few states with possible switching between them (e.g. ).

We conclude that climate changes in the time range of 1-100kyr are long-range correlated confirming the major role of stochasticity in climate . Moreover, our results suggest that the underlying dynamics in the time scales of 1-100kyr is nonlinear. This nonlinearity is specified and quantified by strong long-range correlations in the magnitudes of temperature changes and in a broad multifractal spectrum. Our simple stochastic model suggests that the nonlinearity can be the result of only two random processes that interact with each other. Glaciation models may be generally categorized into two main alternatives: (i) linear mechanisms that are driven by stochastic forcing (e.g. ), and (ii) nonlinear mechanisms without stochastic forcing (e.g. ). Our results and model suggests a third alternative — nonlinear mechanism that inherently involves stochastic forcing. Our results raise a new challenge for the many climate models, and may help guide development of better climate models, which include both periodic and stochastic elements of climate change.

YA and HG thank the Bikura fellowship for financial support. DRB thanks H.E. Stanley for his generous hospitality. We thank P. Cerlini, P. Hybers, V. Schulte-Frohlinde, P.H. Stone, E. Tziperman, and C. Wunsch for helpful discussions.