Effects of counter-rotating-wave terms on the noisy frequency estimation

Quantum parameter estimation is a rapidly developing research field, which has many potential applications from gravitational wave detection Tse et al. (2019); Acernese et al. (2019), quantum radar Maccone and Ren (2020); Huang et al. (2021), quantum illumination Tan et al. (2008); Xu et al. (2021) to various ultrasensitive quantum thermometries Correa et al. (2015); Hovhannisyan et al. (2021); Candeloro et al. (2021) and magnetometers Bouton et al. (2020); Razzoli et al. (2019). As demonstrated in many previous theoretical and experimental studies, using quantum coherence Smirne et al. (2019), quantum entanglement Huelga et al. (1997); Nagata et al. (2007); Pezzè et al. (2018), quantum squeezing Caves (1981); Bai et al. (2019); Chalopin et al. (2018) and quantum criticality Invernizzi et al. (2008); Wang et al. (2014); Frérot and Roscilde (2018); Chu et al. (2021), the performance of quantum metrology can surpass the so-called shot-noise limit or standard quantum limit, which is usually achieved in classical metrological schemes.

Unfortunately, the superiority of quantum metrology is commonly destroyed by decoherence Bai et al. (2019); Chin et al. (2012); Zwick et al. (2016a, b); Matsuzaki et al. (2018); Haase et al. (2018); Tamascelli et al. (2020); Wu and Shi (2020, 2021), which is induced by the inevitable interaction with the surrounding environment. For example, Ref. Bai et al. (2019) reported that the Zeno-limit-type scaling relation, which is generated by the resource of quantum squeezing in the noiseless case, is degraded into the shot-noise limit under the influence of decoherence. In this sense, in any practical scheme, the decoherence should be carefully taken into account, which gives rise to the development of the so-called noisy quantum metrology Haase et al. (2018).

However, almost all the present studies of noisy quantum metrology restricted their attentions to the bosonic environment situation with the pure dephasing mechanism Chin et al. (2012); Matsuzaki et al. (2018); Rossi and Paris (2015); Tatsuta et al. (2019) or the rotating-wave approximation (RWA) Berrada (2013); Li et al. (2017); Wu et al. (2021). On the other hand, the counter-rotating-wave terms play a significant role in the atomic spontaneous emission Yang et al. (2013), the quantum Zeno and anti-Zeno effect Zheng et al. (2008); Ai et al. (2010) as well as the non-Markovianity in open quantum systems Mäkelä and Möttönen (2013); Wu and Liu (2017). A question naturally raised here: what is the influence of the counter-rotating-wave terms on the estimation precision of a noisy quantum metrology? To address the above concern, with the help of the non-perturbative hierarchical equations of motion (HEOM) method Tanimura (2020); an Yan et al. (2004); Xu and Yan (2007); Jin et al. (2007, 2008); Wu (2018a, b), we go beyond the usual RWA treatment and investigate the performance of estimating the frequency of a two-level system which interacts with a bosonic or a fermionic environment.

This paper is organized as follows. In Sec. II, we first recall some basic concepts, mainly about the classical and the quantum Fisher information, in the quantum parameter estimation theory. In Sec. III, we propose our noisy quantum metrology scheme and outline our methodology. The effect of the counter-rotating-wave terms on our metrological precision is analyzed in Sec. IV. Some concerned discussions and the main conclusions of this paper are drawn in Sec. V. In several appendixes, we provide some additional details about the HEOM method as well as the other two additional dynamical approaches. Throughout the paper, we set $\hbar=k_{\mathrm{B}}=1$.

Generally, in a typical quantum metrology scheme, the parameter of interest, say $\theta$ in this section, is encoded into the state of a quantum system (acting as a probe) via a unitary or a nonunitary dynamical process. Then, the information about $\theta$ can be extracted from the $\theta$-dependent state $\varrho_{\theta}$ via repeated quantum measurements. In the above scenario, the metrological precision with respect to a given measurement scheme is constrained by the famous Cram$\mathrm{\acute{e}}$r-Rao bound Kay (1993); Braunstein and Caves (1994); Zhang and Yang (2018)

$$\delta^{2}\theta\geq\frac{1}{\upsilon\mathcal{F}_{\text{C}}},$$

(1)

