# Incoherent witnessing of quantum coherence

Sahar Basiri-Esfahani  Electronic address<EMAIL_ADDRESS>Department of Physics, Swansea University, Singleton Park, Swansea SA2 8PP,
United Kingdom  Farid Shahandeh  Electronic address<EMAIL_ADDRESS>Department of Physics, Swansea University, Singleton Park, Swansea SA2 8PP,
United Kingdom

###### Abstract

Theoretical and experimental studies have suggested the relevance of quantum
coherence to the performance of photovoltaic and light-harvesting complex
molecular systems. However, there are ambiguities regarding the validity of
statements we can make about the coherence in such systems. Here we analyze
the general procedure for coherence detection in quantum systems and show the
counterintuitive phenomenon of detecting a quantum system’s initial coherence
when both the input and output probe states are completely incoherent. Our
analysis yields the necessary and sufficient conditions for valid claims
regarding the coherence of directly inaccessible systems. We further provide a
proof-of-principle protocol that uses entangled probes to detect quantum
coherence satisfying these conditions, and discuss its potency for detecting
coherence.

Quantum coherence—that quantum systems occupy multiple states simultaneously
and hence exhibit interference—is the distinctive feature of quantum systems
compared to classical ones. In recent years, quantum coherence has become a
critical element of developing quantum technologies that aim at improving over
classical approaches Hayes  _et al._ (2013); Scholes  _et al._ (2017); Scully
(2010); Kimble (2008); Streltsov  _et al._ (2017). Of particular importance is
the potential role of quantum coherence in the light-harvesting efficiency of
biochemical processes, e.g., photosynthesis Scholes (2011); Baghbanzadeh and
Kassal (2016); Scholes (2010); Chenu and Scholes (2015), as well as enhancing
the performance of molecular systems such as organic solar cells Tomasi  _et
al._ (2019); Tomasi and Kassal (2020); Fruchtman  _et al._ (2016); Tomasi  _et
al._ (2021); Creatore  _et al._ (2013). This has stimulated studies of the
effects of quantum coherence in such systems, at the heart of which lies
schemes to certify the presence of quantum coherent mechanisms in molecular
systems Engel  _et al._ (2007); Collini  _et al._ (2010); Panitchayangkoon
_et al._ (2011); Abramavicius and Mukamel (2010); Panitchayangkoon  _et al._
(2010).

The main tool to examine quantum coherence in complex systems including
photosynthetic complexes is the ultrafast spectroscopy Engel  _et al._ (2007);
Collini  _et al._ (2010); Panitchayangkoon  _et al._ (2011); Abramavicius and
Mukamel (2010); Panitchayangkoon  _et al._ (2010). Spectroscopic observations,
however, have led to debates Kassal  _et al._ (2013); Wang  _et al._ (2019);
Brumer and Shapiro (2012) mainly because the dynamics of systems in nature, as
opposed to the spectroscopic techniques, is initiated by incoherent inputs
such as sunlight Brumer and Shapiro (2012); Mančal and Valkunas (2010); Jiang
and Brumer (1991). These arguments indicate the need for further investigating
the existence of quantum coherence in systems operating under ambient
conditions and proposals of new protocols to detect quantum coherence using
incoherent light sources rather than coherent lasers Muñoz and Schlawin
(2020). Furthermore, it is intuitive to assume that when quantum channels
suffer from significant decohereing noise the coherence of the system is
untraceable. We approach these arguments from a quantum informational
perspective and pose the following question (Fig. 1): Is it possible to make
deductions about the coherence properties of a system when both the input and
output probe states are incoherent?

Here, we answer to the above question in the affirmative. We show via a
counter-example that quantum coherence of a system can be observed even if
both input and output to the process are fully incoherent. We also provide
rigorous necessary and sufficient conditions for this observation to be
possible. We propose a proof-of-principle protocol to detect the existence of
quantum coherence within a system using incoherent light sources and discuss
the physics behind this counter-intuitive phenomenon.

Figure 1: A coherence-generating quantum channel with incoherent input and
output. In such circumstances, the quantum coherence generated in a directly
inaccessible system seems undetectable.

Let us begin with the description of quantum coherence from the perspective of
quantum information science Streltsov  _et al._ (2017). In this picture,
quantum coherence comprises of the following ingredients. First, given a
system described in a finite-dimensional Hilbert space, an orthonormal basis
${\cal B}_{\rm inc}{=}\\{|i\rangle\\}$ is defined as the incoherent basis.
These vectors represent the basis with respect to which we would like to
consider the quantum coherence properties of the system. Since basis vectors
$|i\rangle$ are also pure quantum states, they define all incoherent quantum
states as probabilistic mixtures of the incoherent basis states. In other
words, every density operator $\varrho_{\rm inc}$ that can be written as
$\varrho_{\rm inc}{=}\sum_{i}p_{i}|i\rangle\\!\langle i|$, for
$|i\rangle{\in}{\cal B}_{\rm inc}$, $p_{i}{\geq}0$, and $\sum_{i}p_{i}{=}1$,
is an incoherent state. The terminology is indeed justified by noting that
every $\varrho_{\rm inc}$ is diagonal in the basis ${\cal B}_{\rm inc}$ with
no off-diagonal matrix elements. We denote the collection of incoherent
density operators by ${\cal S}_{\rm inc}$.

The second ingredient of the theory is the set of incoherent channels denoted
by ${\cal C}_{\rm inc}$. A quantum channel is a transformation that converts
any input density operator to an output density operator. Most generally, a
channel is called incoherent if upon receiving an incoherent input quantum
state outputs an incoherent quantum state. Several classes of such
transformations have been studied within the quantum information literature so
far Chitambar and Gour (2016a); Streltsov  _et al._ (2017).

In order to better understand how these channels look like, it is useful to
introduce a particular transformation called the fully dephasing channel
$\Delta$ that acts on any quantum state as
$\Delta[\varrho]:=\sum_{i}\mbox{$\langle
i|\varrho|i\rangle$}|i\rangle\\!\langle
i|=\sum_{i}\varrho_{i}|i\rangle\\!\langle i|$. Here, we have introduced the
shorthand $\varrho_{i}{:=}\mbox{$\langle i|\varrho|i\rangle$}$ for the
diagonal matrix elements of the density operator $\varrho$ in the incoherent
basis ${\cal B}_{\rm inc}$, or simply the population of each incoherent state
of the system. It is thus evident that the action of the fully dephasing
channel is to completely suppress the coherence of the input quantum state.
Incoherent channels can now be defined as all channels $\Lambda$ that satisfy
Chitambar and Gour (2016b, a)

$\Delta\circ\Lambda[\varrho_{\rm inc}]=\Lambda[\varrho_{\rm inc}]$ (1)

for all incoherent input quantum states $\varrho$. Here, $\circ$ denotes the
composition of quantum channels, meaning that each channel is consecutively
applied to the output of the previous channel to its right. Equation (1) has a
very intuitive physical interpretation: a channel is incoherent if and only if
the further application of a fully dephasing channel on its output is
redundant for every incoherent input. Equation (1) can also be rewritten in
terms of all input quantum states as

$\Delta\circ\Lambda\circ\Delta[\varrho]=\Lambda\circ\Delta[\varrho],$ (2)

wherein a fully dephasing channel has been used to transform the arbitrary
input quantum state into an incoherent one.

The third and final ingredient of the theory of quantum coherence are
incoherent measurements. It is known that every quantum measurement is
described by a set of effects $\mathsf{M}{=}\\{M_{k}\\}$. An effect $M_{k}$ is
positive, i.e., it is Hermitian with nonnegative eigenvalues, and corresponds
to the outcome $k$ of the measurement $\mathsf{M}$. The collection
$\mathsf{M}$ such that $\sum_{k}M_{k}{=}\mathds{1}$ is called a positive
operator-valued measure (POVM) Nielsen and Chuang (2010). Furthermore,
according to the Born rule, given the quantum state $\varrho$ the probability
of outcome $k$ in measurement $\mathsf{M}$ is given by
$p(k|\varrho,\mathsf{M})=\operatorname{Tr}M_{k}\varrho$.

Now, a measurement $\mathsf{M}_{\rm inc}$ is called incoherent if the matrix
representation of each of its effects is diagonal with respect to the
incoherent basis ${\cal B}_{\rm inc}$ Streltsov  _et al._ (2017). The simplest
example of an incoherent measurement is indeed the projective measurement onto
the incoherent basis, i.e. $\mathsf{M}_{\rm
inc}{=}\\{P_{i}{=}|i\rangle\\!\langle i|\\}{:=}\Pi_{\rm inc}$. We denote the
collection of all incoherent measurements by ${\cal M}_{\rm inc}$.

Channels and the probe-system interaction.— The concept of a channel is very
versatile. Consider the scenario in which a probe interacts with a system and
is then measured. It is implicit that the system’s degree of freedom of
interest cannot directly be accessed without the mediation of the probe. As a
result, all we can speak of is the initial and the final quantum state of the
probe. In other words, the system—including its initial and final quantum
states—and its interaction with the probe are subsumed by the quantum channel
that transforms the probe.

A generic channel $\Lambda$ transforming an input state to the output can be
written as Nielsen and Chuang (2010)

$\varrho^{(\rm out)}=\Lambda[\varrho^{(\rm in)}]=\sum_{i}K_{i}\varrho^{(\rm
in)}K_{i}^{\dagger}.$ (3)

Here, the operators $K_{i}$ are called Kraus operators of the channel and
satisfy $\sum_{i}K_{i}^{\dagger}K_{i}{=}\mathds{1}$. This way, it is
guaranteed that the output is a valid normalized density operator; see
Supporting Information-1 (SI-1). The simplest example of a quantum channel is
indeed a unitary one that acts on the input as

$\varrho^{(\rm out)}=\Upsilon[\varrho^{(\rm in)}]=U\varrho^{(\rm
in)}U^{\dagger}.$ (4)

An important question following Eq. (3) is that what happens to the unitarity
and reversibility of quantum mechanics. A fundamental theorem of quantum
mechanics implies that every generic quantum channel as in Eq. (3) can be
purified to a unitary channel acting on the input state and some inaccessible
auxiliary system Nielsen and Chuang (2010), that is,

$\varrho^{(\rm out)}=\Lambda[\varrho^{(\rm in)}]=\operatorname{Tr}_{\rm
a}\Upsilon[\varrho^{(\rm in)}\otimes\pi_{\rm a}^{(\rm in)}].$ (5)

In Eq. (5), $\pi_{\rm a}$ is a suitably chosen state of the ancillary
subsystem ${\rm a}$, $\Upsilon$ is a suitably chosen interaction unitary
channel acting on the joint input-ancilla compound, and
$\operatorname{Tr}_{\rm a}$ is the partial trace taken over the ancilla. We
note that the partial trace represents the inaccessibility of the auxilliary
subsystem to the experimenter.

Equation (5) can be brought to the context of probe-system interaction by
interpreting the input subsystem as the probe (from now on denoted as
$\varrho_{\rm p}$) and the ancilla subsystem as the system of interest (from
now on denoted as $\pi_{\rm s}$). Neglecting the interactions with the
environment for the moment, $\Upsilon$ can be interpreted as the unitary
interaction between the probe and the system, denoted by $\Upsilon_{\rm ps}$.
The inaccessibility of the ancilla thus translates to the fact that direct
measurements on the system of interest are out of rich. Eq. (5) can thus be
rewritten as

$\varrho_{\rm p}^{(\rm out)}=\Lambda[\varrho_{\rm p}^{(\rm
in)}]=\operatorname{Tr}_{\rm s}\Upsilon_{\rm ps}[\varrho_{\rm p}^{(\rm
in)}\otimes\pi_{\rm s}^{(\rm in)}],$ (6)

that allows us to connect the properties of the inaccessible system to the
properties of the channel acting on the probe alone.

In the present work we are interested in the initial coherence properties of
the system $\pi_{\rm s}^{(\rm in)}$ by merely observing the output probe state
$\varrho_{\rm p}^{(\rm out)}$. It is evident from Eq. (6) that there are three
elements that can cause observable coherent effects in the probe, namely, the
initial system state $\pi_{\rm s}^{(\rm in)}$, the initial probe state
$\varrho_{\rm p}^{(\rm in)}$, and the probe-system interaction $\Upsilon_{\rm
ps}$. Thus, to make valid statements about the coherence of $\pi_{\rm s}^{(\rm
in)}$, we must make sure the latter two potential causes are not in effect by
making two assumptions:

  1. (i)

The input state of the probe is incoherent;

  2. (ii)

The channel $\Lambda$ is incoherent for all incoherent input states of the
system and probe.

