text:
# On the photometric variability of blue supergiants in NGC 300 and its impact on the Flux-weighted Gravity–Luminosity Relationship1

## 1 Introduction

The photometric and spectral variability has been established as a characterizing feature of supergiant stars already a few decades ago (Abt 1957, Rosendhal & Wegner 1970, Maeder & Rufener 1972, Rosendhal 1973, Sterken 1977), affecting virtually all luminous stars during their post-main sequence evolution. The light variations appear to be cyclic, characterized by semi-regular, rather than strictly periodic, patterns, obeying a semi-period–luminosity–color relationship (Burki 1978, Maeder 1980) analogous to the relation valid for Cepheids, with timescales ranging from a few days (early B supergiants) to hundreds of days (red supergiants). The maximum light amplitudes, on the order of 0.1 mag, increase with stellar luminosity, and reach in the hot star domain a local maximum for early-B supergiants (Appenzeller 1972, Maeder 1980, van Genderen 1989). The presence of photometric micro-variations on scales $`\mathrm{\Delta }m<0.2`$ mag typifies variable massive OBA supergiants, or $`\alpha `$ Cyg variables following the General Catalogue of Variable Stars (Kholopov et al. 1985; see also reviews by Sterken 1989 and van Genderen 1991, 2001).

The existence of a semi-period–luminosity–color relationship suggests pulsational instabilities as the driving mechanism for the observed variability, where stochastic processes determine the quasi-regularity. Several clues, including the large discrepancy between the theoretical periods for radial pulsations in the fundamental mode and the observed ones (Percy & Welch 1983, van Genderen 1985), together with color variations smaller than predicted (van Genderen 1986), indicate that non-radial oscillations are at work in blue supergiants (Maeder 1986, Baade 1992), possibly appearing in connection with evolutionary stages past the red supergiant phase (Lovey et al. 1984, Schaller 1990). Alternative or additional mechanisms have been proposed to explain the observed photometric and spectroscopic variability, such as orbital motions in binary systems (Harmanec 1987). In the case of late-B and early-A supergiants the variability of spectral line profiles originating in the photosphere is consistent with non-radial pulsation modes, while rotational modulation appears to be affecting the extended stellar wind regions (Kaufer et al. 1996, 1997).

Regardless of the physical reasons for the observed variability, which concerns both the photospheric layers (as deduced from the photometric variations) and the extended winds of blue supergiants (from time series spectroscopic investigations, e.g. Kaufer et al. 1996, Rivinius et al. 1997, Stahl et al. 2003), it is interesting and necessary to determine its relevance when considering the use of blue supergiant stars as stellar candles (Bresolin 2003). Recently Kudritzki, Bresolin & Przybilla (2003) have discussed the preliminary calibration of the Flux-weighted Gravity–Luminosity Relationship (fglr) and its potential use as an extragalactic distance indicator. As originally proposed, such a relation applies to late-B to early-A supergiants, which belong to the visually brightest stars in galaxies, thus allowing tests of its validity to be carried out in galaxies located well beyond the Local Group. Both the fglr and the Wind Momentum–Luminosity Relationship (wlr, Kudritzki et al. 1999) rely on the measurement of the apparent visual magnitude of blue supergiants, which are then transformed to bolometric luminosities from knowledge of the spectral type. An assessment of the reliability of such methods should involve tests concerning their sensitivity to time-dependent stellar properties, such as the mass-loss rate (for the wlr), and to the above mentioned photometric micro-variations (fglr and wlr).

‘Normal’ (non-LBV) blue supergiants have been rarely monitored for variability beyond the Magellanic Clouds. In this paper we report on a medium-term (five months) V monitoring program of blue supergiants in the galaxy NGC 300. The main motivation for this work is to demonstrate that the fglr is insensitive for all practical purposes to the micro-variations normally observed in B and A supergiants.

## 2 Observational data

The photometric measurements for the NGC 300 blue supergiants included in this study have been extracted from images obtained at the ESO/MPI 2.2 m telescope on La Silla, equipped with its Wide Field Imager 8 K $`\times `$ 8 K mosaic camera. These multi-epoch observations are part of the BVRI Cepheid monitoring program described by Pietrzyński et al. (2002), where details on the calibration and the data reduction can be found. The data discussed in the current paper span the period 1999 July 31–2000 January 8, covered by observations on 28 different nights, mostly under photometric conditions, roughly grouped in five blocks of data separated by temporal gaps up to 43 days long. Only the V-band temporal sequence will be discussed here (B-band magnitudes have a somewhat less reliable calibration, but still show variations similar to those observed in V).

Spectra of blue supergiants in NGC 300 have been obtained with the FORS1 multi-object spectrograph at the ESO Very Large Telescope by Bresolin et al. (2002a), and we restrict the following analysis to the spectroscopically confirmed B- and A-type supergiants. For the individual stars we adopt the nomenclature introduced in the spectral catalog presented in that paper.

## 3 Variability of blue supergiants: detection and light curves

The standard deviation from the mean magnitude and the maximum light amplitude were taken as a measure of the stellar variability. We selected stars having small dispersion in magnitude to serve as comparison stars in order to define the observational uncertainty, $`\sigma _0`$, as a function of magnitude. We then arbitrarily retained for further analysis only those objects for which $`\sigma 2\sigma _0`$ (Fig. 1). A selection based on the maximum light amplitude generates a similar sample of stars. This procedure is justified since we are not aiming at a complete compilation of variable supergiants in NGC 300, but instead at an analysis of those spectroscopically confirmed supergiants displaying the largest variability, since these will have a maximum effect on the fglr.

The 16 stars satisfying the selection criterion are labeled in Fig. 1 with the identification number introduced by Bresolin et al. (2002a). Their mean apparent V and absolute $`M_V`$ magnitudes, $`BV`$ colors, standard deviations in the $`V`$ magnitude $`\sigma _\text{V}`$, maximum light amplitudes A$`_\text{V}`$ and spectral types are summarized in Table 1. $`M_V`$’s have been determined assuming a distance modulus $`mM=26.53`$ (Freedman et al. 2001), and estimating the reddening (and the extinction with $`R_V=3.1`$) based on the difference between the observed $`BV`$ and the expected color index for the given spectral type. The photometric measurements differ slightly (mostly at the 0.01-0.02 mag level) from those in Table 2 of Bresolin et al. (2002a), as a result of improved zero-points and photometry. All stars in Table 1 are blue supergiants (early-B to mid-A), with the exception of A5, a red supergiant, which shows the largest light variation among the detected variables, apparently occuring over a period longer than our observing window. The WN11 (or Ofpe/WN9) star B16 analyzed by Bresolin et al. (2002b) also appears in Table 1.

The V light curves of the selected variable stars are displayed in Fig. 2. Despite the gaps in the temporal sequence, photometric variations occuring on timescales of tens of days can be easily detected in most cases, often superimposed on variations at higher frequency.

## 4 Period search

The periodicity of the photometric sequences has been analyzed utilizing the Phase Dispersion Minimization algorithm introduced by Stellingwerf (1978), as implemented in the iraf astronomical data reduction package. Our goal was simply to check whether semi-periods comparable in length to those of Galactic and Magellanic Cloud supergiants could be found also in NGC 300, aware of the fact that the irregular sampling of the data within our observing window is far from being optimal for an accurate frequency analysis.

Only for two stars a clear isolated minimum in the $`\mathrm{\Theta }`$ statistics (which is defined to vary between 0 and 1, see Stellingwerf 1978) is unambiguously identified, leading to a reasonable determination of the semi-periods. Such is the case of D12 ($`\mathrm{\Theta }=0.06`$, $`\mathrm{P}=72`$ days) and A10 ($`\mathrm{\Theta }=0.40`$, $`\mathrm{P}=96`$ days), both of spectral type A2 Ia, shown in the phase diagrams of Fig. 3 (the uncertainty in the quoted periods is on the order of 10%). For the remaining stars in the sample $`\mathrm{\Theta }`$ lies between 0.4 and 0.75, and a unique semi-period determination was not possible, with multiple phase dispersion minima occuring at periods roughly between a few days and 30-40 days. The phased light curves for these stars display a large amount of scatter, but this is typical of the light variations of blue supergiants, as multiple pulsation modes can be excited in blue supergiants (Lucy 1976, van Gent & Lamers 1986, van Genderen, Sterken & de Groot 1998). We note that the semi-periods found for the two A2 supergiants D12 and A10, of absolute magnitude $`M_V=8.3`$ and $`7.8`$, respectively, are roughly twice as long as those listed for similar stars in the Milky Way by Burki (1978). Fundamental periods calculated from the respective stellar parameters are 42 days (D12) and 27 days (A10) (Schaller 1990).

## 5 Impact on the Flux-weighted gravity–Luminosity Relationship

The findings of the previous sections confirm that the typical photometric micro-variability of blue supergiants, known from studies in the Milky Way and the Magellanic Clouds, is also found in our NGC 300 sample. The advantage of these data from the extragalactic point of view is that they allow us to carry out a simple but important test concerning the effect of the stellar variability on the fglr. Given the typical amplitude of the broad-band photometric variations, $``$0.1 mag, it is reasonable not to expect a significant effect, but a quantification of its size is desirable.

The fglr, introduced by Kudritzki, Bresolin & Przybilla (2003), expresses the proportionality between the bolometric magnitude of a blue supergiant and its flux-weighted gravity during its evolution across the Hertzprung-Russell diagram at roughly constant luminosity:

$$M_{\text{bol}}=3.71\mathrm{log}(g/T_{\text{eff,4}}^4)13.49$$
(1)
where $`T_{\text{eff,4}}=T_{\text{eff}}/10^4K`$, with slope and zero-point fixed by a preliminary calibration, based on the stellar parameters determined for stars in different galaxies: Milky Way, Magellanic Clouds, NGC 6822, M31, M33, NGC 300 and NGC 3621.

It is not clear whether during the semi-periodic light variations the bolometric magnitude remains rigorously constant. If that were the case, then the fglr would be unaffected by the blue supergiant variability. On the other hand, one could imagine that the detected light variations are symptomatic of small changes in the total luminosity (as in the pulsating massive star models by Dorfi & Gautschy 2000). In this case, and lacking additional constraints, we can test the corresponding effect on the fglr by determining the relationship at different epochs, letting $`M_{\text{bol}}`$ vary by an amount equal to the broad-band variability. The flux-weighted gravity $`g/T_{\text{eff,4}}^4`$ is kept constant, since spectra for the derivation of gravity and $`T_{\text{eff}}`$ are available for one single epoch only, and therefore we have no way of estimating their possible variations. Only spectral time series would allow us to verify if and how the stellar parameters change during the cyclic changes.

