Matrix Model for Riemann Zeta via its Local Factors

We propose the construction of an ensemble of unitary random matrices (UMM) for the Riemann zeta function. Our approach to this problem is ‘piecemeal’, in the sense that we consider each factor in the Euler product representation of the zeta function to first construct a UMM for each prime p. We are able to use its phase space description to write the partition function as the trace of an operator that acts on a subspace of square-integrable functions on the p-adic line. This suggests a Berry-Keating type Hamiltonian. We combine the data from all primes to propose a Hamiltonian and a matrix model for the Riemann zeta function. arghya.chattopadhyay@gmail.com parikshitdutta@yahoo.co.in suvankar.hri@gmail.com d.ghoshal@gmail.com

apparent lack of pattern in the occurrence of the prime numbers. Montgomery [1] was the first to quantify this randomness through the two-point correlation function of the (non-trivial) zeroes of the zeta function. His serendipitous encounter with Dyson that helped in relating this to the behaviour of a random Gaussian ensemble of large unitary matrices is the stuff of a well known anecdote [2]. Since then extensive numerical work [3] and analysis of higher correlators [4][5][6][7] (see also [8]) have established this connection on a firm footing, and even extended to the closely related Dirichlet L-functions [9]. All these are consistent with the suggestion of Hilbert and Pólya that these zeroes are in the spectrum of an unbounded, self-adjoint operator.
Berry and Keating observed that the fluctuating part of the prime counting function, which is also related to the 'height' of the non-trivial zeroes, resembles the Gutzwiller trace formula that gives, asymptotically, the fluctuating part of the energy distribution of a quantum Hamiltonian in terms of a sum over classically periodic orbits, along with corrections that depend on the deviations around them [10,11]. They, however, need to use the Euler product form of the zeta function in a domain where it is really not applicable. Be that as it may, the prime numbers correspond to the primitive periodic orbits and their powers to the repetitions. They propose that the non-trivial zeroes are the eigenvalues of a quantum dynamical system the classical Hamiltonian of which is H = xp. The authors caution that this cannot quite be, since the spectrum of this Hamiltonian, being related by a canonical transformation to the inverted harmonic oscillator, is not discrete. Moreover, the system is integrable, and hence cannot exhibit chaotic behaviour, which is an expected property of the Hamiltonian. Several possible resolutions were explored to overcome this difficulty. Among them are putting boundary or periodicity conditions on the coordinates and/or momenta. The dilatation symmetry, under which x and p scale so as to cancel each other, was also discussed. For further work on this model, with different boundary conditions see Ref. [12] Connes has also explored a Hamiltonian description for a much wider class of zeta functions [13,14]. This approach involve noncommutative geometry, trace formulas and p-adic number fields. It is, however, rather difficult to access for physicists. Ref. [15] has a comprehensive review and references of various physics related approaches to the Riemann zeta function.
Given an infinite set of distinguished points on a curve, e.g., the non-trivial zeroes of the Riemann zeta on the critical line, one can construct a random ensemble of unitary matrices (UMM) after conformally mapping the points on the unit circle. Using this approach, a unitary one-plaquette model was constructed for the symmetric zeta function (which is an entire function that has zeroes coinciding with the Riemann zeta) by two of us (PD and SD) in Ref. [16]. The parameters of the model are determined in terms of the Li (also known as the Keiper-Li) coefficients. Extending on earlier work [17], a phase space description of the model was developed. In this picture, an eigenvalue is the coordinate and the number of boxes in the Young diagram corresponding to a representation is the momentum. The phase space distribution function was computed and found to have a 'broom-like' structure with spikes at the zeroes of the symmetric zeta function. While the dynamical system is expected to be a many fermion system, we could not get any further insight from this phase space.
In this paper, we shall follow a similar line of approach, however, we shall break it in parts. We shall take up each prime factor in the Euler representation of the Riemann zeta function (the so called local zeta function at a fixed prime) separately and construct an ensemble of unitary matrices (UMM), the eigenvalues of which correspond to its poles that lie on the imaginary axis. It turns out that the phase space density, or more accurately, the fluctuating part of it, can be expressed as the trace of an operator. The operator in question is a generalised Vladimirov derivative that acts on a subspace of square integrable functions supported on the p-adic integers. This suggests that the natural notion of distance in the momentum space is ultrametric, defined in terms of the p-adic norm. The phase space picture also leads us to arrive at a Hamiltonian, which has the form H = xp. Interestingly its action is restricted to a subspace of the full domain -this is reminiscent of a resolution to the shortcomings of this type of Hamiltonian in [10,11]. In this approach, even the classical phase space description is in terms of operators on a Hilbert space, as in the Koopman-von Neumann formalism [18,19]. More specifically, the Hilbert space is a subspace spanned by the Kozyrev wavelets [20] (an analogue of the generalised Haar wavelets) on the p-adic line. The restriction is to the subspace supported on the compact subset Z p of Q p with the translation parameter set to zero and scaling restricted so as to stay within the support.
Finally, we attempt to combine the UMM of all primes (local data) to construct a random matrix model for the Riemann zeta function. (The part of the data that correspond to the trivial zeroes is related to GUE for gamma function, which could be thought of as a local zeta function at infinite place.) Not unexpectedly, this leads to a divergent expression for the coefficients of the model. We propose a renormalisation procedure to deal with this divergence-this requires the introduction of one real parameter, which may be taken to be half to match with the critical line. This shift from the imaginary axis (relevant to the local zeta) to the critical line can be seen to arise from a similarity transformation of the Hamiltonian, although that too does not seem to fix the value of the renormalisation parameter. We end with a discussion on the nature of the phase space, and the Wigner function on it.
It should be emphasised that our objective is not to try to prove the Riemann hypothesis. Rather, our effort is to find a Hamiltonian through the GUE, by first constructing one for each prime factor, and then combining them. There are issues that will require better understanding, and we shall remark on these at appropriate places. The organisation in the following is as outlined in the table of contents. We begin in Section 2 by recalling a few facts about the zeta function and p-adic analysis that we shall use subsequently. The GUE for the local zeta factors are constructed in Section 3, and their phase space picture developed in Section 4. These results are combined in Section 5 in an attempt to build a GUE for the Riemann zeta function. We also discuss our proposal for a renormalisation and attempt to understand this large phase space. The resulting renormalised coefficients for the UMM have been compared In Appendix A with that for the symmetric zeta function constructed earlier in Ref. [16].