We note that assumption (i) necessarily prevents objections of the kind
associated with spectroscopic techniques Kassal  _et al._ (2013); Wang  _et
al._ (2019); Brumer and Shapiro (2012); Mančal and Valkunas (2010); Jiang and
Brumer (1991). Now, let us posit the appropriate incoherent bases for the
probe and the system to be ${\cal B}_{\rm inc;p}$ and ${\cal B}_{\rm inc;s}$,
respectively. Assumption (ii) can then be expressed as

$\sum_{\phi}{{}_{\rm ps}\\!}\langle k,\phi|U_{\rm
ps}|i^{\prime},j^{\prime}\rangle_{\rm ps}\langle i^{\prime},j^{\prime}|U_{\rm
ps}^{\dagger}|l,\phi\rangle_{\rm ps}=0,$ (7)

for all $|i^{\prime}\rangle_{\rm p},|k\rangle_{\rm p},|l\rangle_{\rm
p}\in{\cal B}_{\rm inc;p}$ with $k\neq l$ and all $|j^{\prime}\rangle_{\rm
s}\in{\cal B}_{\rm inc;s}$. Furthermore, $\\{|\phi\rangle_{\rm s}\\}$ is an
arbitrary basis for the Hilbert space of the system. Equation (7) gives us a
general restriction on the unitaries for which we can safely draw conclusions
about the coherence properties of the system based merely on the observation
of the probe.

Coherent scenario.—We can easily verify that Eq. (7) ensures that the output
in Eq. (6) remains incoherent for initially incoherent states of both system
and probe. It is now evident that whenever the output probe state is verified
to be in a coherent superposition of states in ${\cal B}_{\rm inc;p}$, that
must be due to the initial state of the system being coherent with respect to
the incoherent basis ${\cal B}_{\rm inc;s}$. In other words, the unitary
interaction $\Upsilon_{\rm ps}$ transfers the coherence of the system to the
probe. Indeed, this scenario is straightforward: it corresponds to a coherence
generating map $\Lambda$ in Eq. (6) that transforms an incoherent input
(probe) state into a coherent output (probe) state. We are thus mainly
interested in the challenging case in which the map $\Lambda$ is incoherent.

Incoherent scenario.—Whenever the output probe state of the process in Eq. (6)
is incoherent, the naive conclusion is that no signature of the system’s
initial coherence survives the incoherent process $\Lambda$. In the following,
we show that this conclusion is not correct.

Schrödinger versus Heisenber pictures.—According to the standard quantum
mechanics, the probability of outcome $k$ in the measurement $\mathsf{M}$ on
the output of the channel $\Lambda$ for the input state $\varrho^{(\rm in)}$
is given by the Born rule, that is,

$p(k|\varrho^{(\rm
in)},\Lambda,\mathsf{M})=\operatorname{Tr}(M_{k}\Lambda[\varrho^{(\rm
in)}])=\operatorname{Tr}(M_{k}\varrho^{(\rm out)}),$ (8)

where the map $\Lambda$ is defined most generally through its Kraus operators
as introduced in Eq.(3). For the special case of a unitary channel as in Eq.
(4), this reduces to

$p(k|\varrho^{(\rm
in)},\Upsilon,\mathsf{M})=\operatorname{Tr}(M_{k}\Upsilon[\varrho^{(\rm
in)}])=\operatorname{Tr}(M_{k}U\varrho^{(\rm in)}U^{\dagger}).$ (9)

We now recall from elementary quantum mechanics that Eq. (9) describes the
Born rule in the Schrödinger picture wherein the state of the system undergoes
the dynamical evolution according to the unitary $\Upsilon$ and the effect
$M_{\rm k}$ remains stationary.

We know, however, that the unitary quantum evolution can also be expressed in
the Heisenberg picture by using the rule of permutation-under-the-trace as

$p(k|\varrho^{(\rm
in)},\Upsilon,\mathsf{M})=\operatorname{Tr}(U^{\dagger}M_{k}U\varrho^{(\rm
in)}),$ (10)

wherein the dynamical evolution is associated with the observable rather than
the input state of the system.

Similarly to Eq. (10), the case of a general channel of Eq. (8) can also be
recast in the Heisenberg picture as

$p(k|\varrho^{(\rm
in)},\Lambda,\mathsf{M})=\operatorname{Tr}(M_{k}\Lambda[\varrho^{(\rm
in)}])=\operatorname{Tr}(\Lambda^{\rm d}[M_{k}]\varrho^{(\rm in)}).$ (11)

Here, $\Lambda^{\rm d}$ is called the dual channel to $\Lambda$. One can
easily workout the relation between $\Lambda$ and its dual $\Lambda^{\rm d}$
in the above equation (see SI-2) to find

$\Lambda^{\rm d}[M_{k}]:=\sum_{i}K_{i}^{\dagger}M_{k}K_{i}.$ (12)

It can be easily checked that Eq. (12) reduces to the usual passage from the
Schrödinger to the Heisenberg picture for a unitary channel.

Despite the similarities between a unitary and a generic channel in
Schrodinger and Heisenberg representations, the two have a very sharp
contrast. Suppose $\Upsilon$ is an incoherent unitary channel with respect to
the incoherent basis ${\cal B}_{\rm inc}$. Because $\Upsilon$ is invertible,
it must transform pure states in ${\cal B}_{\rm inc}$ to only pure states in
${\cal B}_{\rm inc}$. The latter must also hold for the inverse channel
$\Upsilon^{-1}[\cdot]=U^{-1}\cdot U$. Using the fact that $U^{-1}=U^{\dagger}$
we find $\Upsilon^{-1}[\cdot]=U^{\dagger}\cdot U=\Upsilon^{\rm d}[\cdot]$.
Thus, $\Upsilon^{\rm d}[\cdot]$ must also transform pure states in ${\cal
B}_{\rm inc}$ to only pure states in ${\cal B}_{\rm inc}$, i.e., the dual of
an incoherent unitary channel is also an incoherent unitary channel. This
simple correspondence between the coherence properties of a unitary channel
and its dual, however, does not hold for a generic channel Chitambar and Gour
(2016a). In other words, there are incoherent channels whose dual is not
incoherent. We now move on to show that this asymmetry can be exploited to
prove the coherence of the initial state of an inaccessible system while the
input and output probe are both incoherent.

Figure 2: Schematic of our coherence detection protocol. The probe and monitor
are initially prepared in the maximally entangled state $\varrho_{\rm
mp}^{(\rm in)}$. The probe mode then interacts with the directly inaccessible
system which is part of the channel $\Lambda$. We observe the coherence
emerging in the state of the monitor conditioned on the outcomes of an
incoherent measurement on the probe.

The protocol.—A twin-mode probe is initially prepared in the maximally
entangled state $\varrho_{\rm mp}^{(\rm in)}=|\Phi^{+}\rangle_{\rm
mp}\langle\Phi^{+}|$ where

$|\Phi^{+}\rangle_{\rm mp}=n^{-1/2}\sum_{i=0}^{n-1}|i\rangle_{\rm
m}|i\rangle_{\rm p}.$ (13)

Here, $|i\rangle_{x}{\in}{\cal B}_{\rm inc;x}$ with $x{=}{\rm m},{\rm p}$,
${\cal B}_{\rm inc;m}{=}{\cal B}_{\rm inc;p}$ and $n$ is any natural number
between two and infinity. The first mode is called the monitor and the second
mode is the probe to be interacted with the system. Then, the probe mode is
sent through the channel $\Lambda$ to interact with the system. It is
important to note that the quantum state of the probe mode given by the
marginal density operator $\operatorname{Tr}_{\rm m}\varrho_{\rm mp}^{(\rm
in)}=\sum_{i=0}^{n-1}|i\rangle_{\rm p}\langle i|/n$ is incoherent, as
required.

After the interaction between the probe and the system, an incoherent
measurement on the probe $\rm p$ is carried out. In the postprocessing stage,
the outcomes of the measurement on the monitor are sorted conditioned on the
outcomes of the measurement on the probe. Any coherence within the conditional
data, i.e., coherence within the conditional output state $\varrho_{{\rm
m}{|}k}^{(\rm out)}$ for the outcome $k$ of the probe, indicates that the
initial state of the system $\pi_{\rm s}^{(\rm in)}$ was in a coherent quantum
superposition; see Fig. 2.

In order to understand the working principles of our protocol, suppose that
the channel $\Lambda_{\rm p}$ induced by the probe-system interaction as in
Eq. (6) is asymmetric with respect to the Schrödinger and Heisenberg pictures,
i.e., its dual $\Lambda^{\rm d}_{\rm p}$ is not incoherent. We proceed step by
step according to the protocol. We have

$\begin{split}\varrho_{\rm mp}^{(\rm out)}=\mathsf{I}_{\rm
m}\otimes\Lambda_{\rm p}[\varrho_{\rm mp}^{(\rm
in)}]&=\frac{1}{n}\sum_{i,j,l}|i\rangle_{\rm m}\langle j|\otimes
K_{l}|i\rangle_{\rm p}\langle j|K_{l}^{\dagger}\\\
&=\frac{1}{n}\sum_{i,j,l}K_{l}^{\mathsf{T}}|i\rangle_{\rm m}\langle
j|K_{l}^{*}\otimes|i\rangle_{\rm p}\langle j|.\end{split}$ (14)

In the last step we have used the trick
$\mathsf{I}\otimes\Lambda[\Phi^{+}]=\Lambda^{\mathsf{T}}\otimes\mathsf{I}[\Phi^{+}]$
where $\Phi^{+}$ is the shorthand for $|\Phi^{+}\rangle\langle\Phi^{+}|$ and
$\mathsf{T}$ is the transposition operation (see SI-3 for the proof). By
postselecting on the outcome $k$ of an incoherent measurement on the probe we
obtain the conditional state

$\begin{split}\varrho_{{\rm m}{|}k}^{(\rm out)}=\mbox{$\langle k|\varrho_{\rm
mp}^{(\rm out)}|k\rangle$}_{\rm
p}=\frac{1}{n}\sum_{l}K_{l}^{\mathsf{T}}|k\rangle_{\rm m}\langle
k|K_{l}^{*}.\end{split}$ (15)

Finally, we examine the off-diagonal elements of the monitor’s conditional
state, that is,

$\begin{split}\mbox{$\langle i|\varrho_{{\rm m}{|}k}^{(\rm
out)}|j\rangle$}_{\rm m}&=\frac{1}{n}\sum_{l}\langle
i|K_{l}^{\mathsf{T}}|k\rangle_{\rm m}\langle k|K_{l}^{*}|j\rangle\\\
&=\frac{1}{n}\sum_{l}\langle j|K_{l}^{\dagger}|k\rangle_{\rm m}\langle
k|K_{l}|i\rangle,\end{split}$ (16)

where we have simply transposed each matrix element and rearranged the terms.
If all the off-diagonal elements $\langle i|\varrho_{{\rm m}{|}k}^{(\rm
out)}|j\rangle$ of the conditional state are zero, that is, we are unable to
observe any conditional coherences within the monitor, this is equivalent to
stating that $\sum_{l}\langle j|K_{l}^{\dagger}|k\rangle_{\rm m}\langle
k|K_{l}|i\rangle=0$ for all $k$ and all $i\neq j$. This readily means that the
dual channel $\Lambda^{\rm d}[\cdot]=\sum_{l}K_{l}^{\dagger}\cdot K_{l}$ is
incoherent (see Ref. Chitambar and Gour (2016a) or SI-4 for a proof) which
contradicts our assumption. Therefore, at least one of the conditional states
$\varrho_{{\rm m}{|}k}^{(\rm out)}$ must be coherent.

As we see, the detection power of our protocol is independent of the channel
being incoherent, rather it depends on the coherence properties of the dual of
the channel. Now, since (i) the probe mode is incoherent—observe that
$\operatorname{Tr}_{m}|\Phi^{+}\rangle\\!\langle\Phi^{+}|=I/n$, (ii) the
channel $\Lambda$ is also incoherent, and (iii) the conditional input monitor
states $\langle k|\Phi^{+}\rangle\langle\Phi^{+}|k\rangle_{\rm
p}=|k\rangle_{\rm m}\langle k|/n$ are incoherent, we must have that the
observed conditional coherence is due to the initial coherence of the system’s
initial state $\pi_{\rm s}^{(\rm in)}$.

Example.— Let us demonstrate our findings via a simple physical example.
Suppose the system, the probe, and the monitor are two-level systems. We
assume there are two reservoirs inducing decoherence, one of which is coupled
to the system and the other couples to the probe, and investigate two regimes.

