text:
# An RXTE Archival Search for Coherent X-ray Pulsations in LMXB 4U 1820-30

## 1 Introduction

One of the great scientific expectations when the Rossi X-ray Timing Explorer (RXTE) was launched in 1995 was the discovery of coherent pulsations from low-mass X-ray binaries (LMXBs). At present, only five accreting millisecond pulsars are known. All five are faint transient sources where the pulsations at the spin period of the pulsar were discovered during an outburst (see the contributions of D. Chakrabarty and C. Markwardt to the proceedings of this conference).

There are many things that we can learn from searching for more examples of direct pulsations from LMXBs. For example, coherent pulsations give the precise rotation rate of the neutron star and test the connection between recycled MSPs and LMXBs. They help constrain the models of kHz QPOs and burst oscillations. Timing of the coherent pulses would also allow measurements of the accretion torques, the orbital parameters, and the companion masses. In addition, along with X-ray and optical spectra, pulse timing can give information about the compactness and the equation of state of the neutron star. Finally, setting stringent upper limits to coherent pulsations in several sources in a variety of spectral states will impose constraints on the possible mechanism for supressing coherent pulsations in these sources. For all these reasons, we have started a large-scale search for coherent pulsations from LMXBs.

## 2 Characteristics of 4U 1820$``$30

The source 4U 1820$``$30 is an atoll LMXB in globular cluster NGC 6624. It has an orbital binary period of 685 s (11.4 min) Stella et al. (1987), the shortest known binary orbital period in an LMXB. 4U 1820$``$30 undergoes a regular $``$176 day accretion cycle Priedhorsky and Terrell (1984) switching between high and low luminosity states. In the low state, regular Type I bursts are seen $`\pm `$23 days around the minimum luminosity (Chou and Grindlay, 2001). In the low state, 4U 1820$``$30 has also shown an extremely energetic superburst, likely due to deep ignition of a carbon layer Strohmayer and Brown (2002). Several low frequency QPOs Wijnands et al. (1999) as well as two peaks of kHz QPOs Smale et al. (1997) have been observed from this source. Faint UV and optical counterparts of 4U 1820$``$30 have also been observed King et al. (1993); Anderson et al. (1997). If the secondary star in the system is a white dwarf, its mass is estimated to be 0.058 to 0.078 $`\text{M}_{}`$ Rappaport et al. (1987). If the secondary is a main sequence star, the upper limit on its mass is 0.3 $`\text{M}_{}`$ Arons and King (1993).

Of all known LMXBs, we searched 4U 1820$``$30 first because a lot of high time resolution long data sets of this source are available in the archives, because the presence of kHz QPOs in this source may indicate favorable conditions for detecting pulsations, and because of its small orbital period: the phase modulation search technique that we are using is most sensitive (and provides the biggest increase in sensitivity over previous methods) when the observations are longer than two complete binary orbits Ransom et al. (2003).

## 3 The RXTE Observations

We searched 34 archival RXTE observations collected between 1996 and 2002. The observations were available in event modes, and had a time resolution of 125 $`\mu `$s or better. 32 of the observations had a total time on source bigger than 10 ks. The longest of these was 25 ks long with a total time on source of 16 ks. Some of the observations were segments of longer observations which allowed us to concatenate them and analyse them together. The longest of the concatenated observations was 77.5 ks long with a total time on source of 46.5 ks.

## 4 The data Analysis

Each observation was downloaded, filtered, barycentered, and then split into four energy bands: Soft (2-5 keV), Medium (5-10 keV), Hard (10-20 keV), and Wide (2-20 keV), in order to attempt to maximize the signal-to-noise ratio of the unknown pulsations. This processing used custom Python scripts that call the standard FTOOLs. The processed events were then binned into high time resolution (either 0.122 ms or 0.244 ms) time series. Fast Fourier Transforms of the four time series were then computed. Acceleration searches of FFTs of various durations were performed. ‘Phase-Modulation’ searches were conducted on each of the full duration FFTs. Finally, candidates above our threshold were examined using brute-force folding techniques to determine if they were true pulsations.

## 5 The Search Techniques

After we prepared and binned the data sets, we searched them for pulsations in the Fourier domain using two types of searches which we summarize below.

### 5.1 Acceleration Searches

This is the method traditionally used to compensate for the effects of orbital motion Middleditch and Kristian (1984); Anderson et al. (1990); Johnston and Kulkarni (1991); Wood et al. (1991); Vaughan et al. (1994): A time series is stretched or compressed appropriately to account for a trial constant frequency derivative (i.e. constant acceleration), Fourier transformed, and then searched for pulsations. We used a Fourier-domain variant of this technique Ransom et al. (2002). Unfortunately acceleration searches can detect only the strongest pulsations from systems like 4U 1820$``$30 where the orbital period is short. So while we did not expect acceleration searches to find pulsations, we still used them because they are included in the standard search pipeline that we are going to use to search the other sources.

### 5.2 Phase Modulation Searches

In order to analyze longer observations we used a new ‘phase-modulation’ or ‘sideband’ search technique Ransom et al. (2003). This technique relies on the fact that if the observation time is longer than the orbital period, the orbit phase-modulates the pulsar’s spin frequency: phase modulation results in a family of evenly spaced sidebands in the frequency domain, around the intrinsic pulsar signal. The constant spacing between the sidebands is related to the binary orbital period.

A phase modulation search is conducted by taking short FFTs of the power spectrum of a full observation. The short FFTs cover overlapping portions of various lengths (to account for the unknown semi-major axis of the LMXB) over all the Fourier frequency range of the original power spectrum (to account for the unknown pulsation period). They are then searched for peaks indicating regularly spaced sideband (to detect the orbital period of the LMXB). A detection provides initial estimates of the pulsar period, the orbital period, and the projected radius of the orbit, which are refined by generating a series of complex-valued template responses to correlate with the original Fourier amplitudes (i.e. matched filtering in the Fourier domain) Ransom et al. (2003).

## 6 Results and Upper Limits

Our searches did not detect pulsations.

To set an upper limit on the pulsed fraction of detectable pulsations from 4U 1820$``$30 we used the following procedure: we ran our searches on several data sets containing simulated pulsations of various pulsed fractions for a range of companion masses (see Figure 1) and for several spin frequencies. For the purpose of the simulations, we assumed an orbital period of 685 seconds, the typical length of a long RXTE observation (20.5 ks), an inclination angle of 60 degrees, and a circular orbit. Of the searches that we ran, phase modulation searches run on the full duration data sets were the most sensitive. Figure 1 shows the results of the simulations for this type of search run on data sets containing fake 3 ms pulsations. Every point in the figure corresponds to 100 searched data sets. The shading of every point corresponds to a sigma detection level. The line of sigma = 5 indicates the value of our upper limit on the pulsed fraction for a signal period of 3 ms. This value varied between 0.55% and 1.0% (background subtracted) over the range of companion masses that we used. The value decreased slightly when the period of the signal was larger by a few milliseconds. This means that if there were a 3 ms coherent pulsation coming from 4U 1820$``$30 with a pulsed fraction equal or higher than the upper limit stated above, our searches would have easily (with a 95% confidence) detected it.

Our upper limit on the pulsed fraction of possible signals from 4U 1820$``$30 was about 0.8% (background subtracted) for spin periods under 10 ms. The pulsed fractions of the signals detected from the 5 LMXBs with known signal period ranged between 3% and 7%, for signal frequencies between 2 ms and 6 ms Markwardt et al. (2002); Galloway et al. (2002); Markwardt et al. (2003); van der Klis (2000). This means that our searches would have detected the pulsations from 4U 1820$``$30 if they were as strong as the pulsations detected from the other 5 sources and in the same period range.

## 7 Potential reasons for the Lack of Observed Pulsations

Several theories exist to explain why pulsations with a higher pulsed fraction than our upper limit are not seen. To begin with, an unfavorable rotational geometry, an unfavorable viewing geometry, or both can make pulsations undetectable. The pulsed fraction may be reduced due to scattering in the surrounding medium Titarchuk et al. (2002). Gravitational lens effects on the emission from the hot polar caps can also greatly reduce the pulsed fraction Wood et al. (1988); Meszaros et al. (1988). Finally, screening of the stellar magnetic field by the accreting matter can lead to a less preferential heating of the surface (i.e. no polar caps) and therefore to the absence of observable pulsations Cumming et al. (2001).

This work was supported by the Natural Sciences and Engineering Research Council (NSERC) Julie Payette research scholarship. It would also have been impossible to do this work without the computing power of the Beowulf computer cluster operated by the McGill University Pulsar Group and funded by the Canada Foundation for Innovation.


format: text

text:
# 1 Introduction.

## 1 Introduction.

The phrase ”Information is physical”, enunciated by Rolf Landauer  is the core of the research field known as physics of information. The author claimed that information is everywhere, and, whether it is in biological systems, in a computer, or in a black hole, inevitably has a physical form. Also information, of whatever kind, is associated with matter, radiation or fields of different type. It is impossible to get away from this sort of physical embodiment. Therefore, manipulation of information is inevitably subject to the laws of physics.

On the other hand there is a particular physical conjecture known as the Holographic Principle, which states that the number of the degrees of freedom describing the physics inside a volume (including gravity) is bounded by the area of the boundary which encloses this volume.

As information is physical we can also try the other way around, interpreting physical theories in terms of the language of the Information Theory. More precisely, in this paper I will attempt to connect the holography with one of the cornerstone of Information Theory, the Shannon noiseless coding theorem.

## 2 Information Theory .

The classical theory of information is due to a Bell Labs mathematician, Claude Shannon, who in two seminal works definitively laid down its principles. The works of Shannon deal with two very important tasks of information processing :(i)file compression and (ii)error correction. The Shannon noiseless coding theorem -also called source coding theorem- addresses the first issue: How much redundancy can be eliminated from a given source text without losing information? The Shannon noisy coding theorem addresses the second issue: How much redundancy must be added to a given source text in order to guarantee error free communication over a noisy channel.

In this paper we will concentrate just on the first theorem and its relation with holography.

### 2.1 Definitions and example

Information is discretized: it comes in irreducible packages. The elementary unit of classical information is the bit(or cbit, for classic bit), a classical system with only two states, $`0`$ and $`1`$. Any text can be coded into a string of bits, .

Definitions

* Let $`S:=\{s_1,s_2,\mathrm{},s_n\}`$ be a finite alphabet, endowed with a probability distribution $`P:s_jp_s(s_j)`$
* Codeword - a string of bits which represents a symbol in the alphabet.
* Code, codebook - the set of all possible codewords, e.g., the 4-bit binary code has 16 code-words
* Let $`\mathrm{\Gamma }=(S,P)`$ be a source. An encoding scheme for $`\mathrm{\Gamma }`$ is an ordered pair $`(C,f)`$, where $`C`$ is a code and $`f:SC`$ is an injective function, called encoding function.

For the purposes of noiseless encoding, the measure of efficiency of an encoding scheme is its average codeword length.

Definition: The average codeword length of an encoding scheme $`(C,f)`$ for a source $`\mathrm{\Gamma }=(S,P)`$ is defined by

$$AveLen(C,f)=\underset{i=1}{\overset{N}{}}p_ilen(f(s_i))$$
(1)
where $`len(f(s_i))`$ is the length of the string that represents $`s_i`$.

The Noiseless Coding Theorem, first proved by Claude Shannon in 1948, asserts that the average codeword length of any instantaneous code is bounded from below by $`H(C)`$, i.e.$`H(C)AveLen(C)`$, where $`H`$ is the so called Shannon’s entropy

$$H(C)=\underset{i=1}{\overset{N}{}}p_ilog_2p_i$$
(2)
The randomness of an information source can be describe by its entropy. The operational meaning of entropy is that it determines the smallest number of bits per symbol that is required to represent the total output.

Example

