text:
# Variability of the H𝛼 emission of Cygnus X-1 and its connection with the soft X-ray radiation

## 1 Introduction

Cyg X-1 has been one of the most well-studied X-ray sources since its discovery in 1965 (Bowyer et al. 1965) and identification with the O9.7Iab supergiant HD226868 (Bolton 1972; Webster & Murdin 1972). It provided the first observational evidence for the existence of black holes when it was found to be a spectroscopic binary with an orbital period of 5.6 days. The ephemeris was revised a number of times and a recent estimation based on spectroscopic and photometric orbital variability is presented by Brocksopp et al. (1999a). The masses of the two components were calculated as 17.5 $`M_{}`$ for the supergiant and 10.1 $`M_{}`$ for the black hole, with adopted system inclination $`i`$=35°(Herrero et al. 1995).

Cyg X-1 is a detached double system but the optical component is very close to filling its Roche lobe. Typically of high mass stars, HD22686 is thought to emit a stellar wind, which appears to be attracted towards the black hole as in the focused wind model of Friend & Castor (1982; Gies & Bolton 1986; Sowers et al. 1998). The optical spectrum of the system has only two well-pronounced emission lines: H$`\alpha `$ and He ii$`\lambda `$4686, but high-resolution spectroscopy detects line-asymmetry in most of the hydrogen and red He i lines (Canalizo et al. 1995) which thus cannot be used for precise determination of the orbital elements. Gies & Bolton (1986) and Ninkov et al. (1987) found that the He ii$`\lambda `$4686 emission component has a radial velocity curve that is best interpreted as emission originating between the stars. The H$`\alpha `$ profile appears to have at least two components; P Cygni emission moving with the orbital motion of the primary star and a second emisison component which follows the quite different radial velocity curve of the He ii$`\lambda `$4686 emission line and which is formed in the focused wind flow between the stars (Ninkov et al. 1987; Sowers et al. 1998).

Cyg X-1 demonstrates sporadic flaring at soft X-ray energies (2-12 keV) on different timescales, from days to months and years (Brocksopp et al. 1999b). It is believed that the source of that variability is located near the black hole companion and is connected with the accretion disc. A supergiant wind is also variable on different timescales from days to years (Kaufer et al. 1996; Rivinius et al. 1997) and the variability of the emission lines and the soft X-rays should be somehow connected because they both depend on mass transfer via the wind. Discovery of any correlation between the stellar wind and the X-ray variability is important for understanding the physical processes of mass transfer through the disc and onto the black hole.

Long-term monitoring of variability of the stellar wind from the supergiant of Cyg X-1 has not been previously attempted. Nonetheless some orbit-independent variability of the equivalent widths of the He ii$`\lambda `$4686 emission line was found by Ninkov et al. (1987) and of the emission component of H$`\alpha `$ by Voloshina et al. (1997). The latter also discovered a decrease of EW accompanying the period of high/soft state behaviour in 1996 (see also Brocksopp et al. 1999b). Thus we have monitored the H$`\alpha `$ line at the Crimean Astrophysical Observatory (CrAO) since the 1996 X-ray spectral state change, when a strong increase in soft X-rays was found by the Rossi X-ray Timing Explorer (RXTE); our observations were supported by long-term photometric monitoring provided by Sternberg Astronomical Institute of Moscow.

## 2 Spectroscopic observations and data reduction

The spectroscopic monitoring of Cyg X-1 took place over seven years from 1996 to 2002, during which time we obtained a total of 142 spectra. We used the Coudé focus of the 2.6-m telescope of the CrAO and the detector was an Electronix CCD array (1024$`\times `$260 pixels). All observations were performed in the second order of a diffraction grating with a reciprocal dispersion 3 Å/mm. The spectral resolution was about 35000 and during each observation we obtained a 60 Å spectrum. The exposure time ranged from 30 to 45 minutes, depending on weather conditions resulting in a signal-to-noise ratio of 100 or better. In most cases two or three spectra were obtained during each night and then added for improved signal-to-noise and cosmic ray subtraction.

The spectra were reduced using standard techniques for dark-current subtraction and flat-field division. The subsequent reduction included removal of spurious spikes produced by cosmic-rays, sky-background subtraction, water vapor line subtraction using non-emission fast-rotating early type stars and filtering of the spectra by the moving-average method using three points. For the wavelength calibration we used Th Ar comparison spectra obtained immediately before and after the target observation. The accuracy of the calibration was 1.5 km/s or better. All spectra were corrected to the baricentre of the solar system, both in wavelength and time of mid-exposure. The continuum normalization was provided by fitting of a fifth degree polynomial.

As a final step, equivalent widths (EW) of each spectrum were measured. Estimated errors of EW measurements were better than 0.1Å. The results are presented in Table 1.

During the period of our observations the system was monitored by RXTE. In order to compare variability of H$`\alpha `$ with the X-ray variability, we use the All Sky Monitor (ASM) public archive data from the web. The ASM observations were made in the three energy bands: A (1.5–3 keV), B (3–5 keV) and C (5–12 keV). A detailed description of the ASM, including calibration and reduction, is published in Levine et al. (1996). The daily mean data were converted from counts to energy units ($`\text{keV}\text{cm}^2\text{s}^1`$) using transformation formulae obtained from Zdziarski et al. (2002).

## 3 H$`\alpha `$ line-profile variability

Previous study of the emission spectrum of Cyg X-1 shows that the H$`\alpha `$ line-profile has strong variability on the 5.6-day orbital period (Gies & Bolton, 1986; Ninkov et al. 1987; Sowers et al. 1998). However our observations show high variability of the emission line intensity that is phase-independent and which was not noticed in the previous observations (which took place over a shorter time range). Since mass loss by the optical component is the only source of accretion, any instability of the stellar wind will lead to a variable mass accretion though the response of the accretion disc and that will have some impact on the X-ray emission.

Such complex variability makes analysis of the line-profile variations much more difficult to interpret because the shape and strength of the H$`\alpha `$ emission profile depends on a number of factors:

1. phase of the orbital period (wind is noticeably collimated in the direction of the secondary)
2. variability of the stellar wind density due to variable mass loss rate by the supergiant
3. clumpy structure of the stellar wind providing discrete absorption components (DAC)
4. time-variable and phase-locked photoionization of the stellar wind by the soft X-ray radiation from the accretion disc of the secondary component

In order to overcome this we have divided our spectra into orbital phase groups spanning 0.05 in phase and analysed each group separately according to processes that may have influenced them. Figs. 1 and 2 show three sample line-profiles for each group, chosen at random in terms of continuum or X-ray activity but so as to include extremely strong, weak and intermediate intensity examples of the H$`\alpha `$ emission line in each case (plotted spectra are indicated in Table 1). Spectra drawn with a thick line are those obtained during periods of X-ray activity (i.e. intensity of the 1.5–3 keV X-rays in the exceed 4 $`\text{keV}\text{cm}^2\text{s}^1`$) and it appears that the X-ray variabilty has little effect on the profiles (see Sect. 4).

Figs. 1 and 2 show that the orbital variability of the line profile is mostly in agreement with the focused wind model; there are two emission components as demonstrated by Sowers et al. (1998). One forms a relatively stable P Cyg structure, which is produced by the stellar wind from the hot supergiant and moves around the spectrum in the phase of the orbital period. The other emission component originates between the stars in a focused wind flow from the supergiant to the unseen companion. The orbital movement of this component of the line-profile is different since it appears in the collimated stellar wind close to the L1 point, such that it is blue-shifted during the phase range 0.10–0.60 and merged with the P Cyg emission component at other times.

Irregular variability of the line-profiles is clearly seen during all phase intervals; sometimes (phases: 0.25-0.35, 0.50-0.55, 0.70-0.75) there is an increase or decrease in the intensity of the whole line but for the majority of phase intervals irregular variation of the blue part of the profile seems to dominate. These variations can be connected with the variability of either the collimated wind component or with the absorption part of the P Cyg line-profile from the stellar wind of the supergiant. However since both emission components reflect variable mass loss from the optical star, it is most probable that they are connected with one another.

Fig. 3 shows long-term light curves for the $`RXTE`$ intensity (Band A), hardness ratio (HR1 = Bands B/A) and H$`\alpha `$ EW. As with the line profiles it is clear that there is complex variability in the EW measurements, on both short ($``$days e.g. MJD 51392–51398) and long ($``$few weeks e.g. MJD 51980–52200) timescales but the sampling is often insufficient to be conclusive. However, it is certain that this variability is not connected with the orbital period and instead seems to reflect some other mechanism(s) within the system.

## 4 Connection with the soft X-ray variability

The ASM data plotted in Fig. 3 demonstrate that after the high/soft state period (i.e. the major outburst in the softer X-ray bands) in 1996 (MJD 50100–50400) the subsequent flares increased in both intensity and softness from year to year, reaching a maximum in 2002 with another high/soft state period. These flares have a semi-regular nature with characteristic occurrence time of $``$400 days.

The variability of the EW and ASM X-ray flux are often loosely correlated. For example the most active X-ray periods MJD 50200–50350 and 52100–52500 were both accompanied by a decrease in EW; similarly the more short-lived flares during MJD 51850–51900 and 51440–51480. However a decrease in EW is not proportional to the amplitude of X-ray flares and is not always associated with episodes of X-ray activity. Furthermore the EW of H$`\alpha `$ decreased monotonically from $`1.1`$ Å on MJD 51397 to 0.3 Å on MJD 51456 (Fig. 3b) prior to the X-ray flare, reaching a minimum simultaneously with the small X-ray peak.

To understand the situation better, we combined the EW of H$`\alpha `$ with the three different ASM spectral bands. Fig. 4 shows that the relationship between the EW and the softer X-ray emission depends strongly on the nature of the spectral state; X-ray flares appear only if the EW drops below $`0.4`$ Å. The softest X-ray spectral band (A) is most sensitive to the H$`\alpha `$ variability such that a decrease in EW will probably, although not necessarily, be associated with an increase in X-ray flux as seen in the lightcurves; on the contrary the hardest X-ray band (C) appears insensitive to the EW.

The hardness of the X-ray spectrum can be a useful parameter for understanding the relation between X-ray emission and emission in the H$`\alpha `$ line – the EW and X-ray hardness ratio appear to be correlated for both X-ray colours (Fig. 5), although that of HR1 (i.e. the softer colour) is more significant. The scattering of the data in Fig. 5 is significantly larger than the estimated errors of the observations, both in EW and X-ray and thus some additional factors must play an important role in the connection between optical and X-ray variability.

Comparing Figs. 4 and 5 it is possible to conclude that the loose correlation between the H$`\alpha `$ and X-ray emission is more significant for the X-ray spectral hardness than for the intensity.

## 5 Discussion

A number of sources could be responsible for the presence of some relationship between the variability of the H$`\alpha `$ EW and low-energy X-ray radiation; these include (i) photoionization of the neutral gas by the soft X-ray radiation, (ii) X-ray flares and state transitions of the accretion disc due to variable mass transfer via the stellar wind and/or (iii) variability of X-ray radiation due to variable column density of neutral gas in the line of sight.

### 5.1 Photoionization vs. mass-loss variability