First, we consider the case where the probe-system interaction is fast enough
so that system’s decoherence during the interaction is negligible. However, we
assume that probe mode fully dephases due to its interaction with its
reservoir before we measure it. In this case, $\Lambda_{\rm p}[\varrho_{\rm
p}^{(\rm in)}]=\Delta_{\rm p}\circ\Omega_{\rm p}[\varrho_{\rm p}^{(\rm in)}]$
in which $\Omega_{\rm p}[\varrho_{\rm p}^{(\rm in)}]=\operatorname{Tr}_{\rm
s}\Upsilon_{\rm ps}[\varrho_{\rm p}^{(\rm in)}\otimes\pi_{\rm s}^{(\rm in)}]$.

Second, we analyze the case where either of the system or the probe can
decohere during the interaction such that $\Lambda_{\rm p}[\varrho_{\rm
p}^{(\rm in)}]=\operatorname{Tr}_{\rm s,r}\Upsilon_{\rm psr}[\varrho_{\rm
p}^{(\rm in)}\otimes\pi_{\rm s}^{(\rm in)}\otimes\pi_{\rm r}^{(\rm in)}]$.
Here, the subscript $r$ stands for the reservoirs acting on the system and
probe.

Figure 3: Schematic of our example. Two two-level systems, p and m, are
initially prepared in the maximally entangled state $\varrho_{\rm mp}^{(\rm
in)}$. The probe mode exchanges excitation with the system, s, which is
believed to be initially prepared in a superposition of its energy
eigenstates, at the rate $g_{\rm ps}$. While the system and probe may suffer
from dephasings with rates $\gamma_{\rm s}$ and $\gamma_{\rm p}$,
respectively, the monitor is isolated and decoherence-free. Observing
coherence within the state of the monitor mode conditioned on the outcomes of
an incoherent measurement on the probe certifies the coherence of system’s
initial state.

Let $|0\rangle$ and $|1\rangle$ be the eigenstates of an arbitrary degree of
freedom for the system, the probe, and the monitor. We note that the degree of
freedom may be different for each of the three subsystems. The total
Hamiltonian $(\hbar=1)$ is given by

$H=H_{\rm m}+H_{\rm p}+H_{\rm s}+H_{\rm int},$ (17)

in which $H_{\rm j}=\omega_{\rm j}\sigma_{\rm j}^{+}\sigma_{\rm j}^{-}$ is the
bare Hamiltonian of the mode ${\rm j}={\rm m,p,s}$. $\sigma_{\rm j}^{+}$
($\sigma_{\rm j}^{-}$) is the raising (lowering) operator for mode ${\rm j}$
and $\omega_{\rm j}$ is the excitation energy. The probe-system interaction is
described by the Hamiltonian Scully and Zubairy (1999); Raimond and Haroche
(2006)

$H_{\rm int}=g_{\rm ps}(\sigma_{\rm p}^{+}\sigma_{\rm s}^{-}+\sigma_{\rm
p}^{-}\sigma_{\rm s}^{+}),$ (18)

where $g_{\rm ps}$ is the rate of the excitation exchange between the system
and the probe (see Fig. 3). The unitary interaction generated by this
Hamiltonian satisfies the condition in Eq. (7) which allows us to make valid
statements regarding the system’s initial coherence. We further assume that
there are no decay processes in the protocol that transfer the excitation of
mode ${\rm j}$ to the environment, that is the collective dynamics of the
probe-system preserves the total excitation number. This is essential if we
post-select excited states of the output probe. However, the probe-system
dynamics is susceptible to a pure dephasing process that eliminates the phase
coherence of local excitations of mode ${\rm i}={\rm p,s}$ at a rate of
$\gamma_{\rm i}$. The total dynamics of the probe-system can be modelled by a
Born-Markov master equation of the form Gardiner and Zoller (2004); Zhang  _et
al._ (2007); Caruso  _et al._ (2012)

$\frac{d\varrho}{dt}=-i[H,\varrho]+\mathcal{L}_{\rm deph}[\sigma_{z{\rm
s}}]\varrho+\mathcal{L}_{\rm deph}[\sigma_{z{\rm p}}]\varrho,$ (19)

where $\sigma_{z{\rm i}}$ is the Pauli operator along the $z$ axis of mode
${\rm i}={\rm s},{\rm p}$. The action of the super-operator describing the
dephasing processes of mode $\rm i$ is

$\mathcal{L}_{\rm deph}[\sigma_{z{\rm i}}]\varrho=\frac{\gamma_{\rm
i}}{2}(\sigma_{z{\rm i}}\varrho\sigma_{z{\rm i}}-\varrho).$ (20)

First regime.— The system is initially in an arbitrary pure state of the form
$|\phi(0)\rangle_{\rm s}=\epsilon|1\rangle_{\rm
s}+\sqrt{1-\epsilon^{2}}|0\rangle_{\rm s}$ where $\epsilon$ determines its
degree of coherence. We simulate the performance of our protocol in this
scenario in two steps. First, we let the probe interact with the system via
the interaction Hamiltonian (18). Next, we allow the probe to fully dephase
after the interaction, followed by a measurement of the coherence in the
monitor conditioned on finding the probe in the excited state $|1\rangle_{\rm
p}$. The details of the calculations for this simulation are provided in SI-5.

Figure 4: The evolution of quantum coherence within the state of the system in
the first regime. The magnitude of the monitor’s off-diagonal density matrix
element conditioned on detecting a single excitation in the probe is depicted
(a) versus the degree of the system’s initial coherence, $\epsilon$, at
$g_{\rm ps}t=0.75$; (b) versus the normalized detection time, $g_{\rm ps}t$,
for different values of $\epsilon$.

Figure 4 (a) shows the magnitude of the monitor’s off-diagonal element of the
density matrix, $|\mbox{$\langle 1|\varrho_{{\rm m}{|}1}^{(\rm
out)}|0\rangle$}_{\rm m}|$, versus system’s degree of coherence $\epsilon$.
For $\epsilon=0$ and $\epsilon=1$, where the system is in an incoherent state,
the monitor remains in an incoherent state as expected. However, for any other
non-zero $\epsilon$, where the system contains initial coherence, a
measurement of the off-diagonal element of the monitor results in a non-zero
value with its maximum occurring at $\epsilon=1/\sqrt{2}$, i.e. for the
maximally-coherent initial state of the system. Figure 4 (b) shows the
magnitude of the monitor’s off-diagonal element of the density matrix versus
the detection time for different values of $\epsilon$. These results clearly
show that in situations where the probe is prone to dephasing before we
measurement it, our protocol successfully detects the initial coherence within
the system.

The effect of probe’s dephasing before the interaction with the system can be
seen as a reduction in the initial quantum correlations between the probe and
monitor. It is both natural and correct that such a reduction in correlations
reduces the power of the monitor to accumulate the coherence exhibited by the
probe. We provide the detail of this analysis in SI-6.

Figure 5: The evolution of quantum coherence within the state of the system in
the second regime. The magnitude of the monitor’s off-diagonal density matrix
element conditioned on detecting a single excitation in the probe is depicted
(a) versus the normalised probe-system interaction duration, $g_{\rm
ps}\tau_{\rm ps}$, for the cases where the probe alone decoheres during the
interaction, i.e., $\gamma_{\rm p}=0.1$ and $\gamma_{\rm s}=0$ (solid blue
line), and where both the probe and the system decohere while interacting,
i.e., $\gamma_{\rm p}=0.1$ and $\gamma_{\rm s}=0.1$ (dashed red line); (b)
versus the normalized evolution time of the monitor for $\gamma_{\rm
s}=\gamma_{\rm p}=0.1$ and two different values of interaction duration. In
both graphs the probe continues to dephase after $\tau_{\rm ps}$ until the
measurement is performed at an arbitrary time $t$.

Second regime.— Let us now consider the more realistic regime in which the
system and probe decohere during their interaction due to their coupling to
the environment. We are particularly interested in the relative time scales of
the coherent and dephasing evolutions which allow for detecting the initial
coherence of the system.

We consider the same initial states of $|\phi(0)\rangle_{\rm s}$ and
$|\Phi^{+}\rangle_{\rm mp}$ for the system and monitor-probe, respectively. We
let the system and probe interact for the duration $\tau_{\rm ps}$ in which
both system and probe are susceptible to decoherence. The dynamical evolution
of the entire system is simulated by solving the master equation (19) for the
total state $\varrho_{\rm mps}(t)$. Then, we let the probe alone decohere
further for the duration $\tau_{\rm meas}$ until it is being measured.
Finally, we trace over the state of the system which is inaccessible and
calculate the monitor’s off-diagonal element $\mbox{$\langle 1|\varrho_{{\rm
m}{|}1}^{(\rm out)}|0\rangle$}_{\rm m}$, i.e., conditioned on finding the
probe in the excited state $|1\rangle_{\rm p}$ .

In Figure 5 (a), we depict monitor’s coherence as a function of the probe-
system interaction time. It can clearly be seen that, for fast enough probe-
system interactions such that $\tau_{\rm ps}<{\rm min}(1/\gamma_{\rm
p},1/\gamma_{\rm s})$, the initial coherence of the system will be
successfully observed in the monitor even if the probe fully dephases during
its free evolution time $\tau_{\rm meas}$. However, if $\tau_{\rm ps}$ is of
the order of the dephasing time of either system or probe, that is $\tau_{\rm
ps}\approx{\rm min}(1/\gamma_{\rm p},1/\gamma_{\rm s})$, then the monitor will
not retain the coherence information. Figure 5 (b) illustrates the time
evolution of the induced coherence within the monitor for this case.

We observe that decoherence occurring during the probe-system interaction
causes loss of coherence information in the monitor manifested in the
reduction of $|\mbox{$\langle 1|\varrho_{{\rm m}{|}1}^{(\rm
out)}|0\rangle$}_{\rm m}|$. This is better understood if we think of the probe
as an infinitesimally small part of the bath coupled to the system that is
accessible to us. Upon the system’s interaction with its surroundings, its
coherence information spreads to the entire environment. As a result, the
longer the probe-system interaction is, the smaller the share of the probe
from that information will be. Hence, for the monitor to preserve the
information on the initial coherence of the system, the probe-system
interaction time should be shorter than the coherence lifetime of the system.
Counter-intuitively, further decoherence of the probe does not affect the
monitoring of the system’s coherence because, with the help of quantum
correlations between the probe and monitor, a copy of this information is
stored in the monitor, which is decoherence-free.

In conclusion, we have shown the counterintuitive phenomenon of detecting a
quantum system’s initial coherence when the input and output probe states are
completely incoherent. We achieved this through rigorously examining the
elements of coherence detection schemes from a quantum information
perspective. Our analysis yields the necessary and sufficient conditions to
enable valid claims regarding the coherence of a directly inaccessible system.
We provided a protocol for witnessing the quantum coherence of such systems
that satisfies these conditions. This protocol is then used in a proof-of-
principle demonstration of our results, highlighting its power and
limitations. Our results confirm the necessity of fast measurements when it
comes to the duration of the probe’s interaction with the system.

We believe our analysis inspires novel protocols to detect coherent channels
using entangled probes. The use of a quantum entangled twin-mode probe also
opens up an avenue for experimental schemes in which the probe might undergo
decoherence before its information is retrieved. In particular, we think that,
in the near future, entangled light will play a significant role in probing
all sorts of quantum phenomena, including quantum coherence, in complex
systems.

###### Acknowledgements.

The authors gratefully acknowledge Eric Chitambar, Eric Bittner, Martin
Ringbauer, Ivan Kassal, and Eric Gauger for helpful discussions. S.B-E.
acknowledges funding from the European Union Horizon 2020 research and
innovation programme under the Marie Skłodowska-Curie grant agreement No
663830. F.S. was supported by the Royal Commission for the Exhibition of 1851
Research Fellowship. This project was also supported in part through the Sêr
SAM Project at Swansea University, an initiative funded through the Welsh
Government’s Sêr Cymru II Program (European Regional Development Fund).

