Time-Varying Engineered Reservoir for the Improved Estimation of Atom-Cavity Coupling Strength

Abstract

In this paper, we consider the application of quantum reservoir engineering in quantum metrology. More precisely, we are concerned with a system setup where a sequence of atoms constructing the “time-varying” quantum reservoir interact, in turn, with the trapped field in a cavity through the Jaynes–Cummings Hamiltonian. In particular, we were able to manipulate the initial states of reservoir atoms in order to enhance estimation precision regarding the coupling strength between each atom and the cavity (the coupling strength between each atom and the cavity was assumed to be identical). The novelty of this work lies in alternately preparing the atoms at two different states in a pairwise manner, such that the cavity could converge into a squeezed state with photonic loss to the environment taken into account. The control scheme proposed here thus leads to higher precision compared to the previous work where reservoir atoms were initialized at the same state, which drove the cavity to a coherent state. Detailed theoretical analysis and numerical simulations are also provided. In addition, this system setup and the associated control scheme are easily implemented for quantum metrology, since no entanglement is required for the preparation of atom states, and the final cavity state can stay steady.

1. Introduction

During the past few decades, quantum metrology has been playing an important role in quantum science and technology. The procedure of quantum parameter estimation in general consists of three steps: preparing the probe state, letting the prepared probe state interact with the system (parameterization process), and measuring the final state [1,2]. These three parts all influence the estimation precision. There is established work concerning the preparation of the probe state and the selection of the measurement for the final state. The standard approach is to prepare the optimal probe state with the utilization of appropriate positive-operator-valued measurement (POVM) after the parameterization process [3,4,5,6,7,8,9]. In terms of optimizing the quantum state possessing the information of the parameter to be estimated, Carlton Morris Caves proposed that the squeezed vacuum state could help in improving the estimation precision in 1981 [10]. Since then, quantum properties in nonclassical states began to attract the attention of more researchers, leading to a boom in the development of quantum metrology. For example, it was proven both theoretically and experimentally that entangled coherent states [3,4], spin squeezed states [5,6] and two-mode squeezed vacuum states [7] can contribute to better estimation precision. On the other hand, in order to further improve estimation precision, control was introduced in the parameterization process, such as the gradient ascent pulse engineering (GRAPE) and the asynchronous advantage actor–critic algorithms aimed at optimally adjusting the parameterization process [11,12,13]. The Heisenberg limit can also not be surpassed by merely adding control terms to the Hamiltonian [14,15,16]. However, most of the above-mentioned works that involve adding control need coherent signals to be directly fed to the target system, and entangled states are usually desired to enhance the estimation, which is indeed not simple to prepare in practice.

By contrast, quantum reservoir engineering provides an alternative way that can avoid adding real-time control directly to the target system. Rather, the target system is controlled indirectly by manipulating the reservoir. In this paper, the reservoir is composed of a sequence of atoms (qubits) that interact with the trapped field in a cavity (quantum harmonic oscillator) one after another through the Jaynes–Cummings Hamiltonian. The experimental background of this paper was first developed in [17]. Then, on the basis of this setup, it was proven that nonclassical states of the quantum harmonic oscillator could be stabilized through the interaction between the oscillator and the reservoir qubits [18,19], where it was demonstrated that the cavity could be stabilized at a coherent state with the atoms initialized at the same state, as shown in Figure 1. In particular, control signals do not need to be added to the oscillator, but the initial states of the reservoir qubits can be altered, such that the final state of the oscillator can be manipulated with the purpose of estimating the coupling strength between each qubit and the oscillator. This setting was used for quantum metrology in [20] where the qubits were initialized at an identical state; thus, the harmonic oscillator converged towards a coherent state with the averaged photon counting calculated to quantify the variance in coupling strength. Estimation precision was also proven to reach the Heisenberg limit in [20], where the root-mean-square fluctuation of the atom–cavity coupling strength was proportional to the square of the number of effective coupling atoms. Other reservoir engineering schemes are utilized to estimate qubit–oscillator coupling strength. For example, in [21], a quantum reservoir composed of N quantum harmonic oscillators interacting with the qubit was applied to estimate the coupling strength by calculating the quantum Fisher information (QFI) to quantify the precision. However, preparing the state of a qubit in [20] tended to be easier than preparing the joint state of N harmonic oscillators.

By preparing the reservoir qubits at different states, which may initially involve entanglement, the harmonic oscillator could be stabilized at a squeezed coherent state with the squeezing strength determined by the initial qubit states [22,23]. Mainly inspired by [20,23], we aim to take advantage of a “time-varying” reservoir to enhance the estimation precision of qubit–oscillator coupling strength without entanglement involved in the reservoir qubits. The term “time-varying” means that the reservoir qubits are alternately initialized at two different states in a pairwise manner, as shown in Figure 2, but not initialized at the same state as that shown in Figure 1. In more concrete terms, reservoir qubits are initialized at $|{\psi }_{qred}\rangle$, $|{\psi }_{qblue}\rangle$, $|{\psi }_{qred}\rangle$, $|{\psi }_{qblue}\rangle$, ⋯ such that the harmonic oscillator can converge into a squeezed coherent state, since squeezed states enable us to significantly improve the estimation precision [24]. Furthermore, we took into account photonic loss to the environment in this paper; in [23], only the ideal case was considered. Qubits initialized under the scheme shown in Figure 2, steering the harmonic oscillator to a squeezed state, help in achieving higher precision than that of qubits initialized under the scheme shown in Figure 1, steering the harmonic oscillator to a coherent state in the absence of entanglement in the reservoir, which is not difficult to be realized in experiments. It is, thus, promising that the squeezed states stabilized by means of quantum reservoir engineering, as shown in Figure 2, can be used to provide benefits to related works in the field of quantum information, such as [25,26,27,28].

The organization of this paper is as follows. In Section 2, we give the mathematical description of the system setup, especially concerning the consecutive pairs of separable “time-varying” reservoir qubits. In Section 3, we take into account photonic loss through the oscillator to the environment, providing the form of steady states together with the stability conditions. Furthermore, we prove that steering the oscillator to a squeezed state can give rise to a larger QFI compared to that of previous work where the oscillator was steered to a coherent state. This is also visualized and verified via simulation examples. Critical issues regarding our theoretical analysis and numerical simulation, together with potential future work are discussed in Section 4. Lastly, in Section 5, we conclude our present work.

2. Mathematical Model of the System Setup and Quantification of the Estimation Precision

In the system setting shown in Figure 2, quantum reservoir R is composed of a sequence of separable qubits that are alternately prepared at two different initial states. From a control-theory point of view, the initial states of qubits play the role of control variables, and the harmonic oscillator is considered the target system. In other words, our aim was to control the final state of the harmonic oscillator with the parameter to be estimated embedded, which could be used for quantum metrology. In this section, we recall the detailed mathematical description of this reservoir engineering system. A stabilized squeezed state with an associated QFI is also provided without taking into account photonic loss through the oscillator to the environment (the ideal case) [22,23,24].