where $\delta^{2}\theta$ denotes the variance of the derived estimator, $\upsilon$ is the number of repeated measurements (we set $\upsilon=1$ for the sake of consentience in this paper) and $\mathcal{F}_{\text{C}}$ is the classical Fisher information (CFI) corresponding to the selected measurement scheme. To perform such a measurement in the theory of quantum mechanism, one needs to construct a set of positive-operator valued measurement operators $\{\hat{\Pi}_{u}\}$, which satisfy $\sum_{u}\hat{\Pi}_{u}^{\dagger}\hat{\Pi}_{u}=\mathbf{\hat{1}}$ with discrete measurement outcomes $\{u\}$. For an arbitrary $\theta$-dependent state $\varrho_{\mathrm{\theta}}$, the measurement operator $\hat{\Pi}_{u}$ yields an outcome $u$ with a corresponding probability distribution $p(u|\theta)=\mathrm{Tr}(\hat{\Pi}_{u}\varrho_{\theta}\hat{\Pi}_{u}^{\dagger})$. With all the probabilities from $\{\hat{\Pi}_{u}\}$ at hand, the CFI can be computed via Kay (1993); Braunstein and Caves (1994); Zhang and Yang (2018)

$$\mathcal{F}_{\text{C}}=\sum_{u}p(u|\theta)\bigg{[}\frac{\partial}{\partial\theta}\ln p(u|\theta)\bigg{]}^{2}.$$

(2)

From Eq. (2), one can find the value of CFI strongly relies on the choice of measurement operators. Running over all the possible measurement schemes, one can prove that the ultimate metrological precision is bounded by the following quantum Cram$\mathrm{\acute{e}}$r-Rao inequality Liu et al. (2019)

$$\delta^{2}\theta\geq\frac{1}{\upsilon\mathcal{F}_{\text{Q}}},$$

(3)

where $\mathcal{F}_{\text{Q}}\equiv\mathrm{Tr}(\hat{\varsigma}^{2}\varrho_{\theta})$ with $\hat{\varsigma}$ determined by $\partial_{\theta}\varrho_{\theta}=\frac{1}{2}(\hat{\varsigma}\varrho_{\theta}+\varrho_{\theta}\hat{\varsigma})$ is the so-called quantum Fisher information (QFI). Specially, if $\varrho_{\theta}$ is a two-dimensional density matrix in the Bloch representation, namely, $\varrho_{\theta}=\frac{1}{2}(\mathbf{1}_{2}+\vec{r}\cdot\vec{\sigma})$ with $\vec{r}$ being the Bloch vector and $\vec{\sigma}\equiv(\hat{\sigma}_{x},\hat{\sigma}_{y},\hat{\sigma}_{z})$ being the vector of Pauli matrices, the QFI can be easily calculated via Liu et al. (2019); Zhong et al. (2013)

$$\mathcal{F}_{\text{Q}}=|\partial_{\theta}\vec{r}|^{2}+\frac{(\vec{r}\cdot\partial_{\theta}\vec{r})^{2}}{1-|\vec{r}|^{2}}.$$

(4)

If $\varrho_{\theta}$ is a pure state, Eq. (4) reduces to $\mathcal{F}_{\text{Q}}=|\partial_{\theta}\vec{r}|^{2}$.

Physically speaking, the QFI describes the most statistical information in $\varrho_{\theta}$, while the CFI denotes the statistical message extracted from the selected measurement scheme. Thus, one can immediately conclude that $\mathcal{F}_{\text{Q}}\geq\mathcal{F}_{\text{C}}$. When the selected measurement scheme is the optimal one, $\mathcal{F}_{\text{C}}$ can be saturated to $\mathcal{F}_{\text{Q}}$. Unfortunately, there is no general way to find the optimal measurement operator, which has an explicit physical meaning. In this sense, designing the physically optimal measurement scheme, which can saturate the best attainable precision bounded by the QFI, is of importance in study of quantum metrology.

As discussed in several previous articles Chin et al. (2012); Matsuzaki et al. (2018); Tatsuta et al. (2019); Berrada (2013); Wang et al. (2017), in the ideal noiseless case, the information of the estimated frequency, denoted by $\Delta$ in this paper, can be encoded to the state of a two-level atomic system, which acts as a quantum probe, via a unitary dynamics: $\varrho_{\Delta}=\varrho(t)=e^{-it\hat{H}_{\mathrm{s}}}\varrho(0)e^{it\hat{H}_{\mathrm{s}}}$ with $\hat{H}_{\text{s}}=\frac{1}{2}\Delta\hat{\sigma}_{x}$ being the Hamiltonian of the probe. However, such a unitary encoding scheme breaks down if the environmental influence is taken into account. In the noisy case, the encoding process is realized by a non-unitary dynamics $\varrho_{\Delta}=\varrho_{\text{s}}(t)=\text{Tr}_{\text{e}}[e^{-it\hat{H}}\varrho_{\text{se}}(0)e^{it\hat{H}}]$, where $\varrho_{\text{s}}(t)$ denotes the reduced density operator of the probe by tracing out the environmental degrees of freedom and $\hat{H}$ is the total Hamiltonian of the probe plus the environment.

