Riemann zeros from Floquet engineering a trapped-ion qubit

The non-trivial zeros of the Riemann zeta function are central objects in number theory. In particular, they enable one to reproduce the prime numbers. They have also attracted the attention of physicists working in random matrix theory and quantum chaos for decades. Here we present an experimental observation of the lowest non-trivial Riemann zeros by using a trapped-ion qubit in a Paul trap, periodically driven with microwave fields. The waveform of the driving is engineered such that the dynamics of the ion is frozen when the driving parameters coincide with a zero of the real component of the zeta function. Scanning over the driving amplitude thus enables the locations of the Riemann zeros to be measured experimentally to a high degree of accuracy, providing a physical embodiment of these fascinating mathematical objects in the quantum realm.


INTRODUCTION
The Riemann zeta function ζ(s) is the Rosseta stone for number theory. The stone, found by Napoleon's troops in Egypt, contains the same text written in three different languages, which enabled the Egyptian hieroglyphics to be deciphered. The ζ-function is also expressed in three different "languages": as the series ∑ n n −s over the positive integers n, as the product ∏ p 1/(1 − p −s ) over the prime numbers p, and as the product / Q n ð1 À s=ρ n Þe s=ρ n over the Riemann zeros ρ n 1 . Note that the first two definitions of ζ(s) are only convergent in the half plane Re½s > 1, while the third definition is valid over the entire complex s-plane. Riemann conjectured in 1859 that these zeros would have a real part equal to a half, ρ n ¼ 1 2 þ iE n , where E n is a real number 2 . This is the famous Riemann Hypothesis (RH), one of the six unsolved Millennium problems (https://www.claymath.org/millenniumproblems), whose solution would amplify our knowledge of the distribution of prime numbers with resulting consequences for number theory and factorization schemes 3 . More poetically, in the words of M. Berry, the proof of the RH would mean that "there is music in the prime numbers" 4 .
One of the most interesting ideas to attack the RH is to show that the E n are the eigenvalues of the Hamiltonian of a quantum system. This idea, suggested by Pólya and Hilbert around 1912 5 , began to be taken seriously in the 70s with Montgomery's observation 6 that the Riemann zeros closely satisfy the statistics of the Gaussian unitary ensemble (GUE). In the 80s Odlyzko 7 tested this prediction numerically for 10 5 zeros around the 10 20 th zero, finding only minor deviations from the GUE. These were explained later by Berry and collaborators [8][9][10] using the theory of quantum chaos, and led him to propose that the E n are the eigenvalues of a quantum chaotic Hamiltonian whose classical version contains isolated periodic orbits whose periods are the logarithm of the prime numbers. Much work has been done [11][12][13][14][15][16][17] to find such a Hamiltonian, but so far without a definitive answer.
We present here an experimental observation of the lowest Riemann zeros, which is quite different from the spectral realization described above. Our intention is not to prove the RH, but rather to provide a physical embodiment of these mathematical objects by using advanced quantum technology. The physical system that we consider is a trapped-ion qubit. The ion is subjected to a time-periodic driving field, and consequently its behavior is described by Floquet theory, in which the familiar energy eigenvalues of static quantum systems are generalized to "quasienergies". These quasienergies can be regulated by the parameters of the driving, in a technique termed Floquet engineering. In particular, when the quasienergies are degenerate (or cross) the ion's dynamics is frozen, which can be observed experimentally. The Riemann zeta function enters into this construction in the design of the driving field, which is engineered to produce the freezing of the dynamics when the real part of ζ(s)/s, with s ¼ 1 2 þ iE, vanishes. Thus observing the freezing of the qubit's dynamics as the driving parameters are varied gives a high-precision experimental measurement of the location of the Riemann zeros.

