## Abstract

It was demonstrated in Atas & Bogomolny (2012 *Phys. Rev. E* **86**, 021104) that the ground-state wave functions for a large variety of one-dimensional spin- models are multi-fractals in the natural spin-*z* basis. We present here the details of analytical derivations and numerical confirmations of this statement.

## 1. Introduction

Many-body quantum systems are a recurring subject of theoretical physics. Today's increase in computer power makes it possible to treat problems with a few tens of particles numerically, which opens up new possibilities for their investigation.

Here, we consider one of the oldest many-body interacting models, namely one-dimensional spin chains. The archetypical example is the XYZ Heisenberg model [1] for *N* spins in external fields,
1.1and its various specifications for different values of the parameters. Here are the usual Pauli matrices at site *n*.

Many different methods were developed to find the spectra of such Hamiltonians [2–5]. The calculation of wave functions is more difficult even for integrable models. In the natural basis of *z*-component of each spin, with *σ*_{j}=±1, any wave function of spin can be represented as a formal sum
1.2where the summation is taken over all *M*=2^{N} configurations, with *N* being the total number of spins.

In such a basis, Hamiltonians like (1.1) are represented by an *M*×*M* matrix. In general, the coefficients are obtained by matrix diagonalization. As the size of the Hilbert space grows exponentially with the number of spins, exact diagonalization quickly becomes intractable, and an iterative method of diagonalization is eventually necessary. Even in integrable cases, the wave functions of spin chains still require knowledge of an exponentially large number of coefficients, which look quite erratic (see later diagrams) and their structure is not well understood.

In Atas & Bogomolny [6], it was shown that the ground-state (GS) wave functions for spin chains are multi-fractals on the spin-*z* basis. The multi-fractality as usual (e.g. [7,8]) means that moments of wave functions scale non-trivially with the Hilbert space dimension *M*,
1.3and *τ*(*q*)=*D*_{q}(*q*−1), where *D*_{q} are called fractal dimensions.

More precisely, let *S*_{R}(*q*,*M*) be the Rényi entropy [9] for an eigenfunction (1.2) of an *M*×*M* matrix,
1.4with normalized coefficients , .

Fractal dimensions, *D*_{q}, are determined from the behaviour of Rényi entropy (1.4) in the limit [10]:
1.5If this limit is zero, *D*_{q}=0, only a small number of coefficients give large contributions and one refers to such functions as localized functions. In the opposite situation, when almost all coefficients are of the same order, , it is plain that *D*_{q}=1. Such functions are called delocalized ones. (For *d*-dimensional systems fully delocalized functions have *D*_{q}=*d*.) In the general case, *D*_{q} has a nonlinear dependence on *q* and the corresponding functions are labelled as multi-fractals.

Multi-fractality is a general notion introduced to characterize quantitatively strong and irregular fluctuations of various quantities [7–11]. It appears in very different physical contexts, from turbulence [12,13] to human heartbeat dynamics [14]. Wave function multi-fractality attracted wide attention when it was recognized that it appears at special points in some condensed matter problems (e.g. at the point of metal–insulator transition in the three-dimensional Anderson model [15] and in two-dimensional quantum Hall effect [10]). Simpler critical random matrix models whose eigenfunctions have multi-fractal properties consist of matrices whose off-diagonal elements decrease as the first power of the distance from the diagonal [16,17]: *M*_{ij}∼|*i*−*j*|^{−1} when |*i*−*j*|≫1.

In all such critical models, multi-fractality is a non-trivial consequence of the concurrence between delocalization due to spreading and localization due to randomness. In spin models (1.1), there are no random parameters and the fact that the GS wave functions are multi-fractal may seem strange. Nevertheless, by combining numerical and analytical calculations, it has been proved in [6].

The notion of multi-fractality, similar to localization, depends on the basis. A function may be localized on one basis and delocalized on another. Spin chains are defined in spin basis (cf. (1.1)) and we investigate multi-fractal properties only on this basis though for certain problems the use of another basis, for example the fermion one [4], may be useful.

