The Prime state and its quantum relatives

The Prime state of $n$ qubits, $|\mathbb{P}_n\rangle$, is defined as the uniform superposition of all the computational basis states corresponding to prime numbers smaller than $2^n$. This state encodes, quantum mechanically, arithmetic properties of the primes. We first show that the Quantum Fourier Transform of the Prime state provides a direct access to Chebyshev-like biases in the distribution of prime numbers. We next study the entanglement entropy of $|\mathbb{P}_n\rangle$ up to $n=30$ qubits, and find a relation between its scaling and the Shannon entropy of the density of square-free integers. This relation also holds when the Prime state is constructed using a qudit basis, showing that this property is intrinsic to the distribution of primes. The same feature is found when considering states built from the superposition of primes in arithmetic progressions. Finally, we explore the properties of other number-theoretical quantum states, such as those defined from odd composite numbers, square-free integers and starry primes. For this study, we have developed an open-source library that diagonalizes matrices using floats of arbitrary precision.


INTRODUCTION
Mathematical sequences of integers are key to Number Theory and other areas of Mathematics. Properties of such sequences can be studied on a quantum computer by creating pure states consisting of superpositions of computational-basis vectors that encode the numbers appearing in the sequence, e.g. in binary format using qubits. We call these states number-theoretical quantum states. This novel approach to sequences using the tools of Quantum Information Theory was introduced in Ref. [1], where it was shown that a uniform superposition of prime numbers, i.e. the Prime state, can be created efficiently on a quantum computer by using Grover's search algorithm [2] with a quantum oracle that codifies a primality test.
An example of the potential borne by the use of quantum methods to address problems related to prime numbers is provided by the possible numerical verification of the Riemann Hypothesis [3] on a quantum computer beyond what can be achieved using classical methods [4]. A Prime state with a minimum of 90 logical qubits would be needed for that. This is achieved by estimating, with a quantum counting algorithm [5], the value of the prime counting function π(x), whose fluctuation around its av-erage value x/ ln(x) is bounded by the Riemann Hypothesis.
Some relevant properties of the prime sequence (and presumably other sequences of numbers and their corresponding quantum states) are encoded into the density matrix and the quantum correlations among partitions of the Prime state, and are therefore amenable to investigation using a quantum computer. This approach does not require the full access to the density matrix of the state in order to extract information about the distribution of primes, for that would cost a number of operations that scales exponentially with the number of qubits. Instead, once the Prime state is constructed, several number-theoretical functions become direct observables that can be simply measured in the computational basis. An example is provided by the Chebyshev bias ∆(x) [6], which is the difference in the number of primes below some x appearing in the arithmetic progressions 4k + 3 and 4k + 1, with k integer. This bias can be obtained by measuring just a single qubit of a Prime state with log 2 (x) qubits [1].
It is specially interesting to note that the Prime state bears an amount of entanglement that is almost maximal across many, if not all, bi-partitions (as characterized by e.g. the von Neumann entropy [7]). As a consequence, the Prime state is genuinely quantum, that is, it encodes correlations that cannot be described in classical terms. This fact is to be expected. Very large entanglement is tantamount to large quantum correlations related to nonlocality, and it brings a quantification of surprise among partitions in the system. Digits in a part of the register are almost maximally surprised to relate to those in the rest of the register. If this were not the case, it would be easier to predict the appearance of primes.
Therefore, the main interest in studying the entanglement traits of the Prime state, and other numbertheoretical quantum states, is that it is not at all unreasonable to think that its quantum correlations may be related to deep facts in Number Theory. For instance, the Hardy-Littlewood constants [8] that characterize pairwise correlations among primes appear in the asymptotic expression of the reduced density matrix of the Prime state, for natural bi-partitions [9]. This fact and others that we shall present in this work can be interpreted as convincing evidence supporting the aforementioned statement. Quantum entanglement may indeed help unravel deep truths in Number Theory.
Let us mention a different line of research that brings together Arithmetics with Quantum Mechanics. The so called coprime chain [10] is a one dimensional strongly correlated system where the local degrees of freedom are labelled by integers, and have a nearest-neighbour interaction when they share a common divisor. Another many-body Hamiltonian that has been constructed recently uses interacting hard-core bosons on a ladder system [11]. Here the Prime state emerges as the ground state in a certain limit of the coupling constants. However, in our study we work directly with numbertheoretical states and not with Hamiltonians as in Refs. [10,11]. Another study that considers quantum states built upon sequences of integers can be found in [12].
In the present work, we report several advances in unveiling properties of the Prime state that relate to Arithmetics. First, we introduce the Quantum Fourier Transform (QFT) of the Prime state, which can be efficiently computed and allows to obtain Chebyshev-like biases in the prime number distribution from simple measurements in the computational basis. Then, we study the entanglement present in the Prime state and several of its quantum relatives, up to n = 30 qubits. We find that the scaling of the entropy for bi-partitions of size n 2 is almost identical, and almost maximal, for Prime states written in different qudit basis and for primes in arithmetic progressions. This actual scaling of the entropy is linear in n 2 with slope ∼ 0.88 . . . . We conjecture an analytical relation of this slope with the Shannon entropy of the density of square-free integers. A related con-jecture is posed regarding the intercept of the entropy. We provide an analytical approximation to the eigenvalues of the reduced density matrix of the Prime state, for natural bi-partitions. We also find a numerical approximation to the entropy of Prime states with any number of qubits and natural bi-partitions of any size. Finally, we compute the von Neumann entropy of novel numbertheoretical quantum states, defined from the sequences of square-free integers, odd composite numbers and starry primes, that we shall define below. For this study, we have developed an open-source library that diagonalizes matrices using floats of arbitrary precision [13].
The rest of the paper is organized as follows. Section II briefly reviews the Prime state; section III introduces its Quantum Fourier Transform; section IV describes the entanglement properties of the Prime state and other novel number-theoretical quantum states; lastly, section V summarizes the conclusions.