In this paper, we assume $\hat{H}$ has the following general form Wu and Liu (2017); Wang et al. (2019); Zhang and Xu (2019); Wang et al. (2020)

$$\begin{split}\hat{H}=&\hat{H}_{\text{s}}+\sum_{k}\omega_{k}\hat{c}_{k}^{{\dagger}}\hat{c}_{k}+\sum_{k}\Big{(}g_{k}^{*}\hat{\sigma}_{-}\hat{c}_{k}^{{\dagger}}+g_{k}\hat{\sigma}_{+}\hat{c}_{k}\Big{)}\\ &+\chi\sum_{k}\Big{(}g_{k}^{*}\hat{\sigma}_{+}\hat{c}_{k}^{{\dagger}}+g_{k}\hat{\sigma}_{-}\hat{c}_{k}\Big{)},\end{split}$$

(5)

where $\hat{c}_{k}$ and $\hat{c}_{k}^{\dagger}$ are the annihilation and creation operators of the $k$th environmental mode with frequency $\omega_{k}$, respectively. If $\hat{c}_{k}$ and $\hat{c}_{k}^{\dagger}$ satisfy the canonical commutation relations $[\hat{c}_{k},\hat{c}_{k^{\prime}}^{\dagger}]=\delta_{kk^{\prime}}$, one can deem that the environment is composed of bosons. On the other hand, if $\hat{c}_{k}$ and $\hat{c}_{k}^{\dagger}$ obey $\{\hat{c}_{k},\hat{c}_{k^{\prime}}^{\dagger}\}=\delta_{kk^{\prime}}$, then the environment is a fermionic environment. Operators $\hat{\sigma}_{\pm}$ are defined as $\hat{\sigma}_{-}=\hat{\sigma}_{+}^{{\dagger}}\equiv|-\rangle\langle+|$ with $|\pm\rangle$ being the eigenstates of Pauli-$x$ operator, i.e. $\hat{\sigma}_{x}|\pm\rangle=\pm|\pm\rangle$, and $g_{k}$ are complex numbers quantifying the coupling strength between the quantum probe and the $k$th environmental mode. The dimensionless parameter $\chi$ is a real number, satisfying $\chi\in[0,1]$, characterizes the weight of counter-rotating-wave terms in the interaction Hamiltonian. When $\chi=0$, all the counter-rotating-wave terms are removed and $\hat{H}$ reduces to the damped Jaynes-Cummings model Berrada (2013); Li et al. (2017); Wu et al. (2021); Wang et al. (2017). When $\chi=1$, all the contributions of the counter-rotating-wave terms are taken into consideration and this Hamiltonian is totally beyond the RWA. By introducing such a parameter, a bridge between the RWA regime and the non-RWA regime can be built, which helps us to get a deeper understanding of the effect of counter-rotating-wave terms on the parameter estimation in a dissipative environment.

To obtain the information of the estimated frequency, one need to monitor the reduced density operator of the probe, namely $\varrho_{\text{s}}(t)$. In this paper, we adopt the HEOM method Tanimura (2020); an Yan et al. (2004); Xu and Yan (2007); Jin et al. (2007, 2008); Wu (2018a, b), which can provide a numerically rigours result, to solve this problem (see Appendix A for more details). As benchmarks, we also employ the other two approaches, the Zwanzig-Nakajima master equation technique (see Refs. DiVincenzo and Loss (2005); Burkard (2009) and Appendix B) and the RWA treatment (see Appendix C), to investigate the reduced dynamics of the probe. Together with the purely numerical results from the HEOM, these two additional methods can make our conclusion more persuasive.

In this paper, we assume the initial state of the whole system is given by

$$\varrho_{\text{se}}(0)=\varrho_{\text{s}}(0)\otimes\varrho_{\text{e}}=|\psi_{\text{s}}(0)\rangle\langle\psi_{\text{s}}(0)|\bigotimes_{k}|0_{k}\rangle\langle 0_{k}|,$$

(6)