The purpose of this paper is twofold. First, we present certain details of the analytical calculations in [6]. Second, by mere definition (1.5), the fractal dimensions are defined by the limiting procedure. As there is practically no case where *D*_{q} are known analytically, their careful determination is sensitive to the method of interpolation used to calculate the limit from the data with relatively small *N*. This question is important in applications but rarely discussed in the literature. We collect data obtained by the direct diagonalization of Hamiltonian matrices for models with up to 13 spins and the ones from the iterative Lanczos algorithm [18] up to 19 spins and compare different methods of interpolation.

The plan of this paper is as follows. We start §2 with an informal introduction to the multi-fractal formalism based on an example of the binomial measure. Though it is the simplest mathematical example of multi-fractality, it appears that there exist particular cases of spin chains (1.1) where the GS has exactly the same structure as the binomial measure. In §3, two standard numerical methods of finding GS functions, the direct diagonalization and the Lanczos method, are briefly discussed. Section 4 is devoted to the investigation of the quantum Ising model (QIM) in a transverse field, which is one of the most studied spin chain models. By combining analytical and numerical calculations, we prove that its GS wave function is multi-fractal. A generalization of the Ising model, namely, the XY model, is discussed in §5. Similar to the Ising model, this model is also integrable, which permits certain fractal dimensions to be found analytically. Special attention is given to the so-called factorizing field where the XY model has exact and simple factorizing GS wave function. It is in this case that the GS wave function can be described by the binomial cascade, thus proving rigorously the multi-fractality of this wave function. In §6, properties of GS wave functions for the XXZ and XYZ models are briefly discussed. Section 7 concludes the paper. For clarity, we choose the parameters such that all terms in the Hamiltonian are non-positive. Owing to the Perron–Frobenius theorem, this implies that expansion coefficients of the GS wave function can be chosen non-negative. For the QIM and the XY models, we impose that their GSs are ferromagnetic but for the XXZ and XYZ models we choose them anti-ferromagnetic.

## 2. Binomial measure

To get an insight into the multi-fractality, we consider first the simplest example of the multi-fractal measures built by iterating a procedure called a (binary) multiplicative cascade [7].

In the first step of the cascade, the unit interval is divided into two equal subintervals and one associates a mass *m*_{0} (a measure) to the left subinterval [0,1/2] and a mass *m*_{1} to the right subinterval [1/2,1]. The two positive numbers *m*_{0,1} are such that *m*_{0}+*m*_{1}=1 (convenient parametrization is and ). At stage 2 of the iteration, we apply the same procedure to the two previous subintervals and the measures associated with the four subintervals are
2.1This process can be visualized graphically in a binary tree (figure 1).

After *N* iterations, one gets *M*=2^{N} subintervals of the form
2.2Let us associate a variable *σ*_{j}=−1 for the branch *j* on the binary tree when it turns to the left and *σ*_{j}=1 if it turns to the right. Now each interval (2.2) corresponds uniquely to with *σ*_{j}=±1 and the measure of such an interval is
2.3where is the number of left moves, i.e. the number of occurrences of −1 in the vector .

Iteration of this procedure generates an infinite sequence of measure, which converges to the binomial measure
2.4Figure 2*a* illustrates the function obtained after 12 iterations for . The abscissa (with ) is related to a binary code of binomial cascade as follows:
2.5

As the mass is conserved during the cascade, i.e. the sum of elements along a horizontal line of the tree equals 1 (owing to the binomial theorem, hence the name), measure (2.3) is a probability measure.

### (a) Fractal dimensions

Let *S*_{R}(*q*,*M*) be Rényi entropy (1.4) for a discrete normalized probability distribution after *N* steps (*M*=2^{N}),
2.6For the binomial measure, one gets
2.7where *n* is the number of left moves along the binary tree (or number of occurrences of −1 in sequence ). Therefore,
2.8Here, are the usual binomial coefficients.

Fractal dimensions, *D*_{q}, are defined from the limiting behaviour of the Rényi entropy in the limit (1.5). For the binomial measure
2.9The form of this function is shown in figure 2*b*. In particular, for large *q*, fractal dimensions tend to the limits
2.10As we shall see below, the curves of *D*_{q} as a function of *q* for all spin chain models considered in this paper have the same characteristic shape.