Specifically, as shown in Figure 2, since the reservoir qubits were initialized at states $|{\psi }_{qred}\rangle$, $|{\psi }_{qblue}\rangle$, $|{\psi }_{qred}\rangle$, $|{\psi }_{qblue}\rangle$, ⋯, we could regard one pair of qubits as a whole. In particular, $|{\psi }_{qred}\rangle =cos{u}_1|g\rangle +sin{u}_1|e\rangle$ and $|{\psi }_{qblue}\rangle =cos{u}_2|g\rangle +sin{u}_2|e\rangle$, where the parameters $u_{1,2}\in [0,\frac{\pi }{4})$. According to the results in [23], we know that when the initial phases of the qubits are the same, it could benefit the precision; thus, the initial phases of the qubits were set to zero. Equivalently, the initial state for one pair of qubits can also be described by

$$\begin{array}{cc}|{\psi }_{q^2}\left(0\right)\rangle =& cos{u}_{1}cos{u}_2|gg\rangle +sin{u}_{1}sin{u}_2|ee\rangle \\ \hfill & +cos{u}_{1}sin{u}_2|ge\rangle +sin{u}_{1}cos{u}_2|eg\rangle .\hfill \end{array}$$

(1)

Each qubit interacts with the harmonic oscillator for a short period of time $t_r$ according to the following Jaynes–Cummings Hamiltonian:

$$\begin{array}{c}\hfill {\mathit{H}}_{JC}=i\frac{\mathrm{\Omega }}{2}\left(\right|g\rangle \langle e|{\mathit{a}}^{†}-|e\rangle \langle g\left|\mathit{a}\right),\end{array}$$

(2)

where $\mathrm{\Omega }$ the Rabi oscillation frequency, $\mathit{a}$ is the oscillator mode’s annihilation operator, and $|g\rangle and|e\rangle$ are the qubit’s ground and excited states, respectively. The unitary propagator describing qubit–oscillator interaction is then

$$\begin{array}{c}\hfill {\mathit{U}}_r=|g\rangle \langle g|cos{\theta }_{\mathit{N}}+|e\rangle \langle e|cos{\theta }_{\mathit{N}+\mathit{I}}-|e\rangle \langle g|\mathit{a}\frac{sin{\theta }_{\mathit{N}}}{\sqrt{\mathit{N}}}+|g\rangle \langle e|\frac{sin{\theta }_{\mathit{N}}}{\sqrt{\mathit{N}}}{\mathit{a}}^{†},\end{array}$$

(3)

where

$$\begin{array}{c}\hfill {\theta }_{\mathit{N}}=\theta \sqrt{\mathit{N}}=\frac{1}{2}\mathrm{\Omega }{t}_r\sum _n\sqrt{n}|n\rangle \langle n|,\end{array}$$

where $\mathit{N}={\mathit{a}}^{†}\mathit{a}$ is the photon number operator, $|n\rangle (n=0,1,2,...)$ are the Fock states of the harmonic oscillator mode, and $\mathit{I}$ is the identity operator. Qubit–oscillator interactions operate in the weakly coupled regime; thus, $\theta =\frac{1}{2}\mathrm{\Omega }{t}_r$ (the effective qubit–oscillator coupling strength to be estimated) is sufficiently small.

We first briefly present the scenario where all qubits are initialized at the same state, i.e., $|\psi \rangle =cosu|g\rangle +sinu|e\rangle$ with $u\in [0,\frac{\pi }{4})$, as shown in Figure 1.

By assuming that $\theta$ is sufficiently small, the oscillator state can be stabilized at coherent state $|\alpha \rangle$ with $\alpha =\frac{tan2u}{\theta }$ (u is required to be sufficiently small in order to guarantee the stabilization). The corresponding convergence rate is ${\kappa }_c={\theta }^{2}cos2u$.

Without considering the dissipation to the environment, we now focus on the scenario where the reservoir qubits are initialized at states $|{\psi }_{qred}\rangle$, $|{\psi }_{qblue}\rangle$, $|{\psi }_{qred}\rangle$, $|{\psi }_{qblue}\rangle$, ⋯, as shown in Figure 2. Following a similar procedure to that in [23], the Lindblad master equation for oscillator state $\rho$ by taking the continuous-time approximation is as follows:

$$\begin{array}{c}\hfill \frac{d}{dt}\rho \left(t\right)=-i[\mathit{H},\rho \left(t\right)]+\sum _{j=1}^3\mathcal{L}\left({\mathit{L}}_j\right)\rho \left(t\right),\end{array}$$

(4)

Superoperator $\mathcal{L}\left(\mathit{L}\right)$ denotes the Lindbladian that is defined by

$$\begin{array}{c}\hfill \mathcal{L}\left(\mathit{L}\right)\rho \left(t\right)=\mathit{L}\rho \left(t\right){\mathit{L}}^{†}-\frac{1}{2}\rho \left(t\right){\mathit{L}}^{†}\mathit{L}-\frac{1}{2}{\mathit{L}}^{†}\mathit{L}\rho \left(t\right),\end{array}$$

with $\mathit{L}$ being the coupling operator. The Hamiltonian and coupling operators in Equation (4) are given by

$$\begin{array}{cc}\hfill \mathit{H}& =-i\theta (\mathit{Q}-{\mathit{Q}}^{†}),\mathit{Q}=cos(u_1-u_2)sin(u_1+u_2)\mathit{a},\hfill \\ \hfill {\mathit{L}}_{\mathit{1}}& =\sqrt{2}\theta cos{u}_{1}cos{u}_2\mathit{a}-\sqrt{2}\theta sin{u}_{1}sin{u}_2{\mathit{a}}^{†},\hfill \\ \hfill {\mathit{L}}_{\mathit{2}}& =\theta sin(u_1+u_2)\mathit{a},{\mathit{L}}_{\mathit{3}}=\theta sin(u_1+u_2){\mathit{a}}^{†}.\hfill \end{array}$$

In order to stabilize the oscillator state, we require

$$\begin{array}{c}\hfill sin(u_1+u_2)\mathrm{is}\mathrm{sufficiently}\mathrm{small},\mathrm{and}cos(u_1+u_2)>0,\end{array}$$

such that ${\mathit{L}}_{\mathit{1}}$ dominates the dynamics. Then, reservoir qubits can drive the oscillator to a squeezed state $|\alpha ,r\rangle$ with the parameters

$$\begin{array}{cc}\hfill & \alpha =\frac{tan(u_1+u_2)}{\theta },\hfill \end{array}$$

(5a)

$$\begin{array}{cc}\hfill & r=\frac{1}{2}ln\frac{cos(u_1-u_2)}{cos(u_1+u_2)}.\hfill \end{array}$$

(5b)

