# The Prime state and its quantum relatives

Diego García-Martín  Barcelona Supercomputing Center (BSC), Barcelona, Spain.
Dept. Física Quàntica i Astrofísica, Universitat de Barcelona, Barcelona,
Spain.  Instituto de Física Teórica, UAM-CSIC, Madrid, Spain.  Eduard Ribas
Dept. Física Quàntica i Astrofísica, Universitat de Barcelona, Barcelona,
Spain.  Stefano Carrazza  TIF Lab, Dipartimento di Fisica, Università degli
Studi di Milano and INFN Milan, Milan, Italy.  Quantum Research Centre,
Technology Innovation Institute, Abu Dhabi, UAE.  José I. Latorre  Quantum
Research Centre, Technology Innovation Institute, Abu Dhabi, UAE.  Dept.
Física Quàntica i Astrofísica, Universitat de Barcelona, Barcelona, Spain.
Center for Quantum Technologies, National University of Singapore, Singapore.
Germán Sierra  Instituto de Física Teórica, UAM-CSIC, Madrid, Spain.

###### Abstract

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.

##  1 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 $\pi(x)$, whose fluctuation around its average 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 $\Delta(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 non-
locality, 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 number-theoretical 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 number-theoretical 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 $\frac{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 $\frac{n}{2}$
with slope $\sim 0.88\dots$. We conjecture an analytical relation of this
slope with the Shannon entropy of the density of square-free integers. A
related conjecture 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 number-theoretical 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 2 briefly reviews the
Prime state; section 3 introduces its Quantum Fourier Transform; section 4
describes the entanglement properties of the Prime state and other novel
number-theoretical quantum states; lastly, section 5 summarizes the
conclusions.

##  2 Review of the Prime state

The Prime state ${|\mathbb{P}_{n}{\rangle}}$ 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\geq 2$), i.e.

${|\mathbb{P}_{n}{\rangle}}\equiv\frac{1}{\sqrt{\pi(2^{n})}}\sum_{p:\,prime}^{2^{n}}|p{\rangle}\,,$
(1)

where $\pi(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}+\cdots+p_{1}2^{1}+p_{0}2^{0}$, and
$|p{\rangle}\equiv|p_{n-1}{\rangle}\otimes\cdots\otimes|p_{1}{\rangle}\otimes|p_{0}{\rangle}$.
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,

