text:
# The 1-point PDF of the Initial Conditions of our Local Universe from the IRAS PSC redshift catalogue

## 1 Introduction

In the gravitational instability scenario, the large-scale structure of the Universe is formed by the growth of small fluctuations superimposed on a homogeneous background. In most inflationary models, these fluctuations are usually predicted to follow Gaussian statistics (see, e.g., Linde, 1990). Non-Gaussian initial conditions are however predicted by models based on topological defects (see, e.g., Brandenberger 1998), as well as by some inflationary models (see, e.g., Salopek 1999). At the present time, perturbations at galactic or smaller scales are highly non-linear, while the large-scale distribution of galaxies is still in the linear or mildly non-linear regime. For these fluctuations, the mapping of the mass distribution from the final to the initial configuration is single-valued. The gravitational evolution of the matter field can be inverted, even though at relatively small scales the evolution is beyond the linear stage and strong non-Gaussian features in the density field are already present.

Under the assumption that galaxies trace the large-scale distribution of mass in a specified way, it is possible to recover the initial perturbation field which has given rise to the observed galaxy distribution. It is interesting to try to recover the initial conditions of our local Universe, not simply for understanding our local cosmography, but for other reasons as well. The inversion itself does not require any assumption concerning the detailed statistical properties of the density field. A reconstruction of initial conditions can therefore be used to test the assumption of Gaussianity. Furthermore, it can be used in N-body simulations to mimic the growth of clustering in our local Universe. This type of simulations have many uses, e.g. to produce predictions of our local peculiar velocity field, to quantify Malmquist biases in peculiar velocity measurements etc. (see, e.g., Kolatt et al, 1996; Narayanan et al. 1999).

The reconstruction of initial conditions is hampered by the well-known facts that galaxies are biased tracers of the matter density, and that the details of this bias are poorly known (see, e.g., Tegmark & Peebles 1998; Catelan et al. 1998; Lahav & Dekel 1999; Pearce et al. 1999; Benson et al. 1999; Somerville et al. 1999; Seljak 2000; Sigad, Branchini & Dekel 2000). A model-independent inclusion of bias in a reconstruction algorithm is not possible at present. Besides, for initial Gaussian perturbations the non-Gaussian statistics of the matter density field is determined by gravitational evolution. Galaxy bias influences the statistics of the galaxy density field. A reconstruction procedure applied to a biased galaxy density field but performed assuming no bias can correct only for gravitationally induced non-Gaussianities, thus letting bias induce non-Gaussian features in the recovered initial conditions. In other words, there is a degeneracy between initial non-Gaussianity and the effects of galaxy bias. The Gaussianity of initial conditions is much better tested on very large scales by the cosmic microwave background (CMB). A careful analysis of the COBE data has lead to a detection of small non-Gaussianity (Ferreira, Magueijo & Gorski 1998), which is likely to be an instrumental effect (Bromley & Tegmark 1999; Banday, Zaroubi & Gorski 1999; Contaldi et al. 1999). Future measurements of Boomerang, MAP and PLANCK will either detect non-Gaussianity or strengthen the limits on it. Such possible non-Gaussian features are anyway small. With plausible extrapolation, they are likely to remain considerably smaller, at the $``$10 $`h^1\mathrm{Mpc}`$ scale, than the signal that galaxy bias can induce in a reconstruction scheme as the one presented here (see Verde et al. 1999). As a consequence, the statistics of the reconstructed initial conditions can be used to give new valuable constraints on galaxy bias models, once the Gaussianity of initial conditions is either assumed or tightly constrained by the future CMB measurements.

The mildly non-linear evolution of the matter field is well modeled by the Zel’dovich (1970) approximation. It describes the evolution of the growing mode of perturbations into the non-linear regime, as long as the flow is still laminar. This restriction to the growing-mode dynamics gives an advantage over a backward N-body simulation, in which noise is amplified as a decaying mode. An inversion of the Zel’dovich mapping was proposed by Nusser & Dekel (1992). Their ‘time machine’ code requires the velocity potential field in real-space as an input. When applied to a galaxy catalogue, the initial conditions are not self-consistently obtained from the observed quantity, i.e. the galaxy density field in redshift space. Other methods for reconstructing the initial conditions, related to that of Nusser & Dekel (1992), were proposed by Gramann (1993) and Taylor & Rowan-Robinson (1993). Croft & Gatzañaga (1997) developed PIZA, a method based on the least-action principle (Peebles 1989) and the Zel’dovich approximation. The same least-action principle has been exploited to construct efficient reconstruction algorithms by Nusser & Branchini (1999) and by Golberg and Spergel (1999). Weinberg (1992) and Narayanan & Weinberg (1998) imposed Gaussianity of the recovered initial conditions to improve the performance of a reconstruction algorithm applied to evolved Gaussian fields. Narayanan & Croft (1999) compared the performance of different reconstruction algorithms to N-body simulations. According to them, the most accurate reconstruction schemes are that of Narayanan & Weinberg (1998) and PIZA.

The correctness of the reconstruction of the initial conditions relies on the accuracy of the input galaxy sample, and on the extent of the sampled volume. In this regard, the uniformity of the galaxy selection and the sky coverage are critical aspects. For this reason, the galaxy samples used to construct the galaxy density field have been selected from the objects contained in the all-sky point source catalogue (PSC) of the IRAS satellite, because they are observed on the whole sky by the same instrument, and IR-selected samples are less affected by dust extinction due to the Galactic Plane (although good all-sky optical catalogues as the ORS (Santiago et al. 1995) or the NOG (Marinoni et al. 1999) are now available). The far-IR selection is biased against spheroidal objects, so that IRAS galaxies tend to miss rich structures (even though at scales larger than $`5`$ $`h^1\mathrm{Mpc}`$ this effect is small, see Baker et al. 1999). The QDOT redshift catalogue (Lawrence et al. 1999) was defined by sparsely sampling the PSC galactic objects with 60$`\mu `$ flux larger than 0.6 Jy (the completeness limit). Fisher et al. (1995) defined their sample by limiting the PSC to a 60$`\mu `$ flux of 1.2 Jy. More recently, a redshift catalogue of all the PSC galaxies (hereafter PSCz; Saunders et al. 2000) has been completed. This catalogue provides the best presently available sample to quantify the galaxy density field of our local Universe.

To reconstruct the initial conditions of our local Universe, Nusser, Dekel & Yahil (1995) used the real-space density field of Yahil et al. (1991), obtained from the IRAS 1.2 Jy redshift catalogue, smoothed with a Gaussian filter of width 10 $`h^1\mathrm{Mpc}`$, using an iterative technique based on linear theory. The peculiar velocity potential, obtained from the real-space density field by means of an average non-linear velocity-density relation, was given as an input to the Nusser & Dekel’s (1992) time machine. The 1-point probability distribution function (hereafter PDF) of their recovered initial conditions was consistent with a Gaussian. A reconstruction of the initial conditions was performed also by Kolatt et al. (1996), using the 1.2 Jy catalogue, and Narayanan et al. (1999), using the PSCz catalogue. In both papers the Gaussianity of the initial conditions was not found from data but forced, following the Gaussianization technique of Weinberg (1992) or Narayanan & Weinberg (1999). The recovered initial conditions were used to run N-body simulations to construct mock peculiar velocity catalogues (Kolatt et al. 1996) or to study the effect of biased galaxy formation in the reconstruction (Narayanan et al. 1999).

Monaco & Efstathiou (1999, hereafter ME99) have recently described an algorithm (ZTRACE) for recovering the real-space density field, the peculiar velocities and the initial conditions from a redshift catalogue with known selection function. ZTRACE is based on a self-consistent solution of the Zel’dovich approximation: it finds, with an iterative scheme, the set of initial conditions which evolve into the observed redshift-space density field, under the assumption that galaxies trace mass. The algorithm has been tested using N-body simulations. When the density field is smoothed with a mass-preserving adaptive scheme, the initial conditions and their PDF are reconstructed in an unbiased way, with the exception of the high peaks whose density contrast, linearly extrapolated to the present time, is $`1`$. These have already gone into the highly non-linear regime. As shown by ME99, the quality of the reconstruction is almost independent of the assumed background cosmology.

In this paper we apply the ZTRACE algorithm to the PSCz catalogue. In this way we recover the initial conditions of our local Universe on scales as small as $``$5 $`h^1\mathrm{Mpc}`$, thus pushing the reconstruction to the limit beyond which high non-linearity prevents any recovery. Here we focus on the reconstructed PDF of the initial conditions, testing its Gaussian shape; future work will address the use of the reconstructed initial conditions to simulate our local Universe. Section 2 describes the PSCz catalogue and the construction of the redshift-space density field. Section 3 gives a brief description of ZTRACE, and describes the results of its application to the PSCz catalogue. In Section 4 we analyze some sources of bias in the reconstructed initial PDF, namely the high-density tail, the effect of galaxy bias and the errors in the selection function. Section 5 describes the reconstructed initial PDF of the PSCz catalogue, and Section 6 gives the conclusions.

Throughout this paper distances are given in $`h^1\mathrm{Mpc}`$, with $`h=H_0`$/(100 km/s/Mpc).

## 2 The PSCz catalogue