In general, $\alpha$ can be complex for a squeezed state of the $|\alpha ,r\rangle$ form. However, $\alpha$ in Equation (5a) is real by initially preparing the qubits at $|{\psi }_{qred}\rangle$, $|{\psi }_{qblue}\rangle$, $|{\psi }_{qred}\rangle$, $|{\psi }_{qblue}\rangle$, ⋯.

This can be proven by applying the displacement and squeezing transformation to $\rho \left(t\right)$ as follows:

$$\begin{array}{cc}\hfill \tilde{\rho }\left(t\right)& ={\mathit{S}}^{†}\left(r\right){\mathit{D}}^{†}\left(\alpha \right)\rho \left(t\right)\mathit{D}\left(\alpha \right)\mathit{S}\left(r\right),\hfill \end{array}$$

(6)

$$\begin{array}{cc}\hfill \tilde{\mathit{H}}& ={\mathit{S}}^{†}\left(r\right){\mathit{D}}^{†}\left(\alpha \right)\mathit{H}\mathit{D}\left(\alpha \right)\mathit{S}\left(r\right),\hfill \end{array}$$

(7)

$$\begin{array}{cc}\hfill \widetilde{{\mathit{L}}_{\mathit{1}}}& ={\mathit{S}}^{†}\left(r\right){\mathit{D}}^{†}\left(\alpha \right){\mathit{L}}_{\mathit{1}}\mathit{D}\left(\alpha \right)\mathit{S}\left(r\right),\hfill \end{array}$$

(8)

which results in the following.

$$\begin{array}{c}\hfill \frac{\mathrm{d}\tilde{\rho }\left(t\right)}{\mathrm{d}t}={\theta }^2\left(\frac{{cos}^3(u_1-u_2)}{2cos(u_1+u_2)}+\frac{{cos}^3(u_1+u_2)}{2cos(u_1-u_2)}+cos(u_1-u_2)(u_1+u_2)\right)\mathcal{L}\left(\mathit{a}\right)\tilde{\rho }\left(t\right),\end{array}$$

(9)

with $\alpha$ and r given in Equation (5). By defining the Lyapunov function:

$$\begin{array}{c}\hfill V\left(t\right)=Tr\left({\mathit{a}}^{†}\mathit{a}{\mathit{S}}^{†}\left(r\right){\mathit{D}}^{†}\left(\alpha \right)\rho \left(t\right)\mathit{D}\left(\alpha \right)\mathit{S}\left(r\right)\right),\end{array}$$

we can further obtain that

$$\begin{array}{c}\hfill \frac{d}{dt}V\left(t\right)=Tr\left({\mathit{a}}^{†}\mathit{a}\dot{\tilde{\rho }}\left(t\right)\right)=-2{\theta }^2({cos}^2{u}_1{cos}^2{u}_2-{sin}^2{u}_1{sin}^2{u}_2)V\left(t\right)=-{\kappa }_{s}V\left(t\right),\end{array}$$

(10)

Thus, ${\kappa }_s$ denotes the convergence rate.

For a squeezed state $|\alpha ,r\rangle$, one can calculate the QFI ${\mathcal{F}}_Q$ by using the following method (see, e.g., [29]),

$$\begin{array}{c}\hfill {\mathcal{F}}_Q=2\frac{\mathrm{d}{\nu }^{†}}{\mathrm{d}\theta }{\sigma }^{-1}\frac{\mathrm{d}\nu }{\mathrm{d}\theta },\end{array}$$

where

$$\begin{array}{c}\hfill \nu ={\left(\alpha ,{\alpha }^{*}\right)}^T,\sigma =\left[\begin{array}{cc}cosh2r& sinh2r\\ sinh2r& cosh2r\end{array}\right].\end{array}$$

This formula can be derived on the basis of the definition of the QFI associated with the Bures distance between ${\rho }_{\theta }$ and ${\rho }_{\theta +\mathrm{d}\theta }$ as follows:

$$\begin{array}{c}\hfill {\mathrm{d}}_{Bures}^2\left({\rho }_{\theta },{\rho }_{\theta +\mathrm{d}\theta }\right)=\frac{1}{4}{\mathcal{F}}_Q\mathrm{d}{\theta }^2.\end{array}$$

In this stabilization scenario, we are dealing with the steady state; therefore, the time cost is supposed to be included to evaluate estimation in a control protocol [30]. The appropriate QFI for the steady state is, thus, given by

$$\begin{array}{c}\hfill {\mathcal{F}}_s={\mathcal{F}}_Q{\kappa }_s^2\end{array}$$

(11)

by incorporating the convergence rate. The QFI associated with the squeezed state

$|\frac{tan(u_1+u_2)}{\theta },\frac{1}{2}ln\frac{cos(u_1-u_2)}{cos(u_1+u_2)}\rangle$ is (subscript s stands for “squeezed state”)

$${\mathcal{F}}_s=\frac{16{sin}^2(u_1+u_2){cos}^3(u_1-u_2)}{cos(u_1+u_2)}.$$

(12)

3. Quantification of the Estimation Precision in the Presence of Photonic Loss to the Environment

In Section 2, we present the mathematical description of the system setup, and calculate QFI on the basis of the results in [23] without considering dissipation to the environment. However, in practice, it is inevitable that the qubit–oscillator quantum reservoir engineering system undergoes photonic loss to the environment during the whole process. Therefore, in this section, the associated QFI is amended due to the noisy dynamics of open quantum systems [22,31]. We included the energy decay of the harmonic oscillator via the annihilation operator in a zero-temperature bath. Over a short period of time $t_r$ compared to the oscillator characteristic lifetime $1/\gamma$, the corresponding effect could be modeled with the Kraus map: $\rho \mapsto {\mathit{M}}_0\rho {\mathit{M}}_0^{†}+{\mathit{M}}_1\rho {\mathit{M}}_1^{†}$, where the propagators

$$\begin{array}{cc}\hfill {\mathit{M}}_0& =\mathit{I}-\frac{\gamma t_r}{2}{\mathit{a}}^{†}\mathit{a},{\mathit{M}}_1=\sqrt{\gamma t_r}\mathit{a}.\hfill \end{array}$$

(13)

3.1. Convergence towards a Coherent State

We first consider the scenario where each qubit is prepared at the same state that interacts with the harmonic oscillator, as shown in Figure 1, taking into account the energy decay from the oscillator to the environment. The following Kraus map could then be used to describe the evolution of the oscillator state:

$$\begin{array}{cc}\hfill \rho (t+1)=& {\mathit{M}}_0({\mathit{M}}_g\rho \left(t\right){\mathit{M}}_g^{†}+{\mathit{M}}_e\rho \left(t\right){\mathit{M}}_e^{†}){\mathit{M}}_0^{†}\hfill \\ \hfill & +{\mathit{M}}_1({\mathit{M}}_g\rho \left(t\right){\mathit{M}}_g^{†}+{\mathit{M}}_e\rho \left(t\right){\mathit{M}}_e^{†}){\mathit{M}}_1^{†}.\hfill \end{array}$$

