Monte Carlo simulations of random non-commutative geometries

Random non-commutative geometries are introduced by integrating over the space of Dirac operators that form a spectral triple with a fixed algebra and Hilbert space. The cases with the simplest types of Clifford algebra are investigated using Monte Carlo simulations to compute the integrals. Various qualitatively different types of behaviour of these random Dirac operators are exhibited. Some features are explained in terms of the theory of random matrices but other phenomena remain mysterious. Some of the models with a quartic action of symmetry-breaking type display a phase transition. Close to the phase transition the spectrum of a typical Dirac operator shows manifold-like behaviour for the eigenvalues below a cut-off scale.


Introduction
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 given by Riemannian manifolds, where the algebra is the algebra of functions on the 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 be non-commutative, leading to a generalisation of the notion of geometry.
A random geometry is a class of geometries G that fluctuates according to a probability measure. In this article, the probability measure is taken to be a constant times e −S(D) dD (1) using a real-valued 'action' S(D) and a standard measure dD on the space G of Dirac operators.
To make this computable, the class of geometries is taken to be the Dirac operators on a fixed finite-dimensional Hilbert space H; thus G is a space of matrices. It turns out that the axioms for D for these finite spectral triples are all linear and so G is a vector space [5]. Therefore one can take dD to be its Lebesgue measure, which is unique up to an overall constant. Thus the object of study is a random matrix model where the matrices are constrained to be Dirac operators.
The algebra A in this construction is also fixed, and is taken to be the algebra of n × n matrices, M (n, C). Spectral triples with this algebra are known as fuzzy spaces [19] and are the simplest type of non-commutative spectral triple. Allowing the algebra to be non-commutative is important because it allows a new type of finite-dimensional approximation to a manifold. Staying within the realm of commutative algebras would lead to the algebra of functions on a finite set of points, which is a lattice approximation to a manifold; simple examples of such random commutative spectral triples are studied in [18,12]. The fuzzy spaces are not lattice approximations and so the study of these is complementary to the study of random lattices. Intuitively one can think of the algebra A as consisting of the functions on a space with a certain minimum wavelength that is determined by n; this picture is known to be accurate for the most-studied example of the fuzzy two-sphere [13].
The purpose of this paper is to study the simplest examples of random fuzzy spaces by computing the statistics of the eigenvalues of D using a Markov Chain Monte Carlo algorithm. The examples are determined by the 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 in which there are p gamma matrices that square to 1 and q gamma matrices that square to −1. The examples studied here are types (1, 0), (0, 1), (2, 0), (1,1), (0, 2) and (0, 3). It would be interesting to go to higher types but these types already exhibit an interesting set of different behaviours, some aspects of which are not yet understood from a theoretical point of view.
The form of the action S(D) has not yet been specified. In this paper it is assumed to be spectral, which means it is of the form for some potential function V , with V ≥ b for some b ∈ R and λ i denoting the eigenvalues of the self-adjoint operator D. The Connes-Chamseddine spectral action [7,8], in which V (x) → 0 as x → ∞, is not suitable because the integral for the partition function does not converge. This is because S(µD), for µ ∈ R, converges to a finite constant as µ → ∞.
In fact, it is necessary that V (x) → ∞ instead. The simplest cases are investigated here, namely S(D) = Tr g 2 D 2 + g 4 D 4 (4) with g 4 > 0, or g 4 = 0 and g 2 > 0. By a simple change of variables in the integral, D → µD, one can assume that either g 4 = 1, or g 4 = 0 and g 2 = 1, so one need only study these cases. Note that one could choose other spectral actions and it is possible that the results obtained here might motivate the study of other choices. When g 4 > 0 and g 2 < 0 the potential has a symmetry-breaking double well form. In random matrix models it is known that this potential leads to a phase transition [9] and so this possibility is investigated here. It is shown here that, at least in some of the random Dirac operator models, there is good numerical evidence for the existence of a phase transition.
The eigenvalue distribution is plotted for several values of the coupling constants and matrix sizes to exhibit the typical behaviours. On a (commutative) Riemannian manifold of dimension m the eigenvalue distribution, or density of states, for the Dirac operator is approximately the density of flat space ρ 0 (λ) = |λ| m−1 (5) when the eigenvalues are large enough. Most of the plots for random noncommutative geometries presented here look nothing like this, except that some of them are approximately constant (m = 1) for some range of eigenvalues.
The exception is close to the phase transition, where the distribution does indeed look very much like (5) for a range of eigenvalues below a highenergy cutoff. Thus as far as the eigenvalues are concerned, the random noncommutative geometries are behaving something like random Riemannian geometries in this regime. The results presented here show that this is a promising area for future investigation. A phase transition to a geometric phase in a multi-matrix model with a rather different Yang-Mills-type action has also been observed in [3,23].
The motivation for the present work comes from the close relation between random geometries and models of quantum gravity, though one does not need to know anything about quantum gravity to understand the results. Most approaches to random geometry have been stimulated by work on quantum gravity but some of them (e.g. dynamical triangulations or Liouville gravity) have found wider application. It is possible that random Dirac operators will also find other applications beside quantum gravity. In quantum gravity the maximum eigenvalue has a ready interpretation as a natural cutoff to gravitational phenomena at the Planck scale. However in a wider context it can be interpreted simply as a finite limit to the resolution to which a geometry is defined.
There are features in common with other models of quantum or random geometry, most notably the existence of a phase transition, which is also evident in dynamical triangulations [2] and lattice simulations [15]. There are however some features of our system that are quite different from those of other models. One point is that the requirement that the action has to have a compact global minimum in the non-compact G is a non-trivial constraint on the model. In other theories of discrete geometry, like causal dynamical triangulations [1] or causal set theory [17], the space of geometries explored in Monte Carlo is a combinatorial space of a finite number of elements and the code is guaranteed to reach equilibrium after a finite, although possibly long, time.
Another interesting feature is the freedom to rescale g 2 and g 4 using the change of variables mentioned above. Monte Carlo simulations show that the rescaling does not change the qualitative features or our system, e.g. relative differences between eigenvalues will remain unchanged. This can be explained by the use of the Lebesgue measure, which does not distinguish any particular scale of energy. This is in marked contrast to systems such as the Ising model, or the causal set model of 2d orders [25].
The non-commutativity also distinguishes this approach from most others. Using a finite-dimensional commutative algebra necessarily leads to a lattice model of quantum geometry defined on a finite set of points. The use of non-commutative geometry allows a more general set of finite-dimensional models where the algebra is an algebra of matrices. Thus one can construct perfectly computable models of random geometry that are not lattice models. Moreover, the standard model of particle physics has a non-commutative geometry using exactly the same framework [4,11], so the hope is it will be easy to combine the two into a unified model of gravity and particle physics.
The technical details of the Dirac operators, observable functions and Monte Carlo method are given in section 2. The results for the action Tr D 2 are given in section 3, where it is explained how the results relate to the standard theory of Gaussian random matrices. Actions including a Tr D 4 term are studied in section 4, with particular attention paid to the symmetrybreaking case which exhibits a phase transition. Section 5 discusses the interpretation of the results. The expansions of the action in terms of the constituent matrices of the Dirac operators are given in detail in appendix A.
2 Technical details