For this analysis we use the recently completed PSCz catalogue, which is described in detail in Saunders et al. (2000). The catalogue is based on the IRAS PSC (Joint IRAS Science Working Group 1988), and comprises all galaxies with 60$`\mu `$ fluxes $`>0.6`$ Jy. As the PSC is not complete to 0.6 Jy in 2HCON areas (Beichman et al. 1988), the 1HCON detections recovered from the PSC Reject file were added to the catalogue. Stars, other galactic sources and AGNs were eliminated by means of colour criteria, while residual contamination was eliminated by a combination of optical identifications, examination of the raw IRAS addscan profiles, examination of the Simbad database and, where still unclear, spectroscopic analysis. Redshifts for $``$ 4600 galaxies not belonging to the QDOT (Lawrence et al. 1999) or the 1.2 Jy (Fisher et al. 1995) surveys, and not available in literature, were observed with several telescopes. Faint galaxies, with $`B_J>19.5`$, were not observed; most of these objects are galaxies with $`z>0.1`$, and their omission is irrelevant for the present analysis.

The final catalogue is 98 percent complete in redshift over 84 percent of the sky. The missing zones are the Galactic Plane (where the extinction is $`A_V>1.42`$ mag), the areas covered by the LMC and SMC, and two small strips not observed by IRAS. Fig. 1 shows the Aitoff projected catalogue, with the uncovered zones masked in grey. The PSCz catalogue has already been used in other works on the topology of the large-scale structure (Canavezes et al. 1998), on the cosmological dipole (Rowan-Robinson et al. 1999; Schmoldt et al. 1999a), on the beta parameter (Schmoldt et al. 1999b), on the power spectrum of fluctuations (Sutherland et al. 1999), on the redshift distortions (Tadros et al. 1999) and on the peculiar velocity field (Branchini et al. 1999). Comparisons with other galaxy and galaxy cluster catalogues were presented by Seaborne et al. (1999) and Plionis et al. (1999). Reconstruction analyses were presented by Sharpe et al. (2000), based on the least-action principle, and by the already mentioned paper by Narayanan et al. (1999).

To calculate the galaxy density field from the PSCz catalogue it is necessary to know the selection function, which gives the number density of objects bright enough to be included in the catalogue, as a function of distance. The selection function is well fitted by the following functional form:

$$\mathrm{\Phi }(r)=\mathrm{\Phi }_{}\frac{(r/r_{})^{1\alpha }}{(1+(r/r_{})^\gamma )^{\beta /\gamma }},$$
(1)
where $`r=cz/H_0`$ is the redshift distance, but the volume is corrected for relativistic effects (assuming an Einstein-de Sitter Universe). As ZTRACE applies in a frame comoving with the background, we use redshifts in the CMB frame. A maximum likehood fit, performed with the method of Mann, Saunders & Taylor (1996) and used in Saunders et al. (2000), gives the following best values for the parameters: $`\mathrm{\Phi }_{}=0.0224h^3Mpc^3`$, $`r_{}=64.10`$ $`h^1\mathrm{Mpc}`$, $`\alpha =1.41`$, $`\beta =4.32`$ and $`\gamma =1.39`$. The resulting selection function is only marginally different from that calculated in the Local Group frame. The effects of uncertainties in the selection function will be addressed in Section 4.3.

## 3 Application of ZTRACE to the PSCz catalogue

ZTRACE is based on a self-consistent solution of the Zel’dovich (1970) approximation. In the Lagrangian picture of fluid dynamics a mass element with initial (Lagrangian) comoving position $`𝐪`$ is displaced to the final (Eulerian) position $`𝐱`$ through a displacement field $`𝐒`$:

$$𝐱(𝐪,t)=𝐪+𝐒(𝐪,t).$$
(2)
The map $`𝐒`$ is the dynamical variable in this approach. For small displacements, and for irrotational fluids, the first-order growing mode gives the Zel’dovich (1970) approximation:

$$𝐱(𝐪,t)=𝐪D(t)\phi (𝐪).$$
(3)
Here $`D(t)`$ is the linear growing mode (see, e.g., Peebles 1980); for an Einstein-de Sitter universe it is equal to the scale factor $`a(t)`$. The growing mode is normalized so that $`D(t_0)=1`$, where $`t_0`$ is the present time. $`\phi (𝐪)`$ is a rescaled version of the peculiar gravitational potential, such that, at an initial time $`t_i`$, $`^2\phi =\delta (𝐪,t_i)/D(t_i)`$. The density contrast $`\delta `$ is defined as usual, $`\delta (𝐱)=(\varrho (𝐱)\overline{\varrho })/\overline{\varrho }`$, where $`\overline{\varrho }`$ is the background density. It is useful to define the quantity:

$$\delta _l(𝐪)\delta (𝐪,t_i)/D(t_i).$$
(4)
$`\delta _l`$ is constant in linear theory, and is equal to the initial density contrast linearly rescaled to the present time. This quantity will be called the linear density contrast in the following.

In a redshift survey, the radial coordinate is the redshift distance, not the real distance. It is possible to generalize the Zel’dovich map, Eq. 3, to the redshift space $`𝐬`$:

$$𝐬(𝐪,t)=𝐪D(t)\left[\phi (𝐪)+f(\mathrm{\Omega })\left(\phi (𝐪)\widehat{𝐱}\right)\widehat{𝐱}\right].$$
(5)
Here $`\widehat{𝐱}`$ is the versor of $`𝐱`$, and $`f(\mathrm{\Omega })=d\mathrm{ln}D/d\mathrm{ln}a`$ is usually approximated by $`f(\mathrm{\Omega })\mathrm{\Omega }^{0.6}`$ (Peebles 1980), or by $`f(\mathrm{\Omega },\mathrm{\Omega }_\mathrm{\Lambda })\mathrm{\Omega }^{0.6}+\mathrm{\Omega }_\mathrm{\Lambda }(1+\mathrm{\Omega }/2)/70`$ (Lahav et al. 1991) if the cosmological constant is non-zero. Notice that the last term in the right-hand-side of Eq. 5 is the peculiar velocity along the line of sight.

The density contrast follows from the map in Eq. 3 or Eq. 5 through:

$$1+\delta (𝐪)=[det(\delta _{ab}^K+S_{a,b})]^1.$$
(6)
Here $`\delta _{ab}^K`$ is the Kronecker tensor, and $`S_{a,b}`$ is the derivative of the map $`𝐒`$. ZTRACE is designed to invert the non-linear and non-local set of equations which connect the initial conditions to the redshift-space density, under the hypothesis of laminar flow, i.e. when the Zel’dovich map in redshift space, Eq. 5, is still single-valued. Multi-valued regions appear as soon as different mass elements which lie along the same line of sight are observed at the same redshift. This happens when the perturbation decouples from the Hubble flow, giving rise to multi-stream regions in the distance-redshift relation (often referred to as “triple-valued regions”). The linear density contrast has a value of order one in such regions; ZTRACE is not able to reconstruct the high peaks of the density distribution.

ME99 describe in detail those techniques used to help the convergence of ZTRACE for smoothed density fields with standard deviations $`\sigma `$ up to 0.7111The tests performed by ME99 are limited to $`\sigma 0.7`$. Anyway, for larger $`\sigma `$-values the distortion on the positive tail of the initial PDF becomes very strong.. The application of ZTRACE to the PSCz redshift catalogue can be summarized as follows.

(i) The masked regions are filled with a synthetic catalogue, so as to produce the galaxy density field over a large and connected spherical volume. We have used the procedure developed by Yahil et al. (1991) as modified by Branchini et al. (1999). According to this procedure, the zone of avoidance of the Galactic disk, with galactic latitude $`|b|<8`$ degrees, is divided in longitude and redshift bins. The galaxy number density in each bin is obtained by interpolating between those of the two 8 degree stripes above and below the bin. A synthetic galaxy catalogue is then obtained by randomly sampling the density field within the bins. In this way, the radial distribution follows the same PSCz selection function. The other masked regions of the sky which are outside the zone of avoidance are filled at random, according to the PSCz selection function. This algorithm is considered a better guess with respect to filling at random the masked regions, as it tends not to break the coherent structures, such as the Perseus-Pisces supercluster, which happen to cross the Galactic Plane. As a further test, we have also applied the lognormal interpolation procedure described in Saunders et al. (1999); the results are very similar and are not described here. The reconstruction method is robust with respect to the interpolation algorithm because only a small fraction of the sky is masked.

(ii) The relaxed groups are collapsed into spheroids. The effects of the highly non-linear dynamics are apparent in the redshift space as elongated structures along the line of sight, often referred to as ‘fingers of God’. As pointed out by Gramann, Cen & Gott (1994), a reconstruction scheme improves if such non-linear structures are collapsed before application of the algorithm. In ME99 the collapsing procedure was found to give a modest increase in the performance of ZTRACE. Following ME99, we have found the groups in the PSCz catalogue by applying a standard friends-of-friends algorithm, with radial and tangential linking lengths of 3 and 0.5 $`h^1\mathrm{Mpc}`$.

(iii) The catalogue is smoothed to generate a density field in redshift space. The performance of ZTRACE improves considerably if the density field is obtained with an adaptive, mass-preserving filtering of the galaxy catalogue. This is similar to a smoothing in the Lagrangian space, and correctly takes into account that overdensities come from larger patches which have contracted, while the opposite is true for underdensities. Following ME99, we have constructed the redshift-space galaxy density field by smoothing the PSCz point distribution (with the addition of the synthetic catalogue to fill the masked regions) with the adaptive smoothing code of Springel et al. (1998). This code requires the definition of a reference smoothing radius, which corresponds to the case of null density contrast; the adaptive refinements are then performed so that the smoothing filter always contains the same number of objects. As in ME99, the reference smoothing radius $`R`$ is held constant within a given distance $`D`$, to a value such that 10 galaxies are on average contained in the filter at that distance. At larger distances, the reference smoothing radius is increased according to the selection function, so that 10 galaxies are contained on average within the filter. In this way, the outer parts that are more and more severely smoothed are used by ZTRACE to give the external tides to the inner region.