Zeta functions and p-adic analysis
In this section, we shall recall a few facts about the Riemann zeta function as well as the p-adic number field and complex valued functions defined on it. We do not attempt to be comprehensive or even list the most important issues, but restrict to only those aspects that will be of importance for us in what follows.

Riemann and local zeta functions
The infinite sum over the natural numbers for positive integer values of s ≥ 2, was introduced by Euler, extended for Re(s) > 1 by Chebyshev and analytically continued to the complex s-plane as a meromorphic function by Riemann. This function has a simple pole at the origin, and vanishes for all negative even integers, which are known as its trivial zeroes. More interestingly, it is believed to have an infinite number of non-trivial zeroes. Riemann also showed that ζ(s) satisfies the reflection identity ζ(s) = 2 s π s−1 sin from which one concludes that the non-trivial zeroes must all be in the critical strip 0 ≤ Re(s) ≤ 1. However, he conjectured a much stronger result, namely that the nontrivial zeroes lie on the critical line Re(s) = 1 2 . This is the celebrated Riemann hypothesis that is among the most famous unsolved problems in mathematics. Associated to it is the symmetric zeta-function which is an entire function that satisfies ξ(s) = ξ(1 − s), and has zeroes only at the non-trivial zeroes of ζ(s), at s = γ m = 1 2 + it m . Li proved that the Riemann hypothesis is equivalent to the non-negativity of the following sequence of real numbers [21] (see also [22,23]) called the Li coefficients (sometimes the Keiper-Li coefficients). Using Eq. (2.6), these can be expressed as These numbers appear in the parameters that define the matrix model associated with the zeta function [16].
Euler also expressed the sum Eq. (2.1) as an infinite product over the prime numbers. It is natural to expect, therefore, that the distribution of the zeroes of the zeta-function, i.e., the locations of the t m , to be related to the distribution of the primes. The factors in Eq. (2.6) are called the local zeta-function at the prime p The local zeta function does not have any zero, rather it has equally spaced simple poles at s = 2πi ln p n for n ∈ Z on the vertical line Re(s) = 0.
The prime counting function J(x) = p ∈ primes n∈Z + 1 n Θ(x − p n ) is a monotonically increasing function that jumps by 1 n at every power of prime numbers p n [24] (see also [25,26] for more accessible introductions for physicists). It can be written as ln t is the logarithmic integral function. Riemann showed that ln ζ(s) and J(x) are related by Mellin transformation as follows.
The function J(x) = J (x) + J(x), has an average part J(x) and a fluctuating part J(x), the latter may be expressed in terms of zeroes of the zeta function [24] as A related function that we shall find use for later is the summatory von Mangoldt function or the Chebyshev counting function In the above we have an defined  p (x) = n∈N Θ(x − p n ), a local counting function at a prime p, for which Notice that  p (x) jumps by 1 at the powers p n of the fixed prime p. The following identity proves very useful.
In deriving this, we change variables as t = 1−s and x = e q in the inverse Mellin transform of Eq. (2.12) and pick up the poles at integer multiples of 2πi/ ln p in the t-plane.

A few results from p-adic analysis
The motivation for the local zeta function at prime p comes from the p-adic field Q p and analysis on it. We shall need to use some results from there. Let us recapitulate a few facts about the field Q p and some aspects of p-adic analysis. More details are available at many places, e.g. Refs. [27][28][29].
First, let us fix a prime p. The p-adic norm of a rational number m/n is m n p = p ord p (n)−ordp(m) , where ord p is the highest power of p that divides its argument. Starting with the rationals, the field Q p is the Cauchy completion obtained by including the limits of all converging sequences (upto an appropriately defined notion of equivalence) with respect to the p-adic norm. This norm has the ultrametric property leading to a stronger than usual triangle inequality |ξ − ξ ′ | p ≤ max (|ξ| p , |ξ ′ | p ). Since the real numbers are obtained by the same procedure, but using the absolute value norm, Q p 's are similar to R. Indeed, upto equivalence, the absolute value norm and the p-adic norms are the only possible norms on the set of rationals. An element ξ ∈ Q p admits a Laurent series expansion in p: where N ∈ Z and ξ n ∈ {0, 1, · · · , p − 1}, but 1 ξ 0 = 0. The series above is convergent in the p-adic norm. The subset Z p = {ξ ∈ Q p : |ξ| p ≤ 1}, i.e., elements with norm less than equal to 1, is the compact subring of p-adic integers. One can define functions on Q p , which could be valued in any field. In particular, we shall be interested in complex valued functions. For example, χ : Q p → C, defined as has the property χ (p) (ξ + ξ ′ ) = χ (p) (ξ)χ (p) (ξ ′ ), and is called an additive character of Q p . The totally disconnected nature of the p-adic space means that the (complex valued) continuous functions on Q p are locally constant. One such function is the indicator function and similarly for other open sets. In fact, the indicator functions on open sets are the p-adic analogue of the Gaussian function on the real line, in the sense that these retain their form after a Fourier transformation. There is a translationally invariant Haar measure dξ on Q p . The discrete valuation of p-adic numbers reduces the problem of integration to evaluation of sums. The measure is normalised by The following integral, which is easy to evaluate, will be important for us in what follows.
On the other hand, the usual notion of derivative does not work because of the totally disconnected topology of Q p . The Vladimirov derivative [28], therefore, is defined as an integral kernel which is in fact well-defined for any α ∈ C [30]. Furthermore, The additive character Eq. (2.15) and the indicator function Eq. (2.16) are the ingredients for a set of functions defined by Kozyrev in [20] for n ∈ Z, m ∈ Q p /Z p and j ∈ {1, 2, 3, · · · , p − 1}. They provide an orthonormal basis for Moreover, they satisfy with eigenvalue p α(1−n) . As shown in Ref. [31], this enhances the scaling symmetry of the p-adic wavelets to a larger SL(2,R) group . It is possible to restrict the wavelets to a compact subset 2 , say Z p , just as the Haar wavelets on the real line may be restricted to the finite interval [0, 1]. This would mean restricting scaling of the mother wavelet by the parameter n only in negative integers (only contractions), and the translation parameter m can take a finite number of values that satisfy |m| p = p −n (the parameters of translation are in Q p /Z p ). Restricted thus, the wavelets span a subspace L 2 (Z p ) ⊂ L 2 (Q p ).
Before we close this section, let us get back to the zeta function. The Euler product formula shows that the local zeta functions Eq. (2.7) are naturally associated with the prime numbers, and by extension to Q p . On the other hand, R is the only other completion of the rationals. The symmetric zeta function Eq. (2.3) motivates the following definition called the adelic zeta function that involves the zeta functions for Q p for all primes and for R. It satisfies ζ A (s) = ζ A (1 − s). It is interesting to note that that is, the LHS is the Mellin transform of the Gaussian function. On the other hand, the local zeta function is also a Mellin transform of the p-adic analogue of the Gaussian [32]. Therefore, in a sense π − s 2 Γ s 2 = ζ R (s) ≡ ζ ∞ (s) is the local factor related to R. It is also called the zeta function at infinite place.
The adelic ring A = R × p Q p consists of elements x = (x ∞ , x 2 , x 3 , x 5 , · · · ), where x ∞ ∈ R, x p ∈ Q p , and all but a finite number of x p are p-adic integers (|x| p ≤ 1). Elements can be added and multiplied component by component. An example is n = (n, n, n, · · · ), n ∈ Z, for which the adelic relation |n| p∈primes |n| p = 1, relating the product of all p-adic norms together with the absolute value norm of an integer, holds. We shall not need the adèle, but would like to point out that that the term in Eq. (2.23) is motivated by it.

