# The time problem in primordial perturbations

Alice Boldrin<EMAIL_ADDRESS>National Centre for Nuclear Research,
Pasteura 7, 02-093 Warszawa, Poland  Przemysław Małkiewicz
<EMAIL_ADDRESS>National Centre for Nuclear Research,
Pasteura 7, 02-093 Warszawa, Poland  Patrick Peter<EMAIL_ADDRESS>${\cal G}\mathbb{R}\varepsilon\mathbb{C}{\cal O}$ – Institut d’Astrophysique
de Paris, CNRS & Sorbonne Université, UMR 7095 98 bis boulevard Arago, 75014
Paris, France

###### Abstract

We study the non-unitary relation between quantum gravitational models defined
using different internal times. We show that despite the non-unitarity, it is
possible to provide a prescription for making unambiguous, though restricted,
physical predictions independent of specific clocks. To illustrate this
result, we employ a model of quantum gravitational waves in a quantum
Friedmann universe.

##  I Introduction

Quantum gravity models suffer from the infamous time problem Kuchar (2011);
Isham (1993); Anderson (2010); Kiefer and Peter (2022) as the external and
absolute time on which non-relativistic physics is based, is absent in
Einstein’s theory of gravity. Therefore, one has to rely on largely arbitrary
physical variables, known as internal time variables or internal clocks, to
follow changes occurring in gravitational systems. By virtue of the principle
of general relativity (time-reparametrization invariance), the free choice of
internal time variable has no physical consequence in the classical theory.
Upon passing to quantum theory, however, different choices of internal time
variables are known to produce unitarily inequivalent quantum models Deser et
al. (1961); Beluardi and Ferraro (1995); Hajicek (2000); Catren and Ferraro
(2001); Hajicek and Kijowski (2000); Małkiewicz (2015); Małkiewicz and
Miroszewski (2017); Bojowald and Halnon (2018). The problem of finding the
correct interpretation of these non-equivalent models is commonly referred to
as the time problem.

In this article, we look for the most plausible interpretation of such non-
equivalent clocks. Our analysis is based on the model of primordial
gravitational waves propagating across the Friedmann universe. It is important
to note that similar models were previously used for making predictions for
the primordial amplitude spectrum of density perturbation, which are greatly
constrained by observations (see, e.g., Refs. Peter et al. (2007); Pinto-Neto
(2008); Peter and Pinto-Neto (2008)). Remarkably, to the best of our
knowledge, the time problem has never been studied for such models, so it is
important to clarify the role and the interpretation of internal time
variables in their dynamics. We expect that the ensuing conclusions should
equally apply to all cosmological models.

The fact that the dynamics are unitarily inequivalent in different clocks is
widely known and well-documented with plenty of examples, see e.g. Refs.
Malkiewicz (2017); Małkiewicz et al. (2020, 2022). In this context, it is
sometimes emphasized that the only measurable quantities in quantum gravity
are gauge-invariant variables that do not depend on the employed clock Rovelli
(2002); Dittrich (2007). They are constants of motion. However, they are said
to encode all the relational dynamics in spite of being non-dynamical
themselves. From this point of view, dynamical quantities are not fundamental
and are ambiguously given by one-parameter families of gauge-invariant
quantities, with each family representing the motion with respect to a
specific internal time. The differences are seen as natural rather than
inconsistencies that should be worried about. From our viewpoint, on the other
hand, the dynamical variables can serve as fundamental variables, and the
differences in their dynamics call for a careful interpretation, before
allowing for physical predictions.

The cosmological system examined in this article exhibits, as expected,
dynamical discrepancies when based on different clocks. The discrepancies
concern both the background and perturbation variables. This leads us to ask a
fundamental question: what are the dynamical predictions of quantum
cosmological models, which do not depend on the employed time variable?

We address the above question within the reduced phase space quantization.
Namely, we solve the Hamiltonian constraint and choose the internal time
variable prior to quantization. An alternative approach would be to first
quantize and then solve the constraint quantum mechanically while promoting
one of the variables as internal time. Both approaches lead to the same time
problem and therefore using the technically less involved reduced phase
approach is well-justified (see, however, Ref. Bojowald and Tsobanjan (2023)
for recent developments in the alternative approach). Most significantly,
within the reduced phase space approach, there exists a theory of clock
transformations, which is completely crucial for the purpose of this work
Małkiewicz (2017). Thanks to these precisely defined transformations, we are
able to explore all possible clocks and quantize them with an assumption of
fixed operator ordering. Hence, any quantum ambiguities found arise from the
differences between clocks rather than the differences between quantization
prescriptions.

The outline of this article is as follows. In Sec. II, we make a brief
introduction to the theory of clock transformations in the reduced phase space
of gravitational models. We explain how this theory allows to remove
irrelevant quantization ambiguities when passing to quantum theory based on
different clocks and with different basic dynamical variables. In Sec. III, we
formulate the reduced phase space description of the Friedmann universe with
gravitational waves with respect to a fluid time and obtain the general clock
transformation. In Sec. IV, we quantize our model and establish a convenient
semi-classical approximation. Sec. V deals with concrete clock transformations
applied to our model and makes comparisons between the resulting dynamics. We
summarize our findings, discuss their plausible interpretation and suggest
some directions to move forward in the final Sec. VI.

##  II Clock transformations in totally constrained systems

One crucial characteristic feature of canonical relativity is the appearance
of the Hamiltonian constraint; it is a consequence of the fact that the
dynamics of three-surfaces is generated by infinitesimal time-like
diffeomorphisms, and the latter leave the full four-dimensional spacetime
invariant. It by no means makes the dynamics of three-surfaces spurious, or
redundant. Indeed, the Hamiltonian constraint dynamics is a feature of any
canonical relativistic theory of gravity, be it Einstein’s or any modified
gravity theory, though their dynamics are different. The correct
interpretation of canonical relativity assumes the lack of an absolute,
external time in which three-surfaces evolve, and replaces it with internal
variables that serve as clocks in which the dynamics of three-surfaces takes
place. None of the internal clocks can play a privileged role as the principle
of relativity states. This picture is certainly self-consistent in the
classical theory. At the quantum level, no spacetime exists and, as we will
see later, the principle of relativity takes a somewhat altered form. In order
to study it, we need to extend the canonical formalism by including clock
transformations that transform a canonical description from one internal clock
to another; only then can we move to the quantum level where these new
transformations become a key to unlock the principle of quantum relativity.

Let us consider a system consisting of a set of $N+1$ canonical variables
$\\{q_{\alpha},p^{\alpha}\\}_{\alpha=0,\cdots,N}$ and assume a Hamiltonian
constraint taking the form

$C(q_{\alpha},p^{\alpha})\approx 0,$

where ‘$\approx$’ is the weak equality in the Dirac sense Dirac (1964).
Suppose that one of the positions, say $q_{0}$, varies monotonically with the
evolution generated by the constraint, i.e. $\forall
q_{0},\\{q_{0},C\\}_{\textsc{pb}}\neq 0$. It is then possible to assign to
$q_{0}$ the role of an internal clock in which the evolution of the remaining
variables occurs. This evolution is then governed by a Hamiltonian that is not
a constraint. At this stage, it may seem that the time variable is fixed once
and for all, which would contradict the principle of relativity; we discuss
below in what sense this is not the case.

The reduced Hamiltonian formalism is obtained from the initial symplectic form
$\Sigma=\mathrm{d}q_{\alpha}\wedge\mathrm{d}p^{\alpha}$ (Einstein convention
assumed), evaluated on the constraining surface, namely

$\begin{gathered}\Sigma\big{|}_{C=0}=\left(\mathrm{d}q_{I}\wedge\mathrm{d}p^{I}+\mathrm{d}q_{0}\wedge\mathrm{d}p^{0}\right)\big{|}_{C=0}\\\
\hbox{}=\mathrm{d}q_{I}\wedge\mathrm{d}p^{I}-\mathrm{d}t\wedge\mathrm{d}H,\end{gathered}$
(1)

where $I=1,\cdots,N$, and $H=H\left(q_{0},q_{I},p^{I}\right)$ is the
nonvanishing reduced Hamiltonian such that $p_{0}+H\approx 0$. Note that both
$q_{0}$ (denoted by $t$ from now on to emphasize its role as a time variable)
and $p^{0}$ are removed from the phase space and the remaining dynamical
variables are no longer constrained. Indeed, their dynamics reads

$\frac{\mathrm{d}q_{I}}{\mathrm{d}t}=\frac{\partial H}{\partial
p^{I}}\qquad\hbox{and}\qquad\frac{\mathrm{d}p^{I}}{\mathrm{d}t}=-\frac{\partial
H}{\partial q_{I}},$

which is entirely solved once an arbitrary initial condition
$(q_{I}^{\text{ini}},p^{I}_{\text{ini}},q_{0}^{\text{ini}})$ is provided.

In order to restore the principle of relativity, we need to allow for any
clock, denoted by $\tilde{t}$, which monotonically varies with the evolution
generated by the constraint, $\\{\tilde{t},C\\}_{\textsc{pb}}\neq 0$. This new
clock must be a function of the old clock and the old canonical variables,
$\tilde{t}=\tilde{t}(q_{I},p^{I},t)$. Thus, it must satisfy

$\frac{\mathrm{d}\tilde{t}}{\mathrm{d}t}=\frac{\partial\tilde{t}}{\partial
t}+\underbrace{\frac{\partial\tilde{t}}{\partial q_{I}}\frac{\partial
H}{\partial p^{I}}-\frac{\partial\tilde{t}}{\partial p^{I}}\frac{\partial
H}{\partial q_{I}}}_{\\{\tilde{t},H\\}_{\textsc{pb}}}\neq 0.$ (2)

The original symplectic form induced on the constraint surface $C=0$ must read
in some new canonical variables:

$\Sigma\big{|}_{C=0}=\mathrm{d}\tilde{q}_{I}\wedge\mathrm{d}\tilde{p}^{I}-\mathrm{d}\tilde{t}\wedge\mathrm{d}\tilde{H},$

so that the new reduced formalism is still canonical. This implies that there
must exist an invertible map between the old and the new variables:

$\tilde{t}=\tilde{t}(q_{I},p^{I},t),\
\tilde{q}_{I}=\tilde{q}_{I}(q_{J},p^{J},t),\
\tilde{p}^{I}=\tilde{p}^{I}(q_{J},p^{J},t),$ (3)

and the natural question to ask is whether these transformations are
canonical. In principle, and in all the relevant cases, they most certainly
are not. It can be shown that clock transformations form a group of generally
noncanonical transformations with canonical transformations as its normal
subgroup Malkiewicz (2017); finding them is in general a difficult task.
However, for an integrable dynamical system, the problem can be reduced to
that of solving a set of algebraic equations.

If a dynamical system is integrable, then we may find a complete set of
canonical constants of motion, denoted by $D_{I}$. Let them be functions of
the old internal time and old canonical variables,
$D_{I}=D_{I}(q_{J},p^{J},t)$. Note that substituting back $t\rightarrow
q_{0}$, they must commute with the original constraint,
$\\{D_{I},C\left(q_{\alpha},p^{\alpha}\right)\\}_{\textsc{pb}}=0$. They are
therefore genuine Dirac observables in the constrained system. The new
internal time $\tilde{t}=\tilde{t}(q_{I},p^{I},t)$ and new canonical variables
can then be found according to the algebraic relations

$\displaystyle\tilde{t}=\tilde{t}(q_{I},p^{I},t),\leavevmode\nobreak\
\leavevmode\nobreak\
D_{I}(q_{J},p^{J},t)=D_{I}(\tilde{q}_{J},\tilde{p}^{J},\tilde{t}),$ (4)

where we formally substitute the canonical variables in the expressions for
Dirac observables $D_{I}$, i.e., we assume the same functional dependence of
$D_{I}$ in both sets of variables. The number of $D_{I}$ is equal to the
number of the new canonical variables $\tilde{q}_{J}$ and $\tilde{p}^{J}$, and
thus, leaving aside singular cases, the above relations determine
$\tilde{q}_{J}$ and $\tilde{p}^{J}$ completely. The result is a new canonical
formalism based on a new internal clock. Let us note that, by virtue of Eq.
(4), if a solution to the dynamics is known in one clock, i.e.
$t\rightarrow\left[q_{I}\left(D_{J},t\right),p^{I}\left(D_{J},t\right)\right]$,
then it is readily known for all other clocks and reads
$\tilde{t}\rightarrow\left[\tilde{q}_{I}=q_{I}\left(D_{J},\tilde{t}\right),\tilde{p}^{I}=p^{I}\left(D_{J},\tilde{t}\right)\right]$.
This makes the choice of the new canonical variables $\tilde{q}_{I}$ and
$\tilde{p}^{I}$ via Eq. (4) very convenient: the formal description of the
system is the same in all clocks, only the physical meaning of the clock and
basic variables changes, which is emphasized by the use of a tilde (
$\tilde{}$ ) over the variable names.