where $|\psi_{\text{s}}(0)\rangle=\cos\phi|+\rangle+\sin\phi|-\rangle$, and $|0_{k}\rangle$ denotes the vacuum state of the $k$th environmental mode. The spectral density of the environment $J(\omega)$, which is defined by $J(\omega)=\sum_{k}|g_{k}|^{2}\delta(\omega-\omega_{k})$, has a Lorentzian form

$$J(\omega)=\frac{1}{2\pi}\frac{\gamma\lambda^{2}}{(\omega-\Delta)^{2}+\lambda^{2}},$$

(7)

where $\lambda$ defines the spectral width, $\gamma$ can be approximately interpreted as the probe-environemt coupling strength. The Lorentzian spectral density has a clear Markovian-non-Markovian boundary Breuer and Petruccione (Oxford University Press, Oxford, 2002); Bellomo et al. (2007); Wu and Lin (2017); Wu and Shi (2021) and can give rise to an Ornstein-Uhlenbeck-type environmental correlation function, which guarantees the feasibility of the HEOM algorithm (see Appendix A for details).

As long as $\varrho_{\text{s}}(t)$ is obtained, the CFI and the QFI can be calculated by making use of Eq. (2) and Eq. (4), respectively. During the purely numerical calculations, one needs to handle the first-order derivative to the parameter $\Delta$, say $\partial_{\Delta}\vec{r}$ in Eq. (4). In this paper, the first-order derivative for an arbitrary $\Delta$-dependent function $f_{\Delta}$ is numerically done by adopting the following finite difference method

$$\partial_{\Delta}f_{\Delta}\simeq\frac{-f_{\Delta+2\delta}+8f_{\Delta+\delta}-8f_{\Delta-\delta}+f_{\Delta-2\delta}}{12\delta}.$$

(8)

We set $\delta/\Delta=10^{-6}$, which provides a very good accuracy.

In this section, we try to address the effect of counter-rotating-wave terms on the noisy quantum frequency estimation. To obtain the CFI, the measurement operators in this paper are chosen as $\hat{\Pi}=\{|e\rangle\langle e|,|g\rangle\langle g|\}$ with $|e,g\rangle$ being the eigenstates of Pauli-$z$ operator. In the ideal noiseless case, the CFI with respect to the selected measurement scheme is $\mathcal{F}_{\text{C}}^{\text{ideal}}=t^{2}$, which is completely saturated to the QFI (see Appendix C). This result means the choice of the measurement scheme is the optimal one. Moreover, our measurement scheme is mathematically equivalent to perform a standard Ramsey spectroscopy Huelga et al. (1997); Chin et al. (2012); Wang et al. (2017), which is commonly employed to estimate the phases (or an external magnetic field) of cold atomic systems. As the detailed expositions of a typical Ramsey interferometer setup for a two-level system Huelga et al. (1997); Chin et al. (2012); Wang et al. (2017), the uncertainty of estimating the parameter $\Delta$ from the error propagation formula is in full accord with that of the CFI obtained from Eq. (2). Thus, it is quite natural to generalize the above measurement scheme to the noisy case and evaluate the environmental influences on its performance. We find, in the noisy situation, although the selected measurement scheme cannot completely saturate the ultimate bound given by the quantum Cram$\mathrm{\acute{e}}$r-Rao theorem (see the numerical results in Fig. 2 (a-c) and Fig. 3 (a-c)), the corresponding CFI is in qualitative agreement with the QFI in the non-Markovin and strong-coupling regimes.

Before concluding our work, some important remarks shall be addressed here.

(1) In our work, the strength of the counter-rotating-wave terms in the interaction Hamiltonian, namely the parameter $\chi$, is continuously tunable for the entire range of $\chi\in[0,1]$, which can build a bridge connecting two particular limits: totally with and without RWA. In fact, such a controllable strength of the counter-rotating-wave terms in the quantum Rabi model is called anisotropy Yang and Wang (2017); Baksic and Ciuti (2014); Liu et al. (2017), which leads to a much richer ground-state phase diagram Baksic and Ciuti (2014); Liu et al. (2017), and can be used to explain the anomalous Bloch-Siegert shift in the ultrastrong-coupling regime Xie et al. (2014). These previous studies of the anisotropic quantum Rabi model motivate us to explore the influence of the counter-rotating-wave terms on the metrological precision by introducing the engineered tunable parameter $\chi$.