Unitary matrix model for local zeta
A unitary matrix model for the Riemann zeta function was constructed in [16]. To be precise, the matrix model constructed there was for the logarithm of the symmetric zeta function. In this section, we shall construct a matrix model corresponding to the local zeta functions. We shall largely follow the procedure developed there, however, we shall work directly with the function ζ p (s) and not ln ζ p (s). More precisely, our objective will be to construct an ensemble of unitary matrices, the eigenvalue distribution of which coincide with the distribution of the poles of a local zeta function. Since the eigenvalues of unitary matrices lie on the unit circle, we first map the poles of Eq. (2.7) on the imaginary s-axis to the unit circle on the z-plane by the conformal map The partition function of a generic one plaquette unitary matrix model is where the real numbers β n are the parameters of the theory. In the following, we shall find the set of parameters {β n } corresponding to the poles of ζ p (s).

Eigenvalue distribution
An eigenvalue analysis of this model has been done in [33,34]. Consider an ensemble of N × N matrices, the partition function can be rewritten by change of variables to the eigenvalues The additional factor, known as the Vandermonde determinant, is due to the Jacobian of the transformation. In the limit N → ∞, let us define the continuous variables x and θ(x) as to arrive at the partition function where, the Jacobian has been incorporated in the effective action S eff [θ].
In the large N limit, the dominant contribution comes from extrema of S eff [θ] which is determined by the saddle point equation where, Thus for a given set of the parameters β n , one has to solve saddle point equation Eq. (3.5) to find the density of eigenvalues ρ(θ).

The resolvent
It is difficult, in general, to solve the integral equation Eq. (3.5) to determine the distribution of the eigenvalues. It is easier to work in terms of the resolvant that is analytic both inside, as well as outside, the unit circle. Expanding in Taylor series around z = 0 and ∞ in the two regions Since 1 N TrU k lies between −1 and +1 for all k, the Taylor series is convergent. The resolvant should be such that R(0) = 1 and R(z → ∞) = 0 as well as All these are satisfied by Eq. (3.8). In the large N limit the resolvent satisfies an algebraic (quadratic) equation, rather than an integral equation, that is far easier to solve. In order to obtain this equation one uses the fact that the Haar measure is invariant under a more general transformation U → Ue tA , where A is an anti-symmetric matrix and t is a real parameter. The partition function Eq. (3.4), being an integral over all unitary matrices, is also invariant under this transformation. This yields the identity (3.10) which is known as the Dyson-Schwinger equation.
A solution to this equation is given by where the choice of the sign + or − is for |z| less, respectively greater, than 1, and the function F (z) is The analyticity properties of R(z) inside the unit circle depend on the form of F (z), which determines the phase structure of the model (see [35] for more details). In the 'no-gapped phase' (in which the eigenvalues are distributed over the entire unit circle), only the first terms in the parentheses survive, therefore, F (z) is a perfect square and from which one can find β n if the resolvant is known. A resolvant function that satisfies all the requirements and carries the information about the poles of the local zeta function, is The coefficients (definining the no-gap phase of the matrix model) may in principle be determined from the above, however, we shall not need their explicit forms.

Eigenvalue density from the resolvent
The eigenvalue density can be obtained from real part of the resolvent [36] where the second expression is equivalent to the first from Eq. (3.11). The imaginary part of R(z), on the other hand, gives the derivative of the potential that appears on the RHS of the saddle point equation Eq. (3.5). Thus the eigenvalue density function and the potential for the matrix model of the local zeta function are in 0 < θ < 2π. These functions extremise the effective action We shall work with this definition of the effective action of the matrix model for the local zeta function at the prime p.
It would appropriate to comment on the eigenvalue distribution Eq. (3.18) at this point. At the point θ = 0, the function ρ(θ) becomes negative, contrary to what is expected of a density. We believe that this is an artefact of the conformal mapping from the s-to the z-plane. The number of poles of the local zeta functions (similarly, the number of zeroes of the adelic or symmetric zeta function) increases as Im(s) → ±∞ with Re(s) fixed to 0 (or 1/2 respectively). Consequently, the number of discretely located peaks in the eigenvalue density increases as we approach θ = 0 on the unit circle in z-plane. It is possible to remove the point θ = 0 from the domain of definition of ρ(θ), but this will be at the cost of the normalisibility condition on ρ(θ), which in turn is related to R(0) = 1. This is a technical issue which needs a better understanding. The solution to it may lie in developing a phase space picture for Hermitian matrix models. However, if we accept it as such, the matrix model technique will allow us a phase space description. Nevertheless the resolvent is a valid solution to the Hilbert transform arising from extremising the action. Therefore, it can be used to calculate the moments of the specific model.

