A classification of spin 1/2 matrix product states with two dimensional auxiliary matrices
aa r X i v : . [ qu a n t - ph ] J u l A classification of spin 1/2 matrix product states withtwo dimensional auxiliary matrices
Marzieh Asoudeh † Department of Physics, Shahid Beheshti University, GC19839-63113, Tehran, Iran
Abstract
We classify the matrix product states having only spin-flip and parity sym-metries, which can be constructed from two dimensional auxiliary matrices. Weshow that there are three distinct classes of such states and in each case, wedetermine the parent Hamiltonian and the points of possible quantum phasetransitions. For two of the models, the interactions are three-body and for onethe interaction is two-body.
PACS Number:75.10.Jm m − asoudeh @ sbu.ac.ir Introduction
The problem of determining the ground state of a given many-body Hamiltonian,is an important problem in condensed matter and mathematical physics. There isalready a rich literature on this subject, which dates back to the work of Hans Betheon the Heisenberg spin chain and has continued since then with the works of manyother people including Yang, Baxter, and Lieb, to name only a few. In particularfor spin systems, the exponential increase in the dimension of Hilbert space of sucha system, as the number of particles rise, turns this problem into a computationallyformidable one, beyond the capability of any classical computer. In fact it has nowbeen established that finding the ground state of a given many-body Hamiltonian isthe analog of NP-complete problems for quantum computers [1]. The lesson that welearn from all this is that it is highly improbable that we be able to find a genericsystem with exactly known ground state. Nevertheless, there are systems with ex-actly known ground states and even if such systems are not exactly what we have innature or in the laboratories, they may be good approximations to real systems, orat least may teach us useful and important concepts and methods for studying morerealistic systems.One of the methods, developed in recent years, for investigating this problem isthe Matrix Product State or Finitely Correlated State [2, 3, 4] formalism. It is alsocalled Optimal Ground State formalism in some references [5, 6, 7]. The main themeis that one starts with a state with prescribed symmetries and properties, and thenconstruct the family of Hamiltonians for which this state is an exact ground state.It is obvious that for any given state | ψ i , the equation H | ψ i = 0 has always manysolutions for the unknown H , since the number of equations is much less than thenumber of unknowns. However the problem becomes interesting and quite non-trivialwhen we put physical constraints on the Hamiltonian. That is we demand that i) H be positive, so that | Ψ i is actually the ground state and not an ordinary eigen-state, ii) that it be a sum of local terms, i.e. H = P k h k,k + l , where h k,k + l acts on ablock of l + 1 spins, and iii) that both the state and the Hamiltonian have some rea-sonable physical symmetries, like parity, spin-flip, and at times rotational symmetries.In the past few years, a lot of interest has been attracted to the subject of matrixproduct states [8, 9, 10, 11, 12, 13], specially after the emergence of the field of quan-tum information [14, 15, 16, 17, 18]. The reason is the complementary role that thefields of condensed matter physics and quantum information play in investigation ofmany body systems. On the one hand quantum information starts with propertiesof states, while condensed matter physics, starts from the properties of the Hamilto-nian which embodies the interactions and energy of the system. The matrix productformalism is one of the subjects which lies at the borderline of these two subjects.As is well known, in this formalism, one starts from proposed states whose ex-pansion coefficients are the trace of product of given matrices. While for numericalinvestigations, i.e. the density matrix renormalization group (DMRG), one usually1tarts from large dimensional matrices, to simulate ground states of given Hamilto-nians, in the approach which is used for finding exactly solvable models, one startsfrom low dimensional matrices and finds Hamiltonians for which these states are exactground states. This is the approach which has been used in our works and in manyother works in the past few years [8, 9, 10, 11, 12, 13, 19, 20, 21, 22, 23, 24, 25]. Inthis article we follow this approach and classify all the matrix product states whichcan be constructed from two dimensional matrices. We restrict ourselves to stateswhich allow one or another of the spin-flip or parity symmetries and find that thereare three classes of such matrix product states. We will study these states and findthe parent Hamiltonians and also the points or lines in the space of control parame-ters where a MPS-quantum phase transition [26] (MPS-QPT) may occur. This is aterm, introduced in [27] to differentiate these kinds of QPT’s (characterized by anydiscontinuity in any macroscopic quantity) from the conventional QPT,s in which anon-analyticity in the ground state energy typically occurs.The structure of this paper is as follows: To make the article self-contained, in thenext section we briefly introduce the basic elements of the formalism. In section (3)we discuss the symmetry properties of MPS and classify the spin 1/2 states with twodimensional matrices. We show that there are three classes, denoted by a model A,model B and model C and studied in subsequent sections. We end the paper with adiscussion. First let us review the basics of matrix product states. Consider a homogeneous ringof N sites, where each site describes a d − level state. The Hilbert space of each siteis spanned by the basis vectors | i i , i = 0 , · · · d −
1. A state | Ψ i = X i ,i , ··· i N ψ i i ··· i N | i , i , · · · , i N i (1)is called a matrix product state if there exists D dimensional complex matrices A i ∈ C D × D , i = 0 · · · d − ψ i ,i , ··· i N = 1 √ Z tr ( A i A i · · · A i N ) , (2)where Z is a normalization constant given by Z = tr ( E N ) (3)and E := d − X i =0 A ∗ i ⊗ A i . (4)Here we are restricting ourselves to translationally invariant states, by taking thematrices to be site-independent. By defining the vector valued matrix A = d X i =1 A i | i i , (5)2ne can write the MPS in a more concise way as | ψ i = tr ( A ⊗ A ⊗ · · · A ) , (6)where we use the convention tr ( A ⊗ A ) := tr ( A i A j ) | i i ⊗ | j i . It is important to note that the MPS representation (2) is not unique and atransformation such as A i −→ µU A i U − (7)where U is an invertible matrix, and µ is a constant, leaves the state invariant. Thesimple structure of the MPS allows also an easy calculation of correlation functions.Let O be any local operator acting on a single site. Then we can obtain the one-pointfunction on site k of the chain h Ψ | O ( k ) | Ψ i as follows: h Ψ | O ( k ) | Ψ i = tr ( E k − E O E N − k ) tr ( E N ) , (8)where E O := d − X i,j =0 h i | O | j i A ∗ i ⊗ A j . (9)In the thermodynamic limit N −→ ∞ , equation (8) gives h Ψ | O | Ψ i = h λ max | E O | λ max i λ max , (10)where we have used the translation invariance of the model and λ max is the eigenvalueof E with the largest absolute value and | λ max i and h λ max | are the right and lefteigenvectors corresponding to this eigenvalue, normalized such that h λ max | λ max i = 1.Here we are assuming that the largest eigenvalue of E is non-degenerate. In case λ max is degenerate with degree equal to g , then Eq. (10) will be modified to h Ψ | O | Ψ i = P gi =1 h λ max,i | E O | λ max,i i λ max , (11)The n-point functions can be obtained in a similar way. For example, the two-pointfunction h Ψ | O ( k ) O ( l ) | Ψ i can be obtained as h Ψ | O ( k ) O ( l ) | Ψ i = tr ( E O ( k ) E O ( l ) E N ) tr ( E N ) (12)where E O ( k ) := E k − E O E − k . Note that this is a formal notation which allows usto write the n-point functions in a uniform way, it does not require that E be aninvertible matrix. In the thermodynamic limit the two point function turns out to be h Ψ | O (1) O ( r ) | Ψ i = 1 λ rmax X i λ r − i h λ max | E O | λ i ih λ i | E O | λ max i . (13)3or large distances r ≫
1, this formula reduces to h Ψ | O (1) O ( r ) | Ψ i − h Ψ | O | Ψ i = λ r − λ rmax h λ max | E O | λ ih λ | E O | λ max i , (14)where λ is the second largest eigenvalue of E for which the matrix element h λ | E O | λ max i is non-vanishing and we have assumed that the eigenvectors of E have been normal-ized, i.e. h λ i | λ j i = δ ij . Thus the correlation length is given by ξ = 1ln λ max λ . (15)Any level crossing in the largest eigenvalue of the matrix E signals a possible MPS-QPT . Here we are using the term quantum phase transition in a broader sense thanusual [27], that is, we call any discontinuity in any macroscopic quantity a quantumphase transition, even if the ground state energy itself is a continuous function of thecoupling constants. Also, due to (15), any level crossing in the second largest eigen-value of E implies the correlation length of the system has undergone a discontinuouschange. Given a matrix product state, the reduced density matrix of k consecutive sites isgiven by ρ i ··· i k ,j ··· j k = tr (( A ∗ i · · · A ∗ i k ⊗ A j · · · A j k ) E N − k ) tr ( E N ) . (16)The null space of this reduced density matrix includes the solutions of the followingsystem of equations d − X j , ··· ,j k =0 c j ··· j k A j · · · A j k = 0 . (17)Given that the matrices A i are of size D × D , there are D equations with d k unknowns.Since there can be at most D independent equations, there are at least d k − D solutions for this system of equations. Thus for the density matrix of k sites to havea null space it is sufficient that the following inequality holds d k > D . (18)Let the null space of the reduced density matrix be spanned by the orthogonal vectors | e α i , ( α = 1 , · · · s, ≥ d k − D ). Then we can construct the local hamiltonian actingon k consecutive sites as h := s X α =1 µ α | e α ih e α | , (19)j where µ α ’s are positive constants. These parameters together with the parametersof the vectors | e i i inherited from those of the original matrices A i , determine the totalnumber of coupling constants of the Hamiltonian. If we call the embedding of this4ocal Hamiltonian into the sites l to l + k by h l,l + k then the full Hamiltonian on thechain is written as H = N X l =1 h l,l + k . (20)The state | Ψ i is then a ground state of this hamiltonian with vanishing energy. Thereason is as follows: h Ψ | H | Ψ i = tr ( H | Ψ ih Ψ | ) = N X l =1 tr ( h l,l + k ρ l,l + k ) = 0 , (21)where ρ l,k + l is the reduced density matrix of sites l to l + k and in the last line we haveused the fact that h is constructed from the null eigenvectors of ρ for k consecutivesites. Given that H is a positive operator, this proves the assertion.In view of the above introduction, we have a clear recipe for constructing matrixproduct states and a family of parent Hamiltonians. First one chooses the matricesthrowing away all spurious degrees of freedom by transformations (7) and reducingfurther the degrees of freedom by imposing symmetries. In this way one ends witha reasonable set of matrix product states, which hopefully may have applications indescription of real physical systems. Imposing a continuous symmetry, like rotationaround an axis, restricts the matrices considerably [8, 9, 10]. In this article we re-strict ourselves to discrete symmetries only which allow a larger variety of modelsto be constructed. For two dimensional auxiliary matrices, this is a simple tractableproblem, which we do in this article. For larger matrices, the problem is not so simpleand we defer it to another work. We now classify all the two dimensional matrices which can be used for constructingspin 1/2 matrix product states. We restrict ourselves to the case where these stateshave spin-flip and left-right symmetries.
Consider now a local symmetry operator R acting on a site as R | i i = R ji | j i wheresummation convention is being used. R is a d dimensional unitary representation ofthe symmetry. A global symmetry operator R := R ⊗ N will then change this state toanother matrix product stateΨ i i ··· i N −→ Ψ ′ := tr ( A ′ i A ′ i · · · A ′ i N ) , (22)where A ′ i := R ij A j . (23)5 sufficient but not necessary condition for the state | Ψ i to be invariant under thissymmetry is that there exist an operator U ( R ) such that R ij A j = U ( R ) A i U − ( R ) . (24)Thus R and U ( R ) are two unitary representations of the symmetry, respectively ofdimensions d and D . Equation (24) will be our guiding lines in defining states withprescribed symmetries. Spin-flip symmetry means that ψ i ,i , ··· i N = ψ i N ,i N − , ··· i , (25)where i = 1 − i. For a matrix product state, this requires that there be a matrix like X , such that XA X − = ǫA , XA X − = ǫA , (26)where ǫ = ±
1. Similarly left-right symmetry means that ψ i ,i , ··· i N = ψ i N ,i N − , ··· i . (27)For a matrix product state, this means that there be a matrix Ω such thatΩ A Ω − = σA T , Ω A Ω − = σA T , (28)where the superscript T stands for the transpose and σ = ±
1. These conditionsare general irrespective of the dimension of matrices. For two dimensional matriceshowever, if we take the trace and determinants of both sides of equations (26) and(28), and comparing them, we find that tr ( A ) = ǫtr ( A ) , (29)and det( A ) = det( A ) . (30)The important point is that for two dimensional matrices, the trace and determi-nant are the only invariants under similarity transformations, and hence these twoequations allow us to classify all the matrices A and A which can be used for con-struction of spin 1/2 matrix product states. We will use the freedom (7) and also theabove two conditions to show that there are three distinct classes of matrix pairs andcorresponding matrix product states. In the next three sections, we introduce thematrix pairs and study the properties of the matrices obtained from them. For easeof distinction, we use different notations for the matrix pairs in each section, namelywe denote the matrix pairs by A i , B i and C i in the following sections. If one of the matrices say A is diagonalizable, we can use freedom in re-scaling A −→ µA to put it in the form A = g
00 1 − g ! , (31)6here g is a free parameter. Now take A as an arbitrary matrix of the form A = a bc d ! . If bc = 0, then we can use the transformation (7) with U = √ c √ b ! and a further re-definition of √ bc −→ b to put it in the form A = a bb d ! . (32)For this type of matrices, the MPS is automatically left-right symmetric, with Ω = I .From equations (29) and (30), we find the following constraints on the parameters,1 − g = ( ad − b ) 2 ǫ = ( a + d ) . (33)To solve the second equation, we put a = ǫ + u, b = ǫ − u, (34)which turns the first equation into g = b + u , (35)which can be solved by the parametrization u = g cos θ and b = g sin θ . Therefore thefinal form of the matrices will be as follows: A = g
00 1 − g ! , A = ǫ + g cos θ g sin θg sin θ ǫ − g cos θ ! . (36)The matrices satisfy the symmetry constraints (26) and (28) with Ω = I and X = ǫ sin θ − ǫ sin θ − ǫ cos θ ǫ cos θ ! . (37)We now restrict ourselves to the case ǫ = 1 (We do this also for other models Band C). To study the properties of the corresponding MPS, we should determine theeigenvalues of the transition matrix E A := A ⊗ A + A ⊗ A . For general values ofthe parameters, the analytical form of these eigenvalues are complicated. However wecan gain insight by looking at them for for generic values of the parameter g . Figure(4) shows the eigenvalues as a function of θ for two values of the parameter g . Thesame pattern repeats for other values of g . We first see that there is no level-crossingin the largest eigenvalue and hence no MPS quantum phase transition in this model.However for every value of g , two points are important. The point θ = Π (or − Π), isthe point where the largest and the next to largest eigenvalues become equal. Thisis a point where according to (15), the correlation length becomes infinite. It is seenfrom (36) that at this point, the matrices become diagonal A = diagonal (1 + g, − g )and A = diagonal (1 − g, g ) and according to (6), the un-normalized state becomesthe sum of two product states, namely | ψ i = | φ + i ⊗ N + | φ − i ⊗ N , (38)7igure 1: (Color Online) The eigenvalues of the transition matrix E A for ǫ = 1 and ageneric value of g (i.e. g=1/2) as a function of θ .where | φ + i = (1 + g ) | i + (1 − g ) | i , | φ − i = (1 − g ) | i + (1 + g ) | i . (39)Let us now find the parent Hamiltonian for a fixed value of the parameter g , say g = 1. To this end, we have to see for which value of k (the range of interaction), thesystem of equations (17) have a non-trivial solution. It is seen that the smallest k forwhich there is such a solution is k = 3.The solution space of this system of equations turn out to be spanned by thefollowing vectors: | e A i = − θ )2 | i + | i , (40) | e A i = | i − | i (41) | e A i = | i − | i (42) | e A i = | i − cos ( θ )2 | i . (43)The symmetries of the parent Hamiltonian now show itself in the form that theabove states, which transform into each other under the action of these symmetries.To have a Hamiltonian which respects these symmetries, we construct it as follows:8 A = J ( | e A ih e A | + | e A ih e A | ) + K ( | e A ih e A | + | e A ih e A | ) . (44)The final form of the the full Hamiltonian in terms of Pauli matrices, after neglectingadditive and multiplicative constants becomes H A = N X i =1 J σ zi σ zi +1 + J σ zi σ zi +2 − uJ σ xi σ xi +2 + uJ σ yi σ yi +2 − K σ xi + K σ zi σ xi +1 σ zi +2 , (45)where u = θ and J = J u − , (46) J = J u + 12 − K . (47)This is a three-body Hamiltonian with two free coupling constants. We will say moreabout this in the discussion. In accordance with our proposed notation, we denote the matrices in this case by B and B . In this case the matrices are the same as in the previous case, except thatone of the parameters, say b is zero. There is no transformation which can put B into symmetric form, and we have B = a c d ! . (48)In this case no similarity transformation can put the matrix B into symmetric form,without destroying the diagonal form of B and hence the MPS will not be parityinvariant or left-right symmetric. From Eqs. (29) and (30), we find that1 − g = ad ǫ = a + d, (49)which lead to the following final parametrization for the matrices, where for definite-ness we will show the matrices by a different letter B = g
00 1 − g ! , B = ǫ + g c ǫ − g ! . (50)In this case again we have two free parameters in the MPS state, namely c and g .Moreover the matrices satisfy the spin flip symmetry condition (26) with9 = 1 + ǫ g c ! + 1 − ǫ g c ! , (51)while they do not have any symmetry under parity (i.e. there is no matrix Ω satisfy-ing (28)).The eigenvalues of the transition matrix E B = B ⊗ B + B ⊗ B , (for the case ǫ = 1) now become λ B = 2(1 + g ) , λ B = 2(1 − g ) , λ B , = 2(1 − g ) , (52)independent of the value of c . Figure (5), shows these eigenvalues as a functionof the parameter g . From this figure, a few features can be recognized. First wenote that at g = 0, there is a crossover between the largest and the second-largesteigenvalues. This points to a possible MPS-quantum phase transition at this point.Furthermore at g = ±
1, there is a discontinuity in the derivative of the second largesteigenvalue which points to a discontinuity in the derivative of the correlation length.In view of the general and rather broad definition of MPS quantum phase transition,as the appearance of any discontinuity of a macroscopic observable, these points arealso points of MPS-QPT’s. Moreover since the eigenvalues do not depend on theparameter c , it appears that the above points are really lines in the space of couplingconstants c and g .From (36), it is clear that at g = 0, the un-normalized MPS turns into the followingstate | ψ B i = | χ + i ⊗ N + | χ − i ⊗ N , (53)where | χ ± i = (1 ± g ) | i + | i . (54)Finally we come to the parent Hamiltonian. For this model we find that thesmallest value of k for which the system of equations (17) has a non-trivial solutionis k = 2 and hence we can have a two-local parent Hamiltonian. The solution spaceof the system of equations (17) is spanned by the vectors | e B i = −
12 (1 + g ) | i + | i + 12 ( g − | i (55) | e B i = −
12 (1 + g ) | i + | i + 12 ( g − | i . (56)Interestingly we note that the above vectors transform into each other under spin-flip, but they do not have any transformation property under parity which is tobe expected, since the original matrices had only spin-flip symmetry. The parentHamiltonian which is symmetric under spin flip will be10igure 2: (Color Online) The absolute values of the eigenvalues of the transitionmatrix E B for ǫ = 1, as a function of g . h B = J ( | e B ih e B | + | e B ih e B | ) (57)and the full Hamiltonian will be (after neglecting additive and multiplicative constantsand collecting all the various terms) H B = N X i =1 (1 − g )( σ xi σ xi +1 − σ yi σ yi +1 ) + 1 + 2 g σ zi σ zi +1 + σ xi . (58)This is the Heisenberg XYZ system with specific couplings, that is we have foundexact solution on a submanifold in the space of couplings. Note that since the Hamil-tonian does not depend on the paramter c , while the MPS does, this means thatthere is a large degeneracy in the ground state. Expansion of the MPS in terms ofthe parameter c , i.e. | Ψ( c, g ) i = P n c n | ψ n ( g ) i , will yield the multitude of degenerateground states | ψ n ( g ) i as in [28]. Denoting the matrices by C and C , this is the only remaining case, where C is notdiagonalizable and can be put only in the Jordan form C = g ! . (59)11e take the general form of C to be C = a bc d ! and impose the conditions (29)and (30), from which we obtain the constraints a + d = 2 ǫ (60) ad − bc = 1 . (61)The first constraint is solved by the parametrization a = ǫ + u and d = ǫ − u , whichwhen inserted into the second equation, gives u + bc = 0, the solution of whichis b = µu and c = − uµ . However the parameter µ can be set to unity by a gaugetransformation C i −→ U C i U − with U = µ ! . Thus the final form of thematrices become C = g ! , C = ǫ + u u − u ǫ − u ! . (62)One can verify the existence of both symmetries, with matrices X = u ug ! and Ω = − ! , (63)such that X − C i X = C i , Ω − C i Ω = C Ti . (64)The eigenvalues of the transfer matrix E C = C ⊗ C + C ⊗ C , are found to be λ C = 2 , λ C = 2 − ug, λ C ± = 2 + ug ± q ug + u g . (65)Figure (6) shows the largest eigenvalue of the transfer matrix in the u − g plane. Itis seen that the lines u = 0 and g = 0 are the crossover lines where the largest eigen-value changes and hence an MPS-QPT (MPS quantum phase transition) is expectedto occur on these lines.Finally we can find the parent Hamiltonian of this model. The system of equations(17) has a nontrivial solution for k = 3 and the solution space is spanned by thevectors: | e C i = | i + | i − | i − | i| e C i = (1 + ug )( | i + | i ) + ( | i + | i + | i + | i ) − | i + | i ) | e C i = 2( | i − | i ) − | i − | i + | i − | i )12igure 3: (Color Online) The largest eigenvalue of transfer matrix for model III. Thereare two crossover lines for the largest magnitude eigenvalue. | e C i = 2( | i − | i ) − (1 + ug )( | i − | i + | i − | i ) (66)The interesting point about these vectors is that all of them are invariant (mod-ulo a sign) under the parity and spin-flip transformations and hence the symmetricHamiltonian can be written in the following form, with four free coupling constants: h C = X i =1 J i | e Ci ih e Ci | . (67)The explicit expression of the total Hamiltonian can be obtained along the samelines as for model A and model B. We will not do it here. In the formalism of matrix product states, there is a large room for constructingstates and parent Hamiltonians. What really constrains this freedom and guides usalong a way which may lead to interesting states and Hamiltonians is considerationof symmetries. The other constraining elements is the dimension of matrices whichwe choose. In this article we have classified all such states which are constructed fromtwo dimensional matrices and have two important symmetries, namely the spin flipsymmetry and the parity symmetry. We have shown that there are three differentmodels, two of which lead to parent Hamiltonians with nearest and next-nearestinteraction (models A and C) while one of them lead to a Hamiltonian with nearestneighbor interaction (model B). Furthermore, by calculating the eigenvalues of thetransfer matrix in each case and determining the points of crossover between thelargest and the next to largest eigenvalues of this matrix, we have identified the pointsof possible MPS quantum phase transitions. While in many of the works which havebeen reported on model building in matrix product states, [8, 9, 10, 19, 20, 21, 22, 23],13otation symmetry has been taken into account, a condition which highly restrictsthe form of matrices, in this article we have relaxed this continuous symmetry inorder to find all models compatible with discrete symmetries in order exhaust all thepossibilities with two dimensional matrices. Any model with these symmetries mustbe equivalent to one of the above three models. For example in [27], the followingmodel was suggested A ′ = ! , A ′ = q ! , (68)which has spin-flip symmetry with X = q ! . The Hamiltonian for this modelis [27] H = N X i =1 q − σ zi σ zi +1 − (1 + q ) σ xi + ( q − σ zi σ xi +1 σ zi +2 . (69)It is easy to see that this model is equivalent to model A above. In fact the transfor-mation SA i ( g = 1) S − = 2 A ′ i with S = θ sin θ − cos θ ! in which q = θ , proves the equivalence. I would like to thank V. Karimipour for valuable discussion and comments.
References [1] M. A. Nielsen, and I. L. Chuang;
Quantum computation and quantum informa-tion ,Cambridge University Press, Cambridge, 2000.[2] M. Fannes, B. Nachtergaele and R. F. Werner, Commun. Math. Phys. , 443(1992).[3] I. Affleck, T. Kennedy, E. H. Lieb, H. Tasaki, Commun.Math. Phys. , 477(1988);[4] A. Klumper, A. Schadschneider and J. Zittartz, J. Phys. A (1991) L293; Z. Phys.B, (1992) 281; Europhys. Lett., (1993) 293.[5] H. Niggemann, and J. Zittartz, J. Phys. A: Math. Gen. 31, 9819-9828 (1998).[6] E. Bartel, A. Schadschneider and J. Zittartz, Eur. Phys. Jour. B, , 2, 209-216(2003). 147] M. A. Ahrens, A. Schadschneider, and J. Zittartz, Europhys. Lett.