Here, we implement a modified version of these algorithms using an artificial atom or qubit in the form of a superconducting transmon circuit.18 We show that the transmon can be operated as a dc flux magnetometer with Heisenberg-limited sensitivity. The sensitivity is boosted by a magnetic moment that is about five orders of magnitude larger than that of natural atoms or ions. The idea of the experiment is to combine the extreme magnetic-field sensitivity of superconducting quantum interference devices (SQUIDs) with an enhanced performance brought about by exploiting quantum coherence. The ‘quantum’ in the name of this device refers to the macroscopic complex wave function of the superconducting electronic state. In the SQUID loop geometry, the relative phase of the superconducting wavefunctions across the Josephson junctions acquires a dependence on magnetic flux Φ via the Aharonov–Bohm effect. However, despite its quantum origin, in standard SQUID measurements this phase is a classical variable. In contrast, for the SQUID loop of a transmon qubit, the phase turns into a fully dynamical quantum observable and the flux Φ dependence is encoded in the energy-level separation ℏω 01 (Φ) between the ground state and the first excited state. Therefore, it is possible to exploit the phase difference ϕ = [ω d − ω 01 (Φ)]τ acquired during a time τ by the qubit when it is prepared into a coherent superposition of the ground and excited energy states and driven by an external microwave field at a frequency ω d . Differently from their “natural” counterparts, where the characteristics of the quantum sensor are sample independent and defined by the atomic structure, for artificial-atom systems, such as the transmon, we need to adapt the algorithms by including device-specific properties in a so-called passport—a sample specific Ramsey interference pattern obtained in advance from characterization measurements, see Fig. 1b. Making use of phase estimation algorithms, we demonstrate an enhanced dc-flux sensitivity of the transmon sensor in an enlarged flux range as compared to standard (classical) measurement schemes. Recently, a standard measurement procedure using a flux qubit has been used for the measurement of an ac-magnetic field signal.19

Fig. 1 Experimental layout. The schematic shows a transmon qubit (in blue) comprised of a capacitor and a SQUID loop with two nearly identical junctions. The qubit charging and total Josephson energies are E C = 299 MHz and E JΣ = 26.2 GHz. The qubit is coupled via a gate capacitor C g to a coplanar waveguide resonator (CPW, in green) with a resonance frequency ω r around 2π × 5.12 GHz. The magnetic flux Φ through the transmon’s SQUID loop is controlled by a dc-current flowing through a flux-bias line (in red). An arbitrary waveform generator (AWG) and a microwave analog signal generator are employed to create a Ramsey sequence of two π/2 microwave pulses at a carrier frequency ω d = 2π × 7.246 GHz separated by a time delay τ. The sequence drives the transmon into a superposition of ground and excited states where the state amplitudes depend on the accumulated phase ϕ = [ω d −ω 01 (Φ)]τ. The qubit state is read out nondestructively using a probe pulse sent to the CPW resonator; the reflected signal is downconverted (not shown in the figure), digitized, and analyzed by a computer. Next, the computer updates a flux distribution function \({\cal P}({\mathrm{\Phi }})\) stored in its memory, determines the next optimal Ramsey delay time, and feeds it back into the AWG. a Qubit transition frequency ω 01 (Φ) as a function of magnetic flux Φ (parabolic curve). The bottom inset shows the CPW resonator’s spectrum. The red circles indicate the bias point of our transmon sensor: we operate far away from the ‘sweet spot’ in a regime where the transmon’s frequency ω 01 (Φ) is an approximately linear function of the flux Φ within the entire flux range ΔΦ. For the fluxes around the point considered here, the frequency ω r of the readout CPW resonator remains approximately constant. b A pre-measured sample-specific Ramsey interference fringes pattern defines the “passport” function of our sensor. This can be regarded as a non-normalized probability function P p (τ, Φ) to observe the qubit in the excited state after a Ramsey sequence with a delay τ for a specific value of the magnetic flux Φ. The largest flux value used to obtain the Ramsey interference fringes pattern Φ = 0.1394 Φ 0 corresponds to a frequency detuning Δω = ω d −ω 01 (Φ) = 2π × 15.8 MHz between the drive and the qubit transition frequencies. The flux range of the “passport” ΔΦ~2.5 × 10−3 Φ 0 corresponds to a range 2π × 13.8 MHz in frequency detuning Full size image