## References

  * Hayes  _et al._ (2013) D. Hayes, G. B. Griffin, and G. S. Engel,  _Engineering coherence among excited states in synthetic heterodimer systems,_Science 340, 1431 (2013).
  * Scholes  _et al._ (2017) G. D. Scholes, G. R. Fleming, L. X. Chen, A. Aspuru-Guzik, A. Buchleitner, D. F. Coker, G. S. Engel, R. Van Grondelle, A. Ishizaki, D. M. Jonas, _et al._ ,  _Using coherence to enhance function in chemical and biophysical systems,_Nature 543, 647 (2017).
  * Scully (2010) M. O. Scully,  _Quantum photocell: Using quantum coherence to reduce radiative recombination and increase efficiency,_Physical review letters 104, 207701 (2010).
  * Kimble (2008) H. J. Kimble,  _The quantum internet,_Nature 453, 1023 (2008).
  * Streltsov  _et al._ (2017) A. Streltsov, G. Adesso, and M. B. Plenio,  _Colloquium : Quantum coherence as a resource,_Reviews of Modern Physics 89, 041003 (2017).
  * Scholes (2011) G. D. Scholes,  _Coherence in photosynthesis,_Nature Physics 7, 448 (2011).
  * Baghbanzadeh and Kassal (2016) S. Baghbanzadeh and I. Kassal,  _Geometry, supertransfer, and optimality in the light harvesting of purple bacteria,_The Journal of Physical Chemistry Letters 7, 3804 (2016).
  * Scholes (2010) G. D. Scholes,  _Quantum-coherent electronic energy transfer: Did nature think of it first?_The Journal of Physical Chemistry Letters 1, 2 (2010).
  * Chenu and Scholes (2015) A. Chenu and G. D. Scholes,  _Coherence in energy transfer and photosynthesis,_Annual Review of Physical Chemistry 66, 69 (2015).
  * Tomasi  _et al._ (2019) S. Tomasi, S. Baghbanzadeh, S. Rahimi-Keshari, and I. Kassal,  _Coherent and controllable enhancement of light-harvesting efficiency,_Physical Review A 100, 043411 (2019).
  * Tomasi and Kassal (2020) S. Tomasi and I. Kassal,  _Classification of coherent enhancements of light-harvesting processes,_The Journal of Physical Chemistry Letters 11, 2348 (2020).
  * Fruchtman  _et al._ (2016) A. Fruchtman, R. Gómez-Bombarelli, B. W. Lovett, and E. M. Gauger,  _Photocell optimization using dark state protection,_Physical review letters 117, 203603 (2016).
  * Tomasi  _et al._ (2021) S. Tomasi, D. M. Rouse, E. M. Gauger, B. W. Lovett, and I. Kassal,  _Environmentally improved coherent light harvesting,_The Journal of Physical Chemistry Letters 12, 6143 (2021).
  * Creatore  _et al._ (2013) C. Creatore, M. A. Parker, S. Emmott, and A. W. Chin,  _Efficient biologically inspired photocell enhanced by delocalized quantum states,_Physical review letters 111, 253601 (2013).
  * Engel  _et al._ (2007) G. S. Engel, T. R. Calhoun, E. L. Read, T.-K. Ahn, T. Mančal, Y.-C. Cheng, R. E. Blankenship, and G. R. Fleming,  _Evidence for wavelike energy transfer through quantum coherence in photosynthetic systems,_Nature 446, 782 (2007).
  * Collini  _et al._ (2010) E. Collini, C. Y. Wong, K. E. Wilk, P. M. Curmi, P. Brumer, and G. D. Scholes,  _Coherently wired light-harvesting in photosynthetic marine algae at ambient temperature,_Nature 463, 644 (2010).
  * Panitchayangkoon  _et al._ (2011) G. Panitchayangkoon, D. V. Voronine, D. Abramavicius, J. R. Caram, N. H. C. Lewis, S. Mukamel, and G. S. Engel,  _Direct evidence of quantum transport in photosynthetic light-harvesting complexes,_Proceedings of the National Academy of Sciences 108, 20908 (2011).
  * Abramavicius and Mukamel (2010) D. Abramavicius and S. Mukamel,  _Quantum oscillatory exciton migration in photosynthetic reaction centers,_The Journal of chemical physics 133, 08B603 (2010).
  * Panitchayangkoon  _et al._ (2010) G. Panitchayangkoon, D. Hayes, K. A. Fransted, J. R. Caram, E. Harel, J. Wen, R. E. Blankenship, and G. S. Engel,  _Long-lived quantum coherence in photosynthetic complexes at physiological temperature,_Proceedings of the National Academy of Sciences 107, 12766 (2010).
  * Kassal  _et al._ (2013) I. Kassal, J. Yuen-Zhou, and S. Rahimi-Keshari,  _Does coherence enhance transport in photosynthesis?_The journal of physical chemistry letters 4, 362 (2013).
  * Wang  _et al._ (2019) L. Wang, M. A. Allodi, and G. S. Engel,  _Quantum coherences reveal excited-state dynamics in biophysical systems,_Nature Reviews Chemistry 3, 477 (2019).
  * Brumer and Shapiro (2012) P. Brumer and M. Shapiro,  _Molecular response in one-photon absorption via natural thermal light vs. pulsed laser excitation,_Proceedings of the National Academy of Sciences 109, 19575 (2012).
  * Mančal and Valkunas (2010) T. Mančal and L. Valkunas,  _Exciton dynamics in photosynthetic complexes: excitation by coherent and incoherent light,_New Journal of Physics 12, 065044 (2010).
  * Jiang and Brumer (1991) X.-P. Jiang and P. Brumer,  _Creation and dynamics of molecular states prepared with coherent vs partially coherent pulsed light,_The Journal of chemical physics 94, 5833 (1991).
  * Muñoz and Schlawin (2020) C. S. Muñoz and F. Schlawin,  _Photon correlation spectroscopy as a witness for quantum coherence,_Physical Review Letters 124, 203601 (2020).
  * Chitambar and Gour (2016a) E. Chitambar and G. Gour,  _Comparison of incoherent operations and measures of coherence,_Phys. Rev. A 94, 052336 (2016a).
  * Chitambar and Gour (2016b) E. Chitambar and G. Gour,  _Critical Examination of Incoherent Operations and a Physically Consistent Resource Theory of Quantum Coherence,_Physical Review Letters 117, 030401 (2016b).
  * Nielsen and Chuang (2010) M. A. Nielsen and I. L. Chuang,  _Quantum Computation and Quantum Information: 10th Anniversary Edition_ (Cambridge University Press, 2010).
  * Scully and Zubairy (1999) M. O. Scully and M. S. Zubairy,  _Quantum optics_ (Cambridge University Press, 1999).
  * Raimond and Haroche (2006) J.-M. Raimond and S. Haroche,  _Exploring the Quantum: Atoms, Cavities, and Photons_ (Oxford University Press, 2006).
  * Gardiner and Zoller (2004) C. Gardiner and P. Zoller,  _Quantum noise: a handbook of Markovian and non-Markovian quantum stochastic methods with applications to quantum optics_ (Springer-Verlag Berlin Heidelberg, 2004).
  * Zhang  _et al._ (2007) M. Zhang, H. Y. Dai, Z. Xi, H. W. Xie, and D. Hu,  _Combating dephasing decoherence by periodically performing tracking control and projective measurement,_Physical Review A 76, 042335 (2007).
  * Caruso  _et al._ (2012) F. Caruso, S. Montangero, T. Calarco, S. F. Huelga, and M. B. Plenio,  _Coherent optimal control of photosynthetic molecules,_Physical Review A 85, 042331 (2012).


# The pollution to the $K\pi$-puzzle from the isospin-breaking
$\pi^{0}\\!\\!-\\!\eta\\!-\\!\eta^{\prime}$ mixing effect

Zhen-Hua Zhang  College of Nuclear Science and Technology, University of South
China, Hengyang, 421001, Hunan, China

###### Abstract

The influence of the isospin-breaking
$\pi^{0}\\!\\!-\\!\eta\\!-\\!\eta^{\prime}$ mixing effect on the
$CP$-asymmetries of $B\to K\pi$ processes is examined for the first time. It
is found that this mixing effect brings a large uncertainty both to the
$CP$-asymmetry sum rule of $B\to K\pi$ processes and the $CP$-asymmetry
difference of the $B^{+}\to K^{+}\pi^{0}$ and $B^{0}\to K^{+}\pi^{-}$,
obscuring the significance of the $K\pi$-puzzle. This uncertainty can be so
large that it is even possible to explain the $K\pi$-puzzle by the
$\pi^{0}\\!\\!-\\!\eta\\!-\\!\eta^{\prime}$ mixing effect alone.

##  I Brief introduction to the $K\pi$-puzzle

It is believed that the decay processes $B\to K\pi$ are good probes for new
physics (NP) beyond the Standard Model (SM), as the tree-level amplitudes are
suppressed, making them more sensitive to the potentially NP contributions.
Based on the isospin consideration, the amplitudes of the aforementioned weak
decays $B\to K\pi$ are related through Nir and Quinn (1991); Lipkin  _et al._
(1991); Gronau (1991)

$\mathcal{A}_{B^{+}\to K^{0}\pi^{+}}+\sqrt{2}\mathcal{A}_{B^{+}\to
K^{+}\pi^{0}}=\mathcal{A}_{{B^{0}}\to
K^{+}\pi^{-}}+\sqrt{2}\mathcal{A}_{{B^{0}}\to K^{0}\pi^{0}},$ (1)

from which a sum rule between $CP$ asymmetries of $B\to K\pi$ processes can be
deduced Gronau (2005):

$A_{CP}^{K^{+}\pi^{-}}+\frac{\mathcal{B}(K^{0}\pi^{+})}{\mathcal{B}(K^{+}\pi^{-})}\frac{\tau_{0}}{\tau_{+}}A_{CP}^{K^{0}\pi^{+}}=\frac{2\mathcal{B}(K^{+}\pi^{0})}{\mathcal{B}(K^{+}\pi^{-})}\frac{\tau_{0}}{\tau_{+}}A_{CP}^{K^{+}\pi^{0}}+\frac{2\mathcal{B}(K^{0}\pi^{0})}{\mathcal{B}(K^{+}\pi^{-})}\frac{\tau_{0}}{\tau_{+}}A_{CP}^{K^{0}\pi^{0}}.$
(2)

This $CP$-asymmetry sum rule can be further simplified into a more crude
relation between the $CP$ asymmetries of $B^{+}\to K^{+}\pi^{0}$ and $B^{0}\to
K^{+}\pi^{-}$ Gronau and Rosner (1999):

$A_{CP}^{K^{+}\pi^{0}}-A_{CP}^{K^{+}\pi^{-}}\approx 0,$ (3)

which is clearly in contradiction to the latest world average of the
$CP$-asymmetry difference between these two aforementioned processes Duh  _et
al._ (2013); Lees  _et al._ (2013); Aaltonen  _et al._ (2014); Aaij  _et al._
(2018); Amhis  _et al._ (2021); Aaij  _et al._ (2021),

$\Delta A_{CP}(K\pi)\equiv
A_{CP}^{K^{+}\pi^{0}}-A_{CP}^{K^{+}\pi^{-}}=0.115\pm 0.014.$ (4)

This is basically a short version of the long-standing “$K\pi$-puzzle”.

The $K\pi$-puzzle has gained constant theoretical attentions, especially those
works on the exploration of the possibility of NP Buras  _et al._ (2003,
2004a, 2004b); Barger  _et al._ (2004); Baek  _et al._ (2005); Chang  _et al._
(2009); Fleischer  _et al._ (2018); Beaudry  _et al._ (2018); Bhattacharya
_et al._ (2021). Meanwhile, there are also works attempting to understand it
within SM Bauer  _et al._ (2006); Cheng and Chua (2009); Cheng  _et al._
(2015); Zhou  _et al._ (2017); Kundu  _et al._ (2021). In this paper, we will
try to examine the role of the isospin breaking effect –the
$\pi^{0}\\!\\!-\\!\eta\\!-\\!\eta^{\prime}$ mixing effect– in the
$K\pi$-puzzle.

##  II the pollution of the $\pi^{0}\\!\\!-\\!\eta\\!-\\!\eta^{\prime}$ mixing
effects to the $K\pi$-puzzle

The basic idea is very simple. Since there are $\pi^{0}$s in the final states
of these $B\to K\pi$ decay processes, the isospin-breaking
$\pi^{0}\\!\\!-\\!\eta\\!-\\!\eta^{\prime}$ mixing effect Gross  _et al._
(1979); Leutwyler (1996); Feldmann  _et al._ (1999); Escribano and Royo
(2020); Kordov  _et al._ (2021) takes place. Although this effect seems to be
negligible at first sight as it is small Gronau  _et al._ (2006) –perhaps that
is why the $\pi^{0}\\!\\!-\\!\eta\\!-\\!\eta^{\prime}$ mixing effect has never
been put on the table dealing with the $K\pi$-puzzle– this kind of isospin-
breaking effects could potentially affect $CP$ asymmetries more badly than
expected Zhang  _et al._ (2020).