The use of Dirac observables in the derivation of clock transformations gives
an invaluable advantage when passing to quantum theory. Our goal is to make a
comparison between quantum theories based on different internal clocks of a
single physical system. Therefore, it is of uttermost importance to make sure
that the theories are different only insofar as their clocks differ, and not
due to other quantization ambiguities such as the well-known factor ordering.
This state of affairs can be achieved by fixing a quantum representation of
the Dirac observables and then defining basic and compound observables as
functions of the quantum Dirac observables, both in the original

$\widehat{q}_{I}=q_{I}(\widehat{D}_{J},t),\qquad\widehat{p}^{I}=p^{I}(\widehat{D}_{J},t),$

and the new variables

$\widehat{\tilde{q}}_{I}={q}_{I}(\widehat{D}_{J},\tilde{t}),\qquad\widehat{\tilde{p}}^{I}={p}^{I}(\widehat{D}_{J},\tilde{t}).$

These definitions imply that ${q}_{I}$ and ${p}^{I}$ are promoted to the same
operators as ${\tilde{q}}_{I}$ and ${\tilde{p}}^{I}$, respectively. We invert
this reasoning and start by assuming the same operators for ${q}_{I}$ and
${\tilde{q}}_{I}$ as well as ${p}^{I}$ and ${\tilde{p}}^{I}$. This implies
that the Dirac observables being the same functions in both sets of basic
variables are promoted to the same operators irrespectively of the choice of
clock. Hence, the quantum descriptions in different clocks are formally the
same, only the physical meaning of the basic operators changes from one clock
to another, which is emphasized by the use of tilde. Obviously, a unique
ordering prescription has to be used in all the above formulas. In principle,
after this step, any physically interesting aspect of the quantum theories can
be compared. In the following section, we introduce the model on which we
discuss such comparisons.

##  III Canonical cosmological model

We consider a flat Friedmann-Lemaître-Robertson-Walker (FLRW) universe filled
with radiation and perturbed by gravitational waves; the line element of the
model reads (in units such that $c=1$)

$\mathrm{d}s^{2}=-N^{2}(t)\mathrm{d}t^{2}+a^{2}(t)\left[\delta_{ij}+h_{ij}(\bm{x},t)\right]\mathrm{d}x^{i}\mathrm{d}x^{j},$

where $h_{ij}$ represent the gravitational waves (tensor perturbations); it
satisfies $h_{ij}\delta^{ij}=0$ and $\partial^{j}h_{ij}=0$. Finally, we assume
a toroidal spatial topology with each comoving coordinate $x^{i}\in[0,1)$.
Setting $N\to a$ means one considers the conformal time; we shall henceforth
denote it by $\eta$ to agree with most of the cosmology literature.

###  III.1 Perturbative Hamiltonian

Let us now build the canonical description of these gravitational waves in an
FLRW universe. The relevant canonical variables are the scale factor $a$ and
its conjugate momentum $p_{a}$ to describe the background, while the tensor
perturbations are represented by the gravitational wave amplitude
$\mu^{(\lambda)}=ah^{(\lambda)}$ and its conjugate momentum $\pi^{(\lambda)}$,
with $\lambda\in\\{+,\times\\}$ and
$h_{ij}=\sum_{\lambda}h^{(\lambda)}\varepsilon_{ij(\lambda)}$ (see, e.g. Refs.
Peter and Uzan (2013) and Micheli and Peter (2022) for details on the helicity
expansion).

The matter component is assumed to be a radiation fluid with energy density
$p_{0}$ conjugate to a timelike variable $q_{0}$. The gravitational constraint
is expanded to second order through

$H_{\text{tot}}=H^{\text{(b)}}+\sum_{\bm{k}}H^{\text{(p)}}_{\bm{k}},$

(recall the spatial sections are compact), with the background Hamiltonian
given by

$H^{\text{(b)}}=-\frac{1}{2}p_{a}^{2}-p_{0}.$ (5)

At this stage, one can identify the internal time $q_{0}$ with the conformal
time $\eta$ as it reduces the zeroth order Hamiltonian into

$\begin{gathered}\Sigma\big{|}_{H^{(0)}=0}=\left(\mathrm{d}a\wedge\mathrm{d}p_{a}+\mathrm{d}q_{0}\wedge\mathrm{d}p_{0}\right)\big{|}_{H^{\text{(b)}}=0}\\\
=\mathrm{d}a\wedge\mathrm{d}p_{a}-\mathrm{d}\eta\wedge\mathrm{d}\left(\frac{1}{2}p_{a}^{2}\right),\end{gathered}$
(6)

leading to the physical zeroth-order Hamiltonian

$H^{(0)}=\frac{1}{2}p_{a}^{2},$ (7)

while preserving the form of the perturbation Hamiltonian
$H^{\text{(p)}}_{\bm{k}}$. The latter reads, at second order

$H^{\text{(p)}}_{\bm{k}}\to
H^{(2)}_{\bm{k}}=-\sum_{\lambda=+,\times}H^{(2)}_{\bm{k},\lambda}$ (8)

with

$H^{(2)}_{\bm{k},\lambda}=\frac{1}{2}\left|\pi^{(\lambda)}_{\bm{k}}\right|^{2}+\frac{1}{2}\left(k^{2}-\frac{a^{\prime\prime}}{a}\right)\left|\mu^{(\lambda)}_{\bm{k}}\right|^{2},$
(9)

where a prime stands for a derivative with respect to the conformal time.
Since the tensor perturbations are real, one has
$\mu^{(\lambda)*}_{\bm{k}}=\mu^{(\lambda)}_{-\bm{k}}$. Moreover, since the
background is isotropic, one can restrict attention to upward directed
wavevectors $\bm{k}$ by merely cancelling the factor $\frac{1}{2}$ in
$H^{(2)}_{\bm{k},\lambda}$. This permits to write the final second-order
Hamiltonian as

$H^{(2)}_{\bm{k},\lambda}=\pi^{(\lambda)}_{\bm{k}}\pi^{(\lambda)}_{-\bm{k}}+\left(k^{2}-\frac{a^{\prime\prime}}{a}\right)\mu^{(\lambda)}_{\bm{k}}\mu^{(\lambda)}_{-\bm{k}}.$
(10)

Note that for the radiation fluid we are concerned with here, the Hamiltonian
(7) yields as equations of motion $p_{a}=a^{\prime}$ and $p_{a}^{\prime}=0$,
thus leading to $a^{\prime\prime}=0$: the potential for producing
gravitational waves is indeed classically vanishing if the universe is
radiation dominated.

###  III.2 Dirac observables

Now we shall find the constants of motion that form canonical pairs. To this
end, we need to solve the partial differential equations

$\frac{\mathrm{d}D}{\mathrm{d}\eta}=\frac{\partial
D}{\partial\eta}+\\{D,H^{(0)}+H^{(2)}\\}_{\textsc{pb}}=0.$ (11)

At zeroth order, this is

$\frac{\partial D}{\partial\eta}+p_{a}\frac{\partial D}{\partial a}=0,$

with solutions

$D_{1}=a-p_{a}\eta\qquad\hbox{and}\qquad D_{2}=p_{a}.$ (12)

At first order, Eq. (11) reads

$\frac{\partial\delta D}{\partial\eta}+p_{a}\frac{\partial\delta D}{\partial
a}=\pi^{(\lambda)}_{\bm{k}}\frac{\partial\delta
D}{\partial\mu^{(\lambda)}_{\bm{k}}}-k^{2}\mu^{(\lambda)}_{\bm{k}}\frac{\partial\delta
D}{\partial\pi^{(\lambda)}_{\bm{k}}},$

where we considered the classical solution $a^{\prime\prime}=0$. Since we are
considering only first order perturbations, we demand that $\delta D$ be
linear in the perturbation variables $\mu^{(\lambda)}_{\bm{k}}$ and
$\pi^{(\lambda)}_{\bm{k}}$. The l.h.s. of the above equation is greatly
simplified if $\delta D$ depends only on the variable $y=\eta+a/p_{a}$, so we
look for a solution of the form $\delta
D^{(\lambda)}=\mu^{(\lambda)}_{\bm{k}}\alpha(y)+\pi^{(\lambda)}_{\bm{k}}\beta(y)$,
leading to

$2\frac{\mathrm{d}\alpha}{\mathrm{d}y}\mu^{(\lambda)}_{\bm{k}}+2\frac{\mathrm{d}\beta}{\mathrm{d}y}\pi^{(\lambda)}_{\bm{k}}=\alpha\pi^{(\lambda)}_{\bm{k}}-k^{2}\beta\mu^{(\lambda)}_{\bm{k}}.$

Assuming independent variations of $\mu^{(\lambda)}_{\bm{k}}$ and
$\pi^{(\lambda)}_{\bm{k}}$, one gets
$2\mathrm{d}\alpha/\mathrm{d}y=-k^{2}\beta$ and
$2\mathrm{d}\beta/\mathrm{d}y=\alpha$, and finally
$4\mathrm{d}^{2}\alpha/\mathrm{d}y^{2}=k^{2}\alpha$, so that, setting

$\Omega_{k}=\frac{k}{2}\left(\eta+\frac{a}{p_{a}}\right),$

one gets two independent solutions for each polarization, or in other words
four first-order constants, reading

$\begin{gathered}\delta
D^{(\lambda)}_{1,\bm{k}}=\sqrt{k}\sin\Omega_{k}\,\mu^{(\lambda)}_{\bm{k}}-\frac{\cos\Omega_{k}}{\sqrt{k}}\,\pi^{(\lambda)}_{\bm{k}},\\\
\delta
D^{(\lambda)}_{2,\bm{k}}=\sqrt{k}\cos\Omega_{k}\,\mu^{(\lambda)}_{\bm{k}}+\frac{\sin\Omega_{k}}{\sqrt{k}}\,\pi^{(\lambda)}_{\bm{k}}.\end{gathered}$
(13)

In Eq. (13), the normalisation has been chosen so as to ensure that all these
Dirac observables indeed form canonical pairs, namely

$\left\\{D_{1},D_{2}\right\\}_{\textsc{pb}}=1\qquad\hbox{and}\qquad\left\\{\delta
D^{(\lambda)}_{1,\bm{k}},\delta
D^{(\bar{\lambda})}_{2,\bm{k}}\right\\}_{\textsc{pb}}=\delta_{\lambda\bar{\lambda}}.$

From now on, we drop the index $\lambda$ and consider just a single
polarization mode $(\mu_{\bm{k}},\pi_{\bm{k}})$.

###  III.3 Clock transformations

Having set the full model, and before moving on to its quantum counterpart,
let us fist consider a general clock transformation

$\eta\rightarrow\tilde{\eta}=\eta+\Delta(a,p_{a},\eta),$ (14)

where $\Delta$ is a delay function that in general varies between the
trajectories as well as along them. At the background level, implementing the
recipe given by Eq. (3), i.e.,
$D_{1,2}(a,p_{a},\eta)=D_{1,2}(\tilde{a},\tilde{p}_{a},\tilde{\eta})$ to the
transformation (14), yields

$a-p_{a}\eta=\tilde{a}-\tilde{p}_{a}\tilde{\eta}\quad\hbox{and}\quad
p_{a}=\tilde{p}_{a},$

leading to

$\tilde{a}=a+p_{a}\Delta\quad\hbox{and}\quad p_{a}=\tilde{p}_{a}.$ (15)

In order that the clock transformation (14) actually defines a new and
physically acceptable clock, the delay function $\Delta$ must be subject to
two conditions. First, the new clock must run forward, that is

$\frac{\mathrm{d}\tilde{\eta}}{\mathrm{d}\eta}=1+\frac{\mathrm{d}\Delta}{\mathrm{d}\eta}=1+\frac{\partial\Delta}{\partial\eta}+p_{a}\frac{\partial\Delta}{\partial
a}>0,$ (16)

where in the second equality we used the zeroth-order Hamiltonian $H^{(0)}$
given by Eq. (7) and the associated equations of motion.

The second condition that a clock transformation must satisfy is that the
ranges of the basic variables $a$ and $p_{a}$ must be preserved, thereby
preventing the appearance of non-trivial ranges that may induce new and
potentially unsolvable quantization issues. This second condition implies

$\displaystyle\lim_{p_{a}\rightarrow\pm\infty}\tilde{p}_{a}(a,p_{a},\eta)$
$\displaystyle=\pm\infty,$ (17a)
$\displaystyle\tilde{a}(a,p_{a},\eta)\big{|}_{a=0}$ $\displaystyle=0.$ (17b)

The first equality (17a) is trivially satisfied in the present case because of
(15). For $\Delta=\Delta(a,p_{a})$, to which we shall restrict attention in
what follows, the second equality (17b) is identical to demanding that the
delay function at vanishing scale factor should also vanish,
$\Delta(0,p_{a})=0$. This condition also ensures that the slow-gauge clock is
transformed into another slow-gauge clock, that is, the boundary is reached
within a finite amount of time (see Ref. Małkiewicz et al. (2020)). Such a
condition (17b), although irrelevant in the classical theory, is crucial for
the existence of a bounce at the quantum level where the clock must smoothly
connect contracting and expanding trajectories. Were (17) violated, the clock
transformations would break the bouncing trajectories.

