# Quantum simulation of Riemann-Hurwitz function

###### Abstract

We propose a simple realization of a quantum simulator of the Riemann-Hurwitz (RH) function based on a truncation of its Dirichlet representation. We synthesize a nearest-neighbour-interaction Hamiltonian, satisfying the property that the temporal evolution of the autocorrelation function of an initial bare state of the Hamiltonian reproduces the RH function along the line of the complex plane, with . The tight-binding Hamiltonian with engineered hopping rates and site energies can be implemented in a variety of physical systems, including trapped ion systems and optical waveguide arrays. The proposed method is scalable, which means that the simulation can be in principle arbitrarily accurate. Practical limitations of the suggested scheme, arising from a finite number of lattice sites and from decoherence, are briefly discussed.

###### pacs:

03.65.Db, 42.50.-p, 03.67.Ac,## I Introduction

Riemann function has attracted the interest of mathematicians and physisists for quite a long time. This function plays a key role in number theory, in particular in the distribution of prime numbers Euler (); Schleichbook (), and it is the basis of one of the most fundamental mathematical conjectures, the Riemann hypothesis. In physics, Riemann’s function is found to be related to a wide variety of different physical areas and phenomena, ranging from classical mechanics to statistical and quantum physics. For an extensive review of the topic, see Riemann-RMP (). Several works have highlighted the close connections among the Riemann hypothesis, random matrix theory and the physics of classical and quantum chaos (see, for instance, Berry1 (); Berry2 (); PNAS (); billard (); Bourgade () and references therein). In statistical physics, the Riemann
function can be seen as the partition function of a quantum gas, called the Bose Riemann Gas
with the prime numbers labeling the eigenstates Julia (); uffa1 (); uffa2 (). In quantum mechanics, several attempts have been
made to introduce a quantum system whose spectrum is associated with the Riemann function. In particular,
following the original Hilbert-Pólya conjecture great efforts have been devoted to propose quantum Hamiltonians whose bound states coincide with the zeros of the function Riemann-RMP (); Berry2 (); kepalle (); kepalle1 (); kepalle2 (). Riemann zeros also enter into phenomena like Bose-Einstein condensation or the cosmic microwave background. For example, a Bose-Einstein condensate could be used, at least in principle, to factorise numbers and to calculate the prime
factors Holthaus ().

In a recent work Schleich (), Mack and collaborators showed that a generalization of the Riemann function, introduced by Hurwitz Hurwitz (), can be retrieved from the autocorrelation function of a quantum state propagating in an anharmonic oscillator potential. The connection to this Riemann-Hurwitz (RH) function is made through its Dirichlet representation, and requires that the quantum system be initially prepared into a thermal phase state, the so-called Riemann state Schleich (). Unfortunately, such a thermal phase state is not the mixed thermal state of the anharmonic oscillator, and system preparation might be a nontrivial issue involving a coherent state superposition.

In this work, we propose a simple method to simulate the RH function by using an dimensional quantum system, in which system preparation in a thermal phase state is readily provided by excitation of a bare state of the system. Interestingly, our finite-dimensional system is described by a nearest-neighbour-interaction Hamiltonian, which can be implemented in a variety of quantum or classical systems, such as trapped ions or optical waveguide arrays with engineered hopping rates.

The paper is organized as follows. In Sec.II, we briefly review the connection between the correlation function of a specially-prepared quantum system and the Dirichlet representation of the RH function, introduced in Ref.Schleich (). In Sec.III we introduce a rather general technique to synthesize a finite-dimensional tridiagonal Hamiltonian satisfying the property that the temporal evolution of the autocorrelation function of an initial bare state of the system reproduces the RH function along the complex axis , with . In Sec.IV we briefly present possible physical implementations of the nearest-neighbor Riemannian Hamiltonian, based on trapped ions or optical waveguide array systems. Finally, in the concluding section V we discuss the main limitations of our quantum simulator scheme and briefly suggest possible extensions of our results.