The experiment employs a superconducting circuit in a transmon configuration, consisting of a capacitively-shunted split Cooper-pair box coupled to a λ/4–wavelength coplanar waveguide (CPW) resonator realized in a 90 nm thick aluminum film deposited on the surface of a silicon substrate, see Fig. 1 and SI 1 for an image of the sensor device. The SQUID loop of the transmon has an area of S ≃ 600 μm2, which is chosen large in comparison with standard transmon qubit designs in order to provide a higher sensitivity to magnetic-field changes. The magnetic moment of this artificial atom is μ = Sℏ|dω 01 /dΦ|, directly proportional to the area S and the rate of change with flux Φ of the transition frequency ω 01 . For our device, we obtain dω 01 /dΦ = −2π × 5.3 GHz/Φ 0 at the bias point, resulting in μ = 1.10 × 105 μ B , where μ B is the Bohr magneton. By comparison, the Zeeman splitting due to the magnetic moment of NV centers is 28 GHz T−1, corresponding to a magnetic moment of 2μ B . The sample is thermally anchored to the mixing chamber plate of a dilution refrigerator and cooled down to a temperature of roughly ~20 mK. The qubit has a separate flux-bias line and a microwave gate line, the former allowing to change the qubit transition frequency, while the latter is used for the qubit’s state manipulation. The qubit state is determined by a nondemolition read-out technique (see Methods and SI 1) measuring the probe signal reflected back from the dispersively-coupled CPW resonator. To increase the magnetic field sensitivity, we bias the qubit away from the “sweet spot”, see Fig. 1a. This follows an opposite strategy as compared to the situation where the phase estimation algorithms are employed for quantum computing and simulations; in the latter cases, the qubit sensitivity to flux noise is maximally suppressed by tuning the device to the “sweet spot” characterized by a vanishing first derivative of the energy with respect to flux. Operating away from the “sweet spot” leads to a reduction of the T 2 time. The decoherence rate \(T_2^{ - 1} = \left( {2T_1} \right)^{ - 1} + T_\phi ^{ - 1}\) is the sum of the relaxation (2T 1 )−1 and dephasing \(T_\phi ^{ - 1}\) rates.20 The dephasing rate appreciably increases at our bias point, which reduces T 2 and thus the number of available steps that can be implemented in the Kitaev and Fourier algorithms.

In the experiment, we apply a Ramsey sequence of two consecutive π/2 pulses separated by a time delay τ, which corresponds to an effective spin-1/2 precession around the z-axis of the Bloch sphere. The precession angle ϕ = Δω(Φ)τ is defined by the frequency mismatch Δω(Φ) = ω d −ω 01 (Φ) between the transition frequency ω 01 (Φ) of the transmon qubit and the fixed drive frequency ω d of the π/2 pulses. The Ramsey sequence drives the transmon from its ground state into a coherent superposition of ground and excited states with relative amplitudes determined by the phase ϕ. The theoretical probability to find the transmon in the first excited state is given by

$$P\left[ {\tau ,{\mathrm{\Delta }}\omega ({\mathrm{\Phi }})} \right] = \frac{1}{2} + \frac{1}{2}{\mathrm{exp}}\left( { - \tau {\mathrm{/}}2T_1} \right)\gamma (\tau ){\mathrm{cos}}\left[ {{\mathrm{\Delta }}\omega ({\mathrm{\Phi }})\tau } \right]$$ (1)

and depends both on the delay time τ and on the magnetic flux Φ through the frequency mismatch Δω(Φ). The decay function γ(τ) accounts for qubit dephasing, typically due to charge or flux noise. By design, the transmon artificial atom is rather insensitive to background charge fluctuations. On the other hand, intrinsic 1/f magnetic-flux noise couples to the SQUID loop and is known to be a relevant source for dephasing in flux qubits;21,22,23 in addition, other decoherence mechanisms can be present, see below for details. The dephasing process can be described through an external classical noise source, see Methods. The particular shape of γ(τ) depends on the noise spectral density at low frequencies. “White” noise with a constant power density results in an exponential decay function γ(τ) = exp(−Γ wn τ), while 1/f-noise produces a Gaussian decay γ(τ) = exp[−(Γ 1/f τ)2]. We fit our experimental curves P(τ, Δω(Φ)) by Eq. (1) using both an exponential and a Gaussian decay, see SI 2. For our sample with a relaxation time T 1 of about 260 ns, we cannot distinguish between these two fits, neither in the ‘sweet spot’ nor in the bias point. Fitting the Ramsey oscillation at different fluxes one finds \({\mathrm{\Gamma }}_{{\mathrm{wn}}}^{ - 1} \approx 1250\) ns and \({\mathrm{\Gamma }}_{1/f}^{ - 1} \approx 780\) ns at the ‘sweet spot’. At the bias point, these pure dephasing times reduce to 520 and 420 ns, respectively. The decay rates Γ wn and Γ 1/f in the bias point then can be translated into equivalent white and 1/f flux noises and we find the spectral densities S wn = (5.9 × 10−8 Φ 0 )2/Hz and S 1/f (f) = (1.9 × 10−5 Φ 0 )2/f[Hz], respectively (see Methods).

