\section{Introduction}
\label{sec:intro}

Quantum walks are used to model the evolution of a quantum particle, or excitation, over a discrete set of positions. 
If the walker dynamics is continuous in time, i.e.
time is a real positive parameter, we talk about continuous-time quantum walk (CTQW) \cite{farhi98,kempe003,venegas12}. In many contexts, quantum walks shows a quantum advantage with respect their classical counterpart, by allowing a speed-up in completing certain tasks. Examples are found in quantum computation \cite{childs09,kendon14comp}, and quantum algorithms \cite{gamble10,ambainis2003,kadian21,venegas2008quantum}.
Since a  CTQW evolves over a graph, it is strongly related to  applications over networks, including
quantum spatial search \cite{childs04, Janmark14, childs2014, Chakraborty20_anyGraph, paris-oracle-21,benedetti21,portugal2013quantum},
quantum routing \cite{alastair, Chudzicki10,  paganelli13},
quantum transport and state transfer\cite{Christandl05,kay10, tama16}.
The ability to redirect or control information over a graph in an efficient way is essential to develop protocols involving quantum networks and to deal with a large amount of structured data. To this aim, we develop a protocol to  guide the walker toward a target node  on a graph by exploiting the tools of quantum control.

\par
The theory of quantum control \cite{WisemanMilburn} addresses the problem of preparing a quantum system in a desired quantum state or with some desired quantum properties. In quantum feedback-control strategies the quantum system under control is measured (typically continuously in time) and the information acquired is exploited in order to optimize a feedback operation on the system itself. This kind of strategies has been studied in great detail, in particular with the aim of generating quantum states with non-classical properties such as squeezing or entanglement \cite{Wiseman1993,Wiseman1994,Doherty1999,Thomsen2002,SerafozziMancini,Szorkovszky2011,Genoni2013PRA,Genoni2015NJP,Hofer2015,Martin2015,martin2017optimal,Martin:17,Brunelli2019PRL,martin2020quantum,Jiang2020,Zhang2020,Rossi2020PRL,DiGiovanni2021} or to cool optomehcanical systems towards their ground states, with the experimental results recently demonstrated in~\cite{Rossi2018,Magrini2021,Tebbenjohanns2021}.
\par
In this article, we propose a novel approach to quantum search on graphs and establish a new and unexplored line of research combining CTQWs with quantum feedback-control protocols. In our system, the walker is interacting with an environment that is continuously monitored. As a result, the system evolution is a quantum stochastic trajectory and, based on the result of the measurement, a feedback protocol is applied. We will prove that it is possible to drive the walker towards a target state by optimizing the corresponding target fidelity at each step. 
Our method differs from  other methods that have been developed since it allows the walker to be guided continuously to the target node. 
In the  standard spatial search protocol the oracle is described as a projector operator onto the target state. With our approach we modify this paradigm. We consider a \emph{dynamical} oracle encoded in the feedback operation. 
This means that the final projective measurement on the walker, to be performed at any time $t$ after a certain threshold time $t_{th}$, has a high probability of success. In particular, we find that once the walker reaches the target node it remains stuck in it thanks to the feedback operation. This lift the burden of  performing a final measurement at a very specific time, by allowing us to measure the walker position at any time $t\geq t_{th}$.

\par
The manuscript is organized as follows: in Sec. \ref{qwalks} and in Sec. \ref{sec:qcontrol} we give a brief introduction  to continuous-time quantum walks, and  to continuously monitored quantum systems and  unitary (measurement-based) feedback strategies respectively. In Sec. \ref{sec:qtargeting}, we introduce our search scheme and we describe the idea behind our feedback protocol. In Sec. \ref{sec:results} we analyze our results for the different control strategies we have considered, while
we conclude the manuscript in Sec. \ref{sec:conclusions} with some concluding remarks and some outlooks.


\section{Quantum walks and spatial search on graphs}
\label{qwalks}

A continuous-time quantum walk  describes the continuous motion of a quantum particle over a discrete set of positions. 
Underlying every walk there is a graph $G$, which is described as
 a pair $G=(V, E)$ where $V=\{0,1,\dots,N-1\}$ is the set of vertices and $E$ is the set of undirected edges, i.e. all the pairs of adjacent vertices in $V$. For a CTQW, the  vertices represent the positions that the particle can occupy while 
the edges encode all the possible paths that a walker can move across.
We denote the order of the graph as the number of nodes $N = \vert V\vert$.
We  restrict our discussion to regular graphs, with no loops 
nor multiple edges.
All this information determines the topology of the graph and is encoded in the adjacency matrix $A$. In the position basis
$\mathcal{B}_p=\{\ket{k}\}_{k=0}^{N-1}$,
this is given by the matrix elements
\begin{equation}
    \langle i \vert A \vert j \rangle = A_{ij} =
    \begin{cases}
        1 & \textup{if } (i,j)\in E \\
        0 & \textup{if } (i,j)\not\in E
    \end{cases}.
    \label{eq:adjmat}
\end{equation}
For example, in the case of a cycle graph the only non-zero elements  are 
$A_{j,k}$ and $A_{k,j}$ satisfying the condition ${k=(j+1)\!\!\! \mod N}$, with $j=0,\dots,N-1$.
Associated to $A$ there is a Laplacian matrix, defined as 
$L = D - A$. Here $D$ is a diagonal  matrix, with 
$\langle k \vert D \vert k \rangle=d_k= \sum_{j=0}^{N-1}A_{kj}$  the degree of the $k$th vertex. The Laplacian matrix is  promoted to be the generator of the quantum dynamics and the physical origin behind this choice is the 
correspondence between $L$ and the discretized kinetic operator for regular lattices 
\cite{wong2016laplacian}. 
Hence, the quantum walker Hamiltonian is $H = \gamma L$, where the parameter $\gamma > 0$ is the hopping rate between the nodes and it accounts for the energy scale of the system.
 The generic state of the CTQW at time $t$ is a superposition of over the vertices $\ket{\psi(t)}=\sum_{k} a_k\ket{k}$ with  $a_k=\bra{k}e^{-iHt}\ket{\psi_0}$ and $\ket{\psi_0}$ the initial state of the walker.
 Throughout the paper we set $\hbar=1$.
 %
