# Greedy versus Map-based Optimized Adaptive Algorithms
for random-telegraph-noise mitigation by spectator qubits

Behnam Tonekaboni<EMAIL_ADDRESS>Centre for Quantum Dynamics, Griffith
University, Yuggera Country, Brisbane, Queensland 4111, Australia  Areeya
Chantasri<EMAIL_ADDRESS>Optical and Quantum Physics Laboratory,
Department of Physics, Faculty of Science, Mahidol University, Bangkok, 10400,
Thailand  Centre for Quantum Computation and Communication Technology
(Australian Research Council),
Centre for Quantum Dynamics, Griffith University, Yuggera Country, Brisbane,
Queensland 4111, Australia  Hongting Song<EMAIL_ADDRESS>Centre for Quantum
Computation and Communication Technology (Australian Research Council),
Centre for Quantum Dynamics, Griffith University, Yuggera Country, Brisbane,
Queensland 4111, Australia  Qian Xuesen Laboratory of Space Technology, China
Academy of Space Technology, Beijing 100094, China  Yanan Liu
<EMAIL_ADDRESS>Centre for Quantum Computation and Communication
Technology (Australian Research Council),
Centre for Quantum Dynamics, Griffith University, Yuggera Country, Brisbane,
Queensland 4111, Australia  Howard M. Wiseman<EMAIL_ADDRESS>Centre for Quantum Computation and Communication Technology (Australian
Research Council),
Centre for Quantum Dynamics, Griffith University, Yuggera Country, Brisbane,
Queensland 4111, Australia

###### Abstract

In a scenario where data-storage qubits are kept in isolation as far as
possible, with minimal measurements and controls, noise mitigation can still
be done using additional noise probes, with corrections applied only when
needed. Motivated by the case of solid-state qubits, we consider dephasing
noise arising from a two-state fluctuator, described by random telegraph
process, and a noise probe which is also a qubit, a so-called spectator qubit
(SQ). We construct the theoretical model assuming projective measurements on
the SQ, and derive the performance of different measurement and control
strategies in the regime where the noise mitigation works well. We start with
the Greedy algorithm; that is, the strategy that always maximizes the data
qubit coherence in the immediate future. We show numerically that this
algorithm works very well, and find that its adaptive strategy can be well
approximated by a simpler algorithm with just a few parameters. Based on this,
and an analytical construction using Bayesian maps, we design a one-parameter
($\Theta$) family of algorithms. In the asymptotic regime of high noise-
sensitivity of the SQ, we show analytically that this $\Theta$-family of
algorithms reduces the data qubit decoherence rate by a divisor scaling as the
square of this sensitivity. Setting $\Theta$ equal to its optimal value,
$\Theta^{\star}$, yields the Map-based Optimized Adaptive Algorithm for
Asymptotic Regime (MOAAAR). We show, analytically and numerically, that MOAAAR
outperforms the Greedy algorithm, especially in the regime of high noise
sensitivity of SQ.

##  I Introduction

Qubit decoherence from environmental noises is one of the major obstacles to
building useful large-scale quantum computers [1, 2, 3, 4]. In the past years,
a wide range of noise-mitigation techniques have been proposed to prolong the
coherence time of computational qubits. Some of the techniques, such as the
dynamical decoupling (DD), aim at removing average effects of noise
continuously in time, via applying carefully designed controls onto qubits [5,
6, 7, 8, 9, 10, 11, 12]. Other techniques, such as those that fit in the
category of quantum error correction (QEC), are designed to monitor effects of
noise and correct errors via information-redundancy encoding [13, 14, 15].

However, both DD and QEC approaches require direct controls or measurements on
the qubits of interest, which could in turn create more sources of noise.
Therefore, in the particular case where it is possible to keep the computing
(_data_) qubits well-isolated from their environment, new techniques, that can
both correct noise effects while minimizing direct contact to the data qubits,
are required. Such techniques would demand an additional entity, _e.g._ , a
sensor or a probe, in proximity of the data qubits, that can sense the
problematic noise and any correction to the data qubits can be applied only
when needed.

In this paper, we are interested in the recently proposed idea of spectator
qubits (SQs) [16, 17, 18], which are qubit-type probes that are assumed to be
much more sensitive to target noise than data qubits and can be easily
measured. This could be applicable in solid-state quantum computing, where
different types of qubits with very different properties can co-exist. For
example, the data qubits could be the spin qubits with low sensitivity to
noise [19, 20, 21] and SQs could be quantum dots nearby with an ability of
high-efficiency readout [22, 23, 24]. Assuming a well-isolated spin-donor
qubit in silicon substrate, there will still be lingering effect from low-
frequency noise, such as charge noise, which can affect qubits at different
locations simultaneously [25, 26] and is considered difficult to remove [27].
Therefore, we are interested in constructing a theoretical model using the
data-spectator setup to investigate a plausible optimal regime for the noise
mitigation. It turns out that the solution is not trivial, even with a simple
regime with single data and spectator qubits.

We have presented our most important findings on this problem in a Companion
Letter (CL) [28]. Specifically, we showed that by performing projective
measurements on the SQ at certain times and in certain bases, that are chosen
adaptively in an optimized manner, the decoherence of the data qubit can be
suppressed by an amount which is quadratic in the SQ sensitivity. This result
holds under ideal conditions, and the optimal suppression is quite non-
trivial, as it is better than that can be achieved by simply maximizing the
data qubit coherence in each, arbitrarily small, time step (Greedy algorithm).
In this paper, we give the detailed derivations of the results in the CL [28],
and explore in much more depth the comparison between the Greedy algorithm and
the optimized algorithm. But to explain more precisely the contents of this
paper, it is necessary first to give some more theoretical background.

Figure 1: Schema for the adaptive measurement strategy on SQ. The orange line
shows the dynamics of RTP which jumps between two values indicated with pale
dashed orange lines. The blue line represents the SQ which is reset (triangle)
after each measurement (purple semicircle). Storing measurement results are
represented by purple arrows toward the measurement record line $Y_{n}$
(purple). These measurement records are used via the adaptive strategy (black
arrows) to decide on the next measurement time $t_{n}$ and next measurement
angle $\theta_{n}$. Lastly, at final time $T$, the data qubit (pink) is
corrected using information $Y_{N}$, (Teal arrow downward).

We model the decoherence of the data qubit as being caused by low-frequency
phase noise arising from the random motion of a trapped charge caused by
impurities, as is common in the experimental solid-state qubit setting [29,
30, 31, 27, 32]. Specifically, we model the noise $z(t)$ as a two-state
fluctuator [33, 34, 35, 36], with values $\pm 1$ as shown in Fig. 1. This can
be mathematically described by the random telegraph process (RTP), with two
flip rates: $\gamma_{\uparrow}$ and $\gamma_{\downarrow}$ [37, 38]. To monitor
the noise, a SQ is placed nearby, experiencing the same phase noise, but with
the higher sensitivity. This is described by a total Hamiltonian of both
qubits,

$\hat{H}_{\rm tot}=\frac{\kappa}{2}\hat{\sigma}_{z}^{\rm
d}\,z(t)+\frac{K}{2}\hat{\sigma}_{z}^{\rm s}\,z(t).$ (1)

Here $\kappa$ and $K$ are the data qubits’ and SQ’s noise sensitivities and
$\hat{\sigma}_{z}^{\rm d,s}$ are their respective $z$-Pauli matrices. Also, we
are working in a rotating frame for both qubits (their frequencies need not be
the same or even similar) and with units where $\hbar=1$.

We assume that, in order to obtain information about the noise $z(t)$, the SQ
can be measured after arbitrary waiting times $\\{\tau_{1},\tau_{2},...\\}$,
with arbitrary angles $\\{\theta_{1},\theta_{2},...\\}$, and re-prepared in an
arbitrary state $|{\theta}\rangle^{\rm s}$. Here $\theta$ always indicates an
angle on the equator of the Bloch sphere, defining either a Pauli observable,
$\hat{\sigma}_{\theta}^{\rm s}\equiv\cos\theta\,\hat{\sigma}_{x}^{\rm
s}+\sin\theta\,\hat{\sigma}_{y}^{\rm s}$, or a state $|{\theta}\rangle^{\rm
s}$ which is the $+1$ eigenstate of $\hat{\sigma}_{\theta}^{\rm s}$. The
resulting string of measurement readouts can then be used in estimating the
data qubit phase correction at some final time $T$, which is when the data
qubit is required. The goal is to maximize the resultant coherence of the data
qubit, averaging over all possible realizations of the noise.

Given this model and goal, the question is how one should choose the
$\tau_{n}$’s, the $\theta_{n}$’s, and the final phase correction. In this
work, we apply the separation principle from control theory, and the structure
of the reward function (which defines the goal), to greatly simplify the
problem. We show that the information necessary for all future controls can be
stored in a complex 2-vector and propagated from one measurement to the next
by a _Bayesian map_ , a complex $2\times 2$ matrix. We then use the maps to
numerically find the locally optimal—also known as _Greedy_ —algorithm.

Locally optimal algorithms, in general, are not always globally optimal, and
we find that to be the case for the problem here. We are interested in the
regime where the SQ approach is most useful, to wit,

$T^{-1},~{}\kappa\ll\gamma_{\uparrow,\downarrow}\ll K.$ (2)

We analyze how the Greedy algorithm works in this regime, and this suggests a
simple, single parameter ($\Theta$) family of models where $\tau_{n}=\Theta/K$
and $\theta_{n}=\pm\Theta$. Here the sign of $\theta_{n}$ is determined
_adaptively_ , and reflects the algorithm’s best estimate of $z(t_{n})$, the
current state of RTP. In the asymptotic regime (2), we can solve analytically
for the performance of this family of algorithms. Choosing the optimal
parameter value $\Theta^{\star}$ yields a measurement and control strategy
that is plausibly optimal, and which we call MOAAAR (Map-based Optimized
Adaptive Algorithm in the Asymptotic Regime). This can reduce the decoherence
rate of the data qubit, relative to the no control rate, by the multiplicative
factor

$\displaystyle H^{\star}\,(\bar{\gamma}/K)^{2},\text{ where }H^{\star}\approx
1.254,$ (3)

where $\bar{\gamma}:=(\gamma_{\uparrow}+\gamma_{\downarrow})/2$. That is, as
mentioned above, the data qubit lifetime enhancement offered by the MOAAAR
algorithm scales like the square of the sensitivity of the SQ. While the same
quadratic scaling is observed for the Greedy algorithm, we find that the
prefactor $H^{\star}$ for the decoherence rate with MOAAAR is smaller than
that from Greedy.

The structure of this paper is as follows. First, in Sec. II, we briefly
discuss the physics and mathematical description of the RTP. We then, in Sec.
III, calculate the dephasing of the data qubit in the absence of any control,
and formulate it as a matrix product, which will be useful later. We end this
section by looking at how, in general, additional information can be used to
improve the coherence in the data qubit. In Sec. IV, we construct the setup
and measurement on the SQ and discuss the notion of adaptive measurement
strategy in the context of the SQ. Then, in Sec. V, we use Bayesian statistics
and formulate the expected coherence using the gathered information via the
SQ. We write the Bayesian map as a matrix product which updates a _coherence
vector_ after each measurement on the SQ. We also introduce a phase space
formalism which includes all the sufficient information. We then optimize our
measurement strategy using the locally optimization algorithm (Greedy) in Sec.
VI and the global optimization approach (MOAAAR) in Sec. VII. In Sec. VIII, we
compare these two optimization strategies by constructing a closed-form
formula for each that applies even outside the asymptotic regime. We conclude
our work with a discussion of future research directions in Sec. IX.

##  II Charge noise and random telegraph process (RTP)

Charge noise is one of the causes of decoherence in a wide range of physical
systems including superconducting qubits [39], quantum dots [40], Nitrogen
vacancy centers [41], trapped ions [42], and semiconductors [33, 31, 43, 44,
27, 35]. The noise is often described by stochastic motion of electrons or
holes trapped at impurities in the devices, which can be mathematically
modeled as a random telegraph process (RTP), whose value jumps between two
values [37, 38]. Consider the RTP noise parameter $z_{t}:=z(t)$, which can be
either $+1$ or $-1$. It makes a jump from $z=-1$ to $z=+1$ with rate
$\gamma_{\uparrow}$, and $z=+1\rightarrow-1$ with rate $\gamma_{\downarrow}$.
Define the probability vector

$\underline{P}_{t}:=\left(\begin{matrix}\wp(z_{t}=+1)\\\
\wp(z_{t}=-1)\end{matrix}\right),$ (4)

where $\wp(z_{t}=\pm 1)$ is the probability of $z_{t}=\pm 1$ at time $t$.
Then, with the jump rates defined above, one may write a (classical) master
equation for vector $\underline{{P}}_{t}$ as

$\underline{\dot{P}}_{t}=\mathbf{J}\underline{P}_{t},$ (5)

where $\mathbf{J}$ is the jump rate matrix,

$\mathbf{J}=\left(\begin{matrix}-\gamma_{\downarrow}&+\gamma_{\uparrow}\\\
+\gamma_{\downarrow}&-\gamma_{\uparrow}\end{matrix}\right).$ (6)

If the RTP evolves for a duration of $\tau$, between times $t^{\prime}$ and
$t^{\prime\prime}=t^{\prime}+\tau$, the master equation (5) has a solution of
the form
$\underline{P}_{t^{\prime\prime}}=e^{\mathbf{J}\tau}\underline{P}_{t^{\prime}}$.
One may rearrange the matrix exponential and write the solution as

$\underline{P}_{t^{\prime\prime}}=\left(\mathbf{I}+\frac{1-e^{-2\bar{\gamma}\tau}}{2\bar{\gamma}}\mathbf{J}\right)\underline{P}_{t^{\prime}}.$
(7)

where $\bar{\gamma}=(\gamma_{\uparrow}+\gamma_{\downarrow})/2$ as before, and
$\mathbf{I}$ is the identity matrix. It is worth mentioning that when
$\tau\rightarrow\infty$, the probability vector
$\underline{P}_{t^{\prime\prime}}$ converges to its steady state given by

$\underline{{P}}_{\rm
ss}=\frac{1}{2\bar{\gamma}}\left(\begin{matrix}\gamma_{\uparrow}\\\
\gamma_{\downarrow}\end{matrix}\right).$ (8)

This will be used throughout the paper as the state of RTP when no information
about the noise is available, such as at the initial time.

We model the effect of the charge noise on the data and spectator qubits by
the Hamiltonian in Eq. (1). From this, it is apparent that, under free
evolution in some time interval $(t^{\prime},t^{\prime\prime})$, the qubits’
accumulated phase (angle of rotation around the $\hat{\sigma}_{z}$ axis) will
be proportional to the integral of the noise in that interval,

$x:=\int_{t^{\prime}}^{t^{\prime\prime}}\\!\\!z(s){\rm d}s.$ (9)

To determine the dephasing caused by this random accumulated phase (more on
this in the next section) it is useful to calculate the average,
$\int\wp(x)e^{ikx}{\rm d}x$, of the exponential factor, $e^{ikx}$, for
arbitrary $k$. Actually, it is even more useful to calculate the integral
$\int e^{ikx}\wp(x,z_{t^{\prime\prime}}|z_{t^{\prime}}){\rm d}x$ for arbitrary
$k$, $z_{t^{\prime\prime}}$ and $z_{t^{\prime}}$. This is, of course, just the
Fourier transform of $\wp(x,z_{t^{\prime\prime}}|z_{t^{\prime}})$, and can be
calculated using Eqs. (7) and (9) as

$\displaystyle\int$ $\displaystyle
e^{ikx}\wp(x,z_{t^{\prime\prime}}|z_{t^{\prime}}){\rm
d}x=\exp(-\bar{\gamma}\tau)$
$\displaystyle\times\begin{cases}\cosh\left(\dfrac{\lambda}{2}\tau\right)-s\dfrac{\eta}{\lambda}\sinh\left(\dfrac{\lambda}{2}\tau\right),&\text{for
}z_{t^{\prime}}=z_{t^{\prime\prime}}=s\\\
\dfrac{2\gamma_{\uparrow}}{\lambda}\sinh\left(\dfrac{\lambda}{2}\tau\right),&\text{for
}z_{t^{\prime\prime}}=-z_{t^{\prime}}=+1\\\
\dfrac{2\gamma_{\downarrow}}{\lambda}\sinh\left(\dfrac{\lambda}{2}\tau\right),&\text{for
}z_{t^{\prime\prime}}=-z_{t^{\prime}}=-1\end{cases}$ (10)

where $s\in\\{-1,+1\\}$ and

$\displaystyle\lambda(k)$ $\displaystyle=$
$\displaystyle\sqrt{(\gamma_{\downarrow}+\gamma_{\uparrow})^{2}-4ik(\gamma_{\downarrow}-\gamma_{\uparrow})-4k^{2}},$
(11a) $\displaystyle\eta(k)$ $\displaystyle=$
$\displaystyle(\gamma_{\downarrow}-\gamma_{\uparrow})-2ik,$ (11b)

are functions of the Fourier variable $k$ (See appendix A and Refs. [35, 45]
for the details of its derivation). In the following sections, we will use
these results to derive the statistics of data and spectator phases,
conditioned on measurement results and/or qubit controls.

##  III Data qubit dephasing

In this section, we define the data qubit’s coherence and analyze its
dephasing due to the RTP noise, both for the cases with no control and with
the control phase correction. For the no-control case, we show that the
qubit’s coherence can be written as a matrix equation involving a complex two-
by-two matrix, ${\bf H}$. The matrix representation for the no-control
coherence will be extended to the version including measurement strategies and
readouts of the spectator qubit in Section V.

###  III.1 Formulation of data qubit coherence

From the total Hamiltonian in Eq. (1), if we include the control (phase
correction) applied at the final time, $t=T$, we can write the Hamiltonian of
the data qubit as

$\hat{H}_{\rm d}(t)=\frac{\kappa}{2}\,\hat{\sigma}_{z}^{\rm
d}\,z(t)+\hat{H}_{\rm ctrl}\delta(t-T).$ (12)

The first term of the Hamiltonian describes a stochastic phase rotation around
the $z$-axis of the data qubit’s Bloch sphere. To evaluate the data qubit’s
decoherence due to this term, we consider a qubit’s density matrix in the
$z$-basis affected by a random noise,