The function γ(τ) determines the optimal delay time τ where the sensitivity of the probability P(τ, Δω) to the changes in Δω and hence to a flux is the highest. In the standard (classical) measurement approach, a minimal delay \(\tau = \tau _0 \ll T_2\) sets the frequency range Δω(Φ) ∈ [0, π/τ 0 ] where the phase ϕ and hence P(τ, Δω) can be unambiguously resolved. This defines the range ΔΦ = π(τ 0 dω 01 /dΦ)−1 where the magnetic flux can be resolved with a precision scaling given by the standard quantum limit (see Methods),

$$\left[ {\delta {\mathrm{\Phi }}} \right]_{{\mathrm{class}}} = \left| {\frac{{d\omega _{01}({\mathrm{\Phi }})}}{{d{\mathrm{\Phi }}}}} \right|^{ - 1}\frac{1}{{\tau _0\sqrt {t{\mathrm{/}}T_{{\mathrm{rep}}}} }} = \frac{{A_{{\mathrm{class}}}}}{{\sqrt t }},$$ (2)

where t is the total measurement or sensing time of the experiment and T rep is the time duration of a single Ramsey measurement. A better flux sensitivity can be attained at larger delays τ, where the probability P[τ, Δω(Φ)] is more sensitive to changes in Δω. We obtain the best sensitivity at τ = τ* defined by the condition (2T 1 )−1 − [ln γ(τ)]′ = τ−1 (see Methods),

$$\left[ {\delta {\mathrm{\Phi }}} \right]_{{\mathrm{quant}}} = \left| {\frac{{d\omega _{01}({\mathrm{\Phi }})}}{{d{\mathrm{\Phi }}}}} \right|^{ - 1}\frac{e}{{\tau ^ \ast \sqrt {t{\mathrm{/}}T_{{\mathrm{rep}}}} }} = \frac{{A_{{\mathrm{quant}}}}}{{\sqrt t }}.$$ (3)

The amplitudes A class and A quant in Eqs. (2) and (3) quantify the magnetic flux sensitivities. Measuring at the optimal delay τ = τ* improves the flux resolution by a factor A class /A quant = τ*/(eτ 0 ), which depends on the qubit’s coherence time, the latter serving as the quantum resource in our algorithms. Another important factor which enhances the flux sensitivity is the slope dω 01 /dΦ of the transmon’s spectrum. At our working point ω 01 = 2π × 7.246 GHz, we have dω 01 /dΦ = −2π × 5.3 GHz/Φ 0 . The minimal delay is given by τ 0 ≈31.6 ns, see SI 1. The repetition time T rep = 6.546 μs involves the maximal time duration of the Ramsey sequence, the duration of the probe pulse (2 μs) and the transmon’s relaxation time back into its ground state (4 μs, which is 15 times longer than the T 1 time). Combining these numbers and setting τ*~2T 1 , we estimate the theoretical value of flux sensitivity for our transmon sensor as A quant ≃ 4 × 10−7 Φ 0 Hz−1/2, see Eq. (3), providing an improvement by a factor A class /A quant ~ 6 over the classical sensitivity. Note, that the best sensitivity is attained at T rep = τ* (i.e., for a very fast control and readout) that gives for our sample \(\left[ {\delta \Phi } \right]_{{\mathrm{quant}}}\) ≃ \(1.1 \times 10^{ - 7}\) \(\Phi _0{\kern 1pt} {\mathrm{Hz}}^{ - 1/2}{\mathrm{/}}\sqrt t\).