REVIEW OF THE PRIME STATE
The Prime state |P n of n qubits is defined as the uniform superposition of all the computational-basis states corresponding to prime numbers written in binary format, less than 2 n (we assume that n ≥ 2), i.e.
where π(x) is the prime-number counting function, which gives the number of primes smaller than or equal to x, p = p n−1 2 n−1 + · · · + p 1 2 1 + p 0 2 0 , and |p ≡ |p n−1 ⊗ · · · ⊗ |p 1 ⊗ |p 0 . This state was introduced in Ref. [1], where it was shown that its construction on a quantum computer is efficient. There are two ways of proceeding. The first simpler method is probabilistic; a second one makes use of Grover's algorithm. Let us review in more detail these two ways of constructing the Prime state.
The fundamental element to create a Prime state in either a probabilistic or deterministic way is to design a unitary operation U prime that discriminates primes from composites. This unitary acts as follows, where |k is a n-qubit computational-basis state, |a is a single ancilla qubit, and the action of the X gate is given by the Pauli matrix σ x . Note that if the state of the ancilla is |a = 1 √ 2 (|0 − |1 ), as done in Grover's subroutine, the action of U prime on a superposition is to introduce relative minus sign for prime numbers. The relevant observation here is that the above unitary amounts to code on a quantum computer a primality test, which is a language that belongs to the complexity class P . As a consequence, the operation U prime only involves a polynomial number of quantum gates, which will depend on the specific primality test chosen. For the sake of illustration, a detailed explicit form of a quantum primality oracle with O(n 6 ) gates, based on the classical Miller-Rabin test [14], was produced in Ref. [1]. Other primality algorithms, such as the AKS primality test [15], could be used as well.
We can then consider a first probabilistic way to create the Prime state. We need to apply the unitary U prime on the uniform superposition of all computational-basis vectors, |φ = 1 √ 2 n 2 n −1 i=0 |e i , plus an ancilla in the |0 state; where |φ is created by applying n Hadamard gates, H ⊗n , onto the initial |0 . . . 0 state. This yields: (3) By measuring the ancillary qubit in the above state, and post-selecting whenever the result of the measurement is |1 , one obtains the Prime state. This will occur with probability ∼ 1 n ln 2 , according to the Prime Number Theorem (see below). Because this probability decreases only polynomially with n, this is an efficient method to create the state.
A more refined, deterministic way to create the Prime state consists in using a Grover's algorithm that uses a primality test as oracle. The efficiency of the algorithm hinges on two facts. The first was mentioned previously, that is, the oracle based on U prime is efficient. The second fact that guarantees the efficiency in the construction of the Prime state is the relatively high abundance of primes, π(x), which is asymptotically given by the logarithmic integral function, This result is known as the Prime Number Theorem (PNT) [16], and it implies that the needed number of calls to Grover's oracle is only O( √ n). Indeed, this estimate is given by π 4 N M [17], where N = 2 n is the dimension of the Hilbert space of n qubits and M is the number of solutions to the oracle, i.e. M = π(2 n ). Therefore, the overall computational complexity of generating a Prime state of n qubits on a quantum computer with this method, using e.g. Miller-Rabin primality test, is O(n 6.5 ) [1].
As mentioned previously, once an oracle for the Prime state has been created, it can be used to numerically verify (or falsify) the Riemann hypothesis, which constitutes one of Clay Mathematics Institute's Millennium Problems [18]. The Riemann hypothesis states, for the distribution of prime numbers, that the deviations of π(x) from Li(x) are bounded asymptotically with x as [19] |π The quantum algorithm proposed in Ref. [1] allows to compute π(x) beyond the limits achievable by means of classical computation -once a fault-tolerant universal quantum computer is built-by using a quantum counting algorithm that provides an estimate of the number of terms in superposition in the Prime state. So far, π(x) have been computed up to 10 27 [4], which implies that a Prime state with a minimum of ∼ 90 logical qubits would be needed to surpass this computation (since 2 90 1.238 × 10 27 ).
To do so, it was suggested to apply a quantum counting algorithm [5] that delivers the number of solutions M to an oracle search problem; in this case, the number of primes that are marked by Grover's oracle, G. This algorithm uses Quantum Phase Estimation (QPE) [20] to obtain the eigenvalues of G, which in turn reveal the number of solutions to the query problem. In order to obtain an estimate of π(x) that is meaningful, the precision in the estimation should be smaller than the fluctuations allowed by the Riemann Hypothesis. To achieve such precision, O( √ 2 n ) calls to the oracle are needed. This is quadratically better than the performance of any classical counterpart in an identical oracular setting of the counting problem. Recently, there have been new proposals for quantum counting the solutions to an oracle G that do not rely on QPE [21]. In practice, any quantum counting algorithm may be applied.
For the sake of completeness, let us recall that the best classical algorithms that compute π(x) unconditionally, i.e. the validity of the estimation not depending on the truthfulness of the Riemann Hypothesis or any other unproven conjecture, do not however rely on enumerating all primes. The main algorithm for computing π(x) is a combinatorial method due to Meissel [22] and Lehmer [23], which has subsequently been improved by Lagarias, Miller and Odlyzko [24], Deléglise and Rivat [25], and Gourdon [26]. The latter version of this algorithm, which has time complexity O(x 2/3 ln −2 x), was used to compute π(10 27 ). Another two analytic algorithms for computing π(x), put forward by Lagarias and Odlyzko [27], have time complexity O(x 1/2+ ) and O(x 3/5+ ) respectively, with > 0.