It turns out that the condition (16) is equivalent to the existence of a one-
to-one map between the reduced phase spaces $(a,p_{a})$ and
$(\tilde{a},\tilde{p}_{a})$, i.e. the determinant

$\frac{\partial\left(\tilde{a},\tilde{p}_{a}\right)}{\partial\left(a,p_{a}\right)}=\left|\begin{array}[]{cc}\displaystyle\frac{\partial\tilde{a}}{\partial
a}&\displaystyle\frac{\partial\tilde{a}}{\partial p_{a}}\\\ &\\\
\displaystyle\frac{\partial\tilde{p}_{a}}{\partial
a}&\displaystyle\frac{\partial\tilde{p}_{a}}{\partial
p_{a}}\end{array}\right|>0,$ (18)

which is indeed Eq. (16) when $\partial\Delta/\partial\eta=0$.

At first order, one must solve

$\delta D_{1}(a,p_{a},\mu_{\bm{k}},\pi_{\bm{k}})=\delta
D_{1}(\tilde{a},\tilde{p}_{a},\tilde{\mu}_{\bm{k}},\tilde{\pi}_{\bm{k}})$

and

$\delta D_{2}(a,p_{a},\mu_{\bm{k}},\pi_{\bm{k}})=\delta
D_{2}(\tilde{a},\tilde{p}_{a},\tilde{\mu}_{\bm{k}},\tilde{\pi}_{\bm{k}})$

in order to determine the clock-transformed perturbation variables.
Explicitely, using (13), one gets

$\begin{gathered}\sqrt{k}\sin\Omega_{k}\,\mu_{\bm{k}}-\frac{\cos\Omega_{k}}{\sqrt{k}}\,\pi_{\bm{k}}=\sqrt{k}\sin\tilde{\Omega}_{k}\,\tilde{\mu}_{\bm{k}}-\frac{\cos\tilde{\Omega}_{k}}{\sqrt{k}}\,\tilde{\pi}_{\bm{k}},\\\
\sqrt{k}\cos\Omega_{k}\,\mu_{\bm{k}}+\frac{\sin\Omega_{k}}{\sqrt{k}}\,\pi_{\bm{k}}=\sqrt{k}\cos\tilde{\Omega}_{k}\,\tilde{\mu}_{\bm{k}}+\frac{\sin\tilde{\Omega}_{k}}{\sqrt{k}}\,\tilde{\pi}_{\bm{k}},\end{gathered}$
(19)

where
$\tilde{\Omega}_{k}=\frac{1}{2}k(\tilde{\eta}+\tilde{a}/\tilde{p}_{a})=\Omega_{k}+k\Delta$.
The above algebraic equations (19) can easily be inverted to yield the new
canonical perturbation variables, namely

$\left(\begin{array}[]{c}\tilde{\mu}_{\bm{k}}\\\ \\\
\displaystyle\frac{\tilde{\pi}_{\bm{k}}}{k}\end{array}\right)=\left(\begin{array}[]{cc}\cos
k\Delta&-\sin k\Delta\\\ \\\ \sin k\Delta&\cos
k\Delta\end{array}\right)\left(\begin{array}[]{c}\mu_{\bm{k}}\\\ \\\
\displaystyle\frac{\pi_{\bm{k}}}{k}\end{array}\right).$ (20)

It is important to note that the above are classical relations between
canonical variables belonging to distinct canonical frameworks based on
distinct internal clocks. Although they are canonically inequivalent, these
two frameworks generate the same physical dynamics of the system, which is
required by the principle of relativity.

##  IV Quantization

Having completed the classical treatment of our system, we now move to the
investigation of the possible differences between the respective quantum
dynamics obtained from the quantization of these two different frameworks.

###  IV.1 Semi-classical background

Since, by definition, the scale factor is positive definite ($a>0$), one needs
to quantize our previous system on the half line. Although the position
operator $\widehat{Q}=a$ is self-adjoint on the half line, this is not the
case for the momentum operator $\widehat{P}=i\hbar\partial_{a}$, so we use
instead the symmetric dilation operator

$\widehat{D}=\\{\widehat{P},\widehat{Q}\\}=\frac{1}{2}\left(\widehat{P}\widehat{Q}+\widehat{Q}\widehat{P}\right)=\frac{1}{2}i\hbar\left(a\partial_{a}+\partial_{a}a\right).$

Classically the dilation variable is $d=ap_{a}$, so that the Hamiltonian,
expressed in term of $d$, is
$H^{(0)}=\frac{1}{2}p_{a}^{2}=\frac{1}{2}d^{2}/a^{2}$, and one can define its
quantum counterpart as a symmetric ordering of
$\frac{1}{2}\widehat{Q}^{-2}\widehat{D}^{2}$. Expanding on the basis
($\widehat{Q},\widehat{P}$), this yields

$\widehat{H}^{(0)}=-\frac{1}{2}\frac{\partial^{2}}{\partial
a^{2}}+\frac{\hbar^{2}K}{a^{2}},$

where the value of $K>0$ depends on the ordering; fixing one ordering such
that $K>\frac{3}{4}$ ensures $\widehat{H}^{(0)}$ is self-adjoint on the half
line Vilenkin (1988).

We can find some approximate solutions to the Schrödinger equation with a
family of coherent states (see, e.g., Refs. Klauder (2015, 2012) and Martin et
al. (2022) for the specific case under study here). We choose the coherent
states to read

$|a(\eta),p_{a}(\eta)\rangle=\mathrm{e}^{ip_{a}(\eta)\widehat{Q}/\hbar}\mathrm{e}^{-i\ln[a(\eta)]\widehat{D}/\hbar}|\xi\rangle,$
(21)

where $|\xi\rangle$ is such that the expectation values of $\widehat{Q}$ and
$\widehat{P}$ in $|a(\eta),p_{a}(\eta)\rangle$ are respectively $a(\eta)$ and
$p_{a}(\eta)$, and otherwise arbitrary (see however Ref. Bergeron et al.
(2023)).

The dynamics confined to the coherent states can be deduced from the quantum
action

$\mathcal{S}_{\textsc{q}}=\int\left\\{a^{\prime}(\eta)p_{a}(\eta)-H_{\text{sem}}\left[a(\eta),p_{a}(\eta)\right]\right\\}\mathrm{d}\eta,$
(22)

with the semiclassical Hamiltonian given by

$H_{\text{sem}}=\langle a,p_{a}|\widehat{H}^{(0)}|a,p_{a}\rangle,$ (23)

from which one derives the ordinary Hamilton equations

$a^{\prime}=\frac{\partial{H}_{\text{sem}}}{\partial
p_{a}}\quad\text{and}\quad
p^{\prime}_{a}=-\frac{\partial{H}_{\text{sem}}}{\partial a}.$ (24)

We find that the semiclassical background Hamiltonian reads Martin et al.
(2022)

$H_{\text{sem}}=\frac{1}{2}\left(p_{a}^{2}+\frac{\hbar^{2}\mathfrak{K}}{a^{2}}\right),$
(25)

where the new constant $\mathfrak{K}$ is positive ($\mathfrak{K}>0$). Its
specific value is related with both $K$ and the fiducial state $|\xi\rangle$.
We find the solution to (24) to read
$a^{2}(\eta)=a_{0}+a_{1}\eta+a^{2}(\eta)$, with
$a_{0}a_{2}-a_{1}^{2}/4=\hbar^{2}\mathfrak{K}>0$, so that the equation
$a(\eta)=0$ has no longer any real solution; the singularity is indeed quantum
mechanically avoided. Choosing the origin of time such that $a^{\prime}=0$ for
$\eta=0$ permits to rewrite this solution in full generality as

$\displaystyle a(\eta)$
$\displaystyle=a_{\textsc{b}}\sqrt{1+(\omega\eta)^{2}},$ (26a) $\displaystyle
p_{a}(\eta)$
$\displaystyle=\frac{a_{\textsc{b}}\omega^{2}\eta}{\sqrt{1+(\omega\eta)^{2}}},$
(26b)

where $a_{\textsc{b}}^{4}\omega^{2}=\hbar^{2}\mathfrak{K}$, which in turn
implies
$H_{\text{sem}}=\frac{1}{2}a_{\textsc{b}}^{2}\omega=\frac{1}{2}\hbar\sqrt{\mathfrak{K}}$;
it is clear that the model contains one and only one free parameter, namely
$\mathfrak{K}$. From now on, we assume that the background evolution is given
by Eqs. (26): this means the semi-classical potential

$V_{\text{sem}}=\frac{a^{\prime\prime}}{a}=\frac{\hbar^{2}\mathfrak{K}}{a^{4}}=\left[\frac{\omega}{1+\left(\omega\eta\right)^{2}}\right]^{2},$
(27)

shown in Fig. 1, never vanishes except in the large scale factor limit ($a\gg
1\Longrightarrow\eta\gg\omega^{-1}$). This is appropriate as this is also the
classical limit for which $a^{\prime\prime}\to 0$. A classical radiation-
dominated universe begins or ends with a singularity and produces no
gravitational waves, whereas our quantum radiation-dominated universe
naturally connects the contracting and expanding phases through a bounce which
is subsequently responsible for a non-vacuum spectrum of tensor perturbations,
to which we now turn.

Figure 1: The semi-classical potential $V_{\text{sem}}$ given by Eq. (27) as a
function of the conformal time $\eta$ for various values of the inverse bounce
duration $\omega$. The potential has to be compared with the relevant value of
$k^{2}$ ($k=0.01)$, indicated as a straight line. The corresponding scale
factor time evolution is shown in the insert.

###  IV.2 Quantum perturbations

For a given mode $\bm{k}$, the Hamiltonian $H^{(2)}_{\bm{k}}$, given by Eq.
(10), is easily quantized using the usual prescriptions. We assume that the
background follows the semi-classical approximation described above, so that
the potential for the perturbation is given by $V_{\text{sem}}$ [Eq. (27)].
The basic variables are replaced by a set of operators

$\begin{gathered}\mu_{\bm{k}}\mapsto\widehat{\mu}_{\bm{k}}=\sqrt{\frac{\hbar}{2}}\left[\widehat{a}_{\bm{k}}\mu^{*}_{k}(\eta)+\widehat{a}^{\dagger}_{-\bm{k}}\mu_{k}(\eta)\right],\\\
\pi_{\bm{k}}\mapsto\widehat{\pi}_{\bm{k}}=\sqrt{\frac{\hbar}{2}}\left[\widehat{a}_{\bm{k}}\mu^{*\prime}_{k}(\eta)+\widehat{a}^{\dagger}_{-\bm{k}}\mu^{\prime}_{k}(\eta)\right],\end{gathered}$
(28)

where we assume the Wronskian normalization condition
$\mu^{\prime}_{k}\mu^{*}_{k}-\mu_{k}\mu^{*\prime}_{k}=2i$ for the complex mode
functions $\mu_{k}$. The creation $\widehat{a}^{\dagger}_{\bm{k}}$ and
annihilation $\widehat{a}_{\bm{k}}$ operators satisfy the commutation
relations
$\left[\widehat{a}_{\bm{k}},\widehat{a}^{\dagger}_{\bm{p}}\right]=\delta_{\bm{k},\bm{p}}$
stemming from the canonical ones between the field operators
$\left[\widehat{\mu}_{\bm{k}},\widehat{\pi}_{-\bm{p}}\right]=i\hbar\delta_{\bm{k},\bm{p}}$.

Figure 2: The new time $\tilde{\eta}$ as a function of the original one $\eta$
for three different shapes of delay functions $\Delta_{1}$, $\Delta_{2}$ and
$\Delta_{3}$ defined through Eq. (31) along the original fixed bouncing
trajectory (26). The parameters are chosen as $A=B=D=1$, $C=4$ and $E=2$ for
$\Delta_{1}$, while we set $A=2$, $B=0.2$, $C=0.5$, $D=3$ and $E=4$ for
$\Delta_{2}$, and finally the set $A=-1$, $B=C=1$, $D=0.5$ and $E=3$ defines
$\Delta_{3}$.

In the Heisenberg picture, the equations of motion take the form

$i\hbar\frac{\mathrm{d}\widehat{\mu}_{\bm{k}}}{\mathrm{d}\eta}=\left[H^{(2)}_{\bm{k}},\widehat{\mu}_{\bm{k}}\right]\quad\hbox{and}\quad
i\hbar\frac{\mathrm{d}\widehat{\pi}_{\bm{k}}}{\mathrm{d}\eta}=\left[H^{(2)}_{\bm{k}},\widehat{\pi}_{\bm{k}}\right],$

which imply that the mode function $\mu_{k}(\eta)$ satisfies

$\frac{\mathrm{d}^{2}\mu_{k}}{\mathrm{d}\eta^{2}}+\left(k^{2}-\frac{\hbar^{2}\mathfrak{K}}{a^{4}}\right)\mu_{k}=0,$
(29)

