Approaching Heisenberg-scalable thermometry with built-in robustness against noise

It is a major goal in quantum thermometry to reach a 1/N scaling of thermometric precision known as Heisenberg scaling but is still in its infancy to date. The main obstacle is that the resources typically required are highly entangled states, which are very difficult to produce and extremely vulnerable to noises. Here, we propose an entanglement-free scheme of thermometry to approach Heisenberg scaling for a wide range of N, which has built-in robustness irrespective of the type of noise in question. Our scheme is amenable to a variety of experimental setups. Moreover, it can be used as a basic building block for promoting previous proposals of thermometry to reach Heisenberg scaling, and its applications are not limited to thermometry but can be straightforwardly extended to other metrological tasks.


INTRODUCTION
Accurately measuring temperature is of universal importance, underpinning many fascinating applications in material science 1,2 , medicine and biology 3,4 , and quantum thermodynamics 5 . The advent of quantum technologies has opened up exciting possibilities of exploiting quantum effects to yield the quantum enhancement of thermometric precision that cannot otherwise be obtained using classical methods 6,7 . The emerging field, known as quantum thermometry nowadays, could have a large impact on quantum platforms demanding precise temperature control, such as cold atoms, trapped ions, and superconducting circuits. This has motivated a vibrant activity on quantum thermometry over the past decade [8][9][10][11][12][13][14][15][16][17][18] .
A major goal in quantum thermometry is to reach a 1/N scaling of precision known as Heisenberg scaling (HS), which represents an important quantum advantage of central interest in quantum metrology 19-21 . Here, N stands for the amount of physical resources which usually refers to the number of probes employed. Conventionally, HS is achieved by exploiting highly entangled states 8,9,15,16 . However, along this line, the HS permitted in theory is typically elusive in reality, due to the vulnerability of highly entangled states to noises as well as the difficulty in producing these states. This point has been theoretically shown in refs. [22][23][24][25] and is also reflected in the fact that less experimental progress has been made in implementing Heisenberg-scalable thermometry to date. Indeed, although the issue of how to attain HS has received much attention since the early days of quantum thermometry 8 , the first experiment demonstrating HS was carried out only recently 26 . This experiment explored NOON states to reach HS and observed that the decoherence effect becomes increasingly severe as the order of NOON states increases. As such, HS was only demonstrated for small N, i.e., N N max ¼ 9, for which the scaling advantage is far from allowing one to beat the best possible classical methods. Besides, several no-go theorems [27][28][29][30] show that HS is forbidden to reach in the asymptotic limit N → ∞ in the presence of most types of noises. For these reasons, it is crucial to find a noise-robust scheme capable of reaching HS in the regime of large N 31 which is of interest from a practical standpoint and allowed by these theorems.
Note that a number of noise-robust schemes have been proposed for other metrological tasks [32][33][34][35][36][37][38][39][40][41][42] . However, there has not been a noise-robust scheme in quantum thermometry so far. In this work, we fill the gap. We find that a temperaturedependent phase can be accumulated coherently through the continuous interplay between the strong Markovian thermalization of N probes and the relatively weak coupling of the N probes to a same external ancilla (see Fig. 1). This mechanism enables us to propose a scheme of thermometry to approach HS for a wide range of N, with robustness irrespective of the type of noise in question but without complicated error correction techniques.
A salient feature of our scheme is that the whole estimation procedure is free of entanglement, unlike in previous works 21 , where highly entangled states are either introduced in the statepreparation stage or generated via some interactions in the interrogation stage. Another feature of our scheme is that the probing time is not increasingly long with N. These two features distinguish our scheme from previous schemes like the parallel scheme and sequential scheme 43 on a fundamental level. As detailed below, our scheme is amenable to a variety of experimental setups. Moreover, it can be used as a basic building block for promoting previous proposals of thermometry to reach HS, and its applications are not limited to thermometry but can be straightforwardly extended to other metrological tasks.