Measuring at large time delays τ~T 2 leaves an uncertainty in Δω(Φ) due to the multiple 2π-winding of the accumulated phase, thereby squeezing the flux range ΔΦ ~ 2.5 × 10−3 Φ 0 by the small factor τ 0 /T 2 . The Kitaev and Fourier phase estimation algorithms, avoid this phase uncertainty by measuring the probability P(τ, Δω) at different delays τ k = 2kτ 0 for \(K\sim {\mathrm{log}}_2\left( {T_2{\mathrm{/}}\tau _0} \right)\) consecutive steps k = 0, …, K − 1. As a result, such a metrological procedure is able to resolve the magnetic flux with the quantum limited resolution [δΦ] quant , see Eq. (3), within the original flux range ΔΦ set by the duration ~ τ 0 of the control rf-pulses. The operation of the Kitaev and Fourier metrological procedures can be viewed as a successive determination of the binary digits of the index \(n = \left[ {b_{K - 1} \ldots b_0} \right]\) ≡ \(\mathop {\sum}

olimits_{k = 0}^{K - 1} b_k{\kern 1pt} 2^k\) in the so-called quantum abacus.24 The Kitaev algorithm starts from a minimal delay τ = τ 0 and determines the most significant bit b K−1 in its first step, further proceeding with the less significant bits b K−2 , …, b 0 . The Fourier algorithm works backwards:25 it starts from the maximal delay τ ~ T 2 and first determines the least significant bit b 0 , then gradually learns more and more significant bits b 1 , b 2 , …, b K−1 .

Modified Kitaev and Fourier metrological algorithms

In the present work, we use modified versions of the phase estimation protocols, which take into account the nonidealities present in actual experiments. For brevity, we will still refer to these protocols as the Kitaev and Fourier phase estimation algorithms. We demonstrate the superiority of these algorithms over the standard technique and show that we can beat the standard quantum limit. Instead of relying on the ideal theoretical probability function P[τ, Δω(Φ)] of Eq. (1) these modified Kitaev and Fourier protocols exploit the empirical probability P p (τ, Φ), the so-called passport, which we measure by a set of Ramsey sequences at various magnetic fluxes Φ, representing the result on a discrete equidistant grid in the form P p (τ j ,Φ i ), see Fig. 1b. Here, Φ i = (i − 1)[δΦ] step + Φ 1 with the index i chosen from the flux-index set I 0 = [1, 161] [δΦ] step ≃ 1.59 × 10−5 Φ 0 , Φ 1 ≃ 0.137 Φ 0 , and discrete time delays τ j = (j − 1) × 2 ns, j = 1, … 241 quantifying the time separation between the two π/2 rf-pulses of the Ramsey sequence. In order to increase the signal-to-noise ratio, we average over 65000 Ramsey experiments at each discrete point (τ j , Φ i ). The resulting pattern is only approximately described by Eq. (1) due to the fact that the resonator frequency changes slightly with the applied flux, thus modifying our calibration (see Methods and SI 2). In principle, one can change the working point to an even more sensitive part of the spectrum at the price of a further distortion of this pattern.

Using the qubit passport P p (τ j , Φ i ), one can pose the following metrological question: given an unknown flux Φ within some pre-chosen range, how can one estimate its value using a minimal number of Ramsey measurements? We design two metrological algorithms where the time delay τ of the Ramsey sequence serves as an adaptive parameter whose value is dynamically adjusted. In the course of operation, both our algorithms return a discrete probability distribution \({\cal P}\)(Φ i ), i ∈ I 0 , which reflects our current knowledge about the flux Φ to be measured. This probability distribution is improved in subsequent steps and shrinks to a narrow interval around the actual flux-value when running the algorithm.

Bayesian learning

The elementary building block for both our metrological algorithms is a Bayesian learning subroutine which updates the discrete flux distribution \({\cal P}\)(Φ i ) after each Ramsey measurement of the qubit state. This subroutine takes the time delay τ j between π/2 pulses as an input parameter and performs a sequence of N = 32 Ramsey measurements. Our readout scheme returns a measured variable h N which, at \(N \gg 1\), is equal to the empirical passport probability P p (τ j , Φ i ). At small values N, the readout variable h N is a normally distributed random variable with a mean value given by P p (τ j , Φ i ),

$$p\left( {h_N|\tau _j,{\mathrm{\Phi }}_i} \right) = \frac{1}{{\sqrt {2\pi } \sigma _N}}{\mathrm{exp}}\left[ { - \frac{{\left( {h_N - P_p\left( {\tau _j,{\mathrm{\Phi }}_i} \right)} \right)^2}}{{2\sigma _N^2}}} \right],$$ (4)

