Collective quantum enhancement in critical quantum sensing

Critical quantum sensing (CQS) is by now an established approach, based on the exploitation of quantum properties spontaneously developed in proximity of phase transitions. In the context of quantum metrology, the performance of parameter estimation protocols¹ is assessed based on the precision that can be achieved using a limited amount of fundamental resources, such as the size of the probe system and/or the protocol duration time. By using non-classical properties such as entanglement and superposition, it is possible to overcome the performance of any sensing strategy based on purely classical resources². Numerous theoretical studies have shown that a quantum-enhanced sensing precision can be achieved exploiting static^{3,4,5,6,7,8,9,10,11,12,13} or dynamical^14,15,16,17 properties of many-body systems in proximity of the critical point. Despite critical quantum systems becoming divergingly slow as they approach the phase transition, the optimal precision scaling can, in principle, be saturated with respect to both probe size and time¹⁸. Furthermore, even if prior information on the parameter to be estimated is required to operate the sensor close to the critical point, CQS protocols work well also in a global-sensing scenario when adaptive strategies are applied^19,20,21. First implementations with Rydberg atoms²², nuclear magnetic resonance²³ and superconducting quantum circuits^24,25 demonstrate the experimental feasibility of the CQS approach.

Recently, it has been shown²⁶ that CQS protocols can also be conceived using finite-component phase transitions (FCPTs), where the thermodynamic limit is replaced with a rescaling of the system parameters^27,28. This class of phase transitions can emerge in quantum resonators with atomic^{27,29,30,31,32,33} or Kerr-like^28,34,35,36 nonlinearities, and it is of high relevance for CQS for two main reasons: (i) It provides a tractable theoretical framework to analyze the fundamental precision bounds and to design optimal parameter estimation protocols. For instance, FCPTs have been used to demonstrate the intrinsic constant-factor advantage of dynamical over statical CQS protocols^37,38 and the presence of apparent super-Heisenberg scalings when focusing on a single resource^39,40,41. Continuous-measurement schemes^42,43 have been designed to efficiently retrieve information in the dissipative case, where the optimal precision bounds are in general not achievable. In the presence of a thermal bath, CQS outperforms passive strategies when preparation and measurement times are not negligible⁴⁴. (ii) As demonstrated by recent experiments, FCPTs can be controllably implemented using atomic⁴⁵, polaritonic^46,47 and circuit-QED ^{48,49,50,51,52,53} platforms. Accordingly, CQS protocols can be implemented with small-scale devices, without requiring the implementation and control of many-body quantum systems. CQS protocols based on FCPTs have been designed for quantum resonators^{24,54,55,56,57,58,59}, single trapped-ions⁶⁰, optomechanical^61,62 or magnomechanical⁶³ devices, spin impurities^64,65 and Rabi-like systems^66,67,68,69. Finally, physical devices undergoing FCPTs can be interconnected to form one- or two-dimensional arrays of critical systems, thereby extending CQS approaches to complex extended systems.

Here, we theoretically show that collective quantum enhancement in CQS can be achieved using an array of coupled critical quantum sensors. In particular, we consider a model composed of a chain of nonlinear quantum resonators coupled via a parametric nearest-neighbor interaction. We focus on the weak-nonlinearity limit, where each resonator of the array can become locally critical and undergo an FCPT. We assess the metrological performance of a static CQS protocol by analytically evaluating the quantum Fisher information (QFI) over the ground state manifold, under Gaussian approximation. In proximity to the critical point, which is the optimal working regime, we find closed-form analytical solutions for the asymptotic scaling of the QFI with respect to fundamental resources, i.e., the number of resonators, the number of excitations and the protocol-duration time. This analysis allows us to benchmark the estimation precision achievable with the coupled-resonator chain against that of an equal number of independent critical sensors. We find that the coupled chain yields a quadratic enhancement of the QFI scaling with respect to the number of resonators, compared to the linear scaling observed in an uncoupled array. We demonstrate that this conceptual result still holds when including physically relevant corrections to the minimal model. We first develop a perturbation theory to estimate the saturation point due to local Kerr nonlinearities, and find that the resonator chain presents a practically relevant advantage for finite values of physical parameters. We then show that, using a dynamical CQS protocol, it is still possible to achieve collective quantum advantage in an open-quantum-system setting. We unveil a non-trivial transition between a short-time unitary regime and the stationary regime achieved for long evolution times.

Summary of results

Critical quantum sensing: adiabatic protocol

Here, we briefly introduce the CQS protocol and the methods used to assess parameter-estimation performance. We consider as a sensor a system whose Hamiltonian depends on an unknown parameter ω to be estimated, while all other parameters are known, and such that one control parameter ϵ can be tuned in real-time. A static CQS protocol can be defined as a three-step algorithm: (i) the system is initialized in its ground state in a trivial region of the parameter space, say ϵ = 0; (ii) the control parameter ϵ is adiabatically tuned close to the critical value ϵ_c, slowly enough so that the system never leaves its instantaneous ground state; (iii) a measurement is performed on the final state, and the outcome is used to infer information about the value of ω. The working principle behind adiabatic CQS protocols is the exploitation of the high susceptibility of the system ground state with respect to small variations of the parameter to be estimated.

The ultimate precision achievable in parameter estimation is given by the quantum Cramér-Rao bound². The mean squared error is lower-bounded by (Delta {omega }^{2}ge {left(M{{{mathcal{I}}}}right)}^{-1}), where M is the number of independent measurements and ({{{mathcal{I}}}}) is the QFI. For our adiabatic protocol, the QFI must be evaluated [as in Eq. (8)] on the manifold defined by the system ground state as the parameter ω is varied. The QFI represents the maximum amount of information that can be extracted on average for a given repetition of the protocol, assuming that the optimal observable is measured after the adiabatic sweep. The performance of a parameter-estimation protocol can be assessed by analyzing the scaling of the QFI with the fundamental resources used during the protocol, e.g., the size of the probe system N. For example, in Ramsey interferometry, the QFI for phase estimation can grow at most linearly as ({{{mathcal{I}}}} sim N) when the spins are in a separable state, while it can achieve the quadratic scaling ({{{mathcal{I}}}} sim {N}^{2}) when the spins are entangled⁷⁰. This quadratic enhancement is dubbed Heisenberg scaling and is achievable only using nonclassical resources such as entanglement or squeezing.

