Dynamical phases and quantum correlations in an emitter-waveguide system with feedback

We investigate the creation and control of emergent collective behavior and quantum correlations using feedback in an emitter-waveguide system using a minimal model. Employing homodyne detection of photons emitted from a laser-driven emitter ensemble into the modes of a waveguide allows to generate intricate dynamical phases. In particular, we show the emergence of a time-crystal phase, the transition to which is controlled by the feedback strength. Feedback enables furthermore the control of many-body quantum correlations, which become manifest in spin squeezing in the emitter ensemble. Developing a theory for the dynamics of fluctuation operators we discuss how the feedback strength controls the squeezing and investigate its temporal dynamics and dependence on system size. The largely analytical results allow to quantify spin squeezing and fluctuations in the limit of large number of emitters, revealing critical scaling of the squeezing close to the transition to the time-crystal. Our study corroborates the potential of integrated emitter-waveguide systems -- which feature highly controllable photon emission channels -- for the exploration of collective quantum phenomena and the generation of resources, such as squeezed states, for quantum enhanced metrology.

Introduction. The development of techniques for the manipulation of matter with light has undergone rapid progress in the past decade [1][2][3][4][5][6][7][8][9]. This has enabled the creation of tailored quantum systems for the purpose of quantum simulation and information processing [10][11][12]. It has also opened an avenue for the investigation of novel phases of matter and of the emergence of collective quantum behavior as it appears in the vicinity of a phase transition [13,14]. A paradigmatic example is the (open) Dicke model, which describes the interaction of an ensemble of atoms with a single-mode light field [15]. It has received substantial attention in recent years, not only because of its fundamental theoretical interest, but also because it has been realized in various experimental platforms, such as atom-cavity setups [16,17]. An example of a novel many-body phase is a so-called time crystal. This is a phase of matter in which time-translation invariance is broken [18][19][20][21]. It has been theoretically predicted and analyzed in various scenarios, including disordered closed systems [22][23][24], dissipative systems in continuous time [25][26][27] as well as periodic driven-dissipative systems [28][29][30][31][32]. Moreover, this state of matter has been studied and characterized in a number of recent experiments [33][34][35][36].
In this work we discuss the many-body phases of a lightmatter system that is composed of emitters coupled to light modes which propagate in a waveguide. This setup, which is subject matter of many current experimental studies [37,38], has the appeal that the many-body dynamics can be reduced to a small set of degrees of freedom, and thus lends itself to a largely analytical treatment [39][40][41][42][43][44][45][46][47][48]. We show that by implementing an instantaneous feedback protocol that relies on the measurement of a quadrature of guided light [34,[49][50][51][52][53][54][55][56][57], a rich variety of dynamical phases can be achieved, among them a continuous time crystal. For a large number of emitters, the dynamical and static phases of the system are described by a set of macroscopic spin variables whose expectation values obey a closed set of mean-field equations. We employ the theory of quantum fluctuations [58][59][60][61] to FIG. 1. Measurement-based feedback. A chain of N emitters (modelled as two-level systems) is placed in the vicinity of a waveguide and driven resonantly with Rabi frequency Ω by an external laser with wave vector perpendicular to the chain. Photons with wavelength λ are emitted into the right-and left-propagating modes of the waveguide with equal emission rate γR = γL. When the separation between nearest neighbors in the chain, a, becomes commensurate with λ, the emission can be described as a collective process, where the chain of emitters is treated as a single macroscopic spin of length N/2. The measurement of the light emitted into the right-propagating mode via homodyne detection is fed back as a modulation of the laser field with strength g.
investigate the quantum correlations in the many-body state of the emitter ensemble. This mathematically rigorous approach, which becomes exact in the thermodynamic limit, has been used for isolated systems [62][63][64][65][66][67][68][69], to explore critical phenomena [70][71][72], as well as in open quantum systems [73][74][75][76][77], for instance to witness dissipative generation of entanglement in mesoscopic systems [78][79][80]. Here, we use it to investigate spin squeezing of the emitter ensemble and to demonstrate critical power-law dynamics of quantum correlations at the boundary between two different nonequilibrium phases. Our results show that emitter-waveguide systems exhibit surprisingly complex emergent behavior and allow to realize collective states of matter on demand. This is not only of fundamental interest, but given the potential realizability of such systems as integrated devices, our work may stimulate applications in technological devices, e.g. for sensing and metrology, which are enhanced by quantum manybody effects.
System and equations of motion. The considered setup is sketched in Fig. 1. Here, a chain of N emitters is placed in close vicinity to a waveguide, parallel to its longitudinal axis. An external laser field drives the emitters (modelled as two-level systems with ground and excited states |g and |e , respectively) resonantly, with Rabi frequency Ω. Photons are emitted from the chain into the guided left-and rightpropagating modes of the waveguide at a rate γ R and γ L , respectively. We focus on the case in which the coupling between the emitters and the waveguide is non-chiral, i.e., the emission of photons into both guided modes occurs at the same rate, γ R = γ L . This can be in practice ensured by an appropriate choice of the laser polarization and transition dipole moment of each emitter. To allow an analytical treatment, our minimal model neglects photon losses into unguided radiation modes, c.f. Refs. [40,41,48]. Choosing the laser wave vector perpendicular to the chain and the nearest neighbor distance a between the emitters commensurate with the wavelength λ of the emitted light, i.e. a = mλ with m = 1, 2, . . . , in fact suppresses such unwanted decay channels [45,47]. Under these conditions and employing the Born-Markov approximation, the dynamics of the emitterwaveguide system is described by a particularly simple master equation: Here we have introduced the collective spin operators J α = α are Pauli matrices, e.g. σ x = |g e| + |e g|, corresponding to the k-th emitter. Furthermore, J = J x − iJ y is the lowering operator of the collective spin. The first term of Eq. (1) describes the resonant laser excitation of the emitters. The second and third term describe the incoherent emission of photons into the left-and right-propagating modes of the waveguide, through the dissipator D(J)ρ = JρJ † − 1 2 J † J, ρ . The feedback scheme that we consider for the purpose of this work is based on the measurement of the phase quadrature P (t) of the light emitted into the right-propagating mode. The measured quadrature of the photocurrent at each time determines the feedback, which consists of an instantaneous modulation of the Rabi frequency corresponding to the laser field driving the system. As shown in, e.g., Refs. [52,55], the action of the feedback leads to the following modification to the master equation, where we have set γ R = γ L = γ/2. The dimensionless parameter g represents the feedback strength. The action of the feedback is accounted for by a modified jump operator corresponding to emission into the right-propagating mode and a two-axis counter-twisting Hamiltonian term [81,82]. Dynamical and stationary phases. In the following we study the thermodynamic limit, i.e. when N → ∞. To this end we consider the dynamics of the expectation value of the magnetization operators m α = J α /N with α = x, y, z. In the thermodynamic limit, these operators permit a classical description of the system, since [m α , m β ] = iε αβδ m δ /N ≈ 0 as N → ∞, and converge to their average value, m α → m α , for clustering states [61,83,84]. This implies that expectation values involving them factorize, e.g., m x m z → m x m z . Under the dynamics given by Eq. (2), the exact time-evolution of these operators, in the large N limit, is given by the mean-field equations [77,83] Here we introduced the parameter κ = 2g + 1 and also the rescaled decay rate Γ = γN , in order to ensure a well-defined thermodynamic limit [15]. Eqs. (3) conserve the length j of the average magnetization vector m = ( m x , m y , m z ), which, for the following considerations, we set to 1 2 : j 2 = m x 2 + m y 2 + m z 2 = 1 4 . Let us first consider the particular case where g = − 1 2 , i.e. κ = 0. Here, the mean-field equations of motion become particularly simple and one can, furthermore, identify a second constant of motion, C 0 = m x exp [Γ m y /(2Ω)], which simplifies the discussion considerably: throughout, we consider the initial conditions m x (t = 0) = m y (t = 0) = 0 and m z (t = 0) = − 1 2 , i.e. the emitters are initially in the ground state N k=1 |g k , implying C 0 = 0. Under this simplification the equations of motion of the average magnetization components m y and m z become that of a harmonic oscillator. It is apparent then that, for κ = 0, the feedback eliminates the influence of dissipation and the magnetization undergoes persistent oscillations at frequency 2Ω. However, as we show below, this is not the case at the level of quantum fluctuations. These operators are indeed affected by dissipative effects which do not manifest in the dynamics of the magnetization operators.
Away from this special point, the integration of the first two equations in (3) yields the constant of motion Again, we set C κ = 0, which corresponds to initial conditions with m x (t = 0) = 0. Hence, the dynamics is reduced the two coupled equations which can be further reduced to the second order differential equation This is Newton's equation of motion for a particle with a nonlinear restoring force, which supports oscillating solutions as long as the product of both brackets on the right hand side In Fig. 2(a) we show the corresponding phase diagram. One clearly recognizes the two boundaries which delimit the time crystal phase from a stationary phase where the z-component of the average magnetization assumes a finite constant value m z ss . Moreover, one also observes that for g = − 1 2 (κ = 0) the time crystal phase persists down to the limit of vanishing Ω. In Fig. 2(b) we show a cut through the phase diagram along g = 1 2 in order to compare the result of the mean-field calculation with the exact numerics for finite N obtained from Eq. (2). One clearly sees that the numerical curves approach the mean-field prediction with increasing number of emitters N . This data confirms the position of the critical point as predicted by the mean-field equation (indicated by the red circle). Here the gap of the master operator (2) scales as Λ ∝ N − 1 2 . Fig. 2(c) shows the time evolution of the z-component of the average magnetization vector. In the stationary phase this approaches quickly a stationary value which converges to the mean-field results as the number of emitters grows. In the time crystal phase the finite size simulations exhibit damped oscillations. As the number of emitters grows, the damping rate decreases and the mean-field solution predicts persistent oscillations in the large N limit.
Fluctuations and quantum correlations. We now turn to the analysis of quantum correlations. In particular, we seek to understand whether the feedback can generate spin squeezing, which is a measure of quantum correlations between the emitters and can quantify entanglement in manybody systems [85,86]. Squeezed spin states may find ap-plications in metrological protocols, where they are used to enhance the precision of measurements beyond the standard quantum limit [87]. Recently, spin squeezing in the steady state of dissipative systems has been explored in cavity QED setups [88] and in the non-equilibrium dynamics of an ensemble of superconduncting qubits [89]. In order to study quantum correlations within the emitter ensemble we will exploit the theory of quantum fluctuation operators [58][59][60][61], applied to open quantum systems [74][75][76][77][78][79][80]. For our model, this allows us to obtain rigorous analytical results in the thermodynamic limit.
The basic idea is to introduce operators describing fluctuations of the collective observables J α around their average value J α . Intuitively, these are given by , which are the quantum analogue of fluctuation variables, subject to central limit theorems [58]. To better understand their nature, we look at the commutators between fluctuations, which are proportional to magnetization operators. As such, these commutators converge to scalar quantities, since m δ → m δ for large N [61,83,84], showing that fluctuation operators are bosonic operators in the thermodynamic limit. In fact, it is also possible to obtain the normal modes for quantum fluctuations, which behave as position and momentum operators. To do so, we rotate the coordinate system in such a way that the z-axis aligns with the direction n = m/| m| of the average magnetization vector. Here, we find thatF z (the operator F z in the rotated frame) is a classical degree of freedom that commutes with the other fluctuations [76,77,83]. For the directions perpendicular to n, we can instead define "position" and "momentum" operators, To show that these bosonic fluctuation operators account for correlations between the emitters, we consider their covariance matrix. For the fluctuations F α , this is defined as Contrary to what happens to magnetization operatorswhich play the role, in this quantum context, of the sample mean variable of the law of large numbers -correlations between fluctuations, as given by the quantities in Eq. (8), are not vanishing for N 1. Moreover, for sufficiently clustering states [58][59][60][61], fluctuations possess a Gaussian distribution and the matrix elements Σ αβ are not divergent. In order to obtain the 2 × 2 covariance matrix of the fluctuationsx andp, we rotate Σ [Eq. (8)] into the new reference frame and extract the 2 × 2 principal minor rescaled by j −1 [83].
The construction so far is associated with the quantum state of the emitter ensemble at a fixed time. In order to find the time-evolution of the canonical fluctuation operatorŝ x andp we need to consider the time-dependence of the matrices s and Σ [Eqs. (7) and (8)] in the large N limit. These can be obtained from the master equation Eq. (2), and by exploiting results from Refs. [76,77]. As shown in the supplemental material [83], ones finds that Here A = (Γ/2) diag 1 + κ 2 , 2, 0 , X t,0 is the timeordered exponential of the matrix and X t,u = X t,0 X −1 u,0 . It is noteworthy that the dynamics of the fluctuations, as described by Eq. (9), cannot be obtained by linearizing the mean-field equations (3). This means that fluctuations are affected by a dynamical building-up of correlations which is not captured by magnetization operators [76]. Moreover, the shape of Eq. (9) shows that the dynamics preserves the Gaussian character of fluctuations [77].
We are now in position to analyze spin squeezing. Several measures exist [90][91][92][93] and all of these are based -in analogy with bosonic squeezing -on the Heisenberg uncertainty principle. Here we consider the quantity [82,93] where ∆ 2 J ⊥ denotes the variance of a collective spin operator, J ⊥ , living in the plane orthogonal to the principal direction n of the spin magnetization and the minimization is taken over all possible J ⊥ . In the large N limit, this measure converges to the bosonic squeezing parameter computed for the fluctuation operators [83], and one finds that where λ i are the eigenvalues of the bosonic covariance matrix ofx andp. The state is spin squeezed if ξ < 1. In Fig. 3(a-b), we show the behavior of the spin squeezing in the large N limit. In the time-crystal phase this parameter does not converge to a stationary value [see also inset in Fig. 3(b) and Fig. 3(c), which displays one component of the covariance matrix]. On the other hand, in the stationary phase, for g > − 1 2 , the entries of the covariance matrix reach a stationary state [ Fig. 3(c)]. Here, the squeezing parameter becomes smaller than 1 [inset in Fig. 3(b)], which witnesses the presence of quantum correlations. Interestingly, at the critical line separating the time-crystal phase from the stationary one, spin squeezing displays a power-law behavior with time ξ ∝ t −1 , suggesting a critical building-up of quantum correlations. Moreover, contrary to what one would expect from the phase diagram in Fig. 2(a), the stationary phases to either side of g = − 1 2 are different: for g > − 1 2 , fluctuations reach a stationary behavior and display non-trivial quantum correlations, while for g < − 1 2 , spin squeezing reaches a stationary value above 1 and no stationary state for fluctuations exists. This qualitative change in behavior is rooted in the dependence of the structure of the matrix G [cf. Eq. (10)] on m z , which switches from positive to negative at g = − 1 2 . Conclusions. We have shown that feedback control of a coupled emitter-waveguide system offers the possibility to realize and manipulate dynamical many-body phases, such as a time crystal. Our minimal model allowed us to perform a largely analytical investigation which not only shed light on the behavior of the order parameter but also permitted the a rigorous analysis of fluctuations. This has shown that feedback can control spin squeezing of the quantum state of the emitters, which is particularly pronounced near the transition to the time crystal phase, where it exhibits critical scaling. In an actual experimental setting the finite number of emitters and also the existence of photon decay channels outside of the waveguide will impact on the long time behavior. The phases of the idealized setting discussed here will then be observed in the transient of the dynamics, i.e. they are expected to become meta-stable [94]. This has to be considered when seeking to design protocols that exploit, e.g., squeezing in technological applications such as quantum metrology.
Acknowledgements. The research leading to these results has received funding from the European Union's H2020 research and innovation programme

RECASTING THE DYNAMICAL GENERATOR
In order to find the contribution to the dynamical equations of the different terms in the generator, we rewrite the latter in a convenient form. In addition, we work in the Heisenberg picture, where only operators are evolved. In particular, the Heisenberg representation of the generator appearing in equation (2) can be written as the sum of four different terms (S1) The first contribution contains the local terms of the Hamiltonian and is given by The second contribution is still a coherent Hamiltonian one, and accounts for the coherent all-to-all interaction between emitters generated by the feedback. It can be written as where h is a 3 × 3 Hermitian matrix. In our specific setup, such a matrix assumes the form The last two contributions are due to the dissipative part of the generator, which is split into two terms: where the matrices A and B are given by

MEAN-FIELD DYNAMICS
In this section we give some details on the dynamics of the magnetization operators for the emitters. Given the collective operators J α , these are defined as and, due to the scaling 1 N , such operators form a commutative algebra in the limit N → ∞. This can be seen from the fact that [m α , m β ] ∝ 1 N → 0, for N 1. In addition, the dynamics considered in the main text, as well as the initial state we focus on, are such that only very weak collective correlations between pairs of emitters [76,77] are present. These correlations are not inherited by the operators m α . As such, under the expectation · over the state, these operators m α obey a sort of law of large numbers and converge, for N 1, to their expectation value This also means that the expectation of any product of magnetization operators converges to the product of the limiting operators, i.e. the expectation value factorizes as Our aim is thus to obtain the equations governing the dynamics of the quantities m α in the thermodynamic limit. These are the equations appearing in the main text in Eq. (3). To this end, we first compute the action of the generator on the collective operators J α (using collective spin commutation relations): We do not explicitly compute A[J α ] since this is not relevant for the dynamics of m α . Using the above relation, we straightforwardly find the action of the generator on the m α : The approximate symbol appears because we have set to zero the commutator of two average operators and because we have dropped the term A[J α ]/N since both these terms converge to zero in norm, for large N . We have further defined In order to obtain the semi-classical equations of motion Eq. (3), which are exact in the thermodynamic limit, we consider the expectation of the above equation with respect to the quantum expectation · recalling that, for such operators, expectation values factorize [see Eq. (S6)]. In this way, we find specializing the above equation for α = x, y, z, we find the equations reported in Eq. (3) in the main text.

QUANTUM FLUCTUATION OPERATORS
In order to quantify quantum correlations in the emitter ensemble, as well as to derive their dynamical behavior, we need to be able to control correlations between different collective spin operators, J α , with α = x, y, z. We are interested in correlations of the form where Σ is nothing but a covariance matrix for collective spin operators. We want to look at the behavior of these correlations in the thermodynamic limit. To this end, it is convenient to exploit the formalism of quantum fluctuation operators [58][59][60][61]. For a given quantum expectation · , we can define the fluctuation operator F α (α = x, y, z) as This is essentially the collective operator J α renormalized with respect to its expectation value in the state, and rescaled by a factor of 1 √ N . This operator is a quantum version of the fluctuation variable studied in central limit theorems. Furthermore, note that As shown in Refs. [58][59][60][61], under appropriate conditions on the quantum state defining the expectation · , which are always satisfied in our analysis, these fluctuation operators converge, for large N , to bosonic field operators. To show this, it is convenient to look at the commutation relations between fluctuation operators. These can be derived from the finite-N commutation In the above relations, the definition of s αβ makes apparent that such quantity has the scaling of a magnetization operator and thus converges to a scalar multiple of the identity [see Eq. (S6)] under any expectation taken with the state · . In particular, the actual value to which this operator converges is This shows that the commutator between two fluctuation operators, in the large N limit, is actually a scalar multiple of the identity. Thus, fluctuation operators obey canonical commutation relations and can naturally be identified with bosonic operators. In our specific setting, we have three bosonic field operators whose commutation relations are encoded in the 3 × 3 antisymmetric matrix In particular, for the setup considered here, we can find a mapping from these three generalized bosonic field operators to the standard bosonic position-like and momentum-like operators. Indeed, one can always find a rotation R which brings s into its canonical forms with R being a real unitary matrix, and j = m x 2 + m y 2 + m z 2 . Basically, the rotation simply consists in aligning the z-axis of the reference frame with the principal direction of the magnetization, n = j −1 ( m x , m y , m z ). This rotation R thus identifies three new fluctuation operatorsF α whose commutation relations are F α ,F β = is αβ , and whose correlations are encoded in the appropriately rotated covariance matrix We can then construct effective position and momentum operatorŝ which are such that [x,p] = i. We also note that the operatorF z plays the role of a scalar (classical) random variable since it commutes with the rest of the fluctuation algebra. For the relevant quantum degrees of freedom,x andp, we can also compute the covariance matrixΣ as the 2 × 2 principal minor of the matrixΣ, rescaled by the factor j. This covariance matrix contains the correlations between collective operators in the directions perpendicular to the principal magnetization vector n. The smallest eigenvalue ofΣ thus corresponds to the minimal variance that a collective operator, perpendicular to the direction n, can have, further rescaled by j −1 . When multiplied by 2, the smallest eigenvalue ofΣ is thus nothing but the spin squeezing parameter ξ, evaluated in the thermodynamic limit of a large number of emitters.

TIME-EVOLUTION OF THE COVARIANCE MATRIX
Now that we have understood the structure of fluctuation operators for a fixed state of the quantum system, we need to be able to recover the covariance matrix Σ, at each time t. Essentially, we have to find how this matrix propagates in time according to the generator of the dynamics from a given initial condition. Before, though, note that, from Eq. (S9), it is clear that the dynamics of the anti-symmetric matrix s is completely specified by the dynamics of the magnetization operators.
To simplify the discussion, we introduce the matrix C = F α F β , whose connection with the covariance matrix is explicitly We start by taking the derivative of C. The time derivative of F α is a scalar quantity since the derivative can only act on the state. This gives recalling also that F α has a vanishing expectation on the state · , the time-derivative of the matrix C is simply determined by the generator L as We now study the action of the different terms of the generator, appearing in Eq. (S1), on the product of two fluctuation operators.
First of all, we notice that Then, considering that the expectation value of a single fluctuation operator is zero under expectation over the state, we can freely subtract terms like F α γ D L βγ Jγ √ N to the ones above, to reconstruct For the second term of the generator, the one with all-to-all interacting Hamiltonian, we also have We can thus focus on the action of H C on a single fluctuation operator. This gives (S12) To obtain the second equality, we have simply added and subtracted the proper expectation values, appearing in the last line of the above equation, in order to reconstruct fluctuation operators. We now divide the terms on the right hand side of the second equality into two contributions. We call H C the ones in the second line of the above equation and H C those in the third line. The first term can be understood by looking at commutation relation between fluctuation operators. Indeed, recalling the fact that s αβ = −s βα we can write where we have We now come to the second term; this can be rewritten as where we have The approximate symbol in Eq. (S14) is due to the fact that, in D C , we have considered m µ instead of m µ : this is only valid in the limit N → ∞ [see also discussion of Eq. (S5)]. In addition, when considering the expectation over the state of the term in Eq. (S11), we have that we can safely add or substract a scalar quantity to the term H C [F α ] using the fact that, in any case, F β has zero expectation. Overall, we have found that This concludes the contribution which comes from the Hamiltonian term of the generator.
We now turn to the dissipative contributions and start with the one encoded in A. We have To understand this contribution, we look at a single term in the above summation. This can be reduced to The first two terms on the right hand side of the above equality are zero when taking the expectation over the state and in the thermodynamic limit. This can be understood as follows. All terms in the above equation actually consist of the product of two operators which have the scaling of magnetization operators. For the first two terms, one of the two is given by a fluctuation operator having a further scaling 1 √ N , which indeed transforms it into an operator with a scaling 1 N . However, fluctuation operators are rescaled with respect to their average over the state so that, under any expectation over the state, the operator F α / √ N → 0. Concerning the remaining two terms [the ones in the second line of Eq. (S16)] we can use the fact that they are magnetization operators to argue that they converge to the product of their expectation. This leads to We are thus left with the second dissipative contribution, the one related to the matrix B, which acts on fluctuation operators as (S17) We divide the above term into two pieces. The first on the right hand side of the second equality we call it B while the second B . We focus on the first and, exploiting commutation relations of fluctuation operators, obtain under expectation over the state this contributes with We are then left with the second part of this term. This gives Under expectation we thus have Putting these results together, we obtain the exact differential equation for the evolution of C which, in thermodynamic limit N → ∞, becomes Finally, using the fact that Σ = [C + C T ]/2 we find [76,77] d dt Σ = −sAs + GΣ + ΣG T , with G = D + Q , and Q = Q C + Q B In our specific setting, the relevant matrices have the following form We recall here that all quantities m α are actually time-dependent and obey the system of differential equations (S7). Whenever they appear in Eq. (S22) they must be considered at the running time t. As such the actual solution of Eq. (S22) is given by where the matrix s(t) is the matrix s in Eq. (S9) evaluated at time t and X t,u is the time-order exponential of the matrix G, such that d dt X t,u = G(t)X t,u , and d du X t,u = −X t,u G(u) .