Critical quantum metrology with fully-connected models: from Heisenberg to Kibble–Zurek scaling

We start by introducing and putting into context fully-connected models and quantum metrological protocols. In fully-connected models such as the QR and LMG models, we can define an effective 'system size', η, which controls the non-linearity of the system; in the case of the Rabi model, it is given by the frequency ratio of the qubit and the field, while in the LMG or Dicke model, it corresponds to the number of qubits. In the so-called thermodynamic or scaling limit, η  →  ∞, all of these models can be effectively described as $\hat{H}={\hat{H}}_{0}+{\hat{H}}_{1}$, with (ℏ = 1)

$$\begin{equation}{\hat{H}}_{0}=\omega \left[\frac{{\hat{p}}^{2}}{2}+(1-{g}^{2})\frac{{\hat{x}}^{2}}{2}\right],\end{equation} \tag{ 1 }$$

$$\begin{equation}{\hat{H}}_{1}=\omega \frac{f(g)}{\eta }{\hat{x}}^{4},\end{equation} \tag{ 2 }$$

up to first order in 1/η, where $\hat{x}$ and $\hat{p}$ are the quadratures of a bosonic field, $\hat{x}=\frac{\hat{a}+{\hat{a}}^{{\dagger}}}{\sqrt{2}}$ and $\hat{p}=\frac{i({\hat{a}}^{{\dagger}}-\hat{a})}{\sqrt{2}}$. Here g is an effective and dimensionless coupling strength, ω is the typical frequency scale of the system, while f(g) is a function of the dimensionless coupling that depends on the considered model (but is typically of order 1). In the thermodynamic limit η  →  ∞, this model undergoes a QPT at g_c = 1 (cf section 2). In this study, we will always remain in the normal phase, for g < g_c. This stands in contrast with other approaches where one actually crosses the critical point, typically to exploit symmetry-breaking effects [8, 10, 13, 17].

The working principle of the considered families of protocols is as follows. We assume that all parameters are known, except the one to be estimated, such as for example g or ω. Notice that this framework is relevant in the design of practical sensors [42]. First, the system is initialized in its ground state, far from the critical point. Then the value of a controllable parameter is changed in order to push the system in proximity of the phase transition. Finally, we perform a measurement on the final state, profiting from the high critical susceptibility to gain information about the parameter we want to estimate.

For the sake of clarity, let us compare the working principle of the protocols here considered with the standard interferometric approach to quantum metrology. In the latter, a given probe system is initialized, then the phase to be estimated is imparted on the probe, and finally a measurement is performed. On the contrary, in the protocols here considered the information about the parameter to be estimated is encoded in the probe system during the sweep or quench of the controllable parameter (which is varied to bring the system close to the critical point). In this sense, the quantum criticality is not used only as a mean to generate a quantum state that is subsequently used in a parameter estimation protocol: the critical nature of the system is exploited to efficiently encode information about the parameter to be estimated onto the probe system itself.

The precision achievable by a given estimation protocol is upper-bounded by QFI I_x = 4[⟨∂_x ψ|∂_x ψ⟩ − |⟨∂_x ψ|ψ⟩|²], where |ψ⟩ is the system state at the end of the protocol, which depends on the unknown x. The QFI is a figure of merit of theoretical relevance, which corresponds to the estimation precision when the optimal measurement is performed, a task which can be highly nontrivial. In practice, one rather considers the squared SNR Q_x = x² I_x , which gives the estimation precision obtained with a specific measurement setup.

1.1.1. Scaling regimes

We consider three different preparation protocols, as illustrated in figure 1. First, sudden quenches in which the coupling is abruptly increased from g = 0 to its final value g_f ∼ g_c = 1, followed by an evolution for a time T. Second, adiabatic ramps in which the coupling varies in time g(t) toward the QPT always fulfilling the adiabatic condition $\dot{{\Delta}}(t)\ll {{\Delta}}^{2}(t)$ [56–59], where Δ denotes the energy gap between ground and first-excited states, which vanishes at the QPT as Δ ∼ |g − g_c|^zν , where here zν = 1/2 [22, 25, 45]. We found that if g evolves according to

$$\begin{equation}g(t)={\left(1-\frac{1}{1+{(t/{\tau }_{Q})}^{2}}\right)}^{1/2},\end{equation} \tag{ 3 }$$

with τ_Q some time constant verifying τ_Q ≫ 1/ω, then the evolution remains adiabatic at all time. However, the critical point is only approached in the long-time limit, but never reached. These two families of protocols (sudden quenches and adiabatic ramps) epitomize the 'dynamic' and 'static' approaches, and constitute the two poles of our analysis. The third family of protocols is given by finite-time ramps in which

$$\begin{equation}g(t)=1-{\left(\frac{T-t}{T}\right)}^{r},\end{equation} \tag{ 4 }$$

where r > 0 is an exponent which describes the non-linearity of the ramp close to the QPT [60]. Contrary to the adiabatic ramp, this protocol allows one to reach the critical point in finite time, but it does not ensure perfect adiabaticity. By tuning r and T, we can make the evolution more or less adiabatic, and thus draw a connection between the two previous protocols.

Zoom In Zoom Out Reset image size

As sketched in figure 1, we found three different scaling regimes for the QFI depending on the total duration of the protocol T. The first one, I, concerns fast evolutions, in which the scaling of the QFI depends on the parameter to be estimated in a non-universal manner. Note that this regime is characterized by a small SNR Q_x = x² I_x ≪ 1. The second and third regimes, II and III, are universal in the sense that we obtain the same QFI scaling for almost any parameter x, and any fully-connected model. The regime II is valid for $\frac{1}{\omega }\ll T\ll 1/{\Delta}$, where here Δ represents the gap at g = g_c for finite η. In this regime, the finite-size effects are not yet relevant, and thus one can ignore the quartic term (2), so that the system behaves as in the thermodynamic limit. For a sudden quench in this regime, we found that the QFI scales as (ωT)⁶. For the adiabatic preparation, the precision is dominated by the ground-state properties of the Hamiltonian. The QFI in proximity of the QPT diverges as ${I}_{x}\sim \vert g-{g}_{\text{c}}{\vert }^{-\tilde{\gamma }}$ with $\tilde{\gamma }=2$ [61]. Plugging in the profile (3), we find that the QFI scales as ${(\omega T)}^{\tilde{\gamma }/(z\nu )}={(\omega T)}^{4}$, as previously reported in [22].

The finite-time ramp leads to a scaling which depends on the behavior of g(t) close to the critical point, i.e. on the non-linear exponent r. For r ≫ 1 and T ≲ r/ω, the finite-time ramp mimics a sudden quench, and hence the characteristic T⁶ scaling is recovered. We refer to this regime as IIa. On the contrary, for r/ω ≪ T ≪ 1/Δ, one enters another domain, which we call IIb. In this region, the QFI obeys a Kibble–Zurek scaling law, such that ${Q}_{x}\!\sim \!{(\omega T)}^{\tilde{\gamma }r/(z\nu r+1)}$, which in this case gives a (ωT)^4r/(2+r) scaling. In the limit r ≫ 1 and $T > \frac{r}{\omega }$, i.e., for a very slow and non-linear evolution, the adiabatic scaling (ωT)⁴ is retrieved. Therefore, the finite-time ramp draws a connection between the two extreme cases studied before, that is, for short enough T the ramp behaves as a sudden quench, while for long T and high r, it becomes a fully adiabatic evolution.

Finally, if the protocol time surpasses 1/Δ, we enter the regime III, where finite-size effects can no longer be neglected. In this regime, we find that a sudden quench leads to (ωT)², restoring an Heisenberg-like behavior. To the contrary, for the adiabatic and finite-time ramps, the QFI saturates to its maximum ground-state value ${I}_{x}\sim {{\Delta}}^{-\tilde{\gamma }/(z\nu )}$, and thus it displays a T⁰ scaling. In the thermodynamic limit, the energy gap at the critical point vanishes, Δ  →  0, and the third regime is pushed to T  →  ∞. It is worth mentioning that comparing the three strategies, we find that a sudden quench to g = 1 always gives the highest QFI for a given T, both with and without finite-size effects.

1.1.2. Precision bound

In addition to this, we have derived an upper bound to the QFI which allowed us to re-derive the various scaling regimes, with little to no calculations. This bound is valid if the system state is a squeezed state evolving under a quadratic, time-dependent Hamiltonian of the form

$$\begin{equation}{\hat{H}}_{x}(t)=(\hat{x},\hat{p}){h}_{x}(t){(\hat{x},\hat{p})}^{\text{T}},\end{equation} \tag{ 5 }$$

with h_x (t) being a (in general time-dependent) two-by-two matrix which encodes the parameter x to be evaluated. This simple model captures most of the metrological protocols with fully-connected systems, in the thermodynamic limit. In general, this evolution does not conserve the average number of excitations; hence it belongs to the category of active interferometry. For such an evolution, we found that the QFI is bounded by

$$\begin{equation}{I}_{x}(T)\leqslant 8{\left[{\int }_{0}^{T}\mathrm{d}t\sqrt{\chi {(t)}^{2}+\phi {(t)}^{2}}\left(2N(t)+1\right)\right]}^{2},\end{equation} \tag{ 6 }$$

with T the duration of the protocol, N(t) the average number of excitations at time t, $N(t)=\langle \psi (t)\vert {\hat{a}}^{{\dagger}}\hat{a}\vert \psi (t)\rangle$, and ϕ(t) and χ(t) are the eigenvalues of the matrix M_x (t) = ∂_x h_x (t). Note that N(t) is a very coarse-grained description of the system, while χ and ϕ can be obtained just from ${\hat{H}}_{x}$. Therefore, this expression allows to bound the QFI with minimal information about the state of the system. Notice that, contrary to previous bounds [62–66], this expression takes explicitly into account the fact that the number of excitations varies in time. We stress that this bound efficiently reproduces the scaling in T and general features of regimes II and III. Furthermore, although this bound has been derived for squeezed states and purely quadratic Hamiltonian, we showed that a similar expression can be found when the state is coherently displaced, which could be used to describe more general protocols.

1.1.3. Saturability of the QFI

1.1.4. Effect of decoherence

The effective model that describes the low-energy physics of fully-connected models consists in a non-linear oscillator, whose Hamiltonian is given by equations (1) and (2). We are particularly interested in the regime where η  →  ∞, i.e. in ${\hat{H}}_{0}$ given in equation (1), which we will call in the following the scaling or thermodynamic limit. In this limit, the Hamiltonian exhibits a second-order QPT at the critical point g_c = 1. The ground state of ${\hat{H}}_{0}$ is a squeezed vacuum state, which can be expressed in two equivalent ways, that is,

$$\begin{align}\hfill \vert {\psi }_{b}\rangle & =\mathrm{exp}\left[\frac{1}{2}\left(\vert z\vert {\text{e}}^{-\text{i}\theta }{\hat{a}}^{2}-\vert z\vert {\text{e}}^{\text{i}\theta }{a}^{{\dagger}2}\right)\right]\vert 0\rangle \hfill \\ \hfill & ={(1-4\vert b{\vert }^{2})}^{1/4}\mathrm{exp}\left[{\text{e}}^{\text{i}\theta }\vert b\vert {a}^{{\dagger}2}\right]\vert 0\rangle ,\hfill \end{align} \tag{ 7 }$$

where we have defined z = |z|e^iθ and $b=\frac{1}{2}\,\mathrm{tanh}(\vert z\vert ){\text{e}}^{\text{i}\theta }$. The direction of squeezing is encoded in θ. The squeezing norm can be expressed either through |z| or |b|. Each of these two conventions can be more or less convenient depending on the calculation to make; in the following, we will use both. The number of excitations is related to the squeezing parameter via N = sinh²(|z|). In the ground-state of (1), we have $b=-1/2+{(1+\sqrt{1-{g}^{2}})}^{-1}$ and θ = 0 for 0 ⩽ g ⩽ 1. The quadrature fluctuations change with the coupling as

$$\begin{align}\hfill \langle {\hat{x}}^{2}\rangle & \propto \frac{1}{\sqrt{1-{g}^{2}}},\hfill \\ \hfill \langle {\hat{p}}^{2}\rangle & \propto \sqrt{1-{g}^{2}}.\hfill \end{align} \tag{ 8 }$$

The spectrum is harmonic, and the energy gap is equal to

$$\begin{equation}{\Delta}=\omega \sqrt{1-{g}^{2}},\end{equation} \tag{ 9 }$$

so that Δ ∼ |g − g_c|^zν with zν = 1/2 for |g − g_c| < 1. When the system approaches the critical point, the squeezing diverges and the number of excitations becomes infinite.

For η < ∞, however, the quartic term (2) is no longer negligible near the critical point, and stabilizes the system. Although the model can then no longer be solved exactly, one can still resort to a variational approach to extract scaling arguments, supported by exact numerical simulations [25]. The quartic term will be non-negligible when $\langle (1-{g}^{2}){\hat{x}}^{2}\rangle \sim \langle \frac{f(g)}{\eta }{\hat{x}}^{4}\rangle$, where the average is taken over the ground state in the η  →  ∞ limit. Since the thermodynamic-limit ground state is Gaussian, Wick's theorem gives $\langle {\hat{x}}^{4}\rangle \sim {\langle {\hat{x}}^{2}\rangle }^{2}$ up to an irrelevant prefactor. Hence, we find that the quartic term starts to play a role for 1 − g² ∼ η^−2/3. When the quartic term becomes important, the quadrature variance can no longer diverge. Instead, it saturates $\langle {\hat{x}}^{2}\rangle \sim {\eta }^{1/3}$, and becomes weakly dependent on g. Similarly, the gap stabilizes around a finite value, Δ ∼ ωη^−1/3. Finally, above the critical point g > 1, the effective $\hat{x}$ potential assumes a double-well structure [25, 34, 35], and the ground state becomes degenerate, which marks the phase transition. This structure corresponds to the usual Landau potential for second-order phase transitions.

To summarize, the region g < 1 can be divided into two regimes: on the one hand, for g² < 1 − η^−2/3, the system remains harmonic, and the quantities describing the system scale with g. On the other hand, for 1 − η^−2/3 < g² < 1, the variance and gap saturate at certain values which scales with η. This region is denoted as the critical region. Its typical width, denoted by Γ, scales as Γ ∼ η^−1/ν , with ν = 3/2. As η increases, the critical region shrinks, and the saturation value of the gap becomes smaller.

This behavior can be retrieved by using the concept of scale-invariance and critical exponents [25, 27, 68–73]. We give here a brief sketch of the argument, and refer the reader to appendix A for more details. Even though the model (2) has no spatial structure, we can still define a renormalization group-like approach [70, 71], in which the field quadratures, rather than the position in space, are rescaled. More precisely, performing a transformation $\hat{p}\;\to \;\alpha \hat{p}$, $\hat{x}\to \frac{1}{\alpha }\hat{x}$, $\eta \to \frac{1}{{\alpha }^{6}}\eta$, $\omega \to \frac{1}{{\alpha }^{2}}\omega$, 1 − g² →  α⁴(1 − g²), the Hamiltonian (1) and (2) remains approximately invariant; the invariance becomes exact at the critical point. Thus, quantities such as the energy gap must also be invariant. Now, in the thermodynamic limit, the gap is given by ${\Delta}=\omega {(1-{g}^{2})}^{1/2}$. For general g and η, we can write ${\Delta}\sim \omega {(1-{g}^{2})}^{1/2}f(g,\eta ,\omega )$, with f some scaling function [70, 71]. Such scaling function must satisfy several constraints, since Δ must be invariant under the scaling transformation and must become independent of η in the limit η  →  ∞, and independent of g in the critical region, in the limit 1 − g² ≪ Γ. The simplest expression satisfying these constraints is ${\Delta}\sim \omega {(1-{g}^{2})}^{1/2}f\left(\frac{1-{g}^{2}}{{\Gamma}}\right)\sim \omega {(1-{g}^{2})}^{1/2}f\left(\frac{1-g}{{\eta }^{1/\nu }}\right)$, with f(x)  →  1 for x ≫ 1, and f(x)  →  x^−1/2 for x ≪ 1. In general, for any quantity A that behaves as A ∝ |g − g_c|^α in the η → ∞ limit, one can write $A\propto \vert g-{g}_{\text{c}}{\vert }^{\alpha }{h}_{A}\left(({g}_{\text{c}}-g){{\Gamma}}^{-1}\right)=\vert g-{g}_{\text{c}}{\vert }^{\alpha }{h}_{A}\left(({g}_{\text{c}}-g){\eta }^{1/\nu }\right)$ for η < ∞, with h_A (x) → 1 for x ≫ 1, and h_A (x)  →  x^−α for x ≪ 1. This means that, within the critical region 1 − g² < Γ, the quantity A will saturate at a value that scales as Γ^α = η^−α/ν . Note that this corresponds to the standard finite-size scaling in spatially extended systems. Therefore, knowing how A scales with g in the thermodynamic limit, we can infer how it scales with η in the critical region. This can be summarized by the following heuristic:

Heuristic 1. The scaling of a quantity in the critical region can be obtained by taking its scaling with g in the thermodynamic limit, and substituting 1 − g² with Γ ∼ η^−2/3. Or said differently, a quantity for η finite and 1 − g² ⩽ η^−2/3 is equal to the same quantity in the thermodynamic limit, at a coupling g^*2 = 1 − η^−2/3: the value inside the critical zone is equal to the value near the edge of the zone.

We can readily verify that this heuristic is satisfied for both Δ and $\langle {\hat{x}}^{2}\rangle$. In the thermodynamic limit, we have ${\Delta}\sim \omega \sqrt{1-{g}^{2}}$ and $\langle {\hat{x}}^{2}\rangle \sim {(1-{g}^{2})}^{-1/2}$; near the critical point, we have ${\Delta}\sim \omega {\eta }^{-1/3}\sim \omega \sqrt{{\Gamma}}$ and $\langle {\hat{x}}^{2}\rangle \sim {\eta }^{1/3}\sim {{\Gamma}}^{\,-1/2}$. Exact numerical simulations [25] confirm that this behavior holds in general.

Let us consider a single two-level system, or a spin, interacting with a bosonic mode according to the QR model Hamiltonian,

$$\begin{equation}{\hat{H}}_{\text{QR}}={\Omega}{\hat{\sigma }}_{z}+\omega {\hat{a}}^{{\dagger}}\hat{a}+2\lambda ({\hat{a}}^{{\dagger}}+\hat{a}){\hat{\sigma }}_{x},\end{equation} \tag{ 10 }$$

where ${\hat{\sigma }}_{i}$ are Pauli operators describing the qubit (we take the convention [σ_x , σ_y ] = iσ_z ), and $[a,{\hat{a}}^{{\dagger}}]=1$. This model can describe a large number of physical systems, such as a single atom interacting with a photonic mode in the context of cavity QED [74], the interaction of internal degrees of freedom of a trapped ion with the vibrational motion [75], or artificial atoms interacting with an LC resonator in superconducting circuits [76]. We are interested in evaluating one of the parameters ω, Ω, or λ, assuming the other two are known. This task can be mapped to several estimation problems using the above-mentioned platform; for example, in trapped ions, the qubit frequency Ω can be sensitive to an external magnetic field. Therefore, the estimation of Ω could be useful for space-resolved quantum magnetometry.

In the limit ω ≪ Ω, this model can be mapped to the non-linear oscillator (1) and (2), by the mean of a perturbative Schrieffer–Wolff transformation [25]. Let us define the frequency ratio $\eta =\frac{{\Omega}}{\omega }$. We apply a unitary operator $\hat{S}={\text{e}}^{\text{i}\frac{\lambda }{{\lambda }_{\text{c}}\sqrt{\eta }}{\hat{\sigma }}_{y}(\hat{a}+{\hat{a}}^{{\dagger}})-\frac{{\lambda }^{3}}{3{\lambda }_{\text{c}}^{3}\eta \sqrt{\eta }}{\hat{\sigma }}_{y}{(\hat{a}+{\hat{a}}^{{\dagger}})}^{3}}$, where we have defined ${\lambda }_{\text{c}}=\frac{\sqrt{\omega {\Omega}}}{2}$. This leads to ${\text{e}}^{\text{i}\hat{S}}{\hat{H}}_{\text{QR}}{\text{e}}^{-\text{i}\hat{S}}={\Omega}{\hat{\sigma }}_{z}+\omega \frac{{\lambda }^{2}}{2{\lambda }_{\text{c}}^{2}}{\hat{\sigma }}_{z}{(\hat{a}+{\hat{a}}^{{\dagger}})}^{2}+\omega {\hat{a}}^{{\dagger}}\hat{a}-\omega \frac{{\lambda }^{4}}{8{\lambda }_{\text{c}}^{4}\eta }{\hat{\sigma }}_{z}{(\hat{a}+{\hat{a}}^{{\dagger}})}^{4}$ plus higher-order terms of order $\frac{1}{\eta \sqrt{\eta }}$ and smaller. Such transformation is valid for 0 ⩽ λ ⩽ λ_c, where the spin and the boson are now decoupled, up to second order in perturbation theory. We can now project the spin in its low-energy subspace |↓_z ⟩, and obtain the effective Hamiltonian

$$\begin{equation}{\hat{H}}_{\text{QR,eff}}=\omega \left\{\frac{{\hat{p}}^{2}}{2}+\left[1-{\left(\frac{\lambda }{{\lambda }_{\text{c}}}\right)}^{2}\right]\frac{{\hat{x}}^{2}}{2}+\frac{{\lambda }^{4}}{4{\lambda }_{\text{c}}^{4}\eta }{\hat{x}}^{4}\right\}-\frac{{\Omega}+\omega }{2},\end{equation} \tag{ 11 }$$

which is in the same form of the Hamiltonian $\hat{H}$ given in equations (1) and (2), up to a constant term, and with $g=\frac{\lambda }{{\lambda }_{\text{c}}}=\frac{2\lambda }{\sqrt{\omega {\Omega}}}$, $\eta =\frac{{\Omega}}{\omega }$, and f(g) = g⁴/4. Therefore, the physics of the QR model can be captured by the non-linear oscillator model, where g is the physical spin–boson coupling, normalized by the frequencies, and η is the ratio between the qubit and boson frequencies. Previous numerical simulations confirm that this is a faithful description of the system close to the critical point [25].

For g ⩽ 1, the system finds itself in the so-called normal phase. To first order in 1/η in the η  →  ∞ limit, the ground state reads as |↓_z ⟩|ψ_b ⟩ where the field is in a vacuum squeezed state. As one gets closer to the critical point g_c = 1, the fluctuations and the number of bosons increase, while the spin remains unperturbed due to the large energy difference. Beyond the critical point, for g > 1, the system enters the so-called superradiant phase, with a doubly-degenerate ground state, which feature a bosonic population ∝η and a coherent state in the spin degree of freedom [25]. A direct correspondence can be established between this phenomenology and the superradiant QPT of the Dicke [34, 77–79] and LMG model (see below, and for example reference [27]), with the frequency ratio $\eta =\frac{{\Omega}}{\omega }$ playing the role of a large number of spins. Therefore, we see here that the parameter η can be interpreted as an effective system size, even in this model with no spatial extension. It is worth mentioning that the existence of such superradiant QPT in light–matter systems has been subject to a vibrant theoretical debate [80–85]. Such fundamental limitation can however be sidestep relying on effective implementations of ultrastrongly-coupled systems [86–90], as demonstrated by recent experimental observations of this QPT in distinct platforms [91–93].

The LMG model [55] describes a system of N spins coupled through an all-to-all interaction, whose Hamiltonian can be written as

$$\begin{equation}{\hat{H}}_{\text{LMG}}=h{\hat{J}}_{z}-\frac{{\Lambda}}{N}{\hat{J}}_{x}^{2}.\end{equation} \tag{ 12 }$$

Here ${\hat{J}}_{z}={\sum }_{i}^{N}{\hat{\sigma }}_{z}^{i}$ is a collective operator describing the excitations of the chain of spins. The term ${\hat{J}}_{x}^{2}$ accounts for the all-to-all interaction, ${\sum }_{i,j}{\hat{\sigma }}_{x}^{i}{\hat{\sigma }}_{x}^{j}$, while the parameter Λ controls its strength and so the nature of the ground state of the LMG. The LMG can be considered as a limiting case of the Ising model with long-range interactions, and thus has been proven very useful to test different aspects of critical quantum dynamics and the role of long-range interactions [94–104], which has been experimentally realized in trapped-ion setup [105–108] and with cold gases [109–111].

In the thermodynamic limit N → ∞, the system undergoes a QPT at Λ_c = h, as shown in [45, 46, 112, 113]. For |Λ| ⩽ Λ_c the system is in the normal paramagnetic phase, while it enters in the symmetry-broken ferromagnetic phase for Λ > Λ_c where the ground state is two-fold degenerate, and ⟨J_x ⟩ ≠ 0 acts as a good order parameter. In the thermodynamic limit, this system can be mapped to the non-linear oscillator model (1) and (2), by performing a Holstein–Primakoff transformation [112, 114]. The LMG Hamiltonian commutes with the total spin operator ${\hat{J}}^{2}={\hat{J}}_{x}^{2}+{\hat{J}}_{y}^{2}+{\hat{J}}_{z}^{2}$, and therefore, the Hilbert space can be split into sectors corresponding to the value of ${\hat{J}}^{2}$. In each spin sector, the interaction term can be decomposed as ${\hat{J}}_{x}=\frac{1}{2}({\hat{J}}_{+}+{\hat{J}}_{-})$, where the operators ${\hat{J}}_{\pm }$ describe raising and lowering within the spin ladder. In the subspace with largest angular momentum J = N/2, and for 0 ⩽ Λ ⩽ Λ_c, the Holstein–Primakoff transformation maps the spin states to a bosonic field according to ${\hat{J}}_{z}=-\frac{N}{2}+{\hat{a}}^{{\dagger}}\hat{a}$, and ${\hat{J}}_{+}={\hat{a}}^{{\dagger}}\sqrt{N-{\hat{a}}^{{\dagger}}\hat{a}}$. Intuitively, the operator ${\hat{J}}_{+}$ is mapped on a bosonic creation operator, plus some extra term which encodes the non-linearity of the spin operator. In the limit N  →  ∞, this non-linearity becomes negligible, and we have ${\hat{J}}_{+}\sim \sqrt{N}{\hat{a}}^{{\dagger}}$. Then we can expand the spin operator as ${\hat{J}}_{+}=\sqrt{N}{\hat{a}}^{{\dagger}}-\frac{{\hat{a}}^{{\dagger}2}\hat{a}}{2\sqrt{N}}$, plus higher-order terms. In this manner, one finds

$$\begin{equation}{\hat{H}}_{\text{LMG,eff}}=h\left(\frac{{\hat{p}}^{2}}{2}+\left(1-\frac{{\Lambda}}{h}\right)\frac{{\hat{x}}^{2}}{2}\right)+\frac{{\Lambda}}{4N}\left({\hat{x}}^{4}+\frac{{\hat{x}}^{2}{\hat{p}}^{2}+{\hat{p}}^{2}{\hat{x}}^{2}}{2}-2{\hat{x}}^{2}+\frac{1}{2}\right)-\frac{h}{2},\end{equation} \tag{ 13 }$$

which is very similar to the effective non-linear oscillator, equations (1) and (2), upon the identification ω = h, ${g}^{2}=\frac{{\Lambda}}{h}$, η = N and f(g) = g²/4. Note that the constant term $-\frac{h}{2}$ can be safely ignored. In the limit N → ∞, the quartic term vanish; the ground state is a squeezed state with $\langle {\hat{x}}^{2}\rangle \propto {(1-{g}^{2})}^{-1/2}$, and $\langle {\hat{p}}^{2}\rangle \propto {(1-{g}^{2})}^{1/2}$, which leads to a diverging entanglement among the N coupled spins as g → 1 [115]. The higher-order terms start to become relevant close to the point g = 1. At this stage, the fluctuations of $\hat{x}$ become dominant over $\hat{p}$. Hence, $\langle {\hat{p}}^{4}\rangle \ll \langle {\hat{p}}^{2}{\hat{x}}^{2}\rangle \ll \langle {\hat{x}}^{4}\rangle$, and $\langle {\hat{x}}^{2}\rangle \ll \langle {\hat{x}}^{4}\rangle$. Hence, finite-size effects appear primarily as a quartic potential ${\hat{x}}^{4}$ as in the non-linear oscillator. Therefore, the phenomenology of the LMG reduces to that of the non-linear oscillator (1) and (2).

We now discuss how a metrological protocol exploiting quantum critical effects can be implemented. For concreteness, we will consider a system described by the QR model (10); the discussion would be exactly the same for other fully-connected models. Let us assume, for instance, that we want to evaluate the frequency ω, assuming that the other parameters are known and controllable. We prepare the system in its ground-state at λ = 0, with the boson and qubit decoupled. Then we change the coupling constant from 0 to a certain target value λ, within a time T. As discussed in section 1.1, we consider three families or protocols, i.e., quenches, adiabatic ramps and finite-time ramps. These different profiles are sketched in the left side of figure 1. In all cases, the system at the end of the evolution can be written as

$$\begin{equation}\vert {\psi }_{\omega }(\lambda ,T)\rangle =\vert {\psi }_{f}\rangle .\end{equation} \tag{ 14 }$$

In the left side, we have indicated explicitly that the final state depends both on the final coupling value λ, the protocol duration T, and the (unknown) bosonic frequency ω. On the right side, we have used a short-hand notation to lighten the equations. Then we measure some observable on the system, such as the number of bosonic excitations, spin state, etc. Finally, the measurement results are used to reconstruct the value of ω. Depending on the choice of observable, the evaluation may be more or less precise. A standard result in quantum metrology [116] states that if the choice of observable is optimized, the maximum achievable precision is bounded by the QFI I_ω according to

$$\begin{equation}\delta \omega \geqslant \frac{1}{\sqrt{{I}_{\omega }(\lambda ,t)}},\end{equation} \tag{ 15 }$$

where the QFI reads as

$$\begin{equation}{I}_{\omega }=4\left[\langle {\partial }_{\omega }{\psi }_{f}\vert {\partial }_{\omega }{\psi }_{f}\rangle -\vert \langle {\partial }_{\omega }{\psi }_{f}\vert {\psi }_{f}\rangle {\vert }^{2}\right].\end{equation} \tag{ 16 }$$

Here δω denotes the standard deviation of the estimated ω, and |∂_ω ψ_f ⟩ is the derivative of the state |ψ_f ⟩ with respect to ω. Under certain conditions [116], it is guaranteed that there exists a choice of observable which allow to saturate this bound. For example, in the case of the evaluation of ω in the critical Rabi model, some of us have shown in [22] that quadrature or photon-number measurements on the bosonic field are optimal. In the case of pure state, the QFI can also be identified with the susceptibility [117], a quantity commonly used in the condensed matter community. The same reasoning can be applied if, instead of ω, one aims at evaluating the coupling λ, or any other parameter x. This reasoning can also be extended to mixed states, for which the QFI has a more involved expression [116]. Here we will only consider pure states.

To compute the QFI, we can now apply the mapping which we discussed in the section 2. For fully-connected models, the system can be described in terms of the non-linear oscillator model, equations (1) and (2). We consider first the thermodynamic limit, when the quartic term is negligible. In this case, whether through the quench or ramp process, the bosonic mode evolves under an effective quadratic Hamiltonian. At the end of the evolution, the system is in a squeezed vacuum state |ψ_b(x,T)⟩, where the squeezing parameter b depends both on the unknown parameter x and on the protocol duration T. For squeezed states, the expression (16) can be rewritten in terms of the derivative of the squeezing parameter (see appendix B for details of the derivation)