The size of the probe system is not the only relevant resource. When dealing with critical systems the protocol time must unavoidably be taken into account as a fundamental resource. When the critical point is approached, the energy gap between the ground and first-excited state tends to zero, and so the time T required to perform an adiabatic sweep unavoidably diverges. Analyzing the scaling of the QFI only with respect to the number of probe systems can then lead to an apparent super-Heisenberg scaling. Once time is taken into account^18,26, under very general assumptions, the optimal Heisenberg scaling is given by ({{{mathcal{I}}}} sim {T}^{2}{N}^{2}). However, this is strictly valid only for closed quantum systems. When decoherence and dissipation are included, stricter asymptotic bounds apply^71,72,73. Here, we limit ourselves to the unitary case, as our focus is on the collective enhancement achievable with quantum-correlated probes, and our results can be generalized to the dissipative case.

Critical parametric quantum sensor

Let us now introduce our physical model of a critical quantum sensor, which is a quantum nonlinear resonator with two-photon pumping terms, described by the Hamiltonian

Here and throughout the manuscript, we set ℏ = 1. This system undergoes an FCPT in the limit χ → 0, when ϵ approaches the critical point ϵ_c = ω. This model is of direct experimental relevance⁵⁰, and it effectively describes the critical scaling of fully-connected models³⁸, a broad class of quantum-optical critical systems. The CQS protocol consists of initializing the system in the vacuum state, and adiabatically tuning ϵ to a value close to ϵ_c. To work in this critical regime, we assume to have good prior information on the parameter to be estimated, namely ω = ω₀ + δω, where ω₀ is known and δω ≪ ω₀ is an unknown frequency shift to be estimated. For global quantum sensing, i.e., when prior information is not available, efficient adaptive strategies can be implemented^19,20.

Under Gaussian approximation (χ = 0), we can obtain analytical results. Indeed, we first compute the ground state manifold, and then evaluate the QFI for the estimation of ω. Consequently, we provide an upper-bound on the estimation precision with respect to the involved resources, i.e., the number of photons ({N}_{{{{rm{loc}}}}}=langle {a}_{1}^{{{dagger}} }{a}_{1}rangle) and the protocol duration time T. Here, the subscript stands for local and will be meaningful in the multipartite case. By applying time-dependent perturbation theory, it can be shown²⁶ that the time required to perform the adiabatic sweep is of the order of T ~ (ηλ)⁻¹, where λ is the energy gap between the ground and first-excited state, and η < 1 is a parameter which controls the adiabaticity of the protocol. As shown in Fig. 1, when approaching the critical point the QFI and the photon number diverge, while the energy gap between the ground and the first-excited state vanishes, as expected. For a single-mode critical sensor, the QFI scales as ({{{{mathcal{I}}}}}_{{{{rm{sm}}}}} sim 2{eta }^{2}{T}^{2}{N}_{{{{rm{loc}}}}}^{2}). Here, the subscript “sm” stands for “single-mode”. We stress that in the limit χ → 0, this result can be analytically derived. Notice that in any physical implementation nonlinear terms are always present, and theoretically they are required to obtain a stable dynamics in the whole parameter space³⁴. In the case of finite nonlinearity, this scaling is valid until the photon number becomes large enough to break the Gaussian approximation and reach saturation. As shown in “Methods”, perturbation theory reveals that the saturation point occurs when ({N}_{{{{rm{loc}}}}} sim {root 3 of {omega /132chi }}).

Collective quantum advantage in CQS

In this work, we perform a comparison between an ensemble of M independent critical sensors, with a chain of parametrically-coupled resonators (see Fig. 1). The QFI is additive when using M independent sensors or, equivalently, M uncorrelated repetitions of the parameter-estimation protocol. Accordingly, when M independent critical sensors are used, we straightforwardly obtain ({{{{mathcal{I}}}}}_{{{{rm{ind}}}}} sim 2M{eta }^{2}{T}^{2}{N}_{{{{rm{loc}}}}}^{2}), with the subscript “ind” standing for “independent”. This result can also be rewritten in terms of the total number of photons N = MN_loc, as ({{{{mathcal{I}}}}}_{{{{rm{ind}}}}} sim 2{eta }^{2}{T}^{2}{N}^{2}/M).

Our collective critical sensor is defined as a coupled-resonator chain with Hamiltonian

We derive analytical solutions for the QFI for the estimation of ω, as well as for the photon number and the energy gap, under Gaussian approximation. The result for the QFI is shown in Eqs. (21) and (22) for even and odd M, respectively. We then obtain approximate closed-form analytical solutions in two relevant limits: in proximity of the phase transition and in the continuous limit (large M). The analysis unveils interesting features of this many-body model composed of locally critical constituents. In this summary, we focus on the most effective regime for quantum sensing, which is the critical scaling of a chain with an odd number of resonators. The QFI achievable in this case scales as ({I}_{O}=2{eta }^{2}{M}^{2}{T}^{2}{N}_{{{{rm{loc}}}}}^{2}=2{eta }^{2}{T}^{2}{N}^{2}). Thus, the critical chain presents a quadratic enhancement in the number of resonators M with respect to the non-interacting case. The critical chain is advantageous even when considering a finite value of the Kerr nonlinearity. Indeed, it can be effectively described by a single-mode resonator with a diluted nonlinearity χ/M, and the saturation point is reached when ({N} sim {root 3 of {{M}omega /132chi }}). Finally, the collective quantum advantage can also be achieved in a dissipative setting, when a dynamical protocol is used. Different regimes can be identified with respect to the timescale set by 3/(4MN_locΓ), where Γ is the dissipation rate. For a very short evolution time, the dynamics is effectively unitary, and a scaling advantage is still achievable. Conversely, for a very long evolution time the collective quantum advantage fades away, in accordance with fundamental bounds^71,72,73. We also find a non-trivial transient regime characterized by a collective constant-factor advantage.

Local models

In this section, we analyze the performance of single-and two-mode critical quantum sensors. In both cases, we compute the QFI and express it in terms of the resources involved, i.e., the average number of photons in each mode and the duration time of the sensing protocol. The results will be used as a benchmark, as well as an intermediate step to solve the interacting chain.

Single-mode critical quantum system

A single-mode critical quantum sensor can be modeled, under Gaussian approximation, by the Hamiltonian

where ϵ is the squeezing parameter, ω is the mode frequency and a₁ is the bosonic mode. The full model presented in Eq. (1) undergoes a second-order phase transition in proximity of the critical point ϵ_c = ω in the limit χ → 0. For ϵ < ϵ_c the model is well approximated by the Hamiltonian (3), which is diagonalizable via a Bogoliubov Transformation (BT), obtaining

We defined (lambda =sqrt{{omega }^{2}-{epsilon }^{2}}), and d₁ is the normal mode related to a₁ via the transformation

where the parameters t and s are

Let us focus now on the estimation of the parameter ω. Given a manifold of pure states (| psi rangle), the QFI can be evaluated as²