Matrix model in the phase space
The partition function (and the effective action) of the unitary matrix model can be expressed as an integral over a phase space, in which the angular coordinates are augmented by their conjuguate momenta. We shall review this phase space description, developed in Refs. [16,17] with the motivation that phase space picture may suggest a natural Hamiltonian.
The starting point is the partition function that can alternatively be expressed as a sum over representations of U(N), which are characterized by Young diagrams. In the large N limit, the contributions from a class of diagrams dominate [37]. It is possible to express the partition function in terms of the phase space of a system of free fermions. In this picture the coordinates and momenta are the eigenvalues θ i and the number of hooks h i of the right-most boxes in i-th row of a Young diagram.

Sum over representations
The exponential in the partition function Eq. (3.2) of a single plaquette model may be expanded as where χ R (C(k)) is the character of the conjugacy class C(k) of the permutation group S K (K = n nk n ) and R denotes an irreducible representation of U(N). Finally, using the orthogonality of the characters as the partition function. The sum over representations of U(N) can be traded for a sum over the Young diagrams, with a maximum of N rows and an arbitrary numbers of boxes in each row, as long as the number of boxes in a particular row does not exceed that in a row preceding it. Let λ j be the number of boxes in j-th row, and j λ j = K.
Note that the total number of boxes is the same as the order of the permutation group S K . The characters of the conjugacy classes are determined recursively by the Frobenius formula (see, e.g. [38][39][40]). In the large limit N → ∞, we introduce the continuous variables (where the decreasing sequence of numbers h j = λ j + N − j denote the hook length) in terms of which the partition function, in the large N limit, is where S eff is an effective action defined in terms of h(x) and the parameters k = N 2 k ′ .
The dominant contribution at large N comes from those representations which extremize S eff . These representations are characterised by the so called Young diagram density as required of a density.
In the generic unitary one plaquette model, it is difficult to find u(h) by extremizing the effective action S. This is due to the fact that explicit expressions for the characters of a permutation group in the presence of non-trivial cycles are not available.
The character is determined from the Frobenius formula as . In order to select these coefficients, let us introduce N complex variables (z 1 , z 2 , · · · z N ) to write the character as an integral where the contours are taken around the origin. We should mention that in the above, z i are auxiliary variables, introduced as variables of integration. The notation, however, is chosen in anticipation of the fact that they will soon be identified with the eigenvalues of the U(N) matrices.

Large N saddle
We now exponentiate the terms ( z i ) k i , take the large N limit where we use the continuous functions of x = i/N, and replace the sum over i by an integral over x, to rewrite the characters as In spite of the notation, S χ [h(x), k ′ ] is not really an action.
In the large N limit, the dominant contribution to the characters arise from a saddle point that correspond to a particular distribution of z(x). The conditions to determine the extremum of S χ are where we have introduced ρ(z) = ∂x/∂z in the second line. Since the contours in Eq. (4.10)) go around the origin, we can take z = e iθ and write the above as where now, with an abuse of notation, ρ(θ) = ∂x ∂θ = ie iθ ρ(z) naturally. The imaginary part of the above is exactly the same as the eigenvalue equation Eq. (3.5) with nk ′ n = β n dθ ′ ρ(θ ′ )e inθ ′ , which is the equation of motion for k ′ . This identification is actually consistent and the auxiliary variables which appear in the Frobenius formula can indeed be thought of as the eigenvalues of the unitary matrices under consideration. Not only do they satisfy the same equation as the eigenvalues, the partition function written in terms of these variables exactly matches with the partition function written in terms of the eigenvalues. Thus, in the large N limit, the numbers of cycles k ′ n are fixed in terms of the parameters β n . The density ρ(θ) above is the same as the density of eigenvalues. For more details, we refer to [16].
The real part of Eq. (4.14) h(θ) = 1 2 + n β n cos nθ (4.16) is an algebraic equation that relates the number of boxes in a Young diagram to the eigenvalues of unitary matrices. In terms of θ and ρ(θ), S χ can be expressed as

Many-body wave function
The results in Sections 4.1 and 4.2, in particular Eqs. (4.5), (4.7) and (4.12), allows one to write the partition function of a unitary matrix model as The above is related to the wavefunction of N fermions The many body wavefunction Ψ(z, h) = e −N 2 Sχ leads to the 'fermion density' function Ω(z, h) = Ψ † Ψ.

Simplifying the local model
Let us try to simplify some of the expressions. The potential term in the matrix model involves an infinite sum n β n Tr(U n ) (and its conjugate). Using the fact that the coefficients β n were determined from the resolvant, and the expression Eq. (3.14) for the latter, we can perform the sum With this we can write where, to get to the partition function in the last line, we have performed the integrals over the momenta h i , which yielded δ(θ i + θ ′ i ), using which we eliminated one set of integrals over the coordinates θ.

Density function in the phase space
Both the eigenvalue distribution and the distribution of the Young diagram (of at least a particular class of unitary matrix models) can be obtained from a single distribution function [16,17]. This function, which is one in some region of the (θ, h) space, and zero outside it, may be identified with the phase space distribution of N free fermions in one dimension.
Let us define a complex phase space distribution function When integrated over h, the real part of Ω gives the eigenvalue density ρ(θ), while the imaginary part gives the derivative of the potential V ′ (θ). Since these are even and odd functions of θ, respectively, the real and imaginary part of Ω(θ, h) are also, respectively, even and odd. Hence, when we integrate Ω(θ, h) over θ only the real part contributes to give the Young diagram distribution function u(h) above. Although the distribution function is simple, the shape (or topology) of the imaginary part of Ω(θ, h), which contains information about the derivative of the potential, could be complicated, but is fixed for a given matrix model. Whereas the topology of its real part will change as we go from one phase to another.