The density field is calculated on a $`64^3`$ grid of size 240 $`h^1\mathrm{Mpc}`$. Fig. 2a shows the obtained density field in Super-Galactic coordinates, in the case $`R`$=7.6 $`h^1\mathrm{Mpc}`$ within $`D`$=80 $`h^1\mathrm{Mpc}`$; the severe smoothing of the outer regions is apparent. The high peaks in the density distribution correspond to well-known structures as the Local Supercluster \[at the position $`(x,y)(0,0)`$\], the Great Attractor region \[the two peaks at $`(x,y)(40,20)`$ and $`(40,20)`$\], the Perseus-Pisces supercluster \[the two peaks at $`(x,y)(20,50)`$ and $`(40,20)`$\] and its appendix beyond the Galactic Plane \[the Camelopardalis supercluster, at $`(x,y)(60,10)`$\], and the Coma supercluster \[at $`(x,y)(0,80)`$\] (see, e.g., Branchini et al. 1999).

(iv) The input density field is padded in the volume outside the largest inscribed sphere (of radius 120 $`h^1\mathrm{Mpc}`$), so as to have periodic boundary conditions. This is done so that we can apply FFT’s to solve the equations. ZTRACE is applied to the field; convergence is achieved within $``$12 iterations. Fig. 2b shows the resulting initial conditions in the comoving real space (in terms of the linear density contrast $`\delta _l`$, Eq. 4). (Notice that the linear density contrast is not bound to be $`1`$). The main superclusters are originated from high peaks in the initial conditions. These peaks appear flattened because of the inability, discussed above, of ZTRACE to reconstruct the high-density regions. The different statistical properties of the initial conditions, with respect to the density field in real and redshift space, are apparent.

(v) The PDF of the initial conditions is calculated. In this analysis only the grid points within the distance $`D`$ are considered. The points in the zone of avoidance ($`|b|<8^{}`$) are excluded, so as the central 10 $`h^1\mathrm{Mpc}`$ sphere. This last exclusion is motivated by the inability of ZTRACE to converge in the neighborhood of the observer; the numerical solution tends to oscillate at the origin, where the map given in Eq. 5 is singular. To avoid this the density field is damped over $``$5 $`h^1\mathrm{Mpc}`$. The exclusion of the central region is not a problem as its volume is small, and the PSCz is known to be incomplete at such small distances (Rowan-Robinson et al. 1991).

## 4 errors and Biases in the reconstruction of the initial PDF.

The recovery of the initial PDF can be influenced by high non-linearities, by galaxy bias and by errors in the selection function. As shown in ME99, a wrong assumption on the background cosmology does not influence significantly the shape of the initial PDF.

### 4.1 The high-density tail

As mentioned above, ZTRACE is not able to recover the high peaks, which have decoupled from the Hubble flow. This induces a bias in the initial PDF, whose high density tail is suppressed. The underestimated high-density points influence the initial PDF also at moderate densities. It is useful to quantify this effect before addressing the Gaussian nature of the initial PDF reconstructed from the PSCz catalogue. To this aim, we have extracted mock PSCz catalogues from an N-body simulation of a standard CDM Universe (with $`\mathrm{\Omega }=1`$ and $`\mathrm{\Gamma }\mathrm{\Omega }h=0.5`$) with normalization $`\sigma _8=0.7`$. The simulation, performed with the Hydra code (Couchman, Thomas & Pearce 1995), has already been presented in ME99. It is performed with $`128^3`$ dark matter particles in a box of length 240 $`h^1\mathrm{Mpc}`$ with a $`256^3`$ base mesh with adaptive refinements. As the volume used in the reconstruction is much smaller than the total volume of the simulation, it is convenient to extract more than one simulated PSCz catalogue from the simulation. We have extracted 10 catalogues, centred on random points, under the assumptions that galaxies trace mass and applying the PSCz selection function (in real space). It has been checked in ME99 that the selection in real space with a selection function determined in redshift space does not influence the result. The smoothing radius has been set to $`R=5.1`$ $`h^1\mathrm{Mpc}`$ (constant within $`D=50`$ $`h^1\mathrm{Mpc}`$) or $`R=7.6`$ $`h^1\mathrm{Mpc}`$ (constant within $`D=80`$ $`h^1\mathrm{Mpc}`$). When the smoothing is held constant within $`D=50`$ $`h^1\mathrm{Mpc}`$, the volume actually used in the reconstruction is 3.8 percent of the total volume of the simulation, so that the different realizations are accurately statistically independent. When the smoothing is held constant within $`D=80`$ $`h^1\mathrm{Mpc}`$, the volume used is 15.5 percent, and then the 10 realizations are oversampling the volume. This is approximately taken into account by calculating the variance among the various realizations as if they were $`1/0.1556.45`$ instead of 10.

For each simulated catalogue, the galaxies in the PSCz masked area are removed, then replaced using the same algorithm of Branchini et al. (1999) discussed above. We have checked that the masking procedure introduces noise but no bias in the reconstruction. The reconstruction proceeds as described in Section 3.

Because of sample variance, the reconstructed initial PDFs differ in general in mean ($`\delta _l`$) and standard deviation ($`\sigma _l`$). To construct an average initial PDF it is then useful to consider the PDF of the quantity $`(\delta _l\delta _l)/\sigma _l`$. All the distributions are binned in the same way, and averages and variances are calculated for each bin. Fig. 3 shows the reconstructed initial PDFs, compared with the true Gaussian curve, for the two smoothing radii used. The upper panels show the reconstructed initial PDFs, the errors are the square root of the variance of the mean for the ten catalogues (calculated as if they were 6.45 in the second case, see above). The lower panels show the significance of the difference between the reconstructed and true initial PDF. It is apparent that the true Gaussian shape is recovered accurately, except in the high-density tail.

### 4.2 The effect of galaxy bias.

The simulated galaxy catalogues are created and reconstructed assuming that galaxies trace mass. This may not be true in practice, as galaxies may be biased tracers of the underlying mass distribution. As mentioned in the Introduction, galaxy bias may induce non-Gaussian distortions in the initial PDF, when this is reconstructed simply assuming that galaxies trace mass. To quantify this effect, it is convenient to apply ZTRACE, under the assumption that galaxies trace mass, to biased galaxy catalogues that are simulated assuming some specific bias model.

We have used a (deterministic, i.e. non-stochastic) power-law local bias scheme, already suggested, e.g., by Nusser et al. (1995) and by Lahav & Dekel (1999):

$$1+\delta _g(𝐱)=(1+b_0)(1+\delta (𝐱))^{b_1}.$$
(7)
Here $`\delta _g`$ is the galaxy density contrast. This scheme respects the condition $`\delta _g>1`$. For small density values this scheme reduces to the linear bias, $`\delta _g=(1+b_0)b_1\delta `$. $`b_0`$ is constrained so that $`\delta _g=0`$. In the following we will consider the anti-biased case $`b_1=0.5`$ and $`b_0=0.084`$, and the biased case $`b_1=1.5`$, and $`b0=0.137`$. Fig. 4 shows the two bias curves.

A more physically motivated bias scheme can be obtained by using semi-analytic galaxy formation models (see, e.g., Benson et al. 1999; Somerville et al. 1999). In these models, the assembly of dark matter halos is described by means of N-body simulations, while the physics of the baryonic component is inserted in the halos through a set of simple analytic rules. Galaxy bias can then be obtained by comparing the density of galaxies with the density of dark matter. Narayanan et al. (1999) give an approximate bias function for the IRAS-like galaxies of Benson et al. (1999):

$$1+\delta _g(𝐱)=A(1+\delta (𝐱))^\alpha \left[C+(1+\delta (𝐱))^{(\alpha \beta )/\gamma }\right]^\gamma .$$
(8)
Here $`\alpha =2.9`$ and $`\beta =0.825`$ are the asymptotic logarithmic slopes at small and large densities, while $`A=1.1`$, $`C=0.08`$ and $`\gamma =0.4`$. This bias function is shown in Fig. 4. Compared to the simple power-law curves, it behaves as in the biased case in underdensities, but it is qualitatively more similar to the unbiased case in overdensities. We use the semi-analytic bias scheme as an exercise to quantify the distortion of the initial PDF expected from a more sophisticated bias function. It must be noticed, however, that this bias function is relative to a flat Universe with $`\mathrm{\Omega }_\mathrm{\Lambda }=0.7`$ and $`\sigma _8=0.9`$, and to a density field smoothed with a top-hat filter with radius 3 $`h^1\mathrm{Mpc}`$, while we are applying it to an Einstein-de Sitter universe and to a density field Gaussian smoothed over 5 $`h^1\mathrm{Mpc}`$.

Following the same procedure as in Section 4.1, we have simulated three sets of 10 catalogues, with the three bias functions discussed above and shown in Fig. 4. The bias schemes are applied to the real-space density field of the simulation, smoothed with a Gaussian filter of width 5 $`h^1\mathrm{Mpc}`$. Fig. 5. shows the obtained average initial PDFs; they are to be compared to those shown in Fig. 3. The Gaussian initial PDF is again well reconstructed in the “anti-biased” power-law case in which $`b_1=0.5`$, even though the low-density tail shows a 2-$`\sigma `$ discrepancy when $`R=5.1`$ $`h^1\mathrm{Mpc}`$. In the “biased” power-law case, that in which $`b_1=1.5`$ and especially the semi-analytic case, significant distortions are apparent. For $`R=5.1`$ $`h^1\mathrm{Mpc}`$, the tail at low densities is overestimated, a dip and a bump are present respectively at $`\delta _l/\sigma 0`$ and 1.5. The same features are still present in the case $`R=7.6`$, although they are less significant.