where $a(\eta)$ is given by the semi-classical solution (26a). Using (27),
this transforms into

$\frac{\mathrm{d}^{2}\mu_{k}}{\mathrm{d}\eta^{2}}+\left\\{k^{2}-\left[\frac{\omega}{1+\left(\omega\eta\right)^{2}}\right]^{2}\right\\}\mu_{k}=0,$
(30)

which can be integrated numerically if initial conditions are provided: we
assume that far in the contracting branch, with $\eta_{\text{ini}}<0$ and
$V_{\text{sem}}(\eta_{\text{ini}})\ll k^{2}$, there was no gravitational wave,
so the field was in a vacuum state. This implies the mode function satisfies
$\mu_{k}(\eta_{\text{ini}})=\mathrm{e}^{-ik\eta_{\text{ini}}}/\sqrt{2k}$ and
$\mu^{\prime}_{k}(\eta_{\text{ini}})=-i\sqrt{k/2}\,\mathrm{e}^{-ik\eta_{\text{ini}}}$.

##  V Quantum ‘clocks’

In what follows, we study the effect of clocks on the quantum and semi-
classical dynamics of selected dynamical variables. First, we obtain the
dynamical trajectories in the reduced phase space $(a,{p},{\mu}_{k},\pi_{k})$
that is associated with the initial clock ${\eta}$; note that from that point
on, since there is no risk of confusion, we shall replace what was previously
denoted as $p_{a}$ simply by $p$. Next, we choose a set of delay functions
$\Delta(a,p)$ to define new clocks $\tilde{\eta}$ and obtain the new reduced
phase spaces $(\tilde{a},\tilde{p},\tilde{\mu}_{k},\tilde{\pi}_{k})$
associated with the new clocks. Then, we make use of Eqs (14), (15) and (20)
to transport the dynamical trajectories to these new phase spaces. We assume
that the latter admit a unique physical interpretation, and so the
trajectories can be meaningfully compared in these new variables. In other
words, there are many clocks denoted by ${\eta}$ and only one denoted by
$\tilde{\eta}$. Note that for $\Delta=0$ the clocks ${\eta}$ and
$\tilde{\eta}$ coincide. For this case, we assume that ${\eta}$ and
$(a,{p},{\mu}_{k},\pi_{k})$ are the variables of Sec. II, which sets the
physical meaning of the phase space
$(\tilde{a},\tilde{p},\tilde{\mu}_{k},\tilde{\pi}_{k})$ and the clock
$\tilde{\eta}$.

Figure 3: Semi-classical trajectories obtained in different clocks and mapped
into the initial reduced phase space $(\tilde{a},\tilde{p})$ to compare with
the original trajectory represented by the full black line.

Figure 4: Evolution of the primordial gravity-wave
$\mathfrak{R}\mathrm{e}(\tilde{\mu})$ for two different wavenumbers, $k=0.1$
(top panel) and $k=0.5$ (bottom), and for different clocks based on the first
class of delay function, $\Delta_{1}$, $\Delta_{2}$ and $\Delta_{3}$,
represented by the dotted blue line, dashed red line and dashed-dotted green
line respectively. The original trajectory is represented by the full black
line. In Fig. A3 the same plot for the second class of delay function is
depicted to show how the choice of delay function affects the time of
convergence. Figure 5: Evolution of the primordial gravity-wave
$\mathfrak{R}\mathrm{e}(\tilde{\mu}_{k})$ plotted for 4 different wavenumber
$k$ values. For each fixed $k$ we changed the clock considering the family of
delay functions $\Delta$, whose value is the same as in Fig. 3.

###  V.1 Clock choices and background

In order to illustrate the clock choice issue, we consider a family of delay
functions, namely

$\Delta(a,p)=A\frac{a^{B}}{(a+C)^{D}}\frac{\sin(Ep)}{p},$ (31)

where $A$, $B$, $C$, $D$ and $E$, are arbitrary coefficients, whose values are
limited to ensure that the conditions presented in Sec. III.3 hold. In the
Appendix, we consider another set of acceptable delay functions to show that
our conclusions are not restricted to the choice (31).

A few ‘clocks’ corresponding to the delay function $\Delta(a,p)$ are
represented along a semi-classical dynamical trajectory for different choices
of the free parameters in Fig. 2. It shows that, contrary to the classical
case where the condition (16) holds, the new clocks, in general, are no longer
monotonic due to quantum corrections.

Applying the clock transformation of Fig. 2 to the background solution (26)
yields Fig. 3 once mapped into the reduced phase space, with the original
trajectory superimposed for comparison. All the trajectories originate in the
same classical regime at large $\tilde{a}$ and negative $\tilde{p}$, i.e. at a
time at which the universe is large and contracting. Close to the
$\tilde{a}=0$ boundary, where the quantum behavior dominates, they all somehow
bounce in the variables $\tilde{a}$ and $\tilde{p}$, diverging from one
another and providing different accounts of the bounce. Finally, they reach
the region of large $\tilde{a}$ and positive $\tilde{p}$ where they converge
again to the unique classical behavior representing a large and expanding
universe.

Possible differences between the trajectories include the values of
$\tilde{a}$ and $\tilde{p}$ at which the bounce occurs, the level of asymmetry
between contracting and expanding branches, or even the number of bounces.
These semi-classical trajectories illustrate the non-unitary relation between
different clocks. Nevertheless, they all originate from a unique contracting
classical universe and end towards a similarly unique expanding classical
universe. Therefore, the semi-classical trajectories in different clocks yield
the same outcome for large and classical universes. Notice that the
trajectories’ convergence before and after the bounce can be delayed as much
as one wants by making use of appropriate delay functions, such as that
discussed in the Appendix, i.e. Eq. (32), whose effects on both background and
perturbation trajectories can be seen in the Appendix.

Let us now move to the perturbation of these homogeneous solutions and compare
the different evolution that can result from using different clocks.

###  V.2 Clocks and perturbations

In Fig. 4, we plot the dynamics of the real part of the perturbation variable
$\tilde{\mu}_{k}$ against the delayed time $\tilde{\eta}$ for the three
different functions of Eq. (31) displayed in Fig. 2, and for two values of the
comoving wavenumber $k$. The figure illustrates our general finding that the
absolute clock effect is more or less equally strong and lasting roughly
equally long for all-wavelength perturbations. This is shown more convincingly
in Fig. 5 in which the evolution of 4 different modes is shown as a close-up
in the quantum-dominated bouncing region. This means that the larger the
wavelength of the perturbation, the larger the relative clock effect, and the
longer it lasts in units of its oscillation period. Thus, the clock effect is
more important for phenomena occurring at small timescales and over short
distances. Moreover, the evolving amplitude $\tilde{\mu}_{k}$, in general, is
not a function of the clock $\tilde{\eta}$ due to quantum effects that disrupt
the monotonicity relation between quantized clocks.

Given that both the background and the perturbation modes evolve in such a way
as to reach a unique configuration, the primordial gravity-wave amplitude
$\tilde{\mu}_{k}/\tilde{a}$, which is the quantity one expects to measure in
practice Micheli and Peter (2022), also converges to a unique solution, making
the model predictive.

All the plots above illustrate the non-unitary relation between different
clocks, as well as the spoiling of the clock monotonicity at the quantum
level, which is illustrated in Fig. 2. Nevertheless, similarly to the semi-
classical background trajectories, the perturbation variable
$\mathfrak{R}\mathrm{e}(\tilde{\mu})$ visibly converges to a unique classical
solutions from a well-defined asymptotic past initial condition to the
asymptotic future. Therefore, one can safely extend the background conclusion
to the perturbations: the time development of the mode
$\mathfrak{R}\mathrm{e}(\tilde{\mu})$ using different clocks yields the same
predictions in the large and classical universe regime. The delay of the
convergence due to different choices of delay functions can be seen in Fig. 5.

As a final illustration of the perturbation behavior through the quantum
bounce, we find it useful to inspect the phase space trajectories in the plane
$[\mathfrak{R}\mathrm{e}(\tilde{\mu}_{k}),\mathfrak{I}\mathrm{m}(\tilde{\mu}_{k})]$
as is displayed in Fig. 6. The initial vacuum state is represented by a circle
that is squeezed into an ellipse during the contraction and bounce, squeezing
that represents the amplification of the amplitude of the perturbation. From
the point of view of the time problem, the initial circle and the final
ellipse respectively represent the asymptotic past and future of the
amplitude: from the point of view of physical prediction, the indeterminacy
occurring near the bounce as may develop through various different times
disappears in the asymptotic regimes, so that the existence of a classical
approximation in our trajectory approach ensures the standard procedure of
treating the perturbations leads to physically meaningfull predictions.

Figure 6: Evolution of the real versus imaginary part of $\tilde{\mu}_{k}$ for
a wavenumber $k=0.5$ and a bounce parameter $\omega=1$. The initial circle
represents the initial vacuum state of the perturbation, while the ellipse
shows the final squeezed state, which happens, in the case at hand, to have a
slight phase shift with respect to the real axis. The transition between these
two asymptotic cases differs for the different delay functions $\Delta_{1}$,
$\Delta_{2}$ and $\Delta_{3}$, whose trajectories are represented by the
dotted blue line, dashed red line and dashed-dotted green line respectively,
the original trajectory being represented by the full black line.

##  VI Discussion and perspectives

In this work, we explored the time problem in the framework of quantum fields
on quantum spacetimes. We considered the specific example of primordial
gravitational waves propagating through a bouncing quantum Friedmann universe.
We pointed to several features that we believe to be universal for such
models.

First, we showed that the dynamical variables, such as the scale factor or the
amplitude of a gravitational wave, obtained from different internal clocks,
evolve differently when compared in a clock-independent manner. Second, these
expectation values (background evolution) and mode functions of operators
(perturbations), irrespective of the clock chosen, converge to a unique
evolution for large classically-behaving universes. This is the phase space
domain in which unambiguous predictions can be made. Third, for different
clocks, the dynamics converges to the classical behavior at different times.
In principle, there is no restriction on how far from the bounce the system
must be in order to display the classical behavior. In practice, however, all
the clocks considered were found to converge very quickly, allowing for
unambiguous predictions shortly after the bounce.

Based on the above findings, we postulate that the physical predictions are
only those predictions provided by any clock, which are not altered upon the
clock’s transformation. The fact that for large universes the semi-classical
background dynamics and the quantum perturbation dynamics do not depend on the
clock implies the following: Despite the fact that the dynamical variables are
not Dirac observables, they provide physical predictions for large universes,
which is precisely the regime in which we observe the actual Universe.

Note, however, that the word large is never precisely defined. One could
expect that, at least in principle, some clocks require times larger than the
present age of the Universe to converge to the classical behavior. This,
however, poses no problem to our interpretation as we simply exclude such
clocks and retain only those that behave classically in the domain for which
we make predictions. This may seem arbitrary and unjustified. We must,
however, remember that, as a matter of fact, any semi-classical description of
ordinary quantum mechanics is necessarily restricted to a limited set of
observables, usually the simple ones, while more compound observables often
display classically incompatible behavior (e.g. $\langle
x\rangle^{2}\not=\langle x^{2}\rangle$). For similar reasons, we are allowed
to choose only those clocks in which the dynamics of the relevant observables
is classically consistent.

On the one hand, we proved that the evolution of the expectation values of
some observables constitute physical predictions of quantum cosmological
models. On the other hand, the expectation values are not all that is measured
in the large Universe. In other words, not all objects are classical in the
large Universe. For instance, the position of an electron is a dynamical
variable that can be measured in a lab. So, could the outcomes of such a
measurement also be unambiguously predicted by a quantum cosmological model?
The answer is affirmative. Note that the mode function $\mu_{k}$, whose
dynamics becomes unambiguous in large universe, determines the evolution of
the operator $\widehat{\mu}_{k}$ via Eq. (28). This implies that the
Heisenberg equation of motion encoded in Eq. (30) becomes unambiguous too.
Obviously, the evolution of perturbation in the Schrödinger picture must
consequently become unique as well. Hence, ordinary quantum mechanics of
perturbation modes is recovered in a large universe. This conclusions must
also apply to electrons and, in general, to all non-gravitational degrees of
freedom.

To better understand the origin of the emergence of ordinary quantum
mechanics, notice that any clock transformation (14) involves, by definition,
only background variables. If the latter behave classically, the clock
transformation is completely classical and amounts to a mere (in general,
nonlinear) change of units of time. In Ref. Małkiewicz and Miroszewski (2017),
it was demonstrated that the relational dynamics of a quantum variable in a
classical clock is unambiguous in the sense that switching to another
classical clock does not induce any clock effect.

Let us put to test our approach and our result by addressing a set of
questions that were proposed in Ref. Isham (1993) for assessing the
completeness of any potential solution to the time problem.

1\. How should the notion of time be re-introduced into the quantum theory of
gravity?

Our approach relies on evolving internal variables called clocks. We express
the dynamics of the dynamical variables in terms of these clocks.