Floquet theory
We consider a two-level system subjected to a time-periodic driving, described by the Hamiltonian HðtÞ ¼ Jðσ x þ f ðtÞ 2 σ z Þ, where σ x,z are the standard Pauli matrices and J represents the bare tunneling between the two energy levels. Henceforth we will set _ = 1, and measure all energies (times) in units of J (J −1 ). As H(t) is time-periodic, H(t) = H(t + T), where T is the period of the driving, the system is naturally described within Floquet theory, using a basis of Floquet modes and quasienergies which can be extracted from the unitary time-evolution operator for one driving-period U ¼ T exp½Ài R T 0 Hðt 0 Þdt 0 (where T denotes the time-ordering operator). The Floquet modes, Φ j ðtÞ , are the eigenstates of U, and the quasienergies, ϵ j , are related to the eigenvalues of U via λ j ¼ exp ÀiTϵ j Â Ã . The Floquet modes provide a complete basis to describe the time-evolution of the system, and the quasienergies play an analogous role to the energy eigenvalues of a time-independent system. The state of the qubit can thus be expressed as ψðtÞ j i¼ P j α j exp Àiϵ j t Â Ã Φ j ðtÞ , where the expansion coefficients α j are time-independent, and the Floquet modes are T-periodic functions of time. From this expression, it is clear that if two quasienergies approach degeneracy, the timescale for tunneling between them will diverge as 1/Δϵ. Although in general it is difficult to obtain explicit forms for the quasienergies, even for the case of a two-level system, excellent approximations can be obtained in the high-frequency limit, when Ω = 2π/T is the largest energy scale of the problem, that is, Ω ≫ J. In that case one can derive 18 an effective static Hamiltonian, H eff = J eff σ x , where the effective tunneling is given by Here F(t) is the primitive of the driving function, FðtÞ ¼ R t 0 dt 0 f ðt 0 Þ, and the quasienergies are given by ϵ ± ¼ ± J eff j j. The eigenvalues thus become degenerate when they are zero, corresponding to the vanishing of J eff and the freezing of the dynamics. This expression is accurate to first order in 1/Ω, and although in principle higher-order terms could be calculated using the Magnus expansion 19 , we will work at sufficiently high frequencies for this expression to give results of excellent accuracy.
Equation (1) is the key to our approach. By altering the form of the driving, f(t), we are able to manipulate the effective tunneling and the quasienergies of the driven system. Our aim is to obtain a driving function such that J eff (E) is proportional to the real part of yielding an effective Hamiltonian whose dynamics is intimately related to the properties of the ζ-function. In particular, the effective tunneling will vanish, an effect termed "coherent destruction of tunneling" (CDT) 20 , when E coincides with one of the Riemann zeros. In the Methods section we give the details of the mathematical derivation of the driving function, which enables us to obtain Fourier series for f(t) and F(t) (see Fig. 4) which can be straightforwardly programmed into a waveform generator to provide the experimental driving. We choose to focus on the function −ζ(s)/s for two fundamental reasons. The first is that it has a remarkably simple Fourier transform. This also motivated van der Pol 21 and Berry 22 to use this function as the basis for physical implementations of the Riemann zeros in diffraction experiments (in Fourier optics and in antenna radiation patterns respectively). The second reason is that this function decays slowly as E increases (see Fig. 4a). In previous work 23, 24 we proposed to use Floquet engineering to simulate the Riemann Ξ-function 1 . Although successful, the extremely rapid decay of the Ξ-function meant that only the lowest two Riemann zeros were resolvable. In contrast, the slower decay of ζ(s)/s should allow many more quasienergy crossings to be detectable, and thus more zeros to be identified.