$\displaystyle\rho^{\rm d}(\Xi)=\,$
$\displaystyle\rho_{++}|{+1}\rangle_{z}^{\rm
d}\langle{+1}|+\rho_{+-}e^{-i\phi(\Xi)}|{+1}\rangle_{z}^{\rm d}\langle{-1}|$
$\displaystyle+\rho_{-+}e^{+i\phi(\Xi)}|{-1}\rangle_{z}^{\rm
d}\langle{+1}|+\rho_{--}|{-1}\rangle_{z}^{\rm d}\langle{-1}|,$ (13)

where $|{+1}\rangle_{z}^{\rm d}$ and $|{-1}\rangle_{z}^{\rm d}$ are the two
eigenstates of $\hat{\sigma}_{z}^{\rm d}$, _i.e._ , $\hat{\sigma}_{z}^{\rm
d}|{\pm 1}\rangle_{z}^{\rm d}=\pm|{\pm 1}\rangle_{z}^{\rm d}$. The terms
$\rho_{++},\rho_{--},\rho_{+-},\rho_{-+}$ are density matrix elements of the
data qubit’s state at the initial time ($t=0$), before affected by the noise.
The Hamiltonian Eq. (12) thus introduces an additional phase, denoted by
$\phi=\phi(\Xi)$, which is a function of a set of variables, $\Xi$, which can
include both the noise that affects the data qubit’s phase, and the control
applied to it.

Since the value of $\Xi$ varies from run-to-run, we are interested in the
average final state of the qubit, i.e., $\langle\rho^{\rm
d}(\Xi)\rangle_{\Xi}$, where the subscript $\Xi$ indicates that this is the
random variable with respect to which the average is taken. The data qubit’s
final purity is completely determined by the magnitude of the off-diagonal
elements of Eq. (III.1), which gives $|\rho_{+-}\left\langle
e^{i\phi(\Xi)}\right\rangle_{\Xi}|$. Since the element $\rho_{+-}$ is not
affected by the noise, we can thus define the quantity to be maximized, which
we call the _coherence_ , as

$\mathcal{C}:=\left|\left\langle e^{i\phi(\Xi)}\right\rangle_{\Xi}\right|.$
(14)

The expected average of random phases typically leads to the value of
coherence being less than one (its maximum value) and can be as low as zero
for the case of uniformly distributed phases in $[0,2\pi)$. The reduction of
coherence from random phases is called qubit dephasing. It is also convenient
later to define a _complex coherence_ ,

$\mathcal{A}:=\left\langle e^{i\phi(\Xi)}\right\rangle_{\Xi},$ (15)

without the absolute value, where $\mathcal{C}=|\mathcal{A}|$. The information
of the argument of $\mathcal{A}$ will be useful later, in estimating the phase
information of the data qubit.

Let us define capital-letter variables for a total accumulated noise, from
time $t=0$ to any time $t$,

$X:=\int_{0}^{t}\\!\\!z(s){\rm d}s,$ (16)

and for a total information,

$Y:=\text{all acquired information about }X,$ (17)

which is from any measurements up to the time $t$. In the following
subsections, we simplify the representation of a data qubit’s state in Eq.
(III.1) to showing only its phase. Let us define a state of the data’s qubit
phase as

$|{\phi}\rangle^{\rm d}:=\frac{1}{\sqrt{2}}\left(|{+1}\rangle_{z}^{\rm
d}+e^{i\phi}|{-1}\rangle_{z}^{\rm d}\right).$ (18)

Then, the initial phase of the data qubit is $|{\phi=0}\rangle^{\rm d}$ (a
zero noisy phase) and we consider two cases: the no-control case and the case
with the noise correction. For the former, $\hat{H}_{\rm ctrl}=0$, the total
data noisy phase will be a function of noise, _i.e._ , $\phi(\Xi)=\phi(X)$.
For the latter, we introduce a phase-correction control, $\hat{H}_{\rm
ctrl}=-ic(Y)\hat{\sigma}_{z}^{\rm d}/2$, which describes an instantaneous
rotation of the data qubit’s state around the $z$-axis with an angle $-c(Y)$
that depends on the information, $Y$. In this case, the data qubit’s phase
will be a function of both the noise and the control, _i.e._ ,
$\phi(\Xi)=\phi(X,Y)$. We will show how to compute the data qubit’s coherence
for both cases in the following subsections.

###  III.2 Decay of no-control (nc) coherence

For the ‘no-control’ (nc) case, $\hat{H}_{\rm ctrl}=0$, the data qubit only
evolves stochastically due to the RTP noise. Given the Hamiltonian, Eq. (12),
the data qubit’s state (represented by its noisy phase) at any time $t$,
becomes

$|{\phi(\Xi)}\rangle^{\rm
d}=\exp\left(-i\tfrac{\kappa}{2}\hat{\sigma}_{z}^{\rm d}\\!\int_{0}^{t}\\!{\rm
d}s\,z(s)\right)|{\phi=0}\rangle^{\rm d}=|{\kappa X}\rangle^{\rm d},$ (19)

which means that the data noisy phase, $\phi(\Xi)=\phi(X)=\kappa X$, is simply
proportional to the total accumulated noise. Following Eq. (14), replacing
$\Xi$ with $X$, we can calculate the no-control coherence as

$\displaystyle\mathcal{C}^{\rm nc}$ $\displaystyle:=\left|\left\langle
e^{i\phi(X)}\right\rangle_{X}\right|=\left|\int e^{i\kappa X}\wp(X){\rm
d}X\right|,$ (20)

and its no-control complex coherence

$\displaystyle\mathcal{A}^{\rm nc}$ $\displaystyle:=\int e^{i\kappa
X}\wp(X){\rm d}X,$ (21)

where $\wp(X)$ is the probability density function of the accumulated noise.

Since the variable $X$ is an accumulated RTP noise, it is convenient to write
the probability function $\wp(X)$ as marginalized over the initial
$z_{0}=z(t=0)$ and the end-point $z_{t}=z(t)$ of RTP, _i.e._ ,

$\wp(X)=\sum_{z_{t}}\sum_{z_{0}}\wp(X,z_{t}|z_{0})\wp(z_{0}).$ (22)

Substituting the above equation into Eq. (20), we find that the no-control
coherence is

$\displaystyle\mathcal{C}^{\rm nc}$
$\displaystyle=\left|\sum_{z_{t}}\sum_{z_{0}}\left[\int e^{i\kappa
X}\wp(X,z_{t}|z_{0}){\rm d}X\right]\wp(z_{0})\right|,$ (23)

where the term inside the square brackets is the Fourier form as we defined in
Eq. (II). Now, we see that Eq. (23) can be written as a map-based equation by
defining a vector

$\underline{A}_{0}:=\underline{P}_{t=0}=\left(\begin{matrix}\wp(z_{0}=+1)\\\
\wp(z_{0}=-1)\end{matrix}\right),$ (24)

using the definition in Eq. (4) and defining matrix elements,

$H_{z_{0}}^{z_{t}}(t,\kappa):=\int e^{i\kappa X}\wp(X,z_{t}|z_{0}){\rm d}X.$
(25)

Since $z_{0},z_{t}\in\\{-1,+1\\}$, this defines a two-by-two matrix,

${\bf
H}(t,\kappa):=\left(\begin{matrix}H^{+1}_{+1}(t,\kappa)&H^{+1}_{-1}(t,\kappa)\\\\[5.0pt]
H^{-1}_{+1}(t,\kappa)&H^{-1}_{-1}(t,\kappa)\end{matrix}\right).$ (26)

Thus, the coherence in Eq. (23) is then written as

$\mathcal{C}^{\rm nc}=\left|\underline{I}^{\top}{\bf
H}(t,\kappa)\underline{A}_{0}\right|,$ (27)

where $\underline{{I}}^{\top}=(1,1)$. We can calculate the coherence Eq. (27)
for the no-control case analytically (see details in the Appendix A). As we
will see later, reformulating the coherence as matrix equations simplifies
analytical and numerical calculations for optimal controls, particularly when
measurements on the spectator qubit are involved.

To simplify the discussion later, let us also define a complex 2-dimensional
coherence vector for the no-control case,

$\underline{A}^{\rm nc}={\bf
H}(t,\kappa)\underline{A}_{0}=\left(\begin{matrix}\wp(z_{t}=+1)\mathcal{A}_{|z_{t}=+1}\\\
\\\ \wp(z_{t}=-1)\mathcal{A}_{|z_{t}=-1}\end{matrix}\right).$ (28)

Here we have introduced a complex conditional coherence,

$\mathcal{A}_{|z_{t}}:=\int e^{i\kappa X}\wp(X|z_{t})dX=\langle e^{i\kappa
X}\rangle_{X|z_{t}},$ (29)

of the data qubit. This quantity is conditioned on (hypothetically) finding
out the value of the RTP $z_{t}$ at time $t$; later we will condition on other
data, that may be hypothetical or actual or a combination of both. We can get
back the no-control complex coherence in Eq. (21) by

$\mathcal{A}^{\rm nc}=\underline{{I}}^{\top}\underline{{A}}^{\rm
nc}=\underline{{I}}^{\top}{\bf
H}(t,\kappa)\underline{{A}}_{0}=\sum_{z_{t}}\wp(z_{t})\mathcal{A}_{|z_{t}},$
(30)

where the sum is over $z_{t}\in\\{-1,1\\}$, showing the connection among
matrix equations and the complex coherence. The absolute value of this complex
coherence is the no-control coherence as shown in Eq. (27).

###  III.3 Decay of coherence in asymptotic regime

Here, we show that, in the long-time limit, $t\rightarrow\infty$, and the
asymptotic regime, where $\kappa$ is much smaller than the RTP jump rate,
_i.e._ , $\kappa\ll\bar{\gamma}$, the no-control coherence $\mathcal{C}^{\rm
nc}$ decays exponentially. We start with the no-control complex coherence,
$\mathcal{A}^{\rm nc}$ in Eq. (21), and calculate the ratio of the change over
an infinitesimal time and its value. That is, the ratio is defined as

$\Lambda_{t}:=\frac{{\rm d}\mathcal{A}^{\rm nc}\big{/}{\rm
d}t}{\mathcal{A}^{\rm nc}},$ (31)

which in general can be a time-dependent function denoted by the subscript
$t$. We use the RTP steady state as the initial vector
$\underline{{A}}_{0}=\underline{{P}}_{t=0}=\underline{{P}}_{\rm ss}$ in (8),
and expand the right-hand side of Eq. (31) to second order in $\kappa$ and
find

$\displaystyle\Lambda_{t}=i\left(\gamma_{\uparrow}-\gamma_{\downarrow}\right)\frac{\kappa}{2\bar{\gamma}}+\left(e^{-2t\bar{\gamma}}-1\right)\frac{\kappa^{2}\breve{\gamma}}{2\bar{\gamma}^{2}}+\mathcal{O}(\kappa^{3}),$
(32)

where $\breve{\gamma}\equiv
2\gamma_{\uparrow}\gamma_{\downarrow}/(\gamma_{\uparrow}+\gamma_{\downarrow})$
is the harmonic mean of $\gamma_{\uparrow}$ and $\gamma_{\downarrow}$. Taking
the long-time limit, _i.e._ , $t\rightarrow\infty$, of above the expansion, we
find that the factor $\Lambda_{t}$ approaches a constant complex value,

$\displaystyle\Lambda=\lim_{t\rightarrow\infty}\Lambda_{t}$
$\displaystyle=-\Gamma^{\rm nc}+i\Omega^{\rm nc},$ (33)

where, in the last line, we have defined the average no-control frequency
(phase drift rate) and decay rate (phase diffusion rate) as

$\Gamma^{\rm
nc}:=\frac{\kappa^{2}\breve{\gamma}}{2\bar{\gamma}^{2}}\quad\text{and}\quad\Omega^{\rm
nc}:=\frac{\kappa(\gamma_{\uparrow}-\gamma_{\downarrow})}{2\bar{\gamma}}.$
(34)

As a result, we substitute $\Lambda_{t}$ in Eq. (31) with the constant
$\Lambda$ and solve the differential equation for the no-control complex
coherence to get

$\mathcal{A}^{\rm nc}(t)=\mathcal{A}^{\rm nc}(0)\exp\left(-\Gamma^{\rm
nc}t+i\Omega^{\rm nc}\,t\right),$ (35)

which has the initial condition $\mathcal{A}^{\rm
nc}(0)=\underline{{I}}^{\top}\underline{{P}}_{\rm ss}=1$. We note that the
complex coherence has its imaginary contribution only when the RTP has
asymmetric jump rates. That is, for the symmetric case,
$\gamma_{\uparrow}=\gamma_{\downarrow}$, we have $\Omega^{\rm nc}=0$. Finally,
since $\mathcal{C}^{\rm nc}=|\mathcal{A}^{\rm nc}|$, we obtain the no-control
coherence in the asymptotic regime,

$\mathcal{C}^{\rm nc}(t)=\exp(-\Gamma^{\rm nc}t).$ (36)

In the CL [28] (Figure 1), we show the comparison between Eq. (36) and the
exact coherence in (27). They differ significantly at short times, this is
expected since Eq. (33) is only valid in the limit $t\gg 1/\bar{\gamma}$.
Moreover, even by $t=1/\bar{\gamma}$ one can see that their decay is almost
identical.

###  III.4 Noise correction and optimal choice of control

The data qubit’s coherence can be improved, from the no-control case, by
utilising additional information gathered by a SQ. Before discussing how to
obtain this information (in the following section), we will show how
additional information, in general, helps increase the coherence of the data
qubit.

Let us assume that we have access to an additional information $Y$ about the
unknown accumulated noise $X$. We can use the information to compute an
appropriate phase correction, $-c(Y)$, and apply it to the data qubit as
described in Sec. III.1, so that $\phi(\Xi)=\phi(X,Y)=\kappa X-c(Y)$. Given
the general form of coherence in Eq. (14), and a control function
$c(\bullet)$, we write a phase-corrected coherence as

$\mathcal{C}_{c(\bullet)}:=\Big{|}\big{\langle}e^{i[\kappa
X-c(Y)]}\big{\rangle}_{X,Y}\Big{|},$ (37)

which is similar to Eq. (20), but now the expected average is over all
possible values of $X$ and $Y$.

If $X$ were known, we could trivially find that $c(Y)=\kappa X$ is a perfect
correction that completely removes the phase error from the noise and the data
qubit’s coherence is maximized at one. However, the noise is actually unknown
and we only have its indirect information from $Y$. Therefore, the main task
is to find the control function, $c(\bullet)$, that maximizes the coherence,
$\mathcal{C}_{c(\bullet)}$ in Eq. (37). To accomplish this, we first write the
expected average in Eq. (37) as an explicit summation over values for $Y$
(assumed discrete) and an integral over the continuous value of $X$ with
appropriate probability weights as

$\displaystyle\mathcal{C}_{c(\bullet)}$ $\displaystyle=\bigg{|}\sum_{Y}\int
e^{i\kappa X}e^{-ic(Y)}\wp(Y)\wp(X|Y){\rm d}X\bigg{|}$
$\displaystyle=\bigg{|}\sum_{Y}e^{-ic(Y)}\wp(Y)\big{\langle}e^{i\kappa
X}\big{\rangle}_{X|Y}\bigg{|}$
$\displaystyle\leq\,\sum_{Y}\wp(Y)\left|\big{\langle}e^{i\kappa
X}\big{\rangle}_{X|Y}\right|,$ (38)

where, in the second line, we have rearranged terms into a conditional
expected average, $\big{\langle}e^{i\kappa X}\big{\rangle}_{X|Y}$. In the last
line, the inequality is, in general, a strict inequality. However, it becomes
an equality when

$c(\bullet)=\arg\big{\langle}e^{i\kappa X}\big{\rangle}_{X|\bullet},$ (39)

which thus maximizes the phase-corrected coherence $\mathcal{C}_{c(\bullet)}$.
Given the formula of the coherence in Eq. (37), we can regard Eq. (39) as
defining the optimal estimator of the real phase, $\phi=\kappa X$, for the
purpose of the final phase correction. Indeed, it is just the mean under
_circular statistics_ [46]. We indicate the maximized _control coherence_
defined as

$\mathcal{C}^{\rm
c}:=\mathcal{C}_{c(\bullet)=\arg\mathcal{A}_{|\bullet}}=\sum_{Y}\wp(Y)\left|\mathcal{A}_{|Y}\right|.$
(40)

Here we are again considering the conditional complex coherence,

$\mathcal{A}_{|Y}:=\langle e^{i\kappa X}\rangle_{X|Y}=\int e^{i\kappa
X}\wp(X|Y){\rm d}X,$ (41)

this time conditioned on the information $Y$.

##  IV Spectator qubit (SQ)
for noise sensing

In our model, the information $Y$ used in the phase correction can be obtained
by measuring the SQ which senses the same RTP noise affecting the data qubit.
The Hamiltonian of the SQ,

$\hat{H}_{\rm s}(t)=\frac{K}{2}\hat{\sigma}_{z}^{\rm s}\,z(t),$ (42)

is identical to that of the data qubit, but with $K$ as the SQ’s noise
sensitivity in place of the much smaller $\kappa$ for the data qubit. Thus the
SQ’s Bloch vector also rotates around the $z$-axis. We assume that we have
complete control over the SQ. To maximize its usefulness as a noise probe, we
take it be reset to an equatorial state after each measurement. Thus its state
can always be written as

$|{\Phi}\rangle^{\rm s}=\frac{1}{\sqrt{2}}\left(|{+1}\rangle_{z}^{\rm
s}+e^{i\Phi}|{-1}\rangle_{z}^{\rm s}\right),$ (43)

given the eigenstates $|{\pm 1}\rangle_{z}^{\rm s}$ of the observable
$\hat{\sigma}_{z}^{\rm s}$, _i.e._ , $\hat{\sigma}_{z}^{\rm s}|{\pm
1}\rangle_{z}^{\rm s}=\pm|{\pm 1}\rangle_{z}^{\rm s}$. The definition in Eq.
(43) should be compared to Eq. (18) for the data qubit.

###  IV.1 SQ’s measurement and likelihood function