2\. In particular, should attempts to identify time be made at the classical
level, i.e., before quantization, or should the theory be quantized first?

In our approach, we first reduce the Hamiltonian formalism based on a selected
clock, next we quantize the reduced formalism as if the clock was an external
and absolute time. However, it is neither external nor absolute. The
instantaneous value of the clock determines the instantaneous physical state
of the system. Switching to another clock entails a change in the physical
interpretation of the clock and the entire state of the system.

3\. Can ‘time’ still be regarded as a fundamental concept in a quantum theory
of gravity, or is its status purely phenomenological? […]

Figure A1: Changes in the time variable $\eta$ for the second family of delay
functions $\Delta^{\prime}_{1}$, $\Delta^{\prime}_{2}$ and
$\Delta^{\prime}_{3}$ given by Eq. (32) along a fixed bouncing trajectory,
with parameters chosen such that $\Delta_{1}=ae^{p_{a}}$,
$\Delta_{2}=ae^{3p_{a}/2}$ and $\Delta_{3}=ae^{2p_{a}}$.  Figure A2: Semi-
classical trajectories mapped into the initial reduced phase space $(a,p)$ for
the second class of delay function Eq. (32), with the same parameters as in
Fig. A1.

Figure A3: Evolution of the real part of the primordial gravity-wave
$\mathfrak{R}\mathrm{e}(\tilde{\mu})$ for two different wavenumbers, $k=0.1$
and $k=0.5$, and for different clocks for the second class of delay functions,
$\Delta^{\prime}_{1}$, $\Delta^{\prime}_{2}$ and $\Delta^{\prime}_{3}$,
respectively represented by the dotted blue line, dashed red line and dashed-
dotted green line. The original trajectory is represented by the full black
line.

In our approach, there is no fundamental ‘time’. The fundamental concept is
‘change’ or ‘evolution’, meaning we merely need to assume that the 3+1 split
of the underlying geometry imposes an ordered set of hypersurfaces. As we
showed in this paper, extracting dynamical predictions from such a formalism
is a subtle issue. The clocks serve as tools for deriving the predictions.
Once a class of clocks converges to a unique dynamics, any one of them can be
treated analogously and deserve the qualification of ‘time’, and any quantum
dynamical variable becomes described in them by a unique Schrödinger equation.
This is how ordinary quantum mechanics emerges.

4\. If ‘time’ is only an approximate concept, how reliable is the rest of the
quantum-mechanical formalism in those regimes where the normal notion of time
is not applicable? In particular, how closely tied to the concept of time is
the idea of probability? […]

The quantum-mechanical description in the regime where different clocks
exhibit different dynamics is an essential part of our theory. It describes
the deterministic evolution of the system. However, this regime does not seem
to allow for any meaningful dynamical interpretation in terms of relational
change. Although we have not explicitly addressed this question in the present
work, our approach permits to do it.

To conclude, one can mention that the chosen ’clock’ degrees of freedom,
although perfectly acceptable as such in the classical framework of general
relativity, are arguably not in the quantum regime. They do not qualify as
actual clocks since, along the quantum trajectory, they yield a non monotonic
change of time variable, in other words, they provide different hypersurface
orderings. This might be cured by adding to the classical clock transformation
(14) a quantum term that need be identified. One may also argue that we are
insisting upon using a trajectory to define the background evolution, while
some might insist upon the fact that there is no such thing as a trajectory in
quantum mechanics.

In any case, it is interesting to note that whichever of the possibilities
above happens to be valid, the critical point that is made here is that even
though the quantum-dominated phase is indeed ill-defined both at the
background and perturbation levels from the point of view of time development,
the asymptotic regimes end up being unique. As a result, setting well-
motivated initial conditions in the classical past, one gets unambiguous
physical predictions for the classical future in which we happen to perform
the ensuing measurements. In other words, we have shown that the lack of
predictability in the quantum regime does not exclude the fact that the theory
permits meaningful physical predictions that can be tested with observations.

###### Acknowledgements.

A.B. and P.M. acknowledge the support of the National Science Centre (NCN,
Poland) under the research Grant 2018/30/E/ST2/00370.

Figure A4: Evolution of the primordial gravitational amplitude for different
clocks obtained from the second class of delay function 32. Convergence
happens at a latter time with respect to the first class of delay functions
(31), as can be seen by comparison with Fig. 5.

*

## Appendix A

In this appendix, we consider the alternative choice of a family of two-
parameter delay functions, namely

$\Delta^{\prime}(a,p)=a^{A}e^{Bp},$ (32)

which define a new set of clocks plotted in Fig. A1 for a few relevant values
of the parameters $A$ and $B$. Fig. A2 depicts the trajectories with different
clocks obtained from $\Delta^{\prime}(a,p)$ for which the convergence happens
much later than in the case discussed in the core of this paper, as can be
seen by comparing with Fig. 3. The extent to which this delay can be
increased, and how the matter content of the Universe can affect this limit,
is not dealt with in the present article and will be the subject of a future
work.

One can note that the delay functions (A1) tend to diverge in time from one
another, all of them growing exponentially with the momentum, the phase space
trajectories however do converge to the undelayed one, but at scales that are
increasingly larger with the amplitude of the exponential behavior of the
relevant delay function.

Moving to the perturbations, we performed the same analysis as in the core of
this paper and show the time development of the real part of the mode function
for different values of the wavenumber in Fig. A3, with a special emphasis at
the near-bounce regime in Fig. A4. As for the other family of delay functions,
we find that whenever the classical approximation for the background holds,
one recovers a unique prediction.

## References

  * Kuchar (2011) K. V. Kuchar, Int. J. Mod. Phys. D 20, 3 (2011).
  * Isham (1993) C. J. Isham, NATO Sci. Ser. C 409, 157 (1993), eprint gr-qc/9210011.
  * Anderson (2010) E. Anderson (2010), eprint 1009.2157.
  * Kiefer and Peter (2022) C. Kiefer and P. Peter, Universe 8, 36 (2022), eprint 2112.05788.
  * Deser et al. (1961) S. Deser, R. Arnowitt, and C. W. Misner, Nuovo Cim. 19, 668 (1961).
  * Beluardi and Ferraro (1995) S. C. Beluardi and R. Ferraro, Phys. Rev. D 52, 1963 (1995), eprint gr-qc/9409032.
  * Hajicek (2000) P. Hajicek, Nucl. Phys. B Proc. Suppl. 80, 12/13 (2000), eprint gr-qc/9903089.
  * Catren and Ferraro (2001) G. Catren and R. Ferraro, Phys. Rev. D 63, 023502 (2001), eprint gr-qc/0006027.
  * Hajicek and Kijowski (2000) P. Hajicek and J. Kijowski, Phys. Rev. D 62, 044025 (2000).
  * Małkiewicz (2015) P. Małkiewicz, Class. Quant. Grav. 32, 135004 (2015), eprint 1407.3457.
  * Małkiewicz and Miroszewski (2017) P. Małkiewicz and A. Miroszewski, Phys. Rev. D96, 046003 (2017), eprint 1706.00743.
  * Bojowald and Halnon (2018) M. Bojowald and T. Halnon, Physical Review D 98 (2018), ISSN 2470-0029, URL http://dx.doi.org/10.1103/PhysRevD.98.066001.
  * Peter et al. (2007) P. Peter, E. J. C. Pinho, and N. Pinto-Neto, Phys. Rev. D 75, 023516 (2007), eprint hep-th/0610205.
  * Pinto-Neto (2008) N. Pinto-Neto, Prog. Theor. Phys. Suppl. 172, 178 (2008).
  * Peter and Pinto-Neto (2008) P. Peter and N. Pinto-Neto, Phys. Rev. D 78, 063506 (2008), eprint 0809.2022.
  * Malkiewicz (2017) P. Malkiewicz, Class. Quant. Grav. 34, 145012 (2017), eprint 1601.04857.
  * Małkiewicz et al. (2020) P. Małkiewicz, P. Peter, and S. D. P. Vitenti, Phys. Rev. D 101, 046012 (2020), eprint 1911.09892.
  * Małkiewicz et al. (2022) P. Małkiewicz, P. Peter, and S. D. P. Vitenti, Universe 8, 71 (2022), eprint 2201.09138.
  * Rovelli (2002) C. Rovelli, Phys. Rev. D 65, 124013 (2002), eprint gr-qc/0110035.
  * Dittrich (2007) B. Dittrich, Gen. Rel. Grav. 39, 1891 (2007), eprint gr-qc/0411013.
  * Bojowald and Tsobanjan (2023) M. Bojowald and A. Tsobanjan, Physical Review D 107 (2023), ISSN 2470-0029, URL http://dx.doi.org/10.1103/PhysRevD.107.024003.
  * Małkiewicz (2017) P. Małkiewicz, Class. Quant. Grav. 34, 205001 (2017), eprint 1505.04730.
  * Dirac (1964) P. Dirac, _Lectures on Quantum Mechanics_ , Belfer Graduate School of Science. Monographs series (Belfer Graduate School of Science, Yeshiva University, 1964), URL https://books.google.fr/books?id=oQNRAAAAMAAJ.
  * Peter and Uzan (2013) P. Peter and J.-P. Uzan, _Primordial Cosmology_ , Oxford Graduate Texts (Oxford University Press, 2013), ISBN 978-0-19-966515-0, 978-0-19-920991-0.
  * Micheli and Peter (2022) A. Micheli and P. Peter (2022), eprint 2211.00182.
  * Vilenkin (1988) A. Vilenkin, Phys. Rev. D 37, 888 (1988).
  * Klauder (2015) J. R. Klauder, _Enhanced quantization: Particles, fields and gravity_ (World Scientific, Hackensack, 2015), ISBN 978-981-4644-62-4, 978-981-4644-64-8.
  * Klauder (2012) J. R. Klauder, J. Phys. A 45, 285304 (2012), eprint 1204.2870.
  * Martin et al. (2022) J. d. C. Martin, P. Małkiewicz, and P. Peter, Phys. Rev. D 105, 023522 (2022), eprint 2111.02963.
  * Bergeron et al. (2023) H. Bergeron, J.-P. Gazeau, P. Małkiewicz, and P. Peter (2023), eprint 2310.16868.


# Massively Parallel Open Modification Spectral Library Searching with
Hyperdimensional Computing

Jaeyoung Kang, Weihong Xu, Wout Bittremieux, Tajana Rosing  University of
California, San Diego j5kang, wexu, wbittremieux<EMAIL_ADDRESS>
###### Abstract.

Mass spectrometry, commonly used for protein identification, generates a
massive number of spectra that need to be matched against a large database. In
reality, most of them remain unidentified or mismatched due to unexpected
post-translational modifications. Open modification search (OMS) has been
proposed as a strategy to improve the identification rate by considering every
possible change in spectra, but it expands the search space exponentially. In
this work, we propose HyperOMS, which redesigns OMS based on hyperdimensional
computing to cope with such challenges. Unlike existing algorithms that
represent spectral data with floating point numbers, HyperOMS encodes them
with high dimensional binary vectors and performs the efficient OMS in high-
dimensional space. With the massive parallelism and simple boolean operations,
HyperOMS can be efficiently handled on parallel computing platforms.
Experimental results show that HyperOMS on GPU is up to $17\times$ faster and
$6.4\times$ more energy efficient than the state-of-the-art GPU-based OMS tool
(Bittremieux et al., 2019) while providing comparable search quality to
competing search tools.

Spectral library search, Mass spectrometry, Proteomics, Hyperdimensional
computing

††ccs: Computing methodologies Massively parallel algorithms††ccs: Applied
computing Computational proteomics

##  1\. Introduction

Proteomics plays an essential role in understanding the molecular mechanisms
of proteins, which are responsible for various tasks in a life of a cell
(Aebersold et al., 2016). Mass spectrometry (MS) is one of the most popular
and reliable approaches to identifying and quantifying proteins and peptides
in complex biological samples (Smith et al., 2014). In a typical MS
experiment, a method called tandem mass spectrometry (MS/MS) generates a
massive amount of MS/MS spectra data. Then researchers determine peptide
annotations of the MS/MS spectra via spectral library searching (Lam, 2011).
Peptide sequences are assigned to experimental MS/MS spectra by matching them
against a spectral library of known peptides (see Figure 1).

The challenge of spectral library searching is that a significant portion of
MS/MS spectra acquired during an experiment cannot be directly identified
(Chick et al., 2015) using searching approaches with conventional
similarities, such as cosine similarity. In reality, proteins undergo one or
more post-translational modifications (PTMs), which change their mass and
MS/MS fragmentation pattern. PTMs can be introduced during sample preparation
as an artifact of MS measurement (Bittremieux et al., 2018b), or biologically
relevant PTMs arise _in vivo_ (Walsh et al., 2005). However, spectral
libraries mainly contain reference spectra for unmodified peptides, so PTMs
make experimental spectra challenging to identify as they no longer exactly
match the reference spectra.

Figure 1. Overview of spectral library searching. Standard searching uses a
narrow precursor $m$/$z$ tolerance, while OMS uses a wide precursor $m$/$z$
tolerance during the searching.