Consider an alphabet $`\{a,b,c,d\}`$ with probabilities $`\{\frac{1}{2},\frac{1}{4},\frac{1}{8},\frac{1}{8}\}`$. With four alternatives, it seems natural to assign two-bit binary codewords to each symbol, eg, $`\{00,01,10,11\}`$. If this coding scheme is used, the average codeword length is

$$AveLen(l)=\frac{1}{2}2+\frac{1}{4}2+\frac{1}{8}2+\frac{1}{8}2=2.$$
(3)
However, we can do better than this. The entropy of this alphabet is

$$H=\left(\frac{1}{2}log_2\frac{1}{2}+\frac{1}{4}log_2\frac{1}{4}+\frac{1}{8}log_2\frac{1}{8}+\frac{1}{8}log_2\frac{1}{8}\right)=\frac{7}{4}$$
(4)
and Shannon’s theorem says that a code exists which achieves this rate. In this case, it is easy. The code $`\{0,10,110,111\}`$ assigns shorter codewords to the more frequent symbols. The average codeword length is

$$AveLen(l)=\frac{1}{2}1+\frac{1}{4}2+\frac{1}{8}3+\frac{1}{8}3=\frac{7}{4}.$$
(5)
In the simple example above, the Shannon entropy and the average code length were both $`7/4`$ bits/symbol. It happened this way because the probabilities were powers of $`1/2`$. We are not so fortunate in reality and there will always be a difference between the Shannon’s entropy and the achieved average code length. This is quantified by two numbers, the $`\mathrm{𝑒𝑓𝑓𝑖𝑐𝑖𝑒𝑛𝑐𝑦}`$ and the $`\mathrm{𝑟𝑒𝑑𝑢𝑛𝑑𝑎𝑛𝑐𝑦}`$

$$\mathrm{Efficiency}=\frac{H}{AveLen}1$$
(6)
$$\mathrm{Redundancy}=1\mathrm{Efficiency}.$$
(7)
## 3 Holography

Physicists have the standard vision that the degrees of freedom of the world consist of fields filling space. Improving this scenario some theorists proposed that a small distance cutoff will be required in order to make sense of quantum gravity. Therefore according to that hypothesis the world could be seen as three dimensional discrete lattice with spacing order of the Plank length, $`l_p`$. But recently, ’t Hooft and Susskind suggested an idea even more audacious. According to them the combination of Quantum Mechanics and Gravity requires the three dimensional world to be an image of data that can be stored on a two dimensional projection much like a holographic image. This description only requires one discrete degree of freedom per Planck area and yet it is rich enough to describe all three dimensional phenomena. Therefore it is obvious that the assumption of these ideas imply a radical decrease in the number of degrees of freedom for describing the Universe.

The inspiration for the holographic principle comes from the black hole physics where the entropy of a black hole is given by the Bekenstein-Hawking formula,

$$S=\frac{A}{4G}=\frac{A}{l_p^2}log2,$$
(8)
where A is the horizon area. This result was obtained when they realized the striking resemblance between the laws of black hole mechanics and the laws of thermodynamics. Later, a great amount of work has been done in order to find a precise statistical mechanical interpretation of black hole entropy. One could derive the Bekenstein-Hawking formula by counting black hole microstates and using the Boltzmann entropy formula

$$S=lnW,$$
(9)
where $`W`$ is the number of microstates of the system.

This counting is not a simple task because is quite difficult to identified the black hole microstates. However great progress in this direction has been reported.

Therefore, the basic idea of holography is that the entropy scales like area, or in other words, the total number of states is of order

$$W_Ae^A.$$
(10)
This is a very dramatic change of view because, from the quantum field theory experience, generally one expects that if the energy density is bounded then the maximum entropy is proportional to the volume of space and this means

$$W_Ve^V.$$
(11)
Therefore, for our purposes, the most important conclusion that we can extract from the holographic principle is that the radical reduction of the degrees of freedom necessary for the description of the universe can be expressed in the following simple way,

$$W_VW_A.$$
(12)
## 4 Holography as the best Code

The tools and concepts developed in the above sections will help to show the main idea of this work. Thus, the language of Information Theory will be used to deal with Holography.

Firstly we should point out the relation between the Boltzmann and Shannon formulas. Despite the fact that they had different origins they are conceptually equivalent and the formal relation is very simple,222The relation can be obtained when in (2) the outcomes are equally likely, that is, when $`p_i=\frac{1}{W}`$.

$$S=\frac{1}{log_2e}H.$$
(13)
Here we ignore the units of the thermodynamical entropy and we consider just the bit units.

This relation allows to establish the proper bridge between the two main topics of this work.

Now we arrive to our point. First we consider a gravitational system with a corresponding ‘physical alphabet‘ $`𝐅`$. This ‘alphabet‘ can be thought as the source of the dynamic of the system. The next step is to consider the codes.

Let us suppose that there exits an ‘holographic code’ $`𝐀`$ for this gravitational system and the corresponding average codeword length is

$$AveLen(𝐀)=\frac{A}{4}log_2e.$$
(14)
Similarly, if instead of that we take a ‘volumetric code’$`𝐕`$, we assume that

$$AveLen(𝐕)=Vlog_2e.$$
(15)
On the other hand, if the Holographic Principle is satisfied, the Boltzmann entropy is $`S=A/4`$. Therefore, using (13) we get

$$H=\frac{A}{4}log_2e.$$
(16)
Obtaining,

$$H=AveLen(𝐀).$$
(17)
Consequently, with the support of the Noiseless Coding Theorem we can conclude that the ’holografic code’ is the best encoding scheme that can be constructed.

In this context a volumetric quantum field theory can be seen as a highly inefficient encoding scheme where

$$\mathrm{Efficiency}=\frac{A}{V}1/R1.$$
(18)
## 5 Conclusions

Fermat pointed out that nature is economical, and that light will therefore travel along the path that takes the least time. This is quite a good rule but there are some instances where it is erroneous and the modern statement is extended: the ray path along which light travels is such that the time taken holds a stationary value (i.e. minimum, maximum or constant). Extending this argument to physics of information we could speculate that Holography is a manifestation of the economical character of nature. The minimum amount of information is used to describe the Universe.

However, that is not the end of the story. We can go on further in the same direction studying also the possible connections of Holography with the Shannon’s Second Theorem. Also the arguments presented here could be refined using the Quantum Source Coding Theorem, introduced by Schumacher

Taking into account that the main idea of this work is encoded on the title I will finish with the hope that I introduced the less possible redundancy.

ACKNOWLEDGMENT: The author is grateful to Elcio Abdalla for reading the manuscript and useful suggestions. This work was supported by Fundação de Amparo à Pesquisa do Estado de São Paulo (FAPESP).


text:
# The Statistical Properties of Galaxies Containing ULXs

## 1. Introduction

ROSAT HRI observations showed that ultraluminous X-ray objects (ULXs, a.k.a. IXOs), which were previously detected with Einstein (see Fabbiano, 1989), are compact X-ray sources, are quite common in the local Universe, and are typically offset from the optical galactic nucleus (Colbert & Mushotzky, 1999; Roberts & Warwick, 2000). Here we define ULXs as X-ray sources in galaxies that are point-like, extra-nuclear, and have X-ray luminosities in excess of $`10^{39}\mathrm{ergs}\mathrm{s}^1`$. There is considerable debate as to nature of these sources (see van der Marel, 2003; Miller & Colbert, 2003, for reviews), with popular explanations including anisotropic emission from stellar-mass X-ray binaries and intermediate mass black holes. A constraint on any explanation for these sources would be the frequency at which they occur, i.e., what fraction of galaxies harbor a ULX as a function of ULX luminosity. For example, this frequency would have implications for the typical relativistic boosting factor for jet models of ULXs (Koerding, Falcke, & Markoff, 2002). In this paper we examine the frequency of ULX occurrence as a function of ULX luminosity based on the HRI ULX catalog published in Colbert & Ptak (2002). While HRI observations suffer from significant biases against the detection of ULXs (as discussed below), we note that the relatively large field of view of the HRI resulted in a large number of galaxies having been observed, and of course large area surveys are crucial for establishing the bright end of a luminosity distribution.

## 2. Methodology

In addition to Colbert & Ptak (2002), we also made use of a subset of the RC3 galaxy catalog (de Vaucouleurs et al., 1991) where the distance is determined, most often based on Tully (1988). Colbert & Ptak (2002) chose a recessional velocity cut-off of 5000 km s-1. In this paper we assume the cosmology $`H_0`$ = 75 km s-1 Mpc-1 and $`q_0`$ = 0.5 (for consistency with Tully 1988), and therefore this velocity corresponds to 66.7 Mpc. This resulted in a galaxy sample containing 9452 galaxies. We next determined which of these galaxies had “coverage” by computing the overlap of the RC3 $`D_{25}`$ ellipse and the HRI FOV for each observation that was used in Colbert & Ptak (2002). We considered a galaxy as having coverage if $`90`$% of the galaxy was in the FOV. We found that 766 galaxies had coverage, or $`8\%`$ of the RC3 subset.

In order to properly assess the frequency of ULXs, we need to assess how many galaxies were observed by the HRI with sufficient sensitivity to detect a ULX. We computed the flux limit at the position of each RC3 galaxy with coverage by computing the background level in the point source region size. The point source region used by XAssist (Ptak & Griffiths, 2003) is a circle with a radius given by $`4[\mathrm{exp}(\theta /9)+1]`$ arcseconds, where $`\theta `$ is the off-axis angle of the source in arcminutes. The total number of background counts, b, in this area was then determined using the median background counts per pixel (also computed as a part of the XAssist processing of each field). The $`3\sigma `$ detection threshold count rate is given by $`3\sqrt{b}/t`$, where t is the exposure time, and the limiting flux is then $`3\sqrt{b}f_{conv}v(\theta )/t`$ where $`v(\theta )`$ gives the vignetting correction and $`f_{conv}`$ is the on-axis count rate to flux conversion factor. $`f_{conv}`$ was computed using the FTOOL PIMMS, assuming a power-law slope of $`1.7`$ and the Galactic column density (Dickey & Lockman, 1990, using the FTOOL “nh”). We computed the limiting flux in both the 0.5-2.0 keV bandpass $`F_{\mathrm{lim},0.52.0\mathrm{keV}}`$ and the 2-10 keV bandpass $`F_{\mathrm{lim},210\mathrm{keV}}`$. We used $`F_{\mathrm{lim},210\mathrm{keV}}`$ to compute the limiting luminosity since the 2-10 keV bandpass was used in Colbert & Ptak (2002), for consistency with current Chandra and XMM-Newton surveys. Any ULX could then have been detected in a galaxy when $`L_{\mathrm{lim},210\mathrm{keV}}`$ is $`10^{39}`$erg s-1, while more generally $`L_{\mathrm{lim},210\mathrm{keV}}`$ is $`>10^{39}`$erg s-1.

Finally, it is necessary to determine how many ULXs are likely to be background sources. For this purpose we used the logN-logS given in Hasinger (1998), which we note is consistent with a simple power-law $`\mathrm{log}N(>S)100(S_{14})^{1.4}`$ sources deg-2, where $`S_{14}`$ is the 0.5-2.0 keV flux in units of $`10^{14}\mathrm{ergs}\mathrm{cm}^2\mathrm{s}^1`$. The number of expected background sources is then $`n_{xrb}=100(\frac{F_{\mathrm{lim},0.52.0\mathrm{keV}}}{10^{14}})^{1.4}\pi r^2`$, where r is the major axis radius of each galaxy in degrees. Colbert & Ptak (2002) used $`r=D_{25}`$ (the RC3 major axis diameter) however here we also consider $`r=0.5D_{25}`$. Also note that when $`L_{\mathrm{lim},210\mathrm{keV}}`$ is $`10^{39}`$erg s-1, the relevant number of background sources is the number of source whose apparent luminosity exceeds $`10^{39}`$erg s-1, and since we are assuming a simple power-law logN-logS this simply entails decreasing the total number of background sources expected by the factor ($`L_{\mathrm{lim},210\mathrm{keV}}/10^{39}`$erg s-1)1.4.