(2) Compared with the widely acceptable interaction Hamiltonian between a two-level spin and a bosonic environment, the form of spin-fermion interaction still remains controversial. The spin-fermion model considered in this paper is the most straightforward generalization of the well-accepted spin-boson model and can be regarded as a toy model. However, it is necessary to emphasize that such a spin-fermion model is not merely of academic interest. It may be used to describe certain physical phenomena, say, the magnetism of the weak-moment heavy-fermion compound $\text{UR}_{2}\text{Si}_{2}$ Sikkema et al. (1996); Yamada et al. (2007).

(3) Though the theoretical model and the numerical tool are the same with these of Ref. Wu and Liu (2017), the centered goal of our work is utterly different from that of Ref. Wu and Liu (2017). Our paper concentrate on a concrete physical problem, while the Ref. Wu and Liu (2017) is much closer to a methodology article discussing the HEOM method. What’s more important, we here investigate an alternative topic, which means the physical conclusion in Ref. Wu and Liu (2017) cannot be straightforwardly applied to our present paper.

(4) Simulating an exact non-Markovian dissipative dynamics of an open quantum system in the strong-coupling regime is difficult. As reported in Refs. Zhang et al. (2012); González-Tudela and Galve (2019); Yang et al. (2014); Liu et al. (2016); Rançon and Bonart (2013); McCutcheon et al. (2015), the so-called negative frequency phenomenon may appear in the strong-coupling regime, which leads to isolated eigenenergies in the energy spectrum (the bound state effect) González-Tudela and Galve (2019); Yang et al. (2014); Liu et al. (2016). The appearance of negative frequencies can lead to a long-lived coherence Zhang et al. (2012); González-Tudela and Galve (2019); Liu et al. (2016) as well as the breakdown of canonical thermalization Yang et al. (2014); McCutcheon et al. (2015), which cannot be predicted by the common quantum master equation approach. However, to the best of our knowledge, such a negative frequency phenomenon usually occurs in the unbiased spin-boson model with RWA González-Tudela and Galve (2019); Liu et al. (2016) or the quantum Brownian motion (the Caldeira-Leggett model) without adding the counter term Rançon and Bonart (2013); McCutcheon et al. (2015). For the model considered in this paper (Lorentzian spin-boson model totally beyond the RWA), the negative frequency phenomenon has not been reported. Excluding the negative frequency problem, we believe the numerical correctness from the HEOM method can be fully guaranteed.

In summary, using the numerically rigorous HEOM method, we investigate the exact reduced dynamics of two-level system coupled to a dissipative bosonic or fermionic environment. With these results, we analyze the influence of counter-rotating-wave terms on the performance of a noisy parameter estimation beyond the usual weak-coupling treatment. It is revealed that the counter-rotating-wave terms can effectively enhance the metrological precision in the non-Markovian and strong-coupling regimes, whether the dissipative environment is composed of bosons or fermions. Our results may have certain applications in quantum metrology and quantum sensing.

The authors wish to thank Dr. S.-Y. Bai, Dr. C. Chen, Prof. H.-G. Luo and Prof. J.-H. An for many useful discussions. The work was supported by the National Natural Science Foundation (Grant No. 11704025 and No. 12047501).

In this appendix, we briefly sketch the HEOM method following the detailed expositions in Refs. Wu and Liu (2017); Wu (2018a, b); Suess et al. (2015). Let us consider a general quantum dissipative system whose Hamiltonian can be described by

$$\hat{H}=\hat{H}_{\text{s}}+\sum_{k}\omega_{k}\hat{c}_{k}^{{\dagger}}\hat{c}_{k}+\sum_{k}\Big{(}g_{k}^{*}\hat{L}\hat{c}_{k}^{{\dagger}}+g_{k}\hat{L}^{\dagger}\hat{c}_{k}\Big{)},$$

(A1)

where $\hat{L}$ is the quantum probe’s operator coupled to its surrounding dissipative environment. Taking $\hat{L}=\hat{\sigma}_{-}+\chi\hat{\sigma}_{+}$, Eq. (5) in the main text can be recovered. The dynamics of the Hamiltonian $\hat{H}$ is governed by the common Schr$\ddot{\mathrm{o}}$dinger equation. By introducing the bosonic or the fermionic coherent-state $|\textbf{z}\rangle\equiv\bigotimes_{k}|z_{k}\rangle$ with $\hat{c}_{k}|z_{k}\rangle=z_{k}|z_{k}\rangle$, one can recast the standard Schr$\ddot{\mathrm{o}}$dinger equation into the following non-Markovian quantum state diffusion Jing and Yu (2010); Diósi and Ferialdi (2014); Suess et al. (2014)