Under these assumptions, we have calculated slope and zero-point of the fglr at each available epoch, using our photometric temporal sequence for the NGC 300 blue supergiants considered by Kudritzki et al. (2003), i.e. 14 stars of spectral type between B8 and A2. The stellar parameters for these stars, as calculated by Kudritzki et al. (2003), are summarized in Table 2, and the whole set of $`V`$ magnitudes at 28 different epochs is given in Table 3 (the uncertainty of the individual measurements is on the order of 0.01–0.02 mag). The zero-point and rms of the linear fit are found to be affected below the 10% level, when compared with the regression shown for NGC 300 stars alone ($`rms`$0.2 mag) by Kudritzki et al. (2003), with variations in the slope only up to 4%. It is clear that the supergiant variability introduces a negligible effect on the resulting fglr, which will remain negligible even when an improved calibration of the relation, including its possible dependence on the stellar metallicity, will become available in the near future.

In conclusion, the fglr appears to provide an extragalactic distance indicator which is robust against supergiant photometric variability, provided that the number of stars considered per galaxy is sufficiently large ($`n`$10). As for the future work, besides improving the calibration in terms of statistics and metallicities, additional efforts are required to establish whether spectral variations are important for the derivation of the stellar parameters.

WG gratefully acknowledges financial support for this work from the Chilean Center for Astrophysics FONDAP 15010003. Support from the Polish KBN grant 2P03D02123 is acknowledged.


format: text

text:
# Anomalous temperature-dependent transport in YbNi2B2C and its correlation to microstructural features

## I Introduction

The transport properties of heavy fermion intermetallic systems often display peculiar behaviorsStewart (1984) whose origin remains as yet unestablished and are subject of current interest. In these systems the effective mass of conduction electrons is enhanced as a result of hybridization between localized electronic orbitals and delocalized conduction bands, and apparently these hybridized states tend to be particularly sensitive to crystalline disorder, manifesting strong sample-to-sample variations in the transport propertiesStewart (1984); Franz et al. (1978) which go beyond trivial differences in residual resistivity.

The quaternary compound YbNi2B2C is a ytterbium-based heavy fermionYatskar et al. (1996); Dhar et al. (1996) with an electronic specific heat coefficient, $`\gamma 500`$ mJ/mol K2, and a Kondo temperature, $`T_K10`$ K: a temperature scale that is conveniently isolated from other characteristic temperatures such as superconducting or magnetic condensates ($`T_c,T_N<0.03`$ K, if they exist at all), and crystal electric field splitting ($`T_{CEF}100`$ K)Yatskar et al. (1996); Dhar et al. (1996); Grasser et al. (1996); Rams et al. (2000); Boothroyd et al. (2003), and therefore favorable for the study of Kondo physics. Furthermore, an investigation of the effects of annealing on the resistive behavior of boro-carbide single crystalsMiao et al. (2002a) found that YbNi2B2C displayed radical changes in temperature-dependent resistance below room temperature, pointing to the possibility that the transport properties of this compound could be “tuned” to a certain extent by annealing, and thus offer a model system to study the relationship between disorder and hybridized states near the Fermi level.

In a previous reportAvila et al. (2002), we demonstrated a significant variation in the temperature dependent electrical resistivity and thermoelectric power between as-grown crystals and crystals that had undergone annealing at $`950^{}`$C, whereas the thermodynamic properties (heat capacity and magnetic susceptibility) remained almost unchanged. We interpreted these results in terms of redistributions of local Kondo temperatures associated with ligandal disorder for a small (on the order of 1%) fraction of the Yb sites. This hypothesis left open two obvious questions: (i) What happens to the electrical resistivity for temperatures greater than the maximum possible $`T_K`$ in this proposed distribution? and (ii) What is the microscopic origin of this distribution of Kondo temperatures?

In the present work, we first extend the temperature range of resistivity measurements both higher (up to 1000 K) and lower (down to 0.4 K), in order to fully characterize the transport behavior of as-grown and optimally annealed crystals and provide further evidence of the proposed model. We then present x-ray single-crystal refinements and transmission electron microscopy experiments, which provide new information on the types of disorder that may be causing the anomalous Kondo temperatures in the surrounding Yb sites.

## II Experimental Details

Single crystals of YbNi2B2C were grown from Ni2B flux at high temperatures as described elsewhereCanfield et al. (1998); Yatskar et al. (1996). For the resistivity measurements, a particularly clean and well-formed single crystal plate measuring approximately $`7\times 3\times 0.2`$ mm3 was selected, polished on both surfaces to remove most of the attached flux, and cut using a wire saw into flat bars of typical dimensions $`2\times 0.5\times 0.13`$ mm3 with the length along the  direction. Electrical contacts were placed on three unannealed samples in standard 4-probe linear geometry, using Pt wires attached to a sample surface with Ablebond 88-1 silver epoxy. For each slab, the sample weight and dimensions were carefully measured and an evaluation of the sample densities ($`8.3\pm 0.2`$ g/cm3) was used to estimate an upper limit of $`\pm 10\%`$ for the overall geometry-related uncertainty in calculating resistivity $`\rho =RA/d`$, where $`A`$ is the cross-sectional area and $`d`$ is the distance between voltage leads.

Two measured samples were then selected to undergo annealing for comparison. The Pt wires were removed and the silver epoxy was polished off, leading to a small ($`10\%`$) decrease in the sample’s original weight due to an intentional, slight overpolishing. Since the polishing occurred on a single flat surface parallel to the current direction, the reduction in the sample cross-section was assumed to be the same as the one in mass. The samples were then annealed in vacuum for 150 hours at $`950^{}`$C (details of the annealing procedure are described in ref. Avila et al., 2002).

Electrical resistance measurements below room temperature were performed on commercial Quantum Design PPMS systems, allowing measurements down to 1.8 K, or down to 0.4 K if the 3He option is installed. Above room temperature, measurements were performed on a custom-built high temperature resistivity system (HTRS) which can be mounted on the same system used for sample annealing under vacuum ($`10^6`$ Torr). A quartz probe is inserted in the system with four leads where the other end of the sample contacts are attached using silver epoxy. Temperature is measured with two independent thermometers: Pt resistance temperature detector (RTD) and type-R thermocouple. A 1 mA current is applied on the sample by a Keithley 220 current source and the voltage is read on a HP 34420A nanovoltmeter. The RTD current is applied by a LakeShore 120 current source and the voltage is read on a HP 34401A multimeter. The thermocouple voltage is read by another HP 34420A nanovoltmeter. All these instruments are GPIB interfaced with a computer running a custom-made software for data acquisition. To reduce noise and avoid diode or thermoelectric effects, five readings are taken with the current in each direction, and the ten absolute voltage values are averaged.

A typical HTRS experiment proceeds as follows. 1) The system is heated to 450 K for about an hour in air to cure the contacts, then cooled back to room temperature. 2) The vacuum system is turned on and allowed to achieve a dynamic vacuum of $`10^5`$ Torr or better. 3) The furnace temperature is ramped at 2 K/min or slower up to about 600 K and back to room temperature, while the software acquires data points every 2 min. This first ramp is used to allow further stabilization of the contacts which often undergo small changes when heated, marked by irreproducible shifts in the resistance behavior. 4) The definitive dataset is acquired by ramping at 2 K/min or slower up to about 1000 K and back to room temperature, while the software acquires data points every 2 min. For the same reason of contact stability and in order to obtain a good match between the curves from both measurement systems, the low temperature measurements in the PPMS were performed only after the high temperature ones, and for each sample the data presented in figures 1 and 2 is the cool down from 1000 K to 300 K in the HTRS, followed by the cool down from 300 K to 1.8 K or 0.4 K in the PPMS. A certain level of mismatch still remains at room temperature, both for the absolute value of the measured resistance and for its slope, but these differences are considered small enough to be neglected in the comparative analysis we present in this work.

For the crystallographic studies, room temperature X-ray diffraction (XRD) data of as-grown and annealed YbNi2B2C single crystals were collected using a Bruker APEX diffractometer with Mo $`K_\alpha `$ radiation ($`\lambda =0.71073`$ Å) and a detector-to-crystal distance of 5.08 cm. Data were collected for the full sphere and were harvested by collecting three sets of frames with 0.3 scans in $`\omega `$ for an exposure time of 10 seconds per frame. The range of $`2\theta `$ extended from 3.0 to 56.0. Data were corrected for Lorentz and polarization effects; absorption corrections were based on fitting a function to the empirical transmission surface as sampled by multiple equivalent reflections. Unit cell parameters were indexed by peaks obtained from 90 frames of reciprocal space images and then refined using all observed diffraction peaks after data integration. Together with systematic absences, the space group $`I4/mmm`$ (N. 139) was selected for subsequent structural analysis. The structure solution was obtained by direct methods and refined by full-matrix least-squares refinement of $`F_o^2`$ using the SHELXTL 6.10 packageSheldrick (2000).

For the transmission electron microscopy (TEM) evaluations, the as-grown and annealed single crystals of YbNi2B2C were characterized using a Philips CM30 transmission electron microscope operating at 300 keV. The single crystals were placed in a mortar with denatured ethyl alcohol and ground to a fine powder. A pipette was used to place a small drop of the suspension onto a holey carbon grid. The suspension was air-dried before insertion into the microscope. The crystalline structure of fragments which were electron transparent was evaluated using selected area electron diffraction pattern (SADP). The defect density was measured on select areas of the bright-field images where the thickness was relatively uniform.

## III Resistivity Measurements

The main motivation for extending the resistivity measurements to high temperatures was based on the assumption that, if the excess scattering seen at room temperature and below in unannealed samples arises from the contribution of Yb sites with Kondo temperatures extending up through room temperature, there should be a temperature $`T_K^{max}`$ above which the samples no longer display anomalous scattering, and the difference between the resistivity of unannealed sample and the annealed sample becomes a “simple” Matthiessen impurity term.

The temperature dependence of the electrical resistivity for three unannealed pieces of YbNi2B2C cut from the same crystal is shown in Fig. 1. The general behavior below room temperature is the same as has been previously reported for as-grown samples of this compound,Miao et al. (2002a); Avila et al. (2002); Bud’ko et al. (1997); Hossain et al. (1997); Yatskar et al. (1999); Rathnayaka et al. (1997) i.e., a high-scattering metallic behavior down to about 30 K, below which the resistivity decreases rapidly. If normalized to their respective resistance values at any given temperature, these three curves essentially collapse into a single one in the entire measured temperature interval, meaning that the selected crystal was quite homogeneous and the differences seen in the three measurements result from the overall $`\pm 10`$% uncertainty in estimating the geometrical factor $`A/d`$ for each piece. The residual resistivity ratio (defined as $`R(300)/R(1.8)`$) is $`RRR=10`$ for all three samples, somewhat higher than $`RRR=67`$ found for several other unannealed crystals grown similarly Miao et al. (2002a); Avila et al. (2002), indicating that by measuring the high temperature region first we are already annealing out a small portion of the disorder, but still maintaining most of the characteristic unannealed sample behavior.