Open modification searching (OMS) is a promising approach to circumvent this
limitation and identify modified spectra (Na et al., 2015). Standard spectral
library searching only compares experimental spectra to reference spectra with
a similar precursor mass, i.e., the mass of the unfragmented peptide, as
matching peptides should have an identical mass. In contrast, OMS performs
spectra matching on a wider range of reference spectra, where OMS compares
modified query spectra to their unmodified reference variants, even when their
precursor mass differs due to PTMs. OMS’ wider searching range empowers higher
identification capability, thus enabling the study of more complex protein
interaction in virus-host (Adams et al., 2022) and proteomics analysis of non-
model organisms (Heck et al., 2020).

OMS faces several challenges in terms of low searching speed, efficiency, and
accuracy. Compared to standard searching, the OMS speed is extremely slow due
to the drastically increased search space since it considers all possible
reference matches for each query spectrum (Bittremieux et al., 2018a). This
problem is further exacerbated by the increasing spectral data due to the cost
reduction in the MS experiment (2$\times$ in recent two years) (Aksenov et
al., 2017; Perez-Riverol et al., 2022). Also, large spectral libraries created
by repository-scale mining of open MS data become available (Griss et al.,
2013; Wang et al., 2018). For example, the size of human HCD (higher energy
collisional dissociation) spectral libraries (mas, 2022a; Wang et al., 2018)
hold $2.15$ million data points, which is $4\times$ larger than the previous
NIST-HCD (Wang et al., 2018). MassIVE contains $5.6$ billion spectra, which
corresponds to $448$TB in size (mas, 2022b).

Several tools have been introduced to efficiently perform OMS (Kong et al.,
2017; Bittremieux et al., 2018a; Chi et al., 2018; Devabhaktuni et al., 2019;
Lam et al., 2007; Ma et al., 2014; Bittremieux et al., 2019). These tools use
various techniques to refine the search space, such as fragment ion indexing
(Kong et al., 2017), nearest neighbor searching (Bittremieux et al., 2018a,
2019), or tag-based filtering (Chi et al., 2018; Devabhaktuni et al., 2019).
For example, the state-of-the-art OMS tool ANN-SoLo performs nearest neighbor
searching using GPU and computes shifted cosine similarities on candidates
(Bittremieux et al., 2019). However, existing solutions involve a complex
execution pipeline and exhibit low data parallelism requiring high-precision
floating-point (FP32) arithmetic for good search quality, e.g., shifted cosine
similarity (Bittremieux et al., 2018a). As such, we redesign an OMS algorithm
that only involves hardware-friendly Boolean operations with a simple
execution pipeline.

In this work, we propose novel hyperdimensional computing (HDC)-inspired OMS
algorithm called HyperOMS. The spectral matching algorithm of HyperOMS is
based on the efficient computing paradigm, HDC (Kanerva, 2009), that has shown
high efficiency for pattern matching tasks, especially on parallel computing
platforms, e.g., GPU (Kang et al., 2022b, a). HDC improves the data
separability and robustness (Poduval et al., 2021; Zhang et al., 2021) by
mapping data into high-dimensional (HD) space, where the information is
distributed to every dimension of the HD vector. We leverage HDC’s robustness
to minimize the effects of PTMs. In particular, our method reflects the
spatial and value locality of peaks in the spectrum, making the encoded data
resilient to peak shifts and intensity changes. HyperOMS addresses the search
space challenges in OMS by approximating possible MS peak changes; spectra can
be identified with a single similarity computation. Furthermore, the proposed
HDC-empowered OMS algorithm is hardware-friendly as it replaces FP32
operations with simple Boolean arithmetic by encoding data into binary HVs,
leading to better computational efficiency. We implement the HyperOMS
algorithm on the high-performance HDC library (Kang et al., 2022b, a) on GPU,
obtaining up to $17\times$ speedup and $6.4\times$ energy efficiency over the
state-of-the-art OMS solution, ANN-SoLo (Bittremieux et al., 2019) while
offering comparable search quality.

##  2\. Background and Related Work

###  2.1. Mass Spectrometry Background

MS is used to study the biological process in proteomics via the analysis of
protein expression or state in cells or tissue (Smith et al., 2014). Proteins
are ubiquitous building blocks of life, and they are composed of peptides,
which are chains of amino acids, which can be described as a string of
letters.

During MS data acquisition, peptides are ionized to receive a charge, and
their mass-over-charge ($m$/$z$) is measured. The first intact ions are
measured in an MS scan using data-dependent acquisition, and the resulting MS
spectrum contains the corresponding $m$/$z$ values. The most intense peaks in
the MS spectrum are selected. It is further analyzed in MS/MS scans, i.e., the
second mass spectrometer. Ions with matching $m$/$z$ are isolated and
fragmented to generate MS/MS spectra. Fragmentation occurs along the peptide
backbone in between its constituent amino acids. Peptides are split into their
possible amino acid subsequences. We record the $m$/$z$ and intensity values
of all fragments, and the measured spectrum forms a unique fingerprint of the
measured peptide. Thus, each MS/MS spectrum consists of fragment $m$/$z$
values, spectrum charge, and its intact $m$/$z$ from the preceding MS scan,
called the precursor $m$/$z$ (Figure
1-$\mathbin{\vphantom{\circledast}\text{\ooalign{\smash{\raisebox{-2.58334pt}{\scalebox{2.3}{$\bullet$}}}\cr\smash{{\color[rgb]{1,1,1}\bf\footnotesize$1$}}\cr}}}$).

Spectral library searching determines which peptide corresponds to the
measured spectra (Lam, 2011). (Figure
1-$\mathbin{\vphantom{\circledast}\text{\ooalign{\smash{\raisebox{-2.58334pt}{\scalebox{2.3}{$\bullet$}}}\cr\smash{{\color[rgb]{1,1,1}\bf\footnotesize$2$}}\cr}}}$).
A spectral library contains reference spectra, each with known peptide labels.
We first select the reference candidates with a similar precursor $m$/$z$ to a
query spectrum. Next, similarities between the query spectrum and all
candidates are computed. Finally, the query spectrum is assigned the same
peptide label as its highest-scoring reference match. Here, we apply a false
discovery rate (FDR) filter on search results (Elias et al., 2007) (Figure
1-$\mathbin{\vphantom{\circledast}\text{\ooalign{\smash{\raisebox{-2.58334pt}{\scalebox{2.3}{$\bullet$}}}\cr\smash{{\color[rgb]{1,1,1}\bf\footnotesize$3$}}\cr}}}$),
which is a popular strategy called target–decoy strategy (Elias et al., 2007)
in MS/MS analysis to reduce false positives. Decoy spectra that cannot exist
are added to the spectral library besides the real (target) spectra. Decoy
spectra selected by the searching tool are filtered out. The number of target
SSMs and decoy SSMs at a specific score can be used to compute the FDR.
Typically, an FDR threshold of 1% is used to minimize the number of incorrect
identifications. The performance of different search tools can be compared by
the number of identified spectra at a fixed FDR threshold.

A standard searching strategy can identify directly matching spectra. It
assumes that precursor $m$/$z$ of query and matched reference spectra are
similar (narrow precursor $m$/$z$ tolerance). However, as spectral libraries
mainly contain unmodified reference spectra, they cannot be used to identify
modified ones. Modified ones have a different intact mass, as the
modifications induce mass shifts. Open modification searching (OMS) addresses
these issues by (i) using a wide precursor $m$/$z$ tolerance that exceeds mass
shifts induced by modifications to select reference candidates (Na et al.,
2015), and (ii) using alternative spectrum similarity measures that take peak
shifts due to modifications into account (Bittremieux et al., 2018a). Using a
wide precursor $m$/$z$ tolerance enables finding (partial) matches between
unmodified reference spectra and their modified variants. However, a large
number of candidates need to be evaluated for each query spectrum, which can
be computationally demanding.

###  2.2. Accelerated Spectral Library Searching

OMS has recently become an increasingly popular search strategy, and there
have been several studies to accelerate searches on different hardware
platforms. On CPU, MSFragger (Kong et al., 2017) uses a fragment ion index to
retrieve matching reference spectra. ANN-SoLo (Bittremieux et al., 2018a,
2019) is a state-of-the-art OMS tool that uses GPU-powered nearest neighbor
searching library, Faiss (Johnson et al., 2017), to prune the search space.
Several studies have focused on accelerating spectral library searching using
GPUs for efficient spectrum–spectrum similarity computation (Baumgardner et
al., 2011; Li et al., 2014). (Yang et al., 2019) used a CPU-FPGA architecture
in which multiple FPGAs are used for scalability and parallelism. However,
none of these studies have addressed the OMS challenges.

Figure 2. Overview of OMS process using HyperOMS.

(a) Data preprocessing.

(b) Encoding.

(c) Hamming similarity searching.

Figure 3. Data preprocessing and stages of HyperOMS algorithm.

##  3\. HyperOMS Algorithm

HyperOMS encodes raw spectral data to a binary HD vector called hypervector
(HV) to capture the position and intensities of peaks while considering the
spatial locality and value locality. Although peaks are shifted due to PTMs,
the similarity between a query spectrum and a matching reference spectrum does
not change dramatically. Furthermore, since binary vector representation is
used, HyperOMS replaces complex similarity metrics in existing OMS tools with
a simple Hamming similarity operation.

###  3.1. Overall Flow

Figure 2 shows a flow of HyperOMS. It starts with a data preprocessing step, a
common step for OMS. It refines and vectorizes raw spectra and compresses
them, resulting in spectrum vectors. HyperOMS encodes spectrum vectors into HV
during the encoding step. Next, the Hamming similarity searching step filters
reference spectra according to the query’s precursor $m$/$z$ and spectrum
charge and computes Hamming similarity between query HVs and candidate
reference HVs.

###  3.2. Data Preprocessing

The preprocessing step (1) refines the raw spectra by removing redundant peaks
and (2) vectorizes refined spectra (see Figure 3(a)). First, raw spectra are
refined to gather meaningful peaks
($\mathbin{\vphantom{\circledast}\text{\ooalign{\smash{\raisebox{-2.58334pt}{\scalebox{2.3}{$\bullet$}}}\cr\smash{{\color[rgb]{1,1,1}\bf\footnotesize$1$}}\cr}}}$).
We remove peaks whose intensity is below 1% of the most intense peaks. Low-
intensity peaks are considered noise. In turn, we retain $50$ to $150$ most
intensive peaks of the spectra. Existing studies (Bittremieux et al., 2019;
Lam et al., 2007) have shown that we can effectively refine spectra in this
manner.

Next, we vectorize the filtered spectra
($\mathbin{\vphantom{\circledast}\text{\ooalign{\smash{\raisebox{-2.58334pt}{\scalebox{2.3}{$\bullet$}}}\cr\smash{{\color[rgb]{1,1,1}\bf\footnotesize$2$}}\cr}}}$).
The peaks are discretized by binning the $m$/$z$ range to represent a spectrum
into a sparse vector of floating-point intensity values, called a spectrum
vector. If multiple peaks are assigned to the same $m$/$z$ bin, we sum their
intensity values. A large bin width can lead to loss of information when peaks
are grouped into a single bin. For example, the mass range between 0 $m$/$z$
and 2000 $m$/$z$ and bin width $0.04$ (based on the resolution of the mass
spectrometer) results in a dimensionality of $50,000$. The resulting spectrum
vectors have sparsity less than $1\%$. There are $50$ to $150$ peaks for each
spectrum vector, and its dimensionality is $20,000$ to $50,000$. As such, we
compress spectrum vectors in a compressed sparse row (CSR) format
($\mathbin{\vphantom{\circledast}\text{\ooalign{\smash{\raisebox{-2.58334pt}{\scalebox{2.3}{$\bullet$}}}\cr\smash{{\color[rgb]{1,1,1}\bf\footnotesize$3$}}\cr}}}$).
The preprocessing step is normally run offline, and the resulting data is
stored as a binary file for future use. In the following, we discuss the
HyperOMS algorithm, which first encodes spectrum vector to HV and performs
searching on HV.

###  3.3. Encoding: Spectrum Vectors to Hypervectors

HyperOMS encodes the data into a binary vector representation, which can
enhance the computation efficiency. There have been several efforts to
represent raw data in an HD binary vector, using Locality Sensitive Hashing
(LSH) (Koga et al., 2006; Norouzi et al., 2012; Imani et al., 2020) or HDC
(Imani et al., 2017; Kim et al., 2018; Salamat et al., 2019; Imani et al.,
2020). However, these strategies do not reflect the characteristics of OMS,
including peak shifts and intensity changes. For example, they treat each
feature position (corresponding to peak indices in spectrum vector) as
orthogonal. Peak shifts can lead to significant changes in similarity.
Conversely, our encoding takes both spatial locality (for peak shift) and
value locality (for peak intensity change from instrument error or noise) of
each feature into account. As a result, we can preserve similarity despite
peak changes, which is essential for OMS.

(a)

(b)