## Ii Basic idea

We consider an -dimensional Hilbert space, in which is an orthonormal basis and the Hamiltonian describing our system in such a basis is the matrix . Following the idea in Schleich (), we assume that the Hamiltonian has a logarithmic energy spectrum,

(1) |

where and . Initially, we prepare our system in state , which is assumed to be easy to prepare experimentally. Then this state evolves according to

(2) |

where is the th eigenstate of with an eigenvalue ,

(3) |

and . Let us now suppose that the amplitudes can be set to

(4) |

where is a normalization factor. If the above assumptions are fulfilled, after some simple algebra, the correlation function, i.e. the probability amplitude of state at time , is given by

(5) |

At this point, we can compare this expression with the Dirichlet representation of the RH function, given by

(6) |

where we have set . The Dirichlet series on the right hand side of Eq. (6) converges for , as shown in Fig. 1. An accurate estimate of the RH function can be obtained by truncating the series up to a certain order . This is precisely what our quantum simulator does, see Eq. (5). A detailed discussion of the error arising from truncation will be discussed in Sec.V. Note that the ordinary Riemann function is obtained for . This simple connection allows for the simulation of the RH function, as far as the following requirements are met: (i) the Hamiltonian has a logarithmic spectrum (1); (ii) the initial amplitudes of the eigenstates of obey Eq. (4); (iii) there exists a simple experimental setup, providing such Hamiltonian. In Fig. 1 we show the allowed part of the complex plane, where our approach is applicable. Unfortunately, the most interesting part, , cannot be accessed. We will return to this point in Sec.V. In the next section we propose a practical approach, which allows to meet the above requirements for the Hamiltonian.

## Iii Construction of the Hamiltonian

In this section we propose a rather general technique to synthesize a tridiagonal Hamiltonian matrix with the required spectrum (1), which should also meet the condition (4). To this aim, we start with the diagonal matrix defined by the logarithmic eigenvalues (1), and we make a sequence of similarity transformations (basis changes), which do not change the eigenvalues but modify the eigenvectors of the Hamiltonian matrix. Let us indicate by the -dimensional diagonal matrix

(7) |

where the energies are assumed to follow the logarithmic dependence (1). Next, let us make a basis transformation, by using an orthogonal matrix ,

(8) |

where we choose the elements of the first column of equal to the amplitudes in Eq. (4). The other columns are constructed as orthonormal vectors of the orthogonal complement. This task has infinitely many solutions, which, however, are identical regarding our purposes. Finally, we want to transform the matrix into a tridiagonal form. This can be achieved, for instance, by using a tridiagonalization procedure, based on Householder reflections HR (). The described algorithm leads to a tridiagonal Hamiltonian , which has the desired logarithmic spectrum, and also fulfils condition (4). For the sake of clarity, we present an example of the described procedure for . We choose the parameters in the simulation to be and . In this case, the matrix , according to Eq. (1), is

(9) |

The transformation matrix has the approximate numerical value

(10) |

where the elements of the first column are given by Eq. (4), and the next columns are constructed such that all vector-columns are orthonormal, which leads to an orthogonal matrix, . Using this transformation matrix, we obtain for the Hamiltonian in the new basis

(11) |

Finally, we proceed to the tridiagonalization of . This is achieved by using a sequence of Householder transformations,

(12) |

where

(13) |

and

(14a) | |||

(14b) | |||

(14c) |

The explicit numerical value for the Hamiltonian is

(15) |

where we have neglected the minus signs of the negative off-diagonal elements, since they can be removed by only a simple phase transformation in the amplitudes. Obtained in this way, the Hamiltonian has the desired logarithmic spectrum, and if we diagonalize it,

(16a) | |||

(16b) |

we notice that the first row of is equal to the first column of , and is equal to the numbers in Eq. (4). The fact that and have the same first row follows from the property that the Householder transformations in Eq. (12), with vectors (14), do not change the first row of the matrix, which diagonalizes the Hamiltonian. Hence the matrices and , which diagonalize respectively and , have the same first row. The matrix is in fact the matrix that connects the two bases and , which leads to the possibility to express state as a coherent superposition of the eigenstates of ,

