Domain wall entropy of the bimodal two-dimensional Ising spin glass
Domain wall entropy of the bimodal two-dimensional Ising spin glass
A. Aromsawa and J. Poulter
Department of Mathematics, Faculty of Science, Mahidol University, Rama 6 Road, Bangkok 10400, Thailand
(Dated: November 15, 2018)
We report calculations of the domain wall entropy for the bimodal two-dimensional Ising spin
glass in the critical ground state. The L× L system sizes are large with L up to 256. We find that
it is possible to fit the variance of the domain wall entropy to a power function of L. However,
the quality of the data distributions are unsatisfactory with large L > 96. Consequently, it is not
possible to reliably determine the fractal dimension of the domain walls.
PACS numbers: 75.10.Hk, 75.10.Nr, 75.40.Mg, 75.60.Ch
I. INTRODUCTION
The short-range Ising spin glass is still a source of con-
troversy in spite of its comparative simplicity. In brief,
the exact mechanisms by which widely separated spins
influence each other are not clearly understood. The sys-
tem Hamiltonian, due to Edwards and Anderson1, is
H = −
Jijσiσj (1)
where the nearest-neighbor exchange interactions Jij are
quenched random variables which, for present purposes,
take random sign with equal probability. The two partic-
ular models of disorder widely studied are the bimodal,
or ±J , model where the interactions have fixed magni-
tude and the Gaussian model where they are taken from
a continuous normal distribution of zero mean. In two di-
mensions it is generally accepted that the spin glass only
exists at zero temperature2,3, that is the critical temper-
ature Tc = 0.
Spin correlations in two dimensions for bimodal disor-
der are expected to decay algebraically according to
[< σ0σR >
2]av ∼ R−η (2)
in the ground state. To date, the best estimates from
Monte Carlo simulation4 and exact calculation5 seem to
agree that the exponent is η = 0.14. For Gaussian dis-
order, the exponent is zero since the ground state, apart
from a global inversion, is unique. In contrast, bimodal
disorder comes with a large ground state degeneracy cor-
responding to an entropy of 0.07k per spin6,7,8,9,10,11. In
spite of this it has been argued12 that η = 0 in the ther-
modynamic limit. Very recently, Hartmann13 has esti-
mated η = 0.22 from scaling arguments. This is in fair
agreement with a number of previous estimates7,9,14,15,16.
Nevertheless, the important issue here is universality;
that is whether or not η is positive.
Clear evidence of the importance of long range influ-
ence in the case of bimodal disorder is apparent in the
current controversy regarding the lowest excited state.
This issue dates back to the work of Wang and Swendsen7
which proposed that the energy gap should be 2J , not
4J , in the thermodynamic limit. Although some work4,9
has supported the simple naive 4J scenario, it is now be-
coming clearer that strong disorder is indeed producing
a different result for the infinite system. The energy gap
has been often reported11,17,18 as 2J . Further, it has even
been suggested, based on a finite size scaling study of the
correlation length and the spin glass susceptibility12, that
the specific heat should obey a power law. This may indi-
cate a feature universal with Gaussian disorder. A good
review of the issues involved here is given in Ref. 19.
Essentially it appears to be the case that there exist low
energy excitations with very long range influence.
A standard technique commonly used to investigate
the long range response of a spin glass is to introduce a
domain wall defect14,20,21. This is done, in two dimen-
sions, by drawing a line across the system and reversing
the signs of all bonds cut by the line. With Gaussian
disorder, the ground and low energy excited states are
unique. Consequently, the term domain wall has a clear
physical meaning and appears as an optimium fractal
path that corresponds to the lowest excited state. The
case of bimodal disorder is quite different as a result of the
large degeneracy of both the ground and excited states.
The degeneracy of the ground state is of the order of
exp(0.07L2) for a L × L lattice and that of the first ex-
cited state is probably much larger still11. In particular,
it is often the case that the defect system is not an ex-
cited state of the reference system. Thus, we do not have
real domain walls in the same sense as with the case of
Gaussian disorder. Nevertheless, we might expect some
useful knowledge to be derived from a study of the bi-
modal system in response to reversing the signs of bonds
along a line as with the Gaussian system.
Thinking of the thermodynamic limit, we can discuss
this issue in the light of droplet theory22. For the case of
a continuous defect distribution, Gaussian for instance,
there is a unique ground state and we define a droplet as
a region bounded by a closed path (or surface in three
dimensions). Below the transition temperature, which is
positive in three dimensions at least23, the scaling prop-
erties of droplets of various sizes can be investigated by
reversing all spins inside the closed path. This creates
an excited state of the reference system and it is known
that there exist low energy droplet excitations of large
spatial extent. The idea here is that a domain wall de-
http://arxiv.org/abs/0704.1186v2
fect is much the same as a large droplet excitation in the
thermodynamic limit. For bimodal disorder this compar-
ison is less clear due to the huge ground state degeneracy.
Certain droplets may not represent excitations at all.
The most important prediction of droplet theory is
that the energy difference Edw between the two systems,
with and without the domain wall, can be fitted to a
power function of the system size L. We have, for the
spin glass,
< |Edw| >∼ Lθ (3)
where θ is known as the spin-glass stiffness exponent.
In three dimensions, this exponent is found to be
positive24,25, at about 0.2 for both bimodal and Gaus-
sian disorder, showing consistency with the existence of
a stable spin glass phase at finite temperature. As a mat-
ter of fact the droplet theory22 was originally developed
for this type of case where both the critical temperature
and stiffness exponent are positive. It is also probable
that the Gaussian and bimodal models fit into the same
universality class23 with respect to their transitions at fi-
nite critical temperature. Degeneracy is not an issue due
to the thermal fluctuations.
In two dimensions it was not immediately obvious that
droplet theory22 is entirely appropriate since the criti-
cal temperature is zero and the stiffness exponent is not
positive. For Gaussian disorder the stiffness exponent is
negative26, θ ≈ −0.28, clearly indicating that the spin
glass is unstable at any finite temperature. It is also re-
markable that the values of θ obtained from domain wall
and droplet calculations are in very good agreement. For
bimodal disorder it is generally accepted that the stiffness
exponent due to domain wall defects is zero2,27,28. Never-
theless the situation for droplets is not clear. Hartmann13
has estimated the droplet stiffness exponent by construct-
ing ground and first excited states of three models of bi-
modal disorder and found good agreement with Gaussian
disorder. Nevertheless, universality is not shown since
the correlation function exponent η is not zero.
For the case of Gaussian disorder, in two dimensions,
it has been suggested29,30 that the domain walls are
stochastic Loewner evolution processes31. This theory
is able to relate the domain wall fractal dimension df to
the stiffness exponent via
df = 1 +
4(3 + θ)
and the result df ≈ 1.27 agrees well with the literature26.
There is also an interesting conjecture due to Fisch32 that
θ = (
6 − 3)/2 exactly. Nevertheless, there is no good
reason to believe that Eq. (4) can be used with bimodal
disorder, possibly due to the degeneracy of the ground
state.
The domain wall entropy Sdw is defined in the same
manner as Edw. Droplet theory
22 predicts that Sdw
should take values with random sign and large variance.
In particular, the variance is predicted to scale as
< S2dw > − < Sdw >2∼ Ldf (5)
which provides a possible means to estimate df . In par-
ticular, we can use this to test the appropriateness of
droplet theory for the case where both the critical tem-
perature and the stiffness exponent tend to zero from
above. This is precisely the situation we have for the
two-dimensional model with bimodal disorder.
Previous estimations of the fractal dimension from di-
rect studies of the domain wall entropy have been pub-
lished. Saul and Kardar9 predict df ≈ 1.0. Fisch33 ar-
gues that df might be an increasing function of |Edw|.
The possibility that df = 1.25 in agreement with Eq.
(4) is not ruled out. Finally, Lukic et al34 have reported
df = 1.03(2). These values should also be compared
with those from topological analysis of the ground state.
Romá et al35 report df = 1.30(1) while Melchert and
Hartmann36 find an interval 1.095(1) ≤ df ≤ 1.395(1).
In this article we report calculations of the domain wall
entropy on sample sizes that are much larger than any-
thing done before. Furthermore, our method is applied
at an arbitrarily small temperature and there is no need
to extrapolate to the ground state. Our main conclusion
is that the domain wall fractal dimension for bimodal
disorder, as predicted by droplet theory, is not a well de-
fined quantity. The reason for this is that there exists
no clear prescription for its estimation if the domain wall
entropy distributions are significantly far from normal. A
brief overview of the method is given in Sec. II. This is
followed by our results in Sec. III and a brief discussion
in Sec. IV.
II. BACKGROUND
The planar Ising model is known to be isomorphic
to a system of noninteracting fermions. One particular
illustration37 has been adapted by Blackman38 for disor-
dered systems. For the square lattice, each site is deco-
rated with four fermions. Equivalently, we can decorate
each bond with two fermions. For a system of N lattice
sites, we have 4N fermions in total. It is useful to think
of the two fermions decorating a bond to be placed one on
either side. In this way a plaquette (square) is decorated
with four fermions; left, right, top and bottom.
The partition function for the Ising model on a square
lattice with any set of exchange interactions takes the
Z = 2N [
cosh(Jij/kT )] (detD)
1/2 (6)
where the product is over all nearest neighbor bonds Jij
on the N site lattice and D is a skew-symmetric (4N ×
4N) matrix. The square root of the determinant of D
is also called the Pfaffian. Essentially, it represents the
sum over all closed lattice polygons and is equal to the
product of all the positive eigenvalues of D.
The calculation of the partition function with bimodal
disorder has been described in much detail previously8
and a simple summary should suffice here. At zero tem-
perature, D is a singular matrix with a set of degener-
ate zero eigenvalues exactly equal in number to the total
number of frustrated plaquettes. In order to extract the
physics of the system, degenerate state perturbation the-
ory is applied at an arbitrarily low temperature. The
defect eigenvalues occur in pairs and approach zero as
some power of exp(−2J/kT )
ǫ = ±1
X exp(−2Jr/kT ) (7)
where r is an integer (an order of perturbation theory)
and X is a real number that is independent of tempera-
ture and depends only on the configuration of frustrated
plaquettes. The ground state energy and entropy can be
expressed exactly as
E = −2NJ + 2J
rd (8)
S = k
lnXd (9)
where the sums are over all defect eigenstate pairs.
To summarise the perturbation theory, we first write
the matrix D exactly as the sum of two terms D =
D0 + δD1, where δ = 1 − t with t = tanh(J/kT ). Of
course t = 1 and δ = 0 in the ground state. Both of
the matrices D0 and D1 are independent of temperature.
The matrixD0 has eigenvectors localised inside each frus-
trated plaquette; expanded in the basis of the four deco-
rating fermions. It is these localised states that form the
defect basis for the perturbation theory. The matrix D1
is 2 × 2 block diagonal in the pairs of fermions decorat-
ing the bonds (one fermion either side). All degeneracy
at first order is lifted by diagonalizing D1 in the defect
basis. For example, we can think of just two neighboring
frustrated plaquettes. The perturbation theory gives one
defect pair with rd = 1 and Xd = 1.
In general, the first order calculation will leave some
zero eigenvalues. The corresponding eigenvectors of D1
form the basis for second order. We can imagine a sys-
tem of two next nearest neighbor frustrated plaquettes.
Clearly rd = 2 and Xd = 1 or 2 depending on the ar-
rangement. In order to show this we need to use the
continuum Green’s function5,8 gc = gc1 + gc2. The ma-
trix gc1 is 4×4 block diagonal in the four fermions inside
each plaquette and clearly allows us to connect two frus-
trated plaquettes, across an unfrustrated plaquette. The
second order calculation is performed by diagonalizing
D2 = D1gc1D1. The matrix gc2 is, just like D1, 2 × 2
block diagonal in the pairs of fermions decorating bonds.
A proof that gc2 is irrelevant for the ground state has
been given in Ref. 5.
For higher orders we require Green’s functions for
states whose degeneracy has already been lifted. We de-
fine, for r ≥ 1,
Gr = −
N(r)∑
| r, i〉(1/ǫir)〈r, i | (10)
where | r, i〉 denotes state i (with eigenvalue ǫir) in the set
of states whose degeneracy was lifted at order r; there
are N(r) of these states. At third order the matrix to be
diagonalized is D3 = D2(1 +G1D1)gc1D1 and, generally
at arbitrary order Dn = Dn−1(1 + Gn−2Dn−2) · · · (1 +
G1D1)gc1D1. The perturbation theory is applied order
by order until all degeneracy is lifted.
The scheme outlined above allows exact calculations
of energy and entropy in the ground state. The method
is fully gauge invariant in that it depends only on the
number and distribution of frustrated plaquettes. Fur-
thermore, there is no requirement to extrapolate to the
ground state. The Pfaffian is not calculated at any par-
ticular numerical value of the temperature. We believe
that this method is the best available for calculating the
ground state entropy of large lattices, although matching
algorithms are better for the energy2.
We have used periodic boundary conditions in one di-
mension. The cylindrically wound frustrated patch was
nested in an unfrustrated system of infinite extent in the
second dimension. In this scheme, the introduction of a
domain wall defect is particularly simple. The two pla-
quettes at the ends of the defect, one on each side of the
patch, change their status; frustrated to unfrustrated or
vice versa. For a perfect ferromagnet this gives a domain
wall energy proportional to L as required. For a fully
frustrated system the defect would make no real differ-
ence since the domain wall energy would be independent
of L. For the spin glass the domain wall energies are all
multiples of 4J since we have L even. The probability
that a plaquette is frustrated is expected to be close to
0.5 and it is conceivable that the system with the domain
wall could be interchanged with its reference system in
another realization of disorder. This is consistent with
the prediction of droplet theory that the domain wall
entropy has random sign.
Since the domain wall entropy is generally a small dif-
ference between two larger numbers, we have taken great
care with the floating point computations. Although our
method is analytically exact, it is subject to numerical
propagation error on the computer. Ill-conditioned dis-
order realizations were detected by calculating the cor-
relation function along a path around the cylinder and
repeated in arbitrary precision arithmetic as necessary.
III. RESULTS
We have calculated the domain wall entropy for bi-
modal disorder on L × L lattices where L = 8, 12, 16,
24, 32, 48, 64, 96, 128, 192 and 256. For sizes L ≤ 128,
105 random samples were taken. We also took 4 × 104
for L = 192 and 104 for L = 256.
 500 100 50 10 5
