Synchronization of Interacting Quantum Dipoles

Macroscopic ensembles of radiating dipoles are ubiquitous in the physical and natural sciences. In the classical limit the dipoles can be described as damped-driven oscillators, which are able to spontaneously synchronize and collectively lock their phases. Here we investigate the correspond- ing phenomenon in the quantum regime with arrays of quantized two-level systems coupled via long-range and anisotropic dipolar interactions. Our calculations demonstrate that the dipoles may overcome the decoherence induced by quantum fluctuations and inhomogeneous couplings and evolve to a synchronized steady-state. This steady-state bears much similarity to that observed in classical systems, and yet also exhibits genuine quantum properties such as quantum correlations and quan- tum phase diffusion (reminiscent of lasing). Our predictions could be relevant for the development of better atomic clocks and a variety of noise tolerant quantum devices.


Introduction
Arrays of synchronized oscillators [1] are ubiquitous in biological [2,3], physical [4] and engineering [5] systems and are a resource for technological advances [6]. Although there has been significant progress in the study of synchronization in classical systems [7], the understanding of the same phenomena in the quantum realm remains limited. A major obstacle so far is the general problem of the exponential scaling of the Hilbert space with system size which makes calculations dealing with quantum arrays very challenging. In fact, current investigations have been limited to the exact treatment of arrays of a small number of coupled quantum oscillators [8][9][10][11][12][13][14][15][16][17][18], and large ensembles at the mean field level or by including quantum corrections perturbatively [19][20][21]. Highly symmetric situations with collective coupling mediated, for example, by a cavity mode [22][23][24], have also been studied.
Ensembles of radiating dipoles are a natural platform to study quantum synchronization, where coherence can be generated from an incoherent source. One might regard laser systems, where radiation is amplified by the stimulated emission of photons, as a prototypical example. However, lasing is fundamentally a distinct phenomenon from quantum synchronization. This can be seen from the fact that lasing is possible even in the absence of coupling between the atomic dipoles, as is clear in the single atom laser [25], or in atomic beam lasers where only one atom is present in the cavity at any given time. A more relevant situation is the quantum synchronization that takes place in the context of superradiance [26,27]. It has recently been understood that, in contrast to lasers, steady-state superradiance can produce spectrally pure light [26][27][28][29][30][31][32][33] without stimulated emission. So far this has been demonstrated using a cavity mode as a communication channel that spatially selects an optical mode and enhances the coupling (through the cavity finesse). A more generic and relevant scenario, with great potential and applicability, is the emergence of spontaneous macroscopic quantum synchronization in radiating dipole arrays without a cavity but naturally coupled by the intrinsic anisotropic and long-range dipolar interactions. This is the situation considered in this paper.
Here we demonstrate that in the presence of an incoherent repumping source, dipole induced cooperative emission can dominate over spatial inhomogeneities and quantum fluctuations and lead to a resilient steadystate that exhibits macroscopic quantum phase coherence and intrinsic quantum correlations. An iconic example of a macroscopic coherent state is a Bose-Einstein condensate, achieved in ultra-cold gases at thermal equilibrium. In our case, however, the macroscopic order is reached in the steady-state of an interacting and driven-dissipative system. Moreover, the cooperative behavior can be detected by measuring the spectral purity of the emitted radiation. We note that in clear distinction to previous studies [8, 15-17, 19, 21], our proposal does not rely on an external coherent source or externally generated nonlinearities to seed the collective phase. In our model synchronization emerges as a spontaneously broken symmetry driven by incoherent processes in naturally coupled dipole arrays. As we show, and somewhat counterintuitively, an incoherent drive is sufficient to generate phase coherence in these systems.
Specifically, the systems we consider are dense arrays of frozen quantum dipoles modeled as quantized twolevel systems. By dense arrays of frozen dipoles we mean arrays separated by a distance much closer than the wavelength of the emitted photons and with motional degrees of freedom evolving at a much slower rate than their internal dynamics (figure 1). These conditions can be readily satisfied in a variety of quantum systems found in atomic, molecular and optical physics (e.g., Rydberg gases [34][35][36], alkali vapors [37], alkaline-Earth atoms [38], and polar molecules [39]), chemistry (e.g., J-aggregates of dye molecules [40][41][42]), and biology (e.g., light-harvesting complexes [43,44]). In cold vapors, one possible way to freeze the motion and tightly trap the particles is via an optical lattice potential (figure 1). In this case a sub-optical-wavelength transition must be used in order to reach the tight-packing regime [38,39].
To fully understand synchronization in the complex dipolar system, we analyze each of the ingredients that compete and affect synchronization in a step-by-step procedure: the interplay between repumping and collective emission, inhomogeneity in the coupling constants, quantum correlations, and the competition between elastic and inelastic interactions. The paper is organized as follows: in section 2 we introduce the system in consideration and the master equation we use to describe the dynamics. In section 3 we first provide a simple mean-field description and discuss connections to the classical Kuramoto model (KM)-the iconic model used to describe synchronization in nonlinear coupled oscillators. In section 4 we discuss the phase diagram for the quantum system assuming collective (all-to-all) coupling in the absence of inhomogeneity and elastic interactions, and compare it with the mean-field solution. For this exactly solvable case we are able to explicitly quantify the entanglement and correlations present in the steady state. In section 5 we study how inhomogeneity in the inelastic couplings affects synchronization and focus on the case of power-law decaying interactions. In section 6 we study the emergence of quantum synchronization in radiating dipoles taking the full long-range and anisotropic dipolar interactions into account. Specifically, we consider a one-dimensional geometry, where both the elastic and inelastic interactions can vary significantly accross the dipoles and can be easily adjusted. In Figure 1. Arrays of quantum dipoles spontaneously emit and absorb photons at rate Γ. The photons mediate dipolar interactions between dipoles separated by a distance r ab with both dissipative, f r ( ) ab , and elastic, g r ( ) ab , components. A repumping source at a rate W provides energy to maintain the oscillations and can be implemented using additional internal states that are not shown. (a) A possible implementation using cold atoms in an optical lattice. (b) The dipolar couplings g r ( ) ab and f r ( ) ab exhibit a complex angular distribution as a function of θ, the angle between the dipole orientation (determined by an external electromagnetic field) and r ab | |. The maximum value of f and g for fixed r ab | | is denoted as f max and g max . The cone illustrates the magic angle, section 7 we discuss experimental implementations of our model, and in section 8 we provide a conclusion and an outlook.