Unlike single hot supergiant stars, Cyg X-1 has an additional source of ionisation in the form of X-ray radiation from the accretion disc of the secondary. However, due to orbital motion, it is phase locked with the orbital period and can be identified relatively easily. This is illustrated in Fig. 6 where we display all spectra in the orbital phase intervals 0.95–1.10 (left) and 0.45-0.55 (right). During the first of these phase intervals the black hole is situated behind the optical component, with the red component of the H$`\alpha `$ profile formed in the part of the envelope which is directed towards the X-ray source (resulting in a superposition of the envelope and focussed wind emission components). During the second phase interval the star is behind the black hole and the red emission is no longer under the influence of the X-ray source.

Those spectra obtained during periods of X-ray activity have been indicated in Fig. 6 by an asterisk and show clearly that the effects of photoionisation by the X-ray emission do not dominate the shape of the emission line. However, comparison of the red wing of the red emission peak at different orbital phases shows that photoionisation does have some influence on the line-profile. At phase 0.0 the red wing has negligible extension and the photospheric line profile is often seen; this can be attributed to photoionization by the X-rays since the red wing of the emission forms bewteen the star and the black hole. At phase 0.5 the red wing forms on the far side of the star where it is unaffected by the X-ray radiation and can be extended up to 700 km/s. During both phase intervals plotted in Fig. 6 the aperiodic part of the wind variability has a strong influence over the shape of the H$`\alpha `$ line. Often this sporadic inhomogeneity of the wind or long-term variability of mass loss is so high that it prevails over the (focused wind model controlled) orbital line-profile variability. This can occur at any phase interval of the orbital period and independently of X-ray activity.

In the top panel of Fig. 7 we plot the EW values against orbital phase, using filled circles for the data obtained during periods of X-ray activity (intensities in the 1.5–3 keV band exceed 4 $`\text{keV}\text{cm}^2\text{s}^1`$). While the quiescent (and combined) data show only scatter, the scatter for the filled circles (i.e. epochs of X-ray flaring) is much more constrained.

The presence of the orbital variability is best seen in the EW of the red emission peak alone. The highly scattered wave in Fig. 7b (non-flaring episodes) is in full agreement with the focused wind model and increasing of EW at phase 0.6–0.9 can be seen as the collimated wind component merges with the P Cyg wind component (see Fig. 2). The photoionization shown in Fig. 6 is relatively small in amplitude and the large scattering of the points in Fig. 7b can be best-described by the short- and long-term variability of the stellar wind. Fig. 7c (flaring episodes) shows that photoionisation dominates during phases 0.6–1.3 when the X-ray source is located close to where the red emisison peak is forming. Both curves on Figs. 7b & 7c are statistycally significant – the full amplitudes of EW variability are exceed more than ten times standard deviations of the data. This is in agreement with the orbital variability of the resonance UV lines studied by Theves et al. (1980) and van Loon et al. (2001). The resonance UV lines form in a more extended region (which includes the Roche lobe of the secondary) than H$`\alpha `$ and are significantly influenced by photoionization even during quiescence.

The complex relationship between the X-rays and H$`\alpha `$ EW that was seen in Figs. 4 and 5 cannot be explained purely in terms of time-variable photoionisation. Plotting EWs of the red emission peak against the 1.5–3 keV X-rays for phase intervals $`0.0\pm 0.2`$ (filled circles) and $`0.5\pm 0.2`$ (open circles) shows a similar plot to Fig. 4 for each phase interval (Fig. 8). If photoionization dominated then we would expect to see correlated behaviour in each phase interval; instead a large scatter and independent variations are still present.

The decay of EW prior to the X-ray flare of MJD 51456 would suggest that (at least in some instances) a decrease in mass loss can trigger an X-ray flare. A reduction in mass-loss would increase the size of the Strömgren (i.e. photoionisation) zone surrounding the X-ray source (van Loon et al. 2001). This would increase the accretion rate through the disc and result in X-ray activity, thus increasing the size of the Strömgren zone further. In the case of the 2001–2002 high/soft state of Cyg X-1, once the EW return to $`1.5`$$`2.0`$ Å the Strömgren zone would be compressed to its normal size, as was the case after the 1996 high/soft state . However, this would appear not to be the case for all flares and optical observations with improved time resolution are needed to investigate this further.

### 5.2 Variability of the column density toward secondary

It is well known that the winds of early type stars are driven by resonance line scattering of radiation in the stellar envelope. This radiative driving mechanism is inherently unstable. Small perturbations of an otherwise monotonic wind velocity structure tend to grow into shocks. As a result the wind of this type of star is a mixture of rarified hot material and colder dense regions, with “cool” and “hot” gas in dynamical equilibrium. It is common in the literature to attribute wind “clumps” to dense formations with temperatures much lower than the temperatures of the hot gas. These clumps are presumably responsible for observed line-profile variability or the formations of DACs (discrete absorption components) observed in the optical and UV spectra of O stars. The estimated quantity of these clumps is of the order $`10^3`$ (Howk et al. 2002); their size hierarchy varies and appears dependant on spectral type and possible structure of local magnetic fields. The largest DACs can be seen in the line spectrum of the star for $``$few days but do not exist longer than a few rotational periods of the star (Kaper et al. 1997)

The whole Cyg X-1 system, including the black hole secondary, is surrounded by the relatively dense clumpy stellar wind. Its large amplitude of variability changes the column density of the neutral gas significantly in the direction of the secondary. Two types of such variability are known: orbital, that is connected with the variable column density due to orbital movement of the secondary in the envelope of the primary component (Brocksopp et al. 1999b; Wen et al. 1999); and dips, which usually last in the X-ray spectrum for several minutes up to 8 hours (Balucińska-Church et al. 2000). Both types of X-ray variability are seen during the low/hard state. Dips in the soft X-ray spectrum are relatively widely spread. During the dip, the spectrum is hardened and an iron K absorption edge may be seen, showing that they are due to photoelectric absorption (Kitamoto et al. 1984). Outside the dips the column density of neutral gas in the direction of secondary is $`6.21\pm 0.22\times 10^{21}\text{cm}^2`$ (Schulz et al. 2002) but can easily reach $`1\times 10^{23}\text{cm}^2`$ during dips (Baluchińska-Church et al. 1997). The size of the clumps that are responsible for short-lived dips is $`10^9`$ cm, i.e. relatively small in size. The quantity of clumps on the line of sight depends on the geometrical length of the envelope, so dips are more frequent at superior conjunction (Wen et al. 1999; Baluchińska-Church et al. 2000).

###### Acknowledgements.

The authors are grateful to the anonymous referee for helping to improve this paper. AET and CB gratefully acknowledge receipt of financial support from the Royal Society for collaborate work with the FSU countries.


format: text

text:
# Representations of Ω in number theory: finitude versus parity

## 1. Recursive Enumerability, Algorithmic Randomness and $`\mathrm{\Omega }`$

One of the most startling recent developments in the theory of computation is the discovery of the number $`\mathrm{\Omega }`$, through the subfield of algorithmic information theory. $`\mathrm{\Omega }`$ is a real number between 0 and 1 which was introduced by G. J. Chaitin  as an example of a number with two conflicting properties: it is both recursively enumerable and algorithmically random. Very roughly, this means that $`\mathrm{\Omega }`$ has a simple definition and can be computed in the limit from below, yet we can determine only finitely many of its digits with certainty—for the rest we can do no better than random.

Understanding the full importance of these properties requires some familiarity with the recursive functions—commonly presented through models of computation such as Turing machines or the lambda calculus. For the purposes of algorithmic information theory, however, it is convenient to abstract some of the details from these models and consider a programming language in which the (partial) recursive functions are represented by finite binary strings.111For more details see Chaitin . These strings are just programs for a universal Turing machine (or universal lambda expression) and they take input in the form of a binary string then output another binary string or diverge (fail to halt). For convenience, we will often consider these inputs and outputs to encode tuples of positive integers.

On top of this simplified picture of computation, we impose one restriction which is necessary for the development of algorithmic information theory (and hence $`\mathrm{\Omega }`$). The set of strings that encode the recursive functions must be prefix-free. This means that no program can be an extension of another, and thus each program is said to be self-delimiting. As algorithmic information theory is intricately linked with communication as well as computation, this is quite a natural constraint—if you wish to use a permanent binary communication channel, then you need to know when the end of a message has been reached and this cannot be done if some messages are extensions of others.

There are many prefix-free sets that one could choose and many recursive mappings between these and the recursive functions. These different choices of ‘programming language’ lead to different values of $`\mathrm{\Omega }`$, but this does not matter much as almost all of its significant properties will remain the same regardless. However, to allow talk of $`\mathrm{\Omega }`$ as a specific real number we will use the same language as Chaitin .

Now that we have explained what we mean by a programming language, we can give a quick overview of computability in terms of programs. A program computes a set of $`n`$-tuples if, when provided with input $`x_1,\mathrm{},x_n`$, it returns 1 if this is a member of the set and 0 otherwise. A program computes an infinite sequence if, when provided with input $`n`$, it returns the value of the $`n`$-th element in the sequence. A program computes a real, $`r`$, if it computes a sequence of rationals $`\left\{r_n\right\}`$ which converges to $`r`$ and $`\left|rr_n\right|<\frac{1}{2^n}`$. These sets, sequences and reals that are computed by programs are said to be recursive.

There are also many sets, sequences and reals that cannot be computed, but can be approximated in an important way. A program semi-computes a set of $`n`$-tuples if, when provided with input $`x_1,\mathrm{},x_n`$, it returns 1 if this is a member of the set and diverges otherwise. A program semi-computes an infinite sequence of bits if, when provided with input $`n`$, it returns 1 if the $`n`$-th bit in the sequence is 1 and diverges otherwise. A program semi-computes a real, $`r`$, if, when provided with input $`n`$, it computes a rational number, $`r_n`$, where $`\left\{r_n\right\}`$ converges to $`r`$ from below. These sets, infinite bitstrings and reals that are semi-computed by programs are said to be recursively enumerable or r.e.

There is an important point that needs to be made concerning reals and their representations. Each real number between 0 and 1 has a binary expansion: a binary point followed by an infinite sequence of bits that represents the real.222For numbers that can be expressed with a representation ending in an infinite string of 0’s, there is another representation ending in an infinite sequence of 1’s, but we shall remove this ambiguity by only using representations with an infinite number of 0’s. This will not affect the important reals in this paper, $`\mathrm{\Omega }`$ and $`\tau `$, as they are irrational and thus have unique representations regardless. Throughout this paper, we shall be making considerable use of the binary expansions of real numbers so it is important to point out an oddity in the definitions above: a real is recursive if and only if its binary expansion is recursive, but a real may be r.e. even if its binary expansion is not r.e. We shall thus take care to distinguish the weaker property of being an *r.e. real* from the stronger one of being a *real whose binary expansion is r.e.*

An example of a real that is r.e. but not recursive is $`\tau `$: the real number between 0 and 1, whose $`k`$-th digit is 1 if the $`k`$-th program (in the usual lexical ordering of finite bitstrings) halts when given the empty string as input and 0 if the $`k`$-th program diverges. Equivalently:

(1.1)
$$\tau =\underset{p_n\text{ halts}}{}2^n$$

$`\tau `$ is an r.e. real because there is a computable sequence of rationals $`\left\{\tau _i\right\}`$, where