(14)

Assuming that $\theta \ll 1$ and $\gamma t_r\ll 1$, we expanded the Kraus map given in Equation (14) to the second order in $\theta$ and the first order in $\gamma t_r$, and the following approximated Lindblad master equation could be obtained:

$$\begin{array}{cc}\hfill \frac{\mathrm{d}\rho \left(t\right)}{\mathrm{d}t}=& {\theta }^2{cos}^{2}u\mathcal{L}\left(\mathit{a}\right)\rho \left(t\right)+{\theta }^2{sin}^{2}u\mathcal{L}\left({\mathit{a}}^{†}\right)\rho \left(t\right)+\theta cosusinu[\rho \left(t\right),\mathit{a}-{\mathit{a}}^{†}]\hfill \\ \hfill & -\frac{\gamma t_r}{2}{\mathit{a}}^{†}\mathit{a}\rho \left(t\right)-\frac{\gamma t_r}{2}\rho \left(t\right){\mathit{a}}^{†}\mathit{a}+\gamma t_r\mathit{a}\rho \left(t\right){\mathit{a}}^{†}.\hfill \end{array}$$

(15)

By imposing the displacement transformation onto $\rho \left(t\right)$, we could obtain the corresponding Lindblad master equation for the transformed $\tilde{\rho }\left(t\right)$:

$$\begin{array}{cc}\hfill \frac{\mathrm{d}\tilde{\rho }\left(t\right)}{\mathrm{d}t}=& ({\theta }^2{cos}^{2}u+\gamma t_r)\mathcal{L}\left(\mathit{a}\right)\tilde{\rho }\left(t\right)+{\theta }^2{sin}^{2}u\mathcal{L}\left({\mathit{a}}^{†}\right)\tilde{\rho }\left(t\right)\hfill \\ \hfill & +(\frac{1}{2}{\theta }^2{cos}^{2}u\alpha -\theta sinucosu-\frac{1}{2}{\theta }^2{sin}^{2}u\alpha +\frac{\alpha \gamma t_r}{2})[\tilde{\rho }\left(t\right),{\mathit{a}}^{†}]\hfill \\ \hfill & +(\frac{1}{2}{\theta }^2{sin}^{2}u{\alpha }^{*}+\theta sinucosu-\frac{1}{2}{\theta }^2{cos}^{2}u{\alpha }^{*}-\frac{{\alpha }^{*}\gamma t_r}{2})[\tilde{\rho }\left(t\right),\mathit{a}].\hfill \end{array}$$

(16)

On the basis of the equation of motion above, the coherent state $|\alpha \rangle$ with $\alpha =\frac{\theta sin2u}{{\theta }^{2}cos2u+\gamma t_r}$ of the oscillator could be stabilized, provided that $sinu$ was sufficiently small (e.g. $sinu\ll 1$). Convergence was guaranteed by the ${\theta }^{2}cos2u+\gamma t_r\ge 0$ condition, which could be obtained by considering Lyapunov equation $V\left(t\right)=Tr\left({\mathit{a}}^{†}\mathit{a}{\mathit{D}}^{†}\left(\alpha \right)\rho \left(t\right)\mathit{D}\left(\alpha \right)\right)$ and its derivative

$$\begin{array}{cc}\hfill \frac{d}{dt}V\left(t\right)& =2({\theta }^2{sin}^{2}u-{\theta }^2{cos}^{2}u-\gamma t_r)Tr\left(\tilde{\rho }\left(t\right)a^{†}a\right)\hfill \\ \hfill & =-2({\theta }^{2}cos2u+\gamma t_r)V\left(t\right)=-{\kappa }_{dc}V\left(t\right).\hfill \end{array}$$

(17)

In order to realize $\frac{d}{dt}V\left(t\right)\le 0$, ${\kappa }_{dc}=2\left({\theta }^{2}cos2u+\gamma t_r\right)\ge 0$ is required, with ${\kappa }_{dc}$ as the convergence rate.

In a similar way, to calculate the QFI, as presented in Section 2, we could obtain the QFI regarding a coherent state when dissipation to the environment is taken into account as follows (subscripts d and c stand for “dissipation” and “coherent”, respectively):

$$\begin{array}{cc}\hfill {\mathcal{F}}_{dc}& =\frac{16{sin}^{2}2u{(\gamma t_r-{\theta }^{2}cos2u)}^2}{{({\theta }^{2}cos2u+\gamma t_r)}^2}.\hfill \end{array}$$

(18)

3.2. Convergence towards a Squeezed State

In this part, we focus on the scenario where the reservoir qubits are initialized at the states $|{\psi }_{qred}\rangle$, $|{\psi }_{qblue}\rangle$, $|{\psi }_{qred}\rangle$, $|{\psi }_{qblue}\rangle$, ⋯, as shown in Figure 2. These qubits interact with the harmonic oscillator in the presence of oscillator’s phton loss to the environment. The evolution of the oscillator state after it interacted with one pair of qubits can be described with the following Kraus map:

$$\begin{array}{cc}\hfill \rho (t+1)=& {\mathit{M}}_0{\mathit{M}}_0({\mathit{M}}_{gg}\rho \left(t\right){\mathit{M}}_{gg}^{†}+{\mathit{M}}_{ge}\rho \left(t\right){\mathit{M}}_{ge}^{†}+{\mathit{M}}_{eg}\rho \left(t\right){\mathit{M}}_{eg}^{†}+{\mathit{M}}_{ee}\rho \left(t\right){\mathit{M}}_{ee}^{†}){\mathit{M}}_0^{†}{\mathit{M}}_0^{†}\hfill \\ \hfill +& {\mathit{M}}_0{\mathit{M}}_1({\mathit{M}}_{gg}\rho \left(t\right){\mathit{M}}_{gg}^{†}+{\mathit{M}}_{ge}\rho \left(t\right){\mathit{M}}_{ge}^{†}+{\mathit{M}}_{eg}\rho \left(t\right){\mathit{M}}_{eg}^{†}+{\mathit{M}}_{ee}\rho \left(t\right){\mathit{M}}_{ee}^{†}){\mathit{M}}_1^{†}{\mathit{M}}_0^{†}\hfill \\ \hfill +& {\mathit{M}}_1{\mathit{M}}_0({\mathit{M}}_{gg}\rho \left(t\right){\mathit{M}}_{gg}^{†}+{\mathit{M}}_{ge}\rho \left(t\right){\mathit{M}}_{ge}^{†}+{\mathit{M}}_{eg}\rho \left(t\right){\mathit{M}}_{eg}^{†}+{\mathit{M}}_{ee}\rho \left(t\right){\mathit{M}}_{ee}^{†}){\mathit{M}}_0^{†}{\mathit{M}}_1^{†}\hfill \\ \hfill +& {\mathit{M}}_1{\mathit{M}}_1({\mathit{M}}_{gg}\rho \left(t\right){\mathit{M}}_{gg}^{†}+{\mathit{M}}_{ge}\rho \left(t\right){\mathit{M}}_{ge}^{†}+{\mathit{M}}_{eg}\rho \left(t\right){\mathit{M}}_{eg}^{†}+{\mathit{M}}_{ee}\rho \left(t\right){\mathit{M}}_{ee}^{†}){\mathit{M}}_1^{†}{\mathit{M}}_1^{†}.\hfill \end{array}$$