$$\begin{split}\frac{\partial}{\partial t}|\psi_{t}(\textbf{z}^{*})\rangle=&-i\hat{H}_{\text{s}}|\psi_{t}(\textbf{z}^{*})\rangle+\hat{L}\textbf{z}_{t}^{*}|\psi_{t}(\textbf{z}^{*})\rangle\\ &-\hat{L}^{\dagger}\int_{0}^{t}d\tau C(t-\tau)\frac{\delta}{\delta\textbf{z}_{\tau}^{*}}|\psi_{t}(\textbf{z}^{*})\rangle,\end{split}$$

(A2)

where $|\psi_{t}(\textbf{z}^{*})\rangle$ is the wave function in the coherent-state representation. The stochastic variable (random noise) $\textbf{z}_{t}\equiv i\sum_{k}g_{k}e^{-i\omega_{k}t}z_{k}$ satisfies $\mathcal{M}\{\textbf{z}_{t}\}=\mathcal{M}\{\textbf{z}_{t}^{*}\}=0$ and $\mathcal{M}\{\textbf{z}_{t}\textbf{z}_{\tau}^{*}\}=C(t-\tau)\equiv\sum_{k}|g_{k}|^{2}e^{-i\omega_{k}(t-\tau)}$, where $\mathcal{M}\{...\}$ denotes the statistical mean over all the possible stochastic processes and $C(t)$ is the environmental correction function. By far, no approximation is invoked, which means the non-Markovian quantum state diffusion equation given by Eq. (A2) is exact. However, Eq. (A2) is difficult to solve due to the time-non-local functional derivative term Wu (2018a, b).

Fortunately, for the Lorentzian spectral density considered in this paper, $C(t)$ is an Ornstein-Uhlenbeck-type correlation function, namely $C(t)$ can be expressed as an exponential function: $C(t)=\alpha e^{-\beta t}$ with $\alpha=\frac{1}{2}\gamma\lambda$ and $\beta=\lambda+i\Delta$. Using this particularity, the HEOM algorithm can be realized Wu and Liu (2017); Wu (2018a, b); Suess et al. (2015). By introducing auxiliary operators $\varrho_{t}^{(m,n)}\equiv\mathcal{M}\{|\psi_{t}^{(m)}(\textbf{z}^{*})\rangle\langle\psi_{t}^{(n)}(\textbf{z}^{*})|\}$ with

$$|\psi_{t}^{(m)}(\textbf{z}^{*})\rangle\equiv\bigg{[}\int_{0}^{t}d\tau C(t-\tau)\frac{\delta}{\delta\textbf{z}_{\tau}^{*}}\bigg{]}^{m}|\psi_{t}(\textbf{z}^{*})\rangle,$$

(A3)

the non-Markovian quantum state diffusion equation can be reexpressed as the following HEOM Wu and Liu (2017); Wu (2018a, b); Suess et al. (2015)

$$\begin{split}\frac{d}{dt}\varrho_{t}^{(m,n)}=&(-i\hat{H}_{\text{s}}^{\times}-m\beta-n\beta^{*})\varrho_{t}^{(m,n)}\\ &+m\alpha\hat{L}\varrho_{t}^{(m-1,n)}+n\alpha^{*}\varrho_{t}^{(m,n-1)}\hat{L}^{{\dagger}}\\ &-\hat{L}^{\dagger\times}\varrho_{t}^{(m+1,n)}+\hat{L}^{\times}\varrho_{t}^{(m,n+1)},\end{split}$$

(A4)

where $\hat{X}^{\times}\hat{Y}\equiv\hat{X}\hat{Y}-\hat{Y}\hat{X}$ and $\varrho_{t}^{(0,0)}$ is the reduced density operator of the probe. Similarly, the HEOM of the spin-fermion model is given by Wu and Liu (2017); Suess et al. (2015)

$$\begin{split}\frac{d}{dt}\varrho^{(m,n)}_{t}&=(-i\hat{H}_{s}^{\times}-m\beta-n\beta^{*})\varrho^{(m,n)}_{t}\\ &+\Theta(m)\alpha\hat{L}\varrho^{(m-1,n)}_{t}+\Theta(n)\alpha^{*}\varrho^{(m,n-1)}_{t}\hat{L}^{{\dagger}}\\ &+\Big{[}(-1)^{n}\varrho^{(m+1,n)}_{t}\hat{L}^{\dagger}-\hat{L}^{\dagger}\varrho^{(m+1,n)}_{t}\Big{]}\\ &+\Big{[}(-1)^{m}\hat{L}\varrho^{(m,n+1)}_{t}-\varrho^{(m,n+1)}_{t}\hat{L}\Big{]},\end{split}$$