Dipole-dipole interaction and master equation
In this work we consider arrays of quantum dipoles with two accessible levels, which we denote as | ↓ 〉 and |↑〉. The interactions between two dipoles a and b are described by the functions g r ( ) ab and f r ( ) ab , which depend on the dipoles' separation, r ab | |, and the angle θ between the mean dipole moment and the vector joining the dipoles [38] (see figure 1(a)): where λ is the characteristic wavelength of the dipole-transition, and f (0) Γ = is the spontaneous photon emission rate from a single dipole. The function g r ( ) ab describes the elastic dipole-dipole interactions 4 , while f r ( ) ab gives rise to inelastic collective photon emission. These terms are similar to those that determine the radiation of classical electric dipoles, and the dependence on r ab | |reflects the propagation of photons from one atom to another. The terms 1 ab ζ ∝ account for retardation effects in the far-field regime and those 1 ab account for instantaneous propagation in the near-field. When 1 ab ζ ≪ , the elastic g interactions with a strong angular variation are dominant except close to the magic angle arccos(1 3 ) m θ = , at which they are greatly suppressed. In contrast, f r ( ) ab is almost isotropic in the near-field regime (see figure 1(b)). The spatially uniform behavior of f r ( ) ab at short distance is what gives rise to cooperative effects and superradiant emission [26]. Under generic conditions, however, superradiance is a transient effect that substantially limits the lifetime of dipole excitations. To compensate for the fast decay here we add an incoherent repumping driving term at a rate W. This term is needed to generate a synchronized steady state where longlasting coherence persists. An incoherent repumping drive is commonly used in laser systems to maintain population inversion. It can be implemented by coherently driving, at a rate Ω, the | ↓ 〉 state to an excited level that spontaneously decays, at a rate γ ≫ Ω, to the state |↑〉 [48]. Due to the fast depletion of the excited state, it can be adiabatically eliminated and thus the net process is just an incoherent transfer of population from | ↓ 〉 to |↑〉 at a rate W 2 γ = Ω (see appendix A for details). The evolution of N dipoles is modeled by a quantum master equation for the reduced density matrixρ of the dipoles [26]: The Hamiltonian Ĥ 0 generates the coherent evolution of the dipole array whereˆa z ( , , ) σ + − are the Pauli spin operators for dipole a, a δ denotes its oscillation frequency in the rotating frame defined by the mean frequency of the dipole ensemble, and ℏ is the reduced Planck constant. The Lindblad operator functionals, f W ,  , describe the inelastic photon emission and incoherent repumping processes, respectively. 4 For generic cases, the pathological divergence of g r ( )at r 0 = can be removed by introducing an additional term into the expression of g r ( ) ab , r ( ) k ab [45,46]. This gives rise to the so called Lorentz-Lorentz shift [47]. For the discrete arrays we consider in this work, there is always a cut-off distance and thus no divergence. For example, in the suggested implementation using Sr atoms that is discussed in section 7, the cut-off distance is determined by the optical lattice spacing.