where the propagators

$$\begin{array}{cc}\hfill {\mathit{M}}_{gg}=& cos{u}_{1}cos{u}_2{cos}^2{\theta }_{\mathit{N}}+cos{u}_{1}sin{u}_{2}cos{\theta }_{\mathit{N}}\frac{sin{\theta }_{\mathit{N}}}{\sqrt{\mathit{N}}}{\mathit{a}}^{†}\hfill \\ \hfill & +sin{u}_{1}cos{u}_2\frac{sin{\theta }_{\mathit{N}}}{\sqrt{\mathit{N}}}{\mathit{a}}^{†}cos{\theta }_{\mathit{N}}+sin{u}_{1}sin{u}_2\frac{sin{\theta }_{\mathit{N}}}{\sqrt{\mathit{N}}}{\mathit{a}}^{†}\frac{sin{\theta }_{\mathit{N}}}{\sqrt{\mathit{N}}}{\mathit{a}}^{†},\hfill \\ \hfill {\mathit{M}}_{ge}=& -cos{u}_{1}cos{u}_{2}cos{\theta }_{\mathit{N}}\mathit{a}\frac{sin{\theta }_{\mathit{N}}}{\sqrt{\mathit{N}}}+cos{u}_{1}sin{u}_{2}cos{\theta }_{\mathit{N}}cos{\theta }_{\mathit{N}+\mathit{I}}\hfill \\ \hfill & -sin{u}_{1}cos{u}_2{sin}^2{\theta }_{\mathit{N}}+sin{u}_{1}sin{u}_2\frac{sin{\theta }_{\mathit{N}}}{\sqrt{\mathit{N}}}{\mathit{a}}^{†}cos{\theta }_{\mathit{N}+\mathit{I}},\hfill \\ \hfill {\mathit{M}}_{eg}=& -cos{u}_{1}cos{u}_2\mathit{a}\frac{sin{\theta }_{\mathit{N}}}{\sqrt{\mathit{N}}}cos{\theta }_{\mathit{N}}-cos{u}_{1}sin{u}_2{sin}^2{\theta }_{\mathit{N}+\mathit{I}}\hfill \\ \hfill & +sin{u}_{1}cos{u}_{2}cos{\theta }_{\mathit{N}+\mathit{I}}cos{\theta }_{\mathit{N}}+sin{u}_{1}sin{u}_{2}cos{\theta }_{\mathit{N}+\mathit{I}}\frac{sin{\theta }_{\mathit{N}}}{\sqrt{\mathit{N}}}{\mathit{a}}^{†},\hfill \\ \hfill {\mathit{M}}_{ee}=& cos{u}_{1}cos{u}_2\mathit{a}\frac{sin{\theta }_{\mathit{N}}}{\sqrt{\mathit{N}}}\mathit{a}\frac{sin{\theta }_{\mathit{N}}}{\sqrt{\mathit{N}}}-cos{u}_{1}sin{u}_2\mathit{a}\frac{sin{\theta }_{\mathit{N}}}{\sqrt{\mathit{N}}}cos{\theta }_{\mathit{N}+\mathit{I}}\hfill \\ \hfill & -sin{u}_{1}cos{u}_{2}cos{\theta }_{\mathit{N}+\mathit{I}}\mathit{a}\frac{sin{\theta }_{\mathit{N}}}{\sqrt{\mathit{N}}}+sin{u}_{1}sin{u}_2{cos}^2{\theta }_{\mathit{N}+\mathit{I}}.\hfill \end{array}$$

In this scenario, we required the number of qubits to be even. Assuming that $\theta \ll 1$ and $\gamma t_r\ll 1$ (in practice, $\gamma t_r$ is supposed to be of the $o\left({\theta }^3\right)$ order), we expanded this Kraus map to the second order in $\theta$, and first order in $\gamma t_r$; then, the Kraus map can be approximated with the following Lindblad master equation:

$$\begin{array}{cc}\hfill \frac{\mathrm{d}\rho \left(t\right)}{\mathrm{d}t}=& -i[\mathit{H},\rho \left(t\right)]+\sum _{j=1}^3\mathcal{L}\left({\mathit{L}}_{\mathit{j}}\right)\rho \left(t\right)\hfill \\ \hfill & +2\gamma t_r\mathcal{L}\left(\mathit{a}\right)({\mathit{M}}_{gg}\rho \left(t\right){\mathit{M}}_{gg}^{†}+{\mathit{M}}_{ge}\rho \left(t\right){\mathit{M}}_{ge}^{†}\hfill \end{array}$$

$$\begin{array}{c}\hfill +{\mathit{M}}_{eg}\rho \left(t\right){\mathit{M}}_{eg}^{†}+{\mathit{M}}_{ee}\rho \left(t\right){\mathit{M}}_{ee}^{†}).\end{array}$$

(19)

Equation (19) stabilizes a squeezed state $|\alpha ,r\rangle$ with

$$\begin{array}{cc}\hfill \alpha & =\frac{tan(u_1+u_2)}{\theta +\gamma t_r},\hfill \end{array}$$

(20a)

$$\begin{array}{cc}\hfill r& =\frac{1}{4}ln\frac{{\theta }^2{cos}^2(u_1-u_2)+\gamma t_r}{{\theta }^2{cos}^2(u_1+u_2)+\gamma t_r},\hfill \end{array}$$

(20b)

Provided that $sin(u_1+u_2)$ is sufficiently small. This can be proven by applying the displacement and squeezing transformations to transform the original density operator $\rho \left(t\right)$ into $\tilde{\rho }\left(t\right)$ as follows:

$$\begin{array}{c}\hfill \tilde{\rho }\left(t\right)={\mathit{S}}^{†}\left(r\right){\mathit{D}}^{†}\left(\alpha \right)\rho \left(t\right)\mathit{D}\left(\alpha \right)\mathit{S}\left(r\right),\end{array}$$

whose evolution can be described by