where (| {partial }_{omega }psi rangle) is the partial derivative with respect to the parameter to be estimated. To assess the performance of the adiabatic CQS protocol, we will compute the QFI associated with the system ground state (| grangle). The ground state of the Hamiltonian (3) is the single-mode squeezed vacuum (| grangle ={S}_{1}(xi )| 0rangle). We define the single-mode squeezing operator ({S}_{i}(xi )={e}^{frac{1}{2}left({xi }^{* }{a}_{i}^{2}-xi {a}_{i}^{{{dagger}} 2}right)}) and the squeezing parameter (xi =| xi | =ln left(s+sqrt{{s}^{2}+1}right)). This state has an average number of photons ({N}_{{{{rm{loc}}}}}={sinh }^{2}(| xi | )={s}^{2}), where “loc” stands for “local”. By inserting in (8) the expression of the system ground state, we obtain the single-mode QFI for the estimation of ω:

The dominant contribution to the duration of the sensing protocol is given by the time T required to perform the adiabatic sweep, which must be larger²⁶ than the inverse of the gap between the ground and first-excited state. From (4), it immediately follows that the energy gap is (lambda =sqrt{{omega }^{2}-{epsilon }^{2}}). Consequently, the protocol duration time can be expressed as T ≈ (ηλ)⁻¹, with η < 1. To analyze the critical scaling approached when ϵ → ϵ_c, we define ϵ = (1 − x)ω with 0 < x < 1, and compute the QFI in the asymptotic regime x → 0. Close to the critical point, we have ({{{{mathcal{I}}}}}_{{{{rm{sm}}}}}approx {(2sqrt{2}omega x)}^{-2}), (Tapprox {(omega sqrt{2x})}^{-1}) and ({N}_{{{{rm{loc}}}}}approx {(8x)}^{-frac{1}{2}}). Using the latter, we express x in terms of the number of photons as (xapprox {(2sqrt{2}{N}_{{{{rm{loc}}}}})}^{-2}), which implies T ≈ 2N_loc/ηω. Finally, in terms of the resources involved the single-mode QFI scales as

Two-mode critical quantum system

Similarly to the single-mode case, a two-mode critical quantum system can be modeled with a two-mode squeezing Hamiltonian. The single-resonator critical sensor can be generalized to a two-mode critical system as follows:

where ϵ is the parametric-coupling strength, ω is the frequency of each mode and a₁ and a₂ are two independent bosonic modes. Here, the subscript “tm” means “two-mode”. Also in this case, the model undergoes a second-order phase transition, when a vanishingly small Kerr nonlinearity is included²⁸. The normalization of the coupling parameter is chosen to have the same critical point ϵ_c = ω of the Hamiltonian (1). When ϵ < ϵ_c, the Hamiltonian (11) can be diagonalized with the following BT:

where d₁,  d₂ are the normal modes. We can then write

The ground state of the system is a two-mode squeezed vacuum state (| grangle ={S}_{1,2}(xi )| 0rangle), with ({S}_{i,j}(xi )={e}^{left({xi }^{* }{a}_{i}{a}_{j}-xi {a}_{i}^{{{dagger}} }{a}_{j}^{{{dagger}} }right)}) a two-mode squeezing operator. Notice that the squeezing parameter ξ = ∣ξ∣ is the same as in the single-mode case. The average number of photons in the ground state is (N=2{sinh }^{2}(| xi | )=2{s}^{2}) since now there are two normal modes with the same eigenvalue λ and the same number of photons ({N}_{{{{rm{loc}}}}}=langle {a}_{1}^{{{dagger}} }{a}_{1}rangle =langle {a}_{2}^{{{dagger}} }{a}_{2}rangle ={s}^{2}). Using (8), we can compute the two-mode QFI for the estimation of ω:

In the asymptotic limit x → 0, we have ({{{{mathcal{I}}}}}_{{{{rm{tm}}}}}approx {(2omega x)}^{-2}) and ({N}_{{{{rm{loc}}}}}approx {(8x)}^{-frac{1}{2}}). Since the energy gap is still λ, the protocol duration time remains unchanged with respect to the single mode case. Indeed, expressing x as (xapprox {(2sqrt{2}{N}_{{{{rm{loc}}}}})}^{-2}) leads to T ≈ 2N_loc/ηω. Finally, in terms of the involved resources the two-mode QFI scales as

Comparing this result with Eq. (10), we see that ({{{{mathcal{I}}}}}_{{{{rm{tm}}}}}) is simply twice what is obtained with a single critical sensor. A two-resonator chain is then equivalent to two independent resonators, suggesting that there is no advantage in considering coupled critical sensors. However, in the following, we demonstrate that this conclusion does not hold for larger chains.

Critical resonator chain

Let us finally consider the critical resonator chain. Under Gaussian approximation, it is described by the Hamiltonian

Assuming periodic boundary conditions a_M+j = a_j, the Hamiltonian (16) is diagonalizable in the reciprocal space. By applying the discrete Fourier transform ({a}_{j}={left(Mright)}^{-frac{1}{2}}sum limits_{k}{a}_{k}{e}^{-ikj}), we obtain

Here, δ_M,2m is a Kronecker delta which is nonzero only for even M. The terms H₀,  H_π are single-resonator critical models defined as in Eq. (3), but replacing a₁ with a₀ and a_π, respectively. The sign of the coupling for the a_π mode must be inverted, a phase difference that has no consequence for the current analysis. The periodic boundary conditions imply that we can restrict ourselves to the first Brillouin zone (FBZ), and k = 2nπ/M, with (nin {mathbb{Z}}). The FBZ is defined as follows: if M is even, n ∈ [ − M/2 + 1,  M/2]; if M is odd, n ∈ [ − (M − 1)/2,  (M − 1)/2]. We point out that the mode a_π is present only in the case in which the number of modes M is even.

In reciprocal space, the multi-mode Hamiltonian (17) consists of M independent critical systems. Specifically, two single-mode systems defined by a₀ and a_π, along with a series of two-mode systems defined by pairs {a_k,  a_−k}. Each element can be independently diagonalized using the single- and two-mode BTs derived previously, by replacing (epsilon to epsilon cos (k)). Notice that the BT can be performed only if ϵ < ϵ_c. In the multi-mode case, we have many critical points defined by the dispersion relation ({epsilon }_{{{{rm{c}}}}}=omega /cos (k)). This means that the multi-mode Hamiltonian is diagonalizable via a BT only if the parameter ϵ is smaller than the lowest among the critical points. Therefore, we restrict ourselves to the case ϵ < ω. If this condition is satisfied, applying the BT to (17) results in