RESULTS
The basic dynamical equation Suppose that we are given a sample of temperature T, with which a probe S keeps in contact. The thermalization process of S can be described by a Markovian master equation ∂ t ρ S ðtÞ ¼ L S ρ S ðtÞ 6,7 . Here, ρ S ðtÞ is the evolving state of S. L S is a Liouville superoperator assumed to have the following two properties: (i) it admits a temperature-dependent state ρ T as the unique steady state; and (ii) the nonzero eigenvalues λ μ of L S have negative real parts, where μ is an index labeling the eigenvalues. Property (i) implies that L S ρ T ¼ 0, and property (ii) means that there is a dissipative gap λ :¼ min λμ≠0 Reðλ μ Þ in the Liouvillian spectrum [44][45][46] . We do not impose any restriction on the explicit forms of L S and ρ T , so that the scheme to be presented is applicable to a variety of physical models. In particular, ρ T may or may not be the Gibbs state, depending on the specific physical model in question. Starting from an arbitrary initial state, ρ S ðtÞ undergoing the thermalization process automatically evolves towards ρ T at an exponentially fast rate [44][45][46] . Therefore, we may assume that S (approximately) reaches ρ T at a certain time t 0 .
We weakly couple S to an ancilla A after the time t 0 . The free Hamiltonian of A is denoted by H A . The interaction between S and A is assumed to be of the form H I = ℏgS ⊗ A with ½A; H A ¼ 0, where g is a positive number describing the coupling strength, and S and A are two Hermitian operators of S and A, respectively. In the frame rotating with H A , the dynamics of SA can be described by ∂ t ρðtÞ ¼ L S ρðtÞ À ig½S A; ρðtÞ ¼: LρðtÞ: (1) Note that we temporarily do not take noise into account, as our purpose here is to show how the continuous interplay between the thermalization L S and the interaction H I leads to the coherent accumulation of a temperature-dependent phase. Let us figure out the reduced dynamics of A. Inspired by the Nakajima-Zwanzig projection operator technique 47 , we introduce a superoperator P, defined as PX ¼ ρ T tr S X, for an operator X acting on the joint Hilbert space of SA. Evidently, PρðtÞ ¼ ρ T ρ A ðtÞ, where ρ A ðtÞ denotes the evolving state of A. The superoperator complementary to P, denoted by Q, is defined as QX ¼ X À PX.
To figure out the explicit expression for the first term on the right-hand side of Eq. (5), we deduce from property (i) that L S ρ T ¼ 0 and hence LPρðtÞ ¼ Àig½S A; ρ T ρ A ðtÞ. Then, from the defining property of P, it follows that where φ T :¼ trðSρ T Þ. On the other hand, using property (ii), the norm of the second term on the right-hand side of Eq. (5) can be bounded as Z t t0 dsPLGðt; sÞQLPρðsÞ where K is the superoperator defined as KX :¼ Ài½S A; X, and ε is a dimensionless constant determined by the damping basis of L S 44,45 . The proof of the bound is given in the "Methods" section. Note that λ characterizes the strength of the thermalization L S ; that is, the larger λ is, the stronger the thermalization is. Likewise, g characterizes the strength of the interaction H I . Under the condition that the thermalization is strong whereas the interaction is relatively weak, i.e., λ/g ≫ 1, we deduce from Eq. (7) that the second term on the right-hand side of Eq. (5), i.e., the memory effect, is negligible. Upon neglecting this term and inserting Eq. (6) into Eq. (5), we reach the basic dynamical equation: representing a unitary dynamics able to coherently accumulate φ T in the state of A superposed by eigenstates of A. Hereafter, in line with the studies on quantum phase estimation 43 , we refer to φ T as a temperature-dependent phase.