### (b) Singularity spectrum

The example of the binomial measure is also convenient for informal discussion of other quantities of interest in multi-fractal formalism [7,8,10]. The most important of them is the singularity spectrum, commonly denoted by *f*(*α*). To clarify physical interpretation of formulae below, it is useful (but, of course, not necessary) to define fictitious ‘energy levels’ as follows:
2.11Then expression (1.3) takes the form
2.12This sum can physically be interpreted as the canonical partition function for a system with total energies and temperature 1/*q*. As usual in thermodynamics, one may introduce the density of energy levels as
2.13where *S*(*E*) is the entropy. Then, the previous sum reduces to
2.14By analogy with thermodynamics, it is natural to assume that the entropy *S*(*E*) as well as the total energy *E* are extensive functions of the number of particles *N*. This leads to the following expression for the entropy,
2.15with a certain function *f*(*α*). (In application to fractal dimensions, one has to insist that , where *M* is the dimension of the Hilbert space. This introduces a constant factor between *N* and the number of spins and will result in redefinition of *f*(*α*). We shall implicitly take into account this fact.)

Using (2.15), the integral in (2.14) at large *N* can be calculated by the saddle point method, and it is clear that
2.16For the binomial cascade *τ*(*q*) is known (cf. (2.9)), , and *f*(*α*) can be calculated implicitly by the Legendre transform
2.17In figure 2*c*, the singularity spectrum *f*(*α*) for the binomial cascade with is presented. For other models, *f*(*α*) has a similar form [10].

In the above thermodynamic language *τ*(*q*) and *f*(*α*) play the role, respectively, of the free energy and the entropy per particle. Of course, the multi-fractal formalism itself [7,8] does not require references to thermodynamics. But we think that the latter, being well known, sheds a particular light on the notion of multi-fractality and clarifies the use and the meaning of different quantities.

In principle, *D*_{q} and *f*(*α*) contain the same information, but in this paper we shall focus on the calculation of fractal dimensions only.

## 3. Numerical approaches

For completeness, we briefly discuss the main numerical methods of calculation of the GS wave function for spin models as in (1.1). A basis state of the chain will be denoted with
3.1and |*σ*_{j}〉 the eigenstate of . The dimension of the Hilbert space is *M*=2^{N}.

### (a) Exact diagonalization

To achieve exact diagonalization, one has to write the Hamiltonian matrix on the basis (3.1), and hence sort these vectors. A way to do that is to use the binary code of a configuration: when a spin is *up* (resp. *down*), we replace it by 1 (resp. −1). For *N*=2, one has, for instance,
3.2Consequently, there is a bijection between all integer numbers between 0 and 2^{N}−1 and configurations of the spin chains. The element of the Hamiltonian can then be written as
3.3The terms in (1.1) that contain and flip two adjacent spins at sites *n* and *n*+1. The term flips the spin at site *n* only. These terms give rise to off-diagonal elements of the Hamiltonian matrix. The terms and give contributions to its diagonal part.

Once the Hamiltonian matrix has been defined, one can use a standard diagonalization library and compute eigenvalues and eigenvectors. Typically, these algorithms take a time proportional to *M*^{3}, so that one can only attain easily spin chains up to 13 spins. We shall see later that this is sufficient to compute fractal dimensions with relatively good precision. However, it is necessary to adopt an iterative (and so approximate) method of diagonalization in order to reach spin chains with a large number of spins.

### (b) Lanczos technique

The Hamiltonian matrix is a very sparse matrix: it contains in each row and each column *K*∼*N*≪*M*=2^{N} non-zero matrix elements. The Lanczos algorithm is an iterative algorithm based on power methods [18]. It permits one to find the lowest (or the largest) eigenvalues and corresponding eigenvectors of a square matrix and is particularly useful for finding decompositions of very large sparse matrices.

Basically, the Lanczos algorithm constructs a special basis, the so-called Krylov space, where the Hamiltonian has a tridiagonal representation. The algorithm works as follows (e.g. [19]). Select an arbitrary vector |*ϕ*_{0}〉 in the Hilbert space of the model being studied. Generally, this vector is chosen randomly to ensure that the overlap with the actual GS is non-zero. If some symmetries of the GS are known, then it is convenient to initiate the iterations with a state already belonging to the subspace having those quantum numbers. We then create a new vector by applying to |*ϕ*_{0}〉 and subtracting its projection over |*ϕ*_{0}〉,
3.4which satisfies 〈*ϕ*_{0}|*ϕ*_{1}〉=0. Doing the same thing with |*ϕ*_{1}〉, we construct a new state that is orthogonal to the two previous ones:
3.5

Generalizing to the order *n*, we have
3.6where and the coefficients are given by
3.7with the condition *b*_{0}=0. It is easy to check that the vector created at step *n* is orthogonal to the *n*−1 previous ones. After *n* iterations, the Hamiltonian matrix has the following tridiagonal form on the basis {|*ϕ*_{k}〉}_{0≤k≤n}:

Once in this form, the matrix can be diagonalized easily using standard library subroutines. A number of iterations equal to the size of the Hilbert space are necessary to diagonalize exactly as the model being studied. However, a number *n*∼100≪*M* iterations is sufficient in practice to have good enough accuracy for the lowest energy states. The GS wave function in (1.2) is expressed in the {|*ϕ*_{k}〉}_{0≤k≤n} basis as and the coefficient *c*_{m} is obtained during the diagonalization of *T*_{n}. By this method, with desktop computers one can easily find the lowest states for spin chains up to 20 spins.

## 4. Quantum Ising model

The QIM in a transverse field [20] is a well-studied model of quantum phase transitions [21]. It is defined by the Hamiltonian (1.1) with *Δ*=*α*=0 and *γ*=1,
4.1We consider the case of ferromagnetic Ising model with periodic boundary conditions, *σ*_{N+1}=*σ*_{1}.

### (a) Analytical results

The spectrum of this model can be found analytically by the Jordan–Wigner transformation [4], which maps the spin problem into a free fermions problem (figure 3).

The eigenenergies of the QIM are given by the fermionic filling of one-particle levels,
4.2where *n*_{k}=0,1 is the number of fermions with energy *e*_{k} given by
4.3with *k* taking *N* values depending on the parity of the excitation
4.4where *n*_{down} is the number of spins down.

One has [4]
4.5The GS energy corresponds to all *n*_{k}=0.

The calculation of the GS eigenfunction is more involved owing to the necessity to perform the Bogoliubov transformation [4,22–24]. For each sequence of *σ* coefficients *Ψ*_{σ} in (1.2) for the GS they are given by the determinant of an *N*×*N* matrix,
4.6where the orthogonal matrix *O*_{mn} [4,22] is
4.7and the angle *θ*_{k} (*θ*_{−k}=*θ*_{k}) is given by
4.8For the GS in the QIM, and the summation over *k* indicates the sum over all *l*=0,1,…,*N*−1.

These expressions are useful for numerical calculations of GS properties as well as the asymptotics of fractal dimensions when [22–24]. For illustration, in figure 4*a*, the coefficients of the GS function for the critical QIM with *λ*=1 are presented.

To find , it is necessary to know the largest coefficient in the expansion (1.2). From physical considerations, it is clear that it corresponds to the pure ferromagnetic configuration where all the spins are up. Equation (4.6) when all *σ*=1 gives
4.9When , one obtains [23,24]
4.10Using similar arguments, one concludes that is related with the minimal coefficient in (1.2). As we consider the ferromagnetic Ising model, it is clear that the minimum is attained for a configuration with the maximum number of spins down. If *N* is even, then the configuration (4.6) when all *σ*=−1 is allowed. For odd *N* there exist *N* configurations with *N*−1 spins down and one spin up. In all cases when , one concludes [23,24] that
4.11Using the Kramers–Wannier duality in [23,24], it was shown that for the QIM it is possible to calculate also *D*_{1/2}. Generalizing slightly the arguments of this work, we get
4.12These formulae prove that fractal dimensions of the QIM do exist and are non-trivial.

The above formulae are qualitatively the same for non-critical and critical (that is, *λ*=1) QIM but their sum
4.13has a singularity at *λ*=1, which means that fractal dimensions are also sensitive to quantum criticality.

### (b) Numerics

As mentioned above, the simplest method to find fractal dimensions is direct numerical calculation of the GS wave function for different number of spins and a subsequent extrapolation of Rényi entropy for large *M*. In definition (1.5), it is implicitly assumed that *q* is positive. For many problems (but not for all), fractal dimensions can also be calculated for negative *q* [25,26]. When certain coefficients in (1.2) are zero owing to an exact symmetry (as in QIM), they are not included in the calculation of Rényi entropy (1.4) for *q*≤0. The curves *S*_{R}(*q*,*M*) as a function of were fitted to the simplest functional form
4.14The coefficients *a*_{i} are, of course, functions of *q*. In practically all models, these curves (cf. (1.3)) as a function of *N* are almost straight lines and the slope *a*_{1} gives fractal dimensions with a reasonable precision. To estimate the precision of the fit, we calculate the sum
4.15where the summation is performed over all available values of matrix dimensions *M* and *P*_{q}(*M*) is defined in (1.3).

For the Ising model and all the other models, we perform the Lanczos algorithm with *n*=150 iterations. As in the QIM parity is conserved, the iteration of the Lanczos method starts with the ferromagnetic configuration with all spins up. The quality of the fit is illustrated in figure 5*a*. In figure 5*b*, the comparison between the direct diagonalization and the Lanczos method is presented. Though exact diagonalization was performed for a relatively small number of spins, the results agree well with the Lanczos results.

As an example of the precision of numerical calculations, we compare the exact values of *D*_{1/2} with our numerics. One has from (4.12)
4.16Here, *K* is the Catalan constant.

Numerical fits give
4.17We see that the accuracy of the calculations at this particular point is of the order of 10^{−3}–10^{−4}, which is enough for all practical purposes. More numerical results are presented in [6].

## 5. XY model

The Hamiltonian of this model differs from the QIM by the anisotropy *γ*,
5.1Similar to the QIM, this model is also integrable by the Jordan–Wigner transformation [4,23,24] but the structure of its GS is more complicated (see below). The eigenenergies are given by the same expression (4.2) but with
5.2where *k* is as in (4.5).

As for the Ising model, the GS wave function coefficients are calculated from (4.6) and (4.7) but with *θ*_{k} determined by the expression
5.3The XY model reduces to the Ising model for *γ*=1. In figure 4*b*, an example of the GS wave function for the XY model is presented.

### (a) Ground-state energy degeneracy

Owing to anisotropy *γ*, the parity of the GS is not fixed *a priori*. Depending on parameters *λ* and *γ*, the lowest energy state may have either odd or even parity. Correspondingly, the GS energy is given by the same formula as (4.2) (with *n*_{k}=0),
5.4but either with *k*=2*πl*/*N* or *k*=*π*(2*l*+1)/*N*, which we called even and odd momenta.

The energy difference between these two states is
5.5If this difference is positive, *E*_{GS}=*E*_{odd}; if it is negative, *E*_{GS}=*E*_{even}. To calculate the sum in (5.5) when , one can proceed as follows.

Denote *z*=e^{ik}. For even modes *z*_{l} are roots of *z*^{N}−1=0; for odd modes *z*_{l} are roots of *z*^{N}+1=0. It is easy to check that
5.6with *P*(*z*)=(*z*^{2}−2*λz*+1)^{2}−*γ*^{2}(*z*^{2}−1)^{2} and contour *C* encircles all poles on the unit circle as in figure 6.

Function *f*(*z*) has four square root singularities when *P*(*z*)=0:
5.7If *λ*^{2}+*γ*^{2}<1 these roots are complex; when *λ*^{2}+*γ*^{2}≥1 they are real (cf. figure 6). Contour *C* can be deformed to encircle the cuts of *f*(*z*).

When , one can use and *P*(*z*) close to *z*=*z*_{1} can be approximated as *P*(*z*)=*P*^{′}(*z*_{1})(*z*−*z*_{1}), where
5.8with similar formulae for the other zeros.

Performing straightforward calculations, one gets
5.9where
5.10When *λ*^{2}+*γ*^{2}<1 difference (5.9) oscillates (it has *N* zeros) and is exponentially small, *E*_{even}−*E*_{odd}∼|*r*|^{N}, where
5.11When *λ*^{2}+*γ*^{2}>1 this difference is positive and *E*_{GS}=*E*_{odd} as for the QIM.

To avoid the difficulty of finding the correct GS, we choose parameters *λ* and *γ* outside or inside the unit circle
5.12and perform calculations of fractal dimensions only in such cases.

### (b) Factorizing field

When *λ*^{2}+*γ*^{2}=1 or *λ*=*λ*_{f}, where
5.13the XY model has two degenerate GSs with different parities.

It is known [27,28] that at that field the XY model has two exact factorized GS wave functions
5.14with
5.15This GS is doubly degenerate as both *θ* and −*θ* obey equation (5.15).

For completeness, we sketch here the proof of this result. To find the factorizing field in the XY model, it is sufficient to consider one term in the Hamiltonian (5.1) (so that ), 5.16and try to satisfy the eigenvalue condition 5.17with 5.18Direct calculations give that equation (5.17) will be satisfied if and only if the parameters obey the equations 5.19These equations are easily solved, 5.20and they give equation (5.15).

As *e*=−1, the full GS energy for such a factorized state is
5.21In the factorizing state, the coefficients of the expansion (1.2) are
5.22where *n* is the number of spins up in the state . Comparing it with (2.7), one concludes that it corresponds exactly to the binomial measure discussed in §2. Therefore, fractal dimensions in this case are (cf. (2.9))
5.23To remove the double degeneracy of the GS in the factorizing field, it is convenient to form two states with different parities (4.4),
5.24where *N*_{±} are normalization constants
5.25For such states with fixed parity,
5.26and the fractal dimensions are the same as in (5.23).

To illustrate this case, we computed numerically by the Lanczos method moments of GS wave function (1.3) for the XY model in the factorizing field *γ*=0.6 and *λ*=0.8 (*λ*^{2}+*γ*^{2}=1) and compared them with exact expression (5.26). In the inset of figure 7*a*, the relative errors
5.27are presented for *q*=2,2.5,3.5. Here, (*P*_{q})_{num} are numerically calculated moments and (*P*_{q})_{exact} are exact moments given by (5.26).

The good agreement observed even at large *N* confirms the precision of numerical calculations. Fractal dimensions at these values of parameters are calculated using fit (4.14). From figure 7*b*, one can estimate the accuracy of the interpolation of such a fit for certain values of parameters in the XY models. Though the fit is quite good (i.e. the fitting points are close to the numerical values), the fractal dimensions deviate from the exact ones (5.23) by values of the order of 10^{−2} as shown in figure 7*a*. The main reason for such discrepancies is non-uniform convergence of exact formula (5.23) when in different intervals of *q*. As the form of corrections to limiting values for fractal dimensions is unknown, it is difficult, in general, to estimate the precision of calculation of fractal dimensions. From our experience, we find that it is reasonable to get the absolute error of the order of 10^{−2}−10^{−3} from the data up to 16 spins, though for certain values of *q* higher accuracy may be obtained.

### (c) Limiting values of fractal dimensions

It seems natural to think that the limiting values and as in the QIM should correspond to configurations with, respectively, all spins up and all spins down, and, consequently, be expressed similar to (4.10) and (4.11) as
5.28But it appears that, at small *λ*, there exists another distribution of spins that gives a contribution smaller than the one with all spins down. It corresponds to the anti-ferromagnetic Néel configuration with alternating spins, *σ*_{n}=(−1)^{n}, whose contribution can be calculated analytically as done below.

For simplicity, we consider even *N*. By a simple transformation from (4.6), one gets that for the Néel configuration it is necessary to calculate the determinant
5.29where the orthogonal matrix *O*_{mn} from (4.7) and (5.3) is a Toeplitz matrix, *O*_{mn}=*O*_{m−n}, i.e.
5.30Here, *k*=2*π*(*l*+1/2)/*N*.

When , one can change the summation over *k* into an integration and get
5.31Matrix *K*_{mn} in equation (5.29) is not a Toeplitz matrix but can be written as a block Toeplitz matrix,
5.32where the *Π*_{k} are 2×2 matrices,
5.33The asymptotics of the determinant of a block Toeplitz matrix under certain conditions is given by the formula [29] (see also [30])
5.34where *Φ*(*θ*) is the symbol of the block Toeplitz matrix,
5.35For matrices *Π*_{k} given by (5.33), one gets
5.36Using
5.37one concludes that
5.38Finally,
5.39where
5.40The above formulae give
5.41The imaginary part of will give zero after integration over *θ*, and the real part is
5.42where .

Consequently, the contribution of the Néel configuration to the fractal dimension is
5.43When *λ*=0
5.44from which it is clear that the Néel configuration dominates at small *λ* when *γ*>1.

In general, if *D*_{Neel}>*D*_{−}, , otherwise . For *γ*=1.4, these curves intersect at *λ*≈0.4982 and is built from two curves. Numerical results presented in [6] agree well with this prediction.

## 6. XXZ and XYZ models

In previous sections, we discussed the QIM and the XY model with ferromagnetic GS. Here, we briefly consider the XXZ and XYZ models, choosing the parameters such that their GSs are anti-ferromagnetic.

The XXZ model in zero field is a particular case of the Heisenberg model (1.1) with *γ*=*λ*=*α*=0 and *Δ*≠0,
6.1Owing to the conservation of the *z* component of the total spin, , this Hamiltonian can be diagonalized in subspace with fixed total spin along the *z*-axis. The model can be solved using the coordinate Bethe anzatz [2,31,32] and has a rich phase diagram (e.g. [5]).

For this model, there exists a special point, *Δ*=−1/2, called the combinatorial point, where more information about the GS wave function is available. From the Razumov–Stroganov conjecture [33] proved in [34], it follows that at this value of *Δ* and odd *N*=2*R*+1 the following statements are valid:

— the GS energy is −3

*N*/4;— the largest coefficient in the normalized GS expansion (1.2) (the one for the Néel configuration) is 6.2