where ({lambda }_{k}=sqrt{{omega }^{2}-{epsilon }^{2}{cos }^{2}(k)}), and d_k are independent normal modes. The modes d₀,  d_π are defined as in (5), while the modes d_±k as in (12), with the substitution (epsilon to epsilon cos (k)). The ground state of the system is then a multi-mode squeezed vacuum state

We recall that ({S}_{i}({xi }_{i})={e}^{frac{1}{2}left({xi }_{i}^{* }{a}_{i}^{2}-{xi }_{i}{a}_{i}^{{{dagger}} 2}right)}) with i = {0, π} is a single-mode squeezing operator associated with the modes a₀,  a_π. Conversely, ({S}_{k,-k}({xi }_{k})={e}^{left({xi }_{k}^{* }{a}_{k}{a}_{-k}-{xi }_{k}{a}_{k}^{{{dagger}} }{a}_{-k}^{{{dagger}} }right)}) is a two-mode squeezing operator related to pairs of coupled modes {a_k,  a_−k}. The squeezing parameters are ({xi }_{k}=| {xi }_{k}| =ln left({s}_{k}+sqrt{{s}_{k}^{2}+1}right)) for k ∈ [0, π/2) and ξ_k = ∣ξ_k∣e^iπ for k ∈ [π/2, π]. The parameter s_k is obtained from (7) with the substitution (epsilon to epsilon cos (k)). The total average number of photons in the multi-mode ground state is given simply by the sum of the photons in each mode a_k. Given that the ground state is a multi-mode squeezed vacuum state, each mode contains ({N}_{k}={sinh }^{2}(| {xi }_{k}| )={s}_{k}^{2}) photons. Since ({s}_{k}^{2}={s}_{-k}^{2}), the total number of photons in the ground state can be expressed as

Moreover, the local number of photons at each site j is determined by ({N}_{{{{rm{loc}}}}}=langle {a}_{j}^{{{dagger}} }{a}_{j}rangle ={sum }_{k}langle {a}_{k}^{{{dagger}} }{a}_{k}rangle /M=N/M). A detailed derivation is provided in “Methods” where we show that N_loc is independent of the specific site, i.e., each resonator has the same number of photons. It is important to stress that in the reciprocal space, when the critical point is approached, the photons are strongly confined only in the mode a₀ when M is odd and in the modes a₀,  a_π when M is even. The fraction of photons populating all the other modes vanishes as ϵ → ϵ_c (see Fig. 2). In this scenario, we can effectively approximate the critical chain as a single-mode system with N = MN_loc photons if M is odd, or as a two-mode system with N = MN_loc/2 photons per mode when M is even.

Having obtained its ground state, we can now use (8) to calculate the QFI of the critical chain. From Eq. (19) we see that the ground state is separable in reciprocal space, hence the QFI of the critical chain is given by the sum of the QFIs of each individual mode a₀ and a_π, and of each pair {a_k, a_−k}. We obtain

where ({{{{mathcal{I}}}}}_{{{{rm{E}}}}}) is the QFI for an even number of modes, while ({{{{mathcal{I}}}}}_{{{{rm{O}}}}}) is the QFI for an odd number of modes. Let us now provide closed-form analytical expressions of the multi-mode QFI in relevant asymptotic regimes. First, we consider a finite number of modes and compute explicitly the QFI close to the critical point. Then, we compute the QFI in the continuous limit M → ∞. In both cases, we express the QFI in terms of the involved resources, i.e., the local number of photons N_loc and the total protocol duration time T.

Asymptotic regimes: Critical scaling

In this section, we provide the scaling behavior of the critical chain QFI in the asymptotic regime where ϵ → ϵ_c, assuming a finite M. To analyze the critical scaling we define ϵ = (1 − x)ω with 0 < x < 1, where the critical point is approached for x → 0. We then expand Eqs. (21) and (22) in power series, up to the first order in x. Following this expansion, the QFI separates into two terms: the first depends solely on x and diverges as the critical point is approached, while the second term is a discrete sum dependent only on the number of modes M. This sum can be analytically computed using established results from trigonometric series, allowing us to derive a closed expression for the asymptotic QFI. See “Methods” for further mathematical details. Thus, in the limit x → 0, up to the terms of order x, the QFI for even and odd critical chains can be expressed as

A similar analysis applies to the average number of photons in the multi-mode ground state. Notice that the expression for the number of photons is given by (20), and the only distinction between odd and even M lies in the range of k values in the discrete sum. Then, expanding (20) up to terms of order (sqrt{x}) yields

where γ is the Euler-Mascheroni constant. See “Methods” for further details on the derivation of Eqs. (25–26). Our aim now is to express the multi-mode QFI in terms of the resources involved, namely the local number of photons N_loc and the protocol duration time T. To achieve this, we first express the parameter x as a function of the local number of photons using (25)–(26). Specifically, we find (xapprox {(sqrt{2}M{N}_{{{{rm{loc}}}}})}^{-2}) when M is even and (xapprox {(2sqrt{2}M{N}_{{{{rm{loc}}}}})}^{-2}) when M is odd. It is important to note that these relations hold true under the condition that the local number of photons is sufficiently large, satisfying (pi {N}_{{{{rm{loc}}}}},gg ,leftvert ln left(M/pi right)+gamma -pi /2rightvert) for M even, and (pi {N}_{{{{rm{loc}}}}},gg ,leftvert ln left(2M/pi right)+gamma -pi /2rightvert) for M odd. This assumption is reasonable in the considered asymptotic regime, where M is finite and near the criticality. On the other hand, from (18) we have that the lowest energy gap is λ₀, which means that the protocol duration time is estimated as (Tapprox {(eta {lambda }_{0})}^{-1}approx {(eta omega sqrt{2x})}^{-1}). In terms of the local number of photons we write T ≈ MN_loc/ηω for M even and T ≈ 2MN_loc/ηω for M odd. Finally, inserting the expression of x in (23–24) and taking into account the protocol duration time, the QFI scales as

These solutions then predict a quadratic scaling of the QFI with the number of resonators M, which implies a collective quantum enhancement of the estimation precision. The difference between the even and odd cases arises from the fact that, for an odd M, only the mode with k = 0 becomes critical whereas, for an even M, both the modes with k = 0 and k = π become critical (see Fig. 2). In Fig. 3(a) and (b), we plot the optimized ratio ({{{mathcal{I}}}}/{T}^{2}) as a function of M and as a function of N_loc, respectively. We show that the asymptotic expressions (27)–(28) align perfectly with the full analytical solutions of Eqs. (21)–(22), in the asymptotic regime where x → 0 and M is finite.

Asymptotic regimes: Continuous limit

As M increases, the separation between adjacent modes decreases. This is due to the periodic boundary conditions, which impose Δk = 2π/M, so that Δk → 0 when M → ∞. Consequently, in this regime we have a continuum of modes a_k. Therefore, we can compute the QFI in (21)–(22) by approximating the discrete sums with definite integrals over this continuum, obtaining