FIG. 1: The domain wall energy (in units of J) as a function
of system size L. The error bars are two standard deviations
of the mean and the curve is a best fit to the form A−BL−p
following Ref. 2.
To establish the credentials of our boundary scheme,
we have calculated the domain wall energy. The data is
shown in Fig. 1 where the error bars are two standard
deviations of the mean. Following Ref. 2, we have fitted
the data to a function A − BL−p and find A = 3.7(2),
B = 3.2(1) and p = 0.23(4). We note that the fit is ap-
proaching saturation from below. We believe that this is
as a consequence of our boundary conditions. The prob-
ability of a zero energy (Edw = 0) domain wall is found
to decrease with L, contrary to the situation with free
boundary conditions in the unwound dimension39. Fur-
thermore, since we only use even values of L, the defect
energies are all multiples of 4J . The quality of the non-
linear fit40 is Q = 0.91. Attempts to fit a power law were
not successful. For instance a fit for L ≥ 96 completely
missed the point at L = 256. We conclude that our
method is reliable although larger system sizes L > 256
are required to conclude more convincingly that θ = 0.
The variance of the domain wall entropy Sdw is shown
in Fig. 2. A power law fit for 24 ≤ L ≤ 96 predicts
that, according to Eq. (5), the fractal dimension of the
domain walls is df = 1.090(8) where the quality of the
fit is Q = 0.42. A similar fit (Q = 0.30) of < |Sdw| >2
gives 1.080(9) which agrees well, indicating good quality
data distributions.
A second power law fit for L ≥ 96, also shown in Fig.
2, reveals a significantly higher value df = 1.30(3) with
Q = 0.16. However, a fit of < |Sdw| >2 gives only
1.23(2). Although the quality Q = 0.05 is lower, the
difference is too large to be disregarded. The reason for
 1000
 500 100 50 10 5