— the sum of all terms is given by 6.3

— the minimal coefficient corresponds to half consecutive spins up and other spins down 6.4where

*A*_{R}is the number of alternating sign*R*×*R*matrices given by the formula [35] 6.5

These formulae imply that, for *Δ*=−1/2, fractal dimensions and *D*_{1/2} can be calculated analytically as
6.6The value of the minimal coefficient determines the limiting behaviour of Rényi entropy (1.4) at large negative *q*. The asymptotics of *A*_{R} when is: . Hence, the smallest coefficient in the XXZ model with *Δ*=−1/2 decreases exponentially not with *N* but with *N*^{2}. This quadratic behaviour is a particular case of the emptiness formation probability of a string of *n* aligned spins with *n*∼*N* [36,37]. Such asymptotic behaviour means that negative moments of the GS wave function in the anti-ferromagnetic case require a scaling different from (1.5) and will not be considered here. In figure 8*a*, the wave function of the XXZ model at the combinatorial point is presented.

The numerical calculations were performed by extrapolation of the Rényi entropy separately for odd and even *N*=3 to 19 (cf. figure 9). The necessity of different formulae for odd and even *N* is physically related to the existence of different total spin projection in the GS. For even *N*, the total GS projection along the *z*-axis is zero, but for odd *N* it is .