Notice that, of course, by making this approximation, a numerical error is introduced. This error becomes negligible when the number of modes is sufficiently large, such that Δk → dk. We estimate the absolute error incurred when approximating a discrete sum with an integral via

By imposing (MDelta {k}^{2}max {{partial }_{k}f(k)}/2ll 1), we can evaluate how large M should be so that (29) is a good approximation of the QFI. In the limit of x → 0, we find that (max {{partial }_{k} , f(k)}approx {x}^{-frac{5}{2}}/5), which implies that, to have a negligible absolute error, we require (M,gg ,2{pi }^{2}{x}^{-frac{5}{2}}/5). This indicates that the closer we are to the critical point, the larger the M needed to ensure the validity of (29).

The same reasoning can be applied to the average number of photons in the ground state (20). However, the resulting integral cannot be computed analytically as we did with the QFI. An approximated analytical solution in the limit x → 0 is

Notice that in this regime, the local number of photons N_loc = N/M is independent of the number of modes but depends only on the proximity to the critical point. In this case, the parameter x can be expressed in terms of N_loc as (xapprox 8{e}^{-(1+2{N}_{{{{rm{loc}}}}})pi }). Furthermore, the protocol duration time is approximately (Tapprox {(eta {lambda }_{0})}^{-1}approx {(eta omega sqrt{2x})}^{-1}), which in terms of N_loc becomes (Tapprox {e}^{left(frac{1}{2}+{N}_{{{{rm{loc}}}}}right)pi }/4eta omega). Using these results, in terms of the involved resources we express the QFI as

By looking at (32), one might conclude that the optimal operating regime of the critical chain is the continuous limit, given that the QFI scales exponentially with the photon number resource. However, the validity of this expression is related to (30). As mentioned above, approximating the discrete sums defining the multi-mode QFI with definite integrals is a good approximation only if (M,gg ,2{pi }^{2}{x}^{-frac{5}{2}}/5). Considering the expression of x in terms of N_loc, we find that M must satisfy (M,gg ,C{e}^{5pi {N}_{{{{rm{loc}}}}}}), with (C={pi }^{2}{e}^{frac{5pi }{2}}/320sqrt{2}). This reveals that even for a low number of photons, the required M to ensure a negligible error is extremely large.

In Fig. 3(c), we plot the ratio ({log }_{10}left({{{mathcal{I}}}}/{T}^{2}right)) as a function of N_loc, and compare the exact optimized ratio obtained numerically from (22) with the quadratic scaling of (28) and the exponential scaling of (32). As expected, the QFI goes exponentially with N_loc only when ({N}_{{{{rm{loc}}}}}={{{mathcal{O}}}}(1)). Indeed, as soon as N_loc becomes comparable with M, the QFI is well approximated by the quadratic regime of (28). The exponential in Eq. (32) is then only an apparent super-Heisenberg scaling.

Performance analyses

We can now compare the performance of the critical resonator chain with that of an ensemble of M independent single-mode critical resonators. In the following, we first thoroughly analyze the collective precision enhancement of the array of critical resonators in the relevant asymptotic regimes. Then, we assess the validity of the Gaussian approximation of the critical models and provide an estimate of the saturation of the QFI due to finite Kerr nonlinearity. Finally, we discuss the effects of coupling with a Markovian bath, identifying the regimes in which the collective advantage can be achieved in the dissipative setting.

Asymptotic scaling

We first examine the case of an ensemble of M uncoupled critical resonators. This scenario corresponds to applying the CQS protocol with a single-mode system M times, independently. For independent measurements, the QFI is additive, meaning that overall ({{{{mathcal{I}}}}}_{{{{rm{ind}}}}}=M{{{{mathcal{I}}}}}_{{{{rm{sm}}}}}). For an array of independent critical sensors, from (10) we obtain

Conversely, for a critical chain of M interacting resonators, the optimal QFI scales as

for an even and odd M, respectively. The critical resonator chain presents a collective enhancement in sensing precision, resulting in a quadratic scaling of the QFI with respect to M. In Fig. 4, we show the growth of the QFI with respect to M for fixed values of N_loc, comparing even and odd critical chains with an equivalent ensemble of uncoupled resonators. The plot is numerically computed from Eq. (9) and Eqs. (21)–(22).

This collective enhancement can be understood by analyzing the photon distribution in the modes k in the reciprocal space. Let us focus first on a chain with an odd M. In proximity of the critical point the photons are strongly confined in the single mode a₀, see Fig. 2. As can be seen in Eq. (18), this mode is effectively the closest to its critical point. Given that the QFI is nonlinear with respect to the photon number, it is more favorable to concentrate all the photons into the most critical mode, rather than distributing them uniformly. This also explains the origin of the 1/2 prefactor in the scaling of the QFI for an even chain: when M is even the photons are confined into the two modes {a₀, a_π}, so the collective enhancement is achieved for two uniformly populated modes. In the following, we demonstrate that confining all photons into a single critical mode is advantageous also when considering saturation.

Saturation for finite nonlinearity

Let us now assess the validity of the Gaussian approximation, which is strictly valid for a vanishing nonlinearity (χ = 0) [see Eqs. (1) and (2)]. In any physical implementation, nonlinear effects will inevitably play a role at high energies, and non-Gaussian properties emerge⁷⁴. For finite values of χ, as the critical point is approached, the Kerr term induces saturation in the photon number and, consequently, in the QFI. For a single critical sensor, the maximum achievable QFI for finite χ can be evaluated numerically⁵⁵. However, numerical simulations of the full quantum model for large chains are beyond the capabilities of classical computers. We then take a perturbative approach to assess the saturation point of the QFI as a function of the ratio between the frequency ω and the nonlinearity χ. For the single-mode system, we applied perturbation theory directly on the QFI, the details can be found in “Methods”. In this case, the validity region of the Gaussian approximation is given by

where N_loc is the number of photons in the unperturbed single-mode ground state. Using the bound in (36), we can establish an upper limit also for the QFI in relation to the physical parameters ω and χ. Indeed, recall that, for a single-mode critical resonator T ≈ 2N_loc/ηω, the scaling of the QFI can be written solely in terms of the number of photons as ({{{{mathcal{I}}}}}_{{{{rm{sm}}}}} sim 8{N}_{{{{rm{loc}}}}}^{4}/{omega }^{2}). Thus, for an ensemble of M independent critical resonators, the saturation of the QFI is determined by the bound