Phase space for the local model
We have seen that in the no-gapped phase, ρ (p) (θ) = 1 + n β (p) n cos nθ gives the density of the eigenvalues on the unit circle of the matrix model for the local zeta function. Together with V (p)′ (θ) = n β (p) n sin nθ (where the coefficients β (p) p are to be determined from the resolvant function Eq. (3.14)) one can define a complexified eigenvalue density ̺ (p) (θ) = 2R (p) (θ) − 1 = ρ (p) (θ) + iV (p)′ (θ). This function carries the information of the poles of Eq. (2.7), the local zeta function.
Recall the resolvant for the local zeta function proposed in Eq. (3.14). This gives where the ability to express it as a sum suggests that we may be able to make sense of it as the trace of some operator. We shall now argue that this is indeed possible, albeit not on a Hilbert space that is very familiar to physicists. Rather the operator is on a subspace of L 2 (Q p ), the space of square integrable complex valued functions on the p-adic line. In fact a comparison of the above with Eq. (2.22) makes it immediately obvious that the operator is related to the (generalized) Vladimirov derivative D −i cot θ 2 defined in Eq. (2.18). We write where the trace is expressed in terms of the (normalised) Kozyrev wavelets ψ (p) −n+1,0,1 defined in Eq. (2.19) for n ∈ N. With this choice, the translation parameter m gets restricted to a finite set, as discussed at the end of Section 2. However, we can consistently restrict to the subspace in which m = 0. The parameter j may also be restricted to unity without any loss of generality. The set of Kozyrev wavelet functions, which are eigenfunctions of the Vladimirov derivative form a basis for complex valued, mean-zero, square integrable functions on the p-adic space Q p [20]. Our restricted set of functions span a subspace H (p) − ⊂ L 2 (Z p ) ⊂ L 2 (Q p ) of complex valued, mean-zero, square integrable functions supported on the compact set Z p ⊂ Q p .
It should be emphasised that notwithstanding the notation, we are dealing with a classical system, therefore, the states are classical states. This is in the sense of Koopmanvon Neuman classical mechanics which uses the Hilbert space formalism [18,19]. As in quantum mechanics, classical states are vectors in a Hilbert space and physical observables are Hermitian operators, however, unlike in quantum theory, the operators commute. The dynamical equation is the Liouville equation satisfied by the 'wavefunction'. For further details, we refer to the articles cited above.
We would like to arrive at the density ̺ (p) (θ) from a density in phase space. This is possible thanks to the integral Eq. (2.17) from which we can read the phase space density Ω (p) (θ, h) where θ ∈ [0, 2π] ⊂ R and h takes values in Z ⊂ Z p ⊂ Q p . The momentum density function u (p) (h) Eq. (4.23) is obtained by integrating Ω (p) (θ, h) over the eigenvalues 3 Thus modulo singularities of the second term, we find that the momentum distribution function is a constant, which means that the momenta are equally spaced. This is consistent with the spectrum of the operator D α (with α = −i cot θ 2 ) in Eq. (2.22).

Hamiltonian for the local model
The fluctuating part of the phase space density Eq. (4.27), without the measure factor, is This satisfies the Liouville equation for the Hamiltonian H (p) = n ln p cot θ 2 , which is a function of generalized coordinate X (p) = cot θ 2 and momentum P (p) ∼ ln D (p) ∼ n ln p. We see that δΩ (p) = e −iH (p) , although in reality any function of the Hamiltonian will satisfy the Liouville equation for the density.
This Hamiltonian is of the form H = xp as prescribed in Berry and Keating [10,11], but expressed in the Koopman-von Neumann operator formalism for a classical dynamical system. Moreover, the metric on the integer valued momenta is the non-archimedean padic metric. The Berry-Keating Hamiltonian has a continuous spectrum and is integrable, while any Hamiltonian related to the Riemann zeta function is expected to be chaotic. Furthermore, the spectrum of the corresponding quantum system is not discrete. In order to address these problems, the authors of Ref. [11] suggested imposing restrictions on the accessible region of the phase space. It is interesting to note that for the trace in Eq. (4.25) to be meaningful, that is, for the operator to be trace class, we have a natural truncation to the subspace H (p) − of L 2 (Q p ). The action of the Vladimirov derivative D (p) naturally restricted to this subspace. In fact, the subspace is defined by the span of a subset of orthonormal eigenvectors of this operator.

Matrix model in the large phase space
We shall now attempt to construct a matrix model for the Riemann zeta function, by combining the data from the local zeta functions for all primes. Naturally this involves taking a product (of the infinite number) of partition functions we obtained for the primes individually. Not unexpectedly, some of the results are divergent. We shall adopt the usual approach employed by physicists in treating divergences and try to extract a finite answer through a process of renormalisation.

Combining the local information
In order to define the model we need to specify the potential, which in turn is determined by the coefficients {β m }. Naturally Z = p Z p and S eff = p S (p) eff . Let us, therefore, consider the sum these coefficients for the local models. However, to conform to the standard convention, we include an additional factor of ln p in the sum.
where, the expression in the last line is schematic: We have defined a large Hilbert space in which we have scaled the vectors in the p-th Hilbert space by a factor of √ ln p. The generalised Vladimirov operator D α on the large Hilbert space is defined by the equation above. We may take it to be where we see that it has a 'block-diagonal' form, although in fact, since the basis consists of its eigenvectors, it is actually diagonal. The eigenvalues of the this operator are p nα (n ∈ Z), and the corresponding eigenvectors are Ψ −n+1,0,1 ≡ ψ i.e., it has the wavelet that scales by p n at the p-th place, and the basic mother wavelets at all other places. These eigenvectors have an adelic flavour. We need only those vectors with n ∈ N ⊂ Z.