$$\begin{align}\hfill {I}_{x}& =2\left({\left(\frac{\partial \vert z\vert }{\partial x}\right)}^{2}+{\mathrm{cosh}}^{2}(\vert z\vert ){\mathrm{sinh}}^{2}(\vert z\vert ){\left(\frac{\partial \theta }{\partial x}\right)}^{2}\right)\hfill \\ \hfill & =\frac{8}{{(1-4\vert b{\vert }^{2})}^{2}}{\left\vert \frac{\partial b}{\partial x}\right\vert }^{2}\hfill \\ \hfill & =\frac{8}{{(1-4\vert b{\vert }^{2})}^{2}}\left({\left(\frac{\partial \vert b\vert }{\partial x}\right)}^{2}+\vert b{\vert }^{2}{\left(\frac{\partial \theta }{\partial x}\right)}^{2}\right).\hfill \end{align} \tag{ 17 }$$

Recall that $b=\vert b\vert {\text{e}}^{\text{i}\theta }=\frac{1}{2}\,\mathrm{tanh}(\vert z\vert ){\text{e}}^{\text{i}\theta }$. Therefore, if we want to evaluate a parameter x using the Rabi or LMG model, we can find analytically the expected precision by mapping the model to the effective bosonic model, express the squeezing as a function of x, then use equation (17) to compute the QFI. Rather than the QFI itself, we will mostly focus on the quantity Q_x = x² I_x , which gives the squared SNR of the estimation protocol.

For finite η, the evolution is no longer quadratic, and the state cannot be expressed analytically. However, we can still obtain the QFI relying on numerical simulations truncating the Fock state basis, which complements the analytical results previously obtained.

Although the achievable precision can be obtained exactly by computing the QFI, this computation is, in general, a difficult task. Even for the very simple non-linear oscillator model equations (1) and (2), numerical simulations have to be used unless η  →  ∞. The computation becomes even more challenging when mixed states are considered. Therefore, there has been considerable efforts to derive bounds which are insensitive to the specifics of the evolution. Several bounds can be found in the literature [118–121], although there is sometimes some confusion about their range of validity. For the sake of clarity, we provide here a short review of these bounds, and when they can or cannot be used. Let us consider a probe system evolving for a time T under an Hamiltonian ${\hat{H}}_{x}(t)$, which depends on the unknown parameter x. $\hat{H}$ may be in general time-dependent and/or depend on x is a non-trivial way. At the end of the evolution, the parameter is now encoded in the final state |ψ_x (T)⟩. By measuring this state, we can now evaluate x with a precision bounded by the QFI: ${I}_{x}=4\left[\langle {\partial }_{x}{\psi }_{x}(T)\vert {\partial }_{x}{\psi }_{x}(T)\rangle +{(\langle {\partial }_{x}{\psi }_{x}(T)\vert {\psi }_{x}(T)\rangle )}^{2}\right]$. Without computing this expression exactly, we can find useful bounds by making several assumptions about the evolution. In particular, let us consider the following set of assumptions:

(1) The probe system is composed of a fixed number of probes N. (2) The Hamiltonian ${\hat{H}}_{x}$ is bounded, and acts independently on each probe ${\hat{H}}_{x}={\sum }_{i}\;{\hat{H}}_{x}^{i}$. (3) ${\hat{H}}_{x}$ depends linearly in the unknown parameter x: ${\hat{H}}_{x}=x\hat{A}$, with $\hat{A}$ some operator independent of x. (4) ${\hat{H}}_{x}$ is time-independent.

Although this list of assumptions may seem long, it is satisfied in the vast majority of current metrological protocols, in particular in atomic interferometry and atomic clocks. When these conditions are satisfied, the achievable QFI scales at most quadratically with the time and the number of probes [119]

$$\begin{equation}{I}_{x}(T)\sim {N}^{2}{T}^{2}.\end{equation} \tag{ 18 }$$

This is the so-called Heisenberg limit, which is the backbone of most works in quantum metrology. Ubiquitous as it is, however, the Heisenberg bound only applies when the above list of conditions is satisfied. Several studies have shown how relaxing one or several of these conditions allows to achieve so-called super-Heisenberg scaling. An early example can be found in the work of Boixo et al [120]. They considered a situation in which the Hamiltonian ${\hat{H}}_{x}$ acts on several probes at once, and cannot be written as a sum of local contributions. This is the case, for instance, if the parameter to be evaluated is a interaction strength between neighboring spins on a lattice. Let us consider, for instance, that the Hamiltonian ${\hat{H}}_{x}$ involves two-body interactions, ${\hat{H}}_{x}=x{\sum }_{i}\;{\sigma }_{x}^{i}{\sigma }_{x}^{i+1}$. In this case, Boixo et al [120] have shown that the QFI can scale as

$$\begin{equation}{I}_{x}(T)\sim \frac{{N}^{4}}{{(2!)}^{2}}{T}^{2},\end{equation} \tag{ 19 }$$

which is indeed a super-Heisenberg scaling in N. In general, for a Hamiltonian involving k-body interactions, the QFI can scale as N^2k T². Another possibility is to look at time-dependent protocols. This was first studied in details by Pang and Jordan in reference [121]. Let us consider that the Hamiltonian is now time-dependent, but we still have ${\hat{H}}_{x}(t)=x\hat{A}(t)$ and $\hat{A}(t)$ is bounded. Then we can define its maximum and minimum eigenvalues, which we will call λ_M (t) and λ_m (t), respectively. In this case, Pang and Jordan showed that the QFI could be bounded by

$$\begin{equation}{I}_{x}(T)\leqslant {\left[{\int }_{0}^{T}\mathrm{d}t\ \vert {\lambda }_{M}(t)-{\lambda }_{m}(t)\vert \right]}^{2}.\end{equation} \tag{ 20 }$$

This limit allows to go beyond quadratic time-scaling. Let us consider, for instance, that the Hamiltonian ${\hat{H}}_{x}(t)$ is local, but increases linearly with time. Then we will have in general λ_M,m (t) ∝ Nt, which may result in

$$\begin{equation}{I}_{x}(T)\sim {N}^{2}{T}^{4}.\end{equation} \tag{ 21 }$$

and more generally, if ${\hat{H}}_{x}(t)$ scales like t^α , we can achieve a scaling T^2α+2. Hence, it is possible to achieve faster-than-quadratic scalings in time by using time-dependent Hamiltonians. Therefore, we see the Heisenberg scaling can be surpassed, provided that the system is non-linear in time, or the Hamiltonian is non-local, which means that the eigenvalues and observables do not scale linearly with the system size. We can summarize this in the following heuristic:

Heuristic 2. To beat the Heisenberg scaling, go non-linear.

For time-dependent protocols, controlling the system state may be difficult. A common solution then is to use quantum control techniques, such as counter-diabatic driving or shortcuts to adiabaticity [54]. At first sight, quantum control would seem to be very promising in the context of critical quantum metrology. Indeed, let $\hat{H}=x\hat{A}$ be a Hamiltonian with a critical point. Near this point, the system will generally be infinitely sensitive to a perturbation. However, preparing the system adiabatically will take infinite time because of the vanishing energy gap. One may want to apply quantum control to quickly bring the system near the critical point, and hence enjoy infinite precision in a finite time. However, the control term must be fully known, which means in particular that it must be independent of the unknown parameter x. Hence, the total operation must be of the form

$$\begin{equation}{\hat{H}}_{x}(t)=x\hat{\mathcal{A}}+\hat{\mathcal{B}}(t),\end{equation} \tag{ 22 }$$

where $\hat{\mathcal{B}}$, which contains all the control terms, is independent of x, and $\hat{\mathcal{A}}$ is time-independent. For such a Hamiltonian, it was shown very recently by Gietka et al [38] that the QFI scales at most like T². Hence, naive control shortcuts does not allow to achieve infinite precision in finite time, or even to achieve super-quadratic scalings in time. However, this does not exclude the possibility of reaching super-Heisenberg scaling; it only shows that, if we want to achieve such a result, the part of the Hamiltonian which encodes the parameter, $\hat{A}$, needs to be itself time-dependent. In the course of this article, we will show how simple adiabatic ramps or quenches can indeed be used to achieve non-trivial time scalings in a quantum critical system, without using quantum control.

Let us first consider that we want to estimate ω, the other parameters being known. We switch instantaneously the physical coupling to a target value λ, and then let the system evolve freely. We can eliminate the spin degree of freedom, and we are left with the bosonic field evolving under (1) and (2), with $g=\frac{2\lambda }{\sqrt{\omega {\Omega}}}$. In the thermodynamic limit η  →  ∞, the Hamiltonian is purely quadratic, and the state at all time will be a squeezed state of the form (7). In particular, the squeezing parameter b adopts the following form (see appendix D for the details of the derivation)

$$\begin{equation}b(t)=\frac{2-{g}^{2}}{2{g}^{2}}+\frac{i\sqrt{1-{g}^{2}}}{{g}^{2}\,\mathrm{tan}\left[\sqrt{1-{g}^{2}}\omega t-i\,\text{arctanh}\left(\frac{2\sqrt{1-{g}^{2}}}{2-{g}^{2}}\right)\right]}.\end{equation} \tag{ 26 }$$

If we quench the system away from the critical point, g² < 1, the system will evolve periodically in time, with a period given by the inverse of the gap, i.e. $\tau =\frac{\pi }{\omega \sqrt{1-{g}^{2}}}$. Moreover, from equation (26) it follows that b(nτ) = 0 with n integer, and $b((n+1/2)\tau )=\frac{{g}^{2}}{2(2-{g}^{2})}$ its maximum value. The photon number $N(t)=\frac{4\vert b(t){\vert }^{2}}{1-4\vert b(t){\vert }^{2}}$ has the same periodicity, and its minimum and maximum values are 0 and $\frac{{g}^{4}}{4(1-{g}^{2})}$, respectively. We stress that we have made no approximation here, aside from setting η  →  ∞. By contrast, if we quench the system directly at the critical point g = 1, the expression above can be analytically continued, and simplifies to

$$\begin{equation}b(t)=\frac{\omega t}{2(\omega t-2i)}.\end{equation} \tag{ 27 }$$

In this case, we find that the number of photons grows indefinitely in time according to $N(t)=\frac{{(\omega t)}^{2}}{4}$. This is only possible at g = 1 and η  →  ∞, when there is no stabilizing quartic term in the Hamiltonian, cf equation (2). The time-evolution for the photon number and the squeezing angle are displayed on figure 2, both at and away from the critical point.

Zoom In Zoom Out Reset image size

By combining (26) and (17), we can now compute the QFI and the SNR Q_ω = ω² I_ω . The formula (17) involves the derivative of b with ω. We have so far expressed b in terms of g and ω. However, if we rewrite everything in terms of the physical parameters, we find that g itself depends on ω, since we have $g=\frac{2\lambda }{\sqrt{\omega {\Omega}}}$. Therefore, the QFI will involve two components, coming from $\frac{\partial g}{\partial \omega }\frac{\partial b}{\partial g}$ and $\frac{\partial b}{{\partial }_{\omega }}$, respectively. As we discuss in appendix E, the first contribution actually dominates in most parameter regimes. On figure 3 we plot Q_ω with respect to the duration of the protocol T and the distance to the critical point, 1 − g². The SNR shows plateaux, separated by intervals $\frac{\pi }{\omega \sqrt{1-{g}^{2}}}$. The log–log plot reveals that, for $T > \frac{1}{\omega \sqrt{1-{g}^{2}}}$, the QFI has a secular increase ${\left(\frac{\omega T}{(1-{g}^{2})}\right)}^{2}$. In the top row of figure 4, we study more systematically the scaling of the SNR with T. We see that the SNR shows actually three different regimes. For short durations ωT ≲ 10, we obtain a quartic scaling Q_ω ∝ T⁴. For $10\lesssim \omega T\lesssim \frac{1}{\sqrt{1-{g}^{2}}}$, the SNR scales like T⁶. Finally, beyond the first plateau, for $\omega T\;\gtrsim \;\frac{1}{\sqrt{1-{g}^{2}}}$, the SNR settles on a quadratic scaling T². If we quench the system to values g closer to 1, the quadratic regime kicks in later. In the limit g  →  1, the first plateau extends to T  →  ∞, and the T⁶ scaling lingers for ever.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The same analysis can be conducted when the unknown parameter is the physical coupling λ. In the bottom row of figure 4, we show the scaling of Q_λ = λ² I_λ with time. For short time, the SNR scales like T², instead of T⁴ for Q_ω . For longer time, the features of Q_λ and Q_ω are exactly the same. As we discuss in details in appendix E, this is because a change in either ω or λ have essentially the same effect, which is a renormalization of the effective coupling variable g. We can actually make this argument more general. Let us consider any model which can be mapped to (1). We want to evaluate a physical parameter x. Then as long as the renormalized coupling g depends on x , the SNR Q_x for long T will have the same behavior, independently of x and of the model. In particular, we will get the same T⁶ and T² scalings if we want to evaluate h or Λ in the LMG model, since we have $g=\sqrt{\frac{{\Lambda}}{h}}$ in this case. A similar dynamical behavior was also obtained recently by Chu et al [37]. Here, we show that these scalings have a broad range of application. We will also show in section 5.2 how these results extend to the regime of finite η as well.

We can gain a better intuition of this complex interplay of scalings by using the general bound (25). Let us start with the case g = 1, when the system is quenched at the critical point, and let us consider that we want to evaluate λ. Then we have (we recall that g = λ/λ_c):

$$\begin{equation*}{\partial }_{\lambda }\hat{H}=-\omega \frac{\lambda }{{\lambda }_{\text{c}}^{2}}{\hat{x}}^{2}=-\omega \frac{g}{{\lambda }_{\text{c}}}{\hat{x}}^{2}.\end{equation*}$$

This operator is clearly of the form (23), with ${M}_{\lambda }=-\omega \frac{g}{{\lambda }_{\text{c}}}\left[\begin{matrix}\hfill 1\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill 0\hfill \end{matrix}\right]$. The matrix is already diagonal, with eigenvalues $\phi =-\omega \frac{g}{{\lambda }_{\text{c}}}$ and χ = 0, which are both time-independent. Taking the expression of N(t) and plugging these elements in (25), we find a bound for the SNR, which reads as

$$\begin{align}\hfill {Q}_{\lambda }& ={\lambda }^{2}{I}_{\lambda }\leqslant 8{\omega }^{2}{g}^{4}{\left[{\int }_{0}^{T}\mathrm{d}t\left(2N(t)+1\right)\right]}^{2}\hfill \\ \hfill & =\left(\frac{2}{9}{(\omega T)}^{6}+\frac{8}{3}{(\omega T)}^{4}+8{(\omega T)}^{2}\right).\hfill \end{align} \tag{ 28 }$$

Similarly, for the bosonic frequency ω, we find ${\partial }_{\omega }\hat{H}=\left(\frac{{\hat{p}}^{2}}{2}+(1-{g}^{2})\frac{{\hat{x}}^{2}}{2}\right)+\frac{{g}^{2}}{2}{\hat{x}}^{2}$, and thus $\phi =\chi =\frac{1}{2}$. This leads to a SNR bound

$$\begin{equation}{Q}_{\omega }\leqslant \frac{1}{2}\left(\frac{2}{9}{(\omega T)}^{6}+\frac{8}{3}{(\omega T)}^{4}+8{(\omega T)}^{2}\right).\end{equation} \tag{ 29 }$$

