text:
# Multifractality and critical fluctuations at the Anderson transition

## I Introduction

As is well known, in $`d>2`$ dimensions a disordered electronic system undergoes, with increasing strength of disorder, a transition from the phase of extended states to that of localized states (Anderson transition). Another important realization of the Anderson critical point is the Quantum Hall plateau transition in a 2D system in strong magnetic field. One of the hallmarks of the metal-insulator transition is represented by strong fluctuations of eigenfunctions. These fluctuations can be characterized by a set of inverse participation ratios (IPR)

$$P_q=d^dr|\psi (𝐫)|^{2q}.$$
(1)
The field theory of the Anderson transition is the matrix non-linear $`\sigma `$-model, in the replica  or supersymmetric  formulation. In $`2+ϵ`$-dimensions with $`ϵ1`$ the transition takes place in the weak-coupling regime, allowing for a systematic renormalization-group (RG) treatment, which yields the critical indices in the form of the $`ϵ`$-expansion. In particular, Wegner  found in this way that the IPR show at criticality an anomalous scaling with respect to the system size $`L`$,

$$P_qL^{\tau (q)},\tau (q)=D_q(q1).$$
(2)
Equation (2) should be contrasted with the behavior of the IPR in a good metal (where eigenfunctions are ergodic), $`P_qL^{d(q1)}`$, and, on the other hand, in the insulator (localized eigenfunctions), $`P_qL^0`$.

The scaling (2) characterized by an infinite set of critical exponents $`D_q`$ implies that the critical eigenfunction represents a multifractal distribution . The notion of a multifractal structure was first introduced by Mandelbrot  and was later found relevant in a variety of physical contexts, such as the energy dissipating set in turbulence, strange attractors in chaotic dynamical systems, and the growth probability distribution in diffusion-limited aggregation; see  for a review. More recently, considerable research interest was attracted by the problem of Dirac fermions in a random vector potential, which allows for an exact calculation of the multifractal spectrum . The multifractal exponents play an important role in recent attempts of identification of the conformal theory describing the Quantum Hall plateau transition .

During the last decade, multifractality of critical eigenfunctions at the Anderson transition has been a subject of intensive numerical studies, see Refs.  and references therein. Among all the multifractal dimensions, $`D_2`$ plays the most prominent role, since it determines the spatial dispersion of the diffusion coefficient at the mobility edge .

It should be stressed, however, that Wegner’s result (2) refers to an ensemble-averaged IPR. On the other hand, it is well known that in disordered systems mesoscopic fluctuations from one realization of disorder to another may be very strong. As a result, an average value of some quantity may not provide sufficient information, and one has to speak about the corresponding distribution function. This poses the question of the statistics of the IPR’s $`P_q`$ at criticality, which is a central issue of the present article.

Let us first remind the reader of the existent analytical results concerning the IPR fluctuations. While the direct analytical study of the Anderson transition in 3D is not feasible because of the lack of a small parameter, statistics of energy levels and eigenfunctions in a metallic mesoscopic sample (dimensionless conductance $`g1`$) can be studied systematically in the framework of the supersymmetry method; see  for a review. Within this approach, the IPR fluctuations were studied recently . In particular, the 2D geometry was considered, which, while not being a true Anderson transition point, shows many features of criticality, in view of the exponentially large value of the localization length. It was found that the distribution function of the IPR $`P_q`$ normalized to its average value $`P_q`$ has a scale invariant form. In particular, the relative variance of this distribution (characterizing its relative width) reads

$$\frac{\mathrm{var}(P_q)}{P_q^2}=\frac{Cq^2(q1)^2}{\beta ^2g^2},$$
(3)
where $`C1`$ is a numerical coefficient determined by the sample shape (and the boundary conditions), and $`\beta =1`$ or 2 for the case of unbroken (resp. broken) time reversal symmetry. It is assumed here that the index $`q`$ is not too large, $`q^2\beta \pi g`$. These findings motivated the conjecture  that the IPR distribution at criticality has in general a universal form, i.e. that the distribution function $`𝒫(P_q/P_q^{\mathrm{typ}})`$ is independent of the size $`L`$ in the limit $`L\mathrm{}`$. Here $`P_q^{\mathrm{typ}}`$ is a typical value of the IPR, which can be defined e.g. as a median  of the distribution $`𝒫(P_q)`$. Normalization of $`P_q`$ by its average value $`P_q`$ (rather than by the typical value $`P_q^{\mathrm{typ}}`$) would restrict generality of the statement; see the discussion below. Practically speaking, the conjecture of Ref.  is that the distribution function of the IPR logarithm, $`𝒫(\mathrm{ln}P_q)`$ simply shifts along the $`x`$-axis with changing $`L`$.

The applicability of these results to a generic Anderson transition point has been questioned recently in Ref. . Indeed, the 2D metal represents only an “almost critical” point, and the consideration is restricted to the weak disorder limit $`g1`$ (weak coupling regime in the field-theoretical language), while all the realistic metal-insulator transitions (conventional Anderson transition in 3D, quantum Hall transition etc.) take place in the regime of strong coupling. It was proposed in  (on the basis of numerical simulations for the 3D tight-binding model) that the fractal dimension $`D_2`$ is not a well defined quantity, but rather shows universal fluctuations characterized by some distribution function $`𝒫(D_2)`$ of a width of order unity.

To explore the IPR fluctuations (and also the level statistics, see below) at criticality in the full range from weak to strong coupling, we consider the power-law random banded matrix (PRBM) ensemble. The model is defined  as the ensemble of random Hermitean $`N\times N`$ matrices $`\widehat{H}`$ (real for $`\beta =1`$ or complex for $`\beta =2`$). The matrix elements $`H_{ij}`$ are independently distributed Gaussian variables with zero mean $`H_{ij}=0`$ and the variance

$$|H_{ij}|^2=a^2(|ij|),$$
(4)
where $`a(r)`$ is given by

$$a^2(r)=\frac{1}{1+(r/b)^{2\alpha }}.$$
(5)
At $`\alpha =1`$ the model undergoes an Anderson transition from the localized ($`\alpha >1`$) to the delocalized ($`\alpha <1`$) phase. We concentrate below on the critical value $`\alpha =1`$, when $`a(r)`$ falls down as $`a(r)1/r`$ at $`rb`$. More precisely, we will study the periodic generalization of (5),

$$a^2(r)=\left[1+\frac{1}{b^2}\frac{\mathrm{sin}^2(\pi r/N)}{(\pi /N)^2}\right]^1$$
(6)
(an analog of the periodic boundary conditions), which allows us to get rid of the boundary effects.

In a straightforward interpretation, the PRBM model describes a 1D sample with random long-range hopping, the hopping amplitude decaying as $`1/r^\alpha `$ with the length of the hop. Also, such an ensemble arises as an effective description in a number of physical contexts, such as the quantum Fermi accelerator , the delocalization of impurity-induced quasiparticle states in a 2D $`d`$-wave superconductor , the scattering by a Coulomb center in an integrable billiard , the motion of two interacting particles in a 1D random potential , and the quantum chaos in a billiard with a non-analytic boundary . Very recently, a connection between the level statistics of the PRBM model at $`b1`$ and the correlations in the Luttinger liquid at finite temperature has been established .

At $`\alpha =1`$ the PRBM model was found to be critical for arbitrary value of $`b`$; it shows all the key features of the Anderson critical point, including multifractality of eigenfunctions and non-trivial spectral compressibility . The important property of the ensemble is the existence of the parameter $`b`$ which labels the critical point: Eqs. (4), (6) define a whole family of critical theories parametrized by $`b`$ . This is in full analogy with the family of the conventional Anderson transition critical points parametrized by the spatial dimensionality $`2<d<\mathrm{}`$. The limit $`b1`$ is analogous to $`d=2+ϵ`$ with $`ϵ1`$; it allows for a systematic analytical treatment via the mapping onto a supermatrix $`\sigma `$-model and the weak-coupling expansion . The opposite limit $`b1`$ corresponds to $`d1`$, where the transition takes place in the strong disorder (strong coupling) regime. As we demonstrate below, it is also accessible to an analytical treatment using the RG method introduced by Levitov . Let us also note a similarity with conformal models proposed recently  as candidate theories of the Quantum Hall critical points, which are also parametrized by a continuously changing marginal coupling constant.

In this paper we will combine the analytical study of the eigenfunction and energy level statistics in the two limits $`b1`$ and $`b1`$ with numerical simulations in the full range of $`b`$. The feasibility of the systematic analytical treatment of both regimes, weak-coupling and strong-coupling, as well as of the numerical simulations in a broad range of couplings, makes the PRBM ensemble a unique laboratory for studying general features of the Anderson transition.

As has been already mentioned, we will study not only the eigenfunction fluctuations, but also the energy level statistics. It has been understood  that a scale-invariant level statistics is an intrinsic feature of the Anderson critical point. In particular, the critical level correlations are characterized by a non-trivial value of the spectral compressibility $`0<\chi <1`$ (intermediate between $`\chi =0`$ in a metal and $`\chi =1`$ in an insulator) . While the level correlation function itself is shape-dependent, the value of $`\chi `$ is a fully universal attribute of the critical theory (like critical indices).

Supplementing again the analytical study at $`b1`$ and $`b1`$ by numerical simulations, we are able to calculate the two-level correlation function $`R_2(s)`$ and the spectral compressibility $`\chi `$ in the whole range of $`b`$. Our interest to the critical spectral statistics was additionally motivated by a recent paper , where it was argued that the following exact relation between $`\chi `$ and $`D_2`$ holds:

$$\chi =\frac{dD_2}{2d}.$$
(7)
According to (7), the spectral compressibility should tend to $`1/2`$ in the limit $`D_20`$ (very sparse multifractal), and not to the Poisson value $`\chi =1`$. The derivation of (7) is based, however, on a certain assumption of the decoupling of the energy level and eigenfunction correlations, which is difficult to verify directly. It is important, therefore, to check the validity of the result (7), and such an opportunity is provided by the PRBM model.