Within the context of the $\pi^{0}\\!\\!-\\!\eta\\!-\\!\eta^{\prime}$ mixing
effect, the $\pi^{0}$ meson can be expressed as an admixture of the isospin
eigenstate $\pi_{3}$ and the mass eigenstates $\eta$ and $\eta^{\prime}$,
$|\pi^{0}\rangle=|\pi_{3}\rangle+\epsilon|\eta\rangle+\epsilon^{\prime}|\eta^{\prime}\rangle$,
where $\epsilon$ and $\epsilon^{\prime}$ are small parameters which account
for the mixing between $\pi^{0}$ and $\eta$ and $\eta^{\prime}$, respectively
Leutwyler (1996); Feldmann  _et al._ (1999). The $CP$ asymmetry parameters of
$B^{+}\to K^{+}\pi^{0}$ and $B^{0}\to K^{0}\pi^{0}$ acquire corrections
because of the $\pi^{0}\\!\\!-\\!\eta\\!-\\!\eta^{\prime}$ mixing. For
example, up to $\mathcal{O}({\epsilon^{(^{\prime})}})$, the $CP$ asymmetry of
$B^{+}\to K^{+}\pi^{0}$ is modified into 111An extra factor
$\left(1-{A_{CP}^{K^{+}\pi_{3}}}^{2}\right)$ in front of
$\Delta_{\text{IB}}^{K^{+}\pi^{0}}$ is omitted as it is numerically very close
to one.

$A_{CP}^{K^{+}\pi^{0}}=A_{CP}^{K^{+}\pi_{3}}+\Delta_{\text{IB}}^{K^{+}\pi^{0}},$
(5)

where $A_{CP}^{K^{+}\pi_{3}}$ is in fact the conventionally defined $CP$
asymmetry of $B^{+}\to K^{+}\pi^{0}$ without the
$\pi^{0}\\!\\!-\\!\eta\\!-\\!\eta^{\prime}$ mixing effect, which takes the
form $A_{CP}^{K^{+}\pi_{3}}\equiv\frac{\left|\mathcal{A}_{B^{-}\to
K^{-}\pi_{3}}\right|^{2}-\left|\mathcal{A}_{B^{+}\to
K^{+}\pi_{3}}\right|^{2}}{\left|\mathcal{A}_{B^{-}\to
K^{-}\pi_{3}}\right|^{2}+\left|\mathcal{A}_{B^{+}\to
K^{+}\pi_{3}}\right|^{2}}$, and $\Delta_{\text{IB}}^{K^{+}\pi^{0}}$ represents
the correction of the $\pi^{0}\\!\\!-\\!\eta\\!-\\!\eta^{\prime}$ mixing
effect to the $CP$ asymmetry of $B^{+}\to K^{+}\pi^{0}$, which, after some
algebra, can be expressed as

$\Delta_{\text{IB}}^{K^{+}\pi^{0}}=\Re\left[\left(\frac{\mathcal{A}_{B^{-}\to
K^{-}\eta}}{\mathcal{A}_{B^{-}\to K^{-}\pi^{0}}}-\frac{\mathcal{A}_{B^{+}\to
K^{+}\eta}}{\mathcal{A}_{B^{+}\to
K^{+}\pi^{0}}}\right)e^{i\theta}\epsilon+\left(\frac{\mathcal{A}_{B^{-}\to
K^{-}\eta^{\prime}}}{\mathcal{A}_{B^{-}\to
K^{-}\pi^{0}}}-\frac{\mathcal{A}_{B^{+}\to
K^{+}\eta^{\prime}}}{\mathcal{A}_{B^{+}\to
K^{+}\pi^{0}}}\right)e^{i\theta^{\prime}}\epsilon^{\prime}\right].$ (6)

Note that two strong phases between the amplitudes of $B^{\pm}\to
K^{\pm}\eta^{(\prime)}$ and $B^{\pm}\to K^{\pm}\pi^{0}$, $\theta$ and
$\theta^{\prime}$, are presented explicitly. From Eq. (5) and (6) one can
roughly see that the $\pi^{0}\\!\\!-\\!\eta\\!-\\!\eta^{\prime}$ mixing
correction term can not be simply neglected, provided that the $CP$ asymmetry
parameter $A_{CP}^{K^{+}\pi^{0}}$ is about the same order with the mixing
parameters $\epsilon$ and/or $\epsilon^{\prime}$. Similarly, there is also a
correction term to the $CP$ asymmetry of $B^{0}\to K^{0}\pi^{0}$, which reads

$\Delta_{\text{IB}}^{K^{0}\pi^{0}}=\Re\left[\left(\frac{\mathcal{A}_{\overline{B^{0}}\to\overline{K^{0}}\eta}}{\mathcal{A}_{\overline{B^{0}}\to\overline{K^{0}}\pi^{0}}}-\frac{\mathcal{A}_{B^{0}\to
K^{0}\eta}}{\mathcal{A}_{B^{0}\to
K^{0}\pi^{0}}}\right)e^{i\tilde{\theta}}\epsilon+\left(\frac{\mathcal{A}_{\overline{B^{0}}\to\overline{K^{0}}\eta^{\prime}}}{\mathcal{A}_{\overline{B^{0}}\to\overline{K^{0}}\pi^{0}}}-\frac{\mathcal{A}_{B^{0}\to
K^{0}\eta^{\prime}}}{\mathcal{A}_{B^{0}\to
K^{0}\pi^{0}}}\right)e^{i\tilde{\theta}^{\prime}}\epsilon^{\prime}\right].$
(7)

Since Eqs. (2) and (3) are obtained under the ignorance of the the
$\pi^{0}\\!\\!-\\!\eta\\!-\\!\eta^{\prime}$ mixing effect, the $\pi^{0}$s in
these two equations are in fact the isospin eigenstates $\pi_{3}$s.
Consequently, Eqs. (2) and (3) should be rewritten as

$A_{CP}^{K^{+}\pi^{-}}+\frac{\mathcal{B}(K^{0}\pi^{+})}{\mathcal{B}(K^{+}\pi^{-})}\frac{\tau_{0}}{\tau_{+}}A_{CP}^{K^{0}\pi^{+}}=\frac{2\mathcal{B}(K^{+}\pi_{3})}{\mathcal{B}(K^{+}\pi^{-})}\frac{\tau_{0}}{\tau_{+}}A_{CP}^{K^{+}\pi_{3}}+\frac{2\mathcal{B}(K^{0}\pi_{3})}{\mathcal{B}(K^{+}\pi^{-})}\frac{\tau_{0}}{\tau_{+}}A_{CP}^{K^{0}\pi_{3}},$
(8)

and

$A_{CP}^{K^{+}\pi_{3}}-A_{CP}^{K^{+}\pi^{-}}\approx 0,$ (9)

respectively. It should be pointed out that all the branching ratios and $CP$
asymmetries containing $\pi_{3}$ in the final state are strictly speaking not
physical observables. This however, can be fixed with the aid of Eq. (5) by
rewriting Eqs. (8) and (9) in terms of the physical branching ratios and $CP$
asymmetries with $\pi^{0}$s containing in the final states. The $CP$-asymmetry
sum rule and the $CP$-asymmetry difference now read

$A_{CP}^{K^{+}\pi^{-}}+\frac{\mathcal{B}(K^{0}\pi^{+})}{\mathcal{B}(K^{+}\pi^{-})}\frac{\tau_{0}}{\tau_{+}}A_{CP}^{K^{0}\pi^{+}}=\frac{2\mathcal{B}(K^{+}\pi^{0})}{\mathcal{B}(K^{+}\pi^{-})}\frac{\tau_{0}}{\tau_{+}}A_{CP}^{K^{+}\pi^{0}}+\frac{2\mathcal{B}(K^{0}\pi^{0})}{\mathcal{B}(K^{+}\pi^{-})}\frac{\tau_{0}}{\tau_{+}}A_{CP}^{K^{0}\pi^{0}}-\Delta_{\text{IB}},$
(10)

and

$A_{CP}^{K^{+}\pi^{0}}-A_{CP}^{K^{+}\pi^{-}}\approx\Delta_{\text{IB}},$ (11)

respectively, where $\Delta_{\text{IB}}$ accommodates the
$\pi^{0}\\!\\!-\\!\eta\\!-\\!\eta^{\prime}$ mixing correction and takes the
form 222Contributions from the branching ratios, which are proportional to
$A_{CP}^{K^{+,0}\pi^{0}}\epsilon^{(\prime)}$, have been neglected in
$\Delta_{\text{IB}}$.

$\displaystyle\Delta_{\text{IB}}=\frac{2\mathcal{B}(K^{+}\pi^{0})}{\mathcal{B}(K^{+}\pi^{-})}\frac{\tau_{0}}{\tau_{+}}\Delta_{\text{IB}}^{K^{+}\pi^{0}}+\frac{2\mathcal{B}(K^{0}\pi^{0})}{\mathcal{B}(K^{+}\pi^{-})}\frac{\tau_{0}}{\tau_{+}}\Delta_{\text{IB}}^{K^{0}\pi^{0}}.$
(12)

One interesting behaviour of this modification is that although Eq. (3)
relates only the $CP$ asymmetries of $B^{0}\to K^{+}\pi^{-}$ and $B^{+}\to
K^{+}\pi^{0}$, the isospin-breaking correction term $\Delta_{\text{IB}}$ in
Eq. (11), however, not only contains the contribution of the process $B^{+}\to
K^{+}\pi^{0}$, but also contains that of the process $B^{0}\to K^{0}\pi^{0}$.
The latter turns out to be numerically even more important in
$\Delta_{\text{IB}}$.

We are now in a position to estimate the impact of the the isospin-breaking
correction term $\Delta_{\text{IB}}$ to the $K\pi$-puzzle. Although the
amplitudes in Eqs. (6) and (7) can be calculated theoretically or extracted
from data, the four relative strong phases between these amplitudes, $\theta$,
$\theta^{\prime}$, $\tilde{\theta}$, and $\tilde{\theta}^{\prime}$, are non-
perturbative, preventing us from an accurate prediction of
$\Delta_{\text{IB}}$. What we can do is to give an rough estimation of the
possible range of $\Delta_{\text{IB}}$ by treating the four strong phases as
free parameters varying from 0 to $2\pi$ independently. Based on this
strategy, with the amplitudes borrowed from Ref. Cheng  _et al._ (2015),
$\Delta_{\text{IB}}$ is estimated to be

$\Delta_{\text{IB}}=(-0.37,+0.37)\times\epsilon+(-1.16,+1.16)\times\epsilon^{\prime}.$
(13)

In TABLE 1, with different values of $\epsilon$ and $\epsilon^{\prime}$ from
different references, the corresponding allowed range of $\Delta_{\text{IB}}$
is calculated via Eq. (13). Note that the mixing parameters $\epsilon$ and
$\epsilon^{\prime}$ take quite different values throughout the literature,
resulting in quite different ranges of the $\Delta_{\text{IB}}$, as is
presented in the last column of this table. From TABLE 1 one can see that the
isospin-breaking-correction term $\Delta_{\text{IB}}$ can be as large as a few
percent, which indicates that the influence of the
$\pi^{0}\\!\\!-\\!\eta\\!-\\!\eta^{\prime}$ mixing effect to the
$CP$-asymmetry sum rule of $B\to K\pi$ and the $CP$-asymmetry difference of
the $B^{+}\to K^{+}\pi^{0}$ and $B^{0}\to K^{+}\pi^{-}$ can not be simply
ignored. It is even possible that the $K\pi$-puzzle can be explained by the
$\pi^{0}\\!\\!-\\!\eta\\!-\\!\eta^{\prime}$ mixing effect alone. Take the
$CP$-asymmetry difference between $B^{+}\to K^{+}\pi^{0}$ and $B^{0}\to
K^{+}\pi^{-}$ as example. If we treat the $\Delta_{\text{IB}}$-term as an
uncertainty, with the mixing parameters from Ref. Escribano and Royo (2020),
the $CP$-asymmetry difference now should read as

$A_{CP}^{K^{+}\pi_{3}}-A_{CP}^{K^{+}\pi^{-}}=0.115\pm 0.014\pm
0.041_{\pi^{0}\\!\\!-\\!\eta\\!-\\!\eta^{\prime}}=0.115\pm 0.043,$ (14)

from which one can clearly see that the uncertainty caused by the
$\pi^{0}\\!\\!-\\!\eta\\!-\\!\eta^{\prime}$ mixing effect can be even larger
than the other experimental uncertainties combined. Hence the significance of
the nonzeroness of this $CP$ difference is considerably reduced to less that 3
standard deviations.

