\section{Introduction}

The dynamics of open quantum systems coupled to non-Markovian reservoirs is at the forefront of current interest, as it is relevant to practical as well as to fundamental problems. In that endeavor, one dimensional many-body systems such as quantum spin chains arise in many contexts, as they offer valuable and versatile ground for quantitative studies. Their relevance extends from quantum information theory to  condensed matter physics, owing to their versatility as basic resources for the implementation of solid-state devices for quantum computing and quantum communication tasks \cite{ref1}. Among these tasks, faithful quantum state transfer \cite{ref2,ref3,ref4,ref5,ref6,ref7,ref8,ref9} and long-distance entanglement \cite{ref10,ref11,ref12,ref13,ref14,ref15} have been investigated in detail for a variety of spin chain configurations,  over the last 20 years or so, with  research on these fields continuing at an undiminished pace.

In that general context, the dynamics of an arbitrary length spin XX chain driven by non-Markovian reservoirs at both ends was studied in a recent paper by two of us \cite{ref16}. Our contribution in that paper was the detailed quantitative study of reservoirs, with Lorentzian as well as Ohmic spectral densities. Through a formulation in terms of the wave-function, combined with Laplace transform, we were able to obtain results for chains of arbitrary length, i.e. number of sites, and study dynamical as well as state transfer aspects of the system. Guided by the behaviour of the excitation dynamics, we assessed the qualitative differences of the system from that driven by a Markovian reservoir which is known to entail monotonic decay. That assessment, indirect as it was, can only be considered qualitative, although it did provide a glimpse of the similarities between these two non-Markovian cases. Yet, there are a number of technical measures aiming at a quantitative assessment of the difference between the states of a quantum system, as it evolves under coupling to different external causes; in this case reservoirs. The magnitude of the difference between the evolution under a given reservoir from that under a Markovian one, should therefore provide a quantitative measure for what, following standard terminology \cite{ref17,ref18}, we shall refer to as “degree of non-Markovianity” (DNM). The term degree of Markovianity (DM) can also be found in the literature, meant to indicate the similarity  to the evolution under a Markovian reservoir. Obviously, a large  DNM corresponds to a small DM. We have adopted the term DNM in this paper, as it  appears to have prevailed in the literature. Since  the terms environment and reservoir appear in the literature as synonyms, we are using them interchangeably in this paper.

Having examined much of the existing literature on the issue, we have come to the realization that our scheme provides excellent territory for the quantitative study of non-Markovianity. The relevant literature is vast, ranging from general formal considerations \cite{ref19,ref20,ref21,ref22,ref23,ref24} to examples in small specific systems \cite{ref25,ref26,ref27,ref28}, including methods for the estimation the DNM using machine learning tools \cite{ref29}. The advantage of our system stems from the combination of a realistic arrangement of an XX chain of qubits, coupled to each other, with one end coupled to a reservoir. In our previous paper \cite{ref16}, we had considered a chain coupled to reservoirs at both ends. Although the present quantitative study could as easily address that arrangement, since our objective is the comparative study of non-Markovian reservoirs, we have chosen to consider the same system with only one end qubit coupled to a variety of reservoirs. We calculate first the time evolution of the system under each reservoir, one of which is Markovian. In the second step, we calculate the difference of the state of the system evolved under a given non-Markovian reservoir from the state evolved under the Markovian one. The time evolution of that difference is what we refer to as the DNM of the given reservoir.

We have chosen to address the issue in terms of perhaps the most direct notion, namely the quantum state distance (QSD). A variety of measures appropriate as tools for our task have been proposed over the years \cite{ref30,ref31}. Their objective is the quantification of the difference between two different quantum states of a system, which is of fundamental importance in quantum information processing \cite{ref32,ref33}, state transfer being a case in point \cite{ref34}. The most common quantity pertaining to that issue is the Fidelity (F), for which there are various expressions in the literature \cite{ref35,ref36,ref37}. It is meant to characterize the similarity or difference between the density matrices representing the two states. Obviously it is related to the notion of distance between two quantum states, which explains why the formulation of some distance measures involve one of the F expressions. There is no unique expression or value for the F, which is also the case for the QSD.  In some sense, it may be matter of taste or perhaps usefulness in a particular context which value of F or QSD is adopted.

Our stated objective then is the quantitative evaluation and calibration of non-Markovianity of several non-Markovian reservoirs by means of QSD. In view of the diversity of expressions for the QSD, inevitably an additional component has been injected into our task, namely, the comparative analysis of various QSD measures. Specifically, since the value of DNM does depend on the particular QSD measure, it is important to know whether and to what extend the classification of various reservoirs in terms of DNM is independent of the QSD measure employed. And this completes the road map of our present project.

