\section{Introduction}

The advances of molecular dynamics  simulations~\cite{Rahman1971,Toukan1985,Marti1996}  provided valuable insights into the  microscopic structural and dynamic properties of liquids. It has been shown that for the application of thermal neutron scattering, accurate neutron incoherent inelastic scattering in water can be calculated using molecular dynamics simulated density of states without any free parameter~\cite{Noguere2021}. Such  progress implies that it is now possible to predict the incoherent process reliably prior to any experiment.

However, there exist numerical challenges when evaluating the cross sections numerically using the Gaussian approximation~\cite{Rahman1962}, for instance, the highly oscillatory and singular integrand.
To tackle the problem, physical approximations are often employed to ease such a computation.
Under these approximations, motions of atoms are strictly categorised as a particular type of diffusion and solid-like vibrations.
Therefore, available numerical routines, such as the phonon expansion method~\cite{Sjolander1958}, which are valid for systems that move around equilibrium positions, can be directly used.
The final numerical results of the combined models are often consistent with experimental cross sections in many liquids~\cite{Noguere2021}.

Apart from that, the main difficulty of the computation lies within the finite numerical dynamic range of Fourier transforms. Discrete Fast Fourier transform is typically used to correlate atomic microscopic structure and dynamics with measurable scattering cross sections, i.e. a function of momentum and energy transfers that can span tens of orders of magnitudes.
However, this numerical procedure is known to have spectral leakage~\cite{Harris1978} from a frequency bin to the others close by.
Therefore, a weak and fast decaying signal can be submerged by the leakage from an intense bin, resulting in a significantly smaller dynamic range than the precision that digitised floating-point number can offer. Even when the evaluation of the cross section is assisted with physical approximations, at large energy transfer, e.g. above \SI{1}{\eV}, the numerical error can no longer be ignored. A short-collision-time (STC) approximation is often employed to replace the noise-contaminated results~\cite{MacFarlane2010}.

In this work, we attempt to solve the numerical difficulties using mathematical techniques solely, a convolutional  discrete Fourier transform (CDFT) method is proposed as an alternative to the three physical approximations introduced earlier.

In section~\ref{sLiquid}, the theory of neutron scattering and Gaussian approximation, along with the physical approximations for easing numerical calculation, are introduced.
The convolutional discrete Fourier transform (CDFT) method is firstly introduced in section~\ref{sHDRFT},
and applied to the Gaussian approximation in section~\ref{ssLiquidHDRFT}.
Using the identical input as the state-of-the-art CAB light water cross section~\cite{Damian2013, Damianpriv}, the CDFT calculated scattering function of light water is benchmarked against differential and integral cross sections in section~\ref{sBenchmark}.
This work is concluded in section~\ref{sConclusion}.

\section{Neutron incoherent Inelastic scattering in liquids}
\label{sLiquid}

\subsection{Gaussian approximation}

In Van Hove's space-time correlations~\cite{VanHove1954}, the scattering function $S(Q,\omega)$ is given by the Fourier transform of the intermediate scattering function $F(Q,t)$, which can be approximated by a Gaussian~\cite{Rahman1962}
\begin{align}
	S(Q,\omega) &= \frac{1}{2\pi} \int e^{-i \omega t} F(Q,t) \mathrm{d} t \nonumber\\
	&= \frac{1}{2\pi} \int e^{-i \omega t}  \exp\left[- \frac {Q^2}{2}\Gamma(t)\right]  \mathrm{d} t \label{fGaussian}
\end{align}

With the fluctuation-dissipation theorem, it has been shown that, $\Gamma(t)$  in liquids can be expressed as ~\cite{Rahman1962}
\begin{equation}\label{gamma}
	\Gamma(t) = \frac{\hbar}{m} \int_0^{\infty} \mathrm{d}\omega \frac{\rho(\omega)}{\omega}
	\left[\coth\left(\frac{\hbar \omega}{2k_B T}\right)
	(1-\cos {\omega t})-i \sin{\omega t}\right]
\end{equation}
Here, $\rho(\omega)$ is the density of states, $k_B$ is the Boltzmann constant.

\subsection{Physical models assisted numerical evaluation}
\label{ssSolid}

Physical models are often employed to facilitate the evaluation of Eq.~\ref{fGaussian}.
This section introduces the conventional formulation for the evaluation briefly,  a more detailed discussion can be found in~\cite{CAB2014}.
The key technique for such evaluation is to break down the cumbersome width function $\Gamma(t)$ into two easier-to-evaluate components, i.e. a diffusion part and a solid-like vibration part.

\begin{equation}\label{vdos}
		\rho(\omega) = w_{\mathrm{d}}\rho_{\mathrm{diff}}(\omega) + w_{\mathrm{v}}\rho_{\mathrm{vib}}(\omega)
\end{equation}
Where $w_{\mathrm{d}}$ and $w_{\mathrm{v}}$ are the weights satisfy $w_{\mathrm{d}}+w_{\mathrm{v}}=1$.
The diffusion component is simply described by the Egelstaff-Schofield model~\cite{Egelstaff1962}
\begin{equation}
	\begin{split}
		\rho_{\mathrm{diff}}(\epsilon) &= \frac{4cw_t}{\pi \hbar kT} \sqrt{c^2+1/4}
		 \sinh{\frac{\epsilon}{2kT}}
		  K_1\left( \frac{\epsilon}{kT} \sqrt{c^2+1/4}\right)
	\end{split}
\end{equation}
where $c$ and $w_t$ are the diffusion constant and translational weight, respectively.
\begin{equation}
	c=\frac{M_{\mathrm{diff}}D}{\hbar}
\end{equation}

The diffusion coefficient $D$ that equals $\pi k T \rho(0)/2m$ is correlated linearly with $\rho(0)$, and $M_{\mathrm{diff}}$ is the diffusion mass.

For the solid-like part, the model that describes neutron-phonon interaction in harmonic crystals are used~\cite{Sjolander1958}.

\begin{align}\label{solidF}
	F(Q,t) &= \exp{\left[-\frac{Q^2}{2} \gamma(0)\right]}\exp{\left[\frac{Q^2}{2} \gamma(t)\right]} \\
	&= \exp(-2W)\exp{\left[\frac{Q^2}{2} \gamma(t)\right]}
\end{align}
where $\exp(-2W)$ is the Debye-Waller factor, and
\begin{equation} \label{fsolidgamma}
	\gamma(t) = \frac{\hbar}{m} \int_0^{\infty} \mathrm{d}\omega \frac{\rho_{\mathrm{vib}}(\omega)}{\omega}
	\left[\coth\left(\frac{\hbar \omega}{2kT}\right)\cos{\omega t}+i\sin{\omega t} \right]
\end{equation}

Towards zero energy, for any three-dimensional solid material, the density of states always decay according to a power law~\cite{Sears1995}.
Therefore, the integrand of Eq.~\ref{fsolidgamma} has no singularity at $\omega = 0$.
Indicating by Eq.~\ref{vdos}, to make this combination of physic approximations work, it is crucial to subtract a suitable portion of the diffusive density of states, i.e. $\rho_{\mathrm{diff}}$, from the original liquid state to obtain an $\rho_{\mathrm{vib}}$ that satisfies the power law at small energies, i.e. $a\omega^2$ .

At large energy and momentum transfers, the scattering function became a Gaussian without any detailed structures. Due to the numerical difficulties, conventional evaluation  is likely to fail in that region, hence the short-collision time approximation is often used.

\begin{equation} \label{eSCT}
	S(Q, \omega) = \sqrt{\frac{1}{4\pi k_b T E_r}} \exp\left[ -\frac{(\hbar\omega - E_r)^2}{4\pi E_r}  \right] \,\,\,\,\, \mathrm{, where} \,\,\,\,\, E_r=\frac{\hbar^2 Q^2}{2m}
\end{equation}

\section{The convolutional discrete Fourier transform (CDFT) method}
\label{sHDRFT}

The concept of the CDFT method is based on the Taylor expansion of the function for Fourier transform.
We discuss the case that the function of interest is an exponential function $\exp[-f(t)]$, and $f(t)>0$.
The Fourier transform of that can be expanded as
\begin{align}\label{fft}
	F(\omega) &= \int e^{-i\omega t} \exp[-f(t)] \mathrm{d} t \\
	&= \sum_{n=0}^{\infty} \frac{1}{n!} \int e^{-i\omega t} [-f(t)]^n \mathrm{d} t \\
\end{align}

From the Convolution theorem, it is straight forward to have
\begin{equation}\label{ConvTheo}
\int e^{-i\omega t}  f_1(t) f_2(t) \mathrm{d} t =  F_1(\omega)  \otimes   F_2(\omega)
\end{equation}
Therefore, it is obvious to see that the higher order terms of the summation can be obtained by convoluting the results of lower order terms.
However, this method is not numerically stable. The integral area of the convoluted result is the produce of the areas of the two input functions. After a great number of convolutions, the results can exceed the limiting value that can be expressed by a floating number. Therefore, it is desired to use normalised functions as input.