A crystal of LuNi2B2C was also measured to provide a reference for the non-magnetic (e.g. standard electron-phonon) contribution to scattering down to 16 K, below which it becomes superconducting. The high temperature region is somewhat noisier for this sample because its cross-section is about 3$`\times `$ larger than the others and therefore the actual measured resistance is proportionally smaller. On the other hand, the uncertainty in this sample’s geometrical factor is equally smaller, so its slope becomes a good reference. It is interesting to note that, assuming $`\rho _0l=412\mu \mathrm{\Omega }`$ cm Å for LuNi2B2C,Cheon et al. (1998) near 1000 K the measured resistivity of order $`100\mu \mathrm{\Omega }`$ cm implies a mean free path $`l4`$ Å, already comparable to the unit cell lattice parameter $`a`$ (the Ioffe-Regel criterion), however the resistivity behavior is still quite linear with a slope $`\mathrm{\Delta }\rho /\mathrm{\Delta }T=0.1\mu \mathrm{\Omega }`$ cm/K between 600 K and 900 K. Above 600 K, all three YbNi2B2C samples change slope and roughly follow that of the LuNi2B2C sample (except for the nearly constant extra magnetic contribution).

Samples 1 and 3 were then annealed for 150 hours at 950C and remeasured. Figures 2a and 2b present the comparison between unannealed and annealed resistivity for each sample. Let us go through the comparison in detail for sample 1. At high temperatures (above about 600 K), the curves for the unannealed and annealed conditions of the sample run essentially parallel to each other (and roughly parallel to LuNi2B2C), meaning that in this region there is simply a temperature-independent resistivity factor $`\rho _0`$, in accordance with Matthiessen’s rule. As we cool below 600 K the unannealed curve maintains a more slowly decreasing resistivity (or, conversely, a $`relatively`$ increased level of scattering), when compared to the other two which continue to run parallel down to about 200 K. Around this temperature the annealed YbNi2B2C curve presents a broad shoulder, most likely related to the thermal depopulation of the CEF multiplet with $`T_{cef}100`$ K, and therefore approaches the LuNi2B2C curve. Between 50 K and 10 K the resistivity of the annealed sample resembles a Kondo-minimum type behavior, and below 10 K its resistivity once again drops rapidly due to the onset of coherent scattering of the Yb ions, making its final approach towards the conventional scattering level of the normal-state, non-magnetic LuNi2B2C. Qualitatively similar results are found for sample 3. For both YbNi2B2C samples $`RRR`$ increased to 17 with annealing, consistent with the previous results from annealed crystalsAvila et al. (2002).

Following the encouraging results above, we remeasured these samples in the PPMS with the 3He option in order to extend the resistivity curves down to our lowest measurable temperatures, so as to check how annealing affects the behavior in this regime and verify whether or not they follow the theoretical predictions for a Fermi-liquid.

In figure 3a we show the temperature dependence of resistivity below 30 K of sample 2 which was left unannealed, and samples 1 and 3 (annealed). The decrease in resistivity below 10 K is quite evident in the annealed samples, whereas the resistivity of the unannealed one is already decreasing rapidly below 30 K. In figure 3b we plot the same data in the region below 2 K as a function of $`T^2`$, and all three curves are slightly sub-linear in this plot, indicating that they cannot be completely described by a Fermi-liquid type relation $`\rho (T)=\rho _0+AT^2`$ in this temperature range. In the graph we have included a fit of this expression to the three datasets below 1 K. It seems that with the decrease in disorder there is a decrease in both $`\rho _0`$ and $`A`$. The only other investigations of this type previously reported on an unannealed crystalYatskar et al. (1997, 1999) resulted in $`\rho (T)=12+1.2T^2`$ $`\mu \mathrm{\Omega }`$ cm, and are consistent with this trend. Such a large reduction in $`A`$ is not expected within the framework of the Fermi-liquid models, since our previous studyAvila et al. (2002) found that the electronic specific heat coefficient $`\gamma `$ changes very little with annealing (no more than 10% for the lowest measured temperatures). Even taking into account this small change in $`\gamma `$, the ratio $`A/\gamma _0^2`$ still decreases from $`0.4\times 10^5`$ to $`0.3\times 10^5\mu \mathrm{\Omega }`$ cm (mol K/mJ)2 with annealing. These values are smaller than the Kadowaki-Woods ratioKadowaki and Woods (1986) of $`1\times 10^5\mu \mathrm{\Omega }`$ cm (mol K/mJ)2 found empirically for many Fermi-liquid systems.

The difference of 2-3 $`\mu \mathrm{\Omega }`$ cm in $`\rho _0`$ between our annealed and unannealed YbNi2B2C samples is similar to the results obtained by annealing studies on non-hybridizing TmNi2B2CMiao et al. (2002b) and also comparable to the results of Lu(Ni1-xCox)2B2C substitution studiesCheon et al. (1998) when $`x1\%`$. These similarities suggest that the optimal annealing of single crystals grown by the Ni2B flux growth method is essentially removing lattice imperfections on the order of 1%, and in the coherent scattering regime YbNi2B2C behaves like all other members of the family. Furthermore, the difference in resistivity between annealed and unannealed samples above 600 K (Fig. 2) was found to be 3-6 $`\mu \mathrm{\Omega }`$ cm, which is remarkably similar as well (considering that the separations between the resistivity curves are strongly affected by geometrical uncertainties at these high temperatures) and therefore once again points to the return of a more conventional scattering regime at very high temperatures.

Given the consistency of the extended transport measurements with the hypothesis that ligandal disorder in the as-grown samples is leading to redistributions of Kondo temperatures for a small fraction of the Yb sites, the question that arises naturally is: what are the types of disorder contributing to this behavior? To address this question, we now present a series of structural and microscopic experiments conducted on as-grown and annealed single crystals.

## IV Single Crystal X-Ray Refinements

The single-crystal XRD technique is useful to check for vacancies, substitutions, superstructures, and other perturbations of periodic atomic sites in the crystal. The structure refinement for the as-grown and annealed YbNi2B2C was carried out under the assumption that no structural disorder was present. By this assumption all systematic absences pointed to the expectedSiegrist et al. (1994a) crystal system (Tetragonal) and space group ($`I4/mmm`$) represented in figure 4. Atomic positions were assigned based on site symmetry and bond length and agree with previously published single crystal data for other $`R`$Ni2B2C series membersSiegrist et al. (1994b). However, subtle problems with the initial structure refinement of the as-grown samples implied that the solution was not complete.

For most crystal structure refinements a parameter is added to account for the accordance of extinction. In the SHELXTL software the extinction parameter accounts for both primary and secondary extinctionSheldrick (2000). Before the extinction parameter was added to the refinement of the as-grown YbNi2B2C the R-index was 0.1224 and the thermal parameters for Yb and B were not convergent (table 1). After the extinction parameter was added the refinement appeared to be stable with the exception of the C thermal parameter (this is typical for light elements in the presence of heavy elements). However, the extinction parameter itself (0.24(2)) is larger than what is normally expected (table 2). In solid state structures a large extinction parameter usually indicates the existence of a super-cell structure. Because super-cell structural solutions were not found in this case, as indicated by the absence of additional reflections, the large extinction parameter is most likely due to the presence of dislocations. Refinement of the annealed YbNi2B2C structure presented no problems (extinction parameter 0.065(8)), suggesting that much of the disorder present in the sample before the annealing process has been removed. There is no evidence of vacancies or substitutions within the resolution limit of this technique.

## V Transmission Electron Microscopy

In order to obtain a more detailed and quantitative estimate of these dislocations we studied the TEM patterns of crushed crystals. Figure 5(a) and (b) show two representative bright field TEM images and the corresponding SADPs (insets) of the crushed as-grown and annealed YbNi2B2C single grains, respectively. Both SADPs are fully consistent with the known space group ($`I4/mmm`$ (139)) for this compound, and they do have a common reciprocal lattice vector, (110), which is approximately equal for these two samples. There are no obvious second phases and the diffraction spots are sharp, indicating well crystallized samples. However, the as-grown crystallites show a surprisingly large number of dislocations (Fig. 5a). In contrast, the annealed samples while exhibiting some regions of moderate defect density, the defect density is qualitatively lower than the as-grown sample. Using an estimated thickness of 800 and 500 nm for the as-grown and annealed samples, we estimate a defect density of $`1.5\times 10^8`$ and $`<6\times 10^7`$ cm-2 for the as-grown and annealed samples, respectively.

Although it is not expected that room temperature crushing of a nominally brittle intermetallic would result in extensive dislocation formation, it is not without precedenceBeraha et al. (1997). The difference in the dislocation density between the as-grown and annealed samples may be due to the presence of Frank-Read sources, small defects which when subjected to high shear stresses are capable of overcoming the Peierls forces resisting dislocation mobility. This would imply that the as-grown crystals may have many small, and hard to detect defects which are eliminated by annealing.

## VI Discussion

In our previous studyAvila et al. (2002) we made a semi-quantitative evaluation of the fraction of Yb sites that would need to be perturbed in order to account for the excess resistivity seen at room temperature and below in the unannealed samples. Using the expression for the increase in resistivity caused by a Kondo impurity for temperatures below the Kondo temperature associated with the impurityAvila et al. (2002); Thompson and Fisk (1985):

$`\mathrm{\Delta }\rho _{max}=\frac{\mathrm{}}{e^2}\frac{4\pi c}{pk_F}(2l+1)`$

where $`c`$ is the concentration of Yb “impurities”, $`k_F`$ is the Fermi momentum, $`p`$ is the number of electrons per atom, and $`l`$ is the ytterbium $`4f`$ angular momentum, we estimated an increase in the room temperature resistivity of about $`30\mu \mathrm{\Omega }cm`$ for every 1% concentration of affected Yb sites with $`T_K`$ greater than 300 K. The extended resistivity measurements and analysis performed in the present work have confirmed and reinforced that initially proposed semi-quantitative model. We have observed that when the sample is well into the coherent scattering regime ($`TT_K10`$ K) or well above the maximum perturbed Kondo temperature ($`TT_K^{max}600`$ K) the Yb ions behave like the more conventional, non-hybridizing rare-earth members in the $`R`$Ni2B2C series and the overall resistivity decrease with annealing is of only 2-3 $`\mu \mathrm{\Omega }cm`$, once again consistent with a density of order 1% of conventional scattering sites.