For the critical chain, we focus exclusively on the case where M is odd, as it yields the most favorable QFI scaling. Applying perturbation theory in this case is much more complex, as the Kerr terms in (2) are spatially local, but highly non-local in the reciprocal space. However, close to the criticality (when saturation effects are relevant), only a single critical mode is effectively populated. In “Methods”, we show that when x → 0 the Hamiltonian of the critical chain can be mapped onto a single critical mode, with a rescaled nonlinearity χ/M. In other words, for what concerns the critical scaling, the nonlinearity of the interacting chain is effectively diluted. Under single-mode approximation, we can calculate perturbations to the QFI in spite of the nonlocality of the Kerr terms in the reciprocal space. In addition, we confirm also the validity of this single-mode approximation by evaluating perturbative corrections to the energy gap between the ground and first-excited states in the multi-mode case. This perturbative analysis reveals that, for the critical chain, the Gaussian approximation holds true if

where now N = MN_loc is the number of photons in the multi-mode ground state. It follows that, given T ≈ 2MN_loc/ηω, the multi-mode QFI (28) is asymptotically upper-bounded by

Note that the asymptotic bound (39) is ({root 3 of{M}}) times larger than the single-mode bound (37). This shows that the collective enhancement of the critical chain is not canceled by saturation effects. On the contrary, the critical chain presents a scaling advantage also with respect to the ratio ω/χ.

Driven-dissipative case

Let us now consider an open-quantum-system setting. Notice that an in-depth analysis would require considering a microscopic derivation of the bath properties, which strongly depend on the specific physical platform and the specific readout mechanism. For example, in a chain of solid-state resonators only one or two sites would be coupled with input-output transmission lines. Furthermore, the multi-mode nature of the output field should be taken into account to assess the precision that can be achieved with continuous-measurement schemes^42,43. Here, our goal is to show the conceptual differences between the performances achievable with unitary dynamics and with an open quantum system. Accordingly, we will assume that each resonator is weakly coupled to independent Markovian environments at zero temperature. In this case, the system time evolution can be written in terms of the Lindblad equation

where Γ is the decay rate, and we defined

Also in the dissipative case, we can switch to the Fourier space. It is easy to show that ({sum }_{j}{{{{mathcal{D}}}}}_{j}[rho ]={sum }_{k}{{{{mathcal{D}}}}}_{k}[rho ]), which means that the dissipation acts locally on the operators defined in the momentum space. For the sake of simplicity, we will restrict the analysis to odd resonator chains and to the critical region. As previously discussed, when x → 0 only a single collective mode gives a relevant contribution to the QFI. Consequently, we will consider the dissipative evolution (40) of the k = 0 mode defined in Eq. (17). This simplification does not affect the main results, as including all modes can only increase the estimation precision.

The single-mode case has been thoroughly analyzed⁴⁴. Here, we are going to use those results to analyze the differences in the multi-mode case. So far in this work, we have considered static CQS protocols, where an adiabatic sweep is used to push the system in proximity of the critical point. In the dissipative case, static CQS protocols consists in exploiting the susceptibility of the system-steady state, achieved after a long evolution time determined by the inverse of the smallest Liouvillian eigenvalue. By computing the QFI over the steady-state manifold, for an odd chain we find ({{{{mathcal{I}}}}}_{O} sim {N}^{2}/{Gamma }^{2}={M}^{2}{N}_{{{{rm{loc}}}}}^{2}/{Gamma }^{2}). Whenever time is not considered as a resource, a quadratic enhancement in M can be found in static dissipative protocols. However, the time required to reach the steady-state scales as T_ss ~ MN_loc/Γ, meaning that ({{{{mathcal{I}}}}}_{O} sim M{N}_{{{{rm{loc}}}}}{T}_{{{{rm{ss}}}}}/Gamma). Consequently, when time is a resource no collective advantage is found over the use of M independent single-mode sensors (recall that ({{{{mathcal{I}}}}}_{{{{rm{ind}}}}}=M{{{{mathcal{I}}}}}_{{{{rm{sm}}}}})).

We now show that a quantum advantage in M can be found in the dissipative case even for finite evolution times, when a dynamical CQS protocol is used. In this case, the scaling of the QFI depends on the protocol duration time itself, so we fix a maximum time ({T}_{max }) in which the protocol can be performed. In this way, time is treated as a resource on the same footing of the local photon number, for which we impose a maximum value N_loc. We assume the system is initialized in the vacuum (rho (0)=leftvert 0rightrangle leftlangle 0rightvert). Then, we evaluate the instantaneous QFI on the state ρ(t) evolved according to Eq. (40). This evaluation can be done analytically at arbitrary time⁴⁴. However, the expression of the instantaneous QFI is too complex to show explicitly. In the following, we present approximate solutions for the unitary regime, and we analyze the general behavior in Fig. 5. For short times, we can expand the QFI in series up to the first order in Γ as

For sufficiently short times, we recover the quadratic scaling in both time and photon number observed in the unitary case (35). Indeed, from Eq. (42) we identify the regime in which the quantum collective advantage of the coupled-resonator chain is present also in the dissipative case:

It is important to notice that the time in which the quantum advantage is not spoiled by dissipation is bounded by the resources themselves, in agreement with fundamental bounds^71,72. Comparing Eqs. (43) and (35), one can also see the constant-factor advantage of dynamical over static CQS protocols³⁷, due to the adiabatic parameter η ≪ 1. Finally, in Fig. 5 we analyze the transition from the unitary to the dissipative regime, for different number of modes M. For sufficiently short ({T}_{max }), the dynamics is effectively unitary (black dots), so that ({{{{mathcal{I}}}}}_{O} sim {M}^{2}) while ({{{{mathcal{I}}}}}_{{{{rm{ind}}}}} sim M). Notice that, as ({T}_{max }) is increased, even though the quadratic advantage in M fades there is a wide transition during which ({{{{mathcal{I}}}}}_{O} sim {M}^{alpha }) with α > 1. This behavior is reminiscent with what has been predicted for uncertainties on the system parameters⁷⁵. Interestingly, from Fig. 5 we see that a constant-factor advantage in the collective case can still be observed for large M.

We have demonstrated that a critical quantum sensor based on a chain of M parametrically coupled Kerr resonators exhibits a quadratic enhancement of the QFI with respect to both photon and resonator numbers. Beyond asymptotic scalings, the coupled chain offers practical advantages in terms of physical parameters. This collective quantum advantage in critical quantum sensing is significant even for a small number of resonators. We have analyzed a variety of parameter regimes, and we summarize the QFI scaling for the most interesting cases relevant to quantum sensing applications in Table 1.