(17) |

where have the desired values (4). In such way we have fulfilled all of the necessary restrictions for the Hamiltonian. The described procedure can be applied for an arbitrary dimension , which means that in principle we can simulate RH’s function with an arbitrary accuracy.

As an example, in Fig. 2 we compare the simulation of the RH function, obtained by numerical integration of the Schrödinger equation, with the properly normalized RH function. We see that even when a small dimension for the Hamiltonian is used, , the simulation performs very well, provided that is not too close to the boundary of convergence of the Dirichlet series. The error of the approximation depends on the values of the parameters and . As approaches unity, if is kept fixed, the error will become larger, as shown in the lower panel of Fig. 2. We discuss more on that issue in Sec.V. In the next section, we propose several implementations in particular physical systems.

## Iv Physical Implementation

The tridiagonal Hamiltonian introduced in the previous section arises in many different areas of physics. In particular, tight-binding Hamiltonians with engineered hopping rates and site energies have been extensively investigated in connection to the problem of perfect quantum transfer in quantum networks, and physical implementations based on e.g. spin chain models have been suggested (see, for instance, SC1 (); SC2 (); SC3 (); SC4 () and references therein). Hence we do not aim here to point out all the possible realizations of the proposed scheme. We will specifically focus on two possible physical implementations, namely in trapped ions and in optics.

A first possible realization uses a linear chain of trapped ions. From an experimental point of view, trapped ions are one of the most advanced systems for quantum computation Trapped Ions Comp () and quantum simulation Trapped Ions Sim (). For the purpose of this work, we need to be able to experimentally construct a tridiagonal symmetric interaction matrix, with a control over each element of the Hamiltonian. This can be achieved, for instance, by implementing a nearest neighbour spin-spin interaction

(18) |

where are the Pauli matrices, are the spin-spin coupling constants, and are individual couplings to an external (fictitious) magnetic field. This Hamiltonian commutes with the excitation number operator, which means that the Hilbert space factorizes into subspaces with different number of excitations. If we consider the single-excitation subspace, it is easy to show that, in this basis, the Hamiltonian can be written as a tridiagonal matrix Kay ()

(19) |

where we neglect a common irrelevant term in all the diagonal elements. Traditionally, Hamiltonian (18) is implemented by using optical spin-dependent forces, induced by laser beams Trapped Ions Sim (). In this approach, off-resonant coupling to the motional sidebands is used in order to produce effective spin-spin interactions. Recently it has been shown that a general coupling matrix can be implemented by using such a setup Monroe (). As an example, for a spin chain made of sites the normalized values of couplings and local magnetic fields , needed to realize the logarithmic spectrum (1), are given by Eq. (15). Note that the corresponding magnetic fields show a slow monotonic increase, whereas the hopping rates change in the range of about a mean value.

A second possible physical implementation of the tridiagonal Hamiltonian is provided by light transport in a chain of evanescently-coupled optical waveguides with engineered coupling constants and propagation constants (see, for instance, cazzo () and references therein). In an array of coupled optical waveguides, it is well-known that transport of discretized light is governed by a set of coupled-mode equations described by a tridiagonal Hamiltonian of the same form as Eq. (19), where is the total number of waveguides, are the propagation constants of guided modes in the various waveguides of the array and are the coupling (hopping) constants between adjacent waveguides cazzo (); reviewuff (). In order to independently engineer the coupling and propagation constants , it is convenient to circularly bend the axis of the waveguide array. A schematic of the optical waveguide setting that realizes the Hamiltonian is shown in Fig. 3. Axis bending basically introduces a mismatch between the propagation constants of adjacent waveguides due to different geometric paths. Such mismatch can be properly tailored by controlling the radius of curvature of waveguides and the waveguide separation in the plane of axis bending (see, for instance, reviewuff ()). On the other hand, the hopping rates between adjacent waveguides can be tailored by a proper choice of the waveguide separation. To appreciate the role of axis bending and to properly design the waveguide array, we recall that propagation of discretized light in the bent waveguide array of Fig. 3 is governed by the coupled-mode equations cazzo (); reviewuff ():