$$\begin{array}{cc}\hfill \frac{\mathrm{d}\tilde{\rho }\left(t\right)}{\mathrm{d}t}=& [2{\theta }^2{(cos{u}_{1}cos{u}_{2}coshr+sin{u}_{1}sin{u}_{2}sinhr)}^2+2\gamma t_r{cosh}^{2}r]\mathcal{L}\left(\mathit{a}\right)\tilde{\rho }\left(t\right)\hfill \\ \hfill & +[2{\theta }^2{(cos{u}_{1}cos{u}_{2}coshr+sin{u}_{1}sin{u}_{2}sinhr)}^2+2\gamma t_r{sinh}^{2}r]\mathcal{L}\left({\mathit{a}}^{†}\right)\tilde{\rho }\left(t\right)\hfill \\ \hfill & +\theta l(sinhr+coshr)[{\mathit{a}}^{†},\tilde{\rho }\left(t\right)]+\theta l(sinhr+coshr)[\tilde{\rho }\left(t\right),\mathit{a}]\hfill \\ \hfill & +{\theta }^2(cos{u}_{1}cos{u}_2\alpha -sin{u}_{1}sin{u}_2{\alpha }^{*})\hfill \\ \hfill & \times (cos{u}_{1}cos{u}_{2}coshr+sin{u}_{1}sin{u}_{2}sinhr)[\tilde{\rho }\left(t\right),{\mathit{a}}^{†}]\hfill \\ \hfill & +{\theta }^2(sin{u}_{1}sin{u}_2\alpha -cos{u}_{1}cos{u}_2{\alpha }^{*})\hfill \\ \hfill & \times (cos{u}_{1}cos{u}_{2}sinhr+sin{u}_{1}sin{u}_{2}coshr)[{\mathit{a}}^{†},\tilde{\rho }\left(t\right)]\hfill \\ \hfill & +{\theta }^2(cos{u}_{1}cos{u}_2{\alpha }^{*}-sin{u}_{1}sin{u}_2\alpha )\hfill \\ \hfill & \times (cos{u}_{1}cos{u}_{2}coshr+sin{u}_{1}sin{u}_{2}sinhr)[\mathit{a},\tilde{\rho }\left(t\right)]\hfill \\ \hfill & +{\theta }^2(sin{u}_{1}sin{u}_2{\alpha }^{*}-cos{u}_{1}cos{u}_2\alpha )\hfill \\ \hfill & \times (cos{u}_{1}cos{u}_{2}sinhr+sin{u}_{1}sin{u}_{2}coshr)[\tilde{\rho }\left(t\right),\mathit{a}]\hfill \\ \hfill & +\gamma t_r(sinhr{\alpha }^{*}+coshr\alpha )[\tilde{\rho }\left(t\right),{\mathit{a}}^{†}]-\gamma t_r(sinhr\alpha +coshr{\alpha }^{*})[\tilde{\rho }\left(t\right),\mathit{a}]\hfill \\ \hfill & +{\theta }^2(cos{u}_{1}cos{u}_{2}coshr+sin{u}_{1}sin{u}_{2}sinhr)\hfill \\ \hfill & \times (cos{u}_{1}cos{u}_{2}sinhr+sin{u}_{1}sin{u}_{2}coshr)([[\tilde{\rho }\left(t\right),{\mathit{a}}^{†}],\mathit{a}]+[[\tilde{\rho }\left(t\right),\mathit{a}],\mathit{a}])\hfill \\ \hfill & +\gamma t_{r}sinhrcoshr([[\tilde{\rho }\left(t\right),{\mathit{a}}^{†}],\mathit{a}]+[[\tilde{\rho }\left(t\right),\mathit{a}],\mathit{a}]).\hfill \end{array}$$

(21)

Therefore, if we chose $4r=ln\frac{{\theta }^2{cos}^2(u_1-u_2)+\gamma t_r}{{\theta }^2{cos}^2(u_1+u_2)+\gamma t_r}$ and $\alpha =\frac{tan(u_1+u_2)}{\theta +\gamma t_r}$, Equation (22) could be equivalently written as follows:

$$\begin{array}{c}\hfill \frac{\mathrm{d}\tilde{\rho }\left(t\right)}{\mathrm{d}t}={ϵ}_{ds}\mathcal{L}\left(\mathit{a}\right)\tilde{\rho }\left(t\right),\end{array}$$

(22)

where

$$\begin{array}{cc}\hfill {ϵ}_{ds}=& \frac{1}{2}({\theta }^{2}co{s}^2(u_1-u_2)+\gamma t_r)\sqrt{\frac{{\theta }^2{cos}^2(u_1-u_2)+\gamma t_r}{{\theta }^2{cos}^2(u_1+u_2)+\gamma t_r}}\hfill \\ \hfill & +\frac{1}{2}({\theta }^{2}co{s}^2(u_1+u_2)+\gamma t_r)\sqrt{\frac{{\theta }^2{cos}^2(u_1+u_2)+\gamma t_r}{{\theta }^2{cos}^2(u_1-u_2)+\gamma t_r}}\hfill \\ \hfill & +{\theta }^{2}cos(u_1-u_2)cos(u_1+u_2)+\gamma t_r.\hfill \end{array}$$

Convergence towards a squeezed state of the oscillator in such a scenario also requires $(1-2\gamma t_r)cos(u_1-u_2)cos(u_1+u_2)\ge 0$. This can be viewed by taking into consideration Lyapunov equation $V\left(t\right)=Tr\left({\mathit{a}}^{†}\mathit{a}{\mathit{S}}^{†}\left(r\right){\mathit{D}}^{†}\left(\alpha \right)\rho \left(t\right)\mathit{D}\left(\alpha \right)\mathit{S}\left(r\right)\right)$. Then, it is not difficult to find that

$$\begin{array}{c}\hfill \frac{d}{dt}V\left(t\right)=-2{\theta }^2(1-2\gamma t_r)cos(u_1-u_2)cos(u_1+u_2)V\left(t\right)=-{\kappa }_{ds}V\left(t\right),\end{array}$$

(23)

which should make $\frac{d}{dt}V\left(t\right)\le 0$ hold; thus, ${\kappa }_{ds}=2{\theta }^2(1-2\gamma t_r)cos(u_1-u_2)cos(u_1+u_2)$ is the corresponding convergence rate.

In a similar way, to calculate the QFI as presented in Section 2, we could obtain the QFI regarding a squeezed state when photonic loss from the oscillator to the environment is taken into account as follows (subscripts d and s stand for “dissipation” and “squeezed”, respectively):

$$\begin{array}{cc}\hfill {\mathcal{F}}_{ds}=& 4{\kappa }_{ds}^2{\left|\frac{d\alpha }{d\theta }\right|}^2(cosh2r+sinh2r),\hfill \\ \hfill =& 16{\theta }^4\frac{{sin}^2{u}_s}{{(\theta +\gamma t_r)}^4}\sqrt{\frac{{\theta }^2{cos}^2{u}_d+\gamma t_r}{{\theta }^2{cos}^2{u}_s+\gamma t_r}}{(1-2\gamma t_r)}^2{cos}^2{u}_d,\hfill \end{array}$$

(24)

where $u_d=u_1-u_2$ and $u_s=u_1+u_2$. Apparently, the QFI ${\mathcal{F}}_{ds}$ given in Equation (24) could degenerate into the QFI ${\mathcal{F}}_s$ given in Equation (12) when $\gamma =0$. This, in turn, verified the consistency and correctness of our derivation.

