Synchronization of interacting quantum dipoles

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–18], and large ensembles at the mean field level or by including quantum corrections perturbatively [19–21]. Highly symmetric situations with collective coupling mediated, for example, by a cavity mode [22–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–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 steady-state 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 two-level 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–36], alkali vapors [37], alkaline-Earth atoms [38], and polar molecules [39]), chemistry (e.g., J-aggregates of dye molecules [40–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].

Zoom In Zoom Out Reset image size

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 section 7 we discuss experimental implementations of our model, and in section 8 we provide a conclusion and an outlook.

In this work we consider arrays of quantum dipoles with two accessible levels, which we denote as $| \downarrow \rangle$ and $| \uparrow \rangle$. The interactions between two dipoles a and b are described by the functions $g({{\bf{r}}}_{{ab}})$ and $f({{\bf{r}}}_{{ab}})$, which depend on the dipoles' separation, $| {{\bf{r}}}_{{ab}}|$, and the angle θ between the mean dipole moment and the vector joining the dipoles [38] (see figure 1(a)):

$$\begin{eqnarray*}g({{\bf{r}}}_{{ab}}) & = & -\displaystyle \frac{3{\rm{\Gamma }}}{2}\{{\mathrm{sin}}^{2}\theta \displaystyle \frac{\mathrm{cos}{\zeta }_{{ab}}}{{\zeta }_{{ab}}}+(3{\mathrm{cos}}^{2}\theta -1)[\displaystyle \frac{\mathrm{cos}{\zeta }_{{ab}}}{{({\zeta }_{{ab}})}^{3}}+\displaystyle \frac{\mathrm{sin}{\zeta }_{{ab}}}{{({\zeta }_{{ab}})}^{2}}]\},\\ f({{\bf{r}}}_{{ab}}) & = & \displaystyle \frac{3{\rm{\Gamma }}}{2}\phantom{\rule[-0000]{0ex}{5ex}}\{{\mathrm{sin}}^{2}\theta \displaystyle \frac{\mathrm{sin}{\zeta }_{{ab}}}{{\zeta }_{{ab}}}+(3{\mathrm{cos}}^{2}\theta -1)[\displaystyle \frac{\mathrm{sin}{\zeta }_{{ab}}}{{({\zeta }_{{ab}})}^{3}}-\displaystyle \frac{\mathrm{cos}{\zeta }_{{ab}}}{{({\zeta }_{{ab}})}^{2}}]\}.\end{eqnarray*}$$

Here, ${\zeta }_{{ab}}=2\pi | {{\bf{r}}}_{{ab}}| /\lambda$, where λ is the characteristic wavelength of the dipole-transition, and ${\rm{\Gamma }}=f(0)$ is the spontaneous photon emission rate from a single dipole. The function $g({{\bf{r}}}_{{ab}})$ describes the elastic dipole–dipole interactions⁴ , while $f({{\bf{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 $| {{\bf{r}}}_{{ab}}|$ reflects the propagation of photons from one atom to another. The terms $\propto 1/{\zeta }_{{ab}}$ account for retardation effects in the far-field regime and those $\propto 1/{\zeta }_{{ab}}^{3}$ account for instantaneous propagation in the near-field. When ${\zeta }_{{ab}}\ll 1$, the elastic g interactions with a strong angular variation are dominant except close to the magic angle ${\theta }_{m}=\mathrm{arccos}(1/\sqrt{3})$, at which they are greatly suppressed. In contrast, $f({{\bf{r}}}_{{ab}})$ is almost isotropic in the near-field regime (see figure 1(b)).

The spatially uniform behavior of $f({{\bf{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 long-lasting 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 $| \downarrow \rangle$ state to an excited level that spontaneously decays, at a rate $\gamma \gg {\rm{\Omega }}$, to the state $| \uparrow \rangle$ [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 $| \downarrow \rangle$ to $| \uparrow \rangle$ at a rate $W={{\rm{\Omega }}}^{2}/\gamma$ (see appendix A for details).

The evolution of N dipoles is modeled by a quantum master equation for the reduced density matrix $\hat{\rho }$ of the dipoles [26]:

$$\begin{eqnarray}&&\displaystyle \frac{{\rm{d}}\hat{\rho }}{{\rm{d}}t}=-\displaystyle \frac{{\rm{i}}}{{\rm{\hslash }}}[{\hat{H}}_{0},\hat{\rho }]+{{\mathcal{L}}}_{f}[\hat{\rho }]+{{\mathcal{L}}}_{W}[\hat{\rho }],\end{eqnarray} \tag{ 1 }$$

$$\begin{eqnarray}&&{\hat{H}}_{0}=\displaystyle \frac{{\rm{\hslash }}}{2}\displaystyle \sum _{a=1}^{N}[-{\delta }_{a}{\hat{\sigma }}_{a}^{z}+\displaystyle \sum _{b=1,b\ne a}^{N}g({{\bf{r}}}_{{ab}}){\hat{\sigma }}_{a}^{+}{\hat{\sigma }}_{b}^{-}],\end{eqnarray} \tag{ 2 }$$

$$\begin{eqnarray}&&{{\mathcal{L}}}_{f}[\hat{\rho }]=\displaystyle \frac{1}{2}\displaystyle \sum _{a,b}f({{\bf{r}}}_{{ab}})(2{\hat{\sigma }}_{b}^{-}\hat{\rho }{\hat{\sigma }}_{a}^{+}-{\hat{\sigma }}_{a}^{+}{\hat{\sigma }}_{b}^{-}\hat{\rho }-\hat{\rho }{\hat{\sigma }}_{a}^{+}{\hat{\sigma }}_{b}^{-}),\end{eqnarray} \tag{ 3 }$$

$$\begin{eqnarray}&&{{\mathcal{L}}}_{W}[\hat{\rho }]=\displaystyle \frac{W}{2}\displaystyle \sum _{a}(2{\hat{\sigma }}_{a}^{+}\hat{\rho }{\hat{\sigma }}_{a}^{-}-{\hat{\sigma }}_{a}^{-}{\hat{\sigma }}_{a}^{+}\hat{\rho }-\hat{\rho }{\hat{\sigma }}_{a}^{-}{\hat{\sigma }}_{a}^{+}).\end{eqnarray} \tag{ 4 }$$

The Hamiltonian ${\hat{H}}_{0}$ generates the coherent evolution of the dipole array where ${\hat{\sigma }}_{a}^{(+,-,z)}$ are the Pauli spin operators for dipole a, ${\delta }_{a}$ denotes its oscillation frequency in the rotating frame defined by the mean frequency of the dipole ensemble, and ${\rm{\hslash }}$ is the reduced Planck constant. The Lindblad operator functionals, ${{\mathcal{L}}}_{f,W}$, describe the inelastic photon emission and incoherent repumping processes, respectively.

To obtain a qualitative picture of how synchronization can happen among the dipoles, we first perform a mean-field 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., $\hat{\rho }={\displaystyle \otimes }_{a}{\hat{\rho }}_{a}$, where each ${\hat{\rho }}_{a}={\displaystyle \sum }_{\sigma ,{\sigma }^{\prime }=\uparrow ,\downarrow }{\rho }_{a}^{\sigma ,{\sigma }^{\prime }}| \sigma \rangle \langle {\sigma }^{\prime }|$ is a $2\times 2$ matrix in the pseudospin $1/2$ basis $\{| \uparrow \rangle ,| \downarrow \rangle \}$. The components of the single-dipole density matrix, ${\hat{\rho }}_{a}$, can be visualized as a Bloch vector $\{{S}_{a}^{\perp }(t)\mathrm{cos}{\varphi }_{a}(t),{S}_{a}^{\perp }(t)\mathrm{sin}{\varphi }_{a}(t),{S}_{a}^{z}(t)\}$ = $(1/2)\{{\rho }_{a}^{\uparrow \downarrow }+{\rho }_{a}^{\downarrow \uparrow },-{\rm{i}}({\rho }_{a}^{\uparrow \downarrow }-{\rho }_{a}^{\downarrow \uparrow }),{\rho }_{a}^{\uparrow \uparrow }-{\rho }_{a}^{\downarrow \downarrow }\}$ (figure 2). The mean-field solution yields a system of coupled nonlinear differential equations for ${\rho }_{a}^{\sigma ,{\sigma }^{\prime }}$. For each $a=1,2,...,N$ the parameters evolve as

$$\begin{eqnarray}\displaystyle \frac{{\rm{d}}{S}_{a}^{z}(t)}{{\rm{d}}t} & = & {-S}_{a}^{\perp }(t)\displaystyle \sum _{b=1,b\ne a}^{N}{S}_{b}^{\perp }(t)[f({{\bf{r}}}_{{ab}})\mathrm{cos}[\delta {\varphi }_{{ba}}]-g({{\bf{r}}}_{{ab}})\mathrm{sin}[\delta {\varphi }_{{ba}}]]\\ & & -{\rm{\Gamma }}[\displaystyle \frac{1}{2}+{S}_{a}^{z}(t)]+W[\displaystyle \frac{1}{2}-{S}_{a}^{z}(t)],\end{eqnarray} \tag{ 5 }$$

$$\begin{eqnarray}\displaystyle \frac{{\rm{d}}{S}_{a}^{\perp }(t)}{{\rm{d}}t} & = & {-S}_{a}^{z}(t)\displaystyle \sum _{b\ne a}{S}_{b}^{\perp }(t)[g({{\bf{r}}}_{{ab}})\mathrm{sin}[\delta {\varphi }_{{ba}}]-f({{\bf{r}}}_{{ab}})\mathrm{cos}[\delta {\varphi }_{{ba}}]]\\ & & -\displaystyle \frac{{\rm{\Gamma }}+W}{2}{S}_{a}^{\perp }(t),\end{eqnarray} \tag{ 6 }$$

$$\begin{eqnarray}&&\displaystyle \frac{{\rm{d}}{\varphi }_{a}(t)}{{\rm{d}}t}={\delta }_{a}+\displaystyle \sum _{b=1,b\ne a}^{N}{S}_{a}^{z}\displaystyle \frac{{S}_{b}^{\perp }}{{S}_{a}^{\perp }}[g({{\bf{r}}}_{{ab}})\mathrm{cos}[\delta {\varphi }_{{ba}}]+f({{\bf{r}}}_{{ab}})\mathrm{sin}[\delta {\varphi }_{{ba}}]],\end{eqnarray} \tag{ 7 }$$

where $\delta {\varphi }_{{ba}}(t)={\varphi }_{b}(t)-{\varphi }_{a}(t)$. The term proportional to $f({{\bf{r}}}_{{ab}})$ in equation (7) that contains the sine function can be identified with a similar term in the KM [49]:

$$\begin{eqnarray}&&\displaystyle \frac{{\rm{d}}{\varphi }_{a}}{{\rm{d}}t}={\delta }_{a}+K\displaystyle \sum _{b=1}^{N}\mathrm{sin}[\delta {\varphi }_{{ba}}],\end{eqnarray} \tag{ 8 }$$

where K, the coupling strength per oscillator, must be large enough and positive for synchronization to occur. The term proportional to $g({{\bf{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) the coupling constants are nonuniform and effectively time-dependent, since ${S}_{a}^{\perp }(t)$ and ${S}_{a}^{z}(t)$ are dynamic variables.

Zoom In Zoom Out Reset image size

To investigate whether the mean-field model admits spontaneous synchronization we consider first the simplified case where ${\delta }_{a}=0$ for all dipoles, impose $g({{\bf{r}}}_{{ab}})=0$ for all pairs, and assume a constant collective decay rate ${Nf}({r}_{{ab}})\equiv {f}_{\mathrm{eff}}$. We define a global order parameter Z as $Z{{\rm{e}}}^{{\rm{i}}{\rm{\Phi }}}=\frac{1}{N}{\displaystyle \sum }_{a}{S}_{a}^{\perp }{{\rm{e}}}^{{\rm{i}}{\varphi }_{a}}$ and look for a solution in which Z is time-independent and synchronized oscillators possess a collective frequency $\bar{\omega }$, and thus a macroscopic phase ${\rm{\Phi }}=\bar{\omega }t$. These conditions lead to two equations for the order parameter Z and the collective frequency $\bar{\omega }$ (see appendix B):

$$\begin{eqnarray}&&\bar{\omega }=0,\end{eqnarray} \tag{ 9 }$$

$$\begin{eqnarray}&&Z=\displaystyle \frac{\sqrt{{f}_{\mathrm{eff}}(W-{\rm{\Gamma }})-{(W+{\rm{\Gamma }})}^{2}}}{\sqrt{2}{f}_{\mathrm{eff}}}.\end{eqnarray} \tag{ 10 }$$

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\gt 0$, which translates to the requirement ${S}_{a}^{z}\gt 0$ on average in our model and thus the need to have sufficiently large repump rate. In the limit ${f}_{\mathrm{eff}}\gg {\rm{\Gamma }}$ (e.g., for large N), maximum synchronization is achieved at ${W}_{\mathrm{opt}}={f}_{\mathrm{eff}}/2$, where the order parameter Z reaches a maximum value ${Z}_{\mathrm{max}}\approx \sqrt{1/8}$. For this optimal condition for synchronization the quantum dipoles are ordered with the same phase and radiate with atomic inversion ${S}_{a}^{z}\approx 1/4$. 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\gt {W}_{\mathrm{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 $| \uparrow \rangle$ state (${S}_{a}^{z}\to 1/2$ and ${S}_{a}^{\perp }\to 0$) and phase coherence between them cannot build up.

The cases of a heterogeneous distribution of ${\delta }_{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 ${\delta }_{a}$ decreases the value of the order parameter Z. For instance, if ${\delta }_{a}$ is sampled from a Lorentzian distribution $p({\delta }_{a})={\rm{\Delta }}/[\pi ({{\rm{\Delta }}}^{2}+{\delta }_{a}^{2})]$,

$$\begin{eqnarray}&&Z=\displaystyle \frac{\sqrt{{f}_{\mathrm{eff}}P-{Q}^{2}+2{{\rm{\Delta }}}^{2}-2{\rm{\Delta }}\sqrt{{{\rm{\Delta }}}^{2}+{f}_{\mathrm{eff}}P}}}{\sqrt{2}{f}_{\mathrm{eff}}},\end{eqnarray} \tag{ 11 }$$

where $Q={\rm{\Gamma }}+W$ and $P=W-{\rm{\Gamma }}$. Optimal synchronization is obtained at a smaller repumping rate, ${W}_{\mathrm{opt}}\approx {f}_{\mathrm{eff}}/2-{\rm{\Delta }}/\sqrt{2}$. We note that for given ${f}_{\mathrm{eff}}$ and W, synchronization is destroyed (that is, Z = 0) at ${{\rm{\Delta }}}_{c}=({Q}^{2}-{f}_{\mathrm{eff}}P)/(2Q)$.

For the case with finite $g({{\bf{r}}}_{{ab}})\ne 0$, equation (7) predicts dephasing caused by different $g({{\bf{r}}}_{{ab}})\ne 0$, and a global rotation that does not affect synchronization if all $g({{\bf{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({{\bf{r}}}_{{ab}})$ will be discussed in section 6 by exactly solving the many-body problem.

In the simplified case where ${\delta }_{a}=0$ for all dipoles, $g({{\bf{r}}}_{{ab}})=0$ for all pairs, and a constant collective decay rate ${Nf}({{\bf{r}}}_{{ab}})\equiv {f}_{\mathrm{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³ [51].

4.1. Phase diagram

4.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}^{2}+2N)/3\geqslant {\bar{F}}_{Q}(\hat{\rho })\gt 2\;N/3$ 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 ${\bar{F}}_{Q}({\hat{\rho }}_{c})$ (with the c in ${\hat{\rho }}_{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}_{Q}\gt 0$ (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}(\hat{\rho })$, 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 ${\mathcal{D}}$ (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 ${\mathcal{D}}$ 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 ${\mathcal{I}}$ (see appendix D) and ${\mathcal{D}}$. We find the mixed steady-state contains nonzero quantum correlations in the synchronized regime figure 3(a). Moreover we observed that although both ${\mathcal{D}}$ and Z_Q exhibit a similar dependence with pumping rate W, they do not exactly peak at the same value [15].

Zoom In Zoom Out Reset image size

Although the existence of a nonzero ${\mathcal{D}}$ 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 ${\mathcal{D}}$ in macroscopic samples. In figure 3(b) we show the dependence of ${\mathcal{D}}$ with system size. Our calculation shows that quantum correlations remain a significant fraction of ${\mathcal{I}}$ even in the thermodynamic limit.

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({{\bf{r}}}_{{ab}})$ and $g({{\bf{r}}}_{{ab}})$ are nontrivial functions of ${{\bf{r}}}_{{ab}}$ and contain terms decaying as a power-law with distance, $\propto 1/{| {{\bf{r}}}_{{ab}}| }^{\alpha }$ with $\alpha =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({{\bf{r}}}_{{ab}})\propto {| {{\bf{r}}}_{{ab}}| }^{-\alpha }$, with the exponent α as a variable parameter and set both $g({{\bf{r}}}_{{ab}})=0$ and ${\delta }_{a}=0$.

In the classical regime [62–65], analytical calculations and numerical simulations considering arrays of oscillators interacting via power law interactions on a one-dimensional lattice had identified ${\alpha }_{c}=3/2$ as the critical value of the power law exponent below which long-range phase order is possible [64]. For $\alpha \lt {\alpha }_{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 $\alpha \leqslant D$, entrainment with long-range phase order for $\alpha \lt 3D/2$ and a crossover to exponential decay of correlations at $\alpha =(3D+1)/2$. Reference [65] has also suggested that in the regime $\alpha \gt 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}^{d})}^{2}\equiv \langle \bar{{\hat{\sigma }}_{a}^{+}{\hat{\sigma }}_{b}^{-}}{\rangle }_{d}$, using a cumulant expansion method as described in appendix C. Here the bar followed by a subscript d indicates an average over the pairs of different dipoles a and b contained in a linear cluster of size d. The linear clusters start at the central spin as shown in figure 4. We have confirmed that the cumulant expansion method reproduces well the correlation functions by performing direct comparisons with the exact solution (see appendix C). Figure 4(a) shows the behavior of ${Z}_{Q}^{d}$ as a function of cluster size d and power-law decay exponent α in arrays of dimension D = 1 and 2. Clear global synchronization with an order parameter ${Z}_{Q}^{d}$ independent of d is observed for $\alpha \lesssim D$. For $0\lt \alpha \lesssim D/2$, 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/2\lesssim \alpha \lesssim D$ synchronization remains global and almost independent of d, but the order parameter slowly decreases with α. For $\alpha \gtrsim 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, $\alpha \sim D$ roughly marks the transition between global and local synchronization, although a more quantitative comparison would require far larger systems.

Zoom In Zoom Out Reset image size

An alternative way to characterize domain formation and the fact that it can persist even when there is a variation in the local detunings, ${\delta }_{a}\ne 0$, is to examine pairwise two-time correlation functions, ${Z}_{a,b}(\tau )\equiv \underset{t\to \infty }{\mathrm{lim}}\langle ({\hat{\sigma }}_{a}^{+}(t+\tau )$ $+{\hat{\sigma }}_{b}^{+}(t+\tau ))({\hat{\sigma }}_{a}^{-}(t)+{\hat{\sigma }}_{b}^{-}(t))\rangle$, which can be related to the emission spectrum of the pair of atoms [55]. The oscillations in ${Z}_{a,b}(\tau )$ encode information about the relative precession rate between different dipoles. By parameterizing ${Z}_{a,b}(\tau )$ as ${Z}_{a,b}(\tau )=A\mathrm{cos}({\nu }_{{ab}}\tau )\mathrm{exp}(-\gamma \tau )$ we can extract the relative precession frequency ${\nu }_{{ab}}$ between dipoles a and b, where entrainment of dipoles a and b corresponds to ${\nu }_{{ab}}=0$. To explore the entrainment of dipole pairs in our system, we assign random detunings distributed uniformly in $[-{\rm{\Gamma }}/2,{\rm{\Gamma }}/2]$ to a linear chain of N = 200 dipoles and calculate ${\nu }_{{ab}}$ for $b=1,2,...,100$ with a = 101 corresponding to the central dipole in the chain. Synchronization regimes similar to those shown in figure 4(a) are observed for this D = 1 system, which we illustrate in figure 4(b) by plotting a histogram (top panel) and the distribution of frequencies ν (bottom panel) for three values of α. For global coupling, $\alpha =0$ (dark blue symbols), all the dipoles become entrained with each other ($\nu =0$), indicating complete synchronization; for $\alpha =0.65$ (red symbols) dipoles split into entrained ($\nu =0$) and drifting ($\nu \ne 0$) groups. While not all dipoles are entrained, the entrained dipoles are distributed along the whole array, and thus synchronization is still global; and for $\alpha =2$ (light blue symbols) the majority of dipoles are not entrained. These observations of relative precession frequencies between pairs of oscillators are consistent with the regimes obtained from the order parameter plotted in figure 4(a).

We now treat the full problem of radiating quantum dipoles incorporating elastic interactions $g({{\bf{r}}}_{{ab}})$ and the intricate competition of spatially-dependent and anisotropic couplings (both $g({{\bf{r}}}_{{ab}})$ and $f({{\bf{r}}}_{{ab}})$ have terms with power law dependence $\alpha =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({{\bf{r}}}_{{ab}})$ and $g({r}_{{ab}})$, and set ${\delta }_{a}=0$. We observe a robust synchronized state that exists in a wide parameter space. As long as ${f}_{\mathrm{eff}}\equiv N\bar{f({r}_{{ab}})}$ is large enough, we find that synchronization takes place and is only weakly affected by substantial differences in $g({{\bf{r}}}_{{ab}})$, e.g., variations in the dipole array that modify $g({{\bf{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].

Zoom In Zoom Out Reset image size

For the orientation ${\theta }_{m}=\mathrm{arccos}(1/\sqrt{3})$, the order parameter reaches a significant fraction of ${Z}_{Q}^{\mathrm{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}_{\mathrm{eff}}$ but experiencing only inelastic interactions [$g({{\bf{r}}}_{{ab}})=0$, 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}_{\mathrm{eff}}$.

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 ${}^{87}\mathrm{Sr}$ atoms prepared in two electronic internal states that form the two-level system. The $| \downarrow \rangle$ could then correspond to the long-lived $5{\rm{s}}5{\rm{p}}{}^{3}{P}_{0}$ state, with an intercombination line narrower than ${10}^{-3}$ ${{\rm{s}}}^{-1}$ . This is the state used to operate the most precise atomic clocks [67]. The $| \uparrow \rangle$ could correspond to the $5{\rm{s}}4{\rm{d}}{}^{3}{D}_{1}$ state with a natural linewith ${\rm{\Gamma }}=290\times {10}^{3}$ ${{\rm{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 $\lambda =2.6$ μm and thus, as shown in figure 5, the ratio ${\rm{a}}/\lambda \lt 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 $5{\rm{s}}5{\rm{p}}{}^{3}{P}_{0}$ population to one or several appropriate intermediate states that decay rapidly to the $5{\rm{s}}4{\rm{d}}{}^{3}{D}_{1}$ state [68]. An example would be the $4{\rm{d}}5{\rm{p}}{}^{3}{P}_{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 ${\delta }_{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 $| \downarrow \rangle$ and the lowest vibronic level of the first excited potential chosen as $| \uparrow \rangle$. Incoherent pumping can be realized by driving an optical transition to a higher excited vibronic level $| \varphi \rangle$ in the first excited potential, which decays on picosecond timescales to the state $| \uparrow \rangle$ via non-radiative transitions [72]. Typical values of the fluorescence wavelength ${\lambda }_{{\rm{f}}}$ and lifetime ${\tau }_{{\rm{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–76] typically form low-dimensional molecular aggregates in liquid solution with $a\approx 0.5-2.0\ \mathring{\rm A}$, and ratios ${\rm{a}}/{\lambda }_{{\rm{f}}}$ on the order of ${10}^{-3}$. The typical fluoresence decay rate for these organic chromophores is ${\rm{\Gamma }}\sim 0.1-1$ GHz [72]. In order to achieve $W/{\rm{\Gamma }}=1$ and enter the synchronized phase, the required pumping laser intensity is ${I}_{W}\sim 1-10$ kW cm⁻², which is lower than the theoretical lasing threshold intensities ${I}_{\mathrm{th}}\sim 0.1-1$ MW cm⁻² of dye lasers [72]. Therefore, it should be feasible to achieve steady state synchronization of organic dipoles via incoherent optical driving.

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({{\bf{r}}}_{{ab}})=0$, 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 one-dimensional 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}_{\mathrm{eff}}$ dominates over the elastic pairwise couplings, $g({{\bf{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.

The authors wish to acknowledge useful discussions with David Nesbitt, James K Thompson, Alexey V Gorshkov, Emanuel Knill, Kaden R A Hazzard, Michael Foss-Feig, Zhe-Xuan Gong, Michael L Wall, Dominic Meiser, Xibo Zhang and Timur V Tscherbul. This work was supported by NIST, the NSF (PIF-1211914 and PFC-1125844), AFOSR, AFOSR-MURI, ARO individual investigator awards, and DARPA QuASAR. Computations utilized the Janus supercomputer, supported by NSF (award number CNS-0821794), NCAR, and CU Boulder/Denver.

The simplest realization of the repumping mechanism can be explained using a three-level system, as depicted in figure A1 . Practically, the repumping process could involve a more complicated setup via multiple intermediate levels, but the three-level model is sufficient to capture the relevant physics. The states $| 1\rangle$ and $| 2\rangle$ correspond to the $| \downarrow \rangle$ and $| \uparrow \rangle$ states of our two-level dipole system, and level $| 3\rangle$ is an auxiliary state. The $| 1\rangle \leftrightarrow | 3\rangle$ transition is coupled to a monochromatic light field ${{\bf{E}}}_{\mathrm{ex}}$ with frequency ${\omega }_{L}$. The coupling strength is determined by the Rabi frequency ${\rm{\Omega }}={{\bf{d}}}_{13}\cdot {{\bf{E}}}_{\mathrm{ex}}$, with ${{\bf{d}}}_{13}$ the dipole moment of transition $| 1\rangle \leftrightarrow | 3\rangle$. For the appropriate light polarization and frequency, the electromagnetic field ${{\bf{E}}}_{\mathrm{ex}}$ only drives the transition $| 1\rangle \leftrightarrow | 3\rangle$. The state $| 3\rangle$ is however chosen such that it spontaneously decays to level $| 2\rangle$ at a rate $\gamma \gg {\rm{\Omega }}\gg {\rm{\Gamma }}$. For simplicity, but without loss of generality, we assume the driving field is on resonance with the transition $| 1\rangle \leftrightarrow | 3\rangle$, thus the atom-field interaction is given by the Hamiltonian term ${\hat{V}}_{{\rm{\Omega }}}=({\rm{\Omega }}/2)({\hat{\sigma }}_{31}+{\hat{\sigma }}_{13})$, where we have used the projectors ${\hat{\sigma }}_{{ij}}=| i\rangle \langle j|$. The dynamics of such a three-level system can be described by the following master equation:

$$\begin{eqnarray}&&\displaystyle \frac{{\rm{d}}\hat{\rho }}{{\rm{d}}t}=-{\rm{i}}[{\hat{V}}_{{\rm{\Omega }}},\hat{\rho }]+\displaystyle \frac{\gamma }{2}(2{\hat{\sigma }}_{23}\hat{\rho }{\hat{\sigma }}_{32}-\{{\hat{\sigma }}_{32}{\hat{\sigma }}_{23},\hat{\rho }\}),\end{eqnarray} \tag{ A.1 }$$

where $\hat{\rho }$ represents the reduced state in a frame rotating at the driving frequency ${\omega }_{L}$. Note that for the present discussion we have not included the decay rate Γ since it is assumed to be much slower than γ and is accounted for in equation (3). From the resulting equations of motion, we can adiabatically eliminate the variables associated with the state $| 3\rangle$. The reason is that in the steady state the following stationary conditions are satisfied: $\frac{{\rm{d}}{\rho }_{32}}{{\rm{d}}t}\to 0$ and $\frac{{\rm{d}}{\rho }_{13}}{{\rm{d}}t}\to 0$, ${\rho }_{{ij}}\equiv \langle j| \hat{\rho }| i\rangle$. In addition, ${\rho }_{33}\sim 0$, due to the fact that $\gamma \gg {\rm{\Omega }}$, and thus in steady state the population of $| 3\rangle$ remains negligible. After the adiabatic elimination, the equations of motion for the states $| 1\rangle$ and $| 2\rangle$ are:

$$\begin{eqnarray}&&\displaystyle \frac{{\rm{d}}{\rho }_{11}}{{\rm{d}}t}=-W{\rho }_{11},\end{eqnarray} \tag{ A.2 }$$

$$\begin{eqnarray}&&\displaystyle \frac{{\rm{d}}{\rho }_{22}}{{\rm{d}}t}=\;+W{\rho }_{22},\end{eqnarray} \tag{ A.3 }$$

$$\begin{eqnarray}&&\displaystyle \frac{{\rm{d}}{\rho }_{12}}{{\rm{d}}t}=-\displaystyle \frac{W}{2}{\rho }_{12}.\end{eqnarray} \tag{ A.4 }$$

where $W={{\rm{\Omega }}}^{2}/\gamma$ is the incoherent repumping rate. These equations of motion define the effective Linblad operator used to describe the repumping process, equation (4):

$$\begin{eqnarray}&&{{\mathcal{L}}}_{W}[\hat{\rho }]=\displaystyle \frac{W}{2}[2{\hat{\sigma }}^{+}\hat{\rho }{\hat{\sigma }}^{-}-\{{\hat{\sigma }}^{-}{\hat{\sigma }}^{+},\hat{\rho }\}],\end{eqnarray} \tag{ A.5 }$$

with ${\hat{\sigma }}^{+}={\hat{\sigma }}_{21}$.

Zoom In Zoom Out Reset image size

Take the specific case of strontium atoms in section 7 for example, the states $5{\rm{s}}5{\rm{p}}{}^{3}{P}_{0}$ and $5{\rm{s}}4{\rm{d}}{}^{3}{D}_{1}$, separated by 2600 nm, realize the $| \downarrow \rangle$ and $| \uparrow \rangle$ dipole levels. The intermediate state can be chosen as $4{\rm{d}}5{\rm{p}}{}^{3}{P}_{1}$, which can be coupled to $5{\rm{s}}5{\rm{p}}{}^{3}{P}_{0}$ via external lasers and decays to $5{\rm{s}}4{\rm{d}}{}^{3}{D}_{1}$ at a rate $3.4\times {10}^{7}$ s⁻¹, much faster than the decay from $5{\rm{s}}4{\rm{d}}{}^{3}{D}_{1}$ to $5{\rm{s}}5{\rm{p}}{}^{3}{P}_{0}$, $2.9\times {10}^{5}$ s⁻¹ [38, 69], allowing an effective population transfer from $5{\rm{s}}5{\rm{p}}{}^{3}{P}_{0}$ to $5{\rm{s}}4{\rm{d}}{}^{3}{D}_{1}$. For the fluorescent organic molecule implementation, we propose a two-level configuration consisting of a vibrational level of the ground electronic potential for $| \downarrow \rangle$ and the lowest vibronic level of the first excited potential for $| \uparrow \rangle$. The typical lifetime for this configuration is ∼ns. In this case the incoherent pumping can be realized by driving an optical transition to a higher excited vibronic level $| 3\rangle =| \varphi \rangle$ in the first excited potential, which decays on a timescale of ps to the state $| \uparrow \rangle$ via non-radiative transitions [72].

In the case when the repumping mechanism is achieved via radiative decay from $| 3\rangle \to | \uparrow \rangle$, 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

$$\begin{eqnarray}{{\mathcal{L}}}_{W}[\hat{\rho }] & = & \displaystyle \frac{W}{2}\displaystyle \sum _{a}(2{\hat{\sigma }}_{a}^{+}\hat{\rho }{\hat{\sigma }}_{a}^{-}-{\hat{\sigma }}_{a}^{-}{\hat{\sigma }}_{a}^{+}\hat{\rho }-\hat{\rho }{\hat{\sigma }}_{a}^{-}{\hat{\sigma }}_{a}^{+})\\ & & +\displaystyle \frac{\eta W}{2}\displaystyle \sum _{a\ne b}(2{\hat{\sigma }}_{a}^{+}\hat{\rho }{\hat{\sigma }}_{b}^{-}-{\hat{\sigma }}_{b}^{-}{\hat{\sigma }}_{a}^{+}\hat{\rho }-\hat{\rho }{\hat{\sigma }}_{b}^{-}{\hat{\sigma }}_{a}^{+}),\end{eqnarray} \tag{ A.6 }$$

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

$$\begin{eqnarray}&&Z=\displaystyle \frac{\sqrt{({f}_{\mathrm{eff}}-\eta {NW})(W-{\rm{\Gamma }})-{(W+{\rm{\Gamma }})}^{2}}}{\sqrt{2}({f}_{\mathrm{eff}}-\eta {NW})},\end{eqnarray} \tag{ A.7 }$$

which is decreased for any $\eta \gt 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}_{\mathrm{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.

Zoom In Zoom Out Reset image size

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 $4{\rm{d}}5{\rm{p}}{}^{3}{P}_{1}$ is chosen as $| 3\rangle$, the wavelength of decay $| 3\rangle \to | \uparrow \rangle$ is ${\lambda }_{23}=522$ nm, much smaller than the dipole transition wavelength $\lambda =2600$ nm . As shown in figure 5(a), synchronization can happen with an inter-particle spacing a ∼0.2 λ = 500 nm, which satisfies ${\rm{a}}/{\lambda }_{23}\sim 1$, therefore the collective effect in the pumping process can be much reduced. On the other hand if the decay from $| \uparrow \rangle -| 3\rangle$ 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.

The mean-field ansatz, $\hat{\rho }={\displaystyle \otimes }_{a}{\hat{\rho }}_{a}$, reduces the dynamics to $3N$ 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:

$$\begin{eqnarray}&&{X}_{a}{{\rm{e}}}^{{\rm{i}}{\rm{\Phi }}}=\displaystyle \sum _{b\ne a}f({{\bf{r}}}_{{ab}}){S}_{b}^{\perp }{{\rm{e}}}^{{\rm{i}}{\varphi }_{b}},{Y}_{a}{{\rm{e}}}^{{\rm{i}}{\rm{\Phi }}}=\displaystyle \sum _{b\ne a}g({{\bf{r}}}_{{ab}}){S}_{b}^{\perp }{{\rm{e}}}^{{\rm{i}}{\varphi }_{b}}.\end{eqnarray} \tag{ B.1 }$$

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 ${X}_{a}\approx {f}_{\mathrm{eff}}Z,{Y}_{a}\approx {g}_{\mathrm{eff}}Z$, where the global order parameter Z is defined as $Z{{\rm{e}}}^{{\rm{i}}{\rm{\Phi }}}=\frac{1}{N}{\displaystyle \sum }_{a}{S}_{a}^{\perp }{{\rm{e}}}^{{\rm{i}}{\varphi }_{a}}$ and the effective couplings are given by ${f}_{\mathrm{eff}}={\displaystyle \sum }_{a}{\displaystyle \sum }_{b\ne a}f({{\bf{r}}}_{{ab}})/(N-1)$ and ${g}_{\mathrm{eff}}={\displaystyle \sum }_{a}{\displaystyle \sum }_{b\ne a}g({{\bf{r}}}_{{ab}})/(N-1)$.

The steady-state solution $\dot{Z}=0$, ${\rm{\Phi }}=\bar{\omega }t$ leads to two self-consistent equations for the order parameter Z and the collective frequency $\bar{\omega }$

$$\begin{eqnarray}&&Z=\displaystyle \sum _{a}^{N}\displaystyle \frac{{ZP}[{f}_{\mathrm{eff}}Q+2{g}_{\mathrm{eff}}({\delta }_{a}+\bar{\omega })]}{{NQ}[4{({\delta }_{a}+\bar{\omega })}^{2}+(2{f}_{\mathrm{eff}}^{2}{Z}^{2}+2{g}_{\mathrm{eff}}^{2}{Z}^{2}+{Q}^{2}]},\end{eqnarray} \tag{ B.2 }$$

$$\begin{eqnarray}&&0=\displaystyle \sum _{a}^{N}\displaystyle \frac{{ZP}[{g}_{\mathrm{eff}}Q-2{f}_{\mathrm{eff}}({\delta }_{a}+\bar{\omega })]}{{NQ}[4{({\delta }_{a}+\bar{\omega })}^{2}+(2{f}_{\mathrm{eff}}^{2}{Z}^{2}+2{g}_{\mathrm{eff}}^{2}{Z}^{2}+{Q}^{2}]},\end{eqnarray} \tag{ B.3 }$$

which can be evaluated in the $N\to \infty$ limit as integrals when the detunnings ${\delta }_{a}$ have a known distribution.

We note that in this treatment, when $g({{\bf{r}}}_{{ab}})\ne 0$, the elastic couplings simply induce a global frequency shift $\bar{\omega }={g}_{\mathrm{eff}}Q/(2{f}_{\mathrm{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.

The cumulant expansion method is a useful theoretical tool for including correlation effects beyond the mean-field approximation [29, 78, 79]. We keep two-point correlations such as $\langle {\hat{\sigma }}_{a}^{+,-,z}{\hat{\sigma }}_{b}^{+,-,z}\rangle$, but factorize three-point correlations and higher [80]:

$$\begin{eqnarray}\langle {\hat{\sigma }}_{a}^{\alpha }{\hat{\sigma }}_{b}^{\beta }{\hat{\sigma }}_{c}^{\gamma }\rangle & = & \langle {\hat{\sigma }}_{a}^{\alpha }{\hat{\sigma }}_{b}^{\beta }\rangle \langle {\hat{\sigma }}_{c}^{\gamma }\rangle +\langle {\hat{\sigma }}_{a}^{\alpha }\rangle \langle {\hat{\sigma }}_{b}^{\beta }{\hat{\sigma }}_{c}^{\gamma }\rangle +\langle {\hat{\sigma }}_{a}^{\alpha }{\hat{\sigma }}_{c}^{\gamma }\rangle \langle {\hat{\sigma }}_{b}^{\beta }\rangle \\ & & -2\langle {\hat{\sigma }}_{a}^{\alpha }\rangle \langle {\hat{\sigma }}_{b}^{\beta }\rangle \langle {\hat{\sigma }}_{c}^{\gamma }\rangle .\end{eqnarray} \tag{ C.1 }$$

This factorization closes the set of dynamical equations of motion for all single particle observables $\langle {\hat{\sigma }}_{a}^{+,-,z}\rangle$ and equal time two-point correlations. Two-time correlation functions can be computed by solving [55]:

$$\begin{eqnarray*}\displaystyle \frac{{\rm{d}}\langle {\hat{\sigma }}_{a}^{+}(t+\tau ){\hat{\sigma }}_{b}^{-}(t)\rangle }{{\rm{d}}\tau } & = & -[{\rm{i}}{\delta }_{a}+\displaystyle \frac{{\rm{\Gamma }}+W}{2}]\langle {\hat{\sigma }}_{a}^{+}(t+\tau ){\hat{\sigma }}_{b}^{-}(t)\rangle \\ & & +\displaystyle \frac{1}{2}{f}_{{ab}}\langle {\hat{\sigma }}_{a}^{z}(t)\rangle \langle {\hat{\sigma }}_{b}^{+}(t+\tau ){\hat{\sigma }}_{b}^{-}(t)\rangle \\ & & +\displaystyle \frac{1}{2}\displaystyle \sum _{j\ne a,b}{f}_{{aj}}\langle {\hat{\sigma }}_{a}^{z}(t)\rangle \langle {\hat{\sigma }}_{j}^{+}(t+\tau ){\hat{\sigma }}_{b}^{-}(t)\rangle ,\end{eqnarray*}$$

where we have introduced the approximation $\langle {\hat{\sigma }}_{a}^{z}(t+\tau ){\hat{\sigma }}_{j}^{+}(t+\tau ){\hat{\sigma }}_{b}^{-}(t)\rangle \approx \langle {\hat{\sigma }}_{a}^{z}(t)\rangle \langle {\hat{\sigma }}_{j}^{+}(t+\tau ){\hat{\sigma }}_{b}^{-}(t)\rangle$. Comparisons with exact numerical solutions show that the cumulant expansion captures well the steady-state behavior for inhomogeneous couplings $f({{\bf{r}}}_{{ab}})$, provided the elastic couplings $g({{\bf{r}}}_{{ab}})$ are sufficiently small. In figure C1 we compute the pair-wise two-time correlation function, ${Z}_{a,b}(\tau )\equiv {\mathrm{lim}}_{t\to \infty }\langle ({\hat{\sigma }}_{a}^{+}(t+\tau )$ $+{\hat{\sigma }}_{b}^{+}(t+\tau ))({\hat{\sigma }}_{a}^{-}(t)+{\hat{\sigma }}_{b}^{-}(t))\rangle$, using both the cumulant expansion and the exact solution. The decay rate of these correlations, ${Z}_{a,b}(\tau )=A{{\rm{e}}}^{-\gamma \tau }$, encodes information about the spectral coherence of the emitted radiation (note that here $\nu =0$). The result shows that ${\rm{\Gamma }}/\gamma$ exhibits the same dependence on $W/{\rm{\Gamma }}$ as Z_Q.

Zoom In Zoom Out Reset image size

An individual experimental realization can be considered as a single trajectory, whose evolution can be quite different from the ensemble averaged solution of the master equation. Tracking the evolution of an individual trajectory is equivalent to performing continuous measurements that collect the record of the emitted photons, for example homodyne measurements. The conditional evolution of the system subject to continuous measurements can be modeled by the method of quantum state diffusion [55, 81]. For a single run the state of the system remains pure, ${\hat{\rho }}_{c}=| \psi \rangle \langle \psi |$, but the average over many trials reduces the system into a mixed state and recovers the density matrix obtained from the master equation.

To probe the entanglement of the dipoles, in figure 2 we calculate the average quantum Fisher information for each individual trajectory [53, 54]

$$\begin{eqnarray*}&&{\bar{F}}_{Q}({\hat{\rho }}_{c})=\displaystyle \frac{1}{3}({F}_{Q}({\hat{\rho }}_{c};{\hat{J}}_{x})+{F}_{Q}({\hat{\rho }}_{c};{\hat{J}}_{y})+{F}_{Q}({\hat{\rho }}_{c};{\hat{J}}_{z})),\end{eqnarray*}$$

where ${F}_{Q}(\hat{{\rho }_{c}};\hat{{\mathcal{H}}})=4(| \psi \rangle {\hat{{\mathcal{H}}}}^{2}\langle \psi | -| \psi \rangle \hat{{\mathcal{H}}}\langle {\psi | }^{2})$, and ${\hat{J}}_{x,y,z}$ are collective angular momentum operators.

States with zero entanglement can still be nonclassical. Two systems are correlated if they share information with each other. The total amount of correlation can be quantified by the quantum mutual information ${\mathcal{I}}={{S}}_{A}+{{S}}_{B}-{{S}}_{{AB}}$, where ${{\mathcal{S}}}_{i}$ is the von Neumann entropy of the subsystem $i\in \{A,B,{AB}\}$ (AB is the total system spanned by A and B together) , ${{\mathcal{S}}}_{i}=-\mathrm{Tr}[{\hat{\rho }}_{i}{}_{}{\mathrm{log}}_{2}{\hat{\rho }}_{i}]$, with ${\hat{\rho }}_{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 ${{\mathcal{J}}}_{B| A}=\mathrm{max}\{{{\mathcal{S}}}_{B}-{{\mathcal{S}}}_{B| A}\}$. Here ${{\mathcal{S}}}_{B| A}$ is the von Neumann entropy of subsystem B conditioned on the measurement performed on A and $\mathrm{max}$ represents maximum value obtainable over all local measurements on A. The quantum part, known as the quantum discord, ${{\mathcal{D}}}_{B| A}={\mathcal{I}}-{{\mathcal{J}}}_{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 ${{\mathcal{J}}}_{B| A}$, we consider a set of von Neumann measurements ${\hat{{\rm{\Pi }}}}_{k=1,2}^{A}=\displaystyle \frac{1}{2}(1\pm {\vec{n}}_{k}\cdot {\vec{\sigma }}^{A})$ with ${| {\vec{n}}_{k}| }^{2}=1$, made on the subsystem A and minimize the corresponding conditional entropy ${{\mathcal{S}}}_{B| A}={\sum }_{k=1}^{2}\mathrm{Tr}[{p}_{k}S({\hat{\rho }}_{B| {{\rm{\Pi }}}_{k}^{A}})]$, where ${p}_{k}=\mathrm{Tr}[{\hat{{\rm{\Pi }}}}_{k}^{A}\hat{\rho }],{\hat{\rho }}_{B| {{\rm{\Pi }}}_{k}^{A}}$ $={\mathrm{Tr}}_{A}[{\hat{{\rm{\Pi }}}}_{k}^{A}\hat{\rho }]/{p}_{k}$ [61]. In figure 2(e) we calculate the mutual information from ${\mathcal{I}}$ and the quantum discord from ${\mathcal{D}}$ using as subsystems A and B a pair of dipoles, a and b respectively.