Figure 4. Position HV generation. (a) Strategy overview. (b) Pairwise
similarity (Hamming similarity normalized by the HV dimension size) between
position HVs.

Figure 3(b) shows the encoding process of HyperOMS. Unique position HVs
$\mathbf{F}$ are assigned for each index in a spectrum vector, i.e.,
$\mathbf{F}_{i}$ corresponds to index $i$, and
$\mathbf{F}\in\\{\mathbf{F}_{1},\mathbf{F}_{2},\ldots,\mathbf{F}_{f}\\}$. Note
that $f$ is the dimensionality of spectrum vector. Similarly, we use level HVs
$\mathbf{L}$ to capture different intensity values in each index. We quantize
intensity range into $Q$ levels, and $\mathbf{L}_{i}$ is assigned to each
quantization level $i$ where $i\in[0,Q]$. Given two sets of HVs, $\mathbf{F}$
and $\mathbf{L}$, a spectrum vector is encoded into a HV $\mathbf{I}$ as
follows. Let $\mathds{P}$ be the set of peaks in the spectrum vector,
consisting of tuples $(i,j)$, with $i$ the peak index and $j$ the step value
of its intensity. $\mathbf{I}$ is computed as follows:

(1) $\mathbf{I}=\sum_{(i,j)\in\mathds{P}}\mathbf{F}_{i}\odot\mathbf{L}_{j},$

where $\odot$ indicates element-wise multiplication. In turn, we binarize the
$\mathbf{I}$ for the computational efficiency on hardware; all positive
elements are mapped to $1$ and $-1$ otherwise. The final representation of the
HV is a binary vector. Spectrum vectors corresponding to the query and
reference spectra are encoded to query HVs and reference HVs, respectively.
Encoding of reference spectra is done only once. The reference HVs is reused
for subsequent runs since they are already identified and unlikely to change.

####  3.3.1. Reflecting Spatial Locality

We devise the new position HV generation method to reflect the spatial
locality. Previous studies (Salamat et al., 2019, 2020) have used a
permutation-based or random generation method, which makes $\mathbf{F}_{i}$
and $\mathbf{F}_{j}$ ($i\neq j$) nearly orthogonal. However, they are
vulnerable to peak shifts as they accompany changes in $i$.

Figure 4(a) shows the proposed position HV generation strategy. We randomly
generate $\mathbf{F}_{1}=\\{+1,-1\\}^{D}$. In turn, we flip $\alpha$
components in random positions. As more flips occur, the similarity between
the original vector and the flipped vector decreases. For example, the
similarity ($\delta_{1}$) between $\mathbf{F}_{1}$ and $\mathbf{F}_{2}$ is
larger than the similarity ($\delta_{2}$) between $\mathbf{F}_{1}$ and
$\mathbf{F}_{f}$. The HyperOMS encoding reflects the characteristics of peak
shifts in OMS well: (1) neighboring positions should have spatial locality to
deal with peak shifts, while (2) distant positions need to have adequate
orthogonality. The peak shift changes the index value corresponding to the
intensity in the spectrum vector. Position HVs do not change significantly
even if peak shifts occur; the resulting representation can be tolerable to
PTMs. As depicted in Figure 4(b), for $\mathbf{F}_{i}$ and $\mathbf{F}_{j}$,
the pairwise similarity has a high value when $i\approx j$ and is maximum when
$i=j$ (diagonal elements). Note that we scaled down the $f$ to $128$ and $D$
to $256$ for better visibility.

####  3.3.2. Reflecting Value Locality

The intensity information of the spectrum vectors is captured. We use the
level HV generation method from (Imani et al., 2017; Kim et al., 2018; Salamat
et al., 2019). We allocate a single bit to each of the HV components, i.e.,
$\mathbf{L}_{i}\in\\{-1,1\\}^{D}$. $\mathbf{L}$ that is assigned to each
quantization level needs to reflect the closeness of the intensity, i.e.,
value locality. The similarity between $\mathbf{L}_{i}$ and $\mathbf{L}_{i+1}$
is higher than the similarity between $\mathbf{L}_{i}$ and
$\mathbf{L}_{i+100}$. For instance, for the target level $p$ in percentage,
$\mathbf{L}_{p}$, we can represent this by flipping $(D/2)\times(p/100)$
elements of $\mathbf{L}_{0}$.

###  3.4. Hamming Similarity Search

After the encoding step, HyperOMS finds the matched reference HV that is most
similar to the query HV. It uses Hamming similarity (defined by the number of
equal components in vector pairs) as a similarity metric. Here, reference
spectra that need to be compared primarily need to satisfy spectrum charge and
precursor $m$/$z$ condition per query as discussed in Section 2.1. Each
spectrum has its spectrum charge (+2, +3, $\ldots$) and precursor $m$/$z$
value. We gather reference spectra that (1) have the same spectrum charge as
the query spectra and (2) satisfy the precursor $m$/$z$ tolerance (precursor
$m$/$z$ difference between query and reference) condition.

OMS assumes that precursor $m$/$z$ of selected reference spectra and query
spectra can have a large difference. In other words, a wide precursor $m$/$z$
tolerance is used to match modified spectra to their unmodified variants.
However, we may miss the case of a reference spectrum with a similar precursor
$m$/$z$ that can pass through the FDR filter with high similarity. To avoid
such misidentifications, we adopt cascade search (Kertesz-Farkas et al.,
2015). A narrow precursor $m$/$z$ tolerance is used for standard search and
FDR filtration is applied (Figure
3(c)-$\mathbin{\vphantom{\circledast}\text{\ooalign{\smash{\raisebox{-2.58334pt}{\scalebox{2.3}{$\bullet$}}}\cr\smash{{\color[rgb]{1,1,1}\bf\footnotesize$1$}}\cr}}}$).
In turn, remaining unidentified spectra are processed with a wide precursor
$m$/$z$ tolerance (Figure
3(c)-$\mathbin{\vphantom{\circledast}\text{\ooalign{\smash{\raisebox{-2.58334pt}{\scalebox{2.3}{$\bullet$}}}\cr\smash{{\color[rgb]{1,1,1}\bf\footnotesize$2$}}\cr}}}$).

##  4\. Evaluation

###  4.1. Methodology

The evaluation was performed on a system equipped with Intel i7-8700K with
64GB RAM and NVIDIA Geforce GTX 1080Ti. We have implemented the HyperOMS
algorithm on GPU based on (Kang et al., 2022b, a). To maximize the computation
efficiency, we represent binary HV as a $32$bit integer array using bit
packing, and similarity score computation is done by CUDA intrinsic (__popc).
Since the GPU has limited memory, we split the reference and the query data
into batches. We set the batch size to use the maximum amount of VRAM for GPU-
based solutions. We measured the energy consumption of the CPU and GPU using
Intel Power Gadget and nvidia-smi, respectively.

Workloads. We evaluated HyperOMS on two real-world datasets. Small scale
dataset combines yeast spectral library (Selevsek et al., 2015) with the human
HCD spectral library (total spectra: $1,162,392$) as the reference libraries,
and the iPRG2012 dataset (Chalkley et al., 2014) (total spectra: $15,867$) as
a query. For a larger-scale evaluation, we used the human spectral library
(mas, 2022a; Wang et al., 2018) (total spectra: $2,992,672$) and a HEK293
(Human Embryonic Kidney 293) dataset (Chick et al., 2015) (total spectra per
query: $46,665$ on average), as reference and query spectra, respectively.
Note that HEK293 consists of multiple query files (b1906$\sim$b1938).

We preprocessed query and reference spectra in a similar fashion to existing
works (Bittremieux et al., 2019, 2018a; Lam et al., 2007), using the widely
used configurations listed in Table 1. We removed the duplicates and added
decoy spectra with the same ratio as the existing spectral libraries using the
shuffle-and-reposition method (Lam et al., 2010) for FDR filtering.

Table 1. Spectrum preprocessing settings. | Dataset
---|---
Parameter Name | Small scale | Large scale
Max peaks in spectra | 50
Min / max $m$/$z$ | 101 / 1500
Bin size | 0.05 | 0.04
Precursor $m$/$z$ tolerance (narrow) | 20ppm | 5ppm
Precursor $m$/$z$ tolerance (wide) | 500Da | 500Da

Benchmarks. We compare the search quality of HyperOMS to existing search
tools, including (1) SpectraST (Lam et al., 2007; Ma et al., 2014) that run on
CPU and (2) the state-of-the-art OMS tool, ANN-SoLo, running on CPU
(Bittremieux et al., 2018a) and GPU (Bittremieux et al., 2019). We count the
number of identifications to compare the search quality. All search results
are evaluated at fixed $1\%$ FDR threshold, using Pyteomics (Goloborodko et
al., 2013).

(a)

(b)

(c)

(d)

Figure 5. Search quality comparison on the small scale dataset. (a) HyperOMS
with different HV dimensionality. (b) HyperOMS with different number of bit
flips. (c) HyperOMS with various encoding strategies. (d) Results of baseline
tools.

###  4.2. Impact of Encoding Configuration

HV dimensionality. The HV dimensionality ($D$) plays a critical role in search
quality. A low $D$ limits separability. Figure 5(a) describes the impact of
$D$. The higher $D$ leads to a higher number of identifications. However, the
excessive $D$ leads to an increase in computation and capacity demand. We set
the $D$ to $8192$ (8k).

Flipped bits. The number of flipped bits ($\alpha$) controls the balance
between orthogonality and correlation between bins. A high $\alpha$ increases
the orthogonality of each position, and a low $\alpha$ helps a more number of
adjacent bins to have a correlation (spatial locality). We measured the number
of identifications according to the ratio of flips to $D$, i.e., $\alpha/D$.
As shown in Figure 5(b), an adequate $\alpha/D$ leads to high search quality.
But if we flip a small number of bits, HyperOMS cannot clearly differentiate
the peak position. Also, since the peak shifts due to PTMs are not
significant, spatial locality for a limited range is required. In the rest of
our evaluation, we use $\alpha=D/2$, which shows the highest performance in
most cases.

Encoding strategy. We compare the search quality of HyperOMS with the
different binary encoding strategies. As discussed in Section 3.3, LSH (Koga
et al., 2006; Norouzi et al., 2012; Imani et al., 2020), ID-Level HDC encoding
(Imani et al., 2017; Kim et al., 2018; Salamat et al., 2019) can be used
alternatively to map raw data to binary vectors. Figure 5(c) compares the
search quality of HyperOMS with the (1) proposed encoding method, (2) ID-Level
HDC encoding that can capture the position of feature and its value, and (3)
random projection-based LSH approach. Our encoding offers the best search
quality compared to baselines.

Quantization level. High quantization levels may not be flexible to the peak
intensity changes due to noise and PTMs. Low $Q$ leads to low resolution in
intensity capturing of the encoder. The quantization level $Q$ did not
significantly affect the search quality unless it falls in $[8,32]$.
Therefore, we use $Q=16$.

(a)

(b)

Figure 6. Search result analysis from HyperOMS and baseline tools on the large
scale dataset. (a) Total number of identifications. (b) Venn diagram of b1927
query search result.

###  4.3. Search Quality

We searched the iPRG2012 dataset against the yeast spectral library. As no
ground truth information is available when analyzing complex biological data,
instead, we compare our search quality with a list of consensus
identifications produced by multiple search tools during the iPRG2012 study
(Chalkley et al., 2014). Figure 5(d) shows the search result of baseline
tools. Among $7841$ identifications in the iPRG2012 consensus result, HyperOMS
is able to correctly identify $4141$ spectra (Figure 5(b)). SpectraST and ANN-
SoLo manage to identify $3891$ and $5327$ identifications, respectively.

We compared the performance of HyperOMS with the results from existing tools,
including SpectraST and ANN-SoLo using the large-scale dataset. We used
similar configurations for all tools, listed in Table 1. Figure 6 shows the
number of identifications from the different search tools.

HyperOMS offers a higher search quality than SpectraST, i.e., HyperOMS
identifies more spectra. ANN-SoLo managed to identify more spectra than our
HyperOMS. Nevertheless, as described in Figure 6(b), HyperOMS can identify
spectra that other tools can find (overlapped area). HyperOMS represents
spectra with an approximated form of the original data, which is robust to
PTMs. Therefore, HyperOMS can use Hamming similarity to perform OMS. Besides,
ANN-SoLo uses shifted cosine similarity metric, which is accurate when finding
the original spectra.

The identification rate of HyperOMS can be improved by increasing the HV
capacity. This can be done by (1) increasing $D$ or (2) increasing the
precision of each component in the HV. For example, increasing $D$ from 8k
(8192) to 16k (16384) can yield up to 10% more identifications. However, it
raises the hardware cost, computational complexity, and energy consumption of
the accelerator. Since our main goal is to maximize the speed and energy
efficiency of the OMS while achieving reasonable quality in a biological
sense, we use a binary HV with 8k dimensionality.