To acquire information about the RTP noise, we probe the SQ at various times
throughout the process. Let us denote the measurement times by

$t_{1},t_{2},\dots,t_{n},\dots,t_{N-1},t_{N}=T,$ (44)

where $T$ is the final time that the phase correction is applied on the data
qubit. Given the SQ’s Hamiltonian (42), the integrated noise over time
intervals of duration

$\tau_{n}:=t_{n}-t_{n-1}$ (45)

can be independently probed. Without loss of generality, we initialize the SQ
at the zero-phase state, $|{\Phi=0}\rangle^{\rm s}$, at the initial time,
$t_{0}=0$, and after the $n$th measurement, for $n\in\\{1,N-1\\}$.

Just before the $n$th measurement, at time $t_{n}$, the SQ’s state is given by
$|{\Phi(t_{n})}\rangle^{\rm s}$, where

$\Phi(t_{n})=K\\!\\!\int_{t_{n-1}}^{t_{n}}\\!{\rm
d}t~{}z(t)=K(X_{n}-X_{n-1}):=Kx_{n}.$ (46)

Here we have defined the integrated noise up to time $t_{n}$,

$X_{n}:=\int_{0}^{t_{n}}\\!\\!z(s){\rm d}s,$ (47)

similar to Eq. (16), and defined $x_{n}:=X_{n}-X_{n-1}$, similar to Eq. (9).
We also define $z_{n}:=z(t_{n})$ for the RTP at the measurement time $t_{n}$.

We assume that the measurement of the SQ at time $t_{n}$ is projective. Given
the parameter of interest, $\Phi(t_{n})$, the optimal observable to measure
will be of the form

$\displaystyle\hat{\theta}_{n}:={\mathbb{1}}-|{\theta_{n}}\rangle^{\rm
s}\langle{\theta_{n}}|.$ (48)

Here $\theta_{n}$ is the chosen measurement angle and
$|{\theta_{n}}\rangle^{\rm s}=(1/\sqrt{2})\left(|{+1}\rangle_{z}^{\rm
s}+e^{i\theta_{n}}|{-1}\rangle_{z}^{\rm s}\right)$, and
${\mathbb{1}}=|{+1}\rangle_{z}^{\rm s}\langle{+1}|+|{-1}\rangle_{z}^{\rm
s}\langle{-1}|$ is the identity operator. The outcomes of the projective
measurement,

$y_{n}\in\\{0,1\\},$ (49)

are associated with projecting the SQ’s state to the two eigenstates,
$\\{|{\theta_{n}}\rangle^{\rm s},|{\theta_{n}+\pi}\rangle^{\rm s}\\}$,
respectively. We define $\hat{\theta}_{n}$, Eq. (48), in this way so that the
outcomes have a significance as null and non-null outcomes, respectively, for
a natural choice of $\theta_{n}$, as we will see in Sec. VI.3.

Given the SQ’s state, $|{\Phi(t_{n})}\rangle^{\rm s}=|{Kx_{n}}\rangle^{\rm
s}$, Born’s rule gives the probability of outcome $y_{n}$ as

$\displaystyle\wp(y_{n}|\theta_{n},x_{n})$ $\displaystyle:=\big{|}\,^{\rm
s}\langle\theta_{n}+\pi\,y_{n}|Kx_{n}\rangle^{\rm s}\,\big{|}^{2}$
$\displaystyle=y_{n}+(-1)^{y_{n}}\cos^{2}\left[{\mbox{$\frac{1}{2}$}(\theta_{n}-Kx_{n})}\right].$
(50)

In the context of Bayesian estimation (see next section), this is the
_likelihood function_ , expressed slightly differently as Eq. (9) of the CL
[28].

###  IV.2 Adaptive strategies and the separation principle

Say one has just measured and reprepared the SQ at time $t_{n-1}$. There are
two real parameters that define the next measurement on it: the measurement
angle, $\theta_{n}$, and the waiting time, $\tau_{n}$. Let us combine them and
define—for any $n$—a _measurement setting_

${\mu}_{n}:=\\{\theta_{n},\tau_{n}\\},$ (51)

and define a _measurement strategy_ as

${\mathfrak{S}}=\\{{\mu}_{n}:n\in\\{1,\cdots,N\\}\\}.$ (52)

The aim of our work is to search for measurement strategies, ${\mathfrak{S}}$,
that maximizes the data qubit’s coherence.

In particular, to be general, we must include adaptive strategies, where, at a
local measurement time $t_{n-1}\in\\{t_{1},...,t_{N-1}\\}$, one can use the
information

$Y_{n-1}:=\\{y_{1},y_{2},...,y_{n-1}\\},$ (53)

obtained in the past up to $t_{n-1}$, to find the optimal choice for the next
measurement at $t_{n}$. That is, we can explicitly write a measurement
setting, ${\mu}_{n}(Y_{n-1})$, as a function,

${\mu}_{n}:\\{0,1\\}^{(n-1)}\mapsto U(1)\times\mathbb{R}_{+},$ (54)

that maps measurement results, $Y_{n-1}\in\\{0,1\\}^{(n-1)}$ (in a set of
$(n-1)$-string of binary numbers), to a measurement angle, $\theta_{n}\in
U(1)$ (in a set of the 1D-circular group), and a waiting time,
$\tau_{n}\in\mathbb{R}_{+}$ (in a set of positive real numbers).

From the above, it is apparent that the space of possible strategies
$\left\\{{{\mathfrak{S}}}\right\\}$ is very large for $N$ large. However, we
can simplify the problem of trying to find an optimal strategy by applying the
separation principle [47, 48], as mentioned in the introduction. First, let us
formally define the _expected reward function_ :

$\mathcal{C}^{\rm
c}(T)=\sum_{Y_{N}}\wp(Y_{N})\left|\mathcal{A}_{|Y_{N}}\right|,$ (55)

which is the control coherence, Eq. (40), for the information, $Y=Y_{N}$, at
the final time $T=t_{N}$. It is implicit that the complex conditional
coherence $\mathcal{A}_{|Y_{N}}$, as in Eq. (41), is defined with $X=X_{N}$,
_i.e._ , at the same time ($T=t_{N}$) as $Y_{N}$. We also note that the
summation, $\sum_{Y_{N}}$, is over $Y_{N}\in\\{0,1\\}^{(N)}$ and that
$\wp(Y_{N})=\wp(y_{1},y_{2},...,y_{N})$ is a joint probability of the string
of readouts. Now, the separation principle says that if the reward function is
additive in time—which ours, Eq. (55), is trivially since it is evaluated at
the final time $T$ only—then the optimization problem separates into two
parts.

The first part is optimal estimation of the system. This means finding, at any
time $t_{n}$, the Bayesian probability for the relevant parameters conditioned
on the results obtained so far. That is, one can construct a conditional
distribution $\wp_{n}\equiv\wp(X_{n},z_{n}|Y_{n})$ for the RTP noise ($z_{n}$)
and its accumulated one ($X_{n}$) at time $t_{n}$ conditioned on the current
information $Y_{n}$. The second part is determining the optimal control based
on this information. That is, instead of considering $\mu_{n+1}(Y_{n})$, we
can consider $\mu_{n+1}(\wp_{n})$. Moreover, the goal of the control in the
future of $t_{n}$ (_e.g._ , at the final time $T$) is to maximize the expected
reward function given $Y_{n}$, that is

$\mathcal{C}^{\rm
c}_{|Y_{n}}(T):=\sum_{Y_{N}|Y_{n}}\wp(y_{n+1},...,y_{N}|Y_{n})\left|{\cal
A}_{|Y_{N}}\right|,$ (56)

where the summation is defined as over possible _future_ measurement outcomes,
_i.e._ , $Y_{N}|Y_{n}:=(y_{n+1},\dots,y_{N})$.

Now, given that $\wp_{n}$ is a function of a real variable $X_{n}$ and a
binary variable $z_{n}$, it is not immediately obvious that replacing
$\mu_{n+1}(Y_{n})$ by $\mu_{n+1}(\wp_{n})$ is progress. However, as we will
see in the next section (after considerable calculation), for the case of our
particular conditional expected reward function Eq. (56), only a couple of
statistics derived from the distribution $\wp(X_{n},z_{n}|Y_{n})$ are relevant
(see Section V.4), which greatly simplifies the problem.

Before turning to the details of that calculation, we here also introduce a
more general definition of the conditional control coherence:

$\mathcal{C}^{\rm
c}_{|D_{n}}(t_{m}):=\sum_{Y_{m}}\wp(Y_{m}|D_{n})\left|\langle e^{i\kappa
X_{m}}\rangle_{X_{m}|Y_{m}}\right|,$ (57)

which will be of use in later sections. Here the generalizations beyond (57)
are two-fold. First, we allow conditioning on an arbitrary set of data,
$D_{n}$, which is available at measurement time $t_{n}$, not necessarily equal
to the full record $Y_{n}$ up to that time. Second, the time at which we
imagine implementing the control on the data qubit is any measurement time
$t_{m}$, with $m\geq n$, rather than necessarily the final time $t_{N}=T$ as
in Eq. (56).

##  V Bayesian maps
for phase estimation

In this section, we formulate a Bayesian map, which is central to our results,
especially the MOAAAR algorithm of Sec. VII. We start with the expected reward
function, Eq. (55), and show that it can be written as a matrix equation,
similar to that of the no-control case, with the map ${\bf H}$ in Eq. (27). In
contrast to the no-control case, where there is only one map for the entire
process, the coherence must be updated whenever new information from the SQ
measurements at $t\in\\{t_{1},t_{2},...,t_{N}\\}$ are obtained. We will see
that the new map, denoted by ${\bf
F}(\mu_{n}=\\{\theta_{n},\tau_{n}\\},y_{n})$, is an explicit function of the
measurement setting, ${\mu}_{n}$, and the readout, $y_{n}$, and can be derived
analytically using Bayes’ theorem.

###  V.1 Bayes’ theorem for the expected reward function

The expected reward function, defined in Eq. (55), can be explicitly written
as

$\displaystyle\mathcal{C}^{\rm c}(T)$
$\displaystyle=\sum_{Y_{N}}\wp(Y_{N})\left|\mathcal{A}_{|Y_{N}}\right|,$
$\displaystyle=\sum_{Y_{N}}\wp(Y_{N})\left|\int{\rm
d}X_{N}~{}\wp(X_{N}|Y_{N})e^{i\kappa X_{N}}\right|,$
$\displaystyle=\sum_{Y_{N}}\left|\int{\rm d}X_{N}~{}\wp(Y_{N},X_{N})e^{i\kappa
X_{N}}\right|,$ (58)

where $X_{N}$ is the final accumulated noise,

$X_{N}=x_{1}+x_{2}+...+x_{N},$ (59)

and $Y_{N}$ is a string of measurement results from the SQ,

$\displaystyle Y_{N}=(y_{1},y_{2},...,y_{N}),$ (60)

up until the final time.

The integral in Eq. (V.1) is not trivial since we do not know an explicit form
of the joint probability function $\wp(Y_{N},X_{N})$. However, we have the
likelihood functions $\wp(y_{n}|\theta_{n},x_{n})$ and
$\wp(x_{n},z_{n}|z_{n^{\prime}})$, via Eqs. (IV.1) and (II) respectively.
Thus, we write Eq. (V.1) in terms of these likelihood functions, by first
defining a noise vector, $\vec{x}:=\\{x_{1},\dots,x_{N}\\}$, and its 1-norm
$||\vec{x}||_{1}=\sum_{n}x_{n}$. Using Eq. (59) allows us to write

$\displaystyle\wp(Y_{N},X_{N})$
$\displaystyle=\int\\!\\!\\!\wp(Y_{N},\vec{x})\delta\big{(}X_{N}-||\vec{x}||_{1}\big{)}{\rm
d}^{N}\vec{x},$
$\displaystyle=\int\\!\\!\\!\wp(Y_{N}|\vec{x})\wp(\vec{x})\delta\Big{(}X_{N}-||\vec{x}||_{1}\Big{)}{\rm
d}^{N}\vec{x}.$ (61)

The integral measure is ${\rm d}^{N}\vec{x}\equiv{\rm d}x_{1}\cdots{\rm
d}x_{N}$. Substituting the above joint distribution back into Eq. (V.1) and
integrating over $X_{N}$, we obtain

$\displaystyle\mathcal{C}^{\rm
c}(T)=\sum_{Y_{N}}\left|\int\wp(Y_{N}|\vec{x})\wp(\vec{x})e^{i\kappa||\vec{x}||_{1}}{\rm
d}^{N}\vec{x}\right|.$ (62)

For the joint distribution, $\wp(Y_{N}|\vec{x})$, since we reset the SQ after
every measurement, each result $y_{n}$ depends only on the accumulated noise
$x_{n}$ during the waiting time $\tau_{n}$. Thus, we can write

$\displaystyle\wp(Y_{N}|\vec{x})$
$\displaystyle\equiv\wp(y_{1},...,y_{N}|x_{1},...,x_{N}),$
$\displaystyle=\prod_{n=1}^{N}\wp_{{\mu}_{n}}(y_{n}|x_{n}),$ (63)

where we have added the measurement setting
${\mu}_{n}=\\{\theta_{n},\tau_{n}\\}$ as a subscript of the likelihood
function (IV.1). Furthermore, we construct a joint probability of the
accumulated noises, $x_{1},...,x_{N}$, with the help of RTP variables and Eq.
(II), and write,

$\displaystyle\wp(\vec{x})$ $\displaystyle\equiv\wp(x_{1},...,x_{N}),$
$\displaystyle=\sum_{z_{0},...,z_{N}}\left(\prod_{n=1}^{N}\wp(x_{n},z_{n}|z_{n-1})\right)\wp(z_{0}),$
(64)

as a marginalized probability function over the RTP noises, $z_{0},...,z_{N}$,
at all times. Here the sum $\sum_{z_{0},...,z_{N}}$ is over all possible
values of the RTP. We then substitute Eqs. (V.1) and (V.1) into the reward
function in Eq. (62) and get

$\displaystyle\mathcal{C}^{\rm
c}(T)=\sum_{Y_{N}}\wp(Y_{N})\left|\mathcal{A}_{|Y_{N}}\right|$
$\displaystyle=\,\,\sum_{Y_{N}}\left|\sum_{z_{0},...,z_{N}}\left(\prod_{n=1}^{N}\int
dx_{n}~{}\wp_{{\mu}_{n}}(y_{n}|x_{n})\wp(x_{n},z_{n}|z_{n-1})e^{i\kappa
x_{n}}\right)\wp(z_{0})\right|,$
$\displaystyle=\,\,\sum_{Y_{N}}\left|\sum_{z_{0},...,z_{N}}\left(\prod_{n=1}^{N}F_{z_{n-1}}^{z_{n}}\big{(}{\mu}_{n},y_{n}\big{)}\right)\wp(z_{0})\right|,$
(65)

where in the second line we have defined

$\displaystyle
F_{z_{n-1}}^{z_{n}}\\!\left(\mu_{n},y_{n}\right):=\\!\\!\int\\!\\!dx_{n}\,\wp_{{\mu}_{n}}(y_{n}|x_{n})\,\wp(x_{n},z_{n}|z_{n-1})e^{i\kappa
x_{n}}.$ (66)

These are the linear maps that are the core elements of our Bayesian-map
formalism. In the next subsection, we will derive an analytic expression for
$F_{z_{n-1}}^{z_{n}}\\!\left(\mu_{n},y_{n}\right)$, and rewrite it explicitly
as a matrix.

###  V.2 Analytical expression for $F$-map

In the similar manner as in the construction of the matrix $\mathbf{H}$ in Eq.
(26), we see that the sums and products in Eq. (65) can be rewritten as matrix
multiplication. Specifically, defining the two-by-two matrix

${\bf
F}\left(\mu_{n},y_{n}\right):=\left(\begin{matrix}F^{+1}_{+1}\left(\mu_{n},y_{n}\right)&F^{+1}_{-1}\left(\mu_{n},y_{n}\right)\\\\[5.0pt]
F^{-1}_{+1}\left(\mu_{n},y_{n}\right)&F^{-1}_{-1}\left(\mu_{n},y_{n}\right)\end{matrix}\right),$
(67)

we can write the expected reward function Eq. (65) as

$\displaystyle\mathcal{C}^{\rm c}(T)$
$\displaystyle=\,\,\sum_{Y_{N}}\big{|}\underline{I}^{\top}\,{\bf
F}({\mu}_{N},y_{N})\,\cdots{\bf F}({\mu}_{2},y_{2}){\bf
F}({\mu}_{1},y_{1})\underline{A}_{0}\big{|}$
$\displaystyle=\sum_{Y_{N}}\big{|}\underline{I}^{\top}\,\underline{A}_{N}\big{|},$
(68)

where $\underline{A}_{0}$ and $\underline{I}$ are defined in Eqs. (24) and
(27). Also, in the second line of Eq. (V.2), we defined a new coherence vector

$\underline{{A}}_{N}\equiv{\bf F}({\mu}_{N},y_{N})\,\cdots{\bf
F}({\mu}_{2},y_{2}){\bf F}({\mu}_{1},y_{1})\underline{A}_{0}.$ (69)

Moreover, we can use convolutions and Fourier transform (see details in the
Appendix B) to find that the matrices, $\mathbf{F}$ and $\mathbf{H}$, are
related via

$\displaystyle\mathbf{F}\big{(}\mu_{n}{=}\\{\theta_{n},\tau_{n}\\},y_{n}\big{)}=$
$\displaystyle\frac{1}{4}\Big{[}2\,\mathbf{H}(\tau_{n},\kappa)$
$\displaystyle+(-1)^{y_{n}}e^{-i\theta_{n}}\mathbf{H}(\tau_{n},\kappa+K)$
$\displaystyle+(-1)^{y_{n}}e^{+i\theta_{n}}\mathbf{H}(\tau_{n},\kappa-K)\Big{]}.$
(70)

One can easily confirm from Eq. (V.2) that

$\displaystyle\sum_{y_{n}}\mathbf{F}({\mu}_{n},y_{n})=\mathbf{H}(\tau_{n},\kappa),$
(71)

which can be understood that the no-control map $\mathbf{H}$ (no measurement
information) is equivalent to ignoring the gained information from the
conditional map $\mathbf{F}({\mu}_{n},y_{n})$.