Experiment
The experimental results were obtained by periodically driving a single trapped ion with microwave fields. The two-level system is encoded in the hyper-fine clock transition 0 iin a single ytterbium ( 171 Yb + ) ion confined in a Paul ion trap 25 , as shown in Fig. 1a. This clock qubit has the advantages of highfidelity quantum operations and long coherence time 26,27 . The tunneling, J, in this system is of the order of 10 kHz, giving a resonant Rabi time of~100 μs.
After 1 ms of Doppler cooling and 50 μs of optical pumping, the ion is initialized in the ground state 0 j i with a probability ≥99.5%. The qubit is then driven by a microwave field for multiple periods. The driving function was generated from a programmable arbitrary waveform generator (AWG), which generates a phase modulated microwave with a carrier frequency of 200 MHz and a phase modulation function F(t)/2 24 . It is then mixed with a 12.443 GHz fixed-frequency signal and filtered by a high-pass filter. The amplified microwave fields near 12.643 GHz were delivered to the trapped ion from an horn antenna located outside the vacuum chamber. At the end of the multiple periods, the state is projected onto the eigenvector of σ y with eigenvalue +1, y þ ¼ 1 ffiffi 2 p ð 0 j i þ i 1 j iÞ, by applying a π/4 rotation and normal fluorescence detection. When more than one photon is detected, the measurement result is noted as 1; otherwise it is noted as 0. The time evolution of the state population is recorded as a function of the number of periods.
In Fig. 1b we show the experimental protocol, plotted on the Bloch sphere. As noted previously, the Floquet modes Φ 1 ðE; Ω; tÞ j i¼ aðtÞ 0 j i þ bðtÞ 1 j i and Φ 2 ðE; Ω; tÞ j i¼ b Ã ðtÞ 0 j i À a Ã ðtÞ 1 j i are the eigenstates of the one-period time-evolution operator U(E, Ω), where Ω is the driving frequency, and E is a Fig. 1 Experimental procedure to measure Floquet dynamics in a trapped-ion system. a The qubit is encoded in the clock transition of a single trapped 171 Yb + ion. The qubit is periodically driven by a microwave field generated by an AWG and a fixed-frequency signal source, which are both referenced to a 10 MHz rubidium frequency standard. b On the right Bloch sphere, the qubit is initialized to 0 j i and evolves to UðTÞ 0 j i after one period of driving with function f(E = 16, Ω = 5). The fast dynamics ("micromotion") of the qubit during this elementary period is shown on the Bloch sphere, and arises from the intrinsic time-dependence of the Floquet states Φ j ðtÞ . On the left Bloch sphere, the qubit starts from 0 j i (the north pole) and evolves to UðnTÞ 0 j i for multiple periods n = 1, 2, 3, . . . , 30. The red dots are the points measured at n = 5, 10, 15, 20, 25, 30. The insert shows the probability P(n) of projecting the qubit onto y þ ¼ 1 ffiffi 2 p ð 0 j i þ i 1 j iÞ after n periods. P(n) is a sinusoidal function of n, with frequency proportional to the quasienergy ϵ j (E). For values of E satisfying ϵ j (E) = 0, the state is frozen at 0 j i and P(n) = 1/2 for all n. Since ϵ j (E) is proportional to the real part of −ζ(1/2 + iE)/(1/2 + iE), we can therefore identify Riemann zeros by observing the freezing of the dynamics, i.e., P = 1/2.
driving parameter related to the argument of the zeta function, ζ(1/2 + iE). Starting from the initial state 0 j i, the population of the state projected onto y þ after n periods of driving is Pðn; E; ΩÞ ¼ 1=2 À A sinð2nTϵ j ðE; ΩÞÞ, where A ¼ Refað0Þb Ã ð0Þg, and ϵ j (E, Ω) is the quasienergy. It is clear that if E is equal to the zeros of ϵ j (E, Ω), ϵ j (E, Ω) vanishes and P(n, E, Ω) = 1/2 for all n. While if ϵ j (E, Ω) ≠ 0, P(n, Ω, E) evolves sinusoidally with a frequency proportional to the quasienergy ϵ j (E, Ω). The Riemann zeros can thus be identified by observing the freezing of the evolution of P(n, E, Ω) produced by CDT.
To give a quantitative characterization of the state evolution, we define the S parameter as S(Ω, E, y + ) = ∑ n [P(n, Ω, E) − 1/2], where n = {5, 10, 15, 20, 25, 30} are the number of driving periods in the experiment. It is straightforward to see that the zeros of the quasienergies are the zeros of S(Ω, E, y + ) as well. Therefore a scan of the parameter E allows us to identify the Riemann zeros by observing S(Ω, E, y + ) = 0. To give a direct comparison, we show in Fig. 2 the experimental results of S(Ω, E, y + ) and the real part of g(E) as a function of E for Ω = 2, 5, 8, and 16, respectively. For Ω = 2, we can see that S(Ω, E, y + ) shows a significant distortion from the theoretical zeta function and does not allow the position of the Riemann zeros to be determined. Increasing Ω, however, substantially improves the results, and gives excellent agreement between data and theory over the range E ≤ 200, allowing the first eighty Riemann zeros to be resolved. The measurement can be extended to higher E without loss of efficiency and accuracy. This improvement occurs because larger values of Ω satisfy the highfrequency approximation better, and thus Eq. (1) becomes more precise. The difference between the measured zeros and the exact Riemann zeros is shown in Fig. 2g. Most of the zeros can be identified with an accuracy of 1% or better. In Table 1 we present the agreement quantitatively for four different driving frequencies over a wide range of E. We give further, and more detailed, comparisons of the agreement in the Supplementary Information.
Since increasing Ω further would satisfy the high-frequency limit better, it might be thought that the accuracy of the results can be improved by increasing the driving frequency to arbitrarily high values. This is not the case however. As we show in the Methods section, the best results will be obtained when the driving period T is small enough to satisfy the high-frequency limit, while at the same time T is sufficiently large for it to replace the upper limit of integration in Eq. (5). As a consequence of these opposing requirements, the best results will actually be obtained for mid-range frequencies. In Figs. 2d, f and h we show the results for Ω = 16. Comparing with the Ω = 8 result reveals that increasing the frequency has not improved the accuracy of the results.   (23)