Now it is not hard to see that
where ψ(x) is the summatory von Mangoldt function defined in Eq. (2.11). This expression, however, is formal due to singularities in the integral over θ.

Renormalization of the parameters
Let us try to make sense of the divergent parameters β m by 'renormalisation'. To this end, we 'renormalise' the Kozyrev wavelet states for each prime p by multiplying with a constant. To be precise, let us replace ψ −n+1,0,1 so that the renormalised parameters β ren m are A convenient choice for the constants is c (p) n = p −nµ/2 for some real µ. This is equivalent to a shift in the exponent of the Valdimirov derivative D (p) . Notice that a multiplicative renormalisation of the eigenstates of D (p) is equivalent to a redefinition of the operator itself because of the special functional form of its eigenvalues. As an aside, we note that this renormalisation is not unique, as one can multiply this choice by an infinite series starting with 1. This is equivalent to multiplying the shifted Vladimirov derivative by an infinite series f D (p) = 1 + a 1 D (p) µ 1 + a 2 D (p) µ 2 + · · · , as long as the additional terms lead to convergent integrals. This will of course change the form of the β ren m from its simplest form above. Therefore, we shall use the simplest renormalizsed form.
As a result of the renormalisation, the parameters now take the following form where we have used Eq. (2.11). As expected this has contribution from the trivial zeroes of the Riemann zeta function (the last set of terms) as well as its pole (the first term 1 in the parentheses). We shall discuss the contribution of the trivial zeroes shorty (in Section 5.3), but for now, in the following let us concentrate on the contributions from the non-trivial zeroes only (the second sum). Splitting these zeroes into real and imaginary parts, γ m = σ m + it m , this piece is The renormalised parameters β ren m should be finite for the matrix model to be defined. The requirement for this is that the integrals be convergent. The condition for that is met if σ m ≤ µ for all σ m . A priori not any value of µ is distinguished, however, saturation of the inequality (namely σ m = µ for all m) is certainly more elegant. In this case, thanks to the reflection symmetry, the unique choice is µ = 1 2 .

Matrix model for gamma function and the trivial zeroes
From the resolvant, one can find the parameters of the UMM where θ parametrise the unit circle |z| = 1, as in Eq. (3.1). In order to simplify this, let us use the Weierstrass product representation Since the Euler-Mascheroni constant is γ = lim k→∞ k n=1 1 k − ln k , we may trade the divergent sum for an infinite factor, to write the parameters as x −2n−1 + 1 2 lim n→∞ ln n π (5.10) It is apparent from this that, modulo an infinite constant, the trivial zeroes in Eq. (5.6) are removed by the matrix model for ζ R (s). The divergent constant is probably an artefact of this construction, as one of would not expect any singularity here.
We can also see this if we start with the product of the local zeta ζ p (s) for all primes, together with ζ R (s), the zeta function at infinite place where γ m are the non-trivial zeroes of the Riemann zeta function. Hence, Using Eq. (5.9) and the results in Section 3.2, together with the identity 1 which is essentially the coefficients found in Section 5.2 corresponding to the non-trivial zeroes, but without the renormalisation. In fact, this is what one would expect from Eq. (5.11).