The scheme of thermometry
To estimate the unknown temperature T of a given sample, we put N probes, S 1 ; Á Á Á ; S N , in contact with the sample. Upon preparing each probe in ρ T at time t 0 , we couple all of the probes to a same ancilla A. The above configuration is schematically shown in Fig. 1a. A natural platform for implementing it could be the magnetic resonance force microscopy 48 (see Fig. 1b), where the spins immersed in the sample and the magnetic tip serve as the probes and the ancilla, respectively. The dynamical equation for the N probes and the ancilla in the rotating frame reads with L SnA ρ ¼ L Sn ρ À ig n ½S n A; ρ. Here, L Sn denotes the Liouville superoperator for S n , g n represents the coupling strength between S n and A, and S n is a Hermitian operator acting on S n . For simplicity, we assume that L Sn ¼ L S , g n = g, and S n = S, for all n. Our following discussion can be straightforwardly extended to the general case. To estimate the temperature of a given sample, we put N probes in contact with the sample and weakly couple all of them to a same external ancilla. The continuous interplay between the strong Markovian thermalization of the probes and the relatively weak coupling to the ancilla ensures information transmission from the probes to the ancilla. This leads to the coherent accumulation of a temperature-dependent phase, and measuring the phase yields an estimate of the temperature. b Illustration of one possible implementation of our scheme using magnetic resonance force microscopy. The spins immersed in the sample serve as the probes and the magnetic tip serves as the ancilla. The force produced by the spins on the magnetic tip affects the mechanical vibrations of the cantilever, which can be detected by optical methods.
Let E N ðt 1 ; t 0 Þ denote the dynamical map associated with Eq. (9), Then, the state of A at time t 1 reads To find the explicit expression for ρ A ðt 1 Þ, we resort to the fact that L SnA 's commute with each other and rewrite Eq. (10) as with E SnA ðt 1 ; t 0 Þ ¼ exp½L SnA ðt 1 À t 0 Þ. Besides, according to the basic dynamical equation (8), Substituting Eqs. (12) and (13) into Eq. (11), we have implying that φ T is intensified N times due to the use of N probes. In deriving Eq. (13), we have neglected the memory effect. Indeed, the memory effect can be made arbitrarily weak by choosing either a large λ or a small g so that λ/g is large enough. Here, N should be confined to a certain appropriate range, as detailed below. It is worth noting that the past two decades have witnessed enormous experimental progresses in engineering strongly dissipative processes (see the review article by Müller et al. 49 and references therein). Nowadays, experimentalists are able to realize arbitrary Markovian processes for a number of experimentally mature platforms such as trapped ions 50 and Rydberg atoms 51 . These developments imply that it is experimentally feasible to obtain a large λ. Besides, the value of g is often fixed once the platform adopted has been calibrated, whereas the value of λ may be still allowed to be tuned. For example, in the magnetic resonance force microscopy, the value of g is determined by the distance between the spins and the magnetic tip 48 . So, g is fixed once the distance has been fixed. However, λ is still tunable through varying the intensity and frequency of an external magnetic field 48 . For the above reasons, we consider the scenario that λ is large whereas g is fixed.
Without loss of generality, we set the probing time t 1 − t 0 to be a "unit" of time 1/g for simplicity. It follows from Eq. (14) that Then, an estimate of T can be obtained by measuring ρ A ðt 1 Þ. It remains to find the optimal state ρ A ðt 0 Þ and observable O. This can be carried out with the aid of the quantum Cramér-Rao theorem 52 , stating that the estimation error δT from measuring any observable is bounded by the inequality δT ! 1= ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi νF N ðTÞ p , where ν denotes the number of measurements, and F N ðTÞ is the quantum Fisher information (QFI) for ρ A ðt 1 Þ. Therefore, the optimal ρ A ðt 0 Þ is the state maximizing the QFI, and the optimal O is the observable saturating the inequality. It is well-known 43 Here, a M j i and a m j i denote the eigenvectors of A corresponding to the maximum eigenvalue a M and minimum eigenvalue a m , respectively.
In light of the above analysis, we may now specify our scheme as follows: (i) Put N probes S 1 ; Á Á Á ; S N in contact with the sample and prepare each of them in ρ T at time t 0 . (ii) Initialize the state of ancilla A to be a M j iþ a m j i ð Þ = ffiffi ffi 2 p and couple all of the probes to A for the time 1/g. (iii) Perform the measurement associated with Repeat steps (ii) and (iii) ν times to build a sufficient statistics for determining T. As the QFI for ρ A ðt 1 Þ reads the precision attained in our scheme approaches the HS Notably, the evolving state of the N probes and the ancilla is approximately ρ N T ρ A ðtÞ. Thus, no highly entangled state is involved in the whole estimation procedure of our scheme. Moreover, the probing time in our scheme is a single "unit" 1/g and does not increase as N increases. In passing, we point out that the precision of the form (17) is termed as the weak Heisenberg limit in ref. 53 .
It may be instructive to give a physical picture for comprehending how our scheme works. To gain some insight into the physical origin of Eq. (8), we may interpret the probe employed in our scheme as an effective transducer, which helps to extract the information about temperature from the sample and transmit it to the ancilla under the condition that the Markovian thermalization is strong enough so that the memory effect is negligible. The ancilla, on the other hand, may be interpreted as an information storage and is responsible for storing the information about temperature, with its initial state determining the storage capacity. Such a continuous interplay between the probe and the ancilla leads to the basic dynamical equation (8) permitting the coherent accumulation of the temperature-dependent phase φ T in the state of the ancilla. The accumulating rate is described by g, which is reminiscent of the transduction parameter in electrical or magnetic field sensing 20 . When N probes are coupled to the same ancilla, the accumulating rate is increased from g to Ng as indicated in Eq. (14), since there are now N effective transducers transmitting the information about temperature to the ancilla. We emphasize that the working principle of our scheme is built upon the very nature of open dynamics, which distinguishes our scheme from previous schemes involving Ramsey and Mech-Zehnder interferometers where the dynamics is unitary.
Built-in robustness of our scheme against noises Taking noises into account, we shall rewrite Eq. (9) as ∂ t ρðtÞ ¼ ½ðL S1A þ L noise S1 Þ þ Á Á Á þ ðL SNA þ L noise SN Þ þ L noise A ρðtÞ; Under the condition that the Markovian thermalization L Sn is strong, the QFI of our scheme in the presence of the noises can be evaluated as with coefficients Here, p k and ϕ k j i are the eigenvalue and associated eigenvector of a density operator δ which is determined by L noise A and ρ A ðt 0 Þ. The proof of Eq. (19) is presented in Supplementary Note 1, where the expression of δ is given. Using Eqs. (16) and (19), we arrive at the result that the ratio F noise N ðTÞ=F N ðTÞ is independent of N. This result means that the detrimental influence of noises on our scheme does not become increasingly severe as N increases.
To see the significance of the above result, we may recall the reason why the HS permitted in many previous schemes is extremely vulnerable to noises. So far, a lot of effort has been devoted to exploring the possibilities of encoding the parameter of interest into the relative phase of a quantum system for reaching HS, which has become one main stream of research on Heisenberg-scalable metrology nowadays. Generally speaking, along this line of development, previous schemes require either highly entangled states or increasingly long probing time. A wellknown example is the parallel scheme 43 , in which N probes are employed and a relative phase Nφ is obtained by preparing the initial state of the probes to be an N-body maximally entangled state. Evidently, the N-body maximally entangled state is increasingly vulnerable to noises as N increases. Hence, in the presence of noises, the parallel scheme usually suffers from an exponential drop in performance as N increases 20,21 . Another wellknown example is the sequential scheme 43 , in which a single probe is employed and a relative phase Nφ is obtained with an N times long probing time. Since the decoherence effect is accumulated with time, the longer the probing time is, the severer the detrimental influence of noises on the sequential scheme becomes. Likewise, the performance of the sequential scheme usually drops exponentially as N increases in the presence of noises 20,21 . As a matter of fact, most of the experiments relying on highly entangled states or increasingly long probing time are confined to the regime of very small N 20,21 .
Our scheme provides a different means of encoding the parameter of interest into the relative phase of a quantum system without appealing to highly entangled states and increasingly long probing time. Here, by saying "different means," we mean that the phase Nφ T stems from the information transmission from the N probes to the same ancilla, which are permitted by the very nature of open dynamics. The key advantage of exploring this means for reaching HS is that the detrimental influence of noises is no longer increasingly severe as N increases. This means that our scheme is robust against noises, which, as demonstrated below, allows for approaching HS for a wide range of N. Notably, the robustness of our scheme is intrinsic, since no error correction technique is involved in our scheme. Moreover, it is worth noting that the robustness of our scheme is irrespective of the specific type of noises in question, since Eq. (19) is derived without restricting the explicit forms of L noise Sn and L noise A . In passing, we would like to point out that it has been shown in several works 8,9 that the problem of temperature estimation can be mapped into the problem of phase estimation. Analogous to the parallel scheme, the schemes in these works are based on highly entangled states like NOON states and have been experimentally shown to be very vulnerable to noises 26 .

