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
n1
2
n1
+ ··· + p
1
2
1
+ p
0
2
0
, and
|pi |p
n1
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
|kiX |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
(1p
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 >>
1p
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,
N1
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 AB 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
m1
, `
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|ij|), 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
p1
p2
if h is even
0 otherwise
, (18)
with C
2
=
Q
p>2
1
1
(p1
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
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 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
1m
µ
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 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))
2s1
= 2
ms
Y
p>2
1 +
|µ(p)|
(φ(p))
2s1
= 2
ms
Y
p>2
1 +
1
(p 1)
2s1
, (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