Table Of ContentRandom L´evy Matrices Revisited
Zdzis law Burda1,2, Jerzy Jurkiewicz1,2, Maciej A. Nowak1,2, Gabor Papp3, and Ismail Zahed4
1Marian Smoluchowski Institute of Physics, Jagiellonian University, 30-059 Krak´ow, Reymonta 4, Poland
2Mark Kac Complex Systems Research Centre, Jagiellonian University, Krak´ow, Poland
7
3Institute of Physics, E¨otv¨os University, P´azm´any P.s.1/a, H-1117 Budapest, Hungary
0
0 4Department of Physics and Astronomy, SUNY Stony Brook, NY 11794, USA
2
We compare eigenvalue densities of Wigner random matrices whose elements are independent
n identically distributed (iid) random numbers with a L´evy distribution and maximally random ma-
a
trices with a rotationally invariant measure exhibiting a power law spectrum given by stable laws
J
of free random variables. Wecompute the eigenvalue density of Wigner-L´evy (WL) matrices using
9 (and correcting) the method by Bouchaud and Cizeau (BC), and of free random L´evy (FRL) ro-
2
tationally invariant matrices by adapting results of free probability calculus. We compare the two
types of eigenvalue spectra. Both ensembles are spectrally stable with respect to the matrix addi-
]
h tion. The discussed ensemble of FRL matrices is maximally random in thesense that it maximizes
c Shannon’s entropy. We find a perfect agreement between thenumerically sampled spectra and the
e analytical results already for matrices of dimension N = 100. The numerical spectra show very
m weakdependenceonthematrixsizeN ascanbenoticedbycomparingspectraforN =400. Aftera
pertinentrescalingspectraofWigner-L´evymatricesandofsymmetricFRLmatriceshavethesame
-
t tailbehavior. Aswediscusstowardstheendofthepaperthecorrelationsoflargeeigenvaluesinthe
a
two ensembles are however different. We illustrate the relation between the two types of stability
t
s and show that the addition of many randomly rotated Wigner-L´evy matrices leads by a matrix
.
t central limit theorem to FRL spectra, providingan explicit realization of themaximal randomness
a
principle.
m
- PACSnumbers: 02.50.-r,02.60.-x,89.90.+n
d
n
o
I. INTRODUCTION
c
[
Applications of random matrix theory cover many branches of physics and cross-disciplinary fields [1] involving
3
multivariate analysis of large and noisy data sets [2]. The standard random matrix formulation belongs to the
v
7 Gaussianbasin, with a measure that is Gaussianor polynomialwith finite secondmoment. The ensuing macroscopic
8 spectral distribution is localized with finite supports on the real axis. The canonical distribution for a Gaussian
0 measure is Wigner’s semi-circle.
2 The class of stable (L´evy) distributions [3] is however much larger (the Gaussian class represents only one fixed
0
point in the stability basin of the L´evy class), and one is tempted to ask, why the theory of randomL´evy matrices is
6
not so well established. The case of L´evy randomness is far from being academic, and many distributions in physics
0
and outside (finance, networks) exhibit power-like behavior referred to as fat or heavy tails [4].
/
at One of the chief reasons for why the theory of random L´evy matrices is not yet well understood is the technical
m difficulty inherent to these distributions. First, even for one-dimensional stable distributions, the explicit form of the
probabilitydistributionfunctions(pdf) isknownanalyticallyonlyinfew cases[5,6]. Second,L´evydistributionshave
-
d divergent moments, and the finiteness of the second moment (condition for the Gaussian stability class) is usually a
n key for many of the techniques established in Random Matrix Theory. Third, numerical studies involving power-like
o
behavior on infinite supports require enormous statistics and are very sensitive to systematic errors.
c
In 1994,Bouchaudand Cizeau [7] (hereafter BC) consideredlarge N N randomsymmetric matrices with entries
v: sampled from one-dimensional, stable distributions. In the large N lim×it they obtained analytical equations for the
i
X entries of the resolvent, and then checked their predictions for the spectra using numerically generated spectra by
random sampling. The agreement was fair, although not perfect. Contrary to the standard Gaussian-like ensembles,
r
a the measure in the BC approach was not rotationally invariant.
In2002,followingthe workin[8],we havesuggestedanotherL´evy-typeensemble[9](hereafterFRL).Byconstruc-
tionitsmeasureisrotationallyinvariant. Theaveragespectraldistributioninthisensembleisstableundermatrixthe
convolutionoftwoindependentbutidenticalensembles. Itissimilartothestabilitypropertyofone-dimensionalL´evy
distributions. The measure is non-analytic in the matrix H and universalat largeH with a potential V(H) lnH2.
≈
This weak logarithmic rise in the asymptotic potential is at the origin of the long tail in the eigenvalue spectra.
Thepresentworkwillcomparethe WLandFRLresultsasadvertisedin[10]. Insection2wereanalyzeandcorrect
the originalargumentsforthe resolventpresentedin[7]. Ourintegralequationsforthe resolventandspectraldensity
are different from the ones in [7]. We carry explicit analytical transformations and expansions to provide insights to
2
the spectrum. We show that there is a perfect agreement between the analytical and numerical results obtained by
sampling large L´evy matrices. We also discuss the relation of ours and BC’s results. In section 3 we recall the key
conceptsbehindFRL ensembles. InlargeN the resolventobeysasimple analyticequation. The resultingspectraare
compared to the spectra following from the corrected BC analysis. The WL and FRL matrices represent two types
of stability under matrix convolution. In both cases we have a power behavior in the tails of the spectrum. By a
pertinent rescaling we may in fact enforce the same tail behavior and compare the spectra. The observed differences
disappear in the Gauss limit and become more pronounced in the Cauchy limit. In sectin 4 we explain the relation
between the two types of stability on a simple example: we construct sums of WL matrices rotatedby randomO(N)
matrices and show that the spectrum of these sums converges by a matrix central limit theorem to the pertinent
symmetric FRL spectrum. Our conclusions are in section 5.
II. WIGNER-LE´VY MATRICES
Definition of the ensemble
In a pioneering study on random L´evy matrices, Bouchaud and Cizeau [7] discussed a Wigner ensemble of N N
×
realsymmetricrandommatriceswithelementsbeingiidrandomvariables: withprobabilitydensityfunctionfollowing
a L´evy distribution P(x) N1/µLC,β N1/µx , with µ being the stability index, β – the asymmetry parameter, and
≡ µ
C – the range of the distribution (see below). We shall call these matrices Wigner-L´evy (WL) or Bouchaud-Cizeau
(cid:0) (cid:1)
(BC) matrices. The probability measure for the ensemble of such matrices is given by:
dµ (H)= P(H )dH (1)
WL ij ij
i≤j
Y
ThescalingfactorN1/µ inpdfmakesthelimiting eigenvaluedensityindependentofthematrixsizeN whenN .
Alternatively one can think of the matrix elements H as if they were calculated as H = h /N1/µ with h →bei∞ng
ij ij ij ij
iid random numbers independent of N: p(x)=LC,β(x).
µ
L´evy distributions are notoriously hard to write explicitly (except in few cases), but their characteristic functions
are more user friendly [5]
1
LC,β(x)= dkLˆ(k)eikx (2)
µ 2π
Z
where the characteristic function is given by
logLˆ(k)= C k µ(1+iβ sign(k) tan(πµ/2)). (3)
− | |
The parameters µ, β and C are related to the asymptotic behavior of LC,β(x)
µ
C(1 β)
lim LC,β(x)=γ(µ) ± (4)
x→±∞ µ xµ+1
| |
with the µ-dependent parameter γ(µ) given by
πµ
γ(µ)=Γ(1+µ)sin( ) (5)
2
Here µ is the stability index defined in the interval (0,2], 1 β 1 measures the asymmetry of the distribution
and the range C > 0 is the analogue of the variance, in a s−ense≤tha≤t a typical value of x is C1/µ. A standard choice
corresponds to C =1.
We shallconsider here only the stability index in the range (1,2), althoughas will be shownlater, results obtained
in this range seem to be valid also for µ=1. We also assume that all random variables have zero mean.
Determination of the eigenvalue density
Amethodofcalculatingthe eigenvaluedensityofthe Wigner-L´evymatriceswasinventedbyBouchaudandCizeau
[7]. Let us in this section briefly recall the main steps of the method. It is convenient to introduce the resolvent,
called also Green’s function:
1
g(z)= TrG(z) (6)
Nh i
3
where elements of the matrix G(z) are
G (z)=(z H)−1 (7)
ij − ij
andtheaveragingiscarriedoutusingthe measure(1). Theresolventcontainsthesameinformationastheeigenvalue
density ρ(λ). Indeed if one approaches the real axis one finds that ρ(λ) = 1/πlim Im g(λ+iǫ). This is how
ǫ→0+
−
one usually calculates ρ(λ) from g(z). For the Wigner-L´evy ensemble, where individual matrix elements have large
scale-free statistical fluctuations a slightly different method turns out to be more practical – a method which allows
one to avoidproblems in taking the double limit (first N and ǫ 0+) in which the fluctuations are suppressed
→∞ →
in an uncontrollable way in the presence of an imaginary part of z. In the BC method z is kept on the real axis and
fluctuations are not suppressed, so one can safely take the large N limit.
The method goes as follows [7]. One first generates a symmetric N N random matrix H using the measure (1)
×
andthenbyinvertingH z onecalculatestheresolventG(z)(7). Nextoneaddsanewrow(andasymmetriccolumn)
−
of independent numbers identically distributed as those in the old matrix H. One obtains a new (N+1) (N+1)
matrixH+1,where+1emphasizesthatithasonemorerowandonemorecolumnthanH. Itisconvenientt×onumber
elements of the original matrix H by indices running over the range i,j = 1,...,N, and assign the index 0 to the
ij
new row and the new column so that now indices of H+1 run over the range 0,...,N. If one inverts the matrix
H+1 z one obtains a new (N+1) (N+1) resolvent G+1(z). One can show that it obeys a recursive relation [7]
− ×
N N
1 h h h G (z)
00 0i 0j ij
z =H + H H G (z)= + . (8)
− G+1(z) 00 0i 0j ij N1/µ N2/µ
00 i,j i,j
X X
which relates the element G+1(z) of the (N+1) (N+1) resolventto the elements of the old N N resolventG(z).
One can also derive similar e00quations for off-diag×onal elements of G+1(z): ×
G+1(z) N h G (z)
0j 0i ij
= (9)
G+1(z) N2/µ
00 i
X
The difference between the probability distribution of the elements of G+1(z) and of G(z) dissapears in the limit
N . The diagonal elements of the matrix G+1(z) are identically distributed as the diagonal elements of the
→ ∞
matrix G(z). The same holds for off-diagonal ones. Moreover in this limit all elements of the G matrix become
independent of each other. In particular all diagonal elements of G(z) become independent identically distributed
(iid) random variables as N . One can use equations (8) and (9) to derive self-consistency equations for the
→ ∞
probabilitydensityfunction(pdf) fordiagonalandthepdfforoff-diagonalelements. Wearehereprimarilyinterested
in the distribution of the diagonal elements, as we shall see below. The self-consistency equation for the probability
distributionofdiagonalelementsfollowsfromtheequation(8)andisindependentofthedistributionoftheoff-diagonal
elements. This can be seen as follows. Let us first define after Bouchaud and Cizeau [7] a quantity:
1
S (z)=z (10)
0 − G+1(z)
00
It is merely a convenient change of variables suited to the left hand side of equation (8). It is clear that if one
determines the pdf for S = S (z), one will also be able to determine the pdf for G = G+1(z) since the two pdfs can
0 00
be obtained from each other by the change of variables (10):
1 1
P (G)= P z (11)
G G2 S − G
(cid:18) (cid:19)
where P andP arepdfs for G+1(z) andS (z) respectively. Since the probability distribution P is identicalfor all
G S 00 0 G
diagonal elements of G, it remains to determine the probability distribution P for the quantity S (z).
S 0
Let us sketch how to do that. First observe that the first term in the equation (8) can be neglected at large N, so
the equation assumes the form:
N h2 G (z) h h G (z)
S (z)= 0i ii + 0i 0j ij . (12)
0 N2/µ N2/µ
i i6=j
X X
Now observe that by construction h are independent of G (z). We shall now show that for large N the first term
0i ij
in (12) dominates over the second one. We shall modify here the argument used in [7], where it was assumed that
the off-diagonal elements G N(z), i = j are suppressed by a factor 1/N1/µ with respect to the diagonal elements.
ij
6
4
Instead we note that by construction the quantities G (z) are statistically independent of h , i = 0,...,N. Since
ij 0i
we shall only be interested by the diagonal elements of the resolvent matrix, we may replace the contribution of the
off-diagonal elements h h G (z), i = j by their averaged values. As a result, the contribution of the off-diagonal
0i 0j ij
6
terms averages out. Note that this is also true in the Gaussian limit µ = 2 where following [7] we may replace all
elements of(12) bytheir respectiveaverages. Takingthis into accountandomitting the subleadingcontributionH ,
00
we get
N h2 G (z)
S (z)= 0i ii . (13)
0 N2/µ
i
X
Thus the problem was simplified to an equation where the left hand side (S = z 1/G+1) and the right hand side
0 − 00
depend only of the diagonal elements of the G-matrix, which as we mentioned before, are identically distributed in
the limit N . Using (13) one canderivea self-consistencyequationfor the probability density function (pdf) P
G
→∞
for diagonal elements of the matrix G.
Generalized central limit theorem
To proceed further we apply with Bouchaud and Cizeau [7] the generalized central limit theorem to derive the
universal behavior of the sum on the right hand side of (13) in the limit N :
→∞
i. If the h ’s are sampled from the L´evy distribution LC,β, the squares t = h2 for large t are distributed
• 0i µ i 0i i
solely along the positive real axis, with a heavy tail distribution:
Cdt
i
γ(µ) (14)
∝ t1+µ/2
i
irrespective of β. The sum
N
t
i
(15)
N2/µ
i
X
is distributed following LC′,1(t ). The range parameters C and C′ are related by (3)
µ/2 i
2C′γ(µ/2)=Cγ(µ) (16)
The factor 2 on the left hand side corresponds to the sum (1+β)+(1 β) appearing as a contribution from
−
positiveandnegativevaluesoftheoriginaldistributionofh . Thisrelationisimportantwhencomparingtothe
0i
numerical results below where C = 1 is used. From now on (and to simplify the equations) we assume instead
that C′ =1.
ii. By virtue of the central limit theorem the following sum
•
N
G (z)t
ii i
(17)
N2/µ
i
X
of iid heavy tailed numbers t is for G (z) = (N0) and N L´evy distributed with the pdf: LC(z),β(z),
i ii O → ∞ µ/2
which has the stability index µ/2 and the effective range C(z) and the asymmetry parameter β(z) calculated
from the equations:
N
1
C(z)= G (z)µ/2 (18)
ii
N | |
i
X
and
1 N G (z)µ/2sign(G (z))
β(z)= N i | ii | ii . (19)
1 N G (z)µ/2
P N i | ii |
which follow from thePcomposition rules for the tail amplitudes of iid heavy tailed numbers t defined above.
i
5
Integral Equations
We sawinthe previoussectionthatthe generalizedcentrallimittheoremimpliesfor largeN thatthe “self-energy”
S = S (z) is distributed according to the L´evy law P (S) = LC(z),β(z)(S) with the stability index µ/2 being a half
0 S µ/2
of the stability index of the L´evy law governing the distribution of individual elements of the matrix H, and with
the effective range parameter C(z) and the asymmetry parameter β(z) which can be calculated from the equations
(18) and (19), respectively. One should note that the effective parameters C(z),β(z) of the distribution P (S) are
S
calculated for C = 1 and that they are independent of β of the probability distribution: LC,β(h ) of the H-matrix
µ ij
elements.
The sums on the right hand side of the two equations (18), (19) for C(z) and β(z) have a common form
1 f(G (z)). Since in the limit N , the diagonal elements become iid, the sums on can be substituted
N i ii → ∞
by integrals over the probability density for G :
ii
P
N
1 dG 1
f(G (z)) = dGP (G)f(G)= P z f(G) (20)
N h ii i G G2 S − G
i Z Z (cid:18) (cid:19)
X
where in the second step we used the equation (11). Since the distribution P (S) = LC(z),β(z)(S) is known up to
S µ/2
the values of two effective parameters C(z) and β(z) the equations (18) and (19) can be written as self-consistency
relations for β(z) and C(z):
∞ dG
C(z) = Gµ/2LC(z),β(z)(z 1/G)
− G2| | µ/2 −
Z−∞
∞ dG Gµ/2sign(G)LC(z),β(z)(z 1/G)
β(z) = −−∞ G2| | µ/2 − . (21)
R ∞ dG Gµ/2LC(z),β(z)(z 1/G)
−−∞ G2| | µ/2 −
R
The symbol stands for principal value of the integral. Notice here the difference between our second equation and
−
that in [7]. In addition to [7] we also note that the resolvent takes the form:
R
1 ∞ dG
g(z)= G (z) g(z)= G LC(z),β(z)(z 1/G) (22)
N ii −→ − G2 µ/2 −
i Z−∞
X
The integrals (21) and (22) can be rewritten using the new integration variable x=1/G as
+∞
C(z) = dxx−µ/2LC(z),β(z)(z x)
| | µ/2 −
Z−∞
+∞dx sign(x)x−µ/2LC(z),β(z)(z x)
β(z) = −∞ | | µ/2 − (23)
R +∞dxx−µ/2LC(z),β(z)(z x)
−∞ | | µ/2 −
and R
∞ dx
g(z)= LC(z),β(z)(z x) (24)
− x µ/2 −
Z−∞
All the steps above require both z and G (z) to be strictly real. The argumentcannot be extended to the complex z
ii
plane. So all integrals above should be interpreted as principal value integrals, wherever it is necessary. Since µ< 2
all integrals in (23) are convergentalso in the usual sense.
The equation for g(z) can be rewritten as
∞ dx
g(z)= LC(z),β(z)(x). (25)
− z x µ/2
Z−∞ −
Notice the nontrivialdependence on z in the parameters of the L´evy distribution. Above equation can be a source of
confusion, since its structure resembles another representationof g(z)
∞ dλ
g(z)= ρ(λ) (26)
− z λ
Z−∞ −
6
whichsuperficiallylooksasifonecouldidentify in(25)xandλ andLC(z),β(z)(x)withρ(λ). This isnotthe case,and
µ/2
one should instead invert (25) using the inverse Hilbert transform:
1 ∞ dz
ρ(λ)= g(z). (27)
π2 − z λ
Z−∞ −
In other words, one first has to reconstructnumerically the real partof the resolvent,and only then compute numer-
ically the spectral function ρ(λ), using the ”dispersive relation” (27). This is a difficult and rather subtle procedure.
In the next section we give some analytical insights to the integral equations that would help solve them and extract
the spectral function.
Analytical properties useful for numerics
One cannot do the integrals from previous section analytically. As mentioned, one cannot even write down an
explicit form of the L´evy distribution. It is a great numerical challenge to solve the problem numerically even if
all the expressions are given. One realizes that already when one tries to compute the Fourier integral (2) of the
characteristicfunction since one immediately sees that the integrandin the form (2) is a strongly oscillating function
making the numerics unstable. Fortunately using the power of the complex analysis one can change this integral to
a form which is numerically stable. So in this section we present some analytic tricks which allow one to reduce the
problem of computing the eigenvalue density as formulated in the previous section to a form which is well suited to
the numerical computation.
To simplify (23) we proceed in steps. First, we make use of the L´evy distribution through its characteristic
LC,β(x)= 1 ∞ dkeikxe−C|k|µ/2(1+iζ sign(k)). (28)
µ/2 2π
Z−∞
with ζ =β tan(πµ/4). By rescaling through
k = C−2/µk′, (29)
x = C2/µx′,
z = C2/µz′,
we can factor out the range LC,β(x)=C−2/µL1,β(x′). Second, we make use of the following integrals,
µ/2 µ/2
∞ dx
eikx = 2ieikzsignk,
− z x −
Z−∞ −
∞ cos(k′x′)
dx′ = k′ µ/2−1Γ(1 µ/2)sin(πµ/4), (30)
x′µ/2 | | −
Z0
∞ sin(k′x′)
dx′ = k′ µ/2−1sign(k′)Γ(1 µ/2)cos(πµ/4).
x′µ/2 | | −
Z0
Last, we make use of the change of variables p=k′µ/2. With this in mind, we obtain
4 µ πµ ∞
C2(z′) = Γ 1 sin dpcos(p2/µz′ ζ(z′) p)e−p, (31)
πµ − 2 4 −
∞d(cid:16)psin(p2(cid:17)/µz′(cid:16) ζ(z(cid:17)′)Z0p)e−p
ζ(z′) = 0 − , (32)
∞dpcos(p2/µz′ ζ(z′) p)e−p
R0 −
with ζ(z′) = tan(Rπµ/4) β(z′). For every z′ we can iteratively solve the equation for ζ(z′), then we determine C(z′)
andusez =C2/µ(z′)z′ toexpresseverythingintermsofz. Thesetransformationssolve(23). Usingthesamemethod
we rewrite the equation for g(z) as
2 ∞
g¯(z′)=C(z′)2/µg(z)= dp p(2−µ)/µsin(p2/µz′ p ζ(z′)) e−p. (33)
µ −
Z0
7
These integral forms are useful to study the small-z′ limit. In this case ζ(z′) is an antisymmetric function of z′
and has an expansion in powers of z′. Using ζ(z′) = k z′+ (z′3) we can recursively obtain the coefficients of this
1
O
expansion. The first term is k =Γ(1+2/µ)/2. Similarly C(z′) is a symmetric function in z′
1
4 µ πµ
C2(z′)= Γ 1 sin + (z′2) (34)
πµ − 2 4 O
(cid:16) (cid:17) (cid:16) (cid:17)
For z′ large(z′ )adifferentapproachisneeded. InthiscasewefollowNolan[11]andtreatthetwointegrals
| | →±∞
(numerator and denominator) of (32) together, i.e. we consider the integral
∞
dpe−h(p), (35)
Z0
where
h(p)=(1 iζ)p+iz′p2/µ. (36)
−
Nolan’s idea is to close the contour of integration in the complex p plane in the following way: at p we add an
arc and afterwards continue until p = 0 along the line where h(p)=0. Using the parameterization→p =∞reiθ we get
ℑ
the parametric equation for r(θ) along this line, valid for z′ >0
µ/(2−µ)
sin(θ θ)
0
r(θ)= − . (37)
z′cos(θ0)cos(µ2θ)!
The angle θ for the curve we need (µ 1) is bounded between θ = arctan(ζ(z′)), where r = 0 and θ = πµ/4,
0 1
where cos(2θ) is zero. Let us now intro≥duce a new variableψ through θ =θ ψ and 0 ψ θ +θ −. In thi−s range
µ 0− ≤ ≤ 0 1
we have
µ/(2−µ)
sin(ψ)
V(ψ) = , (38)
cos(θ0)cos(µ2(ψ−θ0))!
h(ψ) = 1 µ/(2−µ)V(ψ) cos(2−µµψ− µ2θ0) .
ℜ z′ cos(θ )cos(2(ψ θ ))
(cid:18) (cid:19) 0 µ − 0
After some manipulations we obtain
θ0+θ1 1 µ/(2−µ)
sin(p2/µz′ ζ(z′) p)e−p = dψ V(ψ)e−ℜh(ψ) (39)
− z′
Z Z0 (cid:18) (cid:19)
2 cos(2−µµ(ψ−θ0)) µ sinθ0
× 2−µ cos(µ2(ψ−θ0) − 2−µ sinψ!
θ0+θ1 1 µ/(2−µ)
cos(p2/µz′ ζ(z′) p)e−p = dψ V(ψ)e−ℜh(ψ)
− z′
Z Z0 (cid:18) (cid:19)
2 sin(2−µµ(ψ−θ0)) + µ cosθ0 .
× 2−µ cos(µ2(ψ−θ0)) 2−µ sinψ !
The resulting integrals look complicated, however they contain both the small-z′ and the large-z′ asymptotics.
For z′ 0 we have ψ = z′p2/µ−1, which reproduces the small-z′ expansion presented above. For z′ we have
ψ =θ →+θ u/z′µ/2. Note that in this limit → ∞
0 1
−
2 2 u
cos( (ψ θ ))=sin( ) (40)
µ − 0 µz′µ/2
and the z′ dependence in Reh(ψ) vanishes in leading order.
The large-z′ asymptotics requires some work. For the leading orders we have
8
0.2
0.15
)
λ
(
ρ
0.1
0.05
0
-4 -2 0 2 4
λ
FIG. 1: Theoretical (black) and numerical (red) eigenvalue distributions for µ=1.95
sin(p2/µz′ ζ(z′) p)e−p Γ(1+µ/2)sinθ (z′)−µ/2 (41)
1
− ≈
Z
cos(p2/µz′ ζ(z′) p)e−p Γ(1+µ/2)cosθ (z′)−µ/2.
1
− ≈
Z
Thus for z′ we have
→∞
ζ(z′)=tan(θ )+ (1/z′µ/2) (42)
1
O
and
C(z′) = z′−µ/4(1+ (1/z′µ/2)) (43)
O
z = √z′(1+ (1/z′µ/2)).
O
All the formulas above apply to the case z′ > 0. One can also derive similar formulas for z′ < 0, it is however more
practical to use the symmetry properties of the functions C(z′) and ζ(z′).
As a final check let us compute g¯(z′). A rerun of the above transformations on the integrals give
2 θ0+θ1 1 2/(2−µ)
g¯(z′) = dψ V(ψ)2/µe−Reh(ψ) (44)
2 µ z′
− Z0 (cid:18) (cid:19)
2 1 + sin(2−µµψ− µ2θ0) .
× µcos(µ2(ψ−θ0)) sinψ !
Notice that both asymptotics follow from this representation. For z′ we have g¯(z′)=1/z′+ . which implies
→∞ ···
g(z)=1/z+ . This can be viewed as a check of the correct normalization of the eigenvalue density distribution.
···
Numerical Comparison
In this section we show perfect agreement between the theoretical analysis for the eigenvalue distribution (27)
and the numerically generated eigenvalue distribution. The latter is obtained by diagonalizing N N random L´evy
×
matricessampledusingthe measure(1). Theformerwasgeneratedbycalculatingnumericallyg(z)asdetailedabove.
We performed numerically the inverse Cauchy transform (27). It is important to note that the integral transforms
9
0.25
0.2
0.15
)
λ
(
ρ
0.1
0.05
0
-7.5 -5 -2.5 0 2.5 5 7.5
λ
FIG. 2: Theoretical (black) and numerical (red) eigenvalue distribution for µ=1.75
0.25
0.2
)
λ( 0.15
ρ
0.1
0.05
0
-10 -5 0 5 10
λ
FIG. 3: Theoretical (black) and numerical (red) eigenvalue distribution for µ=1.50
0.35
0.3
0.25
) 0.2
λ
(
ρ
0.15
0.1
0.05
0
-15 -10 -5 0 5 10 15
λ
FIG. 4: Theoretical (black) and numerical (red) eigenvalue distribution for µ=1.25
10
0.5
0.4
) 0.3
λ
(
ρ
0.2
0.1
0
-10 -5 0 5 10
λ
FIG. 5: Theoretical (black) and numerical (red) eigenvalue distribution for µ=1.00
entering into the definition of the eigenvalue density ρ(λ) in (27) converge slowly for z . For that, we have used
→∞
the asymptotic expansion of g(z) to perform the large-z part of the integrals.
As noted above, all the above analytical results were obtained using a specific choice of the scale factor (16). The
comparisonwiththenumericallygeneratedeigenvaluedistributiongeneratedusingL1,0 distributionrequiresrescaling
µ
through λ φλ and ρ(λ) ρ(λ)/φ with
→ →
Γ(1+µ)cos(πµ) 1/µ
φ= 4 (45)
Γ(1+ µ)
(cid:18) 2 (cid:19)
Below we show a sequence of results for µ = 1.95, 1.75, 1.50, 1.25 and 1.00 with this rescaling. The comparison
is for high statistics 400 400 samples (red). We have checked that the convergence is good already for 100 100
× ×
samples, with no significant difference between N =100,200,400. The numerical results are also not sensitive to the
choiceβ =0. Theagreementbetweenthe resultsfollowingfromthe integralequationsandthe numericallygenerated
6
spectra is perfect. This is true even for µ = 1, where in principle the arguments used in the derivation may not be
valid.
Numerical observation
In the mean-field approximation[7] one can assume that there are no correlations between large eigenvalues of the
Wigner-L´evy random matrix. In this case the eigenvalue density takes the form:
ρ (λ)=LC(λ),β(λ)(λ). (46)
µ µ/2
It is natural to ask how good this mean-field approximation is. This can be done by comparing the mean-field
b
eigenvalue distribution ρ(λ) (46) to the eigenvalue distribution ρ(λ) calculated by the inverse Hilbert transform (27)
of the resolvent g(z) (25) as we did in previous section. We made this comparison numerically. The result of this
numericalexperimentwasthatwithinthenumericalaccuracywhichweachievedthetwocurvesrepresentingρ(λ)and
b
ρ(λ) lied on top of each other in the whole studied range of λ. Since our numericalcodes are written in Mathematica
we could push the numerical accuracy very far, being only limited by the execution time of the code. We have not
b
seen any sign of deviation between the shapes of the two curves. This provides us with a strong numerical evidence
that the mean-field argument [7] gives an exact result but so far we have not managed to prove it. The value of the
eigenvaluedensityρ (0)forλ=0canbecalculatedanalyticallyforthemean-fielddensity(46). Rescalingthedensity
µ
ρ (λ) ρ (φλ)/φ by the factor φ (45) we eventually obtain
µ µ
→
b
b b Γ(1+2/µ) Γ(1+µ/2)2 1/µ
ρ (0)= (47)
µ
π Γ(1+µ)
(cid:18) (cid:19)
b