In the case of regular graphs, 
the matrix $D$ is proportional to the identity operator, thus making the evolutions generated by $L$ and $A$ equivalent.
Before proceeding, we  want to remark that using $L$ as the dynamics generator is only one of the possible (infinitely many) choices for a  Hamiltonian. Indeed any Hermitian operator can be used as a legit CTQW Hamiltonian \cite{wong2016laplacian, frigerio21}.
\\
\par
One  application of CTQWs is  the  spatial search algorithm. The goal is to exploit the coherent evolution of a quantum walker to find a marked vertex on a graph faster than its classical counterpart. 
In the quantum spatial search of a marked node $w$, the walker evolves under the Hamiltonian:
\begin{equation}
    H_S=\gamma L - \ketbra{w}{w} 
\end{equation}
where $\gamma$ is a real parameter and the operator $\ketbra{w}{w}$ is the oracle Hamiltonian i.e. a projector onto the target state. 
The walker is usually initialized in the uniform superposition of all nodes \begin{equation}
    \ket{\psi_0}=\frac{1}{\sqrt{N}}\sum_{k=0}^{N-1}\ket{k},
    \label{psi0}
\end{equation}
with no bias toward  the target state. 
The algorithm is successful if the probability of finding the target node $p_w(t_s)=|\bra{w}e^{-iH_S t_s}\ket{\psi_0}|^2 $ is as close as 1 as possible in a time $t_s=\mathcal{O}(\sqrt{N})$.
It was shown  that a $\sqrt{N}$ speedup can be obtained for specific topologies, such as the complete and hypercube graphs and $(d>4)$-dimensional lattices \cite{farhi981,childs04}. Later studies proved fast search for different kinds of graphs \cite{Chakraborty16,philipp16,wong16,wong18} and
a comprehensive analysis  of the algorithm's performances was carried out in \cite{chakraborty20}, which recovers
previous  graph-dependent results as special cases.
 It is worth mentioning here that considering different oracle operators, such as those which modify the edges connected to the target node, allows to reach
 a search time $t_s=\mathcal{O}(\sqrt{N \ln N})$  in two-dimensional ($d=2$) lattices, by building Hamiltonians that exhibit Dirac points in their dispersion relation  \cite{tanner14, childs2014}. 



\begin{figure*}
    \centering
    \includegraphics[width=\textwidth]{figures/cycle.pdf}
    \caption{Cycle graph embedded in a plane on a circle of unit radius. Left panel: $N=5$; Central panel: $N=11$; Right panel: $N=15$. 
   
    }
    \label{fig:cycleembedded}
\end{figure*}


\section{Continuous monitoring and Feedback control}
\label{sec:qcontrol}

\subsection{Continuously monitored quantum systems}
We here provide a very basic introduction to continuously monitored quantum systems. We refer to the following references \cite{WisemanMilburn,SteckJacobs,Brun2002} for a more detailed introduction and for the derivation of the formulas provided in this Section.
We assume that the quantum system under exam interacts with a large Markovian environment described by a train of input bosonic operators $\hat{a}_{j}(t)$ satisfying the canonical commutation relation $[\hat{a}_{j}(t^\prime),\hat{a}^\dagger_{k}(t)]=\delta_{jk} \delta(t-t^\prime)$. The interaction with the system is then given in terms of the time-dependent interaction Hamiltonian 
\begin{align}
\hat{H}^{int} (t) = i \sum_j \sqrt{\kappa_j}(\hat{c}_j\hat{a}^\dagger_{j}(t) + \hat{c}_j^\dagger \hat{a}_{j}(t)),
\label{eq:intHamiltonian}
\end{align}
in which $\kappa_j$ represent the coupling strengths, while $\hat{c}_j$ are operators acting on the system Hilbert space (one should also notice that the parameter $t$ in the operators $\hat{a}_j(t)$ is just a label denoting the time at which each operator is interacting with the quantum system via the Hamiltonian). 
We  also assume that the environmental modes $\hat{a}_j(t)$ can be measured continuously in time, just after the interaction, in order to gain information on the state of the system itself. Notice that the interaction with the environment can be either considered already present, and that some degree of control on this environment is achievable in order to perform such a measurement, or that such an interaction can be effectively engineered with the purpose of weakly monitoring the system. Both these approaches are nowadays pursued efficiently in different physical platforms, in particular in circuit QED \cite{Murch:2013aa,Naghiloo:2016aa,Hacohen-Gourgy2016,Ficheux2018,Minev:2019aa} and in optomechanical systems \cite{Wieczorek2015,Rossi:2018aa,Rossi2019,Magrini:2021aa,Tebbenjohanns:2021aa}. It is known that in quantum mechanics a measurement modifies the state of the quantum system that is (directly or indirectly) measured, and that the corresponding {\em conditional state} will depend both on the kind of measurement performed and on the outcome of the measurement. We will focus on continuous homodyne detection of the environmental modes with monitoring efficiencies $\eta_j$, corresponding to a set of continuous photocurrents
\begin{align}
dy_t^{(j)} = \sqrt{\eta_j \kappa_j} \Tr[(\hat{c}_j + \hat{c}^\dagger_j) \varrho^c ] \,dt + dW^{(j)}_t \,,
\end{align}
where $dW^{(j)}_t$ denotes the {\em innovation}, that is the difference between the result of the measurement $dy_t^{(j)}$ and the expected results, and mathematically correspond to independent Wiener increments satisfying $dW^{(j)}_t dW^{(k)}_t= \delta_{jk} dt$. The evolution of the quantum state $\varrho^c(t)$ conditioned on the photocurrents $dy_t^{(j)}$ is then given by the following stochastic master equation (SME)
\begin{align}
   d\varrho^c &= -i[\hat{H}_s,\varrho^{c}(t)]dt +\sum_{j}\kappa_j \mathcal{D}[\hat{c}_j] \varrho^{c}(t) dt \nonumber \\
   &\,\, +\sum_{j}\sqrt{\eta_j \kappa_j} \mathcal{H}[\hat{c}_j] \varrho^{c}(t) \,dW_t^{(j)} \,
   \label{eq:sme}