(1.2)
$$\tau _i=\underset{\stackrel{ni}{p_n\text{ halts in }i\text{ steps}}}{}2^n$$

such that $`\left\{\tau _i\right\}`$ converges to $`\tau `$ from below.

Furthermore, it is clear that the binary representation of $`\tau `$ is also r.e. because there is a program that simulates the $`k`$-th program, halting if and only if it does. This program is a slightly modified universal program that first determines the bits of the $`k`$-th program and then simulates it.

$`\tau `$ is not recursive, however, because if a program could compute it to arbitrary accuracy, it would determine whether each program halts or not when given the empty string as input. This is known as the *blank tape problem* and is easily shown to be equivalent to the more general *halting problem*—‘does a given program halt on a given input?’. The halting problem is fundamental to the theory of computation and is the most famous problem that cannot be recursively solved. $`\tau `$ merely encodes the information necessary to solve the halting problem into the binary expansion of a real number and thus provides a very simple example of a non-computable real to which we can contrast the more exotic properties possessed by $`\mathrm{\Omega }`$.

$`\mathrm{\Omega }`$ encodes the halting problem in a more subtle way: it is the *halting probability*. We could, theoretically, generate a random program one bit at a time, by flipping a fair coin and writing down a 1 when it comes up heads and a 0 for tails—stopping if we reach a valid program. The chance of generating any given $`n`$ bit program is therefore $`\frac{1}{2^n}`$. $`\mathrm{\Omega }`$ is the chance that this method of random program construction generates a program that halts. Letting $`\left|p\right|`$ represent the size of $`p`$ in bits, we can also express $`\mathrm{\Omega }`$ as

(1.3)
$$\mathrm{\Omega }=\underset{p\text{ halts}}{}2^{\left|p\right|}$$

As was the case for $`\tau `$, there is a computable sequence of rationals $`\left\{\mathrm{\Omega }_i\right\}`$, where

(1.4)
$$\mathrm{\Omega }_i=\underset{\stackrel{\left|p\right|i}{p\text{ halts in }i\text{ steps}}}{}2^{\left|p\right|}$$

which converges to $`\mathrm{\Omega }`$ from below, showing it to be an r.e. real. However, we shall see shortly that the binary representation of $`\mathrm{\Omega }`$ is *not* r.e.

A real is said to be *algorithmically random*  if and only if the ‘algorithmic complexity’ of each $`n`$-bit initial segment of its binary expansion becomes and remains arbitrarily greater than $`n`$.333This is only one of four common definitions of algorithmic randomness, however, all have been shown to be equivalent. In other words a real, $`r`$, is algorithmically random if and only if any program that has access to outside advice in the form of binary messages requires more than $`n`$ bits of advice to compute the first $`n`$ bits of $`r`$’s binary expansion (for all values of $`n`$ above some threshold).444The reason that slightly more than $`n`$ bits of advice are needed is because in algorithmic information theory the advice comes in self-delimiting messages (which are actually programs that generate the advice—like self-extracting archives) and in order to be self-delimiting, these messages need slightly more bits than they would otherwise. In general, an $`n`$ bit string requires about $`(n+\mathrm{log}n)`$ bits. Chaitin  provides further details. Thus a random real is one for which only finitely many prefixes of its binary expansion can be compressed.

It is easy to see that a random real cannot have an r.e. binary expansion. Let $`x`$ be an arbitrary real whose binary expansion is r.e. By definition, there must be a program, $`p_x`$, that takes a positive integer, $`k`$, and halts if and only if the $`k`$-th bit of $`x`$ is 1. To determine $`n`$ bits of $`x`$, we just need to know how many of these $`n`$ values of $`k`$ make $`p_x`$ halt. We could then simply run $`p_x`$ on all the values of $`k`$ and stop when this many have halted, knowing that no more will halt and thus determining the $`n`$ bits of $`x`$. Since all positive integers less than $`n`$ can be encoded in $`\mathrm{log}n`$ bits (rounding up), we only need to send a message of about $`(\mathrm{log}n+\mathrm{log}\mathrm{log}n)`$ bits. In this manner, any prefix of $`x`$ can be significantly compressed, so $`x`$ cannot be random.

Because of this, we can see that $`\tau `$ too is not random. However, Chaitin  has proven that $`\mathrm{\Omega }`$ *is* random and so cannot be compressed in this manner.555Indeed, it has since been shown through the work of R. Solovay, C. S. Calude, P. Hertling, B. Khoussainov, Y. Wang and T. A. Slaman that the only r.e. random reals are $`\mathrm{\Omega }`$’s for different programming languages. See Calude  for more details. For sufficiently high values of $`n`$, $`n`$ bits of $`\mathrm{\Omega }`$ provide $`n`$ bits of algorithmically incompressible information.

In addition to recursive incompressibility, random reals are also characterised by *recursive unpredictability* . Consider a ‘predictive’ program that takes a finite initial segment of an infinite bitstring and returns a value indicating either ‘the next bit is 1’, ‘the next bit is 0’ or ‘no prediction’. If any such program is run on all finite prefixes of the binary expansion of a random real and makes an infinite amount of predictions, the limiting relative frequency of correct predictions approaches $`\frac{1}{2}`$. In other words when any program is used to predict infinitely many bits of a random real, such as $`\mathrm{\Omega }`$, it does no better than random—even with information about all the prior bits.

The power of this unpredictability can be seen when compare the predictability of $`\tau `$. In this case, the predictive program can easily predict an infinite amount of bits with no errors. This is because infinitely many bits of $`\tau `$ are ’easy’ to compute. For example, consider the halting behaviour of Turing machines: there are infinitely many Turing machines which have no loops in their transition graphs and thus cannot possibly diverge. When the predictive program is asked to predict the $`n`$-th bit of $`\tau `$, it can just check to see if the $`n`$-th program corresponds to such a machine, returning ‘the next bit is 1’ if it does and ‘no prediction’ otherwise.666From the definition of binary programs in algorithmic information theory, there must be a recursive mapping between programs and Turing machines (or any such model).

With its inherent incompressibility and unpredictability, $`\mathrm{\Omega }`$ really does go beyond the type of uncomputability present in a more typical non-recursive real such as $`\tau `$. However, its contrasting property of being an r.e. real makes $`\mathrm{\Omega }`$ seem to be just beyond our reach. In the next section, we will introduce Diophantine equations and show how these can be used to bring uncomputability into the more classical field of number theory. Then, in Section 3, we will show two ways of using Diophantine equations to bring $`\mathrm{\Omega }`$ and randomness to number theory—Chaitin’s original method and our new technique.

## 2. Diophantine Equations and Hilbert’s Tenth Problem

A Diophantine equation is a polynomial equation in which all of the coefficients and variables take only positive integer values. Many natural phenomena with discrete quantities are modelled well by Diophantine equations and they occur frequently in number theory. It is often convenient to express a Diophantine equation with all terms on the left hand side:

(2.1)
$$D(x_1,\mathrm{},x_m)=0$$

Here $`D`$ is a polynomial of $`x_1,\mathrm{},x_m`$ in which the coefficients can take both positive and negative integer values.

The number of solutions for a Diophantine equation varies widely. For example, $`3x_1+6=0`$ has one solution, while $`x_1x_22=0`$ has two and $`x_1x_2x_2=0`$ has infinitely many. Some however, such as $`23x_1=0`$, have no solutions at all. There are many different methods for deciding whether Diophantine equations of certain forms have solutions and determining what these solutions are, but there has been a great desire for a single method that takes an arbitrary Diophantine equation and determines whether or not it has solutions. In 1900, David Hilbert  gave the problem of finding such a method as the tenth in his famous list of important problems to be addressed by mathematicians in the 20th Century. Since then, the task of finding this method has become known simply as Hilbert’s Tenth Problem.

Another area of research concerns families of Diophantine equations. A family of Diophantine equations is a relation of the form:

(2.2)
$$D(a_1,\mathrm{},a_n,x_1,\mathrm{},x_m)=0$$

in which we distinguish between two types of variable. The variables $`x_1,\mathrm{},x_m`$ are called unknowns, while $`a_1,\mathrm{},a_n`$ and called parameters. By assigning values to each of the parameters (and treating them as constants), we pick out an individual Diophantine equation from the family. For example, the family $`a_13x_1=0`$ consists of the equations: $`13x_1=0`$, $`23x_1=0`$, $`33x_1=0`$ and so on.

Each family of Diophantine equations is naturally associated with a certain set of $`n`$-tuples of positive integers, $`𝔇`$, in the following manner:

(2.3)
$$a_1,\mathrm{},a_n𝔇x_1\mathrm{}x_mD(a_1,\mathrm{},a_n,x_1,\mathrm{},x_m)=0$$

In other words, a tuple is in the set if the equation it corresponds to has a solution. Such sets are said to be Diophantine or to have a Diophantine representation. For example, the set of all multiples of 3 is Diophantine because it is represented by the family $`a_13x_1=0`$.

Over the 1950’s and 1960’s, M. Davis, H. Putnam and J. Robinson established several important results regarding which sets are Diophantine. Their key result concerned a characterisation, not of Diophantine sets, but their close relation: *exponential* Diophantine sets.

A family of exponential Diophantine equations is a relation of the form:

(2.4)
$$D(a_1,\mathrm{},a_n,x_1,\mathrm{},x_m,2^{x_1},\mathrm{},2^{x_m})=0$$

where $`D`$ is once again a polynomial, but now some of its variables are exponential functions of others. Davis, Putnam and Robinson  used this additional flexibility to show that all r.e. sets are exponential Diophantine. It had long been known that all exponential (and standard) Diophantine sets are r.e. because it is trivial to write a program that searches for a solution to a given equation and halts if and only if it finds one. Therefore, the new result meant that the exponential Diophantine sets were precisely the r.e. sets.

In 1970, Yu. Matiyasevich  completed the final step, proving that all exponential Diophantine sets are also Diophantine and thus that the Diophantine sets are exactly the r.e. sets—a result now known as the dprm Theorem.

The dprm Theorem provides an intimate link between Diophantine equations and computability, reducing the task of determining whether a set has a Diophantine representation to a matter of programming. For instance, there is a program that takes a single input $`k`$ and halts if and only if the $`k`$-th bit of $`\tau `$ is 1. Thus, the set of positive integers that includes $`k`$ if and only if the $`k`$-th program halts is an r.e. set and via the dprm Theorem, there is a family of Diophantine equations with a parameter $`k`$, that has solutions if and only if the $`k`$-th program halts.

This family of equations provides an example of uncomputability in number theory and shows that Hilbert’s Tenth Problem must be recursively undecidable because a program that finds whether arbitrary Diophantine equations have solutions could be used to determine the bits of $`\tau `$ and thus to solve the halting problem. Indeed, it was long known that the recursive undecidability of Hilbert’s Tenth Problem would follow immediately from the dprm Theorem and this was the main motivation for its proof—the Diophantine representations for all other r.e. sets being largely a bonus.

## 3. Expressing Omega Through Diophantine Equations

