Continuous sensing and parameter estimation with the boundary time-crystal

A boundary time-crystal is a quantum many-body system whose dynamics is governed by the competition between coherent driving and collective dissipation. It is composed of N two-level systems and features a transition between a stationary phase and an oscillatory one. The fact that the system is open allows to continuously monitor its quantum trajectories and to analyze their dependence on parameter changes. This enables the realization of a sensing device whose performance we investigate as a function of the monitoring time T and of the system size N. We find that the best achievable sensitivity is proportional to $\sqrt{T}N$, i.e., it follows the standard quantum limit in time and Heisenberg scaling in the particle number. This theoretical scaling can be achieved in the oscillatory time-crystal phase and it is rooted in emergent quantum correlations. The main challenge is, however, to tap this capability in a measurement protocol that is experimentally feasible. We demonstrate that the standard quantum limit can be surpassed by cascading two time-crystals, where the quantum trajectories of one time-crystal are used as input for the other one.

Dissipative time-crystals provide an example of collective phenomena in which dissipation plays a constructive role.This contrasts its more archetypal role in which dissipation simply destroys quantum correlations, hence being usually detrimental for, e.g., quantum metrology protocols [36,37].Alternative quantum metrology approaches try to harness nonequilibrium phenomena in order to enhance the sensitivity of parameter estimation [38].For instance, this is the case of protocols exploiting dissipative phase transitions [38][39][40][41][42][43][44][45][46].In this regard, another key idea is to exploit the information contained in the emissions of open quantum systems via continuous monitoring protocols [38,[47][48][49][50][51][52][53][54], see a sketch in Fig. 1.Bounds to the fundamental sensitivity achievable by these continuous sensors have been derived [55,56].[14,34], as depicted here, but also free space implementations are possible [35].The two-level systems are driven with a Rabi frequency ω, which represents in our setting the parameter to be estimated.We consider two protocols, a single time-crystal (protocol I) and a two cascaded time-crystals (protocol II), where one system is driven with Rabi frequency ω and the other with ωD.In both cases the parameter ω is estimated by measuring the photon output intensity IT over a time interval of length T or in the stationary limit (T → ∞).The estimation error δω(T ) -and whether it saturates the theoretical sensitivity bound Sω -depends on the observation time and the protocol.
In general, they are difficult to saturate with bare photocounting or homodyne detection protocols [cf.Protocol I in Fig. 1].However, strategies have been developed in order to improve the sensitivity of continuous sensors, based on Bayesian methods [57,58], machine learning techniques [59], or the use of auxiliary systems [60,61].Connecting to the latter idea, Ref. [60] provides a general protocol based on cascading the output [62][63][64] of the sensor to a replica system that is continuously monitored [cf.Protocol II in Fig. 1], which may enable the saturation of the fundamental bound.
In this work we show how dissipative time-crystals can be exploited for sensing applications, reaching a sensitivity which can surpass the standard quantum limit.Recent works have studied BTCs from the perspective of critically enhanced sensing [44,45], focusing on properties of the system alone, while Ref. [65] considered a discrete time-crystal for sensing time-dependent fields.Here, we instead assess the performance of BTCs as continuous sensors.The rationale is that time-crystal oscillations clearly manifest in the photocounting and homodyne detection signals even for finite sizes [66] and that this output is readily accessible in experiments.Firstly, we analyze the fundamental (theoretical) bound on the achievable sensitivity with these devices.Subsequently, we consider two different sensing protocols, based on photocounting experiments, see Fig. 1.Protocol I entails the direct monitoring of the signal while Protocol II relies on indirect monitoring of the photocurrent through a cascaded replica of the BTC, acting as a decoder, in the spirit of Ref. [60].The time-crystal phase offers an enhanced sensitivity bound which scales linearly with the particle number N .This theoretical sensitivity cannot be achieved by Protocol I. Protocol II, on the other hand, allows to achieve the scaling ∼ N 0.78 , thus surpassing the standard quantum limit [38].
The model.-The BTC is composed of N spin-1/2 particles and described by the master equation with D[x]ρ = xρx † − {x † x, ρ}/2 and ρ being the state of the system.We further defined Ŝα = 1 2 N j=1 σ(j) α (α = x, y, z), with σ(j) α being the Pauli matrices and Ŝ± = Ŝx ± i Ŝy .Eq. (1) thus encodes collective spin decay with rate κ and a (resonant) driving with Rabi frequency ω.It preserves the total angular momentum, and we focus throughout on the fully symmetric sector.The above model was introduced in the context of cooperative resonance fluorescence, see, e.g., Refs.[67][68][69][70], and was recently recognized as a simple dissipative model displaying a time-crystal phase [18].
For collective systems such as that of Eq. ( 1), it is customary to rescale the collective decay rate κ by the system size in order to enforce a well-defined thermodynamic limit [22,71].However, we focus on finite-size systems and thus consider a N -independent rate, which further allows for a closer connection with experiments [35].The system displays a stationary regime, characterized by fast relaxation to the stationary state, and an oscillatory regime, featuring long-lived oscillations [70].The quality factor of the oscillations increases with system size [70] and diverges in the thermodynamic limit, where the system approaches the time-crystal phase [18,22].The two regimes are sharply separated by a critical Rabi frequency, ω c = κN/2, as system size increases [70], the long-lived oscillations emerging for ω > ω c .
Continuous sensing.-Our goal is to exploit the above system as a continuous sensor for estimating the Rabi frequency ω.To this end, we consider sensing protocols based on photocounting, so that the quantity that is measured is the time-integrated photon count or output intensity I T up to a measurement time T .The latter is defined as I T = 1 T T 0 dN (t), where dN (t) is a random variable, taking the value 1 when a photon is detected at time t and 0 otherwise [72], with average value given by E[dN (t)] = κdtTr[ Ŝ+ Ŝ− ρ].Here, E[•] represents the average over all possible realizations of the photocounting process [72].A relevant figure of merit for the sensitivity of a protocol is the estimation error, or error propagation formula, which, for the above-introduced output intensity, reads The first term on the right-hand side of the above equation is the standard deviation of the measurement, while the second one represents the susceptibility of the intensity on ω.In practice, the estimation error is the inverse of the signal to noise ratio.The quantum Fisher information of the emission field (QFI), F E , provides a lower bound for continuous sensing protocols [39], which applies to any protocol exploiting the information provided by the emitted field [39,55,60].
Here, we focus on a long measurement time limit, in which the system is described by its stationary state ρss .Since I T obeys a large deviation principle [73,74], the time-integrated intensity tends to its mean [72], lim T →∞ I T = κTr[ Ŝ+ Ŝ− ρss ], and the standard deviation scales, away from phase transitions, as As a consequence, the estimation error asymptotically behaves as Both the prefactor for the standard deviation σ I T and the one for the estimation error δω are time-independent quantities (see also Supplemental Material [75]).The QFI scales linearly with time [55,75], and thus The quantity S ω in the above equation is thus the (theoretical) sensitivity bound for the estimation, i.e., (δω) −1 ≤ S ω .
Fundamental bound on sensitivity.-In Fig. 2 we analyze the theoretically achievable bound on the sensitivity, S ω .Deep inside the stationary phase this bound is constant, S ω ≈ 2/ √ κ, i.e., it does not depend on ω.This result is obtained via a Holstein-Primakoff approach (HP) (see [75]), which shows that spins organize in a large displaced state.Its fluctuations are annihilated by the dissipator and hence they do not contribute to the properties of the emitted light at leading order.The QFI is in this case dominated by the amplitude of the coherent displacement.
In the vicinity of ω c we observe a sharp crossover with the sensitivity bound becoming larger in the oscillatory phase (ω/ω c > 1).This feature becomes more and more pronounced the larger the system size.The sensitivity in the time-crystal phase indeed grows linearly, S ω ∝ N [see Fig. 2(b)], well into the time-crystal regime (while at criticality the scaling is sublinear).From Fig. 2(a) it is also clear that for fixed system size N the sensitivity bound saturates upon increasing ω/ω c .This shows that increasing the system size and increasing the Rabi frequency do not have the same effect on the sensitivity.Therefore, the N -dependency of the bound in the time-crystal phase is a genuine many-body effect.Notice that S ω ∝ N corresponds to a N 2 -scaling for the QFI, meaning that the time-crystal phase theoretically offers a Heisenberg limited sensitivity in the number of particles [38].
The system size enhancement of the sensitivity can be understood from the properties of the emergent manybody oscillations.As shown in Ref. [66], these oscillations translate directly into the photocounting signal, manifesting as an oscillatory detection signal.Increasing system size with fixed ω/ω c has a twofold effect.First, the quality factor of the oscillations increases linearly with N [20,22,70], and thus the oscillations in the emitted field have an increasingly better defined frequency.Second, the amplitude of the oscillations also increases with system size (as more atoms emit synchronously), and thus the signal stands out more clearly from background noise [66].The combination of these effects creates stronger correlations between the system and the emitted field, which are a known source for enhancement of the QFI of the emission field [39,55,60].An indirect signature of such strong correlations is the fact that the reduced state of the system in the time-crystal phase is close to the maximally mixed state [66,76], which witnesses that the total state of system and emission field is highly correlated.In the following, we explore whether sensitivities close to the bound can be achieved using sensing protocols based on photocounting.
Protocol I. -The first protocol we consider is based on counting the number of emitted photons in a time window T .For long measurement times, a single (ideal) measurement run yields the stationary state intensity, while its standard deviation can be systematically studied using large deviations [75].In Fig. 3(a) we present the estimation error for this protocol.The smallest estimation errors are attained in the stationary phase.Indeed at the transition the estimation error increases steeply and in the time-crystal phase it assumes values which are much larger than the bound S ω .The "bad" performance under this measurement protocol can be understood from the data shown in Fig. 3(b).There we observe that in the time-crystal phase the derivative of the intensity with ω diminishes significantly while the standard deviation displays the opposite trend.
Interestingly, simple photocounting saturates the sensitivity bound in the stationary phase.This can be understood by applying the HP approximation to methods from large deviation theory, which allows us to obtain an analytical approximation for the error [75]: δω ≈ 0.5 √ κ.This value agrees with the numerical results of Fig. 3(b) (see black-broken lines).The physical interpretation is the same as for the QFI, i.e., to leading order the emitted light displays the same properties as a state with large coherent amplitude.
Protocol II.-In the second protocol, inspired by Ref. [60], the emission of the system of interest, the sensor, is cascaded [62,63] into an auxiliary system with same degrees of freedom, referred to as decoder (see Fig. 1).This protocol builds on the emergence of dark states in the cascaded system, a situation explored in detail in Ref. [64].Optimal sensing is then achieved close to the parameters in which the dark state is present [60].We consider two BTCs in cascaded configuration: all output light of the sensor (system 1) is fed into the decoder (system 2), while no light of decoder returns to the sensor [62][63][64].Protocol II is described by the master equation: where Ĥc = −iκ( Ŝ Before analyzing the sensitivity of the second protocol we briefly summarize the most important features of the cascaded dynamical system.A crucial result is that the stationary state of Eq. ( 6) is a dark state -for any ratio ω/κ -as long as ω D = ω (see [75] for its analytical expression).Interestingly, well into the stationary phase, it is separable, while for ω/κ 1 it tends to the highly correlated singlet state of the total angular momentum.Once we break the condition ∆ω = ω − ω D = 0, the cascaded system emits light.In Fig. 4(a) we show its stationary emission intensity: lim For ω D /ω c > 1 and ∆ω > 0, the system emits much more intensely than for ω D /ω c < 1.Looking at the purity of the stationary state [Fig.4(b)], it is possible to see a clear separation between these two regions, the brighter one being highly mixed, while the darker one displaying an almost pure stationary state.Notice that the reduced state of the sensor is exactly the same as for Eq. ( 1) [64].
Hence, the mixed bright phase corresponds to the sensor being in the time-crystal phase, while the pure phase corresponds to the sensor being stationary.Using the HP approximation [75], we find that the stationary phase is assumed when ω, ω D < ω c and ∆ω < (ω c − ω D )/2.This result is accurate even for small sizes [see red-dashed line in Fig. 4(a-b)].In this phase, each collective spin organizes itself in a coherent state with large amplitude plus Gaussian bosonic fluctuations.Their stationary state is a two mode separable state, in which each mode displays analogous properties [75].
We proceed by analyzing the sensitivity of Protocol II, which estimates ω from the output intensity of the cascaded system [see Fig. 1].In Fig. 4(c) we provide the estimation error as function of ∆ω for different points of the phase diagram.A first observation is that in the stationary phase there is no improvement with respect to Protocol I.This can be understood again via HP in combination with methods from large deviations (see [75]) which predicts the estimation error to be: δω ≈ 0.5 √ κ.The physical interpretation of this result is the same as for the individual system: the emitted light has essentially the same properties as a coherent state.In contrast, in the cascaded time-crystal phase the sensitivity increases.The estimation error is smaller than in the stationary phase and it improves further with increasing system size, as shown in Fig. 4(c-d).In panel (d) we analyze the system size dependence considering two values of ∆ω close to the dark state condition.By fitting the points for the largest system sizes N (dashed lines), we obtain δω ∝ N −α .The exponents are α = 0.70 for ∆ω/κ = 0.01 and α = 0.78 for ∆ω/κ = 0.005.The sensitivity thus surpasses the standard quantum limit in the number of particles.For large system sizes the system performs better closer to the dark state point ∆ω = 0.The estimation error displays a non-monotonous behavior with respect to ∆ω, which reflects the presence of an optimal value for δω.This can be understood by analyzing separately the two quantities that contribute to the estimation error [see Fig. 4(e)]: as system size increases, the derivative of the stationary intensity develops an increasingly sharper peak closer to zero.In contrast, the (time rescaled) standard deviation of the intensity increases monotonously with system size and ∆ω.The interplay of these two effects makes the optimal estimation error to be found for an intermediate ∆ω that depends on N .
Conclusions.-We have investigated parameter estimation through continuously monitoring quantum trajectories of a BTC.We have shown that the time-crystal phase in principle offers enhanced sensitivity, which manifests through the Heisenberg scaling of the QFI (∝ N 2 ) of the emitted field in the number of particles, N .Investigating two protocols, we have shown that a sensitivity surpassing the standard quantum limit can indeed be practically achieved in a cascaded setup of two timecrystals.Here one time-crystal acts as a sensor and the other one as a decoder.In the future, it would be also interesting to apply Bayesian inference methods [57] in order to improve sensitivity or even saturate the theoretical sensitivity bound [60].In view of recent experiments [35], it would furthermore be interesting to consider other continuous monitoring protocols, such as homodyne detection and/or to include finite detection efficiency.

CALCULATION OF THE QFI OF THE EMISSION FIELD
The QFI of the emission field for an arbitrary measurement time T is given by [55]: Here ρg1,g2 (T ) is the solution of the following deformed master equation [39,55]: with and initial condition ρg1,g2 (0) = ρ(0).Here we are assuming that the parameter g can be encoded both in the Hamiltonian and jump operators.In this work, we will focus on sensing in the stationary state, i.e. in the long time limit, for which the procedure to obtain the Fisher information simplifies to [39,55]: where λ E (g 1 , g 2 ) is the dominant eigenvalue of L(g 1 , g 2 ), i.e., the one with the largest real part.As long as the stationary state is unique, for large times the Fisher information of the emission field increases linearly.Thus we can see the right hand side of Eq. (S4) as the rate at which this Fisher information increases with time in the asymptotic regime: i.e., the maximum possible rate at which sensitivity can be increased by increasing the measurement time.

LARGE DEVIATION APPROACH TO PHOTOCOUNTING
In this work we consider systems described by the standard Markovian master equation in Lindblad form [72]: where Ĥ is the Hamiltonian of the system and L the jump operator.This master equation can be unraveled in a non-linear stochastic Schrödinger equation (SSE), that describes an ideal (unit detection efficiency) photocounting monitoring process [72]: The SSE describes the conditioned (by the measurement protocol) state of the system |Ψ(t) , and expected values with the subscript 'pc' are taken with this state.The quantity dN (t) is a random variable taking the value one if a photon is detected at time t and zero otherwise.Its average value with respect to the instantaneous state in a trajectory is given by: Collecting the number of detections for a measurement time window T , we can define the emission or output intensity as The time-integrated count record is the central quantity of our sensing protocols.In the long-time measurement limit, it approaches the emission intensity in the stationary state of the master equation: lim In the limit T → ∞, the average over realizations in Eq. (S9) can be dropped out, i.e. the time-integrated quantity tends to its mean, as in a law of large numbers.The way in which this occurs can be understood within the framework of large deviations theory.This is of upmost importance in our sensing approach, as it quantifies how the estimation error changes with measurement time.In Ref. [73] a large deviation approach for open quantum systems monitored by ideal photocounting is presented.The central quantity is the partition function: from which we can obtain the moments and cumulants of I T as partial derivatives with respect to the conjugate field s around s = 0. Of particular interest for us are the mean and the variance of the emitted intensity: The partition function can be calculated as Z T (s) = Tr[ρ s (T )], where ρs (T ) is the solution of the tilted master equation [73,74]: This defines the tilted Liouvilian, ∂ T ρs = L(s)ρ s .This generator is not trace preserving, though it preserves positivity.Therefore, its dominant eigenvalue, denoted by θ(s), is real and generally non-zero.This eigenvalue corresponds to the (stationary) scaled cumulant generating function, and one can obtain all cumulants from its relation with the partition function: Thus, for large measurement times T : lim T →∞ As θ(s) does not depend on T , in the long-time limit the variance diminishes as T −1 .This means that the estimation error diminishes as √ T , as stated in the main text.Finally, we find that the prefactor for the estimation error can be written entirely in terms of the scaled cummulant generating function: This formula combined with a Holstein-Primakoff approach can be used to obtain (approximate) analytical expressions for the estimation error in the stationary phase.
equation in the dark state condition ω = ω D for various points in paramater space, finding that it always displays a unique stationary state.Therefore, we conclude that the dark state obtained from Eqs. (S20) (S22) and (S25) is the stationary state of the system.
Benchmarking with known results.-As a first check we find that for S = 1/2, i.e., two spins, a 0 = √ 2, and thus: This is the same result as found in Ref. [64] (notice a differing factor 2 due to the fact that they are defining the Hamiltonian with Pauli matrices).A second interesting scenario is that of the limits ω/κ = 0 and ω/κ 1.In the former case |Ψ = |2S, −2S , i.e., all the spins are in the down state.In the latter, we obtain |Ψ ≈ |0, 0 , which when tracing out the decoder yields the identity matrix for the sensor.Both of these limits coincide with what is known for the individual BTC [66].
Magnetization of each collective spin.-The dark state of Eq. (S20) has some interesting properties with regards to the individual spin expectation values.In particular, we find that: where . . .Ψ denotes an expected value taken with respect to the dark state.The first property [Eq.(S28)] comes from the fact that the application of Ŝ(j) does not change the value m 1 + m 2 = −J and thus it does not connect different sectors of total z-component.The second property [Eq.(S29)] comes from the fact that the application of Ŝ(j) − changes by one the sum m 1 + m 2 and thus it can only give non-zero contributions when sandwiched by the state of the neighbouring sector J + 1.
The states |J, −J are either symmetric or antisymmetric with respect to the exchange of the labels 1 and 2, i.e. with respect to which is the sensor and which is the decoder.This is because of the following property of the Clebsch-Gordan coefficients: C J,M S,m1;S,m2 = (−1) J−2S C J,M S,m2;S,m1 . (S30) By extension, we find the following two properties: C J,M S,m1;S,m2 = (−1) 2(J−2S) C J,M S,m 2 ;S,m 1 C J,M S,m2;S,m1 , (S31) and From Eq. (S31) it follows that the matrix elements J, −J| Ŝ(j) z |J, −J are invariant under exchange of the label j.Instead, Eq. (S32) tells us that the matrix elements J, −J| Ŝ(j) − |J − 1, −J + 1 change sign under exchange of the label j.Therefore, going back to Eqs. (S28) and (S29) we find that: Moreover, since A * J A J−1 is purely imaginary [see Eq. (S25)], we obtain: These results allow us to apply the mean-field theory for the individual BTC to the cascaded system in the stationary phase.In particular, we know that the reduced state of the sensor is the same of the individual BTC.Thus, the mean-field approximation for ω/ω c < 1 yields (see e.g.[22]): For the decoder, using the derived properties, we find: These results are benchmarked in Fig. S1 by solving the master equation numerically.We observe that the ycomponents invert while the z-ones are identical.We also find that the mean-field results are reasonably accurate given the fact that we are considering small sizes.

HOLSTEIN-PRIMAKOFF APPROACH FOR THE INDIVIDUAL SYSTEM
The main results of this section are approximate analytical expressions for the estimation error [see Eqs.(S56) and (S57)] and the QFI [see Eqs.(S78) and (S79)] in the stationary phase.Moreover, we numerically benchmark these expressions in Figs.S2 and S3, respectively, finding good agreement.In order to obtain such expressions, we make use of the Holstein-Primakoff (HP) approximation and we apply it to the calculation of dominant eigenvalue of the tilted master equation [Eq.(S13)] and of the deformed master equation [Eq.( S3)].This approach also provides us important physical insights on the interpretation of the results, as discussed in the main text.

Holstein-Primakoff approximation
The starting point of the HP approximation is to express the spin operators in terms of a bosonic mode (see also Refs.[44,79]): The bosonic mode is assumed to be in a large displaced state around which we study the fluctuations: where β is a complex field to be determined.We proceed to expand the spin operators in powers of the small parameter = 1/ √ S: For our purposes it is only relevant to write down the following explicit expressions: This error is defined as the difference between the exact result and the analytical result of Eq. (S56), both divided by the exact result.We observe that the error increases near the transition, however, it is generally much smaller than 1/N .
We make the following ansatz for the dominant eigenstate of the deformed master equation: We use |E ωj 0 to denote the stationary state of the fluctuations [Eq.(S50)] at the phase diagram point ωj .For the condition ω1 = ω2 = ω this state is the stationary state of the fluctuations.We are interested in obtaining an expression around this point, as the QFI corresponds to the curvature of this eigenvalue around it [see Eq. (S4)].
When applying the approximated spin operators to this state we must be careful, as the Bra and the Ket in Eq. (S70) describe the fluctuations at different points of the phase diagram i.e. at ω1 or ω2 .This implies that, e.
where mωj −,1 denote the HP approximated operators with a β determined for the point in parameter space with ω = ωj .Taking into account this, we can introduce our ansatz into the deformed master equation (S69) and expand order by order.Then, the right hand side of Eq. (S69) gives the following terms.The dominant term: The subdominant term: This term is zero because of the self-consistency condition [Eq.(S45)].We do not go beyond this order.Setting the eigenvalue problem in an analogous way as for the scaled cumulant generating function, we obtain that the eigenvalue associated to We now argue that since |E ω1 0 E ω2 0 | corresponds to the stationary state when ω1 = ω2 , then very close to this point it must correspond to the dominant eigenmode of the deformed master equation [Eq.(S69)] (if the eigenspectrum is gapped, as it is the case in the stationary phase).Based on this argument we make the association: λE (ω 1 , ω2 ) = λ0 (ω 1 , ω2 ), (S76)

FIG. 1 .
FIG. 1. Continuous time-crystal sensor.The boundary time-crystal consists of N two-level systems subject to collective dissipation.It can be realized through atoms confined in a cavity[14,34], as depicted here, but also free space implementations are possible[35].The two-level systems are driven with a Rabi frequency ω, which represents in our setting the parameter to be estimated.We consider two protocols, a single time-crystal (protocol I) and a two cascaded time-crystals (protocol II), where one system is driven with Rabi frequency ω and the other with ωD.In both cases the parameter ω is estimated by measuring the photon output intensity IT over a time interval of length T or in the stationary limit (T → ∞).The estimation error δω(T ) -and whether it saturates the theoretical sensitivity bound Sω -depends on the observation time and the protocol.

FIG. 2 .
FIG. 2. Sensitivity bound.(a) Sensitivity bound Sω [cf.Eq. (5)], as a function of the Rabi frequency ω.(b) Sensitivity bound Sω as a function of the system size N for different values of the Rabi frequency.The black dashed line is a fit of the six points with largest size, yielding Sω ∝ N 1.0 .

FIG. 3 .
FIG. 3. Protocol I. (a) Red solid line: estimation error δω [cf.Eq. (4)], as a function of the Rabi frequency ω for N = 30.The black-dashed line corresponds to the inverse sensitivity bound S −1 ω [Eq.(5)].These quantities are plotted in units of √ κ.(b) Left axis (red): absolute value of the derivative of the stationary intensity, as a function of the Rabi frequency.Right axis (blue): prefactor of the scaling of the standard deviation σI T , as a function of the Rabi frequency.Black broken lines correspond to analytical results for the stationary phase [75].

FIG. 4 .
FIG. 4. Protocol II.(a) Stationary intensity of the cascaded system varying ∆ω = ω − ωD and ωD.At the line ∆ω = 0 the emitted intensity is zero.This quantity is plotted in units of κ.(b) Purity of the stationary state of the cascaded system.In panels (a), (b) the red dashed line corresponds to the mean-field transition line.In both cases N = 10.(c) Estimation error δω [cf.Eq. (4)] varying the Rabi frequency difference between the spin systems ∆ω, and for three different points in the phase diagram.Color lines correspond to the inverse of the sensitivity bound S −1 ω [cf.Eq. (5)] displayed in Fig. 2(a).In this panel we have considered N = 6.(d) Color points: estimation error δω varying system size N for ωD/ωc = 2.The black dotted line corresponds to the sensitivity bound S −1 ω displayed in Fig. 2(b).The color dashed lines correspond to fits for the six data points with largest size, yielding the exponential laws depicted in the figure.The data set contains N ∈ [6, 22], where for N > 11 we have used quantum trajectories methods [75].(e) Left (red) axis: derivative of the (stationary) emitted intensity with respect to ω in the stationary state and as a function of ∆ω.Right (blue) axis: prefactor for the standard deviation of the output intensity, as a function of ∆ω.In this panel we have fixed ωD/ω2 = 2.
α are the total angular momentum operators.Here, Ŝ(1,2) α are macroscopic spin operators, each one associated with N two-level systems.
FIG. S1.Stationary state observables for the cascaded system.Single spin expectation values evaluated in the stationary state of the cascaded system: m(j) α ss = Ŝ(j) α ss/S.The parameters are fixed to the dark state condition ω = ωD with N = 10.The color lines correspond to the exact results solving the master equation.The black broken lines correspond to the mean-field prediction of Eqs.(S35) and (S36).

50 FIG. S2 .
FIG. S2.Benchmark approximate expression for the scaled cumulant generating function.(a) Scaled cumulant generating function [Eq.(S14)] varying s for ω/ωc = 0.5 and N = 30.The blue points are obtained diagonalizing the tilted master equation.The red dashed line corresponds to the analytical result of Eq. (S56).(b) Relative error on θ(s) as a function of ω/ωc, for different system sizes and s = −0.05.This error is defined as the difference between the exact result and the analytical result of Eq. (S56), both divided by the exact result.We observe that the error increases near the transition, however, it is generally much smaller than 1/N .