(20) |

), where is the amplitude of light waves trapped in the -th waveguide, is the hopping rate between waveguides and (with ), is the effective index of the fundamental mode of the individual waveguide with a straight optical axis, and is the propagation constant detuning of waveguide induced by bending of the optical axis , with the -coordinate of waveguide [see Fig. 3]. In the above Eq. (20), is the refractive index of the substrate and with the optical wavelength. For waveguide separation , the hopping rate is given to an excellent accuracy by the exponential law , with and some constants depending on waveguide fabrication parameters that can be experimentally determined. The detuning between waveguide and is determined by the -coordinate difference [see Fig. 3], and thus can be geometrically controlled by selecting the proper value of in the interval . To obtain the correlation function in the photonic structure of Fig. 3, the left-boundary waveguide should be excited at the input plane by a light beam, e.g. using butt-coupling by an optical fiber. The correlation RH function is then simply retrieved by monitoring the amplitude of light that remains trapped in the waveguide of the lattice along the propagation distance . It is worth noting that, with waveguide arrays manufactured by the femtosecond laser writing technology writing (), a careful control of hopping rates and waveguide bending is nowadays possible. For example, tridiagonal Hamiltonian matrices with inhomogeneous hopping rates with up to a few tens have been recently demonstrated in femtosecond laser written wavegude chains palle (). Hence the waveguide design suggested in Fig. 3 is expected to be feasible with current waveguide fabrication technologies.

## V Conclusions and discussion

In this paper we have proposed an experimentally feasible method to simulate the Riemann-Hurwitz function by the autocorrelation function of a finite-dimensional quantum system initially prepared in a bare state. The proposal is based on an engineered one-dimensional tight-binding Hamiltonian, with specific conditions for the eigenvalues and eigenvectors, which can be implemented by e.g. trapped ion systems or evanescently-coupled optical waveguide arrays.
As compared to the recent proposal of Ref. Schleich () and based on the realization of an anharmonic quantum oscillator, our quantum analog simulator uses a more feasible tight-binding Hamiltonian and enables a very simple preparation of the system in a thermal (Riemannian) phase state. Like the proposal of Ref. Schleich (), our analog simulator uses the Dirichlet representation of the RH function, and so it is regrettably not suited to simulate the more interesting behavior of the function along the critical line , where the nontrivial zeros of the Riemann function are located. A desirable extension of our results would be to synthesize a tight-binding Hamiltonian, which can simulate the behavior of the Riemann function along such a line. The starting point of such goal could be provided by the Riemann-Siegel formula Schleichbook (), which gives an approximation of the function by a sum of two finite Dirichlet series which is valid along the critical line. It is envisaged that such two finite sums could be implemented by means of a finite-dimensional block-diagonal Hamiltonian, where each block simulates independently each of the two series in the Riemann-Siegel formula. In this proposal, however, the measurement of the autocorrelation function might be a more challenging task, and further investigation is required.

On a practical level, it should be mentioned that our proposal suffers from some limitations, which are similar to those discussed for other systems (see e.g. Ref. Schleich ()). A first issue is related to the truncation of the Dirichlet series, which makes our approach unsuited to simulate the RH function near the line , where the convergence of the Dirichlet series becomes very slow and thus extremely large values of would be required. In principle, the algorithm presented in Sec.II works well for small as well as large values of (e.g. up to 1000), however in any physical implementation is limited to a relatively small value because of technological issues. For example, in the photonic realization discussed in Sec.IV an upper limit to can be roughly estimated to be of a few tens (e.g. ). As is increased, tolerances in the values of and that ensure the logarithmic spectrum (1) of eigenvalues become more stringent and unrealistic. The limited number of available in any realistic implementation of poses a restriction on the domain of the RH function that we can access in an experiment owing to the slow convergence of the series (6) as . For a given value of , the error introduced by truncation of the series (6) remains at an acceptable small level provided that is larger than a minimum value , that rapidly increases as . To estimate , let us consider as an example the Riemann function, i.e. the case , and let us consider the asymptotic behavior of the error