The structure of the article is as follows. Sect. II is devoted to the IPR statistics and the multifractal spectrum of the PRBM model. In Sect. III we study the two-level correlation function and the spectral compressibility. Sect. IV summarizes our findings. Some of the results of this work have been presented in a brief form in the Letter .

## II Multifractal spectrum and the IPR statistics

We find it convenient to organize this section in the following way. We begin by formulating in Sect. II A our main results concerning the IPR statistics. Then we present, in Sect. II B and Sect. II C, the analytical calculations in the limits of $`b1`$ and $`b1`$, respectively. The numerical simulations (which have been performed for $`\beta =1`$) fully support the analytical findings, ascertaining that the approximations made in the RG treatment are justified. Also, they allow us to explore the intermediate range of $`b1`$. A discussion of finite-size effects playing an important role in the analysis of the scaling of the IPR distribution, is given in Sec. II D.

### A General considerations and a summary of the results

For further needs we define two sets of fractal exponents characterizing the scaling of the average IPR $`P_q`$ and of the typical IPR $`P_q^{\mathrm{typ}}`$, respectively:

$`P_qN^{\stackrel{~}{\tau }(q)},`$ $`\stackrel{~}{\tau }(q)\stackrel{~}{D}_q(q1);`$ (8)
$`P_q^{\mathrm{typ}}N^{\tau (q)},`$ $`\tau (q)D_q(q1).`$ (9)
Note that we consider $`q>0`$ only; negative $`q`$ are outside the range of applicability of our analytical methods. As has been already mentioned, $`P_q^{\mathrm{typ}}`$ can be defined as a median of the distribution $`𝒫(P_q)`$ ; an alternative definition can be $`P_q^{\mathrm{typ}}=\mathrm{exp}\mathrm{ln}P_q`$. Obviously, an information about the IPR distribution function $`𝒫(P_q)`$ is needed in order to judge whether the exponents $`\stackrel{~}{\tau }(q)`$ and $`\tau (q)`$ are equal to each other or not. As we will demonstrate below, in the limit of large system size $`N`$ the distribution $`𝒫(P_q/P_q^{\mathrm{typ}})`$ becomes independent of $`N`$. An important property of this scale-invariant IPR distribution is its power-law “tail” at large $`P_q`$,

$$𝒫(P_q/P_q^{\mathrm{typ}})(P_q/P_q^{\mathrm{typ}})^{1x_q},P_qP_q^{\mathrm{typ}}.$$
(10)
Of course, the far tail of this distribution becomes increasingly better developed with increasing $`N`$. In other words, the point where the distribution deviates from its limiting scale-invariant form moves to infinity as $`N`$ increases.

It is clear that the relation between $`\tau (q)`$ and $`\stackrel{~}{\tau }(q)`$ depends crucially on whether the power-law exponent $`x_q`$ is smaller or larger than unity. If $`x_q>1`$, the two definitions of the fractal exponents are identical, $`\tau (q)=\stackrel{~}{\tau }(q)`$. This situation will be shown to occur at not too large values of $`q`$; in particular, $`x_2>1`$ at any $`b`$. However, with increasing $`q`$ the value of $`x_q`$ decreases. Once it drops below unity, the average $`P_q`$ starts to be determined by the upper cut-off of the power-law “tail”, which depends on the system size. As a result, $`P_q`$ shows scaling with an exponent $`\stackrel{~}{\tau }_q`$ different from $`\tau _q`$. In this situation the average value $`P_q`$ is not representative and is determined by rare realizations of disorder.

The connection between $`\stackrel{~}{\tau }(q)`$ and $`\tau (q)`$ in the regime $`x_q<1`$ can be elucidated best via introducing the singularity spectrum $`f(\alpha )`$, which is the conventional way of analyzing multifractal distributions . To this end, let us note that the average IPR’s $`P_q`$ are (up to a multiplication by $`N`$) the moments of the distribution $`𝒫(|\psi ^2|)`$ of the eigenfunction intensities. The behavior (8) of the moments corresponds to the intensity distribution function of the form

$$𝒫(|\psi ^2|)\frac{1}{|\psi ^2|}N^{1+f(\frac{\mathrm{ln}|\psi ^2|}{\mathrm{ln}N})}$$
(11)
Indeed, calculating the moments $`|\psi ^{2q}|`$ with the distribution function (11), one finds

$$P_q=N|\psi ^{2q}|𝑑\alpha N^{q\alpha +f(\alpha )},$$
(12)
where we have introduced $`\alpha =\mathrm{ln}|\psi ^2|/\mathrm{ln}N`$. Evaluation of the integral by the saddle-point method reproduces the result (8), with the exponent $`\stackrel{~}{\tau }(q)`$ related to the singularity spectrum $`f(\alpha )`$ via the Legendre transform

$$\stackrel{~}{\tau }(q)=q\alpha f(\alpha );q=f^{}(\alpha ).$$
(13)
It is not difficult to see that the condition $`x_q=1`$ is equivalent to $`f(\alpha )=0`$. Indeed, both conditions $`x_q<1`$ and $`f(\alpha )<0`$ characterize the situation when the average value $`P_q`$ is not representative and is determined by rare realizations of disorder. On a more formal level, this can be derived from the formula relating $`x_q`$ and the fractal exponents, see Eq. (16) below and Sec. II C.

We further denote the value of $`\alpha `$ determined by $`f(\alpha )=0`$ as $`\alpha _{}`$, and the corresponding value of $`q`$ as $`q_c`$ (clearly, both $`\alpha _{}`$ and $`q_c`$ depend on $`b`$). The value of $`\tau _q`$ in the region $`q>q_c`$ can be found by observing that $`P_q^{\mathrm{typ}}`$ can be written in the form similar to (12),

$$P_q^{\mathrm{typ}}_{f(\alpha )0}𝑑\alpha N^{q\alpha +f(\alpha )}.$$
(14)
The restriction on the integration range removes from consideration the rare events of such large values of $`|\psi |^2`$ which can be found only in a small fraction ($`N^{f(\alpha )}`$ with $`f(\alpha )<0`$) of all eigenfunctions. Since for $`q>q_c`$ the saddle-point $`\alpha <\alpha _{}`$ is outside the integration domain, the integral (14) is determined in this case by the boundary $`\alpha _{}`$ of the integration range, yielding (see a related discussion in )

$$\tau (q)=q\alpha _{},q>q_c.$$
(15)
The value of the power-law-tail index $`x_q`$ is related to the fractal exponents as follows:

$$x_q\tau (q)=\stackrel{~}{\tau }(qx_q).$$
(16)
To be precise, we were able to derive Eq. (16) for all $`q`$ in the limit $`b1`$, as well as for integer values of $`x_q=1,2,\mathrm{}`$ at arbitrary $`b`$. We expect, however, that this relation is generally valid.

According to what has been said above, the curve $`q_c(b)`$ separates the regions with the two different types of the multifractal behavior: at $`q<q_c(b)`$ we have $`x_q>1`$ and $`\tau (q)=\stackrel{~}{\tau }(q)`$, while at $`q>q_c(b)`$ the tail index $`x_q<1`$ and $`\tau (q)`$ is different from $`\stackrel{~}{\tau }(q)`$ and given by (15). We have calculated the asymptotic form of the “phase boundary” $`q_c(b)`$ in both limits $`b1`$ and $`b1`$,