Reconstruction of prime numbers
In 1859 Riemann found a formula that gives the number of primes π (x) below or equal to x in terms of the non-trivial zeros ρ n ¼ 1 2 þ iE n 2 . A consequence of this result is a theorem due to E. Landau in 1911 who gave the asymptotic behavior of the sum 3,28 where Λ(x) is the Mangoldt function, that is equal to log p when x = p n and is zero otherwise. The function h(x) therefore has peaks at the primes p and their powers p n . We plot the function h(x) in Fig. 3, together with the function J(x) JðxÞ ¼ which jumps by 1 at every prime and by 1/n at the powers p n . Notice that the experimental error in the zeros does not affect appreciably the location of the peaks.  Fig. 3 Primes from zeros. The function in Eq. (2) with the sum restricted to |E n | < 100 (blue) and the function 5J(x) given in Eq. (3) (red) are calculated (a) using the exact values of E n and (b) using the values of E n given in Table II (Extended Data) for Ω = 16. In both cases one can identify the first eight primes and their powers.  4)). c The blue curve shows FðtÞ ¼ cos À1g ðtÞ cos Et ð Þ , the primitive of the driving function, for driving parameter E = 1. The discontinuities ingðtÞ give rise to discontinuities in this function as well. The red curve shows the Fourier expansion of F(t), truncated at 500 terms. We can note how the fine detail is progressively blurred out as t increases. d The driving function, f(t), for E = 1, obtained as f(t) = ∂ t F(t), plotted over 0 ≤ t < T. The discontinuities in F(t) produce δ-function spikes in the driving function. By construction, f(t) is an even function of t, and so the full periodicity of this function is 2T.

DISCUSSION
We have presented an experimental method for measuring the location of the zeros of the Riemann ζ-function, by using Floquet engineering to control the quasienergy levels of a periodically driven trapped ion. The experimentally measured values of the zeros are in excellent agreement with their theoretical values, and we have demonstrated how they can be used to reconstruct the prime numbers. The high level of experimental control over this system, and the implementation of a driving function derived from the complex function g(E) = − ζ(s)/s (where s = 1/2 + iE), allows as many as the first 80 zeros to be resolved. Our analysis indicates that there is a "sweet spot" for the driving frequency, in which Ω is sufficiently large for the system to be in the highfrequency regime, while its period is large in comparison to the width of the Fourier transform of g(s). Using the experimentally measured zeros we have also obtained a good approximation of the lowest prime numbers. This reconstruction suggests the possibility of a direct experimental realization of the primes. The successful realization of the Riemann zeros in a quantum mechanical system represents an important step along a route inspired by the Hilbert and Pólya proposal, and may lead to further insights into the Riemann hypothesis. Another interesting line of research will be to employ more general modulations of the zeta function that would enable the exploration of different intervals along the critical line 29 .