Example
Let us furnish an analytically solvable model to demonstrate the usefulness of our scheme. Consider the physical model of N qubits independently contacting with a Bosonic thermal reservoir of temperature T 12 . We take the qubits to be S 1 ; Á Á Á ; S N . The thermalization process of S n is governed by the Liouvillian Here, H Sn ¼ _Ωσ Sn z =2 is the free Hamiltonian of S n , with σ Sn i , i = x, y, z, denoting the Pauli matrices acting on S n . The last two terms stand for the process of energy exchange with the reservoir, with γ > 0 describing the exchange rate. Note that the dissipative gap λ is proportional to γ, i.e., λ ∝ γ. n ¼ ½exp _Ω=k B T ð ÞÀ1 À1 and D½σ Sn ± ρ ¼ σ Sn ± ρσ Sn ∓ À fσ Sn ∓ σ Sn ± ; ρg=2, with σ Sn ± ¼ ðσ Sn x ∓ iσ Sn y Þ=2. Note that the unique steady state of S n happens to be the Gibbs state ρ T ¼ expðÀH Sn =k B TÞ=Z T , where k B is the Boltzmann's constant and Z T ¼ tr expðÀH Sn =k B TÞ the partition function. A is also taken to be a qubit. The interaction between S n and A is set to be _gσ Sn z σ A z . So, We examine a major source of noises, dephasing. That is, the Liouville superoperator L noise Sn=A describing the noise acting on S n =A reads L noise where κ Sn=A denotes the dephasing rate. The initial state of the probes and the ancilla is set to be ρ N So far, we have specified all the Liouville superoperators entering the full dynamical Eq. (18) as well as the initial state assumed in our scheme. Now, let us figure out the performance of our scheme in the presence of the noises. To this end, we need to find out the output state of A at time t 1 = t 0 + 1/g, which is given by Here, 1/g is the probing time which is independent of N, and where ΔðξÞ ¼ ð2n þ 1Þ 2 ξ 2 À 8iξ À 16 and ω ± ðξÞ ¼ ½Àð2n þ 1Þξ ± ffiffiffiffiffiffiffiffiffi ΔðξÞ p =2, with ξ = γ/g and η ¼ κ A =g characterizing (in units of g) the strength of the Markovian thermalization and that of the dephasing noise acting on the ancilla A, respectively. The proof of Eq. (24) is presented in Supplementary Note 2. It is interesting to note that the noises acting on the probes S 1 ; Á Á Á ; S N do not affect the output state ρ A ðt 1 Þ. The QFI for ρ A ðt 1 Þ reads Equations (24) and (26) are analytical results without numerical approximation.
We first examine the limiting case of ξ → ∞, which corresponds to the ideal situation that the Markovian thermalization is infinitely strong and there is no memory effect. Using Eq. (26) and noting that lim ξ!1 Γ ¼ expðÀ2iφ T Þ, where φ T :¼ trðσ z ρ T Þ ¼ 1=ð2n þ 1Þ, we have that Equation (27) represents a HS of precision. Since we see that the HS in Eq. (27) stems from the relative phase Nφ T acquired in ρ A ðt 1 Þ. Notably, this relative phase is obtained without appealing to highly entangled states and increasingly long probing time. A consequence of this fact is that the detrimental influence of noises does not lead to an exponential drop of the HS. Indeed, using Eq. (27) and noting that which is consistent with our general analysis that F noise N ðTÞ=F N ðTÞ is independent of N.
We now examine the situation that the strength of the Markovian thermalization is large and within current experimental reach. To do this, we numerically compute F noise N ðTÞ in units of F th ðTÞ for different ξ. Here, F th ðTÞ ¼ ð_Ω=2k B T 2 Þ 2 sec h 2 _Ω=2k B T ð Þ is the QFI for the Gibbs state ρ T ¼ expðÀH Sn =k B TÞ=Z T 12 . We choose four values of ξ, namely, 100, 200, 300, and 400. To reach these values, we can set, for instance, g = 100 kHz and γ = 10, 20, 30, 40 MHz, which are experimentally feasible for a number of quantum platforms such as Rydberg atoms and trapped ions 49 . We examine an exemplary temperature satisfying k B T/ℏΩ = 2, which is  Fig. 2, where η is set to be 1/10 in view of the convention that the lifetime of an ancilla is usually assumed to be long in quantum metrology (see, e.g., refs. [33][34][35]37,40 ). Here, we also show the QFI given by Eq. (27) (see the black solid line), which serves as a benchmark for the QFIs associated with the four values of ξ.
The four dashed curves in Fig. 2 correspond to the QFIs associated with the four values of ξ (see the figure legend). As can be seen from this figure, these dashed curves approach the black solid line when N is no larger than a certain N max . This means that the HS described by Eq. (27) is attained in the regime of N N max . It is easy to see that N max depends on ξ; more precisely, the larger ξ is, the greater N max can be. We numerically find that N max may be, respectively, taken to be 109, 217, 326, 435 for the four values of ξ. Note that the HS achieved in current experiments is usually confined to the regime of N N max ¼ 10 31 due to the detrimental effect of noises. The above numerical results demonstrate that our scheme allows for approaching HS in the presence of noises for a wide range of N. Moreover, it is worth noting that strongly dissipative processes have been experimentally realized in a number of quantum platforms (see the review article by Müller et al. 49 ). This means that the HS permitted in the above wide ranges of N may be achieved using a variety of experimental setups (to be listed in the "Discussion" section), which makes our scheme appealing from a practical perspective.  Fig. 3, F noise N ðTÞ decreases as η increases, which is expected from a physical point of view. It is interesting to note that N max does not change much in the course of varying η, indicating that the range in which the HS is attained is irrespective of the value of η. Moreover, in Supplementary Note 3, we examine the influence of the initial state of the N probes and the ancilla on F noise N ðTÞ. We consider the initial state of the form ρ ⊗N ⊗ σ, where ρ and σ are arbitrarily given 2 × 2 density matrices. Both analytical and numerical results show that F noise N ðTÞ does not depend much on ρ but heavily relies on σ. More precisely, F noise N ðTÞ grows quadratically as a function of the l 1 norm of coherence of σ 54 .
We remark that the HS given by Eq. (27) is lost when N is too large. The underlying reason is that there is some memory effect for each E SnA ðt 1 ; t 0 Þ appearing in decomposition (12) of E N ðt 1 ; t 0 Þ and these memory effects may add linearly so that the total memory effect for E N ðt 1 ; t 0 Þ becomes non-negligible for a too large N. Nevertheless, as long as the Markovian thermalization is strong enough, these memory effects can be made very weak and any large value of N max can be obtained in our scheme. By the way, considering that a stronger thermalization allows for producing a larger number of Gibbs states per unit time, a direct approach of taking advantage of strong thermalization for thermometry might be to repeatedly produce and measure Gibbs states with individual probes. Yet, this approach is ineffective in practice. Indeed, real-word quantum measurements suffer from correlated background noises if the measurement time is shorter than the noise correlation time 55,56 . As these noises cannot be averaged out by repetitive measurements 55,56 , the direct approach may not even be able to beat the standard quantum limit in the presence of these noises.