## 3. Results

### 3.1. ULX Frequency

In Table 1 we show the number of galaxies ($`N_{L_{lim}L_X}`$) with HRI coverage and a limiting 2-10 keV luminosity less than or equal to $`L_X`$ ranging from $`10^{39}`$erg s-1to $`10^{41}`$erg s-1(i.e., the number of galaxies where a ULX as bright as $`L_X`$ or brighter could have been detected). We also show the number of galaxies with at least one ULX ($`N_{wULX}`$), and the number of galaxies with at least one ULX whose luminosity exceeds the limiting luminosity ($`N_{L_{ULX}L_X}`$). This latter quantity, $`N_{L_{ULX}L_X}`$, is of interest since it represents the number of galaxies with a ULX at least as bright as the list $`L_X`$, and we then computed the resultant number of galaxies after subtracting the expected number of background sources ($`N_{net,L_{ULX}L_X}`$, i.e., the number of galaxies where the 1-$`\sigma `$ lower limit on the number of ULXs is greater than 0.; the lower limit was calculated using the mehology described in Kraft, Burrows, & Nousek (1991), with the background being the expected number of background sources). Finally we list the values $`N_{wULX}`$, $`N_{L_{ULX}L_X}`$, and $`N_{net,L_{ULX}L_X}`$ using both $`r=D_{25}`$ and $`r=0.5D_{25}`$. Note that in the latter case the total number of ULXs in the Colbert & Ptak (2002) sample is reduced from 87 to 49.

In Table 1 shows that a large fraction of ULXs detected at radii larger $`0.5D_{25}`$ are consistent with background sources, i.e., the net number of galaxies with ULXs is often larger when $`0.5D_{25}`$ is chosen as the limiting radius. Also, in some cases, the number of ULXs detected within the $`0.5D_{25}`$ radius is above the background level when $`r=0.5D_{25}`$ but is not above the background level when $`r=D_{25}`$ (i.e., since in the latter case the expected number of background sources is a factor of four larger). In Figure 1 we plot the fraction of galaxies harboring a ULX as a function of ULX luminosity ($`N_{net,L_{ULX}L_X}/N_{L_{lim}L_X}`$) for both $`r=D_{25}`$ and $`r=0.5D_{25}`$. Here the errors were computed according to the methodology described in Kraft, Burrows, & Nousek (1991), with the background taken as the number of galaxies at each luminosity that are consistent with all ULXs being background sources (i.e., $`N_{L_{ULX}L_X}N_{net,L_{ULX}L_X}`$). The choice of effective radius for the galaxies does not significantly impact the estimate of ULX frequency once background sources are taken into account. The frequency of ULXs observed by the HRI evidently is $`12\%`$ and $`1\%`$ of galaxies harboring a ULX with $`L_X>=10^{39}`$erg s-1and $`L_X>=10^{40}`$erg s-1, respectively. We also plot the galaxy fractions for spirals only (i.e., Hubble type in RC3 $``$ 0). The Hubble type distribution of galaxies with and without a ULX is discussed below, however we note that the fraction of galaxies harboring a ULX is slightly higher in spirals.

### 3.2. Hubble Type Distribution

In Figure 2 we show the Hubble type distribution of the RC3 subsample (discussed in §2 ), the Hubble type distribution of galaxies that had HRI coverage (hereafter HRI galaxies) and the HRI galaxies that contain at least one ULX. The ULX galaxy type distribution is shown both before and after background subtraction, and here we only consider $`r<0.5D_{25}`$ ULXs. There is a slight excess of elliptical ULX galaxies relative to the RC3 sample however the larger excess in the HRI galaxy sample suggests that this bias due in part to a bias toward early type galaxies in the HRI observing patterns relative to RC3 galaxies in general. The most conspicuous feature of the ULX galaxy distribution is the excess of T=5-7 galaxies, corresponding to Sc and Scd galaxies. However, in general the ULX galaxy type distributions are not significantly different statistically from the HRI RC3 galaxy type distributions, with the $`\chi ^2`$ tests giving $`\chi ^2/dof15/8`$, which can only be rejected at the 95% level.

We also investigated the possible connection between star-formation rate (SFR) and ULX frequency by tabulating the far-infrared luminosities (a well-known SFR indicator; Kennicutt 1998) of the galaxy samples discussed above. These distributions are shown in Figure 3, where we used NED IRAS fluxes densities at 60 and 100 $`\mu m`$ to compute $`F_{FIR}=1.26\times 10^{23}[2.58\times 10^{12}f_{60\mu m}+10^{12}f_{100\mu m}]`$. The peak of the distributions with a ULX is slightly skewed toward higher FIR (i.e., near $`\mathrm{log}L_{FIR}=43.5`$) however the means of the distributions do not differ significantly (all are in the range log $`L_{FIR}`$ = 42.7-43.0). The HRI galaxy FIR luminosity distribution differs from the “net” ULX galaxy FIR luminosity distribution at less than the 90% level based on the $`\chi ^2`$ test. Restricting the sample to spiral galaxies (T $`>=`$ 0) shifted the means by only $`\mathrm{log}L_{FIR}0.1`$, reflecting the predominance of spirals in these samples.

Similarly we plot the distribution of $`L_{FIR}/L_K`$ where $`L_K`$ is the K band luminosity computed using the 2mass Extended Source catalog K band magnitude (see Jarrett et al., 2003) and assuming a zero-point flux density of $`6.7\times 10^{21}`$erg cm-2 s-1$`Hz^1`$ (Zombeck, 1990)111see also
http://www.ipac.caltech.edu/2mass/releases/allsky/doc/sec6\_4a.html. The K band luminosity is a proxy for stellar mass since the variation in mass-to-light ratio for K band is only a factor of $`2`$ for various galaxy types (Bell & de Jong, 2001). These distributions are shown in Figure 4 and again we do not see any significant differences between galaxies with ULXs and galaxies observed by the HRI in general.

## 4. Summary & Conclusions

We have assessed the frequency of ULX occurrence in galaxies as a function of ULX luminosity, including proper assessment of the contamination from background sources. We have found that assuming an effective radius of $`D_{25}`$ tends to result in a significant contamination from background sources. We found frequencies of $`12\%`$ and $`1\%`$ for ULX luminosities of at least $`10^{39}`$erg s-1and $`10^{40}`$erg s-1. The Hubble type distribution of galaxies containing at least one ULX is significantly different from those that do not contain a ULX in general, however it is not significantly different from galaxies observed by the HRI. This implies that any bias in the Hubble type distribution of galaxies observed by HRI with ULXs is mostly due to the observing patterns of the HRI. Recently Irwin et al. (2004) reported on an analysis of ULXs in early-type galaxies and found that “very” luminous ULXs (with $`L_X>2\times 10^{39}`$erg s-1) are consistent with the expected number of background sources. This result is based in part on Chandra data and uses luminosities in the 0.3-10.0 keV bandpass as opposed to our choice of the 2-10 keV band. However, as can be seen in Figure 2b, even when considering all ULXs in our sample (i.e., with $`L_X1\times 10^{39}`$erg s-1), the frequency of early-type galaxies containing ULXs is only marginally in excess of zero (with ULXs occurring in globular clusters as in the case of NGC 1399 contributing to the early-type galaxy ULX frequency).

The main caveat with an analysis of the HRI observations is the softness of its bandpass (i.e., 0.1-2.2 keV). Objects absorbed by column densities in excess of $`10^{21}\mathrm{cm}^2`$ are therefore difficult with the HRI. Since columns of this order are typically observed through the disks of galaxies, it is likely that ULXs are being missed that are close to the nucleus of the galaxies and/or are observed in edge-on disks. Chandra and XMM-Newton observations will therefore be much more adept at detecting ULXs in these types of environments. This is particularly true for Chandra whose much smaller PSF relative to the HRI and XMM-Newton allows nuclear ULX to be discerned from AGN at smaller offsets (i.e., 1” compared with 10”). Sipior (2003) finds that 8/41 $`20`$% of nearby galaxies in a relatively unbiased Chandra survey harbor a ULX. This suggests that our HRI survey is missing of order $`40\%`$ of ULXs (albeit with large uncertainty due to the size of the Chandra sample). Hornschemeier et al. (2003) reports that the frequency of ULXs with $`L_X>=2\times 10^{39}`$erg s-1is $`8\%`$ locally and $`36\%`$ at z $`0.1`$ (Hornschemeier et al., 2003) (suggesting that the frequency of ULXs may evolve). We are similarly finding $`8\%`$ (see Table 1) which implies that any biases in our survey are not particularly large for $`L_X>=2\times 10^{39}`$erg s-1ULXs, within the errors.

The lack of any obvious connection between FIR luminosity (and hence SFR) or the FIR/K band luminosity ratio distribution and the presence of a ULX suggests that if such a connection exists, then the biases discussed above are particularly acute in star-forming galaxies. This may not be surprising given that high-SFR galaxies tend to contain significant amounts of molecular clouds and dust which would obscure soft X-rays. In addition the relatively small size of our sample may also be impacting this result. For example, only 23 galaxies have $`\mathrm{log}L_{FIR}>`$43.5 erg s-1and enough HRI sensitivity to detect a ULX, and 4 harbor a ULX (after “background subtraction” and using a galaxy radius of $`0.5D_{25}`$). This amounts to a larger ULX frequency of $`17_8^{+12}`$%, although still statistically consistent with the full sample result. Recent Chandra surveys of the point source luminosity functions in galaxies suggest that, particularly at high star-formation rates, the point source luminosity function normalization is correlated with the SFR (Grimm et al. 2003, Colbert et al. 2004), and that the luminosity function is flatter in spirals (Colbert et al., 2004) or (similarly) in higher SFR galaxies (Kilgard et al., 2002). Having both a flatter luminosity function slope and larger normalization at higher SFR strongly implies that higher SFR galaxies (i.e., with SFR greater than several to ten $`M_{}\mathrm{yr}^1`$ based on Colbert et al. 2004) would tend to harbor multiple ULXs. This is also suggested anecdotally by the large numbers of ULXs found in merging galaxies such as the Antennae (Fabbiano, Zezas, & Murray, 2001). Accordingly the distribution of $`L_{FIR}`$ and $`L_{FIR}/L_K`$ for galaxies with at least one ULX may not be significantly impacted by an increase in the total number of ULXs in galaxies with large SFR. However the resolution of these issues will ultimately require more detailed statistical analysis based on Chandra and XMM-Newton wide-area surveys.

This work was supported by NASA grant NAG 5-11670. We made use of the High Energy Archive and Research Center at NASA/GSFC and the NASA Extragalactic Database. We thanks the anonymous referee for useful comments.


text:
# Surface modes and multi-power law structure in the early-time electromagnetic response of magnetic targets

## Abstract

It was recently demonstrated \[P. B. Weichman, Phys. Rev. Lett. 91, 143908 (2003)\] that the scattered electric field from highly conducting targets following a rapidly terminated electromagnetic pulse displays a universal $`t^{1/2}`$ power law divergence at early time. It is now shown that for strongly permeable targets, $`\mu _c/\mu _b1`$, where $`\mu _b`$ is the background magnetic permeability, the early time regime separates into two distinct power law regimes, with the early-early time $`t^{1/2}`$ behavior crossing over to $`t^{3/2}`$ at late-early time, reflecting a spectrum of magnetic surface modes. The latter is confirmed by data from ferrous targets where $`\mu _c/\mu _b=𝒪(10^2)`$, and for which the early-early time regime is invisibly narrow.

Remote characterization of buried targets is a key goal in many environmental geophysical applications, such as landmine and unexploded ordnance (UXO) remediation serdp . A tool of choice is the time-domain electromagnetic (TDEM) method, in which an inductive coil transmits low frequency \[typically $`𝒪(100\mathrm{Hz})`$\] EM pulses into the ground. Following each pulse \[terminated rapidly on a ramp timescale $`\tau _r=𝒪(10^2\mu \mathrm{s})`$\], the voltage $`V(t)`$ induced by the scattered field is detected by a receiver coil. Standard TDEM sensors are capable of resolving signals from very small (of order 1 gram) metal targets NH01 , and are therefore well suited to UXO detection.