The Dirac operators
The spectral triples considered here are 'real spectral triples', which consist of a finite-dimensional Hilbert space H together with some operators acting in H. These are an algebra A, a chirality operator Γ, an antilinear 'real structure' J and a self-adjoint Dirac operator D. For a given random geometry model H, A, Γ and J are fixed but D is allowed to vary, subject to the axioms of non-commutative geometry.
The axioms are solved in [5] to give explicit forms of the Dirac operator in terms of n × n Hermitian matrices H and n × n anti-Hermitian traceless matrices L according to the formulas below. There are no other constraints on these n × n matrices, so these are the freely-specifiable data for the Dirac operator.
The Dirac operator acts on H = V ⊗ M (n, C), with V = C k the space on which the gamma matrices act. The gamma matrices are assumed to form an irreducible representation of the Clifford algebra, which implies that the chirality operator is trivial for d = p + q odd. The dimension of V is k = 2 d/2 for d even and k = 2 (d−1)/2 for d odd. In the first two cases the sole gamma matrix is just 1 or −i respectively. In the remaining cases the gamma matrices are 2×2 matrices, distinct gamma matrices anti-commuting. As usual, [ · , · ] denotes the commutator and { · , · } the anti-commutator of matrices.
A type (p, q) geometry has a signature s = (q − p) mod 8 which determines some of the characteristics of the spectrum of D. These properties are well-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. The properties can be seen in the Monte Carlo simulations below.
Symmetry For s = 3, 7, if λ is an eigenvalue then so is −λ.
Doubling For s = 2, 3 or 4, each eigenvalue appears with an even multiplicity.
The proof of these is given briefly here. For even s, the chirality operator Γ is non-trivial. It is Hermitian and has eigenvalues ±1. The Dirac operator changes the chirality, DΓ = −ΓD. If v is an eigenvector of eigenvalue λ then Γv is an eigenvector with eigenvalue −λ. As a result, the spectrum of the Dirac operator is symmetric around 0. A similar argument holds for s = 1, 5 using the fact that DJ = −JD so that v and Jv have opposite eigenvalues.
For the doubling property, if s = 2, 3 or 4 then J 2 = −1 (it is 'quaternionic'). Since DJ = JD in these cases, if v is an eigenvector then so is Jv with the same eigenvalue. Moreover, one can check that v and Jv must be linearly independent: suppose the eigenvectors are proportional to each other, i.e., Jv = cv, with c ∈ C, then which is a contradiction.

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 The integral can be approximated as a sum over a discrete ensemble {D j , j = 1, . . . , N }.
so that in the limit taking N → ∞, the average obtained through this discrete sum converges towards the continuum value f . This convergence can be improved by using a Markov Chain Monte Carlo algorithm. In such an algorithm the Dirac operators D j are generated with a probability distribution such that This simplifies the expression for the average and improves convergence by concentrating the sampling on regions which contribute strongly. To generate such an ensemble of Dirac operators the Metropolis-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 next subsection. The proposed operator D p will be accepted as a new part of the chain, If this was the only rule to add new operators to the Markov chain the code would terminate in any sufficiently deep local minimum. To make it possible to escape local minima, the new operator is also accepted if exp(S(D j ) − S(D p )) > p, with p a uniformly distributed random number in [0, 1]. If D p is rejected in both tests then D j+1 = D j . This algorithm ensures a Markov Chain satisfying detailed balance, 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 the initial state D 0 . The states generated before this convergence are not representative of the probability distribution and can not be used to measure observables. We checked that this burn-in process terminated by starting from different initial configurations and monitoring the convergence of the action.
The code is implemented in C++ and all matrix algebra operations use the open source software library Eigen [14].