Mean-field treatment and connection to the KM
To obtain a qualitative picture of how synchronization can happen among the dipoles, we first perform a meanfield treatment and show the close connection between our quantum model and the prototype models for classical synchronization. The mean-field approach assumes uncorrelated dipoles, i.e.,ˆâ × matrix in the pseudospin 1 2 basis { , } |↑〉 | ↓ 〉 . The components of the single-dipole density matrix,ˆa ρ , can be visualized as a Bloch vector S t t S t figure 2). The mean-field solution yields a system of coupled nonlinear differential equations for a , ρ σ σ′ . For each a N 1, 2 ,..., = the parameters evolve as The term proportional to f r ( ) ab in equation (7) that contains the sine function can be identified with a similar term in the KM [49]: where K, the coupling strength per oscillator, must be large enough and positive for synchronization to occur. The term proportional to g r ( ) ab that contains the cosine function appears in the Sakaguchi-Kuramoto model [50]-a more general but similar synchronization model to the KM. Compared to the basic KM, the situation here is more complex. This is due to the fact that in equation (7) and look for a solution in which Z is time-independent and synchronized oscillators possess a collective frequency ω , and thus a macroscopic phase t ω Φ = . These conditions lead to two equations for the order parameter Z and the collective frequency ω (see appendix B): The solution is shown in figure 2(a). The insets show the phase distribution in the steady-state for an array of oscillators initially prepared with a random distribution of phases for two different values of the repump rate W. For a slow repumping rate (bottom inset), the system remains unsynchronized. As the repumping rate is increased beyond a threshold value, the system enters a synchronized state, as can be seen by the appearance of phase locking and the resulting narrow phase spread (top inset). This can be explained by the fact that one necessary condition for synchronization in the KM is K 0 > , which translates to the requirement S 0 a z > on average in our model and thus the need to have sufficiently large repump rate. In the limit f eff ≫ Γ (e.g., for large N), maximum synchronization is achieved at W f 2 opt eff = , where the order parameter Z reaches a maximum value Z 1 8 max ≈ . For this optimal condition for synchronization the quantum dipoles are ordered with the same phase and radiate with atomic inversion S 1 4 a z ≈ . Note that the maximum order parameter is smaller than 1 2 even when fully synchronized because of this required finite value of the atomic inversion. One intriguing aspect is that repumping, which is the process that builds up synchronization, is itself an incoherent process. It is crucial that repumping does not preserve the norm of the collective Bloch vector, allowing it to extend or contract. For large W W opt > , Z decreases again reflecting a suppression of synchronization. In this limit the dipoles are repumped so fast that they are all driven to the |↑〉 state (S 1 2 a z → and S 0 a → ⊥ ) and phase coherence between them cannot build up.
The cases of a heterogeneous distribution of a δ ʼs can be treated at the mean-field level in a simple way. The results, summarized in the appendix B, are qualitatively similar. In general, the inclusion of a finite spread Δ in a δ decreases the value of the order parameter Z.
We note that for given f eff and W, synchronization is destroyed (that is, . For the case with finite g r ( ) 0 ab ≠ , equation (7) predicts dephasing caused by different g r ( ) 0 ab ≠ , and a global rotation that does not affect synchronization if all g r ( ) ab are identical (see appendix B). However, the mean-field ansatz is not appropriate for treating the case where elastic interactions are present since it neglects quantum entanglement between dipoles. The case with finite g r ( ) ab will be discussed in section 6 by exactly solving the many-body problem.

Quantum synchronization for the collective system
In the simplified case where 0 a δ = for all dipoles, g r ( ) 0 ab = for all pairs, and a constant collective decay rate Nf f r ( ) ab eff ≡ , it is possible to exactly solve equation (1), i.e., the full quantum dynamics, even for many particles, allowing us to benchmark the validity of the mean-field solution. This is due to the invariance of the master equation under individual dipole permutations that reduces the scaling of the Liouville space from exponential, 4 N , to polynomial, of order N 3 [51].

Phase diagram
Quantum fluctuations can lead to phase diffusion and to decay of single particle coherences in the steady-state (it is possible forˆ0 a σ 〈 〉 → + even in a synchronized state), so Z cannot be used as a measure of synchronization in a beyond mean-field treatment. However, phase locking in quantum mechanics can be quantified by the degree of spin-spin correlations Z Q , defined by ZˆQ  28,29]. The corresponding phase diagram, shown in figure 2(b), closely resembles the mean-field one.
To demonstrate that Z Q can be used to quantify the emergence of quantum synchronization, regardless of the inherent non-equilibrium and dissipative character of our system, we have also computed pairwise two-time correlation functions (see appendix B). The decay rate of these correlations encodes information about the spectral coherence of the emitted radiation. The range of W Γ values where the emitted light is maximally coherent agrees with the regime where the system is optimally synchronized according to Z Q . Moreover, we have also confirmed the moderate importance of higher order correlations in the synchronized steady-state by comparing the exact solution with a cumulant expansion calculation (which includes lowest order corrections to the mean-field result). We find the cumulant expansion agrees well with the exact solution (see appendix B). The only limit where there are important deviations is at very weak pumping W ≪ Γ where the system favors subradiant emission arising from strong atom-atom correlations (indicated by the purple region in figure 2

Quantum correlations and entanglement
The robust macroscopic quantum coherence exhibited by the synchronized state leads to the natural question of whether or not entanglement can be present in the steady-state even in this dissipative environment. Most previous studies that attempted to address this question have been limited to small systems [13,18,21,52] and focused on the entanglement between a pair of synchronized oscillators. Here, to determine the non-separability of the many-body steady-state, we compute the average of the quantum Fisher information (QFI) and use it as an entanglement witness [53,54] (see appendix D). Any N particle state with N N F N ( 2 ) 3¯(ˆ) 2 3 Q 2 ρ + ⩾ > is entangled (non-separable) and a quantum resource for phase estimation (see appendix C for details).
Due to dissipation, the density matrix of the system is reduced to a mixed state, as obtained from equation (1), which describes the dynamics of the system after a statistical average over many experimental trials. However, the evolution of the system for an individual experimental realization can be quite different. We consider a Gedanken experiment in which one monitors the system evolution and keeps a measurement record of the emitted photons. The evolution of the system is then conditioned on the measurement record [55,56]. This type of conditional evolution has been widely studied in quantum optics and utilized together with quantum feedback control in examples such as the optimal generation of spin squeezed states (see [57,58]). It should be emphasized that the conditional evolution based on the measurement record gives a quantum trajectory that should not be regarded simply as a numerical tool to allow the efficient assembly of ensemble averages. Each quantum trajectory is a potentially realizable physical outcome (even if hard to perform in practice) as allowed by the quantum dynamics of the open quantum system under consideration.
We calculate F (ˆ) Q c ρ (with the c inˆc ρ meaning conditional) for each conditional trajectory and in figure 2(d) we show its average over a sufficiently large set of trajectories at steady state (see appendix C). For this conditional evolution we observe entanglement in a parameter regime that correlates with Z 0 Q > (see figures 2(c) and (d)). On the other hand, if we discard the information present in the measurement record, by using the ensemble averaged ρ obtained from directly solving equation (1), and then computing F (ˆ) Q ρ , the QFI falls below the entanglement witness threshold (see figure 2(d)).
To differentiate quantum effects from classical ones, we further calculate the quantum discord  (see appendix D), which can be considered as a measure of quantum correlation more general than entanglement and more robust in a dissipative environment [59,60]. Separable states with nonzero  are intrinsically nonclassical, since local measurements performed on a subsystem inevitably disturb the whole system [59,61]. We measure classical correlations of the steady state by the difference between the mutual information  (see appendix D) and . We find the mixed steady-state contains nonzero quantum correlations in the synchronized regime figure 3(a). Moreover we observed that although both  and Z Q exhibit a similar dependence with pumping rate W, they do not exactly peak at the same value [15].
Although the existence of a nonzero  has been reported to exist in several quantum synchronization studies [18,61], we want to emphasize that those have been always limited to small systems. To our knowledge our calculations are the first to consider  in macroscopic samples. In figure 3(b) we show the dependence of  with system size. Our calculation shows that quantum correlations remain a significant fraction of  even in the thermodynamic limit.

Synchronization with finite-range interactions
Up to this point we have only considered all-to-all interactions; now we consider the effect of finite range interactions on synchronization. In the dipole array both f r ( ) ab and g r ( ) ab are nontrivial functions of r ab and contain terms decaying as a power-law with distance, r 1 ab ∝ | | α with 1, 2, 3 α = . Instead of dealing with all these terms together, to gain insight on how spatial inhomogeneities affect quantum synchronization, we first study a simpler case assuming a power-law cooperative decay f r r ( ) ab ab ∝ | | α − , with the exponent α as a variable parameter and set both g r ( ) 0 ab = and 0 a δ = . In the classical regime [62][63][64][65], analytical calculations and numerical simulations considering arrays of oscillators interacting via power law interactions on a one-dimensional lattice had identified 3 2 c α = as the critical value of the power law exponent below which long-range phase order is possible [64]. For c α α < , a transition to a state in which a finite fraction of the oscillators is entrained takes place for a sufficiently strong but finite coupling strength in the large system limit. Generalizations of these results to oscillators in D dimensions [64] have also identified three different regimes for synchronization: perfect phase ordering for D α ⩽ , entrainment with long-range phase order for D 3 2 α < and a crossover to exponential decay of correlations at D (3 1) 2 α = + . Reference [65] has also suggested that in the regime D α > global synchronization is absent but local synchronization persists for arbitrary weak coupling with a slowly decaying order parameter.
To quantify the effect of finite-range interactions on synchronization in the quantum regime we compute spin-spin correlations within linear clusters that contain d dipoles, Z ( )ˆQ , the local order parameter Z Q d is almost independent of α and d and the system behaves almost like the all-to-all system. For D D 2 α ≲ ≲ synchronization remains global and almost independent of d, but the order parameter slowly decreases with α. For D α ≳ , synchronization becomes local and correlations quickly decrease with cluster size. The magenta contour provides an indicative scale of the boundary between global and local synchronization. The white contour lines also provide information about the decrease of the order parameter with increasing α and d. We observe that, as in the classical case, D α ∼ roughly marks the transition between global and local synchronization, although a more quantitative comparison would require far larger systems. = .
An alternative way to characterize domain formation and the fact that it can persist even when there is a variation in the local detunings, 0 a δ ≠ , is to examine pairwise two-time correlation functions, , which can be related to the emission spectrum of the pair of atoms [55]. The oscillations in Z ( ) a b , τ encode information about the relative precession rate between different dipoles. By parameterizing Z ( ) , τ ν τ γ τ = − we can extract the relative precession frequency ab ν between dipoles a and b, where entrainment of dipoles a and b corresponds to 0 ab ν = . To explore the entrainment of dipole pairs in our system, we assign random detunings distributed uniformly in

Synchronization of dipoles with elastic interactions
We now treat the full problem of radiating quantum dipoles incorporating elastic interactions g r ( ) ab and the intricate competition of spatially-dependent and anisotropic couplings (both g r ( ) ab and f r ( ) ab have terms with power law dependence 1, 2, 3 α = on distance) (figure 1). We solve the full master equation without any approximation [55] for systems of up to twenty dipoles in a chain using the actual spatial dependence of both f r ( ) ab and g r ( ) ab , and set 0 a δ = . We observe a robust synchronized state that exists in a wide parameter space. As long as f Nf r ( ) ab eff ≡ is large enough, we find that synchronization takes place and is only weakly affected by substantial differences in g r ( ) ab , e.g., variations in the dipole array that modify g r ( ) ab by two orders of magnitude only decrease the order parameter by a factor of two (figure 5) in the steady state. This is in striking contrast to the situation in a system without dissipation, where the elastic interaction is known to generate entanglement between spins and to cause a decay of the order parameter during time evolution [66].
For the orientation arccos(1 3 ) m θ = , the order parameter reaches a significant fraction of Z Q max , indicating the emergence of macroscopic spontaneous synchronization of the radiating quantum dipole array (figure 5). To further emphasize the relevant role played by the inelastic term, in figure 5(b) we compare a solution of the master equation (equation (1)) for two cases: a system of coupled dipoles arranged in a 1D chain and oriented at the magic angle (symbols) and an array of identically coupled dipoles with the same f eff but experiencing only inelastic interactions [ g r ( ) 0 ab = , dashed lines]. The calculated order parameters agree well for the two different cases. The similar behavior demonstrates that in spite of the complex geometry of the dipolar interactions, the capability of the dipole system to synchronize can be characterized to great extent by the quantity f eff .

Experimental implementation
Our calculations above demonstrate the potential for synchronization in a dense array of dipoles. The flexible and precise control exhibited by ultracold atomic systems make them ideal platforms to experimentally investigate the synchronization phenomenon predicted here. Atomic systems operate with a large number of quantum oscillators and also allow for the tunability of the interaction parameters over a broad range.
One possible set-up to observe synchronization consists of arrays of ultracold Sr 87 atoms prepared in two electronic internal states that form the two-level system. The | ↓ 〉 could then correspond to the long-lived P 5s5p 3 0 state, with an intercombination line narrower than 10 3 − s 1 − . This is the state used to operate the most precise atomic clocks [67]. The |↑〉 could correspond to the D 5s4d 3 1 state with a natural linewith 290 10 3 Γ = × s 1 − . Both states can be trapped in an optical lattice at the magic wavelength a 0.2 = μm [38], that generates the same trapping potential for both states minimizing Stark shifts and inhomogeneities in the coupling constants. The dipole-dipole interactions are mediated by photons at the wavelength 2.6 λ = μm and thus, as shown in figure 5, the ratio a 0.08 λ < falls in the parameter regime where dipoles can be synchronized. By changing the angle between the laser beams used to form the lattice potential, the lattice spacing can be varied allowing tunability of the interaction strength between dipoles. The incoherent pumping can be realized by coherently transferring the P 5s5p 3 0 population to one or several appropriate intermediate states that decay rapidly to the D 5s4d 3 1 state [68]. An example would be the P 4d5p 3 1 state [69]. The polarization of the dipoles can be oriented in an arbitrary direction by an electromagnetic field. Although all the dipoles cannot be oriented at the magic angle in a 3D geometry, one may still suppress the elastic interactions by dynamical decoupling techniques adopted from NMR [70]. Those have been already demonstrated in ultracold polar molecule systems [39]. Another possibility is to use a spatial configuration of external fields that induces an 'averaging out' of the dominant elastic interactions [71]. Moreover, by slightly departing from the magic-wavelength condition, the dipoles can be subjected to onsite inhomogeneities that generate different detunings a δ . The phase synchronization can be probed by measuring Z Q , which experimentally can be directly obtained from the fluorescence intensity. As suggested in section 5, phase locking can also be extracted from two-point correlations which can be determined by analyzing the fluorescence spectrum [68]. An intriguing but also more speculative and less controllable realization of our quantum dipole model is the case of fluorescent organic molecules. A possible two-level configuration in those systems consists of a vibrational level of the ground electronic potential chosen as | ↓ 〉 and the lowest vibronic level of the first excited potential chosen as |↑〉. Incoherent pumping can be realized by driving an optical transition to a higher excited vibronic level φ | 〉 in the first excited potential, which decays on picosecond timescales to the state |↑〉 via nonradiative transitions [72]. Typical values of the fluorescence wavelength f λ and lifetime f τ for organic chromophores under a variety of environmental conditions put these systems in a regime of near optimal synchronization. For instance, pseudoisocyanine chloride (PIC) and merocyanine derivatives commonly used in organic light-emitting diodes [73][74][75][76] typically form low-dimensional molecular aggregates in liquid solution with a 0.5 2.0 ≈ − Å, and ratios a f λ on the order of 10 3 − . The typical fluoresence decay rate for these organic chromophores is 0.1 1 Γ ∼ − GHz [72]. In order to achieve W 1 Γ = and enter the synchronized phase, the required pumping laser intensity is I 1 10 W ∼ − kW cm −2 , which is lower than the theoretical lasing threshold intensities I 0.1 1 th ∼ − MW cm −2 of dye lasers [72]. Therefore, it should be feasible to achieve steady state synchronization of organic dipoles via incoherent optical driving.

Conclusion
We have demonstrated that a system of radiating quantum dipoles can be synchronized in the presence of repumping. Our analytic mean-field approach provides a direct analogy between synchronization of quantum dipoles and synchronization of classical phase oscillators. Using exact solutions of the master equation and a cumulant expansion approach, we determined the necessary conditions for synchronization, and the entanglement properties in the steady state of macroscopic ensembles under different measurement protocols. We also analyzed the effect of finite-range interactions in large arrays. To our knowledge those have been previously explored only in the classical regime. These calculations, although restricted to g r ( ) 0 ab = , emphasize how the interplay between incoherent repumping and cooperative dipole decay can give rise to robust quantum synchronization. For treating the general case of dense packed dipoles, we numerically solved the master equation exactly for up to twenty dipoles, and studied the effect of elastic interactions. The onedimensional geometry chosen allows us to tune the dipolar interactions accross a wide range of parameters by simply adjusting the spacing and polarization angle, and investigate the competing roles of the elastic and inelastic interactions. One-dimensional systems are amenable for theoretical investigation since they help us to minimize finite size effects. They can be implemented in current experiments using ultra-cold gases and are relevant for the suggested implementation using organic molecules. There are several organic species such as PIC molecules [74] that are known to form quasi-one-dimensional arrays both as single crystals and in liquid solution. These systems are typically modelled as 1D arrays of two-level dipoles with energetic disorder to match spectroscopic data [77]. For this configuration, we found synchronization in the presence of a large variation of elastic interactions. We expect that in more general geometries, similar degree of synchronization could be achieved provided the collective dissipative coupling f eff dominates over the elastic pairwise couplings, g r ( ) ab . These conditions might be easier to be satisfied in higher dimensional dense packed dipolar arrays given the anisotropic character of the elastic terms which fluctuate both in sign and magnitude across the array. In higher dimensions, geometric effect becomes more important. Nevertheless, whether or not spontaneous synchronization in dipole arrays can be observed under more generic high dimensional geometries remains an open question.
Our results show that the intrinsic macroscopic coherence of the superradiant steady state is inherently resilient to single particle decoherence, spatial inhomogeneities, and noisy environmental effects. This observation could have relevant application to the development of low-threshold organic lasers, highly efficient solar cells, materials with enhanced chemical reactivity, as well as ultra-precise quantum devices, where these effects are anticipated to play an important role. Moreover, since quantum synchronization is imprinted in the spectral purity of the emitted radiation [29], the generated light may potentially serve as a direct diagnostic tool of quantum coherences in generic systems beyond cold gases such as organic molecules.
Take the specific case of strontium atoms in section 7 for example, the states P 5s5p 3 0 and D 5s4d 3 1 , separated by 2600 nm, realize the | ↓ 〉 and | ↑ 〉 dipole levels. The intermediate state can be chosen as P 4d5p 3 1 , which can be coupled to P 5s5p 3 0 via external lasers and decays to D 5s4d 3 1 at a rate 3.4 10 7 × s −1 , much faster than the decay from D 5s4d 3 1 to P 5s5p 3 0 , 2.9 10 5 × s −1 [38,69], allowing an effective population transfer from P 5s5p 3 0 to D 5s4d 3 1 . For the fluorescent organic molecule implementation, we propose a two-level configuration consisting of a vibrational level of the ground electronic potential for | ↓ 〉 and the lowest vibronic level of the first excited potential for |↑〉. The typical lifetime for this configuration is ∼ns. In this case the Figure A1. Three level configuration for the realization of the incoherent pumping process. States 1 | 〉 and 2 | 〉 correspond to the | ↓ 〉 and | ↑ 〉 levels of the two-level dipole system considered in the main text. State 3 | 〉 is a short-lived excited state. Ω is the Rabi frequency of the coherent driving laser, γ is the decay rate from 3 | 〉 to 2 | 〉, and Γ is the decay rate of level 2 | 〉.
incoherent pumping can be realized by driving an optical transition to a higher excited vibronic level 3 φ | 〉 = | 〉 in the first excited potential, which decays on a timescale of ps to the state | ↑ 〉 via non-radiative transitions [72].
In the case when the repumping mechanism is achieved via radiative decay from 3 | 〉 → | ↑ 〉, a relevant question to consider is the effect on synchronization if cooperative decay processes need to be accounted for. To give a quantitative analysis of the role of collective incoherent pumping, here we consider the pumping described by with η describing the effect of the collective pumping, which is 0 for individual pumping, and 1 if all processes are collective. The mean-field order parameter is which is decreased for any 0 η > . As seen from the above expression and shown in figure A2 , although the collective effect in incoherent pumping is to reduce coherences, and disrupt the collective behavior induced by f eff , it is still possible to maintain a finite order parameter. Individual repumping is therefore the most favorable case but some degree of collectiveness is tolerable.
Practically, the collective pumping can be suppressed by choosing intermediate states appropriately. Consider the strontium atoms discussed in section 7 for example, if the state P 4d5p 3 1 is chosen as 3 | 〉, the wavelength of decay 3 | 〉 → | ↑ 〉 is 522 23 λ = nm, much smaller than the dipole transition wavelength 2600 λ = nm . As shown in figure 5(a), synchronization can happen with an inter-particle spacing a ∼0.2 λ = 500 nm, which satisfies a 1 23 λ ∼ , therefore the collective effect in the pumping process can be much reduced. On the other hand if the decay from 3 |↑〉 − | 〉 is non-radiative, as it is the case for organic molecular aggregates, even for arrays with average dipole-spacing of only tens of nanometers the off-diagonal terms can be neglected. There the decay is through vibrational relaxation or charge transfer, so there is no photon exchange between dipoles that can lead to collective incoherent repumping.

Appendix B. Mean-field approach
The mean-field ansatz,ˆâ a ρ ρ = ⊗ , reduces the dynamics to N 3 coupled nonlinear differential equations presented in the main text. In the most generic case we define local order parameters to take into account the effect of the inhomogeneous couplings: If the local order parameters vary slowly over the system size, and can be approximated to be the same for all dipoles one can define where the global order parameter Z is defined as  which can be evaluated in the N → ∞ limit as integrals when the detunnings a δ have a known distribution. We note that in this treatment, when g r ( ) 0 ab ≠ , the elastic couplings simply induce a global frequency shift g Q f (2 ) eff eff ω = that can be eliminated by moving to a rotating frame. However, as discussed in the main text the role of elastic interactions is not just to introduce a simple rotation. In fact the elastic terms tend to generate entanglement and correlations which can not be capture by a simple mean-field treatment. Therefore the applicability of the mean-field solution is restricted to regimes where elastic interactions are highly suppressed. (AB is the total system spanned by A and B together) , , withˆi ρ the reduced density matrix of the subsystem i. A value varying between 0 and 2 is obtained when A and B are pure states or maximally correlated respectively. The mutual information can be separated into a classical and a quantum part. The classical part is max{ } Here B A  is the von Neumann entropy of subsystem B conditioned on the measurement performed on A and max represents maximum value obtainable over all local measurements on A. The quantum part, known as the quantum discord, B A B A = −    , measures the amount of correlations that exceed the classical part and characterizes the 'quantumness' of the system [61]. A state with nonzero quantum discord behaves in a way intrinsically non-classical, since a local measurement performed on one of its subsystems can disturb the whole system. In order to calculate B A  , we consider a set of von Neumann  figure 2(e) we calculate the mutual information from  and the quantum discord from  using as subsystems A and B a pair of dipoles, a and b respectively. The cumulant expansion solution agrees with the full calculation except at W ≪ Γ where subradiant behavior dominates (purple region) [29].