The TEM experiments clearly showed the presence of defects in both unannealed and annealed samples most likely associated with lattice dislocations, and the defect density was estimated to be at least 2.5 times larger in the unannealed samples. Similar qualitative changes were seen in the XRD data. With these results we can propose a scenario where the dislocation sites act as a source of strain fields which disturb nearby Yb sites with intensity inversely proportional to their distance to the defect origin, which would then lead to the appearance of a distribution of local Kondo temperatures along the strain field decay. Estimating the extension of a strain field around a dislocation core is complicated, and will vary according to the nature of the dislocation, the modulii of the material and the length of the Burger’s vector. However, we can once again attempt a semi-quantitative analysis. Using the defect density of $`1.5\times 10^8`$ cm-2 for the unannealed sample and assuming that 1% of the Yb sites are subject to a strain field, we can estimate that all Yb sites within a radius of $`460`$ Å of each dislocation core would be feeling the effect of the strain field. This value corresponds to about 130 unit cells along the plane and 40 unit cells along the $`c`$-axis, and seems rather large given that, according to Eshelby and othersHirsch (1975), the region of elastic limit should be on the order of 5 times the Burger’s vector. For our samples the dislocations appear to be of the $`1/2110\{001\}`$ type with corresponding Burger’s vector $`b=a/2[110]2.47`$ Å. We can still argue an enhanced sensitivity of the hybridized Yb ions to the strain, but even so it is likely that the extension of the strain field (or the fraction of affected Yb sites) is being overestimated.

Whereas there are very clear changes in the microstructural data (TEM and single-crystal XRD) with annealing, the link between these changes and the changes in the transport measurements remains only semi-quantitative.

## VII Conclusion

In this work we have fully characterized the resistivity behavior of annealed and unannealed YbNi2B2C single crystals in the range of $`0.4<T<1000`$ K. Whereas the variability of resistivity behavior for compounds with hybridizing moments has been long suspected, our studies on YbNi2B2C finally allow a clear and controlled demonstration of this effect. We were also able to use single crystal XRD and TEM data to gain confidence in the claim that the changes are indeed intrinsic to the compound which is clearly single phase in both annealed and as-grown conditions. For this particular material and growth process, lattice dislocations seem to be the dominant defect type found in as-grown crystals, and is most likely being responsible for environment changes in the nearby Yb3+ ions which lead to deviations of local Kondo temperatures from their intrinsic value of $`10`$ K. As a consequence, it is important to realize that any attempt at detailed analysis of transport properties such as those found for YbNi2B2C (and these are not uncommon) are highly suspect, and any comparison of experiments with theoretical models must either take into account the sensitivity of the hybridizing moments on the modelling side or guarantee a nearly defect-free sample, a highly non-trivial experimental requirement.

###### Acknowledgements.

We are thankful to F. Borsa for urging us to measure the high-temperature of YbNi2B2C, to J. Schmalian for stimulating discussions, to I. R. Fisher and T. Wiener for the development of the high-temperature resistivity system, and to Y. Mozharivskyj for assistance with the x-ray refinements. Ames Laboratory is operated for the US Department of Energy by Iowa State University under Contract No. W-7405-Eng-82. This work was supported by the Director for Energy Research, Office of Basic Energy Sciences.


text:
# Momentum-position realization of the Einstein-Podolsky-Rosen paradox

## Abstract

We report on a momentum-position realization of the EPR paradox using direct detection in the near and far fields of the photons emitted by collinear type-II phase-matched parametric downconversion. Using this approach we achieved a measured two-photon momentum-position variance product of $`0.01\mathrm{}^2`$, which dramatically violates the bounds for the EPR and separability criteria.

In 1935 Einstein, Podolsky and Rosen (EPR) EPR  wrote one of the most controversial and influential papers of the twentieth century. They proposed a *gedanken* experiment involving two particles entangled simultaneously over a continuum of position and momentum states. By measuring either the position or momentum of one of the particles, the position or momentum of the other could be inferred with complete certainty. Under the assumptions of EPR, the ability to make such an inference meant that the position and momentum of the unmeasured particle were simultaneous realities, in violation of Heisenberg’s uncertainty relation. This thought experiment became known as the EPR paradox. In 1951 Bohm Bohm  cast the EPR paradox into a simpler, discrete form involving spin entanglement of two spin-1/2 particles, such as those produced in the dissociation of a diatomic molecule of zero spin. From Bohm’s analysis sprang Bell’s inequalities Bell ; CHSH  and much of what is now the field of discrete quantum information Aspect81 ; Ekert91 ; Naik00 ; Jennewein00 ; Tittel00 ; Bouwmeester97 .

In recent years, however, there has been a movement toward the study of entanglement of continuous variables as originally discussed by EPR Franson89 ; Kwiat90 ; Rarity90 ; Reid88 ; Ou92 ; Duan00 ; Simon00 ; Julsgaard01 ; Silberhorn01 ; Mancini02 ; Bowen02 ; Korolkova02 ; Braunstein98 ; Gatti03 . Of particular interest was the early work of Reid and Drummond Reid88 . They derived an EPR criterion and showed how it could be implemented with momentum-like and position-like quadrature observables of squeezed light fields. Shortly thereafter, the experiment was realized by Ou et al. Ou92 . Later Duan et al. Duan00  and Simon Simon00  derived necessary and sufficient conditions for the inseparability (entanglement) of continuous variable states. A flurry of experimental activity ensued in both atomic ensembles Julsgaard01  and squeezed light fields Silberhorn01 ; Bowen02 ; Korolkova02 .

Here we report on a demonstration of the EPR paradox using position- and momentum-entangled photon pairs produced by spontaneous parametric downconversion. We find that the position and momentum correlations are strong enough to allow the position or momentum of a photon to be inferred from that of its partner with a product of variances $`0.01\mathrm{}^2`$, which violates the separability bound by two orders of magnitude.

In the idealized entangled state proposed by EPR,

$$|\mathrm{EPR}_{\mathrm{}}^{\mathrm{}}|x,x𝑑x=_{\mathrm{}}^{\mathrm{}}|p,p𝑑p,$$
(1)
the positions and momenta of the two particles are perfectly correlated. This state is non-normalizable and cannot be realized in the laboratory. However, the state of the light produced in parametric downconversion can be made to approximate the EPR state under suitable conditions. In parametric downconversion, a pump photon is absorbed by a nonlinear medium and re-emitted as two photons (conventionally called signal and idler photons), each with approximately half the energy of the pump photon. Considering only the transverse components, the momentum conservation of the downconversion process requires $`𝐩_1+𝐩_2=𝐩_p`$ where $`1,2,p`$ refer to the signal, idler, and pump photons. Provided the uncertainty in the pump transverse momentum is small, the transverse momenta of the signal and idler photons are highly anti-correlated. The exact degree of correlation depends on the structure of the signal+idler state. In the regime of weak generation, this state has the form

$$|\psi _{1,2}=|\mathrm{vac}+A(𝐩_1,𝐩_2)|𝐩_1,𝐩_2𝑑𝐩_1𝑑𝐩_2$$
(2)
where $`|\mathrm{vac}`$ denotes the vacuum state and the two-photon amplitude $`A`$ is

$$A(𝐩_1,𝐩_2)=\chi E_p(𝐩_1+𝐩_2)\frac{\mathrm{exp}(i\mathrm{\Delta }k_zL)1}{i\mathrm{\Delta }k_z}.$$
(3)
Here $`\chi `$ is the coefficient of the nonlinear interaction, $`E_p`$ is the amplitude of the plane-wave component of the pump with transverse momentum $`𝐩_1+𝐩_2`$, $`L`$ is the length of the nonlinear medium, and $`\mathrm{\Delta }k_z=k_{p,z}k_{1,z}k_{2,z}`$ (where $`𝐤=𝐩/\mathrm{}`$) is the longitudinal wavevector mismatch, which generally increases with transverse momentum and limits the angular spread of signal and idler photons. The vacuum component of the state makes no contribution to photon counting measurements and may be ignored. Also, there is no inherent difference between different transverse components; so without loss of generality, we consider scalar position and momentum. The narrower the angular spectrum of the pump field and the wider the angular spectrum of the generated light, the more closely the integral (2) approximates $`_{\mathrm{}}^{\mathrm{}}\delta (p_1+p_2)|p_1,p_2𝑑p_1𝑑p_2=|\mathrm{EPR}`$ and the stronger the correlations in position and momentum become.

The experimental setup used to determine position and momentum correlations is portrayed in Fig. 1. The idea is to measure the positions and momenta by measuring the downconverted photons in the near and far fields respectively Gatti03 . The source of entangled photons is spontaneous parametric downconversion generated by pumping a 2mm thick Type-II BBO crystal with a 30 mW, cw, 390 nm laser beam. A prism separates the pump light from the downconverted light. The signal and idler photons have orthogonal polarizations and are separated by a polarizing beamsplitter. In each arm, the light passes through a narrow 40 $`\mu `$m vertical slit, a 10nm spectral filter, and a microscope objective. The objective focuses the transmitted light onto a multimode fiber which is coupled to an avalanche photodiode single-photon counting module. The spectral filter ensures that only photons with nearly equal energies are detected. To measure correlations in the positions of the photons, a lens of focal length 100 mm (placed prior to the beamsplitter) is used to image the exit face of the crystal onto the planes of the two slits (Fig. 1a). One slit is fixed at the location of peak signal intensity. The other slit is mounted on a translation stage. The photon coincidence rate is then recorded as a function of the displacement of the second slit. To measure correlations in the transverse momenta of the photons, the imaging lens is replaced by two lenses of focal length 100 mm, one in each arm, a distance $`f`$ from the planes of the two slits (Fig. 1b). These lenses map transverse momenta to transverse positions, such that a photon with transverse momentum $`\mathrm{}k_{}`$ comes to a focus at the point $`x=fk_{}/k`$ in the plane of the slit. Again, one slit is fixed at the location of peak count rate while the other is translated to obtain the coincidence distribution.

By normalizing the coincidence distributions, we obtain the conditional probability density functions $`P(x_1|x_2)`$ and $`P(p_1|p_2)`$ (Fig. 2). These probability densities are then used to calculate the uncertainty in the inferred position or momentum of photon 1 given the position or momentum of photon 2:

$`\mathrm{\Delta }x_{\mathrm{inf}}^2`$ $`={\displaystyle (x_1x_2)^2P(x_1|x_2)𝑑x_1}`$
$`\left({\displaystyle (x_1x_2)P(x_1|x_2)𝑑x_1}\right)^2`$ (4)
$`\mathrm{\Delta }p_{\mathrm{inf}}^2`$ $`={\displaystyle (p_1+p_2)^2P(p_1|p_2)𝑑p_1}`$
$`\left({\displaystyle (p_1+p_2)P(p_1|p_2)𝑑p_1}\right)^2.`$ (5)
Because of the finite width of the slits, the raw data in Fig. 2 describe a slightly broader distribution than is associated with the downconversion process itself. By adjusting the computed values of $`\mathrm{\Delta }x_{\mathrm{inf}}`$ and $`\mathrm{\Delta }p_{\mathrm{inf}}`$ to account for this broadening (an adjustment smaller than 10%), we obtain the correlation uncertainties $`\mathrm{\Delta }x_{\mathrm{inf}}=0.027\mathrm{mm}`$ and $`\mathrm{\Delta }p_{\mathrm{inf}}=3.7\mathrm{}\mathrm{mm}^1`$. The measured variance product of the inferred state is