FIG. 2: The variance of the domain wall entropy as a function
of system size L. The error bars are two standard deviations
of the mean. Two power law fits are shown; for 24 ≤ L ≤
96 (dashed line) and L ≥ 96 (dotted line). The powers are
respectively 1.090(8) and 1.30(3).
this discrepancy must lie in the quality of the distribu-
tions for Sdw. For L = 256, for example, the distribu-
tion has skewness 0.64 and kurtosis 2.06. Although the
mean < Sdw >= 3.53 is still much less than the variance
(≈ 350) it reflects a significant sign disparity. The bias
most likely arises due to correlations in the distribution of
frustrated plaquettes. Incidently, the corresponding dis-
tributions for Edw are of excellent quality. For L = 256,
the skewness and kurtosis are respectively −0.006 and
0.07.
IV. DISCUSSION
In summary, the distributions of the domain wall en-
tropy for large L, with bimodal disorder in two dimen-
sions, are found to deviate significantly from normal. In
consequence, even if we assume that droplet theory22 is
appropriate, it is unable to prescribe exactly how to get
the domain wall fractal dimension. This does not neces-
sarily mean that droplet theory is entirely wrong. It just
does not give the whole story, only an approximation,
for this system; having a large ground state degeneracy,
a zero critical temperature and a zero stiffness exponent.
Of course, it is quite probable that corrections to scal-
ing are large and difficult to manage. This is actually a
rather likely scenario in view of the poor results for fitting
the ground state energy with cylindrical winding in one
dimension28. Scaling corrections are an issue probably
related to strong correlations in the distribution of frus-
trated plaquettes for large L. All gauge invariant quanti-
ties like entropy depend only on the frustrated plaquette
distribution; nothing else. The prediction of droplet the-
ory that the domain wall entropy is normally distributed
with zero, or very small, mean and large variance proba-
bly relates to an assumption that a defect system occurs
as a reference system in another realization of disorder.
This assumption may not be true if the frustrated pla-
quette distributions are strongly correlated, as is likely
in view of the anamolous behaviour of the degeneracy of
the first excitations mentioned earlier. A further scenario
is that due to Hartmann13, which proposes that droplet
theory is actually inappropriate for estimating the stiff-
ness exponent of domain wall defects. If this is correct,
it is also unlikely that the fractal dimension of domain
walls has anything to do with droplet theory. Neverthe-
less, our results do indicate an approximate appropriate-
ness for droplet theory in the sense that all the possible
estimates for df do not seem unreasonable.
Previous studies of the domain wall entropy have
worked with much smaller system sizes. Saul and
Kardar9 had sizes up to L = 36 and found df ≈ 1.0,
while Lukic et al34 used sizes up to L = 50 and fitted
< |Sdw| > to find df = 1.03(2). Fisch33 used sizes up to
L = 48 and has introduced the idea that the domain wall
entropy may be significantly correlated with energy. It
is argued that an effective df increases as a function of
|Edw| and convergence to the value 1.25 consistent with
Eq. (4) is not ruled out. We have tested these predictions
and find that, for Edw = 0, the response is in fact much
stronger in both the intermediate and large size regimes.
For L ≥ 96 we find df = 1.43(4) while, for Edw 6= 0,
df = 1.22(3). Also, the probability of finding a disor-
der realization with Edw = 0 is under 0.5 for L > 24,
much less than 0.75. The cause of these discrepancies is
most probably due to boundary conditions and system
size. We do not have any evidence from our work that
the particular values of Edw are significant for the droplet
theory.
We also note that a topological analysis of ground
states35 predicts df = 1.30(1). Essentially, the technique
measured the average length of domain walls. These
lengths respond faster than just L since the domain walls
bend to avoid the rigid lattice. Nevertheless, only one
ground state configuration was studied for each disor-
der realization, completely ignoring the entropy issue. A
study in a similar vein36 looks at the properties of mini-
mum energy domain walls and places the fractal dimen-
sion in an interval 1.095(1) ≤ df ≤ 1.395(1). This may
agree to some extent with the point that it is not pos-
sible, or very difficult, to actually pin down the value of
Acknowledgments
We would like to thank A. K. Hartmann for an en-
lightening correspondence. A. Aromsawa acknowledges
support from the Thailand Research Fund in the form
of a scholarship from the Royal Golden Jubilee Ph.D.
Programme.
1 S. F. Edwards and P. W. Anderson, J. Phys. F: Met. Phys.
5, 965 (1975).
2 A. K. Hartmann and A. P. Young, Phys. Rev. B 64,
180404(R) (2001).
3 J. Houdayer, Eur. Phys. J. B 22, 479 (2001).
4 H. G. Katzgraber and L. W. Lee, Phys. Rev. B 71, 134404
(2005).
5 J. Poulter and J. A. Blackman, Phys. Rev. B 72, 104422
(2005).
6 H.-F. Cheung and W. L. McMillan, J. Phys. C 16, 7027
(1983).
7 J.-S. Wang and R. H. Swendsen, Phys. Rev. B 38, 4840
(1988).
8 J. A. Blackman and J. Poulter, Phys. Rev. B 44, 4374
(1991).
9 L. Saul and M. Kardar, Phys. Rev. E 48, R3221 (1993);
Nucl. Phys. B 432, 641 (1994).
10 J. A. Blackman, J. R. Gonçalves, and J. Poulter, Phys.
Rev. E 58, 1502 (1998).
11 J. Lukic, A. Galluccio, E. Marinari, O. C. Martin, and G.
Rinaldi, Phys. Rev. Lett. 92, 117202 (2004).
12 T. Jörg, J. Lukic, E. Marinari, and O. C. Martin, Phys.
Rev. Lett. 96, 237205 (2006).
13 A. K. Hartmann, arXiv:0704.2748 (2007).
14 A. J. Bray and M. A. Moore, in Heidelberg Colloquium
on Glassy Dynamics, edited by J. L. van Hemmen and I.
Morgenstern, Lecture Notes in Physics Vol. 275 (Springer-
Verlag, Berlin, 1987), p.121.
15 R. N. Bhatt and A. P. Young, Phys. Rev. B 37, 5606
(1988).
16 Y. Ozeki, J. Phys. Soc. Jpn. 59, 3531 (1990).
17 J.-S. Wang, Phys. Rev. E 72, 036706 (2005).
18 H. G. Katzgraber, L. W. Lee, and I. A. Campbell,
cond-mat/0510668 (2005).
19 H. G. Katzgraber, L. W. Lee, and I. A. Campbell, Phys.
Rev. B 75, 014412 (2007).
20 W. L. McMillan, Phys. Rev. B 29, 4026 (1984); Phys. Rev.
B 30, 476 (1984); Phys. Rev. B 31, 340 (1985).
21 A. J. Bray and M. A. Moore, J. Phys. C 17, L463 (1984);
Phys. Rev. B 31, 631 (1985); Phys. Rev. Lett. 58, 57
(1987).
22 D. S. Fisher and D. A. Huse, Phys. Rev. Lett. 56, 1601
(1986); Phys. Rev. B 38, 386 (1988).
23 H. G. Katzgraber, M. Körner, and A. P. Young, Phys. Rev.
B 73, 224432 (2006).
24 A. K. Hartmann, Phys. Rev. E 59, 84 (1999).
25 M. Palassini and A. P. Young, J. Phys. Soc. Jpn. 69, 165
(2000).
26 A. K. Hartmann and A. P. Young, Phys. Rev. B 66, 094419
(2002).
27 C. Amoruso, E. Marinari, O. C. Martin, and A. Pagnani,
Phys. Rev. Lett. 91, 087201 (2003).
28 I. A. Campbell, A. K. Hartmann, and H. G. Katzgraber,
Phys. Rev. B 70, 054429 (2004).
29 C. Amoruso, A. K. Hartmann, M. B. Hastings, and M. A.
Moore, Phys. Rev. Lett. 97, 267202 (2006).
30 D. Bernard, P. Le Doussal, and A. A. Middleton,
cond-mat/0611433 (2006).
http://arxiv.org/abs/0704.2748
http://arxiv.org/abs/cond-mat/0510668
http://arxiv.org/abs/cond-mat/0611433
31 J. Cardy, Ann. Phys. (N.Y.) 318, 81 (2005).
32 R. Fisch, arXiv:0705.0046 (2007).
33 R. Fisch, J. Stat. Phys. 125, 793 (2006).
34 J. Lukic, E. Marinari, O. C. Martin, and S. Sabatini, J.
Stat. Mech. L10001 (2006).
35 F. Romá, S. Risau-Gusman, A. J. Ramirez-Pastor, F. Ni-
eto, and E. E. Vogel, Phys. Rev. B 75, 020402(R) (2007).
36 O. Melchert and A. K. Hartmann, arXiv:0704.2004 (2007).
37 H. S. Green and C. A. Hurst, Order-Disorder Phenomena
(Interscience, London, 1964).
38 J. A. Blackman, Phys. Rev. B 26, 4987 (1982).
39 R. Fisch and A. K. Hartmann, Phys. Rev. B 75, 174415
(2007).
40 W. H. Press, S. A. Teukolsky, W. T. Vettering, and B. P.
Flannery, Numerical Recipes in Fortran (Cambridge Uni-
versity Press, Cambridge, 1992).
http://arxiv.org/abs/0705.0046
http://arxiv.org/abs/0704.2004