To find fractal dimensions, we use three different kinds of fits to the available data with the same number of unknowns:
6.7The point is that the exact asymptotics of contains a logarithmic term, which corresponds to the term *a*_{3} in *f*_{2}(*z*) and *f*_{3}(*z*). From figure 9, it is clear that the second and third fits give better approximation to numerical data than the first one. With the available precision, the fractal dimensions for odd and even *N* are the same but subleading terms in Rényi entropy (1.4) are different. The numerical calculation of *D*_{1/2} gives *D*_{1/2}=0.7949, which agrees well with exact result (6.6).

The XYZ model is similar to the XXZ model but with the anisotropy *γ* present,
6.8Its GS wave function can be found by the algebraic Bethe anzatz [3].

This model also has a special value of the field, 6.9where the GS wave function has the factorized form (5.14) with 6.10As for the XY model at this field, all fractal dimensions are the same as for the binomial measure (5.23).

The combinatorial point for the XYZ model at zero field is *Δ*=(*γ*^{2}−1)/2. At this point, certain properties of the GS are known (or conjectured) [38] but they do not permit fractal dimensions to be found analytically.

As an example, we present in figure 8*b* the GS wave function for the XYZ model with *Δ*=−0.5, *γ*=0.6 and *N*=12. Qualitatively, fractal dimensions in the XYZ model are similar to the XXZ model and are presented in [6].