Fig. 5 shows that the application of ZTRACE with the no-bias assumption can lead to small but significant distortions in the reconstructed initial PDF. However, these distortions are revealed by averaging over ten realizations of the PSCz catalogue, but they are not statistically significant for single biased realizations. The PSCz survey is not large enough to reveal such deviations, but larger surveys such as the 2dF and SDSS surveys may be able to provide constraints on non-linear bias schemes. Of course this conclusion may change if stronger bias functions are considered; however, evidence on IRAS galaxies suggests that they are weakly biased with respect to the matter density (see, e.g., Saunders et al. 1999).

### 4.3 Errors in the selection function

Another potential source of distortion on the initial PDF comes from errors in the selection function. Given the highly correlated nature of the density field, a systematic error in the density estimate in a spherical shell which happens to contain some large voids or clusters can influence the shape of the initial PDF. To test for this, we have collected the selection functions (in the Local Group frame) used by ME99, Sutherland et al. (1999), Canavezes et al. (1998), Branchini et al. (1999), and Springel (1996). We have added to these our own determination of the selection function in the Local Group frame, performed as described in Section 2. The differences between the various functions, which are of order of 10 percent, are due to the slightly different versions of the catalogue and mask used by the various authors, and to differences in the fitting procedure. Fig. 6 shows the residuals of these selection functions with respect to the one used in this paper (Eq. 1). Notice that all the residuals show the same systematics because the selection functions are relative to the Local Group frame. Then, the effect of errors in the selection function is overestimated in this analysis.

The first simulated PSCz catalogue (with no bias) has been reconstructed assuming either the correct or one of the six other selection functions listed above. Fig. 7 shows the seven reconstructed initial PDFs, compared with the average reconstructed initial PDF already shown in Fig. 3. Differently from Fig. 3, the errors shown are the (square root of the) variance between the various realizations (instead of the variance of the mean). This is the correct error to use when assessing the difference between a single realization and the average initial PDF. Again, the lower panels show the significance of the disagreement. The various initial PDFs are hardly distinguishable, no significant distortion is induced. We conclude that the uncertainty in the selection function does not influence significantly the shape of the reconstructed initial PDF.

## 5 The reconstructed initial PDF of the PSCz catalogue

The PSCz catalogue has been adaptively smoothed according to the smoothing strategy described in Section 3. The reference radius has been held constant, within $`D=50`$ and 80 $`h^1\mathrm{Mpc}`$, to the values of $`R=5.1`$ and 7.6 $`h^1\mathrm{Mpc}`$. Fig. 8 shows the reconstructed initial PDFs of the PSCz catalogue compared with the results of the simulated catalogues described in Section 3.1 (with no bias and CMB selection function) and shown in Fig. 3. The error associated to the simulated initial PDFs is the square root of the variance between the realizations.

The value of $`\sigma _8`$ for the PSCz density field turns out to be about 0.8 (Sutherland et al. 1999; Seaborne et al. 1999), slightly larger than the value 0.7 used in our test simulation. As a consequence, the variance of the reconstructed initial PDF is slightly larger than the test one smoothed with the same radius. This has some effect on the spurious high-density cutoff of the PDF induced by ZTRACE. To compare reconstructed PDFs that have roughly the same variance, the PSCz has been smoothed over $`R=7.2`$ and 10.1 $`h^1\mathrm{Mpc}`$ within $`D=75`$ and 105 $`h^1\mathrm{Mpc}`$. Fig. 9 shows the comparison with the same PDFs of simulated catalogues as in Fig. 3 and 8.

The reconstructed initial PDF of the PSCz catalogue disagrees significantly with the simulated one in Fig. 8a (smoothing radius $`R=5.1`$ $`h^1\mathrm{Mpc}`$), where the standard deviation $`\sigma `$ of the two PDFs that are compared are different. Even in this case, significant distortions are visible only on the positive-density side, for $`(\delta _l\delta _l)/\sigma 1`$. In this range of densities the influence of the high-density cut (which is different in the two cases due to the difference in $`\sigma `$) is strong. This discrepancy vanishes in Fig. 9a, when comparing the reconstructed PDFs of simulated and PSCz catalogues with similar variances. Besides, the distortions on the negative-density side, where we expect to find the signature of galaxy bias (see Fig. 4), are never significant, except a 3-$`\sigma `$ discrepancy which is limited to one single bin. Figs. 8b and 9b show the same comparison for the larger smoothing radii. In this case, no strongly significant distortions are visible, although some modest but persistent discrepancy at the 3-$`\sigma `$ level is present both on the positive- and on the negative-density sides.

We conclude that there is no convincing evidence for non-Gaussianity in the initial PDF reconstructed from the PSCz catalogue. Notably, the overall pattern of the residuals is similar to the distortions induced by biased power-law or semi-analytic bias. Anyway, the statistical significance of these distortions is too marginal to draw any conclusion.

The Gaussianity of the reconstructed inital PDF gives support to the hypotheses of Gaussianity of initial conditions and of modest bias of IRAS objects (as described in Section 4.2). In particular, the low-density tail of the initial PDF remains Gaussian down to fluctuations of $``$2.5 times the standard deviation. The fact that ‘the voids are too devoid’ of bright galaxies has often been seen as a possible problem for CDM-like models of structure formation (see, e.g., Peebles 1999). The result presented here shows that the voids, as defined by IRAS galaxies, are as underdense as they are expected to be as an effect of gravitation, at least on scales as small as 5 $`h^1\mathrm{Mpc}`$, where the highly non-linear regime is still not dominant.

## 6 Summary and Conclusions

The ZTRACE algorithm has been applied to the new IRAS PSCz redshift catalogue to reconstruct the initial conditions of our local Universe. The ZTRACE algorithm, recently proposed by ME99, is known to give an unbiased estimate of the initial density which generates, through gravitational collapse in the Zel’dovich (1970) approximation, the observed galaxy density field in the redshift space under the assumption that galaxies trace mass. The reconstructed initial conditions of our local Universe can be used as input to N-body simulations, with the small-scale modes restored with a Monte-Carlo procedure (Kolatt et al. 1996; Narayanan et al. 1999).

The results presented in this paper push the reconstruction of initial conditions to scales of 5 $`h^1\mathrm{Mpc}`$. This is close to the limit at which highly non-linear dynamics becomes dominant, hampering any detailed inversion of the gravitational evolution. This paper improves by a factor of two in scale over the previous determination of the initial PDF, by Nusser et al. (1995). This improvement is due to the increased sampling density and depth of the new PSCz catalogue relative to previous IRAS-based catalogues such as QDOT (Lawrence et al. 1999) or 1.2 Jy Fisher et al. (1995), and to the superior performance of ZTRACE, compared to previous reconstruction algorithms (see ME99 and Narayanan & Croft 1999).

The reconstructed initial PDF of the local Universe is Gaussian to within the accuracy of the method. The lack of detection of non-Gaussianity in the initial PDF supports the hypothesis that initial conditions are Gaussian and IRAS galaxies are not strongly biased with respect to the matter density. Moreover, it shows that the voids are as devoid of bright galaxies as they should be as an effect of gravity (see, e.g., Peebles 1999). The Gaussian nature of the primordial perturbations has been tested on much larger scales and with greater accuracy by measurements of the CMB. Then this result, besides offering a complementary support for Gaussianity on the scale of $``$5–10 $`h^1\mathrm{Mpc}`$, gives original constraints on the nature of galaxy bias, which in general induces distortions on the initial PDF reconstructed under the no-bias hypothesis. However, our tests show that the distortion induced by realistic power-law galaxy bias models, especially on the low-density tail, is just beyond the reach of the PSCz catalogue. Future surveys, like the 2dF or the SDSS, covering larger volume than PSCz, will provide a way to constrain galaxy bias beyond the linear level from the reconstruction of the initial PDF.

## Acknowledgments

The Hydra team (H. Couchman, P. Thomas, F. Pearce) has provided the code for N-body simulations. Volker Springel has kindly provided a C version of the code for adaptive smoothing. We thank the referee, David Weinberg, for his helpful comments. P.M. thanks Christian Marinoni and Inga Schmoldt for discussions. P.M. has been supported by the EC TMR Marie Curie grant ERB FMB ICT961709.


format: text

text:
# KUCP0154 Pareto’s Law for Income of Individuals and Debt of Bankrupt Companies

## 1 Introduction

More than a century ago the Italian sociologist Pareto studied the distribution of personal incomes for the purpose of characterizing a whole country’s economic status . He found power law cumulative distributions with exponents close to $`1.5`$ for several countries, which turned out to be a classic example of fractal distributions . In 1922 Gini checked the same statistics and reported that power laws actually hold, but the values of the exponents vary from country to country . Montroll and Shlesinger analyzed the USA’s personal income data for the year 1935–36 and found that the top 1% of incomes follow a power law with an exponent $`1.63`$, while the rest, who are expected to be salaried, follow a log-normal distribution . Although these results are interesting and suggestive, the data are all old, and a much more precise analysis using contemporary high quality data in digital form is desirable. One of this paper’s aims is to answer this requirement by showing the results of personal income distribution analysis on the latest high precision data reported to the tax office in Japan. As we will show in section 2 we find a clear power law with an exponent very close to $`2`$, the first time it has been reported in this kind of study.

There is another topic in this paper that is related to personal incomes, but in a juridical sense, namely, company incomes. It was recently reported by Okuyama and Takayasu, who analyzed high quality Japanese company data for the year 1997, that the income distribution of companies accurately follows Zipf’s law, that is, a cumulative distribution with exponent $`1`$ . It is anticipated that this law holds rather universally for different years, for different countries, and even for each job category, with some exceptions. The second aim of this paper is to show a complementary result relating to the company income distributions, that is, the distribution of debts of bankrupt companies, which may correspond to the distribution of negative incomes. Although the data available for this purpose is limited we can apparently demonstrate that the debts distribution also follows the same statistics, Zipf’s law.