Hence, for short time, the general bound predicts a quadratic scaling. However, as soon as T is large enough (typically for ωT > 10), the T⁶ term dominates. Therefore, our bound correctly predicts that, for a quench at the critical point g = 1, the SNR features a T⁶ for long enough time. Figures 4(a) and (b) also show the general bound prediction. For Q_λ there is an excellent agreement for all times, i.e., the bound is saturated. For Q_ω , the bound fails to predict the correct scaling for shorter times, but a qualitative agreement is recovered for ωT > 10 (in this regime, the bound actually overestimates the actual result by a factor 2, but captures the scaling). Note that one can retrieve the T⁶ scaling with hardly any calculation by noting that for long enough times, the photon number is large N(t) ≫ 1, and the bound becomes essentially proportional to ${\left[{\int }_{0}^{T}\mathrm{d}t\,N(t)\right]}^{2}$. The dynamics (27) gives a N(t) which scales quadratically with t; therefore, the integral ${\int }_{0}^{T}N(t)\mathrm{d}t$ scales like T³, and the bound is T⁶. Therefore, by simply looking at the scaling of N with time, we can understand the behavior of the SNR for a quench at the critical point.

For a quench away from the critical point, the quadratic scaling Q_ω ∝ T² can also be retrieved from a simple argument. First, let us note that the integral ${\int }_{0}^{T}N(t)\mathrm{d}t$ is monotonic in T. Therefore, the precision predicted by the general bound can only grow with the protocol duration. Second, although its exact expression is a bit involved, the photon number predicted by (26) is periodic in time, of periodicity $\tau =\frac{\pi }{\omega \sqrt{1-{g}^{2}}}$. Let us assume that the average photon number oscillates between zero and its maximum value N_M. Let us define $\alpha =\frac{1}{{N}_{\text{M}}\tau }{\int }_{0}^{\tau }N(t)\mathrm{d}t$. Except in very specific cases (for instance, if N increases in very short bursts), α will be typically of order 1. Then we have ${\int }_{0}^{t}N({t}^{\prime })\mathrm{d}{t}^{\prime }={N}_{\text{M}}(\alpha t+F(t))$, where $F(t)={\int }_{0}^{t}\frac{N(t)}{{N}_{\text{M}}}-\alpha \,\mathrm{d}t$ is a periodic function, with |F| ⩽ 1, and F(nτ) = 0 for integer n. For long t, we will have F(t) ≪ αt. Then, if the QFI saturates the general bound, it follows

$$\begin{equation}{I}_{x}(T)\sim 32({\chi }^{2}+{\phi }^{2}){\alpha }^{2}{N}_{\text{M}}^{2}{T}^{2}\sim {N}_{\text{M}}^{2}{T}^{2},\end{equation} \tag{ 30 }$$

for long T, and

$$\begin{align}\hfill {I}_{x}(T=n\tau )& =8({\chi }^{2}+{\phi }^{2}){\left[2n\tau \alpha {N}_{\text{M}}+n\tau \right]}^{2}\hfill \\ \hfill & ={n}^{2}{I}_{x}(T=\tau ),\hfill \end{align} \tag{ 31 }$$

for every n = 1, 2, .... Therefore, without even computing the general bound, we can deduce from these simple arguments that it must be monotonic, show a secular quadratic increase in T, and be self-similar at intervals nτ. These qualitative features are exactly those of the SNR on figure 3. Furthermore, the maximum number of bosons, as shown in figure 2(c), scales like ${N}_{\text{M}}\propto \frac{1}{1-{g}^{2}}$. Therefore, the general bound prediction, including the ${N}_{\text{M}}^{2}$ prefactor, gives ${Q}_{\omega }\propto {\left(\frac{T}{1-{g}^{2}}\right)}^{2}$, which is what we observe in figure 3. Combining everything together, we now have the following picture: for short times, the QFI is dominated by non-universal terms, and the general bound is generally not saturated. For times 1 ≪ ωT ≲ ωτ, the effect of the gap is not yet relevant. The system behaves essentially as if it were quenched at the critical point, i.e. the photon number increases like N(t) ∼ t², and the bound scales as ${\left[{\int }_{0}^{T}N(t)\mathrm{d}t\right]}^{2}\propto {T}^{6}$. Finally, for durations T larger than τ, the system undergoes periodic oscillations, and the general bound becomes ${N}_{\text{M}}^{2}{T}^{2}\sim {\left(\frac{T}{1-{g}^{2}}\right)}^{2}$. Hence, we have shown that our bound allows to accurately grasp the various scaling regimes of the SNR for ωT  ≳  1.

The previous results have been obtained in the thermodynamic limit η  →  ∞. For η finite, the system evolves under the equations (1) and (2), and the evolution can no longer be exactly solved. Therefore, we resort to numerical simulations with a converged number of Fock basis to find the QFI. We consider a sudden quench to the critical point with η finite. The results are plotted on figure 5. In the limit η  →  ∞, we have a T⁶ scaling in the long-time limit, as already discussed. For η finite, however, we observe a transition from T⁶ to T², which takes place when ωT ∼ η^1/3.

Zoom In Zoom Out Reset image size

This behavior can be understood with the following argument. For finite η, a gap stabilizes around the critical point, of order ωη^−1/3. For short times, this finite gap is not relevant, and the system evolves as it would in the thermodynamic limit. Another formulation would be to say that, for small time, the number of photons is still small, and the quartic term in equation (2) is therefore negligible. However, for longer times, the finite gap will create a periodic revival behavior. This resembles the phenomenology in the thermodynamic limit when quenched away from the critical point. We can here invoke heuristic 1: the behavior for a quench at g = 1 for finite η is similar to the behavior one would obtain in the thermodynamic limit by quenching the system at ${g}^{\ast }=\sqrt{1-{\eta }^{-2/3}}$. Hence, we see that, although the state is now non-Gaussian and the bound (25) is no longer applicable, we have the same essential features as in the thermodynamic limit, with a periodic behavior for N, and a SNR showing a secular T² increase. The difference is that the period of the oscillations, and the boundaries between different scaling regimes, is now given by the parameter η, instead of the effective coupling g.

To evolve the state in a controlled way, a possible solution is to tune the parameters very slowly in time, in order to keep the evolution adiabatic. We will consider the thermodynamic limit η  →  ∞. We take (1), and slowly increase g toward the critical point. We will define = 1 − g² for convenience of notation. As long as the time scale introduced by the external driving $({\Delta}/\dot{{\Delta}})$ is much larger than the typical time scale of the system (1/Δ), i.e., as long as

$$\begin{equation}\frac{\dot{{\Delta}}}{{\Delta}}\ll {\Delta},\end{equation} \tag{ 32 }$$

the evolution will be adiabatic to very good approximation (see references [57–59] for time-dependent perturbation theory in this context, as well as [50]). Here Δ is the energy gap during the evolution, and $\dot{{\Delta}}$ its time derivative. Hence, when initialized in the ground state, the system will remain in it during the evolution. Note that ${\Delta}\propto \omega \sqrt{{\epsilon}}$, and $\dot{{\Delta}}/{\Delta}\sim \dot{{\epsilon}}/{\epsilon}$. In reference [22], a time-profile fulfilling these criteria was derived for , which reads as

$$\begin{equation}{\epsilon}(t)=\frac{1}{1+{(t/{\tau }_{Q})}^{2}},\end{equation} \tag{ 33 }$$

with τ_Q some time constant. During this non-linear ramp, we start from (0) = 1, i.e. g(0) = 0, and then we gradually approach the critical point = 0. The evolution speed $\dot{{\epsilon}}$ is high at first, then gradually decreases to keep up with the closure of the energy gap. More precisely, we have $\dot{{\Delta}}/{{\Delta}}^{2}=(t/\omega {\tau }_{Q}^{2})/\sqrt{1+{(t/{\tau }_{Q})}^{2}}$. Therefore, as long as ${\tau }_{Q}=\frac{1}{\varphi \omega }$, with φ ≪ 1, the criterion (32) will be satisfied, even in the thermodynamic limit, when the gap becomes exactly zero at the critical point. Note, however, that we only approach asymptotically the critical point, but we never reach it. The corresponding time-profile is schematically plotted in green in the left-hand side of figure 1.

Therefore, if we let the system evolve under this ramp for a time T, we expect the system to evolve adiabatically, and to be prepared in the ground-state of the Rabi Hamiltonian. That is, it will be in a squeezed state (7) with $\vert z\vert =-\frac{1}{4}\,\mathrm{log}({\epsilon}(T))$ and θ = 0. Using (17), the QFI can then be computed exactly as

$$\begin{equation}{I}_{\omega }=\frac{1}{8{\omega }^{2}}{\left(\frac{1}{1-{g}^{2}}\right)}^{2}=\frac{1}{8{\omega }^{2}{\epsilon}{(T)}^{2}}.\end{equation} \tag{ 34 }$$

Restoring the dependency of on T, and in the limit of large T, we find

$$\begin{equation}{Q}_{\omega }={\omega }^{2}{I}_{\omega }\sim {\left(\frac{T}{{\tau }_{Q}}\right)}^{4}={\left(\varphi \omega T\right)}^{4}.\end{equation} \tag{ 35 }$$

Therefore, the SNR in this case can scale quadratically in time. This scaling can also be understood through the general bound (25). For long t, $N\sim \frac{1}{\sqrt{{\epsilon}(t)}}\sim t/{\tau }_{Q}$. Since N increases roughly linearly in time, the squared integral in (25) scales as T⁴.

Figure 6 shows the scaling with T of the SNR, computed with exact simulation (see appendix D). We find that the T⁴ scaling is indeed recovered for ωT  ≳  10. The adiabatic behavior, however, is broken for very small T. Indeed, the fulfillment of condition (32) means that the population of excited states oscillates quickly, with an oscillation rate given by the energy gap, which here is of order ω. For T ≫ 1/ω, the system evolves over several periods, these population will be averaged out, and the system will indeed remain in its ground state. However, if we let the system evolve for a time shorter than the oscillation period 1/ω, the system can become excited, and the adiabatic prediction (35) is not satisfied anymore. In this regime, we observe that the SNR scales instead as T⁸. Since this behavior holds only for short times, and is associated with extremely small SNR, it should be irrelevant in practice.

Zoom In Zoom Out Reset image size

Let us stress that these results can also be expressed in terms of critical exponents. Let us assume that the gap scales like Δ ∼ ^zν , and the QFI scales like ${I}_{x}\sim {{\epsilon}}^{-\tilde{\gamma }}$. Then we can design a ramp ${\epsilon}(t)=\frac{1}{1+{(t/{\tau }_{Q})}^{1/(z\nu )}}$, which will satisfy the adiabaticity condition as long as τ_Q ≫ 1/ω. Plugging this expression in the QFI, we find ${I}_{x}\sim {T}^{\tilde{\gamma }/(z\nu )}$. In our case, we have zν = 1/2 and $\tilde{\gamma }=2$, which gives the T⁴ scaling. Finally, we also studied the evaluation of λ; we found that, for ωT  ≳  10, we recover the T⁴ scaling. As in the sudden quench case, the SNR becomes independent of the parameter being evaluated.

The previous ramp, although it optimally keeps the system in the ground state, may be challenging to implement in practice. In general, it might easier to implement ramps according to

$$\begin{equation}g(t)={g}_{\text{c}}\left(1-{\left(\frac{T-t}{T}\right)}^{r}\right),\end{equation} \tag{ 36 }$$

with 0 < r < ∞. The corresponding time-profile is plotted in red in the left-hand side of figure 1. In this subsection, we still assume that we are in the thermodynamic limit. Contrary to the previous case, the critical point g_c = 1 is reached in finite time, g(T) = g_c. However, the dynamics will cease to be adiabatic in the proximity of the QPT, which is at the core of the Kibble–Zurek mechanism [11, 49–53]. In particular, when $\dot{{\Delta}}\sim {{\Delta}}^{2}$ [50], the adiabaticity will be broken, which defines the so-called freeze-out time. In our case, this takes place at a time ${t}_{\text{f}}=T\left[1-{\left(\frac{\omega T}{r}\right)}^{-\frac{2}{r+2}}\right]$, which corresponds to a value

$$\begin{equation}1-{g}_{\text{f}}={\left(\frac{\omega T}{r}\right)}^{-2r/(2+r)}.\end{equation} \tag{ 37 }$$

In this case, the standard Kibble–Zurek argument [49–53] states that the evolution can be decomposed into two parts. First, an adiabatic evolution with the coupling moving from g(0) = 0 to g(t_f) = g_f. Second, an impulse regime in which the system cannot react to external changes imposed by g(t), and thus the system is effectively quenched from g_f to the final value g(T) = g_c = 1. Moreover, the QFI at the end of the evolution becomes approximately equal to the QFI at the freeze-out instant.

Combining this with equation (34), the KZ mechanism predicts a QFI proportional to $\frac{1}{{\omega }^{2}{(1-{g}_{\text{f}}^{2})}^{2}}$, which results in

$$\begin{equation}{Q}_{\omega }={\left(\frac{\omega T}{r}\right)}^{4r/(2+r)}.\end{equation} \tag{ 38 }$$

We computed the SNR under (36) (see appendix D) and compared it to this prediction. The results are plotted in figure 7. On the top panel, we show the scaling of the SNR for relatively small r. We observe that, for 10 ≲ ωT, the SNR scales indeed according to the KZ value. Note that the bound (25) is able to capture this scaling behavior. A closer examination, with higher values of r (bottom row of figure 7) reveals that there are actually three scaling regimes. For very short times ωT ≲ 10, Q_ω scales as T⁴. For 10 < ωT < r, one obtains Q_ω ∝ T⁶, while for ωT > r, we recover the Kibble–Zurek scaling (cf equation (38)). This can be interpreted as follows. For ωT < r, the freeze-out value 1 − g_f is larger than 1. Since g must be positive, this is not possible, which means that the adiabaticity is actually broken from the very beginning of the evolution. Therefore, the entire evolution can be deemed as a sudden quench, in which the coupling is instantaneously brought from g(0) = 0 to g_c = 1, and left to evolve for a time T (cf section 5); we then recover the T⁶ scaling we observed in section 5. For ωT > r, 0 ⩽ g_f ⩽ 1, so that the evolution can be decomposed according to the adiabatic-impulse approximation and the KZ argument holds (cf figure 7).

Zoom In Zoom Out Reset image size

In the limit r  →  ∞, the KZ prediction (38) leads to Q_ω ∝ T⁴. However, this scaling can only be attained for very long protocols, for ωT > r  →  ∞. Furthermore, the SNR carries a very small prefactor $\frac{1}{{r}^{4}}$. We can connect this to the results for the fully adiabatic ramp which we introduced in section 6.1, where Q_ω ∝ T⁴, with a prefactor φ⁴ ≪ 1. The critical point is not reached in the adiabatic ramp, but asymptotically approached for very long T. Therefore, we see that in the limit r ≫ 1 and ωT ≫ 1, the finite-time ramp and the fully adiabatic ones become similar. By contrast, for r ≫ 1 and 1 ≪ ωT ≪ r, the finite-time ramp has the same behavior as the sudden quench, as shown in the lower panel of figure 7. Hence, depending on the tuning of r and T, the non-linear ramp provides an interpolation between the simple linear ramp, the fully adiabatic evolution, and the sudden quench.

Finally, we study the adiabatic (33) and finite-time ramps (36) for finite-size system. The results are displayed in figure 8. For the adiabatic ramp at short times ωT ≪ η^1/3, we observe the same behavior as in the adiabatic limit, with a scaling going from T⁸ to T⁴. This can understood as follows: for finite η, the gap saturates around its minimum value in the critical region, of width Γ ∼ η^−2/3. Now, if we apply the evolution (33) for a total time ωT ≪ η^1/3, we will obtain at the end of the evolution ${\epsilon}(T)\sim {\left(\frac{{\tau }_{Q}}{T}\right)}^{2}\gg {\eta }^{-2/3}={\Gamma}$. In other words, at the end of the evolution, we will have 1 − g² ≫ Γ. Therefore, we remain safely out of the critical zone; the finite-size effects play a negligible role, and we recover the η  →  ∞ results. To the contrary, if ωT ≫ η^1/3, the adiabatic ramp brings us within the critical zone. In this region, the QFI saturates at a value ${I}_{\omega }\sim \frac{1}{{\omega }^{2}{{\Gamma}}^{2}}=\frac{1}{{\omega }^{2}}{\eta }^{4/3}$. This result can be retrieved using (34) and applying heuristic 1. This is depicted in the upper panels of figure 8. For ωT ≫ η^1/3, the SNR saturates at a value ∝η^4/3, and becomes independent of T. Note that this equivalent to perform the adiabatic ramp to get only up to ${g}^{\ast }=\sqrt{1-{\Gamma}}\sim \sqrt{1-{\eta }^{-2/3}}$. We stress that the same results hold for Q_λ .