Low frequency yields increased sensitivity to conducting targets, and increased exploration depth (below 5 m) but leads to nearly complete loss in spatial resolution. Lacking direct spatial imaging capability (enabling straightforward identification), one is reduced to seeking such information indirectly via a careful analysis of the full time dependence of $`V(t)`$. As described in Refs. WL03 ; W03 , this signal is affected by both *intrinsic* (target size, shape, geometry, and other physical characteristics) and *extrinsic* (relative target-sensor position and orientation, transmitter and receiver coil geometries, pulse waveform, etc.) properties, and the key to discrimination is the extraction of the former from the “background” of the latter. The aim of this letter is to further develop such formalism for the early time part of $`V(t)`$ W03 .

Conductor electrodynamics are essentially diffusive, and the basic time scale $`\tau _c=L_c^2/D_c`$ is determined by the target diameter $`L_c`$, and diffusion constant $`D_c=c^2/4\pi \mu _c\sigma _c`$ (in Gaussian units) depending on the (relative) target permeability $`\mu _c`$ and conductivity $`\sigma _c`$. Ferrous targets (e.g., steel) are typically modelled as paramagnetic targets with very large permeability $`\mu _c=𝒪(10^2)`$, and (MKS) conductivity $`\sigma _c=𝒪(10^7\mathrm{S}/\mathrm{m})`$ . Thus, even for targets as small as $`L_c=1`$ cm, one finds decay times $`\tau _c=𝒪(10^2\mathrm{ms})`$, much larger than typical pulse periods. Larger UXO-like targets, which are usually ferrous, have even larger $`\tau _c`$. It will often be the case, therefore, that the *full measured range* of $`V(t)`$ will lie in the early time regime, $`t\tau _c`$. On the other hand, for a nonmagnetic (e.g., aluminum) target of similar size one obtains $`\tau _c=𝒪(1\mathrm{ms})`$, and the measured signal will cover a much broader dynamical range.

The early time response is governed by the dynamics of the screening currents, induced in response to the rapid quenching of the transmitted magnetic field immediately following pulse termination (see Fig. 1). The initial diffusion of these currents inward from the target surface generates a $`t^{1/2}`$ power law in $`V(t)`$, with coefficient reflecting the surface properties of the target W03 . However, underlying this power law is the assumption that $`t/\tau _c`$ is the smallest parameter in the problem. It will now be shown that the background-target permeability contrast $`\mu _b/\mu _c`$ can generate a new small parameter that greatly limits the its range of validity. Specifically, an extended calculation is described that divides the early time regime $`t<\tau _e`$, where $`\tau _e\tau _c`$ is the point at which bulk effects first begin to enter (see further below), into an early-early time regime, $`0<t\tau _{\mathrm{mag}}`$, where the $`V(t)t^{1/2}`$ remains valid, and a late-early time regime $`\tau _{\mathrm{mag}}t\tau _e`$, where a new power law $`V(t)t^{3/2}`$ obtains. In the neighborhood of the *magnetic crossover time*,

$$\tau _{\mathrm{mag}}=\tau _c(\mu _b/\mu _c)^2,$$
(1)
an interpolation between the two power laws occurs (that is visible only if $`\tau _{\mathrm{mag}}\tau _e`$), for which the full functional form is provided. For ferrous targets $`\tau _{\mathrm{mag}}/\tau _e=𝒪(10^4)`$, and the early-early time regime is essentially invisible, and only the $`t^{3/2}`$ behavior is seen, consistent with measured data: see Fig. 2.

At low frequencies the dielectric function in the ground and in the target is dominated by the its imaginary part, $`ϵ=4\pi i\sigma /\omega `$, where $`\sigma (𝐱)`$ is the dc conductivity, and the Maxwell equations may be reduced to a single equation for the vector potential,

$$\times \left(\frac{1}{\mu }\times 𝐀\right)+\frac{4\pi \sigma }{c^2}_t𝐀=\frac{4\pi }{c}𝐣_S,$$
(2)
with magnetic induction $`𝐁=\times 𝐀`$, field $`𝐇=𝐁/\mu `$, and gauge chosen so that the electric field is $`𝐄=(1/c)_t𝐀`$. The transmitter loop generates the source current density $`𝐣_S(𝐱,t)`$. The conductivity and permeability are separated into background \[$`\sigma _b(𝐱)`$, $`\mu _b(𝐱)`$\] and conducting target \[$`\sigma _c(𝐱)`$, $`\mu _c(𝐱)`$\] components, where $`\sigma _c`$, $`\mu _c`$ vanish outside the target volume $`V_c`$. High conductivity contrast, $`\sigma _b/\sigma _c1`$, is required, but $`\mu _b/\mu _c`$ is arbitrary.

Equation (2) is a vector diffusion equation, with contrast $`D_b/D_c=𝒪(10^7\text{}10^9)`$. The “background communication time” between instrument and target separated by distance $`R`$ is $`\tau _b=R^2/D_b`$. For representative values $`\sigma _b=0.1`$ S/m, $`\mu _b=1`$, $`R=3`$ m, one obtains $`\tau _b=𝒪(10^6\mathrm{s})`$, instantaneous even on the scale of $`\tau _r`$. Thus, the electrodynamics of the background may be treated as *quasistatic*: outside the target one may drop the $`_t𝐀`$ term in (2). Further, following pulse termination ($`t>0`$) one has $`𝐣_S0`$, hence $`\times 𝐇=0`$, and one may derive $`𝐇=\mathrm{\Phi }`$ from a scalar potential $`\mathrm{\Phi }`$, which in turn satisfies $`(\mu _b\mathrm{\Phi })=0`$. The external field is therefore entirely determined by the internal field through the boundary condition Jackson  $`\widehat{𝐧}[\mu _b𝐇_b\mu _c𝐇_c]=0`$, where $`𝐇_b`$, $`𝐇_c`$ are the fields just outside and inside the target surface, and $`\widehat{𝐧}`$ is the local surface normal. Thus, one obtains a Neumann-type boundary condition, $`\widehat{𝐧}\mathrm{\Phi }=(\mu _c/\mu _b)\widehat{𝐧}𝐇_c`$, and hence the formal solution

$$\mathrm{\Phi }(𝐱)=_{V_c}d^2rg_N(𝐱,𝐫)\frac{\mu _c}{\mu _b}\widehat{𝐧}𝐇_c(𝐫),$$
(3)
where $`V_c`$ is the target surface, and $`g_N`$ is the Neumann Green function satisfying $`\mu _b^1(\mu _bg_N)=\delta (𝐱𝐱^{})`$ with boundary condition $`\widehat{𝐧}g_N=0`$.

The initial surface screening current appears in the transverse magnetic boundary condition Jackson ; W03 : $`𝐊(𝐫)=(c/4\pi )\widehat{𝐧}\times [𝐇_b(𝐫)𝐇_c(𝐫)]`$. As indicated in Fig. 1, $`𝐇_c`$ is the same as that prior to pulse termination, while $`𝐇_b`$ is given by (3):

$$𝐊=\frac{c}{4\pi }\widehat{𝐧}\times (\mathrm{\Phi }𝐇_c)_{t=0},$$
(4)
which, via (3), determines $`𝐊`$ entirely in terms of $`𝐇_c`$.

To investigate the evolution of $`𝐊`$ for $`t>0`$, we take advantage of the rapid variation of the fields near the surface with the local vertical coordinate $`z`$. Thus, $`z`$-derivatives dominate (2), and to leading order in the small parameter $`ϵ(t)=\sqrt{D_ct/L_c^2}`$ foot1 , one obtains

$$_Z^2𝐀^{}+_t𝐀=0,z<0,$$
(5)
with initial condition $`𝐄(0)=(1/c)_t𝐀(0)=\sigma _c^1𝐊\delta (z)`$. Here $`Z=z/\sqrt{D_c}`$, $`𝐀^{}=𝐀\widehat{𝐧}(\widehat{𝐧}𝐀)`$ is the tangential part of $`𝐀`$, and $`\mu ,\sigma ,D`$ are treated as constants on either side of the boundary. One may choose $`\mathrm{\Delta }𝐀=0`$ for $`t=0`$, and it follows immediately that $`\widehat{𝐧}\mathrm{\Delta }𝐀0`$: $`\mathrm{\Delta }𝐀`$ is purely transverse.

One requires now the boundary condition for (5) at $`Z=0^{}`$. For $`t>0`$ let $`𝐀=𝐀(0)+\mathrm{\Delta }𝐀(t)`$ and $`\mathrm{\Phi }=\mathrm{\Phi }(0)+\mathrm{\Delta }\mathrm{\Phi }(t)`$, where one treats $`|\mathrm{\Delta }𝐀|/|𝐀(0)|,|\mathrm{\Delta }\mathrm{\Phi }|/|\mathrm{\Phi }(0)|=𝒪(ϵ)`$ and $`|_Z\mathrm{\Delta }𝐀|/|\mathrm{\Delta }𝐀|=𝒪(1/ϵ)`$. Keeping only leading terms foot1 , continuity of $`\widehat{𝐧}\times 𝐇`$ (the surface current sheet now has finite thickness) and (4) imply that

$$\sqrt{\frac{4\pi \mu _c}{\sigma _c}}𝐊=_Z\mathrm{\Delta }𝐀^{}|_{Z=0^{}}+\mu _c\sqrt{D_c}\widehat{𝐧}\times \mathrm{\Delta }\mathrm{\Phi }|_{Z=0^+},$$
(6)
in which $`\mathrm{\Delta }\mathrm{\Phi }`$ is given by $`\mathrm{\Delta }𝐀(Z=0^{})`$ via (3). In fact, estimating $`_Z=𝒪[D_c^{1/2}/ϵ(t)L_c]`$, $`g_N=𝒪(1/L_c)`$, and $`\widehat{𝐧}\times =𝒪(1/L_c)`$, the ratio of the second term on the right hand side to the first is $`𝒪[ϵ(t)\mu _c/\mu _b]`$. Their relative order is therefore *time-dependent* foot2 . In particular, at early-early time $`ϵ(t)\mu _b/\mu _c`$ (i.e., $`t\tau _{\mathrm{mag}}`$) one may drop the second term. Lack of time dependence in $`𝐊`$ then leads to the simple homogeneous Neumann boundary condition $`_z𝐄=0`$ for the electric field. This is the limit in which the early time analysis in Ref. W03  was carried out and the $`t^{1/2}`$ behavior of $`V(t)`$ derived.

Thus, for $`\mu _c/\mu _b=𝒪(1)`$ one gains nothing by including the $`\mathrm{\Delta }\mathrm{\Phi }`$ term: it is of the same order as other $`𝒪[ϵ(t)]`$ terms previously dropped from (6) foot1 . However, if $`\mu _c/\mu _b1`$ this term dominates for $`\tau _{\mathrm{mag}}<t<\tau _e`$, which now comprises a large fraction of the early time regime. The remainder of this paper is concerned with extending the theory into this regime.

With the $`\mathrm{\Delta }\mathrm{\Phi }`$ term, (6) is nonlocal, coupling $`\mathrm{\Delta }𝐀(Z=0^{})`$ over the surface. In order to decouple (6) we *diagonalize* it by seeking transverse vector eigenfunctions $`𝜶_n`$, $`n=1,2,3,\mathrm{}`$, satisfying

$$\kappa _n𝜶_n(𝐫)=\mu _c\sqrt{D_c}\widehat{𝐧}\times \psi _n(𝐫)$$
(7)
where $`\kappa _n`$ is the corresponding eigenvalue, and $`𝐫`$ is the surface coordinate. The scalar $`\psi _n`$ is defined by