In the following section we give the details of our personal income distribution analysis. We briefly describe the result of the debt distribution in section 3, and the final section is devoted to discussions of the implications and significance of our results.

## 2 The income and income-tax distributions

The data we have for income and income tax are in many senses complementary. Both data relate to Japan, and the income data is for the fiscal year 1997 and 1998, while the income-tax data is only for 1998. For the fiscal year 1998, the income data contains all 6,224,254 workers who filed tax returns, but it is a coarsely tabulated data, while the income-tax data lists the income tax of individuals who paid tax of ten million yen or more in the same year. In this section we study the distribution of each data set and then combine them to obtain an over-all picture of income distribution in high income range.

The tabulated data for the income distribution of individuals in Japan is publicly available on the web pages of the Japanese Tax Administration . Since the relevant pages are in Japanese only, we quote the data for the fiscal year 1998 in Table 1, where and hereafter we denote the income by $`s`$ in units of million yen. It should be noted that this data is for individuals who filed a tax return, and therefore does not include all individuals with income. Under the Japanese tax system, if a worker has only one source of income (salary) and the income is less than 20,000,000 yen, that person does not have to file a tax return (the tax deducted from the monthly salary is adjusted at the end of the year by the employer and that becomes the final amount of tax paid). In some circumstances, even if the income is below this amount, a person may have to file a tax return if some other conditions are met. Consequently, a large number of individuals in fact submit returns. Therefore for $`s<20`$ this data represents only a portion of the actual number of individuals, and we expect that the deviation between the number of people in this table and the actual number of people in that income range becomes larger in lower income ranges. Rigorously speaking, it is safe to trust only the entries for income $`s20`$ (the last three columns) in this data as a faithful representation of the number of individuals in that income range.

In Fig. 1 we give the rank-size plot, which is the log-log plot of the rank ($`r`$) as a function of income ($`s`$). The vertical axis may be converted to the accumulated probability $`P(s)`$, the probability that a given individual has income equal to or greater than $`s`$. This is done by simple rescaling $`P(s)=r(s)/r(0)`$ where $`r(0)=6,224,254`$, and therefore it merely induces a translation of the vertical axis in the log plot and does not change the shape of the distribution. The ticks for the accumulated probability $`P(s)`$ are placed on the right-hand side of the plot. It is apparent that this plot tends to a linear function in the high income range, say $`s>10^1`$. As we have noted above, this distribution should be equal to that of all the workers only for $`s20`$. Therefore we have calculated the best-fit linear function from the last three data points and have found that it is given by the following;

$$rs^{1.98}.$$
(1)
The broken line in Fig.1 represents this function. The distribution $`P(s)`$ fits the power law (1) very well in the high income range and it gradually deviates from (1) as the income becomes lower. Since this deviation starts to grow as the income $`s`$ becomes lower than 20, we may assume that the power law (1) applies to a much wider range of incomes.

Another set of publicly available data is a list of all taxpayers who paid income-tax of ten million yen or more in 1998. Hereafter we denote the tax by $`t`$ in million yen. There are 84,515 such individuals listed for fiscal year 1998. The rank-size plot of this data is given in Fig.2. This distribution is very close to being linear in this log-log plot. The best-fit result is given by the following;

$$rt^{2.05}.$$
(2)
This function is given by the dotted line in Fig.2.

The fact that the income-tax distribution satisfies Eq.(2), which is very close to the income distribution Eq.(1), implies that income-tax is roughly proportional to the income. This proportionality constant can be fixed by assuming that the income-tax is a monotonically increasing function of the income, which we expect to be true in an averaged sense. By this assumption a given person has the same rank in both the income listing and the income-tax listing. Therefore the relation between income and income-tax can be found in the overlapping region of those two data. The data of incomes, however, is given only in approximate form by Table 1, while the listing of the income-tax covers the bottom one and a half columns. Therefore, the comparison between the income and the income-tax is possible only for the rank $`r=43,455`$ whose income was $`s=50`$. From our listing of the income-tax we find that the 43,455th ranked person has paid income tax of $`t=15.13`$. Therefore we find that

$$t=0.30s.$$
(3)
In the fiscal year 1998, the tax ($`t`$) was obtained from the income ($`s`$) in three steps: (1) various deductions, such as the basic (default) deduction, deductions for the spouse and the dependents, and insurance deductions, were made to calculate the taxable income ($`s^{}`$), (2) the basic tax ($`t^{}`$) is calculated according to the formula $`t^{}=0.5s^{}6.03`$ in this income range, and (3) further deductions, such as the deduction for donation to political parties, deduction for housing loans and further special deductions are made to arrive at the actual tax $`t`$. Although it is difficult to assess the validity of the result (3), it may be understood as an averaged result of these steps, especially since the proportionality constant is meaningfully less than the constant 0.5 of the tax formula.

Using the result (3), we translate the income-tax data to the income data. This way, the shortcoming of the original income data, that it lacks the information in the high income range, can be overcome. The resulting rank-size plot is given in Fig.3.

Again, the power-law behavior is apparent in Fig.3. Caution is needed in applying the best-fit method; The data points of the tax-induced income becomes very dense at lower $`s`$, therefore they would dominate the result if all data points are treated with equal weight in a simple best-fit method. The data of original income table would certainly be neglected in such an analysis. For this reason, we have calculated the best fit to the power law from each of the two columns in Table 1 and each of the five sections $`[50\times 2^n,50\times 2^{n+1}](n=1\mathrm{}5)`$ of the tax-induced income data and have taken the average of the seven linear functions (in the log-log scale). The result is the following;

$$rs^{2.06}.$$
(4)
The broken line in Fig.3 denotes this function. We find that the power-law (4) fits the data over three magnitude of income, $`10<S<10^4`$.

As mentioned earlier, the deviation of the income data points from the power law in the range $`s<20`$ may be explained as the result of the exclusion of the majority of the workers of single salary source. On the other hand, there is a possibility that the distribution in the lower ranges follow the log-normal distribution as was suggested by Montroll and Shlesinger . Therefore it is highly desirable to investigate these lower ranges of income to obtain the overall profile of the income distribution. One such study toward this direction is currently in progress by one of the authors (W.S.) and will be published in near future .

## 3 Debts of bankrupt companies

The data we have gathered covers some 100 big bankruptcies in Japan from 1997 to the end of March, 2000. Unfortunately the exact criteria for the choice of the “big” cases is not known. The resulting rank-size plot for the accumulated probability $`P(d)`$ is given in Fig.4, where we have denoted the debt in units of million yen by $`d`$. The data actually contains much smaller cases, but from the fact that the number of data starts to decrease for $`d<10^5`$ we assume that the data starts to loose reliability in this range. The broken line represents the Zipf law,

$$rd^1,$$
(5)
which were found for the company income previously . This data is in reasonable agreement with this power-law.

## 4 Discussion

As described in section 2 we have discovered a new type of income distribution for individuals, the power law distribution with an exponent close to $`2.0`$. If we assume that all the pioneering results together with our result, then the only consistent interpretation may be as follows: Income distribution generally follows a power law for large incomes but its exponent varies from country to country and year to year. As noted, Pareto and Gini give different exponents using data from different years. However, there is a possibility that such old results may based on imperfect data, and we may still expect a kind of universal behavior for personal income distributions. In our analysis the observed scale range satisfying the power law in Fig.3 is about three decades in the income rank range ($`10`$ to $`10^4`$), therefore, the estimation of the exponent is quite accurate. We need comparably high precision data for other countries or years for further discussion about such universality.

There has been no established theory for the income distributions. A classical theory was proposed by Gibrat in 1932 . He assumed the time evolution of each person’s income to be approximated by a multiplicative stochastic process considering that a random process proportional to the amount of the present income can approximate the increase and decrease of income. By this assumption the resulting distribution of income follows a log-normal distribution which is consistent with the result of Montroll and Shlesinger for salaried people. However, this theory apparently fails to explain the more interesting part of the distribution, the power law tails.

Explaining the existence of power law tails with the level of simplicity in Gibrat’s theory can easily be done as follows by introducing a simple modification to his model. Let us introduce an additive random term together with the multiplicative term, such as

$$I(t+1)=a(t)I(t)+R(t).$$
(6)
If we assume that $`a(t)`$ and $`R(t)`$ are independent random variables then it is mathematically shown that the distribution of $`I(t)`$ follows a general power law even in a case where the fluctuation of $`a(t)`$ is correlated . In the simplest case when both $`a(t)`$ and $`R(t)`$ are white noises, the power law exponent, $`\beta `$, of the distribution of $`I(t)`$ is simply given by the following formula ,

$$a^\beta =1.$$
(7)
Based on this relation the power law distribution with an anticipated exponent can be obtained by only tuning the statistics of growth rate, $`a(t)`$.

There can be many different approaches other than this simple growth model. One straightforward model can be given by relating the company’s incomes to personal incomes by assuming that the high income individuals are mostly company owners. If we can assume that company owner’s income is proportional to the square root of the company’s income, then the established empirical Zipf’s law for the company’s income reduces the squared power law for individuals.

Nonlinear interaction models of incomes may also be important such as the model introduced by Levy and Solomon . An obvious difficulty in such modeling is the mathematical form of interaction among people working in real economics. For example economical interaction such as buying and selling generally occur in a discrete manner both in time and price, therefore, it is even problematic whether we can introduce an analytical function as an interaction term.

