Entanglement and analytical continuation: an intimate relation told by the Riemann zeta function

We propose measurements on a quantum system to realize the Riemann zeta function ζ. A single system, that is classical interference, suffices to create the Dirichlet representation of ζ. In contrast, we need measurements performed on two entangled quantum systems to extend ζ into the critical strip of complex space where the non-trivial zeros of ζ are located. As a consequence, we can view these zeros as a result of a Schrödinger cat which is by its very construction similar to, but in its details very different from, the superposition formed by two coherent states of identical amplitudes but opposite phases. This interpretation suggests that entanglement in quantum mechanics is the analogue of analytic continuation of complex analysis.


Introduction
The analytic continuation of a function [1] is an important tool in complex analysis. According to Schrödinger, entanglement of two quantum systems is the trade mark [2] of quantum theory. In this paper we connect these key ideas of two seemingly distinct branches of science. We use an entangled quantum system consisting of a harmonic oscillator and a two-level atom to construct an analogue computer whose output is the Riemann zeta function [3].

A brief history of the Riemann zeta function
In analytic number theory [3,4] the Riemann zeta function ζ = ζ (s) of complex argument s ≡ σ + iτ defined for 1 < σ by the Dirichlet sum plays a central role. Indeed, its non-trivial zeros contain information about the distribution of prime numbers. Already in 1737 Leonhard Euler had established a connection [5,6] between the primes and ζ . He could cast the Dirichlet sum representing ζ into an infinite product which apart from the argument s of ζ involves only prime numbers. Most importantly, the Euler product formula contains all primes. The Dirichlet sum of ζ given by (1) converges, provided the real part σ of the argument s is larger than unity. Unfortunately, the product expression of ζ suffers from the very same restriction.
It was not till 1859 that Bernhard Riemann presented the first analytical continuation of ζ to the complete complex plane. This paper [7], written on the occasion of Riemann's election as a corresponding member to the Berlin Academy, contains the by now famous Riemann conjecture: All non-trivial zeros of the zeta function have real part one-half.
It is interesting to note that in this paper Riemann did not present a derivation of his claim but only remarked: 'One would, of course, like to have a rigorous proof of this, but I have put aside the search for such a proof after some fleeting vain attempts because it is not necessary for the immediate objective of my investigation'.
For over 70 years mathematicians searched in vain for any hint on how Riemann came up with this far-reaching hypothesis. However, in 1932 Carl Ludwig Siegel [8], who had extensively studied the sketches of calculations left behind by Riemann, recognized that he must have applied the saddle point method to cast the Dirichlet sum into two sums of rather similar form. This asymptotic representation carries the name Riemann-Siegel formula. In lowest order there is a cancellation of terms provided σ = 1/2, which is the conjectured location of the zeros. However, no rigorous proof of the Riemann hypothesis exists even today [5,6] although there have been many attempts to attack this problem by physical methods [9].

Riemann zeta function from time evolution
Our approach toward the Riemann zeta function takes advantage of the time evolution of two interacting quantum systems followed by a joint measurement [10]. Here, the fact that the interaction has entangled the two systems is of utmost importance. Indeed, a single quantum system does not [11] allow us to obtain the behavior of ζ along the critical line σ = 1/2. Conservation of probability restricts us to a domain of the complex plane to the right of σ = 1 where ζ is defined by the Dirichlet sum (1). This fundamental law cannot be violated for a single quantum system. However, in the context of joint measurements on two quantum systems we can circumvent it.
Three steps encapsulate our operational approach to the Riemann zeta function: (i) take two elementary quantum systems and prepare [12,13] them in an appropriate state, (ii) let them interact for a given time according to a well-designed interaction Hamiltonian and (iii) make a joint measurement on the two systems. The probability amplitude to find the reference state in the time-evolved state of the combined system is then given by the Riemann zeta function.
Our quantum system consists of a harmonic oscillator that represents either a single mode of a cavity field [14] or the motion of a trapped ion [15], together with a single two-level atom. Throughout the paper we assume that the oscillator is a field mode but emphasize that it might equally well be the center-of-mass motion of an ion. Oscillator and atom interact with each other in a non-resonant way.
Hamiltonians of this type have been studied theoretically as well as experimentally in great detail [14,15]. In the context of cavity QED, the dependence of the interaction on the photon number n is often linear. However, for the present purpose of constructing an analogue computer, which calculates the Riemann zeta function, we need a logarithmic dependence. This special choice of the interaction Hamiltonian creates due to the time evolution dictated by the Schrödinger equation the logarithmic phases characteristic of ζ .
For this reason it is useful to introduce a quantum state of the field oscillator whose probability amplitude to find n − 1 photons is given by the nth term of the Dirichlet sum. The overlap of the time-evolved state with the initial state provides us with the Dirichlet representation of the zeta function. For this reason we call this state the zeta state. The superposition of a zeta state evolved forward, and the same state evolved backward in time is at the very heart of the Riemann-Siegel formula. For this reason we have dubbed this superposition state the Riemann-Siegel state.
This state is reminiscent of the Schrödinger cat [16] formed by the superposition of two coherent states of a harmonic oscillator of identical amplitudes but opposite phases. Such states have been realized in experiments involving high-Q cavities [17,18] or ions stored in a trap [19].
We recall [17] that it is necessary to employ a second quantum system, such as a twolevel atom, to create this Schrödinger cat. Indeed, the atom is first prepared in a coherent superposition of its internal states. It subsequently interacts non-resonantly with a single mode of the electromagnetic field in a cavity prepared initially in a coherent state. Due to the non-resonant interaction, the atomic states shift the phase of the coherent state in opposite directions [10]. The resulting state of field plus atom is then an entangled state that is a superposition of a coherent state with a given phase and an atomic state, together with a coherent state of opposite phase accompanied by the other atomic state.
Our proposal for engineering a Riemann-Siegel state is motivated by this technique but differs from it in three important points: (i) We do not start from a coherent state but from a zeta state, (ii) whereas in the case of the Schrödinger cat the interaction is linear in the photon number, we employ for the production of the Riemann-Siegel state a Hamiltonian in which n appears logarithmically, and (iii) in our approach the final atomic measurement is slightly more complicated.
In this context it is interesting to note that Riemann in his original paper [7] had found his analytical continuation of ζ by expressing it in terms of the Mellin transform of the Jacobi theta function. A recent paper [20] has used this representation to depict the Riemann zeta function in phase space using parabolic coordinates as appropriate for the quantum mechanical description of the inverted harmonic oscillator. Needless to say, our approach based on a joint measurement of two interacting quantum systems is distinctly different from this work.
We conclude by noting that much progress has been made [21] with regard to the Hilbert-Polya conjecture [5,6,9], which states that the non-trivial zeros of ζ correspond to the eigenvalues of some Hermitian operator. However, as of today no such operator has been found.