Zoom In Zoom Out Reset image size

The results for the finite-time ramp are very similar. Let us consider again the freeze-out value (37). For $\frac{\omega T}{r}\ll {\eta }^{(2+r)/3r}$, it follows 1 − g_f ≫ Γ, i.e. the freeze-out occurs outside of the critical zone, and thus the thermodynamic limit results are recovered. For $\frac{\omega T}{r}\gg {\eta }^{(2+r)/3r}$, the SNR saturates at η^4/3. Those are the features we observe in the lower panels of figure 8, where we plotted the evolution for r = 1. First, we obtain the short-time scaling T⁴, then the KZ scaling T^4r/(2+r) = T^4/3, and finally the saturation with a T⁰ scaling. For larger r, there exists an intermediate region for 1 ≪ ωT ≪ r, where the SNR scales like T⁶.

To sum up, finite-size effects for a sudden quench (cf figure 5) and in adiabatic and finite-time ramps (cf figure 8) have a similar impact. For short durations, the finite-size effect are negligible, and we recover the thermodynamic limit scalings. For long time, the finite-size effects are dominant. The transition between these two regimes is governed by the minimum gap value, ωη^1/3. These various regimes are summarized in figure 1.

To conclude this section, let us discuss how the performances of the different strategies compare against each other. From figures 4 and 7, on the one hand, and figures 5 and 8, on the other hand, we can see that a sudden quench at g = 1 is always the optimal strategy, in the sense that it always gives the highest QFI for a given duration T [124]. This is particularly visible when η becomes large, in which case the precision achievable with the ramp quickly saturates when T increases, while the performances of the quench keep improving. Fully exploiting these improved performances, however, requires a somewhat more complex measurement strategy, as we will see in the next section.

Let us consider first the adiabatic ramp, in the thermodynamic limit. The system is prepared in its ground state, with a coupling value g very close to one. The natural observables in a bosonic system are the quadratures, which can be accessed by homodyne measurement. However, the quadratures always have zero average value in our case. Therefore, measuring $\hat{O}=\hat{x}$ or $\hat{O}=\hat{p}$ will always lead to $\langle \hat{O}\rangle =0$, and hence zero precision. Instead, we need to look at the fluctuations of the quadratures, by measuring ${\hat{x}}^{2}$ and ${\hat{p}}^{2}$. In practice, this could be made by performing an homodyne detection of one quadrature, and integrating the noise of the homodyne current. A natural choice is to measure the squeezed quadrature, by setting $\hat{O}={\hat{p}}^{2}$, this gives $\langle \hat{O}\rangle =\sqrt{1-{g}^{2}}$ (see equation (8)). We then get $\omega \,\frac{\mathrm{d}\langle \hat{O}\rangle }{\mathrm{d}\omega }=(\omega \,\frac{\mathrm{d}g}{\mathrm{d}\omega })\frac{\mathrm{d}\sqrt{1-{g}^{2}}}{\mathrm{d}g}\sim \frac{1}{\omega \sqrt{1-{g}^{2}}}$ [125]. The variance is derived using Wick's theorem, we get $\mathrm{Var}(\hat{O})=3{\langle {\hat{p}}^{2}\rangle }^{2}\sim 1-{g}^{2}$. Putting everything together, we get ${F}_{\omega }({\hat{p}}^{2})\sim {\left(\frac{1}{\omega (1-{g}^{2})}\right)}^{2}$. The precision has then exactly the same scaling in 1 − g² than the QFI obtained in equation (34), eventually, restoring the dependency of 1 − g² = (T) on T, as per (36), we find the same (ωT)⁴ scaling as the QFI. This argument indicates that measuring the noise in one quadrature allows to saturate the QFI. Our exact numerical simulations confirm this insight, as shown in figure 9. Here, the high precision has two origins: on the one hand, the derivative $\frac{\mathrm{d}\langle \hat{O}\rangle }{\mathrm{d}\omega }$ scales like ${(1-{g}^{2})}^{-1/2}$ and diverges as the critical point is approached, thus the observable becomes very sensitive. On the other hand, the noise $\mathrm{Var}(\hat{O})\sim 1-{g}^{2}$ is suppressed. However, $\hat{O}={\hat{p}}^{2}$ is not the only observable which saturates the QFI. Indeed, let us see what happens if we measure the anti-squeezed quadrature instead, and set $\hat{O}={\hat{x}}^{2}$. Then we find $\langle \hat{O}\rangle ={(1-{g}^{2})}^{-1/2}$, $\omega \,\frac{\mathrm{d}\langle \hat{O}\rangle }{\mathrm{d}\omega }\sim {(1-{g}^{2})}^{-3/2}$, and $\mathrm{Var}(\hat{O})\sim {(1-{g}^{2})}^{-1}$. We then obtain ${F}_{\omega }({\hat{x}}^{2})\sim {\left(\frac{1}{\omega (1-{g}^{2})}\right)}^{2}$. Hence, the measurement along the anti-squeezed quadrature gives the same scaling of precision as the measurement along the squeezed one, and also saturates the QFI. This is because, although the variance is much larger (and actually divergent) in this case, the derivative $\frac{\mathrm{d}\langle \hat{O}\rangle }{\mathrm{d}\omega }$ is correspondingly higher. Instead of homodyne detection, we may also measure the number of photons ${\hat{a}}^{{\dagger}}\hat{a}=\frac{{\hat{x}}^{2}+{\hat{p}}^{2}-1}{2}$. The photon number is essentially dominated by the anti-squeezed quadrature, and give the same FI. This is actually a more general property, coming from the critical scaling behavior. If an observable $\hat{O}$ scales like some power ${(1-{g}^{2})}^{\alpha }$ near the critical point, the derivative $\omega \,\frac{\mathrm{d}\langle \hat{O}\rangle }{\mathrm{d}\omega }\sim \frac{\mathrm{d}\langle \hat{O}\rangle }{\mathrm{d}g}$ will scale like ${(1-{g}^{2})}^{\alpha -1}$, and if the dynamics is (at least approximately) Gaussian, we will also have $\mathrm{Var}(\hat{O})\sim {\langle \hat{O}\rangle }^{2}\sim {(1-{g}^{2})}^{2\alpha }$. Combining these quantities, one finds a SNR which scales always like ${(1-{g}^{2})}^{-2}$, independently of the value of α.

Zoom In Zoom Out Reset image size

If we now consider the finite-time ramp, the above arguments remain essentially unchanged. As it can be seen in figure 9, homodyne measurement along both $\hat{x}$ and $\hat{p}$ quadratures gives a FI which increases with the ramp time T, and saturates the QFI in almost any case. We checked other quadratures and obtained every time the same scaling, with only a change in the prefactor. We found a single exception to that rule: for a linear ramp r = 1, measuring the squeezed quadrature ${\hat{p}}^{2}$ gives a precision which saturates to a finite value, and fail to reach the QFI. By contrast, measuring the anti-squeezed ${\hat{x}}^{2}$ always allows to saturate the QFI. Hence, we have a situation in which measuring the noisy quadrature is not just a good strategy, but actually a much better one than measuring the quadrature where noise is suppressed. Note that this is a very isolated case; if a measurement is performed along a slightly tilted direction, or if the ramp is slightly non-linear, the QFI scaling will be restored very fast. These results indicate, on the one end, that homodyne detection is almost always optimal when the ramp is used. On the other end, and counter-intuitively, they show that it is not necessary, or even counter-productive, to measure observables whose fluctuations are suppressed. Several works have considered using critical point as a tool to produce states with reduced quantum noise [12], which can then be used in standard phase-shift sensing protocols. Here, instead, the relevant parameters are encoded in the noise signature itself. Our protocol amplifies the noise, but with an amplification coefficient which depends critically on the parameter to be estimated.

To understand the FI achievable in the sudden quench case, we need to describe the dynamics of the quantity to be measured. In appendix F we provide the analytical expression for the time-evolution of ${\hat{x}}^{2}$, ${\hat{p}}^{2}$ and ${\hat{a}}^{{\dagger}}\hat{a}$, and the associated FI. The key findings are the following: measuring the photon number, or the noise along any quadrature, always yields a precision scaling with the quench duration like (ωT)⁴. We recall that the QFI scales like (ωT)⁶, and therefore, homodyne and photon-number measurement are always suboptimal. However, we also identified a family of observables which allows saturating the QFI. An example of observable saturating the QFI can be expressed as

$$\begin{equation}\hat{A}=-\frac{{\hat{x}}^{2}}{{(\omega T)}^{2}}+{\hat{p}}^{2}.\end{equation} \tag{ 39 }$$

This observable can be accessed via Gaussian operations and homodyne measurements, but it clearly requires a non-trivial measurement setup. These analytical results are confirmed by our numerical simulations. For example, in figure 9 we show the comparison between one- and multi-quadratures measurement strategies, and how the latter can saturate the QFI.

Let us conclude this section with two comments. First, although the sudden quench gives the best QFI, fully reaching this precision comes at the cost of a more complex measurement procedure. If, experimentally, only single-quadrature measurement or photon counting are available, the precision will fall short of the QFI. However, the sudden quench will still give a quartic scaling, while the finite-time ramp will give a KZ scaling ${T}^{\frac{2r}{r+2}}\leqslant {T}^{4}$. Therefore, even when the measurements are constrained, the quench is still the optimal strategy. The advantage will only be less important than in the absence of measurement constraints.

Second, we have focused the discussion on the measurement of ω. The results are unchanged if we want to measure λ instead. Once again, this is because a change in either ω or λ have the same effect, namely a renormalization of the effective coupling g (see appendix F for more details).

The Hamiltonian ${\hat{H}}_{x}$ generates a unitary evolution ${\hat{U}}_{x}(0\;\to \;t)$. The state evolves toward $\vert {\psi }_{t}\rangle ={\hat{U}}_{x}(0\;\to \;t)\vert {\psi }_{0}\rangle$.

The QFI can be rewritten as

$$\begin{equation}{I}_{x}(t)=4\left(\langle {\psi }_{0}\vert {F}_{x}{(t)}^{2}\vert {\psi }_{0}\rangle -\langle {\psi }_{0}\vert {F}_{x}(t){\vert {\psi }_{0}\rangle }^{2}\right),\end{equation} \tag{ C1 }$$

where we have defined

$$\begin{equation}{F}_{x}(t)=i{{\hat{U}}_{x}(0\;\to \;t)}^{{\dagger}}{\partial }_{x}\left({\hat{U}}_{x}(0\;\to \;t)\right).\end{equation} \tag{ C2 }$$

Now, it is straightforward to show that the derivative of F_x can be expressed as: ${\partial }_{t}{F}_{x}(t)={\hat{U}}_{x}{(0\;\to \;t)}^{{\dagger}}({\partial }_{x}{H}_{x}(t)){\hat{U}}_{x}(0\;\to \;t)$. Therefore,

$$\begin{equation}{F}_{x}(T)={\int }_{0}^{T}{\hat{U}}_{x}{(0\;\to \;t)}^{{\dagger}}({\partial }_{x}{H}_{x}(t)){\hat{U}}_{x}(0\;\to \;t)\mathrm{d}t.\end{equation} \tag{ C3 }$$

We want to compute the variance of F_x on the initial state |ψ₀⟩. For this, we can use the convexity of the standard deviation, i.e. $\sqrt{\mathrm{Var}(A+B)}\leqslant \sqrt{\mathrm{Var}(A)}+\sqrt{\mathrm{Var}(B)}$. From this, we can deduce

$$\begin{align}\hfill {I}_{x}(T)& \leqslant 4{\left({\int }_{0}^{T}\mathrm{d}t\sqrt{{\mathrm{Var}}_{{\psi }_{0}}\left({\hat{U}}_{x}^{{\dagger}}(0\;\to \;t)({\partial }_{x}{H}_{x}(t)){\hat{U}}_{x}(0\;\to \;t)\right)}\right)}^{2}\hfill \\ \hfill & =4{\left({\int }_{0}^{T}\mathrm{d}t\sqrt{{\mathrm{Var}}_{{\psi }_{t}}\left({\partial }_{x}{H}_{x}(t)\right)}\right)}^{2}.\hfill \end{align} \tag{ C4 }$$

In the first line, the variance is taken over the initial state |ψ₀⟩; in the second, it is taken over the time-evolved state |ψ_t ⟩. This equation is the most general bound we can set on the precision for a Hamiltonian evolution with a pure initial state. Again, we emphasize that this expression is valid for any Hamiltonian ${\hat{H}}_{x}$, even time-dependent. Note that the term under the integral is time-dependent in two different ways. First, the derivative of the Hamiltonian, ${\partial }_{x}{\hat{H}}_{x}(t)$, can be intrinsically time-dependent. Second, the variance is taken over the state |ψ_t ⟩, which evolves in time.

This expression provides a unified view of all the bounds discussed in the main text. Let us assume that we have ${\hat{H}}_{x}(t)=x\hat{A}(t)$, with $\hat{A}(t)={\partial }_{x}{\hat{H}}_{x}(t)$ independent of x, and ${\hat{H}}_{x}(t)$ is bounded. Then we can define the (time-dependent) extremal eigenvalues λ_M (t) and λ_m (t) of $\hat{A}(t)$. The variance can always be bounded by the difference between these eigenvalues: ${\mathrm{Var}}_{\chi }(\hat{A}(t))\leqslant \frac{1}{4}\vert {\lambda }_{M}(t)-{\lambda }_{m}(t)\vert$, for any state |χ⟩, and all t. Plugging this into (C4), we retrieve the bound (20)

$$\begin{equation}{I}_{x}(T)\leqslant {\left[{\int }_{0}^{T}\vert {\lambda }_{M}(t)-{\lambda }_{m}(t)\vert \right]}^{2}.\end{equation} \tag{ C5 }$$

This bound can be saturated when the state is in a coherent superposition of the two extremal values: $\vert \chi (t)\rangle =\frac{\vert {\lambda }_{M}(t)\rangle +\vert {\lambda }_{m}(t)\rangle }{\sqrt{2}}$ [119, 121]. This condition can be challenging to implement in practice, especially if the Hamiltonian is time-dependent; however, it was shown that this could be done for small systems, using quantum control [121].

Let us now assume that all the conditions (1)–(4) discussed in the main text are satisfied. Then the extremal eigenvalues become time-independent; and because ${\hat{H}}_{x}$ is local, they can scale at most with the number of components. This means we have |λ_M (t) − λ_m (t)| = αN, with α a time-independent, non-universal constant which depends on the precise expression of ${\hat{H}}_{x}$. This finally gives us back the Heisenberg limit

$$\begin{align}\hfill {I}_{x}(T)& \leqslant {\left[{\int }_{0}^{T}\alpha N\,\mathrm{d}t\right]}^{2}\hfill \\ \hfill & \leqslant {\alpha }^{2}{N}^{2}{T}^{2}.\hfill \end{align} \tag{ C6 }$$

Conversely, if ${\hat{H}}_{x}$ acts on k particles at the same time, then the eigenvalues can scale like N^k , and we retrieve the bounds (19). Finally, if we let ${\hat{H}}_{x}$ be time-dependent, but still bounded, we fall back to the discussion of Pang and Jordan [121].