Related to the issue of the income distribution and the dynamics behind it is the distribution of wealth of individuals. One such data set, from the Forbes 400, is analyzed in Ref., and is shown to have the Pareto exponent close to $`1.36`$. We have recently carried out a similar analysis with 1996–1999 data of approximately 400 of the richest people of U.S.A. and 100 of the richest people in Japan and found that the results are consistent with this exponent.

As there is no steady model of personal incomes it is important to approach from many different aspects. Personal income might be related to a kind of popularity such as the citation numbers in scientific publications or access number distribution for Internet homepages. Citation frequency distribution has been examined for publications in the field of physics and a power law tail with exponent about $`3`$ is reported . As for homepage popularity ranking it is known that we have another Zipf’s law for this quantity, namely the power law exponent is about $`1`$ .

Our other finding concerning the distribution of bankrupt companies is considered to be suggestive for the understanding of money flow statistics among companies. The amount of debt of a company is roughly given by accumulated negative incomes that exceeded the whole asset of the company. Therefore, combining our result with the known Zipf’s law for positive incomes we may conjecture that the distributions of money flows among companies both positive and negative generally tend to follow Zipf’s law.

For the purpose of construction of a numerical model of company income statistics, we can summarize the known conditions to be satisfied as follows:

1. The amount of annual summation of outgoing money from each company is of the same order of its assets, also the summation of incoming money is of the same order of the assets .
2. By definition a company’s income is given by the difference of these moneys, that is, incoming money minus outgoing money. The averaged incomes is nearly proportional to a fractional power of asset with exponent estimated to be about 0.85 .
3. The distribution of income follows Zipf’s law for both positive and negative regions.
4. The variance of growth fluctuation of assets is proportional to a fractional power of assets with exponent about $`0.15`$ .
5. The distribution of assets seems to follow a log-normal law , but we need better quality real data to assess this empirical law.

Although no model has successfully satisfied all of these conditions, there are two promising approaches towards this goal. One is a company interaction model based on a zero-sum game  introduced to explain the asset’s variance law 4. In this model companies play a kind of stochastic game in which money flow from losers to winners according to a certain rule. In this model the condition 4 and 5 are satisfied in a rough sense, however, condition 3 is not fulfilled. The condition 1 depends on the definition of the annual year, which can be tuned, while condition 2 has not been checked.

The other approach is a territory occupation model that is based on the assumption that a company’s income may be proportional to the share of customers . Companies under certain stochastic rules in this model share the customers in a new field of a market. From this approach it is easy to satisfy Zipf’s law for positive incomes. However, other conditions are out of the range at present.

Acknowledgments
The authors would like to thank Mr. Nakano of the Japanese Tax Administration for conversation on the properties of the data on the web pages. Numerical computation in this work was in part supported by the computing facility at the Yukawa Institute for Theoretical Physics. The authors would also like to thank Dr. John Constable (Magdalene College, U.K.) for careful reading of the manuscript.


text:
# Chiral non-linear 𝜎-models as models for topological superconductivity

\[

## Abstract

We study the mechanism of topological superconductivity in a hierarchical chain of chiral non-linear $`\sigma `$-models (models of current algebra) in one, two, and three spatial dimensions. The models have roots in the 1D Peierls-Fröhlich model and illustrate how the 1D Fröhlich’s ideal conductivity extends to a genuine superconductivity in dimensions higher than one. The mechanism is based on the fact that a point-like topological soliton carries an electric charge. We discuss a flux quantization mechanism and show that it is essentially a generalization of the persistent current phenomenon, known in quantum wires. We also discuss why the superconducting state is stable in the presence of a weak disorder.

\]

Introduction. In this paper we discuss models of topological superconductivity, where fermionic currents are dynamically coupled with topological textures of chiral bosonic fields . The superconducting state in these models emerges due to an unbroken chiral symmetry. Electronic off-diagonal correlations in topological superconductors are different from those of conventional BCS-type superconductors, but rather close to those of Luttinger liquids . (i) Topological superconductivity takes place in standard models of current algebras (2-3) (also called chiral non-linear $`\sigma `$-models), and (ii) the hydrodynamic action (and wave functions of low energy states) possesses a geometric phase (or $`\vartheta `$-term) – an imaginary part of the Euclidian action

$$𝒮_\mathrm{E}=2\pi \lambda _L^2𝑑x𝑑tj_\mu ^2+i\pi N_f\vartheta .$$
(1)
Here $`j_\mu =(v_0\rho ,\stackrel{}{j})`$ is an electronic current and density deviation from the equilibrium (we set $`v_0=1`$ in the following), $`\lambda _L`$ is the London’s penetration depth, and $`N_f`$ is a degeneracy of an electronic state (e.g., spin).

The first term in (1) gives the London equations. $`\vartheta `$ is an integer valued topological invariant. The theta term of (1) does not affect the equations of motion but the monodromy of the wave function. It is the geometric phase which carries an information about spacetime transformation properties of excitations. In particular, it changes by $`1`$ in the processes of a $`2\pi `$-rotation of a group of $`N_f`$ electrons around some axis or an interchange of positions of two groups of $`N_f`$ electrons. In a multi-connected sample, it also changes by $`1`$ under the translation of $`N_f`$ electrons along a noncontractable path. We show that this term provides a mechanism for flux quantization in a multi-connected geometry. The geometric phase is, perhaps, the most important difference from the s-wave BCS theory, where the Euclidian action is real. It emphasizes a quantum interference aspect of the mechanism. While discussing the $`\vartheta `$-term, we rely on the results of Ref..

Models of current algebras. The models of topological superconductivity in two and three spatial dimensions are the standard non-linear $`\sigma `$-models of current algebras (see e.g. ). These models describe Dirac fermions interacting with a chiral boson field

$`_2`$ $`=`$ $`\overline{\psi }\left(i\widehat{D}m\stackrel{}{n}\stackrel{}{\tau }\right)\psi +W_2[\stackrel{}{n}],`$ (2)
$`_3`$ $`=`$ $`\overline{\psi }\left[i\widehat{D}m(\pi _0+i\gamma _5\stackrel{}{\pi }\stackrel{}{\tau })\right]\psi +W_3[\pi ],`$ (3)
where $`i\widehat{D}=\gamma _\mu (i_\mu +A_\mu )`$ and $`\gamma _\mu `$ are $`2\times 2`$ and $`4\times 4`$ Dirac matrices in spatial dimensions two and three. The fermions are $`SU(2)`$ doublets and $`\stackrel{}{\tau }=(\tau _1,\tau _2,\tau _3)`$ are the Pauli matrices. In addition, a fermion variable $`\psi `$ carries a (suppressed) flavor index running from 1 to $`N_f`$. The parity of the number of flavors $`N_f`$ plays an important role in what follows. The chiral fields take values on two and three-dimensional spheres thus matching the number of spatial dimensions: $`\stackrel{}{n}^2=1`$ and $`\pi _0^2+\stackrel{}{\pi }^2=1`$. We test the system by an external electromagnetic potential $`A_\mu `$. Finally, $`W_{2,\mathrm{\hspace{0.17em}3}}`$ describe a bare action of the chiral field, other than the dynamics induced by fermions. In 2D, e.g., it is $`W_2=\frac{M}{8\pi }(\dot{\stackrel{}{n}}^2s^2(_j\stackrel{}{n})^2)+\mathrm{}`$. The models have an electronic origin and can be viewed as an effective description of certain interacting electronic systems. We do not discuss that here.

We emphasize that the essential properties of these models are inherited from the model of Fröhlich’s ideal conductivity in one spatial dimension . We, therefore, add the latter to the hierarchy. In this model, 1D fermions interact with a chiral field taking values on a circle:

$$_1=\overline{\psi }\left(i\widehat{D}me^{i\gamma _5\varphi }\right)\psi .$$
(4)
It is well known  that $`2\pi `$-phase solitons of the Fröhlich model are charged and move freely through the system, making it an ideal conductor. It is not a superconductor, however, since the restricted spatial geometry does not allow for a transverse motion of a charge flow. As a result, a weak disorder or commensurability effects destroy the ideal conductivity. We will see that in dimensions greater than one (models (2,3)) a phenomenon similar to Fröhlich’s ideal conductivity results in superconductivity even in the presence of weak disorder.

The physics of the models (2-4) is controlled by the chiral symmetry. In 1D it is: $`\varphi \varphi +\alpha ,\psi e^{\frac{i}{2}\gamma _5\alpha }\psi `$. In 2D it is: $`\stackrel{}{n}\stackrel{}{\tau }U^1\stackrel{}{n}\stackrel{}{\tau }U`$ and $`\psi U^1\psi `$, where $`U`$ is an $`SU(2)`$ matrix. In 3D it is $`gU^1gU`$ and $`\psi U^{\gamma _5}\psi `$, where $`g=\pi _0+i\stackrel{}{\pi }\stackrel{}{\tau }`$.

We argue that if the global chiral symmetry is unbroken, the ground states of the 2D and 3D models are superconducting. These models then exhibit a Meissner effect and a flux quantization with period $`\frac{\mathrm{\Phi }_0}{N_f}`$ if $`N_f`$ is even and $`\frac{\mathrm{\Phi }_0}{2N_f}`$ if $`N_f`$ is odd, where $`\mathrm{\Phi }_0=hc/e`$ is the flux quantum. Notice that the period is always a fraction of $`\mathrm{\Phi }_0`$ with an even denominator regardless of whether the number of fermionic species is even or odd. For the most interesting cases of polarized ($`N_f=1`$) and unpolarized ($`N_f=2`$) electrons, the flux is quantized with the period $`\mathrm{\Phi }_0/2`$.