Outline of the paper
Our paper is organized as follows. We introduce in section 2 the truncated zeta state and its time evolution needed in section 3 to create the Dirichlet sum and the Riemann-Siegel representation of ζ by joint measurements. Here we especially emphasize the role of entanglement. Indeed, it is the entanglement between the field oscillator and the atom that allows us to prepare the superposition state of two zeta states of opposite phases, that is the Riemann-Siegel state, and go beyond that domain of the complex plane where the Dirichlet representation of ζ reigns. In this sense entanglement of quantum theory is analogous to analytical continuation of complex analysis.
We dedicate the following sections to a detailed discussion of the zeta as well as the Riemann-Siegel state in phase space. Here we focus on the time evolution of these states, since this feature governs the behavior of the Riemann zeta function along a line in the complex plane parallel to the imaginary axis. Our analysis in section 4 employs the Wigner function [10,22], which brings out most clearly the interference nature of quantum mechanics. Although this distribution is real, it cannot be interpreted as a genuine probability distribution due to the fact that in general the Wigner function assumes negative values. Indeed, the Wigner function of the initial zeta state consists of a positive-valued wedge along the positive x-axis. However, as a function of time this wedge breaks up into several negative-valued wave fronts that curl around the origin of phase space.
Since the Wigner formulation only provides us with the absolute value |ζ | of the zeta function, we analyze in section 5 the Moyal functions [23,24] of these states. When integrated over phase space, they produce the complex value of ζ .
We conclude in section 6 with a brief summary of our results and an outlook. In particular, we show that our method of creating the Riemann zeta function can also be applied to the representation of ζ derived in [25] which not only includes the lowest approximation of the Riemann-Siegel formula employed in the present paper but also the next order correction.
In order to keep the paper self-contained, we briefly review key features of Wigner and Moyal functions in appendix A. Moreover, we dedicate appendix B to a brief discussion of a proposal for realizing an approximate interaction Hamiltonian that displays the logarithmic dependence on the photon number needed to obtain the appropriate phases of the zeta states. Although this Hamiltonian is only an approximation, the resulting scalar product determining ζ leads to a satisfactory agreement with the exact values of ζ .

Zeta states and their properties
Before we discuss the description of the zeta function in phase space, it is useful to introduce the truncated zeta state |σ, τ n 0 ≡ N n 0 (σ ) n 0 n=1 1 n σ/2 e −iτ ln n |n − 1 (2) expressed in terms of a finite number n 0 of photon number states |n of a single mode of the electromagnetic field. The normalization constant N n 0 follows from the condition and depends on σ . Without loss of generality, we choose all normalization constants to be real. We recall that the methods of state engineering [12] and state synthesis [13] are especially tailored to create such truncated states of a single mode of the radiation field. In the limit of n 0 = ∞, we define the zeta state |σ, τ ∞ ≡ |σ, τ = N ∞ n=1 1 n σ/2 e −iτ ln n |n − 1 (4) and the normalization condition (3) reduces with the abbreviation N ∞ ≡ N to In the last step we have used the Dirichlet representation (1) of the zeta function. Since this sum is only convergent for 1 < σ , the existence of the zeta state is restricted to the same domain of the complex plane. The so-defined zeta states are generalizations of the thermal phase state [26] Indeed, both are coherent superpositions of photon number states with real-valued expansion coefficients that decay with n. However, for the thermal phase state (6) we find an exponential dependence of the photon distribution on the photon number n with a maximum at the vacuum, that is at n = 0, whereas for the zeta states (4) we obtain only a polynomial dependence Hence we have replaced in (7) the quantum number n by ln(n + 1). We conclude by noting that the Hamiltonian H R ≡h ln(n + 1) gives rise to the time evolution of the truncated zeta state with τ ≡ t.