We will now show derive our bound (25). The proof contains essentially two parts: the first part is to use (C4), thus we reduce the estimation of the QFI to evaluating the variance of ∂_x H_x at each time. The second step is to notice, using Wick's theorem, that the variance of a quadratic Hamiltonian on a Gaussian state is essentially equal to ${N}^{2}={\langle \hat{N}\rangle }^{2}$. Combining the two, we find that the QFI should be bounded by a quantity like ${\left({\int }_{0}^{T}N(t)\right)}^{2}$. More precisely, we assume that both ${\hat{H}}_{x}$ and its derivative ${\partial }_{x}{\hat{H}}_{x}$ are purely quadratic, with no linear part,

$$\begin{equation}{\partial }_{x}{H}_{x}=(\hat{x}\hat{p}){M}_{x}{(\hat{x}\hat{p})}^{\text{T}}.\end{equation} \tag{ C7 }$$

We also assume that, in the initial state, we have $\langle \hat{x}\rangle =\langle \hat{p}\rangle =0$, which is the case for a vacuum state. We will first assume that ${\hat{H}}_{x}$ is time-independent, and relax this assumption at the end. Starting with (C4), we need to compute the variance of $\hat{O}(t)={\hat{U}}_{x}{(0\;\to \;t)}^{{\dagger}}{\partial }_{x}{H}_{x}{\hat{U}}_{x}(0\;\to \;t)$ over the initial state. We can rewrite:

$$\begin{equation*}\hat{O}(t)=[\hat{x}(t)\hat{p}(t)]{M}_{x}\left[\begin{matrix}\hfill \hat{x}(t)\hfill \\ \hfill \hat{p}(t)\hfill \end{matrix}\right].\end{equation*}$$

Here M_x is a two-by-two Hermitian matrix, and $\hat{x}(t)={\hat{U}}_{x}{(0\;\to \;t)}^{{\dagger}}\hat{x}{\hat{U}}_{x}(0\;\to \;t)$ is the time-evolved quadrature in Heisenberg picture. Because the Hamiltonian is purely quadratic, the quadratures stay centered around zero at all times: $\langle \hat{x}(t)\rangle =\langle \hat{p}(t)\rangle =0$ for all t. Now, at any given t, we can define two quadratures ${({\hat{x}}_{m}(t),{\hat{p}}_{m}(t))}^{\text{T}}=R{(\hat{x}(t),\hat{p}(t))}^{\text{T}}$, which are obtained from the original quadratures by a (time-dependent) rotation R(t), and which satisfy: $\langle :\;{\hat{x}}_{m}{\hat{p}}_{m}\;:\rangle =0$ (where :: means we are taking the symmetric combination $:\;\hat{a}\hat{b}{:=}\frac{\hat{a}\hat{b}+\hat{b}\hat{a}}{2}$). ${\hat{x}}_{m}$ and ${\hat{p}}_{m}$ simply correspond to the directions of maximal and minimum squeezing, respectively. Note that these two quadratures can always be defined, even when the state is mixed. Then we can absorb the rotation in the matrix M_x , and we obtain

$$\begin{equation*}\hat{O}(t)=[\hat{x}(t)\hat{p}(t)]{M}_{x}\left[\begin{matrix}\hfill \hat{x}(t)\hfill \\ \hfill \hat{p}(t)\hfill \end{matrix}\right]=[{\hat{x}}_{m}(t){\hat{p}}_{m}(t)]R(t){M}_{x}R{(t)}^{\text{T}}\left[\begin{matrix}\hfill {\hat{x}}_{m}(t)\hfill \\ \hfill {\hat{p}}_{m}(t),\hfill \end{matrix}\right]\end{equation*}$$

with the rotated matrix

$$\begin{equation*}R(t){M}_{x}R{(t)}^{\text{T}}=\left[\begin{matrix}\hfill {A}_{1}(t)\hfill & \hfill b(t)+ic(t)\hfill \\ \hfill b(t)-ic(t)\hfill & \hfill {A}_{2}(t).\hfill \end{matrix}\right]\end{equation*}$$

In other words, the time-evolution has two different components. It rotates the direction of squeezing, and changes its amount. The second operation is now encoded in the quadrature ${\hat{x}}_{m}$, while the first is contained in the matrix R(t)M_x R(t)^T. Now, the coefficients A_i , b and c can have a very complicated time-dependence in general. However, since R is just a rotation, the eigenstates of M_x , which we call ϕ and χ, remain the same at all time. The goal is now to find a bound on $\mathrm{Var}(\hat{O})$ which depends only on this time-independent factors. Let us expand $\hat{O}(t)$

$$\begin{equation}\hat{O}(t)={A}_{1}{\hat{x}}_{m}^{2}+{A}_{2}{\hat{p}}_{m}^{2}+2b(:\;{\hat{x}}_{m}{\hat{p}}_{m}\;:)-c.\end{equation} \tag{ C8 }$$

The coefficient c is a constant with no effect on the variance, and can be dropped. We can develop

$$\begin{align}\hfill {\hat{O}}^{2}& ={A}_{1}^{2}{\hat{x}}_{m}^{4}+{A}_{2}^{2}{\hat{p}}_{m}^{4}+{A}_{1}{A}_{2}({\hat{x}}_{m}^{2}{\hat{p}}_{m}^{2}+{\hat{p}}_{m}^{2}{\hat{x}}_{m}^{2})+{b}^{2}{({\hat{x}}_{m}{\hat{p}}_{m}+{\hat{p}}_{m}{\hat{x}}_{m})}^{2}+2b(2{A}_{1}:{\hat{x}}_{m}^{3}{\hat{p}}_{m}:+\;2{A}_{2}:{\hat{x}}_{m}{\hat{p}}_{m}^{3}:)\hfill \\ \hfill & ={A}_{1}^{2}{\hat{x}}_{m}^{4}+{A}_{2}^{2}{\hat{p}}_{m}^{4}+(2{A}_{1}{A}_{2}+4{b}^{2}):{\hat{x}}_{m}^{2}{\hat{p}}_{m}^{2}:+({b}^{2}-{A}_{1}{A}_{2})+2b(2{A}_{1}:{\hat{x}}_{m}^{3}{\hat{p}}_{m}:+\;2{A}_{2}:{\hat{x}}_{m}{\hat{p}}_{m}^{3}:),\hfill \end{align} \tag{ C9 }$$

where we have dropped the explicit time-dependence to lighten notation. In the second line, we have used the handy relation: ${\hat{x}}_{m}^{2}{\hat{p}}_{m}^{2}+{\hat{p}}_{m}^{2}{\hat{x}}_{m}^{2}=2:{\hat{x}}_{m}^{2}{\hat{p}}_{m}^{2}:-1$. Now, when we take the average over the initial state, we may use Wick's theorem $\langle {\hat{x}}_{m}{(t)}^{4}\rangle =3{\langle {\hat{x}}_{m}{(t)}^{2}\rangle }^{2}$, $\langle :\;{\hat{x}}_{m}{(t)}^{2}{\hat{p}}_{m}{(t)}^{2}\;:\rangle =\langle {\hat{x}}_{m}{(t)}^{2}\rangle \langle {\hat{p}}_{m}{(t)}^{2}\rangle +2{\langle :\;{\hat{x}}_{m}(t){\hat{p}}_{m}(t)\;:\rangle }^{2}=\langle {\hat{x}}_{m}{(t)}^{2}\rangle \langle {\hat{p}}_{m}{(t)}^{2}\rangle$, and $\langle :\;{\hat{x}}_{m}{(t)}^{3}{\hat{p}}_{m}(t)\;:\rangle =\langle :\;{\hat{x}}_{m}(t){\hat{p}}_{m}{(t)}^{3}\;:\rangle =0$, at all times. Then after straightforward manipulation, we find

$$\begin{equation}\langle {\hat{O}}^{2}(t)\rangle -{\langle \hat{O}(t)\rangle }^{2}=2({A}_{1}^{2}{\langle {\hat{x}}_{m}^{2}\rangle }^{2}+{A}_{2}^{2}{\langle {\hat{p}}_{m}^{2}\rangle }^{2})+4\langle {\hat{x}}_{m}^{2}\rangle \langle {\hat{p}}_{m}^{2}\rangle {b}^{2}-{A}_{1}{A}_{2}+{b}^{2},\end{equation} \tag{ C10 }$$

where all the averages are taken over the initial state. Now, we want to replace the time-dependent components A_i , b and c by the time-independent χ, ϕ. We have the following relations:

$$\begin{align*}\hfill \chi \phi & ={A}_{1}{A}_{2}-({b}^{2}+{c}^{2})\leqslant {A}_{1}{A}_{2}-{b}^{2}\hfill \\ \hfill \chi +\phi & ={A}_{1}+{A}_{2}\hfill \\ \hfill {(\phi -\chi )}^{2}& ={({A}_{1}-{A}_{2})}^{2}+4({b}^{2}+{c}^{2})\geqslant {({A}_{1}-{A}_{2})}^{2}\hfill \\ \hfill {\phi }^{2}+{\chi }^{2}& \geqslant {A}_{1}^{2}+{A}_{2}^{2}.\hfill \end{align*}$$

which gives

$$\begin{align*}\hfill \langle {\hat{O}}^{2}(t)\rangle -{\langle \hat{O}(t)\rangle }^{2}& =2({A}_{1}^{2}{\langle {\hat{x}}_{m}^{2}\rangle }^{2}+{A}_{2}^{2}{\langle {\hat{p}}_{m}^{2}\rangle }^{2})+4\langle {\hat{x}}_{m}^{2}\rangle \langle {\hat{p}}_{m}^{2}\rangle {b}^{2}-{A}_{1}{A}_{2}+{b}^{2}\hfill \\ \hfill & \leqslant 2\left({A}_{1}^{2}{\langle {\hat{x}}_{m}^{2}\rangle }^{2}+{A}_{2}^{2}{\langle {\hat{p}}_{m}^{2}\rangle }^{2}\right)+4\langle {\hat{x}}_{m}^{2}\rangle \langle {\hat{p}}_{m}^{2}\rangle ({b}^{2}+{c}^{2})-{A}_{1}{A}_{2}+({b}^{2}+{c}^{2})\hfill \\ \hfill & =2{\left({A}_{1}\langle {\hat{x}}_{m}^{2}\rangle +{A}_{2}\langle {\hat{p}}_{m}^{2}\rangle \right)}^{2}+(1+4\langle {\hat{x}}_{m}^{2}\rangle \langle {\hat{p}}_{m}^{2}\rangle )({b}^{2}+{c}^{2}-{A}_{1}{A}_{2})\hfill \\ \hfill & =2{\left({A}_{1}\langle {\hat{x}}_{m}^{2}\rangle +{A}_{2}\langle {\hat{p}}_{m}^{2}\rangle \right)}^{2}-\chi \phi (1+4\langle {\hat{x}}_{m}^{2}\rangle \langle {\hat{p}}_{m}^{2}\rangle ).\hfill \end{align*}$$

This bound is saturated whenever c = 0; that is, when ∂_x H_x has no constant term. Now, we can write in general

$$\begin{equation*}{(ac+bd)}^{2}-{(ad+bc)}^{2}\leqslant {(ac+bd)}^{2}\leqslant {(ac+bd)}^{2}+{(ad+bc)}^{2},\end{equation*}$$

so that

$$\begin{equation*}({a}^{2}-{b}^{2})({c}^{2}-{d}^{2})\leqslant {(ac+bd)}^{2}\leqslant \frac{{(a-b)}^{2}{(c-d)}^{2}+{(a+b)}^{2}{(c+d)}^{2}}{2},\end{equation*}$$

for all $a,b,c,\mathrm{d}\in \mathbb{R}$. We can apply this to the first term in the previous equation, and we get

$$\begin{align*}\hfill \langle {\hat{O}}^{2}(t)\rangle -{\langle \hat{O}(t)\rangle }^{2}& \leqslant {({A}_{1}-{A}_{2})}^{2}{\left(\langle {\hat{x}}_{m}^{2}\rangle -\langle {\hat{p}}_{m}^{2}\rangle \right)}^{2}+{({A}_{1}+{A}_{2})}^{2}{\left(\langle {\hat{x}}_{m}^{2}\rangle +\langle {\hat{p}}_{m}^{2}\rangle \right)}^{2}-\chi \phi \left(1+4\langle {\hat{x}}_{m}^{2}\rangle \langle {\hat{p}}_{m}^{2}\rangle \right)\hfill \\ \hfill & ={({A}_{1}-{A}_{2})}^{2}{\left(\langle {\hat{x}}_{m}^{2}\rangle -\langle {\hat{p}}_{m}^{2}\rangle \right)}^{2}+{(\phi +\chi )}^{2}{\left(\langle {\hat{x}}_{m}^{2}\rangle +\langle {\hat{p}}_{m}^{2}\rangle \right)}^{2}-\chi \phi \left(1+4\langle {\hat{x}}_{m}^{2}\rangle \langle {\hat{p}}_{m}^{2}\rangle \right)\hfill \\ \hfill & \leqslant {(\chi -\phi )}^{2}{\left(\langle {\hat{x}}_{m}^{2}\rangle -\langle {\hat{p}}_{m}^{2}\rangle \right)}^{2}+{(\phi +\chi )}^{2}{\left(\langle {\hat{x}}_{m}^{2}\rangle +\langle {\hat{p}}_{m}^{2}\rangle \right)}^{2}-\chi \phi \left(1+4\langle {\hat{x}}_{m}^{2}\rangle \langle {\hat{p}}_{m}^{2}\rangle \right).\hfill \end{align*}$$

Now, we have successfully eliminated all the time-dependent coefficients A_i , b and c, and we are left only with the time-independent eigenvalues ϕ and χ. There is still, however, some time-dependence left in the expectation values of the quadratures. We will now find how this bound can be amended to a more elegant form, involving only the photon number at a given time. The expression above can be massaged as

$$\begin{align*}\hfill \langle {O}^{2}(t)\rangle -{\langle O(t)\rangle }^{2}& \leqslant ({\chi }^{2}+{\phi }^{2})\left[{(\langle {\hat{x}}_{m}^{2}\rangle -\langle {\hat{p}}_{m}^{2}\rangle )}^{2}+{(\langle {\hat{x}}_{m}^{2}\rangle +\langle {\hat{p}}_{m}^{2}\rangle )}^{2}\right]\hfill \\ \hfill & \quad +\chi \phi \left(2{(\langle {\hat{x}}_{m}^{2}\rangle +\langle {\hat{p}}_{m}^{2}\rangle )}^{2}-2{(\langle {\hat{x}}_{m}^{2}\rangle -\langle {\hat{p}}_{m}^{2}\rangle )}^{2}-1-4\langle {\hat{x}}_{m}^{2}\rangle \langle {\hat{p}}_{m}^{2}\rangle \right)\hfill \\ \hfill & =2({\chi }^{2}+{\phi }^{2})[{\langle {\hat{x}}_{m}^{2}\rangle }^{2}+{\langle {\hat{p}}_{m}^{2}\rangle }^{2}]+\chi \phi [4\langle {\hat{x}}_{m}^{2}\rangle \langle {\hat{p}}_{m}^{2}\rangle -1].\hfill \end{align*}$$

The term $4\langle {\hat{x}}_{m}^{2}\rangle \langle {\hat{p}}_{m}^{2}\rangle -1$ is always positive, due to Heisenberg's inequality. Now, we can distinguish two cases: if χϕ ⩾ 0, we use the triangle inequality $\phi \chi \leqslant \frac{{\chi }^{2}+{\phi }^{2}}{2}$, and we find a new bound