Electronic and topological currents. In every model (2-4) a coordinate space matches the target space of chiral fields. Therefore, there exist spatial configurations of a chiral field which wrap over the target space with an arbitrary integer winding number $`Q`$ – topological charge. They are solitons: phase solitons in 1D and skyrmions in $`D=2,3`$. The topological charge $`Q`$ may be written as a spatial integral $`Q=d^DxJ_0`$ of the zeroth component of topological current $`J_\mu `$:

$`J_\mu `$ $`=`$ $`{\displaystyle \frac{1}{2\pi }}ϵ_{\mu \nu }_\nu \varphi ,`$ (5)
$`J_\mu `$ $`=`$ $`{\displaystyle \frac{1}{8\pi }}ϵ_{\mu \nu \lambda }\stackrel{}{n}_\nu \stackrel{}{n}\times _\lambda \stackrel{}{n},`$ (6)
$`J_\mu `$ $`=`$ $`{\displaystyle \frac{1}{12\pi ^2}}ϵ_{\mu \nu \lambda \rho }ϵ_{abcd}\pi _a_\nu \pi _b_\lambda \pi _c_\rho \pi _d.`$ (7)
In configurations without singularities, a topological current conserves $`_\mu J_\mu =0`$. This reflects the fact that topological charge cannot be changed by smooth variations of chiral fields.

It is known that the energy spectra of the models (2-4) in a smooth static soliton background with a topological charge $`Q`$ have an additional $`N_fQ`$ levels with negative energy, compared to the spectrum in a topologically trivial background. Being filled, these levels create an electric charge localized on the solitons. The motion of the solitons then results in an electric current since the electric charge moves together with a soliton as long as the motion is adiabatic . This is summarized by

$$j_0d^Dx=N_fQ;j_\mu (x)=N_fJ_\mu (x)+\mathrm{},$$
(8)
where $`j_\mu =\overline{\psi }\gamma _\mu \psi `$ is an electronic current in a background chiral field and $`J_\mu `$ is a topological current of a chiral field (5-7). The density of the electronic current (8) has corrections (suppressed by the gap $`m`$) due to terms with higher gradients. However, the corrections do not have a topological nature. The change of total charge in a soliton background is exact.

Spin and momentum of a soliton and the $`\vartheta `$-term. Intuitively, it is appealing that a soliton carrying a fermionic charge $`N_fQ`$ has a spin $`N_fQ/2`$ and statistics that depend on whether $`N_fQ`$ is even or odd. It has been shown in Refs. that this is indeed true – an intrinsic quantum number of a soliton – a fermionic number (8) converts into rotational and translational quantum numbers. Spin and statistics of solitons are described by a geometric phase of the wave function. The geometric phase changes by $`\pi N_fQ`$ under an adiabatic $`2\pi `$-rotation of a soliton or an adiabatic exchange of the positions of two solitons. In a multi-connected geometry, the geometric phase also changes by the same amount under an adiabatic translation of a soliton along a noncontractable cycle. Along such a path, the wave function acquires a phase $`e^{i\pi N_fQ}`$, which indicates the spin and statistics of a soliton. The geometric phase also results in a shift in quantization of the orbital momentum $`L=N_f\frac{Q}{2}+l`$ and momentum $`\mathrm{\Pi }_i=\frac{2\pi }{L_i}(N_f\frac{Q}{2}+n_i)`$ along a noncontractable cycle (here $`l`$ and $`n_i`$ are integers and $`L_i`$ is the length of the cycle).

We write fermionic parts of our models as $`=\overline{\psi }𝒟\psi `$ with $`𝒟`$ an operator, dependent on chiral fields. Then the induced action $`𝒲`$ of the non-linear $`\sigma `$-model is a result of an integration over fermions in a fixed background configuration of a chiral field $`e^𝒲=e^{{\scriptscriptstyle 𝑑x𝑑t}}𝑑\psi 𝑑\overline{\psi }=(\text{Det}𝒟)^{N_f}`$. Then the geometric phase appears as an imaginary part of the Euclidian action of a non-linear $`\sigma `$-model. Indeed, in Refs. it has been shown that the geometric phase in the models (2-4) is equal to $`N_f\pi `$ times a topological invariant $`\vartheta `$ that characterizes a topology of spacetime configurations of the chiral fields $`\varphi ,\stackrel{}{n},(\pi _0,\stackrel{}{\pi })`$. The $`\vartheta `$-term is also referred to as a global chiral anomaly, distinct from the local anomaly given by (8).

The $`\vartheta `$-term depends on the compactification of the Euclidian spacetime. If it is a sphere, the $`\vartheta `$-invariant is a representation of the homotopy groups $`\pi _{D+1}(S^D)`$. It is zero in 1D, any integer (the Hopf number) in 2D, and 0 or 1 in 3D.

In a more physical case, the Euclidian time is a circle and the coordinate space may also have some noncontractable cycles. The topological $`\vartheta `$-term then is corrected by terms describing the motion of solitons along the cycles. The explicit expression for the $`\vartheta `$-term requires a higher profile of the topology , but its physical meaning is simple. E.g., an adiabatic rotation of a soliton by an angle $`\varphi `$ results in the $`\vartheta `$-term $`\vartheta [\varphi ]=Q\dot{\varphi }\frac{dt}{2\pi }`$ (see  and references therein). If a soliton is adiabatically translated along a noncontractable spatial cycle, the $`\vartheta `$-term is $`\vartheta [\varphi ]=Q\dot{u}\frac{dt}{2\pi }`$, where $`u`$ is a displacement of the soliton. The latter is the only term which remains in 1D, where it is known to be responsible for the parity effect in the persistent current phenomenon observed in quantum wires .

Non-linear $`\sigma `$-models. Fermionic number (8) and the $`\vartheta `$-term determine the imaginary part of the Euclidian action of a non-linear $`\sigma `$-model. Its real part is given by a regular gradient expansion. To leading order in $`1/m`$, we have (see e.g.,  and references therein)

$`{\displaystyle \frac{𝒲_1}{N_f}}`$ $`=`$ $`i\pi \vartheta [\varphi ]+{\displaystyle d^2x\left(iA_\mu J_\mu +\frac{1}{8\pi }(_\mu \varphi )^2\right)},`$ (9)
$`{\displaystyle \frac{𝒲_2}{N_f}}`$ $`=`$ $`i\pi \vartheta [\stackrel{}{n}]+{\displaystyle d^3x\left(iA_\mu J_\mu +\frac{m}{8\pi }(_\mu \stackrel{}{n})^2\right)},`$ (10)
$`{\displaystyle \frac{𝒲_3}{N_f}}`$ $`=`$ $`i\pi \vartheta [g]+{\displaystyle d^4x\left(iA_\mu J_\mu +\frac{F_\pi ^2}{2}\mathrm{tr}(_\mu g^1_\mu g)\right)},`$ (11)
where $`F_\pi ^2=\frac{1}{2\pi ^2}m^2\mathrm{ln}\frac{\mathrm{\Lambda }}{m}`$ with $`\mathrm{\Lambda }`$ being an ultraviolet cutoff and $`g=\pi _0+i\stackrel{}{\pi }\stackrel{}{\tau }`$ a unitary matrix. The first two terms represent global and local chiral anomalies.

London’s hydrodynamics. The hydrodynamics of a superconductor is defined by a rigidity with respect to transverse charge currents $`j_\mu (k),j_\nu (k)=(4\pi \lambda _L^2e^2)^1(\delta _{\mu \nu }\frac{k_\mu k_\nu }{k^2})`$ at $`k0`$ and can be summarized by the London Lagrangian: $`=2\pi \lambda _L^2j_\mu ^2`$.

In topological liquids, charge current is equal (at small $`k`$) to the topological current (8), so that a criterion for superconductivity translates to a transverse rigidity of topological currents

$$J_\mu (k),J_\nu (k)=(4\pi N_f^2\lambda _L^2e^2)^1\left(\delta _{\mu \nu }\frac{k_\mu k_\nu }{k^2}\right).$$
(12)
In 1D the first terms of the gradient expansion (9) already give (12). We will argue that the transverse rigidity of topological currents takes place in any dimension if the global chiral symmetry is unbroken, i.e., if the chiral fields are in a quantum disordered state: $`\stackrel{}{n}(x),\stackrel{}{n}(0)0`$ at large $`x`$ in 2D, or $`g(x),g^1(0)0`$ in 3D (it is always true in 1D model).

The situation in higher dimensions is obscured by the fact that the models (2-3) are not renormalizable and that the symmetry of the ground state depends on the effect of small distances not included in the continuum effective model. To implement the state with unbroken chiral symmetry (a quantum disordered state), one might choose a non-electronic part $`W_{2,3}`$ of the models (2-3). The latter does not contribute to the anomalous terms. In addition, we assume that the correlation radius of the quantum disordered state $`r_c`$ is large so that the chiral field is smooth and the electronic density follows the topological density according to eq. (8).

The arguments that a quantum disordered state of a non-linear $`\sigma `$-model, regardless of its dimension, is rigid with respect to topological excitations, are general. Here we adopt them for the 2D model following Ref.. We first write the non-linear $`\sigma `$-model (10) in the $`\mathrm{CP}^1`$ representation by writing $`(_\mu \stackrel{}{n})^2`$ as $`4z^{}(i_\mu a_\mu )^2z`$, where $`z=(z_1,z_2)`$ is a complex doublet with constraint $`z^{}z=1`$ related to $`\stackrel{}{n}`$ as $`\stackrel{}{n}=z^{}\stackrel{}{\tau }z`$. The field strength of the gauge field $`a_\mu `$ is proportional to the topological current of $`\stackrel{}{n}`$ according to $`J_\mu =\frac{1}{2\pi }ϵ^{\mu \nu \lambda }_\nu a_\lambda `$ – it describes the fluctuations of the solid angle formed by the $`\stackrel{}{n}`$-field. The rigidity of topological currents (12) means that the propagator of the gauge field $`a_\mu `$ induced by $`z`$-bosons is massless.