(A5)

where $\Theta(x)\equiv x~{}\text{mod}~{}2$ and the auxiliary operators are defined as $\varrho_{t}^{(m,n)}=\mathcal{M}\{|\psi_{t}^{(m)}(\textbf{z}^{*})\rangle\langle\psi_{t}^{(n)}(-\textbf{z}^{*})|\}$.

The initial-state conditions of the auxiliary operators are $\varrho^{(0,0)}_{t}=\varrho_{\text{s}}(0)$ and $\varrho^{(m>0,n>0)}_{t}=0$. In numerical simulations, we need to truncate the hierarchical equations. In practice, we choose a sufficiently large integer $N$ and set $\varrho^{(m,n)}_{t}$ with $m+n>N$ to be zero. Then, $\varrho^{(m,n)}_{t}$ with $m+n\leq N$ form a closed set of ordinary differential equations which can be numerically solved by using the well-developed Runge-Kutta algorithm.

In the special case $\chi=1$, $\hat{H}$ reduces to the standard spin-boson model, whose reduced dynamics can be captured by the famous Zwanzing-Nakajima master equation Nakajima (1958); Zwanzig (1960)

$$\frac{\partial}{\partial t}\varrho_{\mathrm{s}}(t)=-i\mathcal{\hat{L}}_{\mathrm{s}}\varrho_{\mathrm{s}}(t)-\int_{0}^{t}d\tau\hat{\Sigma}(t-\tau)\varrho_{\mathrm{s}}(\tau),$$

(B1)

where we have used the Born approximation, $\mathcal{\hat{L}}_{\mathrm{x}}\mathcal{\hat{O}}=\hat{H}_{\mathrm{x}}^{\times}\mathcal{\hat{O}}$ with $\mathrm{x}=\mathrm{s},\mathrm{e},\mathrm{i}$ are the Liouvillian operators corresponding to the probe, the environment and the probe-environment interaction Hamiltonian, respectively. Here, $\hat{\Sigma}(t)$ is the self-energy super-operator. If the probe-environment coupling is weak, $\hat{\Sigma}(t)$ can be approximately written as DiVincenzo and Loss (2005); Burkard (2009)

$$\hat{\Sigma}(t)\simeq\mathrm{Tr}_{\mathrm{e}}\Big{[}\mathcal{\hat{L}}_{\mathrm{i}}e^{-it(\mathcal{\hat{L}}_{\mathrm{s}}+\mathcal{\hat{L}}_{\mathrm{e}})}\mathcal{\hat{L}}_{\mathrm{i}}\varrho_{\mathrm{e}}\Big{]}.$$

(B2)

In the Bloch representation, the Zwanzig-Nakajima master equation can be expressed as $\partial_{t}\vec{r}(t)=\mathfrak{S}(t)\diamond\vec{r}(t)$, where $\diamond$ denotes the convolution and

$$\mathfrak{S}(t)=\left[\begin{array}[]{ccc}-\mathcal{A}(t)&0&0\\ 0&-\mathcal{B}(t)&-\Delta\delta(t)\\ 0&\Delta\delta(t)&0\\ \end{array}\right],$$

(B3)

with $\mathcal{A}(t)=4\cos(\Delta t)C(t)$, $\mathcal{B}(t)=4C(t)$. By means of the Laplace transform, one can find

$$\begin{split}\tilde{r}_{i}(\zeta)\equiv&\int_{0}^{\infty}dtr_{i}(t)e^{-\zeta t}=\sum_{j}\mathfrak{F}_{ij}(\zeta)r_{j}(0),\end{split}$$

(B4)

where $i,j=x,y,z$. The non-vanishing terms of $\mathfrak{F}_{ij}(\lambda)$ are given by

$$\mathfrak{F}_{xx}(\zeta)=\Big{[}\zeta+\mathcal{\tilde{A}}(\zeta)\Big{]}^{-1},$$

(B5)

$$\mathfrak{F}_{yy}(\zeta)=\bigg{[}\zeta+\mathcal{\tilde{B}}(\zeta)+\frac{\Delta^{2}}{\zeta}\bigg{]}^{-1},$$

(B6)