This work paves the way to the exploration of collective quantum advantage in critical quantum sensing, both at experimental and theoretical levels. The considered class of finite-component phase transitions has already been experimentally observed using atomic⁴⁵, polaritonic^46,47 and circuit-QED^{48,49,50,51,52} devices. Parametric resonators with extremely low static nonlinearity can be implemented using high kinetic inductance cavity arrays⁵³. More generally, nonlinear coupled-resonator chains can be implemented using a broad variety of solid-state platforms, including semiconducting⁷⁶, opto- and electro-mechanical systems^77,78 and superconducting quantum circuits⁷⁹, among others. Further theoretical work is needed to properly support the development of collective CQS protocols. The driven-dissipative case should be addressed considering platform-specific couplings with the environment, which will also depend on the chosen readout mechanism. The need for adiabatic processes, responsible for the η < 1 prefactor in the scaling of the QFI, can be lifted by designing dynamical protocols^37,38. Continuous-measurement schemes^42,43 can be designed to optimize information retrieval from the output signals. Finally, topological properties have been shown⁸⁰ to enhance the sensing precision in sensors based on non-Hermitian lattice dynamics. The inclusion of topological properties in critical chains could potentially enhance the precision and resilience to parameter fluctuations of collective CQS protocols. Besides quantum sensing and metrology, this work presents a phenomenological study of the ground-state properties of an unconventional quantum many-body model, which is critical both at the local and global levels. We analyzed the regime of weak nonlinearity before symmetry breaking occurs. Numerical simulations of the full quantum model are challenging due to large photon numbers and long-range correlations. Phase-space methods appear to be the most suitable for analyzing lattices of critical resonators with weak but finite nonlinearity, where the ground state of the system is approximately Gaussian.

Local number of photons

In this section, we provide details on the calculation of the local number of photons in the critical resonator chain. The local number of photons in the ground state is given by

where (| gleft.rightrangle) is defined in (19). We express the operators a_j in the reciprocal space as

and rewrite the expectation value in (44):

Thus, to obtain N_loc, we have now to compute the expectation values (langle g| {a}_{k}^{{{dagger}} }{a}_{{k}^{{prime} }}| grangle), knowing that the ground state is the multi-mode squeezed vacuum state of Eq. (19). Let us consider first the case ((k,,{k}^{{prime} }),ne ,0) and ((k,,{k}^{{prime} })ne pi). Then, the expectation value becomes

which is valid since, in the case in which (kne {k}^{{prime} }ne {k}^{{primeprime} }), all the operators commute, and the expectation values are null. Following the same reasoning, if k or ({k}^{{prime} }) are equal to 0 or π, the expectation value is always zero unless (k={k}^{{prime} }). It follows that

which does not depend on the particular site location j. The photons in the system are uniformly distributed across each cavity. We can further simplify (49) by expanding the discrete sum over k > 0:

From Eq. (19), the photon number in each collective mode a_k is ({N}_{k}={sinh }^{2}left(| {xi }_{k}| right)={s}_{k}^{2}), where s_k is defined as in (7) with the substitution (epsilon to epsilon cos (k)). Finally, since ({s}_{k}^{2}={s}_{-k}^{2}), we write the local number of photons as

Series expansions

In this section, we present the mathematical details concerning the trigonometric series that appear when computing the QFI and N_loc in the asymptotic limit x → 0. Let us first consider the QFI when M is odd. Expanding the QFI (22) up to the first order in x yields

A closed analytical expression for the QFI in this limit can be obtained using the following identities on trigonometric series⁸¹:

Thus, in our case we have

which then leads to (24). On the other hand, when M is even, expanding the QFI (21) leads to

From the symmetry properties of the cosecant functions and their periodicity, it follows

which then leads to (23).

Let us now focus on the average number of photons. For an even M, expanding (20) up to the first order in x yields

The summation of the cosecant functions can be approximated as⁸²

where γ is the Euler-Mascheroni constant. Thus, we obtain Eq. (25). Conversely, for an odd M, expanding in series (20) results in

Due to the symmetry properties of the cosecant function and its periodicity, we have

from which we retrieve Eq. (26).

Finite non-linearity: single-mode

In this section, we go through the details of the perturbation theory applied to assess the validity of the Gaussian approximation.

The squeezing Hamiltonian (3) is an approximation of the general model:

with χ representing the Kerr non-linearity. The Gaussian approximation is obtained by setting χ = 0. However, when the number of photons in the system state increases beyond a certain limit, dependent on χ and ω, the nonlinear term becomes significant and cannot be neglected. We take a perturbative approach to identify an upper bound on the allowed number of photons, and, consequently, estimate the saturation of the QFI.

Let us start by analysing the energy gap ΔE between the ground and the first excited state. Up to the first order in χ, it is defined as

where (Delta {E}^{(0)}=lambda =sqrt{{omega }^{2}-{epsilon }^{2}}) and (Delta {E}^{(1)}=langle e| {H}_{{{{rm{NL}}}}}| erangle -langle g| {H}_{{{{rm{NL}}}}}| grangle). Here, (| grangle =S(xi )| 0rangle) is the ground state, (| erangle =S(xi )| 1rangle) is the first excited state and H_NL = a^†a^†aa. We recall that the squeezing parameter is defined as (xi =ln left(s+sqrt{{s}^{2}+1}right)), with s given by (7). Since (| grangle) and (| erangle) are single-mode squeezed Fock states, ΔE⁽¹⁾ can be easily computed as

with (N={sinh }^{2}left(| xi | right)) the number of photons in the ground state. In the Gaussian approximation, we have ΔE = ΔE⁽⁰⁾. Thus, when χ > 0, it must hold that

By expressing λ in terms of N as λ = ω/(1 + 2N), we obtain an estimate of the desired upper bound on the number of photons in the ground state, i.e.,

This means that, as far as N is well below the upper bound (66), the Gaussian approximation is valid and the non-linearity can be safely neglected.

The same procedure may be applied directly to the QFI. In the following, we show that this leads to an upper bound with the same scaling in ω/χ as in (66) but slightly tighter. To apply the perturbation theory to the QFI, we use (8) with the ground state defined as

Here, (vert {g}^{(0)}rangle =S(xi )| 0rangle) is the unperturbed ground state, while (vert {g}^{(1)}rangle) is the first-order correction to the unperturbed ground state. Up to the first order in χ, the QFI can be expressed as

where ({{{{mathcal{I}}}}}^{(0)}) is the unperturbed QFI defined in (9), while ({{{{mathcal{I}}}}}^{(1)}=8{rm{Re}} left[langle {partial }_{omega }{g}^{(0)} | {partial }_{omega }{g}^{(1)}rangle right]) is the first order correction. Perturbation theory allows us to easily compute (vert {g}^{(1)}rangle). Indeed, in our case