where the variance \(\sigma _N^2 = \sigma _1^2{\mathrm{/}}N\) can be directly measured, \(\sigma _1^2 \approx 3.5\) (see SI 1 for further explanations on the readout variable h N ). Next, the algorithm makes use of the measurement outcomes h N and updates the flux probability distribution with the help of Bayes’ rule, \({\cal P}\left( {{\mathrm{\Phi }}_i} \right)\) → \(p\left( {h_N|\tau _j,{\mathrm{\Phi }}_i} \right){\cal P}\left( {{\mathrm{\Phi }}_i} \right)\)/\(\mathop {\sum}

olimits_i {\kern 1pt} p\left( {h_N|\tau _j,{\mathrm{\Phi }}_i} \right){\cal P}\left( {\Phi _i} \right)\).

Kitaev algorithm

The Kitaev-type metrological algorithm has been introduced earlier in ref.26. The algorithm involves K steps k = 0, …, K − 1 with optimized Ramsey times τ k , tolerances \(\epsilon _k\), and flux index sets I k ; below, \({\cal N}\)(I) denotes the size of a discrete set I. It is initialized with a uniform discrete distribution \({\cal P}_{\mathrm{0}}\)(Φ i ) which reflects our prior ignorance of the flux to be measured. In the first step k = 0, the algorithm repeats the Bayesian learning subroutine at a zero time delay τ(0) = 0 between π/2 pulses until the probability distribution shrinks to a twice narrower interval I 1 ⊂ I 0 , i.e., \({\cal N}\)(I 1 ) = \({\cal N}\)(I 0 )/2, satisfying \(\mathop {\sum}

olimits_{i \in I_1} {\kern 1pt} {\cal P}_0\left( {{\mathrm{\Phi }}_i} \right) \ge 1 - \epsilon _0\). The flux values Φ i , i ∉ I 1 are discarded. After completing the first step, the algorithm searches for the optimal delay τ j for the next step. The next optimal Ramsey measurement requires a larger delay τ(1) > 0 such that the passport P p (τ(1), Φ i ), i ∈ I 1 , has the largest range: \(\tau ^{(1)} = {\mathrm{argmax}}_{\tau _j}{\kern 1pt} {\mathrm{\Delta }}P\left( {\tau _j} \right)\) where ΔP(τ j ) = \({\mathrm{max}}_{i \in I_1}P_p\left( {\tau _j,{\mathrm{\Phi }}_i} \right)\) − \({\mathrm{min}}_{i \in I_1}P_p\left( {\tau _j,{\mathrm{\Phi }}_i} \right)\). The algorithm thus sweeps over the passport data P p (τ j , Φ i ) to find the optimal delay τ(1) with maximal range ΔP(τ(1)). Subsequently, a new distribution \({\cal P}_1\left( {{\mathrm{\Phi }}_{i \in I_1}} \right)\) = \({\cal N}^{ - 1}\left( {I_1} \right)\) and \({\cal P}_1\left( {{\mathrm{\Phi }}_{i

otin I_1}} \right) = 0\) is initialized and the algorithm proceeds to the next step by running the Bayesian learning with the new optimal delay τ 1 . After K steps, the algorithm localizes Φ within a 2K times narrower interval I K , \({\cal N}\)(I K ) = \({\cal N}\)(I 0 )/2K, with an error probability \(\epsilon = 1 - \mathop {\prod}

olimits_{k = 0}^{K - 1} \left( {1 - \epsilon _k} \right)\).

Quantum Fourier algorithm