DISCUSSION
The key to the implementation of our scheme is to realize the weak coupling of N probes to a same ancilla. Apart from the magnetic resonance force microscopy, this kind of coupling is achievable in a variety of experimental setups, including superconducting qubits coupled to a microwave resonator 57 , Rydberg atoms coupled to a microwave cavity 58 , trapped ions coupled to a common motional mode 59 or an optical cavity mode 60 , and nitrogen-vacancy (NV) centers in diamond coupled to a microwave mode in a superconducting transmission line cavity 61 .
Thanks to the freedom in choosing explicit forms of L S and ρ T , various kinds of probes can be used in our scheme. Particularly, the probes in previous proposals of thermometry may be used in our scheme, and our scheme may be exploited to promote these proposals to reach HS. To this end, one only needs to find an ancilla that is able to couple with the N probes given in a previous proposal. The ancilla can be a NV center if hybridization proposals 62 are under consideration, or a measuring apparatus when quantum nondemolition measurements are used 63,64 , or even an environment to which we have access 20,65 .
The applications of our scheme are not limited to thermometry but can be straightforwardly extended to other metrological tasks. Indeed, our scheme is applicable to any metrological task in which the parameter of interest can be encoded in the unique steady state of a Markovian system. The metrological tasks of this kind are far more than quantum thermometry and are frequently encountered in quantum sensing, especially in noisy quantum metrology 65 . Two examples are the detection of rotating fields 66 and the sensing of low-frequency signals 67 .
In summary, we have proposed a scheme of thermometry to approach HS for a wide range of N, based on the finding that an Nfold increase of the temperature-dependent phase can be obtained from the continuous interplay between the strong thermalization of N probes and the weak coupling to the same external ancilla. As opposed to conventional ones, our scheme gets rid of a number of experimentally demanding requirements,  e.g., the preparation of highly entangled states, and is implementable in a variety of experimental setups. More importantly, it offers the key advantage of robustness irrespective of the type of noise in question, without resorting to complicated error correction techniques. Therefore, our scheme provides a feasible and robust pathway to the HS, with the entanglement-free feature which is noteworthy in view of previous proposals of Heisenbergscalable thermometry. Two directions for future work are to exploit our scheme to promote previous proposals of thermometry to reach HS and to apply our scheme to other metrological tasks. Besides, throughout this paper, we have assumed the Markovian dissipative process described by the Lindblad equation, which is an approximative description relying on a number of simplifications. A more detailed analysis of our scheme starting from a fully microscopic description, while going beyond the scope of the present work, would be an interesting topic for future studies.