$$\mathrm{\Delta }x_{\mathrm{inf}}^2\mathrm{\Delta }p_{\mathrm{inf}}^2=0.01\mathrm{}^2.$$
(6)
Also shown in Fig. 2 are the predicted probability densities. These curves contain no free parameters and are obtained directly from the two-photon amplitude $`A(p_1,p_2)`$, which is determined by the optical properties of BBO and the measured profile of the pump beam. Fig. 2 indicates that the correlation widths we obtained are intrinsic to the downconversion process and are limited only by the degree to which it deviates from the idealized EPR state (1). The value of $`\mathrm{\Delta }p_{\mathrm{inf}}`$ is limited by the finite width of the pump beam. The pump photons in a Gaussian beam of width $`w`$ have an uncertainty $`\mathrm{}/2w`$ in transverse momentum which, due to conservation of momentum, is imparted to the total momentum $`p_1+p_2`$ of the signal and idler photons. The value of $`\mathrm{\Delta }x_{\mathrm{inf}}`$ is limited by the range of angles over which the crystal generates signal and idler photons. If the angular width of emission is $`\mathrm{\Delta }\varphi `$, then the principle of diffraction indicates that the photons cannot have a smaller transverse dimension than $`(k_{s,i}\mathrm{\Delta }\varphi )^1`$. Careful analysis based on the angular distribution of emission yields $`\mathrm{\Delta }x_{\mathrm{inf}}=1.88(k_{s,i}\mathrm{\Delta }\varphi )^1`$. With the measured beam width of $`w=0.17\mathrm{mm}`$ and predicted angular width $`0.012\mathrm{rad}`$, the theory predicts $`\mathrm{\Delta }x_{\mathrm{inf}}^2\mathrm{\Delta }p_{\mathrm{inf}}^2=0.0036\mathrm{}^2.`$ This is somewhat smaller than the experimentally calculated value of $`0.01\mathrm{}^2`$, even though the data appears to closely match the theoretical curves. The reason for this discrepancy is that the experimental distributions have small ($`1`$% of the peak) but very broad wings. The origin of these coincidence counts is unknown; they are perhaps due to scattering from optical components. If these counts are treated as a noise background and subtracted, the experimentally obtained uncertainties come into agreement with the theoretically predicted values, yielding an uncertainty product of $`0.004\mathrm{}^2`$.

To interpret these results, it is helpful to consider the relationship between the original EPR paradox and the issues of entanglement, non-locality, and quantum signatures, which have been the subjects of more modern studies. Because of lingering doubts about quantum theory, raised in part by the EPR paradox, there has been much interest in “quantum signatures”, i.e. phenomena that confirm the predictions of quantum mechanics and cannot be explained by classical mechanics. However, the intent of EPR was not to reveal a discrepancy between classical and quantum theory, but to reveal an inconsistency, or incompleteness, within quantum theory. This inconsistency was revealed by showing that measurements of one particle could be used (seemingly) to infer the state of an unmeasured particle with greater certainty than is allowed by the uncertainty relation, i.e. $`\mathrm{\Delta }x_{\mathrm{inf}}^2\mathrm{\Delta }p_{\mathrm{inf}}^2\mathrm{}^2/4`$ EPR ; Reid88 . A key component of EPR’s argument was the assumption of “local reality.” Under this assumption, the statistics of each particle depend only on a (hidden) state parameter $`s`$ which is determined while the particles are close enough to interact. This assumption is commonly taken to mean that the joint probablity of any pair of observables $`a_1`$ and $`b_2`$ must be expressible as

$$P(a_1,b_2)=\underset{s}{}P_sP(a_1|s)P(b_2|s)$$
(7)
where $`P_s`$ is the probability of state $`s`$ and $`P(a_1|s)`$ and $`P(b_2|s)`$ are the conditional probabilities for particle 1 and particle 2 respectively. This constraint led to the Bell tests of local realism Bell ; CHSH  and their subsequent experimental realizations (e.g., Aspect81 ). Another consequence of EPR  has been the development of the concept of entanglement, which is closely connected to non-locality. A bipartite system is entangled (inseparable) if its density matrix *cannot* be written in the form

$$\rho =\underset{s}{}P_s\rho _{1,s}\rho _{2,s}.$$
(8)
Measurement of either member of an entangled system projects both members into a mixture of states consistent with the result of the measurement. It is now generally accepted that this mutual projection occurs even when the particles are widely separated; hence, within the framework of quantum theory, entanglement and non-local behavior have a one-to-one relationship. This conclusion is confirmed by the fact that any system which satisfies eqn. (8) also satisfies eqn. (7), and vice versa. However, entanglement does not always rule out local realism in the context of a theory other than quantum mechanics. For example, the strong position and momentum correlations of the entangled EPR state can be readily explained by local-realistic classical mechanics, which does not impose an uncertainty relation on position and momentum.

The close connection between entanglement and non-locality has prompted refinement of EPR’s criterion, resulting in a number of different separability tests and entanglement measures. Of relevance to the original EPR paradox are the tests for separability of continuous- variable systems developed by Duan *et al.* Duan00 , Simon Simon00 , and Mancini *et al.* Mancini02 . The tests in Duan00  and Simon00  involve sums of dimensionless variances and are not scale invariant; hence it is not clear how, or even if, they may be applied to the present EPR experiment, which involves dimensional position and momentum. The criterion derived by Mancini *et al.* Mancini02  is more useful here. It states that separable systems must satisfy the joint uncertainty relation

$$(\mathrm{\Delta }x_{12})^2(\mathrm{\Delta }p_{12})^2\mathrm{}^2$$
(9)
where $`(\mathrm{\Delta }x_{12})^2=(x_1x_2)^2(x_1x_2)^2`$ and $`(\mathrm{\Delta }p_{12})^2=(p_1+p_2)^2(p_1+p_2)^2.`$ The angle brackets denote expected values over the respective joint probability distributions. In our experiment the widths of the conditional probability distributions $`P(x_1|x_2)`$ and $`P(p_1|p_2)`$ are essentially independent of $`x_2`$ and $`p_2`$ over most of their ranges, so that $`\mathrm{\Delta }x_{\mathrm{inf}}^2`$ and $`\mathrm{\Delta }p_{\mathrm{inf}}^2`$ are nearly equal to $`(\mathrm{\Delta }x_{12})^2`$ and $`(\mathrm{\Delta }p_{12})^2`$. Hence our results constitute a two order-of-magnitude violation of the joint uncertainty relation (9), which in this case is both a separability criterion and a local-realism criterion.

Finally, we note that the EPR paradox does not represent a true inconsistency. It is generally accepted that the EPR argument fails because the assumption of local realism is invalid: The position or momentum of the unmeasured particle becomes a reality when, and only when, the corresponding quantity of the other particle is measured. As the measurement involves only one quantity or the other, the position and momentum of the unmeasured particle need not be simultaneous realities.

In conclusion, we have reported the experimental realization of Einstein, Podolsky and Rosen’s paradox using momentum-position entangled photons. We have measured position and momentum correlations resulting in a variance product which dramatically violates the original EPR criterion and a modern inseparability criterion. Compared to squeezed-light realizations of the EPR paradox, the momentum-position realization has several attractive features which make it promising for further development. For one, the entanglement is observed using direct photon detection, which is experimentally simpler than homodyne detection. Secondly, the entanglement does not reside in the photon count, which frees this quantity to be used for postselection. Since the position and momentum measurements involve only those photons that are detected, the measured entanglement is not degraded by optical loss which inevitably occurs in real systems. For both of these reasons, systems with very small values of the joint uncertainty product can be readily achieved in practice. This capability has already been used to achieve near- diffractionless coincidence imaging Bennink . We believe that the work presented here sets the stage for many more interesting applications to come.

JCH acknowledges support from the NSF and University of Rochester. RWB gratefully acknowlwedges support by ARO under award DAAD19-01-1-0623 and ONR under award N00014-02-1-0797

Sean Bentley is now at Adelphi University.


text:
# Multistability, noise and attractor-hopping: The crucial role of chaotic saddles

## Abstract

We investigate the hopping dynamics between different attractors in a multistable system under the influence of noise. Using symbolic dynamics we find a sudden increase of dynamical entropies, when a system parameter is varied. This effect is explained by a novel bifurcation involving two chaotic saddles. We also demonstrate that the transient lifetimes on the saddle obey a scaling law in analogy to crisis.

PACS number 05.45+b

Systems with a large number of coexisting stable states have been the subject of increasing interest recently. This multistable behavior occurs in many different fields like optics Wiesenfeld:1990 , chemistry Marmillot:1991 , neuroscience schiff:1994 , semiconductor physics Prengel:1994 , plasma physics Chern:1991  and coupled oscillators kim:1997b . Moreover, addition of noise to multistable systems has led to several interesting phenomena, like noise-induced preference of attractors Kaneko:1997 ; Kraut:1999 , directed diffusion hondou:1995 , noise-enhanced multistability kim:1997a  and chaotic itinerancy ikeda:1989 . The latter effect, which consists of hopping between different attractors caused by the noise, has also been observed experimentally in an optical system Arecchi:1990 .
In this Letter we focus on the dynamics of the attractor-hopping process by taking a simple model as a paradigm. Due to the noise the attractors become ’metastable’ and the trajectory starts hopping between the different attractors. In the case of fractal basin boundaries, this hopping process consists of two steps. In the first one the trajectory leaves the open neighborhood about the attractor according to Arrhenius law grassberger:1989 , in the second one, the trajectory bounces around on the chaotic saddle before cascading again into the open neighborhood of an attractor. The general mechanism for the first step of the hopping process is the same for bi- and multistable systems. In both cases we find Arrhenius law; quanitative differences are due to different relative sizes of basins of attraction Kraut:1999  However, for the attractor-hopping dynamics, our study points out a major difference between bistable and multistable systems: While in bistable systems the structure of the saddle separating the two attractors (saddle point or chaotic saddle) does not play a role for the hopping characteristics, it is essential for the attractor-hopping in a system possessing a multitude of attractors. The structure of the chaotic saddles influences the nature of the hopping process, in particular it determines which transitions between attractors are possible. Bifurcations in the chaotic saddles lead to changes in the ’accessibility’ of attractors in the hopping process.
Employing symbolic dynamics, we assign one symbol to each attractor, thus transforming our hopping time series onto a symbolic string. As a next step we compute both, the Shannon and the topological entropy, which can be regarded as measures of complexity for the hopping dynamics. As a system parameter is varied the complexity of the hopping process changes beyond a certain threshold value. We show, that this change can be related to a novel bifurcation, namely a merging of two chaotic saddles accompanied by the emergence of additional points filling the gap between the formely separated saddles. This bifurcation is mediated by a snapback repellor marotto:1978 , whose eigenvalues determine the scaling law for the transient lifetimes on one chaotic saddle close to the bifurcation point.
Our basic prototype model of multistablity, which captures the main features of highly multistable systems, namely more than two final states (attractors) and a complexly interwoven fractal basin boundary separating these states, is given by the 10-fold iterate of two coupled logistic maps:

$`x_{k+1}`$ $`=`$ $`1.0\alpha x_k^2+\gamma (y_kx_k)`$
$`y_{k+1}`$ $`=`$ $`1.0\alpha y_k^2+\gamma (x_ky_k).`$ (1)
We fix $`\gamma `$, the coupling strength, at $`\gamma =0.29`$, and vary $`\alpha `$, the nonlinearity, in the range $`\alpha [0.72,0.755]`$. For this parameter set, there exists a stable period $`10`$-orbit for the two coupled logistic maps themselves. Since we consider the 10-fold iterate, our map exhibits 10 coexisting fixed point attractors, where each 5 lie above and below the symmetry axis $`x=y`$, respectively. Using only initial conditions below the symmetry axis, we can restrict our study to a system possessing 5 coexisting fixed point attractors, 5 repellors with two unstable directions and 10 saddle points. Each of the 5 attractors is surrounded by an open neighborhood of different size, in which all initial conditions converge to the corresponding attractor inside this neighborhood. The basins of attraction as a whole have a complex fractal structure due to a homoclinic bifurcation which occurs already for smaller values of the nonlinearity $`\alpha `$. As a next step we apply noise to the system, using Gaussian, white noise with standarddeviation $`\sigma `$ added to the $`x`$ and $`y`$ component of the 10-fold iterate of (Multistability, noise and attractor-hopping:  The crucial role of chaotic saddles). Small noise causes the system now to alternate between the different states, where long periods close to a fixed point, comparable to laminar motion, are interrupted by short, sudden bursts, which are reminiscent to intermittent behavior. During the bursts the motion takes place on the chaotic saddles separating the attractors, until the trajectory is again being injected in the neighborhood of one of the $`5`$ fixed points. This complex dynamics is depicted in Fig. 1.

According to our aim we are now interested in the characteristic properties of the hopping dynamics. Thus we address the question, how the complexity of the hopping dynamics changes with the nonlinearity parameter $`\alpha `$, namely whether every fixed point has a positive transition probability to every other fixed point, via a transient on a chaotic saddle. For this purpose, we employ the concept of symbolic dynamics. We use an encoding scheme in which we assign one symbol to each attractor poon:1995 : Neglecting the number of iterations the trajectory spends close to a fixed point, a symbol is given for every fixed point and only after a jump out of an attractor a new symbol is bestowed, according to the attractor where the trajectory lands at. Using this scheme, we focus only on the structural properties of the jumps, not taking fully into account the complete temporal evolution in each iteration. As a quantitative measure of the complexity of the symbol string we use the Shannon entropy, in analogy to the Kolmogorov-Sinai entropy kolmogorov:1958 , given by

$`h_S`$ $`=`$ $`\underset{n\mathrm{}}{lim}{\displaystyle \frac{H_n}{n}}=\underset{n\mathrm{}}{lim}(H_{n+1}H_n)`$ (2)
$`=`$ $`\underset{n\mathrm{}}{lim}{\displaystyle \frac{1}{n}}\left({\displaystyle \underset{|S|=n}{}}p(S)\mathrm{log}p(S)\right),`$
where $`S=s_1s_2\mathrm{}s_n`$ denotes a finite symbol sequence consisting of $`n`$ elements $`s_i=1,2,\mathrm{},5`$, $`p(S)`$ its probability of occurrence and $`H_n`$ the block entropy of block length $`n`$. Numerically, the quantity $`lim_n\mathrm{}(H_{n+1}H_n)`$ converges already for $`n=1`$, indicating that the hopping dynamics is a Markov process of first order. This corresponds to the intuitive expectation, that the dynamics possesses only memory of the last state it was dwelling on. Higher correlations are suppressed by the long laminar-like motion. Consequently a transition matrix between the different fixed points can be constructed, where the entries are the probabilities for a transition between two states.
Another interesting quantity to investigate is the topological entropy. One way of computing it is to built the Stefan transition matrix derrida:1978 , where the entries of the $`5\times 5`$ matrix are $`1`$ and $`0`$, depending on whether a transition took place or not, respectively. The logarithm of the largest eigenvalue of this matrix yields the topological entropy. In our numerical implementation, we use a certain, very small cutoff limit of $`0.001`$ for the transition probabilities, below which we regard the value as zero. Thus we neglect transition, having extremely small probabilities compared to the other transitions. However, the exact value of this cutoff does not change the results significantly, as long as it is small enough. The entropies shown in Fig. 2 remain unchanged for cutoff limits between $`10^3`$ to $`10^6`$. It is important to note that, although the escape times out of the stable states depend on the noise level, the transition probabilities remain constant for a wide range of noise values. Thus the considered entropies are not affected by the noise amplitude for a rather large interval of noise levels. In Fig. 2 we present the evolution of these two quantities as the nonlinearity $`\alpha `$ changes. The curve for the topological entropy appears to be a monotonic curve with a devil’s staircase like behavior biham:1989 .

Every time it increases, at least one new transition between hitherto unconnected fixed points is created. Below the homoclinic bifurcation $`\alpha <0.72`$, where the basin boundaries are smooth, only one transition per fixed point is possible, as the fixed points are located on a closed curve made up by the unstable manifolds of the saddlepoints. These unstable manifolds impose a direction. Consequently, the topological entropy is zero. Above a value of $`\alpha 0.75`$, all transitions can occur and the maximum of the topological entropy is reached. The Shannon entropy is always less than the topological entropy. Besides the overall increase of both entropies, the plateau between $`0.74<\alpha <0.744`$ suggests, that there is a change in the dynamical behavior beyond it. This change is due to a change in the topological structure of the system and relates to a new bifurcation as we show next.
We compute the chaotic saddle, responsible for the chaotic dynamics, by the PIM-triple (proper interior maximum) method nusse:1989 . The result is depicted in Fig. 3, and 4, where the chaotic saddles are shown for two different values of $`\alpha `$ together with the $`4`$-fold preimages $`(+)`$ of one repellor $`R`$ marked with a $`\mathrm{}`$.

Slightly above the homoclinic bifurcation ($`\alpha =0.725`$) the saddles are very thin, consisting of two separate rings. At a value of $`\alpha =0.755`$, there exists only one large piece of the nonattracting chaotic set (not shown). This results from the saddle merging, which took place at $`\alpha 0.7430`$ (Fig. 4). The two bands of chaotic saddles merge and the repellor, formerly located ouside of the saddles, is now embedded in the large saddle, including some of its preimages. The merging is accompanied by the successive filling of the formerly empty gap between the two saddles. A similar effect has been observed in chin:1992  where it is called ’spilling’, in chaotic scattering lai:1993 , and in the gap filling crisis szabo:1996 . In contrast to the latter this happens smoothly in the saddle merging.

However, the actual bifurcation takes place at a slightly lower value of $`\alpha 0.7428`$ (not shown). At this parameter value some of the preimages of the marked repellor touch its unstable manifold, which constitutes the border of the chaotic saddle. Through this mechanism, a snapback repellor has developed. A fixed point p is called a snapback repellor, if i) all eigenvalues of p have absolute values larger than 1, ii) there exists a $`qW_{loc}^u(p)`$, the unstable manifold of p, such that the M-fold iterate of the map $`^M(q)=p`$ for some positive integer M and iii) $`detD^M(q)0`$marotto:1978 . While conditions i) and iii) are true for all $`\alpha `$ considered here, it is condition ii) that becomes fulfilled during the bifurcation, as the unstable manifold of the repellor touches its preimages. This bifurcation is responsible for the emergence of the additional points in the gaps, (see Fig. 4). At a slightly higher value of $`\alpha 0.7430`$ the preimages touch not only the unstable manifold, constituting the border of the chaotic saddle, but the saddle itself. As the chaotic saddle is an invariant set, the repellor becomes thus embedded in the saddle and the merging of the two rings has taken place. This merging has an important consequence for the dynamics of the attractor-hopping. Not only jumps between neighbors, but arbitrary jumps become, in principle, possible. Above the critical value, the connected pieces of the chaotic saddle enable the trajectory to jump between the formerly separated parts, thereby causing a sharp increase in the number of the possible transitions, as can be seen from the topological entropy above $`\alpha =0.7430`$ in Fig. 2.
Since the bifurcation described above can be considered as a kind of crisis, it is interesting to see, whether the scaling laws for the characteristic transient lifetimes ott:1993  apply in this case as well. We compute the average escape times out of one (still separated) chaotic saddle ring slightly below the bifurcation. In Fig. 5 the data are plotted, resulting in a scaling relation:

$$\tau (\alpha \alpha _c)^ϵ,$$
(3)
where $`\alpha _c=0.7428`$. The best fit to the data yields an exponent $`ϵ=1.86`$. That value can be reproduced by

$`ϵ`$ $`=`$ $`1/2+(\mathrm{ln}(\beta _1)/\mathrm{ln}(\beta _2))`$ (4)
$`ϵ`$ $`=`$ $`1/2+(\mathrm{ln}(2.45)/\mathrm{ln}(1.9))=1.85,`$ (5)
where $`\beta _1`$ and $`\beta _2`$ are the expanding eigenvalues of the mediating repellor.

This coincidence has been confirmed also for other values of the coupling strength $`\gamma =0.285,0.287,0.292`$. It is worth mentioning that in the original derivation of this formula, which applies to the heteroclinic tangency crisis, $`\beta _1`$ and $`\beta _2`$ are the expanding and contracting eigenvalues of the mediating saddle point grebogi:1986 . However, in our case both, $`\beta _1`$ and $`\beta _2`$, are expanding eigenvalues of the mediating repellor. We conjecture this relation to hold for the general case of a crisis caused by a snapback repellor. The reason, that a law for the heteroclinic tangency is supposed to hold in a homoclinic case of a snapback repellor is, that the mechanism for approaching the fixed point is different here. The fixed point is not reached by successive iterations, but suddenly with a jump of noninfinitesimal size from a finite distance.
In conclusion, we have studied a simple paradigmatic map possessing a multitude of coexisting stable states under the influence of small noise. We have shown that the nature of the attractor-hopping process depends, in contrast to the bistable case, crucially on the structure of the chaotic saddles separating the attractors. Bifurcations of the saddles lead to a change in the hopping dynamics which is manifested in a sudden increase of dynamical entropies. In particular we described a novel bifurcation, a merging of two chaotic saddles involving a snapback repellor, which results in a change of the ’accessibility’ of the attractors. Finally, a conjecture for the scaling law of the transient lifetimes on the chaotic saddle has been provided. The hopping dynamics, the bifurcation and the scaling law should be observable in experimental realizations of such systems. We expect this bifurcation also to occur in systems with three or more dimensions, where the mediating saddle has at least two unstable directions.
We acknowledge Y. Maistrenko for valuable discussions and the DFG for financial support.