In the next session, we review various QSD measures, their properties and the modifications necessary in order to accommodate issues pertaining to open systems. In section III we describe the theoretical approach to our problem in terms of the time-dependent Schr{\"o}dinger’s equations of motion in Laplace space, from which one can obtain even closed form solutions for the judiciously transformed amplitudes, for chains of arbitrary  length, as well as arbitrary initial conditions, within the single-excitation subspace. Using the amplitudes of the chain sites, we construct the density matrices necessary for the calculation of the various QSD measures. In section IV we present the results of our study while in section V we provide an overview as well as some concluding remarks related to our work.

\section{Brief summary and necessary modifications of QSD measures}

QSD measures need to be formulated in terms of density operators, so as to allow their applicability to mixed states. In books and papers, more often than not, expressions for F or QSD are assumed, either tacitly or explicitly, to involve Hermitian density operators of trace equal to unity. Coupling of an initially closed system to any reservoir, however, leads to decay. As a consequence, the reduced density operator of the system becomes non-Hermitian, while the trace decays usually to zero, as time tends to infinity. This implies that in expressions involving the density operator, its Hermitian adjoint, as well as the trace must be retained explicitly, even though typically at $t=0$ the trace is normalized to unity.

To illustrate the issue, let us consider the most common and straightforward QSD measure, referred to as the trace distance. It is defined by
\begin{equation}
D_T \equiv \frac{1}{2}  \Tr \abs{ρ-σ}
\end{equation}
where $\abs{A}$, for a matrix $A$, stands for $\sqrt{A^{\dagger}A}$ which stands for the positive definite square root of the matrix under the radical symbol, with $A^{\dagger}$ being the Hermitian adjoint of $A$.
Applying the above definition to the expression for the trace distance, we obtain

\begin{equation}
D_T = \frac{1}{2}  \Tr \sqrt{ \left( ρ^{\dagger} ρ + σ^{\dagger} σ - ρ^{\dagger} σ - σ^{\dagger} ρ \right)}
\label{D_T_extended}
\end{equation}

Now, if both density operators under the radical are Hermitian, the expression for the trace distance shown in Eqn. (\ref{D_T_extended}) reduces to

\begin{equation}
D_T = \frac{1}{2}  \Tr \left( ρ-σ \right)
\end{equation}
This expression may be valid for closed systems. Incidentally, if the two density operators are normalized to unity, this expression will always be zero, which seems like a tautology.  However, since for a system coupled to a reservoir, the density operators become non-Hermitian,  the correct expression for the trace distance is given by Eqn. (\ref{D_T_extended}).

Let us consider now another expression for QSD, namely the Hellinger measure $D_H$, defined as

\begin{equation}
D_H^2 \equiv \Tr \left( \sqrt{ρ} - \sqrt{σ} \right)^2
\end{equation}
which upon expanding the square becomes
\begin{equation}
D_H^2 = \Tr ρ + \Tr σ - 2\Tr \left( \sqrt{ρ} \sqrt{σ} \right)
\label{Hellinger_squared}
\end{equation}
where the invariance of the trace under cyclic permutation of the factors has been used. If the system is closed, with density operators normalized to unity, Eqn. (\ref{Hellinger_squared}) reduces to $D_H^2 = 2 \left[ 1- \Tr \left( \sqrt{ρ} \sqrt{σ} \right)  \right]$, leading to the expression
\begin{equation}
D_H = \sqrt{2} \left[ 1- \Tr \left( \sqrt{ρ} \sqrt{σ} \right)  \right]^{1/2}
\label{Hellinger_closed}
\end{equation}
for the Hellinger distance measure, which typically is the expression cited for this QSD measure. However, in the general case of open systems, following Eqn. (\ref{Hellinger_squared}), we have:

\begin{equation}
D_H = \left[  \Tr ρ + \Tr σ - 2\Tr \left( \sqrt{ρ} \sqrt{σ} \right) \right]^{1/2}
\label{Hellinger_open}
\end{equation}

When we are dealing with open systems, which implies coupling to reservoirs, the density operators tend to zero as $ t \rightarrow \infty$. In that limit, Eqn. (\ref{Hellinger_open}) correctly yields zero for the distance between the two states, as they both decay. On the other hand, Eqn. (\ref{Hellinger_closed}) would lead to the value $\sqrt{2}$ which on physical grounds, is at best problematic, as it does not account for the decay of the traces of the density operators.

There is one more measure that we employ in our calculations. It is known as the Bures distance measure, usually denoted by $D_B$. Since it is typically defined in terms of one of the expressions for F, we list below the three most common expressions, labelled for our convenience $F_1$, $F_2$  and $F_3$. They are:

\begin{subequations}
\begin{equation}
F_1 \left( ρ, σ \right) = \left( \Tr \sqrt{ \sqrt{ρ} σ \sqrt{ρ}}  \right)^2
\end{equation}

\begin{equation}
F_2 \left( ρ, σ \right) = \Tr \sqrt{\sqrt{ρ} σ \sqrt{ρ}} = \sqrt{ F_1 \left( ρ, σ \right)}
\label{F_2}
\end{equation}

\begin{equation}
F_3 \left( ρ, σ \right) = \Tr \left( ρ σ \right)
\end{equation}
\end{subequations}

The Bures QSD is usually defined as $D_B^2 \left( ρ, σ \right) \equiv 2 \left[ 1 - \sqrt{ F_1 \left( ρ, σ \right)} \right]$, which, in view of Eqn. (\ref{F_2}) above, can be simply written as

\begin{equation}
D_B^2 \left( ρ, σ \right) = 2 \left[ 1 -  F_2 \left( ρ, σ \right) \right]
\label{Bures_QSD}
\end{equation}
from which we obtain

\begin{equation}
D_B \left( ρ, σ \right) = \sqrt{2} \left[  1 - \Tr \sqrt{\sqrt{ρ} σ \sqrt{ρ}} \right]^{1/2}
\end{equation}
Again, this expression is valid as long as the traces of both density operators remain equal to one, because if in any distance measure we set $ρ=σ$, we should obtain zero. This condition is satisfied as long as  $\Tr ρ = \Tr σ =1$ for all times. In the presence of dissipation, however, the traces decay to zero, in which case for $ρ = σ$ we would obtain the nonphysical value $\sqrt{2}$. This contradiction is amended if the above expression is modified as shown in the following equation:

\begin{equation}
D_B \left( ρ, σ \right) = \sqrt{2} \left[ \frac{1}{2}\left(\Tr ρ+ \Tr σ \right) - \Tr \sqrt{\sqrt{ρ} σ \sqrt{ρ}} \right]^{1/2}
\label{Bures_open}
\end{equation}
This is the expression we are using for the Bures QSD in this paper. In the long time limit and in the presence of dissipation, the modified Bures distance tends to zero as it should. Part of our investigation addresses the rate with which the various QSD measures tend to their final value. We are in addition interested in the consistency of the values of the DNM obtained by different measures. For whatever it is worth, it seems to us that one could as well define a QSD measure through an extension of Eqn. (\ref{Bures_QSD}), by inserting any of the expressions for F.

\section{Theory}

\begin{figure*}[t]
	\centering
	\includegraphics[width=17.8cm]{Figure_1.pdf}
		\caption[Fig1]{Schematic representation of the system at study and our methodology. A chain of an arbitrary number of coupled qubits is interacting with an external reservoir at its one end. The system is studied for two different cases of external reservoirs, namely for a Markovian and a non-Markovian reservoir with spectral density $J(ω)$. Through our methodology we calculate the qubit amplitudes in the single-excitation subspace via which we construct the density matrix of the chain for each case. The QSD between the two density matrices is then calculated using various measures, allowing us to assess the degree of non-Markovianity of the system for different types of spectral densities $J(ω)$.}
		\label{Fig1}
\end{figure*}

As discussed in section II, all expressions for QSD measures involve the density matrices $ρ$ and $σ$ of the quantum states whose distance is to be calculated. In general, the temporal evolution of the density matrix of a system coupled to an external environment is given by the solution of the master equation, obtained after tracing out the environmental degrees of freedom from the density matrix of the compound system (system + environment) \cite{ref38}. Despite the ongoing progress in the derivation of closed form expressions for master equations of quantum systems coupled to non-Markovian reservoirs \cite{ref39,ref40,ref41,ref42,ref43}, such expressions are usually obtained for special cases of environmental spectral densities under certain approximations and often are too complicated to be handled even numerically. Consequently, an attempt to quantify the degree of non-Markovianity in systems coupled to environments with a variety of possible spectral densities using numerical solutions of the respective master equations would be pointless -if possible at all- for certain forms of environmental spectral densities, such as those studied in this work.

To address this difficulty, we have developed  a formulation enabling the calculation of the amplitudes of the wavefunction of the open system, from which we then construct its density matrix. This formulation has been described in detail in a previous paper \cite{ref16}, where it was used to calculate aspects of the dynamics of a XX chain boundary-driven by non-Markovian environments. As demonstrated below, one of the advantages of this approach is its applicability to any possible form of environmental spectral density, as well as an arbitrary number of qubits in the chain.

The system under consideration consists of $N$ identical qubits in a one-dimensional configuration coupled to each other via a nearest-neighbor coupling $\mathcal{J}$. The $N^{th}$ qubit of the chain is coupled, via a coupling strength $g$, to an external environment of an arbitrary spectral density $J(ω)$. Our goal is to assess the degree of non-Markovianity of that system by means of the QSD between its density matrix and the density matrix of the system consisting again by $N$ qubits but with the $N^{th}$ qubit coupled to a Markovian environment. The Markovian environment is known to induce a shift of the energy of the last qubit, as well as a decay \cite{ref38}. We neglect the shift, as it is of no significance in the issue of dissipation and denote the decay by $\gamma_M$. The problem is studied using several QSD measures, for various forms of non-Markovian spectral densities $J(ω)$. Hereafter, we denote the density matrix of the Markovian open system by $σ(t)$ and the density matrix of the non-Markovian open system by $ρ(t)$. A schematic presentation of the systems under study and our methodology is depicted in Fig. 1.

The Hamiltonian of the compound system $\mathcal{H} = \mathcal{H}_S + \mathcal{H}_E + \mathcal{H}_I $, consists of three parts; namely, the Hamiltonian of the XX chain $\mathcal{H}_S$, the Hamiltonian of the environment which is coupled to $N^{th}$ (last) qubit of the chain $\mathcal{H}_E$, and the chain-environment interaction Hamiltonian $\mathcal{H}_I$, given by the relations $(\hbar = 1)$:

\begin{subequations}
\begin{equation}
\begin{split}
\mathcal{H}_S   =  & \omega_e \sum_{i=1}^N \ket{e}_i \prescript{}{i}{\bra{e}}  + \omega_g \sum_{i=1}^N \ket{g}_i \prescript{}{i}{\bra{g}}  \\
& + \sum_{i=1}^{N-1} \frac{\mathcal{J}}{2} \left( \sigma_i^{+} \sigma_{i+1}^{-} + \sigma_i^{-} \sigma_{i+1}^{+} \right),
\end{split}
\end{equation}

\begin{equation}
\mathcal{H}_E = \sum_\lambda \omega_{\lambda} {a^{E}_{\lambda}}^{\dagger} a_{\lambda}^{E} ,
\end{equation}

\begin{equation}
\mathcal{H}_I =  \sum_\lambda g ( \omega_{\lambda} ) \left( a_{\lambda}^{E} \sigma_N^{+} + {a_{\lambda}^{E}}^{\dagger} \sigma_N^{-} \right),
\end{equation}
\end{subequations}
where $\omega_g$ and $\omega_e$ are, respectively, the energies of the ground and excited state of each spin (all spins are assumed identical),  $\omega_{\lambda}$  is the energy of the $\lambda$-mode  of the environment consisting of an infinite number of bosons, $a_{\lambda}^{E}$  and ${a^{E}_{\lambda}}^{\dagger}$, are the bosonic annihilation and creation operators of the environment,  and $\sigma_i^{+}= \ket{e}_i \prescript{}{i}{\bra{g}}$ and $\sigma_i^{-}= \ket{g}_i \prescript{}{i}{\bra{e}}$, $i=1, \dots, N$ are the qubit raising and lowering operators, respectively.

The wavefunction of the compound system in the single-excitation space is expressed as:

\begin{equation}
\ket{\Psi (t)} = \sum_{i=1}^N c_i(t) \ket{\psi_i}  + \sum_{\lambda} c_{\lambda}^{E}(t) \ket{\psi_{\lambda}^{E}},
\end{equation}
where,

\begin{subequations}
\begin{equation}
\ket{\psi_i} \equiv \ket{g}_1 \ket{g}_2 \dots \ket{g}_{i-1} \ket{e}_i \ket{g}_{i+1} \dots \ket{g}_N \ket{0}_{E},
\end{equation}
and
\begin{equation}
\ket{\psi_{\lambda}^{E}} \equiv \ket{g}_1 \ket{g}_2 \dots \ket{g}_N \ket{0 0 \dots 0 1_{\lambda} 0 \dots 0 0}_{E}.
\end{equation}
\end{subequations}

\subsection{Construction of $σ(t)$ (Markovian system)}
As already noted above, a Markovian environment (i.e. an environment with slowly varying spectral density) interacting with a two-level quantum system (qubit) within the Born approximation, induces a shift as well as a decay in the energy of the latter. Formally, this result is obtained by tracing  the environmental degrees of freedom out of the density matrix of the compound system, a procedure found in many quantum optics and open quantum systems textbooks \cite{ref38, ref44}.  In our problem, the Markovian environment will only affect  the energy of the last qubit ($N^{th}$ qubit) of the chain. Neglecting the shift and keeping the decay rate $\gamma_M$, the time dependent Schr\"{o}dinger equation (TDSE) for the chain amplitudes reads:

\begin{subequations}
\begin{equation}
\begin{split}
i \frac{dc_i(t)}{dt} = & \left[ \omega_e + (N-1) \omega_g \right] c_i(t) + \frac{\mathcal{J}}{2} \left[ c_{i-1} (t) + c_{i+1} (t) \right], \\
& \quad i = 1, \dots, N-1,
\end{split}
\label{Ampl_j_Markovian}
\end{equation}
\begin{equation}
i \frac{dc_N(t)}{dt}= \left[ \omega_e + (N-1) \omega_g -i γ_M \right] c_N(t) + \frac{\mathcal{J}}{2} c_{N-1} (t)
\end{equation}
\label{Markovian_Eqns_of_motion}
\end{subequations}
where, in order to compress the notation in our equations, we allowed $i$ to run from 1 to $N-1$ in Eqn. (\ref{Ampl_j_Markovian}), noting that $c_{i-1} (t)=0$ for $i=1$, because qubit $1$ is coupled only to qubit $2$. To simplify the set of Eqns. (\ref{Markovian_Eqns_of_motion}), it is useful to introduce the amplitude transformations $c_i(t) = e^{-i\left[  \omega_e + (N-1) \omega_g \right]t} \Tilde{c}_i(t)$, $i= 1, \dots, N-1$ and $c_N(t) = e^{-i\left[  \omega_e + (N-1) \omega_g -i \gamma_M \right]t} \Tilde{c}_N(t)$, in terms of which our equations become:

\begin{subequations}
\begin{equation}
\frac{d \Tilde{c}_i(t)}{dt} = -i \frac{\mathcal{J}}{2} \left[ \Tilde{c}_{i-1} (t) + \Tilde{c}_{i+1} (t) \right], \quad i = 1, \dots, N-1,
\label{Ampl_j_Markovian_tilde}
\end{equation}
\begin{equation}
\frac{d \Tilde{c}_N(t)}{dt}= -i \frac{\mathcal{J}}{2}  e^{\gamma_M t}  \Tilde{c}_{N-1} (t)
\end{equation}
\label{Markovian_Eqns_of_motion_Tilde}
\end{subequations}
where, as above, $\Tilde{c}_{i-1} (t)=0$ for $i=1$.  The above set of $N$ coupled linear differential equations can be readily solved numerically to yield the time dependence of the tilde amplitudes of the chain sites. The density matrix $\sigma(t)$ of the chain can then be constructed in terms of the site amplitudes as shown the following equation:

\begin{equation}
\sigma(t) = \sum_{i,j}^N c_i(t) c_j^*(t) \ket{\psi_i} \bra{\psi_j}
\end{equation}

\subsection{Construction of $ρ(t)$ (non-Markovian system)}
The construction of $ρ(t)$, i.e. the density matrix of the chain with the last qubit interacting with a non-Markovian reservoir, is a  subtler task because the memory effects of the reservoir make the equation of the last qubit amplitude more complicated. In two of our recent papers, we examined the dynamics of XX chain boundary driven by non-Markovian reservoirs \cite{ref16}, as well as the connection between exceptional points and the Quantum Zeno effect in a simpler quantum system consisting of two interacting qubits one of which is coupled to a non-Markovian reservoir \cite{ref45}. On the basis of the formulation developed in those papers, as applied to the present system,  we outline here the basic procedure to obtain the amplitudes of the chain sites. In order to discriminate the amplitudes from those of the Markovian damped chain, we denote them by $c'_i (t)$, $i = 1, \dots N $.

Adopting the amplitude transformations $c_i' (t) = e^{-i \left[ \omega_e + (N-1) \omega_g \right] t} {\Tilde{c}'}_i (t)$ , $i = 1, \dots N $ and $ {c'_{\lambda}}^{E}(t) = e^ {-i \left( N\omega_g + \omega_{\lambda} \right)t}  {{\Tilde{c}}_{\lambda}}^{'E}(t)$, the TDSE yields the following equations of motion:
\begin{subequations}
\begin{equation}
i \frac{ d \Tilde{c}'_i (t)}{dt} = \frac{\mathcal{J}}{2} \left[ \Tilde{c}'_{i-1} (t) + \Tilde{c}'_{i+1} (t) \right], \quad i = 1, \dots, N-1,
 \label{Amplitudej}
\end{equation}

\begin{equation}
i \frac{ d \Tilde{c}'_N (t)}{dt} = \frac{\mathcal{J}}{2} \Tilde{c}'_{N-1} (t) + \sum_{\lambda} g(\omega_{\lambda}) e^{-i \Delta_{\lambda} t}  \Tilde{c}_{\lambda}^{'E}(t),
 \label{AmplitudeN}
\end{equation}

\begin{equation}
i \frac{ d \Tilde{c}_{\lambda}^{'E}(t)}{dt} = g(\omega_{\lambda}) e^{i \Delta_{\lambda} t} \Tilde{c}'_N (t),
\label{Env1}
\end{equation}
\end{subequations}
where $\Delta_{\lambda} \equiv \omega_{\lambda} - (\omega_e - \omega_g) \equiv \omega_{\lambda} - \omega_{eg}$ and $\Tilde{c}'_{i-1} (t)=0$ for $i=1$.

Formal integration of Eqn. (\ref{Env1}), under the initial condition $\Tilde{c}_{\lambda}^{'E}(0)=0$,  and substitution into Eqn. (\ref{AmplitudeN}), yields:

\begin{equation}
\begin{split}
\frac{ d \Tilde{c}'_N (t)}{d t} = & -i \frac{\mathcal{J}}{2} \Tilde{c}'_{N-1} (t) \\
& - \int_0^t \sum_{\lambda} e^{- i \Delta_{\lambda} \left(  t - t' \right)} \left[  g(\omega_{\lambda}) \right]^2   \Tilde{c}'_N (t')  dt'.
\end{split}
\label{Amplitude_c_N_new}
\end{equation}

The summation over the environmental modes is at this point replaced by an integration which requires the specification of the environment's spectral density $J(\omega)$, according to the relation $\sum_{\lambda} \left[  g(\omega_{\lambda}) \right]^2 \rightarrow \int d \omega J (\omega)$. The resulting set of differential equations then becomes:

\begin{subequations}
\begin{equation}
\frac{ d \Tilde{c}'_i (t)}{dt} = -i \frac{\mathcal{J}}{2} \left[ \Tilde{c}'_{i-1} (t) + \Tilde{c}'_{i+1} (t) \right], \quad i = 1, \dots, N-1,
\end{equation}

\begin{equation}
\frac{ d \Tilde{c}'_N (t)}{d t}  =  -i  \frac{\mathcal{J}}{2} \Tilde{c}'_{N-1} (t) -  \int_0^t R (t-t') \Tilde{c}'_N (t') dt',
\label{Non-Markovian_Diff_Eqn_Nth_site}
\end{equation}
\label{Non-Markovian_Diff_Eqns}
\end{subequations}
where
\begin{equation}
R (t) \equiv \int_0^{\infty} J (\omega) e^{- i \left( \omega - \omega_{eg} \right)   t } d\omega.
\label{R(t)}
\end{equation}
The second term in Eqn. (\ref{Non-Markovian_Diff_Eqn_Nth_site}) reflects the possibility for the excitation to be transferred from the $N^{th}$ qubit to the non-Markovian reservoir and vice versa.

Taking the Laplace transform of Eqns. (\ref{Non-Markovian_Diff_Eqns}) and following the procedure presented in \cite{ref16} we find that the Laplace transform $F_1(s)$ of the tilde amplitude $\Tilde{c}'_1 (t)$ is given by the expression:
\begin{widetext}
\begin{equation}
F_1(s) = \frac{ik c_N (0) - k^2 \left[s + B (s) \right] A_1(s) c_{N-1} (0)  + (ik) \sum_{m=1}^{N-2} \Big\{ ik  \left[ s + B (s) \right] A_{N-m}(s) - A_{N-1-m}(s) \Big\} c_m (0)}{ik \left[ s+ B (s) \right] A_N (s) - \left[ 1+  s+ B(s)  \right] A_{N-1} (s)}.
\label{F1}
\end{equation}
where $k \equiv \frac{2}{\mathcal{J}}$, $B(s)$ is the Laplace transform of the function $R(t)$ and $A_m (s)$ is given by the following expression:
\begin{equation}
A_m(s) \equiv \frac{\left[ (iks) +  i \sqrt{k^2 s^2 +4} \right]^m - \left[ (iks) - i \sqrt{k^2 s^2 +4} \right]^m}{2^m i \sqrt{k^2 s^2 +4}},\quad m = 1, \dots, N.
\label{A_m}
\end{equation}
\end{widetext}
The Laplace transforms $F_i(s)$ of the remaining qubit amplitudes are connected to $F_1(s)$ via the relation:
\begin{equation}
\begin{split}
F_i (s) & = A_i(s) F_1 (s) - (ik) \sum_{n=1}^{i-1} A_{i-n}(s) c_n (0), \\
& \quad i = 2, \dots, N,
\label{F_i-F_1}
\end{split}
\end{equation}
The effectiveness of this method rests on the fact that it allows the derivation of closed form expressions for the Laplace transforms of the tilde amplitudes for an arbitrary number of qubits $N$, an arbitrary spectral density $J(\omega)$ for the environment that interacts with the $N^{th}$ qubit of the chain, as well as arbitrary initial conditions. The Laplace inversion integrals can be readily calculated numerically (or even analytically in some cases) upon specification of these parameters to yield the tilde amplitudes of the chain sites, in terms of which we can calculate $c_i' (t)$ , $i = 1, \dots N $, and thus construct the density matrix for the chain via the relation:

\begin{equation}
\rho(t) = \sum_{i,j}^N c'_i(t) {c_j'}^{*}(t) \ket{\psi_i} \bra{\psi_j}
\end{equation}

Specific types of non-Markovian reservoirs, such as Lorentzian \cite{ref46,ref47}, Lorentzian squared \cite{ref48} and Ohmic \cite{ref49,ref50} are employed in order to assess the non-Markovianity of the chain in each case. The analytical derivations of the functions $R(t)$ and $B(s)$ for these types of reservoirs are given in Appendix A.

\section{Results $\&$ Discussion}

\begin{figure}[t]
	\centering
	\includegraphics[width=8.6cm]{Figure_2}
		\caption[Fig2]{Dynamics of the excitation survival probability of a single qubit coupled to various types of reservoirs: (a) Markovian reservoir with $γ_M =0.01$, (b) Lorentzian reservoir with $g=1$, $γ=0.03$ and $Δ_c=0$, (c) Lorentzian squared reservoir with $g=1$, $γ=0.3$ and $Δ_c=0$, (d) Ohmic reservoir with $g=1$, $S=1.5$, $ω_c=8$ and qubit frequency $ω_{eg}=10$.  }
		\label{Fig2}
\end{figure}

\begin{figure*}[t]
	\centering
	\includegraphics[width=17.6cm]{Figure_3}
		\caption[Fig3]{QSD measures between the Markovian-damped and the non-Markovian-damped system using various types of non-Markovian reservoirs in the single-qubit case. The types of reservoirs chosen and their parameters are: Markovian reservoir with $γ_M =0.01$,  Lorentzian reservoir with $g=1$, $γ=0.03$ and $Δ_c=0$, Lorentzian squared reservoir with $g=1$, $γ=0.3$ and $Δ_c=0$, Ohmic reservoir with $g=1$, $S=1.5$, $ω_c=8$ and qubit frequency $ω_{eg}=10$. (a) Markovian-Lorentzian Trace Distance, (b) Markovian-Lorentzian-squared Trace Distance, (c) Markovian-Ohmic Trace Distance, (d) Markovian-Lorentzian Hellinger Distance, (e) Markovian-Lorentzian-squared Hellinger Distance, (f) Markovian-Ohmic Hellinger Distance, (g) Markovian-Lorentzian Bures Distance, (h) Markovian-Lorentzian-squared Bures Distance, (i) Markovian-Ohmic Bures Distance.}
		\label{Fig3}
\end{figure*}

The stated objective of this paper is the quantitative assessment of the DNM of non-Markovian reservoirs in a realistic context provided by a spin chain. As outlined in the previous sections, for the task we have chosen to employ measures of QSD between the state of the chain coupled to various reservoirs. An obviously suitable for the purpose quantity is the QSD between the density operator of the chain evolved under a Markovian from that under a non-Markovian reservoir. From a physical standpoint,  it stands to reason that the larger the values of the QSD, irrespective of the measure under consideration, the less Markovian the system's evolution will be, which means higher DNM. Given that reservoirs by definition involve certain parameters, characterising their properties, in addition to the system-reservoir coupling constant, a decision has to be made as to the choice of parameters that make the value of QSD most meaningful. It is rather obvious that a direct comparison between, for example, an evolution under a Markovian reservoir with damping rate $\gamma_M$ and a Lorentzian one with arbitrarily chosen parameters $g$, $\gamma$ and $\Delta_c$ would be pointless. We suggest that, for the comparison to be meaningful, a feature common to the effect of all reservoirs entering the comparison should be adopted. With that in mind, we note that for all reservoirs the excitation of all qubits will eventually go to zero. Since initially the excitation is in the first qubit, without claim to uniqueness, we have chosen the parameters of the various reservoirs such that the half life of the excitation of the first qubit is approximately the same to its half life under a Markovian reservoir with a given $\gamma_M$.

\subsection{Single qubit}
We begin our analysis by considering the simplest case of a single qubit coupled to a variety of non-Markovian reservoirs. These  results serve as a point of reference in the comparative analysis of the interplay between non-Markovianity  and qubit-qubit interaction, as reflected in the behaviour of the QSD measures as a function of the time, analysed in section IVB.

In Fig. 2 we plot the time dynamics of the single qubit excitation survival probability for various types of reservoirs. In particular, we examine the cases of Markovian, Lorentzian, Lorentzian squared and Ohmic reservoirs, with the parameters of each non-Markovian reservoir chosen as described above. The Lorentzian and Lorentzian squared spectral densities result in general to fast oscillations in the single-qubit survival of excitation dynamics, indicating the exchange of the excitation between the qubit and the environment within finite times. The Ohmic spectral distribution tends in general to be much broader than the Lorentzian and Lorentzian squared distributions, leading to a qubit excitation survival probability that resembles more the respective survival probability of the excitation of the Markovian-damped qubit.

The QSD measures between the Markovian-damped single-qubit system and the non-Markovian ones are presented in Fig. 3. In all three of the considered non-Markovian reservoirs, the coupling strength $g$ between the single-qubit and the reservoir has been kept the same. As expected, all QSD measures are zero for $t=0$, as well as in the long time limit, as  there is practically no excitation left in the open system.

We focus first on the comparison  of the degrees of non-Markovianity, for different types of non-Markovian reservoirs, resulting from the same QSD measure. The trace distance between the Markovian and the Lorentzian damped system (Fig. 3(a)) exhibits rapid oscillations between two bounds, whose values change over time. The upper bound is initially at its maximum value and over time decreases to zero, whereas the lower bound initially increases, reaches a maximum, after which it decreases merging with the upper bound,  tending eventually to zero for long times. The  dynamics of the trace distance between the Markovian and the Lorentzian-squared damped system (Fig. 3(b)) also exhibit rapid oscillations, following a trend similar to the Markovian-Lorentzian case, with the exception that the lower bound can be zero within finite times, indicating that, at those instants, the Lorentzian-squared damped system resembles the Markovian-damped one. The oscillations in both Figs. 3(a) and 3(b) are indicative of the fast excitation exchange between the qubit and the non-Markovian reservoirs. In Fig. 3(c) we plot the trace distance between the Markovian-damped and the Ohmic-damped qubit state as a function of  time. Clearly, the behaviour of the trace distance dynamics in this case differs substantially from the dynamics of the Markovian-Lorentzian and the Markovian-Lorentzian-squared cases, in that it does not exhibit any oscillations, reaches a maximum after some finite time before decreasing, with the overall values of the trace distance remaining much smaller than the respective values in panels (a) and (b) of Fig. 3. We may therefore safely conclude that, under the prescribed conditions on parameters, the Ohmic-damped single qubit system displays the lowest DNM.

In addition to the trace distance, we have calculated QSD's in terms of two other relevant measures, namely the Hellinger and Bures measures, with the results shown in Fig. 3. The overall behavior bears significant similarity to that of the trace distance, with some non-negligible quantitative differences. The oscillations reflecting the exchange of excitation between qubit and reservoir are still present in both measures for the Markovian-Lorentzian and Markovian-Lorentzian-squared cases, whereas they are absent for the Markovian-Ohmic one. However both of these measures produce generally slower decrease of the QSD from the respective trace distance measure. For example, in terms of the trace distance, the value of the QSD between Lorentzian and Markovian in the vicinity of $t = 100$ hovers around $0.15$, whereas the Hellinger and Bures measures indicate values around a mean of $0.4$. Moreover, the oscillations of the QSD for the Markovian obtained from these two measures are damped much more slowly than the oscillations in the trace distance. The Ohmic case seems to follow the same trend, in that these two measures produce larger values for the QSD and a much slower evolution towards the expected value of $0$. In summary, for all three non-Markovian reservoirs in the single qubit case, the trace distance is found to show faster decay of non-Markovianity than the Hellinger and Bures measures do. It is worth noting that, despite the non-trivial formal difference of the expressions for the Hellinger and Bures QSD (Eqns. (\ref{Hellinger_open}) and (\ref{Bures_open})), they lead to practically the same dynamics, indicating the same DNM.

\subsection{Five mutually interacting qubits}

\begin{figure}[H]
	\centering
	\includegraphics[width=8.6cm]{Figure_4}
		\caption[Fig4]{Dynamics of the excitation survival probability of the first qubit of a chain consisting of $N=5$ qubits with qubit-qubit coupling $\mathcal{J}=1$. The initial excitation is chosen to be on the first qubit while the last qubit of the chain is coupled to various types of reservoirs: (a) Markovian reservoir with $γ_M =0.01$, (b) Lorentzian reservoir with $g=1$, $γ=0.03$ and $Δ_c=0$, (c) Lorentzian squared reservoir with $g=1$, $γ=0.3$ and $Δ_c=0$, (d) Ohmic reservoir with $g=1$, $S=1.5$, $ω_c=8$ and qubit frequency $ω_{eg}=10$.}
		\label{Fig4}
\end{figure}

\begin{figure*}[t]
	\centering
	\includegraphics[width=17.6cm]{Figure_5}
		\caption[Fig5]{QSD measures between the Markovian-damped and the non-Markovian-damped system using various types of non-Markovian reservoirs for the case of $N=5$ qubits with qubit-qubit coupling $\mathcal{J}=1$. The types of reservoirs chosen and their parameters are: Markovian reservoir with $γ_M =0.01$,  Lorentzian reservoir with $g=1$, $γ=0.03$ and $Δ_c=0$, Lorentzian squared reservoir with $g=1$, $γ=0.3$ and $Δ_c=0$, Ohmic reservoir with $g=1$, $S=1.5$, $ω_c=8$ and qubit frequency $ω_{eg}=10$. (a) Markovian-Lorentzian Trace Distance, (b) Markovian-Lorentzian-squared Trace Distance, (c) Markovian-Ohmic Trace Distance, (d) Markovian-Lorentzian Hellinger Distance, (e) Markovian-Lorentzian-squared Hellinger Distance, (f) Markovian-Ohmic Hellinger Distance, (g) Markovian-Lorentzian Bures Distance, (h) Markovian-Lorentzian-squared Bures Distance, (i) Markovian-Ohmic Bures Distance. In the inset of panel (f) we show the dynamics of the probability of finding the excitation in the Ohmic environment for the case of $N=5$ qubits. The time window where the probability dynamics abruptly become oscillatory coincides with the time window where the Markovian-Ohmic QSD measures follow the same behaviour.}
		\label{Fig5}
\end{figure*}

The merit of our approach rests with the possibility to obtain straightforwardly the Laplace transforms of the qubits tilde state amplitudes for an arbitrary number of qubits $N$, because $N$ appears in our equations as a modifiable parameter. Having calibrated the measures of QSD in characterising non-Markovianity in the simplest set-up of a single qubit, in this section we explore the role of non-Markovianity in the more realistic situation of a chain involving qubit-qubit interaction. For a concrete quantitative analysis, we have chosen five qubits in an XX chain with nearest neighbor coupling $\mathcal{J}$. The calculation follows the formulation presented in section III, with the initial excitation in qubit $1$ and the last qubit $5$ coupled to the reservoir. The question now is how the non-Markovianity of the reservoir affects the exchange of excitation between  the qubits, as manifested in the evolution of QSD between the Markovian and the various non-Markovian reservoirs.

In a set-up of this type, the crucial parameter is $\mathcal{J}$ as it controls the communication between qubits, which typically is the dominant aspect in applications. In order to remain within a realistic scenario, we have chosen the Markovian damping rate $\gamma_M$ much smaller than $\mathcal{J}$, keeping the choice of the parameters of the non-Markovian reservoirs as in the previous section, so that the initial excitation of the first qubit decays on a time scale approximately the same for all reservoirs. At first,  we calculate the time evolution of the excitation probability of the first qubit for all of the reservoirs under consideration,  with the results shown in Fig. 4. The glaring difference between the single-qubit dynamics in Fig. 2 and the $N=5$ dynamics of the first qubit excitation in  Fig. 4, is the appearance of oscillations in the dynamics of the latter, indicating the exchange of population between the qubits of the chain. In the cases of Lorentzian and Lorentzian squared  reservoirs, such oscillations are superimposed on the oscillations between the last qubit and the reservoir, resulting to the characteristic oscillatory dynamics of $P_1(t)$ depicted in Figs. 4(b) and 4(c). For the Markovian case  (Fig. 4(a)) the oscillations are due solely to the coupling between the qubits, because as we have seen in the case of the single qubit (Fig. 2(a)) the Markovian damping is monotonic. Again, for the Ohmic-damped case, the dynamics of the first qubit population (Fig. 4(d)) bears strong similarity to the Markovian-damped chain, with a rather unexpected exception, in that the oscillations are found to collapse abruptly after a finite time, a behaviour that has been also reported in previous works involving qubits coupled to Ohmic reservoirs \cite{ref16,ref51}.

We consider now the dynamics of the QSD measures exploring the non-Markovianity of the three reservoirs in the context of a chain of $N=5$ interacting qubits, with the results shown in Fig. 5. First we note that the time scales of decay for all QSD measures are longer than the respective time scales in the single-qubit case of Fig. 3. This behaviour was to be expected, since the excitation is now spread over the whole chain and naturally it takes more time to be lost in the reservoir. At the same time, the values of the QSD measures are in general larger than the respective values of Fig. 3, indicating that the DNM of the open system is affected by the number of qubits in the chain. The Markovian-Lorentzian and Markovian-Lorentzian-squared QSD measures are seen to display an overall similar behaviour. We should, however,  keep in mind that the values of the parameters for the Lorentzian and Lorentzian-squared spectral densities have been  chosen according to the criterion explained in the discussion of Fig. 1. In contrast, an open system damped by a Lorentzian squared environment would display a much larger DNM compared to the same system damped by a Lorentzian spectral environment characterized by the same parameters ($g$, $γ$ and $Δ_c$). In fact, under those conditions, the excitation throughout the chain in the Lorentzian squared case remains equal to one for a very long time such that it resembles a practically steady state. All QSD measures we examined (Trace distance, Hellinger distance and Bures distance) exhibit dynamics following the same trend, although the exact distance values may differ in general between different measures. As clearly seen in Fig. 5, the Hellinger and Bures distances display minor differences in their dynamical behaviour for all of the reservoirs studied.

A rather interesting phenomenon appears in the dynamical behaviour of the QSD measures in the case of the Ohmic-damped chain. In contrast to the single-qubit case, the interplay between qubit-qubit and qubit- Ohmic reservoir interactions gives rise to an intricate behavior in the dynamics for all of the Markovian-Ohmic QSD measures considered. The dynamics are initially smooth for all QSD measures with the distance being zero at $t=0$ and slowly increasing afterwards. However, there is a time when they suddenly change and exhibit rapid oscillations. These rapid oscillations appear only in a time window which is the same for all QSD measures considered. Beyond that time window, the dynamics are smoothed out with the distance slowly decaying back to zero. The time window of that sudden modification of the behaviour of the QSD measures coincides with the time window within which the total probability of finding the excitation in the Ohmic reservoir, exhibits similar behaviour (see inset of Fig. 5(f)). Changes in the dynamical behaviour of the system can be also observed through abrupt changes in the dynamics of the qubit excitations of the chain, some of which were reported in one of our previous papers \cite{ref16}. In summary, the qubit-qubit interaction is found to enhance the DNM obtained in all QSD measures. Conversely, this may be interpreted as an enhancement of the excitation storage in the chain due to the non-Markovian character of the reservoir.

\section{Summary and Outlook}

The work in this paper was undertaken with the dual purpose of evaluating the DNM of some typical non-Markovian reservoirs, as well as the generalization of the most common QSD measures employed for the purpose. In order to offer an analysis connected with a realistic system typical in quantum information processing, we chose an XX chain of interacting qubits with a reservoir connected to one end of the chain. QSD measures require the density matrices of the systems under consideration. Since non-Markovian reservoirs do not in general lend themselves to descriptions in terms of a master equation, we employed an approach developed in earlier work \cite{ref16} based on the solution of the Schr{\"o}dinger equation, in which the amplitudes of the excitation of the sites as a function of time are obtained by means of Laplace transform. With the amplitudes in hand we construct the density matrix of the system which is necessary for the calculation of the desired QSD measure.

The expressions for QSD measures found in most of the literature explicitly or implicitly assume density matrices of  trace unity. However the density matrix of an open system by definition will eventually decay to zero, which means the trace cannot remain equal to unity. That led us to a generalization of the expressions of certain QSD measures so that they are valid for open as well as closed systems.

For the DNM we have chosen the distance of the quantum state of our system as evolved under a given non-Markovian reservoir from the state evolved under a Markovian one. A meaningful comparison required the judicious choice of the parameters entering the expressions of the spectral densities of the reservoirs involved in each QSD calculation. Given the flexibility of our approach as to the number of qubits it can handle, we studied first the case of a single-qubit coupled to a reservoir, so as to have a frame for the evaluation of the interplay between qubit-qubit interaction and qubit-reservoir. The results from the single qubit analysis showed that the Hellinger and Bures measures indicated slightly higher DNM for all three reservoirs than the trace distance measure did. On the other hand, all three measures produced by far much lower DNM for the Ohmic reservoir compared to the Lorentzian or Lorentzian squared. Moreover, the absence of any noticeable back-flow of excitation from the Ohmic reservoir to the qubit suggests that the Ohmic reservoir behaves practically in a Markovian fashion.

The Lorentzian and Lorentzian squared reservoirs, on the basis of all three measures were found to exhibit essentially the same DNM. This might seem counter intuitive, because inspection of the relevant formal expressions in the appendix suggests that the squared Lorentzian is more peaked than the Lorentzian with the same parameters, which would imply higher DNM. However, the Lorentzian squared employed in all of our comparative calculations (Figs. 2-5) was not the square of the accompanying Lorentzian, owing to the adjustment of the parameters explained in section III. If we calculated evolution of the qubit excitation dynamics and DNM for a Lorentzian and a Lorentzian squared spectral density characterized by the same parameters ($g$, $γ$ and $Δ_c$), we would indeed find more back-flow and higher DNM for the Lorentzian squared, which brings up an interesting issue in reservoir engineering. If in any open system subjected to dissipation through a Lorentzian reservoir, it were possible to modify the reservoir to the square of the Lorentzian, dissipation would slow down significantly. Our consideration of the Lorentzian squared spectral density is not a mere mathematical exercise. It was inspired by the observation of non-Lorentzian line shapes \cite{ref48,ref52}, including the Lorentzian squared, albeit in different physical context.

The results for the case of five interacting qubits chain bear similarities to as well as differences from the single-qubit case. The similarity is seen in the DNM by the three QSD measures which produce more or less similar values, with the value of the trace distance somewhat lower than the values of the Hellinger and Bures measures. In other words, the three measures provide mutually compatible results. However, the values of the DNM for all three reservoirs for the chain are noticeably higher than those for the single-qubit. Moreover, for the Ohmic reservoir, the QSD obtained by all three measures is not only quantitatively but even qualitatively different from that of the single qubit. The smooth dynamics obtained in the single qubit case, is in this case interrupted by a region of abrupt oscillations, within a certain time window in which the dynamics of the probability of finding the excitation in the Ohmic environment was found to exhibit the same behaviour.
However, the stunning difference between the single and five qubit cases lies in the retention of population distributed among the interacting qubits and the much slower decay of the excitation in all of the reservoirs in the latter case, which is compatible with the consistently higher DNMs reported in our work. It can then be argued that the exchange of excitation between the qubits slows down the dissipation which implies that the larger the number of the qubit the slower the loss into the reservoir. Perhaps the most significant message conveyed by the comparison of the single-qubit and chain of interacting qubits results is that the DNM is not independent of the quantum system in the context of which it is evaluated.

In closing, we believe that the appropriate modifications of the QSD measures, that extend their use to non-trace-preserving processes, which are inevitable in open systems, and the detailed analysis of their applicability to a qubit chain system, have revealed subtle aspects of non-Markovianity. Extension of our approach and analysis to other open systems can be expected to offer new insights into the performance of protocols pertaining to quantum information technology tasks as well as to the statistical mechanics of spin chains. Even for the single chain coupled to non-Markovian reservoirs, issues such as the dynamics of entanglement, and the DNM beyond single excitation represent potentially fruitful territory for exploration.

\section*{Acknowledgments}
GM would like to acknowledge the Hellenic Foundation for Research and Innovation (HFRI) for financially supporting this work under the 3rd Call for HFRI PhD Fellowships (Fellowship Number: 5525).

\begin{figure}[H]
	\centering
	\includegraphics[width=5cm]{HFRI_Logo}
\end{figure}


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

\textit{Introduction.---}Ever since the establishment of quantum mechanics a century ago, the collapse of the wave function has been its central topic much debated. Numerous interpretations and theories have been proposed, including the Copenhagen interpretation \cite{bohr1928quantum, neumann1955mathematical}, many-worlds interpretation \cite{RevModPhys.29.454}, decoherence theories \cite{adler2003decoherence}, de Broglie--Bohm theory \cite{PhysRev.85.166, PhysRev.85.180}, quantum Bayesianism \cite{PhysRevA.65.022305, fuchs2010qbism, RevModPhys.85.1693}, and many others \cite{RevModPhys.85.471}. It is intriguing that if there exists a theory that is fully compatible with continuous quantum mechanics yet is capable of describing the dynamics of wave function collapse. After all, during the measurement process, the atoms and molecules in the measuring devices and the probed quantum objects are both governed by the Schr\"{o}dinger equation, hence their interaction and dynamics during the collapse of the wave function should also in principle be governed by the Schr\"{o}dinger equation.

In this Letter, based on the recent work \cite{hu2022quantum} of ours, we are going to show, starting with pure state cases, how the continuous collapse of the wave function process can actually be described by the Schr\"{o}dinger equation with a stochastic time-dependent Hamiltonian. We then solve analytically for the state-dependent Hamiltonian responsible for measurements on an arbitrary $n$-level system and position measurements on an harmonic oscillator in the ground state, and propose a few experimental schemes to verify and utilize the conclusions. To address the uniqueness of the Hamiltonian, we discuss how one can find all possible Hamiltonians that may be responsible for the wave function collapse such that the true Hamiltonian must be one of them. We then discuss and show how the above formalism can be used also to describe the wave function collapse for mixed quantum states. Our results show that the wave function collapse could be fully describable by a modified Schr\"{o}dinger equation where the Hamiltonian is state-dependent, which gives a new approach to applying control to measured systems.

\textit{Wave function collapse and continuous quantum measurement.}---While the projective measurement is usually treated for convenience to be happened instantaneously, it is reasonable to expect that actual physical processes always take some finite time to complete \cite{minev2019catch, mandelstam1945uncertainty, margolus1998maximum, deffner2013energy}. If we continuously measure a quantum object such that the density matrix $\hat{\rho}(t)$ of the object change gradually rather than all in a sudden, then this constitutes a continuous quantum measurement \cite{caves1987quantum, diosi1988continuous, brun2002simple, korotkov1999continuous, jacobs2006straightforward, chantasri2013action}. This process, when written in the infinitesimal form, can be described by some measurement operator $\hat{\mathcal{M}}_{\delta t}$ such that the density matrix $\hat{\rho}(t)$ at later time $t+\delta t$ changes to
\begin{equation} \label{eq1}
\hat{\rho}(t+\delta t)=\frac{\hat{\mathcal{U}}_{\delta t} \hat{\rho}(t) \hat{\mathcal{U}}_{\delta t}^{\dagger}}{\operatorname{Tr}{[\hat{\mathcal{U}}_{\delta t} \hat{\rho}(t) \hat{\mathcal{U}}_{\delta t}^{\dagger}]}}, \quad \hat{\mathcal{U}}_{\delta t} \equiv e^{-i \hat{H}_s(t) \delta t} \hat{\mathcal{M}}_{\delta t},\end{equation}
where $\hat{H}_s(t)$ denotes the Hamiltonian of the system, $\hat{\mathcal{U}}_{\delta t}$ denotes the total evolution contributed by both $\hat{H}_s(t)$ and the measurement effect, and $\operatorname{Tr}{[\hat{\mathcal{U}}_{\delta t} \hat{\rho}(t) \hat{\mathcal{U}}_{\delta t}^{\dagger}}]$ serves as a normalization factor. Now if we consider the limit where the continuous quantum measurement disturbs the system strongly and happens extremely fast, then the whole process essentially becomes the projective measurement \cite{caves1987quantum} and can therefore be regarded as the continuous collapse of the wave function. In this extreme case, the effect from the measurement dominates such that we can take $\hat{H}_s(t) \approx 0,$ which will be assumed to be the case for the rest of the paper. Eq.~(\ref{eq1}) can thereby be rewritten as $\hat{\rho}(t+\delta t)=\frac{\hat{\mathcal{M}}_{\delta t} \hat{\rho}(t) \hat{\mathcal{M}}_{\delta t}^{\dagger}}{\operatorname{Tr}{[\hat{\mathcal{M}}_{\delta t} \hat{\rho}(t) \hat{\mathcal{M}}_{\delta t}^{\dagger}]}}$. For a given $\mathcal{\hat{M}}_{\delta t}$, we can then calculate the evolution of $\hat{\rho}(t)$ and record its evolution history as quantum trajectories \cite{Carmichael1993, wiseman1996quantum, brun2002simple, Bassi_2005, daley2014quantum}.

\textit{Reconstruct the Hamiltonian from the quantum trajectories.}---For simplicity, let us first consider the case where $\hat{\rho}(0)$ is initially pure. By requiring that during the projective measurement process, the quantum state does not become more mixed (since information is being extracted from the system), we can see that $\hat{\rho}(t)$ will remain pure at any later time $t$. This indicates that the evolution is unitary and can therefore be described by the Schr\"{o}dinger equation, at least effectively or phenomenologically. As a result, the problem becomes that with a given history of quantum trajectory $\ket{\psi(t)}$, which could either be calculated theoretically from certain $\hat{\mathcal{M}}_{\delta t}$, or reconstructed experimentally from measurement records (\textit{e.g}. via quantum tomography \cite{murch2013observing, weber2014mapping, jullien2014quantum, PhysRevX.6.041052}), how one can find an effective measuring Hamiltonian $\hat{H}(t)$ acting on the system alone that is responsible for the process. This problem, which is equivalent to solving the Schr\"{o}dinger equation inverse problem, has been solved in our recent work \cite{hu2022quantum}. We showed that for a given continuous quantum trajectory $\ket{\psi(t)}$, one can reconstruct its corresponding Hamiltonian $\hat{H}(t)$ by
\begin{equation} \label{eq2}
	\hat{H}(t) = i (\ketz{\p_t \tilde{\psi}(t)}\braz{\tilde{\psi}({t})}-\ketz{\tilde{\psi}(t)}\braz{\p_t\tilde{\psi}(t)})+\dot{\phi}(t) \mathbbm{1},
\end{equation}
where $\phi(t)\equiv\int-i\left\langle\partial_{t} \psi(t)|\psi(t)\right\rangle d t \in \mathbb{R}$ and $\ketz{\tilde{\psi}(t)}\equiv e^{i \phi(t)}\ket{\psi(t)}$, such that $\hat{H}(t)\ket{\psi(t)}=i \ket{\p_t \psi(t)}$. Then the general procedures for finding the ``measuring Hamiltonian'' responsible for the continuous wave function collapse of pure states are as follows:
\begin{enumerate}
	\item Assuming the process to be time-continuous, calculate theoretically (or measure experimentally) how the quantum state $\ket{\psi(t)}$ changes infinitesimally to $\ket{\psi(t+\delta t)}$, for instance, via certain  measurement operator $\hat{\mathcal{M}}_{\delta t}$,
\begin{equation} \label{eq3}
	\ket{\psi(t+\delta t)}=\f{\hat{\mathcal{M}}_{\delta t}}{\sqrt{\bra{\psi(t)}\hat{\mathcal{M}}_{\delta t}^\dagger \hat{\mathcal{M}}_{\delta t} \ket{\psi(t)}}} \ket{\psi(t)}.
\end{equation}
where ${\sqrt{\bra{\psi(t)}\hat{\mathcal{M}}_{\delta t}^\dagger \hat{\mathcal{M}}_{\delta t} \ket{\psi(t)}}}$ serves merely as a normalization factor.
	\item Calculate $\ket{\p_t \psi(t)}$ using its limit definition (or approximate with experimental data),
	\begin{equation} \label{eq4}
		\ket{\p_t \psi(t)}=\lim_{\delta t\to0} \f{\ket{\psi(t+\delta t)}-\ket{\psi(t)}}{\delta t}.
	\end{equation}
	\item Make use of Eq.~(\ref{eq2}) to reconstruct the measuring Hamiltonian $\hat{H}(t)$.
\end{enumerate}

Note that the measuring Hamiltonian is expected to generally be time-dependent and stochastic, such that it is compatible with the random nature of quantum mechanics and the Born rule. More interestingly, it is also state-dependent, meaning it will generally be different for different state. Such time- and state-dependency is determined by nature from the real-time interaction between the classical probe and the quantum object. In the above calculation. the stochasticity of the Hamiltonian comes from the stochasticity of quantum trajectories $\ket{\psi(t)}$. Physically, the stochasticity in the Hamiltonian may come from various sources, such as the vacuum fluctuation \cite{murch2013observing}, or more predominantly, just from the enormous number of degrees of freedom in the interaction between classical measuring devices and quantum objects. This interpretation differs from that of objective collapse theories, which holds that the quantum mechanics is intrinsically stochastic, even for the evolution of an isolated system \cite{RevModPhys.85.471}. It also differs from that of all sorts of hidden variable theories, since our formalism assumes it is the measuring device itself that causes the random nature of quantum mechanics.

Below we will show a few examples demonstrating our formalism.

\textit{Example I -- Two level system.}---To begin with, let us first consider the wave function of a two level system initially being at $\ket{\psi(0)}=a(0)\ket{0}+b(0)\ket{1}$ where $|a(0)|^2+|b(0)|^2=1$. We take our measuring operator $\hat{\mathcal{M}}_{\delta t}$ to be \cite{jacobs2014quantum}
\begin{equation} \label{eq5}
	\hat{\mathcal{M}}_{\delta t}=(\delta t / 2 \pi \tau)^{1 / 4} \exp \left[-\frac{\delta t}{4 \tau}\left(r_t \mathbbm{1}-\hat{\sigma}_{z}\right)^{2}\right],
\end{equation}
the physical meaning of which is that during the short time interval $[t,t+\delta t]$, we are measuring the $z$ component of the qubit with the measuring strength $1/\tau$, and the instantaneous readout $r_t \in (-\infty, +\infty)$ \cite{PhysRevLett.60.1351} at time $t$, which is a random variable, follows the conditional probability distribution $P(r_t\,|\,\rho(t))=\operatorname{Tr}[\hat{\rho}(t) \hat{\mathcal{M}}_{\delta t}\hat{\mathcal{M}}_{\delta t}^\dagger]$. The Gaussian part of the above equation indicates that the measurement effect is weighted towards the two eigenvalues $+1$ and $-1$ of the Pauli matrix $\hat{\sigma}_z$, so that the outputs of sequential applications of $\hat{\mathcal{M}}_{\delta t}$ will at last conform with the empirical facts that quantum states are always found to be at their eigenstates after von Neumann measurements. Note that different measurement approaches may give rise to different choices of $\hat{\mathcal{M}}_{\delta t}$. Assuming the measurement is strong so that the Hamiltonian of the system itself can be ignored, one then may calculate how $\ket{\psi(0)}$ will evolve under sequential applications of $\hat{\mathcal{M}}_{\delta t}$, until it stabilizes around its eigenstate $\ket{0}$ or $\ket{1}$.

By applying the aforementioned formalism we developed (Eq.~(\ref{eq2})-(\ref{eq4})),  we obtain the measuring Hamiltonian
\begin{equation} \label{eq6}
\begin{aligned}
	\hat{H}(t)&=\begin{pmatrix}
	0&ia(t) b^*(t) \f{r_t}{\tau}\\
	-ia^*(t) b(t) \f{r_t}{\tau}&0
\end{pmatrix}\\
&=-\frac{r_t}{2 \tau}\left(-y(t) \hat{\sigma}_{x}+x(t) \hat{\sigma}_{y}\right),
\end{aligned}
\end{equation}
where $x(t)=-\operatorname{Im}[a(t)b^*(t)]/2$ and $y(t)=\operatorname{Re}[a(t)b^*(t)]/2$ are the components of the Bloch vector $\vec{v}_\rho$ of the density matrix. It can be easily verified that the $\hat{\mathcal{M}}_{\delta t}$ introduced in Eq.\,(\ref{eq5}) and the Hamiltonian we just obtained result in exactly the same evolution of $\ket{\psi(t)}$. The random variable $r_t$ is chosen from the distribution $P(r_t\,|\,\rho(t))$ depending on $a(t)$ and $b(t)$, both of which go into $H(t)$, which in turn determines the value of $a(t+\delta t)$ and $b(t+\delta t)$.

\begin{figure}[!htp]
\subfloat{\includegraphics[width=.36\textwidth,scale=1]{plot01.eps}\hfill}
\caption{\label{fig:epsart1} The Bloch vector $\vec{v}_\rho$ of qubit is always orthogonal to and rotates around the Bloch vector $\vec{v}_H$ of the Hamiltonian given by Eq.\,(\ref{eq6}) during the measurement process, as visualized for two-level case. This is also the case in $n$-level system. The qubit then moves stochastically on the displayed longitudinal arc, until it stabilizes at $\ket{0}$ or $\ket{1}$. When one simultaneously measure the qubit along different axes (e.g. $\hat{x}$ and $\hat{z}$), it will make the qubit subject to the rotations of two different axes, $\vec{v}_{H_1}$ and $\vec{v}_{H_2}$, which is highly chaotic, and also makes the original eigenstates no longer be stabilization points.}
\end{figure}

There are a few comments worth mentioning regarding the Hamiltonian. First, it has instantaneous eigenvalues $E_{\pm}=\pm|a(t)||b(t)| \frac{r_t}{\tau}$ and instantaneous energy variance $[\Delta H(t)]^2= |a(t)|^2 |b(t)|^2 \f{r_t^2}{\tau^2}$ such that when we make a very strong measurement (\textit{i.e.} $1/\tau \to \infty$), its eigenvalues and energy variance diverge (or more physically, become very large). Not only does this validate our earlier treatment $H_s(t) \approx 0$, but also it is consistent with the Heisenberg uncertainty principle $\Delta H \Delta t \gtrsim \hbar/2$, and the idea of quantum speed limit \cite{mandelstam1945uncertainty, margolus1998maximum, deffner2013energy, hu2022quantum}, that a rapidly changing quantum state inevitably has a large energy variance. The fact that $[\Delta H(t)]^2 \to \infty$ as $1/\tau \to \infty$ also echoes our earlier claim that all physical processes take finite amount of time to happen, since in reality the energy supplied during measurement cannot be infinite.

Second, when the quantum state reaches one of its eigenstate, $\ket{0}$ or $\ket{1}$, which means $a(t)=0$ or $b(t)=0$, the measuring Hamiltonian vanishes, so that the quantum state is stabilized and no longer evolves, had we assumed $H_s(t)\approx 0$. In addition, while the mean energy of the Hamiltonian is zero, an arbitrary real function $\dot{\phi}(t)$ can be added to the diagonal term such that $H(t) \to H(t)+\dot{\phi}(t) \mathbbm{1}$, which only changes the global $U(1)$ phase of $\ket{\psi(t)}$ and therefore can also be used to describe the measuring process. Finally, if we plot on the Bloch sphere the evolution of the qubit during the measurement, we notice that the qubit rotates exactly around the axis $-y(t)\hat{x}+x(t)\hat{y}$ pointed by the Bloch vector $\vec{v}_H$ of the Hamiltonian, and that the Bloch vector $\vec{v}_\rho=x(t)\hat{x}+y(t)\hat{y}$ of the qubit is orthogonal to $\vec{v}_H$ (Fig.\,\ref{fig:epsart1}). It can be shown that this conclusion actually also applies to the measurement on $n$-level system, if we reconstruct the Hamiltonian via Eq.\,(\ref{eq2}) (Supplemental Materials I). Moreover, this observation illustrates why we cannot make precise simultaneous measurements in the different bases, such as $\hat{x}$ and $\hat{z}$, as this will force the qubit to stochastically rotate around ``two'' different axes, $\vec{v}_{H_1}$ and $\vec{v}_{H_2}$, making the total evolution highly chaotic. It also makes the measuring Hamiltonian no longer vanish at the eigenstates of $\hat{\sigma}_x$ and $\hat{\sigma}_z$, so the qubit will continue to evolve and escape from its original eigenstates, making the readouts spread over a large region \cite{ruskov2010qubit, chantasri2018simultaneous}.

To experimentally verify and utilize the conclusions drawn from the above discussions, one can implement real-time quantum feedback control to ``freeze'' the qubit under the continuous quantum measurement. For instance, for a given known initial state $\ket{\psi(0)}$, one can apply a real-time feedback control Hamiltonian $\hat{H}_c(t)=-(\hat{H}_s(t)+\hat{H}_m(t))$, where $\hat{H}_m(t)$ denotes the measuring Hamiltonian defined by Eq.\,(\ref{eq6}). The total Hamiltonian acting on the qubit during the measurement then is $\hat{H}_s(t)+\hat{H}_m(t)+\hat{H}_c(t)=0$, meaning the qubit will be ``frozen''. Experimentally such scheme can be more easily implemented via kicked quantum nondemolition measurements, e.g., by applying periodic weak measurements with certain time interval $T$ such that $\hat{U}(T)=\exp(-i \hat{H}_s T)=I$, where we have assumed $\hat{H}_s$ to be time-independent. In this way, the contribution of the evolution from $\hat{H}_s$ is effectively zero in each measurement snapshot, and
the waiting time between two successive measurements provides sufficient room for feedback controls to be implemented \cite{PhysRevB.74.085307}.

Another way to experimentally verify and utilize the conclusions is related to the most probable path. In \cite{PhysRevA.88.042110}, it is shown by using the stochastic path integral and the variational principle, and later verified experimentally in \cite{weber2014mapping}, that there exists a most probable path of a continuously monitored quantum system, the evolution of which are subject to the effects from both the Hamiltonian of system and the measuring processes. We can thereby solve for an effective Hamiltonian $\hat{H}_\text{eff}(t)$ responsible for such most probable path via Eq.\,(\ref{eq2}-\ref{eq4}), and by applying the counter-Hamiltonian $-\hat{H}_\text{eff}(t)$, to cancel out the most probable path. Experiments can then be designed to check whether the most probable path vanishes under the application of $-\hat{H}_\text{eff}(t)$ so as to verify if $\hat{H}_\text{eff}(t)$ describes well the dynamics of a continuously monitored quantum system. In Supplemental Materials \rom{2}, we provide an example to show the calculation of $\hat{H}_\text{eff}(t)$ of a continuously monitored qubit with its Hamiltonian $\hat{H}=\f{\epsilon}{2}\hat{\sigma}_z$.

\textit{Example II -- $n$-level system.}---Now let us consider the measurement on an arbitrary $n$-level system initially at $\ket{\psi(0)}=\sum_i^n a_i(0)\ket{i}$ where $\sum_i^n|a_i(0)|^2=1$. Similarly, let the measurement operator be $\hat{\mathcal{M}}_{\delta t}=\sum_{i}^{n} M_{i}(t)|i\rangle\langle i|$ where $
	M_{i}(t)=(\delta t / 2 \pi \tau)^{1 / 4} \exp [-\frac{\delta t}{4 \tau}\left(r_t-\lambda_{i}\right)^{2}]$. The $\lambda_i$ represents the $i$-th eigenvalue corresponding to the $i$-th eigenstate of all possible outcomes. By using again Eq.\,(\ref{eq2})-(\ref{eq4}), it can be shown (Supplemental Materials \rom{3}) that the measuring Hamiltonian responsible for the wave function collapse is given by (do not confuse the imaginary number $i$ with the summation index $i$)
	\begin{equation}\label{eq7}
		[H(t)]_{i j}=i \frac{a_{i}(t) a^*_{j}(t)}{4 \tau} \sum_{k}\left|a_{k}(t)\right|^{2}\left(\eta_{i k}+\eta_{k j}\right),
	\end{equation}
	where $\eta_{i j} \equiv\left(\lambda_{i}-\lambda_{j}\right)\left(2 r_t-\lambda_{i}-\lambda_{j}\right)$. The instantaneous energy variance $[\Delta H(t)]^2$ is calculated to be
\begin{equation} \label{eq8}
	[\Delta H(t)]^2=\sum_{j}^{n}\left[\frac{\left|a_{j}(t)\right|^{2}}{16 \tau^{2}}\left(\sum_{i}^{n}\left|a_{i}(t)\right|^{2} \eta_{ij}\right)^{2}\right].
\end{equation}

\textit{Example III -- harmonic oscillator at the ground state.}---Let us consider the ground state wave function of an harmonic oscillator peaked around $x_0$ and at time $t=0$,
\begin{equation}
	\psi(x, 0)=\left(\frac{m \omega}{\pi \hbar}\right)^{1 / 4} \exp \left(-\frac{m \omega}{2 \hbar}\left(x-x_{0}\right)^{2}\right).
\end{equation}
The measurement operator for a weak position measurement can be given by \cite{caves1987quantum} $\hat{\mathcal{M}}_{\delta t}=\int(\delta t / 2 \pi \tau)^{1 / 4} \exp \left[-\frac{\delta t}{4 \tau}(r_t-x)^{2}\right]|x\rangle\langle x| d x$. Expressing the wave function $\psi(x, 0)$ by the ket state $\ket{\psi(0)}=\int \psi(x,0)\ket{x}dx$, and applying Eq.\,(\ref{eq2})-(\ref{eq4}), we obtain (Supplemental Materials \rom{4}) the measuring Hamiltonian at time $t=0$,
\begin{widetext}
\begin{equation}\label{eq10}
	\hat{H}(0)=i\hbar\int \f{\eta_{x x^\prime}}{4\tau}\sqrt{\f{m \omega}{\pi \hbar}}\exp{\left( -\f{m \omega ((x-x_0)^2+{(x^\prime}-x_0)^2)}{2\hbar}\right)} \ket{x}\bra{x^\prime}dxdx^\prime,
\end{equation}
\end{widetext}
where $\eta_{xx^\prime} \equiv (x-x^\prime)(2r_t-x-x^\prime)$, which has instantaneous energy uncertainty
\begin{equation}\label{eq11}
	\Delta H(0)=\sqrt{\frac{\hbar \left(\hbar+4 m \omega (r_t-x_0)^{2} \right)}{32 m^{2} \omega^{2}\tau^{2} }} \hbar.
\end{equation}
It can be shown that at time $t=\delta t$, the wave function is still a Gaussian, so we can apply the same procedure repeatedly, which in turn indicates that the wave function will always be a Gaussian at any later time. The corresponding measuring Hamiltonian $\hat{H}(t)$ will therefore also share a similar form as $\hat{H}(0)$, except that the $x_0$ and $\omega$ in $\hat{H}(0)$ will need to be replaced by some other numbers. More specifically, it can be shown by solving a differential equation that the position uncertainty $\Delta X(t)$ at time $t$ can be related to $\Delta X(0)$ by $[\Delta X(t)]^{2}=\frac{[\Delta X(0)]^{2} \tau}{[\Delta X(0)]^{2} t+\tau}$. By replacing $\omega$ with $\Delta X(t)$ and using Eq.\,(\ref{eq11}), we obtain
\begin{equation} \label{eq12}
	 \Delta H(t)=\frac{\sqrt{2 (r_t-x_t)^{2} [\Delta X(t)]^{2}+[\Delta X(t)]^{4}}}{2 \sqrt{2} \tau} \hbar,
	 \end{equation}
which approximates to $|r_t-x_t|\hbar/(2\sqrt{rt})$ as $\tau \to 0.$
Such expression of $\Delta H(t)$ gives the energetic resources needed to quantify a given measurement dynamics. It is interesting to see that if the readout $r_t$ is too far away from the mean position $x_t$, then $\Delta H(t)$ becomes large. It is also worth noting that the Hamiltonian $\hat{H}(0)$ is non-local, since $\bra{x}\hat{H}(0)\ket{x^\prime}$ is generally nonzero for $x \neq x^\prime$. This seems to imply that position measurements are intrinsically non-local in quantum mechanics --- imagine the extreme case that the position measurement on a free particle in a very large room (e.g. $\sim$100 km$^2$) and the particle is found almost instantaneously to be somewhere in the room, the process of which would require detectors placed all over the room with good timing coordination such that the measuring Hamiltonian has to be regarded effectively non-local.

\textit{Discussions on the non-uniqueness of the measuring Hamiltonian.}---The measuring Hamiltonian may not be unique. Two totally different Hamiltonians with different spectra may result in exactly the same quantum evolution. For example, both $\hat{H}=-\hat{\sigma}_z$ and $\hat{H}(t)=\f{1}{3}\mathbbm{1}-\f{2}{9}\sqrt{2}\cos{(2t)}\hat{\sigma}_x+\f{2}{9}\sqrt{2}\sin{(2t)}\hat{\sigma}_y-\f{8}{9}\hat{\sigma}_z$ drive $\ket{\psi(t)}=\f{e^{it}}{\sqrt{3}}\ket{0}+\sqrt{\f{2}{3}}e^{-i t}\ket{1}$, and there is no way to tell which one is the actual Hamiltonian responsible for the driving just from (hypothetically) observing the quantum trajectory of $\ket{\psi(t)}$. This implies that the measuring Hamiltonian we obtained from Eq.~(\ref{eq2})-(\ref{eq4}) may not be the actual Hamiltonian responsible for the wave function collapse. Is there a way to find the true measuring Hamiltonian? With only the information from the quantum trajectories $\ket{\psi(t)}$, the answer is negative. Nevertheless, we prove the following theorem in the Supplemental Materials \rom{5}.

\textit{Theorem}. \textit{Let $\ket{\psi(t)}\neq0$ be an arbitrary pure quantum state at time $t$. Let $V$ be the linear space formed by all Hermitian operators $\hat{T}(t)$ satisfying $\hat{T}(t)\ket{\psi(t)}=0.$ If one can find at least one Hamiltonian $\hat{H}(t)$ satisfying $\hat{H}(t)\ket{\psi(t)}=i \ket{\p_t \psi(t)}$, then one can find all possible Hamiltonians $\hat{H}^\prime(t)$ satisfying $\hat{H}^\prime(t)\ket{\psi(t)}=i \ket{\p_t \psi(t)}$. In particular, $\hat{H}^\prime(t)$ can be given by $\hat{H}^\prime(t)=\hat{H}(t)+\hat{T}(t)$ where $\hat{T}(t)$ is arbitrarily chosen from $V$.}

The above theorem implies that even though we may never be able to know the actual measuring Hamiltonian, we can be assured that it must be one of the possible $\hat{H}^\prime(t)$, if we are able to solve for all $\hat{T}(t)$. Moreover, all such Hamiltonians $H^\prime(t)$ should have the same energy variance equal to $\braket{\p_t \psi(t)}{\p_t \psi(t)}-|\braket{\psi(t)}{\p_t\psi(t)}|^2$.

\textit{Measurement on a mixed quantum state.}---Here we briefly discuss the generalization of the measurement on a mixed quantum state. Consider an ensemble of (proper) mixed quantum states $\rho(0)=\sum_i^n p_i(0) \ket{\psi_i(0)}\bra{\psi_i(0)}$ at $t=0$, which is prepared by mixing up $n$ orthonormal pure states $\ket{\psi_i(0)}$ with relative population $p_i(0)$, such that for a quantum object randomly selected from the ensemble, there is a classical probability $p_i(0)$ for the selected object to be in $\ket{\psi_i(0)}$. Now suppose one makes a measurement on the selected object, then according to the earlier theory we developed, one concludes that there should be a Hamiltonian responsible for the measuring process, since the selected quantum object is actually pure --- it is just that one has lost track (information) of the state. In this case, one can describe the measuring process by telling that at time $t=0$, there is a probability $p_i(0)$ that the measuring Hamiltonian is given by $\hat{H}_{\ket{\psi_i(0)}}$. Such description $\{ p_i(0)\}$ needs to be Bayesian updated at later time $t$ to $\{ p_i(t)\}$ once a series of readouts $\{r_t\}$ are accumulated, indicating the process of purification. Note that in this scheme, the selected quantum object actually remains pure all the time, corresponding to a unitary process. It is the observer's updating knowledge of the state that makes the evolution appear non-unitary for the observer.

\textit{Conclusion.}---In this Letter, we demonstrate how one can describe the wave function collapse process by using the Schr\"{o}dinger equation and the formalism we developed in \cite{hu2022quantum} and Eq.~(\ref{eq2})-(\ref{eq4}). We solve analytically the measuring Hamiltonian (Eq.(\ref{eq6}),(\ref{eq7}), and (\ref{eq10})) responsible for the process in varied cases, and calculate its instantaneous energy variance (Eq.(\ref{eq8}) and (\ref{eq12})). We then discuss the non-uniqueness of the measuring Hamiltonian, and how our results can also be applied to describe the wave function collapse process for mixed quantum states.

\textit{Acknowledgement.}---We are grateful to the support from the Army Research Office (ARO) under the grant W911NF-22-1-0258.