METHODS
To obtain Eq. (7), it is convenient to express L as L ¼ L S þ gK with KX :¼ Ài½S A; X. Using the equalities L S P ¼ PL S ¼ 0 44,45 , we have PL ¼ gPK; LP ¼ gKP: Note that E 1 E 2 k k E 1 k k E 2 k k for two superoperators E 1 and E 2 . Here and throughout this paper, the Hilbert-Schmidt norm is adopted. That is, for an operator X, the norm reads X k k :¼ ffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi trðX y XÞ q , and for a superoperator E, it is the induced norm defined as E k k :¼ sup X k k 1 EðXÞ k k. It follows from Eq. (30) that PL k k g P k k K k k and LP k k g K k k P k k. Using these two inequalities and noting that ρðsÞ k k¼ ffiffiffiffiffiffiffiffiffiffiffiffiffi ffi trρ 2 ðsÞ p 1, we have R t t0 dsPLGðt; sÞQLPρðsÞ R t t0 ds PLGðt; sÞQLPρðsÞ k k PL k k R t t0 ds Gðt; sÞQ k k LP k k αg 2 ; where α :¼ P k k 2 K k k 2 R t t0 ds Gðt; sÞQ k k . It can be shown that Gðt; sÞQ k k ε exp½ εg K k k À λ ð Þ ð t À sÞ: The proof of Eq. (32) is very technical and presented in Supplementary Note 4, where the expression for ε is given. Under the condition λ=g > ε K k k; Gðt; sÞQ k kdecreases exponentially as t − s increases. Then, α P k k 2 K k k 2 ε=ðλ À εg K k kÞ ½ : Inserting Eq. (34) into Eq. (31), we can immediately obtain Eq. (7).

DATA AVAILABILITY
Numerical data from the plots presented are available from D.-J.Z. upon request.

CODE AVAILABILITY
The codes used for this study are available from D.-J.Z. upon request.