(21) |

as . The asymptotic form of for the Riemann function is given by Theorem 4.11 in Titchmarsh Tich (). If we consider for the sake of clearness the case, one has Tich ()

(22) |

Note that, as , the first term on the right hand side in Eq. (22) gets singular, and to keep the error small one should take large enough to counteract the singular term. For a given value of close to (but larger than) one, an estimate of the truncation index that ensures a small error is readily obtained from Eq. (22) and reads , where

(23) |

The behavior of versus shows a steep increase below . For example, one has for , respectively. This shows that an accurate estimate of the Riemann function can be realized in practice for larger than , and that for the estimation is very accurate even for few lattice sites . Another limitation of our proposed scheme in the estimation of the RH correlation function is decoherence. Obviously, our analysis assumes that the quantum evolution of the system is coherent, i.e. we neglected interaction with the environment. This requires that the observation time is smaller than the coherence time . This limitation will depend on the specific physical implementation of the quantum simulator. In trapped ions coherence time varies in the range of milliseconds to minutes, depending on the qubit realization and the ions in use. On the other hand, in the case of the optical-waveguides realization decoherence is irrelevant and can be ignored, though the maximum value of is still limited by the sample length. The two limitations for the parameters and discussed above, arising from decoherence and from series truncation, are illustrated in Fig. 1.

###### Acknowledgements.

This work was supported by the Fondazione Cariplo (Grant No. 2011-0338).## References