While the dprm Theorem demonstrates the existence of $`\tau `$ and uncomputability in number theory, it also denies the possibility of finding a similar family of Diophantine equations expressing $`\mathrm{\Omega }`$ and randomness. This is due to the fact discussed in Section 1 that, while $`\mathrm{\Omega }`$ is an r.e. real, its sequence of bits is *not* r.e. However, the dprm Theorem only prohibits a direct Diophantine representation of $`\mathrm{\Omega }`$ and says nothing about the more subtle properties of Diophantine equations in which these bits could perhaps be encoded.

Chaitin  takes such an approach. While there is no program of one variable, $`k`$, that halts if and only if the $`k`$-th bit of $`\mathrm{\Omega }`$ is 1, Chaitin provides a program, $`P`$, that takes two variables, $`k`$ and $`N`$, and computes $`\mathrm{\Omega }`$ somewhat less directly. For a given value of $`k`$, $`P`$ can be thought of as making an infinite series of ‘guesses’ as to the value of the $`k`$-th bit of $`\mathrm{\Omega }`$—when $`P`$ is run on $`k`$ and $`N`$, it gives the $`N`$-th guess as to the $`k`$-th bit of $`\mathrm{\Omega }`$. What is impressive is that $`P`$ gets infinitely many of these guesses right and only finitely many wrong.

How does $`P`$ do this? It simply computes the sequence $`\{\mathrm{\Omega }_i\}`$ discussed in Section 1 until it gets to $`\mathrm{\Omega }_N`$ and then returns the $`k`$-th bit of $`\mathrm{\Omega }_N`$. Just as $`\{\mathrm{\Omega }_i\}`$ forms a sequence of approximations to $`\mathrm{\Omega }`$, so the $`k`$-th bit of each $`\{\mathrm{\Omega }_i\}`$ forms a sequence of approximations to the $`k`$-th bit of $`\mathrm{\Omega }`$.

Consider this $`k`$-th bit of each $`\{\mathrm{\Omega }_i\}`$ as $`i`$ is increased. This bit could change between 0 and 1 many times, but since $`\{\mathrm{\Omega }_i\}`$ approaches $`\mathrm{\Omega }`$, it must eventually remain fixed, at which point it must have the same value as the $`k`$-th bit of $`\mathrm{\Omega }`$. Therefore, if the $`k`$-th bit of $`\mathrm{\Omega }`$ is 1, the $`k`$-th bit of $`\{\mathrm{\Omega }_i\}`$ must be 0 for only finitely many values of $`i`$, and so $`P`$ must return 0 for finitely many values of $`N`$ and 1 for infinitely many. On the other hand, if the $`k`$-th bit of $`\mathrm{\Omega }`$ is 0, then the $`k`$-th bit of $`\{\mathrm{\Omega }_i\}`$ must be 1 for only a finite number of values of $`i`$ and $`P`$ must return 1 for finitely many values of $`N`$ and 0 for infinitely many. Either way, as $`N`$ increases, the output of $`P`$ applied to $`k`$ and $`N`$ limits to the $`k`$-th bit of $`\mathrm{\Omega }`$.

It may seem as though this program is computing the bits of $`\mathrm{\Omega }`$ but this is not quite the case. $`P`$ just computes the $`N`$-th ‘guess’ of the $`k`$-th bit. From the infinite sequence of such guesses, the $`k`$-th bit could be determined but $`P`$ does not and cannot put the guesses together like that—it just returns one of them.

Since recursive functions are just a special type of r.e. function, we can apply the dprm Theorem and see that there must be a family of Diophantine equations

(3.1)
$$\chi _1(k,N,x_1,\mathrm{},x_m)=0$$

that has solutions for given values of $`k`$ and $`N`$ if and only if $`P`$ returns 1 when provided with these as input. For a given value of $`k`$, there are solutions for infinitely many values of $`N`$ if and only if the $`k`$-th bit of $`\mathrm{\Omega }`$ is 1.

Thus, by using a more subtle property of the family of Diophantine equations, Chaitin was able to show that algorithmic randomness occurs in number theory: as $`k`$ is varied, there is simply no recursive pattern to whether this family of equations has solutions for finitely or infinitely many values of $`N`$.

By modifying Chaitin’s method slightly, we can find a new way of expressing the bits of $`\mathrm{\Omega }`$ through a family of Diophantine equations . Consider a new program, $`Q`$, that also takes inputs $`k`$ and $`N`$, and begins to compute the sequence $`\{\mathrm{\Omega }_i\}`$. For each value of $`\mathrm{\Omega }_i`$, $`Q`$ checks to see if it is greater than $`\frac{N}{2^k}`$, halting if this is so, and continuing through the sequence otherwise. Since $`\{\mathrm{\Omega }_i\}`$ approaches $`\mathrm{\Omega }`$ from below, we can see that $`\mathrm{\Omega }_i>\frac{N}{2^k}`$ implies that $`\mathrm{\Omega }>\frac{N}{2^k}`$ and conversely, if $`\mathrm{\Omega }>\frac{N}{2^k}`$ there must be some value of $`i`$ such that $`\mathrm{\Omega }_i>\frac{N}{2^k}`$. Therefore, $`Q`$ will halt on $`k`$ and $`N`$ if and only if $`\mathrm{\Omega }>\frac{N}{2^k}`$. Alternatively, we could say that $`Q`$ recursively enumerates the pairs $`k,N`$ such that $`\mathrm{\Omega }>\frac{N}{2^k}`$.

Just as we could determine the $`k`$-th bit of $`\mathrm{\Omega }`$ from the number of values of $`N`$ that make $`P`$ return 1, so we can determine it from the number of values of $`N`$ for which $`Q`$ halts. In what follows, we shall refer to these quantities as as $`p_k`$ and $`q_k`$ respectively.

Unlike $`p_k`$, $`q_k`$ is always finite. Indeed, an upper bound is easily found. Since $`\mathrm{\Omega }<1`$, only values of $`k`$ and $`N`$ such that $`\frac{N}{2^k}<1`$ can possibly be less than $`\mathrm{\Omega }`$ and thus make $`Q`$ halt. Since both $`k`$ and $`N`$ take only values from the positive integers we also know that $`\frac{N}{2^k}>0`$ and thus for a given $`k`$, there are less than $`2^k`$ values of $`N`$ for which $`Q`$ halts and $`q_k\{0,1,\mathrm{},2^k1\}`$.

From the value of $`q_k`$, it is quite easy to derive the first $`k`$ bits of $`\mathrm{\Omega }`$. Firstly, note that $`q_k`$ is equal to the largest value of $`N`$ such that $`\frac{N}{2^k}<\mathrm{\Omega }`$—unless there is no such $`N`$, in which case it equals 0. Either way, its value can be used to provide a very tight bound on the value of $`\mathrm{\Omega }`$: $`\frac{q_k}{2^k}<\mathrm{\Omega }\frac{q_k+1}{2^k}`$. Since $`\mathrm{\Omega }`$ is irrational, we can strengthen this to $`\frac{q_k}{2^k}<\mathrm{\Omega }<\frac{q_k+1}{2^k}`$, which means that the first $`k`$ bits of $`\frac{q_k}{2^k}`$ are exactly the first $`k`$ bits of $`\mathrm{\Omega }`$.

This gives some nice results connecting $`q_k`$ and $`\mathrm{\Omega }`$. The first $`k`$ bits of $`\frac{q_k}{2^k}`$ are just the bits of $`q_k`$ when written with enough leading zeros to make $`k`$ digits in total. Thus $`q_k`$, when written in this manner, provides the first $`k`$ bits of $`\mathrm{\Omega }`$. Additionally, we can see that $`q_k`$ is odd if and only if the $`k`$-th bit of $`\mathrm{\Omega }`$ is 1.

Now that we know the power and flexibility of $`q_k`$, it is a simple matter to follow Chaitin in bringing these results to number theory. The function computed by $`Q`$ is r.e. so, by the dprm Theorem, there must be a family of Diophantine equations

(3.2)
$$\chi _2(k,N,x_1,\mathrm{},x_m)=0$$

that has a solution for specified values of $`k`$ and $`N`$ if and only if $`Q`$ halts when given these values as inputs. Therefore, for a particular value of $`k`$, this equation only has solutions for values of $`N`$ between 0 and $`2^k1`$ with the number of solutions, $`q_k`$, being odd if and only if the $`k`$-th bit of $`\mathrm{\Omega }`$ is 1.

This new family of Diophantine equations improves upon the original one in a couple of ways. Whereas the first method expressed the bits of $`\mathrm{\Omega }`$ in the fluctuations between a finite and infinite amount of values of $`N`$ that give solutions, the second keeps this value finite and bounded, with the bits of $`\mathrm{\Omega }`$ expressed through the more mundane property of parity. It is the fact that this quantity is always finite that leads to many of the new features of this family of Diophantine equations. $`p_k`$ is infinite when the $`k`$-th bit of $`\mathrm{\Omega }`$ is 1 and, since there is only one way in which it can be infinite, it can provide no more than this one bit of information. On the other hand, $`q_k`$ can be odd (or even) in $`2^{k1}`$ ways, which is enough to give $`k1`$ additional bits of information, allowing the first $`k`$ bits of $`\mathrm{\Omega }`$ to be determined.

The fact that $`q_k`$ is always finite also provides a direct reduction of the problem of determining the bits of $`\mathrm{\Omega }`$ to Hilbert’s Tenth Problem. To find the first $`k`$ bits of $`\mathrm{\Omega }`$, one need only determine for how many values of $`N`$ the new family of Diophantine equations has solutions. Since we know that there can be no solutions for values of $`N`$ greater than or equal to $`2^k`$, we could determine the first $`k`$ bits of $`\mathrm{\Omega }`$ from the solutions to $`2^k`$ instances of Hilbert’s Tenth Problem. In fact, we can lower this number by taking advantage of the fact that if there is a solution for a given value of $`N`$ then there are solutions for all lower values. All we need is to find the highest value of $`N`$ for which there is a solution and we can do this with a bisection search, requiring the solution of only $`k`$ instances of Hilbert’s Tenth Problem.777For details see .

Finally, the fact that $`q_k`$ is always finite allows the generalisation of these results from binary to any other base, $`b`$. If we replace all above references to $`2^k`$ with $`b^k`$ we get a new program, $`Q_b`$, with its associated family of Diophantine equations. For this family, the value of $`q_k`$ now gives us the first $`k`$ digits of the base $`b`$ expansion of $`\mathrm{\Omega }`$: it is simply the base $`b`$ representation of $`q_k`$ with enough leading zeroes to give $`k`$ digits. The value of the $`k`$-th digit of $`\mathrm{\Omega }`$ is simply $`q_k`$ mod $`b`$.

Chaitin  did not stop with his Diophantine representation of $`\mathrm{\Omega }`$, but instead moved to exponential Diophantine equations where his result could be presented more clearly. He made this move to take advantage of the theorem that all r.e. sets have *singlefold* exponential Diophantine representations, where a representation is singlefold if each equation in the family has at most one solution.

We can denote the singlefold family of exponential Diophantine equations for the program $`P`$ by

(3.3)
$$\chi _1^e(k,N,x_1,\mathrm{},x_m^{})=0$$

For a given $`k`$, this equation will have exactly one solution for each of infinitely many values of $`N`$ if the $`k`$-th bit of $`\mathrm{\Omega }`$ is 1 and exactly one solution for each of finitely many values of $`N`$ if the $`k`$-th bit of $`\mathrm{\Omega }`$ is 0. We can make use of this to express the bits of $`\mathrm{\Omega }`$ through a more intuitive property.