text:
# Velocity Statistics in the Two-Dimensional Granular Turbulence

## Abstract

We studied the macroscopic statistical properties on the freely evolving quasi-elastic hard disk (granular) system by performing a large-scale (up to a few million particles) event-driven molecular dynamics systematically and found that remarkably analogous to an enstrophy cascade process in the decaying two-dimensional fluid turbulence. There are four typical stages in the freely evolving inelastic hard disk system, which are homogeneous, shearing (vortex), clustering and final state. In the shearing stage, the self-organized macroscopic coherent vortices become dominant. In the clustering stage, the energy spectra are close to the expectation of Kraichnan-Batchelor theory and the squared two-particle separation strictly obeys Richardson law.

The dynamics of granular materials becomes one of the most important topics in the studies of nonlinear, dissipative, and non-equilibrium statistical physics. Granular media are collections of macroscopic particles with rough surfaces, dissipative, and frictional interactions. The granular systems require an energy source in order to be in a steady state and the external gravitational force much affect their dynamics.

To focus on the dissipative features, smooth inelastic hard sphere (IHS) model are often used as an ideal model. The freely cooling granular fluid has been studied as an ideal dissipative particle system in the absence of external force. Since the system is only composed of an inelastic hard spheres and no relevant energy scale exists, the restitution coefficient between collision particles is the only parameter in terms of non-equilibrium. The assumption of an inelastic hard sphere potential is also employed in kinetic theory, which facilitates comparisons between theory and simulation. In order to construct the theory of the macroscopic phenomenology in non-equilibrium dissipative particle system, an IHS model is the most promising as a microscopic model. A linear stability analysis of hydrodynamic equations for IHS model has revealed that the initial spatially homogeneous cooling state is unstable to the formation of vortices and clusters. The shearing (vortex) and cluster instabilities were theoretically predicted and were tested by molecular dynamics (MD) simulations .

Since the total energy is monotonically decreasing in the freely evolving process, a steady state in terms of energy fluctuation can be realized by scaling the velocity of the entire particle to the total energy remaining constant . Here, we introduce the new-scaled time $`t_s`$ which is described by

$$t_s=_0^t\frac{dt}{\beta (t)},\beta (t)=\sqrt{T(0)/T(t)},$$
(1)
where $`T(t)`$ is the average kinetic energy per particle as a function of usual time $`t`$. The great advantage of this scaling operation is that the trajectories of particles don’t change compared with the non-scaled case in case of hard sphere system. Therefore, one simply replaces the usual time $`t`$ with the new-scaled time $`t_s`$. This operation is the same as the well-known velocity scaling method  in the usual MD simulation, in which the velocities of all particle $`𝐯_i(t)`$ are scaled following each collision by the factor $`\beta (t)`$ (i.e. $`\beta (t)𝐯_i(t)`$) and the total energy is kept fixed strictly all over the time.

The 2D turbulence in nature is a large-scale fluid motion in the atmosphere or ocean dynamics on earth. The most remarkable feature of 2D turbulence is described by the enstrophy cascade dynamics, which is completely different from that of 3D turbulence represented by the K41 theory. The existence of enstrophy cascade process was originally proposed by Kraichnan and Batchelor. The theory expected that the enstrophy injected at a prescribed scale is dissipated at smaller scales, undergoing a cascading process at a constant enstrophy transfer rate; this led to predicting a $`k^3`$ spectrum for the energy, in a range of scales extending from the injection to the dissipative scale. The granular kinetic energy spectrum in the connection with the fluid turbulence was firstly pointed out by Taguchi  in 2D granular vibrated beds. He obtained the results of the $`k^{5/3}`$ spectrum in his simulation with a few hundred particles. Another important nature of 2D turbulence is the self-organized coherent vortices, which develop into larger ones through the merging process between vortices with same sign of circulation.

In this letter, we especially focus on the velocity statistics and statistical laws on the fluid turbulence. To specify what is the universal characters in dissipative system both macroscopic equation and microscopic dissipative particle, we performed extensive event-driven molecular dynamics simulation systematically on a freely cooling process in 2D IHS model with velocity scaling thermostat. We found a strong similarity between 2D IHS model and 2D NS fluid turbulence.

The freely 2D IHS (granular) model is so simple that the system is completely characterized with only three dimensionless parameters: the restitution coefficient $`r`$, the total number of disks $`N`$, and the packing fraction $`\nu `$. The system size in the unit of disk diameter $`d`$ is $`L/d=\sqrt{\pi N/\nu }/2`$. All the disks are identical, namely the system is monodispersed. The collision is instantaneous and only binary collisions occur. When two disks, $`i`$ and $`j`$, with respective velocities $`𝐯_i`$ and $`𝐯_j`$ collide, the velocities after the collision, $`𝐯_i^{}`$ and $`𝐯_j^{}`$, are given by

$`𝐯_i^{}`$ $`=`$ $`𝐯_i{\displaystyle \frac{1}{2}}(1+r)[𝐧(𝐯_i𝐯_j)]𝐧`$ (2)
$`𝐯_j^{}`$ $`=`$ $`𝐯_j+{\displaystyle \frac{1}{2}}(1+r)[𝐧(𝐯_i𝐯_j)]𝐧,`$ (3)
where $`𝐧`$ is the unit vector parallel to the relative position of the two colliding disks in contact. Our system consists of more than 250 thousands hard disks (up to a few million) placed in a square box with periodic boundaries without no external force. To perform such a large-scale simulation, we implemented the simple and efficient event-driven algorithm, which can actually simulate more than a few million particles even in the personal computer . Initially, the system is prepared as the equilibrium state by the long enough preliminary run with the restitution coefficient $`r=1`$, in which the density is uniform and the disk velocities are Maxwell-Boltzmann distribution. The packing fraction and the restitution coefficient ($`\nu ,r`$) were varied from dilute to dense and within shearing regime, which is estimated by the criterion of McNamara and Young , respectively. In case of $`(N,\nu )=(1024,0.25)`$, McNamara and Young have found the final states have three typical states, which are kinetic, shearing, and collapse regime. However, the spatio-temporal structure of shearing regime have not known yet especially at the macroscopic level.

The criterion of kinetic-shearing boundary in Ref.  is based on the results of Jenkins and Richman , in which the high wave number cutoff for the unstable shear modes were derived. On the contrary, the shearing-collapse boundary is estimated by 1D theoretical analyses for inelastic collapse, which are known as the phenomenon on the number of collisions during a finite time diverges. By using the the regime criterion described by McNamara and Young , one can find both regime boundaries become close to the unity in the thermodynamic limit. These theoretical expectation indicate that the system always unstable even in the quasi-elastic limit. This fact implies the important conjecture discussed later when we consider a large-scale IHS model as the macroscopic fluid model.

In the large-scale simulations, the restitution coefficients within shearing regime become quasi-elastic ($`r1`$). All our simulations by changing various parameters within the shearing regime, the system evolves to the final steady state after several stages. As described in McNamara and Young , we can calculate the packing fraction $`\nu (𝐱)`$, velocity $`𝐮(𝐱)=(u_x(𝐱),u_y(𝐱))`$, and temperature $`T(𝐱)`$ at any point $`𝐱`$. Figure 1 represents typical evolving process for four normalized properties as a function of new-scaled time $`t_s`$ in 2D IHS model. During relaxation process, four stages can be distinguished. After the homogeneous cooling state (HCS) continues for a certain time from initial thermal equilibrium state(first stage), the short-range velocity correlation for each time (i.e. both pre- and postcollisional velocity correlation) within the distance $`_𝐱𝐮(𝐱)𝐮(𝐱^{})`$ sharply deviate from zero (solid line in Fig.1), where $`𝐱^{}`$ is location around $`𝐱`$ with a distance of disk diameter $`d`$ (second stage). Note that there are several works on the existence for the short-range velocity correlations even in HCS . Our simulation also seems to show that velocity correlation ’gradually’ increase from the beginning of simulation in first stage. Therefore, it might to be difficult to determine the exact time between first and second stage. In this stage, coherent vortices self-organize and the coherent vortices develop into larger ones through the mutual confluence and the merging process among vortices with same sign of circulation (Fig. 2 (a)). These self-organized vortices are found firstly by McWilliams  in the direct numerical simulation (DNS) of 2D NS fluid turbulence. In the third stage, the density fluctuation $`\nu _{rms}`$ (dotted line in Fig.1), which is calculated the square root of the space average $`(\nu (𝐱)\nu )^2`$, exhibit instability compressive flow (Fig. 2 (b)). Finally, steady state is realized (fourth stage). In the final steady states, the spatial correlation of inelastic hard disks, which gradually increases from the beginning of simulation, might be reached beyond the system size and begin to interfere with each other though the periodic boundary condition. Our simulations, in the shearing regime, there is no sign for the inelastic hard disks assembling to one cluster during the simulation time. This is because the shear mode expanded to the whole system might be stable though the periodic boundary condition. Therefore, we call them ‘final steady states’. We found there are four characteristic final steady states, which are shear (laminar, oscillating, and turbulent) and vortex (one pair of vortices with opposite sign of circulation) flows, by changing both packing fraction and the restitution coefficient within shearing regime, systematically. Vortex flows of final pattern are also observed in the DNS simulation for 2D incompressible turbulence with the periodic boundary condition.

The velocity anisotropy $`A=(u_x(𝐱)^2u_y(𝐱)^2)/(u_x(𝐱)^2+u_y(𝐱)^2)`$, in which one can distinguish the final state as a shear ($`A=1`$ or $`1`$) or a vortex ($`A=0`$) flows, and the enstrophy $`Z=|\omega (𝐱)|^2`$, where $`\omega (𝐱)=rot𝐮(𝐱)`$, are also plotted by dot dashed line and dashed line in Fig. 1, respectively. We confirmed that the total vorticity $`_𝐱\omega (𝐱)`$ is zero throughout the simulation.

The 2D fluid turbulence has different characters on the statistical law between forced and freely decaying case. However, in the granular case, no systematic consideration seems to exist yet. In the previous studies related to energy spectrum in granular material, energy injection (thermostat) are driven by the vibration cycle  and the periodical-stochastic thermostat . These energy injections resemble in the studies of forced fluid turbulence. On the other hand, velocity scaling thermostat , in which the system is driven continuously, is thought as corresponding to a freely decaying case, because the statistics itself don’t change by introducing the new-scaled time $`t_s`$. Actually, we found the energy spectra (power spectra of velocity field) in 2D IHS model with the velocity scaling thermostat are close to the expectation of Kraichnan-Batchelor theory ($`E(k)k^3`$) after the third stage (Fig. 3). Therefore, our simulations show the enstrophy cascade, as is expected by the theory for freely decaying 2D fluid turbulence. We have confirmed this power law by changing several different system size. In Fig. 3, we can estimate the characteristic spatial scale ($`k_d0.3`$) for minimal dissipative domain (such as Kolmogorov scale in the fluid turbulence), which is composed of about a thousand disks.