Table 1: The $\pi^{0}\\!\\!-\\!\eta\\!-\\!\eta^{\prime}$ mixing parameters quoted from different references and the corresponding range of $\Delta_{\text{IB}}$. Only the central values of $\epsilon$ and $\epsilon^{\prime}$ are used when obtaining the range of $\Delta_{\text{IB}}$. reference | $\epsilon(\%)$ | $\epsilon^{\prime}(\%)$ |  $\Delta_{\text{IB}}(\%$)
---|---|---|---
Kroll Kroll (2005) | $1.7\pm 0.2$ | $0.4\pm 0.1$ | $(-1.1,+1.1)$
Escribano & RoyoEscribano and Royo (2020) | $0.1\pm 0.9$ | $3.5\pm 0.9$ | $(-4.1,+4.1)$
Benayoun et. al. Benayoun  _et al._ (2021) | $4.16\pm 0.20$ | $1.05\pm 0.05$ | $(-2.8,+2.8)$

##  III Conclusion

In conclusion, the contribution of the isospin-breaking
$\pi^{0}\\!\\!-\\!\eta\\!-\\!\eta^{\prime}$ mixing effect to the $CP$
asymmetries of $B\to K\pi$ is investigated for the first time in this paper.
It is found that the $\pi^{0}\\!\\!-\\!\eta\\!-\\!\eta^{\prime}$ mixing effect
pollutes the $K\pi$-puzzle more dramatically than the naive expectation. This
pollution is imbedded in the parameter $\Delta_{\text{IB}}$, which brings
quite a large uncertainty to the $CP$-asymmetry sum rule of $B\to K\pi$ and
the $CP$-asymmetry difference of the $B^{+}\to K^{+}\pi^{0}$ and $B^{0}\to
K^{+}\pi^{-}$. The analysis of this paper shows that it is even possible to
explain the $K\pi$-puzzle by the $\pi^{0}\\!\\!-\\!\eta\\!-\\!\eta^{\prime}$
mixing effect alone. Consequently, the implications of the $K\pi$-puzzle
should be reconsidered.

###### Acknowledgements.

This work was supported by National Natural Science Foundation of China under
Contracts Nos. 11705081 and 12192261.

## References

  * Nir and Quinn (1991) Y. Nir and H. R. Quinn, Phys. Rev. Lett. 67, 541 (1991).
  * Lipkin  _et al._ (1991) H. J. Lipkin, Y. Nir, H. R. Quinn, and A. Snyder, Phys. Rev. D 44, 1454 (1991).
  * Gronau (1991) M. Gronau, Phys. Lett. B 265, 389 (1991).
  * Gronau (2005) M. Gronau, Phys. Lett. B 627, 82 (2005), arXiv:hep-ph/0508047 .
  * Gronau and Rosner (1999) M. Gronau and J. L. Rosner, Phys. Rev. D 59, 113002 (1999), arXiv:hep-ph/9809384 .
  * Duh  _et al._ (2013) Y. T. Duh _et al._ (Belle), Phys. Rev. D 87, 031103 (2013), arXiv:1210.1348 [hep-ex] .
  * Lees  _et al._ (2013) J. P. Lees _et al._ (BaBar), Phys. Rev. D 87, 052009 (2013), arXiv:1206.3525 [hep-ex] .
  * Aaltonen  _et al._ (2014) T. A. Aaltonen _et al._ (CDF), Phys. Rev. Lett. 113, 242001 (2014), arXiv:1403.5586 [hep-ex] .
  * Aaij  _et al._ (2018) R. Aaij _et al._ (LHCb), Phys. Rev. D 98, 032004 (2018), arXiv:1805.06759 [hep-ex] .
  * Amhis  _et al._ (2021) Y. S. Amhis _et al._ (HFLAV), Eur. Phys. J. C 81, 226 (2021), arXiv:1909.12524 [hep-ex] .
  * Aaij  _et al._ (2021) R. Aaij _et al._ (LHCb), Phys. Rev. Lett. 126, 091802 (2021), arXiv:2012.12789 [hep-ex] .
  * Buras  _et al._ (2003) A. J. Buras, R. Fleischer, S. Recksiegel, and F. Schwab, Eur. Phys. J. C 32, 45 (2003), arXiv:hep-ph/0309012 .
  * Buras  _et al._ (2004a) A. J. Buras, R. Fleischer, S. Recksiegel, and F. Schwab, Phys. Rev. Lett. 92, 101804 (2004a), arXiv:hep-ph/0312259 .
  * Buras  _et al._ (2004b) A. J. Buras, R. Fleischer, S. Recksiegel, and F. Schwab, Nucl. Phys. B 697, 133 (2004b), arXiv:hep-ph/0402112 .
  * Barger  _et al._ (2004) V. Barger, C.-W. Chiang, P. Langacker, and H.-S. Lee, Phys. Lett. B 598, 218 (2004), arXiv:hep-ph/0406126 .
  * Baek  _et al._ (2005) S. Baek, P. Hamel, D. London, A. Datta, and D. A. Suprun, Phys. Rev. D 71, 057502 (2005), arXiv:hep-ph/0412086 .
  * Chang  _et al._ (2009) Q. Chang, X.-Q. Li, and Y.-D. Yang, JHEP 05, 056 (2009), arXiv:0903.0275 [hep-ph] .
  * Fleischer  _et al._ (2018) R. Fleischer, R. Jaarsma, and K. K. Vos, Phys. Lett. B 785, 525 (2018), arXiv:1712.02323 [hep-ph] .
  * Beaudry  _et al._ (2018) N. B. Beaudry, A. Datta, D. London, A. Rashed, and J.-S. Roux, JHEP 01, 074 (2018), arXiv:1709.07142 [hep-ph] .
  * Bhattacharya  _et al._ (2021) B. Bhattacharya, A. Datta, D. Marfatia, S. Nandi, and J. Waite, Phys. Rev. D 104, 051701 (2021), arXiv:2104.03947 [hep-ph] .
  * Bauer  _et al._ (2006) C. W. Bauer, I. Z. Rothstein, and I. W. Stewart, Phys. Rev. D 74, 034010 (2006), arXiv:hep-ph/0510241 .
  * Cheng and Chua (2009) H.-Y. Cheng and C.-K. Chua, Phys. Rev. D 80, 114008 (2009), arXiv:0909.5229 [hep-ph] .
  * Cheng  _et al._ (2015) H.-Y. Cheng, C.-W. Chiang, and A.-L. Kuo, Phys. Rev. D 91, 014011 (2015), arXiv:1409.5026 [hep-ph] .
  * Zhou  _et al._ (2017) S.-H. Zhou, Q.-A. Zhang, W.-R. Lyu, and C.-D. Lü, Eur. Phys. J. C 77, 125 (2017), arXiv:1608.02819 [hep-ph] .
  * Kundu  _et al._ (2021) A. Kundu, S. K. Patra, and S. Roy, Phys. Rev. D 104, 095025 (2021), arXiv:2106.15633 [hep-ph] .
  * Gross  _et al._ (1979) D. J. Gross, S. B. Treiman, and F. Wilczek, Phys. Rev. D 19, 2188 (1979).
  * Leutwyler (1996) H. Leutwyler, Phys. Lett. B 374, 181 (1996), arXiv:hep-ph/9601236 .
  * Feldmann  _et al._ (1999) T. Feldmann, P. Kroll, and B. Stech, Phys. Lett. B 449, 339 (1999), arXiv:hep-ph/9812269 .
  * Escribano and Royo (2020) R. Escribano and E. Royo, Phys. Lett. B 807, 135534 (2020), arXiv:2003.08379 [hep-ph] .
  * Kordov  _et al._ (2021) Z. R. Kordov _et al._ , (2021), arXiv:2110.11533 [hep-lat] .
  * Gronau  _et al._ (2006) M. Gronau, Y. Grossman, G. Raz, and J. L. Rosner, Phys. Lett. B 635, 207 (2006), arXiv:hep-ph/0601129 .
  * Zhang  _et al._ (2020) Z.-H. Zhang, S. Yang, and X.-H. Guo, JHEP 10, 020 (2020), arXiv:2005.09157 [hep-ph] .
  * Note (1) An extra factor $\left(1-{A_{CP}^{K^{+}\pi_{3}}}^{2}\right)$ in front of $\Delta_{\REV@text{IB}}^{K^{+}\pi^{0}}$ is omitted as it is numerically very close to one.
  * Note (2) Contributions from the branching ratios, which are proportional to $A_{CP}^{K^{+,0}\pi^{0}}\epsilon^{(\prime)}$, have been neglected in $\Delta_{\REV@text{IB}}$.
  * Kroll (2005) P. Kroll, Mod. Phys. Lett. A 20, 2667 (2005), arXiv:hep-ph/0509031 .
  * Benayoun  _et al._ (2021) M. Benayoun, L. DelBuono, and F. Jegerlehner, (2021), arXiv:2105.13018 [hep-ph] .


# Statistical Characterization of a 1D Random Potential Problem – with
applications in score statistics of MS-based peptide sequencing

Gelio Alves and Yi-Kuo Yu111To whom correspondence should be addressed. E-mail
address<EMAIL_ADDRESS>National Center for Biotechnology Information,
National Library of Medicine, National Institutes of Health, Bethesda, MD
20894

###### Abstract

We provide a complete thermodynamic solution of a 1D hopping model in the
presence of a random potential by obtaining the density of states. Since the
partition function is related to the density of states by a Laplace transform,
the density of states determines completely the thermodynamic behavior of the
system. We have also shown that the transfer matrix technique, or the so-
called dynamic programming, used to obtain the density of states in the 1D
hopping model may be generalized to tackle a long-standing problem in
statistical significance assessment for one of the most important proteomic
tasks – peptide sequencing using tandem mass spectrometry data.

###### keywords:

Statistical Significance , Dynamic Programming , Mass Spectrometry , Directed
Paths in Random Media , Peptide Identification

##  1 Introduction

Important in both fundamental science and numerous applications, optimization
problems of various degrees of complexity are challenging (see [1] for an
excellent introduction). Optimization conditioned by constraints that may vary
from event to event is of especial theoretical and practical importance. As a
first example, when dealing with a system under a random potential, each
realization of the random potential demands a separate optimization resulting
in a different ground state. The thermodynamic behavior of such a system in a
quenched random potential crucially depends on the random potential realized.
A similar but practical problem may arise in routing passengers at various
cities to reach their destinations. In the latter case, the optimal routing
depends on the number of passengers at various locations, the costs from one
location to the others, which likely to vary from time to time. This type of
conditional optimization also occurs in modern proteomics problem, that is, in
the mass spectrometry (MS) based peptide sequencing. In this case, each tandem
MS (MS2) spectrum constitute a different condition for optimization which aims
to find a database peptide or a de novo peptide to best explain the given MS2
spectrum.

When the cost function of an optimization problem can be expressed as a sum of
independent local contributions, the problem usually can be solved using the
transfer matrix method that is commonly employed in statistical physics. A
well-studied example of this sort in statistical physics is the directed
polymer/path in a random medium (DPRM) [2, 3, 4]. Even when a small non-local
energetics is involved, the transfer matrix approach still proves useful [5].
As an example, the close relationship between the DPRM problem and MS-based
peptide sequencing, where a small nonlocal energetics is necessary to enhance
the peptide identifications, was sketched in an earlier publication [5] and
the cost value distribution from many possible solutions other than the
optimal one is explored. Indeed, obtaining the cost value distribution from
all possible solutions in many cases is harder than finding the optimal
solution alone. In this paper, we will provide the solution to a generic
problem that enables a full characterization of the peptide sequencing score
statistics, instead of just the optimal peptide. The 1D problem considered is
essentially a hopping model in the presence of a random potential. The
solution to this problem may also be useful in other applications such as in
routing of passengers and even internet traffic.

In what follows, we will first introduce the generic 1D hopping model in a
random potential, followed by its transfer matrix (or dynamic programming)
solution. We then discuss the utility of this solution in the context of MS-
based peptide sequencing, and demonstrate with real example from mass spectrum
in real MS-based proteomics experiments. In the discussion section, we will
sketch the utility of the transfer matrix solution in other context and then
conclude with a few relevant remarks.

##  2 1D hopping in random potential