- (1) L. Euler, Comm. Acad. Sci. Imp. Petropol. 9, 160 (1737); B. Riemann, Monatsb. der Berliner Akad., 671 (1859).
- (2) H. Maier and W.P. Schleich, Prime Numbers 101: A Primer on Number Theory (Wiley-Interscience, New York, 2012).
- (3) Dániel Schumayer and David A. W. Hutchinson, Rev. Mod. Phys. 83, 307 (2011).
- (4) M.V. Berry, RiemannÕs zeta function: A model fro quantum chaos, in Quantum Chaos and Statistical Nuclear Physics, T.H. Seligman and H. Nishioka eds., Lecture Notes in Phys. 263 (Springer-Verlag, New York, 1986), pp. 1-17.
- (5) M.V. Berry and J.P. Keating, The Riemann zeros and eigenvalue asymptotics, SIAM Review, 41, 236 (1999).
- (6) T. Kriecherbauer, J. Marklof, and A. Soshnikov, PNAS 98, 10531 (2001).
- (7) L. A. Bunimovich, C. P. Dettmann, Phys. Rev. Lett. 94, 100201 (2005).
- (8) P. Bourgade and J. P. Keating, Seminaire Poincare XIV, 115 (2010).
- (9) B. Julia, in Number Theory and Physics, edited by M. W. J. M. Luck, P. Moussa (Springer, Berlin), p. 276.
- (10) B.L. Julia, Physica A 203, 425 (1994).
- (11) I. Bakas and M.J. Bowick, J. Math. Phys. 32, 1881 (1991).
- (12) D. Schumayer, B.P. van Zyl, and D.A.W. Hutchinson, Phys. Rev. E 78, 056215 (2008).
- (13) G. Sierra and J. Rodriguez-Laguna, Phys. Rev. Lett. 106, 200201 (2011).
- (14) M.V. Berry and J.P. Keating, J. Phys. A: Math. Theor. 44, 285203 (2011).
- (15) C. Weiss, S. Page, and M. Holthaus, Physica A 341, 586 (2004).
- (16) R. Mack, J. P. Dahl, H. Moya-Cessa, W. T. Strunz, R. Walser, and W. P. Schleich, Phys. Rev. A 82, 032119 (2010).
- (17) A. Hurwitz, Z. für Math. und Physik 27, 86 (1882).
- (18) J. H. Wilkinson, Numer. Math 4, 354 (1962); R. S. Martin, C. Reinsch, and J. H. Wilkinson, Numer. Math. 11, 181 (1968).
- (19) M. Christandl, N. Datta, A. Ekert, and A.J. Landahl, Phys. Rev. Lett. 92, 187902 (2004); M. Christandl, N. Datta, T.C. Dorlas, A. Ekert, A. Kay, and A.J. Landahl, Phys. Rev. A 71, 032312 (2005).
- (20) G. M. Nikolopoulos, D. Petrosyan, and P. Lambropoulos, Europhys. Lett. 65, 297 (2004); J. Phys. Condens. Matter 16, 4991 (2004).
- (21) S. Bose, Contemp. Phys. 48, 13 (2007); D. Burgarth, Eur. Phys. J. Special Topics 151, 147 (2007).
- (22) M.-H. Yung, Phys. Rev. A 74, 030303 (2006); G.M. Nikolopoulos, Phys. Rev. Lett. 101, 200502 (2008); T. Caneva, M. Murphy, T. Calarco, R. Fazio, S. Montangero, V. Giovannetti, and G. E. Santoro, Phys. Rev. Lett. 103, 240501 (2009).
- (23) J. I. Cirac and P. Zoller, Phys. Rev. Lett. 74, 4091 (1995); D. Kielpinski, C. Monroe and D. J. Wineland, Nature 417, 709 (2002); H. Häffner, C. F. Roos, R. Blatt, Phys. Rep. 469, 155 (2008).
- (24) D. Porras and J. I. Cirac, Phys. Rev. Lett. 92, 207901 (2004); M. Johanning, A. F. Varòn, and C. Wunderlich, J. Phys. B 42, 154009 (2009); K. Kim, M.-S. Chang, S. Korenblit, R. Islam, E. E. Edwards, J. K. Freericks, G.-D. Lin, L.-M. Duan and C. Monroe, Nature 465, 590 (2010); R. Blatt and C. F. Roos, Nature Phys. 8, 277 (2012).
- (25) A. Kay, Int. J. Quantum Inform. 8, 641 (2010).
- (26) S. Korenblit, D. Kafri, W. C. Campbell, R. Islam, E. E. Edwards, Z.-X. Gong, G.-D. Lin, L.-M. Duan, J. Kim, K. Kim and C. Monroe, New J. Phys. 14, 095024 (2012).
- (27) D.N. Christodoulides, F. Lederer, and Y. Silberberg, Nature (London) 424, 817 (2003); I. Garanovich, S. Longhi, A.A. Sukhorukov, and Y.S. Kivshar, Phys. Rep. 518, 1 (2012).
- (28) S. Longhi, Laser Photonics Rev. 3, 243 (2009).
- (29) R.R. Gattass and E. Mazur, Nature Photon.2, 219 (2008); G. Della Valle, R. Osellame, and P. Laporta, J. Opt. A 11, 013001 (2009); A. Szameit and S. Nolte, J. Phys. B: At. Mol. Opt. Phys. 43, 163001 (2010).
- (30) R. Keil, A. Perez-Leija, F. Dreisow, M. Heinrich, H. Moya-Cessa, S. Nolte, D.N. Christodoulides, and A. Szameit, Phys. Rev. Lett. 107, 103601 (2011); A. Crespi, S. Longhi, and R. Osellame, Phys. Rev. Lett. 108, 163601 (2012); M. Bellec, G.M. Nikolopoulos, and S. Tzortzakis, Opt. Lett. 37, 4504 (2012).
- (31) E. C. Titchmarsh, The theory of the Riemann zeta-function, 2nd edition, edited by D.R. Heath-Brown (The Clarendon Press, Oxford University Press, New York, 1986).