It is obvious that if the chiral field $`\stackrel{}{n}`$ is ordered, than $`z`$-bosons are condensed and $`J_\mu (k),J_\nu (k)\kappa k^2`$, where the stiffness $`\kappa `$ is proportional to the amount of the condensate of z-bosons. This combined with (8) leads to an insulating behavior.

On the contrary, if the chiral symmetry is unbroken the system does not respond to the constant transverse field $`a_\mu `$. This means that the gauge field is massless $`a_\mu (k),a_\nu (k)(r_ck^2)^1`$, which leads to (12) with $`r_c\lambda _L^2e^2`$.

Combining (12), the fermionic number of solitons, and the geometric phase (9-11), we obtain the hydrodynamic description of a topological superconductor (1) in D=2,3 and Fröhlich’s ideal conductor in D=1. In terms of topological currents, it is

$$𝒮=N_f(2\pi N_f\lambda _L^2e^2J_\mu ^2+ieJ_\mu A_\mu )𝑑x𝑑t+i\pi N_f\vartheta .$$
(13)
Disorder. It is generally correct that if the Meissner effect takes place in a system with no disorder, then a superconducting state remains intact in the presence of weak disorder, regardless of the mechanism that led to superconductivity. In 1D there are no transverse directions and no true Meissner effect. As a result the Fröhlich ideal conductivity is destroyed by the backscattering of a single impurity . A $`2k_f`$-harmonic of the impurity potential $`V_{2k_f}`$ adds the term $`|V_{2k_f}|\mathrm{cos}(\varphi \alpha )`$ to the effective Lagrangian (9), where $`\alpha `$ is the phase of the potential. It pins a soft translational mode of the charge density wave at $`\varphi =\alpha `$, makes fluctuations of $`\varphi `$ massive, and suppresses correlations of topological currents at small momenta.

In higher dimensions, the existence of transverse currents eviscerates this effect. In 2D, similar to the 1D case, we add a backscattering term $`\stackrel{}{V}_{2k_f}\stackrel{}{n}`$ to (10), where $`\stackrel{}{V}_{2k_f}`$ is a vector, representing a $`2k_f`$ harmonic (scattering between Dirac points of the spectrum) of an impurity potential. This term has a tendency to align the $`\stackrel{}{n}`$-field along $`\stackrel{}{V}_{2k_f}`$, but does not change the correlation of topological currents (solid angles of $`\stackrel{}{n}`$) as long as the potential is smaller than a scale provided by the correlation radius $`r_c`$. The system remains in the quantum disordered phase, while the superfluid density is reduced by the amount of the order $`\stackrel{}{V}_{2k_f}^2`$. One can see this formally by writing the impurity term in $`\mathrm{CP}^1`$ representation $`\stackrel{}{V}_{2k_f}z^{}\stackrel{}{\tau }z`$. This term shifts the local chemical potentials of the $`z`$ bosons in opposite directions for the bosons with isospin along and opposite to $`\stackrel{}{V}_{2k_f}`$. If $`z`$-bosons are far from the condensation point (ordering of $`\stackrel{}{n}`$), the small change of chemical potential cannot drive the condensation, and fluctuations of an internal gauge field $`a_\mu `$ remain massless. The superconducting state cannot be destroyed by an infinitesimal amount of impurities. The important difference with 1D is that the impurity potential does not couple to a topological density of solitons (i.e. to solid angle of $`\stackrel{}{n}`$-field in this case). Similar arguments are valid for the 3D model.

Flux quantization. The geometric phase $`i\pi N_f\vartheta `$ in the low energy action (1) provides the mechanism of the flux quantization in topological superconductors in a multi-connected geometry. Consider a sample in the shape of a torus, threaded by magnetic flux $`\mathrm{\Phi }`$ through some particular torus cycle, which we label by 1. Due to the gauge invariance all energy levels are $`\mathrm{\Phi }_0`$-periodic functions of $`\mathrm{\Phi }`$. We will show that the energy of a macroscopic sample at a fixed chemical potential is periodic with a smaller period: $`\mathrm{\Phi }_0/N_f`$ if $`N_f`$ is even and $`\mathrm{\Phi }_0/2N_f`$ if $`N_f`$ is odd.

Starting from the low energy action (13), we pass to the Hamiltonian in a sector with a given topological charge $`Q`$. We introduce displacements $`\stackrel{}{u}`$ according to $`J_0=\stackrel{}{}\stackrel{}{u},\stackrel{}{J}=\dot{\stackrel{}{u}}`$. The relevant part of the $`\vartheta `$-term can be written in terms of zero modes (the homogeneous part) of displacements $`\overline{u}_i=u_id^Dx`$. E.g, a soliton moving around the cycle 1 contributes $`\vartheta =Q\frac{\dot{\overline{u}}_1}{L_1}𝑑t`$ (similar to the one in 1D, where $`\overline{u}_1/L_1=\varphi /2\pi `$).

Neglecting all excitations transverse to the torus cycle pierced by the flux, we set all currents except for $`J_1`$ to zero. Then we neglect all harmonics of the displacement $`u_1`$ except for the zero mode. Their contribution to the energy does not depend on the flux. Setting the gauge potential to be a constant along the cycle 1, $`\stackrel{}{A}=(2\pi \mathrm{\Phi }/L_1\mathrm{\Phi }_0,\mathrm{\hspace{0.17em}0},\mathrm{\hspace{0.17em}0})`$, we obtain the flux dependent part of the Lagrangian (we return to the real time)

$$L=\frac{2\pi N_f^2\lambda _L^2e^2}{V}\left(\dot{\overline{u}}_{1}^{}{}_{}{}^{2}Q^2\right)2\pi N_f\left(\frac{Q}{2}+\frac{\mathrm{\Phi }}{\mathrm{\Phi }_0}\right)\frac{\dot{\overline{u}}_1}{L_1}.$$

Varying it with respect to $`\dot{\overline{u}}_1`$ we find canonical momentum $`\mathrm{\Pi }_1=4\pi N_f^2\lambda _L^2e^2\frac{\dot{\overline{u}}_1}{V}\frac{2\pi N_f}{L_1}(\frac{Q}{2}+\frac{\mathrm{\Phi }}{\mathrm{\Phi }_0})`$, where we observe a shift of the momentum by the $`\vartheta `$-term. Representing the canonical momentum of the zero mode as $`\mathrm{\Pi }_1=i\frac{}{\dot{\overline{u}}_1}`$, we obtain the wave function of the zero mode $`\mathrm{exp}(i\frac{2\pi n_1}{L_1}\overline{u}_1)`$, where $`n_1`$ is an integer, and the flux dependent part of the energy of the sliding charge density mode

$$E_{n_1}=\frac{\pi }{2\lambda _L^2e^2}\frac{V}{L_1^2}\left(\frac{n_1}{N_f}+\frac{\mathrm{\Phi }}{\mathrm{\Phi }_0}+\frac{Q}{2}\right)^2+2\pi N_f^2\lambda _L^2e^2\frac{Q^2}{V}.$$
(14)
If $`N_f`$ is even, it is obvious that the energy is a $`\mathrm{\Phi }_0/N_f`$\- periodic function of the flux (c.f., ). However, the energy remains $`\mathrm{\Phi }_0/2`$-periodic even at $`N_f=1`$. In fact, if $`N_f`$ is odd the energy is $`\mathrm{\Phi }_0/2N_f`$-periodic. In macroscopic samples, the energy of vacuum polarization by a soliton – the last term in (14) – is negligible compared to the kinetic energy. If $`N_f`$ is odd, the ground state possesses no solitons $`Q=0`$ as long as $`\mathrm{\Phi }<\mathrm{\Phi }_0/4N_f`$. If, however, $`\mathrm{\Phi }_0/4N_f<\mathrm{\Phi }<\mathrm{\Phi }_0/2N_f`$, the ground state at the fixed chemical potential acquires a topological charge $`Q=1`$ and $`N_f`$ additional particles. The freedom to create a soliton, when the flux exceeds $`\mathrm{\Phi }_0/4N_f`$, halves the naive period $`\mathrm{\Phi }_0/N_f`$ of the kinetic energy.

In the 1D setup (a quantum wire) $`V=L_1`$, $`\lambda _L^2e^2=1/4N_f`$, and (14) gives the persistent current $`j=E/\mathrm{\Phi }`$ for the Peierls-Fröhlich model (4). In this case the vacuum polarization energy caused by a static soliton (the last term in eq.(14)) destroys periodicity for odd $`N_f`$ (c.f., ).

This mechanism is to be compared with the familiar flux quantization in the BCS theory. In that case, a flux-dependent part of the energy is given by zero modes of the phase $`\overline{\phi }=2\pi n\frac{x}{L_1}`$ of the order parameter along a noncontractable cycle: $`E_n=\frac{\pi V}{2\lambda _L^2e^2L^2}\left(\frac{n}{2}+\frac{\mathrm{\Phi }}{\mathrm{\Phi }_0}\right)^2`$. It is $`\mathrm{\Phi }_0/2`$-periodic and has deep minima at fluxes which are multiples of $`\mathrm{\Phi }_0/2`$. Although the topological mechanism leads to a similar result, the state counting and the underlying physical picture is different.

We would like to thank R.B. Laughlin, P.A. Lee, A.M. Polyakov, and X.-G. Wen for discussions. A.A. was supported by NSF DMR 9813764. P.W. was supported by grants NSF DMR 9971332 and MRSEC NSF DMR 9808595.