## 7. Conclusion

Wave functions for *N*-spin problems can be represented only as a collection of an exponentially large (as ) number of coefficients corresponding to all possible basis states. This proliferation makes difficult not only the analysis but even the representation of such functions. There exists no ‘natural’ ordering of wave function coefficients and simple attempts similar to the use of binary codes as in the examples above lead to irregular and complex structures even in one-dimensional models.

The multi-fractal formalism has been developed to measure quantitatively how irregular different objects are. The main result of this paper and of [6] is the demonstration that the multi-fractal formalism is an adequate and useful language to describe GS wave functions of spin chains. This formalism can informally be considered as an analogue of the usual thermodynamic formalism but applied not to true energies of interacting particles but to ‘pseudo-energies’ associated with wave function coefficients (cf. (2.11)).

We stress that the models considered have no random parameters. The origin of their multi-fractality is not an interplay between the localization in random relief and the spreading as in critical models such as the three-dimensional Anderson model at the metal–insulator point but the complexity of the internal structure of many-body Hilbert space. In this respect, the simple example of the binomial measure is a characteristic one. It clearly shows how the necessity of choosing one path from an exponentially large number of other possibilities results in irregular multi-fractal structure.

The central object of the multi-fractal formalism is fractal dimensions, *D*_{q}, which play the role similar to the free energy in the usual thermodynamics with *q* being the inverse temperature. It is well established that, for classical Hamiltonians with local interactions, the free energy exists in thermodynamic limit . Much less is known rigorously about fractal dimensions. Our results strongly suggest that for quantum *N*-body Hamiltonians with local interactions they do exist (i.e. the limit gives a well-defined answer) but we are unaware of formal proof even in simple models like QIM.