This algorithm starts from the Ramsey measurement with an optimal time delay τ(s)~T 2 . The starting delay τ(s) is a free input parameter of the algorithm. Similarly to the Kitaev algorithm, the quantum Fourier algorithm runs the Bayesian learning subroutine until the flux probability distribution \({\cal P}_{\mathrm{0}}\)(Φ i ), i ∈ I 0 , squeezes to a twice narrower subset S 1 ⊂ I 0 such that \(\mathop {\sum}

olimits_{i \in S_1} {\kern 1pt} {\cal P}_0\left( {{\mathrm{\Phi }}_i} \right) \ge 1 - \epsilon _0\). However, in contrast to the Kitaev algorithm, the passport function P p (τ(s), Φ i ) is an ambiguous function of Φ i at the large delay τ(s). As a result, S 1 is not a single interval but rather a set of n ~ τ(s)/τ 0 disjoint narrow intervals S 1 = I 1 ∪ … ∪ I n of almost equal lengths, see Fig. 2. Hence, after completing the first step, the flux value is distributed among n equiprobable alternatives I i . The Fourier algorithm discriminates between these n alternatives in the next steps. First, it searches for the next optimal delay τ j , where it is possible to rule out half of the remaining alternatives in the most efficient way. At each delay τ j the algorithm splits the remaining intervals I i , i = 1, …, n into two approximately equal-in-size groups \(A = I_{i_1} \cup \ldots \cup I_{i_{[n/2]}}\) and \(B = I_{i_{[n/2] + 1}} \cup \ldots \cup I_{i_n}\) which are ordered by the passport function, P p (τ j , Φ ∈ A) > P p (τ j , Φ ∈ B). Then it finds the probability distance ΔP(τ j ) = \({\mathrm{min}}_{i \in A}P_p\left( {\tau _j,{\mathrm{\Phi }}_i} \right)\) − \({\mathrm{max}}_{i \in B}P_p\left( {\tau _j,\Phi _i} \right) > 0\) separating the two sets A and B. Repeating this procedure at all available delays τ j , the algorithm finds the optimal delay τ(1) with maximal ΔP(τ j ) over the discrete set of delays τ j . In the next step, the algorithm discriminates between A and B by repeating the Bayesian learning subroutine approximately [ΔP(τ(1))]−2 times and sets S 2 = A or B. Continuing in this way, the algorithm returns a single interval I out where the actual value of the flux Φ(Φ i ), i ∈ I out , is located. Figure 2 shows how the flux distribution function \({\cal P}\)(Φ i ) develops in time during the execution of the Kitaev and Fourier algorithms.

Fig. 2 Evolution of the flux probability distribution \({\cal P}_k\)(Φ i −Φ), during the run of the first k = 1, …, 4 steps of the Fourier (red panels) and Kitaev (blue curves) estimation algorithms. The magnetic flux is measured with respect to a reference flux (as explained in the Methods). The actual flux value is shown by a thick black line in the flux-step plane. The Kitaev algorithm starts at a zero delay τ = 0 and a first step returns a broad probability distribution with a single peak centered near the actual flux value. During the run of the Kitaev algorithm this peak narrows down. The Fourier algorithm starts from the Ramsey measurement at large delay τ(s) = 360 ns, with the first step returning a probability distribution with six out of twelve flux intervals assuming a non-vanishing value. Hence, this first step selects half of the n = τ(s)/τ 0 ~12 different flux intervals ΔΦ m given by Δω m = Δω 0 + 2πi/τ s , m = 0, …, n−1, where Δω m ≡ ω 01 (ΔΦ m ) − ω d is the frequency interval corresponding to the flux interval ΔΦ m , and determines the parity of the yet unknown index m ∈ [0, n − 1] associated with the true flux interval. In the second step, the Fourier algorithm proceeds to a shorter delay and rules out another half of the remaining six intervals. In the next two steps the algorithm discriminates between the remaining three alternatives and ends up with the correct flux interval. The green line at the fourth step displays the probability distribution learned by the standard (classical) procedure during the same number of Ramsey measurements as was required by the quantum procedures. The distributions obtained at the step number 4 for the Kitaev and Fourier estimation algorithms and in the standard (classical) measurement are shown in the inset Full size image

Results

The superiority of our quantum metrological algorithms is clearly demonstrated by the scaling behavior of the magnetic flux resolution with the total sensing time of the flux measurement, see Fig. 3. We run each algorithm n = 25 times at every flux value Φ = Φ i , i ∈ I 0 , within the entire flux range, and find the corresponding arrays of estimated values \({\hat{\mathrm \Phi }}_{ji}\), j = 1, …, n. The estimate \({\hat{\mathrm \Phi }} = {\hat{\mathrm \Phi }}({\cal P}({\mathrm{\Phi }}))\) is defined as the most likely value derived from the observed probability distribution \({\cal P}\)(Φ). For a probability distribution \({\cal P}_i\)(Φ) measured at a known flux value Φ = Φ i the corresponding estimate \({\hat{\mathrm \Phi }}_{ji}\) is a random quantity due to statistical nature of the measurement procedure. We define an aggregated resolution δΦ as an ensemble standard deviation of the random variables \({\hat{\mathrm \Phi }}_{ji} - {\mathrm{\Phi }}_i\),

$$\delta {\mathrm{\Phi }}^2 = \frac{1}{{{\cal N}\left( {I_0} \right)}}\mathop {\sum}\limits_{i \in I_0} \frac{1}{{n - 1}}\mathop {\sum}\limits_{j = 1}^n \left[ {{\hat{\mathrm \Phi }}_{ji} - {\mathrm{\Phi }}_i} \right]^2.$$ (5)