The Monte Carlo move
To construct a Markov Chain on the space of Dirac operators G = {D}, a move 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 proposed operator can only depend on the current operator D j . The space G is a vector space, so a simple additive move with δD a Dirac operator, will always be ergodic, and as long as δD does not depend on past states the Markov property is also satisfied. As shown in section 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 Dirac operator composed from δH i , δL i . Generate a random n × n matrix R with matrix elements in the complex range [−1 − i, 1 + i] and define 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 acceptance rate of proposed moves is a r = (#accepted moves)/(#proposed moves) 0.5 (where #proposed moves counts all D p generated). At the beginning of a simulation the acceptance rate is tested and l adjusted, larger if the acceptance rate is too large, smaller if the acceptance rate is too small, such that a r 0.5 is satisfied within a tolerance of 1%. The number of attempted Monte Carlo moves is called the Monte Carlo time τ MC . Note that the move for the L i does not preserve the condition that it is trace-free. However since the L i appear only in commutators, the trace decouples and its value does not affect the Dirac operator.

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 of the left and right actions of M . These are written as matrices using the tensor product The Dirac operator D can then be written as a kn 2 × kn 2 matrix that acts on the tensor space V ⊗ C n ⊗ C n . The matrix operations needed in the computer code are matrix multiplication, addition and calculation of eigenvalues. The run time of these grows like O(m 2 ), O(m 2 ) and O(m 3 ) respectively for m × m matrices. Therefore it makes 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 the geometries investigated are collected in appendix A.

Observables
Given a Dirac operator D, the eigenvalues {λ i } can be computed. The two main observables of interest are the j-th eigenvalue, ordering the eigenvalues from lowest to highest f j (D) = λ j (20) and the distribution of eigenvalues at eigenvalue λ Since eigenvalue calculations are computationally expensive, the eigenvalues are only measured every 4n attempted Monte Carlo moves. This improves run time, and reduces the correlation of the measurements. The action S and the acceptance rate of moves are recorded at every step to monitor the algorithm.
At later points it will become useful to measure some additional observables that are computed from the eigenvalues, for example, Tr D 2 . For certain cases it has also proven instructive to examine the non-physical degrees of freedom of the matrices H and L via their eigenvalues.
The average of any observable can be calculated directly from the set of measurements. However to estimate the statistical error on our measurements it is necessary to take the correlation between successive states in with σ λ i the variance of the eigenvalue, τ A,λ i the integrated autocorrelation time of the eigenvalue and M the number of measurements performed [22]. In figure 1, autocorrelations for the simulations of a type (1, 0) geometry with S = Tr D 4 for size n = 5 and n = 15 are shown. The figures show the autocorrelation for both the action and the smallest eigenvalue of D.
The autocorrelation time is determined on the data after the burn-in is completed. In practice the burn-in phase was combined with the adjustment of l. Then τ M C was counted from 5000 moves after the last adjustment to l.
This burn-in and adjustment period takes up most of the simulations. We found that for the eigenvalue distribution and the average eigenvalues, 200 measurements (corresponding to 4n·200 attempted Monte Carlo moves) lead to very good results. To determine the phase transition, 2000 measurements were used to ensure that statistical fluctuations were not mistaken for a phase transition.

Results for D 2 action
In this section the Monte Carlo simulations for the simplest possible action S = Tr D 2 are examined. The one-dimensional Clifford algebras, type (1, 0) and (0, 1) are examined first and the results understood using the standard theory of Gaussian matrix models. After this, some numerical results for the two-and three-dimensional types are shown.
3.1 The simplest cases: type (1, 0) and (0, 1) The n 2 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, 0) case one has eigenvalues while for the (0, 1) case 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 Dirac operator: 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 will be seen later that a peak in the eigenvalue distribution at, or near, 0 is a feature of some other random fuzzy spaces.
The second point is that the spectrum of the (0, 1) case is symmetric about the origin, as λ jk = −λ kj . This is in accordance with its signature s = 1, which means that each Dirac operator has symmetric spectrum. The spectrum of a (1, 0) Dirac operator is typically not symmetric since s = 7 in this case. This means that our simulation gives an eigenvalue distribution that is not exactly symmetric, though it will eventually converge to a symmetric distribution as the Monte Carlo time increases.
For the (1, 0) case, using the simplified action (47) one can transform the integral over the Dirac operator into an integral over the Hermitian matrix H.
The (0, 1) case is similar, but one has to take into account the fact that the integration over Dirac operators is an integration over traceless matrices L. Using (51) gives An example of average ordered eigenvalues generated by the Monte Carlo simulation is shown in figure 2.
These random matrix models are close to the Gaussian Hermitian matrix model [20,6], which has the similar action with integration over all Hermitian matrices. A standard technique in random matrix models is to calculate the joint probability density for the eigenvalues µ j . The formula is [21] P (µ 1 , . . . , µ n ) = C exp − S(µ 1 , . . . , µ n ) The terms with the differences of eigenvalues result from the Jacobian for the change of variables from the matrix elements to the eigenvalues. Since this term is small when two eigenvalues are close, this results in the phenomenon of the repulsion of eigenvalues.
The matrix M can be split into traceless and trace parts, and these are statistically independent. It follows that expectation values in the (0, 1) model can be calculated as expectation values of observables in the Gaussian Hermitian matrix model that are independent of the trace of M . This is done by writing M = −iL. The action in the (1, 0) model transforms directly to the Gaussian Hermitian matrix model by rescaling the trace by a factor √ 2. The transformation is A standard result (the Wigner semicircle law [26]) is that the analogue of the eigenvalue distribution (21) for the Gaussian Hermitian matrix model converges as n → ∞ to the density of states with A = 1. In our simulations using actions S (1,0) and S (0,1) we find that the semicircle law is also a good approximation for the eigenvalues of H and L. It is already well-satisfied for n = 5 and improves at higher n, as shown in figure 3. The reason for this is that in the Gaussian Hermitian matrix model, the variable 1 n Tr M is normally-distributed with variance 1/(4n 2 ), and so adjusting the eigenvalues with a fixed multiple of this makes no difference to the density of states in the limit n → ∞.
Another standard result from random matrix theory is that the correlation between different fixed eigenvalues of M vanishes as n → ∞. Thus for large n, the joint probability distribution away from the diagonal µ 1 = µ 2 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 the semicircle law as a convolution, with a correction for the behaviour of the correlations on the diagonal. This is shown as follows. Let f (λ) be an observable for a random fuzzy space, with λ = µ 1 ± µ 2 . Then assuming a product probability distribution, one has with density of states for the Dirac operator the convolution   which is an elliptic integral. This integral is the same for type (1, 0) and (0, 1), since σ(µ) = σ(−µ). The Monte Carlo simulation of the eigenvalue density 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 for the (0, 1) case. The correction to the product probability density gives a contribution only near the diagonal µ 1 = µ 2 . The approximate form is an additional contribution to f of [24] − sin 2 πn(µ 1 − µ 2 )σ((µ 1 + µ 2 )/2) 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 scales as 1/n.