with (vert {j}^{(0)}rangle =S(xi )| jleft.rightrangle). The sum is performed for j ∈ [1,  4], since when j > 4 all the expectation values are null due to the orthogonality of the Fock states. The transition elements in (69) can be then analytically evaluated as (vert {j}^{(0)}rangle) are single-mode squeezed Fock states. The first order correction to the ground state is then

The amplitude probability functions v_ω,ϵ and h_ω,ϵ are defined as

and depend on ω and ϵ according to N = N(ω, ϵ) and λ = λ(ω, ϵ). If we now rewrite the squeezing operator as (S(xi )={e}^{,{f}_{omega ,epsilon }O}), where f_ω,ϵ = ∣ξ∣/2 and O = a² − a^†2, the derivatives with respect to the parameter ω read

where the apex stands for ∂_ω. Finally, using (72), we are able to compute the first-order correction to the QFI:

From the definition N = s², we obtain an expression for ϵ in terms of the unperturbed ground state number of photons, i.e., (epsilon =2omega sqrt{N+{N}^{2}}/left(1+2Nright)). With this in mind, the leading terms of the zeroth-order QFI and its first-order correction are

In the Gaussian approximation, it holds that ({{{mathcal{I}}}}={{{{mathcal{I}}}}}^{(0)}). Thus, when χ > 0, we have to impose (| chi {{{{mathcal{I}}}}}^{(1)}| ,,ll, ,| {{{{mathcal{I}}}}}^{(0)}|), from which we obtain a new estimate for the upper bound on the number of photons in the ground state:

This bound has the same scaling in ω/χ as that in (66), but is tighter by roughly a factor of two.

Finite non-linearity: resonator chain

For the multi-mode case, the Hamiltonian in (16) can be thought of as the Gaussian approximation of the general model

where χ represents the local Kerr non-linearity. To assess the validity of the Gaussian approximation, we can apply perturbation theory, as we did for the single-mode case. In the following, we consider M odd, since this choice leads to the most favorable scaling of the QFI. Close to the criticality, the majority of the photons are then confined in the mode a₀, leaving the other modes a_±k unpopulated. Thus, we approximate the unperturbed ground state defined in (19) as

This allows us to map the multi-mode system to a single-mode system. The low-energy excited states can be written as (vert {j}^{(0)}rangle approx {left(j!right)}^{-frac{1}{2}}{S}_{0}({xi }_{0}),{({a}_{0}^{{{dagger}} })}^{j}| 0rangle). The single-mode system has now N = MN_loc photons in the mode a₀, where N_loc is the local number of photons in each site j. Let us now focus on the effect of the nonlinear terms, which can be written in the reciprocal space as

Given the structure of the low-excited states, the expected value and transition elements of the nonlinear term are nonvanishing only for k₁ = k₂ = k₃ = k₄ = 0. In this single-mode approximation, the Hamiltonian (76) can be mapped onto a single critical resonator with a rescaled nonlinear term, such that

We can then perform the same calculations we did for the perturbation of the QFI in the single-mode case, but with the ground state (77) and a rescaled nonlinear Hamiltonian. It follows that the bounds found for the single-mode case are still valid for the multi-mode one, provided that we replace N → MN_loc and χ → χ/M. Finally, to ensure the validity of the Gaussian approximation in the multi-mode scenario, we require

So far, we have assumed that the unperturbed multi-mode ground state can be written as (77). This allowed us to evaluate the first-order correction to the QFI. In this section, we analytically compute the first-order correction to the multi-mode energy gap, without resorting to Eq. (77). In the limit where N ≈ N₀, with N the number of photons in the ground state and N₀ the number of photons in the mode a₀, the single-mode bound (66) is retrieved, thus endorsing the validity of the single-mode approximation. The multi-mode energy gap between the ground and the first excited state is defined as in (63), with (Delta {E}^{(0)}={lambda }_{0}=sqrt{{omega }^{2}-{epsilon }^{2}}) and (Delta {E}^{(1)}={E}_{e}^{(1)}-{E}_{g}^{(1)}), where ({E}_{e}^{(1)}) and ({E}_{g}^{(1)}) are the first order corrections to the eigenvalues of the first-excited and ground state, respectively. Let us start by computing the first-order correction to the ground state ({E}_{g}^{(1)}). By using (78) and the relabeling ({k}_{1}to {k}^{{prime} }+q), k₂ → k − q, ({k}_{3}to {k}^{{prime} }), we write ({E}_{g}^{(1)}) as

where we have used ([{a}_{{k}^{{prime} }},{a}_{k-q}^{{{dagger}} }]={delta }_{{k}^{{prime} },k-q}) and (MN={sum }_{k,q}langle g| {a}_{k}^{{{dagger}} }{a}_{k}| grangle). Since the ground state (| grangle) is defined as in (19), any expectation value (langle g| {a}_{k}| grangle) is null. This means that the non-zero terms in (81) are those in which the ladder operators refer to the same mode in pairs, leading to the constraints: (i) ({k}^{{prime} }+q={k}^{{prime} }) and k − q = k, meaning q = 0; (ii) ({k}^{{prime} }+q=k-q) and ({k}^{{prime} }=k), meaning q = 0 and (k={k}^{{prime} }); (iii) ({k}^{{prime} }+q=k), meaning (q=k-{k}^{{prime} }). Implementing these constraints in (81) allows us to obtain a much simpler expression:

On the other hand, the first excited state in the multi-mode system is

Thus, the expression (82) holds also for ({E}_{e}^{(1)}), provided that (| grangle to | erangle) and N → N + 1 + 2N₀ (it is straightforward to verify that the number of photons in (| erangle) is indeed N + 1 + 2N₀). The first-order correction to the multi-mode energy gap is then

Since the first excited state (| erangle) differs from ground state (| grangle) only in the terms related to the mode a₀, the differences between expectation values in (84) will always give zero unless k = 0 or ({k}^{{prime} }=0). Then, we can further simplify ΔE⁽¹⁾:

The remaining expectation values in (85) can be computed easily, being (| erangle) and (| grangle) single-mode squeezed Fock states for what concerns a₀. We have that (langle e| {a}_{0}^{{{dagger}} }{a}_{0}{a}_{0}^{{{dagger}} }{a}_{0}| erangle =1+12{N}_{0}+15{N}_{0}^{2}) and (langle g| {a}_{0}^{{{dagger}} }{a}_{0}{a}_{0}^{{{dagger}} }{a}_{0}| grangle =2{N}_{0}+3{N}_{0}^{2}). Finally, we write

and, up to the first order in χ, the multi-mode energy gap is given by

Close to the critical point, we have N = MN_loc ≈ N₀, so that

We retrieve the same result obtained in the single-mode case in (66), with a larger number of photons N = MN_loc and a rescaled non-linearity χ/M.