The Prime state and its quantum relatives
Diego Garc
´
ıa-Mart
´
ın
1,2,3
, Eduard Ribas
2
, Stefano Carrazza
4,5
, Jos
´
e I. Latorre
5,2,6
, and
Germ
´
an Sierra
3
1
Barcelona Supercomputing Center (BSC), Barcelona, Spain.
2
Dept. F
´
ısica Qu
`
antica i Astrof
´
ısica, Universitat de Barcelona, Barcelona, Spain.
3
Instituto de F
´
ısica Te
´
orica, UAM-CSIC, Madrid, Spain.
4
TIF Lab, Dipartimento di Fisica, Universit
`
a degli Studi di Milano and INFN Milan, Milan, Italy.
5
Quantum Research Centre, Technology Innovation Institute, Abu Dhabi, UAE.
6
Center for Quantum Technologies, National University of Singapore, Singapore.
The Prime state of n qubits, |P
n
i, is de-
fined as the uniform superposition of all the
computational-basis states corresponding to
prime numbers smaller than 2
n
. This state en-
codes, quantum mechanically, arithmetic prop-
erties of the primes. We first show that the
Quantum Fourier Transform of the Prime state
provides a direct access to Chebyshev-like bi-
ases in the distribution of prime numbers. We
next study the entanglement entropy of |P
n
i up
to n = 30 qubits, and find a relation between
its scaling and the Shannon entropy of the den-
sity of square-free integers. This relation also
holds when the Prime state is constructed us-
ing 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 arith-
metic progressions. Finally, we explore the
properties of other number-theoretical quan-
tum states, such as those defined from odd
composite numbers, square-free integers and
starry primes. For this study, we have devel-
oped an open-source library that diagonalizes
matrices using floats of arbitrary precision.
1 Introduction
Mathematical sequences of integers are key to Num-
ber Theory and other areas of Mathematics. Proper-
ties of such sequences can be studied on a quantum
computer by creating pure states consisting of super-
positions 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 quan-
tum 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 quan-
tum 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 func-
tion π(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 informa-
tion 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 func-
tions 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].
Accepted in Quantum 2020-11-28, click title to verify. Published under CC-BY 4.0. 1
arXiv:2005.02422v3 [quant-ph] 3 Dec 2020
It is specially interesting to note that the Prime
state bears an amount of entanglement that is al-
most 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 can-
not 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 par-
titions 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 entan-
glement traits of the Prime state, and other number-
theoretical quantum states, is that it is not at all un-
reasonable to think that its quantum correlations may
be related to deep facts in Number Theory. For in-
stance, the Hardy-Littlewood constants [8] that char-
acterize pairwise correlations among primes appear in
the asymptotic expression of the reduced density ma-
trix 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 entangle-
ment may indeed help unravel deep truths in Number
Theory.
Let us mention a different line of research that
brings together Arithmetics with Quantum Mechan-
ics. The so called coprime chain [10] is a one dimen-
sional strongly correlated system where the local de-
grees of freedom are labelled by integers, and have a
nearest-neighbour interaction when they share a com-
mon 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 re-
late to Arithmetics. First, we introduce the Quan-
tum Fourier Transform (QFT) of the Prime state,
which can be efficiently computed and allows to ob-
tain Chebyshev-like biases in the prime number dis-
tribution from simple measurements in the computa-
tional basis. Then, we study the entanglement present
in the Prime state and several of its quantum rela-
tives, up to n = 30 qubits. We find that the scaling of
the entropy for bi-partitions of size
n
2
is almost iden-
tical, 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 en-
tropy of the density of square-free integers. A related
conjecture is posed regarding the intercept of the en-
tropy. 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 nu-
merical 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. Sec-
tion 2 briefly reviews the Prime state; section 3 intro-
duces its Quantum Fourier Transform; section 4 de-
scribes 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 |P
n
i of n qubits is defined as the
uniform superposition of all the computational-basis
states corresponding to prime numbers written in bi-
nary format, less than 2
n
(we assume that n ≥ 2),
i.e.
|P
n
i ≡
1
p
π(2
n
)
2
n
X
p: prime
|pi, (1)
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
|pi ≡ |p
n−1
i ⊗ ··· ⊗ |p
1
i ⊗ |p
0
i. This state was intro-
duced in Ref. [1], where it was shown that its con-
struction 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 al-
gorithm. 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
|ki⊗|ai ≡
|ki⊗X |ai if k is prime
|ki⊗|ai if k is not prime
,
(2)
Accepted in Quantum 2020-11-28, click title to verify. Published under CC-BY 4.0. 2
where |ki is a n-qubit computational-basis state, |ai 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 |ai =
1
√
2
(|0i − |1i), 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 uni-
tary amounts to code on a quantum computer a pri-
mality 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 pro-
duced 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, |φi =
1
√
2
n
P
2
n
−1
i=0
|e
i
i,
plus an ancilla in the |0i state; where |φi is created
by applying n Hadamard gates, H
⊗n
, onto the initial
|0 . . . 0i state. This yields:
r
2
n
− π(2
n
)
2
n
X
p: not prime
|pi|0i+
r
π(2
n
)
2
n
X
p: prime
|pi|1i.
(3)
By measuring the ancillary qubit in the above state,
and post-selecting whenever the result of the mea-
surement is |1i, 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 men-
tioned previously, that is, the oracle based on U
prime
is efficient. The second fact that guarantees the ef-
ficiency in the construction of the Prime state is the
relatively high abundance of primes, π(x), which is
asymptotically given by the logarithmic integral func-
tion, Li(x) ≡
R
x
2
dt
ln t
, i.e.
π(x) ∼ Li(x)
x→∞
−−−−→
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(
√
n). Indeed, this
estimate is given by
π
4
q
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 nu-
merically 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 asymptot-
ically with x as [19]
|π(x) −Li(x)| ≤
1
8π
√
x ln x for x > 2657 . (5)
The quantum algorithm proposed in Ref. [1] allows to
compute π(x) beyond the limits achievable by means
of classical computation –once a fault-tolerant univer-
sal 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 com-
putation (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 so-
lutions M to an oracle search problem; in this case,
the number of primes that are marked by Grover’s or-
acle, G. This algorithm uses Quantum Phase Estima-
tion (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 Rie-
mann Hypothesis. To achieve such precision, O(
√
2
n
)
calls to the oracle are needed. This is quadratically
better than the performance of any classical counter-
part 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) uncon-
ditionally, i.e. the validity of the estimation not de-
pending on the truthfulness of the Riemann Hypoth-
esis 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´eglise 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
). An-
other two analytic algorithms for computing π(x), put
Accepted in Quantum 2020-11-28, click title to verify. Published under CC-BY 4.0. 3
forward by Lagarias and Odlyzko [27], have time com-
plexity O(x
1/2+
) and O(x
3/5+
) respectively, with
> 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 |ψi =
P
j
x
j
|e
j
i of n qubits, expressed in the
computational basis, {|e
j
i}, is given by
QF T |ψi =
2
n
−1
X
k=0
y
k
|e
k
i , y
k
≡
1
√
2
n
2
n
−1
X
j=0
x
j
e
2π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 expo-
nent is the imaginary unit. The QFT is an efficient
subroutine, requiring O(n
2
) gates on a quantum com-
puter [17].
In the present work, we calculate the QFT of the
Prime state, |
ˆ
P
n
i ≡ QF T |P
n
i, that is,
|
ˆ
P
n
i =
1
p
2
n
π(2
n
)
2
n
−1
X
k=0
X
p: prime
e
2πi pk/2
n
|ki. (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 |ki and |2
n
− ki is the same. Indeed:
P (2
n
− k) =
1
2
n
π(2
n
)
X
p: prime
e
2πi p
e
−2πi pk/2
n
2
(8)
=
1
2
n
π(2
n
)
X
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
i is that
the QFT applied to the Prime state provides a method
for estimating Chebyshev-like biases in the distribu-
tion of prime numbers. These biases reflect the unbal-
ance in the number of primes, below a certain value,
appearing in different arithmetic progressions.
Figure 1: Numerical simulation of the Quantum Fourier
Transform of a Prime state with n = 15 qubits. The prob-
ability 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
i is not
included in the Hilbert space of n qubits.
Let us therefore consider arithmetic progressions of
the form αk + β, for k = 0, 1, 2, . . . , with α, β co-
primes (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 count-
ing function, π(x). The asymptotic behaviour of these
modular counting functions is determined by
π
α,β
(x) ∼
1
φ(α)
Li(x)
x→∞
−−−−→
1
φ(α)
x
ln x
, (9)
where φ(α) is the Euler’s totient function, which gives
the number of coprimes to α smaller than α. For ex-
ample, 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, asymptoti-
cally 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 func-
tions
∆
α; β
1
,β
2
(x) ≡ π
α,β
1
(x) − π
α,β
2
(x) . (10)
These biases can be numerically estimated, to a de-
sired precision , from the probability peaks appear-
ing in |
ˆ
P
n
i. 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 deriva-
tion of these results). The expressions for these peaks
Accepted in Quantum 2020-11-28, click title to verify. Published under CC-BY 4.0. 4
are,
P (N/3) =
∆
3;2,1
(N)
2
+ π
3,1
(N)π
3,2
(N) − π(N) + 2
Nπ(N)
,
P (N/4) =
1 + ∆(N)
2
Nπ(N)
,
P (N/6) =
∆
6;5,1
(N)
2
+ π
6,1
(N)π
6,5
(N) − 3 π
6,1
(N) + 3
Nπ(N)
,
(11)
where ∆(N) ≡ π
4,3
(N) − π
4,1
(N) is the Chebyshev
bias. There exists a direct relation between the func-
tions π
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
0
), and p = 6k
0
+ 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
0
+ 1), so
p = 6k
0
+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,
π
3,1
(x) = π
6,1
(x) , π
3,2
(x) = π
6,5
(x) + 1 , (12)
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 equa-
tions 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 can-
not be directly measured in the computational basis.
Instead, one needs to resort to ancillary qubits to rep-
resent binary fractions and hence rational numbers;
the QFT acts upon the two registers, the one repre-
senting 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 mea-
sure 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). Ac-
cording to Eq. (4) and Eq. (9), these probabilities
will then decrease as ∼
1
φ(α)
2
1
n ln 2
, where α = 3, 6 re-
spectively. 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 prob-
abilities, {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 measure-
ments performed, #M, as
i
=
q
p
i
(1−p
i
)
#M
. Hence,
in order to achieve a small multiplicative error ε
i
,
say ε
i
≈ 0 and
i
<< p
i
, then the number of mea-
surements, or runs of the quantum circuit, must be
#M >>
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(
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
i 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/3i,
N−1
X
k=0
k6=N/3
y
k
|e
k
i|0i +
p
P (N/3) |N/3i|1i. (13)
States like the one in Eq. (13) are ready for Amplitude
Estimation. Classically, the best algorithm for com-
puting modular prime counting functions, π
α,β
(x), is
a variation of the Meisser-Lehmer method for comput-
ing π(x), and has time complexity O(x
2/3
ln
−2
x) [34].
To the best of our knowledge, modular prime count-
ing 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 denomi-
nator Nπ(N) is expected to cause an exponential de-
crease 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
i, all
the information contained in the probability peaks
is extracted simultaneously, rather than sequentially,
by accumulation of statistical knowledge on measure-
ment 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-
Accepted in Quantum 2020-11-28, click title to verify. Published under CC-BY 4.0. 5
weighted superposition provides a means to count
the number of terms in the computational basis. In-
deed, the probability of finding the |0 . . . 0i 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 ba-
sis). A derivation of this result can be found in Ap-
pendix B. For |
ˆ
P
n
i, this probability reads
P (0) =
π(N)
N
, (14)
which decreases as ∼
1
n ln 2
in the limit n 1, ac-
cording 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 expres-
sion for this peak is
P (N/2) =
π(N)
2
− 4π(N) + 4
Nπ(N)
, (15)
which also decreases linearly in n
−1
.
Because the precision required to meaningfully test
Riemann’s hypothesis, Eq. (5), implies that the mul-
tiplicative error must become smaller and smaller
for increasing x (as
√
x ln x
π(x)
decreases), direct sam-
pling from |
ˆ
P
n
i is not competitive with previously
mentioned quantum counting algorithms, requiring
O(2
n
n
−3
) repetitions. On the other hand, it is possi-
ble to prepare a state like that in Eq. (13), but mark-
ing 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
i for estimating π(2
n
).
The QFT may nonetheless be useful to estimate
the number of terms in an equally-weighted super-
position, for states for which one does not know of
any oracle; for instance, a ground state of a Hamilto-
nian obtained from a variational algorithm such as the
Variational Quantum Eigensolver [29] or the Quan-
tum 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 entangle-
ment, 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 avail-
able 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 superposi-
tion. In turn, these observables may reveal informa-
tion 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
ρ
A
=
1
d
( + `
N
C
m
) , (16)
where d = 2
m−1
, `
N
≡
Li
2
(N)
Li(N)
N→∞
−−−−→
1
n ln 2
, and C
m
is a d × d symmetric Toeplitz matrix defined by
(C
m
)
i,j
≡ (1−δ
i,j
) C(2|i−j|), i, j = 1, . . . , d , (17)
where C(h) are the Hardy-Littlewood constants, de-
fined by
C(h) =
2 C
2
Q
p>2; p|h
p−1
p−2
if h is even
0 otherwise
, (18)
with C
2
=
Q
p>2
1 −
1
(p−1
2
)
the twin prime con-
stant 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 con-
jecture [8] assures that the number of prime pairs
(p, p + k) up to a certain number x, i.e. π
k
(x), is
given by
π
k
(x) ∼ C(k) Li
2
(x)
x→∞
−−−−→ C(k)
x
(ln x)
2
, (19)
where Li
2
(x) =
R
x
2
dt
(ln t)
2
.
An analytical approximation to the positive eigen-
values of the matrix C
m
was conjectured in Ref.
[7]. That approximation was suitable for the high-
est eigenvalues, but failed to describe the lowest ones.
Moreover, the degeneracies of the former failed be-
yond some point [37]. A recent article gives a heuris-
tic derivation of the Hardy-Littlewood conjecture (19)
starting from the pair correlation formula for the Rie-
mann zeros, including lower order terms [38]. The
Accepted in Quantum 2020-11-28, click title to verify. Published under CC-BY 4.0. 6
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 & 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.
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,
C(h) =
∞
X
k=1
µ(k)
φ(k)
2
c
k
(h) (20)
where φ(x) is the Euler’s totient function, µ(x) is the
M¨obius function (see below), and
c
k
(h) =
k
X
l=1
gcd(l,k)=1
e
2π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 pro-
cesses, 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
γ
k
' 2
m
µ
2
(k)
1
φ
2
(k)
−
1
φ
m
, k = 1, 3, 5, . . . , k
m
,
(22)
where the M¨obius function, µ(x), takes values ±1 for
square-free numbers and 0 otherwise. The values of
k
m
and φ
m
are fixed by the dimension and normal-
ization 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 dif-
ferent values of k, k
0
, in which case the degeneracy is
higher. An example of such an accidental degener-
acy is bγ
13
= bγ
21
=
1
144
, which has a total degeneracy
φ(13) + φ(21) = 24 (we denote by bγ
k
= µ
2
(k)/φ
2
(k)).
Another example is bγ
35
= bγ
39
=
1
576
, that has a de-
generacy φ(35) + φ(39) = 48.
Using Eq. (22), the eigenvalues of the density ma-
trix (16) are approximated by
λ
k
' 2
1−m
µ
2
(k)
1 +
2
m
n ln 2
1
φ
2
(k)
−
1
φ
m
,
(23)
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 entangle-
ment spectrum are plotted. These correspond to the
highest eigenvalues of the density matrix of the Prime
state, since the entanglement spectrum {ε
k
} of a den-
sity matrix ρ with eigenvalues {λ
k
} is defined as
ε
k
= −ln λ
k
. (24)
Using again Eq. (22), we can estimate the contri-
bution of the positive eigenvalues to the trace of the
powers of C
m
, from which all R´enyi entropies can be
readily derived [43]. The computation is as follows,
Tr C
s
m
'
k
m
X
k=1 (odd)
φ(k) γ
s
k
m1
−−−−→ 2
ms
∞
X
k=1 (odd)
|µ(k)|
(φ(k))
2s−1
= 2
ms
Y
p>2
1 +
|µ(p)|
(φ(p))
2s−1
= 2
ms
Y
p>2
1 +
1
(p − 1)
2s−1
, (25)
where we have used that for any function f(k) which
is multiplicative, it holds that
∞
X
k=1 (odd)
|µ(k)|f(k) =
Y
p>2
(1 + f(p)) , (26)
Accepted in Quantum 2020-11-28, click title to verify. Published under CC-BY 4.0. 7
Figure 3: top) Ratio of Tr (C
s
m
) to Eq. (25), for s = 2, 3, 4, 5
and m =
n
2
, up to n = 32 qubits. bottom) Natural logarithm
of (1− ratio of Tr (C
s
m
) to Eq. (25)), up to n = 32 qubits.
where the product runs over the odd prime num-
bers. 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
N, 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 eigen-
values 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 contri-
bution of the negative eigenvalues of C
m
to the trace
of the powers of this matrix is asymptotically negligi-
ble, 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 ρ,
Figure 4: top) Value of c
π
(n) computed according to Eq.
(29), compared to the conjectured value H
3
π
2
, up to n =
30 qubits; bottom) Value of γ(n) computed according to
Eq. (30), up to n = 30 qubits, compared to the conjectured
value 1 +
3
π
2
.
with eigenvalues {λ
i
}, is defined by
S = −Tr (ρ log
2
ρ) = −
X
i
λ
i
log
2
λ
i
, (27)
and it constitutes a relevant measure of bipartite en-
tanglement. The von Neumann entropy of the re-
duced 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 indi-
cates 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. Random states, i.e. those with random complex
coefficients, follow a “volume law”, with an entangle-
ment 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/π) n/2 = 0.363 . . . n/2 [45]. So the
Accepted in Quantum 2020-11-28, click title to verify. Published under CC-BY 4.0. 8
Figure 5: Natural logarithm of the absolute value of the
difference δ, between the observed and predicted values of c
π
and γ, up to n = 30 qubits. That is, δ ≡
c
π
(n) − H
3
π
2
and δ ≡
γ(n) −
1 +
3
π
2
, respectively. An exponential
decrease in δ is observed before reaching a plateau.
Prime state is not a typical random positive state ei-
ther. 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].
We herein conjecture that the entropy of the nat-
ural equal-sized bi-partition of the Prime state of n
qubits is asymptotically given by
S(n) = c
π
n
2
−γ , c
π
= H
3
π
2
, γ = 1+
3
π
2
,
(28)
where H(p) ≡ −p log
2
(p) − (1 − p) log
2
(1 − p) 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 partic-
ular, it has been shown that it gives half the topo-
logical 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 |1i, be-
cause 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
π
. 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 preci-
sion. The code is publicly available on GitHub, to-
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 b
n
2
c.
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.
gether with all the numerical results presented in this
paper [13].
Figure 4-top shows the difference in entropy be-
tween two consecutive values of S(n), for even n. That
is, the dots in the plot correspond to
c
π
(n) = S(n) − S(n − 2) , (29)
up to n = 30 qubits. Oscillatory convergence towards
the predicted value H
3
π
2
= 0.886082085 . . . is ob-
served, in good agreement with Eq. (28). The value
of the intercept γ, obtained as
γ(n) = (S(n) − S(n − 2))
n
2
− S(n) , (30)
is shown in Fig. 4-bottom, up to n = 30 qubits. Os-
cillatory 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 con-
vergence is almost identical in both cases, with an ex-
ponential approximation toward the conjectured val-
ues as the number of qubits grows, as shown in Fig.
Accepted in Quantum 2020-11-28, click title to verify. Published under CC-BY 4.0. 9
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).
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 nat-
ural equal-sized bi-partition of size b
n
2
c, 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 Neu-
mann 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 indepen-
dent 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-
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¨obius 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
¨
obius 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.
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
S(m, n)
S(n/2, n)
=
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.
4.2 Arithmetic Prime states
Other number-theoretical quantum states apart from
the Prime state will be defined below, and their en-
tanglement traits studied. Table 1 summarizes our
findings. We define arithmetic Prime states as the
uniform superpositions of primes belonging to arith-
metic progressions, i.e.
|P
n; α,β
i ≡
1
p
π
α,β
(2
n
)
2
n
X
p: prime
p=β mod α
|pi, (32)
Accepted in Quantum 2020-11-28, click title to verify. Published under CC-BY 4.0. 10
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.
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,
ρ
A
(α) =
φ(α)
d
( + l
N
C
m
(α)) , (33)
where again d = 2
m−1
, l
N
≡
Li
2
(N)
Li(N)
N→∞
−−−−→
1
n ln 2
, and
C
m
(α) is a symmetric Toeplitz matrix defined by
(C
m
(α))
i,j
≡ (1 − δ
i,j
) C(α |i − j|) . (34)
Notice that Eq. (33) generalizes the asymptotic ex-
pression of the reduced density matrix of the Prime
state. Indeed, all primes but 2, which does not con-
tribute in the asymptotic limit, belong to the arith-
metic progression 2k + 1. Thus, φ(2) = 1, C
m
(2) =
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 π
α; β,β
0
(x),
π
α; β,β
0
(x)
x→∞
−−−−→
C(|β − β
0
|)
φ(α)
Li
2
(x) , (35)
where π
α; β,β
0
is the number of prime pairs (p, p
0
) such
that p, p
0
are primes smaller than, or equal to, x of the
form p = αk + β, p = αk + β
0
, with (α, β) and (α, β
0
)
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, |P
n; α,β
i, when α is a power
of 2, is once again dictated by a straight line with
slope ∼ 0.884, see Fig. 8 and Table 1. 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 . . . , s =
5.1830 . . . , s = 7.6977 . . . and s = 10.4260 . . . , for
α = 4, 8, 16, 32 respectively.
Taking into account all the results presented in this
section, we generalize the conjecture in Eq. (28), stat-
ing 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
π
b
n
2
c − γ
α; d
, c
π
= H
3
π
2
, (37)
where bxc 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.
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 com-
posite (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
i ≡
1
p
C
c,n
2
n
−1
X
c: odd
c: composite
|ci , 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 en-
tropy remains approximately constant around a mean
Accepted in Quantum 2020-11-28, click title to verify. Published under CC-BY 4.0. 11
Figure 9: top) Value of γ
4
(n) computed according to Eq.
(30), up to n = 30 qubits, compared to the value 1 + 3 ×
3
π
2
in Eq. (36); bottom) Value of γ
8
(n) computed according
to Eq. (30), up to n = 30 qubits, compared to the value
1 + 5 ×
3
π
2
in Eq. (36).
value of ∼ 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 → ∞, by
ρ
A
=
1
d
( + P
m
) , (39)
where
(P
m
)
ij
=
(
0 if i = j
1 if i 6= j
. (40)
From Eq. (39), the von Neumann entropy S in the
asymptotic limit is given by
S ' 0 , (41)
in agreement with the observed numerical results.
One might hope that, because the entanglement en-
tropy of the odd composite state is small, there could
Figure 10: Von Neumann entropy S of the natural equal-
sized bi-partition of the odd composite, odd square-free and
M
¨
obius states, up to n = 28 qubits.
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 num-
bers are composite, and the composite state carries
their equal superposition. This very high density of
equally-weighted superpositions approaches the prod-
uct state |φi =
1
√
2
n
P
2
n
−1
i=0
|e
i
i. Thus, tensor net-
works 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 de-
scription.
Another number-theoretical quantum state, rela-
tive 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
|S
n
i ≡
1
p
C
s,n
2
n
X
s: odd
|µ(s)| |si , C
s,n
=
2
n
X
s: odd
|µ(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 ∼ 20 qubits.
Another possibility is to define a uniform superpo-
sition of square-free numbers and let the signs of the
amplitudes be given by the M¨obius function. This
state, that we shall call M¨obius state [7], is then de-
Accepted in Quantum 2020-11-28, click title to verify. Published under CC-BY 4.0. 12
fined by
|µ
n
i ≡
1
p
C
µ,n
2
n
X
s=1
µ(s) |si , C
µ,n
=
2
n
X
s=1
|µ(s)| .
(43)
The von Neumann entropy of the natural equal-sized
bi-partition of the M¨obius state scales linearly with a
slope close to 1. This means that the M¨obius 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 nat-
ural 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
≤ p
∗
n
≤ p
n+1
, (44)
where p
n
is the n-th prime. These numbers are re-
markable in that a ζ
∗
function, defined via the Euler
product formula, i.e.
ζ
∗
(s) ≡
Y
p
∗
1
1 − p
∗
−s
, (45)
has an analytic extension with exactly the same ze-
ros, and the same multiplicities, than the original Rie-
mann zeta function [50]. The starry-Prime state is
then defined by
|P
∗
n
i ≡
1
p
π(2
n
)
2
n
−1
X
p
∗
: starry
prime
|p
∗
i, (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
n
2
with
slope ∼ 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 ef-
ficiently placed in a linear superposition, the Prime
state. This cast of a classical sequence onto a quan-
tum 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 mea-
surements. We have also delved into the entangle-
ment properties of the Prime state, which turn out
to be intimately related to the correlations between
pairs of prime numbers. These correlations were pos-
tulated by Hardy and Littlewood almost a century
ago, and played a critical role in the connection be-
tween Number Theory and the Theory of Quantum
Chaos based on heuristic analogies. According to it,
prime numbers are classical objects (labelling peri-
odic orbits in a chaotic dynamical system), while the
Riemann’s zeros are quantum objects (eigenenergies
of a hypothetical Riemann Hamiltonian). In our ap-
proach, the Hardy and Littlewood conjecture provides
us with the entanglement Hamiltonian of the Prime
state, whose spectrum we have found to be charac-
terized 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 ap-
plies 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 en-
tropy of the Prime state does not depend on the basis
Accepted in Quantum 2020-11-28, click title to verify. Published under CC-BY 4.0. 13
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 diagonal-
izes 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´erez-Garc´ıa and S. Ramos-Calderer for useful dis-
cussions. DGM and JIL acknowledge
CaixaBank
for its support of this work through Barcelona
Supercomputing Center’s project CaixaBank Com-
putaci´on Cu´antica. DGM and JIL are supported by
Project PGC2018-095862-B-C22, and Project Quan-
tum CAT (001- P-001644) co-funded by Generali-
tat de Catalunya and the European Union Regional
Development Fund. SC is supported by the Euro-
pean 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´on y Univer-
sidades (grant PGC2018-095862-B-C21), the Comu-
nidad de Madrid (grant QUITEMAD+ S2013/ICE-
2801), the Centro de Excelencia Severo Ochoa Pro-
gramme (grant SEV-2016-0597) and the CSIC Re-
search Platform on Quantum Technologies PTI-001.
References
[1] J.I. Latorre and G. Sierra, “Quantum compu-
tation of prime number functions”, Quant. Inf.
Comput. 14, 577 (2014).
[2] L.K. Grover, “A fast quantum mechanical algo-
rithm 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 con-
firmed π(10
27
) prime counting func-
tion record”, Mersenne Forum (2015),
https://github.com/kimwalisch/primecount.
[5] G. Brassard, P. Høyer and A. Tapp, “Quan-
tum 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 entangle-
ment in the primes”, Quant. Inf. Comput. 15,
622 (2015).
[8] G.H. Hardy and J.E. Littlewood, “Some Prob-
lems of ’Partitio Numerorum.’ III. On the Ex-
pression of a Number as a Sum of Primes”, Acta
Math. 44, 1 (1923).
[9] We consider natural bi-partitions those that sep-
arate the first m qubits and the remainder n −m
qubits, without any reordering.
[10] G. Mussardo, H. Giudici and J. Viti, “The Co-
prime 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 Se-
quence 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, Ox-
ford (1938).
[17] M.A. Nielsen and I.L. Chuang, Quantum Com-
putation and Quantum Information: 10th An-
niversary Edition, Cambridge University Press,
New York (2011).
[18] https://www.claymath.org/millennium-
problems.
[19] L. Schoenfeld, “Sharper bounds for the Cheby-
shev functions θ(x) and ψ(x). II”, Math. Com-
put. 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 algo-
rithms revisited”, Proc. R. Soc. London A 454,
339 (1998).
[21] S. Aaronson, “Quantum Lower Bound for Ap-
proximate Counting Via Laurent Polynomials”,
arXiv:1808.02420 (2018); S. Aaronson and P.
Rall, “Quantum Approximate Counting, Simpli-
fied”, SOSA20, 24 (2020).
[22] Meissel, “
¨
Uber die Bestimmung der Primzahlen-
menge 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).
Accepted in Quantum 2020-11-28, click title to verify. Published under CC-BY 4.0. 14
[24] J.C. Lagarias, V.S. Miller and A.M. Odlyzko,
“Computing π(x): the Meissel-Lehmer method”,
Math. Comp. 44, 573 (1985).
[25] M. Del´eglise and J. Rivat, “Computing π(x):
the Meissel, Lehmer, Lagarias, Miller, Odlyzko
method”, Math. Comp. 65, 235 (1996).
[26] X. Gourdon, “Computation of π(x): improve-
ments to the Meissel, Lehmer, Lagarias, Miller,
Odllyzko, Del´eglise and Rivat method”, unpub-
lished (2001). As cited in Ref. [4].
[27] J.C. Lagarias and A.M. Odlyzko, “Computing
π(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 Algo-
rithm”, arXiv:1411.4028 (2014).
[31] If p
i
is the exact probability, and ˆp
i
the estimate
of it, the statistical additive error is defined as
p
i
= ˆp
i
+
i
, and the multiplicative error as ˆp
i
=
(1 + ε
i
) p
i
.
[32] G. Brassard, P. Høyer, M. Mosca and A. Tapp,
“Quantum Amplitude Amplification and Estima-
tion”, Contemp. Math. 305, 53 (2002).
[33] A. Montanaro, “Quantum speedup of Monte
Carlo methods”, Proc. R. Soc. A 471, 20150301
(2015).
[34] M. Del´eglise, P. Dusart and X.F. Roblot, “Count-
ing 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 en-
tanglement, 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 corre-
lations from the pair correlation of Riemann ze-
ros”, 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 cir-
cle and the sieve: Ramanujan-Fourier series”,
arXiv:math/0601574 (2006).
[43] A. R´enyi, “On measures of information and en-
tropy”, Proceedings of the fourth Berkeley Sym-
posium on Mathematics, Statistics and Probabil-
ity 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¨obius function, randomness, and
dynamics”, online notes (2011),
http://publications.ias.edu/sarnak/paper/512.
[48] F. Cellarosi and Y. G. Sinai, “Ergodic Proper-
ties 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 ζ 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, “
¨
Uber die
Summe
P
n
µ
2
(n)/ϕ(n)”, Nederl. Akad. Weten-
sch.‘ Proc. Ser. A 67, Indag. Math. 26, 582
(1964).
Accepted in Quantum 2020-11-28, click title to verify. Published under CC-BY 4.0. 15
Appendices
A Computation of relevant peaks in |
ˆ
P
n
i
We compute here the most relevant peaks of the QFT of the Prime state, |
ˆ
P
n
i. We use Eq. (6), where
x
j
=
1
√
π(N)
for prime j and 0 otherwise. Recall that N = 2
n
is the dimension of the Hilbert space.
P (N/2) =
1
Nπ(N)
X
p: prime
e
pπi
2
=
1
Nπ(N)
X
p: prime
p=0 mod2
e
pπi
+
X
p: prime
p=1 mod2
e
pπi
2
=
1
Nπ(N)
|π
2,0
(N) − π
2,1
(N)|
2
=
1
Nπ(N)
|1 − (π(N) − 1)|
2
=
π(N)
2
− 4π(N) + 4
Nπ(N)
(47)
P (N/3) =
1
Nπ(N)
X
p: prime
e
2
3
πi p
2
=
1
Nπ(N)
X
p: prime
p=0 mod3
e
2
3
pπi
+
X
p: prime
p=1 mod3
e
2
3
pπi
+
X
p: prime
p=2 mod3
e
2
3
pπi
2
=
1
Nπ(N)
π
3,0
(N) + π
3,1
(N) e
2
3
πi
+ π
3,2
(N) e
4
3
πi
2
=
1
Nπ(N)
1 +
−
1
2
+ i
√
3
2
π
3,1
(N) +
−
1
2
− i
√
3
2
π
3,2
(N)
2
=
1
Nπ(N)
1 −
1
2
(π
3,1
(N) + π
3,2
(N)) − i
√
3
2
∆
3;2,1
(N)
2
=
1
Nπ(N)
1 +
1
4
(π
3,1
(N) + π
3,2
(N))
2
− π
3,1
(N) − π
3,2
(N) +
3
4
∆
3;2,1
(N)
2
=
∆
3;2,1
(N)
2
+ π
3,1
(N)π
3,2
(N) − π(N) + 2
Nπ(N)
(48)
P (N/4) =
1
Nπ(N)
X
p: prime
e
1
2
πi p
2
=
1
Nπ(N)
X
p: prime
p=2 mod4
e
1
2
pπi
+
X
p: prime
p=1 mod4
e
2
4
pπi
+
X
p: prime
p=3 mod4
e
2
4
pπi
2
=
1
Nπ(N)
−π
4,0
(N) + π
4,1
(N) e
1
2
πi
+ π
4,3
(N) e
3
2
πi
2
=
1
Nπ(N)
|−1 + i (π
4,1
(N) − π
4,3
(N))|
2
=
1 + ∆(N)
2
Nπ(N)
(49)
Accepted in Quantum 2020-11-28, click title to verify. Published under CC-BY 4.0. 16
P (N/6) =
1
Nπ(N)
X
p: prime
e
2
6
πi p
2
=
1
Nπ(N)
X
p: prime
p=2 mod6
e
1
3
pπi
+
X
p: prime
p=3 mod6
e
1
3
pπi
+
X
p: prime
p=1 mod6
e
2
6
pπi
+
X
p: prime
p=5 mod6
e
2
6
pπi
2
=
1
Nπ(N)
π
6,0
(N) e
2
3
πi
+ π
6,3
(N)e
πi
+ π
6,1
(N) e
1
3
πi
+ π
6,5
(N) e
5
3
πi
2
=
1
Nπ(N)
−
1
2
+ i
√
3
2
− 1 +
1
2
+ i
√
3
2
π
6,1
(N) +
1
2
− i
√
3
2
π
6,5
(N)
2
=
1
Nπ(N)
π
6,1
(N) + π
6,5
(N) − 3
2
− i
√
3 (∆
6;5,1
(N) − 1)
2
2
=
1
Nπ(N)
9 + π
6,1
(N)
2
+ π
6,5
(N)
2
+ 2 π
6,1
(N)π
6,5
(N) − 6 π
6,1
(N) − 6 π
6,5
(N)
4
+
3 ∆
6;5,1
(N)
2
+ 3 − 6 ∆
6;5,1
(N)
4
=
∆
6;5,1
(N)
2
+ π
6,1
(N)π
6,5
(N) − 3 π
6,5
(N) + 3
Nπ(N)
(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.
Accepted in Quantum 2020-11-28, click title to verify. Published under CC-BY 4.0. 17
B Quantum Fourier Transform of a uniform superposition
Let us consider the definition of the QFT of a general quantum state |ψi of n qubits, expressed in the compu-
tational basis, {|e
j
},
|ψi =
2
n
−1
X
j=0
x
j
|e
j
i
QF T
−−−−−−−→
2
n
−1
X
k=0
y
k
|e
k
i, (51)
where
y
k
≡
1
√
2
n
2
n
−1
X
j=0
x
j
e
2πi jk/2
n
, (52)
and the i in the exponent is the imaginary unit. The probability of measuring the |0 . . . 0i state in Eq. (51) in
the computational basis thus is:
P (0) = |y
0
|
2
=
1
√
2
n
2
n
−1
X
j=0
x
j
e
0
2
=
1
2
n
2
n
−1
X
j=0
x
j
2
. (53)
If the state |ψi is a uniform superposition of M computational-basis vectors, i.e. x
j
= 1/
√
M for all j such
that x
j
6= 0, and N = 2
n
is the dimension of the Hilbert space, then
P (0) =
1
N
M−1
X
j=0
1
√
M
2
=
1
N
M
√
M
2
=
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.
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
¯ρ
A
=
1
d
(1 + `
N
C
m
), `
N
N→∞
−−−−→
1
n ln 2
, (55)
where C
m
is a matrix constructed with the even values of the Hardy-Littlewood constants C(2k),
(C
m
)
i,j
= (1 − δ
i,j
) C(2|i − j|), i, j = 1, 2, . . . , d = 2
m−1
. (56)
In order to justify the conjecture (22) for the eigenvalues of C
m
, we define the density matrix
ρ
D
=
1
D
(1 + `
D
C
D
), `
D
=
1
2 ln(2D)
, (57)
with
(C
D
)
i,j
= (1 − δ
i,j
) C(2|i − j|), i, j = 1, 2, . . . , 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 [38, 41, 42]
C(h) =
∞
X
k=1
µ(k)
φ(k)
2
c
k
(h) , (59)
Accepted in Quantum 2020-11-28, click title to verify. Published under CC-BY 4.0. 18
where
c
k
(h) =
k
X
l=1
gcd(l,k)=1
e
2πi hl/k
, (60)
are the Ramanujan’s sums. Replacing Eq. (59) into Eq. (58) gives
(C
D
)
i,j
= (1 − δ
i,j
)
∞
X
k=1
µ(k)
φ(k)
2
c
k
(2|i − j|)), i, j = 1, 2, . . . , D , (61)
which leads us to write C
D
as a sum of matrices,
C
D
=
∞
X
k=1
µ(k)
φ(k)
2
R
k,D
, (62)
where
(R
k,D
)
i,j
= (1 − δ
i,j
) c
k
(2|i − j|), i, j = 1, 2, . . . , 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), gcd(k
1
, k
2
) = 1 . (64)
In particular,
c
2k
(2h) = c
2
(2h) c
k
(2h) = c
k
(2h) , if k odd , (65)
where we have used that c
2
(h) = cos(πh), for integer h. Equation (65) implies that
R
2k,D
= R
k,D
, if k odd . (66)
The M¨oebius and Euler’s totient functions are also multiplicative, hence
µ(2k) = −µ(k), φ(2k) = φ(k), if k odd , (67)
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
C
D
= 2
∞
X
k=1(odd)
µ(k)
φ(k)
2
R
k,D
. (68)
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
D = p
2
. . . p
L+1
, (69)
where p
2
= 3, p
3
= 5 . . . . Second, we truncate the sum in Eq. (68) to the integers k that divide D, obtaining
the matrix
˜
C
D
= 2
X
k|D
µ(k)
φ(k)
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
=
˜
R
k,D
− φ(k) I
D
, (71)
where I
D
is the identity matrix of dimension D, and
(
˜
R
k,D
)
j,j
0
= c
k
(2(j − j
0
)) =
k
X
l=1
gcd(l,k)=1
e
4πi l(j−j
0
)/k
. (72)
Accepted in Quantum 2020-11-28, click title to verify. Published under CC-BY 4.0. 19
The negative term in Eq. (71) guarantees that (
˜
R
k,D
)
j,j
= 0, because
k
X
l=1
gcd(l,k)=1
1 = φ(k) . (73)
Replacing Eq. (71) into Eq. (70) gives
˜
C
D
= 2
X
k|D
µ(k)
φ(k)
2
˜
R
k,D
− 2
X
k|D
µ
2
(k)
φ(k)
1
D
(74)
= 2
X
k|D
µ(k)
φ(k)
2
˜
R
k,D
− 2
D
φ(D)
1
D
,
where we used the relation [16]
X
d|n
µ
2
(d)
φ(d)
=
n
φ(n)
. (75)
It turns out that the matrices
˜
R
k,D
commute among themselves, so they can be diagonalized simultaneously
and, using Eq. (74), we can obtain the spectrum of
˜
C
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) , (76)
that can be easily proved using Eq. (60). Hence, the eigenvalues, r
k,D
(a) (a = 1, . . . , D), of
˜
R
k,D
are given by
the Fourier transform of its entries,
r
k,D
(a) =
D
X
h=1
c
k
(2h)e
−2πiha/D
=
k
X
l=1
gcd(l,k)=1
D
X
h=1
ω
h
l,a
, ω
l,a
= e
2πi
(
2l
k
−
a
D
)
. (77)
Let us first notice that ω
l,a
is a D root of unity (use that k|D),
ω
D
l,a
= e
2πi
(
2lD
k
−a
)
= 1, ∀l, ∀a . (78)
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
ω
l,a
= 1 ⇐⇒ a =
2lD
k
(mod D) =⇒ r
k,D
(a) = D . (79)
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
P
D
h
ω
h
l,a
= 0, and then
ω
l,a
6= 1 ⇐⇒ a 6=
2lD
k
(mod D) =⇒ r
k,D
(a) = 0 . (80)
The number of vanishing eigenvalues is therefore D − φ(k). This proves that the matrices
˜
R
k,D
are diagonal
in the Fourier variables a. This result, together with Eq. (74), will allow us to obtain the eigenvalues of
˜
C
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
.
The previous computation is simplified by the following remarkable fact: for each a = 1, . . . , D, there is only
one pair (l, k) where condition (79) holds, while for the remaining pairs (l
0
, k
0
) 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 =
2l
1
D
k
1
(mod D) =
2l
2
D
k
2
(mod D) . (81)
Accepted in Quantum 2020-11-28, click title to verify. Published under CC-BY 4.0. 20
k \ 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).
2 4 6 8 10 12 14
5
10
γ
i,15
, γ
˜
k,15
Figure 13: Plot of the eigenvalues of the matrix C
D
: γ
i,D
for i = 2, . . . , 15 (red dots), and those of
˜
C
D
: ˜γ
k,D
(blue dots) for
D = 15 and k = 3, 5, 15. The highest eigenvalue (not plotted) is given by γ
1,D
= 25.0742 and ˜γ
1,D
= 26.25.
This equation implies
2l
1
k
1
−
2l
2
k
2
D = 0 (mod D) =⇒ 2
l
1
k
1
−
l
2
k
2
∈ Z . (82)
Since 1 ≤ l
i
< k
i
(i = 1, 2), the term
l
1
k
1
−
l
2
k
2
belongs to the interval (−1, 1). Hence Eq. (82) holds in three
cases,
2
l
1
k
1
−
l
2
k
2
= 0, ±1 =⇒ 2(l
1
k
2
− l
2
k
1
) = 0, ±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
=⇒ x
1
k
2
= k
1
& k
2
= k
1
x
1
=⇒ x
1
x
2
= 1 =⇒ x
1
= x
2
= 1 =⇒ l
1
= l
2
& k
1
= k
2
, (84)
which means that the decomposition a = 2lD/k (modD) is unique, as stated above. Table 2 shows the
eigenvalues r
k,D
(a) for the case D = 15.
In summary, the spectrum of the matrix
˜
R
k,D
is given by
Spec
˜
R
k,D
=
(
D, φ(k)
0, D − φ(k)
. (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
,
Spec R
k,D
=
(
D −φ(k), φ(k)
−φ(k), D −φ(k)
. (86)
which shows explicitly that it is a traceless matrix. Finally, using Eq. (74), together with the fact that only
one matrix
˜
R
k,D
contributes with D to the total sum, we obtain the eigenvalues of the matrix
˜
C
D
,
˜γ
k,D
= 2D
1
φ
2
(k)
−
1
φ(D)
, k|D (87)
Accepted in Quantum 2020-11-28, click title to verify. Published under CC-BY 4.0. 21
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 matrix
˜
C
D
,
using the relation [16]
X
k|D
φ(k) = D . (88)
The matrix
˜
C
D
was introduced in order to obtain a simple analytic approximation to the eigenvalues γ
i,D
of
the matrix 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 φ(k), which in this case are 1, 2, 4, 8. Equation
(87) gives the positive and negative values of ˜γ
k,D
. The change of sign occurs at a value k
D
satisfying
φ
2
(k
D
) ' φ(D) . (89)
The Prime state
Equation (87) leads us to propose the approximation given in Eq. (22) for the eigenvalues of the matrix C
m
γ
k
' 2
m
µ
2
(k)
1
φ
2
(k)
−
1
φ
m
, k = 1, 3, . . . , k
m
, (90)
with a degeneracy φ(k). The parameter k
m
is fixed imposing that the sum of the degeneracies is the dimension
of the matrix C
m
, that is,
A(k
m
) ≡
k
m
X
k=1(odd)
µ
2
(k)φ(k) ' 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) ' αk
2
+ fluctuations (k 1) , (92)
with
α =
2
5
Y
p
1 +
p − 1
p
2
1 −
1
p
= 0.171299873 . . . (93)
In this equation the product runs over all the prime numbers {p}. This result follows from the equation
[1]
n
X
k=1
µ
2
(k)φ(k) '
1
2
Y
p
1 +
p − 1
p
2
1 −
1
p
+ fluctuations (n 1) . (94)
Replacing Eq. (92) into Eq. (91) gives the asymptotic behaviour of k
m
,
k
m
'
2
m/2
√
2α
(m 1) , (95)
5000 6000 7000 8000 9000 10 000
k
0.1706
0.1708
0.1710
0.1712
0.1714
0.1716
0.1718
A(k)/k
2
6000 7000 8000 9000 10000
k
-0.0010
-0.0005
0.0005
B(k)-
1
2
(log k + β)
Figure 14: Left: plot of A(k)/k
2
. The straight line is α = 0.173304. Right: plot of B(k) −
1
2
(ln k + β) (see Eq. (99)). Both
plots are in the interval k ∈ (5000, 10000).
1
We thank A. C´ordoba and F. Chamizo for providing us this result. See also [51, 52].
Accepted in Quantum 2020-11-28, click title to verify. Published under CC-BY 4.0. 22
15 20 25 30 35
n/2
15
20
25
30
S(n)
30 40 50 60 70
n
0.80
0.85
0.90
0.95
1.00
1.05
1.10
S(n)-S(n-2)
Figure 15: Left: plot of S(n) using Eq. (104) in the interval n/2 ∈ [13, 35]. The straight line is the linear fit of the points,
given by 0.886793
n
2
− 0.972872. Righ: plot of S(n) − S(n − 2) in the interval n ∈ [28, 70]. The dotted line is H(3/π
2
).
where [x] denotes the integer part of x. To find φ
m
we impose the vanishing of the trace of C
m
,
0 = tr C
m
=
k
m
X
k=1(odd)
φ(k)γ
k
' 2
m
k
m
X
k=1(odd)
µ
2
(k)
1
φ(k)
−
φ(k)
φ
m
, (96)
that yields
φ
m
=
A(k
m
)
B(k
m
)
, (97)
where
B(k
m
) ≡
k
m
X
k=1(odd)
µ
2
(k)
φ(k)
. (98)
The asymptotic behaviour of this sum is given by (see Fig. 14-right)
B(k) '
1
2
(ln k + β) + o(k
−1/2
) , (99)
with
β = γ +
1
2
ln 2 +
X
p
ln p
p(p − 1)
= 1.679135304 . . . , (100)
where γ is the Euler’s constant. This result can be derived from the equation
[1]
n
X
k=1
µ
2
(k)
φ(k)
' ln k + γ +
X
p
ln p
p(p − 1)
+ o(n
−1/2
) (n 1) . (101)
Using Eqs. (91), (95) and (99) we find
φ
m
'
2
m+1
m ln 2 + δ
, δ = 2β − ln(2α) = 4.42946304048 . . . (102)
The eigenvalues of the density matrix in Eq. (55) can be approximated using Eq. (90)
λ
k
=
1
d
|µ(k)|(1 + `
N
γ
k
), k = 1, 3, . . . , k
m
(103)
= 2
1−m
|µ(k)|
1 +
2
m−1
m ln 2
1
φ
2
(k)
−
1
φ
m
,
that yields the von Neumann entropy,
S(n) = −
k
m
X
k=1(odd)
φ(k) λ
k
log
2
λ
k
, (104)
Accepted in Quantum 2020-11-28, click title to verify. Published under CC-BY 4.0. 23
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/π
2
), as shown in Fig. 15-right.
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, |P
n
i, its
von Neumann entropy is given by a function of the form
S(m, n) = n f
m
n
, (105)
where f
m
n
can be expressed as a Fourier series,
f
m
n
=
∞
X
k=0
a
k
sin
2πk
m
n
+ b
k
cos
2πk
m
n
. (106)
Alternatively, Eq. (105) can be normalized by the entropy of the half chain, S(
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 (
1
2
) = 0.4430 . . . Such
latter coefficients, {a
0
k
, b
0
k
}, obtained numerically for k = 0, 1, . . . , 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
0
k
b
0
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 |P
n
i, normalized by the entropy of half chain,
for natural bi-partitions of size m, and for k = 0, 1, . . . , 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.
Accepted in Quantum 2020-11-28, click title to verify. Published under CC-BY 4.0. 24
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,
ρ
ii
=
2
n−m
− π
M,i
(N)
N/2 − π(N)
, i odd , (107)
where N = 2
n
and M = 2
m
. The off-diagonal elements are given by
ρ
ij
=
2
n−m
− π
M; i,j
(N)
N/2 − π(N)
, i, j odd . (108)
In the limit n, m → ∞, Eq. (107) and Eq. (108) become
ρ
ii
=
2
n−m
− N/(φ(M) ln N)
N/2 − N/ ln N
, i odd , (109)
and
ρ
ij
=
2
n−m
− NC(|i − j|)/(φ(M) ln
2
N)
N/2 − N/ ln N
, i, j odd , (110)
respectively. Substituting φ(M) = 2
m−1
and rearranging terms, we obtain
ρ
ii
=
2
−m
(1 − 2/(n ln 2))
1/2 − 1/(n ln 2)
, i odd , (111)
and
ρ
ij
=
2
−m
(1 − 2C(|i − j|)/(n
2
ln
2
2))
1/2 − 1/(n ln 2)
, i, j odd . (112)
Combining Eq. (111) and Eq. (112), we obtain, in the limit n → ∞, the expression in the main text,
ρ
A
=
1
d
( + P
m
) . (113)
Accepted in Quantum 2020-11-28, click title to verify. Published under CC-BY 4.0. 25