###  V.3 Interpretation of complex 2-vector $\underline{A}_{n}$

Following the definition of $\underline{{A}}_{N}$ in Eq. (69), we define a
coherence vector $\underline{A}_{n}=(A_{n}^{z=+1},A_{n}^{z=-1})^{\top}$ for an
intermediate time $t_{n}$. This can be computed from the vector at the
previous time step using the measurement result $y_{n}$ via

$\underline{A}_{n}={\bf F}({\mu}_{n},y_{n})\underline{A}_{n-1}.$ (72)

We can interpret the two elements of $\underline{{A}}_{n}$, similar to Eq.
(28), as

$\underline{A}_{n}=\left(\begin{matrix}\underline{{A}}_{n}^{+1}\\\
\underline{{A}}_{n}^{-1}\end{matrix}\right)=\left(\begin{matrix}\wp(Y_{n},z_{n}=+1){\cal
A}_{|Y_{n},z_{n}=+1}\\\ \wp(Y_{n},z_{n}=-1){\cal
A}_{|Y_{n},z_{n}=-1}\end{matrix}\right).$ (73)

That is, the coherence vector $\underline{{A}}_{n}$ contains information about
the measurement results $Y_{n}=(y_{1},...,y_{n})$ and the complex conditional
coherence defined as

${\cal A}_{|Y_{n},z_{n}}:=\langle e^{i\kappa
X_{n}}\rangle_{X_{n}|Y_{n},z_{n}},$ (74)

conditioned on both the results and the hypothetically knowable $z_{n}$, at
time $t_{n}$. With this, we can see a relationship between two types of
complex conditional coherences, in Eq. (41) and Eq. (74), via

$\wp(Y_{n})\mathcal{A}_{|Y_{n}}=\underline{{I}}^{\top}\underline{{A}}_{n}=\sum_{z_{n}}\wp(Y_{n},z_{n})\mathcal{A}_{|Y_{n},z_{n}},$
(75)

where the summation is over the RTP state $z_{n}\in\\{-1,+1\\}$.

The vector $\underline{{A}}_{n}$ does not let one calculate exactly the
probability of the measurement results, $\wp(Y_{n})$. However, if one
considers $\kappa$ (the sensitivity of the data qubit to noise) to be a
variable rather than a fixed number, then one can take the limit $\kappa\to 0$
in Eq. (74). This gives $\lim_{\kappa\rightarrow 0}{\cal A}_{|Y_{n},z_{n}}=1$
and, from Eq. (73), we see that

$\wp(Y_{n})=\lim_{\kappa\rightarrow
0}\underline{{I}}^{\top}\underline{{A}}_{n}.$ (76)

This motivates defining a probability vector,

$\underline{{\check{A}}}_{n}=\lim_{\kappa\rightarrow
0}\underline{{A}}_{n}=\left(\begin{matrix}\wp(Y_{n},z_{n}=+1)\\\
\wp(Y_{n},z_{n}=-1)\end{matrix}\right),$ (77)

which should not be confused with the unconditioned RTP probability vector
$\underline{{P}}_{t}$.

Similar to the coherence vector $\underline{{A}}_{n}$ in Eq. (72), the
probability vector can be updated after every SQ measurement via

$\underline{{\check{A}}}_{n}=\check{\bf
F}({\mu}_{n},y_{n})\underline{{\check{A}}}_{n-1},$ (78)

where we have defined a new probability map with $\kappa\rightarrow 0$,

$\check{\bf F}({\mu}_{n},y_{n}):=\lim_{\kappa\to 0}{\bf F}({\mu}_{n},y_{n}).$
(79)

This $\check{\bf F}$ matrix contains only real numbers and thus can map
between two real-number vectors containing only probability functions.
Interestingly, the initial state of the probability vector is the same as
$\underline{{A}}_{0}$ and the RTP initial state, $\underline{{P}}_{0}$, _i.e._
,

$\displaystyle\underline{{\check{A}}}_{0}=\underline{{A}}_{0}=\underline{{P}}_{0}=\underline{{P}}_{\rm
ss}.$ (80)

Here we have also set them to be equal to the RTP steady state as per Eq. (8).
We will see later that $\underline{{\check{A}}}_{n}$ and $\wp(Y_{n})$ are
useful to both understanding dynamics of the SQ and simulating the data
qubit’s coherence numerically.

###  V.4 Two real sufficient statistics: $(\alpha,\zeta)$

As mentioned in Section IV.2, the separation principle allows us to separate
the problem at any time $t_{n}$ into the estimation part, using the
information $Y_{n}$ obtained up to that time, and the control part, which uses
that estimate to maximize the expected future cost function, conditioned on
the present $Y_{n}$. Here we show that, for our cost function, the estimation
problem reduces to calculating the complex 2-vector $\underline{{A}}_{n}$.
Furthermore, except at the final time $T$, just two real numbers derived from
$\underline{{A}}_{n}$ are sufficient for the control problem.

We start with the expected reward function conditioned on $Y_{n}$ given in Eq.
(56), which we rewrite here for convenience:

$\mathcal{C}^{\rm
c}_{|Y_{n}}(T):=\sum_{Y_{N}|Y_{n}}\wp(y_{n+1},...,y_{N}|Y_{n})\left|{\cal
A}_{|Y_{N}}\right|.$ (81)

We then use Bayes’ rule to write
$\wp(y_{n+1},...,y_{N}|Y_{n})=\wp(Y_{N})/\wp(Y_{n})$ and use the relationship
$\wp(Y_{N})|{\cal A}_{|Y_{N}}|=|\underline{I}^{\top}\,\underline{A}_{N}|$ from
Eq. (75) to obtain

$\displaystyle\mathcal{C}^{\rm c}_{|Y_{n}}(T)$ $\displaystyle=$
$\displaystyle[\wp(Y_{n})]^{-1}\sum_{Y_{N}|Y_{n}}\big{|}\underline{I}^{\top}\,\underline{A}_{N}\big{|},$
(82a) $\displaystyle=$
$\displaystyle[\wp(Y_{n})]^{-1}\sum_{Y_{N}|Y_{n}}\big{|}\underline{I}^{\top}\,{\bf
F}({\mu}_{N},y_{N})\,\cdots$ $\displaystyle\hskip 28.45274pt{\bf
F}({\mu}_{n+2},y_{n+2}){\bf F}({\mu}_{n+1},y_{n+1})\underline{A}_{n}\big{|},$

where in the second line we used Eq (69), only expanding out the $\bf F$ maps
for the future outcomes $y_{n+1},...,y_{N}$.

Now, to maximize this reward function, the only things we have the ability to
control, in the future of $t_{n}$, in Eq. (82) is the SQ measurement strategy,
${\mu}_{n+1},...,{\mu}_{N}$. What current information is relevant to this
strategy? First, since $Y_{n}$ is already known, $\wp(Y_{n})$ is a fixed
scalar multiple, and so is not relevant. That leaves only the vector
$\underline{{A}}_{n}$. In other words, if we know $\underline{{A}}_{n}$, we
know everything relevant for choosing the SQ measurement strategy to maximize
the final objective $\mathcal{C}^{\rm c}_{|Y_{n}}(T)$. From Eq. (40), it also
contains all the relevant present information for calculating, at the final
time, the optimal control to apply to arrive at Eq. (82). Here we are
concerned with the question ‘what parts of $\underline{{A}}_{n}$ are relevant
for the SQ measurement strategy?’ The complex-2 vector $\underline{{A}}_{n}$
defines four real parameters. It turns out that only two of them are relevant
as we now show.

Since the two (complex) elements of $\underline{{A}}_{n}$ are associated with
two plausible values of the unknown RTP, _i.e._ , $z_{n}\in\\{-1,+1\\}$ as in
Eq. (73), one might guess that these elements should tell us something
relevant about the data qubit’s phase and coherence if one could,
hypothetically, find out which state $z_{n}$ the RTP is actually in. One such
quantity is the scaled ratio of the two elements of $\underline{{A}}_{n}$,

$\displaystyle\alpha_{n}:$
$\displaystyle=\frac{K}{\kappa}\arg\frac{A_{n}^{z=+1}}{A_{n}^{z=-1}},$ (83a)
$\displaystyle=\frac{K}{\kappa}\left({\arg\mathcal{A}_{|Y_{n},+1}-\arg\mathcal{A}_{|Y_{n},-1}}\right),$
(83b)

using the definition in Eq. (73). This quantity $\alpha_{n}$ is the first of
the two sufficient statistics and can be interpreted as the _maximum
difference of an optimal control phase_ from finding out if the RTP state was
$z_{n}=+1$ or $z_{n}=-1$, scaled so as to be of order unity.

To understand Eq. (83a), we recall that the optimal control (phase correction)
is given by $c(Y)=\arg\mathcal{A}_{|Y}$ from Eq. (39). If, at a current time
$t_{n}$, a phase correction, $c(Y_{n})$, were to be applied to the data qubit,
by further assuming the hypothetical $z_{n}$, one could consider refining the
control to

$c^{\prime}(Y_{n},z_{n})=\arg\mathcal{A}_{|Y_{n},z_{n}}=\arg\left\langle
e^{i\kappa X_{n}}\right\rangle_{{X_{n}}|Y_{n},z_{n}},$ (84)

where $z_{n}$ is included in the condition of the expected average. The
bracketed term in Eq. (83b) is therefore the difference of the control phases
for the two hypothetical $z_{n}$. To understand the scaling factor that makes
this phase difference of order unity, $O(1)$, it is necessary to consider the
dynamics of the RTP and the measurements, as we do in the following paragraph.

The value of $z_{n}$ is relevant to the data qubit phase only for the recent
period of evolution, of order $1/\breve{\gamma}$. For any time earlier than
that, the RTP state is pretty much uncorrelated with its current state. Thus
the maximum difference that finding out $z_{n}$ could make on the data qubit’s
phase estimate is the difference in the cumulative phase, $\kappa X_{n}$, over
that time scale, for $z_{n}=+1$ versus $z_{n}=-1$. That is, of order
$\kappa/\breve{\gamma}$, which is small. In actuality, when we are gathering
information about $z$ through a good adaptive protocol (see later sections)
the difference is even smaller, because from the preceding measurement result
we have a very good idea of the value of $z$. Hence if $z$ is unknown that is
substantially because it may have changed its value since the last
measurement. In other words, the preceding time scale for which finding out
the value of $z$ is relevant is of order $\tau$, not of order
$1/\breve{\gamma}$. As we will see later, the optimal value for $\tau$ is of
order $1/K$. Thus, the difference that finding out $z_{n}$ could make for the
phase estimate is of order $\kappa/K$. Hence, to scale the phase difference to
the order unity, we divide it by the factor $\kappa/K$ and get $\alpha_{n}$ as
in Eq. (83a).

The second of the two sufficient statistics is related to the modulus of the
two elements of $\underline{{A}}_{n}$, defined as

$\zeta_{n}:=\frac{|A_{n}^{z=+1}|-|A_{n}^{z=-1}|}{|A_{n}^{z=+1}|+|A_{n}^{z=-1}|}.$
(85)

This can be interpreted as an _approximated mean of $z_{n}$ conditioned on
$Y_{n}$_. To see this, we consider the conditional complex coherence,
$\mathcal{A}_{|Y_{n},\pm 1}$. As we have already seen above, their arguments
are very close, differing by only $O(\kappa/K)$ as per Eq. (83b). But this
implies that their moduli must also be very close. In fact, the relative
difference in their moduli is quadratically smaller than the difference in
their arguments. That is because the relative change in the data qubit
coherence $|\mathcal{A}|$ over some short time interval is quadratic in the
uncertainty in the data qubit’s phase. (One can see this by Taylor expanding
the decoherence to the second order of the phase.) But the difference in the
uncertainty in the data qubit’s phase, for the different values of $z_{n}$,
will not be larger than the corresponding difference in their mean, because
the difference between the two RTP values is the maximum uncertainty the RTP
value can have.

Therefore, in the relevant time for the refined best-estimate phase shift,
$\tau=O(1/K)$, the difference in the relative decoherence resulting from
refining one’s knowledge to $z=+1$ versus $z=-1$ would be at most of order
$O((\kappa/K)^{2})$. That is, for a good measurement protocol, we can be
confident that
$|\mathcal{A}_{|Y_{n},+1}|/|\mathcal{A}_{|Y_{n},-1}|=1+O\big{(}(\kappa/K)^{2}\big{)}$.
With relative errors of the same magnitude, we can thus replace
$|\mathcal{A}_{|Y_{n},z_{n}}|$ by $|\mathcal{A}_{|Y_{n}}|$ in Eq. (73), to get

$|A_{n}^{z_{n}}|\approx|{\cal A}_{|Y_{n}}|\,\wp(Y_{n},z_{n}).$ (86)

Hence, $\zeta_{n}$ in Eq. (85) can be approximated as

$\zeta_{n}\approx\frac{\sum_{z_{n}}z_{n}\,\wp(Y_{n},z_{n})}{\sum_{z_{n}}\wp(Y_{n},z_{n})},$
(87)

which is exactly the mean of $z_{n}$ conditioned on the current measurement
record, $Y_{n}$.

There are two more (third and fourth) parameters from the vector
$\underline{{A}}_{n}$ that are not relevant to finding the optimal SQ control.
The third parameter is the argument of the sum of elements, _i.e._ ,

$\varphi_{n}:=\arg\mathcal{A}_{|Y_{n}}=\arg\left(A_{n}^{z=+1}+A_{n}^{z=-1}\right).$
(88)

This would be relevant as the final control on the data qubit, $c(Y_{n})$, if
$t_{n}=T$ were the final time. However, because $t_{n}$ is not the end of the
protocol, the quantity simply expresses how far the mean phase has drifted so
far. Since the random phase driven by the RTP is a true $U(1)$ process, having
an estimate of the phase at any particular time is not relevant to future
decoherence. One can see this from Eq. (82) that a global phase change to the
vector $\underline{{A}}_{n}$, _i.e._ , changing the value of $\varphi_{n}$,
does not affect the coherence on the left of the equation.

Fourthly, we can define the final real parameter from the sum of moduli of the
elements as

$r_{n}:=|A_{n}^{z=+1}|+|A_{n}^{z=-1}|\approx\wp(Y_{n})|\mathcal{A}_{|Y_{n}}|,$
(89)

where we have used Eq. (86) with the same degree of approximation as discussed
around that equation. As already argued, $\wp(Y_{n})$ is irrelevant to
optimizing the future control, and so is any loss of coherence already
suffered such that $|\mathcal{A}_{|Y_{n}}|<1$, by a similar argument used for
$\varphi_{n}$ above.

Parameter |  Description
---|---
$r_{n}$, Eq. (89) |  $\approx$ $Y_{n}$-probability $\times$ conditional coherence.
$\varphi_{n}$, Eq. (88) |  $=$ optimal control phase, not knowing $z_{n}$.
$\zeta_{n}$, Eq. (85) |  $\approx$ mean of $z_{n}$ conditioned on $Y_{n}$.
$\alpha_{n}$, Eq. (83a) |  $=$ maximum difference that finding out $z_{n}$ could make to optimal control phase, scaled.
Table 1: We present four (real) statistics parameters along with their
physical explanations in this table. The vector $\underline{{A}}_{n}$ is
described by these parameters.

We summarize the meaning of the four real parameters that make up the vector
$\underline{{A}}_{n}$ in table 1. Another way to summarize them is to note
that in the regime of interest it is possible to approximate
$\underline{{A}}_{n}$ by a simple expression using the above four parameters,

${A}_{n}^{z=\pm 1}\approx
r_{n}\exp(i\varphi_{n})\times\frac{1\pm\zeta_{n}}{2}\exp\left(\pm
i\frac{1\mp\zeta_{n}}{2}\frac{\alpha_{n}\kappa}{K}\right).$ (90)

Here the approximation sign allows for a relative error in $r_{n}$ of
$O\left((\kappa/K)^{2}\right)$ and an absolute error in $\varphi_{n}$ of
$O\left((\kappa/K)^{3}\right)$. This expression is for interest only; we do
not use it in the remainder of the paper.

Returning to the two sufficient statistics, $\alpha_{n}\in\mathbb{R}$ in Eq.
(83a) and $\zeta_{n}\in[-1,1]$ in Eq. (85), these can thus be regarded as the
only two moments of $\wp(X_{n},z_{n}|Y_{n})$ (the solution to the estimation
part of the problem) that are relevant to the SQ measurement part of the
problem. Thus, as mooted in Sec. IV.2, we can vastly reduce in size the set of
control protocols $\mathfrak{S}$ Eq. (52) that we consider by simplifying the
functional form of $\mu_{n+1}$ from that in Eq. (54) to

${\mu}_{n+1}:\mathbb{R}\times[-1,1]\mapsto U(1)\times\mathbb{R}_{+}.$ (91)

That is, $\mu_{n+1}$ is now a function of $(\alpha_{n},\zeta_{n})$ rather than
of $Y_{n}$ directly. Finding the exactly optimal control strategy
$\left\\{{\mu_{n}:n\in\\{1,\cdots,N\\}}\right\\}$ is still a difficult task,
especially as, in general, $N$ will not, in fact, be fixed for a fixed $T$.
However, if we seek a close-to-optimal, rather than exactly optimal, control
strategy ${\mathfrak{S}}$, then this difficulty can be avoided. Because we are
interested in the long time limit, we will have $N\gg 1$, and for the vast
majority of the control choices we can ignore the _edge effects_ (where $n$ is
small or $n$ is close to $N$), and choose the $\mu_{n+1}$ to all be the same
function of their two arguments. That is, we have only to optimize a single
function,

${\mu}:\mathbb{R}\times[-1,1]\mapsto U(1)\times\mathbb{R}_{+}.$ (92)

In the future sections, we will make further assumptions to restrict the
function $\mu$, and we will see how the stochastic behaviour of $\alpha_{n}$
and $\zeta_{n}$ is greatly useful in developing the intuition to do this.

##  VI Local optimisation algorithm (Greedy)