$$\psi _n=\widehat{}_g^N\mu _b^1\widehat{𝐧}\times 𝜶_n,$$
(8)
where, to condense the notation, we define the Neumann Green function operator via $`\widehat{}_g^N\varphi (𝐫)=_{V_c}d^2r^{}g_N(𝐫,𝐫^{})\varphi (𝐫^{})`$. Substituting (7) into (8), one obtains a *scalar* eigenvalue problem,

$$\kappa _n\psi _n=\widehat{}\psi _n,\widehat{}\widehat{}_g^N\mu _b^1\widehat{}_\mathrm{\Delta },$$
(9)
where the generalized surface Laplacian is $`\widehat{}_\mathrm{\Delta }\varphi (𝐫)=\widehat{𝐧}\times \left[\mu _c\sqrt{D_c}\widehat{𝐧}\times \varphi (𝐫)\right]`$. On a sphere one obtains $`\widehat{}_\mathrm{\Delta }=\mu _c\sqrt{D_c}L_c^2𝐋^2`$, where $`L_c`$ is now the radius, and $`𝐋=i𝐱\times `$ is the angular momentum operator Jackson .

As the $`𝜶_n`$ are derived from the scalar $`\psi _n`$ via (7), they do not form a complete set of 2D vector fields. The missing fields consist of the kernel of (7) (the space of functions $`\kappa _n=0`$), spanned by all fields $`𝜷`$ with vanishing normal magnetic field, $`\widehat{𝐧}\times 𝜷=0`$. It follows that $`𝜷=^{}\varphi \widehat{𝐧}\times \widehat{𝐧}\times \varphi `$ for some other scalar $`\varphi `$. Let $`\varphi _n`$, $`𝜷_n`$ be chosen as eigenfunctions of the generalized transverse Laplace equation,

$`\lambda _n\varphi _n(𝐫)`$ $`=`$ $`\widehat{𝐧}\times \left[\mu _c^1D_c^{1/2}\widehat{𝐧}\times \varphi _n(𝐫)\right]`$
$`𝜷_n(𝐫)`$ $`=`$ $`^{}\varphi _n,`$ (10)
with a new set of eigenvalues $`\lambda _n`$.

The set $`\{𝜶_n,𝜷_n\}`$ forms a complete basis, and we perform the *magnetic surface mode expansion*

$`\mathrm{\Delta }𝐀(𝐱,t)`$ $`=`$ $`{\displaystyle \underset{n=1}{\overset{\mathrm{}}{}}}\left[A_n^{(1)}(Z,t)𝜶(𝐫)+A_n^{(2)}(Z,t)𝜷(𝐫)\right]`$
$`\sqrt{{\displaystyle \frac{4\pi \mu _c}{\sigma _c}}}𝐊(𝐫)`$ $`=`$ $`{\displaystyle \underset{n=1}{\overset{\mathrm{}}{}}}\left[K_n^{(1)}𝜶(𝐫)+K_n^{(2)}𝜷(𝐫)\right].`$ (11)
From (7), $`\widehat{}_g^N\mu _b^1\widehat{𝐧}\times \mathrm{\Delta }𝑨=_nA_n^{(1)}\psi _n`$, and from (10), $`\widehat{𝐧}\times (\mu _c^1D_c^{1/2}\mathrm{\Delta }𝐀)=_n\lambda _nA_n^{(2)}\varphi _n`$ (and similarly for $`𝐊`$), hence appropriate orthogonality relations for $`\psi _n,\varphi _n`$ may be used to determine $`A_n^{(i)},K_n^{(i)}`$.

Substituting (11) into (5) and (6), complete separation of the surface modes is achieved:

$$(_t_Z^2)A_n^{(i)}=0,Z<0,$$
(12)
with initial and boundary conditions

$`_tA_n^{(i)}|_{t=0^+}`$ $`=`$ $`K_n^{(i)}\delta (Z)`$
$`(_Z+\kappa _n\delta _{i1})A_n^{(i)}|_{Z=0^{}}`$ $`=`$ $`K_n^{(i)}.`$ (13)
The solutions are $`A_n^{(i)}(Z,t)=K_n^{(i)}H(Z,t;\kappa _n\delta _{i1})`$, where

$`H(Z,t;\kappa )`$ $`=`$ $`\frac{1}{\kappa }\left[\mathrm{erfc}\left(\frac{|Z|}{\sqrt{4t}}\right)e^{\kappa ^2t\kappa Z}\mathrm{erfc}\left(\frac{2\kappa tZ}{\sqrt{4t}}\right)\right]`$
$`H(Z,t;0)`$ $`=`$ $`\sqrt{\frac{4t}{\pi }}e^{Z^2/4t}|Z|\mathrm{erfc}\left(\frac{|Z|}{\sqrt{4t}}\right),`$ (14)
where $`\mathrm{erfc}(x)`$ is the complementary error function.

Fields *external* to the target are obtained by extending $`𝜶_n,𝜷_n`$ into the exterior space. First, let

$`\psi _n(𝐱)`$ $`=`$ $`{\displaystyle \frac{1}{\kappa _n}}{\displaystyle _{V_c}}d^2r^{}g_N(𝐱,𝐫^{}){\displaystyle \frac{1}{\mu _b}}\widehat{}_\mathrm{\Delta }\psi _n(𝐫^{})`$
$`\varphi _n(𝐱)`$ $`=`$ $`{\displaystyle _{V_c}}d^2r^{}g_D(𝐱,𝐫^{})\varphi _n(𝐫^{}),`$ (15)
in which $`g_D(𝐱,𝐱^{})`$ is the Dirichlet Green function satisfying $`\sigma _b^1(\sigma _bg_D)=\delta (𝐱𝐱^{})`$ and vanishing on $`V_c`$, and $`𝐱`$ is no longer restricted to the surface. These definitions guaranteed continuity at the boundary. The vector eigenfunctions are now extended by solving

$`\mu _b\times 𝜶_n(𝐱)=\psi _n(𝐱),(\sigma _b𝜶_n)=0,`$ (16)
while imposing continuity of $`𝜶_n^{}`$ at the surface, and

$$𝜷_n(𝐱)=\varphi _n(𝐱).$$
(17)
Since $`\times 𝜷_n=0`$ it does not contribute to the external magnetic field, and hence to any inductive measurement. The external field is now simply

$$\mathrm{\Delta }𝐀(𝐱,t)=\underset{n=1}{\overset{\mathrm{}}{}}\left[A_n^{(1)}(0^{},t)𝜶_n(𝐱)+A_n^{(2)}(0^{},t)𝜷_n(𝐱)\right],$$
(18)
with time-dependence given by the internal field on the boundary, hence governed by the function

$`H(0,t;\kappa )={\displaystyle \frac{1}{\kappa }}\left[1e^{\kappa ^2t}\mathrm{erfc}\left(\kappa \sqrt{t}\right)\right]`$ (19)
$`\{\begin{array}{cc}\sqrt{\frac{4t}{\pi }}\left[1\frac{1}{2}(\pi \kappa ^2t)^{1/2}+𝒪(\kappa ^2t)\right],\hfill & \kappa ^2t1\hfill \\ \frac{1}{\kappa }\left\{1(\pi \kappa ^2t)^{1/2}+𝒪[(\kappa ^2t)^{3/2}]\right\}\hfill & \kappa ^2t1.\hfill \end{array}`$
The time derivative, entering the electric field and voltage, displays the promised $`t^{1/2}`$ and $`t^{3/2}`$ power laws in the two limits. If $`\psi _n`$ varies on scale $`L_nL_c`$, one may estimate from (9) $`\kappa _n=𝒪[\tau _{\mathrm{mag}}^{1/2}(L_c/L_n)]`$, hence crossover point $`\tau _{\mathrm{mag},n}1/\kappa _n^2=𝒪[\tau _{\mathrm{mag}}(L_n/L_c)^2]`$. One expects $`L_n=L_c`$ only for the fundamental mode, hence (1) actually represents an *upper bound* on the spectrum of crossover times $`\tau _{\mathrm{mag},n}`$.

The exact solution for a homogeneous sphere in a homogeneous background serves to clarify all of the above concepts. The operators $`\widehat{}_g^N`$ and $`\widehat{}_\mathrm{\Delta }`$ now commute and may be simultaneously diagonalized using spherical harmonics: $`\widehat{}_g^NY_{lm}=Y_{lm}/(l+1)L_c`$, $`\widehat{}_\mathrm{\Delta }Y_{lm}=[\mu _c\sqrt{D_c}l(l+1)/L_c^2]Y_{lm}`$. Thus, using $`\psi _{lm}=\varphi _{lm}=L_c^2Y_{lm}`$ one obtains $`\kappa _{lm}=l\tau _{\mathrm{mag}}^{1/2}`$ (hence mode length scale $`L_{lm}=L_c/l`$), and $`\lambda _{lm}=l(l+1)/L_c^2\mu _c\sqrt{D_c}`$. The vector functions are $`𝜶_{lm}(𝐱)=i\mu _b\sqrt{(l+1)/l}(L_c^l/x^{l+1})𝐗_{lm}`$, $`𝜷_{lm}(𝐱)=i\sqrt{(l+1)/l}\times [(L_c^l/x^{l+1})𝐗_{lm}]`$, where $`𝐗_{lm}=[l(l+1)]^{1/2}\widehat{𝐋}Y_{lm}`$ are the vector harmonics Jackson . The exact solution associated with the magnetic modes $`𝜶_{lm}`$ is governed by the usual bulk exponentially decaying mode expansion WL03 ; W03 , from which one obtains $`A_{lm}^{(1)}(t)=K_{lm}^{(1)}\sqrt{4\tau _c}[H_l(0)H_l(t/\tau _c)]`$, where

$$H_l(\tau )=\underset{n=1}{\overset{\mathrm{}}{}}\frac{j_l(\zeta _{ln})^2}{j_l(\zeta _{ln})^2j_{l+1}(\zeta _{ln})j_{l1}(\zeta _{ln})}\frac{e^{\zeta _{ln}^2\tau }}{\zeta _{ln}^2},$$
(20)
where $`j_l(x)`$ are the spherical Bessel functions, and the scaled decay rates are given by the roots of

$$0=\frac{\mu _b}{\mu _c}\zeta _{ln}j_{l1}(\zeta _{ln})+l\left(1\frac{\mu _c}{\mu _b}\right)j_l(\zeta _{ln}).$$
(21)
In Fig. 3 $`H_l(t/\tau _c)`$ is compared to the early time prediction $`H_l(t/\tau _c)H_l(0)(4\tau _c)^{1/2}H(0^{},t;\kappa _{lm})`$, with $`H_l(0)=(\mu _b/2\mu _c)/[l+(l+1)\mu _b/\mu _c]`$, for various $`l`$, $`\mu _c/\mu _b`$.

The author is indebted to E. M. Lavely for numerous discussions. The support of SERDP, through contract No. DACA 72-02-C-0029, is gratefully acknowledged.


text:
# Exponential distribution of financial returns at mesoscopic time lags: a new stylized fact

## 1 Introduction

The empirical probability distribution functions (EDFs) for different assets have been extensively studied by the econophysics community in recent years . Stock and stock-index returns have received special attention. We focus here on the EDFs of the returns of individual large American companies from 1993 to 1999, a period without major market disturbances. By ‘return’ we always mean ‘log-return’, the difference of the logarithms of prices at two times separated by a time lag $`t`$.

The time lag $`t`$ is an important parameter: the EDFs evolve with this parameter. At micro lags (typically shorter than one hour), effects such as the discreteness of prices and transaction times, correlations between successive transactions, and fluctuations in trading rates become important. Power-law tails of EDFs in this regime have been much discussed in the literature before . At ‘meso’ time lags (typically from an hour to a month), continuum approximations can be made, and some sort of diffusion process is plausible, eventually leading to a normal Gaussian distribution. On the other hand, at ‘macro’ time lags, the changes in the mean market drifts and macroeconomic ‘convection’ effects can become important, so simple results are less likely to be obtained. The boundaries between these domains to an extent depend on the stock, the market where it is traded, and the epoch. The micro-meso boundary can be defined as the time lag above which power-law tails constitute a very small part of the EDF. The meso-macro boundary is more tentative, since statistical data at long time lags become sparse.