According to Bochner's theorem, the Fourier transform of a probability measure on $f(t)$ is necessarily a normalized continuous positive-definite function.
For a finite discrete function, to make the function continuous in the overall domain of $t$, it is required to change the function to make it zero at both boundaries.
For an event function, it is obvious to get $r(t)=f_{max}-f(t)$. Eq.~\ref{fft} can then be rearranged as
\begin{equation}\label{eConvFirstStep}
	\begin{split}
		F(\omega) &= e^{-f_{max}} \int e^{-i\omega t} e^{r(t)} \mathrm{d} t \\
		&= e^{-f_{max}} \sum_{n=0}^{\infty} \frac{r^n(0)}{n!}  \int e^{-i\omega t} \left[ \frac{r(t)}{r(0)} \right]^n \mathrm{d} t 	\end{split}
\end{equation}

Denoting $g_n(\omega)$ function as
\begin{equation}\label{gnw}
 	g_n(\omega)=\int e^{-i\omega t}  \left[ \frac{r(t)}{r(0)} \right]^n  \mathrm{d} t
 \end{equation}

It can be see immediately that $g_0$ is a $\mathrm{sinc}$ function and not involved in the convolution procedure.
On the other hand, the integral of $g_1$ is normalised (see Eq.~\ref{g1normal}) and is the fundamental input of the convolutions.
 \begin{equation}\label{g1normal}
 \begin{split}
	\int g_1(\omega) d \omega &= \iint e^{-i \omega t} \frac{r(t)}{r(0)} \mathrm{d} t  \mathrm{d}\omega \\
	&= \int \delta (t) \frac{r(t)}{r(0)} \mathrm{d} t \\
	&= 1
 \end{split}
 \end{equation}

Therefore, the integral of higher order $g$ functions are also unity.

As $r(0)=f_{max}$, Eq.~\ref{eConvFirstStep} becomes
 \begin{equation}
	 F(\omega) = e^{-f_{max}} \sum_0^{\infty} \frac{(f_{max})^n}{n!} g_n(\omega)

\end{equation}

Until now, the numerical procedure is straight forward. Here we introduced the distortion factor, which is used to maximise the dynamic range of the convolution.

We defined an distort function as $e^{aw}$. It has the relation with the convolution as
\begin{equation}\label{eRelation}
	\begin{split}
	F_1(\omega) \otimes F_2(\omega) &= \int F_1(\mu) \cdot F_2(\omega-\mu) d\mu \\
		&= \int e^{a\mu} \hat{F_1}(\mu) \cdot e^{a(\omega-\mu)} \hat{F_2}(\omega-\mu) d\mu \\
		&= e^{a\omega} \int \hat{F_1}(\mu) \cdot \hat{F_2}(\omega-\mu) d\mu \\
		&= e^{a\omega}\hat{F_1}(\mu) \otimes \hat{F_2}(\omega-\mu)
	\end{split}
\end{equation}

Denoting
\begin{equation}
	\hat{g}_n(\omega)  = e^{-a\omega} g_n(\omega)
\end{equation}

Using Eq.~\ref{eRelation}, we finally have,
\begin{equation}\label{fft_final}
	F(\omega) = \exp(-f_{max}+a\omega) \sum_{n=0}^{\infty} \frac{(f_{max})^n}{n!} \hat{g}_n(\omega)

\end{equation}
where $a$ is a constant, which only affects the dynamic range of the final convolution. To optimise the result for a higher dynamic range, the CDFT method evaluates the Fourier transform by convoluting the distorted function $\hat{g}_n(\omega)$, instead of the original function $g_n(\omega)$.
The significant improvement of such procedure is introduced in section~\ref{sBenchmark}.

\subsection{Applying the CDFT Method to the cross section calculation}
\label{ssLiquidHDRFT}
	Recall that
\begin{equation}\label{esqwmaster}
	S(Q,\omega) = \frac{1}{2\pi} \int e^{-i \omega t}
	\exp\left(- \frac {Q^2}{2}\Gamma(t)\right) \mathrm{d} t
\end{equation}

Using Eq.~\ref{eConvFirstStep}, the integration becomes
\begin{equation}\label{eFFTequv}
		S(Q,\omega) = \exp\left({-\frac{\Gamma_{max}Q^2}{2}}\right)
		\sum_{n=0}^{\infty} \left(\frac{\Gamma_{max}Q^2}{2}\right)^n \frac{1}{n!}
		\int e^{-i\omega t}\left[\frac{r(t)}{r(0)}\right]^n \mathrm{d} t
\end{equation}
where $r(t)=\Gamma_{max}-\Gamma(t)$. Numerically, this equation is equivalent to the direct fast Fourier transform.

Based on Eq.~\ref{fft_final}, the integration can also be expressed as
\begin{equation}\label{eCDFT}
	S(Q, \omega) = \exp{\left(-\frac{\Gamma_{max}Q^2}{2}+a \omega\right)}
	\sum_{n=0}^{\infty} \left(\frac{\Gamma_{max}Q^2}{2}\right)^n \frac{1}{n!} \hat{g}_n(\omega)
\end{equation}
Note that for our typical calculations are room temperature for light water, we have found that the distortion coefficient $a=\omega/3k_b T$ works quite well.

At large momentum transfer, the factor  $\exp (\Gamma_{max} Q^2)$ in Eq.~\ref{eCDFT} may overflow. To prevent that, we apply the convolutional method again.
Denoting $h(\omega)=S(Q,\omega)$, at a given $Q$ where the scattering function is obtained by Eq.~\ref{eCDFT}, according to Eq.~\ref{esqwmaster}, the scattering function at higher Q values can be calculated as
\begin{equation}\label{finalsqw}
	 S(\sqrt{2^{n-1}} Q ,\omega) = \underbrace{h(\omega) \otimes h(\omega) \dots \otimes h(\omega)}_{n}

\end{equation}

Notice that the self-scattering function satisfies the sum-rule $\int S(Q,\omega) \mathrm{d} \omega=1$, hence a large number of $h(\omega)$ functions can be convoluted without encountering the numerical overflowing problem.
The distortion factor, similar to Eq.~\ref{eCDFT}, are also applied in our numerical implementation.

\section{Results and discussion}
\label{sBenchmark}

\subsection{Evaluation of the self-scattering function}
The density of states (DOS) for both hydrogen and oxygen in light water for the CAB model are used as the input of our calculation~\cite{Marquez2013}.
The data, shown in Fig.~\ref{fig:vdos}, are provided by the one of the authors of the CAB model~\cite{Damianpriv}.

\begin{figure}[H]
\centering
\includegraphics[width={0.5\textwidth},page=1]{vdos_gamma.pdf}
\caption{The density of states for hydrogen and oxygen in light water for CAB model at \SI{297}{\K}~\cite{Marquez2013}. }
\label{fig:vdos}
\end{figure}

The $\Gamma(t)$ functions defined in Eq.~\ref{gamma} are evaluated using the Filon method~\cite{Filon1968} in the range between \si{-200}\si{\pico\second} and \si{200}\si{\pico\second},
of which the detailed implementation can be found in~\ref{sFilonMethod}.
For light water, a molecular dynamic simulation of \si{100}\si{\pico\second} is able to capture the important slow dynamics for scattering function~\cite{CAB2014}, the time range chose in this work is considered to be adequate.
The result of $\Gamma(t)$ for hydrogen at \SI{297}{\K} is shown in Fig.~\ref{fig:gamma}.
Complex fast dynamics are  concentrated in the region less than \si{1}\si{\pico\second}.
The real part of $\Gamma(t)$ increases linearly when the physics is mainly governed by the diffusion process at large times.
\begin{figure}[H]
\centering
\includegraphics[width={0.5\textwidth},page=2]{vdos_gamma.pdf}
\caption{Width function $\Gamma(t)$ for hydrogen at \SI{297}{\K}.}
\label{fig:gamma}
\end{figure}

Fig.~\ref{fig:sqw_smallQ} compares the numerical results of three different methods for the Fourier transform. The reduced momentum transfer is \si{1}\si{\AA^{-1}}.
At small energy transfers, results from three methods are all agreed. When the absolute value of the momentum transfer is greater than \num{2}\si{\eV}, the scattering function calculated by the CDFT method showed a much higher dynamic range than the others.

\begin{figure}[H]
\centering
\includegraphics[width={0.5\textwidth}]{sqw_smallQ.pdf}
\caption{Comparison of calculated scattering function from direct fast Fourier transform (i.e. Eq.~\ref{esqwmaster}), CDFT without distortion (i.e. Eq.~\ref{eFFTequv}) and with distortion (i.e. Eq.~\ref{eCDFT}), $Q=\si{1}\si{\AA^{-1}}$  }
\label{fig:sqw_smallQ}
\end{figure}

When the momentum transfer is very large, the short-collision-time (SCT) approximation becomes valid.
Fig.~\ref{fig:sqw_largeQ} compares the numerical results from direct fast Fourier transform, CDFT method and the prediction of the SCT model, which is only valid when moment transfer is high.
The result from the CDFT method showed good agreement with the prediction of SCT method. On the other hand, the results from direct fast Fourier transform agree with neither of the other results, suggesting such a method could suffer from significant error when the momentum transfer is high.
\begin{figure}[H]
\centering
\begin{minipage}[t]{0.45\linewidth}
	\includegraphics[width=\textwidth, page=2]{sqw_largeQ.pdf}
\end{minipage}
\begin{minipage}[t]{0.45\linewidth}
	\includegraphics[width=\textwidth, page=3]{sqw_largeQ.pdf}
\end{minipage}
\caption{Comparison of calculated scattering function from direct fast Fourier transform (i.e.~\ref{esqwmaster}), CDFT with distortion (i.e. Eq.~\ref{finalsqw}) and the prediction of the SCT approximation (i.e. Eq.~\ref{eSCT}), at larger Q.}
\label{fig:sqw_largeQ}
\end{figure}