The first strategy we explore is a ‘greedy’ one [49]. That is, it based on the
natural idea of maximizing the expected reward function $\mathcal{C}^{\rm c}$
locally in time. This is easier than maximizing the actual reward function,
which is $\mathcal{C}^{\rm c}$ at the final time $T$. Greedy algorithms are
not optimal for typical hard problems, but even if they are not, they may be
relatively close to optimal. That appears to be the situation for our problem,
as we will find. We call this strategy the Greedy algorithm, or simply
‘Greedy’ for short.

We conceptualize the Greedy algorithm as follows. Say that we are at some time
$t$ such that $t_{n}<t\leq t_{n+1}$. By this we mean that exactly $n$
measurements on the SQ have already happened and their measurement results,
$Y_{n}$, are known, and it is still possible to decide to make the $(n+1)$th
measurement at time $t$. Imagine now that the final time is the immediate
future,

${\mathfrak{T}}=t+{\rm d}t.$ (93)

Then the Greedy algorithm aims to maximize the conditional expected reward
function $\mathcal{C}^{\rm c}_{|Y_{n}}({\mathfrak{T}})$ as in Eq. (82a) where
we have replaced $T$ with ${\mathfrak{T}}$.

We emphasize that any—even random—measurement on the SQ provides some
information that can be used to increase the coherence of the data qubit. Thus
to maximize $\mathcal{C}^{\rm c}_{|Y_{n}}({\mathfrak{T}})$ at least one
measurement on the SQ at time ${\mathfrak{T}}=t+{\rm d}t$ should be performed.
Greedy maximizes the expected reward function at the immediate future time, by
deciding whether or not to measure at time $t$ and optimizing the measurement
angles at both times. In the following subsections, we present the details of
Greedy in Sec. VI.1 and then demonstrate its numerical implementation in Sec.
VI.2. Studying the patterns in the numerically found Greedy algorithm in Sec.
VI.3 then motivate us to propose and analyze various ‘experimentally’ simpler
greedy algorithms in Secs. VI.4 and VI.6, where in Sec.VI.5 we analyze the
phase space dynamics of the simpler greedy algorithms.

###  VI.1 The Greedy algorithm to maximize $\mathcal{C}^{\rm
c}_{|Y_{n}}({\mathfrak{T}})$

Let the current time be denoted $t=t_{n}+\tau$, with $\tau>0$. The goal is to
determine the next measurement setting, ${\mu}_{n+1}$, based on the current
information, $Y_{n}$. To do that, the Greedy algorithm compares the reward
$\mathcal{C}^{\rm c}_{|Y_{n}}({\mathfrak{T}})$ from two prospective scenarios:

  1. (i)

The SQ is measured only at time ${\mathfrak{T}}=(t_{n}+\tau)+{\rm d}t$,

  2. (ii)

The SQ is measured now, at $t=t_{n}+\tau$, and again at
${\mathfrak{T}}=(t_{n}+\tau)+{\rm d}t$.

We then calculate the coherence vectors associated with the two scenarios at
time ${\mathfrak{T}}$ using equation (72):

$\displaystyle\underline{A}^{\rm(i)}$ $\displaystyle:=\mathbf{F}({\mu}_{(\rm
i)},y_{n+1})\underline{A}_{n},$ (94a) $\displaystyle\underline{A}^{\rm(ii)}$
$\displaystyle:=\mathbf{F}({\mu}_{(\rm
ii)}^{\prime},y_{n+2})\mathbf{F}({\mu}_{(\rm ii)},y_{n+1})\underline{A}_{n},$
(94b)

where $y_{n+1}$ and $y_{n+2}$ are measurement results at time $t$ and
${\mathfrak{T}}=t+{\rm d}t$ respectively. Scenario (i) consists of only one
measurement with a setting ${\mu}_{(\rm i)}:=\\{\theta_{\rm(i)},\tau+{\rm
d}t\\}$ whereas scenario (ii) consists of two consecutive measurements with
settings ${\mu}_{(\rm ii)}:=\\{\theta_{\rm(ii)},\tau\\}$ and ${\mu}_{(\rm
ii)}^{\prime}:=\\{\theta_{\rm(ii)}^{\prime},{\rm d}t\\}$ for the first and
second measurement respectively. Note that these angles should be thought of
not as numbers but as functions of $\underline{A}_{n}$ and, for
$\theta_{\rm(ii)}^{\prime}$, also of $y_{n+1}$.

Despite the fact that scenario (ii) has one more measurement than scenario
(i), the coherence vectors are both calculated at the hypothetical final time
${\mathfrak{T}}=t+{\rm d}t$. Thus, the reward functions—as it is defined in
Eq. (82a)— for the two cases at time ${\mathfrak{T}}$ are

$\displaystyle\mathcal{C}^{\rm c}_{\rm(i)}(\theta_{\rm(i)}):$
$\displaystyle=\mathcal{C}^{\rm
c}_{|Y_{n}}({\mathfrak{T}})\Big{|}_{\text{(i)}}$
$\displaystyle=[\wp(Y_{n})]^{-1}\sum_{y_{n+1}}\left|\underline{I}^{\top}\underline{A}^{\rm(i)}\right|,$
(95a) $\displaystyle\mathcal{C}^{\rm
c}_{\rm(ii)}(\theta_{\rm(ii)},\theta_{\rm(ii)}^{\prime}):$
$\displaystyle=\mathcal{C}^{\rm
c}_{|Y_{n}}({\mathfrak{T}})\Big{|}_{\text{(ii)}}$
$\displaystyle=[\wp(Y_{n})]^{-1}\sum_{y_{n+2},y_{n+1}}\left|\underline{I}^{\top}\underline{A}^{\rm(ii)}\right|,$
(95b)

where on the left hand sides of the above equations, we defined
$\mathcal{C}^{\rm c}_{\rm(i)}$ and $\mathcal{C}^{\rm c}_{\rm(ii)}$ as
functions of the prospective angles. We then find the optimal angles by
maximizing both $\mathcal{C}^{\rm c}_{\rm(i)}$ and $\mathcal{C}^{\rm
c}_{\rm(ii)}$ over all possible angles. Formally, the maximized reward
functions are given by

$\displaystyle\mathcal{C}^{\rm op}_{\rm(i)}$
$\displaystyle:=\max_{\theta_{\rm(i)}}~{}\mathcal{C}^{\rm
c}_{\rm(i)}(\theta_{\rm(i)}),$ (96a) $\displaystyle\mathcal{C}^{\rm
op}_{\rm(ii)}$
$\displaystyle:=\max_{\theta_{\rm(ii)},\theta_{\rm(ii)}^{\prime}}\mathcal{C}^{\rm
c}_{\rm(ii)}(\theta_{\rm(ii)},\theta_{\rm(ii)}^{\prime}).$ (96b)

We compare the optimum values in Eqs. (96) to decide whether or not to measure
at time $t$. The decision-making procedure is as follows: if $\mathcal{C}^{\rm
op}_{\rm(i)}>\mathcal{C}^{\rm op}_{\rm(ii)}$, it means that no measurement at
time $t$ is required to maximize the reward function at ${\mathfrak{T}}$.
Thus, Greedy decides to wait and moves to the next infinitesimal time to
repeat the procedure again. If, on the other hand, $\mathcal{C}^{\rm
op}_{\rm(i)}\leq\mathcal{C}^{\rm op}_{\rm(ii)}$, Greedy chooses to measure the
SQ at $t$, so that $\tau$ is the actual waiting time $\tau_{n+1}$ from $t_{n}$
to $t_{n+1}=t$ and the measurement angle $\theta_{\rm(ii)}$ is the actual
$\theta_{n+1}$. Thus, we write

$\displaystyle{\mu}_{n+1}=\\{\theta_{\rm(ii)},\tau\\},$ (97)

as the next measurement setting chosen by the Greedy. We note that, the pre-
factor $[\wp(Y_{n})]^{-1}$ in Eq. (95) does not influence the Greedy decision-
making procedure and knowing only $\underline{{A}}_{n}$ suffices. Indeed,
following the argument of Sec. V.4, knowing only the two real parameters
$\alpha_{n}$ and $\zeta_{n}$ suffices.

Now that we have identified ${\mu}_{n+1}$, we repeat the procedure, starting
at time $t_{n+1}$, to find the measurement settings for the next step, _i.e._
, ${\mu}_{n+2}$, using $\underline{{A}}_{n+1}$. Repeating in this fashion till
$t_{N}=T$ (where we always measure), equips us with a single realization of
$Y_{N}$ that corresponds to a specific sequence of settings,
$\\{{\mu}_{1},\dots,{\mu}_{N}\\}$ and a trajectory of the reward function,
$\mathcal{C}^{\rm c}_{|Y_{n}}(t_{n})$ for $n\in\\{1,\dots,N\\}$. Then, to
calculate the expected reward function, similar to what we defined in Eq.
(V.1), we average over the information $Y_{n}$ as

$\displaystyle\mathcal{C}^{\rm c}(t_{n})$
$\displaystyle=\left\langle\mathcal{C}^{\rm
c}_{|Y_{n}}(t_{n})\right\rangle_{Y_{n}}=\sum_{Y_{n}}\,\wp(Y_{n})\,\mathcal{C}^{\rm
c}_{|Y_{n}}(t_{n}).$ (98)

The hope, then, is that by locally maximizing $\mathcal{C}^{\rm c}$ in time,
that the final $\mathcal{C}^{\rm c}(T)$ will be much larger than the
unconditioned coherence, $\mathcal{C}^{\rm nc}$, and perhaps even close to
optimal.

In the following subsection, we will implement the “full-Greedy” algorithm as
presented here, and evaluate Eq. (98) numerically. However, instead of
numerically generating all possible trajectories associated with all possible
measurement results $Y_{N}$, we will perform a stochastic simulation where the
measurement output at each step, _i.e._ , $y_{n}$, is generated randomly
according to its probability. Then, by generating enough stochastic
trajectories, we can evaluate the expected reward, $\mathcal{C}^{\rm
c}(t_{n})$ in Eq. (98) as an average. But in later subsections we will see how
it can be possible to simplify the algorithm so as to perform the exact sum.

Figure 2: Decoherence in data qubit for $\kappa=0.2$ and $K=20$. The thick
brown curve represents no-control decoherence $(1-\mathcal{C}^{\rm nc})$ in
asymptotic regime as in Eq. (36), whereas the green lines represent $100$
trajectories from numerical simulations of the Greedy algorithm. We plotted
the trajectories in light green so that the saturation in the colour shows the
relative probability of a particular decoherence value at a particular time.
Note that Greedy selects a different waiting time at each step, but the
difference between them in this case is of order $10^{-3}/\bar{\gamma}$, which
is too small to see in the figure, giving the appearance that all branches
occur at the same time. Three typical trajectories are bolded with thick
dashed lines and one atypical trajectory, which happened only once, is shown
in dotted red. The thick dark green line is the average of the $100$
trajectories, while the black dashed line represents the exact mean of
Greedy4, introduced in Sec. VI.4.

###  VI.2 Numerical implementation of Greedy

For our numerical calculations, we choose
$\gamma_{\uparrow}=\gamma_{\downarrow}=\bar{\gamma}=1$ for simplicity and the
other dynamical parameters ($\kappa$ and $K$) are in the units of
$\bar{\gamma}$. Equating RTP up and down jump rates yields symmetric dynamics
that are easier to analyze. Also, we assume that, at $t=0$, the RTP is in the
steady state, _i.e._ , $\underline{A}_{0}=\underline{P}_{\rm ss}$, where
$\underline{P}_{\rm ss}$ is defined in Eq. (8). We choose $\kappa=0.2$ and a
range of $K$ values from $2$ to $100$. At the higher end of ranges for $K$,
condition (2) will be satisfied. For the (theoretically infinitesimal) time
step ${\rm d}t$, it is most critical that it is small compared to $1/K$. We
choose ${\rm d}t=0.001/K$.

The Greedy algorithm requires optimization over the three measurement angles
as in Eqs. (96). Fortunately, our locally optimal strategy has the same
structure as that introduced in Ref. [50]. This means we can use analytical
solutions described in [51] to reduce the range of each possible angle from a
real interval $[-\pi/2,\pi/2]$ to only three discrete choices. This is
described in Appendix C.

As we discussed above, after finding an optimal ${\mu}_{n+1}$, we need to
simulate the measurement outcome $y_{n+1}$ in order to calculate
$\underline{{A}}_{n+1}$. For our stochastic simulation we randomly generate
$y_{n+1}$ using a probability function $\wp(y_{n+1}|Y_{n})$, where $Y_{n}$ is
the measurement record assumed known up to that time. This function is encoded
in vector $\underline{{\check{A}}}_{n+1}$, following Eq. (76), where we write

$\wp(y_{n+1}|Y_{n})=\frac{\wp(Y_{n+1})}{\wp(Y_{n})}=\frac{\underline{{I}}^{\top}\underline{{\check{A}}}_{n+1}}{\underline{{I}}^{\top}\underline{{\check{A}}}_{n}}.$
(99)

To be able to use the above equation, we update the vector
$\underline{{\check{A}}}_{n}$ after each measurement using Eq. (78), with the
initial probability vector $\underline{{\check{A}}}_{0}=\underline{{P}}_{\rm
ss}$, same as $\underline{{A}}_{0}$ as in Eq. (80).

We show the conditional coherence $\mathcal{C}^{\rm c}_{|Y_{n}}(t_{n})$ in
Fig. 2 from numerical implementation of the Greedy algorithm using
aforementioned parameters and $K=20$. The light green lines are $100$
trajectories of $\mathcal{C}^{\rm c}_{|Y_{n}}(t_{n})$ where three of them are
exaggerated (bold-dashed lines) as examples. To calculate the expected reward
function $\mathcal{C}^{\rm c}(t_{n})$, we average over these $100$
trajectories. It should be noted that the waiting times, $\tau_{n}$, in
different trajectories do not match, so we could not simply calculate the
average decoherence at each $t_{n}$. Instead, for each trajectory we
interpolated a line between consecutive $(t_{n},\mathcal{C}^{\rm
c}_{|Y_{n}}(t_{n}))$ points, giving rise to the broken lines shown in Fig. 2.
Then we took the average of the broken lines, which is shown in dark green.
Comparing the average with the no-control coherence $\mathcal{C}^{\rm nc}(t)$,
we see that the Greedy algorithm can massively improve the data qubit’s
coherence.

Figure 3: Greedy chosen values of next measurement parameters,
$\\{\theta_{n+1},\tau_{n+1}\\}$, together with the information gained via
$n$th measurement, $\zeta_{n}$, $\alpha_{n}$ and $y_{n}$, in a single
trajectory of the Greedy algorithm for $\kappa=0.2$ and $K=20$. For
$\theta_{n+1}$ and $\zeta_{n}$, their absolute values are shown in the
vertical axes, while their signs are depicted via filled and hollow markers
for positive and negative signs, respectively.

###  VI.3 Greedy choices of measurement settings

We would like to understand how the greedy algorithm chooses the measurement
settings ${\mu}_{n+1}=\\{\theta_{n+1},\tau_{n+1}\\}$, based on the measurement
result $y_{n}$ and relevant statistics $\alpha_{n}$ and $\zeta_{n}$. To do
that, we look at individual trajectories of the Greedy chosen settings. A
typical one is presented in Fig. 3. For $\zeta_{n}$ and $\theta_{n+1}$, their
values can change sign, so we plotted the absolute values of them in order to
see the pattern better, while the sign of them are depicted by filled and
hollow markers for positive and negative respectively.

We observe in Fig. 3 that, after an initial transient, the statistical
parameters, $\alpha_{n}$ and $|\zeta_{n}|$, jump between only two values each,
and these jumps are correlated with the measurement result $y_{n}$. We also
observe that whenever a non-null measurement $(y_{n}=1)$ occurs, the sign of
$\zeta_{n}$ changes, that is depicted by the filled markers switched to hollow
markers. Moreover, from Fig. 3(d), we see that the Greedy choices for the
waiting time $\tau_{n+1}$ and the absolute value of measurement angle
$|\theta_{n+1}|$ are, apart from transients, again confined to two values and
their jumps are correlated with jumps in $\\{\alpha_{n},\zeta_{n}\\}$.
Finally, the sign of $\theta_{n+1}$ is the same as the sign of $\zeta_{n}$.
These observations lead us to proposing a new approach to implement a simpler
version of Greedy algorithm that we explain next.

###  VI.4 “Greedy4” with four parameters

We can capture the Greedy choices of measurement parameters as just explored
by setting

$\displaystyle\theta_{n+1}$ $\displaystyle:=s_{n}\Theta_{\rm
G}(\alpha_{n},|\zeta_{n}|),$ (100a) $\displaystyle\tau_{n+1}$
$\displaystyle:=\frac{\Delta_{\rm G}(\alpha_{n},|\zeta_{n}|)}{K}.$ (100b)

Here the function $\Theta_{\rm G}$ sets the absolute value of the Greedy
choice of measurement angle, while

$s_{n}:={\rm sign}(\zeta_{n}),$ (101)

determines its sign, and the function $\Delta_{\rm G}$ sets the waiting time
$\tau_{n+1}$. This $\Delta_{\rm G}$ is dimensionless and of order unity,
since, as we have discussed, the waiting time scales with $1/K$.

Although, in general, $\Delta_{\rm G}$ and $\Theta_{\rm G}$ may be functions
of both $\alpha_{n}$ and $\zeta_{n}$, we see from Fig. 3 that knowing
$\alpha_{n}$ is sufficient. In fact, since (after some transient steps)
$\alpha_{n}$ jumps between two values, denoted by
$\\{\alpha^{>},\alpha^{<}\\}$, it is sufficient to know the binary variable
$a_{n}\in\\{>,<\\}$, defined as

$a_{n}=\begin{cases}>&\text{if }\alpha_{n}>\bar{\alpha}\ ,\\\ <&\text{if
}\alpha_{n}<\bar{\alpha}\ ,\end{cases}$ (102)

where $\bar{\alpha}:=(\alpha^{>}+\alpha^{<})/2$. Note that, most of the time,
$a_{n}=\ >$. Thus, we simplify $\Theta_{\rm G}$ and $\Delta_{\rm G}$ as binary
functions, each with two values as