The first quantitative phenomenological observation in developed turbulence was shown by Richardson , in which the two fluid particle separation $`R=|𝐫_i𝐫_j|`$ obeys power law ($`R^2t^3`$). We also show that the time ($`t_s`$) dependence of the space-averaged squared two-disk separation in 2D IHS model strictly obeys Richardson dispersion law in the third stage (Fig. 4). In the inset of Fig. 4, the enstrophy evolution is also plotted in terms of the new-scaled time $`t_s`$. The enstrophy seems to decay power-law behavior in the second stage, in which the coherent vortices self-organize. However, since the second stage itself is relatively short, this behavior needs further confirmation. As the intermittency of vorticity, McWilliams found that the probability distribution function of vorticity significantly deviate from Gaussian . By calculating flatness of vorticity ($`f_\omega =\omega ^4/\omega ^2^2`$) in 2D IHS model, our simulations also show the deviation from Gaussian after third stage.

How should we understand the obtained results? The different points between fluid turbulence and granular turbulence are compressibility, the origin of dissipation (i.e. viscosity and inelasticity between particles) and the ratio of particle size and system size. Is the granular turbulence close to fluid turbulence in the thermodynamic limit? Let us assume the extreme condition that is dense, thermodynamic and elastic limit. In this situation, a little amount of dissipation in the system always results in unstable state. We also obtained the fact that the velocity correlation length (Kolmogorov scale) becomes larger in dense (i.e., quasi-incompressible) system, but less than system size when we consider the thermodynamic limit. Therefore, this extreme condition seems to really correspond to NS fluid turbulence.

In this letter, we showed the remarkable similar aspects on the statistical law of vorticity between 2D IHS (granular) turbulence and 2D NS fluid turbulence. These results were obtained by only solving a simple Newton’s equation system for inelastic hard disks in terms of event-driven scheme. From the microscopic dynamics of inelastic hard disk to the macroscopic fluid, it is important to study the origins of statistical law for turbulence at the microscopic level, but there are very few studies so far from this point of view. The discussion for three limits (dense, thermodynamic, and elastic) might make a connection between 2D granular turbulence and 2D fluid turbulence.

I would like to thank to Profs. Y. Hiwatari, H. Nakanishi, H. Hayakawa, and S. Nosé for making helpful comments. I also acknowledge helpful discussion with Prof. T. Gotoh and Dr. T. Watanabe on the two dimensional fluid turbulence. This work was partially supported by the Ministry of Education, Science, Sports and Culture, Grant-in-Aid for Scientific Research (C), 15560042, 2003.


text:
# A comparison of methods for confidence intervals

## I P.d.f. for life time estimators

For a given value $`\tau `$ of the true life time, the P.D.F. of a measurement is

$$\frac{\mathrm{d}W}{\mathrm{d}t}=\frac{1}{\tau }\mathrm{exp}\left(\frac{t}{\tau }\right),$$

and so for an experiment with $`n`$ measurements

$$\mathrm{d}W=\frac{1}{\tau ^n}\underset{k=1}{\overset{n}{}}\mathrm{d}t_k\mathrm{exp}\left(\frac{t_k}{\tau }\right).$$
(1)
The negative log likelihood function is

$$L=n\mathrm{ln}\tau +\frac{1}{\tau }\underset{k=1}{\overset{n}{}}t_k.$$
(2)
The maximum likelihood estimator of the lifetime can easily be found minimizing $`L`$

$$\widehat{\tau }=\frac{1}{n}\underset{k=1}{\overset{n}{}}t_k;\mathrm{min}L=L_0=n+n\mathrm{ln}\widehat{\tau },$$
(3)
so the probability (1) can be transformed to

$$\frac{\mathrm{d}W}{\mathrm{d}\widehat{\tau }}=\frac{1}{(n1)!}\left(\frac{n\widehat{\tau }}{\tau }\right)^{n1}\frac{n}{\tau }\mathrm{exp}\left(\frac{n\widehat{\tau }}{\tau }\right).$$
(4)
Given some true value $`\tau `$ then for any algorithm that defines a confidence interval $`\widehat{\tau }_{\mathrm{\Delta }\tau ^{()}}^{+\mathrm{\Delta }\tau ^{(+)}}`$ we can evaluate the coverage $`P`$:

$$P=\frac{1}{(n1)!}\frac{n}{\tau }\underset{\widehat{\tau _1}}{\overset{\widehat{\tau _2}}{}}\left(\frac{n\widehat{\tau }}{\tau }\right)^{n1}\mathrm{exp}\left(\frac{n\widehat{\tau }}{\tau }\right)d\widehat{\tau },$$
(5)
where $`\widehat{\tau _1}+\mathrm{\Delta }\tau ^{(+)}=\tau ;\widehat{\tau _2}\mathrm{\Delta }\tau ^{()}=\tau .`$

## II Likelihood Function confidence interval

The conventional Likelihood function method for finding a 68% confidence interval Hudson ; LF\_method1  is to find the values of $`\tau `$ for which

$$\mathrm{\Delta }L=LL_0=\frac{1}{2}.$$

In our case

$$\mathrm{\Delta }L=n\left(\frac{\widehat{\tau }}{\tau }1+\mathrm{ln}\frac{\tau }{\widehat{\tau }}\right)=\frac{1}{2}.$$
(6)
For example, for $`n=5`$ the limits are

$$[\tau _1,\tau _2]=[0.6595\widehat{\tau },1.6212\widehat{\tau }].$$
(7)
The coverage of this interval, from Equation (5), is

$$\frac{1}{4!}\frac{5}{\tau }\underset{\widehat{\tau _1}}{\overset{\widehat{\tau _2}}{}}\left(\frac{5\widehat{\tau }}{\tau }\right)^4\mathrm{exp}\left(\frac{5\widehat{\tau }}{\tau }\right)d\widehat{\tau }=0.6747,$$

where the integration limits, corresponding to (7), are

$$[\widehat{\tau _1},\widehat{\tau _2}]=[0.6168\tau ,1.5163\tau ].$$

The coverage is close to, but significally different from, the nominal value of 0.6827.

Examples of confidence intervals obtained by this means are shown in Table 1, as the values in parantheses. The 95% confidence interval was obtained using the rule $`\mathrm{\Delta }L=2`$, and 90% upper limit using a one side interval for which

$$\mathrm{\Delta }L=\left[\text{erf}^1\left(20.91\right)\right]^20.821.$$

The coverage given by such intervals is shown in Fig. 1, evaluated using a Monte Carlo method.

## III Bayesian confidence interval

For comparison we can estimate a Bayesian confidence interval for the same example of $`n=5`$. In the Bayesian approach Jeffreys ; BayesAll ; Neyman2 , the likelihood function is considered to be a probability density for the true parameter $`\tau `$. Assuming a flat prior distribution for $`\tau `$ this is

$$\frac{\mathrm{d}W}{\mathrm{d}\tau }=P(\tau )=\frac{1}{\tau ^n}\mathrm{exp}\left(\frac{n\widehat{\tau }}{\tau }\right).$$

After normalization (for $`n2`$) this becomes:

$$P(\tau )=\frac{\left(n\widehat{\tau }\right)^{n1}}{(n2)!\tau ^n}\mathrm{exp}\left(\frac{n\widehat{\tau }}{\tau }\right),$$

which for $`n=5`$ gives

$$\underset{0}{\overset{\tau }{}}dW=\left[1+\frac{5\widehat{\tau }}{\tau }+\frac{1}{2}\left(\frac{5\widehat{\tau }}{\tau }\right)^2+\frac{1}{6}\left(\frac{5\widehat{\tau }}{\tau }\right)^3\right]e^{\frac{5\widehat{\tau }}{\tau }}.$$
(8)
The 68% central confidence region for this distribution is (see Fig. 2):

$$\tau =\widehat{\tau }\left(1_{0.1552}^{+1.3974}\right)$$

The coverage of this region is actually not 68.27% but 64.31%.

## IV Neyman’s confidence interval

Neyman Neyman1 ; Neyman2 ; Barlow  proposed a frequentist construction of a confidence zone (or confidence belt) as follows (see Figure 3):

1. One obtains functions $`\widehat{\tau }_1(\tau )`$ and $`\widehat{\tau }_2(\tau )`$ of the true parameter $`\tau `$ such that

$$\begin{array}{c}\underset{0}{\overset{\widehat{\tau }_1(\tau )}{}}\frac{\mathrm{d}W(\widehat{\tau };\tau )}{\mathrm{d}\widehat{\tau }}d\widehat{\tau }=\frac{1\beta }{2};\hfill \\ \underset{\widehat{\tau }_2(\tau )}{\overset{\mathrm{}}{}}\frac{\mathrm{d}W(\widehat{\tau };\tau )}{\mathrm{d}\widehat{\tau }}d\widehat{\tau }=\frac{1\beta }{2},\hfill \end{array}$$

where $`\beta `$ is the confidence level required, here $`\beta =0.6827`$. For $`n=5`$ these are simply $`\widehat{\tau }_1(\tau )=0.568\tau `$, $`\widehat{\tau }_2(\tau )=1.433\tau `$, as shown in Figure 3.
2. One defines the inverse functions

$$\tau _1(\widehat{\tau })=\widehat{\tau }_2^1(\widehat{\tau });\tau _2(\widehat{\tau })=\widehat{\tau }_1^1(\widehat{\tau })$$

which, for a given value of $`\widehat{\tau }`$, define the borders of the confidence interval for $`\tau `$, with coverage $`\beta `$.

In our example, there are $`\tau _1(\widehat{\tau })=0.698\widehat{\tau }`$, $`\tau _2(\widehat{\tau })=1.760\widehat{\tau }`$.

Thus the result of a lifetime experiment of this type can be written

$$\tau =\widehat{\tau }\left(1_{0.302}^{+0.760}\right).$$

The coverage evaluated is 0.6826 — the difference of 0.0001 is purely due to rounding errors.

Table 1 shows these intervals for several values of $`n`$, with the likelihood approximation shown in parantheses for comparison.

Table 2 compares the coverage of all three methods for the $`n=5`$ case.

## V  Conclusion

* Neyman’s method for confidence intervals provides exact coverage, by construction.
* The intervals from $`\mathrm{\Delta }L=1/2`$ agree well with the Neyman intervals for large $`n`$, but differ for small $`n`$, as can be seen in Table 1. In such cases they undercover, i.e. the interval is smaller than the true one.
* Bayesian confidence intervals give very different results, and can undercover or overcover.