A ramification of lower search quality could be missing potentially relevant
biomarker proteins in the context of a healthy versus diseased study, or
missing data similarly impacting other downstream biological interpretations.
Nevertheless, the HyperOMS identification rate is within the range of the
state-of-the-art in MS identification. For example, we can typically expect an
identification rate of 33–66% currently for human samples that we have used,
and HyperOMS satisfies the expected range criterion. One advantage of HyperOMS
is that there is a search quality–efficiency trade-off that can be tuned using
the hyperparameters. Furthermore, the user can decide between different search
engines based on their requirements. For example, HyperOMS runs much faster
with superior energy efficiency compared to existing OMS tools (Section 4.4).
It could be used to efficiently process extremely large proteomics datasets
consisting of tens of thousands of query files, which are being generated
increasingly often recently.

Figure 7. Performance and energy efficiency comparison.

###  4.4. Speed and Energy Efficiency Improvement

We compare the execution time and the energy consumption of HyperOMS on GPU to
the state-of-the-art OMS tool ANN-SoLo, which offers a faster search speed
with the GPU acceleration than other baseline (Bittremieux et al., 2019). Note
that we measured the second run of the ANN-SoLo since the reference data is
likely to be pre-encoded in reality. ANN-SoLo saves the pre-indexed
information of the reference library on the first run and reuses it in the
subsequent run. For the large-scale dataset, we averaged the measurements from
multiple queries (b1906 $\sim$ b1938).

Figure 7 compares the end-to-end runtime. ANN-SoLo builds the index on the CPU
while the encoding of HyperOMS is done on the GPU. The HyperOMS encoding is
parallelized over HV dimensions and datapoints. The encoding stage of
HyperOMS, which corresponds to the index build of ANN-SoLo, is up to
$8.6\times$ faster than ANN-SoLo. HyperOMS uses HD binary vector and easily
parallelizable Hamming similarity computation, while ANN-SoLo uses FP32
vector. The search process of HyperOMS GPU achieves on average $82\times$
speedup over ANN-SoLo on CPU and $11.2\times$ speedup ANN-SoLo on GPU.
Ultimately, HyperOMS GPU offers significant speedup in all stages of the OMS.
For the end-to-end execution, HyperOMS GPU gains an average speedup of
$15.7\times$ over the state-of-the-art OMS tool running on the same GPU.

Besides, HyperOMS running on the GPU requires more power than the ANN-SoLo, as
it has high parallelism. However, the increased power consumption is
compensated by reduced execution time, improving energy efficiency. Overall,
HyperOMS results in $7.8\times$ and $5\times$ energy efficiency improvement
over CPU and GPU on average, respectively, as shown in Figure 7.

##  5\. Conclusion

We proposed HyperOMS, HDC-inspired massively parallel algorithm for OMS of MS-
based proteomics. The proposed algorithm encodes spectra into binary HVs,
considering the spatial and value locality of peaks. Therefore, HyperOMS
simplifies the execution pipeline and maximizes the computation efficiency and
parallelism by using a binary vector with boolean operations. Our evaluation
results show that HyperOMS offers comparable search quality to existing OMS
tools. Furthermore, HyperOMS on NVIDIA Geforce GTX 1080Ti yields up to
$17\times$ speedup and $6.4\times$ improved energy efficiency over the state-
of-the-art GPU-based OMS solution.

###### Acknowledgements.

This work was supported in part by CRISP, one of six centers in JUMP (an SRC
program sponsored by DARPA), SRC Global Research Collaboration (GRC) grant,
and NSF grants #1826967, #1911095, #2052809, #2112665, #2112167, and #2100237.

## References

  * (1)
  * mas (2022a) 2022a.  _MassIVE-KB_.  https://massive.ucsd.edu/ProteoSAFe/static/massive.jsp
  * mas (2022b) 2022b.  _MassIVE: Mass Spectrometry Interactive Virtual Environment_.  https://massive.ucsd.edu/
  * Adams et al. (2022) C. Adams et al. 2022\.  Open modification searching of SARS-CoV-2–human protein interaction data reveals novel viral modification sites.  (March 2022).
  * Aebersold et al. (2016) R. Aebersold et al. 2016\.  Mass-Spectrometric Exploration of Proteome Structure and Function.  _Nature_ 537, 7620 (Sept. 2016), 347–355.
  * Aksenov et al. (2017) A. A Aksenov et al. 2017\.  Global chemical analysis of biology by mass spectrometry.  _Nature Reviews Chemistry_ 1, 7 (2017), 1–20.
  * Baumgardner et al. (2011) L. A. Baumgardner et al. 2011\.  Fast Parallel Tandem Mass Spectral Library Searching Using GPU Hardware Acceleration.  _J. Proteome Res._ 10, 6 (June 2011), 2882–2888.
  * Bittremieux et al. (2018a) W. Bittremieux et al. 2018a.  Fast Open Modification Spectral Library Searching through Approximate Nearest Neighbor Indexing.  _J. Proteome Res._ 17, 10 (Sept. 2018), 3463–3474.
  * Bittremieux et al. (2018b) W. Bittremieux et al. 2018b.  Quality Control in Mass Spectrometry-Based Proteomics.  _Mass Spectrom. Rev._ 37, 5 (Sept. 2018), 697–711.
  * Bittremieux et al. (2019) W. Bittremieux et al. 2019\.  Extremely Fast and Accurate Open Modification Spectral Library Searching of High-Resolution Mass Spectra Using Feature Hashing and Graphics Processing Units.  _J. Proteome Res._ 18, 10 (Aug. 2019), 3792–3799.
  * Chalkley et al. (2014) R. J Chalkley et al. 2014\.  Proteome Informatics Research Group (iPRG)_2012: A Study on Detecting Modified Peptides in a Complex Mixture.  _MCP_ 13, 1 (Jan. 2014), 360–371.
  * Chi et al. (2018) H. Chi et al. 2018\.  Comprehensive Identification of Peptides in Tandem Mass Spectra Using an Efficient Open Search Engine.  _Nat. Biotechnol._ (Oct. 2018).  https://doi.org/10.1038/nbt.4236
  * Chick et al. (2015) J. M Chick et al. 2015\.  A Mass-Tolerant Database Search Identifies a Large Proportion of Unassigned Spectra in Shotgun Proteomics as Modified Peptides.  _Nat. Biotechnol._ 33, 7 (June 2015), 743–749.
  * Devabhaktuni et al. (2019) A. Devabhaktuni et al. 2019\.  TagGraph Reveals Vast Protein Modification Landscapes from Large Tandem Mass Spectrometry Datasets.  _Nat. Biotechnol._ 37, 4 (April 2019), 469–479.
  * Elias et al. (2007) Joshua E Elias et al. 2007\.  Target-Decoy Search Strategy for Increased Confidence in Large-Scale Protein Identifications by Mass Spectrometry.  _Nature Methods_ 4, 3 (Feb. 2007), 207–214.
  * Goloborodko et al. (2013) A. A. Goloborodko et al. 2013\.  Pyteomics—a Python Framework for Exploratory Data Analysis and Rapid Software Prototyping in Proteomics.  _J. Am. Soc. Mass Spectrom._ 24, 2 (Jan. 2013), 301–304.
  * Griss et al. (2013) J. Griss et al. 2013\.  PRIDE Cluster: Building a Consensus of Proteomics Data.  _Nat. Methods_ 10, 2 (Jan. 2013), 95–96.
  * Heck et al. (2020) M. Heck et al. 2020\.  Proteomics in Non-model Organisms: A New Analytical Frontier.  _J. Proteome Res._ 19, 9 (2020), 3595–3606.
  * Imani et al. (2017) M. Imani et al. 2017\.  Voicehd: Hyperdimensional computing for efficient speech recognition. In _ICRC_.
  * Imani et al. (2020) M. Imani et al. 2020\.  DUAL: Acceleration of Clustering Algorithms using Digital-based Processing In-Memory. In _MICRO_.
  * Johnson et al. (2017) J. Johnson et al. 2017\.  Billion-scale similarity search with GPUs.  _arXiv preprint arXiv:1702.08734_ (2017).
  * Kanerva (2009) P. Kanerva. 2009\.  Hyperdimensional Computing: An Introduction to Computing in Distributed Representation with High-Dimensional Random Vectors.  _Cognitive Computation_ 1, 2 (Jan. 2009), 139–159.
  * Kang et al. (2022a) J. Kang et al. 2022a.  OpenHD: A GPU-Powered Framework for Hyperdimensional Computing.  _IEEE TC_ (2022).
  * Kang et al. (2022b) J. Kang et al. 2022b.  XCelHD: An Efficient GPU-Powered Hyperdimensional Computing with Parallelized Training. In _ASP-DAC_. 220–225.
  * Kertesz-Farkas et al. (2015) Attila Kertesz-Farkas et al. 2015\.  Tandem Mass Spectrum Identification via Cascaded Search.  _J. Proteome Res._ 14, 8 (Aug. 2015), 3027–3038.
  * Kim et al. (2018) Y. Kim et al. 2018\.  Efficient human activity recognition using hyperdimensional computing. In _IoT_. 1–6.
  * Koga et al. (2006) H. Koga et al. 2006\.  Fast agglomerative hierarchical clustering algorithm using Locality-Sensitive Hashing.  _KAIS_ 12, 1 (July 2006), 25–53.
  * Kong et al. (2017) A. T Kong et al. 2017\.  MSFragger: ultrafast and comprehensive peptide identification in mass spectrometry–based proteomics.  _Nat. Methods_ 14, 5 (April 2017), 513–520.
  * Lam (2011) H. Lam. 2011.  Building and searching tandem mass spectral libraries for peptide identification.  _MCP_ 10, 12 (2011).
  * Lam et al. (2007) H. Lam et al. 2007\.  Development and validation of a spectral library searching method for peptide identification from MS/MS.  _PROTEOMICS_ 7, 5 (March 2007), 655–667.
  * Lam et al. (2010) H. Lam et al. 2010\.  Artificial decoy spectral libraries for false discovery rate estimation in spectral library searching in proteomics.  _J. Proteome Res._ 9, 1 (2010), 605–610.
  * Li et al. (2014) Y. Li et al. 2014\.  Accelerating the scoring module of mass spectrometry-based peptide identification using GPUs.  _BMC Bioinform._ 15, 1 (2014), 121.
  * Ma et al. (2014) C. W. M. Ma et al. 2014\.  Hunting for Unexpected Post-Translational Modifications by Spectral Library Searching with Tier-Wise Scoring.  _J. Proteome Res._ 13, 5 (April 2014), 2262–2271.  https://doi.org/10.1021/pr401006g
  * Na et al. (2015) S. Na et al. 2015\.  Software Eyes for Protein Post-Translational Modifications.  _Mass Spectrom. Rev._ 34, 2 (April 2015), 133–147.
  * Norouzi et al. (2012) M. Norouzi et al. 2012\.  Hamming Distance Metric Learning. In _NeurIPS_ , Vol. 25.
  * Perez-Riverol et al. (2022) Y. Perez-Riverol et al. 2022\.  The PRIDE database resources in 2022: a hub for mass spectrometry-based proteomics evidences.  _Nucleic acids research_ 50, D1 (2022).
  * Poduval et al. (2021) P. Poduval et al. 2021\.  Robust In-Memory Computing with Hyperdimensional Stochastic Representation. In _NANOARCH_. 1–6.
  * Salamat et al. (2019) S. Salamat et al. 2019\.  F5-hd: Fast flexible fpga-based framework for refreshing hyperdimensional computing. In _FPGA_. ACM, 53–62.
  * Salamat et al. (2020) S. Salamat et al. 2020\.  Accelerating Hyperdimensional Computing on FPGAs by Exploiting Computational Reuse.  _IEEE TC_ 69, 8 (2020), 1159–1171.
  * Selevsek et al. (2015) N. Selevsek et al. 2015\.  Reproducible and Consistent Quantification of the Saccharomyces cerevisiae Proteome by SWATH-mass spectrometry.  _MCP_ 14, 3 (March 2015), 739–749.
  * Smith et al. (2014) R. Smith et al. 2014\.  Proteomics, Lipidomics, Metabolomics: A Mass Spectrometry Tutorial from a Computer Scientist’s Point of View.  _BMC Bioinform._ 15, Suppl 7 (2014), S9.
  * Walsh et al. (2005) C. T Walsh et al. 2005\.  Protein Posttranslational Modifications: The Chemistry of Proteome Diversifications.  _Angew. Chem. Int. Ed._ 44, 45 (Nov. 2005), 7342–7372.
  * Wang et al. (2018) M. Wang et al. 2018\.  Assembling the Community-Scale Discoverable Human Proteome.  _Cell Systems_ 7, 4 (Oct. 2018), 412–421.e5.
  * Yang et al. (2019) M. Yang et al. 2019\.  A Complete CPU-FPGA Architecture for Protein Identification with Tandem Mass Spectrometry. In _ICFPT_. 295–298.
  * Zhang et al. (2021) S. Zhang et al. 2021\.  Assessing Robustness of Hyperdimensional Computing Against Errors in Associative Memory. In _ASAP_. 211–217.