Fig. 3 Observed scaling behavior of the flux resolution versus total sensing time for the three different metrological procedures, Kitaev (colored circles), Fourier (black diamonds) and standard (red crosses). The Kitaev algortihm has been run with constant tolerances \(\epsilon _k\) = \(\epsilon\) for each step k = 1, …, 5 and for five different values of \(\epsilon\) as indicated by different colors. The Fourier algorithm has been performed with the step-dependent tolerances \(\epsilon _k\) = 0.182, 0.076, 0.039, 0.02, 0.01 for k = 1, …, 5. We show the result of the Fourier algorithm only for the final two steps, k = 4 (filled diamonds) and k = 5 (empty diamonds), running the algorithm with four different starting delays, τ(s) = 300, 320, 340, and 360 ns (all collapsed to the same data points). The phase estimation algorithms lag behind in precision at short times when compared to the standard procedure, but rapidly gain precision at longer times. The black solid line represents the scaling law for a numerical simulation of the standard procedure with a regular passport function given by Eq. (1). The crossover to the red solid line is due to the irregularity of the passport function Full size image

In case of the Fourier algorithm, such a definition is meaningful only at the final step of the algorithm where \({\cal P}\)(Φ) becomes a single-peaked function. The sensing time t is defined through the total number of calls of the Bayesian learning subroutine m, t = NT rep m. The scaling behavior of the measured flux resolution δΦ(t) with sensing time t is shown in Fig. 3 for both our algorithms and is compared with the scaling δΦ std (t) of the standard (classical) procedure, where all Ramsey measurements are done at a zero delay τ = 0.

Both quantum algorithms clearly outperform the standard procedure, with the Kitaev algorithm appearing slightly more efficient than the Fourier one. We explain this by the fact that the Fourier algorithm strongly relies on the periodicity of the Ramsey interference pattern given by Eq. (1), whereas our readout scheme produces a slightly distorted pattern. On the other hand, the Kitaev algorithm turns out to be more stable to the irregularities in the measured passport function P p (τ, Φ). The magnetic flux sensitivities A quant range within 5.6–7.1 × 10−6 Φ 0 Hz−1/2 for the Kitaev algorithm and within 6.5–8.5 × 10−6 Φ 0 Hz−1/2 for the Fourier procedure. These sensitivites are an order of magnitude worse than the theoretical bound 4.0 × 10−7 Φ 0 Hz−1/2 set by Eq. (3). The discrepancy has two main reasons. First, our readout scheme is not a single-shot measurement, which leads to a factor 32 increase of the T rep time. Second, we spend part of the time resource for the intermediate steps with τ k ≤ T 2 during the run of the phase estimation procedure. Finally, for our transmon, the SQUID area S ≃ 20 × 30 μm2 results in a magnetic filed sensitivity in the range 19.3–29.3 pT Hz−1/2.

Decoherence processes define the most important factor limiting the sensitivity of our device. E.g., the intrinsic 1/f flux noise22,23 caused by magnetic impurities constitutes a relevant source of decoherence. At short times 0 < τ < 2T 1 , τ the duration of a single Ramsey sequence, the presence of 1/f noise can be accounted for by a finite coherence time T ϕ of the qubit. Assuming that dephasing originates exclusively from intrinsic flux noise results in an upper limit S 1/f (f) = (1.9 × 10−5 Φ 0 )2/f[Hz] for the noise spectral function. At much larger time scales, as defined by the entire duration t of the Kitaev or Fourier procedure, 1/f noise causes low-frequency flux fluctuations \(\left\langle {\delta {\mathrm{\Phi }}^2} \right\rangle \sim {\int}_{1/t}^{1/\tau ^ \ast } {\kern 1pt} S_{1/f}(f)df\). As follows from Fig. 3, the 5-step Kitaev procedure takes ≈0.05 s, which provides a value \(\left\langle {\delta {\mathrm{\Phi }}^2} \right\rangle\)~(6.4 × 10−5 Φ 0 )2 for the flux fluctuations, about twice larger than the actually achieved flux resolution δΦ~3 × 10−5 Φ 0 . This suggests that 1/f flux noise has a smaller weight and another, non-magnetic decoherence mechanism is present in our device. One of the potential candidates derives from electron tunneling at defects inside the dielectric layer of the qubit’s Josephson junctions. These fluctuating charges produce 1/f noise in the critical current and hence affect the transition frequency of the transmon atom.27 The 1/f flux noise may become more pronounced at a larger size L of the transmon’s SQUID loop as the flux-noise spectral density grows linearly with the loop size.23 Consequently, the flux resolution [δΦ] degrades \(\propto \sqrt L\) when increasing the loop size, while the corresponding magnetic field resolution [δB] ∝ [δΦ]/L2 still improves as L growths, see Methods.