Higher types
While the spectra for geometries with one-dimensional Clifford algebra are easy to understand, those with a two-dimensional Clifford algebra are less straightforward. The average eigenvalues and the eigenvalue distributions are shown 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 three types are symmetric about the origin and the third one exhibits eigenvalue doubling, all in accordance with the properties for s = 6, 0 and 2 derived in section (2.1).
The action Tr D 2 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) actions previously analysed. In particular, these matrices are statistically independent. The eigenvalues of the H i , L i are still approximated well by the semicircle law (32) with A = 1/2. However, the main difference in analysing the two-dimensional cases is that the eigenvalues of the Dirac operator are not simply related to the eigenvalues of L i , H i .
For the case (0, 2) the multiplicity of eigenvalue 0 is at least 2n, as shown by examining the Dirac operator using a basis so that γ 1 γ 2 = diag(i, −i). The Dirac operator acts on the space C 2 ⊗ M (n, C). All v ⊗ m in this space for which have eigenvalue 0. Picking a basis on v one can choose v 1 = 1, v 2 = 0 and v 1 = 0, v 2 = 1. There will then be n linearly independent matrices m that commute with L 1 + iL 2 and n that commute with L 1 − iL 2 , hence 2n 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 shows there is eigenvalue repulsion for this Dirac operator also, though we do not have a theoretical understanding of this phenomenon. The types (2, 0) and (1, 1) also have a feature at the origin. The density of eigenvalues is sharply lower in a narrow dip at the origin and there is a somewhat wider upward spike around this. This is shown for the (2, 0) case in figure 7, which zooms in on a region around eigenvalue 0. The gap in the middle is further evidence of eigenvalue repulsion, this time between eigenvalue λ and the opposite eigenvalue −λ that is required by the symmetry of the spectrum of D about 0.
The numerical results also indicate that the range of the eigenvalues remains unchanged under the change of matrix size, and the distribution becomes smoother, appearing to converge to a smooth limiting distribution in the same way as for random matrices. Another similar feature is that for        larger n the fluctuation of each individual eigenvalue becomes smaller. This can be seen in figure 8. The leftmost bump in each plot is the smallest eigenvalue while the rightmost bump is the largest, the eigenvalues in between were chosen to be symmetric, and include the central most eigenvalues. The eigenvalue distribution for the type (0, 3) case is plotted in figure  9. This appears to be smooth at the origin, similar to the (1, 0) case. The common property of these cases is that the Dirac operators do not have a symmetric spectrum. Thus a small eigenvalue λ does not have to be close to any other eigenvalue. There is nothing special about the origin, and in particular, the eigenvalue repulsion hypothesis does not lead to any special behaviour here.