$$\begin{align*}\hfill \langle {\hat{O}}^{2}(t)\rangle -{\langle \hat{O}(t)\rangle }^{2}& \leqslant 2({\phi }^{2}+{\chi }^{2})\left[{\left(\langle {\hat{x}}_{m}^{2}\rangle +\langle {\hat{p}}_{m}^{2}\rangle \right)}^{2}-\langle {\hat{x}}_{m}^{2}\rangle \langle {\hat{p}}_{m}^{2}\rangle -\frac{1}{4}\right]\hfill \\ \hfill & \leqslant 2({\phi }^{2}+{\chi }^{2})\left[{\left(\langle {\hat{x}}_{m}^{2}\rangle +\langle {\hat{p}}_{m}^{2}\rangle \right)}^{2}-\frac{1}{2}\right].\hfill \end{align*}$$

where, in the last line, we have used Heisenberg's inequality $\langle {\hat{x}}_{m}^{2}\rangle \langle {\hat{p}}_{m}^{2}\rangle \geqslant \frac{1}{4}$. Recall that $\langle \hat{x}\rangle =\langle \hat{p}\rangle =0$. If χϕ ⩽ 0 instead, we have directly

$$\begin{align*}\hfill \langle {\hat{O}}^{2}(t)\rangle -{\langle \hat{O}(t)\rangle }^{2}& \leqslant 2({\chi }^{2}+{\phi }^{2})[{\langle {x}_{m}^{2}\rangle }^{2}+{\langle {p}_{m}^{2}\rangle }^{2}]\hfill \\ \hfill & =2({\phi }^{2}+{\chi }^{2})\left[{\left(\langle {x}_{m}^{2}\rangle +\langle {p}_{m}^{2}\rangle \right)}^{2}-2\langle {x}_{m}^{2}\rangle \langle {p}_{m}^{2}\rangle \right]\hfill \\ \hfill & \leqslant 2({\phi }^{2}+{\chi }^{2})\left[{\left(\langle {x}_{m}^{2}\rangle +\langle {p}_{m}^{2}\rangle \right)}^{2}-\frac{1}{2}\right].\hfill \end{align*}$$

Hence, in both cases, we find the same bound. It is now easy to rewrite everything in terms of the photon number; at each time, we have $N(t)=\langle \frac{{\hat{x}}_{m}^{2}(t)+{\hat{p}}_{m}^{2}(t)-1}{2}\rangle$, and therefore

$$\begin{equation*}\langle {\hat{O}}^{2}(t)\rangle -{\langle \hat{O}(t)\rangle }^{2}\leqslant 2({\phi }^{2}+{\chi }^{2})\left[{\left(2N(t)+1\right)}^{2}-\frac{1}{2}\right].\end{equation*}$$

Finally, combining this with the expression (C4), we find

$$\begin{align}\hfill {I}_{x}(t)& \leqslant 8({\chi }^{2}+{\phi }^{2}){\left({\int }_{0}^{T}\mathrm{d}t\sqrt{{(2N(t)+1)}^{2}-\frac{1}{2}}\right)}^{2}\hfill \\ \hfill & \leqslant 8({\chi }^{2}+{\phi }^{2}){\left[{\int }_{0}^{T}\mathrm{d}t\left(2N(t)+1\right)\right]}^{2}.\hfill \end{align} \tag{ C11 }$$

This derivation was performed assuming that ${\hat{H}}_{x}$ was time-independent. Now, if ${\hat{H}}_{x}$ is explicitly time-dependent, we can still define a matrix M(t) and its two eigenvalues ϕ(t) and χ(t). Those are also the eigenvalues of the rotated matrix R(t)M(t)R(t)^T, at any given time t. The steps of the derivation remain exactly the same, and we find instead

$$\begin{equation}{I}_{x}(T)\leqslant 8{\left[{\int }_{0}^{T}\mathrm{d}t\sqrt{\chi {(t)}^{2}+\phi {(t)}^{2}}\left(2N(t)+1\right)\right]}^{2}.\end{equation} \tag{ C12 }$$

Note that although χ and ϕ are now time-dependent, there are still entirely determined by ${\hat{H}}_{x}$. Thus, in this case too, we have successfully separated the QFI into a contribution which depends on the Hamiltonian only, and one which depends only on the average photon number.

We will now show how the bound above can be extended when the state is still Gaussian, but with a non-zero displacement. We now assume that the initial state is an arbitrary pure Gaussian state, and the Hamiltonian now has the form

$$\begin{equation}{\hat{H}}_{x}=(\hat{x}\hat{p}){h}_{x}{(\hat{x}\hat{p})}^{\text{T}}+u{(\hat{x}\hat{p})}^{\text{T}}.\end{equation} \tag{ C13 }$$

For now, we will assume that u is independent of x. Then the derivative ${\partial }_{x}{\hat{H}}_{x}$ is still purely quadratic, but the state has now a non-zero displacement which may change in time. We can apply once again apply some time-dependent rotation ${({\hat{x}}_{m},{\hat{p}}_{m})}^{\text{T}}=R(t){(\hat{x},\hat{p})}^{\text{T}}$, so that the rotated quadratures verify the following property: ${\hat{x}}_{m}=\hat{X}+\alpha$ and ${\hat{p}}_{m}=\hat{P}+\beta$, with α and β scalars, and we have the following properties: $\langle {\hat{X}}^{n}\rangle =\langle {\hat{P}}^{n}\rangle =0$ for odd n, and $\langle :\hat{X}\hat{P}\;:\rangle =0$. In other words, what we have done is apply first a rotation, then a displacement, to obtain quadratures centered around zero and aligned with the squeezing. Then we can again absorb the rotation inside M_x , and we obtain again equation (C8). If we now develop ${\hat{x}}_{m}=\hat{X}+\alpha$, we find

$$\begin{equation}\hat{O}(t)={A}_{1}{\hat{X}}^{2}+{A}_{2}{\hat{P}}^{2}+2b(:\hat{X}\hat{P}\;:)+2({A}_{1}\alpha +b\beta )\hat{X}+2({A}_{2}\beta +b\alpha )\hat{P}+({A}_{1}{\alpha }^{2}+{A}_{2}{\beta }^{2}+2b\alpha \beta -c).\end{equation} \tag{ C14 }$$

The last term is again a scalar (albeit time-dependent) with no effect on the variance, so it can be safely neglected. Let us write μ₁ = 2(A₁ α + bβ) and μ₂ = 2(A₂ β + bα). Then we can develop

$$\begin{align}\hfill {\hat{O}}^{2}& ={A}_{1}^{2}{\hat{X}}^{4}+{A}_{2}^{2}{\hat{P}}^{4}+(2{A}_{1}{A}_{2}+4{b}^{2}):{\hat{X}}^{2}{\hat{P}}^{2}:+({b}^{2}-{A}_{1}{A}_{2})+2b(2{A}_{1}:{\hat{X}}^{3}\hat{P}:+2{A}_{2}:\hat{X}{\hat{P}}^{3}:)+{\mu }_{1}^{2}{\hat{X}}^{2}+{\mu }_{2}^{2}{\hat{P}}^{2}\hfill \\ \hfill & \quad +2{\mu }_{1}{\mu }_{2}:\hat{X}\hat{P}:+2{\mu }_{1}{A}_{1}{\hat{X}}^{3}+2{\mu }_{1}{A}_{2}:\hat{X}{\hat{P}}^{2}:+4{\mu }_{1}b:{\hat{X}}^{2}\hat{P}:+2{\mu }_{2}{A}_{2}{\hat{P}}^{3}+2{\mu }_{2}{A}_{1}:{\hat{X}}^{2}\hat{P}:+4{\mu }_{2}b:\hat{X}{\hat{P}}^{2}:.\hfill \end{align} \tag{ C15 }$$

The first line is just the expression (C9) where we replaced ${\hat{x}}_{m}$ and ${\hat{p}}_{m}$ by $\hat{X}$ and $\hat{P}$, and the second line is an extra term which disappears when α = β = 0. Now we can apply Wick's theorem and exploit the relations $:\hat{X}\hat{P}{:=}0$ etc. Restoring the expression of μ₁ and μ₂, we finally obtain the expression of the variance

$$\begin{align}\hfill \langle {\hat{O}}^{2}(t)\rangle -{\langle \hat{O}(t)\rangle }^{2}& =2({A}_{1}^{2}{\langle {\hat{X}}^{2}\rangle }^{2}+{A}_{2}^{2}{\langle {\hat{P}}^{2}\rangle }^{2})+4\langle {\hat{X}}^{2}\rangle \langle {\hat{P}}^{2}\rangle {b}^{2}-{A}_{1}{A}_{2}+{b}^{2}\hfill \\ \hfill & \quad +4{({A}_{1}\alpha +b\beta )}^{2}\langle {\hat{X}}^{2}\rangle +4{({A}_{2}\beta +b\alpha )}^{2}\langle {\hat{P}}^{2}\rangle .\hfill \end{align} \tag{ C16 }$$

This expression is just (C10), plus an extra term which couples the variance of the quadratures and the displacement. The first term captures the effect of squeezing, the second the interplay between squeezing and displacement. We now make the following manipulation: we integrate again the displacement in the first term (noting that $\langle {\hat{x}}_{m}^{2}\rangle =\langle {\hat{X}}^{2}\rangle +{\alpha }^{2}$, and the same for ${\hat{p}}_{m}$), but not in the second term. We obtain

$$\begin{align}\hfill \langle {\hat{O}}^{2}(t)\rangle -{\langle \hat{O}(t)\rangle }^{2}& =2({A}_{1}^{2}{\langle {\hat{x}}_{m}^{2}\rangle }^{2}+{A}_{2}^{2}{\langle {\hat{p}}_{m}^{2}\rangle }^{2})+4\langle {\hat{x}}_{m}^{2}\rangle \langle {\hat{p}}_{m}^{2}\rangle {b}^{2}-{A}_{1}{A}_{2}+{b}^{2}\hfill \\ \hfill & \quad +8b\alpha \beta ({A}_{1}\langle {\hat{X}}^{2}\rangle +{A}_{2}\langle {\hat{P}}^{2}\rangle )-(4{\alpha }^{2}{\beta }^{2}{b}^{2}+2{A}_{1}^{2}{\alpha }^{4}+2{A}_{2}^{2}{\beta }^{4}).\hfill \end{align} \tag{ C17 }$$

The first term can now be bounded, in exactly the same way as in the previous subsection, by $2({\phi }^{2}+{\chi }^{2})[{(2N(t)+1)}^{2}-\frac{1}{2}]$ noting that $N=\frac{\langle {\hat{x}}_{m}^{2}+{\hat{p}}_{m}^{2}\rangle -1}{2}$. For the second term, we can bound it by applying the Cauchy–Schwarz and triangle inequalities

$$\begin{align*}\hfill 8b\alpha \beta ({A}_{1}\langle {\hat{X}}^{2}\rangle +{A}_{2}\langle {\hat{P}}^{2}\rangle )& \leqslant 8\alpha \beta b\sqrt{{A}_{1}^{2}+{A}_{2}^{2}}\sqrt{{\langle {\hat{X}}^{2}\rangle }^{2}+{\langle {\hat{P}}^{2}\rangle }^{2}}\hfill \\ \hfill & \leqslant 8\alpha \beta b\sqrt{{A}_{1}^{2}+{A}_{2}^{2}}(\langle {\hat{X}}^{2}\rangle +\langle {\hat{P}}^{2}\rangle )\hfill \\ \hfill & \leqslant 2({b}^{2}+{A}_{1}^{2}+{A}_{2}^{2})(\langle {\hat{X}}^{2}\rangle +\langle {\hat{P}}^{2}\rangle )({\alpha }^{2}+{\beta }^{2})\hfill \\ \hfill & \leqslant \frac{1}{2}({b}^{2}+{A}_{1}^{2}+{A}_{2}^{2}){(\langle {\hat{X}}^{2}\rangle +\langle {\hat{P}}^{2}\rangle +{\alpha }^{2}+{\beta }^{2})}^{2}\hfill \\ \hfill & =\frac{1}{2}({b}^{2}+{A}_{1}^{2}+{A}_{2}^{2}){(\langle {\hat{x}}_{m}^{2}\rangle +\langle {\hat{p}}_{m}^{2}\rangle )}^{2}.\hfill \end{align*}$$

And finally, noting that ${b}^{2}+{A}_{1}^{2}+{A}_{2}^{2}\leqslant 2({b}^{2}+{c}^{2})+{A}_{1}^{2}+{A}_{2}^{2}={\chi }^{2}+{\phi }^{2}$, we can bound the variance by

$$\begin{align}\hfill \langle {\hat{O}}^{2}(t)\rangle -{\langle \hat{O}(t)\rangle }^{2}& \leqslant 2({\phi }^{2}+{\chi }^{2})\left[{(2N(t)+1)}^{2}-\frac{1}{2}\right]+\frac{1}{2}({\phi }^{2}+{\chi }^{2})[{(2N(t)+1)}^{2}]\hfill \\ \hfill & \quad -(4{\alpha }^{2}{\beta }^{2}{b}^{2}+2{A}_{1}^{2}{\alpha }^{4}+2{A}_{2}^{2}{\beta }^{4})\hfill \\ \hfill & \leqslant \frac{5}{2}({\phi }^{2}+{\chi }^{2})[{(2N(t)+1)}^{2}],\hfill \end{align} \tag{ C18 }$$

which is very similar to our previous bound, but with a different prefactor.

Finally, let us consider an even more general case, when the linear term in ${\hat{H}}_{x}$, u, can now depend on x. We define v = ∂_x u. We can again apply a quadrature rotation, and we get ${\partial }_{x}{\hat{H}}_{x}=({\hat{x}}_{m}{\hat{p}}_{m}){M}_{x}{({\hat{x}}_{m}{\hat{p}}_{m})}^{\text{T}}+v{R}^{-1}{({\hat{x}}_{m}{\hat{p}}_{m})}^{\text{T}}$. The rotated vector vR⁻¹ can be written as (d₁, d₂); although the expression of these two terms is unknown, we know we have ${d}_{1}^{2}+{d}_{2}^{2}=\vert v{\vert }^{2}$, since the rotation preserves the norm. This term adds additional linear $\hat{X}$ and $\hat{P}$ term. We find that the expressions (C14) and (C15) will remain the same, but we now have μ₁ = 2(A₁ α + bβ + d₁) and μ₂ = 2(A₂ β + bα + d₂). Integrating the displacement in the first, but not the second term, we get

$$\begin{align}\hfill \langle {\hat{O}}^{2}(t)\rangle -{\langle \hat{O}(t)\rangle }^{2}& =2({A}_{1}^{2}{\langle {\hat{x}}_{m}^{2}\rangle }^{2}+{A}_{2}^{2}{\langle {\hat{p}}_{m}^{2}\rangle }^{2})+4\langle {\hat{x}}_{m}^{2}\rangle \langle {\hat{p}}_{m}^{2}\rangle {b}^{2}-{A}_{1}{A}_{2}+{b}^{2}+8b\alpha \beta ({A}_{1}\langle {\hat{X}}^{2}\rangle +{A}_{2}\langle {\hat{P}}^{2}\rangle )\hfill \\ \hfill & \quad +4{d}_{1}({A}_{1}\alpha +b\beta )\langle {\hat{X}}^{2}\rangle +4{d}_{2}({A}_{2}\beta +b\alpha )\langle {\hat{P}}^{2}\rangle +{d}_{1}^{2}\langle {\hat{X}}^{2}\rangle +{d}_{2}^{2}\langle {\hat{P}}^{2}\rangle \hfill \\ \hfill & \quad -(4{\alpha }^{2}{\beta }^{2}{b}^{2}+2{A}_{1}^{2}{\alpha }^{4}+2{A}_{2}^{2}{\beta }^{4}).\hfill \end{align} \tag{ C19 }$$

The first terms can be bounded as before. For the other terms, we find using Cauchy–Schwarz and triangle inequalities