$\Theta_{\rm G}(a_{n})=\begin{cases}\Theta_{\rm G}^{>}&\text{if }a_{n}=\ >\
,\\\ \Theta_{\rm G}^{<}&\text{if }a_{n}=\ <\ ,\end{cases}$ (103a) and
$\Delta_{\rm G}(a_{n})=\begin{cases}\Delta_{\rm G}^{>}&\text{if }a_{n}=\ >\
,\\\ \Delta_{\rm G}^{<}&\text{if }a_{n}=\ <.\end{cases}$ (103b)

It should be noted that the above equations could have been defined in terms
of $|\zeta_{n}|$ rather than $\alpha_{n}$. However, because the difference
between $\alpha^{>}$ and $\alpha^{<}$ is more pronounced, it is more natural
to use this as the variable on which $\Theta_{\rm G}$ and $\Delta_{\rm G}$
depend.

We calculate numerical values of $\\{\Delta_{\rm G}^{>},\Delta_{\rm
G}^{<},\Theta_{\rm G}^{>},\Theta_{\rm G}^{<}\\}$ as follows. Each Greedy
trajectory provides us with a sequence of Greedy chosen measurement angles
$\\{\theta_{n}\\}$, a sequence of waiting times $\\{\tau_{n}\\}$, and an
associated sequence of $\\{\alpha_{n}\\}$. First, we discard the transient
data (the first few time steps) and then separate the chosen angles into two
groups by checking if the associated $\alpha_{n}$ is larger or smaller than
the threshold $\bar{\alpha}$. That is, we have two sets of angles defined as

$\displaystyle\vartheta^{>}$ $\displaystyle=\\{|\theta_{n}|:a_{n}=\ >\\},$
(104a) $\displaystyle\vartheta^{<}$ $\displaystyle=\\{|\theta_{n}|:a_{n}=\
<\\}.$ (104b)

Similarly for the waiting times, scaled by $K$, we define two sets as

$\displaystyle\Upsilon^{>}$ $\displaystyle=\\{\tau_{n}/K:a_{n}=\ >\\},$ (105a)
$\displaystyle\Upsilon^{<}$ $\displaystyle=\\{\tau_{n}/K:a_{n}=\ <\\},$ (105b)

The above procedure is repeated for all 100 trajectories. Then, we take all
members of $\vartheta^{>}$s and average them to be the numerical value of
$\Theta_{\rm G}^{>}$. The numerical value of $\Theta_{\rm G}^{<}$,
$\Delta_{\rm G}^{>}$ and $\Delta_{\rm G}^{<}$ are calculated similarly,
following the conditions of $a_{n}$ in Eqs. (103). Fig. 4 shows the values of
$\Delta_{\rm G}^{>},\Delta_{\rm G}^{<},\Theta_{\rm G}^{>},\Theta_{\rm G}^{<}$
for a range of the sensitivity $K$ from $2$ to $100$. All four of these
parameters are within 2.5% of $\pi/2$ across the whole range.

Figure 4: Calculated values of the Greedy4 parameters, from top to bottom,
$\Theta_{\rm G}^{<},\Delta_{\rm G}^{<},\Theta_{\rm G}^{>},\text{ and
}\Delta_{\rm G}^{>}$, for a range of $K$. See text for details of calculation.
The horizontal dashed line shows the value of $\pi/2$; see Eq. (112).

Now that we know the binary values of $\Theta_{\rm G}$ and $\Delta_{\rm G}$,
we no longer need the full optimization part of the Greedy algorithm. As a
result, we obtain a new adaptive algorithm in which an experimenter can simply
use (103) with the knowledge of $\alpha_{n}$ to decide measurement settings
for the next step. We call this new adaptive algorithm “Greedy4” (Greedy with
sub-index $4$), since it is defined by four parameters, namely $\Delta_{\rm
G}^{>},\Delta_{\rm G}^{<},\Theta_{\rm G}^{>}$ and $\Theta_{\rm G}^{<}$. We
emphasize again that in order to determine these four parameters, the full
Greedy algorithm, including the optimization part, must be run. Once the four
parameters are found, the experimenter could use them for future experiments.

To calculate the expected reward $\mathcal{C}^{\rm c}(t_{n})$ as in Eq. (98)
for Greedy4, we numerically generate all possible $Y_{N}$ and average the
associated reward functions with the probability weight $\wp(Y_{N})$ from Eq.
(76). Each realisation of $Y_{N}$, like full-Greedy, is associated with a
different sequence of waiting times $\tau_{n}$, which take two values.
However, as shown in Fig. 4, the difference between $\Delta_{\rm G}^{>}$ and
$\Delta_{\rm G}^{<}$ is not great. Moreover, as we have already seen in Fig.
3, most of the time the value $\Delta_{\rm G}^{>}$ is chosen. The other value
is chosen only when the system detects that a jump has occurred, which happens
at approximately the same rate as the jumps, $\breve{\gamma}$. Thus it is not
a bad approximation to take the average as if every $\tau_{n}$ were the same,
an effective waiting time is given by

$\tau_{\rm eff}=\frac{\Delta^{>}}{K}(1-\breve{\gamma}\,\tau_{\rm
eff})+\frac{\Delta^{<}}{K}\breve{\gamma}\,\tau_{\rm eff}.$ (106)

Here $\tau_{\rm eff}$ appears on the right hand side because
$\breve{\gamma}\tau_{\rm eff}$ is the probability, in one waiting time, that
the RTP jumps. Solving for $\tau_{\rm eff}$ we obtain

$\tau_{\rm eff}=\frac{\Delta^{>}}{K-\breve{\gamma}(\Delta^{<}-\Delta^{>})}.$
(107)

Using this $\tau_{\rm eff}$, we find the expected reward function for Greedy4,
which is shown with a black dashed line in Fig. 2. This line is an _exact_
average in the sense that it is not stochastic. Moreover, this line agrees
with the stochastic average using the full Greedy, within the stochastic
error, and confirms that the expected reward function calculated from the full
Greedy is in agreement with Greedy4.

Figure 5:  Dynamics of the Greedy4 algorithm in phase space (top row) and
probability of $\alpha$ (bottom row). We used $\kappa=0.2$ in all plots, while
$K=2,20,\text{and}~{}100$, from left to right. The value of $\bar{\alpha}$ is
indicated by a dashed line in each plot. The blue arrows represent maps due to
null results ($y=0$), while the red arrows represent non-null outcomes
($y=1$). Less frequent transitions are shown by red dashed arrows. We have
only shown half of the arrows in the instance of $K=2$; the reader can
visualize the other half by mirroring the image around the $\zeta=0$ axis.

###  VI.5 Dynamics of Greedy4 in $(\alpha,\zeta)$ phase space

In this subsection, we use the sufficient statistics $(\alpha,\zeta)$ as a way
to visualize the dynamics of the Greedy4 algorithm in “phase space”. The plots
in Fig. 5 show phase space dynamics of 100 trajectories for $K=2,20\text{ and
}100$. We use light green dots to represent all $(\alpha,\zeta)$ states and
thus the green color saturation indicates the relative probability of being in
that state. Blue arrows in the figure show the state transition due to the
null result, _i.e._ , $y=0$, while the red arrows show that of the non-null
result, _i.e._ , $y=1$. The maps used for the transitions are shown as labels.
Since Greedy4 chooses between the binary values of the measurement angle and
the waiting time, we can simplify the notation for the $\mathbf{F}$ map
introduced in Sec. V.2 as

$\mathbf{F}_{s,a}^{y}:=\mathbf{F}\left(s\Theta_{\rm G}({a}),\Delta_{\rm
G}({a})/K,y\right),$ (108)

where $\Theta_{\rm G}({a})$ and $\Delta_{\rm G}({a})$ are the functions
defined in Eqs. (103). Note that, as a consequence of taking
$\gamma_{\uparrow}=\gamma_{\downarrow}$, the phase diagrams for
$(\alpha,\zeta)$ are symmetric in reflection around the $\zeta=0$ line.
Furthermore, the bottom panel of each sub-figure depicts the normalized
histograms for different values of $\alpha$, which can be interpreted as the
probability of the system being on that specific $\alpha$.

In these phase diagrams, we see that the statistical states $(\alpha,\zeta)$
have a tendency to be in only a few states, instead of spreading over the
entire phase space. This was already noted in the discussion of Fig. 3 for
$K=20$, but here we visualize the correlated changes in $\alpha,\zeta$ more
easily, as well as the transition from $K$ of order $1$ to $K\gg 1$. In all
cases, the two most occupied states turn out to correspond to the eigenstates
with the largest absolute eigenvalues of $\mathbf{F}_{+,>}^{0}$ and
$\mathbf{F}_{-,>}^{0}$ (Recall that these are $2\times 2$ complex linear maps
acting on $\underline{{A}}$, and it turns out that they each have two
eigenstates). We denote these stable eigenstates (the ones with largest
absolute eigenvalues) by $\underline{{E}}_{\,+,>}^{0}$ and
$\underline{{E}}_{\,-,>}^{0}$, respectively. The $(\alpha,\zeta)$-values for
these eigenstates are depicted as the dashed circles in the phase space plot.

To understand the state dynamics, let us first consider a state being in
$\underline{{E}}_{\,+,>}^{0}$, the stable eigenstate of
$\mathbf{F}_{\,+,>}^{0}$. In this case, any null-result $(y=0)$ measurement
will map the state to itself and no effect on the system’s phase space
configuration. On the other hand, a non-null measurement ($y=1$) will map the
state to a point on the left side, with $\alpha<\bar{\alpha}$ while switching
the sign of $\zeta$ from positive to negative. Any subsequent non-null result,
which is less common, will cause the sign of $\zeta$ to switch again, while
$\alpha$ remains in the same region $\alpha<\bar{\alpha}$. However, since the
probability of two consecutive non-null results is low, the subsequent
measurement most likely yields a null result, mapping the system state to a
point on the right side with $\alpha>\bar{\alpha}$, while keeping the sign of
$\zeta$ unchanged. The following null results, as they are more likely, will
move the state towards the other stable eigenstate,
$\underline{{E}}_{\,-,>}^{0}$.

The above description holds for any value of $K>\bar{\gamma}$. But there are
differences as $K$ changes. First consider the case for a relatively low
sensitivity, $K=2$, _i.e._ , of the same order as $\bar{\gamma}=1$. In this
case we observe in Fig. 5(a) that the dynamics in the phase space can be
grouped into eight locations. We understand the dynamics as follows. For small
$K$, we find that the likelihood of consecutive non-null measurement results
is higher compared to that of the larger $K$. So, if the system configuration
has moved to the left side of the phase space, as a consequence of a
measurement result, any further non-null measurements yield maps that changes
the sign of $\zeta$ but keeps the $\alpha$ on the left side. These maps as we
showed them with red dashed arrows, reveals four main locations on the left
side of the phase space. Furthermore, we observe that a map due to a null
result, $(y=0)$, brings the system from the left to the right side of the
phase space, but the resulting configuration is distanced from the
eigenstates. Then, further null results bring the system to the eigenstate. As
a result, there are four main locations on the right side. Note that we have
not shown all the transitions in the $K=2$ case to avoid messiness (see
caption).

As $K$ increases, the first thing that happens is that the right-most points
on the phase diagram all converge to the eigenstates
$\underline{{E}}^{0}_{\pm,>}$. That is, for $\alpha>\bar{\alpha}$ we have only
these two eigenstates to consider. This is seen already in Fig. 5(b). As $K$
increases even further, reaching the asymptotic regime $K\gg\bar{\gamma}$, the
chance of two consecutive non-null results is negligible. This leads to there
also being only two points on the left side, as seen in Fig. 5(c). Thus, the
dynamics are simplified and can be grasped by concentrating on just four
points: $\underline{{E}}_{\,\pm,>}^{0}$, and
$\mathbf{F}_{\pm,>}^{1}\underline{{E}}_{\,\pm,>}^{0}$. Moreover, the system
spends the majority of its time in one of the eigenstates and occasionally,
with a non-null ($y=1$) measurement result, goes to the $a=\,<$ state while
altering $s={\rm sign}(\zeta)$. The next measurement almost always yields
$y=0$, causing a jump back to an eigenstate while keeping $s$ the same.

###  VI.6 “Greedy2” in the asymptotic regime

In the asymptotic regime $K\gg\bar{\gamma}$, we show that the Greedy4
algorithm in Eqs.(103) can be simplified even further. By inspecting Fig. 4
carefully, we observe that in this regime the four parameters converge to only
two, because we have the identities $\Theta_{\rm G}^{>}=\Delta_{\rm G}^{>}$
and $\Theta_{\rm G}^{<}=\Delta_{\rm G}^{<}$. From the original definitions in
Eq. (100), it means that the Greedy-chosen measurement angles and waiting
times are related via a simple relation:

$\theta_{n+1}=s_{n}K\tau_{n+1}.$ (109)

As a result, we have a simpler algorithm with just two parameters:
$\Theta_{\rm G}^{>}$ and $\Theta_{\rm G}^{<}$. We refer to this algorithm as
“Greedy2”, which is valid only for large $K$ or the asymptotic regime. The
algorithm can be summarized as

$\displaystyle\theta_{n+1}=$ $\displaystyle\,\,s_{n}\Theta_{\rm G}(a_{n}),$
(110a) $\displaystyle\tau_{n+1}=$ $\displaystyle\,\,\Theta_{\rm G}(a_{n})/K,$
(110b)

where $\Theta_{\rm G}(a_{n})$ is given by Eq. (103a).

This choice has an intuitive explanation, as follows. The term $K\tau_{n+1}$
is simply a SQ’s phase evolved (after the SQ’s phase reset to zero) during the
time interval $(t_{n},t_{n+1})$, where $\tau_{n+1}=t_{n+1}-t_{n}$, given that
the RTP stays unchanged with a value of $s_{n}$ during that time. That is,
under the approximation that $z(t)$ stays equal to our best estimate of its
value, $s_{n}$ at the start of the interval, the evolved SQ’s phase is given
by

$\Phi(t_{n+1})=K\int_{t_{n}}^{t_{n+1}}\\!\\!z(t){\rm
d}t\,\approx\,{s}_{n}K\tau_{n+1}.$ (111)

This approximation can be justified in the asymptotic regime, as our
confidence in the value of $z_{n}$ is high ($|\zeta_{n}|\approx 1$) and the
probability of a jump in the interval is low, $O(\breve{\gamma}/K)$.
Therefore, Eq. (109) simply says that the measurement angle should be chosen
to be exactly the most likely SQ phase at the time of measurement, _i.e._ ,
$\theta_{n+1}=\Phi(t_{n+1})$. In other words, the measurement probes whether
the SQ is in the most likely state. A null result ($y_{n+1}=0$) is the answer
‘yes’ while the other result ($y_{n+1}=1$) is the answer ‘no’. That is, in
this aspect, Greedy2 performs an optimal choice for minimum-error binary state
discrimination [52, 53], where one state (which is pure) arises if and only if
$z(t)=s_{n}$ for the entire interval, and the other state (which is mixed)
arises if and only if $z(t)=-s_{n}$ for some of the interval.

Note that the above simple intuition does not explain the value chosen by
Greedy2 for its two parameters, $\Theta_{\rm G}^{>}$ and $\Theta_{\rm G}^{<}$.
However, inspecting Fig. 4 again, we see that one of these values does have a
simple interpretation in the asymptotic limit:

$\displaystyle\lim_{K/\bar{\gamma}\to\infty}\Theta^{>}_{\rm G}=$
$\displaystyle\,\,\pi/2.$ (112)

Note that this is the more important of the two parameters, as in the
asymptotic regime, the measurement result is almost always a null result, the
system state is almost always in the region $\alpha_{n}>\bar{\alpha}$, and the
measurement choice is almost always $\Theta_{\rm G}(a_{n})=\Theta_{\rm
G}^{>}$. The choice of $\pi/2$ for $s_{n}\theta_{n+1}$ and $K\tau_{n+1}$
implies that the most likely pure state (from $z(t)=s_{n}$ for the entire
interval) and the second-most likely pure state (from $z(t)=-s_{n}$ for the
entire interval) are orthogonal at the measurement time, and can be perfectly
distinguished by the measurement.

##  VII Map-based Optimized Adaptive Algorithm for Asymptotic Regime (MOAAAR)

While using the Greedy algorithm for the SQ offers a great reduction in the
decoherence of the data qubit, it is only locally optimal. There is no reason
to expect it to be globally optimal. In this section we introduce and
analytically study a new algorithm that is plausibly optimal in the asymptotic
regime, $K\gg\gamma_{\uparrow,\downarrow}\gg\kappa$. However, to design the
algorithm we use the knowledge and intuitions gained from the numerical
investigation of Greedy in this regime.

First, recall from Sec. VI.6 that in asymptotic regime, Greedy chooses
measurement settings that are perfectly matched according to Eqs. (110).
Second, recall from Sec. VI.5 that out of the two values of $\Theta_{\rm
G}(a_{n})$, the choice is almost always $\Theta_{\rm G}(>)$, which we called
$\Theta_{\rm G}^{>}$. Considering just the single parameter, $\Theta_{\rm
G}^{>}$, one could use it to define what one might call a Greedylike1
algorithm, in which $\Theta_{\rm G}(a)$ in Greedy2 is replaced by the constant
$\Theta_{\rm G}^{>}$. This would reproduce the performance of the Greedy
algorithm in the asymptotic regime, although it would not actually approximate
the Greedy algorithm itself, which always uses two distinctly different
angles, even in the asymptotic regime, as per Greedy2.

Like the Greedy algorithm, this Greedylike1 algorithm would still be adaptive,
in general, because of the $s_{n}$ in Eq. (110a). However, we already know the
value of $\Theta_{\rm G}^{>}$ in the asymptotic regime, Eq. (112). While this
choice is intuitive, it in fact corresponds to a non-adaptive measurement, as
$\theta=s\pi/2$ is measuring in the same basis for $s=\pm 1$. It would be
surprising if the globally optimal performance in the asymptotic regime could
be attained by a non-adaptive algorithm, and indeed we find that this is not
the case. We do this by cleaving to the well-motivated relations between the
measurement angles and the waiting times in Eqs. (110), and to the
understanding that the asymptotic performance can be expected to be reproduced
with $\Theta$ fixed, but without cleaving to the particular choice of $\Theta$
made by the Greedy algorithm. That is, as expected, we can do better by
choosing $\Theta$ in a globally optimal manner rather than a locally optimal
(greedy) manner.