Driving function derivation
Our starting point is the function gðEÞ ¼ Àζ 1=2þiE ð Þ 12þiE , where ζ(s) is the standard Riemann zeta function. We plot the behavior of this function in Fig. 4a. As van der Pol showed in 1947, its Fourier transform (Fig. 4b) can be written in the surprisingly simple form 21 gðtÞ ¼ e t=2 À e Àt=2 ½e t : where [x] is the integer part of x. As we can see from Fig. 4b, this function is localized around the origin, with an envelope of the form exp Àjtj=2 ð Þ . By dividing the range of integration for the Fourier transform into two halves, it is straightforward to show that the real component of g(E) is given by ðtÞ cos Et dt: In order to observe the location of the Riemann zeros, our interest is focused on values of E > 10. Accordingly we can simply discard the first term, as over this range its magnitude is smaller than the experimental uncertainty in the measurements. Our aim is to obtain a driving function f(t) such that the effective tunneling is proportional to the real component of g(E), that is, where α is the constant of proportionality. Comparing Eq. (5) with Eq. (1), and assuming that the driving period T is sufficiently large to replace the upper limit of integration in Eq. (5) . This choice of α also imposes the condition that the argument of the inverse cosine function is bounded within ± 1 as required, sincegð0Þ is the global maximum ofgðtÞ. We can note that replacing the upper limit of integration with T represents an important restriction on the value of Ω. This replacement means that T must be large in comparison with the width ofgðtÞ, and thus the driving frequency Ω must correspondingly be low. However, for Eq. (1) to be an accurate description of the system's dynamics requires a high value of Ω, so that the system is in the high-frequency regime. Therefore, good results will be obtained in an intermediate range of frequency, when both of these conditions can be adequately satisfied.
We show the form of the F(t) and the driving function for a particular value of E in Fig. 4c and d. The finite discontinuities present in g(E) also produce discontinuities in F(t), and thus δ-function spikes in f(t). A convenient way to obtain f(t) numerically is to expand F(t) in a Fourier series, differentiate the series term by term, and then to re-sum it. As in ref. 23 , we want the driving function to be of definite parity, so that the two Floquet states will be of opposite parity, and so can cross as the driving parameter E is varied. If this parity condition were not satisfied, the von Neumann-Wigner theorem would prevent the quasienergies becoming degenerate, and they could only form broader avoided crossings instead. For this reason we choose to expand F(t) as a Fourier sine series, so that its derivative, f(t) is a cosine series, and is thus an even function of time. Sufficient terms must be included in the series to ensure that the fine structure in f(t) is reproduced with sufficient resolution. Typically in the experiment the series was truncated at 500 terms.
In an analogous way, we could construct a driving function related to the imaginary component of g(E). Unfortunately, however, Im½gðEÞ does not cross the x-axis sharply in the way that Re½gðEÞ does, making the zeros more difficult to locate accurately.

Experimental details
A long coherence time of the system is vital in the experiment. The hyperfine splitting of the ( 171 Yb + ) ion, ω hf ¼ 12642812118:5 þ ω B ð Þ Hz, has a second-order Zeeman shift ω B = 310.8B 2 , where B is the magnetic field in Gauss (G). We used Sm 2 Co 17 permanent magnets to generate a static magnetic field of around B = 9.15 G to reduce the 50 Hz ac-line noise. The whole platform is shielded in a 2-mm-thick μ-metal enclosure to reduce the residual fluctuating magnetic fields 30 . During the experiment, we still observed a slow drift of~±30 Hz of the clock transition in 10 h. This corresponds to ΔB∼ ± 0.005 G, which is mainly due to the temperature drift in the laboratory. This drift is not negligible. Therefore the clock transition frequency was frequently measured by Ramsey type measurements and calibrated by updating the AWG wave frequency during the experiment every half hour.

DATA AVAILABILITY
Source data and all other data that support the plots within this paper and other findings of this study are available from the corresponding authors upon reasonable request.