The first result is that we extend to meso time lags a stylized fact known since the 19th century  (quoted in ): with a careful definition of time lag $`t`$, the variance of returns is proportional to $`t`$.

The second result is that log-linear plots of the EDFs show prominent straight-line (tent-shape) character, i.e. the bulk (about 99%) of the probability distribution of log-return follows an exponential law. The exponential law applies to the central part of EDFs, i.e. not too big log-returns. For the far tails of EDFs, usually associated with power laws at micro time lags, we do not have enough statistically reliable data points at meso lags to make a definite conclusion. Exponential distributions have been reported for some world markets  and briefly mentioned in the book  (see Fig. 2.12). However, the exponential law has not yet achieved the status of a stylized fact. Perhaps this is because influential work  has been interpreted as finding that the individual returns of all the major US stocks for micro to macro time lags have the same power law EDFs, if they are rescaled by the volatility.

The Heston model is a plausible diffusion model with stochastic volatility, which reproduces the timelag-variance proportionality and the crossover from exponential distribution to Gaussian. This model was first introduced by Heston, who studied option prices . Later Drăgulescu and Yakovenko (DY) derived a convenient closed-form expression for the probability distribution of returns in this model and applied it to stock indexes from 1 day to 1 year . The third result is that the DY formula with three lag-independent parameters reasonably fits the time evolution of EDFs at meso lags.

## 2 Probability distribution of log-returns in the Heston model

In this section, the Heston model  and the DY formula  are briefly summarized. The price $`S_t`$ of a model stock obeys the stochastic differential equation of multiplicative Brownian motion: $`dS_t=\mu S_tdt+\sqrt{v_t}S_tdW_t^{(1)}`$. Here the subscript $`t`$ indicates time dependence, $`\mu `$ is the drift parameter, $`W_t^{(1)}`$ is a standard random Wiener process, and $`v_t`$ is the fluctuating time-dependent variance. The detrended log-return is defined as $`x_t=\mathrm{ln}(S_t/S_0)\mu t`$, although detrending is a minor correction at meso lags. In the Heston model, the variance $`v_t`$ obeys the mean-reverting stochastic differential equation:

$$dv_t=\gamma (v_t\theta )dt+\kappa \sqrt{v_t}dW_t^{(2)}.$$
(1)
Here $`\theta `$ is the long-time mean of $`v`$, $`\gamma `$ is the rate of relaxation to this mean, and $`\kappa `$ is the variance noise. We take the Wiener processes $`W_t^{(1,2)}`$ to be uncorrelated.

The DY formula  for the probability density function (PDF) $`P_t(x)`$ is:

$`P_t(x)={\displaystyle \underset{\mathrm{}}{\overset{+\mathrm{}}{}}}{\displaystyle \frac{dk}{2\pi }}e^{ikx+F_{\stackrel{~}{t}}(k)},F_{\stackrel{~}{t}}(k)={\displaystyle \frac{\alpha \stackrel{~}{t}}{2}}\alpha \mathrm{ln}\left[\mathrm{cosh}{\displaystyle \frac{\mathrm{\Omega }\stackrel{~}{t}}{2}}+{\displaystyle \frac{\mathrm{\Omega }^2+1}{2\mathrm{\Omega }}}\mathrm{sinh}{\displaystyle \frac{\mathrm{\Omega }\stackrel{~}{t}}{2}}\right],`$ (2)
$`\stackrel{~}{t}=\gamma t,\alpha =2\gamma \theta /\kappa ^2,\mathrm{\Omega }=\sqrt{1+(k\kappa /\gamma )^2},\sigma _t^2x_t^2=\theta t.`$ (3)
The variance $`\sigma _t^2x_t^2`$ (3) of the PDF (2) increases linearly in time, while $`x_t=0`$. The three parameters of the model are $`\gamma `$, $`\theta `$ and $`\alpha `$. At short and long time lags, the PDF (2) reduces to exponential (if $`\alpha =1`$) and Gaussian :