State engineering with the Riemann zeta function
Our goal is to build an analogue computer based on a quantum system whose output is ζ = ζ (s). Here we take advantage of the time evolution of an appropriate quantum state due to a specifically designed Hamiltonian. The crucial ingredient of our approach is the interaction HamiltonianĤ RS ≡Ĥ R (|e e| − |g g|) =h ln(n + 1) (|e e| − |g g|) (11) describing the interaction of a two-level atom consisting of a ground state |g and an excited state |e with a single mode of the radiation field. Here is the Rabi frequency that establishes the coupling between the atom and the field. This interaction is motivated by the Hamiltonian of the strongly detuned Jaynes-Cummings model frequently used in cavity QED [10,14] and ion trap physics [15]. In our case the photon number operatorn is replaced by ln(n + 1).
The Riemann zeta function arises as an appropriate projection of the time-evolved quantum state of the field and the atom. Indeed, the real part σ of the complex-valued argument s of ζ = ζ (s) determines the initial occupation probability of the field oscillator whereas the imaginary part τ is proportional to time.

Dirichlet representation
With the help of the zeta states (4) we can summarize the engineering of the Dirichlet sum (1) in a compact way. Indeed, when we choose the initial state | (0) ≡ |σ, 0 |e of the combined system consisting of field plus atom, the time-evolved state following from (11) reads | (t) = |σ, t |e . Therefore, the probability amplitude to find the initial state in the time-evolved state takes the form When we recall the normalization constant (5) we arrive at Hence, the probability amplitude C representing the overlap between the initial state |σ, 0 and the time-evolved state |σ, t follows the complex-valued Riemann zeta function ζ = ζ (s) = ζ (σ, τ ) in its dependence on the imaginary part τ ≡ t of the argument s. As a consequence, C(t) as a function of the interaction time allows us to read out ζ = ζ (σ, t) along an axis in complex space parallel to the imaginary axis. The photon statistics (8) depends on the parameter σ and therefore defines the location of this line.
So far, the atomic states have not entered into this analysis. Indeed, the atom has started in the excited state and remains in this state throughout the dynamics. Hence, a single quantum system suffices to obtain the Dirichlet representation of the Riemann zeta function.
In [11] we have followed this line of thought and have represented the Riemann zeta function by wave packet dynamics. In particular, we have analyzed the motion of a particle in a potential whose energy eigenvalues E n satisfy the law where , γ and a are parameters. Moreover, we have constructed such a potential 1 .
A wave packet given initially by a superposition of energy eigenstates of this nonlinear oscillator in the form of the zeta state |σ, 0 develops logarithmic phase factors leading to |σ, t . As a result the auto-correlation function of this wave packet is also determined by (14).
We emphasize that there is really no need for a quantum system. For example, an interferometer with classical light can implement the Dirichlet representation of the Riemann zeta function provided we introduce logarithmic phase shifts in the individual arms numbered by the index n. This realization shows that the convergence condition of the Dirichlet sum corresponds to the requirement of conservation of intensity. Indeed, as σ approaches unity from above, more and more arms are necessary to implement ζ . Still, the intensities in the individual arms have to add up to the intensity entering the interferometer. This classical example translates the mathematical convergence criterion of the Dirichlet sum into the physical conservation law of intensity and clearly demonstrates that a single-system approach cannot lead us beyond the Dirichlet wall at σ = 1.

Riemann-Siegel formula
There exist many possibilities [3] to perform an analytical continuation of the Riemann zeta function and overcome this wall. For the present discussion the most convenient one is the one based on the Riemann-Siegel formula [8]. Indeed, we can approximate ζ by the expression [4] ζ (s) ∼ = where Here denotes the Gamma function and [x] is the nearest integer smaller than or equal to x. Equation (15) represents only the lowest approximation of the complete Riemann-Siegel formula discussed in section 6. When we express the Riemann-Siegel formula (15) in the real and imaginary parts of the argument s and shift the summation index by unity, that is the two truncated Dirichlet sums have opposite phases ±iτ ln(n + 1). This behavior is reminiscent of the generation [10,14,15] of Schrödinger cats [16] with the help of the Hamiltonian (12) of the strongly detuned Jaynes-Cummings model. Indeed, there each part of the entangled state accumulates a phase ±i t n.
Hence, we start with the initial state consisting of the truncated zeta state |σ, 0 n 0 defined by (2) and a superposition of the two atomic states with equal weight. The time evolution of the state (18) with the HamiltonianĤ RS given by (11) yields the Riemann-Siegel state Here we have made use of the relation (10). The reference state with the normalization provides us with the probability amplitude When we recall the definition (2) of the truncated zeta states, we find the explicit expression that is, the overlap is proportional to the Riemann-Siegel formula (17) provided and Since β ≡ β(σ, τ ) the normalization constant N n 0 given by (21) depends on τ ≡ t as well. Therefore, the reference state | defined by (20) seems to be time dependent. But we have to keep in mind that in the definition of the reference state τ ≡ s is merely a parameter chosen to match the time-evolved state | (t) producing in the overlap (23) the Riemann-Siegel formula ζ RS .
Moreover, we find from (19) that | (t) becomes entangled for t > 0 independent of the value of σ . In contrast, the reference state | is always entangled, except for σ = 1/2, that is on the critical line. From the definition (24) of β we find β 1 2 , τ = χ 1 2 + iτ , which together with the expression following from (16) yields the product state of field oscillator plus atom.