If we treat $`N`$ in this equation as an unknown instead of a parameter, we get a new (very similar) family of exponential Diophantine equations with only one parameter

(3.4)
$$\chi _1^e(k,x_0,x_1,\mathrm{},x_m^{})=0$$

Since the previous family was singlefold and $`N`$ has become another unknown, there will be exactly one solution to this single parameter family for each value of $`N`$ that gave a solution to the double parameter family. Thus, (3.4) has infinitely many solutions if and only if the $`k`$-th bit of $`\mathrm{\Omega }`$ is 1.

This same approach can be used with our method . There is a two-parameter singlefold family of exponential Diophantine equations for $`Q`$ and this can be converted to a single parameter family of exponential Diophantine equations

(3.5)
$$\chi _2^e(k,x_0,x_1,\mathrm{},x_m^{})=0$$

with between 0 and $`2^k1`$ solutions, the quantity being odd if and only if the $`k`$-th bit of $`\mathrm{\Omega }`$ is 1.

Finally, we have also shown  that both Chaitin’s finitude-based method and our parity-based method can be used to generate polynomials for $`\mathrm{\Omega }`$. For a given family of Diophantine equations with two parameters,

(3.6)
$$D(k,N,x_1,\mathrm{},x_m)=0$$

we can construct a polynomial, $`W`$, where

(3.7)
$$W(k,x_0,x_1,\mathrm{},x_m)x_0\left(1(D(k,x_0,x_1,\mathrm{},x_m))^2\right).$$

Note that the parameter, $`N`$, is again treated as an unknown and thus denoted $`x_0`$.

If we restrict the values of the variables to positive integers then, for a given $`k`$, this polynomial takes on exactly the set of all values of $`N`$ for which (3.6) has solutions. We can thus use this method on $`\chi _1=0`$ and $`\chi _2=0`$, generating polynomials that express $`p_k`$ and $`q_k`$ in the number of distinct positive integer values they take on for different values of $`k`$. We therefore have a polynomial whose number of distinct positive integer values fluctuates from odd to even and back in an algorithmically random manner as a parameter $`k`$ is increased.

## Acknowledgements

We would like to thank Gregory Chaitin for helpful discussion and Cristian Calude for kindly offering to present this paper for us at DMTCS’03. TDK wishes to acknowledge the continuing support of Peter Hannaford.


text:
# Scaling Law for the Distribution of Fluctuations of Share Volume

## 1 Intoroduction

There is a long tradition of studies on power-law scaling in financial markets. Scaling of market prices was first reported by Mandelbrot in 1963 in his seminal work on cotton prices . He presented empirical evidence that the cumulative distribution of daily fluctuations in cotton prices followed a power-law asymptotic behavior, characterized by an exponent $`\alpha 1.7`$; he also proposed stable L$`\stackrel{´}{e}`$vy distributions  as candidates for the probability density function of price changes in financial assets. Fama  analyzed quantitatively the daily data on the share prices of 30 companies over a 5-year period, and showed a power-law distribution of the price changes. Recent studies have confirmed the presence of scaled distributions of price changes in various financial markets: exchange rate markets , German stocks , the stock price index , ,and individual stocks . These results suggest the possibility of the existence of universal properties of price changes.

Although trading volume is, like price, an important quantity that characterizes the activities of financial markets, only a few attempts have so far been made to understand the statistical properties of trading volumes111Gopikrishnan, et.al.  are among those who first studied the statistical properties of share volume. They showed that the distribution of share volume follows a power law with an exponent close to $`1.7`$.. In this paper we study the statistical properties of fluctuations of share volume. By analyzing a database of the daily transactions for all securities listed on the Tokyo Stock Exchange, we found that the cumulative distributions of daily fluctuations of share volumes can be well described by a power-law decay and characterized by an exponent within the stable L$`\stackrel{´}{e}`$vy domain $`0<\alpha <2`$. Furthermore, we found that the distributions for a number of companies follow Zipf fs law, that is, the cumulative distributions of fluctuations of share volume have exponents very close to unity 222Zipf fs law is found in various fields, e.g. ,     ..

The paper is organized as follows. The following section presents the details of our empirical results on the distribution of fluctuations in share volume, and the final section presents some concluding remarks.

## 2 The Distribution of Fluctuations of Share Volume

We investigated the statistical properties of the daily data of share volume, using a database of the daily volumes for all securities listed on the Tokyo Stock Exchange. We selected 1050 companies that each had an unbroken series of daily trading activity for the entire 27-year period from January 1975 to January 2002. Each time series of daily trading volume had approximately 7000 data points, corresponding to the number of trading days in the 27-year period.

The basic quantity studied for individual companies was the volume $`Q(t)`$, defined as the number of shares traded in a trading day. For each company, we analyzed increments of volume, $`V=Q(t)Q(t1)`$. Nissin Steel, a large Japanese steelmaker, serves as a typical example of the 1050 companies we analyzed. Fig. 1 shows the probability distribution of the increments in share volume $`V`$ for Nissin Steel. At a glance, the contrast with the Gaussian distribution curve is striking. The mean value is close to zero. Skewness is positive, and kurtosis is significantly positive. The important result is that the probability distribution of share volume exhibits apparently greater probability mass in the tails and in the center than dose the standard Normal, which is characterized by the empirical distribution that has been repeatedly observed in various market data.

Fig. 2 (a) and (b) show the cumulative probability distributions of increments of share volume $`V`$, i.e. the probability for an increment larger or equal to a threshold $`x`$, for Nissin Steel in log-log scale. The cumulative distribution $`P(Vx)`$ shows a power-law decay:

$$P(Vx)\frac{1}{x^\alpha }.$$
(1)
Regression fits in the region $`4\times 10^2x2\times 10^4`$ yield