$U_{prime}\,|k{\rangle}\otimes|a{\rangle}\equiv\left\\{\begin{array}[]{ll}|k{\rangle}\otimes
X\,|a{\rangle}&\;\,\textrm{if $k$ is prime}\\\ \\\
|k{\rangle}\otimes|a{\rangle}&\;\,\textrm{if $k$ is not
prime}\end{array}\right.\,,$ (2)

where $|k{\rangle}$ is a $n$-qubit computational-basis state, $|a{\rangle}$ is
a single ancilla qubit, and the action of the $X$ gate is given by the Pauli
matrix $\sigma_{x}$. Note that if the state of the ancilla is
$|a{\rangle}=\frac{1}{\sqrt{2}}(|0{\rangle}-|1{\rangle})$, 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,
$|\phi{\rangle}=\frac{1}{\sqrt{2^{n}}}\,\sum_{i=0}^{2^{n}-1}|e_{i}{\rangle}$,
plus an ancilla in the $|0{\rangle}$ state; where $|\phi{\rangle}$ is created
by applying $n$ Hadamard gates, $H^{\otimes n}$, onto the initial $|0\dots
0{\rangle}$ state. This yields:

$\sqrt{\frac{2^{n}-\pi(2^{n})}{2^{n}}}\sum_{\overline{p}:\,not\,prime}|\overline{p}{\rangle}|0{\rangle}+\sqrt{\frac{\pi(2^{n})}{2^{n}}}\sum_{p:\,prime}|p{\rangle}|1{\rangle}\,.$
(3)

By measuring the ancillary qubit in the above state, and post-selecting
whenever the result of the measurement is $|1{\rangle}$, one obtains the Prime
state. This will occur with probability $\sim\frac{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, $\pi(x)$, which is asymptotically given
by the logarithmic integral function,
$Li(x)\equiv\int_{2}^{x}\frac{dt}{\ln\,t}$, i.e.

$\pi(x)\sim Li(x)\;\xrightarrow{x\rightarrow\infty}\;\frac{x}{\ln\,x}\,.$ (4)

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(\sqrt{n})$.
Indeed, this estimate is given by $\frac{\pi}{4}\,\sqrt{\frac{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=\pi(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 $\pi(x)$ from $Li(x)$ are bounded asymptotically with
$x$ as [19]

$|\pi(x)-Li(x)|\leq\frac{1}{8\pi}\sqrt{x}\ln{x}\qquad\textrm{for $x>2657$}\,.$
(5)

The quantum algorithm proposed in Ref. [1] allows to compute $\pi(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, $\pi(x)$ have been computed up to $10^{27}$ [4], which implies
that a Prime state with a minimum of $\sim 90$ logical qubits would be needed
to surpass this computation (since $2^{90}\simeq 1.238\times 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 $\pi(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(\sqrt{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 $\pi(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 $\pi(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 $\pi(10^{27})$. Another two analytic algorithms for computing
$\pi(x)$, put forward by Lagarias and Odlyzko [27], have time complexity
$O(x^{1/2+\epsilon})$ and $O(x^{3/5+\epsilon})$ respectively, with
$\epsilon>0$.

##  3 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
$|\psi{\rangle}=\sum_{j}x_{j}\,|e_{j}{\rangle}$ of $n$ qubits, expressed in
the computational basis, $\\{|e_{j}{\rangle}\\}$, is given by

$QFT\,|\psi{\rangle}=\sum_{k=0}^{2^{n}-1}y_{k}\,|e_{k}{\rangle}\;,\quad
y_{k}\equiv\frac{1}{\sqrt{2^{n}}}\sum_{j=0}^{2^{n}-1}x_{j}\,e^{2\pi
i\,jk/2^{n}}\,,$ (6)

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,
$|\hat{\mathbb{P}}_{n}{\rangle}\equiv QFT\,{|\mathbb{P}_{n}{\rangle}}$, that
is,

$|\hat{\mathbb{P}}_{n}{\rangle}=\frac{1}{\sqrt{2^{n}\pi(2^{n})}}\sum_{k=0}^{2^{n}-1}\left(\sum_{p:\,prime}e^{2\pi
i\,pk/2^{n}}\right)|k{\rangle}\,.$ (7)

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{\rangle}$ and $|2^{n}-k{\rangle}$ is the same. Indeed:

$\displaystyle P(2^{n}-k)$ $\displaystyle=$
$\displaystyle\frac{1}{2^{n}\,\pi(2^{n})}\left|\sum_{p:\,prime}e^{2\pi
i\,p}\,e^{-2\pi i\,pk/2^{n}}\right|^{2}$ $\displaystyle=$
$\displaystyle\frac{1}{2^{n}\,\pi(2^{n})}\left|\sum_{p:\,prime}e^{-2\pi
i\,pk/2^{n}}\right|^{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.

Figure 1: Numerical simulation of the Quantum Fourier Transform of a Prime
state with $n=15$ qubits. The probability peaks appear symmetrically
distributed with respect to the central value $k=2^{n}/2$. Note as well that
there is no peak for the value $k=2^{n}$, because the state $|2^{n}{\rangle}$
is not included in the Hilbert space of $n$ qubits.

The main result of the computation of $|\hat{\mathbb{P}}_{n}{\rangle}$ 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 $\alpha
k+\beta$, for $k=0,1,2,\dots$, with $\alpha,\beta$ coprimes (i.e.
$\gcd(\alpha,\beta)=1$). Dirichlet proved that there exists an infinite number
of primes in any of these sequences [16]. The number of primes in the sequence
$\alpha k+\beta$ (with $\alpha,\beta$ coprimes), below a certain value $x$, is
denoted by the modular prime counting function $\pi_{\alpha,\beta}(x)$, in
close analogy to the prime counting function, $\pi(x)$. The asymptotic
behaviour of these modular counting functions is determined by

$\pi_{\alpha,\beta}(x)\sim\frac{1}{\phi(\alpha)}\,Li(x)\;\xrightarrow{x\rightarrow\infty}\;\frac{1}{\phi(\alpha)}\frac{x}{\ln\,x}\,,$
(9)

where $\phi(\alpha)$ is the Euler’s totient function, which gives the number
of coprimes to $\alpha$ smaller than $\alpha$. For example, if $p$ is a prime,
one has $\phi(p)=p-1$, because all the integers below $p$ do not divide it.
Notice as well that $\phi(\alpha)$ is the number of arithmetic progressions of
the form $\alpha k+\beta$ that can contain primes, which are those in which
$\beta$ is coprime to $\alpha$. Thus, Eq. (9) means that the number of primes
is, on average, asymptotically equi-distributed among the existing
progressions $\alpha k+\beta$ ($\alpha,\beta$ coprimes), for fixed $\alpha$.
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

$\Delta_{\alpha;\,\beta_{1},\beta_{2}}(x)\equiv\pi_{\alpha,\beta_{1}}(x)-\pi_{\alpha,\beta_{2}}(x)\,.$
(10)

These biases can be numerically estimated, to a desired precision $\epsilon$,
from the probability peaks appearing in $|\hat{\mathbb{P}}_{n}{\rangle}$. 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,

$\begin{split}P(N/3)&=\frac{\Delta_{3;2,1}(N)^{2}+\pi_{3,1}(N)\pi_{3,2}(N)-\pi(N)+2}{N\pi(N)}\,,\\\
\\\ P(N/4)&=\frac{1+\Delta(N)^{2}}{N\pi(N)}\,,\\\ \\\
P(N/6)&=\frac{\Delta_{6;5,1}(N)^{2}+\pi_{6,1}(N)\pi_{6,5}(N)-3\,\pi_{6,1}(N)+3}{N\pi(N)}\,,\end{split}$
(11)

where $\Delta(N)\equiv\pi_{4,3}(N)-\pi_{4,1}(N)$ is the Chebyshev bias. There
exists a direct relation between the functions $\pi_{3,1}(x)$, $\pi_{3,2}(x)$
and $\pi_{6,1}(N)$, $\pi_{6,5}(N)$. In fact, $\pi_{3,1}(x)=\pi_{6,1}(x)$ and
$\pi_{3,2}(x)=\pi_{6,5}(x)+1$. This is so because every prime $p$ counted by
$\pi_{3,1}(x)$ is of the form $p=3k+1$, but since $p$ is odd, $k$ must be even
($k=2k^{\prime}$), and $p=6k^{\prime}+1$, so $\pi_{3,1}(x)=\pi_{6,1}(x)$. On
the other hand, every prime $p$ counted by $\pi_{3,2}(x)$ is of the form
$p=3k+2$, with $k$ odd ($k=2k^{\prime}+1$), so $p=6k^{\prime}+5$; in this
case, 2 is the only prime counted by $\pi_{3,2}(x)$ and not by $\pi_{6,5}$,
and thus $\pi_{3,2}(x)=\pi_{6,5}(x)+1$. These relations,

$\pi_{3,1}(x)=\pi_{6,1}(x)\quad,\quad\pi_{3,2}(x)=\pi_{6,5}(x)+1\,,$ (12)

allow to obtain $\pi_{3,1}(x)$, $\pi_{3,2}(x)$, $\pi_{6,1}(N)$,
$\pi_{6,5}(N)$, $\Delta_{3;\,2,1}(N)$ and $\Delta_{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, $\pi(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 ancillary qubits to
represent binary fractions and hence rational numbers; the QFT acts upon the
two registers, the one representing the integer part and the one representing
the fraction (see Appendix A for a numerical check of the error incurred using
nine ancillary qubits).

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\pi(N)$. According to Eq. (4) and Eq. (9), these
probabilities will then decrease as
$\sim\frac{1}{\phi(\alpha)^{2}}\frac{1}{n\ln\,2}$, where $\alpha=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],
$\\{\epsilon_{i}\\}$, arise from finite sampling, and depend on the number of
measurements performed, $\\#M$, as
$\epsilon_{i}=\sqrt{\frac{p_{i}(1-p_{i})}{\\#M}}$. Hence, in order to achieve
a small multiplicative error $\varepsilon_{i}$, say $\varepsilon_{i}\approx 0$
and $\epsilon_{i}<<p_{i}$, then the number of measurements, or runs of the
quantum circuit, must be $\\#M>>\frac{1-p_{i}}{p_{i}}$. 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(\frac{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 $|\hat{\mathbb{P}}_{n}{\rangle}$ 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{\rangle}$,

$\sum_{\begin{subarray}{c}k=0\\\ k\neq
N/3\end{subarray}}^{N-1}y_{k}\,|e_{k}{\rangle}|0{\rangle}+\sqrt{P(N/3)}\,|N/3{\rangle}|1{\rangle}\,.$
(13)

States like the one in Eq. (13) are ready for Amplitude Estimation.
Classically, the best algorithm for computing modular prime counting
functions, $\pi_{\alpha,\beta}(x)$, is a variation of the Meisser-Lehmer
method for computing $\pi(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, $\Delta(N)$, is unlikely to be of practical use because the
denominator $N\pi(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 $|\hat{\mathbb{P}}_{n}{\rangle}$, 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 $\epsilon$.

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\dots 0{\rangle}$ 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 $|\hat{\mathbb{P}}_{n}{\rangle}$, this
probability reads

$P(0)=\frac{\pi(N)}{N}\,,$ (14)

which decreases as $\sim\frac{1}{n\ln\,2}$ in the limit $n\gg 1$, according to
the PNT, Eq. (4). The prime counting function, $\pi(N)$, appears as well in
the central peak, where $k=N/2$, as shown in Appendix A. The expression for
this peak is

$P(N/2)=\frac{\pi(N)^{2}-4\pi(N)+4}{N\pi(N)}\,,$ (15)

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 $\frac{\sqrt{x}\,\ln\,x}{\pi(x)}$ decreases), direct
sampling from $|\hat{P}_{n}{\rangle}$ 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 $\pi(N)=M/N$, but it is
easier to prepare than the one in Eq. (13). So there is no point in using
$|\hat{\mathbb{P}}_{n}{\rangle}$ for estimating $\pi(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].

##  4 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].

###  4.1 The Prime state

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], $\rho_{A}$, is asymptotically given
by $\overline{\rho}_{A}$ [7], whose expression is

$\overline{\rho}_{A}=\frac{1}{d}\left({\bf\mathbbm{1}}+\ell_{N}{\bf
C}_{m}\right)\,,$ (16)

where $d=2^{m-1}$,
$\ell_{N}\equiv\frac{Li_{2}(N)}{Li(N)}\xrightarrow{N\rightarrow\infty}\frac{1}{n\ln\,2}$,
and ${\bf C}_{m}$ is a $d\times d$ symmetric Toeplitz matrix defined by

$({\bf C}_{m})_{i,j}\equiv(1-\delta_{i,j})\,C(2|i-j|),\quad i,j=1,\dots,d\,,$
(17)

where $C(h)$ are the Hardy-Littlewood constants, defined by

$C(h)=\left\\{\begin{array}[]{lc}2\,C_{2}\,\prod_{p>2;\,p|h}\frac{p-1}{p-2}&\textrm{
if $h$ is even}\\\ \\\ 0&\textrm{otherwise}\end{array}\right.\,,$ (18)

with $C_{2}=\prod_{p>2}\left(1-\frac{1}{(p-1^{2})}\right)$ 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 $\overline{\rho}_{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. $\pi_{k}(x)$, is given by

$\pi_{k}(x)\sim
C(k)\,Li_{2}(x)\;\xrightarrow{\,x\rightarrow\infty\,}\;C(k)\,\frac{x}{(\ln\,x)^{2}}\,,$
(19)

where $Li_{2}(x)=\int_{2}^{x}\frac{dt}{(\ln\,t)^{2}}$.

An analytical approximation to the positive eigenvalues of the matrix ${\bf
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,

$\displaystyle C(h)$ $\displaystyle=$
$\displaystyle\sum_{k=1}^{\infty}\left(\frac{\mu(k)}{\phi(k)}\right)^{2}c_{k}(h)$
(20)

where $\phi(x)$ is the Euler’s totient function, $\mu(x)$ is the Möbius
function (see below), and

$\displaystyle c_{k}(h)$ $\displaystyle=$
$\displaystyle\sum_{\begin{subarray}{c}l=1\\\ {\rm
gcd}(l,k)=1\end{subarray}}^{k}e^{2\pi i\,hl/k}$ (21)

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].

Figure 2: Comparison of the exact lowest 72 entanglement energies of the
entanglement spectrum of the Prime state to those computed using Eq. (23), for
$n=30$ qubits. We expect this approximation to improve for $n\gtrsim 64$,
based on the very good approximation obtained for the slope of the
entanglement entropy, Eq. (29), shown in Fig. 15-right in Appendix C.

Based on Eq. (20), a more precise description of the eigenvalues of ${\bf
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 $\gamma_{k}$ of ${\bf C}_{m}$ with degeneracy
$\phi(k)$. These eigenvalues are given by

$\gamma_{k}\simeq
2^{m}\,\mu^{2}(k)\left(\frac{1}{\phi^{2}(k)}-\frac{1}{\phi_{m}}\right),\quad
k=1,3,5,\dots,k_{m}\,,$ (22)

where the Möbius function, $\mu(x)$, takes values $\pm 1$ for square-free
numbers and 0 otherwise. The values of $k_{m}$ and $\phi_{m}$ are fixed by the
dimension and normalization of the density matrix in Eq. (16), and $k_{m}$
scales as $2^{m/2}$ (see Appendix C for details). There are cases where the
eigenvalue is the same for two different values of $k,k^{\prime}$, in which
case the degeneracy is higher. An example of such an accidental degeneracy is
$\widehat{\gamma}_{13}=\widehat{\gamma}_{21}=\frac{1}{144}$, which has a total
degeneracy $\phi(13)+\phi(21)=24$ (we denote by
$\widehat{\gamma}_{k}=\mu^{2}(k)/\phi^{2}(k)$). Another example is
$\widehat{\gamma}_{35}=\widehat{\gamma}_{39}=\frac{1}{576}$, that has a
degeneracy $\phi(35)+\phi(39)=48$.

Using Eq. (22), the eigenvalues of the density matrix (16) are approximated by

$\lambda_{k}\simeq
2^{1-m}\,\mu^{2}(k)\left(1+\frac{2^{m}}{n\ln\,2}\left(\frac{1}{\phi^{2}(k)}-\frac{1}{\phi_{m}}\right)\right)\,,$
(23)

with a degeneracy $\phi(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 $\\{\varepsilon_{k}\\}$ of a density
matrix $\rho$ with eigenvalues $\\{\lambda_{k}\\}$ is defined as

$\varepsilon_{k}=-\ln\,\lambda_{k}\,.$ (24)

Figure 3: top) Ratio of ${\rm Tr}\left({\bf C}_{m}^{s}\right)$ to Eq. (25),
for $s=2,3,4,5$ and $m=\frac{n}{2}$, up to $n=32$ qubits. bottom) Natural
logarithm of ($1-\,$ratio of ${\rm Tr}\left({\bf C}_{m}^{s}\right)$ to Eq.
(25)), up to $n=32$ qubits.

Using again Eq. (22), we can estimate the contribution of the positive
eigenvalues to the trace of the powers of ${\bf C}_{m}$, from which all Rényi
entropies can be readily derived [43]. The computation is as follows,

$\displaystyle{\rm Tr}\;{\bf C}_{m}^{s}$ $\displaystyle\simeq$
$\displaystyle\sum_{k=1\,({\rm odd})}^{k_{m}}\phi(k)\,\gamma_{k}^{s}$ (25)
$\displaystyle\xrightarrow{\,m\gg 1\,}$ $\displaystyle 2^{ms}\sum_{k=1\,({\rm
odd})}^{\infty}\frac{|\mu(k)|}{(\phi(k))^{2s-1}}$ $\displaystyle=$
$\displaystyle
2^{ms}\prod_{p>2}\left(1+\frac{|\mu(p)|}{(\phi(p))^{2s-1}}\right)$
$\displaystyle=$ $\displaystyle
2^{ms}\prod_{p>2}\left(1+\frac{1}{(p-1)^{2s-1}}\right)\,,$

where we have used that for any function $f(k)$ which is multiplicative, it
holds that

$\sum_{k=1\,({\rm
odd})}^{\infty}|\mu(k)|\,f(k)=\prod_{p>2}\left(1+f(p)\right)\,,$ (26)

where the product runs over the odd prime numbers. We have also used that
$\phi(p)=p-1$ for prime numbers, and that $\mu(k)^{2}=|\mu(k)|$. In Eq. (25)
we dropped the negative term proportional to $1/\phi_{m}$ of Eq. (22) because
its contribution vanishes in the limit $N,m\gg 1$ for $s>1$. For $s=2$, Eq.
(25) coincides with the heuristically-derived asymptotic expansion of ${\rm
Tr}\;{\bf C}_{m}^{2}$, and was conjectured to hold for any $s\geq 2$ in Ref.
[7]. A direct connection thus exists between Eq. (22) and the latter
conjecture. For $s=1$, the Riemann zeta function $\zeta(1)$ arises and the
product diverges because we have dropped the negative eigenvalues of the
matrix ${\bf C}_{m}$. Figure 3 shows the ratio of the exact value of the trace
of the matrix ${\bf C}_{m}^{s}$ 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
${\bf C}_{m}$ to the trace of the powers of this matrix is asymptotically
negligible, for $s\geq 2$. In contrast, their contribution to the von Neumann
entropy does not appear to be so.

Figure 4: top) Value of $c_{\pi}(n)$ computed according to Eq. (29), compared
to the conjectured value $H\left(\frac{3}{\pi^{2}}\right)$, up to $n=30$
qubits; bottom) Value of $\gamma(n)$ computed according to Eq. (30), up to
$n=30$ qubits, compared to the conjectured value $1+\frac{3}{\pi^{2}}$.

The von Neumann entropy $S$ of a density matrix $\rho$, with eigenvalues
$\\{\lambda_{i}\\}$, is defined by

$S=-\textrm{Tr}\,(\rho\log_{2}\,\rho)=-\sum_{i}\lambda_{i}\log_{2}\lambda_{i}\,,$
(27)

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\,\frac{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 $\frac{n}{2}$ with a slope
equal to 1. Random states, i.e. those with random complex coefficients, follow
a “volume law”, with an entanglement entropy of half chain that scales as
$n/2-1/2$ for big $n$, where $n$ is the number of qubits [44]. If we restrict
to real positive random coefficients, the scaling is
$(1-2/\pi)\,n/2=0.363\dots n/2$ [45]. So the Prime state is not a typical
random positive state either. Notice as well that the high entanglement of the
Prime state implies that it is not possible to apply standard classical
techniques such as Matrix Product States (MPS), or other more-refined tensor
networks, to efficiently simulate the Prime state on a classical computer
[46].

Figure 5: Natural logarithm of the absolute value of the difference $\delta$,
between the observed and predicted values of $c_{\pi}$ and $\gamma$, up to
$n=30$ qubits. That is,
$\delta\equiv\left|c_{\pi}(n)-H\left(\frac{3}{\pi^{2}}\right)\right|$ and
$\delta\equiv\left|\gamma(n)-\left(1+\frac{3}{\pi^{2}}\right)\right|$,
respectively. An exponential decrease in $\delta$ is observed before reaching
a plateau.

We herein conjecture that the entropy of the natural equal-sized bi-partition
of the Prime state of $n$ qubits is asymptotically given by

$S(n)=c_{\pi}\,\frac{n}{2}-\gamma\quad,\quad
c_{\pi}=H\left(\frac{3}{\pi^{2}}\right)\quad,\quad\gamma=1+\frac{3}{\pi^{2}}\,,$
(28)

where $H(p)\equiv-p\log_{2}(p)-(1-p)\log_{2}(1-p)$ is the Shannon entropy and
$3/\pi^{2}=1/(2\,\zeta(2))$ is equal to a half of the asymptotic density of
odd square-free integers [16]. The constant $3/\pi^{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 [47, 48]. The
term -1 in the intercept comes from the fact that the least relevant qubit is
basically in the state $|1{\rangle}$, 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 of $n$, by the observed oscillatory behaviour of the slope
$c_{\pi}$. To the best of our knowledge, this is the first open-source 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].

Figure 6: top) Von Neumann entropy $S$ of the Prime state of $n$ qubits (up to
$n=29$) and its Quantum Fourier Transform (up to $n=22$), for the natural bi-
partition of size $\lfloor\frac{n}{2}\rfloor$. Notice that the entropies of
the Prime state and its QFT are almost identical; bottom) Von Neumann entropy
$S$ for the natural equal-sized bi-partition of Prime states expressed in
$q$-nary bases using $n$ qudits.

Figure 4-top shows the difference in entropy between two consecutive values of
$S(n)$, for even $n$. That is, the dots in the plot correspond to

$c_{\pi}(n)=S(n)-S(n-2)\,,$ (29)

up to $n=30$ qubits. Oscillatory convergence towards the predicted value
$H\left(\frac{3}{\pi^{2}}\right)=0.886082085\dots$ is observed, in good
agreement with Eq. (28). The value of the intercept $\gamma$, obtained as

$\gamma(n)=\left(S(n)-S(n-2)\right)\frac{n}{2}-S(n)\,,$ (30)

is shown in Fig. 4-bottom, up to $n=30$ qubits. Oscillatory convergence and
good agreement with the predicted value $1+\frac{3}{\pi^{2}}=1.30396355\dots$
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_{\pi}$ and $\gamma$ converge may be
seen as an indication that both conjectures are closely related.

Figure 7: top) Von Neumann entropy $S$ of Prime states with different number
of qubits, up to $n=28$, for the natural bi-partition that separates the $m$
first (i.e. least significant) qubits; bottom) Same von Neumann entropy $S$,
rescaled by $S(n/2)$.

We have also extended the computations to the natural equal-sized bi-partition
of size $\lfloor\frac{n}{2}\rfloor$, with odd $n$ up to $n=29$ qubits, see
Fig. 6-top. The scaling of the entropy is very similar in this case to that of
the Prime state with an even number of qubits, as it could be expected.
Besides, we have computed the von Neumann entropy of Prime states for qudit
systems, with local dimensions $d=3,4,5$, for natural equal-sized bi-
partitions, see Fig. 6-bottom and Table 1. 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 $\sim 0.88\dots$, 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 1.

Furthermore, we have computed the entropy for natural bi-partitions of size
$m$ of the Prime state, for different number of qubits $n$, up to $n=28$
qubits, see Fig. 7-top. Assuming a scaling behaviour of the form
$S(m,n)=n\,f(m/n)$, it follows that

$\frac{S(m,n)}{S(n/2,n)}=\frac{f(m/n)}{f(1/2)}\,.$ (31)

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$.

$State$ | $Slope$ | $Intercept$
---|---|---
Prime (qubits) | 0.885791 | -1.299313
Prime (qutrits) | 0.887102 | -0.733640
Prime (quatrits) | 0.885822 | -0.649876
Prime (5-qudits) | 0.886266 | -0.421618
QFT Prime | 0.882952 | -0.382451
Arithmetic Prime (mod 4) | 0.884033 | -1.885681
Arithmetic Prime (mod 8) | 0.883698 | -2.575452
Arithmetic Prime (mod 16) | 0.884932 | -3.339827
Arithmetic Prime (mod 32) | 0.888128 | -4.169144
Möbius | 0.999843 | -0.953885
Starry Prime | 0.985580 | -0.948598
Table 1: Characterization of the linear scaling of the von Neumann entropy,
for several number-theoretical quantum states. The entropy corresponds to the
equal-sized natural bi-partition, and the slopes and intercepts have been
calculated using Eq. (29) and Eq. (30), for $n=30$; except for the Möbius
state ($n=28$), the QFT of the Prime state ($n=22$), and for qudit systems
($n=18,14,12$ for $d=3,4,5$ respectively). In the case of arithmetic Prime
states, averages over the different arithmetic progressions are shown.

###  4.2 Arithmetic Prime states

Figure 8: Von Neumann entropy $S$ for the natural equal-sized bi-partition of
Prime states in arithmetic sequences, compared to that of the Prime state, up
to $n=30$ qubits. Observe the collapse of data for the sequences with the same
moding, when $n$ is large.

Other number-theoretical quantum states apart from the Prime state will be
defined below, and their entanglement traits studied. Table 1 summarizes our
findings. We define arithmetic Prime states as the uniform superpositions of
primes belonging to arithmetic progressions, i.e.

$|\mathbb{P}_{n;\,\alpha,\beta}{\rangle}\equiv\frac{1}{\sqrt{\pi_{\alpha,\beta}(2^{n})}}\sum_{\begin{subarray}{c}p:\,prime\\\
p=\beta\,\textrm{mod}\,\alpha\end{subarray}}^{2^{n}}|p{\rangle}\,,$ (32)

with $\alpha,\beta$ coprime numbers. When $\alpha$ 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,

$\overline{\rho}_{A}(\alpha)=\frac{\phi(\alpha)}{d}\left(\mathbbm{1}+l_{N}{\bf
C}_{m}(\alpha)\right)\,,$ (33)

where again $d=2^{m-1}$,
$l_{N}\equiv\frac{Li_{2}(N)}{Li(N)}\xrightarrow{N\rightarrow\infty}\frac{1}{n\ln\,2}$,
and ${\bf C}_{m}(\alpha)$ is a symmetric Toeplitz matrix defined by

$({\bf C}_{m}(\alpha))_{i,j}\equiv(1-\delta_{i,j})\,C(\alpha\,|i-j|)\,.$ (34)

Notice that Eq. (33) generalizes the asymptotic expression of the reduced
density matrix of the Prime state. Indeed, all primes but 2, which does not
contribute in the asymptotic limit, belong to the arithmetic progression
$2k+1$. Thus, $\phi(2)=1$, ${\bf C}_{m}(2)={\bf C}_{m}$, and we recover Eq.
(16). The validity of Eq. (33) and Eq. (16) is conditional to the first Hardy-
Littlewood conjecture and to the following asymptotic behaviour for the prime-
number function $\pi_{\alpha;\,\beta,\beta^{\prime}}(x)$,

$\pi_{\alpha;\,\beta,\beta^{\prime}}(x)\,\xrightarrow{\;x\rightarrow\infty\;}\,\frac{C(|\beta-\beta^{\prime}|)}{\phi(\alpha)}\,Li_{2}(x)\,,$
(35)

where $\pi_{\alpha;\,\beta,\beta^{\prime}}$ is the number of prime pairs
$(p,p^{\prime})$ such that $p,p^{\prime}$ are primes smaller than, or equal
to, $x$ of the form $p=\alpha k+\beta$, $p=\alpha k+\beta^{\prime}$, with
$(\alpha,\beta)$ and $(\alpha,\beta^{\prime})$ coprimes pairs. This asymptotic
behaviour, Eq. (35), was numerically observed in Ref. [7].

We have found that the scaling of the entropy of arithmetic Prime states,
$|\mathbb{P}_{n;\,\alpha,\beta}{\rangle}$, when $\alpha$ is a power of 2, is
once again dictated by a straight line with slope $\sim 0.884$, see Fig. 8 and
Table 1. The intercept, however, is different for distinct values of $\alpha$.
It is also independent of $\beta$. The latter fact is easily understood from
the asymptotic expression of the reduced density matrix
$\overline{\rho}_{A}(\alpha)$, which does not depend on $\beta$. This is
ultimately so because prime numbers are equally distributed on average among
the different arithmetic progressions, for a given $\alpha$, as expressed in
Eq. (9). Hence, we label this constant as $\gamma_{\alpha}$. We have found the
following behaviour for $\gamma_{\alpha}$ when $\alpha$ is a power of 2,

$\gamma_{\alpha}\simeq 1+s\frac{3}{\pi^{2}}\quad,\quad\alpha=2^{k}\quad,\quad
s=1,3,5,\dots\,,$ (36)

where $s$ is the $k$-th odd square-free number. Figure 9 shows the
independence of $\gamma_{4}$ and $\gamma_{8}$ from $\beta$, for a large number
of qubits, and the convergence towards a value close to that in Eq. (36).
Using the results obtained for $\gamma_{4}$, $\gamma_{8}$, $\gamma_{16}$ and
$\gamma_{32}$ for $n=30$ qubits, averaged among the different progressions,
the value $s=\left(\gamma_{\alpha}-1\right)\frac{\pi^{2}}{3}$ is found to be
$s=2.9137\dots$, $s=5.1830\dots$, $s=7.6977\dots$ and $s=10.4260\dots$, for
$\alpha=4,8,16,32$ respectively.

Taking into account all the results presented in this section, we generalize
the conjecture in Eq. (28), stating that, for Prime and arithmetic Prime
states of $n$ qudits, the von Neumann entropy $S$ of equal-sized bi-partitions
is asymptotically given by

$S(n)=c_{\pi}\,\lfloor\frac{n}{2}\rfloor-\gamma_{\alpha;\,d}\quad,\quad
c_{\pi}=H\left(\frac{3}{\pi^{2}}\right)\,,$ (37)

where $\lfloor x\rfloor$ is the floor function and $\gamma_{\alpha;\,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.

Figure 9:  top) Value of $\gamma_{4}(n)$ computed according to Eq. (30), up to
$n=30$ qubits, compared to the value $1+3\times\frac{3}{\pi^{2}}$ in Eq. (36);
bottom) Value of $\gamma_{8}(n)$ computed according to Eq. (30), up to $n=30$
qubits, compared to the value $1+5\times\frac{3}{\pi^{2}}$ in Eq. (36).

###  4.3 Quantum relatives of the Prime state

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

$|\not{\mathbb{P}}_{n}{\rangle}\equiv\frac{1}{\sqrt{C_{c,n}}}\sum_{\begin{subarray}{c}c:\,odd\\\
c:\,composite\end{subarray}}^{2^{n}-1}|c{\rangle}\quad,\quad
C_{c,n}=2^{n-1}-\pi(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 $\sim 2$,
irrespective of the number of qubits $n$, see Fig. 10. This behaviour is
reminiscent of finitely-correlated systems away from criticality [49]. 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\rightarrow\infty$, by

$\overline{\rho}_{A}=\frac{1}{d}\left(\mathbbm{1}+{\bf P}_{m}\right)\,,$ (39)

where

$\left({\bf P}_{m}\right)_{ij}=\left\\{\begin{array}[]{lc}0&{\rm if}\;i=j\\\
1&{\rm if}\;i\neq j\end{array}\right.\,.$ (40)

From Eq. (39), the von Neumann entropy $S$ in the asymptotic limit is given by

$S\simeq 0\,,$ (41)

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
$|\phi{\rangle}=\frac{1}{\sqrt{2^{n}}}\,\sum_{i=0}^{2^{n}-1}|e_{i}{\rangle}$.
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

$|\mathbb{S}_{n}{\rangle}\equiv\frac{1}{\sqrt{C_{s,n}}}\,\sum_{s:\,odd}^{2^{n}}\left|\mu(s)\right|\,|s{\rangle}\quad,\quad
C_{s,n}=\sum_{s:\,odd}^{2^{n}}|\mu(s)|\,.$ (42)

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\sim
20$ qubits.

Figure 10: Von Neumann entropy $S$ of the natural equal-sized bi-partition of
the odd composite, odd square-free and Möbius states, up to $n=28$ 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

$|\mu_{n}{\rangle}\equiv\frac{1}{\sqrt{C_{\mu,n}}}\,\sum_{s=1}^{2^{n}}\mu(s)\,|s{\rangle}\quad,\quad
C_{\mu,n}=\sum_{s=1}^{2^{n}}\left|\mu(s)\right|\,.$ (43)

The von Neumann entropy of the natural equal-sized bi-partition 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 1). 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

$p_{n}\leq p^{*}_{n}\leq p_{n+1}\,,$ (44)

where $p_{n}$ is the $n$-th prime. These numbers are remarkable in that a
$\zeta^{*}$ function, defined via the Euler product formula, i.e.

$\zeta^{*}(s)\equiv\prod_{p^{*}}\frac{1}{1-{p^{*}}^{-s}}\,,$ (45)

has an analytic extension with exactly the same zeros, and the same
multiplicities, than the original Riemann zeta function [50]. The starry-Prime
state is then defined by

$|\mathbb{P^{*}}_{n}{\rangle}\equiv\frac{1}{\sqrt{\pi(2^{n})}}\,\sum_{\begin{subarray}{c}p^{*}:\,starry\\\
\quad prime\end{subarray}}^{2^{n}-1}|p^{*}{\rangle}\,,$ (46)

where $\\{p^{*}\\}$ is any set of starry primes up to $2^{n}-1$. We have found
that the entropy of the natural equal-sized bi-partition 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 $\frac{n}{2}$ with slope $\sim 1$ (see Fig. 11 and Table 1). 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.

Figure 11: Von Neumann entropy $S$ of random starry-Prime states for the
natural equal-sized bi-partition, compared to that of the Prime state, up to
$n=30$ qubits. Each red dot is an average over 10 random starry-Prime states.

##  5 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.

Acknowledgements. We thank J. Andrade, A. Botero, J. Keating, A. LeClair, G.
Mussardo, D. Pérez-García and S. Ramos-Calderer for useful discussions. DGM
and JIL acknowledge CaixaBank for its support of this work through Barcelona
Supercomputing Center’s project CaixaBank Computación Cuántica. DGM and JIL
are supported by Project PGC2018-095862-B-C22, and Project Quantum CAT (001-
P-001644) co-funded by Generalitat de Catalunya and the European Union
Regional Development Fund. SC is supported by the European Research Council
under the European Unions Horizon 2020 research and innovation Programme
(grant agreement number 740006). GS is supported by the Ministerio de Ciencia,
Innovación y Universidades (grant PGC2018-095862-B-C21), the Comunidad de
Madrid (grant QUITEMAD+ S2013/ICE-2801), the Centro de Excelencia Severo Ochoa
Programme (grant SEV-2016-0597) and the CSIC Research Platform on Quantum
Technologies PTI-001.

## References

  * [1] J.I. Latorre and G. Sierra, “Quantum computation of prime number functions”, Quant. Inf. Comput. 14, 577 (2014).
  * [2] L.K. Grover, “A fast quantum mechanical algorithm for database search”, Proc. STOC. May 1996, 212 (1996).
  * [3] B. Riemann, “On the Number of Prime Numbers less than a Given Quantity”, Monatsberichte der Berliner Akademie November 1859, 671 (1859), English version.
  * [4] K. Walisch and D. Baugh, “New confirmed $\pi(10^{27}$) prime counting function record”, Mersenne Forum (2015), https://github.com/kimwalisch/primecount.
  * [5] G. Brassard, P. Høyer and A. Tapp, “Quantum Counting”, Proc. 25th ICALP, LNCS 1443, Springer-Verlag, 820 (1998).
  * [6] M. Rubinstein and P. Sarnak, “Chebyshev’s Bias”, Exp. Math. 3, 173 (1994).
  * [7] J.I. Latorre and G. Sierra, “There is entanglement in the primes”, Quant. Inf. Comput. 15, 622 (2015).
  * [8] G.H. Hardy and J.E. Littlewood, “Some Problems of ’Partitio Numerorum.’ III. On the Expression of a Number as a Sum of Primes”, Acta Math. 44, 1 (1923).
  * [9] We consider natural bi-partitions those that separate the first $m$ qubits and the remainder $n-m$ qubits, without any reordering.
  * [10] G. Mussardo, H. Giudici and J. Viti, “The Coprime Quantum Chain”, J. Stat. Mech. Theory Exp. 2017, 033104 (2017).
  * [11] G. Mussardo, A. Trombettoni and Zhaon Zhang, “Prime suspects in a Quantum Ladder”, arXiv:2005.01758 (2020).
  * [12] F.V. Mendes and R.V. Ramos, “Quantum Sequence States”, arXiv:1408.4838 (2014).
  * [13] S. Carrazza, D. García-Martín and J.I. Latorre, Quantum-TII/qprime: Qprime v1.0.0 (May 2020). https://doi.org/10.5281/zenodo.3787043.
  * [14] L. Miller. “Riemann’s hypothesis and tests for primality”, J. Comput. Sys. Sci. 13, 300 (1976); M. O. Rabin. “Probabilistic algorithm for testing primality”, J. Number Theory 12, 128 (1980).
  * [15] M. Agrawal, N. Kayal and N. Saxena, “Primes is in P”, Ann. Math. 160, 781 (2004).
  * [16] T.M. Apostol, Introduction to Analytic Number Theory, Springer-Verlag, New York (1976); G.H. Hardy and E.M. Wright, An Introduction to the Theory of Numbers, Oxford Universiy Press, Oxford (1938).
  * [17] M.A. Nielsen and I.L. Chuang, Quantum Computation and Quantum Information: 10th Anniversary Edition, Cambridge University Press, New York (2011).
  * [18] https://www.claymath.org/millennium-problems.
  * [19] L. Schoenfeld, “Sharper bounds for the Chebyshev functions $\theta(x)$ and $\psi(x)$. II”, Math. Comput. 30, 337 (1976).
  * [20] A. Kitaev, “Quantum measurements and the Abelian stabilizer problem”, arXiv:quant-ph/9511026 (1995); R. Cleve, A. Ekert, C. Macchiavello and M.Mosca, “Quantum algorithms revisited”, Proc. R. Soc. London A 454, 339 (1998).
  * [21] S. Aaronson, “Quantum Lower Bound for Approximate Counting Via Laurent Polynomials”, arXiv:1808.02420 (2018); S. Aaronson and P. Rall, “Quantum Approximate Counting, Simplified”, SOSA20, 24 (2020).
  * [22] Meissel, “Über die Bestimmung der Primzahlenmenge innerhalb gegebener Grenzen”, Math. Ann. 2, 636 (1870).
  * [23] D.H. Lehmer, “On the exact number of primes less than a given limit”, Illinois J. Math. 3, 381 (1959).
  * [24] J.C. Lagarias, V.S. Miller and A.M. Odlyzko, “Computing $\pi(x)$: the Meissel-Lehmer method”, Math. Comp. 44, 573 (1985).
  * [25] M. Deléglise and J. Rivat, “Computing $\pi(x)$: the Meissel, Lehmer, Lagarias, Miller, Odlyzko method”, Math. Comp. 65, 235 (1996).
  * [26] X. Gourdon, “Computation of $\pi(x)$: improvements to the Meissel, Lehmer, Lagarias, Miller, Odllyzko, Deléglise and Rivat method”, unpublished (2001). As cited in Ref. [4].
  * [27] J.C. Lagarias and A.M. Odlyzko, “Computing $\pi(x)$: an Analytic Method”, J. Algorithms 8, 173 (1987).
  * [28] P. Shor, “Polynomial-time Algorithms for Prime Factorization and Discrete Logarithms on a Quantum Computer”, SIAM J. Sci. Statist. Comput. 26, 1484 (1997).
  * [29] A. Peruzzo, J. McClean, P. Shadbolt, M.H. Yung, X.Q. Zhou, P.J. Love, A. Aspuru-Guzik and J. L. O’Brien, “A variational eigenvalue solver on a quantum processor”, Nat. Commun. 5, 4213 (2013).
  * [30] E. Farhi, J. Goldstone and S. Guttman, “A Quantum Approximate Optimization Algorithm”, arXiv:1411.4028 (2014).
  * [31] If $p_{i}$ is the exact probability, and $\hat{p}_{i}$ the estimate of it, the statistical additive error is defined as $p_{i}=\hat{p}_{i}+\epsilon_{i}$, and the multiplicative error as $\hat{p}_{i}=(1+\varepsilon_{i})\,p_{i}$.
  * [32] G. Brassard, P. Høyer, M. Mosca and A. Tapp, “Quantum Amplitude Amplification and Estimation”, Contemp. Math. 305, 53 (2002).
  * [33] A. Montanaro, “Quantum speedup of Monte Carlo methods”, Proc. R. Soc. A 471, 20150301 (2015).
  * [34] M. Deléglise, P. Dusart and X.F. Roblot, “Counting primes in residue classes”, Math. Comp. 73, 1565 (2004).
  * [35] http://mathworld.wolfram.com/
ModularPrimeCountingFunction.html.

  * [36] I. Bengtsson and K. Zyczkowski, Geometry of quantum states: an introduction to quantum entanglement, Cambridge University Press, New York (2018).
  * [37] We thank Alonso Botero for noticing this fact.
  * [38] J.P. Keating and D.J. Smith, “Twin prime correlations from the pair correlation of Riemann zeros”, J. Phys. A: Math. Theor. 52, 365201 (2019).
  * [39] E. B. Bogomolny and J. P. Keating, “Gutzwiller’s trace formula and spectral statistics: beyond the diagonal approximation”, Phys. Rev. Lett. 77, 1472 (1996).
  * [40] M. V. Berry and J. P. Keating, “The Riemann zeros and eigenvalue asymptotics”, SIAM Rev. 41, 236 (1999).
  * [41] H.G. Gadiyar and R. Padma, “Ramanujan-Fourier series, the Wiener-Khintchine formula and the distribution of prime pairs”, Phys. A 269, 503 (1999).
  * [42] H. G. Gadiyar and R. Padma, “Linking the circle and the sieve: Ramanujan-Fourier series”, arXiv:math/0601574 (2006).
  * [43] A. Rényi, “On measures of information and entropy”, Proceedings of the fourth Berkeley Symposium on Mathematics, Statistics and Probability 1960, 547 (1961).
  * [44] D.N. Page, “Average entropy of a subsystem”, Phys. Rev. Lett. 71, 129 (1993).
  * [45] T. Grover and M.P.A. Fisher, “Entanglement and the Sign Structure of Quantum States”, Phys. Rev. A 92, 042308 (2015).
  * [46] J.I. Latorre, “Entanglement entropy and the simulation of quantum mechanics”, J. Phys. A: Math. Theor. 40, 6689 (2007).
  * [47] P. Sarnak, “Three lectures on the Möbius function, randomness, and dynamics”, online notes (2011), http://publications.ias.edu/sarnak/paper/512.
  * [48] F. Cellarosi and Y. G. Sinai, “Ergodic Properties of Square-Free Numbers”, arXiv:1112.4691 (2011).
  * [49] L. Amico, R. Fazio, A. Osterloh and V. Vedral, “Entanglement in Many-Body Systems”, Rev. Mod. Phys. 80, 517 (2008).
  * [50] E. Grosswald and F. J. Schnitzer, “A class of modified $\zeta$ and L-functions”, Pacific. Jour. Math. 74, 357 (1978).
  * [51] D. R. Ward, “Some Series Involving Euler’s Function”, J. London Math. Soc. 2, 210 (1927).
  * [52] J. H. van Lint and H.E. Richert, “Über die Summe $\sum_{n}\mu^{2}(n)/\varphi(n)$”, Nederl. Akad. Wetensch.‘ Proc. Ser. A 67, Indag. Math. 26, 582 (1964).

##  Appendices

##  Appendix A  Computation of relevant peaks in
$|\hat{\mathbb{P}}_{n}{\rangle}$

We compute here the most relevant peaks of the QFT of the Prime state,
$|\hat{\mathbb{P}}_{n}{\rangle}$. We use Eq. (6), where
$x_{j}=\frac{1}{\sqrt{\pi(N)}}$ for prime $j$ and 0 otherwise. Recall that
$N=2^{n}$ is the dimension of the Hilbert space.

$\begin{split}P(N/2)=\quad&\frac{1}{N\pi(N)}\left|\sum_{p:\,prime}e^{p\pi
i}\right|^{2}=\;\frac{1}{N\pi(N)}\left|\sum_{\begin{subarray}{c}p:\,prime\\\
p=0\,\textrm{mod2}\end{subarray}}e^{p\pi
i}+\sum_{\begin{subarray}{c}p:\,prime\\\
p=1\,\textrm{mod2}\end{subarray}}e^{p\pi i}\right|^{2}=\\\ \\\
&\frac{1}{N\pi(N)}\left|\pi_{2,0}(N)-\pi_{2,1}(N)\right|^{2}=\;\frac{1}{N\pi(N)}\left|1-(\pi(N)-1)\right|^{2}=\quad\frac{\pi(N)^{2}-4\pi(N)+4}{N\pi(N)}\end{split}$
(47)
$\begin{split}P(N/3)=\quad&\frac{1}{N\pi(N)}\left|\sum_{p:\,prime}e^{\frac{2}{3}\pi
i\,p}\right|^{2}=\;\frac{1}{N\pi(N)}\left|\sum_{\begin{subarray}{c}p:\,prime\\\
p=0\,\textrm{mod3}\end{subarray}}e^{\frac{2}{3}p\pi
i}+\sum_{\begin{subarray}{c}p:\,prime\\\
p=1\,\textrm{mod3}\end{subarray}}e^{\frac{2}{3}p\pi
i}+\sum_{\begin{subarray}{c}p:\,prime\\\
p=2\,\textrm{mod3}\end{subarray}}e^{\frac{2}{3}p\pi i}\right|^{2}=\\\ \\\
&\frac{1}{N\pi(N)}\left|\pi_{3,0}(N)+\pi_{3,1}(N)\,e^{\frac{2}{3}\pi
i}+\pi_{3,2}(N)\,e^{\frac{4}{3}\pi i}\right|^{2}=\\\ \\\
&\frac{1}{N\pi(N)}\left|1+\left(-\frac{1}{2}+i\frac{\sqrt{3}}{2}\right)\pi_{3,1}(N)+\left(-\frac{1}{2}-i\frac{\sqrt{3}}{2}\right)\pi_{3,2}(N)\right|^{2}=\\\
\\\
&\frac{1}{N\pi(N)}\left|1-\frac{1}{2}\left(\pi_{3,1}(N)+\pi_{3,2}(N)\right)-i\frac{\sqrt{3}}{2}\,\Delta_{3;2,1}(N)\right|^{2}=\\\
\\\
&\frac{1}{N\pi(N)}\left(1+\frac{1}{4}\left(\pi_{3,1}(N)+\pi_{3,2}(N)\right)^{2}-\pi_{3,1}(N)-\pi_{3,2}(N)+\frac{3}{4}\Delta_{3;2,1}(N)^{2}\right)=\\\
\\\
&\frac{\Delta_{3;2,1}(N)^{2}+\pi_{3,1}(N)\pi_{3,2}(N)-\pi(N)+2}{N\pi(N)}\end{split}$
(48)
$\begin{split}P(N/4)=\quad&\frac{1}{N\pi(N)}\left|\sum_{p:\,prime}e^{\frac{1}{2}\pi
i\,p}\right|^{2}=\;\frac{1}{N\pi(N)}\left|\sum_{\begin{subarray}{c}p:\,prime\\\
p=2\,\textrm{mod4}\end{subarray}}e^{\frac{1}{2}p\pi
i}+\sum_{\begin{subarray}{c}p:\,prime\\\
p=1\,\textrm{mod4}\end{subarray}}e^{\frac{2}{4}p\pi
i}+\sum_{\begin{subarray}{c}p:\,prime\\\
p=3\,\textrm{mod4}\end{subarray}}e^{\frac{2}{4}p\pi i}\right|^{2}=\\\ \\\
&\frac{1}{N\pi(N)}\left|-\pi_{4,0}(N)+\pi_{4,1}(N)\,e^{\frac{1}{2}\pi
i}+\pi_{4,3}(N)\,e^{\frac{3}{2}\pi
i}\right|^{2}=\;\frac{1}{N\pi(N)}\left|-1+i\left(\pi_{4,1}(N)-\pi_{4,3}(N)\right)\right|^{2}=\\\
\\\ &\frac{1+\Delta(N)^{2}}{N\pi(N)}\end{split}$ (49)
$\begin{split}P(N/6)=&\quad\frac{1}{N\pi(N)}\left|\sum_{p:\,prime}e^{\frac{2}{6}\pi
i\,p}\right|^{2}=\;\frac{1}{N\pi(N)}\left|\sum_{\begin{subarray}{c}p:\,prime\\\
p=2\,\textrm{mod6}\end{subarray}}e^{\frac{1}{3}p\pi
i}+\sum_{\begin{subarray}{c}p:\,prime\\\
p=3\,\textrm{mod6}\end{subarray}}e^{\frac{1}{3}p\pi
i}+\sum_{\begin{subarray}{c}p:\,prime\\\
p=1\,\textrm{mod6}\end{subarray}}e^{\frac{2}{6}p\pi
i}+\sum_{\begin{subarray}{c}p:\,prime\\\
p=5\,\textrm{mod6}\end{subarray}}e^{\frac{2}{6}p\pi i}\right|^{2}=\\\ \\\
&\frac{1}{N\pi(N)}\left|\pi_{6,0}(N)\,e^{\frac{2}{3}\pi i}+\pi_{6,3}(N)e^{\pi
i}+\pi_{6,1}(N)\,e^{\frac{1}{3}\pi i}+\pi_{6,5}(N)\,e^{\frac{5}{3}\pi
i}\right|^{2}=\\\ \\\
&\frac{1}{N\pi(N)}\left|\left(-\frac{1}{2}+i\frac{\sqrt{3}}{2}\right)-1+\left(\frac{1}{2}+i\frac{\sqrt{3}}{2}\right)\pi_{6,1}(N)+\left(\frac{1}{2}-i\frac{\sqrt{3}}{2}\right)\pi_{6,5}(N)\right|^{2}=\\\
\\\
&\frac{1}{N\pi(N)}\left|\frac{\pi_{6,1}(N)+\pi_{6,5}(N)-3}{2}-i\frac{\sqrt{3}\,\left(\Delta_{6;5,1}(N)-1\right)}{2}\right|^{2}=\\\
\\\
&\frac{1}{N\pi(N)}\left(\frac{9+\pi_{6,1}(N)^{2}+\pi_{6,5}(N)^{2}+2\,\pi_{6,1}(N)\pi_{6,5}(N)-6\,\pi_{6,1}(N)-6\,\pi_{6,5}(N)}{4}+\frac{3\,\Delta_{6;5,1}(N)^{2}+3-6\,\Delta_{6;5,1}(N)}{4}\right)=\\\
\\\
&\frac{\Delta_{6;5,1}(N)^{2}+\pi_{6,1}(N)\pi_{6,5}(N)-3\,\pi_{6,5}(N)+3}{N\pi(N)}\end{split}$
(50)

Given that $N/3$ and $N/6$ are not integers when $N=2^{n}$, ancillary qubits
are required to represent rational numbers as binary fractions. A numerical
check of the error incurred in the estimation of the peak when using nine
ancillary qubits to represent the non-integer part of $N/3$ or $N/6$ is shown
in Fig. 12, for $n=10$-30 qubits. The relative error is around $10^{-6}$ for
$n=30$ qubits. If more precision is required, one can simply add more
ancillas.

Figure 12: Relative percentage error in the estimation of Eqs. (48), (50) due
to the limited precision in the representation of $N/3$ and $N/6$, using nine
ancillary qubits, for $n=10-30$ qubits.

##  Appendix B Quantum Fourier Transform of a uniform superposition

Let us consider the definition of the QFT of a general quantum state
$|\psi{\rangle}$ of $n$ qubits, expressed in the computational basis,
$\\{|e_{j}\\}$,

$|\psi{\rangle}=\sum_{j=0}^{2^{n}-1}x_{j}\,|e_{j}{\rangle}\xrightarrow{\quad
QFT\quad}\sum_{k=0}^{2^{n}-1}y_{k}\,|e_{k}{\rangle}\,,$ (51)

where

$y_{k}\equiv\frac{1}{\sqrt{2^{n}}}\,\sum_{j=0}^{2^{n}-1}x_{j}\,e^{2\pi
i\,jk/2^{n}}\,,$ (52)

and the $i$ in the exponent is the imaginary unit. The probability of
measuring the $|0\dots 0{\rangle}$ state in Eq. (51) in the computational
basis thus is:

$P(0)=|y_{0}|^{\,2}=\left|\frac{1}{\sqrt{2^{n}}}\,\sum_{j=0}^{2^{n}-1}x_{j}\,e^{0}\,\right|^{2}=\frac{1}{\,2^{n}\,}\,\left|\sum_{j=0}^{2^{n}-1}x_{j}\,\right|^{2}\,.$
(53)

If the state $|\psi{\rangle}$ is a uniform superposition of $M$ computational-
basis vectors, i.e. $x_{j}=1/\sqrt{M}$ for all $j$ such that $x_{j}\neq 0$,
and $N=2^{n}$ is the dimension of the Hilbert space, then

$P(0)=\frac{1}{\,N\,}\,\left|\sum_{j=0}^{M-1}\frac{1}{\sqrt{M}}\,\right|^{2}=\frac{1}{\,N\,}\,\left|\frac{M}{\sqrt{M}}\,\right|^{2}=\frac{M}{\,N\,}\;.$
(54)

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.

##  Appendix 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 by

$\bar{\rho}_{A}=\frac{1}{d}({\bf 1}+\ell_{N}{\bf
C}_{m}),\qquad\ell_{N}\xrightarrow{N\rightarrow\infty}\frac{1}{n\ln 2}\,,$
(55)

where ${\bf C}_{m}$ is a matrix constructed with the even values of the Hardy-
Littlewood constants $C(2k)$,

$({\bf C}_{m})_{i,j}=(1-\delta_{i,j})\,C(2|i-j|),\quad
i,j=1,2,\dots,d=2^{m-1}.$ (56)

In order to justify the conjecture (22) for the eigenvalues of ${\bf C}_{m}$,
we define the density matrix

${\rho}_{D}=\frac{1}{D}({\bf 1}+\ell_{D}{\bf
C}_{D}),\qquad\ell_{D}=\frac{1}{2\ln(2D)}\,,$ (57)

with

$({\bf C}_{D})_{i,j}=(1-\delta_{i,j})\,C(2|i-j|),\quad i,j=1,2,\dots,D.$ (58)

If $D=2^{m-1}$ and $n=2m$, Eq. (57) and Eq. (58) become Eq. (55) and Eq. (56)
respectively. The Hardy-Littlewood constants can be expressed as [41, 42, 38]

$\displaystyle C(h)$ $\displaystyle=$
$\displaystyle\sum_{k=1}^{\infty}\left(\frac{\mu(k)}{\phi(k)}\right)^{2}c_{k}(h)\,,$
(59)

where

$\displaystyle c_{k}(h)$ $\displaystyle=$
$\displaystyle\sum_{\begin{subarray}{c}l=1\\\ {\rm
gcd}(l,k)=1\end{subarray}}^{k}e^{2\pi i\,hl/k}\,,$ (60)

are the Ramanujan’s sums. Replacing Eq. (59) into Eq. (58) gives

$({\bf
C}_{D})_{i,j}=(1-\delta_{i,j})\sum_{k=1}^{\infty}\left(\frac{\mu(k)}{\phi(k)}\right)^{2}c_{k}(2|i-j|)),\quad
i,j=1,2,\dots,D\,,$ (61)

which leads us to write ${\bf C}_{D}$ as a sum of matrices,

${\bf
C}_{D}=\sum_{k=1}^{\infty}\left(\frac{\mu(k)}{\phi(k)}\right)^{2}R_{k,D}\,,$
(62)

where

$({R}_{k,D})_{i,j}=(1-\delta_{i,j})\,c_{k}(2|i-j|),\quad i,j=1,2,\dots,D\,.$
(63)

The sums $c_{k}(n)$ are multiplicative functions in $k$, namely

$c_{k_{1}k_{2}}(h)=c_{k_{1}}(h)\,c_{k_{2}}(h),\qquad{\rm
gcd}(k_{1},k_{2})=1\,.$ (64)

In particular,

$c_{2k}(2h)=c_{2}(2h)\,c_{k}(2h)=c_{k}(2h)\,,\qquad{\rm if}\;k\;{\rm odd}\,,$
(65)

where we have used that $c_{2}(h)=\cos(\pi h)$, for integer $h$. Equation (65)
implies that

${R}_{2k,D}={R}_{k,D},\qquad{\rm if}\;k\;{\rm odd}\,.$ (66)

The Möebius and Euler’s totient functions are also multiplicative, hence

$\mu(2k)=-\mu(k),\qquad\phi(2k)=\phi(k),\qquad{\rm if}\;k\;{\rm odd}\,,$ (67)

where we have used that $\mu(2)=-1,\phi(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

${\bf C}_{D}=2\sum_{k=1({\rm
odd})}^{\infty}\left(\frac{\mu(k)}{\phi(k)}\right)^{2}R_{k,D}\,.$ (68)

Our aim is to find the eigenvalues of the matrix ${\bf 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

$D=p_{2}\dots p_{L+1}\,,$ (69)

where $p_{2}=3,p_{3}=5\dots$ . Second, we truncate the sum in Eq. (68) to the
integers $k$ that divide $D$, obtaining the matrix

$\tilde{{\bf
C}}_{D}=2\sum_{k|D}\left(\frac{\mu(k)}{\phi(k)}\right)^{2}R_{k,D}\,.$ (70)

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

$R_{k,D}=\tilde{R}_{k,D}-\phi(k)\;{\bf I}_{D}\,,$ (71)

where ${\bf I}_{D}$ is the identity matrix of dimension $D$, and

$(\tilde{R}_{k,D})_{j,j^{\prime}}=c_{k}(2(j-j^{\prime}))=\sum_{\begin{subarray}{c}l=1\\\
{\rm gcd}(l,k)=1\end{subarray}}^{k}e^{4\pi i\,l(j-j^{\prime})/k}\;.$ (72)

The negative term in Eq. (71) guarantees that $(\tilde{R}_{k,D})_{j,j}=0$,
because

$\sum_{\begin{subarray}{c}l=1\\\ {\rm
gcd}(l,k)=1\end{subarray}}^{k}1=\phi(k)\,.$ (73)

Replacing Eq. (71) into Eq. (70) gives

$\displaystyle\tilde{{\bf C}}_{D}$ $\displaystyle=$ $\displaystyle
2\sum_{k|D}\left(\frac{\mu(k)}{\phi(k)}\right)^{2}\tilde{R}_{k,D}-2\,\sum_{k|D}\frac{\mu^{2}(k)}{\phi(k)}{\bf
1}_{D}$ $\displaystyle=$ $\displaystyle
2\sum_{k|D}\left(\frac{\mu(k)}{\phi(k)}\right)^{2}\tilde{R}_{k,D}-2\,\frac{D}{\phi(D)}{\bf
1}_{D}\,,$

where we used the relation [16]

$\sum_{d|n}\frac{\mu^{2}(d)}{\phi(d)}=\frac{n}{\phi(n)}\,.$ (75)

It turns out that the matrices $\tilde{R}_{k,D}$ commute among themselves, so
they can be diagonalized simultaneously and, using Eq. (C), we can obtain the
spectrum of $\tilde{{\bf C}}_{D}$. This fact follows from the following
periodicity property of the Ramanujan sums, i.e.

${\rm if}\;k|D\Longrightarrow c_{k}(2(h+D))=c_{k}(2h)\,,$ (76)

that can be easily proved using Eq. (60). Hence, the eigenvalues,
$r_{k,D}(a)\;(a=1,\dots,D)$, of $\tilde{R}_{k,D}$ are given by the Fourier
transform of its entries,

$\displaystyle r_{k,D}(a)$ $\displaystyle=$
$\displaystyle\sum_{h=1}^{D}c_{k}(2h)e^{-2\pi
iha/D}=\sum_{\begin{subarray}{c}l=1\\\ {\rm
gcd}(l,k)=1\end{subarray}}^{k}\sum_{h=1}^{D}\omega_{l,a}^{h},\qquad\omega_{l,a}=e^{2\pi
i\,\left(\frac{2l}{k}-\frac{a}{D}\right)}\,.$ (77)

Let us first notice that $\omega_{l,a}$ is a $D$ root of unity (use that
$k|D$),

$\omega_{l,a}^{D}=e^{2\pi i\left(\frac{2lD}{k}-a\right)}=1,\quad\forall
l,\;\forall a\,.$ (78)

We have to consider two cases, when $\omega_{l,a}=1$ and when $\omega_{l,a}$
is a primitive root of unity. The former case occurs if $a$ is given in terms
of $l$ as

$\omega_{l,a}=1\Longleftrightarrow a=\frac{2lD}{k}\;({\rm
mod}\;D)\Longrightarrow r_{k,D}(a)=D\,.$ (79)

Taking into account that ${\rm gcd}(l,k)=1$, there are $\phi(k)$ eigenvalues
of this type. In the case where $\omega_{l,a}$ is a primitive root of unity,
one has $\sum_{h}^{D}\omega_{l,a}^{h}=0$, and then

$\omega_{l,a}\neq 1\Longleftrightarrow a\neq\frac{2lD}{k}\;({\rm
mod}\;D)\Longrightarrow r_{k,D}(a)=0\,.$ (80)

The number of vanishing eigenvalues is therefore $D-\phi(k)$. This proves that
the matrices $\tilde{R}_{k,D}$ are diagonal in the Fourier variables $a$. This
result, together with Eq. (C), will allow us to obtain the eigenvalues of
$\tilde{{\bf C}}_{D}$. Indeed, for a given $a=1,\dots,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. (C) with the corresponding
weight $(\mu(k)/\phi(k))^{2}$.

The previous computation is simplified by the following remarkable fact: for
each $a=1,\dots,D$, there is only one pair $(l,k)$ where condition (79) holds,
while for the remaining pairs $(l^{\prime},k^{\prime})$ it is the condition
(80) that is satisfied. The proof is as follows. Let $a$ be an integer
satisfying (79) for two pairs $(l_{1},k_{1})$ and $(l_{2},k_{2})$, that is,

$a=\frac{2l_{1}D}{k_{1}}\;({\rm mod}\,D)=\frac{2l_{2}D}{k_{2}}\;({\rm
mod}\,D)\,.$ (81)

This equation implies

$\left(\frac{2l_{1}}{k_{1}}-\frac{2l_{2}}{k_{2}}\right)D=0\;({\rm
mod}\,D)\Longrightarrow
2\left(\frac{l_{1}}{k_{1}}-\frac{l_{2}}{k_{2}}\right)\in{\mathbb{Z}}\,.$ (82)

Since $1\leq l_{i}<k_{i}\;(i=1,2)$, the term
$\frac{l_{1}}{k_{1}}-\frac{l_{2}}{k_{2}}$ belongs to the interval $(-1,1)$.
Hence Eq. (82) holds in three cases,

$2\left(\frac{l_{1}}{k_{1}}-\frac{l_{2}}{k_{2}}\right)=0,\pm
1\,\Longrightarrow 2(l_{1}k_{2}-l_{2}k_{1})=0,\pm k_{1}k_{2}\,.$ (83)

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

$l_{1}k_{2}=l_{2}k_{1}\Longrightarrow
x_{1}k_{2}=k_{1}\;\&\;k_{2}=k_{1}x_{1}\Longrightarrow
x_{1}x_{2}=1\Longrightarrow x_{1}=x_{2}=1\Longrightarrow
l_{1}=l_{2}\;\&\;k_{1}=k_{2}\,,$ (84)

which means that the decomposition $a=2lD/k\;({\rm mod}D)$ is unique, as
stated above. Table 2 shows the eigenvalues $r_{k,D}(a)$ for the case $D=15$.

$k\setminus a$ | 1 | 2 | 3 | 4 | 5 | 6 | 7 | 8 | 9 | 10 | 11 | 12 | 13 | 14 | 15
---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---
$1$ | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 1
$3$ | 0 | 0 | 0 | 0 | 1 | 0 | 0 | 0 | 0 | 1 | 0 | 0 | 0 | 0 | 0
$5$ | 0 | 0 | 1 | 0 | 0 | 1 | 0 | 0 | 1 | 0 | 0 | 1 | 0 | 0 | 0
$15$ | 1 | 1 | 0 | 1 | 0 | 0 | 1 | 1 | 0 | 0 | 1 | 0 | 1 | 1 | 0
Table 2: The eigenvalues $r_{k,D}(a)/D$ for $D=15$. The four rows corresponds
to the odd integers $k=1,3,5,15$ that divide $D$. The entry is 1 when $a$ is
given by Eq. (79). Figure 13: Plot of the eigenvalues of the matrix ${\bf
C}_{D}$: $\gamma_{i,D}$ for $i=2,\dots,15$ (red dots), and those of
$\tilde{\bf C}_{D}$: $\tilde{\gamma}_{k,D}$ (blue dots) for $D=15$ and
$k=3,5,15$. The highest eigenvalue (not plotted) is given by
$\gamma_{1,D}=25.0742$ and $\tilde{\gamma}_{1,D}=26.25$.

In summary, the spectrum of the matrix $\tilde{R}_{k,D}$ is given by

${\rm Spec}\;\tilde{R}_{k,D}=\left\\{\begin{array}[]{cc}D,&\phi(k)\\\
0,&D-\phi(k)\\\ \end{array}\right.\,.$ (85)

where the first column denotes the eigenvalue and the second one its
degeneracy. Using Eq. (71) we obtain the spectrum of ${R}_{k,D}$,

${\rm Spec}\;{R}_{k,D}=\left\\{\begin{array}[]{cc}D-\phi(k),&\phi(k)\\\
-\phi(k),&D-\phi(k)\\\ \end{array}\right.\,.$ (86)

which shows explicitly that it is a traceless matrix. Finally, using Eq. (C),
together with the fact that only one matrix $\tilde{R}_{k,D}$ contributes with
$D$ to the total sum, we obtain the eigenvalues of the matrix $\tilde{\bf
C}_{D}$,

$\tilde{\gamma}_{k,D}=2D\left[\frac{1}{\phi^{2}(k)}-\frac{1}{\phi(D)}\right],\quad
k|D$ (87)

which has degeneracy $\phi(k)$. We have simplified the notation using the fact
that $\mu(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 matrix
$\tilde{\bf C}_{D}$, using the relation [16]

$\sum_{k|D}\phi(k)=D\,.$ (88)

The matrix $\tilde{\bf C}_{D}$ was introduced in order to obtain a simple
analytic approximation to the eigenvalues $\gamma_{i,D}$ of the matrix ${\bf
C}_{D}$. Fig. 13 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 $\phi(k)$, which in this case are $1,2,4,8$. Equation
(87) gives the positive and negative values of $\tilde{\gamma}_{k,D}$. The
change of sign occurs at a value $k_{D}$ satisfying

$\phi^{2}(k_{D})\simeq\phi(D)\,.$ (89)

### The Prime state

Figure 14: Left: plot of $A(k)/k^{2}$. The straight line is $\alpha=0.173304$.
Right: plot of $B(k)-\frac{1}{2}(\ln k+\beta)$ (see Eq. (99)). Both plots are
in the interval $k\in(5000,10000)$.

Equation (87) leads us to propose the approximation given in Eq. (22) for the
eigenvalues of the matrix ${\bf C}_{m}$

$\gamma_{k}\simeq
2^{m}\mu^{2}(k)\left(\frac{1}{\phi^{2}(k)}-\frac{1}{\phi_{m}}\right),\;k=1,3,\dots,k_{m}\,,$
(90)

with a degeneracy $\phi(k)$. The parameter $k_{m}$ is fixed imposing that the
sum of the degeneracies is the dimension of the matrix ${\bf C}_{m}$, that is,

$A(k_{m})\equiv\sum_{k=1({\rm odd})}^{k_{m}}\mu^{2}(k)\phi(k)\simeq
2^{m-1}\,,$ (91)

where $k_{m}$ is chosen as the highest square free odd integer. The asymptotic
behaviour of the sum is given by (see Fig. 14-left)

$A(k)\simeq\;\alpha k^{2}+{\rm fluctuations}\quad(k\gg 1)\,,$ (92)

with

$\alpha=\frac{2}{5}\prod_{p}\left(1+\frac{p-1}{p^{2}}\right)\left(1-\frac{1}{p}\right)=0.171299873\dots$
(93)

In this equation the product runs over all the prime numbers $\\{p\\}$. This
result follows from the equation[1] 11footnotetext: We thank A. Córdoba and F.
Chamizo for providing us this result. See also [51, 52].

$\sum_{k=1}^{n}\mu^{2}(k)\phi(k)\simeq\;\;\frac{1}{2}\prod_{p}\left(1+\frac{p-1}{p^{2}}\right)\left(1-\frac{1}{p}\right)+{\rm
fluctuations}\quad(n\gg 1)\,.$ (94)

Replacing Eq. (92) into Eq. (91) gives the asymptotic behaviour of $k_{m}$,

$k_{m}\simeq\left[\frac{2^{m/2}}{\sqrt{2\alpha}}\right]\quad(m\gg 1)\,,$ (95)

where $[x]$ denotes the integer part of $x$. To find $\phi_{m}$ we impose the
vanishing of the trace of ${\bf C}_{m}$,

$0={\rm tr}\,{\bf C}_{m}=\sum_{k=1({\rm odd})}^{k_{m}}\phi(k)\gamma_{k}\simeq
2^{m}\sum_{k=1({\rm
odd})}^{k_{m}}\mu^{2}(k)\left(\frac{1}{\phi(k)}-\frac{\phi(k)}{\phi_{m}}\right)\,,$
(96)

that yields

$\phi_{m}=\frac{A(k_{m})}{B(k_{m})},$ (97)

where

$B(k_{m})\equiv\sum_{k=1({\rm odd})}^{k_{m}}\frac{\mu^{2}(k)}{\phi(k)}\,.$
(98)

The asymptotic behaviour of this sum is given by (see Fig. 14-right)

$B(k)\simeq\frac{1}{2}(\ln k+\beta)+o(k^{-1/2})\,,$ (99)

with

$\beta=\gamma+\frac{1}{2}\ln 2+\sum_{p}\frac{\ln p}{p(p-1)}=1.679135304\dots,$
(100)

where $\gamma$ is the Euler’s constant. This result can be derived from the
equation[1]

$\sum_{k=1}^{n}\frac{\mu^{2}(k)}{\phi(k)}\simeq\;\;\ln
k+\gamma+\sum_{p}\frac{\ln p}{p(p-1)}+o(n^{-1/2})\quad(n\gg 1)\,.$ (101)

Using Eqs. (91), (95) and (99) we find

$\phi_{m}\simeq\frac{2^{m+1}}{m\ln
2+\delta},\qquad\delta=2\beta-\ln(2\alpha)=4.42946304048\dots$ (102)

Figure 15: Left: plot of $S(n)$ using Eq. (104) in the interval
$n/2\in[13,35]$. The straight line is the linear fit of the points, given by
$0.886793\frac{n}{2}-0.972872$. Righ: plot of $S(n)-S(n-2)$ in the interval
$n\in[28,70]$. The dotted line is $H(3/\pi^{2})$.

The eigenvalues of the density matrix in Eq. (55) can be approximated using
Eq. (90)

$\displaystyle\lambda_{k}$ $\displaystyle=$
$\displaystyle\frac{1}{d}\,|\mu(k)|(1+\ell_{N}\gamma_{k}),\qquad\qquad
k=1,3,\dots,k_{m}\,$ $\displaystyle=$ $\displaystyle
2^{1-m}|\mu(k)|\left(1+\frac{2^{m-1}}{m\ln
2}\left(\frac{1}{\phi^{2}(k)}-\frac{1}{\phi_{m}}\right)\right)\,,$

that yields the von Neumann entropy,

$\displaystyle S(n)$ $\displaystyle=$ $\displaystyle-\sum_{k=1({\rm
odd)}}^{k_{m}}\,\phi(k)\,\lambda_{k}\log_{2}\lambda_{k}\,,$ (104)

where $n=2m$. Fig. 15-left shows the linear behaviour of $S(n)$ with a slope
that seems to converge to the value $H(3/\pi^{2})$, as shown in Fig. 15-right.

##  Appendix D Fourier expansion of the entropy of the Prime state for natural
bi-partitions

Given a natural bi-partition with the first (i.e. least significant) $m$
qubits of the $n$-qubit Prime state, $|\mathbb{P}_{n}{\rangle}$, its von
Neumann entropy is given by a function of the form

$S(m,n)=n\,f\left(\frac{m}{n}\right)\,,$ (105)

where $f\left(\frac{m}{n}\right)$ can be expressed as a Fourier series,

$f\left(\frac{m}{n}\right)=\sum_{k=0}^{\infty}a_{k}\sin\left(2\pi
k\,\frac{m}{n}\right)+b_{k}\cos\left(2\pi k\,\frac{m}{n}\right)\,.$ (106)

Alternatively, Eq. (105) can be normalized by the entropy of the half chain,
$S(\frac{n}{2})$, as in Eq. (31). This results in a Fourier series analogous
to that in Eq. (106), but with the coefficients divided by
$f(\frac{1}{2})=0.4430\dots$ Such latter coefficients,
$\\{a^{\prime}_{k},b^{\prime}_{k}\\}$, obtained numerically for
$k=0,1,\dots,7$, are presented in Table.3. The resulting truncated Fourier
series, compared to the exact values appearing in Fig. 7, is shown in Fig. 16.

$k$ | $a^{\prime}_{k}$ | $b^{\prime}_{k}$
---|---|---
0 | 0 | 5.38987832e-01
1 | -6.41605273e-02 | -4.46006651e-01
2 | -8.36285927e-03 | -2.72896597e-02
3 | -1.80403548e-02 | -3.91333186e-02
4 | -4.62380210e-03 | -6.47395477e-03
5 | -8.37937907e-03 | -9.22366149e-03
6 | -3.03787769e-03 | -2.01306807e-03
7 | -4.11451014e-03 | -2.29208566e-03
Table 3: Coefficients of the Fourier expansion of the entropy of the Prime
state $|\mathbb{P}_{n}{\rangle}$, normalized by the entropy of half chain, for
natural bi-partitions of size $m$, and for $k=0,1,\dots,7$. Figure 16:
Comparison of the truncated Fourier series obtained with the coefficients
given in Table 3, and the exact entropies plotted in Fig. 7.

##  Appendix E Reduced density matrix of the odd composite state

We provide here the derivation of the asymptotic expression for the reduced
density matrix of the odd composite state upon tracing out the last (i.e. most
significant) $n-m$ qubits. We begin with the diagonal elements,

$\rho_{ii}=\frac{2^{n-m}-\pi_{M,i}(N)}{N/2-\pi(N)}\quad,\quad i\;\;odd\,,$
(107)

where $N=2^{n}$ and $M=2^{m}$. The off-diagonal elements are given by

$\rho_{ij}=\frac{2^{n-m}-\pi_{M;\,i,j}(N)}{N/2-\pi(N)}\quad,\quad
i,j\;\;odd\,.$ (108)

In the limit $n,m\rightarrow\infty$, Eq. (107) and Eq. (108) become

$\rho_{ii}=\frac{2^{n-m}-N/(\phi(M)\ln\,N)}{N/2-N/\ln\,N}\quad,\quad i\;\;odd\
\,,$ (109)

and

$\rho_{ij}=\frac{2^{n-m}-NC(|i-j|)/(\phi(M)\ln^{2}\,N)}{N/2-N/\ln\,N}\quad,\quad
i,j\;\;odd\,,$ (110)

respectively. Substituting $\phi(M)=2^{m-1}$ and rearranging terms, we obtain

$\rho_{ii}=\frac{2^{-m}\,(1-2/(n\ln\,2))}{1/2-1/(n\ln\,2)}\quad,\quad
i\;\;odd\ \,,$ (111)

and

$\rho_{ij}=\frac{2^{-m}(1-2C(|i-j|)/(n^{2}\ln^{2}\,2))}{1/2-1/(n\ln\,2)}\quad,\quad
i,j\;\;odd\,.$ (112)

Combining Eq. (111) and Eq. (112), we obtain, in the limit
$n\rightarrow\infty$, the expression in the main text,

$\overline{\rho}_{A}=\frac{1}{d}\left(\mathbbm{1}+{\bf P}_{m}\right)\,.$ (113)