Only in very rare models (like spin chains in a factorizing field) can all fractal dimensions be calculated analytically. In some models, one can find an exact expression for fractal dimensions at certain values of *q*. Here, we present details of the calculations of limiting values and *D*_{1/2} in QIM, the XY model and XXZ model at *Δ*=−1/2. These results are important as they prove the existence of fractal dimensions at least for special values of *q* and can be used to control the precision of numerical calculations.

In general, fractal dimensions have to be estimated numerically and the main question is how carefully the limit can be found from the numerical data with *N* of the order of 10–20. In all models discussed in this paper, it appears that the corrections to the limit are reasonably small but the form of interpolation may depend on *q*. In many cases, the precision of 10^{−2}–10^{−3} can be achieved without considerable numerical efforts.

Though in this paper only spin chain models are considered, preliminary results on one-dimensional bosonic and fermionic models demonstrate that their GS wave functions are also multi-fractals. Based on these calculations, we conjecture that the multi-fractality of the GS wave function is a common feature of generic *N*-body quantum Hamiltonians with local interactions. The discussion of multi-fractal properties of excited states is postponed to a future publication.

## Funding statement

Y.Y.A is supported by the CFM Foundation.

## Acknowledgements

The authors are greatly indebted to G. Roux for many useful discussions and especially for permission to use his code to implement the Lanczos algorithm. Numerous conversations with O. Giraud are much appreciated.

## Footnotes

One contribution of 13 to a Theo Murphy Meeting Issue ‘Complex patterns in wave functions: drums, graphs and disorder’.

- © 2013 The Author(s) Published by the Royal Society. All rights reserved.