3.3. Analysis of Enhanced Estimation Precision via the Time-Varying Reservoir

In this subsection, we show the significant increase in QFI in the presence of photonic loss from the oscillator to the environment while the oscillator is converging into a squeezed state instead of a coherent state.

Referring to the physical parameters in [22], we set dissipation rate $\gamma =$ 20 Hz, and the interaction time $t_r$ to $1.25\times {10}^{-2}$ ms, with the coupling strength mostly chosen to be $\theta =\frac{\pi }{8}$ in simulations. For example, if u was set to $u=0.1$ when all the reservoir qubits were initialized at the same state $|\psi \rangle =cos0.1|g\rangle +sin0.1|e\rangle$, the final state of the oscillator could be visualized through the Wigner distribution as depicted in Figure 3a. The state in this case was coherent.

By contrast, we chose $u_1=0.395$ and $u_2=0.39$, as shown in Figure 3b. Namely, the reservoir qubits were initially prepared at states $|{\psi }_{qred}\rangle =cos0.395|g\rangle +sin0.395|e\rangle$ and $|{\psi }_{qblue}\rangle =cos0.39|g\rangle +sin0.39|e\rangle$ pair by pair. Specifically, the Wigner distribution and the corresponding x–p phase for the final state of the oscillator are shown in Figure 3b, which shows that the state was squeezed. The Wigner distributions depicted in Figure 3 indicate that the form of steady states numerically obtained with the discrete-time Kraus map was consistent with the form of steady states theoretically derived on the basis of the approximated continuous-time Lindblad master equation.

In Figure 4, we plotted the values of QFI, as the initial states of qubits varied. For example, the blue dotted line in Figure 4b corresponds to the scenario where the oscillator was moving towards a coherent state in our simulation in the presence of photonic loss, with dissipation rate $\gamma =20$ Hz. It is obvious that ${\mathcal{F}}_{dc}\le 5.1$.

In fact, according to Equation (18), the maximal ${\mathcal{F}}_{dc}$ could be obtained theoretically. By calculating the derivative of ${\mathcal{F}}_{dc}$ with respect to $cos2u$ ($u\in (0,u_{max})$, $u_{max}<\frac{\pi }{4}$), we have that

$$\begin{array}{c}\hfill \frac{\partial {\mathcal{F}}_{dc}}{\partial (cos2u)}=-32\frac{(\frac{\gamma t_r}{{\theta }^2}-cos2u)[2\frac{\gamma t_r}{{\theta }^2}(1-cos{2u}^2)+cos2u(\frac{{\gamma }^2{t}_r^2}{{\theta }^4}-{cos}^{2}2u)]}{{(\frac{\gamma t_r}{{\theta }^2}+cos2u)}^3}.\end{array}$$

The stability condition also requires u to be sufficiently small. According to both numerical and theoretical analyses, $u_{max}\approx \frac{\pi }{10}$. Then, due to $\frac{\gamma t_r}{{\theta }^2}\ll 1$ and $0<\frac{\gamma t_r}{{\theta }^2}<\frac{{cos}^{2}2u-1+\sqrt{2{cos}^{4}2u-2{cos}^{2}2u+1}}{cos2u}$, we have $\frac{\partial {\mathcal{F}}_{dc}}{\partial (cos2u)}<0$ in this region, which indicates that ${\mathcal{F}}_{dc}$ increased when u increased ($cos2u$ decreased). Therefore, $max\left({\mathcal{F}}_{dc}\right)\approx 16(1-{cos}^2(2\times \frac{\pi }{10}))=5.53$.

The maximal ${\mathcal{F}}_{ds}$ can be obtained in a similar way. By calculating the derivative of ${\mathcal{F}}_{ds}$ with respect to $u_s$, i.e.,

$$\begin{array}{cc}\hfill \frac{\partial {\mathcal{F}}_{ds}}{\partial u_s}=& 16\frac{{\theta }^4{cos}^2{u}_d}{{(\theta +\gamma t_r)}^4}{(1-2\gamma t_r)}^2\left[2sin{u}_{s}cos{u}_s\sqrt{\frac{{\theta }^2{cos}^2{u}_d+\gamma t_r}{{\theta }^2{cos}^2{u}_s+\gamma t_r}}\right.\hfill \\ \hfill +& \left.\frac{{sin}^2{u}_s}{2\sqrt{\frac{{\theta }^2{cos}^2{u}_d+\gamma t_r}{{\theta }^2{cos}^2{u}_s+\gamma t_r}}}\frac{{\theta }^2{cos}^2{u}_d+\gamma t_r}{{({\theta }^2{cos}^2{u}_s+\gamma t_r)}^2}\right],\hfill \end{array}$$

For $\frac{\gamma t_r}{{\theta }^2}\ll 1$ and $u_s\in (0,\frac{\pi }{4}]$, $\frac{\partial {\mathcal{F}}_{ds}}{\partial u_s}>0$ thus holds in this region. Additionally, $u_d$ should be reduced in order to gain a larger value for ${\mathcal{F}}_{ds}$. Hence, for $u_s\in (0,u_{s_{max}}]$ ($u_{s_{max}}$ was around $\frac{\pi }{4}$ on the basis of our numerical and theoretical analysis), we have that $max{\mathcal{F}}_{dc}\approx 16{sin}^2\left(\frac{\pi }{4}\right)/cos\left(\frac{\pi }{4}\right)=11.31$.

In contrast to the blue dotted line in Figure 4b, the red and green solid lines show the values of QFI concerning the steady state of the oscillator when the oscillator was stabilized at a squeezed state. Specifically, the red line corresponds to the case where the reservoir qubits were initialized in pairs with $u_1=0.4$, as $u_2$ varied from 0 to 0.3, while the green line corresponds to the case where the reservoir qubits were initialized in pairs with $u_1=0.3$, as $u_2$ varied from 0 to 0.3. In both cases, larger values of QFI were obtained compared to the case where all the qubits had been initialized at the same state (u varied from 0 to 0.3). Therefore, steering the oscillator towards a squeezed state enabled us to obtain a much larger value of QFI. If we chose $u_1=0.395$, and $u_2=0.39$, then ${\mathcal{F}}_{ds}=11.25$ was almost double $max{\mathcal{F}}_{dc}$. In addition, the red line presents larger values of QFI than the green line does. This is consistent with the theoretical analysis above that QFI increased as $u_s$ increased when the oscillator was converging towards a squeezed state.

We also took into account the cases of $\gamma =200$ Hz and $\gamma =2$ KHz for comparison, with $\gamma$ = 0 Hz included in Figure 4 in order to see the effect of dissipation. In particular, the cases of $\gamma =0$ Hz, $\gamma =200$ Hz and $\gamma =2$ KHz are demonstrated in Figure 4a,c,d, respectively. Apparently, when $\gamma$ grew, the value of QFI was reduced. It is not surprising that the values of QFI were almost the same for the two cases of $\gamma$ = 20 Hz and $\gamma$ = 0 Hz because, in principle, the steady state of the harmonic oscillator in the presence of photonic loss was not very different from the steady state in the ideal case when $\gamma$ was sufficiently small.