Along the $x$-axis, let us consider a particle that can hop with a set of
prescribed distances $\\{m_{i}\\}_{i=1}^{K}$ towards the positive $\hat{x}$
direction. That is, if the particle is currently at location $x_{0}$, it can
move to location $x_{0}+m_{1}$, $x_{0}+m_{2}$, …$x_{0}+m_{K}$ in the next time
step. At each hopping step, the particle will accumulate an energy $-s(x)$
from location $x$ that it just visited. The score $s(x)$ (negative of the on-
site potential energy) is assumed positive and may only exist at a limited
number of locations. For locations that $s(x)$ do not exist, we simply set
$s(x)=0$ there. The energy of a path starting from the origin specified by the
sequential hopping events $p\equiv\\{m_{h_{1}},m_{h_{2}},\ldots,m_{h_{L}}\\}$
would have visited locations $\\{x_{1},x_{2},\ldots,x_{L}\\}$ with
$x_{i}\equiv\sum_{j=1}^{i}m_{h_{j}}$ and has energy

$E_{p}(x=x_{L})\equiv-\sum_{i=1}^{L-1}s\left(x_{i}\right)\equiv-S_{p}(x)\;.$

In general, there can be more than one path terminated at the same point.
Treating each path as a state with energy given by $E_{p}$, one ends up having
the following recursion relation for the partition function
$Z(x)\equiv\sum_{p}e^{-\beta E_{p}(x)}$

$Z(x)=\sum_{i=1}^{K}e^{\beta s(x-m_{i})}Z(x-m_{i})\;,$ (1)

where $\beta=1/T$ plays the role of inverse temperature (with $k_{B}=1$
chosen). If one were only interested in the best score terminated at point
$x$, it will be given by the zero temperature limit $\beta\to\infty$ and the
recursion relation may be obtained by taking the logarithm on both sides of
(1) and divided by $\beta$ then taking $\beta\to\infty$ limit to reach

$S_{\rm best}(x)=\max_{1\leq i\leq K}\\{s(x-m_{i})+S_{\rm
best}(x-m_{i})\\}\;,$ (2)

where $S_{\rm best}(x)$ records the best path score among all paths reaching
position $x$. This update method, also termed dynamic programming, records the
lowest energy and lowest energy path reaching a given point $x$. The lowest
energy among all possible at position $x$ is simply $-S_{\rm best}(x)$ and the
associated path can be obtained by tracing backwards the incoming steps. It is
interesting to observe that one can also obtain the worst score at each
position via dynamical programming

$S_{\rm worst}(x)=\min_{1\leq i\leq K}\\{s(x-m_{i})+S_{\rm
worst}(x-m_{i})\\}\;.$ (3)

The full thermodynamic characterization demands more information than the
ground state energy. In principle, one may obtain the full partition function
using eq. (1) evaluated at various temperatures. This procedure, however,
hinders analytical property such as determination of the average energy

$\langle E\rangle\equiv-\frac{\partial\ln Z}{\partial\beta}\;.$

A better starting point may be achieved if one can obtain the density of
states $D(E)$. In this case, we have

$\displaystyle Z$ $\displaystyle\equiv$ $\displaystyle\int dEe^{-\beta E}D(E)$
$\displaystyle\langle E\rangle$ $\displaystyle=$ $\displaystyle\frac{\int
dEe^{-\beta E}ED(E)}{\int dEe^{-\beta E}D(E)}\;.$

Note that if the ground energy $E_{\rm grd}$ of the system is bounded from
below, the partition function is simply a Laplace transform of a modified
density of states given by

$Z=e^{-\beta E_{\rm grd}}\int_{0}^{\infty}dEe^{-\beta E}\tilde{D}(E)$

where $\tilde{D}(E)\equiv D(E-E_{\rm grd})$ and

$\langle E\rangle=E_{\rm grd}+\frac{\int_{0}^{\infty}dEe^{-\beta
E}E\tilde{D}(E)}{\int_{0}^{\infty}dEe^{-\beta E}\tilde{D}(E)}$

This implies that the density of states $D(E)$ together with the ground state
energy $E_{\rm grd}$ determine all the thermodynamic behavior of the system.
In the next section, we will explain how to obtain the density of states using
the dynamical programming technique as well as how to extend this approach to
more complicated situations that will be useful in characterizing the score
statistics in MS-based peptide sequencing.

##  3 Obtaining the Density of States

The density of states is related to the energy histogram in a simple way. The
number of states between energies $E$ and $E+\eta$ (with $\eta\ll 1$) is given
by $D(E)\eta$. If we happen to use $\eta$ as the energy bin size for energy
histogram, the count $C(E)$ in the bin with energy $E$ is simply $D(E)\eta$
and the density of states $D(E)=C(E)/\eta$. For simplicity, we will assume
that the all the on-site energies $-s(x)$ are integral multiple of $\eta$.
This implies that each path energy/score is also an integral multiple of
$\eta$. In the following subsections, we will use score density of states
instead of energy density of states.

###  3.1 The Simplest Case and its Application

We denote by $C(x,N)$ the number of paths reaching position $x$ with score
$N\eta$. With this notation, we can easily write down the recursion relation
for $C(x,N)$ as follows

$C(x,N)=\sum_{i=1}^{K}C(x-m_{i},N-\frac{s(x-m_{i})}{\eta})\;.$ (4)

This recursion relation allows us to compute the density of states in the same
manner as computing the partition function (1) except that we need to have an
additional dimension for score at each position $x$. As an even simpler
application of this recursion relation, suppose that one is only interested in
the number of paths reaching position $x$, one may sum over the energy part on
both side of (4) and arrives at

$C(x)=\sum_{i}C(x-m_{i}),$ (5)

which enables a very speedy way to compute the total number of paths reaching
position $x$. In the context of de novo peptide sequencing [6], this number
corresponds to the total number of a͡ll possible de novo peptides within a
given small mass range. Although simply obtained, this number may be useful
for providing rough statistical assessment in de novo peptide sequencing.

###  3.2 The More Realistic Case

In general, one may wish to associate with each hop an energy $h$ or one may
wish to introduce some kind of score normalization based on the number of
hopping steps. This is indeed the case when applying this framework to MS-
based peptide sequencing where a peptide length factor adding or multiplying
to the overall raw score is a common practice. In this case, it becomes
important to keep track the number of hops made in each path. We may further
categorize the counter $C(x,N)$ into $\sum_{L}C(x,N,L)$. That is, we may
separate the paths with different number of steps from one another and arrive
at a finer counter $C(x,N,L)$ which records the number of paths reaching
position $x$ with score $N\eta$ and with $L$ hopping steps.

It is rather easy to write down the recursion relation obeyed by this fine
counter

$C(x,N,L)=\sum_{i=1}^{K}C(x-m_{i},N-\frac{s(x-m_{i})}{\eta},L-1)\;.$ (6)

This recursion relation allows us to renormalize the raw score based on the
number of steps taken. For example, for RAId_DbS [7], a database search method
we developed, we divide the raw score obtained by $2(L-1)$ for any peptide
(path) of $L$ amino acids (hopping steps) to get better sensitivity in peptide
identification.

In principle, the recursion relations given by (4-6) are all one-dimensional
updates. The only difference is the internal structure of counters at each
position $x$. For (5), the counter is just an integer and has no further
structure. For (4), the counter at each position has a 1D structure indexed by
the score. For (6), the counter at each position $x$ has a 2D structure
indexed by both the score and the number of hopping steps. This means that in
terms of solving the problem using dynamical programming, it is always a 1D
dynamical programming with different degrees of internal structure that may
lengthen the execution time when shifting from the simplest case (5) to the
more complicated case (6). Obviously at each position $x$, there is an upper
bound and a lower bound for score and also for the number of hopping steps
accumulated. We shall call them $S_{\rm best}(x)$, $S_{\rm worst}(x)$, $L_{\rm
max}(x)$ and $L_{\rm min}(x)$ respectively. The first two quantities may be
obtained by eqs. (2) and (3) respectively. We provide the recursions for the
two latter quantities below

$\displaystyle L_{\rm max}(x)$ $\displaystyle=$ $\displaystyle\max_{1\leq
i\leq K}\\{L_{\rm max}(x-m_{i})\\}+1\;,$ (7) $\displaystyle L_{\rm min}(x)$
$\displaystyle=$ $\displaystyle\min_{1\leq i\leq K}\\{L_{\rm
min}(x-m_{i})\\}+1\;.$ (8)

Eqs. (2-3) and (7-8) provides the ranges for both the scores and the number of
cumulative hopping steps at each position $x$ via simple dynamic programming.
As we will discuss later, this information enables a memory-efficient
computations of score histograms.

##  4 Application in MS-based Peptide Sequencing

In this section, we focus on an important subject in modern biology – using MS
data to identify the numerous peptides/proteins involved in any given
biological process. Because of the peptide mass degeneracies and the limited
measurement accuracy for the peptide mass-to-charge ratio, using MS2 spectra
is more effective in peptide identifications. In a MS2 setup, a selected
peptide with its mass identified by the first spectrometer is fragmented by
noble gas, and the resulting fragments are analyzed by a second mass
spectrometer. Although such MS2-based proteomics approaches promise high
throughput analysis, the confidence level assignment for any peptide/protein
identified is challenging.

The majority of peptide identification methods are so-called database search
approaches. The main idea is to theoretically fragment each peptide in a
database to obtain the corresponding theoretical spectra. One then decides the
degree of similarity between each theoretical spectrum and the input query
spectrum using a scoring function. The candidate peptides from the database
are ranked/chosen according to their similarity scores to the query spectrum.
Although one may assign relative confidence levels among the candidate
peptides via various (empirical) means, an objective, standardized calibration
exists only recently [8]. In our earlier publications [5, 9], we proposed to
tackle this difficulty by using a de novo sequencing method to provide an
objective confidence measure that is both database-independent and takes into
account spectrum-specific noise. In this paper, we will provide concrete
algorithms for such purpose.

To begin, consider a spectrum $\sigma$ with parent ion mass range
$[w-\delta,w+\delta]$, we denote by $\Pi(w,\delta)$ the set of all “possible”
peptides with masses in this range. Given a peptide $\pi$ from
$\Pi(w,\delta)$, the associated quality score $S(\pi,\sigma)$ is defined by a
prescribed scoring system. The score distribution of $S(\pi,\sigma)$ within
$\Pi(w,\delta)$ provides naturally a likelihood measure for any given peptide
$\pi$ to the the correct one.

However, as described earlier [5], this seemingly straightforward idea faces
two difficulties in terms of implementations. First, unlike the DPRM problem
for which the function to be optimized is defined without ambiguity, the
choice of the scoring function is somewhat empirical because the parameters
used in the scoring must be trained using a training data set. Further,
because of different instruments and experimental setups, it seems impossible
to design a scoring system such that the correct peptide for each spectrum has
the highest score among all possible peptides; the application of a given
scoring function to general cases may require a leap of faith. Second, even
after the scoring function is chosen, it is not known how to find the peptide
$\pi_{o}$ that maximizes $S(\pi,\sigma)$ as well as the score distribution
${\rm pdf}(S)$ within $\Pi(w,\delta)$ other than by the generally impractical
procedure of examining all members of $\Pi(w,\delta)$.

The first difficulty can be alleviated by validating high scoring de novo
peptides via database searches [9] and is not the main focus of the current
paper. Note that a partial solution to the second problem via iterative
mapping when nonlocal score contributions exist is provided earlier [5]. Here
we tackle the second problem head on when the scoring function used does not
contain nonlocal contribution other than a final renormalization with respect
to the peptide length. Our algorithms contains two parts: computer memory
allocations and dynamical programming update. Prior to discussing these two
parts, however, we first address the important issue of choosing a good mass
unit.

###  4.1 Choosing a Good Mass Unit

The goal here is to choose a mass unit $\Delta$ and expresses the molecular
mass of each amino acid as an integral multiple of this unit. For example, one
may choose $\Delta$ to be $0.1$ Dalton (Da), and round the molecular mass of
each amino acid to be an integral multiple of $0.1$ Da. Once a mass unit is
chosen, all the masses under consideration are integral multiples of this
unit. It turns out that different choices of the mass unit leads to different
maximum cumulative mass error. As a specific example, consider using
$\Delta=0.1$ Da as the mass unit. The mass of Alanine, with true mass
$71.03711538$ Da, is now represented as $710\Delta$. This molecular mass
expression is $0.03711538$ Da smaller than the true molecular mass of Alanine.
When this happens, the integral mass representation has a mass smaller than
the true mass, and we call such type of mass error a down-error. Now the amino
acid Tryptophan with molecular mass $186.07931613$ Da will be assigned an
integral mass of $1861\Delta$, which has an extra of $(0.1-0.07931613)$ Da
compared to the true mass. We call this type of mass error the up-error.