\end{align}
where $\hat{H}_s$ is the Hamiltonian describing the evolution of the quantum system only, and where we have introduced the two following superoperators
\begin{align}
    \label{eq:superop1}
    \mathcal{D}[\hat{c}]\bullet &= \hat{c}\bullet \hat{c}^\dag - (\hat{c}^\dag\hat{c}\bullet+\bullet\hat{c}^\dag\hat{c})/2 \,, \\
    \label{eq:superop2}
    \mathcal{H}[\hat{c}]\bullet &= \hat{c}\bullet + \bullet \hat{c}^\dag - \Tr[(\hat{c} + \hat{c}^\dag)\bullet]\bullet \,.
\end{align}
The continuous outcomes of the photocurrents $\{dy_t^{(j)}\}$ thus define a particular conditional trajectory for the conditional state of the quantum system. By averaging over all the possible trajectories, i.e. over all the possible outcomes of the measurements, we obtain the evolution of the unconditional state $\varrho^u = \mathbbm{E}_{\sf traj}[\varrho^c]$ that, by exploiting the property $\mathbbm{E}_{\sf traj}[dW_t^{(j)}]=0$, is a Markovian master equation in the Lindblad form
\begin{align}
    \frac{d\varrho^u}{dt} &= -i[\hat{H}_s,\varrho^{u}(t)] +\sum_{j}\kappa_j \mathcal{D}[\hat{c}_j] \varrho^{u}(t). \label{eq:uncondMQE}
\end{align}

The evolution of the conditional states described by the SME (\ref{eq:sme}) can be equivalently described via the formula \cite{rouchon2015arxiv,Rouchon2015} 
\begin{widetext}
\begin{align}
    \varrho^c(t+dt) = \frac{\hat{M}_{{\bf dy}_t} \varrho^c(t) \hat{M}_{{\bf dy}_t}^\dag + \sum_j (1-\eta_j) \hat{c}_j \varrho^c(t) \hat{c}_j^\dag \,dt }{\Tr[\hat{M}_{{\bf dy}_t} \varrho^c(t) \hat{M}_{{\bf dy}_t}^\dag + \sum_j (1-\eta_j) \hat{c}_j \varrho^c(t) \hat{c}_j^\dag \,dt]} \,,
\end{align}
\end{widetext}
where we have introduced the family of Kraus operators
\begin{equation}
    \hat{M}_{{\bf dy}_t} = \mathbb{I} -i \hat{H}_s dt - \sum_{j}\left( \frac{\kappa_j}{2} \hat{c}^\dagger_j\hat{c}_j dt - \sqrt{\eta_j \kappa_j}\, \hat{c}_j\, dy_t^{(j)}\right) \,,
    \label{eq:Kraus}