$$q_c(b)\{\begin{array}{cc}(2\pi \beta b)^{1/2},\hfill & b1\hfill \\ 2.4056,\hfill & b1\hfill \end{array}$$
(17)
Notice that $`q=2`$ always belongs to the low-$`q`$ phase, i.e. $`\tau (2)=\stackrel{~}{\tau }(2)`$ for all $`b`$. For $`q>q_c(b)`$ we find from Eqs. (15), (16) that $`x_q=q_c(b)/q`$, while in the opposite regime $`q<q_c(b)`$ the value of $`x_q`$ is determined by the form of the function $`\stackrel{~}{\tau }(q)`$. In particular, at $`b1`$ we have

$$\stackrel{~}{\tau }(q)(q1)\stackrel{~}{D}_q(q1)(1q/2\beta \pi b),$$
(18)
yielding $`x_q2\beta \pi b/q^2`$ for $`q<(2\beta \pi b)^{1/2}`$. In the other limit, $`b1`$, the function $`\stackrel{~}{\tau }(q)`$ has a somewhat more complicated form

$$\stackrel{~}{\tau }(q)\frac{4b}{\sqrt{\pi }}\frac{\mathrm{\Gamma }(q1/2)}{\mathrm{\Gamma }(q1)},$$
(19)
and Eq. (16) does not seem to have an analytical solution for $`x_q`$. However, for the particularly important case $`q=2`$ we find $`x_2=3/2`$, while all higher integer $`q=3,4,\mathrm{}`$ are already above the phase boundary $`q_c(b1)2.4`$.

### B Regime $`b1`$

The regime $`b1`$ can be studied via the mapping onto the supermatrix $`\sigma `$-model . The $`\sigma `$-model action has in momentum space the form

$$S[Q]=\beta \mathrm{Str}\left[\frac{1}{t}\frac{dk}{2\pi }|k|Q_kQ_k\frac{i\pi \nu \omega }{4}Q_0\mathrm{\Lambda }\right],$$
(20)
where $`Q_k=_re^{ikr}Q(r)`$ and $`Q(r)`$ is a $`4\times 4`$ ($`\beta =2`$) or $`8\times 8`$ ($`\beta =1`$) supermatrix field constrained by $`Q^2(r)=1`$ and belonging to a coset space with the origin $`\mathrm{\Lambda }=\mathrm{diag}(\mathrm{𝟏},\mathrm{𝟏})`$. Furthermore, $`\nu `$ is the density of states given by the Wigner semicircle law

$$\nu (E)=\frac{1}{2\pi ^2b}(4\pi bE^2)^{1/2},|E|<2\sqrt{\pi b},$$
(21)
and $`t1`$ is the coupling constant,

$$\frac{1}{t}=\frac{\pi }{4}(\pi \nu )^2b^2=\frac{b}{4}\left(1\frac{E^2}{4\pi b}\right).$$
(22)
For a system of finite size $`N`$ with the periodic generalization (6) of the $`1/r`$ decay law of $`a(r)`$ the $`k`$–integration in (20) is replaced by summation in the usual way:

$$\frac{dk}{2\pi }\frac{1}{N}\underset{k=2\pi n/N;n=0,\pm 1,\pm 2,\mathrm{}}{}.$$

The eigenfunction statistics can now be studied via the same methods as for conventional metallic samples. The main difference between the action (6) and that of the diffusive $`\sigma `$-model is in the replacement of the diffusion operator $`\frac{\pi \nu }{8}Dk^2`$ by $`\frac{1}{t}|k|`$. Consequently, all calculations within the weak coupling expansion of the $`\sigma `$-model are generalized to the PRBM case by substituting $`\mathrm{\Pi }(k)=t/8|k|`$ for the diffusion propagator $`\mathrm{\Pi }(k)=1/\pi \nu Dk^2`$. In particular, calculating the average IPR $`P_q`$, one finds the following result for the fractal dimensions $`\stackrel{~}{D}_q`$ :

$$\stackrel{~}{D}_q1q\frac{t}{8\pi \beta },q<\frac{4\pi \beta }{t}.$$
(23)
Fig. 1 shows that this result is in good agreement with numerical simulations.

Similar to the case of 2D  or $`(2+ϵ)`$ dimensions , Eq. (23) describes weak multifractality: the deviation of the fractal exponents from the “normal” value $`d=1`$ is proportional to the small parameter $`t`$. The Legendre transform of $`\stackrel{~}{\tau }(q)=(q1)\stackrel{~}{D}_q`$ yields the singularity spectrum (for definiteness, we concentrate on the band center $`E=0`$, where $`t=4/b`$)

$$f(\alpha )1\frac{(\alpha _0\alpha )^2}{4(\alpha _01)};\alpha _0=1+\frac{1}{2\beta \pi b},$$
(24)
which crosses the $`x`$-axis at the point

$$\alpha _{}\left[1\frac{1}{(2\beta \pi b)^{1/2}}\right]^2,$$
(25)
corresponding to $`q_c(b)=(2\beta \pi b)^{1/2}`$.

In Fig. 2 we confront our analytical findings with data from numerical simulations. At $`b=4`$ the parabola represents the numerical data well up to $`q8`$. The deviations from the asymptotic (parabolic) form are much more pronounced at $`b=1`$. These deviations are a precursor of the crossover to the small-$`b`$ regime (Sec. II C), where the parabolic approximation breaks down completely. The sign of the deviations (downwards) is consistent with the fact that at $`b=1/2\pi `$ the parabolic approximation would predict $`\alpha _{}=0`$, while we expect $`\alpha _{}>0`$ for all $`b`$, in view of the absence of localized states.

Now we turn to the IPR fluctuations, which are found by generalizing the results obtained for metallic samples . In particular, the IPR variance is given for $`qq_c(b)`$ by

$`{\displaystyle \frac{\mathrm{var}(P_q)}{P_q^2}}`$ $`=`$ $`{\displaystyle \frac{2}{\beta ^2}}q^2(q1)^2{\displaystyle \frac{1}{N^2}}{\displaystyle \underset{k}{}}\mathrm{\Pi }^2(k)`$ (26)
$`=`$ $`{\displaystyle \frac{1}{24\beta ^2}}{\displaystyle \frac{q^2(q1)^2}{b^2}},`$ (27)
where the $`k`$-summation goes over the non-zero harmonics $`k=2\pi j/N`$ with $`j=\pm 1,\pm 2,\mathrm{}`$. Equation (27) is the PRBM counterpart of formula (3) for 2D metallic systems. The higher moments of the IPR distribution were studied by Prigodin and Altshuler  (see also ); generalizing these results, we find that the irreducible moments (cumulants) of the order $`2n2\pi \beta b/q^2`$ are given by

$`{\displaystyle \frac{P_q^n}{P_q^n}}`$ $`=`$ $`{\displaystyle \frac{(n1)!}{2}}\left[{\displaystyle \frac{2}{\beta }}q(q1)\right]^n{\displaystyle \frac{1}{N^n}}{\displaystyle \underset{k}{}}\mathrm{\Pi }^n(k)`$ (28)
$`=`$ $`(n1)!\left({\displaystyle \frac{q(q1)}{2\pi \beta b}}\right)^n\zeta (n),`$ (29)
where $`\zeta (n)`$ is the Riemann $`\zeta `$-function. Defining in analogy with

$$\stackrel{~}{P}=\left[\frac{P_q}{P_q}1\right]\frac{2\pi \beta b}{q(q1)},$$
(30)
we have for the cumulants of $`\stackrel{~}{P}`$

$`\stackrel{~}{P}=0;`$ (31)
$`\stackrel{~}{P}^n=(n1)!\zeta (n)K_n,n=2,3,\mathrm{}.`$ (32)
This allows us to restore the corresponding distribution function:

$`𝒫(\stackrel{~}{P})`$ $`=`$ $`{\displaystyle _{\mathrm{}}^{\mathrm{}}}{\displaystyle \frac{ds}{2\pi }}e^{is\stackrel{~}{P}}\mathrm{exp}\left[{\displaystyle \underset{n=2}{\overset{\mathrm{}}{}}}K_n{\displaystyle \frac{(is)^n}{n!}}\right]`$ (33)
$`=`$ $`{\displaystyle _{\mathrm{}}^{\mathrm{}}}{\displaystyle \frac{ds}{2\pi }}e^{is(\stackrel{~}{P}+𝐂)}\mathrm{\Gamma }(1+is)`$ (34)
$`=`$ $`e^{\stackrel{~}{P}𝐂}\mathrm{exp}(e^{\stackrel{~}{P}𝐂}),`$ (35)
where $`𝐂0.5772`$ is the Euler constant. The restriction on $`n`$ given above (29) implies that Eq. (35) is valid for $`P_q/P_q11`$.

The similarity with the 2D metallic regime extends also to the asymptotic behavior of the distribution. Specifically, at $`P_q/P_q11`$ the exponential falloff (35) crosses over to a power-law tail (see  for the discussion of this tail in 2D)

$$𝒫(P_q)(P_q/P_q)^{1x_q}.$$
(36)
To calculate $`x_q`$, we consider the moments

$$P_q^n=\underset{r_1,\mathrm{},r_n}{}|\psi (r_1)|^{2q}\mathrm{}|\psi (r_n)|^{2q}.$$
(37)
For $`n`$ not too large the sum is dominated by the main part of the total configuration space, with all points $`r_i`$ lying far from each other, so that $`P_q^nP_q^nN^{n\tau (q)}`$. In contrast, for sufficiently large $`n`$ the integral is dominated by the contribution from $`r_1r_2\mathrm{}r_n`$, yielding $`P_q^nN^{\stackrel{~}{\tau }(nq)}`$. Therefore, if

$$n\tau (q)=\stackrel{~}{\tau }(nq),$$
(38)
we have the marginal situation, which implies that $`𝒫(P_q)P_q^{1n}`$. This completes the derivation of Eq. (16) in the range $`q<q_c(b)`$. Using now (16) in combination with Eq. (23), we find

$$x_q=\frac{2\pi \beta b}{q^2},q^2<2\pi \beta b.$$
(39)
Note that an analogous consideration in the 2D case yields $`x_q=2\pi \beta g/q^2`$, in full agreement with the result of the optimum-fluctuation method .

We have therefore presented an explicit calculation of the IPR distribution function at $`b1`$ and $`q<q_c(b)`$. The distribution function is scale-invariant and has the form (35), (30) at $`P_q/P_q11`$ and the power-law tail (36), (39) at $`P_q/P_q11`$. Fig. 3 shows results of the numerical simulations for the distribution of the IPR’s $`P_q`$ with $`q=2`$, 4, and 6 at $`b=4`$ (the corresponding value of $`q_c`$ being $`q_c=(8\pi )^{1/2}5`$). It is seen that at $`q=2`$ the analytical formula (35) nicely describes the “main body” of the distribution, with the upward deviations at large $`\stackrel{~}{P}`$ indicating the crossover to the power-law tail (36). The asymptotic behavior (36) is outside the reach of our numerical simulations for $`q=2`$, however, since the condition of its validity $`\stackrel{~}{P}2\pi \beta /q(q1)12.5`$ corresponds to very small values of the distribution function $`𝒫(\stackrel{~}{P})10^5`$, and its clear resolution would require a much larger statistical ensemble. The situation changes, however, with increasing $`q`$ (see the data for $`q=4`$, 6 in Fig. 3). Equation (35) gets inapplicable (since the condition of its validity $`qq_c`$ is not met anymore), and the power-law asymptotic behavior (36) becomes clearly seen. In particular, the inset of Fig. 3 shows the tail for $`q=4`$; the extracted value of the index $`x_41.7`$ is in good agreement with the prediction of the $`b1`$ theory, $`x_4=\pi /2`$.

In conclusion of this subsection, we comment on the notion of the termination of the multifractal spectrum (24), which has been discussed in the literature on disordered Dirac fermions . It is important to realize that, in the present context, there are two types of such termination, depending on whether one studies $`P_q^{\mathrm{typ}}`$ or $`P_q`$. In the former case, the relevant values of $`\alpha `$ are those where $`f(\alpha )0`$, so that the singularity spectrum $`f(\alpha )`$ effectively terminates at $`\alpha _{}`$ (which corresponds to $`q=q_c(b)(2\pi \beta b)^{1/2}`$). For $`q>q_c(b)`$ $`\tau (q)`$ is given by Eq. (15), so that the fractal exponent $`D_q=q\alpha _{}/(q1)`$ saturates as $`D_q\alpha _{}`$ in the limit $`q\mathrm{}`$. In contrast, if the exponent $`\stackrel{~}{\tau }(q)`$ describing the scaling of the average $`P_q`$ is studied, then the behavior (23) continues up to $`q\pi \beta b`$, which corresponds to $`\alpha =0`$. This type of termination (which takes place at $`q`$ parametrically much larger than $`q_c(b)`$) has a physically transparent origin: since $`\alpha =\mathrm{ln}|\psi ^2|/\mathrm{ln}N`$ and in view of the wave function normalization, the allowed values of $`\alpha `$ are restricted by $`\alpha 0`$. More detailed discussion of the behavior of $`\stackrel{~}{\tau }(q)`$ in the vicinity of $`q=\pi \beta b`$ (i.e. the precise form of this termination) is outside the scope of the present article.

### C Regime $`b1`$

In the case $`b1`$ the problem can be studied via the renormalization-group method of Levitov . The idea of the method is as follows. One starts from the diagonal part of the matrix $`\widehat{H}`$, each eigenstate being localized on a single site. Then one includes into consideration non-diagonal matrix elements $`H_{ij}`$ with $`d(i,j)=1`$, where $`d(i,j)`$ is the distance between the sites $`i`$ and $`j`$ with periodic boundary conditions taken into account,

$$d(i,j)=\mathrm{min}\{|ij|,N|ij|\}.$$
(40)
Now one argues that most of these matrix elements are essentially irrelevant, since their typical value is $`b`$, while the energy difference $`|E_iE_j|`$ is typically of order unity. Only with a small probability ($`b`$) is $`|E_iE_j|`$ also of the order of $`b`$, so that the matrix element mixes strongly the two states. Following Levitov, we will say that these two states are in resonance. In this case one is led to consider a two-level problem

$$\widehat{H}_{\mathrm{two}\mathrm{level}}=\left(\begin{array}{cc}E_1& V\\ V& E_2\end{array}\right).$$
(41)
The corresponding eigenfunctions and eigenenergies are

$$\psi ^{(+)}=\left(\begin{array}{c}\mathrm{cos}\theta \\ \mathrm{sin}\theta \end{array}\right);\psi ^{()}=\left(\begin{array}{c}\mathrm{sin}\theta \\ \mathrm{cos}\theta \end{array}\right);$$
(42)
$$E_\pm =\frac{E_1+E_2}{2}\pm |V|\sqrt{1+\tau ^2},$$
(43)
where

$`\mathrm{tan}\theta =\tau +\sqrt{1+\tau ^2},`$ (44)
$`\tau ={\displaystyle \frac{\omega }{2V}},\omega =E_1E_2.`$ (45)
In the next RG step the matrix elements $`H_{ij}`$ with $`d(i,j)=2`$ are taken into account, then those with $`d(i,j)=3`$, and so forth until $`d(i,j)=N/2`$. Each time a resonance is encountered, the Hamiltonian is re-expressed in terms of the new states. Since the probability of a resonance at a distance $`r=d(i,j)`$ is $`b/r`$, the typical scale $`r_2`$ at which a resonance state formed at a scale $`r_1`$ will be again in resonance satisfies

$$\mathrm{ln}\frac{r_2}{r_1}\frac{1}{b},$$
(46)
so that $`r_2`$ is much larger than $`r_1`$. Therefore, when considering the resonant two-level system at the scale $`r_2`$, one can treat the $`r_1`$-resonance state as point-like. Furthermore, it is easy to see that the Gaussian statistics of the matrix element coupling the states on the scale $`r_2`$ is not affected by the transformation to the new basis induced by the $`r_1`$-resonance.

Now we consider the evolution of the IPR distribution with the distance $`r`$; we will denote the corresponding distribution function as $`f(P_q,r)`$. When a resonance occurs, two states with IPR’s $`P_q^{(1)}`$ and $`P_q^{(2)}`$ are replaced by two new states with the IPR’s

$`P_q^{(+)}=P_q^{(1)}\mathrm{cos}^{2q}\theta +P_q^{(2)}\mathrm{sin}^{2q}\theta ,`$ (47)
$`P_q^{()}=P_q^{(1)}\mathrm{sin}^{2q}\theta +P_q^{(2)}\mathrm{cos}^{2q}\theta .`$ (48)
We thus have for real matrices ($`\beta =1`$)

$`{\displaystyle \frac{}{r}}f(P_q,r)=2\nu {\displaystyle _{\mathrm{}}^{\mathrm{}}}𝑑\omega {\displaystyle _{\mathrm{}}^{\mathrm{}}}𝑑V{\displaystyle \frac{1}{\sqrt{2\pi }}}{\displaystyle \frac{\stackrel{~}{r}}{b}}e^{V^2\stackrel{~}{r}^2/2b^2}`$ (49)
$`\times [f(P_q,r)+{\displaystyle }dP_q^{(1)}dP_q^{(2)}f(P_q^{(1)},r)f(P_q^{(2)},r)`$ (50)
$`\times \delta (P_qP_q^{(1)}\mathrm{cos}^{2q}\theta P_q^{(2)}\mathrm{sin}^{2q}\theta )],`$ (51)
where $`\nu =\frac{1}{\sqrt{2\pi }}e^{E^2/2}`$ is the density of states and

$$\stackrel{~}{r}=\frac{N}{\pi }\mathrm{sin}\frac{\pi r}{N}.$$
(52)
The difference between $`\stackrel{~}{r}`$ and $`r`$ is irrelevant for the present calculation, since the $`r`$-integral will be of logarithmic nature and thus dominated by $`rN`$. However, this difference is important for the calculation of the level correlation function below (Sect. III). Transforming the integration measure according to

$$d\omega =2Vd\tau ,d\tau =\frac{1}{2\mathrm{sin}^2\theta \mathrm{cos}^2\theta }d\theta ,$$
(53)
calculating the $`V`$-integral, and specializing on the band center ($`E=0`$) for definiteness, we reduce the evolution equation (49) to the form

$`{\displaystyle \frac{}{\mathrm{ln}r}}f(P_q,r)={\displaystyle \frac{2b}{\pi }}{\displaystyle _0^{\pi /2}}{\displaystyle \frac{d\theta }{\mathrm{sin}^2\theta \mathrm{cos}^2\theta }}`$ (54)
$`\times [f(P_q,r)+{\displaystyle }dP_q^{(1)}dP_q^{(2)}f(P_q^{(1)},r)f(P_q^{(2)},r)`$ (55)
$`\times \delta (P_qP_q^{(1)}\mathrm{cos}^{2q}\theta P_q^{(2)}\mathrm{sin}^{2q}\theta )].`$ (56)
Eq. (54) is a kind of kinetic equation (in the fictitious time $`t=b\mathrm{ln}r`$), with the two terms in the square brackets describing the scattering-out and scattering-in processes, respectively.

Figure 4 shows the results of the numerical integration of Eq. (54) for $`q=2`$ with the initial condition $`f(P_2)=\delta (P_21)`$ at $`t=0`$. It is seen that at sufficiently large $`t`$ the distribution of $`\mathrm{ln}P_2`$ acquires a limiting form, shifting with $`t`$ without changing its shape. This conclusion of scale-invariance of the IPR distribution will be corroborated below by analytical arguments.

Turning to the theoretical analysis, we consider first the average value $`P_q`$. Multiplying Eq. (54) by $`P_q`$ and then integrating over $`P_q`$, we get the evolution equation for $`P_q`$

$$\frac{P_q}{\mathrm{ln}r}=2b\stackrel{~}{T}(q)P_q$$
(57)
with

$`\stackrel{~}{T}(q)`$ $`=`$ $`{\displaystyle \frac{1}{\pi }}{\displaystyle _0^{\pi /2}}{\displaystyle \frac{d\theta }{\mathrm{sin}^2\theta \mathrm{cos}^2\theta }}(1\mathrm{cos}^{2q}\theta \mathrm{sin}^{2q}\theta )`$ (58)
$`=`$ $`{\displaystyle \frac{2}{\sqrt{\pi }}}{\displaystyle \frac{\mathrm{\Gamma }(q1/2)}{\mathrm{\Gamma }(q1)}}`$ (59)
$`=`$ $`{\displaystyle \frac{1}{2^{2q3}}}{\displaystyle \frac{\mathrm{\Gamma }(2q1)}{\mathrm{\Gamma }(q)\mathrm{\Gamma }(q1)}}.`$ (60)
We assumed here that $`q>1/2`$, which is the condition of the existence of the integral in Eq. (58). For smaller $`q`$ the resonance approximation breaks down. Integrating (57) from $`r=1`$ to $`rN`$, we find the multifractal behavior $`P_qN^{\stackrel{~}{\tau }(q)}`$ with the exponents

$$\stackrel{~}{\tau }(q)=2b\stackrel{~}{T}(q).$$
(61)
The function $`\stackrel{~}{T}(q)`$ is shown in Fig. 5. Its asymptotics are

$`\stackrel{~}{T}(q){\displaystyle \frac{1}{\pi (q1/2)}},`$ $`q1/2;`$ (62)
$`\stackrel{~}{T}(q){\displaystyle \frac{2}{\sqrt{\pi }}}q^{1/2},`$ $`q1.`$ (63)
Legendre transformation of (61) produces the $`f(\alpha )`$-spectrum of the form

$$f(\alpha )=2bF(A);A=\alpha /2b,$$
(64)
where $`F(A)`$ is the Legendre transform of $`\stackrel{~}{T}(q)`$. The function $`F(A)`$ is shown in the inset of Fig. 5, its asymptotics are

$`F(A){\displaystyle \frac{1}{\pi A}},`$ $`A0;`$ (65)
$`F(A){\displaystyle \frac{A}{2}},`$ $`A\mathrm{}.`$ (66)
Furthermore, it changes sign at $`A_{}0.5104`$, corresponding to $`q_c2.4056`$.

These analytical findings are fully supported by numerical simulations as can be seen from Fig. 6.

We return now to the IPR distribution function. The scale invariance of the limiting distribution has been already demonstrated via the numerical solution of Eq. (54), see Fig. 4. To show this also analytically, we make the Ansatz

$$f(P_q,r)=r^{\tau (q)}f_0(P_qr^{\tau (q)}).$$
(67)
Substituting (67) into (54), we get the equation

$`\tau (q)[f_0(\stackrel{~}{P}_q)+\stackrel{~}{P}_qf_0^{}(\stackrel{~}{P}_q)]={\displaystyle \frac{2b}{\pi }}{\displaystyle _0^{\pi /2}}{\displaystyle \frac{d\theta }{\mathrm{sin}^2\theta \mathrm{cos}^2\theta }}`$ (68)
$`\times [f_0(\stackrel{~}{P}_q)+{\displaystyle }d\stackrel{~}{P}_q^{(1)}d\stackrel{~}{P}_q^{(2)}f_0(\stackrel{~}{P}_q^{(1)})f_0(\stackrel{~}{P}_q^{(2)})`$ (69)
$`\times \delta (\stackrel{~}{P}_q\stackrel{~}{P}_q^{(1)}\mathrm{cos}^{2q}\theta \stackrel{~}{P}_q^{(2)}\mathrm{sin}^{2q}\theta )].`$ (70)
The fact that the scale $`r`$ has dropped out from Eq. (68) implies the consistency of the Ansatz (67) for the fixed-point distribution. To demonstrate the significance of this statement, we make a more general Ansatz for the limiting distribution

$$f(P_q,r)=\frac{1}{P_q\sigma (r)}g_0\left(\frac{\mathrm{ln}(P_qr^{\tau (q)})}{\sigma (r)}\right),$$
(71)
which allows for a change of the width of the distribution of $`\mathrm{ln}P_q`$ with $`r`$, characterized by a function $`\sigma (r)`$ (note that $`\sigma (r)`$ is defined up to a constant factor, which can be absorbed into the definition of the function $`g_0`$). At $`\sigma (r)=1`$ this reduces to our earlier Ansatz (53), while at $`\sigma (r)=\mathrm{ln}r`$ we get the form proposed in Ref. . Substituting (71) in (54), we find that the $`r`$-dependence drops out of the resulting equation for $`g_0`$ if and only if $`\sigma (r)=\mathrm{const}`$. This means inconsistency of the Ansatz (71) with a non-constant $`\sigma (r)`$ and, in particular, excludes the possibility of a universal distribution of fractal exponents \[$`\sigma (r)=\mathrm{ln}r`$\] advocated in Ref. .

We turn now to the power-law tail of this scale-invariant distribution, $`f_0(\stackrel{~}{P}_q)\stackrel{~}{P}^{x_q1}`$. In order to calculate the index $`x_q`$, we consider Eq. (68) in the limit $`\stackrel{~}{P}_q1`$. It is easy to see that the integral $`𝑑\stackrel{~}{P}_q^{(1)}𝑑\stackrel{~}{P}_q^{(2)}\mathrm{}`$ in the r.h.s. of (68) is dominated by the region $`\stackrel{~}{P}_q^{(1)}\stackrel{~}{P}_q`$, $`\stackrel{~}{P}_q^{(2)}1`$ (or vice versa), the contribution of the region $`\stackrel{~}{P}_q^{(1)}\stackrel{~}{P}_q^{(2)}\stackrel{~}{P}_q`$ being suppressed by an additional factor of $`\stackrel{~}{P}_q^{x_q}`$. Furthermore, when $`\stackrel{~}{P}_q^{(1)}\stackrel{~}{P}_q`$ and $`\stackrel{~}{P}_q^{(2)}1`$, we can neglect $`\stackrel{~}{P}_q^{(2)}`$ in the argument of the $`\delta `$-function. The integrals over $`\stackrel{~}{P}_q^{(1)}`$ and $`\stackrel{~}{P}_q^{(2)}`$ become then trivial, and Eq. (68) reduces to

$$\tau (q)x_q=\frac{2b}{\pi }_0^{\pi /2}\frac{d\theta }{\mathrm{sin}^2\theta \mathrm{cos}^2\theta }(1\mathrm{sin}^{2qx_q}\theta \mathrm{cos}^{2qx_q}\theta ).$$
(72)
Comparing this with (58), (61), we see that the r.h.s. of (72) is nothing but $`\stackrel{~}{\tau }(qx_q)`$, so that Eq. (72) can be rewritten in the form (16).

We analyze now Eq. (16) in the regimes $`q<q_c`$ and $`q>q_c`$. In the case $`q<q_c`$ we expect $`x_q>1`$ and $`\tau (q)=\stackrel{~}{\tau }(q)`$. The latter statement can be directly proven by applying the operation $`𝑑\stackrel{~}{P}_q\stackrel{~}{P}_q\mathrm{}`$ to Eq. (68). \[Clearly, this proof breaks down for $`q>q_c`$ because of the divergence of the integral $`𝑑\stackrel{~}{P}_q\stackrel{~}{P}_qf_0(\stackrel{~}{P}_q)`$.\] Graphical interpretation of Eq. (16) for $`q<q_c`$ is shown in Fig. 7; its solution $`x_q>1`$ decreases with increasing $`q`$, reaching unity at $`q=q_c`$, as expected. For the most frequently studied case $`q=2`$ (“conventional” IPR) we find $`x_2=3/2`$. As to the $`q>q_c`$ regime, we have then $`\tau (q)=q\alpha _{}`$, and the solution of (16) has a very simple form

$$x_q=\frac{q_c}{q},q>q_c.$$
(73)
Let us remind the reader that up to now we considered in this subsection the ensemble of real matrices ($`\beta =1`$). However, all the above considerations are also applicable to the case $`\beta =2`$, with a minor modification. Specifically, the measure of the $`V`$-integration in Eq. (49) should be modified:

$`{\displaystyle _{\mathrm{}}^{\mathrm{}}}𝑑V{\displaystyle \frac{\stackrel{~}{r}}{b\sqrt{2\pi }}}e^{V^2\stackrel{~}{r}^2/2b^2}\mathrm{}`$ (74)
$`{\displaystyle _0^{\mathrm{}}}𝑑VV{\displaystyle \frac{2\stackrel{~}{r}^2}{b^2}}e^{V^2\stackrel{~}{r}^2/b^2}\mathrm{}.`$ (75)
This leads, after the $`V`$-integration, to the replacement of $`b`$ by $`\frac{\pi }{2\sqrt{2}}b`$. With this substitution, all results of this subsection remain valid for $`\beta =2`$.

### D Finite-size effects in the scaling of the IPR distribution

Since in reality one always has to deal with systems of a finite size, the understanding of finite-size effects is important for an accurate interpretation of numerical data. In Figs. 8 we show the evolution of the distribution $`𝒫(\mathrm{ln}P_2)`$ with $`N`$ for three values of $`b`$, representative of the small-$`b`$, the large-$`b`$, and the crossover regimes.

In the case $`b1`$, the evolution of the IPR distribution with $`N`$ is governed by the “time” $`t=b\mathrm{ln}N`$, so that $`t1`$ is the condition of the proximity to the fixed point. Therefore, at small $`b`$ one needs exponentially large values of $`N`$ in order to reach the limiting distribution. (Note that this is not true for the average $`P_q`$, the evolution of which is governed by Eq. (57), implying a much weaker condition $`N1`$ for the scaling regime.) The logarithmically slow approach to the limiting distribution is clearly seen in Figs. 8c, 9.

At $`b1`$ the convergence to the fixed-point distribution is much faster (Fig. 8b). In this regime the condition for the scaling of the IPR distribution is $`N^{1/2}/b1`$, as can be seen by comparing the relative magnitude of the IPR fluctuations at the critical point, $`[\mathrm{var}(P_2)]^{1/2}/P_20.41/b`$ \[see Eq. (27)\] with that in the Gaussian orthogonal ensemble, $`[\mathrm{var}(P_2)]^{1/2}/P_21.64/N^{1/2}`$. Extrapolating the numerically found values of the relative variance to $`1/N=0`$, we find good agreement with the theoretical prediction, as shown in Fig. 10. Note a qualitative difference in the approach to the fixed-point distribution at small and large $`b`$: while at $`b1`$ the height of the distribution $`𝒫(\mathrm{ln}P_2)`$ decreases with $`N`$, approaching the limiting value from above, the behavior is opposite at $`b1`$.

Extrapolating the small-$`b`$ and large-$`b`$ results to the crossover range $`b1`$, we find simply the condition $`N1`$ for the proximity to the fixed-point distribution. Therefore, in the crossover regime $`b1`$ the limiting distribution can be reached most easily. This is precisely what we observe in our numerical simulations. Figure 8a, representing the evolution of the IPR distribution at $`b=1`$ demonstrates the almost perfect scale-invariance of the distribution with only very small deviations (less than $`3\%`$) over a range of system sizes from $`N=256\mathrm{}4096`$.

Let us also comment on other types of ”finite size effects” that appear in numerical simulations. Numerically, it is impossible to perform an average at precisely a given value of the energy. Instead, one averages over an energy interval $`\mathrm{\Delta }E`$ that one would like to choose as big as possible in order to improve statistics. On the other hand, it is clear that $`\mathrm{\Delta }E`$ should not be too big in order to avoid mixing of different critical theories. In our simulations we have chosen $`\mathrm{\Delta }E`$ to be about $`10\%`$ of the bandwidth. This value is still small enough, the corresponding variation of the density of states $`\nu (E)`$ being of the order of $`1\%`$.

Furthermore, the size $`S`$ of the matrix ensemble the average was taken over is an important parameter in the simulations. Typical values we have used are: $`N=256,S=30000;N=512,S=10000;N=1024,S=5000;N=2048,S=1000;N=4096,S=100`$. In some cases, like for the two-level correlation function $`R_2(s)`$ or the full IPR distribution function, the ensemble has to be much larger. In these cases we give the numbers explicitly in the figure caption.

## III Level statistics

The two-level correlation function is defined in the usual way,

$$R_2(\omega )=\frac{1}{\nu (E)^2}\nu (E+\omega /2)\nu (E\omega /2),$$
(76)
where $`\nu (E)=N^1\mathrm{Tr}\delta (E\widehat{H})`$ is the fluctuating density of states. At the critical point $`R_2(\omega )`$ acquires a scale-invariant form (if considered as a function of $`s=\omega /\mathrm{\Delta }`$, the frequency normalized to the mean level spacing $`\mathrm{\Delta }=1/N\nu (E)`$) . The distinct feature of the critical level statistics is a non-trivial value of the spectral compressibility $`0<\chi <1`$ characterizing the linear behavior of the variance of the number $`n()`$ of levels in an energy window $``$ ,

$$\mathrm{var}[n()]=\chi n(),n()\frac{}{\mathrm{\Delta }}1.$$
(77)
The compressibility $`\chi `$ can be expressed through the connected part $`R_2^{(c)}(s)=R_2(s)1`$ of the critical level correlation function as follows:

$$\chi =_{\mathrm{}}^{\mathrm{}}𝑑sR_2^{(c)}(s).$$
(78)
Recently, it was argued in Ref.  that Eq. (7) constitutes an exact relation between the spectral compressibility $`\chi `$ and the fractal dimension $`D_2`$. The derivation of (7) in  is based on Dyson’s idea of Brownian motion through the ensemble of Hamiltonians combined with some assumption of the decoupling of the energy-level and wave function correlations previously proposed in . While this decoupling has been proven to work up to three-loop order in the $`1/g`$-expansion in 2D , its applicability in the strong-coupling regime remained in the status of a conjecture. The PRBM model allows us to check the validity of the relation (7). Similarly to the IPR distribution function, the level correlation function can be calculated analytically in the two limits $`b1`$ and $`b1`$ and numerically in the full range of $`b`$.

In the $`b1`$ regime the two-level correlation function is obtained by an appropriate generalization of the earlier findings for the diffusive samples ; the results can be found in . In particular, considering for simplicity the $`\beta =2`$ ensemble at the band center, the level correlation function has the form

$$R_2^{(c)}(s)=\delta (s)\frac{\mathrm{sin}^2(\pi s)}{(\pi s)^2}\frac{(\pi s/4b)^2}{\mathrm{sinh}^2(\pi s/4b)}.$$
(79)
The correlation function (79) follows the RMT result $`R_2^{(c)}(s)=\delta (s)\mathrm{sin}^2(\pi s)/(\pi s)^2`$ up to the scale $`sb`$ (playing the role of the Thouless energy here), and then begins to decay exponentially. The spectral compressibility at $`b1`$ is given by

$$\chi \frac{1}{2\pi \beta b},b1.$$
(80)
Comparing this with (23), one finds  that the formula (7) is indeed satisfied to leading order in $`1/b`$.

We now turn to the opposite limit $`b1`$. The evolution equation for $`R_2(\omega ,r)`$ can be written down in analogy with Eq. (49):

$`{\displaystyle \frac{R_2(\omega ,r)}{r}}`$ $`=`$ $`{\displaystyle \frac{2}{N}}{\displaystyle _{\mathrm{}}^{\mathrm{}}}𝑑\tau {\displaystyle _{\mathrm{}}^{\mathrm{}}}𝑑V{\displaystyle \frac{1}{\sqrt{2\pi }}}{\displaystyle \frac{\stackrel{~}{r}}{b}}e^{V^2\stackrel{~}{r}^2/2b^2}2|V|`$ (81)
$`\times `$ $`[\delta (\omega 2V\sqrt{1+\tau ^2})\delta (\omega 2V\tau )].`$ (82)
Equation (81) should be integrated over $`r`$ from $`r=0`$ to $`N/2`$ with the boundary condition $`R_2(\omega ,0)=1`$; the result $`R_2(\omega ,N/2)`$ will then give the sought level correlation function. Evaluating the $`V`$-integral in (81) and changing the variables to $`z=2r/N`$ and $`x=2V/\omega `$, we get at $`s=\omega /\mathrm{\Delta }>0`$

$`R_2^{(c)}(s)`$ $`=`$ $`1+{\displaystyle _0^1}𝑑z{\displaystyle \frac{s}{\pi b}}\mathrm{sin}{\displaystyle \frac{\pi z}{2}}{\displaystyle _0^1}{\displaystyle \frac{dx}{\sqrt{1x^2}}}`$ (83)
$`\times `$ $`\mathrm{exp}\left({\displaystyle \frac{s^2x^2}{4\pi b^2}}\mathrm{sin}^2{\displaystyle \frac{\pi z}{2}}\right).`$ (84)
After some algebra, we find the level correlation function to be given by

$$R_2^{(c)}(s)=\delta (s)\mathrm{erfc}\left(\frac{|s|}{2\sqrt{\pi }b}\right),$$
(85)
where $`\mathrm{erfc}(x)=(2/\sqrt{\pi })_x^{\mathrm{}}\mathrm{exp}(t^2)𝑑t`$ is the error function and we have included the $`\delta (s)`$ contribution due to the self-correlation of the energy levels. Substitution of (85) in (78) yields the spectral compressibility

$$\chi 14b,b1.$$
(86)
We see, therefore, that in the limit of small $`b`$ the level repulsion is efficient in a narrow region $`|s|b`$ only, and the spectral compressibility tends to the Poisson value $`\chi =1`$. Hence, formula (7), which would predict $`\chi 1/2`$ at $`b0`$, is violated. Similar violation of (7) is indicated by numerical data for the tight-binding model in dimensions $`d4`$ . Most likely it is never an exact relation, but rather an approximation valid in the weak-multifractality limit only. Strictly speaking, our results do not rule out the possibility that Eq. (7) is exact at $`b`$ exceeding a certain value $`b_c`$ and breaks down at $`b<b_c`$. However, we do not see a physical reason for such a qualitative change being induced by the variation of $`b`$.

These results are fully supported by our numerical data. In particular, Fig. 11 represents the level correlation function $`R_2(s)`$ at $`b=0.1`$ showing a nice agreement with Eq. (85). Note that the finite size effects in the level correlation function at $`b1`$ are much weaker than in the IPR distribution (Sec. II D). Indeed, the only assumption (besides $`b1`$) used in the derivation of the evolution equation (81) is $`N1`$, and then Eq. (85) is obtained without any further approximations. Therefore, in contrast to the IPR statistics, which reaches its fixed-point form at exponentially large $`N`$ (the condition being $`b\mathrm{ln}N1`$), the level statistics acquires the $`N`$-invariant form already at $`N1`$, see Figs. 1112. At $`b1`$ the fixed-point condition is $`Nb^2`$ (the same as for the IPR distribution, see Sec. II D).

To find numerically the spectral compressibility $`\chi `$, we plot the level number variance $`\mathrm{var}[n()]`$ versus the average $`n()`$ (Fig. 12). The data show an extended plateau region in $`\mathrm{var}[n()]/n()`$, determining $`\chi `$. The upper bound for this region is set by the matrix size $`N`$, while the lower bound is $`b`$ (the value of the upper limit at which the integral (78) saturates). We see that the data traces are independent of the system size $`N`$ (with exception of the large-$`n`$ cutoff determined by $`N`$) within the statistical errors. The numerically obtained spectral compressibility in the broad range of $`b`$ is shown in Fig. 13; in the large-$`b`$ and small-$`b`$ regions it agrees well with the corresponding analytical asymptotics.

The above calculation is easily modified to the $`\beta =2`$ case, by replacing the measure of the $`V`$-integration in (81) according to (74). Performing the $`V`$-integral, we now get

$`R_2^{(c)}(s)`$ $`=`$ $`1+{\displaystyle _0^1}𝑑z{\displaystyle \frac{s^2}{\pi b^2}}\mathrm{sin}^2{\displaystyle \frac{\pi z}{2}}{\displaystyle _0^1}𝑑u`$ (87)
$`\times `$ $`\mathrm{exp}\left[(1u^2){\displaystyle \frac{s^2}{2\pi b^2}}\mathrm{sin}^2{\displaystyle \frac{\pi z}{2}}\right],`$ (88)
which yields the result

$$R_2^{(c)}(s)=\delta (s)\mathrm{exp}\left(\frac{s^2}{2\pi b^2}\right).$$
(89)
The spectral compressibility is thus equal to

$$\chi 1\pi \sqrt{2}b.$$
(90)
Again we see that the region of level repulsion shrinks in the limit $`b0`$, with the compressibility approaching the Poisson value $`\chi =1`$.

It is worth mentioning that the above results for the case $`\beta =2`$ can also be obtained by exploiting the relation  between the PRBM model and a random matrix ensemble introduced by Moshe, Neuberger, and Shapiro . This mapping, which becomes exact in both limits $`b1`$ and $`b1`$, relates the level statistics (but not the eigenfunction statistics!) of the two ensembles. On the other hand, the level correlation function of the ensemble of Moshe et al. can be calculated exactly in the case $`\beta =2`$ and is in fact identical to the density correlation function of a 1D non-interacting Fermi gas at a finite temperature . Applying the results of , one obtains for the $`\beta =2`$ PRBM ensemble precisely the results (79), (80) at $`b1`$  and (89), (90) at $`b1`$ .

## IV Conclusions

In this paper, we have presented a detailed study of the statistics of eigenfunctions and energy levels in the family of the critical PRBM models. We have obtained analytical results for the IPR distribution function, the multifractal spectrum and the level correlation function in the two limits of weak and strong multifractality ($`b1`$ and $`b1`$). The analytical results are fully supported by numerical simulations, which also have allowed us to explore the crossover region ($`b1`$).

On the qualitative level, our main findings can be summarized as follows:

1. The distribution function of the IPR (normalized to its typical value $`P_q^{\mathrm{typ}}`$) is scale-invariant in the limit of large system size $`N`$. In other words, the distribution function of the IPR logarithm, $`𝒫(\mathrm{ln}P_q)`$ shifts along the $`x`$-axis with increasing $`N`$, without changing its form and width.
2. The scaling of $`P_q^{\mathrm{typ}}`$ with the system size defines the fractal exponent $`D_q`$, which is a non-fluctuating quantity.
3. The scale-invariant distribution $`𝒫(zP_q/P_q^{\mathrm{typ}})`$ has a power-law tail $`z^{1x_q}`$. At sufficiently large $`q`$ one finds $`x_q<1`$, and the average value $`P_q`$ becomes non-representative and scales with a different exponent $`\stackrel{~}{D}_qD_q`$.
4. The critical spectral statistics shows a crossover from a “quasi-metallic” (close-to-RMT) behavior at $`b1`$ to a “quasi-insulating” (close-to-Poisson) one at $`b1`$. In particular, the spectral compressibility changes from 0 to 1, thus violating the relation (7) in the strong-multifractality regime.

Finally, it is worthwhile to comment on the extent of universality of the IPR distribution. Like the conductance distribution or the level statistics , the IPR distribution at criticality does depend on the system geometry (i.e., on the shape and on the boundary conditions). However, for a given geometry it is independent of the system size and of microscopic details of the model, and is an attribute of the relevant critical theory.

Discussions with Y.V. Fyodorov, V.E. Kravtsov, L.S. Levitov, D.G. Polyakov, I. Varga and I.Kh. Zharekeshev are gratefully acknowledged. This work was supported by the SFB 195 der Deutschen Forschungsgemeinschaft.


format: text

text:
# Kilohertz quasi-periodic oscillations in low mass X–ray binary sources and their relation with the neutron star magnetic field

## 1. Introduction

Color-color diagrams (CDs) have been widely used to gain insight on low mass X–ray binaries (LMXRBs) hosting a neutron star (e.g. Hasinger & van der Klis 1989). These sources are usually classified as atoll or Z sources depending on the path that they trace in the CDs as their luminosity vary. Atoll sources are characterised by an upward curved branch (lower and upper banana state) and one or more harder island states (cf. Fig. 1). Z sources produce a Z-shaped path in the CD with three branches, from hard to soft: horizontal branch (HB), normal branch (NB) and flaring branch (FB). Type I X–ray burst properties correlate well with position in the CD (van der Klis et al. 1990). This led to argue that position on the CD along the atoll was a good indicator of $`\dot{M}`$ (but X–ray flux was not). For Z sources (a subset of which show sporadic X–ray bursts too), optical and UV emission as well as quasi periodic oscillations (QPOs) models indicate that $`\dot{M}`$ increases as HB$``$NB$``$FB (cf. van der Klis 1995).

RossiXTE observations have shed new light on these sources, revealing QPOs in the kHz frequency range (for a review see van der Klis 2000). In atoll sources, kHz QPOs occur in the island and lower banana part of the CD (Méndez 1998 and references therein). In the case of Z sources kHz QPOs are present all the way down the HB and disappear in the lower part of the NB (e.g. Wijnands & van der Klis 1998b and references therein). KHz QPOs are therefore detected in a narrow range of fluxes and in similar positions in the CDs. Whereas the disappearance of the kHz QPO activity at the low luminosity end may be related in some atoll sources to a lack of sensitivity, this is certainly not the case at the bright end.

In this letter we try to link the source behavior in the CD to the physics of the accretion process. We mainly concentrate on atoll sources, since for Z sources data are more sparse and they still lack an indication of the spin period from e.g. type I X–ray bursts.

## 2. Accretion regimes

Accretion phenomena onto a fast spinning, weakly magnetic neutron star three radii which do not vary with the mass inflow rate can be identified:

* the neutron star radius, $`R`$, which depends on the neutron star mass and equation of state. We scale here $`R=10R_6`$ km for a neutron star of mass $`M=2M_2M_{}`$;
* the marginally stable orbit radius, $`r_{\mathrm{ms}}`$. In the case of a rapidly rotating neutron star a small correction to the Schwarzschild value has to be included. As a first order in the specific angular momentum $`j`$ \[$`j=2\pi cI/(GM^2P)`$, with $`c`$ the speed of light, $`I`$ and $`P`$ the neutron star moment of inertia and spin period and $`G`$ the gravitational constant\], $`r_{\mathrm{ms}}(62.31j)GM/c^2`$ (Kluźniak 1998). In the case of ‘soft’ equations of state this radius is larger than $`R`$, whereas for stiff equation of states it is smaller and does not play any role. In the case of a spin frequency of 1.8 ms (as in the case of Aql X-1; Zhang et al. 1998a), and for a $`2M_{}`$ neutron star with $`I=2\times 10^{45}`$ g cm-2, $`r_{\mathrm{ms}}16M_2`$ km;
* the corotation radius, $`r_{\mathrm{cor}}`$, at which a test particle in Keplerian orbit corotates with the neutron star, $`r_{\mathrm{cor}}=(\frac{GMP^2}{4\pi ^2})^{1/3}`$. In the case of a spin frequency of 1.8 ms ($`P=1.8P_{1.8\mathrm{ms}}`$ ms) we have $`r_{\mathrm{cor}}=28M_2^{1/3}P_{1.8\mathrm{ms}}^{2/3}`$ km.

As the mass inflow rate changes the limiting radius at which the neutron star magnetic field pressure balances the ram pressure of the mass inflow towards the neutron star (i.e. the magnetospheric radius, $`r_\mathrm{m}`$) changes too. In the case of spherical accretion with a dipolar magnetic field $`r_\mathrm{m}(\frac{\mu ^4}{2GM\dot{M}^2})^{1/7}`$ (where $`\dot{M}`$ is the mass inflow rate onto the neutron star magnetosphere and $`\mu =BR^3`$ is the magnetic moment). A parameter $`\xi `$ is usually considered to account for the case of disk accretion, with $`\xi 0.51`$ (Ghosh & Lamb 1979, 1992). Wang (1995) pointed out that $`\xi `$ is in the range 0.87–0.95. In the following, we assume $`\xi =1`$. For an accretion luminosity of $`L=GM\dot{M}/R=10^{37}L_{37}\mathrm{erg}\mathrm{s}^1`$, we derive a magnetospheric radius of $`r_\mathrm{m}=20B_8^{4/7}R_6^{10/7}M_2^{1/7}L_{37}^{2/7}\xi `$ km (where $`B=B_8\mathrm{\hspace{0.17em}10}^8`$ G is the magnetic field). For higher luminosities the radiation pressure starts playing a role (region A in the standard disk model) and the magnetospheric radius changes to $`r_\mathrm{m}^\mathrm{A}=33B_9^{20/39}L_{38}^{2/13}M_2^{1/39}R_6^{18/13}\xi `$ km (cf. Ghosh & Lamb 1992; here we neglect the weak dependence on the viscosity parameter).

## 3. Application to Atoll sources

As a working example we can consider the case of the transient atoll source Aql X-1. This is the prototype of soft X–ray transient sources (Campana et al. 1998a) and, together with the accreting X–ray binary pulsar SAX J1808.4–3658 (Wijnands & van der Klis 1998a; Chakrabarty & Morgan 1998), is the only LMXRB which showed signs of centrifugal inhibition of accretion (Campana et al. 1998b; Zhang, Yu & Zhang 1998; Gilfanov et al. 1998). During a type I burst of Aql X-1 observed with RossiXTE, a coherent periodicity was discovered at 549 Hz (Zhang et al. 1998a). This frequency may correspond to the spin frequency of the neutron star or, like the case of 4U 1636–53, to twice its value (Miller 1999). In the following we adopt a spin period of 1.8 ms.

In Fig. 1 we show the CDs of Aql X-1 based on RossiXTE. Data fall into three main distinct groups, which can be readily identified with the extreme island, island and banana states based on the hard color. By adopting an absorbed blackbody plus power law model, Reig et al. (2000) derived the 2–10 keV X–ray luminosity in the different states of the CD diagram (for a distance of 2.5 kpc). This varied from $`2.0\times 10^{34}`$ erg s-1 in the extreme island state (not shown in Fig. 1) to $`3.2\times 10^{35}`$ erg s-1 in the island state. In the banana state, the luminosity increased from left to right, from $`2.4\times 10^{36}`$ erg s-1 to $`3.9\times 10^{36}`$ erg s-1 (Reig et al. 2000). Assuming the spectrum at a comparable luminosity observed with BeppoSAX (Campana et al. 1998b), we derive a bolometric correction of $`4`$ in the 0.1–40 keV energy band. We remark, however, that it has long been recognised that the observed flux does not provide a confident measure of the accretion luminosity. Different contributions may come from e.g. the neutron star rotational energy which is extracted by the disk (Priedhorsky 1986).

Despite the large variations in X–ray luminosity, kHz QPOs are observed in selected regions of the CDs. In the transient sources Aql X-1 and 4U 1608–52 (Méndez et al. 1999) there exist RossiXTE observations in which kHz QPO activity was not detected. The corresponding limiting RossiXTE PCA rates for the two sources are $`600`$ and $`400`$ c s-1, respectively (2–16 keV; Méndez 1998). Moreover, the persistent atoll source 4U 0641+091 showed kHz QPOs in the range 250–1800 c s-1 (2–60 keV; van Straaten et al. 2000). We conclude that, at least in the case of Aql X-1, the disappearance of the kHz QPO activity at the faint end is likely not due to poor statistics. For what concerns the maximum flux at which kHz QPOs have been detected, no problems of statistics exists ($`1500`$ and $`3000`$ c s-1 for Aql X-1 and 4U 1608–52, respectively, cf. Méndez 1998). The stopping occurs for all the atoll sources in the lower banana state.

The existence of selected regions in the CD (even if slightly different from source to source) in which kHz QPOs are detected is noteworthy. We now try to identify the onset and the disappearance of kHz QPO activity with changes in the magnetospheric radius.

($`i`$) As the accretion rate increases, the onset of the kHz activity can be related to the exit of the source from the propeller regime, when most of the matter is inhibited to accrete onto the neutron star surface by centrifugal forces (i.e. $`r_\mathrm{m}r_{\mathrm{cor}}`$) and therefore kHz QPOs can be hardly detected. The luminosity jump expected for the transition from the propeller to the accretion regime in the case of a ms spinning neutron star is just a factor of a few (Corbet 1996; Campana & Stella 2000). This flux reduction is independent on the mechanism invoked to explain the kHz QPOs, like beat-frequency models (BFMs; Strohmeyer et al. 1996; Miller, Lamb & Psaltis 1998) or general relativistic models (GRMs, Stella & Vietri 1998, 1999).

The minimum luminosity for which kHz QPOs have been detected during the fading end of an Aql X-1 outburst is $`1.1\times 10^{36}\mathrm{erg}\mathrm{s}^1`$ (2–10 keV; Zhang et al. 1998a). This luminosity is a factor of $`2`$ higher than the luminosity corresponding at the steep increase in the hardness ratio and associated to the centrifugal barrier onset as interpreted by Zhang, Yu & Zhang (1998) and Campana et al. (1998b). Considering the bolometric correction given above we assume a 0.1–40 keV luminosity of $`L_{\mathrm{min}}=4.2\times 10^{36}\mathrm{erg}\mathrm{s}^1`$. This luminosity represents an upper limit to the accretion luminosity (which defines the magnetospheric radius), but at the onset of the centrifugal barrier external contribution such as the luminosity extracted from the neutron star rotational energy should weaken considerably.

Equating the magnetospheric to the corotation radius, we derive an upper limit to the magnetic field value of:

$$B\stackrel{<}{}1\times 10^8L_{\mathrm{min}}^{1/2}R_6^{5/2}M_2^{1/3}P_{1.8\mathrm{ms}}^{7/6}\xi ^{7/4}\mathrm{G}.$$
(1)
In the case of Miller’s et al. BFM the relevant radius for kHz QPOs is the sonic radius, $`r_\mathrm{s}`$, which lies inside the magnetospheric radius. Since, $`r_\mathrm{s}<r_\mathrm{m}<r_{\mathrm{cor}}`$, the constraint above applies as well.

($`ii`$) The disappearance of the kHz QPOs at high luminosities may be related to the disappearance of the magnetosphere, when $`r_\mathrm{m}`$ becomes comparable to the neutron star radius, $`R`$. A similar line of reasoning has been applied by Cui et al. (1998) interpreting the transition from island and banana and the disappearance of the kHz QPO, as the disk becoming disangaged from the magnetic field. In particular, by comparing $`r_\mathrm{m}`$ with $`r_{\mathrm{ms}}`$ Cui et al. (1998) were able to limit the neutron star magnetic field to $`0.60.9\times 10^8`$ G.

The disappearance of the kHz QPOs can be interpreted in the framework of BFMs as the drop out of the magnetosphere without which the beating mechanism cannot work. In the Miller’s et al. BFM the relevant radius is the sonic point radius which, by definition, lies always inside the magnetospheric radius. On the other hand, the disappearance of the magnetosphere may leave no space for an orbiting blob to acquire the necessary eccentricity to generate kHz QPOs as predicted by GRMs (Stella & Vietri 1999) or if the magnetosphere itself is needed to generate the small eccentricity in a resonant interaction (Vietri & Stella 1998). In general both the sonic point model and the GRMs do not require the presence of a magnetosphere to work. However, we expect that the neutron stars in LMXRBs do have a magnetic field if they are the progenitors of millisecond radio pulsars. This link has been recently confirmed by the detection of coherent pulsations in the transient LMXRB SAX J1808–3658 (Wijnands & van der Klis 1998a).

In the case of Aql X-1, kHz QPOs stop to be detected at a luminosity higher than $`3.6\times 10^{36}\mathrm{erg}\mathrm{s}^1`$ (2–10 keV; Zhang et al. 1998a). Given the small range of luminosity we adopt the same bolometric correction, resulting in 0.1–40 keV luminosity of $`L_{\mathrm{max}}=1.4\times 10^{37}\mathrm{erg}\mathrm{s}^1`$. In this case, we are going to derive a lower limit on the magnetic field and no problems are posed by taking a higher luminosity than the true accretion one. We can obtain a conservative estimate of the magnetic field from $`r_\mathrm{m}\stackrel{<}{}R`$:

$$B\stackrel{>}{}3\times 10^7L_{\mathrm{max}}^{1/2}R_6^{3/4}M_2^{1/4}\xi ^{7/4}\mathrm{G}.$$
(2)
Different parameters can play a role in changing the exact value of these constraints as well as the uncertainties in the theory underlying the location of the magnetospheric radius. Despite these uncertainties, the two limits are well consistent and constrain the magnetic field of Aql X-1 to $`B0.31\times 10^8`$ G.

### 3.1. 4U 1820–30

A further case is represented by the atoll source 4U 1820–30 for which indication of a saturation in the kHz QPO frequencies versus the flux intensity has been observed (Zhang et al. 1998b; Kaaret et al. 1999 and references therein). 4U 1820–30 follows the general behaviour of atoll sources described above (Zhang et al. 1998b). The saturation starts near the lower banana onset (Zhang et al. 1998b). KHz QPOs appear at a 0.3–40 keV luminosity of $`5\times 10^{37}\mathrm{erg}\mathrm{s}^1`$ (for a distance of 8 kpc) and disappear at $`8\times 10^{37}\mathrm{erg}\mathrm{s}^1`$, i.e. spanning a vary small range of luminosities. We adopt these as bolometric luminosities and derive, as described above, $`B\stackrel{<}{}9\times 10^8`$ G (for a neutron star spin period of 3.6 ms, as inferred from the kHz QPO separation) and $`B\stackrel{>}{}0.8\times 10^8`$ G, respectively. A further constraint comes from the interpretation of the saturation of the kHz QPO frequencies as due to the magnetospheric radius reaching the last stable orbit. If this is the case, the upper kHz QPO frequency is related to the Keplerian frequency at the disk inner edge. In the case of GRMs, the ceasing of the kHz QPOs derives from the lack of space between the point at which the blob leaves the accretion disk and the neutron star.

A saturation value of 1065 Hz has been measured in 4U 1820–30, leading to an estimate of the neutron star mass of $`2.1M_{}`$. Alternatively, this implies a radius of the marginally stable orbit of $`r_{\mathrm{ms}}=17.4`$ km. The saturation occurs at a luminosity of $`6.6\times 10^{37}\mathrm{erg}\mathrm{s}^1`$ (0.3–40 keV), which provides a direct estimate of the magnetic field value of $`B2\times 10^8`$ G (close to the lower limit derived above as expected due to the fact that $`r_{\mathrm{ms}}R`$). In the case of Miller’s et al. BFM the estimate of the magnetic field gives only a lower limit, since the sonic radius is smaller than the magnetospheric radius.

The saturation in the kHz QPO frequencies has been observed only in 4U 1820–30. A possible explanation relies on the large neutron star mass inferred for this object, which in turn implies a smaller radius than for a $`1.61.8M_{}`$ neutron star by a factor of 10–20%. This might be enough to prevent the observation of the kHz QPO saturation in lighter neutron stars.

## 4. Application to Z sources

In Z sources kHz QPOs are usually observed in the upper part of the NB down to the end of the HB (e.g. Cyg X-2, Wijnands et al. 1998; GX 17+2, Wijnands et al. 1997). Flux determinations are very rare. Only in the case of GX 349+2 (Zhang, Strohmeyer & Swank 1997) the onset of kHz QPO occurs at a 2–10 keV luminosity of $`1.2\times 10^{38}\mathrm{erg}\mathrm{s}^1`$ (for a distance of 8.5 kpc). This is consistent with the idea that the lower apex of the Z path occurs when the source attains the Eddington luminosity. Despite the variations observed in the CDs, Z sources vary very little in luminosity, with typical RossiXTE PCA rates in excess of 1000 c s-1 (Wijnands & van der Klis 1998b) such that problems with a too poor statistics do not exist.

We consider the case of Cyg X-2, which showed kHz QPO down to apex NB/FB. The 0.1–20 keV luminosity extrapolated by EXOSAT observations at the NB/FB apex is $`L_{\mathrm{max}}2\times 10^{38}\mathrm{erg}\mathrm{s}^1`$ (at 8 kpc; Chiappetti et al. 1990). The luminosity is high and the inner disk regions are likely in the radiation pressure regime such that the magnetospheric radius $`r_\mathrm{m}^A`$ should be adopted. Using this luminosity to infer the lower bound on the magnetic field we obtain $`B\stackrel{>}{}1\times 10^8`$ G. This is a conservative lower limit since at these luminosities the coupling between disk and neutron star should be high. KHz QPOs have been observed at the end of the HB, for which the minimum luminosity observed with EXOSAT is just a factor of 2 smaller. From this we derive $`B\stackrel{<}{}8\times 10^8`$ G, assuming a neutron star spin period of 2 ms, still waiting for the discovery of a coherent periodicity from bursting Z sources. Higher fields may be obtained by assuming faster spin periods.

## 5. Discussion

We interpret the onset and the end of kHz QPOs in terms of the variable extent of the neutron star magnetosphere. The turn-on of kHz QPOs at low fluxes is interpreted as the starting of accretion onto the neutron star surface, overcoming the effects of the centrifugal barrier (here we do not claim that this barrier is working perfectly, but if it can stop the large fraction of the infalling matter, the flux reduction will prevent to detect kHz QPOs). The disappearance of the kHz QPOs corresponds instead to the ‘disruption’ of the magnetosphere on the neutron star (see also Cui et al. 1998). Given these assumptions, we estimate in the case of the atoll source Aql X-1 a magnetic field of $`B0.31\times 10^8`$ G. Similar results can be obtained for other atoll sources. Interpreting the saturation of the kHz QPO frequency for increasing luminosities in 4U 1820–30, leads to a magnetic field estimate of $`2\times 10^8`$. In the case of Z sources we derive a rough estimate based on luminosity obtained with EXOSAT. In the case of Cyg X-2, we obtain a very conservative lower limit of $`B18\times 10^8`$ G. This is consistent with the idea that Z sources have higher magnetic field than atoll sources.

This scenario allows for clear predictions on the luminosities at which the kHz QPOs set in and out, which can be verified with RossiXTE monitoring of low mass X–ray binaries.

I thank P. Reig for permission of publishing Fig. 1. I acknowledge useful discussions with T. Belloni, M. Méndez & L. Stella.