QUANTUM FOURIER TRANSFORM OF THE PRIME STATE
The Quantum Fourier Transform (QFT) is a unitary operation that plays a key role in several important algorithms, such as Quantum Phase Estimation [20] or Shor's algorithm [28]. Its action on a quantum state |ψ = j x j |e j of n qubits, expressed in the computational basis, {|e j }, is given by x j e 2πi jk/2 n , where the {y k } are the discrete Fourier transform of the original amplitudes {x j }, and the i in the exponent is the imaginary unit. The QFT is an efficient subroutine, requiring O(n 2 ) gates on a quantum computer [17].
In the present work, we calculate the QFT of the Prime state, that is, Notice that the amplitudes in Eq. (7) are symmetric with respect to the central value k = N/2, where N = 2 n . This means that the probability of measuring the state |k and |2 n − k is the same. Indeed: P (2 n − k) = 1 2 n π(2 n ) p: prime e 2πi p e −2πi pk/2 n 2 (8) = 1 2 n π(2 n ) p: prime e −2πi pk/2 n 2 = P (k) . This is a general property of the QFT of a state with real amplitudes, that is verified by the Prime state, see Fig.  1.
The main result of the computation of |P n is that the QFT applied to the Prime state provides a method for estimating Chebyshev-like biases in the distribution of prime numbers. These biases reflect the unbalance in the number of primes, below a certain value, appearing in different arithmetic progressions.
Let us therefore consider arithmetic progressions of the form αk + β, for k = 0, 1, 2, . . . , with α, β coprimes (i.e. gcd(α, β) = 1). Dirichlet proved that there exists an infinite number of primes in any of these sequences [16]. The number of primes in the sequence αk + β (with α, β coprimes), below a certain value x, is denoted by the modular prime counting function π α,β (x), in close analogy to the prime counting function, π(x). The asymptotic behaviour of these modular counting functions is determined by where φ(α) is the Euler's totient function, which gives the number of coprimes to α smaller than α. For example, if p is a prime, one has φ(p) = p−1, because all the integers below p do not divide it. Notice as well that φ(α) is the number of arithmetic progressions of the form αk+β that can contain primes, which are those in which β is coprime to α. Thus, Eq. (9) means that the number of primes is, on average, asymptotically equi-distributed among the existing progressions αk + β (α, β coprimes), for fixed α. Nevertheless, when considering a finite natural sequence of primes, there exists biases in the distribution of the primes among different progressions that are quantified by the functions These biases can be numerically estimated, to a desired precision , from the probability peaks appearing in |P n . In particular, relevant peaks appear at values k = N/3, N/4, N/6, and their mirror images about the central peak (see Appendix A for a derivation of these results). The expressions for these peaks are, where ∆(N ) ≡ π 4,3 (N ) − π 4,1 (N ) is the Chebyshev bias. There exists a direct relation between the functions π 3,1 (x), π 3,2 (x) and π 6,1 (N ), π 6,5 (N ). In fact, π 3,1 (x) = π 6,1 (x) and π 3,2 (x) = π 6,5 (x) + 1. This is so because every prime p counted by π 3,1 (x) is of the form p = 3k + 1, but since p is odd, k must be even (k = 2k ), and p = 6k + 1, so π 3,1 (x) = π 6,1 (x). On the other hand, every prime p counted by π 3,2 (x) is of the form p = 3k + 2, with k odd (k = 2k + 1), so p = 6k + 5; in this case, 2 is the only prime counted by π 3,2 (x) and not by π 6,5 , and thus π 3,2 (x) = π 6,5 (x) + 1. These relations, allow to obtain π 3,1 (x), π 3,2 (x), π 6,1 (N ), π 6,5 (N ), ∆ 3; 2,1 (N ) and ∆ 6; 5,1 (N ), by substituting Eq. (12) into the probability peaks P (N/3) and P (N/6) in Eq. (11), and solving the resulting system of two equations and two unknowns. Of course, π(N ) must be known in advance, but this is already provided by the Prime state combined with some quantum counting algorithm. Notice as well that N/3 and N/6 do not correspond to integer values for qubit systems, where N = 2 n . This means that the associated peaks cannot be directly measured in the computational basis. Instead, one needs to resort to the peak at the closest integer value, which will be N/3 ± 1 and N/6 ± 1, depending on N .
It is important to consider how efficient it is to measure the peaks in Eq. (11), to a desired precision, when n grows large. Note that P (N/3) and P (N/6) are dominated by the cross product of the modular prime counting functions, divided by N π(N ). According to Eq. (4) and Eq. (9), these probabilities will then decrease as ∼ 1 φ(α) 2 1 n ln 2 , where α = 3, 6 respectively. This decrease is only polynomial in n, and therefore, the method is efficient. This is so because the number of measurements required to achieve a target precision in the estimation of outcome probabilities, {p i }, is dictated by the sampling process of a multinomial distribution. In this context, the statistical additive errors [31], { i }, arise from finite sampling, and depend on the number of measurements performed, #M , as i = pi(1−pi) #M . Hence, in order to achieve a small multiplicative error ε i , say ε i ≈ 0 and i << p i , then the number of measurements, or runs of the quantum circuit, must be #M >> 1−pi pi . Therefore, it is clear that estimating probabilities with a small multiplicative error remains efficient as long as these probabilities show at most a polynomial decrease with increasing n.
Moreover, direct application of the QFT on the Prime state, followed by the procedure of Amplitude Estimation [32], allows to obtain a dependency of the additive error on the number of runs of the quantum circuit that is O( 1 #M ), i.e. a quadratic improvement compared to direct sampling from the multinomial distribution. This is the same technique employed for instance in the quantum Monte Carlo algorithm [33]. In this case, all it takes to prepare |P n for Amplitude Estimation is to apply a single n + 1-qubit controlled-X gate on an ancilla qubit, that singles out the desired peak, e.g. |N/3 , States like the one in Eq. (13) are ready for Amplitude Estimation. Classically, the best algorithm for computing modular prime counting functions, π α,β (x), is a variation of the Meisser-Lehmer method for computing π(x), and has time complexity O(x 2/3 ln −2 x) [34]. To the best of our knowledge, modular prime counting functions have been computed for values only up to x = 10 9 [35]. Hence, a quantum computer could be helpful in this sort of computations.
In contrast to P (N/3) and P (N/6), the peak P (N/4), which delivers the Chebyshev bias, ∆(N ), is unlikely to be of practical use because the denominator N π(N ) is expected to cause an exponential decrease in this probability, as the number of qubits n increases. As already explained, this would demand an exponentially large number of measurements. Yet, it is interesting to appreciate how the amplitudes of the Prime state interfere under the QFT to give some remarkable meaning, in terms of number-theoretical functions, to certain values of k. Note as well that if one is willing to directly sample from |P n , all the information contained in the probability peaks is extracted simultaneously, rather than sequentially, by accumulation of statistical knowledge on measurement outcomes. This means that the values of all peaks are estimated altogether with a fixed number of samples, to a desired precision .
It also happens that the QFT of an equally-weighted superposition provides a means to count the number of terms in the computational basis. Indeed, the probability of finding the |0 . . . 0 state on a measurement after the application of the QFT equals M/N , where again N = 2 n is the dimension of the Hilbert space of n qubits and M is the number of terms in the superposition (in the computational basis). A derivation of this result can be found in Appendix B. For |P n , this probability reads which decreases as ∼ 1 n ln 2 in the limit n 1, according to the PNT, Eq. (4). The prime counting function, π(N ), appears as well in the central peak, where k = N/2, as shown in Appendix A. The expression for this peak is which also decreases linearly in n −1 .
Because the precision required to meaningfully test Riemann's hypothesis, Eq. (5), implies that the multiplicative error must become smaller and smaller for increasing x (as √ x ln x π(x) decreases), direct sampling from |P n is not competitive with previously mentioned quantum counting algorithms, requiring O(2 n n −3 ) repetitions. On the other hand, it is possible to prepare a state like that in Eq. (13), but marking the peak P (0) instead, and then apply Amplitude Estimation. However, note that if one has an oracle, it is possible to generate, with a single query to that oracle, the state in Eq. (3). Such state is also ready for AE, and can be used to estimate π(N ) = M/N , but it is easier to prepare than the one in Eq. (13). So there is no point in using |P n for estimating π(2 n ).
The QFT may nonetheless be useful to estimate the number of terms in an equally-weighted superposition, for states for which one does not know of any oracle; for instance, a ground state of a Hamiltonian obtained from a variational algorithm such as the Variational Quantum Eigensolver [29] or the Quantum Approximate Optimization Algorithm [30].

ENTANGLEMENT OF NUMBER-THEORETICAL QUANTUM STATES
The entanglement properties of number-theoretical quantum states can be studied by taking different bi-partitions of the states and quantifying their entanglement, with figures of merit such as the von Neumann entropy, the purity or the entanglement spectrum, among many possibilities. Unfortunately, a practical way to quantify genuine multipartite entanglement, not related to partitions of the system, is not available for states with a large number of qubits [36].
The entanglement present in number-theoretical quantum states is the result of correlations among the digits of the numbers appearing in the superposition. In turn, these observables may reveal information about pairwise correlations between the numbers in the sequence. For instance, the reduced density matrix of the Prime state upon taking a natural bi-partition A-B with the first (i.e. least significant) m qubits [9], ρ A , is asymptotically given by ρ A [7], whose expression is (17) where C(h) are the Hardy-Littlewood constants, defined by the twin prime constant and {p} the prime numbers. The constants C(h) characterize the pairwise correlations among prime numbers, and it is precisely due to these correlations that the reduced density matrix ρ A is not maximally mixed. To be precise, the first Hardy-Littlewood conjecture [8] assures that the number of prime pairs (p, p + k) up to a certain number x, i.e. π k (x), is given by where An analytical approximation to the positive eigenvalues of the matrix C m was conjectured in Ref. [7]. That approximation was suitable for the highest eigenvalues, but failed to describe the lowest ones. Moreover, the degeneracies of the former failed beyond some point [37]. A recent article gives a heuristic derivation of the Hardy-Littlewood conjecture (19) starting from the pair correlation formula for the Riemann zeros, including lower order terms [38]. The reverse statement was known since long and played an important role in the connection between Number Theory and Quantum Chaos [39,40]. The derivation of Ref. [38] employs the following expression for the Hardy-Littlewood constants, where φ(x) is the Euler's totient function, µ(x) is the Möbius function (see below), and are the Ramanujan's sums. Equation (20) is the Ramanujan-Fourier series of the constants C(h) and has also appeared in the context of random processes, more concretely, in the relation between the power spectrum and the correlation function using the Wiener-Khintchine formula [41,42].
Based on Eq. (20), a more precise description of the eigenvalues of C m can be obtained (for a justification see Appendix C). To every square-free odd integer k, that is, to every odd integer k such that it does not contain any prime raised to a power larger than 1 in its prime decomposition, we associate an eigenvalue γ k of C m with degeneracy φ(k). These eigenvalues are given by such an accidental degeneracy is γ 13 = γ 21 = 1 144 , which has a total degeneracy φ(13) + φ(21) = 24 (we denote by γ k = µ 2 (k)/φ 2 (k)). Another example is γ 35 = γ 39 = 1 576 , that has a degeneracy φ(35) + φ(39) = 48.
Using Eq. (22), the eigenvalues of the density matrix (16) are approximated by with a degeneracy φ(k). The good agreement between the exact eigenvalues and those computed from Eq. (23) for n = 30 qubits, is shown in Fig. 2, where the lowest 72 entanglement energies of the entanglement spectrum are plotted. These correspond to the highest eigenvalues of the density matrix of the Prime state, since the entanglement spectrum {ε k } of a density matrix ρ with eigenvalues {λ k } is defined as Using again Eq. (22), we can estimate the contribution of the positive eigenvalues to the trace of the powers of C m , from which all Rényi entropies can be readily derived [43]. The computation is as follows, where we have used that for any function f (k) which is multiplicative, it holds that where the product runs over the odd prime numbers. We have also used that φ(p) = p − 1 for prime numbers, and that µ(k) 2 = |µ(k)|. In Eq. (25) we dropped the negative term proportional to 1/φ m of Eq. (22) because its contribution vanishes in the limit m 1 for s > 1. For s = 2, Eq. (25) coincides with the heuristically-derived asymptotic expansion of Tr C 2 m , and was conjectured to hold for any s ≥ 2 in Ref. [7]. A direct connection thus exists between Eq. (22) and the latter conjecture. For s = 1, the Riemann zeta function ζ(1) arises and the product diverges because we have dropped the negative eigenvalues of the matrix C m . Figure 3 shows the ratio of the exact value of the trace of the matrix C s m to the asymptotic analytical formula Eq. (25), and the exponential convergence of this ratio towards 1, up to n = 32 qubits. These results indicate that the contribution of the negative eigenvalues of C m to the trace of the powers of this matrix is asymptotically negligible, for s ≥ 2. In contrast, their contribution to the von Neumann entropy does not appear to be so.
The von Neumann entropy S of a density matrix ρ, with eigenvalues {λ i }, is defined by and it constitutes a relevant measure of bipartite entanglement. The von Neumann entropy of the reduced density matrices of equal-sized bi-partitions of the Prime state, with an even number of qubits, was numerically calculated in Ref. [7], and found to scale linearly with the size of the bi-partition. To be precise, the best fit to a line, for n = 20 -30 qubits and for the natural equal-sized bi-partition, is S(n) = 0.88612902 n 2 − 1.30405956. This result indicates that the Prime state is highly entangled, but not maximally so, because the maximal possible scaling for the entropy would be linear in n 2 with a slope equal to 1. Notice as well that this implies that it is not possible to apply standard techniques such as Matrix Product States (MPS), or other more-refined tensor networks, to efficiently simulate the Prime state on a classical computer [44].
We herein conjecture that the entropy of the natural equal-sized bi-partition of the Prime state of n qubits is asymptotically given by is the Shannon entropy and 3/π 2 = 1/(2 ζ(2)) is equal to a half of the asymptotic density of odd square-free integers [16]. The constant 3/π 2 also appears in the study of topological dynamical systems. In particular, it has been shown that it gives half the topological entropy of the square-free flow [45,46]. The term -1 in the intercept comes from the fact that the least relevant qubit is basically in the state |1 , because all primes but 2 are odd, and hence this qubit does not contribute to the entropy in the asymptotic limit. In order to test the conjectured Eq. (28), we have computed the entropy of the Prime state up to n = 30 qubits. Diagonalization using quadruple-precision floating-point numbers was implemented to meet the precision demanded, for large values , respectively. An exponential decrease in δ is observed before reaching a plateau. of n, by the observed oscillatory behaviour of the slope c π . To the best of our knowledge, this is the first opensource library that diagonalizes matrices using float128; as a matter of fact, this library works in arbitrary precision. The code is publicly available on GitHub, together with all the numerical results presented in this paper [13].
up to n = 30 qubits. Oscillatory convergence towards the predicted value H 3 π 2 = 0.886082085 . . . is observed, in good agreement with Eq. (28). The value of the intercept γ, obtained as is shown in Fig. 4-bottom, up to n = 30 qubits. Oscillatory convergence and good agreement with the predicted value 1 + 3 π 2 = 1.30396355 . . . is observed as well. As a matter of fact, the behaviour of the convergence is almost identical in both cases, with an exponential approximation toward the conjectured values as the number of qubits grows, as shown in Fig. 5. The remarkably-similar fashion in which c π and γ converge may be seen as an indication that both conjectures are closely related.
We have also extended the computations to the natural equal-sized bi-partition of size n 2 , with odd n up to n = 29 qubits, see Fig. 6-top. The scaling of the entropy  Table I. This is equivalent to expanding the prime numbers in the ternary basis, quaternary basis, etc. and putting them in superposition. We find that the scaling is independent of the basis, i.e. linear with slope ∼ 0.88 . . . , so in this respect it is universal. This result should be expected, because otherwise there would exist in some sense a 'preferred' basis for expressing prime numbers. It is also truly noteworthy that we have found the same scaling (and almost the same values) for the Von Neumann entropy of the natural equal-sized bi-partitions of the QFT of the Prime state, see Fig. 6-top and Table  I.
This is precisely the behaviour observed when plotting the l.h.s of Eq. (31) against m/n, see Fig. 7-bottom. This realization allows for a numerical approximation to f (m/n) as a Fourier series, see Appendix D for the precise values of the first few coefficients. It follows that an estimate of S(m, n) may be obtained from this series for any value of n, m.
Other number-theoretical quantum states apart from the Prime state will be defined below, and their entanglement traits studied. Table I   i.e.
with α, β coprime numbers. When α is a power of 2, a completely analogous derivation to that of Eq. (16) leads to an asymptotic expression for the reduced density matrix of arithmetic Prime states upon taking a natural bi-partition A-B with the first (i.e. least significant) m qubits, where again d = 2 m−1 , l N ≡ Li2(N )
We have found that the scaling of the entropy of arithmetic Prime states, |P n; α,β , when α is a power of 2, is once again dictated by a straight line with slope ∼ 0.884, see Fig. 8 and Table I. The intercept, however, is different for distinct values of α. It is also independent of β. The latter fact is easily understood from the asymptotic expression of the reduced density matrix ρ A (α), which does not depend on β. This is ultimately so because prime numbers are equally distributed on average among the different arithmetic progressions, for a given α, as expressed in Eq. (9). Hence, we label this constant as γ α . We have found the following behaviour for γ α when α is a power of 2, γ α 1 + s 3 π 2 , α = 2 k , s = 1, 3, 5, . . . , (36) where s is the k-th odd square-free number. Figure 9 shows the independence of γ 4 and γ 8 from β, for a large number of qubits, and the convergence towards a value close to that in Eq. (36). Using the results obtained for γ 4 , γ 8 , γ 16 and γ 32 for n = 30 qubits, averaged among the different progressions, the value s = (γ α − 1) π 2 3 is found to be s = 2.9137 . . . that, for Prime and arithmetic Prime states of n qudits, the von Neumann entropy S of equal-sized bi-partitions is asymptotically given by where x is the floor function and γ α; d is a constant that depends upon the specific arithmetic progression used to define the state, as well as the qudit basis in which the numbers are expressed; in addition, it also depends on whether n is even or odd.
One may ask whether it could be possible to construct a state that is the uniform superposition of odd composite (i.e. non-prime) numbers, and whether or not this state holds large entanglement, in the hope of finding a shortcut to the distribution of primes by looking at the distribution of composites. This state is defined by |/ P n ≡ 1 C c,n 2 n −1 c: odd c: composite |c , C c,n = 2 n−1 − π(2 n ) . (38) We have computed the von Neumann entropy for the natural equal-sized bi-partition of such odd composite state, up to n = 28 qubits, and found that this entropy remains approximately constant around a mean value of ∼ 2, irrespective of the number of qubits n, see Fig.  10. This behaviour is reminiscent of finitely-correlated systems away from criticality [47]. A slight decline in the entropy can be observed as n increases, though. An asymptotic expression for the reduced density matrix of the odd composite state upon taking a bi-partition A-B with the first (i.e. least significant) m qubits, can be obtained, as shown in Appendix E. It is given, when n, m → ∞, by where From Eq. (39), the von Neumann entropy S in the asymptotic limit is given by in agreement with the observed numerical results. One might hope that, because the entanglement entropy of the odd composite state is small, there could indeed be a shortcut to efficiently compute the prime number distribution by using tensor networks. But this is incorrect, as can be argued readily. Most numbers are composite, and the composite state carries their equal superposition. This very high density of equally-weighted superpositions approaches the product state |φ = 1 √ 2 n 2 n −1 i=0 |e i . Thus, tensor networks of low ancillary dimension only describe trivial properties of the composite state. To learn something about primes from the composite state, we need to go to the subleading terms in any figure of merit, which implies exponential effort in the tensor network description.
Another number-theoretical quantum state, relative of the Prime state, is the odd square-free state, i.e. the uniform superposition of odd square-free numbers, which are odd integers that do not contain any prime raised to a power larger than 1 in its prime decomposition. This state is defined by  We have computed the von Neumann entropy of this state, up to n = 28 qubits, and found that it follows a similar behaviour to that of the odd composite state, although now the entropy shows a slight increase as the number of qubit grows (see Fig. 10). The crossing point of both entropies is located at n ∼ 20 qubits.
Another possibility is to define a uniform superposition of square-free numbers and let the signs of the amplitudes be given by the Möbius function. This state, that we shall call Möbius state [7], is then defined by (43) The von Neumann entropy of the natural equal-sized bipartition of the Möbius state scales linearly with a slope close to 1. This means that the Möbius state is maximally entangled, up to a constant (see Fig. 10 and Table  I). Note that the appearance of negative signs in the coefficients of the superposition allows for the appearance of almost maximal entanglement. This shows that there are two ways of getting maximal entanglement. One is to play with random coefficients in the quantum state. Natural cancellations make the reduced density matrix very mixed. Alternatively, we may impose that all coefficients are zero or one (up to normalization). Then, maximal entanglement is related to the thin distribution of numbers present in the state.
Finally, we have computed the entropy of the natural equal-sized bi-partition of starry -Prime states. Starry primes are integer numbers, such that the n-th starry prime, p * n , satisfies where p n is the n-th prime. These numbers are remarkable in that a ζ * function, defined via the Euler product formula, i.e.
has an analytic extension with exactly the same zeros, and the same multiplicities, than the original Riemann zeta function [48]. The starry-Prime state is then defined by where {p * } is any set of starry primes up to 2 n − 1. We have found that the entropy of the natural equal-sized bipartition of random starry-prime states -where the terms entering into the superposition Eq. (46) were taken at random according to Eq. (44)-, scales almost in a maximal way, i.e. linearly in n 2 with slope ∼ 1 (see Fig. 11 and Table I). This result shows that primes are less random than starry primes. The extra element of randomness in the choice of numbers within intervals defined by primes is sufficient to raise the slope of the entanglement entropy to a value near 1.

CONCLUSIONS
Prime numbers are classical mathematical objects. However, using a quantum computer, they can be efficiently placed in a linear superposition, the Prime state. This cast of a classical sequence onto a quantum state brings a novel quantum perspective to Number Theory that we are just beginning to explore.
In this work, we have shown how Chebyshev-like biases arise from the Quantum Fourier Transform of the Prime state, and how to estimate them from measurements. We have also delved into the entanglement properties of the Prime state, which turn out to be intimately related to the correlations between pairs of prime numbers. These correlations were postulated by Hardy and Littlewood almost a century ago, and played a critical role in the connection between Number Theory and the Theory of Quantum Chaos based on heuristic analogies. According to it, prime numbers are classical objects (labelling periodic orbits in a chaotic dynamical system), while the Riemann's zeros are quantum objects (eigenenergies of a hypothetical Riemann Hamiltonian). In our approach, the Hardy and Littlewood conjecture provides us with the entanglement Hamiltonian of the Prime state, whose spectrum we have found to be characterized by the odd square free integers. Furthermore, the entanglement entropy of the Prime state is closely related to the Shannon entropy of half the density of the square free integers. The previous result also applies to the quantum superposition of prime numbers in arithmetic progressions, which is a kind of quantum version of Dirichlet's theorem about these sequences. Likewise, we have verified that the entanglement entropy of the Prime state does not depend on the basis used, i.e. qubits, qutrits, etc., which is an expected result but that confirms its universality. We have furthermore characterized entanglement properties of other quantum states built upon different sequences of numbers. From a more practical point of view, we have developed an open-source library that diagonalizes matrices using floats of arbitrary precision.
We expect this novel quantum approach to Number Theory to be fruitful in the future, as it may help us deepen into the profound mysteries of numbers.

B. Quantum Fourier Transform of a uniform superposition
Let us consider the definition of the QFT of a general quantum state |ψ of n qubits, expressed in the computational basis, {|e j }, where and the i in the exponent is the imaginary unit. The probability of measuring the |0 . . . 0 state in Eq. (51) in the computational basis thus is: If the state |ψ is a uniform superposition of M computational-basis vectors, i.e. x j = 1/ √ M for all j such that x j = 0, and N = 2 n is the dimension of the Hilbert space, then Equation (54) implies that the QFT can be used as an efficient method for counting the number of terms in a uniform superposition -for it takes only O(n 2 ) gates to implement it on a quantum computer [17]-, as long as M/N and 1 − M/N remain large enough, i.e. do not decrease exponentially with the number of qubits, for this would require in practice an exponential precision in the estimation and therefore an exponential number of samples/repetitions of the protocol.
C. The Hardy-Littlewood-Ramanujan density matrices As explained in the main text, the reduced density matrix of the Prime state with n qubits can be approximated bȳ where C m is a matrix constructed with the even values of the Hardy-Littlewood constants C(2k), In order to justify the conjecture (22) for the eigenvalues of C m , we define the density matrix with where are the Ramanujan's sums. Replacing Eq. (59) into Eq. (58) gives which leads us to write C D as a sum of matrices, where The sums c k (n) are multiplicative functions in k, namely In particular, where we have used that c 2 (h) = cos(πh). Equation (65) implies that The Möebius and Euler's totient functions are also multiplicative, hence where we have used that µ(2) = −1, φ(2) = 1. From the latter equations we find that the sum over the even integers k in Eq. (62) is equal to the sum over the odd integers, and thus Our aim is to find the eigenvalues of the matrix C D . This is a difficult problem in general. However, we can find a good approximation performing some simplifications. First of all, we choose D as the product of the first L odd prime numbers, that is where p 2 = 3, p 3 = 5 . . . . Second, we truncate the sum in Eq. (68) to the integers k that divide D, obtaining the matrixC We shall show below that the eigenvalues of this matrix provide a good approximation to those in Eq. (68). Using Eq.(63) and Eq. (60) one can write R k,D as where I D is the identity matrix of dimension D, and The negative term in Eq. (71) guarantees that (R k,D ) j,j = 0, because k l=1 gcd(l,k)=1 Replacing Eq. (71) into Eq. (70) gives where we used the relation [16] d|n It turns out that the matricesR k,D commute among themselves, so they can be diagonalized simultaneously and, using Eq. (74), we can obtain the spectrum ofC D . This fact follows from the following periodicity property of the Ramanujan sums, i.e.
if k|D =⇒ c k (2(h + D)) = c k (2h) , that can be easily proved using Eq. (60). Hence, the eigenvalues, r k,D (a) (a = 1, . . . , D), ofR k,D are given by the Fourier transform of its entries, Let us first notice that ω l,a is a D root of unity (use that k|D), We have to consider two cases, when ω l,a = 1 and when ω l,a is a primitive root of unity. The former case occurs if a is given in terms of l as Taking into account that gcd(l, k) = 1, there are φ(k) eigenvalues of this type. In the case where ω l,a is a primitive root of unity, one has D h ω h l,a = 0, and then The number of vanishing eigenvalues is therefore D − φ(k). This proves that the matricesR k,D are diagonal in the Fourier variables a. This result, together with Eq. (74), will allow us to obtain the eigenvalues ofC D . Indeed, for a given a = 1, . . . , D, the eigenvalue r k,D (a) can be either D or 0 depending on whether a satisfies condition (79) or (80). The latter conditions have to be analyzed for all values of k|D and their contribution added in Eq. (74) with the corresponding weight (µ(k)/φ(k)) 2 .
Recall that k i (i = 1, 2) are odd integers and so is k 1 k 2 . On the other hand, the l.h.s. of that equation is even, which shows that this condition cannot be fulfilled. We are left with the case l 1 k 2 = l 2 k 1 . Since l 1 does not divide k 1 there must exists an integer x 1 such that l 1 = l 2 x 1 . Similarly, l 2 = l 1 x 2 for some x 2 . Using these equations we find which means that the decomposition a = 2lD/k (modD) is unique, as stated above. Table II shows   In summary, the spectrum of the matrixR k,D is given by where the first column denotes the eigenvalue and the second one its degeneracy. Using Eq. (71) we obtain the spectrum of R k,D , .
which shows explicitly that it is a traceless matrix. Finally, using Eq. (74), together with the fact that only one matrixR k,D contributes with D to the total sum, we obtain the eigenvalues of the matrixC D , , k|D which has degeneracy φ(k). We have simplified the notation using the fact that µ(k) 2 = 1 when k|D. Adding these degeneracies over all the values of k that divide D, we recover the total dimension of the matrixC D , using the relation [16] k|D φ(k) = D .
The matrixC D was introduced in order to obtain a simple analytic approximation to the eigenvalues γ i,D of the matrix C D . Fig. 12 shows the eigenvalues of both matrices for D = 15. Even for this small value of D, they are numerically close and follow the degeneracies given by φ(k), which in this case are 1, 2, 4, 8. Equation (87)  that yields where B(k m ) ≡ km k=1(odd) The asymptotic behaviour of this sum is given by (see Fig. 13-right) with β = γ + 1 2 ln 2 + p ln p p(p − 1) = 1.679135304 . . . , where γ is the Euler's constant. This result can be derived from the equation [1] n k=1 µ 2 (k) φ(k) ln k + γ + where n = 2m. Fig. 14-left shows the linear behaviour of S(n) with a slope that seems to converge to the value H(3/π 2 ), as shown in Fig. 14-right.