The complete scattering functions for light water molecule is shown in Fig.~\ref{fig:hotSqw}. For the generation of the data, Eq.~\ref{eCDFT} and Eq.~\ref{finalsqw} are used to calculate the function below and above $Q=\SI{0.3}{\AA^{-1}}$, respectively.
The peak position of the SCT approximation (i.e. Eq.~\ref{eSCT}) at given Q are shown as broken lines. These lines are well aligned with the peaked intensities  contributed by hydrogen and oxygen.

\begin{figure}[H]
\centering
\includegraphics[width={0.5\textwidth}]{sqw_hot.png}
\caption{The scattering function of light water at \SI{297}{\K} calculated using Eq.~\ref{eCDFT} and Eq.~\ref{finalsqw} }
\label{fig:hotSqw}
\end{figure}

\subsection{Quasi-elastic scattering}

Fig.~\ref{fig:hwhm} compares the half width at half maximum (HWHM) of the scattering function at small momentum transfers. Our calculated results are slightly greater than that from a simple isotropic three-dimensional diffusion model (or the Fick's law, see for example chapter 5.4 in~\cite{Squires2012}).
According to that model, the scattering function at small energy and momentum transfers can be described as a Lorentzian
\begin{align}\label{Lorentzian}
S(Q,\omega) = \frac{1}{\pi \hbar}\frac{DQ^2}{(DQ^2)^2+\omega^2}
\end{align}

\begin{figure}[H]
\centering
\includegraphics[width={0.5\textwidth}]{hwhm.pdf}
\caption{Comparison of the HWHM from CDFT and Fick's law.}
\label{fig:hwhm}
\end{figure}

The discrepancies are because, unlike the Fick's law, the HWHM from our calculation is not contributed by the diffusive process alone~\cite{Levesque1970}. Fig.~\ref{fig:sqwOrder} shows the contributions to the scattering function from different scattering orders. At small energies, multiple scatterings of diffusive motion and energetic phonons can also contribute to the energy region that the diffusive peak allocated. For example, in a two phonon scattering case, a neutron absorbs and emits phonons of very similar energies.
As can be observed in the figure, the contribution from higher order scattering is with a broader HWHM.

The intensity of those peaks grow quickly at greater momentum transfer and the shape diffusion peak will eventually be submerged by the multiple scattering contributions. In that case, the short-collision-time approximation becomes valid.
\begin{figure}[H]
\centering
\includegraphics[width={0.5\textwidth}]{sqw_difforder.pdf}
\caption{Comparison of scattering function with different scattering orders, $Q=1$.}
\label{fig:sqwOrder}
\end{figure}

\subsection{Benchmarking against experimental data}
The double differential scattering cross section for light water is computed and compared with experiment results from Esch~\cite{Esch1971}.
A comparison for $E_0 = \SI{0.154}{\eV}$ and $E_0 = \SI{0.632}{\eV}$  with different angles were shown in Fig.~\ref{fig:ddxs_154} and Fig.~\ref{fig:ddxs_632}. Excellent agreements are observed.

\begin{figure}[H]
\centering
\begin{minipage}[t]{0.45\linewidth}
	\includegraphics[width=\textwidth, page=3]{ddcs.pdf}
\end{minipage}
\begin{minipage}[t]{0.45\linewidth}
	\includegraphics[width=\textwidth, page=4]{ddcs.pdf}
\end{minipage}
\caption{Double differential scattering cross section for light water, $E_0 = \SI{0.154}{\eV}$.}
\label{fig:ddxs_154}
\end{figure}

\begin{figure}[H]
\centering
\begin{minipage}[t]{0.45\linewidth}
	\includegraphics[width=\textwidth, page=1]{ddcs.pdf}
\end{minipage}
\begin{minipage}[t]{0.45\linewidth}
	\includegraphics[width=\textwidth, page=2]{ddcs.pdf}
\end{minipage}
\caption{Double differential scattering cross section for light water, $E_0 = \SI{0.632}{\eV}$ .}
\label{fig:ddxs_632}
\end{figure}

	Total cross section is calculated from \SI{1}{\mu\eV} to \SI{10}{\eV} at \SI{297}{\K} to assess  the accuracy of the CDFT method.
Fig.~\ref{fig:totxs} shows the comparison of several results from both theoretical calculation and experiment, which are from Heinloth (1961)~\cite{Heinloth1961} and russell (1966)~\cite{Russell1966}.
Along with the CAB model cross section, the total cross section from the CDFT method are in good agreement with experimental data.

\begin{figure}[H]
\centering
\includegraphics[width={0.5\textwidth}]{xs.pdf}
\caption{Total cross section for light water at \SI{297}{\K}.}
\label{fig:totxs}
\end{figure}

\section{Conclusion}
\label{sConclusion}
A numerical procedure, called convolutional discrete Fourier transform (CDFT) method, for evaluating the Fourier transform of an exponential function is proposed. We showed that this method can accurately evaluate neutron incoherent inelastic scattering function in light water. The computed scattering function is consistent with the state-of-the-art CAB model. Instead of physical approximations assisted evaluation of the scattering function, the CDFT method relies on no additional physical approximations. The results from the CDFT method should provide more confidence when working with a new liquid, on with the validity of the commonly used physical approximations are unknown.

The CDFT method will be released in an open source package, along side with the NCrystal project~\cite{Cai2020}.

\section*{Acknowledgements}
This research is supported by the National Natural Science Foundation of China (Grant No.12075266). The authors are grateful to Dr. Damian, J. I. Marquez for valuable discussions and substantial support.


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

\IEEEPARstart{Q}{uantum} computation and quantum information~\cite{nielsen2002quantum} are no longer promising research fields but current realities whose applicability is and will be increasing in the next decade, in particular, that related to computational intelligence~\cite{manju2014,Nguyen2020}.  Quantum computation is driven by the so-called quantum bits (qubits), a quantum generalization of the classical bit. The two basic states of a qubit are $\ket{0}$ and $\ket{1}$,  corresponding with the states zero and one, respectively, of a classical bit. A qubit $\ket{\Psi}$ generalizes its classical counterpart because it allows states formed by the superposition of $\ket{0}$ and $\ket{1}$, namely, $\ket{\Psi}=\alpha \ket{0} + \beta \ket{1}$, where $\alpha$ and $\beta$ are complex coefficients. The measurement of a qubit in superposition state involves that it will collapse to one of its basic states, but there is no way to determine in advance which one; the unique available information is that the probability of $\ket{0}$ is $|\alpha|^2$ and the probability of $\ket{1}$ is $|\beta|^2$, hence, $|\alpha|^2 + |\beta|^2=1$. The main operation when dealing with qubits is the unitary transformation $U$. When $U$ is applied to a superposition state, the result is another superposition state that superposes all basis vectors. This is an appealing characteristic of unitary transformations, which is called quantum parallelism because it can be employed to evaluate the different values of a function f(x) for a given input x at the same time. The unitary transformation $U(t,t_0)$ evolves a qubit state $\ket{\Psi(t_0)}$ to $\ket{\Psi(t)}=U(t,t_0)\ket{\Psi(t_0)}=\mathcal{T}\exp[-(i/\hbar)\int_{t_0}^tH(t')dt']\ket{\Psi(t_0)}$, where $H(t')$ is its Hamiltonian. Consequently, quantum control arises as the most critical problem in realizing quantum computation. It aims at designing a time-dependent Hamiltonian $H(t)$, which drives the qubit to its target by a unitary transform. It can be indeed trivial, e.g., one may have the quantum NOT gate by a resonant pulse $H=\hbar\Omega\sigma_x/2$, flipping a qubit within operation time $T=\hbar\pi/2\Omega$. However, the solution is far from optimal since it is not robust, i.e., any slight systematic error $T\rightarrow T+\delta T$ or equivalently $\Omega\rightarrow\Omega+\delta\Omega$ will lead to fidelity loss. Meanwhile, qubits can never be perfectly isolated from the external environment, where quantum noises induce decoherence. In this way, optimal quantum control enables the implementations of high-fidelity and high-robustness gate operations, which will be the milestones in fault-tolerant universal quantum computation~\cite{shor1996fault,preskill1998fault,gottesman2010introduction}.

Physicists suggest many protocols for such goals, including adiabatic quantum evolution~\cite{RevModPhys.79.53}, composite pulses~\cite{PhysRevA.70.052318,PhysRevLett.106.233001,rong2015experimental},  pulse-shaping engineering~\cite{PhysRevA.75.062326,PhysRevLett.109.060401,PhysRevLett.111.050404},  shortcuts to adiabaticity~\cite{RevModPhys.91.045001,torrontegui2013shortcuts,PhysRevLett.104.063002}, etc. Particularly, machine learning algorithms (ML) are combined with them for further optimizations~\cite{PhysRevApplied.6.054005,PhysRevLett.123.100501,PhysRevA.103.L040401,ai2021experimentally,PhysRevX.11.031070}. It is also natural to consider applying reinforcement learning (RL) individually for quantum control tasks~\cite{PhysRevX.8.031086,porotti2019coherent,niu2019universal,dalgaard2020global,PhysRevA.97.052333,PhysRevA.99.042327,ostaszewski2019approximation,PhysRevX.11.031070,PhysRevLett.127.190403,chunlin2014,MARTINGUERRERO2022457,app11188589}. In the past few years, deep reinforcement learning (DRL) has solved pulse design for fast and robust quantum state preparation~\cite{henson2018approaching,zhang2019does,haug2020classifying}, gate operation~\cite{an2019deep}, and quantum Szilard engine~\cite{PhysRevA.100.042314}. However, as we highlighted in our previous researches~\cite{PhysRevA.103.L040401,ai2021experimentally}, we are still far from exploiting the power of RL for quantum control because of quantum measurement. RL requires the observation of states for outputting an action, conflicting with quantum mechanics' fundamental feature that the state is destroyed after direct quantum measurement. Most RL models for quantum control are trained in numerical environments instead of real quantum devices for saving resources, followed by fixed pulses after observation and evaluation. These fixed pulses hardly prevail over gradient-based optimization methods. Another approach combines the model with a quantum environment for evaluation, outputting instant action after observing a state. Even though the quantum state is discarded after extraction of quantum information by direct measurement,  the RL model can still realize a quantum control with lack of efficiency if historical actions are stored for repetitive operation of $n(n+1)/2$ steps, retrieving the last destroyed state, where $n$ is the maximum number of timesteps in each episode.

In this work, we present a stream learning approach to quantum control by employing the RL algorithm for closed-loop quantum control. In this paradigm, qubit's wave functions are no longer destroyed but slightly perturbed after information extraction via weak measurement. The model observes the state, which contains weak values as the partial information of the qubit with less confidence, resulting in an action to evolve the quantum environment to the next timestep. Our scheme reflects the spirit of stream learning once the length of each timestep is sufficiently small, resembling the dynamics of continuous measurement. It also enables transfer learning by adapting the model to the environment during the evaluation while external noises patterns are changing. We reckon that stream learning in the quantum realm enhances the performance of quantum computing and quantum information processing in real-time experiments, accelerating its development from noisy intermediate-scale devices to the next level.

\begin{figure*}[t]
\centering
\includegraphics[width=1\linewidth]{pipulse.pdf}
\caption{The expectations on $Y$ (red) and $Z$ (blue) directions of repetitively measured qubits, when a weak measurement takes place per (a) $\Delta t/T=0.1$, (b) $\Delta t/T=0.01$, and $\Delta t/T=0.001$. The qubit is driven by a resonant $\pi$-pulse, whose measurement-free dynamics is plotted in dashed curves. By averaging over the trajectories in the ensemble of $N=20$ qubits (shaded curves), we have the average trajectories plotted in solid curves. For the closed-loop quantum control, the characteristics of dynamics vary from each other by the scale of the measurement interval. The corresponding weak value feedbacks $q_0$ from measuring the inaccurate Gaussian pointer $\sigma=10$ are also recorded and averaged, showing a low signal-to-noise ratio.  Additionally, we provide the dynamics of a qubit without continuous measurement, being driven by $\pi$-pulse, for comparison purposes. The corresponding expectations are plotted by dashed curves.}
\label{fig:pipulse}
\end{figure*}

\section{Physics models}
\label{sec:models}

\subsection{Open quantum system}
The dynamics of isolated quantum systems are always governed by Schr\"odinger equation $i\hbar\partial_t\ket{\psi(t)}=H(t)\ket{\psi(t)}$, where the operator $H(t)$ is the Hamiltonian of the quantum system, with its expectation be in the unit of energy. It is equivalent to the so-called von Neumann equation $\dot{\rho}(t)=-(i/\hbar)[H(t),\rho(t)]$, where the pure state wave function $\ket{\psi(t)}$ is extended to density matrix $\rho(t)=\ket{\psi(t)}\bra{\psi(t)}$. However, the isolated quantum system never exists in real world. The external environment affects the quantum system by coupling itself to the system, inducing undesired dynamics, such as decoherence. Generally speaking, one can always write down the total Hamiltonian $H_{T}=H_{S}+H_{E}+H_{I}$, including the system Hamiltonian $H_S$, the environment Hamiltonian $H_E$, and the coupling interaction Hamiltonian $H_I$. The new system is assumed to follow the von Neumann equation. One retrieves the information of the original system by tracing out the environmental subsystem $\rho=\text{Tr}_{E}(\rho_T)$, resulting in the Lindblad master equation
\begin{eqnarray}
\label{eq:lindblad}
\dot{\rho}(t)&=&-\frac{i}{\hbar}[H(t),\rho(t)]\nonumber\\
&+&\sum_n\frac{1}{2}\left[2C_n\rho(t)C^\dag_n-\rho(t)C^\dag_nC_n-C^\dag_nC_n\rho(t)\right],
\end{eqnarray}
where $C_n=\sqrt{\gamma_n}A_n$ are the collapse operators, $A_n$ are the operators that couples the system to environment in $H_I$, and $\gamma_n$ are the corresponding rates.  The density matrix is assumed to be initially in the product state as $\rho_T(0)=\rho(0)\otimes\rho_E(0)$, i.e., the original system and the environment are not correlated at $t=0$. They still remains separable $\rho_T(t)\approx\rho(t)\otimes\rho_E$ during the evolution since the environment does not evolve significantly. The environment is considered to be Markovian, requiring the fast decays of its correlation functions than those of the system.  it is also worthwhile to mention that the evolution of an open quantum system is no longer unitary. Therefore, the density matrix of the quantum system is not a pure state, but a mixed state $\rho=\sum_np_n\ket{\psi_n}\bra{\psi_n}$ instead, with $p_n$ be the classical probability of being in $\ket{\psi_n}$ state.

\subsection{Weak measurement and continuous measurement}
A major difficulty of applying ML algorithms in the quantum regime is caused by measurement, which is usually cost-less in the classical realm. The quantum system is destroyed once measurement happens, being projected to an eigenstate once quantum information is extracted. Measuring a wave function by operator $\hat{A}$ outputs eigenvalues, whose expectation follows $\langle\hat{A}\rangle=\bra{\psi}\hat{A}\ket{\psi}$. It can also be expressed in the language of density matrix as $\langle\hat{A}\rangle=\text{Tr}(\rho\hat{A})$. Aharonov's work~\cite{PhysRevLett.60.1351} proposed its extension, which extracts partial information from the quantum system without destroying it. The weak value $A_w=\bra{\psi_f}\hat{A}\ket{\psi_i}/\langle\psi_f\ket{\psi_i}$ is no longer real eigenvalues of the operator, but exotic values instead or even complex, where $\ket{\psi_i}$ and $\ket{\psi_f}$ are pre/post-selected states. The post-selection operation does not always succeed, and the wave function is discarded once the operation fails. Thus, we couple the quantum system to a pointer for entanglement and measure the pointer projectively for a weak value, which is actually the original framework proposed by Aharonov, instead of the later developed pre/post-selection formalism. To be more specific, a Gaussian pointer $\ket{\Phi}=\int(2\pi\sigma^2)^{-1/4}\exp(-q^2/4\sigma^2)\ket{q}dq$ is coupled to the qubit $\ket{\Psi}=[\cos(\alpha/2),\sin(\alpha/2)]^{\text{T}}$, following the interaction Hamiltonian $H_{int}=g(t)p\otimes\hat{A}$, where $\sigma$ is the standard deviation of the pointer's position, $p$ is its conjugate momentum operator, and $g(t)$ is the coupling strength. A non-correlated initial state $\ket{\Phi(q)}\otimes\ket{\Psi}$ is evolved by the Hamiltonian, entangling as $\cos(\alpha/2)\ket{\Phi(q-a_1)}\otimes\ket{a_1}+\sin(\alpha/2)\ket{\Phi(q-a_2)}\otimes\ket{a_2}$, where $a_i$ and $\ket{a_i}$ are the eigenvalues and eigenstates of the operator $\hat{A}$ to be weakly measured, respectively, if $\int_0^{t_0} g(t)dt=1$. If one aims at performing a weak measurement on the Z direction, i.e., the Pauli-Z operator $\hat{A}=\hat{\sigma_z}$, the measurement outputs of the pointer's position follow the probability distribution
\begin{equation}
P(q)\approx\frac{1}{\sqrt{2\pi\sigma^2}}\exp\left[-\frac{(q-\cos\alpha)^2}{2\sigma^2}\right],
\end{equation}
shifting a displacement of the expectation $\bra{\Psi}\hat{\sigma_z}\ket{\Psi}=\cos\alpha$. Correspondingly, the wave function of the qubit is slightly perturbed as
\begin{eqnarray}
\label{eq:perturbed}
\ket{\tilde{\Psi}}\propto\frac{1}{(2\pi\sigma^2)^{1/4}}\left\{\cos\frac{\alpha}{2}\exp\left[-\frac{(q_0-1)^2}{4\sigma^2}\right]\ket{0}\right. \nonumber\\
\left.+\sin\frac{\alpha}{2}\exp\left[-\frac{(q_0+1)^2}{4\sigma^2}\right]\ket{1}\right\},
\end{eqnarray}
if the weak value $q_0$ is the measurement feedback of the projective measurement.

Furthermore, quantum information can be continuously extracted from the quantum system, realizing a continuous measurement as the information obtained per measurement goes to zero. In this framework, the total operation time is divided into intervals of timestep $\Delta t$,  so that a weak measurement is performed in each interval. The limit $\Delta t\rightarrow 0$ results in continuous measurement, in which stochastic differential equations govern its dynamics~\cite{gross2018qubit,jacobs2006straightforward}. In Fig.~\ref{fig:pipulse}, we illustrate the dynamics of stochastic Schr\"odinger equations, aiming at flipping a qubit by fixed resonant $\pi$-pulse, with different scales of time interval $\Delta t$.

\section{Numerical experiments}
\label{sec:ne}
\subsection{Physical system and task}

As we introduced in Sec.~\ref{sec:models}, systems under quantum noises are no longer governed by Schr\"odinger equation, but Lindblad master equation instead. For pure dephasing, the diagonal Lindblad operators reads $C_n=\sqrt{\gamma_n}\ket{n}\bra{n}$, giving the master equation
\begin{equation}
\dot{\rho} = -\frac{i}{\hbar}[H(t),\rho(t)] + \Gamma\left[\text{diag}(\rho)-\rho\right],
\end{equation}
if $\gamma_0=\gamma_1$, affecting the coherence by reducing the off-diagnol elements of the density matrix. For relaxation, we consider the energy dissipation from the qubit to the external environment on the $X$ direction, i.e., $C=\sqrt{\gamma}\hat{\sigma_x}$. The non-unitary evolution governed by master equation with Lindblad terms converts the density matrix into mixed state,  in which the classical probability $p_n$ of being in $\ket{\psi_n}$ cannot be retrieved. Hence, the perturbed system after weak measurement cannot be analytically calculated by Eq.~\eqref{eq:perturbed}~\cite{PhysRevLett.124.140504}. Here we extend it to the case of the density matrix. Similarly, a Gaussian pointer of pure state $\rho_p=\ket{\Phi}\bra{\Phi}$ is coupled to a two-level system of mixed state $\rho$ by the interaction Hamiltonian $H_{int}=g\delta(t-t')p\otimes\hat{\sigma_z}$. The collective system is evolved from the initial state $\rho_{ini}=\rho_p\otimes\rho$ to $\rho_{fin}$ after the coupling by
\begin{equation}
\rho_{fin}=\exp(-igp\otimes\hat{\sigma_z})\rho_{ini}\exp(igp\otimes\hat{\sigma_z}),
\end{equation}
shifting the pointer by $\langle\hat{\sigma_z}\rangle=\text{Tr}(\rho\hat{\sigma_z})$ when $g=1$. One retrieves the wave pointer after the coupling by tracing out the qubit. The measurement of the pointer's position projects the pointer to its eigenstate $\ket{q_0}$, where the projection operator of the collective system reads $\hat{P}=\ket{q_0}\bra{q_0}\otimes{I}$. In this way, we have the qubit's density matrix after the weak value feedback of $q_0$ by the projection operator and tracing out the pointer.

After clarifying the calculation of state perturbation in the language of the density matrix, we formulate the specific task to be studied by stream learning. We consider the optimal control of a continuously measured qubit within operation time $T$ by ML algorithm. We aim at flipping it from $|0\rangle$ to $|1\rangle$ by a sequence of pulses on the $X$ direction. Each pulse lasts a small interval of $\Delta t$, being described by the driving Hamiltonian $H=\Omega\hat{\sigma_x}$, followed by a weak measurement on the $Z$ direction. We assume that the measurement process is impulsive, i.e., the coupling and projective measurement on the pointer are instant, being independent of the dynamical evolution. Meanwhile, the control pulses may also be imprecise, including slight detuning $H=\Omega\hat{\sigma_x}+\Delta\hat{\sigma_z}$ and amplitude error $\Omega\rightarrow\Omega+\delta\Omega$. The weak value and the last pulse amplitude are fed to the ML model as streaming data. Accordingly, instant feedback from the model controls the quantum system for the next timestep.

\begin{figure*}[h!]
	\centering
	\includegraphics[width=0.8\linewidth]{transfer.pdf}
	\caption{(a) The expectations on $X$ (green), $Y$ (red), and $Z$ (blue) directions of repetitively measured qubits once the agent in Fig.~\ref{fig:threemodels}(c) is employed for driving the qubits in the presence of hybrid detuning errors and quantum noises. The solid curves in Fig.~\ref{fig:threemodels}(c) are plotted in dashed curves as the benchmark. The average fidelity decreases from $0.9636$ to $0.8805$ with the standard deviation of $0.1874$. Imprecise control pulses outputted by the agent as step functions are plotted in different colors. (b) The expectations after training for about 2000 episodes. The solid curves in Fig.~\ref{fig:transfer}(a) are plotted in dashed curves for the benchmark. The agent fits the new environment, retrieving an average fidelity of $0.9513$ with the standard deviation of $0.0353$. Corrected control pulses outputted by the agent as step functions are plotted in different colors.}
	\label{fig:transfer}
\end{figure*}

\subsection{Numerical setup}
We apply the DRL method to our task for the stream learning approach. The environment includes a qubit, which is continuously measured and perturbed for weak values, controlled by the agent's pulses. An artificial neural network (ANN), as the agent, observes a state from the qubit, outputting an action for the control problem, which is trained by a deep learning algorithm for approximating the optimal policy function $\pi(a|s)$. The environment receives the action from the agent, evolving to the next timestep. It also computes the new RL state and the corresponding reward. It is worth noting that the environment in the quantum realm is different from other physical environments. The RL environment requires the numerical simulation if quantum information, e.g., density matrix elements or fidelity, is encoded in the RL state. Instead of a black box in other physical environments, one has to compute the density matrix based on the weak value and the control pulses since elements of the density matrix cannot be obtained from the qubit without destroying it.

In our practice, we set the tunable range of the Rabi frequency (pulse amplitude) as action $\Omega\in[0,3\pi]$ in dimensionless unit, which is renormalized to $\tilde{\Omega}\in[0,1]$ for fitting the neuron. Total operation time $T=1$ is uniformly separated by $n=100$ control pulses. Each pulse drives the qubit for a time interval of $\Delta t=0.01$. For saving the computational resources, the position space of the pointer is within $q\in[-50,50]$, being uniformly separated by $\Delta x=1$, which is large enough for $g=1$. Accordingly, the momentum operator $p$ is constructed by $[q,p]=i\hbar$ with boundary conditions. The density matrix of the collective quantum matrix has a size of $202\times202$. The coarse grained position space leads to weak values $q_0$ of integer number, which is renormalized to $\tilde{q_0}=(q_0+50)/100\in[0,1]$. Thereby, the state is defined as
$
s(t_i)=\{\tilde{\Omega}(t_{i-1}),\tilde{q}_0(t_i),i/n,|\rho_{11}(t_i)|,|\rho_{12}(t_i)|, |\rho_{21}(t_i)|,|\rho_{22}(t_i)|\}$,
including the last action as renormalized pulse, renormalized weak value, current system time, and elements of the density matrix. The RL state is observed by the agent, an ANN with three fully connected hidden layers of $64$ neurons activated by ReLU, evolving to the next state by the numerical simulation part of the environment, which receives an action from the agent.

\begin{figure*}[h!]
	\centering
	\includegraphics[width=1\linewidth]{threemodels.pdf}
	\caption{The expectations on $X$ (green), $Y$ (red), and $Z$ (blue) directions of repetitively measured qubits, which are driven by trained DRL agents under (a) detuning, (b) dephasing, and (c) $\hat{\sigma_x}$ relaxation, respectively. The dynamics of each qubit are plotted by shaded curves, being averaged for the solid curve. Fidelities are calculated by $F=\left(\text{Tr}\sqrt{\sqrt{\rho}\ket{1}\bra{1}\sqrt{\rho}}\right)^2$, leading to average fidelities of $0.9922$, $0.9919$, and $0.9636$ for each case, with the standard deviation of $0.0012$, $0.0012$, and $0.0361$. Control pulses provided by the agent as step functions are also plotted in different colors. We set the detuning strength $\Delta=0.05\Omega$, the dephasing rate $\Gamma=0.05$, and the relaxation rate $\gamma=0.05$. Other parameters are the same as those in previous figures, which are listed in the main text. Dashed curves, as baselines,  are derived by averaging the expectations of qubits under $\pi$-pulse control.}
	\label{fig:threemodels}
\end{figure*}

\subsection{Training of the agent and results}
We train three individual models for driving the qubit at the presence of detuning, dephasing, and relaxation on the $X$ direction, respectively. The agents approximate the optimal policy, which maximizes accumulated artificial rewards. We keep the generality in the design of reward functions since we have no specific preference for the shape of the pulses. For the $\ket{0}\rightarrow\ket{1}$ task, we reward the agent by $r(t_i)=|\rho_{22}(t_i)-1|$ per timestep as a negative value, aiming at a fast flipping operation. The agent gets an extra reward of $1000$ if $\rho_{22}$ meets the threshold $|\rho_{22}|>0.99$, and terminates the episode for calculating the total reward early. We also notice that punishment of $100$ if $|\rho_{11}|>0.05$ at the final timestep helps the convergence of the model.

Fig.~\ref{fig:threemodels} shows the high-fidelity closed-loop quantum control at different errors or noises. For relaxation on the $X$ direction, we change the terminal condition by $|\rho_{22}(t_i)|>0.99$ for four neighboring timesteps, which prevents the model from converging on trivial resonant $\pi$-pulses. We used the Proximal Policy Optimization (PPO) method~\cite{schulman2017proximal} to train the agent, with the learning rate being 1e-3 and a batch size of $20$. Other hyperparameters are the default setting by Tensorforce v0.5.3~\cite{tensorforce}. Moreover, we set a random error on the action, characterized by a centered Gaussian distribution with a standard deviation of $0.02$, equivalently implementing the time-varying systematic error in the quantum system. The models give control pulses that are robust against systematic errors. It is worth mentioning that a trade-off usually exists between fidelity and robustness. We obtain the models in Fig.~\ref{fig:threemodels} after about 2000, 3000, and 8000 episodes for controlling the system under detuning, dephasing, and relaxation on the $X$ direction, respectively.

\subsection{Transfer of the agent}
A ML model is online for service after being trained for a particular task, such as flipping a qubit under $\hat{\sigma_x}$ relaxation as Fig.~\ref{fig:threemodels}(c) does. One can evaluate the model by querying the information from the environment to check its validity after it is online. The flipped qubit can go for further tasks, which are independent of the model's duty. In this way, one can also evaluate it by checking the results of additional tasks without querying the environment and the model. Once the performance of a well-trained model deviates from its expectation, one concludes that the qubit in the environment changes, i.e., the quantum errors or noises are shifted to other patterns, if further tasks are reliable or not happening at all. Hence, we want another model for precise control in the new environment. Instead of discarding the current model and training a new one, which is apparently inefficient, we transfer the agent to a new environment for the same task, exploring its capability to fit the new environment without much effort.

We test the proposal by starting from the trained agent in Fig.~\ref{fig:threemodels}(c). By directly evaluating the agent in a new environment in the presence of detuning $\Delta=0.1\Omega$, dephasing rate $\Gamma=0.05$, and $\hat{\sigma_x}$ relaxation rate $\gamma=0.05$, the average final state deviates significantly from the previous result [c.f. Fig.~\ref{fig:transfer}(a)], resulting in a fidelity decrease as well. Furthermore, we train the agent with the same setting for about $2000$ episodes. By additional $20\%$ of the total episodes, the agent retrieves its performance before the environment shift happens [c.f. Fig.~\ref{fig:transfer}(b)].

\section{Discussion and Outlook}
Based on the numerical experiment in Sec.~\ref{sec:ne}, we reckon that DRL can be employed as a stream learning approach to closed-loop quantum control. The fidelity can be further improved by fine-tuning in another training environment, with different thresholds and the design of reward function. Interestingly, we found out that the optimized policy from the agent is interpretable to some extent, as Fig.~\ref{fig:threemodels}(c) and Fig.~\ref{fig:transfer} shown. The maximal tunable Rabi frequency is $3\pi$, which is $1.5$ times the $\pi$-pulse for an operation time of $T=1$. The agent drives the qubit with a relatively high frequency for reaching a large $\rho_{22}$ sooner. It is understandable since continuous measurement can be described in the language of superoperators, affecting the dynamics like quantum noises, which effectively can be suppressed by reduced operation time. Later, the pulse strength decreases significantly once $\rho_{22}$ is large enough, converging to a small constant value for more precise operations. Accordingly, it is the weak measurement predominant in the state evolution instead of the control pulse. It is also similar to the quantum Zeno effect that locks a wave function on its eigenstate by repeating projective measurements.

Now we further discuss this topic after analyzing the results above. In Sec.~\ref{sec:ne}, we explained that the RL environment consists of a qubit and a numerical simulation part. The qubit can be physical, e.g., constructed in superconducting circuits, trapped ions, photonics, etc., or simulated by classical computers as we performed in numerical experiments. Here we emphasize again that the numerical simulation for calculating the qubit dynamics is compulsory if we include quantum information $\rho_{ii}$ or fidelity in the RL state and reward. Although we can perform the weak measurement, extracting partial quantum information and converting it to weak value $q_0$ without destroying the quantum state, it is still impossible to retrieve the total information of the density matrix by a single shot of measurement. We cannot treat the qubit as a black box, as we usually do in other classical scenarios, where the observation of the RL state is instant and cost-less. By contrast, we have to calculate the qubit dynamics based on the actions and feedback, deducing $\rho_{ii}$ without operating on the qubit. It becomes a setback when one performs stream learning in the quantum realm since simulating the quantum dynamics is time-consuming, e.g., about $7$ seconds for an episode in our numerical experiment. However, stream learning for real quantum devices requires the simulation speedup of about $10^5$ times (compared to the T1 time of state-of-the-art superconducting qubit). A possible solution is to train another ANN to mimic the dynamics of the quantum system, with available information as input, outputting the quantum information to be deduced without measurements. The training of such ANN needs plenty of training data and adequate training methods, which goes beyond the scope of this work.

Another method to avoid the black-box problem is to exclude the quantum information in the RL state. The RL state may contain the weak value $q_0$ and other classical information such as last action, the system time, etc.  It should be reminded that the threshold criteria for early termination are no longer available in this paradigm. Thus, the training environment only rewards the agent by a constant at the end of each episode once a projective measurement on the target state succeeds. However, the agent struggles to learn the precise control since the signal-to-noise ratio of $q_0$ is too low. The reward criteria also need a large ensemble (batch size) to evaluate the fidelities of final  quantum states. In this way, the problem becomes more challenging, which can be applied for evaluating RL algorithms.

\section*{Acknowledgments}
This work is partially supported from NSFC (12075145), STCSM (2019SHZDZX01-ZX04 and 20DZ2290900), SMAMR (2021-40), Program for Eastern Scholar, QMiCS (820505) and OpenSuperQ (820363) of the EU Flagship on Quantum Technologies, EU FET Open Grant Quromorphic (828826) and EPIQUS (899368), Spanish Government PGC2018-095113-BI00 (MCIU/AEI/FEDER, UE), and Basque Government IT986-16. X.C. acknowledges the Ram\'on y Cajal program (RYC-2017-22482).

\bibliographystyle{ieeetran}


\section{Introduction}
Multiple object tracking with spawning events has become ubiquitous in recent years. However, due to the nature of observations that may or may not be generated from multiple time-varying objects, the association of the observations to the corresponding objects as well as tracking the time-varying number of objects are difficult and computationally demanding tasks. In this paper, we exploit the Bayesian nonparametric modeling over multiple objects as well as nonparametric modeling for observation association. Bayesian nonparametric frameworks have shown promising results in recent years. Due to the efficiency of Bayesian nonparametric modeling, in  \cite{ caron2012, moraffah2019random, moraffah2018,  moraffah2019PY}, authors employed nonparametric Bayesian modeling to improve the tracking of multiple objects. Hierarchical Bayesian nonparametric modeling has been exploited for switching in linear dynamic systems \cite{Fox2011} and for tracking an object with multiple observation sources \cite{moraffah2019-mdm, moraffah2020bayesian}.  In \cite{moraffah2019bayesian_thesis, moraffah2019inference-arxiv}, authors provide mathematical guarantees as well as Gibbs sampling formulations for the dependent nonparametric models. Moreover, Bayesian nonparametric frameworks are introduced for tracking objects in high clutter \cite{moraffah_clutter}. It is worth pointing out that even though all these frameworks provide auspicious results that are shown to outperform the previously introduced methods, they all assume each measurement is originated from one and only one object. This assumption introduces extensive challenges in applications that may violate this premise for which each measurement may or may not be originated from more than one object. For instance, the sudden creation of many objects due to fragmentation may potentially cause huge risks for tracking the objects. In this case, a simple tracking model cannot provide an accurate estimate of the object's trajectory.  In recent years, authors proposed solutions for spawning through the random finite set (RFS) theory \cite{bar1990tracking, Faber2016}. Zero-Inflated Poisson spawn model is introduced which is derived from cardinalized probability hypothesis density (CPHD) filter \cite{Bryant17}. This model only accounts for survival and birth. In \cite{Bryant18}, the generalized labeled multi-Bernoulli (GLMB)-based filter is used to address the birth and death issue posed by the CPHD-based technique. These models utilize Bayesian analysis on random sets and introduce Bayesian filtering techniques. However, It is shown that RFS-based solutions are computationally demanding. In this paper, we introduce Bayesian nonparametric techniques based on the Dependent Dirichlet process (DDP) for multiple object tracking and the Indian Buffet Process (IBP) to account for both tracking time-varying number of objects and spawning. The nature of the IBP prior allows each measurement to come from multiple objects. Our nonparametric spawning model accounts for the measurement association to multiple objects where each measurement may be originated from more than one object. Moreover, the number of objects that originate a measurement is assumed to be unknown.   In this paper, we propose a novel Bayesian nonparametric approach that not only enables the joint estimation of a spawned object's state as well as the object cardinality but also addresses the measurement to objects association.  We incorporate a dependent process to account for multiple object tracking introduced in \cite{moraffah2018} and the IBP-based technique to address spawning.

The rest of the paper is organized as follows. In section \ref{sec:background}, we provide a quick review of the Dependent Dirichlet process (DDP) and the Indian Buffet Process (IBP). Section \ref{sec:PF} and section \ref{sec:proposed_method} state the problem formulations and our proposed framework, respectively. In section \ref{sec:ex}, we then demonstrate through experiments that our proposed framework outperforms the GLMB-based spawning technique introduced in \cite{Bryant18}.

\textbf{Notation:} We denote the the $j$th observation and the $\ell$th object state at time $k$ by $\by_{j,k}$ and $\bx_{\ell,k}$, respectively. The set of all possible states at time $k$ is shown by $\mathbb{X}_k$. Parameters associated with states and observations are $\theta$ and $\phi$, respectively. We indicate the $j$th element of vector $V$ by $[V]_j$. We assume exchangeability throughout this paper.

\section{Background}
\label{sec:background}
\subsection{Dependent Dirichlet Process (DDP)}

In recent years, the use of Dirichlet Processes (DPs) for modeling complex data has been well-established. However, DPs may not capture the dependency in data. To address the dependency issue, MacEachern introduced dependent nonparametric Bayesian processes as  \cite{maceachern1999dependent}
\begin{eqnarray}
G_x(\cdot) = \sum_{i=1}^\infty \pi_{i,x}\delta_{\theta{i,x}}
\end{eqnarray}
where $\pi_{i,x} = v_{i,x}\prod_{j=1}^{i-1}(1-v_{j,x})$. In this setup,  $v_{i,\bX}$ is considered to be a stochastic process. It is worth emphasizing that unlike vanilla Dirichlet process, in this formulation both $v_{i,\bX}$ and $\theta_{i,\bX}$ are stochastic processes. Moraffah in \cite{moraffah2018} introduced a class of Dependent Dirichlet process prior that can be employed in tracking multiple time-varying objects. This formulation captures birth, death, and remaining of time-varying objects in the field of view. Taking advantage of the Chinese Restaurant Process (CRP) representation of DP, the $\text{DDP}(\alpha)$ prior is constructed analogously as follows:
\begin{itemize}
\item \textbf{Case 1:}  Assigning the $\ell$th object at time $k$  to one of the survived and transitioned to time $k$ which have been occupied by at least one of the first $(\ell - 1)$ w.p.
\begin{eqnarray}
\begin{split}
\Pi_1 (\text{$j$th cluster}|\theta_{1,k},\dots,&\theta_{\ell-1,k}) = \\
&\frac{[V_{k|k-1}]_j + [V_k]_j}{\alpha+\sum_i[V_{k|k-1}]_i + [V_k]_i }
\end{split}
\end{eqnarray}
where $V_{k|k-1}$ and $V_k$ are the vector of the number of objects in each cluster after transitioning to time $k$ and at time $k$, respectively.
\item \textbf{Case 2:} Assigning the $\ell$th object to one of the survived and transitioned cluster at time $k$ where it has not been assigned to any of the previous objects w.p.
\begin{eqnarray}
\begin{split}
\Pi_1 (\text{$j$th cluster}|\theta_{1,k},\dots,&\theta_{\ell-1,k}) = \\
&\frac{[V_{k|k-1}]_j }{\alpha+ \sum_i  [V_{k|k-1}]_i + [V_k]_i }
\end{split}
\end{eqnarray}
\item \textbf{Case 3:} Assigning a new cluster to the $\ell$ object w.p.
\begin{eqnarray}
\begin{split}
\Pi_1 (\text{new cluster}|\theta_{1,k},\dots,&\theta_{\ell-1,k}) = \\
&\frac{\alpha}{\alpha+ \sum_i[V_{k|k-1}]_i + [V_k]_i }.
\end{split}
\end{eqnarray}
\end{itemize}

\subsection{Indian Buffet Process (IBP)}
Indian Buffet Process (IBP) is a class of nonparametric models that can be considered as a prior on the sparse infinite binary matrices \cite{JMLR:v12:griffiths11a}. Unlike CRP which is the distribution on the partitions of data induced by DP-- implying each data point belongs to one and only one cluster-- the IBP provides distributions over infinite binary matrices where each data may comprise more than one feature. Using CRP analogy, the IBP prior is constructed as follows \cite{JMLR:v12:griffiths11a}:
\begin{itemize}
\item The first customer starts at the left of the buffet, and takes food from each dish, stopping after a Poisson($\gamma$) number of dishes.
\item The $j$th customer samples dishes proportional to their popularity w.p. $n_i/j$ and tries Poisson($\gamma/j$) new dishes, where $n_i$ is the number of previous customers who have tried $i$th dish.
\end{itemize}
The aforementioned procedure produces a customer-dish infinite binary matrix $A$ which is known as the feature matrix.

\section{Problem Formulation}
\label{sec:PF}
We consider time-varying multi-object tracking problem where each measurement may or may not be generated from the unknown number of objects. Additionally, not only is the total number of objects unknown but since it is unknown which objects generate each measurement, the number of objects that generate a specific measurement is also assumed to be unknown.

Let us consider a multiple object tracking model with time-varying object cardinality $N_k$ ($N_k$ is unknown) and time-varying observation cardinality $M_k$, at time step $k$. Let $\mathbb{X}_k = \{\bx_{1,k},\dots, \bx_{N_k, k}\}$ and $\mathbb{Y}_k = \{\by_{1,k},\dots, \by_{M_k,k}\}$ be the time varying collection of states and measurements  at time step $k$, respectively. We assume that each measurement may or may not be generated from more than one objects, however, the measurement-to-object association is not known in advance. Transitioning from time $(k-1)$ to $k$, objects may enter (new objects may appear), leave with probability $(1- \text{P}_{\cdot, k \mid k-1})$, or remain in the field of view with probability $\text{P}_{\cdot, k \mid k-1}$ and then transition according to the (probability) transition kernel $p_{{\underline{\theta}}}(\bx_{\cdot,k} | \bx_{\cdot,k-1})$. The objective is to
\begin{enumerate}
\item provide measurement-to-object association, and
\item jointly estimate each object's trajectory and object cardinality.
\end{enumerate}
\section{Tracking with spawning framework}
\label{sec:proposed_method}
 We utilize Bayesian nonparametric modeling to place a Dependent Dirichlet process mixture model (DDPMM) introduced in \cite{moraffah2018} on the object states and stated in Section \ref{sec:background}. It is worth noting that our problem setup gives rise to a Bayesian nonparametric setting, where an infinite number of objects may generate each measurement. In particular,  we employ an Indian Buffet Process (IBP) as a prior over the allocation matrix $A$. To do so, we define the infinite binary allocation matrix $A = [a_{i,j}]$ as follows

\begin{eqnarray}
a_{i,j} =
\left\{
	\begin{array}{ll}
		1  & \mbox{if } \text{$i$th measurement is originated} \\
		&\mbox{} \text{ from the $j$th object} \\
		0 & \mbox{ } \text{otherwise}.
	\end{array}
\right.
\end{eqnarray}
\normalsize
Note that we assume matrix $A$ is the infinite latent matrix that determines from which objects each measurement is originated as the number of objects generating each measurement is unknown. It is also worth mentioning that as all Bayesian nonparametric models suggest, only a finite number of elements in matrix $A$ is non-zero \cite{JMLR:v12:griffiths11a}. We use an Indian buffet process (IBP) as a prior over such infinite binary matrices. We then infer the object identities and, hence, successfully track the multiple objects using its corresponding measurements in a robust manner.

\noindent Our proposed model consists of two main parts:
\begin{enumerate}
\item[(A)] One is a dependent Dirichlet Process (DDP) component which creates a multi-state prior by learning multiple object clusters or labels over related information;
\item[(B)] an Indian buffet process (IBP) which accounts for the measurement-to-object association.
\end{enumerate}
Integrating (A) with (B), we define a framework that can track multiple time-varying objects using their associated measurements.

\subsection{Dependent Dirichlet Process Prior on State Transition  }

In this section, we exploit the properties of the Dependent Dirichlet process (DDP) introduced in \cite{moraffah2018, maceachern1999dependent}. Let $G$ be a probability random measure drawn from a DDP($\alpha$). Note that it is shown that DDP is discrete with probability one. Therefore, using DDP as a prior over states is restrictive. We thus employ the dependent Dirichlet mixture model. In particular, let $\mathbb{X}^{\ell-1}_{k}$ and $\Theta^*_{k|k-1}$ be the configuration at time $k$ up to the $\ell$th object and the set of unique parameters survived and evolved to time $k$, respectively, then the generative model is defined as follows:
\begin{eqnarray}
\label{eq: mainmodel1}
&G\sim \text{DDP}(\alpha)\notag\\
&\Theta_k | G \sim G\\
&{\bf{x_{\ell,k}}} | {\mathbb{X}}^{\ell-1}_{k}, {\mathbb{X}}^*_{k|k-1}, \Theta^*_{k|k-1}, \Theta_k\sim F\notag
\end{eqnarray}
where $F$ is a distribution whose probability density function depends on the physical model.

\subsection{Measurement-to-Object Association }
Measurements can be originated from multiple different objects. However, objects as well as  the number of objects from which a measurement is generated are not known. In this section, we propose a Bayesian nonparametric model to address the association problem. to this end, let $A_k$ be the association matrix where each entry of the matrix at time $k$, where $a^k_{i,j} = 1$ if the $j$th measurement is generated from the $i$th object and $a^k_{i,j} = 0$ otherwise. Note that matrix $A$ is a latent matrix. We define an Indian buffet process as a prior over $A$. It is worth noting that this prior similar to Chinese restaurant process imposes a \textit{"rich gets richer"} effect and thus the matrix is sparse even though it has infinite dimension. To generate the measurements, we define $W_k = \mathbb{X}_k \odot A_k$ where $\odot$ indicates the Hadamard product (element-wise product). The generative model is defined as follows
\begin{eqnarray}
\label{eq: mainmodel2}
&A_k\sim \text{IBP}(\gamma)\notag\\
&\phi_{j,k}| {G_m}\sim G\\
&\by_{j,k}| \phi_{j,k}, \mathbb{X}_k, A\sim R(\cdot|\phi_{j,k}, \{\bx_{\ell,k}\}_{\{\ell:a_{j,\ell}\neq 0\}}, W_k)\notag
\end{eqnarray}
where $G$, $F$, and $R$ are some distribution following the physical model. One should determine the hyperparameters $\alpha$ and $\gamma$, in order to accurately estimate the trajectory of each object. We take a Bayesian approach and use Gamma  distributions as priors over the hyperparameters.

\section{Inference}
In this section, we aim to estimate the latent variables as well as the object trajectory, i.e., the state variables at each time $k$, conditioned on the observations up to time $k$. It is worth mentioning we are only interested in the state estimation and the rest of parameters and latent variables are nuisance. Nonetheless, for the sake of completeness, we provide the full posterior distribution. Marginalizing out is simple due to IBP representation of the model. The posterior distribution of interest is $\mathbb{P}(\{{\bf{x_{\ell,k}}}\}_{k, \ell}, \Theta_k, \Phi_k, W_k |\{\mathbb{Y}_k\}_k, \alpha, \gamma)$, where $\Theta_k, \Phi_k$ are the set of parameters for the model at time $k$. In particular, we are interested in
\begin{multline}
\label{eq:posterior}
\mathbb{P}(\{{\bf{x_{\ell,k}}}\}_{k, \ell}, \Theta_k, \Phi_k, W_k |\{\mathbb{Y}_k\}_k, \alpha, \gamma) \propto \\
\mathbb{P}({\bf{x_{\ell,k}}}, \Theta_k, \Phi_k, A_k , \{\mathbb{Y}_k\}_k| \alpha, \gamma)\\
 = \mathbb{P}(A_k|\gamma)\prod_{k=1}^{K}\bigg[\prod_{\ell = 1}^{D_k}\mathbb{P}(\theta_{\ell,k}|\theta_{\ell, k-1}, \alpha)\\
 \prod_{j=1}^{M_k}\prod_{\ell=1}^{N_k}\mathbb{P}({\bf{x_{\ell,k}}}|{\bf{x_{\ell,k-1}}}, \theta_{\ell,k})\mathbb{P}(\phi_{j,k}) \\
 \big[\mathbb{P}(y_{j,k}|\phi_{j,k}, \{{\bf{x_{\ell,k}}}\})\big]^{\mathbbm{1}(\ell:a^k_{j,\ell}\neq 0)}\bigg]
\end{multline}
where $\mathbbm{1}(\cdot)$ is an indicator function; i.e., the indicator function is $1$ if the argument is true and otherwise it is zero. Note that all terms in the posterior equation (\ref{eq:posterior}) are defined through equations (\ref{eq: mainmodel1}) and (\ref{eq: mainmodel2}). Computing the exact posterior mean is impossible as there is no closed form solution for the integral. Thus, we take advantage of Markov chain Monte Carlo methods. In this paper, we employ the collapsed Gibbs sampling technique to sample from the posterior distribution.

\section{Experiment}
\label{sec:ex}
\begin{figure*}[th]
	\centering
	\hspace*{-.4in}
	\subfloat[ ]{\label{figa}
		\resizebox{2.6in}{!}{\includegraphics{figs/x-y.eps}}}
	%
	\hspace*{-.31in}
	\subfloat[ ]{\label{figb}
		\resizebox{2.6in}{!}{\includegraphics{figs/card.eps}}}
	%
	\hspace*{-.3in}
	\subfloat[ ]{ \label{figc}
		\resizebox{2.6in}{!}{\includegraphics{figs/ospa.eps}}}
	%
	\vspace*{-3mm}
	\caption{(a) True and estimated x-coordinate and y-coordinate using the proposed method for 4 objects; (b) Cardinality estimation through the proposed method; (c) OSPA metric for the location and cardinality for order $p=1$ and cut-off $c=100$. The results are by 1000 Monte Carlo runs.}
	\label{fig}
\end{figure*}
\subsection{Experiment I: Performance of Proposed Method}
\label{subsec:ex1}
The performance of the proposed framework is demonstrated through experiments.  To this end, we use the same model as \cite{Bryant18}.  In this model, number of objects varies due to birth and death of the objects. To this end, the state model is an autonomous linear Gaussian with state vector $\bx = [x, y, v_x, v_y]^T$ where the first two coordinates denote the $x$ and $y$- coordinate of the object position and the last two coordinates are the object velocity following:
\begin{eqnarray}
\bx_k = A \bx_{k-1} + \zeta_k, \hspace{1cm} \zeta_k\sim\mathcal{N}(0, \Sigma_s)
\end{eqnarray}
where
\[
A = \begin{bmatrix}
I_2 & \Delta I_2\\
0_2 &I_2
\end{bmatrix} \hspace{0.75cm}\text{and} \hspace{0.75cm}
\Sigma_s = \sigma_s^2\begin{bmatrix}
\frac{\Delta^4 }{4} I_2& \frac{\Delta^3 }{2} I_2\\
\frac{\Delta^3 }{2} I_2 & \Delta^2  I_2
\end{bmatrix}.
\]
We set $\Delta = 1$s, $\sigma_s = 1$m$\text{s}^{-2}$, $I_2$, and $0_2$ are $2\times 2$ identity and zero matrices, respectively. The observations are also assumed to follow a linear Gaussian equation:
\begin{eqnarray}
\by_k = B \by_{k-1} + \eta_k, \hspace{1cm} \eta_k\sim\mathcal{N}(0, \Sigma_o),
\end{eqnarray}
where
\[
B = \begin{bmatrix}
I_2 & 0
\end{bmatrix} \hspace{0.75cm}\text{and} \hspace{0.75cm}
\Sigma_s = \sigma_o^2 I_2, \hspace{0.2cm} \sigma_o^2 = 10 m. \]
For the sake of this simulation, we assume all the densities in our framework are too Gaussian. We also placed a Normal-inverse-Wishart on the parameters of Gaussian distribution as $\text{NIW}(\mu_0, 0, \nu, I)$ for $\mu_0 = 0.001$ and $\nu = 50$ and Gamma distribution priors  on our hyperparameters-- concentration parameter $\alpha$-- as Gamma$(\alpha;  1, 0.1)$. In  Figure (\ref{figa}), we display the $x$ and $y$- coordinate range as well as the estimation of them through the proposed method.

The OSPA metric of order $p=1$ with cut-off $c=100$ for the proposed method is depicted in Figure (\ref{figb}). We also successfully estimate the time-varying cardinality of objects which is shown in Figure (\ref{figc}). To estimate the trajectory we the associations we employ collapsed Gibbs sampling techniques and run 1000 Monte Carlo runs.

\subsection{Experiment II: Comparison to GLMB with Spawning}

In the previous section, we demonstrated the performance of our proposed method. In this section, we provide a comparison between our proposed work and the GLMB-based spawning technique. In order to provide a fair comparison, we used the same setup introduced in \cite{Bryant18}. This experiment is performed to demonstrate the advantage of our proposed framework over the GLMB-based technique. In Figure \ref{fig:traj_comp}, we show the $x$ and $y$ coordinates of the range as well as an estimation of them through the proposed method. Figure \ref{fig:card_comp}, however, compare the cardinality estimation through our proposed method versus the cardinality estimation using GLMB-based algorithms. In Figure \ref{fig:traj_comp}, the OSPA comparison between the two methods is provided. As depicted, under the same conditions, the proposed method outperforms the existing GLMB-based spawning technique. The OSPA metric for comparison is provided for order $p=1$ with cut-off $c=100$. The comparison between the Bayesian nonparametric modeling and the GLMB-based technique verifies the advantage of the former.

\begin{figure}[h]
    \centering
    \includegraphics[width=0.5\textwidth]{figs/traj_comp.eps}
    \caption{True and estimated object trajectory for 10 objects using our proposed method based on DDP.}
    \label{fig:traj_comp}
\end{figure}

\begin{figure}[h]
    \centering
    \includegraphics[width=0.5\textwidth]{figs/card_comp.eps}
    \caption{True and estimated object cardinality at each time step using DDP-bases and GLMB-based algorithms for 1000 Monte Carlo runs.}
    \label{fig:card_comp}
\end{figure}

\begin{figure}[h]
    \centering
    \includegraphics[width=0.5\textwidth]{figs/ospa_comp.eps}
    \caption{OSPA metric for both cardinality and the location for p = 1, c = 100.}
    \label{fig:traj_comp}
\end{figure}

\section{Conclusion}
In this paper, we introduced a class of robust Bayesian nonparametric models to address multi-object tracking with spawning. Our proposed model has two folds. The first fold exploits the Dependent Dirichlet process to account for the birth and death of objects and the second fold utilizes the Indian Buffet Process as a prior over the association matrix for which one can associate each measurement to the corresponding objects. Integrating with a Markov chain Monte Carlo sampling technique -- in this paper, we used collapsed Gibbs sampling-- to sample from the posterior distribution. This framework accurately and robustly provides objects' trajectory estimation, object cardinality as well as measurement-to-object association. We demonstrated the performance of our proposed framework through simulations and showed that this framework, under the same conditions,  may outperform the RFS-based techniques.

\bibliographystyle{IEEEtran}