\end{equation}
with ${\bf dy}_t = \{dy_t^{(1)},...,dy_t^{({\tiny K})}\}$ denoting the vector of the outcomes of the $K$ mesaurement channels.\\
In our protocol we  consider an initial pure state and perfect monitoring efficiency, i.e. $\eta_j = 1$ for all channels. 
Under these assumptions the conditional evolution is described by a stochastic Schr{\"o}dinger equation, or equivalently via the Kraus operators as follows:
\begin{align}
    |\psi^c(t+dt)\rangle = \frac{\hat{M}_{{\bf dy}_t} |\psi^c(t)\rangle}
    {\sqrt{\langle \psi^c(t) | \hat{M}_{{\bf dy}_t}^\dag \hat{M}_{{\bf dy}_t} |\psi^c(t)\rangle}} \,.
    \label{eq:sse}
\end{align}

\subsection{Unitary quantum feedback}
In addition to conditioning the evolution of the quantum state, the outcomes of the measurement performed can  in principle be exploited to further modify the dynamics of the system. 
In this respect, here we  briefly introduce unitary measurement-based quantum feedback. The idea is that, once the measurement outcomes ${\bf dy}_t$ are obtained, one preforms a unitary operation $\hat{U}_{fb}(t)$ on the quantum state, typically optimized in order to achieve a certain goal as, for example, the preparation of a certain target quantum state. 
This unitary operation may depend only on the last measurement outcomes, or on the whole history of outcomes, and thus on the whole trajectory of the conditional state. In the first instance one talks about Markovian quantum feedback and one can derive a corresponding Markovian feedback master equation \cite{Wiseman1994,WisemanMilburn}. In this work we will focus on the second kind of feedback, and thus our feedback strategy will be optimized by knowing both the last measurement outcomes and the conditional state $|\psi^c(t)\rangle$ (and as a consequence the whole measurement history). In order to obtain the corresponding evolution, we exploit the formulas involving the Kraus operators. 
In particular, if the feedback operation is performed after the measurement, by assuming unit measuring efficiency, initial pure states and no-delay between measurement and feedback, the conditional state at each instant is described via the formula
\begin{align}
       |\psi^{fb}(t+dt)\rangle =  \frac{\hat{U}_{fb} \hat{M}_{{\bf dy}_t} |\psi^{fb}(t)\rangle}
    {\sqrt{\langle \psi^{fb}(t) | {\hat{M}_{{\bf dy}_t}}^\dag \hat{M}_{{\bf dy}_t} |\psi^{fb}(t)\rangle}} \,.
    \label{eq:purefeedback} 
\end{align}
This formula is particularly useful for our numerical approach
where one needs to substitute the time differential $dt$ with a finite but small time increment $\Delta t$, while the Wiener increment $dy_t^{(j)}$ must be replaced by a Gaussian random variable $\Delta y_t^{(j)}$ with zero mean and variance $\Delta t$. The finite increments to the measurement records are
\begin{equation}
   \Delta y_t^{(j)} = \sqrt{\kappa_j} \langle \psi^{c}(t) \vert (\hat{c}^\dagger_j + \hat{c}_j) \vert \psi^{c}(t) \rangle \Delta t + \Delta W_t^{(j)} .
\end{equation}
Due to the finite nature of $\Delta t$, the deterministic identity $\Delta W_t^{(j) 2} = \Delta t$ is no longer satisfied, thus corrections must be considered. This is accomplished by adding an extra term, known as  Euler-Millstein correction \cite{Rouchon2015}, in the Kraus operators that now read 
\begin{align}
    \hat{M}_{{\bf \Delta y}_t} & = \mathbb{I} - i \hat{H}_s \Delta t - \sum_{j=x,y} \bigg(\frac{\kappa_j}{2} \hat{c}_j^\dagger \hat{c}_j \Delta t - \nonumber \\
    & \quad -\sqrt{\kappa_j} \hat{c}_j \Delta y_t^{(j)} - \frac{\kappa_j}{2}\hat{c}_j^2({\Delta y_t^{(j)}}^2-\Delta t) \bigg) \,.
\end{align}



\section{Quantum search assisted by feedback}
\label{sec:qtargeting}
Our idea is to continuously monitor the position of a quantum walker during its evolution on a cycle graph, and then to use this information to apply feedback unitary operations as a dynamical oracle with the aim of finding a particular target node. 
The walker is initially prepared in the uniform superposition of all nodes of the graph as in Eq. \eqref{psi0}.
The first step is to describe the continuous monitoring. 
In particular we assume to be able to couple  our environment to the following jump operators 
\begin{align}
    \hat{c}_1 = \hat{x} = \sum_{k=0}^{N-1} \cos(\frac{2\pi k}{N}) \vert k \rangle \langle k \vert\,,\label{eq:jumpoperators1}  \\
	\hat{c}_2 = \hat{y} = \sum_{k=0}^{N-1} \sin(\frac{2\pi k}{N}) \vert k \rangle \langle k \vert \,,
	\label{eq:jumpoperators2}
\end{align}
whose eigenvalues exactly correspond to the coordinates of the position of the $N$ nodes of the graph, corresponding to 
equally spaced points 
on a unit radius ring centered on $(0,0)$ in the $(x,y)$ plane. See Figure \ref{fig:cycleembedded} for cycle graphs with $N=5, 11, 15$. 
We show in Appendix \ref{app:uncondME} that the unconditional evolution, corresponding to the master equation (\ref{eq:uncondMQE}) with the cycle graph Hamiltonian $H_S = \gamma L$ and 
jump operators \eqref{eq:jumpoperators1}-\eqref{eq:jumpoperators2}, leads to a symmetric dephasing-like evolution in the position basis, {thus reflecting the translation invariance of the  graph's nodes and} further validating our choice.
\\
The second  step is to define our feedback strategy. We parametrize the  unitary feedback operator as 
\begin{align}
 
   &\hat{U}_{fb}( {\boldsymbol \theta})= e^{-i \hat{H}_{fb}( {\boldsymbol \theta}) dt} \,,
\end{align}
where, in general, we  assume to have a finite number of control parameters $ {\boldsymbol \theta}=\{\theta_k\}$ corresponding to a set of control operators $\{\hat{h}_k\}$, such that the feedback Hamiltonian reads
\begin{equation}
\hat{H}_{fb}( {\boldsymbol \theta}) = \sum_k \theta_k \hat{h}_k \,.
\label{eq:feedbHamiltoniansOK}
\end{equation}
The definition of the control operators is indeed crucial. Two natural choices can be considered: the first one is to choose the on-site projectors, i.e. $\hat{h}^{(os)}_k = |k\rangle\langle k|$. However preliminary numerical simulations show that this choice is in general not efficient. 
In order to understand why this is the case, one can for example observe that, being $\{|k\rangle\}$  eigenstates of the control operators above, if the walker during the evolution happens to be in a node $|\bar{k}\rangle$ different from the target, the unitary operation will not be able to change its state and thus the feedback is useless for our purposes. The second natural choice is to consider the hopping operators 
\begin{align}
    \hat{h}^{(hop)}_k = |k\rangle\langle k+1| + |k+1\rangle\langle k| \,, \label{eq:hopping}
\end{align}
with the usual boundary condition $|N\rangle =|0\rangle$.
This set of feedback control operations represents the ability of controlling each coupling between adjacent nodes individually. 
Finally, one needs to decide how to optimize the feedback operation, that is the set of control parameters $ {\boldsymbol \theta}$, in order to find the target node on the graph.
This is typically done by defining a reward function $\Lambda(|\psi^{fb}\rangle)$, that in our case naturally corresponds to the fidelity between the conditional state after the feedback operation $|\psi^{fb}\rangle$ in Eq. (\ref{eq:purefeedback}), and the target state (that we will hereafter denote as $|0\rangle$), i.e.
\begin{equation}
\label{eq:rewardfun}
    \Lambda(|\psi^{fb}(t) \rangle) = \mathcal{F}_{\vert 0 \rangle}(|\psi^{fb}(t) \rangle) = |\langle 0 \vert \psi^{fb}(t)\rangle |^2.
\end{equation}
We will thus choose the parameters ${\boldsymbol \theta}$ as the ones maximizing the fidelity at each step of the trajectory. The same figure of merit will be then used in order to assess the performance of our protocol. We will indeed numerically evaluate
\begin{align}
    \overline{\mathcal{F}}_{\vert 0 \rangle} = \mathbbm{E}_{\sf traj}[\mathcal{F}_{\vert 0 \rangle}(|\psi^{fb}(t) \rangle) ] \,,
\end{align}
that is the fidelity averaged over all the possible trajectories conditioned by the continuous monitoring.
\section{Results}
\label{sec:results}

In the following, we  present our main results, dividing them in three different settings: in Sec. \ref{sec:unbounded} we address the numerical optimization of the feedback operation with unbounded control parameters, that is without posing any bound on the search domain for the parameters $\{\theta_i\}$. Then, in Sec. \ref{sec:multfeedbHnumOpt} we consider numerical optimization of the feedback but with a bounded domain, that from the physical point of view may represent constraints on the physical implementation of the feedback operations. In Sec. \ref{sec:multfeedbHswitch} we study the case of {\em digital feedback} \cite{riste2015}, in which the value of the feedback couplings are not only bounded, but can take values only from a discrete set (one should notice that unlike  previous examples of digital feedback \cite{riste2015}, we still have a continuous measurement output, and only the feedback operations are discrete). 
In all these examples we  consider the system initially prepared in the quantum  state defined in Eq. (\ref{psi0}), corresponding to the uniform superposition over the $N$ nodes.
\par
The algorithm we used to numerically optimize the feedback couplings is provided by the \texttt{SciPi} library, and in particular the \texttt{scipy.optimize.minimize} function \cite{scipy.opt.min}. The method used for the different strategies are the following: for the unbounded optimization (Sec. \ref{sec:unbounded}) we use the \texttt{BFGS} method; for the bounded optimization (Sec. \ref{sec:multfeedbHnumOpt}) we used the \texttt{L-BFGS-B} method. Differently, the method used for digital feedback (Sec. \ref{sec:multfeedbHswitch}) is a brute force one, i.e. we consider all the possible combinations of the finite discrete value and we select the optimal one.\\
We remark that we have also investigated the scenario with a single feedback Hamiltonian controlling collectively all the couplings via a single parameter $\theta$, i.e. via the Hamiltonian $\hat{H}_{fb} = \theta \sum_k  \hat{h}^{(hop)}_k$. However, as we  show in Appendix \ref{app:singlefeedH}, this kind of feedback is not particularly efficient for our purposes.
\subsection{Numerical optimization with unbounded controls}
\label{sec:unbounded}
We here show the results of our protocol when considering a feedback operation via the control operators introduced in Eq. (\ref{eq:hopping}) and optimized control parameters $\{\theta_k\}$ with unbounded domain. A remark is in order here: the first attempt we have pursued was to follow the approach described by Martin et al. in \cite{Zhang2020,martin2020quantum}, where the feedback operation is assumed to be infinitesimal, i.e. via couplings written as
\begin{align}
    {\boldsymbol \theta}\,dt = A \,{\bf dW} + {\bf B} \,dt \,,
\end{align}
where ${\boldsymbol \theta} = (\theta_1, \dots, \theta_K)^{\sf T}$ is the vector of the feedback couplings, ${\bf dW} = (dW_t^{(x)}, dW_t^{(y)})$ is the vector of the Wiener increments describing the measurement, while $A$ and ${\bf B}$ are respectively a $(N\times 2)$-dimensional matrix and a $2$-dimensional vector describing the feedback strategy. The optimization in this scenario could be done analytically if some conditions are fulfilled, as described in \cite{martin2020quantum} (see Appendix \ref{app:singlefeedH} for further details on this method). However for our problem we verify that these conditions are never satisfied and thus a numerical optimization with unbounded couplings has to be performed.  
\begin{figure*}
\centering
\includegraphics[width=\textwidth]{fig_noBound.pdf}
\caption{Results for multiple feedback control Hamiltonians $\hat{h}_{k}^{(hop)}$ defined in Eq. \eqref{eq:hopping}. $(a)$: average reward function $\overline{\mathcal{F}}_{\vert 0 \rangle}$ as a function of time (black dashed line: threshold $\mathcal{F}^{th} = 0.95$). $(b_{1})-(b_{4})$: averaged feedback couplings $\bar{\theta}_{k}$, corresponding to the $k$th Hamiltonian in \eqref{eq:feedbHamiltoniansOK}: $(b_{1})$ and $(b_{3})$: short time behavior; $(b_{2})$ and $(b_{4})$ time-asymptotic behavior. 
Red line: $N=5$; blue line: $N=11$; orange line: $N=15$. Other parameters: $\eta =1$, $\gamma\Delta t=0.01$, number of stochastic trajectory $N_{tj}=5000$.
}
\label{fig:multFeedUnbound}
\end{figure*}
\par
The results are depicted in Fig. \ref{fig:multFeedUnbound} for $N=5,11,15$, where we observe that the protocol is particularly efficient in reaching the target state.
This efficiency decreases as the size of the graph increases, although the threshold value is reached on rather small time scales $\gamma t <1$.
In Fig. \ref{fig:multFeedUnbound} we also consider the average feedback coupling $\bar{\theta}_0$ between  nodes $\ket{0}$ and $\ket{1}$ and the coupling $\bar{\theta}_1$ between  nodes $\ket{1}$ and $\ket{2}$ (see Fig. \ref{fig:cycleembedded} for reference): after  an initial transient, the average value of the feedback couplings reaches an asymptotic value, meaning that the feedback operation is stable after having reached the target node. In particular, we notice that the average coupling $\bar{\theta}_0$,
between nodes $\ket{0}$ and $\ket{1}$, tends to the asymptotic value $\bar{\theta}_0 \to -\gamma$. In fact,  when the protocol has  {\em almost} localized the walker in the desired node, the role of the feedback is  to try to stop the dynamics by nullifying the corresponding couplings in the Hamiltonian. We have also numerical evidence that, within numerical noise, the average feedback couplings are symmetric with respect to the $x$-axis ,i.e. $\bar{\theta}_0 = \bar{\theta}_{N-1}$, $\bar{\theta}_1 = \bar{\theta}_{N-2}$, etc., in the configuration in which the target node is placed in $(1,0)$, see Fig. \ref{fig:cycleembedded}.
The values for $\bar{\theta}_k$ reported in Fig. \ref{fig:multFeedUnbound} are  noisy, which suggests that in some trajectories larger values of the optimal couplings are present. This is confirmed by inspecting the values of the feedback couplings for a single trajectory, in which at a certain time $t$ larger feedback couplings might be needed to maintain the state with a high fidelity with respect to $\vert 0 \rangle$. We believe that this behaviour is due to the 
 stochasticity of the single random trajectory, i.e. there might be a time-step in which the measurement project the state far away from the target node, and thus a large correction is needed.

\subsection{Numerical optimization with bounded controls}
\label{sec:multfeedbHnumOpt}
As we have seen in the previous Section, not only numerical optimization of the feedback strategy is necessary, but also large absolute values of the feedback couplings might be needed to perform an efficient search. 
Hence, to test the limits of our protocol, we consider the case where the feedback couplings $\boldsymbol{\theta}$ belong to a bounded domain. In this case  each $\theta_{k}$ can take values from the interval $[-\xi \gamma, \xi \gamma]$, where we introduced the bounding (dimensionless) parameter $\xi$ that quantifies the range of  values admitted for  the feedback couplings $\boldsymbol{\theta}$. 
We  consider values of the bounding parameter $\xi\geq 1$. We have indeed numerical evidence that for $\xi<1$, that is for feedback couplings smaller than the Laplacian parameter $\gamma$, the protocol fails. As we noticed in the example above, once the walker has been localized over the target the role of the feedback is to stop the dynamics and this effect cannot be achieved efficiently if in general $|\theta_0| < \gamma$.
\par
The numerical results are provided in Fig. \ref{fig:multipleFeedbackBound} for $N=11$ (the results are qualitatively similar also for $N=5$ and $N=15$), and for different values of $\xi$ ranging between $\xi=1$ and $\xi=100$. 
As expected, as $\xi$ grows, the efficiency of the protocol improves, i.e. the minimum time $t_{th}$ required
to reach the
 threshold $\mathcal{F}_{th}$  decrease on average.
 \\
The decreasing of the optimal time $\gamma t_{th}$  for increasing values of $\xi$ can also be seen 
from Fig. \ref{fig:tthfuncxi}, where we plot the ratio of $\gamma t_{th}/\mathcal{F}_{th}$, a quantity that  corresponds to the effective time necessary to reach the target on average \cite{cattaneo18}.
We observe that this quantity in general may quantitatively depend on the chosen threshold value $\mathcal{F}_{th}$  but still it gives  always the same qualitative behaviour. 
We  see that above a certain value of $\xi$, the ratio reaches a minimum asymptotic value. 
Larger values of $\xi$ are necessary to reduce the effective time $\gamma t_{th}/\mathcal{F}_{th}$  when the size $N$ is increased. 
From Fig. \ref{fig:multipleFeedbackBound} and \ref{fig:tthfuncxi},  we notice that the order of magnitude of the $t_{th}$ is  larger if compared with the one obtained for the unbounded feedback (see Fig. \ref{fig:multFeedUnbound}).
\begin{figure*}
\centering
\includegraphics[width=\textwidth]{boundN112v}
\includegraphics[width=\textwidth]{boundN11ang}
\caption{Results for multiple feedback Hamiltonians $\hat{h}_k^{(hop)}$ with bounded control parameters $\boldsymbol\theta$ and for a graph with $N=11$ nodes. $(a)$: average reward function $\overline{\mathcal{F}}_{\vert 0 \rangle}$ as a function of time (black dashed line: threshold $\mathcal{F}^{th} = 0.95)$; $(b)$: average mean value of the position operator $\langle\hat{x}\rangle$ and $\langle\hat{y}\rangle$; $(c)$ and $(d)$: average feedback coupling $\bar{\theta}_{0}$ and $\bar{\theta}_{1}$ respectively, as a function of time (short time behavior). Insets: average feedback coupling $\bar{\theta}_{0}$ and $\bar{\theta}_{1}$ for a larger time $t$ (asymptotic behavior). 
In all  plots we used  $\xi=1$ (red line),  $\xi=5$ (blue line)  $\xi=50$ (orange line),  $\xi=100$ (purple line),  $\eta=1$, $dt=0.01$, number of stochastic trajectories $N_{tj}=5000$.
}
\label{fig:multipleFeedbackBound}
\end{figure*}
\begin{figure}
\centering
\includegraphics[width=0.48\textwidth]{compxibound}
\caption{Plot of $\gamma t_{th}/\mathcal{F}^{th}$ as a function of the bounding parameter $\xi$ for different graph sizes: $N=5$ (red line)  $N=11$ (blue line) and $N=15$ (orange line). The value $t_{th}$ is the time at which the $\mathcal{F}^{th}=0.95$ is reached on average. By changing the value of the  threshold parameter $\mathcal{F}^{th}$ the qualitative behavior of the curves does not change and the values obtained for different $\xi$ have the same order of magnitude.    We used  the same data and  the same set of parameters of Fig. \ref{fig:multipleFeedbackBound}.
}
\label{fig:tthfuncxi}
\end{figure}
\par
In Fig. \ref{fig:multipleFeedbackBound} we also report the average values of the feedback couplings $\bar{\theta}_{0}$ and $\bar{\theta}_{1}$.
Excluded a small noisy contribution,  we have numerical evidence that the time behaviour of $\bar{\theta}_{0}$ and $\bar{\theta}_{1}$  is equal to their symmetric counterpart $\bar{\theta}_{10}$ and $\bar{\theta}_{9}$. 
The qualitative behavior of $\bar{\theta}_{0}$ is the same also for the other values of $N$ and $\xi$ considered: after a first positive peak, it follows a minimum and then a second maximum, which is smaller than the first, and eventually it tends to the finite asymptotically value $\bar{\theta}_{0}\to - \gamma$, confirming our previous intuition. The situation is slightly different for $\bar{\theta}_{1}$, where for $N=11$ and $N=15$ there is a sequence of minima and maxima which asymptotically tend to a value close to $0$  while for $N=5$ the asymptotic values remains positive. %
This behaviour leads to the observation that, for smaller graphs, the couplings  $\theta_1$ and $\theta_{N-2}$ are more relevant with respect to graphs with a larger size.  The other feedback couplings are not particularly interesting, since their average is approximately $0$ everywhere, so we decide to not report them.
\begin{figure*}
\centering
\includegraphics[width=\textwidth]{singtrj1}
\includegraphics[width=\textwidth]{singtrj2}
\caption{Single trajectory for the multiple-feedback Hamiltonian protocol with bounded domain. $(a)$: reward function $\mathcal{F}_{\vert 0\rangle}(\bullet)$ as a function of time for $N=5$, $\xi =1$ (red line) and $N=11$, $\xi=5$ (blue line). In this case, we have respectively $\gamma t_{th}=1.80$ and $\gamma t_{th} = 2.23$, highlighted by the red circles. $(b)$: expectation value of the operator $\hat{c}_{x}$ and $\hat{c}_{y}$ for the same trajectory of $(a)$. Big blue dot:  nodes of the cycle graph in the plane as can be see in Fig. \ref{fig:cycleembedded}. $(c)$: feedback couplings $\boldsymbol{\theta}\,dt$ as a function of time for the $N=5$ trajectory: Red line: $\theta_{0}$; Blue line: $\theta_{1}$; Orange line: $\theta_{2}$; Purple line: $\theta_{3}$; Green line: $\theta_{4}$. $(d)$: feedback couplings $\boldsymbol{\theta} \,dt$ as a function of time for the $N=11$ trajectory: Red line: $\theta_{0}$; Blue line: $\theta_{1}$; Orange line: $\theta_{2}$. All the other feedback couplings are null (apart for some numerical noise of order $10^{-8}$) and we do not report them here. In $(c)$ and $(d)$ the vertical dashed line corresponds to the threshold time $\gamma t_{th}$ of the single trajectory. The parameters considered are $\eta=1$ and $dt=0.01$.\\
}
\label{fig:singtrjboundedxi}
\end{figure*}
\par
So far we have discussed only averaged results on a large number of trajectories. 
To understand in detail the behavior of the protocol, in Fig. \ref{fig:singtrjboundedxi} we report the results for a single stochastic trajectory with bounded domain for $N=5$  with $\xi=1$ and $N=11$  with $\xi=5$. During the transient evolution, when the feedback is driving the walker to the target node $\vert 0 \rangle$,  the value of the reward function is particularly affected by the measurement. After reaching the threshold value, the feedback operation tries to keep the walker into the target vertex. However, the stabilization procedure is not given by constant values of the feedback couplings. Similar to what we have discussed in the previous Section, this means that corrections are necessary also after having reached the target state with high fidelity.  These corrections are responsible for the noise we see in the averaged feedback couplings reported in Fig. \ref{fig:multipleFeedbackBound}.  By inspecting  the single trajectories, we also notice that most of the values taken by $\boldsymbol{\theta}$ during the evolution are either equal to $0$ or to  $\pm \gamma$.  
So it is worth exploring a scenario where the
 feedback couplings $\boldsymbol{\theta}$ belong to  discrete set of possible values.
\begin{figure*}
\centering
\includegraphics[width=\textwidth]{stepFeed.pdf}
\caption{Digital feedback control for $N=5$. Left panel: three possible value of the control parameters: $\{0,\pm \gamma \}$. Right panel: five possible value of the control parameters: $\{0,\pm \gamma,\pm \xi\gamma \}$ with $\xi=5$. Main plot: averaged fidelity $\overline{\mathcal{F}}_{\vert 0\rangle}$ as a function time. Inset: averaged feedback couplings $\overline{\boldsymbol{\theta}}$ as a function of time: Red line: $\theta_{0}$; Blue line: $\theta_{1}$; Orange line: $\theta_{2}$; Purple line: $\theta_{3}$; Green line: $\theta_{4}$. The red and the green line are superposed, as well as the the blue and the purple. The number of trajectories is $N_{tj}=5000$, and the parameters considered are $\eta=1$ and $dt=0.01$.
}
\label{fig:stepfeedback}
\end{figure*}
\subsection{Numerical optimization with digital feedback control}
\label{sec:multfeedbHswitch}
 We now explore a digital feedback protocol, where the feedback control parameters $\boldsymbol{\theta}$ are picked, at each time step, from a discrete number of values. We study this strategy only for a cycle graph of order $N=5$ since the numerical algorithm we employed is particularly demanding. \\
First, we consider only three possible values for the couplings $\theta_{k}$, belonging to the set $\{ 0, \pm \gamma \}$. 
In this case the optimization algorithm explores all the values of the reward function for all the possible combinations that the five feedback couplings $\boldsymbol{\theta}$ may realize, i.e. $5^{3}=125$ possible combination, and select the one that realizes the maximum $\mathcal{F}_{\vert 0 \rangle}(\bullet)$. In the left panel of Fig. \ref{fig:stepfeedback} we report the results obtained by repeating this algorithm at each step and averaging over $N_{tj}=5000$ trajectories. We see that the threshold value of $\mathcal{F}_{th}$ is reached for a time $\gamma t_{th}=6.40$, which is slightly larger than the value obtained via the continuous bounded protocol $\gamma t_{th} = 6.39$. Regarding the feedback couplings $\boldsymbol{\theta}$, their average values oscillate in the transient time, and after the threshold time $t_{th}$, they stabilize around asymptotic values. We found that approximately  $\bar{\theta}_0 = \bar{\theta}_4 = -\gamma$, 
while   $\bar{\theta}_{1},\bar{\theta}_{2}$ and $\bar{\theta}_{3}$ correspond in general to positive values.
\par
Then, we consider five possible values for the feedback couplings $\theta_{k}$, belonging to  the set $\{0, \pm \gamma, \pm \xi \gamma\}$, with $\xi$ playing the same role of the bounding factor we introduced before. Here, we may consider the two extra switchers as a boosted feedback operation, i.e. a larger coupling strength than the standard $\gamma$. In this case, the number of possible combination increase and it is equal to $5^{5}=3125$. The averaged results are reported in Fig. \ref{fig:stepfeedback}(b) for $\xi=5$ and for $N_{tj}=5000$ trajectories. 
The threshold value is reached, on average, by the time $\gamma t_{th}=1.12$, which in this case is slightly smaller than the continuous-bounded protocol threshold time $\gamma t_{th}=1.24$. This result is indeed unexpected since the digital feedback is an instance of  the strategies allowed by the continuous bounded domain. 
This suggests that all the algorithms we employed for the optimization of the continuous bounded domain are not particularly efficient in finding the optimal values of the parameters $\boldsymbol{\theta}$, while the brute-force spanning algorithm we considered in this Section cannot fail, as all the possible combinations are tested. 
\\

The averaged values of the feedback couplings reported in the inset of Fig. \ref{fig:stepfeedback} show that the extra switchers are considerably used in the initial stage of the evolution. As time increases, also in this case the values of $\bar{\theta}_{0}$ and $\bar{\theta}_{4}$  reach the  asymptotic value of $\bar{\theta}_0 = \bar{\theta}_4 = - \gamma$. The other couplings, instead, reach an asymptotic value larger than one and are particularly noisy. This is a sign that, in each trajectory, the walker dynamics must  often be corrected by the boosted positive feedback couplings.
\section{Conclusions}
\label{sec:conclusions}
The ability to control or manipulate the dynamics of a quantum walker over a network is important for the development of quantum computation, quantum algorithms and simulations. In this work we proposed a new protocol for searching a target node over a cycle graph by means of a continuous-time quantum walk. The CTQW  interacts with environmental bosonic modes that are continuously monitored and then a proper feedback operation is applied to drive the walker toward the target state. 
The feedback thus plays the role of a dynamic oracle, able to recognize the marked vertex and to change the values of the couplings between the nodes. In this work we analyzed and compared the performances of three different feedback strategies. In the first one, we optimized the feedback couplings without posing any bound on their values; then we considered the case of bounded control, by introducing a bounding parameter $\xi$; finally, we studied the case of digital feedback, where the optimal couplings were picked from a discrete set of values. We show how all the three strategies are able to localize the walker on the target node. In particular, as expected, the minimum time  necessary to reach a threshold  target fidelity is lower in the unbounded case, while the continuous bounded control and the digital feedback strategies achieve similar results. Furthermore, for all considered strategies, we show that once the target vertex is reached, the feedback operates to keep the walker in this position. This is an important difference with respect to standard spatial search protocols \cite{childs04,portugal2013quantum}, where the target is found, with high probability, at a specific time or in a very narrow time window. The implications are relevant, especially at the experimental and operational level, as in our protocol one does not need to perform the final position measurement at a specific time, but rather at any time larger than the known threshold.

Although in this paper we focus on the cycle graph, our scheme can be, in principle, generalized to any graph with appropriate adaptations both in the feedback operations and system dynamics, i.e. respectively by changing the feedback Hamiltonian and the system-environment coupling.



\acknowledgments
The authors acknowledge useful discussion with Francesco Albarelli. MGAP is member of INdAM-GNFM.