The ratio of the mass error to the real molecular mass when multiplied by
$3000$ Da provides the cumulative maximum error that can be induced by a
single amino acid at $3,000$ Da mass. For a fixed mass unit, we went over this
mass error analysis for each of the twenty amino acids and documented the
largest up-error and down-error. The larger one between the maximum up-error
and the maximum down-error is called the max-error. To search for best mass
units that minimize the max-error at $3,000$ Da, we went over all possible
mass unit ranging from $0.005$ Da to $1.005$ Da in step of $10^{-6}$ Da.
Interestingly enough, we found a discrete list of mass units that have smaller
max-error compared to their nearby mass units. These numerically found magic
mass units are summarized in table 1.

Table 1: A list of best mass units in Da. The abbreviation “m.u.e.” stands for “maximum up-error,” while “m.d.e.” stands for “maximum down-error.” The maximum up-error, maximum down-error, and max-error are evaluated in extrapolation to $3,000$ Da as described in the text. The abbreviation “a.a.w.m.u.” stands for “amino acid with maximum up-error,” while “a.a.w.m.d.” stands for “amino acid with maximum down-error.”  mass unit | m.u.e. | a.a.w.m.u. | m.d.e. | a.a.w.m.d. | max-error
---|---|---|---|---|---
0.006070 | 0.041980 | Tryptophan | 0.037455 | Cysteine | 0.041980
0.007300 | 0.041495 | Methionine | 0.061276 | Asparagine | 0.061276
0.017540 | 0.094183 | Cysteine | 0.121977 | Proline | 0.121977
0.021500 | 0.199585 | Arginine | 0.182283 | Asparagine | 0.199585
0.054470 | 0.453793 | Asparagine | 0.347792 | Alanine | 0.453793
0.065400 | 0.553492 | Lysine | 0.536989 | Alanine | 0.553492
0.109450 | 0.908287 | Proline | 0.900898 | Lysine | 0.908287
0.110300 | 0.962781 | Histidine | 0.858742 | Lysine | 0.962781
0.110320 | 0.960176 | Aspartate | 0.907801 | Histidine | 0.960176
0.500208 | 0.980357 | Cysteine | 0.983149 | (Iso)Leucine | 0.983149
1.000416 | 0.980357 | Cysteine | 0.983149 | (Iso)Leucine | 0.983149

Once a mass unit is chosen, all the amino acid masses are effectively
integers. To obtain the score histogram of all de novo peptides when queried
by a spectrum $\sigma$ with parent molecular mass $w$ (with N- and C- terminal
groups of the peptide stripped away), we first construct a mass array where
index $k$ corresponds a molecular mass $k\Delta$. To encode all possible
peptides with molecular mass up to $w$, we need to have an array of size
$w/\Delta+1$. Apparently, when a larger mass unit is used, the size of the
mass array is smaller and thus reduces computation time. However, as one may
see from table 1, the larger mass unit is also accompanied by a larger max-
error and might not be preferred when high mass accuracy is the first
priority.

###  4.2 Efficient Memory Allocation

The basic idea of our algorithm is to encode all possible peptides in the mass
array by linking pointers, analogous to the consecutive hopping steps in the
1D hopping model. For an amino acid $a$, let $n(a)$ represents its
corresponding integer mass in unit of $\Delta$. For a peptide made of
$[a_{1},a_{2},\ldots,a_{M}]$, it will have a hopping trajectory in the
molecular array given by $[0,x_{1},x_{2},\ldots,x_{M}]$ with $x_{i\geq
1}\equiv\sum_{j=1}^{i}n(a_{j})$. Let us also denote $x_{M}$ by $x_{F}$ to
indicate that it is the terminating point of the path. Apparently, all
possible peptides with molecular masses equal to $x_{F}\Delta$ will all have
corresponding hopping paths starting at the origin and terminating at $x_{F}$.
Through appropriate pointer linking, one may therefore encode all possible
peptides with molecular mass $x_{F}\Delta$ in a one-dimensional mass array.

For a given spectrum $\sigma$, depending on the score function used, one may
calculate local score contributions at each mass index. This step is done once
only for the whole mass array, and need not be repeated for each candidate
peptide. In a typical MS2 experimental spectrum, there always exists some
level of parent ion mass uncertainty. Once the size of the mass uncertainty is
specified, we only need to examine de novo peptides whose corresponding
hopping paths terminating at a few consecutive mass indices. This indicates
that some of the mass indices of the aforementioned mass array may not even be
used in this context. Below we describe how to efficiently obtain relevant
mass indices and only allocate computer memories for those masses.

Assume that the possible terminating points are $F_{1},F_{2},\ldots,F_{k}$
with $F_{j+1}=F_{j}+1$. The update rules described in Eqs. (2-3), (5), and
(7-8) will also be used at this stage. The following pseudocode describes our
algorithm.

Initialize the mass_index = 0 entry
$S_{\rm best}=S_{\rm worst}=L_{\rm max}=L_{\rm min}=0$; $C$=1;
REMARK: Max_aa is the maximum number of amino acids considered
for (aa_index = 0; aa_index $<$ Max_aa; aa_index ++) {
label occupancy of n(aa_index);
at n(aa_index) attach a pointer back to 0;
update $S_{\rm best}$, $S_{\rm worst}$, $L_{\rm max}$, $L_{\rm min}$, $C$ at
n(aa_index);
}
for (mass_index = 1; mass_index $<=F_{k}$; mass_index ++){
if (mass_index occupied ?) {
for (aa_index = 0; aa_index $<$ Max_aa; aa_index ++) {
label occupancy of (mass_index + n(aa_index));
at mass_index+n(aa_index) attach a pointer to mass_index;
update $S_{\rm best}$, $S_{\rm worst}$, $L_{\rm max}$, $L_{\rm min}$, $C$ at
(mass_index + n(aa_index));
}
}
}
for (mass_index = $F_{k}$ ; mass_index $>=F_{1}$; mass_index --){
backtrack all possible paths $\to$ final occupied entries;
}

The last step in the algorithm above identifies relevant mass indices,
mass_indices that will be traversed by the hopping paths of all peptides with
molecular masses in the range $[F_{1}\Delta,F_{k}\Delta]$. We only need to
allocate computer memory associated with those sites. For each of these
relevant sites, we also know the values of $S_{\rm best}$, $S_{\rm worst}$,
$L_{\rm max}$, $L_{\rm min}$, and the total number of peptides reaching that
site through the algorithm above. One may therefore allocate a 2D array of
size $(S_{\rm best}(i)-S_{\rm worst}(i))/\eta\times(L_{\rm max}(i)-L_{\rm
min}(i))$ for each relevant mass_index $i$ for later use.

###  4.3 Main Algorithm and some Results

Once memory allocation for relevant mass_indices is done, we can efficiently
go through those relevant sites to obtain the 2D score histogram that we
mentioned. In the pseudocode below, update is performed using eq. (6). We now
demonstrate the very simple main algorithm

Initialize all the fine counters $C(x,N,L)=0$
except $C(x=0,N=0,L=0)$=1;
for (aa_index = 0; aa_index $<$ Max_aa; aa_index ++) {
update $C(x,N,L)$ at $x=$n(aa_index);
}
for (mass_index in ascendingly ordered relevant mass_indices){
for (aa_index = 0; aa_index $<$ Max_aa; aa_index ++) {
update $C(x,N,L)$ at $x=$(mass_index + n(aa_index));
}
}

We now define the final 2D counter

$Y(N,L)\equiv\sum_{i=1}^{k}C(F_{i},N,L)\;.$ (9)

Apparently, in the 1D hopping model when allowing $k$ consecutive terminating
points, the resulting density of states $D(E)$ can now be expressed as
$D(E=-N\eta)=\sum_{L}Y(N,L)/\eta$. If one were interested in normalizing the
final score in a path-length dependent manner, one will has the following
generic transformation

$H(E)=\sum_{L}\int
dE^{\prime}\frac{Y(E^{\prime}=-N\eta,L)}{\eta}\delta\left(E-f(E^{\prime},L)\right)$
(10)

where $f(E^{\prime},L)$ is a generic length-normalized energetic function that
takes the raw energy $E^{\prime}$ with $L$ hopping steps and turn them into a
new energy $f(E^{\prime},L)$, and $\int dE^{\prime}\to\eta\sum_{N}$ is
understood.

Using a real experimental MS2 spectrum of parent ion mass $2254.7\pm 3.0$ Da
and a raw scoring function (RAId_DbS [7] raw score without divided by $2(L-1)$
with $L$ being the peptide length), we obtained a 2D score histogram. From
this 2D score histogram, we can compute the average peptide length $\langle
L\rangle$ as well. We then transform the 2D score histogram using two
different $f$ functions. In the first case,
$f(E^{\prime},L)=E^{\prime}/2(\langle L\rangle-1)$, meaning that one just
divides the score by a constant given by $2(\langle L\rangle-1)$. In the
second case, we use the RAId_DbS scoring function where
$f(E^{\prime},L)=E^{\prime}/2(L-1)$. In Figure 1, we show the two resulting
score histograms along with the fits to theoretical distribution function [7].
As one may see from the figure, both histograms are well fitted by the
theoretical distribution function over at least 15 order of magnitudes. There
is difference, however, in the histograms obtained. In the first case, where
the score is merely divided by the average length, we have a wider score
distribution than that of the second case. This implies that a high scoring
hit out of the first type of scoring function will have a larger $P$-value
than that of the second type. This is perfectly reasonable because when using
the first type of raw scoring, very long peptides which by random chances are
more likely to hit on fragment peaks in the mass spectrum are less penalized
than the shorter peptides. As a consequence, one anticipates more false long
peptides out of the first type of scoring method than that of the second
scoring method. Therefore, one should assign a larger $P$-value to the former
case and a smaller $P$-value to the latter case. It is apparently important to
be able to obtain score histograms of the second scoring method. However, this
can only be achieved if one keeps the length information in the dynamical
programming update, see eq. (6).

Figure 1: Score histograms for raw score and RAId_DbS score. Note that the two
histogram cross each other at large score regime, indicating that the raw
score function might not be as effective as the RAId_DbS score, see text for
details.

##  5 Discussion, Summary, and Outlook

Our method may also be extended to other applications. In the case of
passenger routings, the $x$-axis actually represents time. The local score may
be viewed as the additional cost that may vary for different stops. Once the
problem is laid out, the 2D histogram obtained from our solution indicates the
number of equivalent routes in terms of additional costs and the total number
of stops. This problem should be interesting in its own right.

In this paper, we developed a new approach to obtain the density of states of
a 1D hopping problem in random potential. We have extended the simplest case
scenario and have shown that we can apply this method to provide a complete
score histogram for MS-based peptide sequencing problem. This important
information may be used for a more objective statistical significance
assignment in peptide identification. Our algorithm may also serve as a speedy
de novo algorithm. If one is only interested in getting the best scoring
peptide with length normalized score, one only needs to keep track of $S_{\rm
best}(x,L)$. Furthermore, it is straightforward to include in our de novo
algorithm post-translationally modified amino acids. The effect is simply an
enlargement of the alphabet. That is, instead of having 20 amino acids, we
will simply have more allowed masses but without needing to change any part of
the algorithm.

In the near future, we would like to build a web application that allows the
users to obtain information of interest. For example, a user might be
interested in knowing: given a parent ion molecular mass and a mass error
tolerance, how many de novo peptides can there be? Furthermore, we plan to
provide users with the full score histogram when a query spectrum is provided
and a scoring method is chosen. Our approach, founded on statistical physics,
can easily address this type of questions to provide useful information for
biological researches.

## Acknowledgement

This work was supported by the Intramural Research Program of the National
Library of Medicine at the National Institutes of Health.

## References

  * [1] M.R. Garey and D.S. Johnson, Computers and intractability, W.H. Freeman and company, New York (1979).
  * [2] D.A. Huse and C.L. Henley, Phys. Rev. Lett. 54, 2708 (1985).
  * [3] M. Kardar, Nucl. Phys. B290, 582 (1987).
  * [4] D.S. Fisher and D.A. Huse, Phys. Rev. B. 43, 10728 (1991).
  * [5] T.P. Doerr, G. Alves and Y.-K. Yu, Physica A 354, 558-570 (2005).
  * [6] J.A. Taylor and R.S. Johnson, Raid Commu. Mass Spect. 11, 1067 (1997) .
  * [7] G. Alves, A.Y. Ogurtsov and Y.-K. Yu, Biology Direct 2, art. no. 25 (2007).
  * [8] G. Alves, A.Y. Ogurtsov, W.W. Wu, G. Wang, R.-F. Shen and Y.-K. Yu, Biology Direct 2, art. no. 26 (2007).
  * [9] G. Alves and Y.-K. Yu, Bioinformatics 21, 3726-3732 (2005).