$$P_t(x)\{\begin{array}{cc}\mathrm{exp}(|x|\sqrt{2/\theta t}),\hfill & \stackrel{~}{t}=\gamma t1,\hfill \\ \mathrm{exp}(x^2/2\theta t),\hfill & \stackrel{~}{t}=\gamma t1.\hfill \end{array}$$
(4)
In both limits, it scales with the volatility: $`P_t(x)=f(x/\sqrt{x_t^2})=f(x/\sqrt{\theta t})`$, where $`f`$ is the exponential or the Gaussian function.

## 3 Data analysis and discussion

We analyzed the data from Jan/1993 to Jan/2000 for $`27`$ Dow companies, but show results only for four large cap companies: Intel (INTC) and Microsoft (MSFT) traded at NASDAQ, and IBM and Merck (MRK) traded at NYSE. We use two databases, TAQ to construct the intraday returns and Yahoo database for the interday returns. The intraday time lags were chosen at multiples of 5 minutes, which divide exactly the 6.5 hours (390 minutes) of the trading day. The interday returns are as described in  for time lags from 1 day to 1 month = 20 trading days.

In order to connect the interday and intraday data, we have to introduce an effective overnight time lag $`T_n`$. Without this correction, the open-to-close and close-to-close variances would have a discontinuous jump at 1 day, as shown in the inset of the left panel of Fig. 1. By taking the open-to-close time to be 6.5 hours, and the close-to-close time to be 6.5 hours + $`T_n`$, we find that variance $`x_t^2`$ is proportional to time $`t`$, as shown in the left panel of Fig. 1. The slope gives us the Heston parameter $`\theta `$ in Eq. (3). $`T_n`$ is about 2 hours (see Table 1).

In the right panel of Fig. 1, we show the log-linear plots of the cumulative distribution functions (CDFs) vs. normalized return $`x/\sqrt{\theta t}`$. The $`\text{CDF}_t(x)`$ is defined as $`_{\mathrm{}}^xP_t(x^{})𝑑x^{}`$, and we show $`\text{CDF}_t(x)`$ for $`x<0`$ and $`1\text{CDF}_t(x)`$ for $`x>0`$. We observe that CDFs for different time lags $`t`$ collapse on a single straight line without any further fitting (the parameter $`\theta `$ is taken from the fit in the left panel). More than 99% of the probability in the central part of the tent-shape distribution function is well described by the exponential function. Moreover, the collapsed CDF curves agree with the DY formula (4) $`P_t(x)\mathrm{exp}(|x|\sqrt{2/\theta t})`$ in the short-time limit for $`\alpha =1`$ , which is shown by the dashed lines.

Because the parameter $`\gamma `$ drops out of the asymptotic Eq. (4), it can be determined only from the crossover regime between short and long times, which is illustrated in the left panel of Fig. 2. We determine $`\gamma `$ by fitting the characteristic function $`\stackrel{~}{P}_t(k)`$, a Fourier transform of $`P_t(x)`$ with respect to $`x`$. The theoretical characteristic function of the Heston model is $`\stackrel{~}{P}_t(k)=e^{F_{\stackrel{~}{t}}(k)}`$ (2). The empirical characteristic functions (ECFs) can be constructed from the data series by taking the sum $`\stackrel{~}{P}_t(k)=\mathrm{Re}_{x_t}\mathrm{exp}(ikx_t)`$ over all returns $`x_t`$ for a given $`t`$ . Fits of ECFs to the DY formula (2) are shown in the right panel of Fig. 2. The parameters determined from the fits are given in Table 1.

In the left panel of Fig. 3 we compare the empirical PDF $`P_t(x)`$ with the DY formula (2). The agreement is quite good, except for the very short time lag of 5 minutes, where the tails are visibly fatter than exponential. In order to make a more detailed comparison, we show the empirical CDFs (points) with the theoretical DY formula (lines) in the right panel of Fig. 3. We see that, for micro time lags of the order of 5 minutes, the power-law tails are significant. However, for meso time lags, the CDFs fall onto straight lines in the log-linear plot, indicating exponential law. For even longer time lags, they evolve into the Gaussian distribution in agreement with the DY formula (2) for the Heston model. To illustrate the point further, we compare empirical and theoretical data for several other companies in Fig. 4.

In the empirical CDF plots, we actually show the ranking plots of log-returns $`x_t`$ for a given $`t`$. So, each point in the plot represents a single instance of price change. Thus, the last one or two dozens of the points at the far tail of each plot constitute a statistically small group and show large amount of noise. Statistically reliable conclusions can be made only about the central part of the distribution, where the points are dense, but not about the far tails.

## 4 Conclusions

We have shown that in the mesoscopic range of time lags, the probability distribution of financial returns interpolates between exponential and Gaussian law. The time range where the distribution is exponential depends on a particular company, but it is typically between an hour and few days. Similar exponential distributions have been reported for the Indian , Japanese , German , and Brazilian markets , as well as for the US market  (see also Fig. 2.12 in ). The DY formula  for the Heston model  captures the main features of the probability distribution of returns from an hour to a month with a single set of parameters. We believe that econophysicists should be aware of the presence of the exponential distribution and recognize it as another “stylized fact” in the set of analytical tools for financial data analysis.

We thank Chuck Lahaie from the Robert H. Smith School of Business at UMD for help with the TAQ database.


text:
# Stylized facts in minority games with memory: a new challenge

## 1 Introduction

Remarkable new phenomenology arises when agents are granted the possibility of not participating in a game . In the context of the Minority Game (MG) , this extension – which takes the name of Grand Canonical MG (GCMG) – has been introduced independently by several authors . The MG captures some aspects of financial market dynamics, namely the interplay between exploitable predictability, price fluctuation and trader behaviour (we refer the interested reader to  for more details). Furthermore, GCMGs reproduce some of the main stylized facts of financial market phenomenology, such as large price changes and volatility clustering . Remarkably, these features emerge close to a phase transition, thus suggesting that financial markets operate in the vicinity of a critical point .

For physicists, this is almost intuitive. Nevertheless, it is hard to convince economists that it is actually true. The MG is able to propose a coherent account of why financial markets are at criticality, based on market (in)efficiency (or predictability) . The idea of explaining anomalous fluctuations as a critical phenomenon is not new. Other models, notably those based on two-dimensional percolation , also relate stylized facts with criticality. There is however little reason to consider the market as bi-dimensional (although the trading floors are), or to believe it should spontaneously reach the critical percolation probability. The fact that all kinds of financial markets generate power-law-distributed price changes implies that a reasonable model of a financial market should not need the fine-tuning of a parameter.

An annoying feature of the GCMG is that, close to the phase transition, its dynamics displays a peculiar dependence on initial conditions: anomalous fluctuations only materialize in some realizations of the dynamics, whereas others exhibit plain Gaussian price fluctuations . The introduction of a finite memory in the learning dynamics of agents, such that the impact that events occurred far in the past have on the agents’ choices in the present vanishes with the time lag, is the natural remedy to this inconvenient. This modification was introduced in  for a model in which agents play the MG strategically. Here we focus on the standard MG, where agents adopt a naïve price-taking behavior .

We shall first show that indeed a finite memory eliminates the dependence on initial conditions. A direct extrapolation of theoretical results  suggests that a finite memory leads to an increase of fluctuations. This is confirmed by numerical simulations. Then we shall see that, however, memory brings about new non-trivial effects in the GCMG. In particular, it causes a build-up of fluctuations in the stationary state that sets in after very large times. Our discussion will be mainly informal and we shall give only the minimal technical details. We will refer to original papers for the detailed definition of the models.

## 2 Finite memory score in the standard MG

In few words, the MG models a system of $`N`$ interacting agents who should repeatedly take a binary decision, such as buy/sell. Each agent aims at taking the minority decision. This captures the fact that it is generally convenient to buy when the majority sells and vice-versa. Agents resort to trading strategies to process market information. The latter takes one of $`P=\alpha N`$ possible values. The performance of each strategy is monitored by a score function. The dynamics of the score of strategy $`s`$ has the form

$$U_{i,s}(t+1)=U_{i,s}(t)a_{i,s}^{\mu (t)}A(t)$$
(1)
where $`a_{i,s}^{\mu (t)}=\pm 1`$ is the action prescribed by strategy $`s`$ when the information takes the value $`\mu (t)`$ and $`A(t)`$ is the sum of all actions ($`\pm 1`$) of agents. At each time step, agents adopt the strategy with the highest score. The last term in (1) is the payoff of agent $`i`$ for using strategy $`s`$ at time $`t`$. As the particular form of the dynamics is not at stake in what follows, we refer to  for more details. We just remind that the key observables are the magnitude of fluctuations $`\sigma ^2=A^2`$ and the information content $`H=(1/P)_\mu A|\mu ^2`$, which details how the outcome $`A(t)`$ is correlated with the information $`\mu (t)`$ . We also remind that the value $`N=P/\alpha _c`$, with $`\alpha _c=0.3374\mathrm{}`$, separates an asymmetric, information rich phase with $`H>0`$ (for $`N<P/\alpha _c`$) from a symmetric phase with $`H=0`$ ($`N>P/\alpha _c`$). In the latter, the stationary state depends on the asymmetry

$$U_0=U_{i,1}(0)U_{i,2}(0)$$

of the initial conditions (we focus on the MG with two strategies per agent). This dependence is due to the fact that, under (1), all market fluctuations since time $`t=0`$ are remembered and contribute with the same weight to $`U_{i,s}(t)`$, irrespective of how far in the past they took place.

This suggests to generalize (1) to

$$U_{i,s}(t+1)=(1\lambda /P)U_{i,s}(t)a_{i,s}^{\mu (t)}A(t).$$
(2)
with $`\lambda >0`$ a constrant. Fig. 1 shows that a time dependent $`\sigma ^2=A^2_t`$ 111Here $`\mathrm{}_t`$ stands for an average over a long but finite time interval around $`t`$. converges, for long times, to a value which is independent of initial conditions $`U_0`$. This implies that the initial asymmetry $`U_0`$ has no influence on the stationary state. This is confirmed by Fig. 2, where we plot $`\sigma ^2/N`$ as a function of $`\alpha =1/n_s`$ for different $`U_0`$ and $`\lambda `$. For a fixed $`\alpha `$ and for a given realization of the game, $`\sigma ^2`$ and $`H`$ are increasing functions of $`\lambda `$, although in the symmetric phase ($`\alpha <\alpha _c=0.3374\mathrm{}`$) $`\sigma ^2`$ varies very slowly with $`\lambda `$ (see Fig. 3). While only small values of $`\lambda `$ make physical sense, for completeness the same figure reports also the queer effects of very large $`\lambda `$. From (2), $`\lambda =P`$ cancels the contribution of $`U_i(t)`$ in $`U_i(t+1)`$; larger $`\lambda `$ implies that this contribution is of opposite sign to the usual MG.

## 3 Finite memory in the Grand Canonical MG

Let us focus now on the simplest GCMG , where the agents (speculators) possess a single strategy each and may abstain from playing if the game is not profitable enough. A measure of profitability is given by a parameter $`ϵ`$, which measures the gains from investing outside of the market. Again, we refer to the original paper  for a detailed account. The important new feature of this model is that the number of participants is not fixed a priori, but is dynamically determined by market gains. The model, which includes a number $`N_p`$ of producers who trade no matter what, according to a given strategy, can be solved exactly in the limit of infinitely many speculators $`N_s\mathrm{}`$ (with $`\alpha =P/N_s=1/n_s`$ and $`n_p=N_p/P`$ finite) via a static replica calculation (this solution can be also recovered with the dynamical generating functional technique ). This allows for a systematic understanding of the model’s properties. It turns out that a phase transition occurs for $`ϵ=0`$ at $`n_s=1/\alpha >n_s^c`$. The transition is discontinuous across the $`ϵ=0`$ line, where the number of active agents jumps abruptly. Stylized facts are observed, for $`ϵ>0`$, close to the point $`ϵ=0`$, $`n_sn_s^{}(P)>n_s^c`$ in finite systems,222$`n_s^{}(P)`$ increases with $`P`$ . and the size of the critical region can also be characterized. More precisely, we observe the occurrence of power-law distribution of returns $`A(t)`$ and volatility clustering , or long-time correlations333Note that other functions for describing the long-range correlation have been proposed  $`|A(t)||A(t+\tau )|\tau ^\gamma `$. These findings are summarized in Fig. 4.

First we see that fat tails or volatility clustering both emerge close to the phase transition $`n_sn_s^{}(P)`$ and $`ϵ0`$. Then, these two features do not appear in all runs. Some runs, even in the critical region, show plain Gaussian returns, sometimes with long range correlations and sometimes without. In the real world, this would imply that some markets show fat-tailed price changes, and some others not, which is at odds with the observed universality.

In the GCMG, this issue signals a non-ergodic behavior and is related to the dependence on initial conditions. Both effects arise because of the presence of infinite memory in the strategy scores dynamics. This suggests to introduce a finite memory $`\lambda >0`$ in the dynamics444A second solution may be to allow for the evolution of agent’s strategies. In that case, as observed in , an infinite memory would make no sense., as done in the previous section for the standard MG (see (1) and (2)). We remark that, technically, finite score memory can also be implemented by fixing a time window $`T`$ such that only the $`T`$ most recent scores are kept by the agents, as in 555Such games with binary payoffs  can be written as a matrix; predicting large price changes amounts to studying its eigenvalues .. The two procedures are equivalent for $`TP/\lambda `$.

The dynamics of the GCMG with finite memory is surprising. At first, fluctuations are reduced and are Gaussian. Many iterations are needed before stylized facts emerge again. Remarkably, the typical time for the emergence of anomalous fluctuations is much larger than characteristic times, which are of order $`P/\lambda `$. When one moves deeper into the crowded market region $`N_s/Pn_s^c`$, one sees even stronger effects, as shown in Fig. 6. Still, one has to wait a quite long time before wild fluctuations set in. But now large price fluctuations are proportional to the number of speculators ($`A(t)N_s`$) and they appear associated to a particular information pattern $`\mu `$. The sample of Fig. 6 shows that the value of $`\mu `$ for which $`A(t)N_s`$ may even change in the course of time. What is striking is that a modification like memory, which naively should make the dynamics smoother by introducing a finite cutoff in time dependencies, actually provokes the boost of wild fluctuations. The fact that this effects sets in for times much longer than the memory decay strongly suggests that its origin lies in collective fluctuation phenomena. Finally, as for the plain MG, a finite memory $`\lambda >0`$ eliminates all dependences on initial conditions and non-ergodicity effects.

## 4 Conclusions

In short, we have shown that finite score memory suppresses non-ergodic behavior and dependence on initial conditions in MGs. At the same time, it brings about new spectacular fluctuation phenomena in the GCMG. Interestingly, all the analytical tools used so far for the study of MG fail when $`\lambda >0`$. On one side, as soon as $`\lambda >0`$ the stationary state can no longer be related to the minimum of a global function, in contrast with the $`\lambda =0`$ case . This seems to rule out replica approaches. On the other, a finite memory destroys the frozen component of the dynamical variables which allows one to extract information from the generating functional approach . This makes finite memory MG a quite challenging problem.


text:
# A method for direct observation of quantum tunneling in a single molecule.

## Abstract

An application of impedance measurement technique (IMT) for a detection of quantum tunneling in molecular structures is investigated. A charged particle which tunnels in a two-well potential is electrically coupled to a high-quality superconducting LC-circuit(tank) that makes possible a measurement of the electric susceptibility of the molecule at the resonant frequency of the tank. The real part of this susceptibility bears information about the tunneling rate through a measurable parameter - a phase angle between the tank voltage and a bias current applied to the tank. It is shown that the present approach is highly sensitive and allows the monitoring of the tunnel motion of charged nuclei in a single molecule.

Tunneling of a particle between minima of a double-well potential represents a pure quantum phenomenon which takes place in many physical and chemical systems Bender1 ; McMahon1 . Splittings of molecular energy levels due to tunneling motion of constituent nuclei have been observed using spectroscopical methods Oppen1 ; Frost1 ; Sekiya1  as well as methods of nuclear magnetic resonance relaxometry Horsewill1 . The measurements have been performed with samples containing a large number of molecules. Recent progress in molecular manipulation with a scanning tunneling microscope tip Wolkow1 ; Hla1  makes possible a creation of single-molecule devices, and, in particular, apparatuses which use tunneling effects QC-NH3 .

Here we propose to employ an impedance measurement technique Greenberg2002 ; APL\_IMT  for contactless characterization of quantum tunneling in single molecules. This technique has provided a simple and reliable method for experimental investigation of tunneling in magnetic systems, especially, in superconducting flux qubits Izmalkov1 ; Grajcar1 ; Izmalkov2 . In the framework of IMT method the qubit’s loop is inductively coupled to a high-quality LC circuit (tank) driven by a time-dependent bias current. An interaction of the tank with the qubit modifies an effective inductance of the tank and produces a shift of its resonant frequency. As a result, a voltage in the LC-circuit, $`V_T(t)=V_T\mathrm{cos}(\omega t+\mathrm{\Theta }),`$ has been displaced in phase from the bias current, $`I_{bias}(t)=I_{ac}\mathrm{cos}\omega t,`$ by the angle $`\mathrm{\Theta }.`$ A dependence of the angle $`\mathrm{\Theta }`$ on the bias applied to the qubit has demonstrated IMT dips which are indicative of quantum tunneling in the qubit. A tunneling rate between wells is determined by the width of the dips. From the theoretical point of view the tangent of the difference between phases of voltage and current, $`\mathrm{tan}\mathrm{\Theta },`$ is proportional to the real part of a magnetic susceptibility of the qubit.

To study intra-molecular tunneling of a charged particle, say, a proton or another light nucleus, in a two-well potential $`U(z)`$, we have to measure an electric susceptibility of such a molecule. To do this, the molecule is placed between plates of the tank’s capacitor $`C_T,`$ so that the charged particle tunnels along lines of the intra-capacitor electric field $`E_z=V_T/d.`$ Here $`V_T`$ is a voltage applied to capacitor’s plates which are separated by the distance $`d`$. We assume that the two-well potential $`U(z)`$ has minima at $`z=\pm z_0.`$ Taking into account the lowest tunneling doublet of energy eigenstates we describe quantum dynamics of the system by the Hamiltonian

$$H=\frac{\mathrm{\Delta }}{2}\sigma _x+\frac{\epsilon }{2}\sigma _z(\lambda V_T+Q_0+f)\sigma _z+H_T.$$
(1)
Here $`\mathrm{\Delta }`$ is a rate related to the particle’s tunneling between the wells, $`\epsilon `$ is a bias describing an asymmetry of the wells. This bias can be changed by applying an additional constant voltage to the capacitor $`C_T.`$ A position of the particle is characterized by the Pauli matrix $`\sigma _z:\widehat{z}=z_0\sigma _z,`$ so that an interaction with the intra-capacitor electric field $`E_z`$ is given by the term $`\lambda V_T\sigma _z,`$ where $`\lambda =e(z_0/d).`$ We add a qubit’s internal heat bath with a variable $`Q_0`$ and an auxiliary external force $`f`$. It should be noted that thermal fluctuations of the tank voltage contribute to decoherence and relaxation of the qubit together with internal mechanisms. We consider the tank, having a resonant frequency $`\omega _T=1/\sqrt{L_TC_T},`$ as a quantum harmonic oscillator with the Hamiltonian ($`\mathrm{}=1`$):

$$H_T=\omega _T(a^+a+1/2)(a+a^+)Q_bL_T\widehat{I}_TI_{bias}+H_{TB}.$$
(2)
The tank’s excitations are described by creation/annihilation operators $`a^+,a([a,a^+]_{}=1);`$ in so doing for the operators of voltage and current we obtain the expressions: $`\widehat{V}_T=i\sqrt{\omega _T/2C_T}(a^+a),\widehat{I}_T=\sqrt{\omega _T/2L_T}(a^++a).`$ A linewidth broadening $`\gamma _T`$ and a finite quality factor of the tank, $`Q_T=\omega _T/(2\gamma _T),`$ are caused by the tank’s own heat bath, which is characterized by the variable $`Q_b`$ and the Hamiltonian $`H_{TB}`$.

Following to Ref.15 we derive the equation for the averaged voltage in the tank coupled to the double dot:

$$\left(\frac{d^2}{dt^2}+\gamma _T\frac{d}{dt}+\omega _T^2\right)\widehat{V}_T=\lambda \frac{\omega _T^2}{C_T}\sigma _z+\frac{1}{C_T}\dot{I}_{bias}.$$
(3)
The particle’s variable $`\sigma _z`$ is functionally-dependent on the tank’s voltage:

$$\sigma _z(t)=\lambda 𝑑t_1\frac{\delta \sigma _z(t)}{\delta f(t_1)}V_T(t_1),$$
(4)
where the functional derivative corresponds to the electric susceptibility $`\chi _z(\omega )`$ of the particle in the double-dot system

$$\frac{\delta \sigma _z(t)}{\delta f(t_1)}=\frac{d\omega }{2\pi }e^{i\omega (tt_1)}\chi _z(\omega ).$$
(5)
The bias current, $`I_{bias}(t)=I_{ac}\mathrm{cos}\omega t,`$ applied to the tank, induces oscillations of the tank’s voltage, $`V_T(t)=V_T\mathrm{cos}(\omega t+\mathrm{\Theta }),`$ with the amplitude $`V_T`$ and the phase $`\mathrm{\Theta }`$ which obey the equation

$$\left[\omega ^2i\gamma _T\omega +\omega _T^2\frac{\lambda ^2\omega _T^2}{C_T}\chi _z(\omega )\right]V_Te^{i\mathrm{\Theta }}=\frac{i\omega }{C_T}I_{ac}.$$
(6)
For the case when the frequency of the bias current is exactly equal to the initial resonant frequency of the tank, $`\omega =\omega _T,`$ the current-voltage phase shift is given by the equation

$$\mathrm{tan}\mathrm{\Theta }=2\lambda ^2\frac{\overline{Q}_T}{C_T}\chi _z^{}(\omega _T),$$
(7)
where $`\overline{Q}_T=\omega _T/(2\overline{\gamma }_T)`$ is an effective quality factor of the tank, modified because of the double dot contribution to the broadening of the tank’s linewidth: $`\overline{\gamma }_T=\gamma _T+(\lambda ^2\omega _T/C_T)\chi _z^{\prime \prime }(\omega _T).`$

The electric susceptibility of the double dot, $`\chi _z(\omega ),`$ can be found in the framework of the theory of open quantum systems ES1981 ; ASRabi2003 . Quantum dynamics of the charged particle in two-well potential is governed by the Heisenberg equations:

$`\dot{\sigma }_x=\epsilon \sigma _y+2(Q+f)\sigma _y,`$
$`\dot{\sigma }_y=\mathrm{\Delta }\sigma _z+\epsilon \sigma _x2(Q+f)\sigma _x,`$
$`\dot{\sigma }_z=\mathrm{\Delta }\sigma _y.`$ (8)
Here $`Q=Q_0+\lambda \stackrel{~}{V}_T`$ is the operator of the total dissipative environment which surrounds the double dot. For the susceptibility of this environment we have the relation: $`\chi (\omega )=\chi _0(\omega )+\chi _T(\omega ),`$ where $`\chi _0(\omega )`$ corresponds to the internal mechanisms of dissipation in the double dot, and

$$\chi _T(\omega )=\frac{e^2}{C_T}\left(\frac{z_0}{d}\right)^2\frac{\omega _T^2}{\omega _T^2\omega ^2i\omega \gamma _T}$$
(9)
describes dissipative properties of the tank which are due to its coupling to the bath $`Q_b.`$ The spectral functions of the heat bath fluctuations, $`S(\omega )=S_0(\omega )+S_T(\omega ),`$ is related to the imaginary part of the heat bath susceptibility by the fluctuation-dissipation theorem: $`S(\omega )=\chi ^{\prime \prime }(\omega )\mathrm{coth}(\omega /2T),`$ with $`T`$ being a temperature of the bath. Averaging the Heisenberg equations (A method for direct observation of quantum tunneling in a single molecule.) over the equilibrium fluctuations of environment followed by an application of Bloch-Redfield approximation allows us to find the electric susceptibility of the charged particle confined in the two-well potential

$$\chi _z(\omega )=2\mathrm{\Delta }\frac{i\omega +S(\omega +\omega _0)+S(\omega \omega _0)}{(i\omega +\gamma _0)[i(\omega \omega _0)+\gamma ][i(\omega +\omega _0)+\gamma ]}\sigma _x^0.$$
(10)
Here $`\omega _0=\sqrt{\mathrm{\Delta }^2+\epsilon ^2}`$ is the ground-state splitting for the tunneling particle, $`\sigma _x^0=(\mathrm{\Delta }/\omega _0)\mathrm{tanh}(\omega _0/2T)`$ is a steady-state value of the matrix $`\sigma _x.`$ Coupling of the double dot to the dissipative environment results in the frequency-dependent decoherence and relaxation rates $`\gamma `$ and $`\gamma _0`$

$`\gamma (\omega )={\displaystyle \frac{\mathrm{\Delta }^2}{\omega _0^2}}S(\omega )+{\displaystyle \frac{\epsilon ^2}{\omega _0^2}}\left(1{\displaystyle \frac{\omega _0}{\omega }}\right)S(\omega +\omega _0)+{\displaystyle \frac{\epsilon ^2}{\omega _0^2}}\left(1+{\displaystyle \frac{\omega _0}{\omega }}\right)S(\omega \omega _0),`$
$`\gamma _0(\omega )={\displaystyle \frac{\mathrm{\Delta }^2}{\omega _0^2}}\left[S(\omega +\omega _0)+S(\omega \omega _0)\right].`$ (11)
In particular, a back-action of the tank on the tunneling particle leads to the small decoherence and relaxation rates

$`\gamma (\omega _0)={\displaystyle \frac{1}{2}}\gamma _0(0)=\left({\displaystyle \frac{\mathrm{\Delta }}{\omega _0}}\right)^2S_T(\omega _0)={\displaystyle \frac{e^2}{2Q_TC_T}}\left({\displaystyle \frac{z_0}{d}}\right)^2\left({\displaystyle \frac{\mathrm{\Delta }}{\omega _0}}\right)^2\left({\displaystyle \frac{\omega _T}{\omega _0}}\right)^3.`$ (12)
It should be emphasized that the current-voltage angle $`\mathrm{\Theta }`$ (7) is proportional to the double-dot susceptibility, $`\chi _z^{}(\omega _T),`$ taken at the resonant frequency of the tank $`\omega _T,`$ which is much lower than the tunneling frequency $`\omega _0.`$ Taking this fact into account we find the real part of the electric susceptibility: $`\chi _z^{}(\omega _T)=2\mathrm{\Delta }\sigma _x^0/\omega _0^2.`$ As a result, the phase angle between the voltage in the tank and the bias current is determined by the following equation

$$\mathrm{tan}\mathrm{\Theta }=\frac{e^2}{C_T}\left(\frac{2z_0}{d}\right)^2Q_T\frac{\mathrm{\Delta }^2}{(\mathrm{\Delta }^2+\epsilon ^2)^{3/2}}\mathrm{tanh}\left(\frac{\sqrt{\mathrm{\Delta }^2+\epsilon ^2}}{2T}\right).$$
(13)
By way of illustration, we consider tunneling in an ammonia molecule NH3 which is described by the two-well potential energy FLP . The existence of two delocalized states with eigenenergies separated by a tunneling splitting $`\mathrm{\Delta }/h=23.8`$ GHz has motivated to harness this molecule as a quantum bit QC-NH3 . According to this proposal the NH3-molecule can be encapsulated into a fullerene $`C_{60}`$ to be addressed individually. A direct chemisorption of the ammonia molecule on the surface of a semiconductor represents another possiblity Brown1 . The fullerene ball incorporating the single ammonia molecule is inserted between the capacitor’s plates having an area about 1 $`\mu m^2`$. If the distance between plates is of order of a diameter of $`C_{60}`$molecule ParkC60 : $`d=1.2\times 10^9m`$ , then for the capacitance $`C_T`$ we get an estimation: $`C_T`$ 10 fF. Together with a value for the tank’s inductance, $`L_T`$ 1 $`\mu H,`$ it gives the resonant frequency of the tank: $`\omega _T/2\pi =`$ 1.6 GHz which is much less than the tunneling frequency $`\mathrm{\Delta }/h:\mathrm{}\omega _T/\mathrm{\Delta }=0.06.`$ The electric dipole moment of the ammonia molecule Brown1 , $`2ez_0,`$ which is involved in the formula (13), is of order of 0.31 $`e\times 10^{10}`$ m . For the case when losses in the tank are mainly due to a shunt resistance $`R_T,`$ which is parallel to the capacitor $`C_T`$, the linewidth of the tank is defined as: $`\gamma _T=1/(R_TC_T),`$ so that the quality factor is $`Q_T=\omega _TR_TC_T/2,`$ and the current-voltage angle $`\mathrm{\Theta }`$ does not depend on the capacitance of the tank $`C_T:`$

$$\mathrm{tan}\mathrm{\Theta }=\frac{\mathrm{}\omega _T}{2\mathrm{\Delta }}\frac{R_T}{(\mathrm{}/e^2)}\left(\frac{2z_0}{d}\right)^2\frac{\mathrm{\Delta }^3}{(\mathrm{\Delta }^2+\epsilon ^2)^{3/2}}\mathrm{tanh}\left(\frac{\sqrt{\mathrm{\Delta }^2+\epsilon ^2}}{2T}\right).$$
(14)
The current-voltage angle $`\mathrm{\Theta }`$ can be measured in the experiment as a function of the bias $`\epsilon `$ applied to the tunneling particle. With a reasonable estimation for the shunt resistance ParkC60 : $`R_T=50\times 10^6`$ Ohm, the single ammonia molecule demonstrates a pronounced low-frequency response: $`\mathrm{tan}\mathrm{\Theta }=0.23,`$ at temperatures $`T<\mathrm{\Delta }/k_B=`$ 1.1 K, and $`\epsilon =0.`$ Measurement-induced decoherence, $`\gamma =\gamma (\mathrm{\Delta })`$ (12), is negligibly small in the process: $`\gamma /\mathrm{\Delta }10^{11},`$ that allows us to characterize the IMT procedure as a weak continuous measurement Averin1 ; Korotkov1 . To be realistic, for the example under discussion - the ammonia molecule confined inside of $`C_{60}`$ \- we have to take into account screening effects due to a fullerene cage QC-NH3 . Because of the screening, a magnitude of the IMT dip (13,14) decreases that can be compensated by a proper increasing of the shunt resistance $`R_T`$ and the quality factor of the tank $`Q_T`$.

In conclusion, we have proposed to use the impedance measurement technique for continuous monitoring of quantum tunneling in molecular systems which are characterized by a two-well potential energy. To this end, the molecule(s) should be inserted between the plates of a capacitor $`C_T`$ comprising a high-quality superconducting tank along with an inductance $`L_T`$. We have shown that the present approach is sensitive enough to detect a tunneling motion in a single molecule.

The author is grateful to M. Grajcar, E. Il’ichev, A. Izmalkov, S. Rashkeev, M. Steininger, and A. Zagoskin for helpful comments.