Results for actions with D term
The Tr D 4 term in the action leads to interactions between the L i , H i that compose the Dirac operator. An extreme case of this is for type (0, 3), in which a four point interaction of all four matrices H, L 1 , L 2 and L 3 is present. These terms make it harder to understand the system analytically, however for the simulations they are no obstacle.
The simple action S(D) = Tr D 4 leads to behavior very similar to that for the action Tr D 2 . This is shown in figure 10. Some characteristics, like the shoulders for type (1, 1) and (2, 0) are more pronounced, but the overall shape is quite similar.
Combining the two terms together gives the action  For positive values of g 2 the behaviour of the numerical simulations is somewhere between the Tr D 2 case and the Tr D 4 and does not show qualitatively new features. However when g 2 is negative this is a symmetry-breaking potential with two minima, shown in figure 11. The question of how the eigenvalues behave in this case is interesting and a variety of behaviours is exhibited depending on the type of the gamma matrices. This is shown in figure 12. The different types are described here for values of g 2 decreasing from 0.
(1, 0) Two peaks start to form at around g 2 = −3 then grow and separate sharply between g 2 = −3 and g 2 = −3.5, leaving the centre of the distribution empty. Since the Dirac operator is not symmetric, the Monte Carlo simulation can and does settle in just one of the peaks, though one expects that the Markov Chain would eventually explore both peaks equally given a long enough run.
(0, 1) Two peaks form at around g 2 = −3 and grow slowly and steadily. The central part of the distribution remains. One can understand this from the fact that the eigenvalues of L settle into two peaks, and since it is traceless, the favoured configuration has the same number of eigenvalues in each peak. The differences between eigenvalues of L in the same minimum 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 2 and grow until the central part of the distribution vanishes suddenly between g 2 = −2.5 and −3.
(1, 1) Two peaks develop at small g 2 and grow until the central part of the distribution vanishes suddenly between g 2 = −2 and −2.5.
(0, 2) This is the most mysterious case. Two slight peaks develop but the eigenvalues do not separate into two peaks for the whole of the range of g 2 tested. Instead some further substructure to the eigenvalue distribution develops.
(0, 3) This is similar to the (0, 1) case, with the sharp change occurring between g 2 = −4 and g 2 = −4.5. In figure 12f the Markov Chain for g 2 = −4.5 has to a certain degree explored both minima.
These descriptions can be compared to the plots of the order parameter ∂ log Z ∂g 2 = Tr D 2 (42) and the autocorrelation time, which is usually expected to increase near a continuous phase transition due to the long-range order ('critical slowing down'). These are plotted in figure 13. One can see that there is good evidence for a phase transition for the types (1, 0), (2, 0), (1, 1) and (0, 3). In these plots, the order parameter changes gradient at around the values of g 2 described above, and the autocorrelation time of Tr D 2 has a peak around this value also. It is difficult to see any clear signal from the plots for types (0, 1) and (0, 2). It is remarkable that the types for which the evidence for a phase transition is clearest are the ones where the Dirac operator contains an anti-commutator with a Hermitian matrix H. Unlike the L matrices, the trace of H appears to play a crucial role. The observable measures the strength of Tr H, calculated as a square so that positive and negative 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, 0), there are two matrices H 1 and H 2 , and the F for both combined is In both cases, F = 1 if the matrices are pure trace. The plots show that Tr H develops a large expectation value at the phase transition. In the case of (2, 0) the Monte Carlo data used for the plots developed a preference for Tr H 1 rather than Tr H 2 but this is of no significance due to the rotational symmetry between the two 1 . The sum of squares is the correct rotationallyinvariant observable.