Towards a Hamiltonian in the large phase space
We saw in Section 4.7 that (the fluctuating part of) the phase space density of the unitary matrix model for the local zeta function at a prime p suggests an xp-type Hamiltonian H (p) = cot θ 2 ln D (p) . For the matrix model corresponding to the Riemann zeta function, it is natural to propose the Hamiltonian where the definition of ln D is similar to that is Eq. (5.3). This is only schematic, as it may lead to divergences. Indeed, we have seen that the parameters β m of the combined matrix model are divergent, and needs a renormalisation before the can make sense. It would not be a surprise if after combining the data for the local models, the full Hamiltonian only makes sense as a quantum operator. Let us start with the quantum HamiltonianĤ =xp, wherex andp satisfy the canonical commutation relation [x,p] = i (in units = 1). One may define another Hamiltonian H ′ by a similarity transformĤ BothĤ andĤ ′ lead to the same classical Hamiltonian. Thus, for µ ∈ R, one gets a class of Hamiltonians that differ by a complex shift inx. The action of this similarity transformation on the operator e −iĤ is as follows We shall now propose a construction for the operators on the Hilbert space spanned by the Kozyrev wavelets. From its spectrum, we see that the momentum operator is the logarithm of the Vladimirov derivative, but we need to identify the generalised coordinate.
To this end, let us note that the recursive construction of the wavelet basis defines the raising and lowering operators J as shown in [31]. Consider a Schwinger realisation of this algebra using a pair of (bosonic) creation and annihilation operators a † (p) , a (p) and b † (p) , b (p) (satisfying [a (p) , a † (p) ] = 1 and Since the momentum operatorP (p) is the difference of the two number operators N, it is natural to identify the coordinateX (p) with the difference in the phase operator Φ (such that [Φ, N] = i) The construction of the phase operator has a long history, and is somewhat indirect (see [41] and references therein, where the construction of analytic functions of Φ has also been discussed). The coordinate assumes the continuous values cot θ 2 , and the momentum has a discrete spectrum. The Hilbert space that is relevant for us, is spanned by only the positive values of ln D (p) , therefore, it should suffice to 'freeze' one of the oscillators.

Wigner function in the phase space
What is the phase space? Although a proper understanding of its geometry and analytical structure will need further investigation, let us take a first step in this direction. The Hilbert space H (p) − (in the Koopman-von Neumann formalism) for the matrix model at a prime p, consists of a suitably truncated subspace of complex valued, mean-zero, square integrable functions on a compact subspace of Q p . The large Hilbert space is a consistent truncation of functions with these properties on a compact subspace of the product space p Q p .
For a local matrix model, the coordinate and momentum take values in the set R × Z, where the first factor is from the eigenvalues and the second corresponds to the number of boxes in a Young diagram in a representation. In expressing (the fluctuating part of) the phase space density as the trace of an operator, however, the topology in the momentum space was determined by the ultrametric norm of the p-adic numbers. Thus, the local phase space should be taken as R × Z p . When we combine these local constituents to define the large phase space, the Hamiltonian is the sum of the prime factors acting on vectors in the large Hilbert space. Hence the proposed large phase space is where in the second expression, we have taken the coordinate as the 'diagonal' in ⊗ p R, as it is common for all p. This space is reminiscent of the adelic space.
In order to understand this large phase space better, let us go back to the small phase space for which our description has been classical. We shall try to understand the corresponding quantum system. From the spectrum of the momentum, we can label the eigenstates of the momentum operatorP (p) = ln D (p) (see Section 4.7) as |n (p) ≡ |ψ  It is worth emphasising that the sum is only over positive integers, and the same is true of the inner product Using the identity Eq. (2.13), this can alternatively be written as In what follows, we shall sometime omit the label p to avoid cluttering up formulas.
Either of the forms is consistent with the periodic nature of (the difference of) the phase operator. Equivalently, with the fact that the eigenvalues of the UMM, which are the poles of the local zeta function, repeat periodically along the imaginary axis. We also note that the above is the Schrödinger wavefunction Ψ (p) x (x ′ ) in the position representation. Let us construct the Wigner distribution function W (p) x ′ (x, q) for the position eigenstate labelled by x ′ , which is defined in terms of the operatorρ (p) where we have used the wavefunction Eq. (5.20), performed the integrals over q and one of the κ's, and used the identity Eq. (2.13) and the definition of the local counting function  p (x). The above form requires that q > κ/2, and since 0 ≤ κ < ∞, this is the expression in q ≥ 0. Before we proceed further it would be prudent to verify that the same Wigner function is obtained starting with its momentum space definition From Eq. (5. 19) we see that, upto a factor of ln p, scalar products between the position and momentum eigenstates are supported only at positive integer points on the real line, that is It is easy to check that this will indeed lead to the same Wigner function as in Eq. (5.22) after we substitute the scalar product above and redefine of the dummy variable of integration in the expression Eq. (5.23). The integral of the Wigner function over momentum should give the probability density in the coordinate space. Let us check that this is the case.
We see the periodicity of the eigenvalue distribution, related to the poles of the local zeta function.
In the above, we have omitted the terms involving the limit Λ → ∞, as the contributions of these terms to an integral with any well-behaved function will vanish. Therefore, the fluctuating part of the eigenvalue density is We see that with x − x ′ = cot θ 2 , Eq. (5.24) matches exactly with the above, provided we identify the divergent factors.
On the other hand, after integrating the Wigner function over position, we find which is related to the momentum distribution (upto a divergent factor). The divergence seems to suggest that the Wigner function ought to be 'renormalised', which could perhaps be achieved by a similarity transform.
In the large phase space, the position eigenstates in H − , when expanded in the momentum basis, after the conventional scaling as in Eq. where n ′ (p ′ ) |n (p) = δ n (p) n ′ (p ′ ) δ pp ′ is the normalization of the momentum states. The inner product of these states that follows from the above is The above can be re-expressed in two different ways. Either we may introduce integrals for each p, as in Eq. (5.21), or we may use only one integral over κ, leading to the two following expressions which ought to be equivalent.
In the above, ψ(x) is the summatory von Mangoldt function Eq. (2.11) that is related to the counting of prime numbers. We also note that since x is chosen to be the same for all primes, that is along the diagonal in ⊗ p R, one would not expect to have any periodicity since the periods along the different directions are incommensurate.
Hence we see that when we combine the information from all primes to get the Wigner function in the large phase space, there are two ways of parametrising it. One way is to think of the coordinates and momenta as infinite component vectors (x, q) = (x, x, · · · ), (q (2) , q (3) , · · · ) . The other is to simply take (x, q) ∈ R 2 . The coordinates for all primes have been identified, thus justifying the variable x. For the canonically conjugate momentum operator P = p ln D (p) on the other hand, its eigenvalues (and hence expectation values) are real numbers. Thus, by the correspondence principle, the momentum may be taken to be real valued. In the following, we shall work with the latter choice, that is, take the classical phase space of the matrix model for the Riemann zeta function to be R 2 . There may be a deeper connection between the two descriptions. In this context, it may be relevant to mention that Ref. [42] shows a canonical bijection between the spaces of locally constant complex valued functions on ⊗ p Q n p and distributions on R n . It will be worth exploring if this result (for n = 1) explains the equivalence of the two ways of parametrising the large phase space.
Let us start with the Wigner function on the large phase space for the operatorρ x ′ +iµ = p ln pρ (p) x ′ +iµ , in which we have included the conventional additional factor of ln p to be consistent with Eq. (5.1), and at the same time have shifted x ′ → x ′ + iµ in anticipation.
When integrated over x, we get which involves the square of the derivative of ψ(x). Since the integrated Wigner function is the probability density, the conclusion to draw from here is that e −µq dψ ( e q ) dq is the wavefunction (in momentum space). On the other hand, when integrating over q, if we take dψ(e q ) dq = p ln p n∈N δ(q − n ln p) (5.31) where we have used γ m = 1 2 + it m . This is an example of what are known as trace formulas [43].
We see from its definition Eq. (5.29) that dq W x ′ (x, q) = |Ψ x ′ (x)| 2 in the square of the wavefunction Eq. (5.34) (with x ′ shifted by iµ). Hence from the trace formula we find that the wavefunction is In the above, a Dirac delta function with a complex argument is to be understood in terms of the residue of the function that appears in the integrand with it. One would expect the terms involving the Γ-function to be accounted for by the contribution from ζ R . Although we have not verified it, we think that it would be natural for the wavefunction for the adelic zeta function Eq. (2.23) to involve only tm δ(x − x ′ − t m ) that correspond to the non-trivial zeroes. Working backwards, one can see that this would be equivalent to the following inner product where, the last expression is true for µ = 1 2 . This may be realised through the following similarity transform (defined with a star-product appropriate to the phase space formalism) W x ′ +iµ (x, q) = e µq ⋆ W x ′ (x, q) ⋆ e −µq (5.35) of the naive Wigner function one would have assumed from the individual prime factors (together with the contributions from the Gamma function). This is consistent with the similarity transform Eq. (5.13) on the Hamiltonian, and the fact that the support of the eigenvalues of the UMM for the Riemann zeta function is shifted from those of the local zeta functions. The square of the wavefunction Eq. (5.34) is again of the same form, upto an infinite factor, as in the case of the Wigner function for a fixed p.
Their approach does not have any explicit connection with prime factorisation or the adelic structure. We also refer to [45] for a critique and [46] for a reply to it. (Let us also mention Ref. [15] and references therein for physicists approaches to various aspects of the Riemann zeta function.) We approach this problem using the phase space description of the unitary matrix model (UMM), a random ensemble of unitary matrices, the eigenvalue distribution of which is known to exhibit the same statistical features as the Riemann zeroes [1,3,5].
In fact a phase space description for a UMM for the Riemann zeta function was developed in [16]. In this paper, we split this analysis into two parts, first we have broken it to its 'prime parts'. We construct UMMs for each prime factor of the Riemann zeta function, that is a UMM for each of the local zeta function at a prime p. Each of these zeta functions has regularly spaced poles on the imaginary axis. We have then combined these UMMs to define a UMM for the Riemann zeta function. We encounter divergences in this process.
In particular, the coefficients that define the UMM involve an infinite sum, and diverge in this process. We propose a 'renormalisation' to extract a finite answer to properly define the UMM. The coefficients that we get through this process match with those found in Ref. [16], which constructed a UMM for the symmetric zeta function. As in the Berry-Keating proposal, the Hamiltonian we arrive at has the xp form, however, a direct comparison is difficult for many reasons. First, we work in the Koopmanvon Neumann formalism in which a classical system has a Hilbert space description. Secondly, our approach involves the local fields Q p of p-adic numbers and complex valued functions over these spaces, since we work with the local zeta functions. Thirdly, in trying to make sense of the divergent answer, we have to do a similarity transform on this basic Hamiltonian. These features could be advantageous, since the xp Hamiltonian cannot be right in its original form. In fact, the truncation to a subspace of the Hilbert space of complex valued square integrable functions on Q p that we have is reminiscent of the truncation of the allowed region of the classical phase space in Ref. [10]. We also recall the work of Connes [13] in which p-adic numbers and the adelic ring play essential roles.
We have only initiated a programme here, but have not resolved all the issues satisfactorily. There are several points that will need further clarity and resolution. First among these is the density function for the eigenvalues, which exhibits an essential singularity at z = 1 (equivalently θ = 0). Formally, this makes the eigenvalue density non-positive. The issue has its root in the fact that there are an infinite number of poles of the local zeta function (respectively, the zeroes of the adelic or symmetric zeta function) beyond any finite value of Im(s) as it approaches ±i∞, with Re(s) fixed at 0 or 1/2, as appropriate. After the conformal mapping of the critical line in the s-plane to the unit circle, needed for description in terms of unitary matrices, this results in an accumulation point in the finite part of the z-plane. This problem would not arise if one could develop a phase space description for the partition function of a random ensemble of Hermitian matrices, however, presently the phase space picture (in terms of eigenvalues and representations) is available only for a unitary ensemble. Secondly, only the fluctuating part of the corresponding phase space density can be related to the Hamiltonian. This, however, is not different from the situation in Ref. [10]. Thirdly, the geometrical and other properties of the phase space need a thorough investigation, so that the nature of the Hamiltonian and its spectrum could be understood better and directly, independent of the matrix model. Related to this are the issues of self-adjointness or otherwise of these operators. It would be desirable to analyse the quantum mechanical problem. This is possible in principle since we have a concrete proposal with an explicit Hamiltonian and the Hilbert space on which it acts as an operator. However, the existing approaches to quantum mechanics on local fields (see e.g., [47,48]) are not directly applicable, since in our proposal the coordinates are real valued while the momenta take values in Z p or its adelic counterpart, or, to be precise the topology (distances and measures) in this part is determined by these ultrametric spaces. This would fit with a symplectic structure, if one exists, in the adelic space. Finally, the local zeta functions at primes p (as well as that at the 'infinite place', i.e., R) have an infinite sequence of poles repeating with regular periodicities. After combining to form the large matrix model, these poles disappear, and the special points are the zeroes of the Riemann zeta function. It would be instructive to understand how this transformation takes place.
We close by emphasising the main ingredients of our approach once again. First, we construct a one-plaquette unitary matrix model given a distribution of special points (poles or zeroes of a function), and make use of its phase space formulation. We apply this to the prime factors in the Euler product representation of the Riemann zeta function, and express the phase space density as the trace of an operator on a suitable Hilbert space of functions on a subspace of the p-adic number field. This operator leads to the suggested Hamiltonian for the zeta function. It is important to remember that this is not an attempt to prove the Riemann hypothesis, but to find a Hamiltonian related to its zeroes.

A Comparison with the matrix model for symmetric zeta
A one-plaquette unitary matrix model for the symmetric zeta function Eq. (2.3) was constructed in Ref. [16]. The parameters β m were found to be given by the Li coefficients Eq. (2.4). Let us compare these with our renormalised β ren m coefficients obtained in Section 5.2. Notice that for obvious reasons the conformal map Eq. (3.1) from the z-to the s-plane that we have used is different from the one in [16]. In order to facilitate comparison one may start by a translation so that the poles of the local zeta function lie on the line Re(s) = s 0 and then map the shifted line to the unit circle. This changes the expression of the coefficients in Eq. This should be compared with the coefficients in [16] β sym m = −  Apart from a constant factor the only difference of this with the renormalized coefficients Appendix A is that the symmetric zeta function appear instead of the Riemann zeta function, however, this discrepancy is due to not incorporating the contribution of the matrix model of ζ R discussed in Section 5.3.