To be precise, we here propose an adaptive (in general) algorithm for
measurement angles and waiting times,

$\displaystyle\theta_{n+1}=$ $\displaystyle\,s_{n}\Theta,$ (113a)
$\displaystyle\tau_{n+1}=$ $\displaystyle\,\Theta/K,$ (113b)

as in Eqs. (110), but with the function $\Theta_{\rm G}$ replaced by a fixed
parameter $\Theta$. This simplifies the calculation of the expected reward
function in the asymptotic regime, as we need consider only the following set
of four single-parameter maps:

$\mathbf{F}_{s}^{y}(\Theta):=\mathbf{F}(s\Theta,\Theta/K,y),$ (114)

where $s\in\\{+1,-1\\}$ and $y\in\\{0,1\\}$. Similar to the preceding section,
we denote the stable eigenstates of these maps by $\underline{{E}}_{\,s}^{y}$.

In this section, we first show how the problem of maximizing the expected
reward function is, in the regime of interest (2), equivalent to minimizing an
expected decoherence rate. We then apply the adaptive strategy, Eqs. (113),
and derive an analytical expression for the decoherence rate as a function of
$\Theta$ in the asymptotic regime. This can be used to search for a globally
optimal parameter value, $\Theta^{\star}$, which defines our Map-based
Optimized Adaptive Algorithm for Asymptotic Regime (MOAAAR). We call it Map-
based because of the central role of the Bayesian maps, as mentioned in the
preceding paragraph. We emphasize that it is adaptive because we will find
$\Theta^{\star}\neq\pi/2$, which means that, via $s_{n}$, previous measurement
results do alter the basis in which the SQ is measured. (Recall that this is
in contrast with the $\Theta_{\rm G}^{>}=\pi/2$ choice of the Greedylike1
algorithm in this regime.)

###  VII.1 Expected reward and decoherence rate

Recall the original expected reward function Eq. (55) defined for the final
time $T$, when the phase-correction control is to be applied. This is a sum
over exponentially many terms, which makes it infeasible to do an analytic
maximization. However, from the numerical results in the previous section
(Greedy), we know that the system’s state $(\alpha,\zeta)$ effectively jumps
between only a finite number of states, and explores all of them over time. As
we will see, this is also true for the globally optimized MOAAAR strategy.
This means that the stochastic evolution of our state of knowledge is ergodic
[37], so the average decay of the coherence over a typical trajectory over
many steps ($T\gg K^{-1}$) is the same as an ensemble average over one step.
There will be deviations from this due to initial [$t=O(1/\breve{\gamma})$]
and final [$T-t=O(1/\breve{\gamma})$] transients, but these can be ignored for
long times as we are interested in. Thus we can then safely use an expected
reward function at an intermediate time $t_{n}$, $\mathcal{C}^{\rm c}(t_{n})$
as a proxy of the expected reward at the actual final time, $\mathcal{C}^{\rm
c}(T)$.

We also saw in the preceding section that in the regime $K\gg\bar{\gamma}$,
the Greedy algorithm makes the system spend almost all of the time close to
the two stable eigenstates $\underline{{E}}_{\,\pm}^{y=0}$. We find this to be
even more true with MOAAAR; we will quantify deviations from this assumption
in Sec. VII.3. Thus we can approximate the coherence (the reward function) at
time $t_{n}$ as

$\displaystyle\mathcal{C}^{\rm c}(t_{n})=$
$\displaystyle\,\sum_{Y_{n}}\wp(Y_{n})|\mathcal{A}_{|Y_{n}}|=\sum_{Y_{n}}\,\big{|}\underline{{I}}^{\top}\,\underline{{A}}_{n}\big{|}$
$\displaystyle=$ $\displaystyle\sum_{\begin{subarray}{c}Y_{n}|\\\
s_{n}=+1\end{subarray}}\\!|\beta_{+}(Y_{n})\,\underline{{I}}^{\top}\,\underline{{E}}_{+}^{0}|+\\!\\!\sum_{\begin{subarray}{c}Y_{n}|\\\
s_{n}=-1\end{subarray}}\\!|\beta_{-}(Y_{n})\,\underline{{I}}^{\top}\,\underline{{E}}_{-}^{0}|.$
(115)

The first line follows from the definitions of the reward function in Eq. (98)
or Eq. (55). In the second line, we have separated the summation,
$\sum_{Y_{n}}\big{|}\underline{{I}}^{\top}\,\underline{{A}}_{n}\big{|}$, into
two contributions associated with $s_{n}=+1$ and $s_{n}=-1$, and used the fact
that the state $\underline{{A}}_{n}$ is almost always proportional to one of
the two eigenstates ($\underline{{E}}_{\,+}^{0}$ and
$\underline{{E}}_{\,-}^{0}$). Which of the two states pertains at time $t_{n}$
is determined solely by $s_{n}$ (which is a function of $Y_{n}$), while
$\beta_{+}(Y_{n})$ and $\beta_{-}(Y_{n})$ are the appropriate proportionality
scalars.

Let us now consider the expected reward at the next SQ’s measurement at
$t_{n+1}=t_{n}+\tau$, which can be computed by applying one additional map,
${\mathbf{F}}_{\pm}^{y}(\Theta)$, to Eq. (115) and summing over its
measurement results, $y_{n+1}\in\\{0,1\\}$, for which we use the dummy
variable $y$. This gives

$\displaystyle\mathcal{C}^{\rm c}(t_{n+1})=$
$\displaystyle\sum_{y}\sum_{\begin{subarray}{c}Y_{n}|\\\
s_{n}=+1\end{subarray}}\\!|\beta_{+}(Y_{n})\,\underline{{I}}^{\top}\,{\mathbf{F}}_{+}^{y}(\Theta)\,\underline{{E}}_{+}^{0}|$
$\displaystyle+\sum_{y}\sum_{\begin{subarray}{c}Y_{n}|\\\
s_{n}=-1\end{subarray}}\\!|\beta_{-}(Y_{n})\,\underline{{I}}^{\top}\,{\mathbf{F}}_{-}^{y}(\Theta)\,\underline{{E}}_{-}^{0}|.$
(116)

Note that the signs $s_{n}$ determine which $\mathbf{F}$-map is applied, and
matches the subscript of the eigenstates.

Next, we rearrange the formula Eq. (116) such that it is written in terms of
the coherence at time $t_{n}$ in Eq. (115), with one (justified)
approximation. We can interpret the two terms in Eq. (115) as two weighted
conditional coherence contributions, namely $P_{\rm
ss}(s_{n}=+1)\,\mathcal{C}^{\rm c}_{|s_{n}=+1}(t_{n})$ and $P_{\rm
ss}(s_{n}=-1)\,\mathcal{C}^{\rm c}_{|s_{n}=-1}(t_{n})$, respectively. [Recall
that we defined a general conditioned coherence in Eq. (57).] Since $s_{n}$ is
our best guess for the RTP $z_{n}$, and since we are not in a transient
regime, we can take these probabilities to be the same as the steady-state
probabilities for the RTP. Also, in this long time steady-state regime, we can
safely say that the value that $s_{n}$ happens to have at this particular time
has almost no bearing on the coherence at this time, which is determined by a
long history of stochastic events, only the most recent of which is relevant
to $s_{n}$. That is, we can make the approximation $\mathcal{C}^{\rm
c}_{|s_{n}=+1}(t_{n})\approx\mathcal{C}^{\rm
c}_{|s_{n}=-1}(t_{n})\approx\mathcal{C}^{\rm c}(t_{n})$. Rewriting the terms
involving $\beta_{\pm}(Y_{n})$ in Eq. (115) in terms of conditional
coherences, we can simplify the expected reward at time $t_{n+1}$ to

$\displaystyle\mathcal{C}^{\rm c}(t_{n+1})=$
$\displaystyle\,\sum_{y}\Bigg{\\{}\frac{P_{\rm
ss}(s_{n}=+1)|\underline{{I}}^{\top}\,{\mathbf{F}}_{+}^{y}(\Theta)\underline{{E}}_{+}^{0}|}{|\underline{{I}}^{\top}\,\underline{{E}}_{+}^{0}|}\mathcal{C}^{\rm
c}_{|s_{n}=+1}(t_{n})$ $\displaystyle\quad+\frac{P_{\rm
ss}(s_{n}=-1)|\underline{{I}}^{\top}\,{\mathbf{F}}_{-}^{y}(\Theta)\underline{{E}}_{-}^{0}|}{|\underline{{I}}^{\top}\,\underline{{E}}_{-}^{0}|}\mathcal{C}^{\rm
c}_{|s_{n}=-1}(t_{n})\Bigg{\\}}$ $\displaystyle\approx$
$\displaystyle\,[1-{\bar{\Gamma}(\Theta)}\tau]\,\mathcal{C}^{\rm c}(t_{n}).$
(117)

This describes an exponential decay with an expected decoherence rate defined
as

$\bar{\Gamma}(\Theta):=\sum_{s=\pm}{P}_{\rm
ss}(s)\frac{|\underline{I}^{\top}\underline{E}_{s}^{0}|-\sum_{y}|\underline{I}^{\top}\,{\bf
F}_{s}^{y}(\Theta)\,\underline{E}_{s}^{0}|}{\tau|\underline{I}^{\top}\underline{E}_{s}^{0}|}.$
(118)

Here we have dropped the $n$ subscript entirely as it is no longer necessary.

The $\bar{\Gamma}(\Theta)\tau$ in Eq. (VII.1) can also be thought of as a
steady-state average of the relative change of the conditional coherence, as
the system’s state evolves from $\underline{{A}}_{n}$ (conditioned on $Y_{n}$)
at time $t_{n}$ to the next time $t_{n+1}=t_{n}+\tau$ with an unknown
$y_{n+1}$. This relative change was called $\delta(\underline{{A}}_{n})$ in CL
[28] (where the vector $\underline{{A}}_{n}$ is its argument, not the scalar
quantity which is changing). In terms of this, we have a conditional rate,

$\displaystyle\Gamma(\Theta,\underline{{A}}_{n})$
$\displaystyle=\frac{\delta(\underline{{A}}_{n})}{\tau},$
$\displaystyle=\frac{1}{\tau}\frac{\left|\mathcal{A}_{|Y_{n}}\right|-\sum_{y_{n+1}}\wp({y_{n+1}}|Y_{n})\left|\mathcal{A}_{|Y_{n+1}}\right|}{\left|\mathcal{A}_{|Y_{n}}\right|},$
(119)

and the average rate in Eq. (118) is obtained by averaging over $Y_{n}$. This
average exponential-decay rate well approximates the coherence decay rate for
a typical trajectory, as argued above. Therefore we can use
$\bar{\Gamma}(\Theta)$ as a proxy for the global reward function (55). That
is, we reformulate the problem of maximizing the expected final coherence as
the problem of minimizing the average decoherence rate. This can be done
essentially analytically, as we see in the following subsection.