On the other hand, it is also very important to estimate the resources that we must utilize. First, the convergence rates for the oscillator in the cases of stabilization at a coherent state and a squeezed state were theoretically compared as follows. Bearing in mind the assumption that $\frac{\gamma t_r}{{\theta }^2}\ll 1$, one has that

$$\begin{array}{cc}\hfill \frac{{\kappa }_{ds}}{{\kappa }_{dc}}& =\frac{\left(1-2\gamma t_r\right)cos\left(u_1-u_2\right)cos\left(u_1+u_2\right)}{cos2u+\frac{\gamma t_r}{{\theta }^2}}\hfill \\ \hfill & =\frac{\left(1-2\gamma t_r\right)cos\left(u_1-u_2\right)cos\left(u_1+u_2\right)}{2{cos}^{2}u-\left(1-\frac{\gamma t_r}{{\theta }^2}\right)}\hfill \\ \hfill & \approx \frac{{cos}^2{u}_1+{cos}^2{u}_2-1}{2{cos}^{2}u-1}.\hfill \end{array}$$

(25)

Therefore, the convergence rates in these two cases are typically of the same speed order. For example, if we chose $u_1<u<\frac{\pi }{4}$ and $u_2<u<\frac{\pi }{4}$, ${\kappa }_{ds}$ was rendered to be larger than ${\kappa }_{dc}$, which indicates that fewer reservoir qubits were needed to harvest a squeezed state of the oscillator under such circumstances.

In Figure 5, we demonstrate the number of consumed reservoir qubits to reach the steady state of the oscillator. Here, $2n$ equals the number of qubits, and fidelity represents the similarity between the current and steady states of the oscillator. Fidelity $\mathrm{F}$ can be calculated as follows [32]:

$$\begin{array}{c}\hfill \mathrm{F}(\rho \left(t\right),\rho (\infty ))=\mathrm{Tr}\left(\sqrt{\sqrt{\rho \left(t\right)}\rho (\infty )\sqrt{\rho \left(t\right)}}\right).\end{array}$$

(26)

In more concrete terms, as shown in Figure 5a, with $u_1=0.18$ and $u_2=0.12$, for $\theta =\frac{\pi }{12}$, $\theta =\frac{\pi }{10}$, and $\theta =\frac{\pi }{8}$, we plot how fidelity ${\mathrm{F}}_{\mathrm{s}}$ varied as the number of qubits increased when the oscillator was converging towards a squeezed state using the red solid, blue dashed, and green dotted lines respectively. Additionally, Figure 5b shows, with $u=0.2$, the evolution of ${\mathrm{F}}_{\mathrm{c}}$ when the oscillator was converging towards a coherent state, where the red solid, blue dashed, and green dotted lines correspond to the cases where $\theta =\frac{\pi }{12}$, $\theta =\frac{\pi }{10}$ and $\theta =\frac{\pi }{8}$, respectively. Typically, fewer than 120 qubits (60 pairs of qubits) were needed to reach the steady state. The smaller the $\theta$ was, the more reservoir qubits were required to obtain the steady state, which could give increase estimation precision.

In our system setup, the oscillator could be stabilized at a squeezed state using the “time-varying” reservoir as specified in Section 3. Since the final state was not temporarily generated, it was easier and more convenient for us to impose or repeat the measurement of such a steady squeezed state in practice.

4. Discussion

In Section 3, we allowed dissipation through the oscillator to the environment. Our theoretical analysis shows that the steady state of the harmonic oscillator in the presence of photonic loss should not be very different from the steady state in the ideal case under the assumption of $\frac{\gamma t_r}{{\theta }^2}\ll 1$. This could also be observed by comparing the simulation results in Figure 4a,b. From another point of view, our system setup that could stabilize the oscillator state presented robustness to photonic loss.

On the other hand, the initial states of the reservoir qubits were chosen depending on the problems that needed explaining in our numerical simulation. In Figure 4, we focus on the comparative analysis regarding the values of QFI between different cases where the oscillator was steered towards a coherent or squeezed state. Therefore, $u_2$ varies from 0 to 0.3 in Figure 4 because u had to take values from 0 to 0.3 to ensure the convergence of the oscillator state. Squeezed states led to larger values of QFI when $u_2$ and u varied in the same range. In Figure 3b, a steady squeezed state of the oscillator is depicted. Parameters $u_1=0.395$ and $u_2=0.39$ were, thus, chosen in this case in order to render the squeezing obvious according to Equation (20). In Figure 5, we chose $u_1=0.18$ and $u_2=0.12$ ($u=0.2$) since the purpose was to show that, by appropriately preparing the initial states of qubits (i.e. $u_1<u=0.2<\frac{\pi }{4}$ and $u_2<u=0.2<\frac{\pi }{4}$), convergence rate ${\kappa }_{ds}$ when the oscillator was steered towards a squeezed state could be made larger than convergence rate ${\kappa }_{dc}$ when the oscillator was steered towards a coherent state according to Equation (26). There was, indeed, a trade-off between fast convergence and strong squeezing strength (large value of QFI).

In the future, we will consider variants in the system setup by adjusting the initial states of reservoir qubits, and the interaction between each qubit and the oscillator with the purpose of seeking other applications in quantum computation and quantum metrology. In terms of other imperfections in practice, dephasing and relaxation errors playing the role of noise to the preparation of qubits states, along with missing qubits, may also be considered.

5. Conclusions

Mainly inspired by [20] where the Serge Haroche experimental setting was utilized to restore the Heisenberg limit regarding estimating the coupling strength between each qubit and the harmonic oscillator, in this paper, we altered the setting by incorporating a time-varying quantum reservoir that enabled us to obtain even higher estimation precision. The contributions of this paper are threefold. First, in contrast to the results in [23], we took into account photonic loss through the harmonic oscillator to the environment and analytically provided the steady state of the oscillator under such circumstances. With consecutive pairs of separable time-varying input qubits, we could stabilize the oscillator at a squeezed state. The corresponding convergence condition and convergence rate were given as well. Second, we proved that stabilizing the oscillator at a squeezed state could improve the estimation precision of qubit–oscillator coupling strength by theoretically calculating the QFI, which was demonstrated in and verified by simulation examples in comparison with the case where all the qubits were prepared at an identical state, as discussed in the previous work. Last but not least, we discuss the approximate number of reservoir qubits that needs to be utilized for the oscillator to reach the final state. By appropriately choosing the initial qubit states, the time-varying reservoir could help us in achieving faster convergence. Since this final state of the oscillator could stay steady by absorbing qubits from the reservoir, it is more convenient for implementing quantum metrology experimentally.