Monte Carlo simulations of random non-commutative geometries
MMonte Carlo simulations of randomnon-commutative geometries
John W. Barrett, Lisa GlaserSchool of Mathematical SciencesUniversity of NottinghamUniversity ParkNottingham NG7 2RD, UK [email protected], [email protected]
May 11th, 2016
Abstract
Random non-commutative geometries are introduced by integrat-ing over the space of Dirac operators that form a spectral triple witha fixed algebra and Hilbert space. The cases with the simplest typesof Clifford algebra are investigated using Monte Carlo simulations tocompute the integrals. Various qualitatively different types of be-haviour of these random Dirac operators are exhibited. Some featuresare explained in terms of the theory of random matrices but otherphenomena remain mysterious. Some of the models with a quarticaction of symmetry-breaking type display a phase transition. Close tothe phase transition the spectrum of a typical Dirac operator showsmanifold-like behaviour for the eigenvalues below a cut-off scale.
A spectral triple is a way of encoding a geometry using a Dirac operator [10].There is a Dirac operator D acting on a Hilbert space H and an algebra A that acts on the same space. Examples with a commutative algebra are given1 a r X i v : . [ g r- q c ] M a y y Riemannian manifolds, where the algebra is the algebra of functions onthe manifold and the Dirac operator is the usual one acting on spinor fields.However, the point of spectral triples is that the algebra is allowed to benon-commutative, leading to a generalisation of the notion of geometry.A random geometry is a class of geometries G that fluctuates accordingto a probability measure. In this article, the probability measure is taken tobe a constant times e − S ( D ) d D (1)using a real-valued ‘action’ S ( D ) and a standard measure d D on the space G of Dirac operators.To make this computable, the class of geometries is taken to be the Diracoperators on a fixed finite-dimensional Hilbert space H ; thus G is a space ofmatrices. It turns out that the axioms for D for these finite spectral triplesare all linear and so G is a vector space [5]. Therefore one can take d D to beits Lebesgue measure, which is unique up to an overall constant. Thus theobject of study is a random matrix model where the matrices are constrainedto be Dirac operators.The algebra A in this construction is also fixed, and is taken to be thealgebra of n × n matrices, M ( n, C ). Spectral triples with this algebra areknown as fuzzy spaces [19] and are the simplest type of non-commutativespectral triple. Allowing the algebra to be non-commutative is importantbecause it allows a new type of finite-dimensional approximation to a man-ifold. Staying within the realm of commutative algebras would lead to thealgebra of functions on a finite set of points, which is a lattice approximationto a manifold; simple examples of such random commutative spectral triplesare studied in [18, 12]. The fuzzy spaces are not lattice approximations andso the study of these is complementary to the study of random lattices. In-tuitively one can think of the algebra A as consisting of the functions ona space with a certain minimum wavelength that is determined by n ; thispicture is known to be accurate for the most-studied example of the fuzzytwo-sphere [13].The purpose of this paper is to study the simplest examples of randomfuzzy spaces by computing the statistics of the eigenvalues of D using aMarkov Chain Monte Carlo algorithm. The examples are determined bythe value of n and the type of gamma matrices used in the Dirac operator(explicit formulas are given in section 2.1). A type ( p, q ) geometry is one inwhich there are p gamma matrices that square to 1 and q gamma matricesthat square to −
1. The examples studied here are types (1 , , , , ,
2) and (0 , S ( D ) has not yet been specified. In this paper itis assumed to be spectral, which means it is of the form S ( D ) = Tr V ( D ) = (cid:88) i V ( λ i ) (2)for some potential function V , with V ≥ b for some b ∈ R and λ i denotingthe eigenvalues of the self-adjoint operator D . The Connes-Chamseddinespectral action [7, 8], in which V ( x ) → x → ∞ , is not suitable becausethe integral for the partition function Z = (cid:90) G e − S ( D ) d D (3)does not converge. This is because S ( µD ), for µ ∈ R , converges to a finiteconstant as µ → ∞ .In fact, it is necessary that V ( x ) → ∞ instead. The simplest cases areinvestigated here, namely S ( D ) = Tr (cid:0) g D + g D (cid:1) (4)with g >
0, or g = 0 and g >
0. By a simple change of variables in theintegral, D → µD , one can assume that either g = 1, or g = 0 and g = 1,so one need only study these cases. Note that one could choose other spectralactions and it is possible that the results obtained here might motivate thestudy of other choices.When g > g < m the eigenvalue distribution, ordensity of states, for the Dirac operator is approximately the density of flatspace ρ ( λ ) = | λ | m − (5)when the eigenvalues are large enough. Most of the plots for random non-commutative geometries presented here look nothing like this, except thatsome of them are approximately constant ( m = 1) for some range of eigen-values. 3he exception is close to the phase transition, where the distributiondoes indeed look very much like (5) for a range of eigenvalues below a high-energy cutoff. Thus as far as the eigenvalues are concerned, the random non-commutative geometries are behaving something like random Riemanniangeometries in this regime. The results presented here show that this is apromising area for future investigation. A phase transition to a geometricphase in a multi-matrix model with a rather different Yang-Mills-type actionhas also been observed in [3, 23].The motivation for the present work comes from the close relation betweenrandom geometries and models of quantum gravity, though one does not needto know anything about quantum gravity to understand the results. Mostapproaches to random geometry have been stimulated by work on quantumgravity but some of them (e.g. dynamical triangulations or Liouville gravity)have found wider application. It is possible that random Dirac operators willalso find other applications beside quantum gravity. In quantum gravitythe maximum eigenvalue has a ready interpretation as a natural cutoff togravitational phenomena at the Planck scale. However in a wider contextit can be interpreted simply as a finite limit to the resolution to which ageometry is defined.There are features in common with other models of quantum or randomgeometry, most notably the existence of a phase transition, which is alsoevident in dynamical triangulations [2] and lattice simulations [15]. Thereare however some features of our system that are quite different from those ofother models. One point is that the requirement that the action has to havea compact global minimum in the non-compact G is a non-trivial constrainton the model. In other theories of discrete geometry, like causal dynamicaltriangulations [1] or causal set theory [17], the space of geometries exploredin Monte Carlo is a combinatorial space of a finite number of elements andthe code is guaranteed to reach equilibrium after a finite, although possiblylong, time.Another interesting feature is the freedom to rescale g and g using thechange of variables mentioned above. Monte Carlo simulations show thatthe rescaling does not change the qualitative features or our system, e.g.relative differences between eigenvalues will remain unchanged. This can beexplained by the use of the Lebesgue measure, which does not distinguishany particular scale of energy. This is in marked contrast to systems such asthe Ising model, or the causal set model of 2d orders [25].The non-commutativity also distinguishes this approach from most oth-ers. Using a finite-dimensional commutative algebra necessarily leads to alattice model of quantum geometry defined on a finite set of points. The useof non-commutative geometry allows a more general set of finite-dimensional4odels where the algebra is an algebra of matrices. Thus one can constructperfectly computable models of random geometry that are not lattice mod-els. Moreover, the standard model of particle physics has a non-commutativegeometry using exactly the same framework [4, 11], so the hope is it will beeasy to combine the two into a unified model of gravity and particle physics.The technical details of the Dirac operators, observable functions andMonte Carlo method are given in section 2. The results for the action Tr D are given in section 3, where it is explained how the results relate to thestandard theory of Gaussian random matrices. Actions including a Tr D term are studied in section 4, with particular attention paid to the symmetry-breaking case which exhibits a phase transition. Section 5 discusses theinterpretation of the results. The expansions of the action in terms of theconstituent matrices of the Dirac operators are given in detail in appendixA. The spectral triples considered here are ‘real spectral triples’, which consistof a finite-dimensional Hilbert space H together with some operators actingin H . These are an algebra A , a chirality operator Γ, an antilinear ‘realstructure’ J and a self-adjoint Dirac operator D . For a given random geom-etry model H , A , Γ and J are fixed but D is allowed to vary, subject to theaxioms of non-commutative geometry.The axioms are solved in [5] to give explicit forms of the Dirac operatorin terms of n × n Hermitian matrices H and n × n anti-Hermitian tracelessmatrices L according to the formulas below. There are no other constraintson these n × n matrices, so these are the freely-specifiable data for the Diracoperator.The Dirac operator acts on H = V ⊗ M ( n, C ), with V = C k the spaceon which the gamma matrices act. The gamma matrices are assumed toform an irreducible representation of the Clifford algebra, which implies thatthe chirality operator is trivial for d = p + q odd. The dimension of V is k = 2 d/ for d even and k = 2 ( d − / for d odd. In the first two cases thesole gamma matrix is just 1 or − i respectively. In the remaining cases thegamma matrices are 2 × · , · ] denotes the commutator and { · , · } the anti-commutator ofmatrices. 5 ype (1,0) D = { H, · } (6) Type (0,1) D = − i [ L, · ] (7) Type (2,0) ( γ ) = ( γ ) = 1. D = γ ⊗ { H , · } + γ ⊗ { H , · } (8) Type (1,1) ( γ ) = 1 , ( γ ) = − D = γ ⊗ { H, · } + γ ⊗ [ L, · ] (9) Type (0,2) ( γ ) = ( γ ) = − D = γ ⊗ [ L , · ] + γ ⊗ [ L , · ] (10) Type (0,3) ( γ ) = ( γ ) = ( γ ) = − D = { H, · } + γ ⊗ [ L , · ] + γ ⊗ [ L , · ] + γ ⊗ [ L , · ] (11)A type ( p, q ) geometry has a signature s = ( q − p ) mod 8 which deter-mines some of the characteristics of the spectrum of D . These properties arewell-known, holding also for the case of a Riemannian geometry in dimension d , which is a type (0 , d ) spectral triple with signature s = d mod 8. Theproperties can be seen in the Monte Carlo simulations below. Symmetry
For s (cid:54) = 3 ,
7, if λ is an eigenvalue then so is − λ . Doubling
For s = 2 , s , the chirality operator Γis non-trivial. It is Hermitian and has eigenvalues ±
1. The Dirac operatorchanges the chirality, D Γ = − Γ D . If v is an eigenvector of eigenvalue λ thenΓ v is an eigenvector with eigenvalue − λ . As a result, the spectrum of theDirac operator is symmetric around 0. A similar argument holds for s = 1 , DJ = − J D so that v and J v have opposite eigenvalues.For the doubling property, if s = 2 , J = − DJ = J D in these cases, if v is an eigenvector then so is J v with the same eigenvalue. Moreover, one can check that v and J v mustbe linearly independent: suppose the eigenvectors are proportional to eachother, i.e.,
J v = cv , with c ∈ C , then − v = J v = J ( cv ) = ¯ cJ v = c ¯ cv , (12)which is a contradiction. 6 .2 A Monte Carlo algorithm for matrix geometries An observable f ( D ) is a real- or complex-valued function of Dirac operators.The expectation value of f is defined to be (cid:104) f (cid:105) = 1 Z (cid:90) f ( D ) e − S ( D ) d D. (13)The integral can be approximated as a sum over a discrete ensemble { D j , j =1 , . . . , N } . (cid:104) f (cid:105) N = (cid:80) j f ( D j ) e − S ( D j ) (cid:80) j e − S ( D j ) (14)so that in the limit taking N → ∞ , the average obtained through this dis-crete sum converges towards the continuum value (cid:104) f (cid:105) . This convergence canbe improved by using a Markov Chain Monte Carlo algorithm. In such an al-gorithm the Dirac operators D j are generated with a probability distributionsuch that P ( D j ) = e − S ( D j ) (cid:80) i e − S ( D i ) . (15)This simplifies the expression for the average (cid:104) f (cid:105) N = 1 N N (cid:88) j =1 f ( D j ) (16)and improves convergence by concentrating the sampling on regions whichcontribute strongly. To generate such an ensemble of Dirac operators theMetropolis-Hastings algorithm is used [16]. In this algorithm a proposed D j +1 is generated from D j by a move which will be defined in the nextsubsection. The proposed operator D p will be accepted as a new part of thechain, D j +1 = D p , if S ( D p ) < S ( D j ). If this was the only rule to add newoperators to the Markov chain the code would terminate in any sufficientlydeep local minimum. To make it possible to escape local minima, the newoperator is also accepted if exp( S ( D j ) − S ( D p )) > p , with p a uniformlydistributed random number in [0 , D p is rejected in both tests then D j +1 = D j . This algorithm ensures a Markov Chain satisfying detailedbalance, which ensures that the transition probability converges [22].After a sufficient number of moves, the probability distribution for D j converges towards the desired configuration and becomes independent of the7nitial state D . The states generated before this convergence are not rep-resentative of the probability distribution and can not be used to measureobservables. We checked that this burn-in process terminated by startingfrom different initial configurations and monitoring the convergence of theaction.The code is implemented in C++ and all matrix algebra operations usethe open source software library Eigen [14]. To construct a Markov Chain on the space of Dirac operators G = { D } , amove that proposes a new Dirac operator D p based on the last Dirac operator D j is needed. The Markov Chain property requires that the next proposedoperator can only depend on the current operator D j . The space G is a vectorspace, so a simple additive move D → D + δD , (17)with δD a Dirac operator, will always be ergodic, and as long as δD doesnot depend on past states the Markov property is also satisfied. As shown insection 2.1 the Dirac operator is defined using a choice of Hermitian matrices H i and anti-Hermitian matrices L i . To construct δD we define it as a Diracoperator composed from δH i , δL i . Generate a random n × n matrix R withmatrix elements in the complex range [ − − i, i ] and define δH i = l ( R i + R ∗ i ) δL i = l ( R i − R ∗ i )where l is a real constant that is determined at the start of each simulation.The value of l determines how ‘long’ the steps in the configuration space are.A Monte Carlo algorithm has the best thermalisation properties if the accep-tance rate of proposed moves is a r = ( / ( (cid:39) . D p generated). At the beginning ofa simulation the acceptance rate is tested and l adjusted, larger if the ac-ceptance rate is too large, smaller if the acceptance rate is too small, suchthat a r (cid:39) . τ MC .Note that the move for the L i does not preserve the condition that itis trace-free. However since the L i appear only in commutators, the tracedecouples and its value does not affect the Dirac operator.8 .4 Calculating the action The expression for the Dirac operator contains terms [ M, · ] and { M, · } for M ∈ M ( n, C ). The commutators and anti-commutators require the use ofthe left and right actions of M . These are written as matrices using thetensor product [ M, · ] = M ⊗ I n − I n ⊗ M T (18) { M, · } = M ⊗ I n + I n ⊗ M T . (19)The Dirac operator D can then be written as a kn × kn matrix that actson the tensor space V ⊗ C n ⊗ C n .The matrix operations needed in the computer code are matrix multipli-cation, addition and calculation of eigenvalues. The run time of these growslike O ( m ), O ( m ) and O ( m ) respectively for m × m matrices. Therefore itmakes sense to write the action in terms of the much smaller n × n matrices H i , L i to accelerate the simulations. The details of this calculation for thegeometries investigated are collected in appendix A. Given a Dirac operator D , the eigenvalues { λ i } can be computed. The twomain observables of interest are the j -th eigenvalue, ordering the eigenvaluesfrom lowest to highest f j ( D ) = λ j (20)and the distribution of eigenvalues at eigenvalue λf λ ( D ) = 1 kn (cid:88) j δ ( λ − λ j ) (21)Since eigenvalue calculations are computationally expensive, the eigen-values are only measured every 4 n attempted Monte Carlo moves. Thisimproves run time, and reduces the correlation of the measurements. Theaction S and the acceptance rate of moves are recorded at every step tomonitor the algorithm.At later points it will become useful to measure some additional observ-ables that are computed from the eigenvalues, for example, Tr D . For certaincases it has also proven instructive to examine the non-physical degrees offreedom of the matrices H and L via their eigenvalues.The average of any observable can be calculated directly from the set ofmeasurements. However to estimate the statistical error on our measure-ments it is necessary to take the correlation between successive states in9
00 1000 1500 Τ Mc - - (a) Autocorrelation of S for n = 5 ,
500 1000 1500 Τ Mc - - (b) Autocorrelation of λ min for n = 5 , Figure 1: Fall-off of the autocorrelation for the action and the minimumeigenvalue for a type (1 ,
0) geometry with S = Tr D . The blue line is n = 5and the yellow line is n = 15. The horizontal axis is Monte Carlo time.the Markov Chain into account. The error bars shown on plots of averageeigenvalues show the statistical error Err( λ i ) calculated asErr( λ i ) = (cid:114) σ λ i τ A,λ i M (22)with σ λ i the variance of the eigenvalue, τ A,λ i the integrated autocorrelationtime of the eigenvalue and M the number of measurements performed [22].In figure 1, autocorrelations for the simulations of a type (1 ,
0) geometrywith S = Tr D for size n = 5 and n = 15 are shown. The figures show theautocorrelation for both the action and the smallest eigenvalue of D .The autocorrelation time is determined on the data after the burn-in iscompleted. In practice the burn-in phase was combined with the adjustmentof l . Then τ MC was counted from 5000 moves after the last adjustment to l .This burn-in and adjustment period takes up most of the simulations. Wefound that for the eigenvalue distribution and the average eigenvalues, 200measurements (corresponding to 4 n ·
200 attempted Monte Carlo moves) leadto very good results. To determine the phase transition, 2000 measurementswere used to ensure that statistical fluctuations were not mistaken for a phasetransition. D action In this section the Monte Carlo simulations for the simplest possible action S = Tr D are examined. The one-dimensional Clifford algebras, type (1 , ,
1) are examined first and the results understood using the standard10heory of Gaussian matrix models. After this, some numerical results for thetwo- and three-dimensional types are shown. (1 , and (0 , The n eigenvalues of Dirac operators (6), (7) can be written in terms of the n eigenvalues µ j of the matrix H or the eigenvalues iµ j of L . For the (1 , λ jk = µ j + µ k (23)while for the (0 ,
1) case λ jk = µ j − µ k (24)This follows from the fact that eigenvectors of D are of the form u j ⊗ u k ,with u j the eigenvectors of H or L .The first point is that the (0 ,
1) case has eigenvalue 0 with multiplicity n given by the terms j = k . This can also be seen directly from the Diracoperator: all matrices in M ( n, C ) that commute with L have eigenvalue 0,and there are always at least n linearly independent such matrices. It willbe seen later that a peak in the eigenvalue distribution at, or near, 0 is afeature of some other random fuzzy spaces.The second point is that the spectrum of the (0 ,
1) case is symmetricabout the origin, as λ jk = − λ kj . This is in accordance with its signature s = 1, which means that each Dirac operator has symmetric spectrum. Thespectrum of a (1 ,
0) Dirac operator is typically not symmetric since s = 7in this case. This means that our simulation gives an eigenvalue distribu-tion that is not exactly symmetric, though it will eventually converge to asymmetric distribution as the Monte Carlo time increases.For the (1 ,
0) case, using the simplified action (47) one can transform theintegral over the Dirac operator into an integral over the Hermitian matrix H . S (1 , ( D ) = Tr D = 2 n Tr H + 2(Tr H ) (25)= 2 n (cid:88) i µ i + 2 (cid:88) i (cid:88) j µ i µ j (26)The (0 ,
1) case is similar, but one has to take into account the fact thatthe integration over Dirac operators is an integration over traceless matrices L . Using (51) gives S (0 , ( D ) = Tr D = − n Tr L (27)= 2 n (cid:88) i µ i , (28)11 - - - X Μ i \ (a) Type (1 ,
0) average eigenvalues of H - - - X Μ i \ (b) Type (0 ,
1) average eigenvalues of − iL - - X Λ i \ (c) Type (1 ,
0) average eigenvalues of D - - X Λ i \ (d) Type (0 ,
1) average eigenvalues of D Figure 2: Average ordered eigenvalues for
H, L and D for the cases (1 ,
0) and(0 ,
1) with n = 5.An example of average ordered eigenvalues generated by the Monte Carlosimulation is shown in figure 2.These random matrix models are close to the Gaussian Hermitian matrixmodel [20, 6], which has the similar action (cid:101) S ( M ) = 2 n Tr M = 2 n (cid:88) j µ j , (29)with integration over all Hermitian matrices.A standard technique in random matrix models is to calculate the jointprobability density for the eigenvalues µ j . The formula is [21] P ( µ , . . . , µ n ) = C exp (cid:16) − (cid:101) S ( µ , . . . , µ n ) (cid:17) (cid:89) j 0) model transforms directly tothe Gaussian Hermitian matrix model by rescaling the trace by a factor √ M = H + 1 n ( √ − 1) Tr H. (31)A standard result (the Wigner semicircle law [26]) is that the analogue ofthe eigenvalue distribution (21) for the Gaussian Hermitian matrix modelconverges as n → ∞ to the density of states σ ( µ ) = lim n →∞ < f µ ( M ) > = (cid:40) πA (cid:112) A − µ for − √ A ≤ µ ≤ √ A A = 1.In our simulations using actions S (1 , and S (0 , we find that the semicirclelaw is also a good approximation for the eigenvalues of H and L . It is alreadywell-satisfied for n = 5 and improves at higher n , as shown in figure 3. Thereason for this is that in the Gaussian Hermitian matrix model, the variable n Tr M is normally-distributed with variance 1 / (4 n ), and so adjusting theeigenvalues with a fixed multiple of this makes no difference to the densityof states in the limit n → ∞ .Another standard result from random matrix theory is that the correla-tion between different fixed eigenvalues of M vanishes as n → ∞ . Thus forlarge n , the joint probability distribution away from the diagonal µ = µ is simply the product of the density of states [24]. Therefore, for large n ,one can calculate the eigenvalue distribution of the Dirac operator from thesemicircle law as a convolution, with a correction for the behaviour of thecorrelations on the diagonal.This is shown as follows. Let f ( λ ) be an observable for a random fuzzyspace, with λ = µ ± µ . Then assuming a product probability distribution,one has (cid:104) f (cid:105) = (cid:90) σ ( µ ) σ ( µ ) f ( µ ± µ )d µ d µ = (cid:90) σ D ( λ ) f ( λ )d λ (33)with density of states for the Dirac operator the convolution σ D ( λ ) = (cid:90) σ ( λ ∓ µ ) σ ( µ )d µ , (34)13 - Μ H Μ L (a) Type (1 , n = 2 - - Μ H Μ L (b) Type (0 , n = 2 - - Μ H Μ L (c) Type (1 , n = 5 - - Μ H Μ L (d) Type (0 , n = 5 - - Μ H Μ L (e) Type (1 , n = 15 - - Μ H Μ L (f) Type (0 , n = 15 Figure 3: The semicircle law is compared with the density of states for H or L . 14 - Λ H Λ L (a) (1 , n = 5 - - - Λ H Λ L (b) (0 , n = 5 Figure 4: The eigenvalue density for the Dirac operator compared with theconvolution of two semicircle functions σ D , with correction applied in the(0 , 1) case.which is an elliptic integral. This integral is the same for type (1 , 0) and(0 , σ ( µ ) = σ ( − µ ). The Monte Carlo simulation of the eigenvaluedensity at finite n is shown in figure 4. The continuous line is the curve σ D ( λ ) for the (1 , 0) case but a significant correction term is added to σ D forthe (0 , 1) case.The correction to the product probability density gives a contributiononly near the diagonal µ = µ . The approximate form is an additionalcontribution to (cid:104) f (cid:105) of [24] − (cid:90) sin (cid:0) πn ( µ − µ ) σ (( µ + µ ) / (cid:1) π n ( µ − µ ) f ( µ ± µ ) d µ d µ (35)For the (0 , 1) case ( − ), this formula contributes significantly near λ = 0,accounting for the gap at the origin in figure 4(b) with a width that scalesas 1 /n . While the spectra for geometries with one-dimensional Clifford algebra areeasy to understand, those with a two-dimensional Clifford algebra are lessstraightforward. The average eigenvalues and the eigenvalue distributions areshown in figure 5 for the case n = 5 and in figure 6 for the larger matrices n = 15. The individual eigenvalues are more easily seen in figure 5. All threetypes are symmetric about the origin and the third one exhibits eigenvaluedoubling, all in accordance with the properties for s = 6, 0 and 2 derived insection (2.1). 15he action Tr D is 2 times the sum of quadratic actions for each matrix L i or H i , these quadratic actions being exactly the (0 , 1) and (1 , 0) actionspreviously analysed. In particular, these matrices are statistically indepen-dent. The eigenvalues of the H i , L i are still approximated well by the semi-circle law (32) with A = 1 / 2. However, the main difference in analysing thetwo-dimensional cases is that the eigenvalues of the Dirac operator are notsimply related to the eigenvalues of L i , H i .For the case (0 , 2) the multiplicity of eigenvalue 0 is at least 2 n , as shownby examining the Dirac operator D (0 , = γ ⊗ [ L , · ] + γ ⊗ [ L , · ] (36)= γ (cid:0) [ L , · ] − γ γ ⊗ [ L , · ] (cid:1) (37)= γ (cid:18) [ L + iL , · ] 00 [ L − iL , · ] (cid:19) , (38)using a basis so that γ γ = diag( i, − i ). The Dirac operator acts on thespace C ⊗ M ( n, C ). All v ⊗ m in this space for which (cid:18) [ L + iL , · ] 00 [ L − iL , · ] (cid:19) (cid:18) v mv m (cid:19) (39)= (cid:18) v [ L + iL , m ] v [ L − iL , m ] (cid:19) = 0 (40)have eigenvalue 0. Picking a basis on v one can choose v = 1 , v = 0and v = 0 , v = 1. There will then be n linearly independent matrices m that commute with L + iL and n that commute with L − iL , hence 2 n eigenvalues equal to 0 for the (0 , 2) type geometry. Just as in the (0 , 1) case,there is a gap in the eigenvalue spectrum around the spike at 0. This showsthere is eigenvalue repulsion for this Dirac operator also, though we do nothave a theoretical understanding of this phenomenon.The types (2 , 0) and (1 , 1) also have a feature at the origin. The densityof eigenvalues is sharply lower in a narrow dip at the origin and there isa somewhat wider upward spike around this. This is shown for the (2 , λ and the opposite eigenvalue − λ that is required by the symmetryof the spectrum of D about 0.The numerical results also indicate that the range of the eigenvalues re-mains unchanged under the change of matrix size, and the distribution be-comes smoother, appearing to converge to a smooth limiting distribution inthe same way as for random matrices. Another similar feature is that for16 - - X Λ i \ (a) Type (2 , - - - Λ H Λ L (b) Type (2 , 10 20 30 40 50 i - - X Λ i \ (c) Type (1 , - - - Λ H Λ L (d) Type (1 , 10 20 30 40 50 i - - X Λ i \ (e) Type (0 , - - - Λ H Λ L (f) Type (0 , Figure 5: The average eigenvalues, and the histograms of the eigenvaluedistribution for the different types of two-dimensional Clifford algebra. Theaction is S ( D ) = Tr D and the matrix size n = 5.17 00 200 300 400 i - - - X Λ i \ (a) Type (2 , - - - Λ H Λ L (b) Type (2 , 100 200 300 400 i - - - X Λ i \ (c) Type (1 , - - - Λ H Λ L (d) Type (1 , 100 200 300 400 i - - - X Λ i \ (e) Type (0 , - - - Λ H Λ L (f) Type (0 , Figure 6: The average eigenvalues, and the histograms of the eigenvaluedistribution for the different types of two-dimensional Clifford algebra. Theaction is S ( D ) = Tr D and the matrix size n = 15.18 10 15 20 25 30 i - - - X Λ i \ (a) Average eigenvalues - - Λ H Λ L (b) Eigenvalue distribution Figure 7: Zooming in to a region near eigenvalue 0. Type (2 , n = 15. - - - Λ H Λ L (a) Type (2 , n = 5 - - - Λ H Λ L (b) Type (2 , n = 10 - - - Λ H Λ L (c) Type (2 , n = 15 - - - Λ H Λ L (d) Type (0 , n = 5 - - - Λ H Λ L (e) Type (0 , n = 10 - - - Λ H Λ L (f) Type (0 , n = 15 Figure 8: The distribution of single eigenvalues at different sizes.19 00 200 300 400 i - - X Λ i \ (a) Type (0 , - - Λ H Λ L (b) Type (0 , Figure 9: The average eigenvalues and the eigenvalue distribution for type(0 , S = Tr D and the matrix size n = 15.larger n the fluctuation of each individual eigenvalue becomes smaller. Thiscan be seen in figure 8. The leftmost bump in each plot is the smallest eigen-value while the rightmost bump is the largest, the eigenvalues in betweenwere chosen to be symmetric, and include the central most eigenvalues.The eigenvalue distribution for the type (0 , 3) case is plotted in figure9. This appears to be smooth at the origin, similar to the (1 , 0) case. Thecommon property of these cases is that the Dirac operators do not have asymmetric spectrum. Thus a small eigenvalue λ does not have to be closeto any other eigenvalue. There is nothing special about the origin, and inparticular, the eigenvalue repulsion hypothesis does not lead to any specialbehaviour here. D term The Tr D term in the action leads to interactions between the L i , H i thatcompose the Dirac operator. An extreme case of this is for type (0 , H , L , L and L is present.These terms make it harder to understand the system analytically, howeverfor the simulations they are no obstacle.The simple action S ( D ) = Tr D leads to behavior very similar to thatfor the action Tr D . This is shown in figure 10. Some characteristics, likethe shoulders for type (1 , 1) and (2 , 0) are more pronounced, but the overallshape is quite similar.Combining the two terms together gives the action S = Tr (cid:0) g D + D (cid:1) . (41)20 - - Λ H Λ L (a) (1 , n = 15 - - - Λ H Λ L (b) (0 , n = 15 - - Λ H Λ L (c) (2 , n = 15 - - Λ H Λ L (d) (1 , n = 15 - - Λ H Λ L (e) (0 , n = 15 - - - Λ H Λ L (f) (0 , n = 15 Figure 10: The eigenvalue distribution for the action S = Tr D .21 - Λ H Λ L Figure 11: The potential V = λ + g λ for g = − − . − − . − − . − − . − 5. The lines are coloured from red ( g = − 1) through toyellow ( g = − g the behaviour of the numerical simulations is some-where between the Tr D case and the Tr D and does not show qualitativelynew features. However when g is negative this is a symmetry-breakingpotential with two minima, shown in figure 11. The question of how theeigenvalues behave in this case is interesting and a variety of behaviours isexhibited depending on the type of the gamma matrices. This is shown infigure 12. The different types are described here for values of g decreasingfrom 0.(1 , 0) Two peaks start to form at around g = − g = − g = − . 5, leaving the centre of thedistribution empty. Since the Dirac operator is not symmetric, theMonte Carlo simulation can and does settle in just one of the peaks,though one expects that the Markov Chain would eventually exploreboth peaks equally given a long enough run.(0 , 1) Two peaks form at around g = − L settle into two peaks, and since it istraceless, the favoured configuration has the same number of eigenval-ues in each peak. The differences between eigenvalues of L in the sameminimum remains small, giving the central peak in the distribution for D . The n eigenvalues exactly 0 are also still present.(2 , 0) Two peaks develop at small g and grow until the central part of thedistribution vanishes suddenly between g = − . − , 1) Two peaks develop at small g and grow until the central part of the22 - Λ H Λ L (a) Type (1 , - - Λ H Λ L (b) Type (0 , - - Λ H Λ L (c) Type (2 , - - Λ H Λ L (d) Type (1 , - - Λ H Λ L (e) Type (0 , - - Λ H Λ L (f) Type (0 , Figure 12: The eigenvalues of S = Tr ( D + g D ) for n = 10 and g = − − . − − . − − . − − . − 5. The lines are coloured from red( g = − 1) through to yellow ( g = − g = − − . , 2) This is the most mysterious case. Two slight peaks develop but theeigenvalues do not separate into two peaks for the whole of the rangeof g tested. Instead some further substructure to the eigenvalue dis-tribution develops.(0 , 3) This is similar to the (0 , 1) case, with the sharp change occurringbetween g = − g = − . 5. In figure 12f the Markov Chainfor g = − . ∂ log Z∂g = (cid:104) Tr D (cid:105) (42)and the autocorrelation time, which is usually expected to increase near acontinuous phase transition due to the long-range order (‘critical slowingdown’). These are plotted in figure 13. One can see that there is goodevidence for a phase transition for the types (1 , , , 1) and (0 , g described above, and the autocorrelation time of Tr D has a peak aroundthis value also. It is difficult to see any clear signal from the plots for types(0 , 1) and (0 , H . Unlike the L matrices, thetrace of H appears to play a crucial role. The observable F = (Tr H ) n Tr H (43)measures the strength of Tr H , calculated as a square so that positive andnegative values do not cancel, and as a fraction of the total strength of H .The averages of this are plotted in figure 14. In the case of (2 , H and H , and the F for both combined is F = (Tr H ) + (Tr H ) n Tr ( H + H ) (44)In both cases, F = 1 if the matrices are pure trace. The plots show thatTr H develops a large expectation value at the phase transition. In the caseof (2 , 0) the Monte Carlo data used for the plots developed a preference forTr H rather than Tr H but this is of no significance due to the rotationalsymmetry between the two . The sum of squares is the correct rotationally-invariant observable. We have checked this with additional simulations and found that the trace degree offreedom is in general distributed randomly between the two matrices. - - - - g X T r D \ Τ H T r D L (a) Type (1 , - - - - - g X T r D \ Τ H T r D L (b) Type (0 , - - - - - g X T r D \ Τ H T r D L (c) Type (2 , - - - - - g X T r D \ Τ H T r D L (d) Type (1 , - - - - - g X T r D \ Τ H T r D L (e) Type (0 , - - - - - g X T r D \ Τ H T r D L (f) Type (0 , Figure 13: The mean of order parameter Tr D and the autocorrelation time τ as g varies. 25 - - - g (a) Type (1 , - - - - g (b) Type (1 , - - - - g (c) Type (2 , - - - - g (d) Type (0 , Figure 14: Fraction F measuring the square of the proportion of H that is inthe trace part of H as g varies. The plot for type (2 , 0) shows the fraction F for H (red), F (green) and combining both F (blue).26 Conclusion A model of random geometry has been presented here as random Dirac op-erators in non-commutative geometry. The integrals can be interpreted asmulti-matrix models but with a new type of observables, namely the eigen-values of the Dirac operator. The one-dimensional cases can be understoodusing theoretical results from random matrices but the higher-dimensionaltypes are not so easy and will require further study to obtain analytic results.Numerical results have been presented showing various phenomena that de-pend strongly on the type of gamma matrices used, particularly whetherthe spectrum of a Dirac operator for that signature is symmetric about theorigin.From the numerical results it is clear that some features are similar tothe properties of the eigenvalues of random matrices: the eigenvalue distribu-tions appear to converge in the large n limit and the dispersion of individualordered eigenvalues decreases; also there is some evidence of a degree ofeigenvalue repulsion at the origin.The most interesting results are for the quartic action with negative g , sothat the potential is of symmetry-breaking type. For some types, the eigen-value spectrum changes suddenly when g reaches a critical negative valueand the observable Tr D is a good order parameter for this transition. Thisis taken as a strong indication that a sharp phase transition would occur inthe large n limit. The types where this transition is clear are those wherethe Dirac operator contains at least one term involving an anti-commutatorwith a random hermitian matrix H . Then the trace of H develops a largeexpectation value, becoming the dominant contribution to H after the tran-sition. This can’t happen with commutator terms as the trace of the randommatrix decouples in this instance.For generic g , the eigenvalue distribution of D is not a good approxi-mation to the behaviour for the Dirac operator on any fixed (commutative)Riemannian manifold (5), except that one could possibly argue that the dis-tribution is approximately constant for some ranges (e.g. figure 10(f), whichlooks like a one-dimensional manifold). The exception to this is near thephase transition, where the curves in figure 12 do appear to have the rightpower-law behaviour. Two of these distributions are highlighted in figure 15,showing the types (1 , 1) and (2 , 0) at values of g just below the value for thephase transition.These are compared with the eigenvalue distribution for the fuzzy spherefrom [5], shown in figure 15(c). This is a type (1 , 3) spectral triple, havingsignature s = 2, and has exactly the same spectrum as the Dirac operatoron the Riemannian round S but with a maximum eigenvalue cut-off and27 - - Λ H Λ L (a) Type (1 , g = − . - - - Λ H Λ L (b) Type (2 , g = − - - Λ H Λ L (c) Type (1,3) Fuzzy S Figure 15: Eigenvalue distributions near the phase transition compared withthe fuzzy sphere. Matrix size n = 10.fermion doubling. The distributions are remarkably similar, the main differ-ences being the gap at the origin, the size of which depends on the distancefrom the phase transition, and the fact that the fuzzy sphere has exactlyinteger eigenvalues with multiplicity, due to its rotational symmetry. Thefeature that is common to the plots is the approximately linear increase ofthe eigenvalue density with eigenvalue that is characteristic of Riemannianmanifolds of dimension two, i.e., m = 2 in (5). The simulations show thatdecreasing g further increases the gap in the middle of the spectrum whereasthe middle of the spectrum fills up for values of g greater than the criticalvalue.These results are somewhat preliminary and we have not yet carried outa systematic study of the phase transition.The study of random non-commutative geometries gives an insight intothe closely-related problem of the quantization of this fascinating theory. Aquantized non-commutative space is a potential candidate for a quantumtheory of gravitational interactions and will allow a better understandingof fundamental interactions. Independent of these physical applications, it28lso is an interesting modification of the well-known matrix models, intro-ducing non-trivial interactions and observables among several matrices. Theresults reported here will be a basis for further investigations into the phasetransition, the continuum limit, and other possible actions on this space ofgeometries. This work of JWB on this project was supported by STFC Particle PhysicsTheory Consolidated Grant ST/L000393/1, and LG was supported by fund-ing from the European Research Council under the European Union SeventhFramework Programme (FP7/2007-2013) / ERC Grant Agreement n.306425“Challenging General Relativity”.This research was supported in part by Perimeter Institute for Theo-retical Physics. Research at Perimeter is supported by the Government ofCanada through Industry Canada and by the Province of Ontario throughthe Ministry of Economic Development and Innovation.The dissemination of this work was aided by EU COST Action MP1405‘Quantum structure of spacetime’. A Dirac operators for the fuzzy geometrieswe examined In this appendix matrices H will be Hermitian and matrices L anti-Hermitian.The L matrices are not assumed to be traceless, though the actions are in-dependent of trace part of these matrices. The bracketing convention forthe trace of an expression is Tr AB ≡ Tr( AB ) and Tr A n ≡ Tr( A n ), butTr A + B ≡ (Tr A ) + B . A.1 The (1 , geometry γ = 1 (45) D = { H, ·} (46)Tr D = 2 n Tr H + 2(Tr H ) (47)Tr D = 2 n Tr H + 8 Tr H Tr H + 6(Tr H ) (48)29 .2 The (0 , geometry γ = − i (49) D = γ ⊗ [ L , · ] (50)Tr D = − (cid:0) n Tr L − L ) (cid:1) (51)Tr D = 2 n Tr L − L Tr L + 6(Tr L ) (52) A.3 The (2 , geometry γ = (cid:18) − (cid:19) γ = (cid:18) (cid:19) (53) D = γ ⊗ { H , ·} + γ ⊗ { H , ·} (54)The gamma matrix trace identities are Tr ( γ i γ j ) = 2 δ ij and Tr (cid:0) γ i γ j γ k (cid:1) = 0and Tr (cid:0) γ i γ j γ k γ l (cid:1) = 2( δ ij δ kl − δ ik δ jl + δ il δ jk ).Tr D = 4 n (Tr H + Tr H ) + 4 (cid:0) (Tr H ) + (Tr H ) (cid:1) (55)Tr D = 2Tr (cid:0) { H , ·} (cid:1) + 2Tr (cid:0) { H , ·} (cid:1) + 8Tr (cid:0) { H , ·} { H , ·} (cid:1) (56) − { H , ·} { H , ·} { H , ·} { H , ·} ) (57)= 4 n (cid:18) Tr H + Tr H + 4 Tr H H − H H H H (cid:19) (58)+ 16 (cid:18) Tr H (cid:0) Tr H + Tr H H (cid:1) (59)+ Tr H (cid:0) Tr H H + Tr H (cid:1) + (Tr H H ) (cid:19) (60)+ 12 (cid:18) (Tr H ) + (Tr H ) (cid:19) + 8 Tr H Tr H (61) A.4 The (1 , geometry γ = (cid:18) − (cid:19) γ = (cid:18) − (cid:19) (62)30 = γ ⊗ { H, ·} + γ ⊗ [ L, · ] (63)The gamma matrix trace identities are Tr ( γ i γ j ) = 2 η ij , where η i,j = diag( − , (cid:0) γ i γ j γ k (cid:1) = 0 and Tr (cid:0) γ i γ j γ k γ l (cid:1) = 2( η ij η kl − η ik η jl + η il η jk ).Tr D = 4 n (Tr H − Tr L ) + 4 (cid:0) (Tr H ) + (Tr L ) (cid:1) (64)Tr D = 2Tr (cid:0) { H, ·} (cid:1) + 2Tr (cid:0) [ L, · ] (cid:1) − (cid:0) { H, ·} [ L, · ] (cid:1) (65)+ 4Tr ( { H, ·} [ L, · ] { H, ·} [ L, · ]) (66)= 4 n (cid:18) Tr H + Tr L − H L + 2 Tr HLHL (cid:19) (67)+ 16 (cid:18) Tr H (cid:0) Tr H − Tr L H (cid:1) (68)+ Tr L (cid:0) − Tr L + Tr H L (cid:1) + (Tr HL ) (cid:19) (69)+ 12 (cid:18) (Tr H ) + (Tr L ) (cid:19) − H Tr L (70) A.5 The (0 , geometry γ = (cid:18) i − i (cid:19) γ = (cid:18) − (cid:19) (71) D = γ ⊗ [ L , · ] + γ ⊗ [ L , · ] (72)The gamma matrix trace identities are Tr ( γ i γ j ) = − δ ij and Tr (cid:0) γ i γ j γ k (cid:1) = 0and Tr (cid:0) γ i γ j γ k γ l (cid:1) = 2( δ ij δ kl − δ ik δ jl + δ il δ jk ).Tr D = − n (Tr L + Tr L ) + 4 (cid:0) (Tr L ) + (Tr L ) (cid:1) (73)Tr D = 2Tr (cid:0) [ L , · ] (cid:1) + 2Tr (cid:0) [ L , · ] (cid:1) + 8Tr (cid:0) [ L , · ] [ L , · ] (cid:1) (74) − L , · ] [ L , · ] [ L , · ] [ L , · ]) (75)= 4 n (cid:18) Tr L + Tr L + 4 Tr L L − L L L L (cid:19) (76)+ 16 (cid:18) − Tr L (Tr L + Tr L L ) (77) − Tr L (Tr L + Tr L L ) + (Tr L L ) (cid:19) (78)+ 12 (cid:18) (Tr L ) + (Tr L ) (cid:19) + 8 Tr L Tr L (79)31 .6 The (0 , geometry γ = (cid:18) i − i (cid:19) γ = (cid:18) − (cid:19) γ = (cid:18) ii (cid:19) γ = (cid:18) − − (cid:19) (80) D = { H, } + γ ⊗ [ L , · ] + γ ⊗ [ L , · ] + γ ⊗ [ L , · ] (81)To calculate the action, the following gamma matrix identities are used:Tr ( γ i γ j ) = − δ ij and Tr (cid:0) γ i γ j γ k γ l (cid:1) = 2( δ ij δ kl − δ ik δ jl + δ il δ jk ), and γ i γ j γ k = ε ijk , the totally antisymmetric tensor with ε = − D = 4 n (Tr H − Tr L − Tr L − Tr L ) (82)+ 4((Tr H ) + (Tr L ) + (Tr L ) + (Tr L ) ) (83)Tr D = 2Tr (cid:0) { H, ·} (cid:1) + 2Tr (cid:0) [ L , · ] (cid:1) + 2Tr (cid:0) [ L , · ] (cid:1) + 2Tr (cid:0) [ L , · ] (cid:1) (84)+ 8Tr (cid:0) [ L , · ] [ L , · ] (cid:1) + 8Tr (cid:0) [ L , · ] [ L , · ] (cid:1) (85)+ 8Tr (cid:0) [ L , · ] [ L , · ] (cid:1) − (cid:0) { H, ·} [ L , · ] (cid:1) (86) − (cid:0) { H, ·} [ L , · ] (cid:1) − (cid:0) { H, ·} [ L , · ] (cid:1) (87) − L , · ] [ L , · ] [ L , · ] [ L , · ]) − L , · ] [ L , · ] [ L , · ] [ L , · ]) (88) − L , · ] [ L , · ] [ L , · ] [ L , · ]) − { H, ·} [ L , · ] { H, ·} [ L , · ]) (89) − { H, ·} [ L , · ] { H, ·} [ L , · ]) − { H, ·} [ L , · ] { H, ·} [ L , · ])(90) − { H, ·} [ L , · ] [ L , · ] [ L , · ]) + 8Tr ( { H, ·} [ L , · ] [ L , · ] [ L , · ]) (91) − { H, ·} [ L , · ] [ L , · ] [ L , · ]) + 8Tr ( { H, ·} [ L , · ] [ L , · ] [ L , · ]) (92) − { H, ·} [ L , · ] [ L , · ] [ L , · ]) + 8Tr ( { H, ·} [ L , · ] [ L , · ] [ L , · ]) (93)= 4 n (cid:18) Tr (cid:0) L + L + L + H (cid:1) (94) − L L L L + L L L L + L L L L ) (95) − HL HL + HL HL + HL HL ) (96)+ 4Tr (cid:0) L L + L L + L L (cid:1) − (cid:0) H L + H L + H L (cid:1) (97) − HL L L + HL L L + HL L L ) (98)+ 4Tr ( HL L L + HL L L + HL L L ) (cid:19) (99) − 16 Tr L (cid:0) Tr L + Tr L L + Tr L L − H L (cid:1) (100) − 16 Tr L (cid:0) Tr L + Tr L L + Tr L L − H L (cid:1) (101) − 16 Tr L (cid:0) Tr L + Tr L L + Tr L L − H L (cid:1) (102)32 16 Tr H (cid:0) Tr H − L H − L H − L H (cid:1) (103)+ 16 (cid:18) Tr L Tr ( H [ L , L ]) − Tr L Tr ( H [ L , L ]) (104)+ Tr L Tr ( H [ L , L ]) − H Tr ( L [ L , L ]) (cid:19) (105)+ 12 (cid:0) (Tr L ) + (Tr L ) + (Tr L ) + (Tr H ) (cid:1) (106)+ 8 (cid:0) Tr L Tr L + Tr L Tr L + Tr L Tr L (cid:1) (107) − 24 Tr H (cid:0) Tr L + Tr L + Tr L (cid:1) (108)+ 16 (cid:0) (Tr L L ) + (Tr L L ) + (Tr L L ) (cid:1) (109)+ 48 (cid:0) (Tr HL ) + (Tr HL ) + (Tr HL ) (cid:1) (110)This expression shows that the complexity of terms rises quickly with in-creasing d = p + q . A.7 Powers of (anti-)commutators For anti commutators we have { H, ·} = H ⊗ I n + I n ⊗ H T (111) { H, { H, ·}} = I n ⊗ H T H T + 2 H ⊗ H T + HH ⊗ I n (112)Tr ( { H, { H, ·}} ) = 2 n Tr H + 2(Tr H ) (113) { H, { H, { H, { H, ·}}}} = (cid:0) I n ⊗ H T H T + 2 H ⊗ H T + HH ⊗ I n (cid:1) (114)= I n ⊗ ( H T ) + 4 H ⊗ ( H T ) + 6 H ⊗ ( H T ) (115)+ 4 H ⊗ H T + H ⊗ I n (116)Tr ( { H, { H, { H, { H, ·}}}} ) = 2 n Tr H + 8 Tr H Tr H + 6(Tr H ) (117)And for commutators[ L, · ] = L ⊗ I n − I n ⊗ L T (118)[ L, [ L, · ]] = I n ⊗ L T L T − L ⊗ L T + LL ⊗ I n (119)Tr ([ L, [ L, · ]]) = 2 n Tr L − L ) (120)[ L, [ L, [ L, [ L, · ]]]] = (cid:0) I n ⊗ L T L T − L ⊗ L T + LL ⊗ I n (cid:1) (121)= I n ⊗ ( L T ) − L ⊗ ( L T ) + 6 L ⊗ ( L T ) (122) − L ⊗ L T + L ⊗ I n (123)Tr ([ L, [ L, [ L, [ L, · ]]]]) = 2 n Tr L − L Tr L + 6(Tr L ) (124)33 eferences [1] J. Ambjørn, A. Goerlich, J. Jurkiewicz, and R. Loll. NonperturbativeQuantum Gravity. arXiv:1203.3591 , March 2012.[2] Jan Ambjørn, Bergfinnur Durhuus, and Thordur Jonsson. Quantum Ge-ometry . Cambridge Monographs on Mathematical Physics. CambridgeUniv. Press, Cambridge, UK, 2005.[3] Takehiro Azuma, Subrata Bal, Keiichi Nagao, and Jun Nishimura. Non-perturbative studies of fuzzy spheres in a matrix model with the Chern-Simons term. JHEP , 05:005, 2004.[4] John W. Barrett. A Lorentzian version of the non-commutative geome-try of the standard model of particle physics. J. Math. Phys. , 48:012303,2007.[5] John W. Barrett. Matrix geometries and fuzzy spaces as finite spectraltriples. J. Math. Phys. , 56:082301, 2015.[6] S. Chadha, G. Mahoux, and M. L. Mehta. A method of integration overmatrix variables: II. Journal of Physics A: Mathematical and General ,14(3):579, March 1981.[7] Ali H. Chamseddine and Alain Connes. The Spectral Action Principle. Communications in Mathematical Physics , 186(3):731–750, July 1997.arXiv: hep-th/9606001.[8] Ali H. Chamseddine and Alain Connes. The Uncanny Precision of theSpectral Action. Communications in Mathematical Physics , 293(3):867–897, February 2010. arXiv: 0812.0165.[9] Giovanni M. Cicuta. Phase transitions and random matrices. In Randommatrix models and their applications , volume 40 of Math. Sci. Res. Inst.Publ. , pages 95–109. Cambridge Univ. Press, Cambridge, 2001.[10] Alain Connes. Noncommutative geometry . Academic Press, Inc., SanDiego, CA, 1994.[11] Alain Connes. Noncommutative geometry and the standard model withneutrino mixing. JHEP , 11:081, 2006.[12] Luiz C. de Albuquerque, Jorge L. deLyra, and Paulo Teotonio-Sobrinho.Fluctuating dimension in a discrete model for quantum gravity basedon the spectral action. Phys. Rev. Lett. , 91:081301, 2003.3413] H. Grosse and P. Presnajder. The Dirac operator on the fuzzy sphere. Lett. Math. Phys. , 33:171–182, 1995.[14] Ga¨el Guennebaud, Benoˆıt Jacob, et al. Eigen v3.http://eigen.tuxfamily.org, 2010.[15] Herbert W. Hamber. Quantum Gravity on the Lattice. Gen. Rel. Grav. ,41:817–876, 2009.[16] W. K. Hastings. Monte Carlo sampling methods using Markov chainsand their applications. Biometrika , 57(1):97–109, April 1970.[17] Joe Henson, David P. Rideout, Rafael D. Sorkin, and Sumati Surya.Onset of the Asymptotic Regime for Finite Orders. arXiv:1504.05902[gr-qc] , April 2015. arXiv: 1504.05902.[18] Alexander Holfter and Mario Paschke. Moduli spaces of discrete gravity.1. A Few points... J. Geom. Phys. , 47:101, 2003.[19] J. Madore. The Fuzzy sphere. Class. Quant. Grav. , 9:69–88, 1992.[20] M. L. Mehta and M. Gaudin. On the density of Eigenvalues of a randommatrix. Nuclear Physics , 18:420–427, August 1960.[21] Madan Lal Mehta. Random Matrices . Academic Press, October 2004.[22] M. E. J. Newman and G. T. Barkema. Monte Carlo Methods in Statis-tical Physics . Clarendon Press, February 1999.[23] Denjoe O’Connor, Brian P. Dolan, and Martin Vachovski. Critical Be-haviour of the Fuzzy Sphere. JHEP , 12:085, 2013.[24] L. Pastur and M. Shcherbina. Universality of the local eigenvalue statis-tics for a class of unitary invariant random matrix ensembles. J. Statist.Phys. , 86(1-2):109–147, 1997.[25] Sumati Surya. Evidence for the continuum in 2d causal set quantumgravity. Classical and Quantum Gravity , 29(13):132001, July 2012.[26] Eugene P. Wigner. On the Distribution of the Roots of Certain Sym-metric Matrices.