Figure 6: The $(\alpha,\zeta)$ phase-space plots for MOAAAR using optimal
measurement angle $\Theta=\Theta^{\star}$ for different values of $K$. Similar
to the Greedy case, we chose $\bar{\gamma}=1$ and $\kappa=0.2$. It is clear
that when $K$ is large ($K=100$), the system’s states are almost always in
only two stable states (corresponding to the stable eigenstates of the two
maps $\mathbf{F}_{\pm}^{0}(\Theta^{\star})$.

###  VII.2 Analytical expression for decoherence

We wish to find the minimum of the function $\bar{\Gamma}(\Theta)$ in Eq.
(118). Recall that $\Theta$ is the single parameter that defines our adaptive
algorithm as per Eqs. (113). This in term defines the four maps in Eq. (114)
(from two choices of $s\in\\{+1,-1\\}$ and two choices of $y\in\\{0,1\\}$), as
a function of $\Theta$, via the original definition in Eq. (66). We are
interested in the stable eigenstates $\underline{{E}}_{\pm}^{0}$ of the map
${\bf F}^{0}_{\pm}(\Theta)$. With the help of Mathematica, we analytically
solve for these eigenvectors and their expansions in the asymptotic regime,
taking $\gamma_{\uparrow},\gamma_{\downarrow}\ll K$. Ignoring an irrelevant
complex scalar multiple, these eigenvectors can be represented by the two
sufficient statistics introduced in Sec. V.4. These are the parameters
$\alpha_{s}^{0}$ and $\zeta_{s}^{0}$ found by substituting
$\underline{{A}}_{n}=\underline{{E}}_{s}^{0}$ into
$\alpha_{n}(\underline{{A}}_{n})$ (83a) and $\zeta_{n}(\underline{{A}}_{n})$
(85) respectively. Keeping terms up to $O(1/K)$, we get

$\displaystyle\alpha_{s}^{0}=$ $\displaystyle\,\,\frac{1}{12(2\Theta+\sin
2\Theta)^{2}}(N_{\Theta}+s\frac{\gamma_{\downarrow}-\gamma_{\uparrow}}{K}M_{\Theta}),$
(120a) $\displaystyle\zeta_{s}^{0}=$
$\displaystyle\,\,s\times\bigg{\\{}1-\frac{\gamma_{s}}{K}\big{[}\csc\Theta(\cos\Theta+\Theta\csc\Theta)\big{]}\bigg{\\}},$
(120b)

where we have used $\gamma_{s}$ with the meaning
$\gamma_{+}=\gamma_{\downarrow}$, $\gamma_{-}=\gamma_{\uparrow}$. The terms
$N_{\Theta}$ and $M_{\Theta}$ are lengthy functions of $\Theta$ alone, and are
shown in Appendix. D.

By substituting $\alpha_{s}^{0}$ and $\zeta_{s}^{0}$ in Eqs. (120) for the
eigenstates $\underline{E}_{s}^{0}$ in Eq. (118) and expanding terms to the
lowest order of $(\kappa/K)$, we obtain, to lowest order,

$\bar{\Gamma}(\Theta)\approx H_{\Theta}\breve{\gamma}\kappa^{2}/(2K^{2}),$
(121)

where

$\displaystyle
H_{\Theta}=3\Theta^{2}\csc^{4}\Theta-[2\Theta(\Theta-\cot\Theta)+1]{\rm
csc}^{2}\Theta+\tfrac{1}{3}\Theta^{2}-1.$ (122)

This is our most important result, also shown in the CL [28], where we plot
$H_{\Theta}$. We find that this function has an absolute minimum of
$H^{\star}\approx 1.254$ at $\Theta=\Theta^{\star}\approx 1.50055$. This
measurement angle $\Theta^{\star}$ defines our MOAAAR, and is plausibly the
lowest possible decoherence rate using projective measurements on a SQ within
the regime of Eq. (2). We show in Fig. 6 the phase space plots for
$K=2,20,100$ using the MOAAAR algorithm. This confirms that, even more so than
for the Greedy algorithm in Fig. 5, as $K/\bar{\gamma}$ grows, the system
spends more and more time in the two stable eigenstates.

###  VII.3 Approximation errors and approaching rates

With the analytical expressions for MOAAAR, we can analyze any errors that
could occur from the approximations we have made. First, let us revisit the
approximation, $z_{n}\approx s_{n}={\rm sign}(\zeta_{n})$, used prior to Eq.
(VII.1) to justify using the steady state probability of $z_{n}$ in place of
those for $s_{n}$. From the analytical expression in Eqs. (120), we can see
that $\zeta_{s}^{0}=s[1-O(\gamma_{s}/K)]$, for any value of $\Theta$. Using
the result (correct to even higher order) for $\zeta_{n}$ in Eq. (87), we can
thus say that,

$\displaystyle\zeta_{n}\approx$
$\displaystyle\,\,\frac{\wp(Y_{n},z_{n}=+1)-\wp(Y_{n},z_{n}=-1)}{\wp(Y_{n})}$
$\displaystyle=$ $\displaystyle\,\,\wp(z_{n}=+1|Y_{n})-\wp(z_{n}=-1|Y_{n}).$
(123)

Combining the above, we have that, in the stable eigenstates,

$\wp(z_{n}\neq s_{n})=O(\gamma_{s}/K).$ (124)

Therefore, not only do the probabilities of $s_{n}$ and $z_{n}$ match on
average, but the variables themselves are almost always identical when the
system is in the stable eigenstates. However, this raises the question of the
probability that the system is in one of these eigenstates, to which we now
turn.

Figure 7: Plots showing the exponential approach to the two stable eigenstates
on the right side of the phase space. Subplots (a) and (b) are phase-space
$(\alpha,\zeta)$ plots under the adaptive algorithm using $\Theta=1.0$ (blue
dots) and $\Theta=\Theta^{\star}$ (MOAAAR, maroon dots) respectively. To
generate these points, we use $\mathbf{F}_{s}^{1}E^{0}_{s}$ for both
$s\in\\{+1,-1\\}$ as a starting point and let the state evolve for the number
of time steps $n=9$. In each case, the eigenstates are shown as dashed
circles. (c) A Log-scale plot of distances from the eienstate’s values,
$|\alpha_{n}-\alpha_{s}^{0}|$, using the same data in (a) and (b). The last
maroon data point (at $n=9$) has reached the numerical machine precision in
Mathematica. Here $K=20$, and $\kappa=0.2$ as usual.

Because the phase space variables, $\alpha$ and $\zeta$, are continuous, we
cannot strictly talk about the probability of the system being in a stable
eigenstate. However, we can calculate the probability that it is arbitrarily
close to such an eigenstate. For a given $s$, if null results ($y=0$) keep
occurring, the system’s state moves towards the eigenstate
$\underline{{E}}_{\,s}^{0}$ exponentially fast. This is shown in Fig. 7 (blue
dots), which takes an arbitrary value $\Theta=1.0$ as an example, as well as
the optimal value $\Theta^{\star}$. (We do this because we want the results in
this subsection, used to justify the arguments of the preceding subsection, to
hold for any $\Theta$, not just the optimal value $\Theta^{\star}$ found via
those arguments.) The sign of $\zeta$ stays fixed, while the value of $\alpha$
moves towards $\alpha_{s}^{0}$ exponentially fast for both $s=\pm 1$.
Moreover, we see that the transient values of $|\zeta|$ are even closer to $1$
than the values in the stable eigenstates. Thus, it is very safe to assume Eq.
(124) for the whole evolution.

The approximation that the system state is almost always in a stable
eigenstate $\underline{{E}}_{\pm}^{0}$ was also used in the first step, Eq.
(115), towards finding the decoherence rate in Eq. (118). Thinking about the
alternate expression for the decoherence rate, Eq. (VII.1), we have to ask how
much difference arises, when averaging over all records, as a result of the
system state’s deviation from the stable eigenstates. After a non-null result,
$\alpha$ jumps a finite distance from $\alpha_{s}^{0}$. After $n$ time steps,
the distance (in $\alpha$) between the system’s state and the fixed point is
an exponentially decreasing function, _i.e._ ,
$|\alpha_{n}-\alpha_{s}^{0}|\sim e^{-kn}$, where $k$ is some dimensionless
approaching constant, the slope of the curves in Fig. 7(c). We can thus
approximate the number of steps after a jump that the system is still more
than $\epsilon$-distanced away from the eigenstates as
$n_{\epsilon}=O[\log(1/\epsilon)/k]$. Now, the proportion of its evolution
that the system spends this far from the eigenstates is $n_{\epsilon}$ times
the probability per time step of a jump away from the eigenstates. But the
latter is just equal to the probability of a non-null result. This is the same
order as Eq. (124), the probability that the best estimate of $z_{n}$ is
wrong, which is $O(\breve{\gamma}/K)$. Since a jump takes a system a finite
distance away in the $\alpha$-direction, we make the pessimistic assumption
that if the system is more than $\epsilon$-distanced away from the eigenstates
then the relative error in the quantity of interest (the decoherence rate)
could be of order unity. If, on the other hand, the system is less than
$\epsilon$-distanced away from the eigenstates (which is almost all the time),
then if we are again pessimistic, we will say that the relative error in the
quantity of interest is $O(\epsilon)$. Adding these two contributions
together, we get the relative error to be
$O[\log(1/\epsilon)\,\breve{\gamma}/Kk]+O(\epsilon)$. The minimum value of
this relative error is achieved when we choose $\epsilon$ to scale as
$\breve{\gamma}/Kk\ll 1$, which gives a relative error that is small in the
asymptotic regime. Thus the errors in the approximations leading to Eq. (121)
are insignificant in the regime (2).

###  VII.4 Numerical comparison of performance of Greedy and MOAAAR and other
algorithms

Having completed the analytical analyses of MOAAAR, we turn to numerics to
compare it with the Greedy algorithm. We also compare it to some other, non-
optimized, algorithms. In Fig. 8, we plot the data qubit’s decoherence as a
function of time. Although the time is not long compared to $\bar{\gamma}=1$,
the decoherence for the optimized algorithms are very close to linear in time,
which is a consequence of the very simple phase-space dynamics of those
algorithms. All of the numerical results were calculated _exactly_ using the
coherence (or the expected reward) definition in Eq. (98), where all possible
trajectories were generated using the ${\bf F}$ maps that correspond to the
algorithm-chosen measurement angles and waiting times, for all possible
realizations of $Y_{n}$. For the adaptive schemes we took the first
measurement angle to be $\pi/2$, since in the absence of any initial
information about $z_{1}$, there is no reason to choose one sign for $\theta$
rather than the other. This is the automatic choice of the full Greedy
algorithm, but not of Greedy4. We also do this for later plots.

From worst-to-best, the curves in Fig. 8 are as follows: (i) the no-control
case, Eq. (36). (ii) a non-adaptive algorithm with $\theta_{n+1}=\pi/2$ and
$\tau_{n+1}=1/K$. Even this completely unoptimized algorithm already does much
better. (iii) the adaptive algorithm of Eqs. (113) with the choice $\Theta=1$,
giving $\theta_{n+1}=s_{n}$ and $\tau_{n+1}=1/K$, this being the same waiting
times as in (ii). Here we see that making the measurement angle adaptive
significantly improves the performance. (iv) the Greedy4 algorithm of Eqs.
(103), which does much better again. (v) MOAAAR, which like (iii) uses Eqs.
(113), but this time with the optimized $\Theta^{\star}=1.50055$. This is
better than Greedy, but only by a small amount.

Figure 8: Plots of data qubit’s decoherence using different SQ’s measurement
and control strategies. Analytical results are shown with solid curves, while
numerical results are shown with points connected via dashed lines. From top
to bottom: the no-control case (solid brown), a non-adaptive algorithm with
$\theta_{n+1}=\pi/2$ and $\tau_{n+1}=1/K$ (red squares), an adaptive algorithm
with $\theta_{n+1}=s_{n}$ and $\tau_{n+1}=1/K$ (blue squares), the analytical
MOAAAR from Eq. (121) (solid black), the Greedy4 algorithm (green dots), and
the MOAAAR with $\Theta^{\star}=1.50055$ (maroon square). In all calculation,
we used $\gamma_{\uparrow}=\gamma_{\downarrow}=1$, $\kappa=0.2$, and $K=20$.

Finally in Fig. 8, we also plot the analytical result for MOAAAR from Eq.
(121). This agrees quite well with the numerics, but not perfectly. Indeed,
the predicted decoherence is worse than what is achieved from the exact
simulations, and in fact is worse even than the decoherence achieved by
Greedy. This may surprise the reader, since we claimed in the preceding
section that all the approximations leading to Eq. (121) were justified in the
asymptotic regime, Eq. (2). The issue is that $K/\bar{\gamma}=20$ is not far
enough into the asymptotic regime to see the accuracy of the asymptotic
analytical result.

We address this in the next plot, Fig. 9, which compares the decay rates of
the two algorithms for various values of $K$. These rates are obtained as the
slopes of straight lines fitted to data such as those shown in Fig. 8,
throwing away the first two data points to avoid transient effects. We see
that, using the rate in Eq. (121) and $H_{\Theta}$ in Eq. (122), the MOAAAR
and Greedy decoherence rates, scaled by $2K^{2}/(\breve{\gamma}\kappa^{2})$,
do appear to approach the analytic asymptotic values of $H^{\star}\approx
1.254$ and $H_{\pi/2}\approx 1.290$, respectively, for large $K$. This plot
better shows that the decoherence rates for MOAAAR are always smaller than
those for Greedy, even for relatively small $K$, where both algorithms give
almost the same decoherence rates, which can be seen in the inset of Fig. 9.

In Fig. 9, we also plot the decoherence rates calculated with an analytical
closed-form expressions, which is derived in Section VIII. This closed-form
expression agrees very well with the numerical results, especially for MOAAAR,
over almost two orders of magnitude of variation in $K$. Moreover, the closed
form expression can be shown to approach the asymptotic values as
$K\to\infty$, with a relative deviation scaling as $\bar{\gamma}/K$.

Figure 9: The decoherence rates of the data qubit for the two competing
strategies: the MOAAAR (maroon) and the Greedy4 algorithm (green), for
different values of $K$, setting $\gamma_{\uparrow}=\gamma_{\downarrow}=1.0$
and $\kappa=0.2$. Numerical results for the rates, shown as maroon disks and
the green circles, are calculated from the slopes of curves in Fig. 8
(throwing away the first two points with transient effects). Here, the solid
maroon and dashed green curves are not simply lines connecting data points,
but rather the closed-form approximation of the rates in Eq. (131) for MOAAAR
and Eq. (130) for Greedy4, respectively. The fainter maroon horizontal dashed
line near the top shows the analytical result of the MOAAAR decoherence rate
in the asymptotic regime Eq. (121) with $\Theta^{\star}$, while the fainter
green solid line above it shows the same for Greedy, with $\Theta=\pi/2$ as in
Eq. (112).

###  VII.5 Further numerical evidence for optimality of MOAAAR

In the previous section, we showed that MOAAAR performs better than other
discussed algorithms. Of course, there are infinitely many possible algorithms
that one could compare to, even restricting to the form of strategy in Eq.
(91), derived from the separation principle (Sec. IV.2) applied to the
sufficient statistics of Sec. V.4. It is because we do not know any way to
make such a comparison that we have claimed MOAAAR only to be plausibly
optimal, in the asymptotic regime. The claim of plausibility is, naturally,
subjective, but is based on our detailed study of the problem. Nevertheless,
the reader may ask whether more evidence for the optimality of MOAAAR can be
given. To address this, in this subsection we consider a more general
algorithm.

As the reader will recall, MOAAAR’s best measurement angle
$\Theta=\Theta^{\star}\approx 1.50055$ was found under the constraint in Eqs.
(113), which assumed the relationship $\tau=\Theta/K$ between the magnitude of
the measurement angle and the waiting time. This was chosen based on an
analysis of Greedy2 in the asymptotic regime, and has an intuitive
interpretation in terms of optimal state discrimination; see Sec. VI.6. In
this subsection, we relax this constraint and numerically calculate the
decoherence rate over a wide range of parameters: $\Theta\in[0,\pi]$ and
$\tau\in[1/K,6/K]$. We keep the natural adaptive strategy for the measurement
angle, $\theta_{n+1}=s_{n}\Theta$, where $s_{n}={\rm sign}(\zeta_{n})$, as any
other choice would break the symmetry between the RTP states.

Figure 10: Numerical decoherence rates for different SQ’s adaptive measurement
strategies defined by measurement angles $\theta_{n}=s_{n}\Theta$ and waiting
times $\tau_{n}=\tau$, with
$\gamma_{\uparrow}=\gamma_{\downarrow}=1,\kappa=0.2$, and $K=100$. Subplot (a)
shows the decoherence rate as a 3D colored surface plot with a colored contour
plot (projection of the surface plot on the $\Theta$–$K\tau$ plane) for
$\Theta\in[0,\pi]$ and $\tau\in[1/K,6/K]$. Different colors (on the color bar)
indicate different values of the decoherence rate divided by the factor
$\breve{\gamma}\kappa^{2}/(2K^{2})$. Note the log scale for this axis and
colour spread. The tiny black area of the plane at $\Theta\approx\pi$,
$K\tau\approx 1$ shows where we do not have extreme confidence in the
calculated rate. Subplot (b) shows a zoomed-in region, $\Theta\in[1.4,1.6]$
and $\tau\in[1.4/K,1.6/K]$, which includes a segment of the line
$\tau=\Theta/K$ (white dashed line), which was the ansatz adopted to find
MOAAAR, including the point
$\\{\Theta^{\star},\tau^{\star}\\}=\\{1.50055,1.50055/K\\}$ (red dot) defining
MOAAAR itself.

We first investigate the phase-space dynamics and find that, outside the
regime where $\tau\approx\Theta/K$, the very simple dynamics as shown in Fig.
6 or Fig. 7 are no longer evident. Instead, these strategies generally result
in highly scattered state dynamics in the phase space. This means that we can
no longer use eigenstates of the ${\bf F}$ maps to calculate the decoherence
rate. We therefore need to use the brute-force calculation following Eq. (98)
to compute the exact coherence for a given strategy from averaging over all
possible realizations of measurement results $Y_{n}$. From the coherence data
over time, similar to the plots in Fig. 8, we then compute a decoherence rate
by fitting the data with a linear function and extracting its slope.

We also note that the further the measurement strategies are from the regime
$\tau\approx\Theta/K$, the longer time it takes for the coherence process to
reach its asymptotic behaviour of exponential decay. (The “long” times we
consider are still much shorter than the reciprocal of the decay rate, which
is why using a linear fit for the long-time decoherence curve is permissible.)
This is expected as, if more points in phase space are accessible, it may take
longer to reach the ergodic regime. Thus, to avoid the transient effects, we
perform the numerical calculations for more steps than previously (in this
case, $T/\tau=15$ steps) and only use the last $5$ steps for the linear-model
regression. Also, in order to be confident about the decoherence rate
extracted from the slope, we calculate the coefficient of determination, known
as the $R^{2}$ value [54], of the linear regression fit. Its value ranges from
$0$ to $1$, with $R^{2}=1$ representing perfect fitting. We find that we can
have extremely high confidence in our linear fits ($R^{2}>0.998$) except in a
tiny region of the parameter space we explore, which is very far from the
optimal regime; see Fig. 10(a).

We present in Fig. 10 our numerical results for the decoherence rate divided
by the factor $\breve{\gamma}\kappa^{2}/(2K^{2})$, the same as in Fig. 9. In
Fig. 10(a), we show this for the full range of parameters. There is one global
minimum and another local minimum, which correspond to
$\\{\Theta,\tau\\}\approx\\{1.5,1.5/K\\}$ and
$\\{\Theta,\tau\\}\approx\\{1.5,(1.5+\pi)/K\\}$, respectively. (The pattern of
local minima is presumably repeated for $\tau\approx(1.5+n\pi)/K$ for all
integers $n$.) In order to see clearly the optimal point, we zoom into the
region $\Theta\in[1.4,1.6]$ and $\tau\in[1.4/K,1.6/K]$ in Fig. 10(b), where
the decoherence rates are calculated with finer grid values of $\Theta$ and
$\tau$. The results show that the smallest values of the decoherence rate all
occur around the line $\Theta=K\tau$ (shown as the white dashed line). We also
find that $\Theta=1.50055$ and $\tau=1.50055/K$ (shown as the red dot) indeed
gives the smallest value of
$\bar{\Gamma}/(\check{\gamma}\kappa^{2}/2K^{2})=1.244$, which agrees perfectly
with the optimal measurement strategy chosen by MOAAAR. These numerical
results thus do offer more evidence that MOAAAR is a plausibly optimal
adaptive algorithm for measuring the SQ.

##  VIII Closed-form expression for decoherence rate in non-asymptotic regime

In Section VII.4, we defined a simple adaptive strategy, MOAAAR, by finding
the optimal value $\Theta^{\star}$ that minimizes the average decoherence
rate, Eq. (118), computed from the stable eigenstates of the two $\bf F$ maps
in the asymptotic limit $K\gg\bar{\gamma}$. In this section, we show that we
can generalize this to find a closed-form expression for the average
decoherence rate even in the non-asymptotic regime. Moreover, it will apply
for Greedy4 and Greedy2 as well as MOAAAR.

We do this by taking into account the small probabilities that the system’s
states can be outside of the eigenstates. In particular, guided by the
numerical results in Fig. 5 and Fig. 6, we assume that as well as the two
stable eigenstates (on the right side of the phase space) the system also
spends time in another two less likely states on the left side of the phase
space. Using all the four states and considering the probability transfers
among them, we can then recalculate the average decoherence rate in closed
form. For simplicity, we consider a symmetric-jump case, _i.e._ ,
$\gamma_{\uparrow}=\gamma_{\downarrow}$ (though all of this can be simply
adapted to the asymmetric case), where there is a symmetry in swapping between
$s=+1$ and $s=-1$.

We first define the four states, using the maps defined in Eq. (108), which we
reproduce here for convenience:

$\mathbf{F}_{s,a}^{y}:=\mathbf{F}\left(s\Theta({a}),\Delta({a})/K,y\right).$
(125)

Recall that $a=\ <$ refers to the left side of phase space and $a=\ >$ to the
right side. The two right states (with $s=\pm 1$) are the stable (largest
absolute eigenvalue) eigenvectors, $\underline{{E}}_{s}^{0}$, of ${\bf
F}_{s,>}^{0}$. The left states are the states that these jump to when a non-
null result occurs, ${\bf F}_{s,>}^{1}\underline{{E}}_{s}^{0}$. Because we are
assuming the symmetry $\gamma_{\uparrow}=\gamma_{\downarrow}$, we need only
one right and one left state, and we define two _normalized complex vectors_
for these as

$\displaystyle\underline{{r}}$
$\displaystyle=\underline{{E}}_{+}^{0}\,/\,\underline{{I}}^{\top}\underline{{E}}_{+}^{0}\,,$
(126a) $\displaystyle\underline{{l}}$ $\displaystyle={\bf
F}_{\\!\\!+,>}^{1}~{}\underline{{r}}\,/\,\underline{{I}}^{\top}{\bf
F}_{\\!\\!+,>}^{1}~{}\underline{{r}}\,.$ (126b)

Note that $\underline{{l}}$ is a state with $s=-1$ and $a=\ <$.

In order to calculate the average decoherence rate, we also need to know the
probability that the system be in each of the four states. In the symmetric
case, we can reduce the problem to only finding probabilities of the system
being either on the right or left sides of the phase space, of which we denote
the probabilities by $p^{>}$ and $p^{<}=1-p^{>}$, respectively. To compute the
probabilities, we need the real matrices defined in Eq. (79). Using the same
compact notation as in Eq. (125), we write these as
$\check{\mathbf{F}}_{s,a}^{y}$. These matrices satisfy

$\forall s,a,\ \sum_{y}\underline{{I}}^{\top}\check{{\bf
F}}_{s,a}^{y}=\underline{{I}}^{\top}.$ (127)

That is, for all four cases $s\in\\{+1,-1\\}$ and $a\in\\{>,<\\}$,
$\sum_{y}\check{\bf F}_{s,a}^{y}$ is a transition matrix, known in stochastic
process theory as a column-stochastic matrix, with $\check{\bf F}_{s,a}^{1}$
and $\check{\bf F}_{s,a}^{0}$ known as sub-stochastic matrices [55, 56]. From
these we define _normalized probability vectors_ , $\check{\underline{{r}}}$
and $\check{\underline{{l}}}$ similar to Eqs. (126) as

$\displaystyle\check{\underline{{r}}}$
$\displaystyle=\underline{{\check{E}}}_{+}^{0}\,/\,\underline{{I}}^{\top}\underline{{\check{E}}}_{+}^{0}\,,$
(128a) $\displaystyle\check{\underline{{l}}}$ $\displaystyle={\check{\bf
F}}_{\\!\\!+,>}^{1}\check{\underline{{r}}}\,/\,\underline{{I}}^{\top}{\check{\bf
F}}_{\\!\\!+,>}^{1}~{}\check{\underline{{r}}}\,,$ (128b)

where $\check{\underline{{E}}}_{+}^{0}$ is the eigenvector of $\check{\bf
F}^{0}_{+,>}$ with the largest absolute eigenvalue. Note that the divisors
here ensure that $\check{\underline{{r}}}$ and $\check{\underline{{l}}}$ have
positive elements, which are the conditioned Bayesian probabilities for $z=\pm
1$ when the system has normalized coherence vector $\underline{{r}}$ and