Riemann zeta function from Wigner phase space
In the preceding section we have introduced the zeta and the Riemann-Siegel states that allow us to obtain the Riemann zeta function in the Dirichlet and the Riemann-Siegel representations. One way to gain deeper insight into the structure of these states is to represent them in quantum phase space. Here the Wigner function [10,[22][23][24] is particularly useful since it reflects the interference properties of these quantum states. Moreover, the quantum mechanical scalar product between two quantum states, which is crucial in our realization of the Riemann zeta function, is given by the overlap of the Wigner functions corresponding to the two states.
In this section we analyze the behavior of the Wigner function, or the Wigner function matrix of the zeta and the Riemann-Siegel states, respectively. For an overview of these phase space representations, we refer to appendix A. We emphasize that there exist many experimental techniques [29][30][31] to reconstruct such Wigner functions from measured data. Therefore, the results of the present section are not only a visualization of the states underlying our quantum system giving rise to ζ , but these visualizations are even experimentally accessible.

Wigner function of the zeta state
Throughout our paper we use the definition of the Wigner function where the position x and the wave vector k as the conjugate variables represent phase space and Planck's constant does not enter. Substitution of the zeta state (4) into this definition yields the Wigner function Here, W |n m| denote the Moyal functions of the photon number states defined by (A.10) and W |n ≡ W |n n| is the corresponding Wigner function. For a detailed evaluation of Wigner functions of a superposition of photon number states, we refer to section A.3 of the appendix.
In figure 1 we represent the Wigner function W |σ,τ of the zeta state |σ, τ by colored density plots for several values of τ and a fixed value of σ . Here and in the remainder of our paper, warmer colors indicate positive values of the Wigner function whereas colder colors denote negative values. The 'thermometer' on the right-hand side of the pictures provides us with the actual values.
We emphasize that the initial zeta state |σ, 0 is represented in phase space by a wedge aligned along the x-axis. It is reminiscent of a thermal phase state (6). in phase space, which is inversely proportional to n, leads to diffusion. Indeed, the parts of the wedge corresponding to small values of n move faster than those where n is large. Moreover, due to the minus sign in (29) the motion in phase space is in the clockwise direction as familiar from unitary time evolution in quantum mechanics. As a result, the original wedge starts to break up and curls around the origin as depicted in figure 1 for different times τ . When we apply the familiar identity that the scalar product between two quantum states can be expressed [10] by an integral over phase space to the probability amplitude C given by (13) and (14), we find Hence, the absolute value squared of the Riemann zeta function ζ (σ, τ ) is determined by the overlap in phase space between the wedge representing the initial Wigner function W |σ,0 and the Wigner function W |σ,τ of the time-evolved state |σ, τ .
In figure 2 we recall in the upper left picture the initial Wigner function W |2,0 for the zeta state |2, 0 . In the middle row we display the Wigner function W |2,τ for the three times τ = 1, 5 and 10. Finally, in the bottom row we show the products of the Wigner functions of the initial state and the time-evolved states which, when integrated over phase space, determine according to (30) the values |ζ (2, τ )| 2 marked in the top right picture by dots.

Wigner matrix of entangled states
As shown in section A.2 of the appendix the Wigner representation of the entangled state | = |ψ e |e + |ψ g |g Figure 2. The absolute value squared of the Riemann zeta function |ζ (2, τ )| 2 in its dependence along a line in complex space parallel to the imaginary axis with σ = 2 represented as overlap in Wigner phase space. According to (13) the Riemann zeta function given by the Dirichlet sum (1) follows from the scalar product between the initial zeta state |σ, 0 and the time-evolved state |σ, τ . In the Wigner formulation of quantum mechanics, this scalar product is given by the overlap in phase space between the corresponding Wigner functions W |2,0 (top row, left picture) and W |2,τ presented in the middle row for three characteristic values of τ . The bottom row shows the product of the two Wigner functions W |2,0 and W |2,τ . The integral of this product over phase space determines according to (30) the values of |ζ (2, τ )| 2 marked in the top right picture by dots. is given by the Hermitian operator W | ≡ W |ψ e |e e| + W |ψ g |g g| + W |ψ e ψ g | |e g| + W |ψ g ψ e | |g e|.
Here W |ψ i and W |ψ i ψ j | with i, j ∈ {e, g} denote the Wigner and Moyal functions of the unnormalized field states |ψ i , respectively. We first describe in this section the Wigner matrix elements of the Riemann-Siegel state (19) and its reference state (20), and then show how their overlap produces the Riemann-Siegel formula (17).  (19) is given by the expression It is crucial to keep in mind that here n 0 is a constant whereas τ = t depends explicitly on time.
In order to gain more insight into the structure of the Wigner matrix (32), we first show in figure 3 for the example of σ = 1/2 the Wigner function of the initial truncated zeta state |σ, 0 n 0 for different values n 0 . Since the Moyal functions W |n m| are symmetric in k, this symmetry is also contained in W |σ,0 n 0 . Moreover, the maximum of the photon statistics of the truncated zeta state at n = 0 makes W |0 the dominant contribution of the first sum in (33). It has a maximum at the origin, as confirmed by the left picture of figure 3, which gets shifted on the x-axis to the right due to the influence of the Moyal functions in the double sum. Here the main contribution originates from W |1 0| . Furthermore, the double sum in (33) creates n 0 − 1 minima around the origin which get blurred due to the interference of both sums. As expected, the Wigner function develops a wedge along the x-axis for n 0 → ∞ since in this limit the truncated state |σ, 0 n 0 approaches the zeta state |σ, 0 , but equality can only be reached for σ > 1.
For τ > 0 the Wigner function W |σ,τ n 0 corresponding to the truncated zeta state |σ, τ n 0 is given by The time dependence in the double sum rotates the wedge clockwise around the origin, again with different velocities for different n, which leads to a break-up of the wedge. In the upper row of figure 4, we show this behavior for 1 2 W | 1 2 ,τ 10 . The Wigner function 1 2 W | 1 2 ,−τ 10 displays the same patterns mirrored at the x-axis since |σ, −τ n 0 rotates counter-clockwise, that is W |σ,τ n 0 (x, k) = W |σ,−τ n 0 (x, −k), as depicted in figure 4 in the row below.
Next we turn to the off-diagonal matrix element Since the Wigner matrix is Hermitian, we find W |σ,−τ n 0 σ,τ | = W * |σ,τ n 0 σ,−τ | . In contrast to the diagonal elements W |σ,τ n 0 given by (35), both sums in (36) have timedependent phase factors. However, they simply act as complex scaling factors. Indeed, we add in the double sum the real parts of the unrotated Moyal functions {W |n−1 m−1| } whereas in W |σ,τ n 0 defined by (35) we sum over {e −iτ ln n m W |n−1 m−1| }, which corresponds to adding rotated Moyal functions. Therefore, the off-diagonal elements do not rotate in phase space and the symmetry in k is preserved for all times, that is Furthermore, the initial off-diagonal elements are equal to the diagonal ones (35), W |σ,0 n 0 σ,0| = W |σ,0 n 0 , depicted in figure 3.
For τ = 0 and n 0 > 1, however, the contributions in the real or imaginary parts of W |σ,τ n 0 σ,−τ | are scaled by the cosines or sines of the complex scaling factors that suppress positive or even produce negative contributions. This fact leads to a different distribution of positive and negative regions in the patterns of the whole Wigner function as shown in the bottom rows of figure 4.

Reference state. The definition (20) of the reference state | yields with (31)
The shapes of the diagonal matrix elements shown in the two upper rows in figure 5 are identical to those of the Wigner function W |σ,0 n 0 depicted in figure 3. However, the scaling factor |β| 2 , defined by (24), decreases in the critical strip for increasing values of σ . Hence, the matrix element with W |2−3σ,0 n 0 gets more pronounced than the corresponding matrix element with W |σ,0 n 0 for σ < 1/2, as shown in the upper pictures of figure 5, although we would assume the contrary by the photon statistics (34) of the truncated zeta state.
In the bottom rows we have depicted the real and imaginary parts of N 2 n 0 β W |σ,0 n 0 2−3σ,0| for τ = 700 (n 0 = 10) and different values of σ . The pattern of the real part is for σ = 1/3 almost symmetric whereas that of the imaginary part is nearly antisymmetric since arg χ( 1 2 + 700 i) ≈ 0.89 π .  (19) for σ = 1/2 and n 0 = 10 giving rise via (32) to the Wigner matrix W | . As expected, the diagonal matrix element W | 1 2 ,τ 10 associated with the excited atomic state, given by (35) and depicted in the upper row, rotates clockwise around the origin, whereas W | 1 2 ,−τ n 0 , connected with the ground state and shown below, rotates counter-clockwise. As in the case of W |σ,τ shown in figures 1 and 2, the Wigner functions W |σ,τ n 0 and W |σ,−τ n 0 get blurred due to the different velocities of the individual photon number states. However, the offdiagonal matrix elements W |σ,τ n 0 σ,−τ | given by (36) do not rotate. The timedependent factors in front of each Moyal function merely act as complex scaling factors. Therefore, the shapes remain symmetric in k as indicated by the two lower rows.  (20) at τ = 700 (n 0 = 10) for σ = 1/3 and 1/2. The diagonal matrix element including W |2−3σ,0 10 shown in the second row is due to its scaling factor |β| 2 for σ < 1/2 more pronounced than the corresponding matrix element with W |σ,0 10 associated with the excited state. On the critical line, the two matrix elements are the same as shown in (39). Moreover, the diagonal matrix elements are, for all values of τ symmetric in k, independent of the value of σ . However, the real and imaginary parts of the off-diagonal matrix element, connected to the Moyal functions W |σ,0 10 2−3σ,0| and depicted in the bottom rows, are only on the critical line symmetric in k. The patterns for σ = 1/3 are nearly symmetric or anti-symmetric, respectively, since arg χ ≈ 0.89 π .
We conclude the analysis of the Wigner matrix of the reference state by a closer look at the special case σ = 1/2. Since the reference state is on the critical line given by the product state (26) the Wigner matrix (38) reduces to where we have used (25). All matrix elements have the shape of W | 1 2 ,0 n 0 . They are merely scaled by the factor 1/2 and in the case of the off-diagonal matrix elements additionally by the phase factor e ±iϕ . Therefore they are again symmetric in k.

Overlap.
We have now all ingredients to represent the overlap (23) determining the Riemann-Siegel formula in terms of combinations of the corresponding Wigner functions. In this way, we show how the values of ζ emerge from the patterns created in phase space by the overlap of the reference and the Riemann-Siegel state.
According to (A.8) and (A.9) we find where the kernel K defined by (A.9) for the Riemann-Siegel representation reads Here we have used the Hermiticity of the Wigner matrix to combine the off-diagonal matrix elements. When we recall the properties of the matrix elements, (41) reveals that K is only on the critical line symmetric in k. In figure 6 we depict K for different values of σ and τ . The right pictures show that for σ = 1/2 and τ 1 (n 0 = 10), chosen near a maximum of |ζ | 2 , areas in phase space with positive weight dominate whereas at τ 2 , being a zero of ζ , negative contributions cancel positive ones quite well. The quantitative analysis of the data in table 1 confirms this observation.
For σ = 1/3 we see that the structures of the left pictures of figure 6 are less pronounced than those for σ = 1/2 shown on the right. Nevertheless, the value of K at τ 1 evaluated in table 1 is for σ = 1/3 nearly the same as on the critical line, that is for σ = 1/2. However, for τ 2 it is by a factor of 30 higher for σ = 1/3. This behavior strongly indicates that ζ (1/3 + iτ 2 ) cannot be zero, as expected, since due to the structure of the main part (17) of the Riemann-Siegel formula the zeros of ζ have to be on the critical line.

Riemann zeta function described by Moyal functions
In the preceding section we have presented a description of the Riemann zeta function based on the Wigner function. Unfortunately, this approach provides us only with the absolute value of ζ . However, according to (A.4) the scalar product of the operator |ψ φ|. In this section we discuss the Moyal functions leading by integration to the Dirichlet sum and the Riemann-Siegel representation of the zeta function.

Zeta state
With the help of (42) we find from (13) and (14) the expression  (44) In contrast to the Wigner function W |σ,τ given by (28) of the zeta state |σ, τ , both sums contain time-dependent factors leading to complex-valued Moyal functions W |σ,τ σ,0| for τ > 0. Figure 7 shows the Moyal function W |2,τ 2,0| represented by real and imaginary parts (two top rows), and by amplitude and phase (two bottom rows) for four different values τ = t. The initial Moyal function is, of course, given by the Wigner function W |2,0 , which is real. At times τ = 0 the wedge along the positive x-axis curls again around the origin but the negative areas that appear in the real part as well as the imaginary part of W |2,τ 2,0| are more pronounced than in the Wigner function W |2,τ shown in figure 1. The absolute value and the argument of W |2,τ 2,0| also reflect the rotation of the wedge.

Riemann-Siegel state
The Riemann-Siegel formula (17) is expressed by the weighted sum (22) of the two scalar products n 0 σ, 0|σ, τ n 0 and n 0 2 − 3σ, 0|σ, −τ n 0 involving the truncated zeta states |σ, 0 n 0 , |σ, ±τ n 0 and |2 − 3σ, 0 n 0 . With the help of (42) we can now represent (22) as the phase space integral In contrast to the Wigner function matrix W | given by (32), the Moyal function M | | ≡ W |σ,τ n 0 σ,0| + β(σ, τ ) W |σ,−τ n 0 2−3σ,0| (45) Figure 7. Complex-valued Riemann zeta function ζ (2, τ ) represented in phase space by the Moyal function W |2,τ 2,0| defined by (44), and depicted in its real and imaginary parts (top rows) and absolute value and argument (bottom rows) for four different values of τ = t. The initial real part is, of course, the Wigner function W |2,0 , which is real. For τ = 1 we find a shorter positive wedge with a negative ridge, whereas the imaginary part displays a wedge with mainly negative parts. For larger times the wedges curl more and more around the origin while the negative parts are more dominant than in the Wigner function W |2,τ . The lower rows, which depict the absolute value and the phase of W |2,τ 2,0| , confirm the curling of the wedge.
of the entangled states | and | is not an operator but the sum of two Moyal functions of the form W |σ 1 ,τ n 0 σ 2 ,0| . The representation W |σ 1 ,τ n 0 σ 2 ,0| = N n 0 (σ 1 ) N n 0 (σ 2 ) n 0 n,m=1 e −iτ ln n n σ 1 /2 m σ 2 /2 W |n−1 m−1| , In the upper pictures of figure 8 we show the real and imaginary parts of M | | and in the lower ones the corresponding absolute value and phase for the same values of σ and τ as in figure 6. Apart from the absolute value |M | | | of the Moyal function for σ = 1/2, which is symmetric in k, the distributions have in general no overall symmetries.

Conclusion and outlook
The quantum dynamics of a two-level atom coupled in a logarithmic way to a harmonic oscillator such as a single mode of the radiation field allows us to construct an analogue computer for the Riemann zeta function. Here joint measurements on both systems, that is on atom plus field, play a central role. For the realization of the Dirichlet representation of ζ , the internal states of the atom are of no importance. Indeed, classical interference suffices to obtain ζ = ζ (s) for arguments s in the domain of the complex plane with 1 < s. However, in the region only accessible by an analytical continuation of ζ the entanglement between the atom and the field is crucial. In this sense the tool of analytical continuation of complex analysis is intimately related to the concept of entanglement of quantum mechanics.
We conclude by briefly discussing three questions that emerge in this context. Obviously, a complete analysis goes beyond the scope of the present paper. For this reason, we only intend to provide some ideas of how to approach these problems.
Throughout our paper we have relied on the Riemann-Siegel representation of ζ . However, the expression (17) with n 2 × 10 −8 τ and δ ≡ √ τ − √ 2π n 0 + 1 2 . Here a k is given by the recurrence formula (n + 1) √ τ a n+1 = −(n + 1 − σ ) a n + i a n−2 , where a −2 ≡ a −1 ≡ 0 and a 0 ≡ 1 and . Hence, we need to generalize our approach and provide a procedure to implement (47). Many ideas offer themselves. Here we only outline one of them. According to [25] we can also represent ζ by a series of convergent sums where the main term with n B n 0 and already contains ζ RS and its first correction term. Here, Erfc denotes the complementary error function [32]. Moreover, this representation is free of discontinuities and, in contrast to the Riemann-Siegel formula, does not jump at τ = 2π n 2 0 . For the sake of simplicity, we use in what follows n B = n 0 and define the reference state Its overlap with the Riemann-Siegel state (19) produces with the Berry-Keating representation ζ BK given by (48). Thus, we can include the first correction term of the Riemann-Siegel formula by using a different reference state. The second question is directed toward the experimental realization of the interaction Hamiltonian (11) that is logarithmic in the photon number. In appendix B we show that, in general, we can simulate the behavior of the interaction by an effective Hamiltonian reminiscent of nonlinear Jaynes-Cummings models [33].
Finally, we wonder whether it is really necessary to employ a second quantum system, that is the atom, in order to prepare a superposition of |σ, τ n 0 and its phase-conjugated image, that is of |σ, −τ n 0 . Would a single system be sufficient after all? The familiar problem of a particle in a box suggests an affirmative answer. Indeed, there unitary time evolution can create at a well-defined moment of time a superposition of a wave packet and its phase-conjugated image.
However, due to the vanishing boundary conditions at the infinitely steep and infinitely tall walls, the particle in the box [34,35] is very different from any other realistic model in quantum mechanics. Indeed, the time evolution of a wave packet in a box is equivalent to the free evolution of a periodic array of a superposition of this wave packet and its phase-conjugated wave. Hence, in the example of the box the appearance of the desired superposition state is not a consequence of the time evolution but of the initial condition. This fact confirms our claim that two quantum systems are needed to obtain the wave and its phase-conjugated counterpart from time evolution.
Obviously, many open questions will remain. Nevertheless, we are convinced that analogue computers for mathematical functions of the type reported in this paper may shine some new light on old and deep questions of number theory. W T Strunz, R Walser and A Wünsche on this topic. This work has been supported by the Landesstiftung Baden-Württemberg in the framework of the Quantum Highway A8 and the Ministerium für Wissenschaft und Kunst Baden-Württemberg.

Appendix A. Wigner representation
In the present appendix we extend the definition of the Wigner function to operatorsÔ ≡ |ψ φ| consisting of two different states |φ and |ψ and show that these Moyal functions [23,24] provide us with a phase space representation of the scalar product ψ|φ . Moreover, we recall [36] a generalization of the Wigner formalism for states of a single quantum system to entangled states | leading to a Wigner matrix W | .

A.1. Moyal functions
In analogy to the Wigner function (A.1), we define the Moyal function [23,24] by where |ψ and |φ are two arbitrary states of a single system that are, in general, not orthogonal to each other. Thus, the Moyal functions are complex, that is Since the entangled state | has to be normalized, the states |ψ i cannot be normalized and therefore their Wigner functions W |ψ i fulfil the normalization condition which is an immediate consequence of (A.7).

A.2. Wigner matrix
Following the example of [36] we define the so-called Wigner function matrix The Wigner function matrix W | contains the complete information about the two entangled systems. However, when we trace over one subsystem, e.g. | j , the Wigner matrix reduces to the Wigner function and information about the other system is erased. Moreover, since this representation only contains the diagonal matrix elements W |ψ i the information about the entanglement between the two systems is contained in the off-diagonal elements of the Wigner matrix. They emerge in the scalar product of the two states | and | where the kernel K of the integral reads Here we have used (A.5) to rewrite the product of two scalar products and the property (A.3) of the Moyal functions. Equation (A.8) clearly shows that not only the diagonal terms of the Wigner matrices contribute to the scalar product but also the off-diagonal terms.

A.3. Wigner representation of superposition of photon number states
Next we briefly recall the form of the Wigner function for an arbitrary superposition |ψ ≡ ∞ n=1 ψ n |n − 1 of photon number states which yields the Wigner function where the complex-valued Moyal function Here L m n denotes the generalized Laguerre polynomial and M is a quantity with the units of mass.
Separation into diagonal and off-diagonal terms yields Finally, we find that the entangled state produces the Wigner matrix The diagonal matrix elements with i = j reduce to However, for i = j only the case ψ ( j) * n = ψ (i) n leads to the simplified expression

Appendix B. Riemann Hamiltonian from a Jaynes-Cummings model
In this appendix we derive a Jaynes-Cummings-like Hamiltonian as an approximation for the Riemann-Siegel HamiltonianĤ RS defined by (11). Although this approximate Hamiltonian does not quite display the logarithmic dependence ofĤ RS , it nevertheless provides us with a surprisingly good approximation of the Riemann zeta function. When we start in the Schrödinger picture with the Hamiltonian where the free evolution is given byĤ 0 ≡hωn +¯h 2σ z and the Hermitian operatorf only acts in the Hilbert space of the field, the transformation into the interaction picture yieldsĤ m|f |n e −iωt (n−m) |m n| (σ e −i t +σ † e i t ).
The Pauli operators are given byσ ≡ |e g| andσ z ≡ |e e| − |g g|, respectively. After a lengthy but straightforward calculation we find that in second-order perturbation theory the time evolution resulting from this Hamiltonian is given by the effective Hamiltonian Here, we have introduced the inverse energies ε (±) j ≡ 1 2h 1 ν j ± 1 ν − j and ε 0 ≡ 1 hν 0 with ν j ≡ + jω. We note that only for the choicef ≡h ln(n + 1) is the effective HamiltonianĤ eff given by (B.2) in complete agreement with the HamiltonianĤ RS defined in (11). Since our realization of the Riemann zeta function is determined by the overlap of two quantum states, an approximation ofĤ RS suffices as shown by the following example of the zeta state.
Indeed, when we compare the effective overlap between the state created by the time evolution of the zeta state (4) due to the effective Hamiltonian (B.2) and the initial state to the overlap (13), we find that (B.3) benefits from three properties: (i) the argument of the exponential function in (B.3) must be equal to ± t ln(n + 1) only modulo 2π , (ii) the probability amplitudes are proportional to n −σ , which reduces for large n-values the importance of deviations ofĤ eff fromĤ RS , and (iii) deviations for one special n can be compensated by the others since only the sum over all n matters. Hence, even if the effective Hamiltonian does not provide a good approximation of the Riemann-Siegel Hamiltonian the overlap can nevertheless approximate the zeta function very accurately.
To bring out this fact most clearly we now, inspired by the familiar Jaynes-Cummings model [10] wheref ∼â +â † , choosê Moreover, we can reduce the number of coefficients f µ by choosing small Rabi frequencies, that is ω. Hence, the largest contribution in the phase of (B.3) is given by shown beneath in blue, which differs slightly from the value 1 depicted by the red line. Including more coefficients f 2µ and optimizing them will lead to much better agreement.
In the bottom frame of figure B.1 we display the ratio of the exponents, which shows that (+) n is larger than . Hence, ζ eff provides an approximation of the zeta function although the exponent F n−1 + C n−1 is not even close to the logarithm. The resulting curve as a function of the dimensionless time, that is, of the imaginary part τ of the argument of ζ , oscillates around the exact curve |ζ (2, τ )| depicted by the red curve. The middle frame depicts in blue the ratio , defined by (B.6), between |ζ eff | and the exact values |ζ |. However, as expected, this ratio is much smaller than the ratio of the exponents (+) n (B.7) shown at the bottom. The red lines indicate perfect agreement in the middle and the bottom pictures.
In summary, this example shows that, in general, an approximation of the Hamiltonian H RS is possible. Since the quality of the approximation depends on many parameters, a detailed analysis of the optimization problem must be the subject of a further investigation which would exceed the scope of this paper.