$$\mathfrak{F}_{zz}(\zeta)=\Big{[}1+\zeta^{-1}\mathcal{\tilde{B}}(\zeta)\Big{]}\mathfrak{F}_{yy}(\zeta),$$

(B7)

$$\mathfrak{F}_{yz}(\zeta)=-\mathfrak{F}_{zy}(\zeta)=-\Delta\zeta^{-1}\mathfrak{F}_{yy}(\zeta).$$

(B8)

Then, for an arbitrary given initial state, the reduced dynamics of Bloch vector $\vec{r}(t)$ can be completely determined by Eq. (B4) with the help of inverse Laplace transform.

In the special case $\chi=0$, the RWA completely removes the counter-rotating-wave terms in $\hat{H}$, resulting in the total excitation number operator $\hat{\mathcal{N}}=\hat{\sigma}_{+}\hat{\sigma}_{-}+\sum_{k}\hat{c}_{k}^{\dagger}\hat{c}_{k}$ being a constant of motion. This character leads to the spin-boson model and the spin-fermion model share the same reduced dynamical behavior in the RWA case Wu and Liu (2017). As displayed in Refs. Wu and Liu (2017); Breuer and Petruccione (Oxford University Press, Oxford, 2002), in the basis of $\{|+\rangle,|-\rangle\}$, the exactly analytical expression of $\varrho_{\mathrm{s}}(t)$ can be expressed as

$$\varrho_{\mathrm{s}}(t)=\left[\begin{array}[]{cc}\varrho_{++}(0)\mathcal{G}_{t}^{2}&\varrho_{+-}(0)\mathcal{G}_{t}e^{-i\Delta t}\\ \varrho_{-+}(0)\mathcal{G}_{t}e^{i\Delta t}&1-\varrho_{++}(0)\mathcal{G}_{t}^{2}\\ \end{array}\right],$$

(C1)

where

$$\mathcal{G}_{t}=\exp\Big{(}-\frac{1}{2}\lambda t\Big{)}\bigg{[}\cosh\Big{(}\frac{1}{2}\Omega t\Big{)}+\frac{\lambda}{\Omega}\sinh\Big{(}\frac{1}{2}\Omega t\Big{)}\bigg{]},$$

(C2)

is the decoherence factor and $\Omega\equiv\sqrt{\gamma^{2}-2\gamma\lambda}$.

With Eq. (C1) at hand, the analytical expressions of CFI and QFI can be easily computed. Assuming the initial state parameter $\phi$ is chosen as $\phi=\pi/4$, we find

$$\mathcal{F}_{\text{C}}^{\text{RWA}}=\frac{t^{2}\mathcal{G}_{t}^{2}\sin^{2}(\Delta t)}{1-\mathcal{G}_{t}^{2}\cos^{2}(\Delta t)},$$

(C3)

and $\mathcal{F}_{\text{Q}}^{\text{RWA}}=t^{2}\mathcal{G}_{t}^{2}$. One can easily check that $\mathcal{F}_{\text{Q}}^{\text{RWA}}\geq\mathcal{F}_{\text{C}}^{\text{RWA}}$, due to the fact that $0\leq\mathcal{G}_{t}\leq 1$. Moreover, in the ideal noiseless case, we have $\mathcal{G}_{t}=1$, which leads to $\mathcal{F}_{\text{C}}^{\text{ideal}}=\mathcal{F}_{\text{Q}}^{\text{ideal}}=t^{2}$. This result means the selected measurement scheme is the optimal one in the noise-free situation.

In Fig. 5 (a), we display the the dynamics of the population difference $\langle\hat{\sigma}_{z}(t)\rangle\equiv\text{Tr}_{\text{s}}[\hat{\sigma}_{z}\varrho_{\text{s}}(t)]$ from the purely numerical HEOM method. In the especial case $\chi=0$, one can see our numerical results are completely consistent with that of the exactly analytical RWA result $\langle\hat{\sigma}_{z}(t)\rangle_{\text{RWA}}=\mathcal{G}_{t}\cos(\Delta t)$. For the case $\chi=1$, as long as the probe-environment coupling is not too strong, our HEOM results are in good agreement with the predictions from the Zwanzig-Nakajima master equation approach (see Fig. 5(b)-(d)). However, the HEOM results are believed to be more reliable, because the Zwanzig-Nakajima master equation approach neglects the higher-order probe-environment coupling terms, which is invalid in non-Markovian and strong-coupling regimes. These results verify the feasibility and the validity of the HEOM method.