$$\alpha =\{\begin{array}{cc}\hfill 0.97\pm 0.007,R^2=0.996& \text{(positive tail)}\hfill \\ \hfill 1.02\pm 0.003,R^2=0.999& \text{(negative tail)}\hfill \end{array}$$
(2)
For both positive and negative tails, a power-law asymptotic behavior well fits the data over the range $`4\times 10^2x2\times 10^4`$. These estimates of the exponent $`\alpha `$ are well inside the stable L$`\stackrel{´}{e}`$vy range, which requires $`0<\alpha <2`$. More noteworthy is that these are very close to unity, i.e., the distribution of share volume for Nissin Steel follows Zipf fs law.

To confirm the robustness of the above analysis, we repeated this analysis for each time series of share volume for each of the 1050 companies. For all of the companies, the asymptotic behavior of the functional form of the cumulative distributions was consistent with a power law (1). The estimates of $`\alpha `$ were sensitive to the bounds of the regression used for fitting. Thus we used the coefficient of determination $`R^2`$ of the linear regression line against a standard and determined the appropriate values of the exponent $`\alpha `$. We chose the results of the regression fit as an appropriate value of the exponents where the coefficient of determination was greater than $`0.995`$. The estimates of the exponent $`\alpha `$ were, for all but $`6`$ of the $`1050`$ companies, within the stable L$`\stackrel{´}{e}`$vy domain, $`0\alpha <2`$. In Figure 3 we show the histogram for $`\alpha `$, obtained from power-law regression fits to the positive tails (a) and negative tails (b) of the individual cumulative distributions of all $`1050`$ companies. Table 1 lists the proportion of observations $`\alpha `$ in each class. As shown in Table 1, the most probable value is the class $`0.91.1`$ for the estimates $`\alpha `$ of the positive tails and the class $`1.11.3`$ for the estimates $`\alpha `$ of the negative tails. Thirty-eight percent of all the estimates $`\alpha `$ for the positive tails were between $`0.9`$ and $`1.1`$, while $`31`$ percent of $`\alpha `$ for the negative tails were between $`0.9`$ and $`1.1`$. More than $`30`$ percent of the companies had a distribution of fluctuations of share volume characterized by an exponent close to unity, ($`0.9<\alpha <1.1`$).

## 3 Concluding Remarks

We have presented the statistical properties of the distribution of fluctuations in share volume for individual companies. We found (1) that the distribution of fluctuations in share volume for almost all of the companies was consistent with a power-law behavior, characterized by an exponent inside the L$`\stackrel{´}{e}`$vy stable regime, $`0<\alpha <2`$. We also found that more than 35 percent of the distributions of fluctuations of share volume followed the Zipf fs law approximately.

Share prices are largely determined largely by the relation between supply and demand. Therefore, in order to get a comprehensive understanding of the market dynamics, we need to investigate the relationship between price and volume. Our primitive work  on the share returns for the same companies showed that the return distributions were characterized by the double exponential distribution. This result was inconsistent with previous works on stock returns , . Two questions must be considered next. The first is, What mechanism generates Zipf fs law of volume? The second question is whether or not there exist statistical laws that present the relationship between volume and price.

Although so far no theory has been established to describe the relationship between price and volume in financial markets333Our recent work  presented a stockastic model of financial markets, which related share returns to trading-volume changes., our empirical results give conditions that any empirically accurate theories of stock markets have to satisfy. These empirical and theoretical studies remain for future work.

## 4 Acknowledgment

We thank Michiyo Kaizoji for classifying data and encouragements.


text:
# Weakly Non-ideal Bose Gas: Comments on Critical Temperature Calculations

## I Introduction

During the recent upsurge of activity in the field of Bose-Einstein condensation (BEC), some older questions concerning the weakly interacting Bose gas have been revived. One such un-resolved question is the effect of a weak repulsion on the critical temperature $`T_\mathrm{c}`$. Calculations on this question has produced widely dissimilar results: increases and decreases of $`T_\mathrm{c}`$ proportional to $`a`$, $`\sqrt{a}`$, $`a\mathrm{ln}a`$, etc. have been reported by various authors. (Here $`a`$ is the scattering length.) It is only during the past 2-3 years that the community has been converging toward a consensus.

Considering the diverse approaches to the $`T_\mathrm{c}`$ problem, and the variety of issues that have arisen in this connection, it seems appropriate at this point to attempt to draw some connections between the different approaches, and point out some misconceptions and fine points. The purpose in this article is to review and comment on existing work, rather than to present new calculations.

The model is described by the Hamiltonian

$$\widehat{H}=\underset{𝐤}{}ϵ_𝐤\widehat{b}_𝐤^{}\widehat{b}_𝐤+\frac{1}{2V}\underset{𝐩,𝐪,𝐤}{}U(𝐤)\widehat{b}_{𝐩+𝐤}^{}\widehat{b}_{𝐪𝐤}^{}\widehat{b}_𝐩\widehat{b}_𝐪,$$

where $`ϵ_𝐤=k^2/2m`$ is the free-gas spectrum, and $`\widehat{b}`$, $`\widehat{b}^{}`$ are bosonic operators. In principle, one is interested in any possible interaction function $`U(𝐤)`$, but the simplest possible form is a delta function in real space, so that $`U(𝐤)=U`$ is momentum-independent. To first order $`U`$ is related to the $`s`$-wave scattering length $`a`$ by $`U=4\pi \mathrm{}^2a/m`$. A dimensionless measure of the interaction is the quantity $`an^{1/3}`$, where $`n=N/V`$ is the density. The program is to examine the shift in the critical temperature $`T_\mathrm{c}`$, as compared to a noninteracting gas of the same density, as a function of $`an^{1/3}`$, for small $`an^{1/3}`$.

In section II, the background is set by summarizing different calculations and basic results. Section III contains the comments on various approaches to the $`T_\mathrm{c}`$ problem. In the concluding section IV, I suggest some directions for additional calculations.

## II Background and Perspective

Early predictions for the shift in transition temperature appear in the work of Lee, Yang and Huang, around forty years ago, as summarized in baym\_bigpaper . Recently, Huang has published huang\_prl  a prediction for an increase in $`T_\mathrm{c}`$ due to the inter-particle repulsion, proportional to the square root of the interaction strength, $`T_\mathrm{c}T_\mathrm{c}^{(0)}\sqrt{an^{1/3}}`$. Schakel in a short note schakel\_comment , and Baym et al in baym\_bigpaper , have both discussed the problem with this work; in particular baym\_bigpaper  contains an analysis showing how Huang’s result is generated spuriously from a truncation of the virial expansion used at mean field level.

A similar spurious result appearing in the literature is that of Toyoda toyoda , who in 1982 performed mean-field-like calculations deriving a Landau-Ginzburg type thermodynamic function. His prediction was a *decrease* of the transition temperature, $`\mathrm{\Delta }T_\mathrm{c}\sqrt{an^{1/3}}`$. Again, this is believed to be an artifact of mean field approximations. The present understanding is that there is no mean field effect on the transition temperature, and that higher-order correlations are required to find the leading shift in $`T_\mathrm{c}`$.

A more recent work that also predicts a negative shift in $`T_\mathrm{c}`$ is the canonical-ensemble work of Wilkens *et al* wilkens . This work has been critiqued in mueller\_baym\_finite . I discuss some related issues in section III.5.

In addition, in view of the continuing appearance of negative-shift predictions, I give in section III.2 a simple argument for the *increase* of the transition temperature due to the addition of a repulsive interaction. This argument makes use of eq (2), which is derived independent of any particular method of treating interaction effects.

During the past several years, a large body of work has been published by Baym, Laloë and collaborators baym\_prl ; baym\_N ; baym\_a2log ; HGL\_ursell ; mueller\_baym\_finite ; baym\_bigpaper ; baymfriends\_numerics , and a number of authors have done followup work arnold\_N ; arnoldmoore ; svistunov . These workers predict a positive shift in $`T_\mathrm{c}`$ that to leading order is *linear* in the scattering length; $`\mathrm{\Delta }T_\mathrm{c}=c_1an^{1/3}+𝒪(a^2n^{2/3})`$. Attempts to calculate the coefficient $`c_1`$, however, have continued to give fluctuating results, and this has prompted the Baym group to predict baym\_a2log  a logarithmic contribution at the *next* order, i.e., a $`a^2\mathrm{ln}a`$ contribution. The presence of such a non-analyticity might explain the difficulty in any numerical estimate of $`c_1`$.

Stoof’s renormalization-group analysis stoof1 ; stoof2  of the three-dimensional bose gas has led to a prediction stoof2 ; stoof\_private  of an increase of the transition temperature proportional to $`a|\mathrm{ln}a|`$, i.e., a non-analyticity at *leading* order. A leading non-analytic behavior would actually explain even better the uncertain results obtained in attempts to calculate the coefficient $`c_1`$, based on the assumption of linear shift $`\mathrm{\Delta }T_\mathrm{c}c_1an^{1/3}`$.

An older attack on the $`T_\mathrm{c}`$ problem is that of Kanno kanno12 ; kanno3 . In a series of papers, Kanno introduced a “quasilinear” canonical transformation approach to the many-body problem, which he used to calculate (among other things) the free energy and hence transition temperature of the bose gas. Since Kanno’s work has mostly been ignored by the current $`T_\mathrm{c}`$ community, I give a description of his work in section III.6. I also show how some of his derivations can be performed in a more transparent way.

Schakel’s work on developing an effective theory schakel\_ijmp  of the weakly non-ideal Bose gas has led to the prediction schakel\_boulevard  of a shift in the transition temperature that is strictly linear in the interaction. The connection of this calculation with other approaches has not yet been explored in detail. Some issues are pointed out in section III.4.

Our recent report our\_f02  addresses the question of $`T_\mathrm{c}`$ within a quasiparticle framework. Other than working out details of a quasiparticle description (ultraviolet divergence, effect of quasiparticle lifetime), we also provide a non-selfconsistent method of regularizing the infrared divergences appearing in the description of the transition point. This results in a shift of $`T_\mathrm{c}`$ approximately $``$ $`a\sqrt{|\mathrm{log}a|}`$. The procedure is uncontrolled; however, this result points strongly to a non-analytic contribution in the dependence of $`T_\mathrm{c}`$ on the interaction.

Other relevant calcultaions include Kleinert’s “five-loop” determination kleinert  of the coefficient $`c_1`$ mentioned above, assuming a puely linear leading shift, and the “linear $`\delta `$-expansion” calculations delta  of de Souza Cruz *et al*. Kleinert has questioned kleinert  the $`\delta `$-expansion method, but a detailed analysis of the relationship of this approach to others seems to be lacking in the literature.

### II.1 Mean Field

For a momentum-independent potential, the two mean field contributions to the self-energy (Hartree and Fock) are equal, $`\mathrm{\Sigma }^{(1)}(𝐤,\omega )=Un+Un=2Un`$. Since the self-energy is momentum-independent, it can simply be absorbed into a redefinition of the chemical potential.

$$E_𝐤^{(1)}=\frac{k^2}{2m}\mu +\mathrm{\Sigma }^{(1)}(𝐤,E_𝐤^{(1)})=\frac{k^2}{2m}\mu ^{}$$
(1)
$$N=N_0+\frac{V}{\lambda ^3}g_{3/2}(e^{\beta \mu ^{}})$$

Here $`\mu ^{}=\mu \mathrm{\Sigma }^{(1)}=\mu 2Un`$. The transition now simply corresponds to $`\mu ^{}=0`$ instead of $`\mu =0`$, and $`T_\mathrm{c}`$ remains unchanged. The lesson here is that a momentum-dependent self-energy is required for a shift in $`T_\mathrm{c}`$.

In second-order calculations, (e.g. in our\_f02 ,) it is usual to have the 1st-order self-energy already absorbed in the chemical potential, i.e., to use $`\mu ^{}`$ as the chemical potential. This allows one to exclude those second-order diagrams which can be obtained by 1st-order self-energy insertions into 1st-order diagrams.

*Spatially extended interaction*. In a real system, of course, the description of the interaction in terms of a single parameter ($`U_0`$ or $`a`$) is valid only for very weak interactions or at very low temperatures/densities; if these conditions are not satisfied one has to take into account the detailed shape of the potential $`V(x)=V_0\sigma (x)`$ which is not a delta function in real space (not a constant in momentum space) any more. In such a case we get a momentum-dependent Fock term, and hence a contribution to $`T_\mathrm{c}`$ from the first-order self-energy. Assuming a Gaussian $`V(x)`$ of width $`a`$ and height $`V_0`$, one gets a correction to the spectrum that can be incorporated as an effective mass $`m^{}=m[1V_0a^2nm/3]^1`$, resulting in a negative shift fetter :

$$\frac{\mathrm{\Delta }T_\mathrm{c}}{T_\mathrm{c}^{(0)}}=\frac{1}{3}mnV_0a^2.$$

The negative shift at higher densities is present in liquid helium and presumably in other Bose fluids. However, for weak interactions, this shift is less important than the leading $`𝒪(a)`$ effect, which is the current theoretical challenge.

### II.2 $`T_\mathrm{c}`$ in Liquid Helium

As a highly correlated liquid, the properties of liquid helium are not expected to be described well by the dilute weakly-interacting Bose gas model. However, the deviation of the lambda transition temperature (2.2K), from the BEC temperature of an ideal gas with the same density (3.1K), can be mostly explained using a simple effective mass description. In other words, the experimental $`T_\mathrm{c}`$ is reproduced quite closely by replacing the helium atomic mass $`m`$ by a quasiparticle effective mass $`m`$* in the ideal gas expression $`T_\mathrm{c}^{(0)}=\left(2\pi /m[\zeta (\frac{3}{2})]^{2/3}\right)n^{2/3}`$.

The effective mass in the liquid exceeds the bare mass, because of the inertia of the medium which has to make way for any atom to move. This results in a lowering of $`T_\mathrm{c}`$. The effect is essentially equivalent to that described in the last section for a momentum-dependent interaction.

Feynman feynman53  calculated a partition function for liquid helium that takes into account the effective mass enhancement described above. Using this partition function, Chapline calculated chapline  a $`T_\mathrm{c}`$ quite close to the actual $`\lambda `$-point of helium. This is equivalent to using the effective mass in the ideal gas expression for $`T_\mathrm{c}^{(0)}`$.

According to current understanding, therefore, one expects the following variation of $`T_\mathrm{c}`$ with $`n^{1/3}a`$: an initial linear-like increase, followed by a maximum and then a decrease due to the effective-mass mechanism described above, until eventually $`\mathrm{\Delta }T_\mathrm{c}`$ becomes negative. This kind of dependence has been seen in the Monte-Carlo calculation of Grüter, Ceperley, Laloë baymfriends\_numerics , and in Reppy’s experiments reppy  described below.

### II.3 Relevant Experiments

One could think of studying the $`T_\mathrm{c}`$ problem experimentally, e.g., in the context of trapped atomic Bose gases. Unfortunately, trap effects dominate in the dependence of $`T_\mathrm{c}`$ on the interaction for these systems. Turning on a repulsion for a trapped Bose gas causes an expansion of the gas, reducing the density and hence decreasing the transition temperature $`T_\mathrm{c}n^{2/3}`$ trap\_tc . In contrast to the intrinsic effect under review here, trapping effects on $`T_\mathrm{c}`$ are found at mean-field level already.

A different experimental realization of the weakly non-ideal bose gas is the 4He-Vycor system explored by Reppy *et al* reppy . Here the small-$`an^{1/3}`$ regime is reached by enormously decreasing the 4He density. Reppy’s investigation of the transition temperature shows a linear-like increase of $`T_\mathrm{c}`$ at low $`an^{1/3}`$, followed by a maximum and a decrease to values lower than $`T_\mathrm{c}^{(0)}`$. The general behavior is satisfying, but for the initial increase, the data is not yet good enough to choose between $`a`$, $`a\mathrm{ln}a`$, or $`a\sqrt{\mathrm{ln}a}`$-like behavior.

## III Comments on Various Approaches

### III.1 Baym and collaborators

The work of Baym, Laloë and collaborators baym\_prl ; baym\_N ; baym\_a2log ; mueller\_baym\_finite ; baym\_bigpaper ; HGL\_ursell ; baymfriends\_numerics  establishes definitely that a repulsive interaction *increases* $`T_\mathrm{c}`$, and that the dependence of $`\mathrm{\Delta }T_\mathrm{c}`$ on the interaction is of linear order, as opposed to quadratic or $`\sqrt{a}`$.

One should note, however, that the Baym group’s work relies heavily on power-counting and dimensional arguments. This kind of argument typically does not catch logarithmic corrections. The possibility of leading non-analytic dependence of $`\mathrm{\Delta }T_\mathrm{c}`$ on the interaction, therefore, needs to be further examined.

It has been suggested schakel\_private  that the large-$`N`$ calculation baym\_N  is a reliable indicator that the leading shift is purely linear, because of intuition gained from other systems where large-$`N`$ calculations capture all logarithms. I would like to point out, however, that the $`N2`$ models belong to a different universality class, and have different excitation spectra as compared to the usual ($`N=2`$) Bose gas. More care is therefore required in interpreting results from the large-$`N`$ approach.

The Baym-Laloë collaboration seems to have been the first to relate the single-particle spectrum at $`T_\mathrm{c}`$ to the shift $`\mathrm{\Delta }T_\mathrm{c}`$. Another noteworthy aspect is the use of Ursell operators HGL\_ursell ; baym\_bigpaper , and a corresponding diagram technique, to describe particle correlations. (As far as I can tell, for this problem, the Ursell operator technique does not give any extra information compared to the usual Green function techniques.)

### III.2 Direction of Shift, Infrared Spectrum

Once the shift in $`T_\mathrm{c}`$ is related to the modification of the single-particle spectrum, the direction of the shift can be argued using well-known facts about the spectrum of the Bose gas.

We will use the following expression for the shift, as appearing in our\_f02 :

$$\frac{\mathrm{\Delta }T_\mathrm{c}}{T_\mathrm{c}^{(0)}}\frac{2}{3n}\frac{d^3k}{(2\pi )^3}\left[f_\mathrm{c}(\xi _𝐤)f_\mathrm{c}(ϵ_𝐤)\right].$$
(2)
Here $`f_\mathrm{c}`$ is the Bose function at the critical point; $`f_\mathrm{c}(x)=(e^{x/kT_\mathrm{c}}1)^1`$; and $`\xi _𝐤`$ is the quasiparticle spectrum at $`T_\mathrm{c}`$, i.e., $`\xi _𝐤=ϵ_𝐤+\mathrm{\Sigma }(𝐤,\xi _𝐤;T=T_\mathrm{c})`$.

The expression used by Baym *et al* baym\_prl ; baym\_bigpaper , $`\mathrm{\Delta }T_\mathrm{c}𝑑kU(k)/[k^2+U(k)]`$, is the high-temperature (zero Matsubara frequency) version of (2). Noting that the major contribution to (2) comes from the infrared, one can write approximately

$$\mathrm{\Delta }T_\mathrm{c}_0^{\mathrm{cutoff}}𝑑kk^2[ϵ_𝐤^1\xi _𝐤^1].$$
(3)
This equation might introduce or leave out logarithmic corrections; however the following argument is unchanged for (2).

If the quasiparticle spectrum $`\xi _𝐤`$ is “harder” than the bare spectrum $`ϵ_𝐤=k^2/2m`$, i.e., if $`\xi _𝐤`$ is sub-quadratic in $`k`$ or quadratic with an effective mass $`m^{}<m`$, then the $`ϵ_𝐤`$ integral dominates and the shift is positive. On the other hand if the spectrum at the transition point is “softened” by the interaction, then eq (3) \[ or (2)\] leads to a decrease in the transition temperature.

Since the weakly interacting bose gas is known to have a linear (Bogoliubov) infrared spectrum below the critical temperature, and a quadratic spectrum above $`T_\mathrm{c}`$, it is extremely reasonable to assume, even without any RG or perturbation-theoretic calculations, that the spectrum at the critical point should be something intermediate between linear and quadratic. As pointed out above, this means a *positive* shift in the transition temperature.

The long-wavelength spectrum at the transition point should be $`\xi _𝐤k^{2\eta }`$, where $`\eta `$ is the anomalous dimension. The present system falls in the same universality class as the $`XY`$ model, or the $`N=2`$ quantum rotor model. For this universality class, large-$`N`$ and $`ϵ`$-expansion calculations give respectively $`\eta 0.14`$ and $`\eta 0.02`$, and a recent Monte Carlo calculation eta\_montecarlo  has produced the value $`\eta 0.038`$.

Second-order treatments of the self-energy (e.g., the Baym group’s baym\_prl ; baym\_bigpaper  and our our\_f02  work) over-estimates the modification of the spectrum, giving a $`k^{3/2}`$ infrared behavior ($`\eta =0.5`$). On the other hand, some other treatments (such as Stoof’s stoof2  and Schakel’s schakel\_ijmp ; schakel\_boulevard ) neglect the modification altogether ($`\eta =0`$).

### III.3 Renormalization group treatment (Stoof)

The RG flow equations derived by Bijlsma & Stoof stoof2  allow one to numerically locate the transition point. Under RG flow, the physical effective chemical potential flows to positive and negative values, if one starts out respectively in the Bose-condensed and the symmetry-unbroken phase. The critical point is therefore identified by finding parameters for which the chemical potential flows to zero. Stoof finds stoof2 ; stoof\_private  that the dependence of $`\mathrm{\Delta }T_\mathrm{c}`$ on $`a/\lambda `$ has more structure than purely linear; $`\mathrm{\Delta }T_\mathrm{c}(a/\lambda )\mathrm{ln}(a/\lambda )`$.

The differences between Stoof’s and the Baym group’s results merit further study. Stoof has suggested stoof\_private  that the source of the difference is not the infrared but the *ultraviolet*. The correction to the spectrum, $`U(k)\mathrm{\Sigma }(𝐤,E_𝐤)\mathrm{\Sigma }(0,0)`$, has an ultraviolet $`\mathrm{ln}k`$ behavior in the Baym group’s treatment baym\_bigpaper ; baym\_prl . In the RG framework, the corresponding quantity is the running coupling constant $`\mu (\mathrm{\Lambda })`$, which shows no ultraviolet logarithms in the calcultation of ref stoof2 . This leads to a logarithmic difference in the two predictions for $`T_\mathrm{c}`$. Stoof also suggests stoof\_private  that an improved RG calculation, including the momentum dependence of the effective two-body interaction at the transition point, would shed light on the difference between the two approaches.

### III.4 Schakel: IR regularization

In treating his effective-action theory near the transition point, Schakel schakel\_ijmp ; schakel\_boulevard  is also faced with infrared (IR) divergences. One of the remarkable features of this work is an innovative analytic continuation procedure used to regulate IR problems.

Baym et al have remarked baym\_bigpaper  that Schakel’s prescription ($`\zeta `$-function regularization) must be incorrect because it gives an incorrect result for the compressibility of the ideal gas. This argument, unfortunately, is derived from an equation that Schakel derived explicitly for the nontrivial broken-symmetry ground state schakel\_boulevard ; schakel\_private  of an interacting gas ($`\mu >0`$). The same equations are not applicable to the noninteracting case ($`\mu 0`$) for which the structure of the vacuum is different (more trivial). The situation is therefore more complicated than indicated in Baym et al’s note. In response, Schakel has provided a more detailed explanation and justification schakel\_rejoinder  for the $`\zeta `$-function continuation procedure.

It seems to me that the regularization procedure itself is quite possibly sound. My concern with Schakel’s derivation is different — just below the transition point (broken-symmetry phase), Schakel uses a quadratic spectrum. A more accurate treatment would use the phonon-like spectrum $`E_𝐤\sqrt{ϵ_𝐤^2+2Un_0(T)ϵ_𝐤}`$. The effect such a spectrum would have on the thermodynamic potential or on $`T_\mathrm{c}`$ is not clear. (Redoing Schakel’s $`T_\mathrm{c}`$ calculation with a phonon-like spectrum seems quite difficult.)

### III.5 Canonical Ensemble: Role of $`N_0(T_\mathrm{c})`$.

In Wilkens et al’s paper wilkens , calculations have been performed using a canonical ensemble of $`N`$ particles, and the shift $`\mathrm{\Delta }T_\mathrm{c}`$ obtained is linear and *negative* in the scattering length. It has been pointed out mueller\_baym\_finite  that the error lies in evaluating quantities (derivatives) at $`N_0=0`$, instead of evaluating them at the expectation value of $`N_0(T_\mathrm{c})`$. Mueller et al mueller\_baym\_finite  have also used finite-size scaling ideas in combination with a first-order calculation in the canonical ensemble, and obtain a positive linear shift.

The idea that the zero-momentum state occupancy $`N_0`$ becomes “microscopic” at and above the transition needs to be clarified. One finds $`lim_{N0}(N_0/N)=0`$ for $`TT_\mathrm{c}`$, where $`N`$ is the total number of particles, but not necessarily $`N_0𝒪(1)`$. In fact, for both the ideal and the interacting Bose gas, $`N_0(T_\mathrm{c})N^{2/3}`$, which is large compared to $`1`$ though small compared to $`N`$.

In view of the continuing appearance of the $`N_0(T_\mathrm{c})𝒪(1)`$ misconception, I give here a simple derivation of $`N_0(T_\mathrm{c})=\gamma N^{2/3}`$, together with an expression for the coefficient $`\gamma `$, for the ideal-gas case.

Let us consider the ideal-gas relationship $`N=N_0+\frac{V}{\lambda ^3}g_{3/2}(e^\alpha )`$, where $`\alpha =\beta \mu `$, and $`g_\nu `$ refers to a Bose-Einstein integral function BE\_functions . Since $`\alpha `$ is small near the transition, one can use the Robinson expansion BE\_functions  for $`g_{3/2}`$. Also, inverting $`N_0=1/(e^{\beta ϵ_0+\alpha }1)`$ one gets $`\alpha 1/N_0`$. Combining these yields, near the transition,

$$N_0^{3/2}+\left[\frac{V}{\lambda ^3}\zeta \left(\frac{3}{2}\right)N\right]\sqrt{N_0}\frac{V}{\lambda ^3}(2\sqrt{\pi })0.$$

Since $`n\lambda ^3=\zeta (\frac{3}{2})`$ at the transition, the coefficient of $`\sqrt{N_0}`$ needs to vanish. Therefore

$$N_0(T_\mathrm{c})\left[\frac{2\sqrt{\pi }}{\zeta (3/2)}\right]^{2/3}N^{2/3}.$$
(4)
The $`𝐤=0`$ occupancy $`N_0`$ can actually be calculated for all $`T`$, using an extra Lagrange multiplier $`\mu _0`$ conjugate to $`N_0`$, in addition to usual chemical potential $`\mu `$ coupling to $`N`$. Using this “two chemical potentials” formalism our\_unpubl , one can show that for the ideal gas $`N_0=1`$ occurs at a temperature $`T2.6T_\mathrm{c}^{(0)}`$, well above the transition point.

For the interacting-gas case, an equation corresponding to (4) is difficult to derive, but finite-size scaling arguments cardy\_book ; mueller\_baym\_finite  show that $`N_0N^{2/3}`$ continues to hold, even though the nonideal Bose gas belongs to a different universality class.

The definition of the transition point used in ref wilkens  is in terms of the probability

$$P_n(T,N)\mathrm{tr}\left[\delta _{\widehat{N}_0,n}e^{\beta \widehat{H}}\delta _{\widehat{N},N}\right]$$
(5)
for having a condensate occupancy $`N_0=n`$. The claim is that $`P_n(T,N)`$ decreases monotonically with $`n`$ for all temperatures above the transition, and has a peak at nonzero $`n`$ for in the condensed phase, and therefore the transition should be defined as the thermodynamic limit of the zero of $`D(T,N)=P_0(T,N)P_1(T,N)`$.

The correspondence between the canonical particle-counting statistics, and the grand-canonical ideas of average occupancies, is not trivial. One important point to note is that a monotonically decreasing $`P_n`$ does not imply a zero expectation value for $`N_0`$, even though $`P_n`$ has a maximum at $`n=0`$. (For example, using $`P_ne^{\alpha n}`$ gives $`N_0=_nnP_n=1/[e^\alpha 1]`$.) The stipulation of a monotonically decreasing $`P_n`$ therefore does not imply a zero or $`𝒪(1)`$ value of $`N_0`$ at the transition.

### III.6 Kanno’s Approach: Using Modified Thermodynamic Potential

Kanno developed an unusual approach kanno12  to study many-body systems, and also used this approach to calculate Bose gas properties kanno3 .

*The formalism*.
Kanno’s calculation of thermodynamic quantities involves a “quasilinear” canonical transformation kanno12  of the free-particle operators $`\widehat{b}`$ ($`\widehat{b}^{}`$):

$$\begin{array}{c}\widehat{b}_i=\widehat{\gamma }_i+\underset{j;k,l}{}A_i(j;k,l)\widehat{\gamma }_j^{}\widehat{\gamma }_k\widehat{\gamma }_l\hfill \\ \hfill +\underset{j,k;l,m,n}{}A_i(j,k;l,m,n)\widehat{\gamma }_j^{}\widehat{\gamma }_k^{}\widehat{\gamma }_l\widehat{\gamma }_m\widehat{\gamma }_n+\mathrm{},\end{array}$$
(6)
resulting in a Hamiltonian having a particular (diagonal-like) form. The formalism involves examining the properties of the coefficients $`A`$, and the coefficients appearing in the quasi-diagonal hamiltonian. Applying this formalism to the weakly repulsive bose gas kanno3 , one can derive the Bogoliubov spectrum, entropy and pressure above and below the transition, and many other quantities.

Kanno does not provide direct connections to the usual language of perturbative many-body theory. I would like to point out that many of his physical results, if not all, are equivalent to results obtained in more mainstream approaches. For example, the quantities $`A(𝐩)`$ and $`B(𝐩)`$, appearing in the second paper of the series kanno12 , are equivalent to the normal and anomalous self-energies that appear in modern descriptions.

*The shift in $`T_\mathrm{c}`$*.
Kanno’s defines kanno3  a modified thermodynamic potential

$$\mathrm{\Omega }^{}(N_0,\mu ^{},T)=\frac{1}{\beta }\mathrm{ln}\mathrm{tr}^{}e^{\beta (\widehat{H}\mu ^{}\widehat{N^{}})},$$

with an unusual Lagrange multiplier $`\mu ^{}`$ conjugate to the number of non-condensate particles, $`N^{}=NN_0`$, as opposed to the usual chemical potential which is conjugate to the total particle number $`N`$. It follows that

$$NN_0=N^{}=\frac{\mathrm{\Omega }^{}(N_0,\mu ^{},T)}{\mu ^{}}.$$
(7)
Once $`\mathrm{\Omega }^{}`$ is calculated in some approximation, comparing this relation to $`NN_0=Vg_{3/2}(e^{\beta \mu })/\lambda ^3`$ for the free bose gas at the transition point yields the shift $`\mathrm{\Delta }T_\mathrm{c}`$, in terms of the quantity $`\mu ^{}(T_\mathrm{c})`$.

The quantity $`\mu ^{}`$ itself is obtained from a consideration of the free energy $`F(N_0,N,T)=\mathrm{\Omega }^{}(N_0,\mu ^{},T)+\mu ^{}N^{}`$. The actual value of $`N_0`$ chosen by the system is the one that minimizes $`F(N_0,N,T)`$, i.e.,

$$0=\frac{F}{N_0}|_{T=T_\mathrm{c}}=\frac{\mathrm{\Omega }^{}}{N_0}|_{T=T_\mathrm{c}}\mu ^{}(T_\mathrm{c}).$$
(8)
This equation allows one to express $`\mu ^{}(T_\mathrm{c})`$ in terms of the interaction strength.

*Computing $`\mathrm{\Omega }^{}`$ in more usual language*.
As advertised in section II, I show here how the modified thermodynamic potential $`\mathrm{\Omega }^{}`$ (which is the quantity relevant to the $`T_\mathrm{c}`$ problem) can be obtained from a diagrammatic second-order perturbative calculation.

The Feynman rules for a diagrammatic perturbative calculation of the modified therodynamic potential $`\mathrm{\Omega }^{}(\mu ^{},N_0,T)`$ are very similar to those for the usual $`\mathrm{\Omega }(\mu ,T)`$. One has to use a momentum-dependent chemical potential $`\mu _𝐤=\mu ^{}(1\delta _{𝐤,0})`$, and separate out $`𝐤=0`$ from the momentum-sums at each stage. Two first order and two second-order diagrams need to be calculated:

$$\mathrm{\Omega }^{}\mathrm{\Omega }^{(0)}+\text{}+\text{}+\text{}+\text{}.$$

Other second-order diagrams are excluded by considering the chemical potential to be renormalized upto first order, as discussed in section II.1.

The zeroth and first order terms can be evaluated exactly. UV divergences appearing at second order can be removed by eliminating the bare potential $`U`$ in favor of the 2-body $`t`$-marix, as done for the self-energy in our\_f02 . The result is

$$\begin{array}{c}\mathrm{\Omega }^{}=\frac{V}{\beta \lambda ^3}g_{5/2}(e^{\beta \mu ^{}})+\frac{t}{V}N^2+4t^2N_0_𝐩_𝐪\frac{f(ϵ_𝐩\mu ^{})f(ϵ_𝐪\mu ^{})}{ϵ_𝐩ϵ_𝐪ϵ_{𝐩𝐪}}\hfill \\ \hfill +2t^2V_𝐩_𝐪_𝐤\frac{f(ϵ_𝐩\mu ^{})f(ϵ_𝐤\mu ^{})f(ϵ_𝐪\mu ^{})}{ϵ_𝐩+ϵ_𝐤ϵ_𝐪ϵ_{𝐩+𝐤𝐪}}.\end{array}$$
(9)
The expression for $`\mathrm{\Omega }^{}`$ obtained in Ref kanno3  is essentially the same as eq (9).

*Comparison & remarks*.
The parameter $`\mu ^{}(T_\mathrm{c})`$ corresponds to the infrared cutoff used in the the non-selfconsistent regulariztion of our\_f02 . Eq (8) gives the same result for the parameter $`\mu ^{}(T_\mathrm{c})`$ as that obtained in our\_f02  for the chemical potential $`\mu `$ using $`\mu =\mathrm{\Sigma }(0,0)`$. The reason for this is that $`F/N_0`$ represents the energy of a zero-momentum quasiparticle with respect to the chemical potential, i.e., $`\mathrm{\Sigma }(0,0)\mu `$. Therefore eq (8) is equivalent to the usual condition $`\mathrm{\Sigma }(0,0)=\mu `$ for the Bose-condensation point.

The shift $`\mathrm{\Delta }T_\mathrm{c}`$ obtained using eqs (7), (8), (9) is the same as that obtained by our non-selfconsistent regularization scheme in our\_f02 .

The distinctive idea in Kanno’s work is the introduction of a modified thermodynamic function $`\mathrm{\Omega }^{}`$ such that $`\mathrm{\Omega }^{}/\mu ^{}`$ replaces the right side in the criticality condition $`N=(V/\lambda ^3)\zeta (\frac{3}{2})`$. This bypasses the issue of self-consistency and the modification of the spectrum at $`T_\mathrm{c}`$, and the use of the chemical potential as a cutoff folllows naturally.

## IV Summary, Further Calculations

To summarise, this article reviews aspects of the calculation of the $`T_\mathrm{c}`$ shift due to a small repulsive interaction in a Bose gas. This question has generated a number of calculations involving diverse techniques; the present discussion is a first step toward developing an organized overview of the different approaches.

The issue that needs to be resolved urgently, in my opinion, is the possibility of a logarithm-like correction to the leading shift in $`T_\mathrm{c}`$. Presumably, a better understanding of the spectrum at the critical point (both at the infrared and the ultraviolet) should help settle this issue. Two directions need to pursued in this regard.

First, as suggested by Stoof, an improved RG calculation, taking into account the momentum dependence of the effective interaction, would improve the ultraviolet behavior of the spectrum in his formalism. Also, it might be worth retaining the wavefunction renormalization (critical index $`\eta `$), which changes the infrared spectrum at $`T_\mathrm{c}`$ from $`k^2`$ to $`k^{2\eta }`$. This is neglected in stoof2  on the grounds that $`\eta `$ is small.

A second approach suggests itself when we note that the quasiparticle treatments our\_f02 ; baym\_bigpaper ; baym\_prl  have all calculated self-energies in terms of the two-body $`t`$-matrix $`t=4\pi a/m`$. This quantity describes all possible collisions between *two* particles, without taking into account the fact that the surrounding medium has an effect on these collisions. Many-body corrections arising from the surrounding gas can be treated by using the many-body $`T`$-matrix, which is the sum of ladder diagrams. The two-body $`t`$-matrix has no temperature-dependence, while the temperature-dependence of the effective potential (vanishing at $`T_\mathrm{c}`$) is captured by the many-body $`T`$-matrix stoof\_variational ; stoof2 ; shi+griffin . Stoof suggests stoof\_private  that this might explain why the quasiparticle analysis (perturbative treatment of spectrum modification) over-adjusts the infrared spectrum.

Determining the self-energies and quasiparticle spectra in terms of the many-body $`T`$-matrix has been attempted our\_unpubl , using the relations provided by Shi & Griffin shi+griffin  between the bare potential $`U`$, the $`t`$-matrix, and the many-body $`T`$-matrix. Our attempt was unsuccessful due to the difficulty of disentangling multiple frequency and momentum inter-dependences at second-order.

Finally, I mention some other work that would be interesting to the $`T_\mathrm{c}`$ community.

* One possibility is to use finite-size scaling ideas, in conjunction with second-order perturbation theory in the grand-canonical ensemble. This would be a complementary to the calculation of Mueller & Baym mueller\_baym\_finite , which uses finite-size scaling with first-order perturbation theory in the canonical ensemble.
* It would be nice to see Schakel’s calculation of the thermodynamic potential, just below $`T_\mathrm{c}`$, done with the spectrum treated more properly, either with a Bogoliubov-like spectrum or with the $`E_𝐤(T_\mathrm{c})k^{2\eta }`$ spectrum.
* As mentioned in section II, a detailed comparison is required between the $`\delta `$-expansion work delta  and other approaches to $`T_\mathrm{c}`$.

###### Acknowledgements.

Informative discussions with H.T.C. Stoof and F. Laloë were made possible by the generosity of ECT\*, Trento (Italy). Adriaan M.J. Schakel answered many e-mail questions. My PhD advisor, A. E. Ruckenstein, introduced me to the problem, and the calculations in our\_f02  were performed under his guidance.