Conclusion
A model of random geometry has been presented here as random Dirac operators in non-commutative geometry. The integrals can be interpreted as multi-matrix models but with a new type of observables, namely the eigenvalues of the Dirac operator. The one-dimensional cases can be understood using theoretical results from random matrices but the higher-dimensional types are not so easy and will require further study to obtain analytic results. Numerical results have been presented showing various phenomena that depend strongly on the type of gamma matrices used, particularly whether the spectrum of a Dirac operator for that signature is symmetric about the origin. From the numerical results it is clear that some features are similar to the properties of the eigenvalues of random matrices: the eigenvalue distributions appear to converge in the large n limit and the dispersion of individual ordered eigenvalues decreases; also there is some evidence of a degree of eigenvalue repulsion at the origin.
The most interesting results are for the quartic action with negative g 2 , so that the potential is of symmetry-breaking type. For some types, the eigenvalue spectrum changes suddenly when g 2 reaches a critical negative value and the observable Tr D 2 is a good order parameter for this transition. This is taken as a strong indication that a sharp phase transition would occur in the large n limit. The types where this transition is clear are those where the Dirac operator contains at least one term involving an anti-commutator with a random hermitian matrix H. Then the trace of H develops a large expectation value, becoming the dominant contribution to H after the transition. This can't happen with commutator terms as the trace of the random matrix decouples in this instance.
For generic g 2 , the eigenvalue distribution of D is not a good approximation to the behaviour for the Dirac operator on any fixed (commutative) Riemannian manifold (5), except that one could possibly argue that the distribution is approximately constant for some ranges (e.g. figure 10(f), which looks like a one-dimensional manifold). The exception to this is near the phase transition, where the curves in figure 12 do appear to have the right power-law behaviour. Two of these distributions are highlighted in figure 15, showing the types (1, 1) and (2, 0) at values of g 2 just below the value for the phase transition.
These are compared with the eigenvalue distribution for the fuzzy sphere from [5], shown in figure 15(c). This is a type (1, 3) spectral triple, having signature s = 2, and has exactly the same spectrum as the Dirac operator on the Riemannian round S 2 but with a maximum eigenvalue cut-off and fermion doubling. The distributions are remarkably similar, the main differences being the gap at the origin, the size of which depends on the distance from the phase transition, and the fact that the fuzzy sphere has exactly integer eigenvalues with multiplicity, due to its rotational symmetry. The feature that is common to the plots is the approximately linear increase of the eigenvalue density with eigenvalue that is characteristic of Riemannian manifolds of dimension two, i.e., m = 2 in (5). The simulations show that decreasing g 2 further increases the gap in the middle of the spectrum whereas the middle of the spectrum fills up for values of g 2 greater than the critical value. These results are somewhat preliminary and we have not yet carried out a systematic study of the phase transition.
The study of random non-commutative geometries gives an insight into the closely-related problem of the quantization of this fascinating theory. A quantized non-commutative space is a potential candidate for a quantum theory of gravitational interactions and will allow a better understanding of fundamental interactions. Independent of these physical applications, it also is an interesting modification of the well-known matrix models, introducing non-trivial interactions and observables among several matrices. The results reported here will be a basis for further investigations into the phase transition, the continuum limit, and other possible actions on this space of geometries.

A Dirac operators for the fuzzy geometries we 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 independent of trace part of these matrices. The bracketing convention for the trace of an expression is Tr AB ≡ Tr(AB) and Tr A n ≡ Tr(A n ), but Tr A + B ≡ (Tr A) + B.