$$\begin{align*}\hfill {d}_{1}({A}_{1}\alpha +b\beta )\langle {\hat{X}}^{2}\rangle +{d}_{2}({A}_{2}\beta +b\alpha )\langle {\hat{P}}^{2}\rangle & \leqslant \sqrt{{d}_{1}^{2}+{d}_{2}^{2}}\sqrt{{({A}_{1}\alpha +b\beta )}^{2}{\langle {\hat{X}}^{2}\rangle }^{2}+{({A}_{2}\beta +b\alpha )}^{2}{\langle {\hat{P}}^{2}\rangle }^{2}}\hfill \\ \hfill & \leqslant \sqrt{{d}_{1}^{2}+{d}_{2}^{2}}\sqrt{{({A}_{1}\alpha +b\beta )}^{2}+{({A}_{2}\beta +b\alpha )}^{2}}\sqrt{{\langle {\hat{X}}^{2}\rangle }^{2}+{\langle {\hat{P}}^{2}\rangle }^{2}}\hfill \\ \hfill & \leqslant \sqrt{{d}_{1}^{2}+{d}_{2}^{2}}\sqrt{{({A}_{1}\alpha +b\beta )}^{2}+{({A}_{2}\beta +b\alpha )}^{2}}(\langle {\hat{X}}^{2}\rangle +\langle {\hat{P}}^{2}\rangle )\hfill \\ \hfill & \leqslant \sqrt{{d}_{1}^{2}+{d}_{2}^{2}}\sqrt{({\alpha }^{2}+{\beta }^{2})}\sqrt{({A}_{1}^{2}+2{b}^{2}+{A}_{2}^{2})}(\langle {\hat{X}}^{2}\rangle +\langle {\hat{P}}^{2}\rangle )\hfill \\ \hfill & \leqslant \sqrt{{d}_{1}^{2}+{d}_{2}^{2}}\sqrt{({\alpha }^{2}+{\beta }^{2})}\sqrt{({\chi }^{2}+{\phi }^{2})}(\langle {\hat{X}}^{2}\rangle +\langle {\hat{P}}^{2}\rangle )\hfill \end{align*}$$

$$\begin{equation*}{d}_{1}^{2}\langle {\hat{X}}^{2}\rangle +{d}_{2}^{2}\langle {\hat{P}}^{2}\rangle \leqslant ({d}_{1}^{2}+{d}_{2}^{2})(\langle {\hat{X}}^{2}\rangle +\langle {\hat{P}}^{2}\rangle )\end{equation*}$$

$$\begin{align*}\hfill 8b\alpha \beta ({A}_{1}\langle {\hat{X}}^{2}\rangle +{A}_{2}\langle {\hat{P}}^{2}\rangle )& \leqslant 2({b}^{2}+{A}_{1}^{2}+{A}_{2}^{2})(\langle {\hat{X}}^{2}\rangle +\langle {\hat{P}}^{2}\rangle )({\alpha }^{2}+{\beta }^{2})\hfill \\ \hfill & \leqslant 2({\chi }^{2}+{\phi }^{2})(\langle {\hat{X}}^{2}\rangle +\langle {\hat{P}}^{2}\rangle )({\alpha }^{2}+{\beta }^{2}).\hfill \end{align*}$$

We have now successfully isolated the d factors, so that the bound only depends on ${d}_{1}^{2}+{d}_{2}^{2}={v}^{2}$. Putting everything together, the second term in (C19) can now be bounded by

$$\begin{align*}\hfill & (\langle {\hat{X}}^{2}\rangle +\langle {\hat{P}}^{2}\rangle )\left[2({\chi }^{2}+{\phi }^{2})({\alpha }^{2}+{\beta }^{2})+4\vert v\vert \sqrt{{\chi }^{2}+{\phi }^{2}}\sqrt{{\alpha }^{2}+{\beta }^{2}}+\vert v{\vert }^{2}\right]\hfill \\ \hfill & \quad \leqslant 2(\langle {\hat{X}}^{2}\rangle +\langle {\hat{P}}^{2}\rangle ){\left[\sqrt{({\chi }^{2}+{\phi }^{2})({\alpha }^{2}+{\beta }^{2})}+\vert v\vert \right]}^{2}\hfill \\ \hfill & \quad \leqslant 2(\langle {\hat{X}}^{2}\rangle +\langle {\hat{P}}^{2}\rangle )({\chi }^{2}+{\phi }^{2}+\vert v{\vert }^{2})({\alpha }^{2}+{\beta }^{2}+1)\hfill \\ \hfill & \quad \leqslant \frac{1}{2}({\chi }^{2}+{\phi }^{2}+\vert v{\vert }^{2}){(\langle {\hat{X}}^{2}\rangle +\langle {\hat{P}}^{2}\rangle +{\alpha }^{2}+{\beta }^{2}+1)}^{2}\hfill \\ \hfill & \quad =\frac{1}{2}({\chi }^{2}+{\phi }^{2}+\vert v{\vert }^{2}){(2N+2)}^{2}.\hfill \end{align*}$$

And finally, we can bound the variance of $\hat{O}$ as

$$\begin{equation}\langle {\hat{O}}^{2}(t)\rangle -{\langle \hat{O}(t)\rangle }^{2}\leqslant 2({\phi }^{2}+{\chi }^{2})\left[{(2N(t)+1)}^{2}-\frac{1}{2}\right]+2({\phi }^{2}+{\chi }^{2}+\vert v{\vert }^{2})[{(N(t)+1)}^{2}].\end{equation} \tag{ C20 }$$

The bounds (C18) and (C20) are most likely not tight, and do not explicitly give back the previous bound (25) in the limit α, β, v → 0. However, they show that, even in the best possible case, adding a non-zero displacement to the state should not allow to achieve better scaling that what we obtained with vacuum squeezed states. Most importantly, these bounds differ from (25) only by some prefactors. This indicates that much of the insight we have obtained about scaling regimes should apply also when we have displacement.

Let us study what happens when κT ≪ 1. In this regime, the quench is too short for the system to reach the steady-state. All the quantities will depend on T. We can develop the expressions above in terms of κT. For instance, we will have ${C}_{1}\sim 1-2{(\kappa T)}^{2}+\frac{{(2\kappa T)}^{3}}{3}+O({(\kappa T)}^{4})$, ${C}_{2}\sim 1-\frac{4}{3}{(\kappa T)}^{3}+\frac{{(2\kappa T)}^{4}}{8}+O({(\kappa T)}^{5})$, and ${C}_{3}\sim 1-\frac{2}{3}{(\kappa T)}^{4}+\frac{{(2\kappa T)}^{5}}{30}+O({(\kappa T)}^{6})$. After a tedious but straightforward development of the quadratures, we arrive to

$$\begin{align}\hfill \langle {\hat{x}}^{2}\rangle & =\frac{1}{2}+\frac{{(\omega T)}^{2}}{2}+\frac{2{\omega }^{2}\kappa {T}^{3}}{3}+(1-{g}^{2})\left(-\frac{{(\omega T)}^{2}}{2}-\frac{{(\omega T)}^{4}}{6}+\frac{2}{3}{\omega }^{2}\kappa {T}^{3}+\frac{4}{15}{\omega }^{4}\kappa {T}^{5}\right),\hfill \\ \hfill \langle {\hat{p}}^{2}\rangle & =\frac{1}{2}+(1-{g}^{2})\left(-\frac{{(\omega T)}^{2}}{2}+\frac{2}{3}{\omega }^{2}\kappa {T}^{3}\right),\hfill \\ \hfill \langle :\;\hat{x}\hat{p}\;:\rangle & =\frac{\omega T}{2}-\frac{\omega \kappa {T}^{2}}{2}+\frac{1-{g}^{2}}{2}\left(-\omega T-\frac{2}{3}{(\omega T)}^{3}+\omega \kappa T+{\omega }^{3}\kappa {T}^{4}\right).\hfill \end{align} \tag{ G3 }$$

Note that if we set κ = 0, we recover the expressions (F7).

Equipped with the expressions above, we can now extract the relevant parameters υ, θ and |z|. Keeping only the largest coefficients, we find

$$\begin{equation*}\mathrm{tan}(2\theta )=\frac{2\langle :\;\hat{x}\hat{p}\;:\rangle }{\langle {\hat{p}}^{2}\rangle -\langle {\hat{x}}^{2}\rangle }\sim -\frac{2}{\omega T}+(1-{g}^{2})\frac{2\omega T}{3},\end{equation*}$$

$$\begin{equation*}{\upsilon }^{2}=4(\langle {\hat{x}}^{2}\rangle \langle {\hat{p}}^{2}\rangle -{\langle :\;\hat{x}\hat{p}\;:\rangle }^{2})\sim 1+\frac{1}{0}3{\omega }^{2}\kappa {T}^{3}-(1-{g}^{2})\left(\frac{7}{10}{\omega }^{4}\kappa {T}^{5}\right),\end{equation*}$$

and

$$\begin{equation*}\mathrm{Tr}[\sigma ]=\upsilon \,\mathrm{cosh}(2\vert z\vert )=\frac{{(\omega T)}^{2}}{2}-(1-{g}^{2})\frac{{(\omega T)}^{4}}{6}.\end{equation*}$$

Now, we must distinguish two important cases. We can define the time-scale T₀ = ω^−2/3 κ^−1/3. As long as T ≪ T₀, we have υ² ∼ 1, and the state is almost pure. In this regime, we can work out the dominant terms for all observables, and their derivative. Up to prefactors, the results are cosh(2|z|) ∼ sinh(2|z|) ∼ (ωT)², υ² − 1 ∼ ω² κT³, $\omega \,\frac{\mathrm{d}\vert z\vert }{\mathrm{d}\omega }\sim {(\omega T)}^{2}$, $\omega \,\frac{\mathrm{d}\theta }{\mathrm{d}\omega }\sim \omega T$, $\omega \,\frac{\mathrm{d}\upsilon }{\mathrm{d}\omega }\sim {\omega }^{4}\kappa {T}^{5}$. The various terms in (G1) scale therefore as ${\mathrm{sinh}}^{2}(2\vert z\vert ){\left(\frac{\mathrm{d}\theta }{\mathrm{d}\omega }\right)}^{2}\propto {(\omega T)}^{6}$, ${\left(\frac{\mathrm{d}\vert z\vert }{\mathrm{d}\omega }\right)}^{2}\propto {(\omega T)}^{4}$, and $\frac{1}{{\upsilon }^{2}-1}{\left(\frac{\mathrm{d}\upsilon }{\mathrm{d}\omega }\right)}^{2}\propto {\omega }^{6}\kappa {T}^{7}$. Recall that the term $\frac{{\upsilon }^{2}}{1+{\upsilon }^{2}}$ is always bounded between 1/2 and 1 and can be ignored. Since we are in the regime $\frac{1}{\omega }\ll T\ll \frac{1}{\kappa }$, the dominant term is the first one, that is, we recover the T⁶ scaling we obtained in the absence of decoherence. Hence, the short-time dynamics for the noisy quench is essentially identical to the quench without boson loss, and gives the same QFI scaling. Note that the QFI arises from the interplay between dθ/dω and sinh(2|z|). During the evolution, the system experiences both phase-shift and amplification in a correlated manner. Hence, to saturate the QFI a measurement sensitive to this correlation is required. By contrast, a measurement of a single quadrature will provide a T⁴ (cf section 7.2). This means that the fluctuations need to be quadrature dependent, although not necessarily below shot noise.

In the case T ≫ T₀, by contrast, υ becomes large, and the thermal noise cannot be neglected anymore. In this case, we obtain modified scalings for the observables $\mathrm{cosh}(2\vert z\vert )\sim \mathrm{sinh}(2\vert z\vert )\sim \frac{\omega }{\sqrt{\kappa }}{T}^{1/2}$, υ² − 1 ∼ ω² κT³, $\omega \,\frac{\mathrm{d}\vert z\vert }{\mathrm{d}\omega }\sim {(\omega T)}^{2}$, $\omega \,\frac{\mathrm{d}\theta }{\omega }\sim \omega T$, $\omega \,\frac{\mathrm{d}\upsilon }{\omega }\sim {(\omega T)}^{2}$. The terms in the QFI scale now as ${\mathrm{sinh}}^{2}(2\vert z\vert ){\left(\frac{\mathrm{d}\theta }{\mathrm{d}\omega }\right)}^{2}\propto \frac{{\omega }^{4}}{\kappa }{T}^{3}=\frac{{\omega }^{2}}{{\kappa }^{2}}\frac{{T}^{3}}{{T}_{0}^{3}}$, ${\left(\frac{\mathrm{d}\vert z\vert }{\mathrm{d}\omega }\right)}^{2}\propto {(\omega T)}^{4}$, and $\frac{1}{{\upsilon }^{2}-1}{\left(\frac{\mathrm{d}\upsilon }{\mathrm{d}\omega }\right)}^{2}\propto \frac{{\omega }^{2}}{\kappa }T$. The dominant term is again the first one (note that the second term gives a larger scaling, but is smaller in absolute value, because $T\ll \frac{1}{\kappa }$). Hence, we see that we now achieve a cubic scaling with T.

In the regime κT ≫ 1, the quench is long enough for the system to reach its steady-state. All observables, as well as the QFI, saturate at a constant value. Working out the expressions (G2) leads to

$$\begin{align*}\hfill \langle {\hat{x}}^{2}\rangle & =\frac{1}{2}+\frac{{\omega }^{2}}{4{\kappa }^{2}}-(1-{g}^{2})\left(\frac{{\omega }^{2}}{4{\kappa }^{2}}+\frac{{\omega }^{4}}{4{\kappa }^{4}}\right),\hfill \\ \hfill \langle {\hat{p}}^{2}\rangle & =\frac{1}{2}-(1-{g}^{2})\left(\frac{{\omega }^{2}}{4{\kappa }^{2}}\right),\hfill \\ \hfill \langle :\;\hat{x}\hat{p}\;:\rangle & =\frac{\omega }{4\kappa }-(1-{g}^{2})\left(\frac{\omega }{2\kappa }\right),\hfill \end{align*}$$

from which we can obtain the dominant terms for the observables and their derivatives, $\mathrm{cosh}(2\vert z\vert )\sim \mathrm{sinh}(2\vert z\vert )\sim \frac{\omega }{\kappa }$, $\mathrm{tan}(2\theta )\sim \frac{\kappa }{\omega }$, ${\upsilon }^{2}-1\sim \frac{{\omega }^{2}}{{\kappa }^{2}}$, $\omega \,\frac{\mathrm{d}\vert z\vert }{\mathrm{d}\omega }\sim \frac{{\omega }^{2}}{{\kappa }^{2}}$, $\omega \,\frac{\mathrm{d}\theta }{\mathrm{d}\omega }\sim \frac{\omega }{\kappa }$, $\omega \,\frac{\mathrm{d}\upsilon }{\mathrm{d}\omega }\sim \frac{{\omega }^{3}}{{T}^{3}}$. Plugging everything in the QFI, we find that ${Q}_{\omega }\sim \frac{{\omega }^{4}}{{\kappa }^{4}}$.

To summarize, the various regimes for the quench in the presence of boson loss are the following: for $T\ll \frac{1}{\omega }$, the dynamics follow non-universal behavior. For $\frac{1}{\omega }\ll T\ll \frac{1}{{\omega }^{2/3}{\kappa }^{1/3}}$, the thermal noise is still negligible; the system behaves as in the absence of dissipation, and the QFI scales like T⁶. For $\frac{1}{{\omega }^{2/3}{\kappa }^{1/3}}\ll T\ll \frac{1}{\kappa }$, the thermal noise becomes important; the system enters a thermal squeezed state, with both squeezing and thermal noise increasing with T. In this regime, the QFI achieves a modified scaling T³. Finally, for $T\gg \frac{1}{\kappa }$, the system reaches its steady-state, and the QFI saturates at a value $\frac{{\omega }^{4}}{{\kappa }^{4}}$. All of these findings are corroborated by the simulations presented in the main text.