Non-Markovian temperature sensing

We investigate the sensing performance of a single-qubit quantum thermometer within a non-Markovian dynamical framework. By employing an exactly numerical hierarchical equations of the motion method, we go beyond traditional paradigms of the Born-Markov theory, the pure dephasing mechanism, and the weak-coupling approximation, which were commonly used in many previous studies of quantum thermometry. We find (i) the non-Markovian characteristics may boost the estimation efficiency, (ii) the sensitivity of quantum thermometry can be effectively optimized by engineering the proportions of different coupling operators in the whole sensor-reservoir interaction Hamiltonian, and (iii) a threshold, above which the strong sensor-reservoir coupling can significantly enhance the sensing precision in the long-encoding-time regime. Our results may have certain applications for high-resolution quantum thermometry.


I. INTRODUCTION
Temperature is the most fundamental cornerstone in the theory of both classical and quantum thermodynamics. How to measure the temperature of a quantum reservoir with a high sensitivity has recently attracted much attention [1][2][3][4][5][6]. Quantum thermometry pursues a highly precise measuring of temperature, which can surpass the standard bound set by classical statistics, with the help of coherence, quantum squeezing, entanglement, or other quantum resources. In a typical thermometry scheme, a two-level system [7][8][9] or a harmonic oscillator [10][11][12], working as the quantum sensor, is coupled to the quantum reservoir of interest. Due to the sensor-reservoir interaction, the information about the reservoir's temperature is encoded into the state of the sensor. Then, by measuring a certain sensor's observables, the knowledge about the temperature can be obtained. A highly sensitive quantum thermometer has wide applications in the fields of physics, biology, and material science.
Roughly speaking, there are two totally different quantum thermometers in the present existent temperature sensing schemes: the fully thermalized thermometer [4,10,11,[13][14][15][16] and the partly thermalized thermometer [6][7][8][9][17][18][19][20][21]. If the sensor-reservoir interaction is weak and the encoding time is sufficiently long, one can regard the sensor as being completely thermalized, namely, the sensor evolves to its thermal equilibrium state which is independent of the encoding time. In many previous treatments [4,13,22], such a thermal equilibrium state is approximately expressed as a canonical Gibbs state at the same temperature as the reservoir. In this case, the temperature can be readout directly from the canonical Gibbs state, and the corresponding temperature uncertainty is related to the heat capacity of the sensor [13,22,23]. On the other scenario, if the sensor-reservoir interaction is too strong, resulting in the so-called bound state effect [24][25][26], or the encoding time * weiwu@lzu.edu.cn is short, the sensor cannot be completely thermalized, which implies the state of the sensor is unstable in the time domain. In this situation, to obtain the information of temperature, one needs to monitor the time evolution of the sensor's reduced density matrix. In comparison with the fully thermalized thermometer, the partly thermalized thermometer may have certain advantages, because the quantum coherence can be utilized as the quantum resource to boost the sensing performance, which means the maximum precision is generally achieved out of the thermal equilibrium regime [9,17,19].
However, many previous studies of the partly thermalized thermometer restricted their attentions to the Born-Markov approximation [7,8,21], the weak-sensorreservoir-coupling regime [27], or the pure dephasing encoding mechanism [6,17,19]. Such treatments can provide an analytical result as well as a intuitionistic picture, but inevitably ignore certain important physical phenomena. For example, as reported in Refs. [28][29][30][31], the authors demonstrated that the non-Markovian effect can severely influence the performance of quantum sensing; and the authors of Refs. [9,29,32] reported that the pure dephasing mechanism is not the optimal sensor-reservoir encoding scenario. Thus, to obtain a global view and more physical insights into the quantum thermometry, going beyond the above traditional paradigms is highly desirable.
To address these concerns, we employe the hierarchical equations of motion (HEOM) approach [33][34][35][36][37][38][39][40], which is a non-perturbatively numerical method, to investigate the effect of the reservoir's characteristic, as well as the form of the sensor-reservoir coupling operator on the performance of a single-qubit quantum thermometer in both partly and fully thermalized situations. This paper is organized as follows. In Sec. II, we briefly outline some basic formalism about the quantum parameter estimation. In Sec. III, we employe the HEOM method to investigate the sensitivity of a single-qubit quantum thermometry. The conclusion of this paper is drawn in Sec. IV. In the two appendixes, we provide some additional materials about the HEOM method, the Born-Markov master equation approach, as well as the varia-tional polaron transformation. Throughout the paper, we set = k B = 1 for the sake of simplicity.

II. QUANTUM PARAMETER ESTIMATION
In this section, we would like to recall some basic concepts as well as some important formalism in quantum parameter estimation theory. Generally speaking, to sense a physical quantity λ labeling on a quantum system, one first needs a quantum sensor, which is initially prepared in a suitable input state in , and coupled to the sensor to the system of interest. Due to the sensorsystem interaction, the message about λ can be encoded into the output state of the sensor via a λ-dependent mapping out = M λ ( in ) = λ . Such an encoding process can be realized by a unitary rotation [41][42][43] or a nonunitary reduced dynamics [29][30][31].
As long as the output state λ is obtained, the information about λ can be extracted by measuring the sensor's observable. More specifically speaking, one can construct a measurement operatorÔ and calculate the expected valueŌ ≡ Tr( λÔ ) as well as the variancě O 2 ≡ O 2 −Ō 2 . Then, the sensing precision with respect to { λ ,Ô} can be evaluated via the standard error propagation formula [44,45] Running over all the possible measurement schemes, one can find the optimal measurement operatorÔ opt corresponding to the ultimate sensing precision. The ultimate sensing precision is constrained by the famous quantum Cramér-Rao bound [23,46] where F λ ≡ Tr(ς 2 λ ) withς determined by ∂ λ λ = 1 2 (ς λ + λς ) is the quantum Fisher information (QFI) [23]. Specifically, if the output state λ is a twodimensional density matrix described in the Bloch representation, namely, λ = 1 2 (1 2 + r r r ·σ σ σ) with r r r being the Bloch vector andσ σ σ ≡ (σ x ,σ y ,σ z ) being the vector of Pauli matrices, the QFI can be easily calculated via [23,47] F λ = |∂ λ r r r| 2 + (r r r · ∂ λ r r r) 2 1 − |r r r| 2 .
When λ is a pure state, Eq. (2) further reduces to F λ = |∂ λ r r r| 2 . Physically speaking, the error propagation formula tells us how to extract the message about λ from the output state λ with respect to the given measurement operatorÔ, while F λ fully quantifies the most statistical information about λ contained in the output state λ . There is no general way to find the optimal measurement scheme saturating the best attainable precision determined by the QFI. In this sense, how to construct the optimal measurement scheme is of importance in the research of quantum sensing.

Sensor
Encoding Measurement Reservoir Readout FIG. 1. Schematic diagram of our quantum thermometry scheme in which a two-level system is used to as the sensor estimate the temperature of a bosonic reservoir.

III. OUR SCHEME
In this section, we propose our scheme and analyze its sensing efficiency. A two-level system or a qubit, acting as the quantum sensor, is employed to estimate the temperature of a bosonic reservoir (see Fig. 1). The total Hamiltonian of the sensor plus the finite-temperature reservoir is described as follows: whereĤ s = 1 2 σ z + 1 2 ∆σ x is the Hamiltonian of the sensor with being the bias and ∆ being the tunneling parameter, operatorsb k andb † k are annihilation and creation operators of the kth bosonic mode with frequency ω k , respectively. The parameter g k labels the coupling strength between the sensor and the kth bosonic mode. Generally, it is very convenient to encode the frequency dependence of the interaction strengths in the so-called spectral density, which is defined by J(ω) ≡ k g 2 k δ(ω − ω k ). In this paper, we consider an Ohmic spectral density with a Drude-type cutoff where χ qualifies the sensor-reservoir coupling strength, and ω c denotes the cutoff frequency. In our scheme, we consider a general coupling operator, i.e.,Ŝ has the form of [48] S =σ θ ≡ sin θσ z + cos θσ x .
By varying the coupling angle θ, all the directions on the x-z plane of the Bloch sphere can be taken into account. In Refs. [49][50][51][52], sin θσ z is named as the diagonal coupling term, while cos θσ x is called the off-diagonal coupling operator. However, such a naming convention may lead to misunderstanding when the diagonal and the off-diagonal terms are included in bothĤ s andŜ at the same time. To avoid such a problem, in this paper, we prefer to use the z-type and x-type coupling terms to indicate sin θσ z and cos θσ x inŜ, respectively. Due to the rotational symmetry ofσ θ , we restrict our study to 0 ≤ θ < π. If θ = π/2, a more z-type coupling operator is included inŜ; when θ = 0, only the x-type coupling is considered [49]. As demonstrated in Refs. [49,50,53], the authors demonstrated the x-type coupling term plays a positive role in restraining coherence and can lead to a much richer ground-state phase diagram in the spin-boson systems. These results inspire us to explore the effect of z-type and x-type coupling on sensitivity of quantum thermometry.
To realize our sensing scheme, we assume the whole sensor-reservoir system is initially prepared in a product state where |ψ s (0) = cos α|e + sin αe −iϕ |g with |e, g being the eigenstates of the Pauli z-operator andĤ b = k ω kb † kb k denotes the Hamiltonian of the reservoir. The parameter β ≡ T −1 denotes the inverse temperature and is the quantity of interest to be estimated in this paper.
During the purely numerical calculations to the QFI, one needs to handle the first-order derivative to the parameter β, say ∂ β r r r in Eq. (2). In this paper, the firstorder derivative for an arbitrary β-dependent function f β is numerically done by adopting the following finite difference method: We set δ/β = 10 −6 , which provides a very good accuracy.

A. Non-Markovian temperature sensing
In this subsection, we first investigate the dynamical behavior of QFI within an exact non-Markovian framework. Depending on the characteristic of the reservoir spectral density, quantum reservoirs can be roughly classified into two different categories: the fast reservoir ω c > max{ , ∆} and the slow reservoir ω c ≤ max{ , ∆} [52,54,55]. An interesting question arises here: what is the influence of the reservoir's characteristic on the sensing performance? In this subsection, we try to address this question.
We first consider the so-called deep fast-reservoir regime, i.e., ω c max{ , ∆}, in which the reservoir is memoryless (see Fig. 6 in Appendix A and Refs. [52,[54][55][56] for more details). In this regime, the dynamical behavior of QFI from the HEOM approach is displayed by the magenta triangles in Fig. 2. One can see the QFI gradually increases from its initial value F β (0) = 0 to the maximum value. Then, the QFI begins to decrease and eventually tends to a steady value in the long-encodingtime limit. Such behavior of F β (t) is physically reasonable because the whole temperature-sensing process includes the following three steps. First, no message about  the temperature is contained in the initial state of the sensor, resulting in F β (0) = 0 at the beginning. Then, the sensor-reservoir interaction, which generates the temperature's information in s (t) with the evolution of encoding time, leads to the increase of QFI. Finally, in the long-time limit, the value of F β (t) remains unchanged because, in this step, the message about the temperature is completely from the steady-state reduced density matrix s (∞), which is independent of the encoding time. The above result means there exists an optimal encoding time which can maximize the value of QFI. The appearance of maximal QFI can be understood as the physical result of the interplay between the encoding and the decoherence [29][30][31]57]. In Fig. 3, we plot the maximal QFI, max t F β , as a function of the coupling angle θ. One can see the maximal QFI is quite sensitive to the coupling angle. Using our above results, one can design the most efficient sensor-reservoir interaction Hamiltonian to obtain the maximum precision estimation by adjusting the weight of the x-type coupling operator in the whole sensor-reservoir interaction Hamiltonian.
Next, we consider the slow-reservoir regime, for example, the blue rectangles with ω c = 0.5∆ in Fig. 2, in which the reservoir-induced decoherence should be strongly non-Markovian [52,54,55]. In this situation, we find the dynamics of QFI presented by the HEOM method exhibits a collapse-and-revival phenomenon resulting in multiple local maximums, which is completely different from that of the fast reservoir case. A similar result was also reported in Refs. [28,29] and can be regarded as evidence of the information's backflow from the reservoir back to the sensor. Moreover, we also display the maximum QFI with respect to the optimal encoding time versus the coupling angle in Fig. 3. Multiple localmost-efficient coupling angles are revealed (see the the blue rectangles in Fig. 3). Additionally, we find the maximal QFI in the slow-reservoir regime can be much larger than that of the fast-reservoir case for the entire range of 0 ≤ θ < π. This result suggests the non-Markovian effect, generated by the characteristics of slow reservoir, may be used to improve the estimation precision.

B. Breakdown of the canonical statistics
In this subsection, we shall discuss the feature of steady-state QFI in the long-encoding-time regime. If the sensor-reservoir coupling is weak, the long-time steady state of a fully thermalized thermometer can be described by the canonical Gibbs state, namely, s (∞) G s = exp(−βĤ s )/Tr[exp(−βĤ s )]. In such an approximate treatment, by diagonalizingĤ s , one can easily find the ratio of gg (∞)/ ee (∞) is a monotonously and exponentially increasing function where Ω = √ 2 + ∆ 2 shall be regarded as the sensor's Rabi frequency. And the corresponding QFI with respect to the canonical Gibbs state is given by From the above expression, one can easily find F G β (Ω) is not a monotonic function. With a fixed temperature, F G β (Ω) is a increasing function for Ω ∈ [0, Ω * ) with Ω * being determined by ∂ Ω F G β (Ω) = 0; while, for Ω ∈ [Ω * , ∞), F G β (Ω) monotonically decreases. The above analysis means the steady-state behaviors of the population ratio and the QFI strongly depend on the sensor's Rabi frequency Ω.
However, as demonstrated in many previous articles [54,[58][59][60], in the deep slow-reservoir regime, which generally induces a strongly non-Markovian effect, the reservoir-induced decoherence can give rise to a drastic frequency renormalization, namely, Ω →Ω. The renormalized frequencyΩ is smaller than the original frequency Ω and their deviation becomes larger as the increase of the sensor-reservoir coupling [54,[58][59][60][61]. This result can be qualitatively understood by performing a polaron transformation. Making use of the polaron transformation, as displayed in Appendix B, the sensor's Rabi frequency is renormalized as where˜ = sin θ + ∆ cos θ,∆ = ∆ sin θ − cos θ, and η is the renormalized factor, whose explicit expression is given by Eq. (B8), induced by the polaron transformation. In Table I, we display the frequency renormalization within the framework of polaron transformation theorỹ Ω P /Ω versus the sensor-reservoir coupling strength χ. It is clear to see a larger χ leads to a smallerΩ P /Ω. Such a frequency renormalization phenomenon is responsible for the so-called localized-delocalized transition in the wellknown spin-boson model [58,59]. Thus, the most simple and straightforward correction to the previous assumption of canonical statistics is replacing the original Rabi frequency Ω by a smaller renormalized frequencyΩ. Under such treatment, one can see the frequency shift Ω →Ω shall lead to a smaller value of gg (∞)/ ee (∞) [see Figs. 4(a) and 5(a)] compared with that of the canonical Gibbs state case, which implies the breakdown of canonical statistics [62][63][64]. Moreover, if Ω > Ω * , the frequency renormalization results in a larger QFI in comparison with that of the canonical statistics [see Fig. 4(b)]. On the contrary, if Ω < Ω * , the frequency renormalization in turn diminishes the sensing precision [see Fig. 5 The above analytical analysis is completely based on the polaron transformation, which only provides a qualitative interpretation as well as a pictorial illustration. To obtain a quantitative result, we display the steadystate population ratio and the QFI as functions of the coupling strength from the HEOM and the Born-Markov master equation methods (see Appendix A for details) in Figs. 4(c,) and 4(d) and Figs. 5(c,) and 5(d). As a result of neglecting the higher-order terms of the sensorreservior interaction, the steady-state solutions from the Born-Markov master equation approach coincide with the canonical statistics, i.e., the steady-state population ratio and the QFI are independent of the coupling  strength. In sharp contrast to this, one can see the numerical results obtained by the HEOM method are in good agreement with our previous discussions about the frequency renormalization within the framework of the polaron transformation. In Table I, we also display the renormalized Rabi frequency predicted by the HEOM method, namelyΩ H ≡ β −1 ln[ gg (∞)/ ee (∞)], as the function of the coupling strength χ. It is clear to see the evolutionary trend ofΩ H /Ω in qualitative agreement with that ofΩ P /Ω. These results demonstrate, in the strong-coupling regime, the noncanonical distribution occurs, which plays a complicated role in the quantum thermometry. The result presented by the red stars in Figs. 4(c,) and 4(d) and Figs. 5(c,)and 5(d) cannot be predicted by using the Gorini-Kossakowski-Sudarshan-Lindblad master equation formula or other common perturbative methods, where the effect of frequency shift is generally washed out [65,66].
Therefore, we uncover a threshold Ω = Ω * , above such critical Rabi frequency, i.e., Ω > Ω * , the sensing precision can be effectively enhanced by strengthening the sensor-reservoir coupling. In this sense, the sensitivity of quantum thermometry can be improved by engineering the sensor's bare frequency. Here, we concentrate the frequency renormalization mechanism on the slow-reservoir regime, the reason is twofold. First, in the deep fastreservoir regime, the decoherence is Markovian which means the canonical statistics assumption is valid with no need for corrections. Second, as ω c / max{ , ∆} → ∞, the renormalized factor η approaches to 1 accordingly, resulting in the disappearance of frequency shift phenomenon in the deep-reservoir regime. Our result is consistent with that of Ref. [10], in which the authors found the sensing precision of a Brownian thermometer can be improved by increasing the coupling at low temperature. Their result can be easily understood by using our threshold mechanism: the critical Rabi frequency Ω * approaches to zero in the low-temperature limit, which means the condition Ω > Ω * becomes very easy to be satisfied. Thus, as long as the Brownian sensor has a finite bare frequency, i.e., Ω > 0, the sensor-reservoir can boost the thermometric precision without doubt. Of course, it is worth mentioning that our scenario is totally different from that found in Ref. [10]. Their quantum thermometry is based on the Caldeira-Leggett Hamiltonian, namely, the sensor in Ref. [10] is a continuous-variable system, rather than the discrete-variable sensor (qubit) in this paper.

IV. CONCLUSION
In summary, going beyond the usual assumptions of the Born-Markov theory, the pure dephasing mechanism and the weak-coupling approximation, we propose a single-qubit thermometer scheme and analyze its sensing performance with the help of the HEOM approach. We find, when the encoding time is short and the sensor is partly thermalized, the non-Markovian effect induced by the slow reservoir's characteristics may enhance the efficiency of the quantum thermometry. By optimizing the encoding time, we investigate the relation between the maximal QFI and the sensor-reservoir coupling angle. It is revealed that the performance of our quantum thermometry can be further boosted by modulating the weight of the x-type coupling operator in the whole sensor-reservoir Hamiltonian. Moreover, in the slowreservoir regime, by strengthening the sensor-reservoir coupling, we find the noncanonical feature appears in the sensor's steady state and the corresponding QFI obtained by the HEOM method can be larger than the value pred-icated by the canonical Gibbs state as well as the Born-Markov master equation approach, under suitable conditions. A possible threshold mechanism, which is based on the frequency renormalization appearing in the strongly non-Markovian regime, is proposed to explain the above result in the long-encoding-time limit. Finally, due to the generality of the qubit-based quantum thermometer, we expect our results to be of interest for certain potential applications in the researches of quantum metrology and quantum sensing.

V. ACKNOWLEDGMENTS
The authors thank Dr. Si-Yuan Bai, Professor Jun-Hong An and Professor Hong-Gang Luo for many fruitful discussions. W. Wu wishes to especially thank Dr. Chong Chen and an anonymous referee for their suggestions on the discussion of frequency renormalization. The work was supported by the National Natural Science Foundation (Grant No. 11704025).

VI. APPENDIX A: HEOM
The HEOM method is a purely numerical technique, which exactly maps the Schrödinger equation or quantum Liouville equation to a set of ordinary differential equations by introducing auxiliary density matrices. These ordinary differential equations can be treated numerically by the Runge-Kutta algorithm. Without invoking the Born-Markov, weak-coupling and rotating-wave approximations, the HEOM can provide a rigorous numerical dynamics for an open quantum system [33][34][35][36][37][38][39][40].
To realize the above HEOM method, it is generally required that the reservoir correlation function C(t), which is defined by can be expressed as a sum of exponential functions [33][34][35][36][37][38][39][40]67]. Here,B ≡ k g k (b † k +b k ). Fortunately, for the Ohmic spectral density considered in this paper, the reservoir correlation function satisfies such requirement. Substituting Eq. (4) into Eq. (A1), one can find where υ = ω c δ 0 + 2 π(1 − δ 0 )β −1 denotes the th Matsubara frequency and are the expansion coefficients. Here, we set a truncation parameter ε to ensure the above series remain finite.
In this paper, we constantly increasing the number of the differential equations as well as the value of ε until the final result converges. As a benchmark, we also employe the most common Born-Markov master equation to compute the reduced dynamics of the sensor. As displayed in Ref. [66], the Born-Markov master equation reads whereΥ Here,Ŝ(t) ≡ e itĤsŜ e −itĤs and C R(I) (t) denotes the real (imaginary) part of the correlation function. By neglecting the imaginary Lamb-shift terms [66], the master equation given by Eq. (A7) can be further simplified, which shall be very convenient for numerical simulations in practice.
To verify the feasibility of the HEOM method, we here make a comparison between the purely numerical result obtained by the numerical HEOM method and that of the Born-Markov master equation approach. In Fig. 6, we display the evolution of population difference Good agreement is found between results from the two different approaches in the deep fast-reservoir regime, see Fig. 6(a). On the contrary, a relatively large deviation is found in the deep slow-reservoir regime, see Figs. 66(d) and 6(e). Moreover, we observe that such a deviation becomes more and more evident as ω c /∆ approaches to zero with fixed = 0.5∆. This result is physically reasonable, because the reservoir is memoryless in the deep fast-reservoir regime. While, with the decrease of ω c /∆, the non-Markovianity becomes non-negligible [56], which means the results from the Born-Marovian master equation method are unreliable. Our result is consistent with that of Ref. [52] in which only the z-type coupling term is taken into consideration.

VII. APPENDIX B: VARIATIONAL POLARON TRANSFORMATION
To obtain the explicit expression ofΩ P , we first need a rotation matrixR ≡ exp(− i 2 φσ y ) with φ = arctan(cot θ) to diagonalizeŜ. By doing so, the original Hamiltonian H can be transformed intoĤ =R †ĤR , wherê withg k = sin(θ + φ)g k . For the special scope 0 ≤ θ < π considered in this paper,g k coincidentally reduces to g k . Employing such rotation,Ĥ now has the standard form of well-known spin-boson model.
Next, we apply the variational polaron transformation toĤ asĤ = exp(P)Ĥ exp(−P). Here, the polaron generatorP is given by [54] where ξ k = g k S(ω k ) with S(ω k ) = 1 + η 2∆2 ω kΩP coth 1 2 βω k tanh 1 2 βΩ P −1 , (B3) being self-consistently determined by minimizing the Bogoliubov upper bound for the free energy [54]. The expression ofĤ can be written asĤ =Ĥ s +Ĥ b +Ĥ sb , whereĤ s denotes the transformed Hamiltonian of the sensorĤ The last term in the above expression is an energy shift induced by the polaron transformation, which has no influence on the sensing accuracy in the assumption of canonical statistics. From the expression ofĤ s , one can immediately find the Rabi frequency isΩ 2 P =˜ 2 + η 2∆2 , which recovers Eq. (10) in the main text. The transformed sensor-reservoir interaction HamiltonianĤ sb can be written asĤ sb =σ xBx +σ yBy +σ zBz , wherê (B8)