Interestingly, the non-ideality of the qubit’s passport strongly affects the performance of the standard procedure as well. Its scaling behavior δΦ std ∝ t−α exhibits a crossover in the scaling exponent α, assuming a value α ≈ 0.39 at short sensing times, while at large times α decreases to a much smaller value ≈0.046. The scaling exponent 0.39 deviates from the shot noise exponent 1/2 due to the cases when the actual flux value is located near the boundaries of the flux interval Φ ∈ [Φ 1 , Φ 161 ], where the passport function P p (0, Φ) has an extremum and the scaling exponent for the standard procedure is reduced to 1/4, δΦ(t) ∝ t−1/4. As a result, the aggregated scaling exponent of Eq. (5) is reduced below 1/2. On the other hand, the crossover to α ≈ 0.046 is a consequence of the irregularity of the passport function set by low-frequency noise fluctuations during the passport measurement. Indeed, at large sensing times, the standard procedure needs to distinguish fluxes within a narrow interval where the passport function P p (τ = 0, Φ) has a non-regular and non-monotonic dependence on Φ. As a result, the Bayesian learning procedure fails to converge to a correct flux value. In contrast, at a larger scale of Φ, the passport function is smooth and monotonic and the standard procedure behaves properly. These arguments are indeed confirmed by a numerical simulation with a regular passport function given by Eq. (1). Importantly, both our quantum metrological algorithms are more stable than the standard procedure with respect to passport imperfections and their scaling behavior at large sensing times coincides with the scaling behavior resulting from a regular passport function. The quantum algorithms suffer, however, from the same irregularity problem at larger sensing times, not shown in Fig. 3.

Finally, we discuss how our metrological algorithms use the quantum resource of qubit coherence in order to acquire information about the measured flux. We quantify the quantum coherence resource spent in a given measurement by the total phase accumulation time \(\tau _\phi = N\mathop {\sum}

olimits_k {\kern 1pt} \tau _km_k\), where m k is the number of calls for the Bayesian learning subroutine with delay τ k . The amount of information ΔI acquired during the measurement is given by a decrease of the Shannon entropy ΔI = H(\({\cal P}_{\mathrm{0}}\)) − H(\({\cal P}\)), where \({\cal P}_{\mathrm{0}}\)(Φ i ) and \({\cal P}\)(Φ i ) are the initial and final probability distributions and H(\({\cal P}\)) = \(- \mathop {\sum}

olimits_i {\kern 1pt} {\cal P}\left( {{\mathrm{\Phi }}_i} \right){\mathrm{log}}_2\left( {{\cal P}\left( {{\mathrm{\Phi }}_i} \right)} \right)\). The scaling behavior \({\mathrm{\Delta }}I\left( {\tau _\phi } \right) \propto \tau _\phi ^\alpha\) separates the classical domain with 0 < α ≤ 0.5 from the quantum domain with 0.5 < α ≤ 1, where α = 1 corresponds to the ultimate Heisenberg limit. Indeed, in the ideal case where no relaxation and decoherence phenomena are present, the quantum algorithms double the flux resolution (squeeze the flux distribution function into a twice narrower range) for each next step of the procedure. This means that the associated Shannon entropy decreases by ln(2) and one learns one bit of information for each doubling of the Ramsey delay time. In contrast, the classical procedure with \(N \gg 1\) repetitions results in a Gaussian probability distribution of the measured quantity where the precision scales as \(\delta {\mathrm{\Phi }} \propto 1{\mathrm{/}}\sqrt N\). Hence the associated Shannon entropy scales as τ0.5 with an invested total phase accumulation time τ = Nτ 0 . We run each quantum algorithm 25 times for every flux value and average over the obtained information gains and phase times. The resulting scaling dependence is shown in Fig. 4 and demonstrates that both Kitaev and Fourier algorithms indeed belong to the quantum domain with a scaling exponent within [0.624, 0.654] (95% confidence interval). The scaling exponent is still below the Heisenberg limit, which is a consequence of the finite dephasing time T 2 : at large time delays τ ~ T 2 the visibility of the Ramsey interference pattern decreases, requiring more Ramsey measurements in order to learn the next bit of information.