Cooperative Effects in Closely Packed Quantum Emitters with Collective Dephasing

TL;DR

This paper demonstrates that strong collective dephasing can restore cooperative effects in closely packed quantum emitters, enhancing dipole forces compared to non-interacting ensembles, with implications for solid state quantum systems.

Contribution

It reveals that collective dephasing can revive cooperative effects in dense quantum emitter ensembles, a phenomenon previously suppressed by dipole-dipole interactions.

Findings

Enhanced dipole force in strongly driven, collectively dephased ensembles.

Restoration of cooperative effects despite dense packing.

Relevance to solid state quantum systems like color centers and superconducting qubits.

Abstract

In a closely packed ensemble of quantum emitters, cooperative effects are typically suppressed due to the dephasing induced by the dipole-dipole interactions. Here, we show that by adding sufficiently strong collective dephasing cooperative effects can be restored. In particular, we show that the dipole force on a closely packed ensemble of strongly driven two-level quantum emitters, which collectively dephase, is enhanced in comparison to the dipole force on an independent non-interacting ensemble. Our results are relevant to solid state systems with embedded quantum emitters such as colour centers in diamond and superconducting qubits in microwave cavities and waveguides.

Figures

b, we show here the values for $\eta$ for each configuration. In the region $\bar{r}/\lambda_{0}\in\left[10^{-2},10^{-1}\right]$, $\langle\eta\rangle$ is not appreciably above $1$ although a vast majority of random configurations exhibit $\eta>1$. Some configuration even show large values of $\eta$, which are likely to be related to geometrical configurations with strong symmetries (see the discussion concerning the equidistant circular arrangement). The most remarkable aspect of this figure concerns the region of maximum CE. As already mentioned, this region is consistent with the optimal den

b, we show here the values for $\eta$ for each configuration. In the region $\bar{r}/\lambda_{0}\in\left[10^{-2},10^{-1}\right]$, $\langle\eta\rangle$ is not appreciably above $1$ although a vast majority of random configurations exhibit $\eta>1$. Some configuration even show large values of $\eta$, which are likely to be related to geometrical configurations with strong symmetries (see the discussion concerning the equidistant circular arrangement). The most remarkable aspect of this figure concerns the region of maximum CE. As already mentioned, this region is consistent with the optimal den

Equations

$${\hat{H}}_A\equiv \frac{\hslash \Omega (x)}{2}{\hat{S}}^x-\frac{\hslash \Delta }{2}{\hat{S}}^z.$$

$${\mathcal{L}}_H\hat{\rho }$$

$${\mathcal{L}}_{\Gamma }\hat{\rho }$$

$${\mathcal{L}}_{\gamma }\hat{\rho }$$

$$\begin{array}{rl}{g}_{mn}\equiv -\frac{3\Gamma }{4}& [\left(1-{\cos }^2{\theta }_{mn}\right)\frac{\cos \xi }{\xi }\\ & -\left(1-3{\cos }^2{\theta }_{mn}\right)\left(\frac{\cos \xi }{{\xi }^3}+\frac{\sin \xi }{{\xi }^2}\right)].\end{array}$$

$$\begin{array}{rl}{\Gamma }_{mn}\equiv \frac{3\Gamma }{4}& [\left(1-{\cos }^2{\theta }_{mn}\right)\frac{\sin (\xi )}{\xi }\\ & -\left(1-3{\cos }^2{\theta }_{mn}\right)\left(\frac{\sin \xi }{{\xi }^3}-\frac{\cos \xi }{{\xi }^2}\right)].\end{array}$$

$${\mathbb{F}}_{\mathrm{dp}}=-\frac{\hslash \nabla \Omega (x){|}_{x_0}}{2}⟨{\hat{S}}^x⟩,$$

$$\begin{array}{rl}⟨{\hat{S}}^x{⟩}_{\text{ind}}& =N\frac{4\Delta {\Omega }_0\Gamma }{\Gamma \left(4{\Delta }^2+{\gamma }_{\perp }^2\right)+2{\Omega }_0^2{\gamma }_{\perp }}\\ & \le N\frac{{\Omega }_0\Gamma }{\sqrt{\Gamma {\gamma }_{\perp }\left(\Gamma {\gamma }_{\perp }+2{\Omega }_0^2\right)}}.\end{array}$$

$${\Delta }_0\equiv -\sqrt{\frac{{\gamma }_{\perp }\left(\Gamma {\gamma }_{\perp }+2{\Omega }_0^2\right)}{4\Gamma }}.$$

$${\mathcal{L}}_H\hat{\rho }\equiv \frac{1}{2\text{i}}[{\Omega }_0{\hat{S}}^x-({\Delta }_0+\bar{g}){\hat{S}}^z+2\bar{g}{\hat{S}}^{+}{\hat{S}}^{-},\hat{\rho }].$$

$$\underset{{\gamma }_c/\Gamma \to \infty }{\lim }\left(\eta -1\right)$$

$$\underset{{\Omega }_0/\Gamma \to \infty ,{\gamma }_c=0}{\lim }\left(\eta -1\right)$$

$$\dot{\hat{\mu }}=\frac{1}{\text{i}\hslash }\left[{\hat{H}}_I,\hat{\mu }\right]+{\mathcal{L}}_{\Gamma }\hat{\mu }-\frac{\kappa }{4}([{\hat{S}}^{-},[{\hat{S}}^{+},\hat{\mu }]]+[{\hat{S}}^{+},[{\hat{S}}^{-},\hat{\mu }]]).$$

$$\frac{d\mathbf{X}}{dt}$$

$$\frac{d\mathbf{Y}}{dt}$$

$$\underset{{\gamma }_c/\Gamma \to \infty }{\lim }\left(\eta -1\right)$$

$$\underset{{\Omega }_0/\Gamma \to \infty ,{\gamma }_c=0}{\lim }\left(\eta -1\right)$$

$${\hat{S}}_z|m,{\alpha }_m⟩=2m|m,{\alpha }_m⟩,$$

$${\partial }_t\hat{\mu }=\mathcal{P}{\mathcal{L}}_1\mathcal{Q}\rho +{\mathcal{L}}_{1c}\hat{\mu }.$$

$$Q\hat{\rho }\approx {\int }_0^{\infty }\mathcal{Q}{e}^{{\mathcal{L}}_{0}t}{\mathcal{L}}_{1d}\mu ,$$

$${\partial }_t\hat{\mu }={\mathcal{L}}_{1c}\hat{\mu }+{\int }_0^{\infty }\mathcal{P}{\mathcal{L}}_{1d}{e}^{{\mathcal{L}}_{0}t}{\mathcal{L}}_{1d}\hat{\mu }.$$

$$\mathcal{P}{\mathcal{L}}_{1d}{e}^{{\mathcal{L}}_{0}t}{\mathcal{L}}_{1d}\hat{\mu }$$

$$-e^{{\mathcal{L}}_{0}t}{\hat{S}}^x\hat{\mu }{\hat{S}}^x+e^{{\mathcal{L}}_{0}t}\hat{\mu }{\hat{S}}^x{\hat{S}}^x].$$

$$\mathcal{P}{\hat{S}}^x{e}^{{\mathcal{L}}_{0}t}{\hat{S}}^x\hat{\mu }=e^{-{\gamma }_{c}t-i{\Delta }_{0}t}{\hat{S}}^{+}{\hat{S}}^{-}\hat{\mu }+e^{-{\gamma }_{c}t+i{\Delta }_{0}t}{\hat{S}}^{-}{\hat{S}}^{+}\hat{\mu }.$$

$${\hat{V}}_{\mathrm{ep}}$$

$${\hat{H}}_0/\hslash =\sum _{\mathbb{k}}{\omega }_k{\hat{a}}_{\mathbb{k}}^{†}{\hat{a}}_{\mathbb{k}}+\frac{{\omega }_0}{2}{S}^z.$$

$${\hat{V}}_{\mathrm{ep}}(t)$$

$${\hat{R}}_l(t)$$

$$k_x,k_y,k_z=0,\pm \frac{2\pi }{L},\pm \frac{4\pi }{L},\dots .$$

$$D({\omega }_k)=\{\begin{array}{ll}\frac{V{\omega }_k^2}{c_s^3(2\pi {)}^2},& \mathrm{when}{\omega }_k\le {\omega }_D.\\ 0,& \mathrm{when}{\omega }_k>{\omega }_D.\end{array}$$

$$\frac{d\rho }{dt}=\sum _{lm}$$

$$-i\left[-S_{ml}{\hat{\sigma }}_m^z{\hat{\sigma }}_l^z,\rho \right],$$

$${\gamma }_{ml}\approx \underset{{\omega }_k\to 0}{\lim }\frac{V{G}_0}{2\pi c_s^3}{\omega }_k^3\frac{\sin ({\omega }_k{\tilde{r}}_{ml})}{{\omega }_k{\tilde{r}}_{ml}}\frac{k_{B}T}{\hslash {\omega }_k},$$

$${\gamma }_{ml}\approx 4{\pi }^2{G}_0\frac{k_{B}T}{\hslash {\omega }_f}\frac{\sin (k_f{r}_{ml})}{k_f{r}_{ml}}.$$

Full text

A collection of two-level quantum emitters (TLEs) with sub-wavelength average separations can show remarkable cooperative behaviour like superradiant emission Agarwal74 ; Gross82 ; Ficek87 . The study of optical response in such systems has predominantly been restricted to the emission properties or the propagation of light within the TLE ensemble. This is because the systems usually considered in the early days Gross82 ; Malcuit87 , as well as in some recent works Bachelard11 ; Guerin16 ; AraujoKemp16 ; Jennewein16 , are gaseous clouds of atoms. With the advent of artificial atoms in solid state systems, e.g. quantum dots Brandes05 , superconducting qubits Nissen13 ; Mlynek14 and colour centers in diamond NVreview ; Juan16 ; Bradac16 , it is now possible to study the impact of cooperative effects on other aspects of the optical response. In particular, a recent experiment Juan16 studied the dipole force on optically trapped nanodiamonds containing a high density of Nitrogen vacancy (NVs) centers. An intriguing result in Juan16 was that the observed dipole force originating from the emitters could not be correctly accounted for by considering the emitters to respond independently.

In this work, we focus on cooperative effects in a small and closely packed ensemble of TLEs subject to strong coherent driving and collective dephasing. In particular, we show that the dipole force on such an ensemble can be larger than on an equivalent one where each TLE spontaneously emits independently. For the emitter separations that we consider here, the dipole-dipole interaction can be larger than the line-width of the individual emitters. Furthermore, spontaneous emission is not perfectly collective. In this situation, one typically expects cooperative effects to be suppressed Gross82 ; Ficek83 . Here, we show that the combination of strong driving and large collective dephasing can restore cooperative effects, even in the presence of dipole interaction shifts and non-collective spontaneous emission. While there have been previous studies of cooperative effects with strong driving fields FicekTanasReview02 ; AminCordes78 ; Agarwal79 ; Agarwal80 ; Drummond80 ; Carmichael82 ; Campagno82 , the role of collective dephasing has received less attention Nissen13 ; Juan16 . In the context of Quantum Information, collective decoherence in general, and collective dephasing in particular, has been studied both theoretically Palma96 ; Duan98 ; Lidar03 and experimentally Kwiat00 ; Kielpinski01 . There, particular attention was paid to the existence and robustness of so called decoherence free sub-spaces (DFS) under collective dephasing Lidar03 . Moreover, recent studies Lanyon13 ; Carnio15 have shown that collective dephasing could also be used as a resource to generate strong but separable correlations. We note that, due to their promise as a general passive strategy to protect quantum resources from noise, the study of DFS continues to be an active area of research, see Lidar12 for a recent review on the theoretical aspects and DFSexpts for experimental implementations. Of late, novel applications of DFS, such as generation of arbitrary photonic states Tudela15 and universal quantum computation in waveguide QED Tudela16 as well as quantum repeaters with trapped ions, Zwerger17 have also been proposed.

Let us consider a collection of $N$ identical TLEs with resonance frequency $\omega_{0}\equiv ck_{0}\equiv 2\pi c/\lambda_{0}$. The matrix element of the dipole moment operator $\hat{\mathbb{d}}$ is given by $\langle{e}|{\hat{\mathbb{d}}}|{g}\rangle\equiv\boldsymbol{\epsilon}_{a}d$, where $|{g}\rangle$ ($|{e}\rangle$) denotes the ground (excited) state and $\left|{\boldsymbol{\epsilon}_{a}}\right|=1$. For simplicity, we assume $\boldsymbol{\epsilon}_{a}$ and $d$ to be real. The TLEs, each with position $\mathbb{r}_{m}$ ($m=1,\ldots,N$), are fixed on a background matrix with center-of-mass position $\mathbb{x}$, namely $\mathbb{r}_{m}^{\prime}\equiv\mathbb{r}_{m}-\mathbb{x}$ for all $m$ is a constant of motion. We consider the TLEs to be driven by a classical electromagnetic field of the form $\mathbb{E}(\mathbb{r},t)=\mathbb{E}(\mathbb{r})\cos\left(\omega_{d}t\right)\equiv\boldsymbol{\epsilon}_{d}E_{0}f(\mathbb{r})\cos\left(\omega_{d}t\right)$, where $\boldsymbol{\epsilon}_{d}$ and $E_{0}$ are real and $\left|{\boldsymbol{\epsilon}_{d}}\right|=1$. We assume that $f(\mathbb{r}_{m})\approx\mathbb{f}(\mathbb{x})$ for all $m$. In a frame rotating with the drive and assuming the rotating wave approximation, the hamiltonian describing the interaction of the identical TLEs with $\mathbb{E}(\mathbb{r},t)$ is given by

[TABLE]

Here, $\Omega(\mathbb{x})\equiv-2dE_{0}(\boldsymbol{\epsilon}_{a}\cdot\boldsymbol{\epsilon}_{d})f(\mathbb{x})/\hbar$ is the Rabi frequency and $\Delta\equiv\omega_{d}-\omega_{0}$ the detuning. Hereafter, we use the following notation for the spin operators: $\hat{\sigma}^{\alpha}_{m}$ denotes the Pauli matrix (for $\alpha={x,y,z}$) and the ladder operator (for $\alpha=\pm$) of the $m$th TLE, collective operators are denoted by $\hat{S}^{\alpha}\equiv\sum_{m=1}^{N}\hat{\sigma}^{\alpha}_{m}$. Apart from the dynamics induced by the interaction with the external driving, we assume the TLEs to experience collective dephasing as well as dipole-dipole interaction and spontaneous emission due to the interaction with free electromagnetic field modes in the vacuum state. The overall dynamics of such an ensemble of identical TLEs can be then described by the master equation $\dot{\hat{\rho}}=\mathcal{L}\hat{\rho}\equiv\left({\mathcal{L}_{H}+\mathcal{L}_{\Gamma}+\mathcal{L}_{\gamma}}\right)\hat{\rho}$, where Agarwal74 ; Gross82 ; Ficek87

[TABLE]

The term $\mathcal{L}_{H}$ describes coherent dynamics given by the interaction with the external field, Eq. (1), and the dipole-dipole interaction given by $\hat{H}_{I}\equiv\sum_{m\neq n}\hbar g_{mn}\hat{\sigma}_{m}^{+}\hat{\sigma}_{n}^{-}$, where Kraemer16

[TABLE]

Here, $\xi=k_{0}r_{mn}\equiv k_{0}\left|{\mathbb{r}_{m}-\mathbb{r}_{n}}\right|$, and $\cos\theta_{mn}\equiv(\mathbb{r}_{m}-\mathbb{r}_{n})\cdot\boldsymbol{\epsilon}_{a}/r_{mn}$. The term $\mathcal{L}_{\Gamma}$ describes the spontaneous emission of the TLEs with correlated emission rates given by

[TABLE]

The diagonal term $\Gamma_{mm}\equiv\Gamma/2=2d^{2}\omega_{0}^{3}/(3\hbar c^{3})$ is the individual spontaneous emission rate of the TLEs. The term $\mathcal{L}_{\gamma}$ describes collective dephasing with a rate given by $\gamma_{c}$. It is convenient to introduce the rate $\bar{g}\equiv\sum_{n\neq 1}\left|{g_{1n}}\right|$, which parameterizes the strength of the dipole-dipole interaction, and $\overline{\Gamma}\equiv\sum_{n\neq 1}\Gamma_{1n}/(N-1)$, which parameterizes the cooperativity of the spontaneous emission. The physical origin of collective dephasing is left unspecified in the theoretical treatment here. Note that it can, for instance, arise via correlated magnetic field fluctuations for ions Lanyon13 or due to interactions with phononic baths in the case of colour centers Palma96 ; Fu2009 ; Albrecht13 (see Supplement for additional details regarding such situations).

We are interested in the closely packed regime defined by $k_{0}r_{mn}\leq 1$ for any pair of TLEs. While in this regime the spontaneous emission is predominantly collective, namely $\bar{\Gamma}\lesssim\Gamma/2$, the strong dipole-dipole interaction $\bar{g}\gg\Gamma$ typically suppresses any cooperative effect. However, we show below that strong collective dephasing $\gamma_{c}\gg\Gamma$ together with strong driving $\Omega(\mathbb{x})\gg\Gamma$ can recover cooperative effects. In particular, we concentrate on the steady state value of $\hat{S}^{x}$, namely $\langle\hat{S}^{x}\rangle\equiv\textrm{tr}[\hat{S}^{x}\hat{\rho}_{s}]$, where $\mathcal{L}\hat{\rho}_{s}=0$. This is related to the dipole force exerted by the driving field to the matrix hosting the TLEs. Indeed, assuming that the motion of the background matrix in the applied field is slow compared to the emitter dynamics, the dipole force is given by

[TABLE]

where $\mathbb{x}_{0}$ is the equilibrium position of the matrix. For an ensemble of $N$ independent TLEs, namely with $g_{mn}\equiv 0$ and $\Gamma_{mn}\equiv 0$ (for $m\neq n$), one has that

[TABLE]

Here, $\Omega_{0}\equiv\Omega(\mathbb{x}_{0})$, $\gamma_{\perp}\equiv(\Gamma+2\gamma_{c})$, and the upper bound is achieved at the optimal detuning

[TABLE]

We are interested in the parameter $\eta\equiv\langle\hat{S}^{x}\rangle/\langle\hat{S}^{x}\rangle_{\text{ind}}$ evaluated at $\Delta=\Delta_{0}$ (note that $\langle\hat{S}^{x}\rangle$ is not maximized at $\Delta_{0}$). In particular, we refer to situations when $\eta>1$ as cooperative enhancement (CE). We remark that since we are interested in closely packed ensembles, we do not consider variations of the Rabi frequency with emitter locations $\mathbb{r}_{m}$, which can also lead to interesting modifications of collective effects Zoubi10 ; AsenjoGarcia17 .

Let us consider $N$ identical TLEs, randomly positioned in a three dimensional volume with an average separation given by $\bar{r}\equiv\sum_{m>n}r_{mn}/N$, see Fig. 1a. We generate multiple random configurations at a fixed $\bar{r}$ using the procedure described in Damanet16 (see Supplement for details) and numerically calculate $\hat{\rho}_{s}$ in each case using Qutip13 . The average $\eta$ over 1000 configurations is plotted in Fig. 1b as a function of $\bar{r}/\lambda_{0}$ for $N=6$ with (solid line) and without (dotted line) collective dephasing. The shaded regions correspond to the interval where a majority, 68$\%$, of the values for $\eta$ lie. In the absence of collective dephasing ($\gamma_{c}=0$), there is no CE ($\eta\leq 1$). However, in the presence of strong collective dephasing ($\gamma_{c}/\Gamma\approx 10^{4}$) and strong driving ($\Omega_{0}/\Gamma=10^{3}$) there is a range of separations $\bar{r}$ where there is CE ($\eta>1$). This is the main finding of this work. In the inset of Fig. 1b, the parameters $\bar{g}/\Gamma$ and $2\bar{\Gamma}/\Gamma$ as a function of $\bar{r}$ (averaged over 1000 configurations) are plotted. Note that CE vanishes both at large average separations due to the non-collective nature of spontaneous emission and at small distances due to the large dipole-dipole interaction.

Let us analytically support these statements for the simplest $N=2$ case Ficek83 . The hamiltonian including the dipole interaction is

[TABLE]

In this case, the hamiltonian is collective and commutes with $\hat{\mathbf{S}}^{2}$. In contrast, in the finite emitter separation regime, where the so-called small sample limit Dicke54 ; Agarwal74 ; AminCordes78 ; Gross82 ; FicekTanasReview02 can not be used, the spontaneous emission Eq. (3) is still non-collective i.e. $\chi\equiv 2\overline{\Gamma}/\Gamma<1$. Following Ficek83 and as shown in Supplement , one can analytically calculate $\eta$. Fig. 2 plots the CE region ($\eta>1$) in the plane $(\Omega_{0}/\Gamma,\gamma_{c}/\Gamma)$ for different $\bar{r}$. Note that CE requires both large dephasing and large driving. Furthermore, from the lengthy analytical expression for $\eta$, one can obtain that in the limit of large dephasing $\eta$ reads

[TABLE]

In this limit, CE ($\eta>1$) requires sufficiently large driving $\Omega_{0}>\sqrt{\Gamma^{2}+2\overline{\Gamma}\Gamma}$. Alternatively, one can also show that in the limit of no dephasing ($\gamma_{c}=0$) and large driving, $\eta$ reads

[TABLE]

and hence there is no CE ($\eta<1$). This is in agreement with previous studies of resonance fluorescence AminCordes78 ; FicekTanasReview02 ; Agarwal79 ; Agarwal80 at strong driving, which results from the increased occupation of the bright Dicke subspace in the small sample limit. It was shown in Ficek83 that when $\chi<1$ such an enhancement is absent. Interestingly, we claim here that large collective dephasing can restore cooperative effects.

Returning to Fig. 1b, $N=6$, we observe that there is an optimal separation distance $\bar{r}/\lambda_{0}$ where $\eta$ reaches a maximum. The dipole-dipole interactions for $N>2$ do not conserve permutation symmetry in general. As a result, large dipole interactions induce strong local dephasing apart from providing energy shifts that prevent the TLEs to be polarized at the chosen detuning $\Delta_{0}$. We see in Fig. 1b that in regions with CE the dipole interactions satisfy $\gamma_{c}\gg\overline{g}\sim\Gamma$. In this manner we can understand the maximum $\eta$ in Fig. 1b as occurring at separations where the dipole interactions are small enough to not dephase the collective behaviour fostered by the cooperative spontaneous emission and collective dephasing. This statement can be quantified by noting that the product $g_{12}(1-2\Gamma_{12}/\Gamma)$, for a pair of TLEs with parallel moments and $\theta_{12}=\pi/2$, is locally minimized at $r_{12}\approx 0.2\lambda_{0}$. This is consistent with the position of the peak in Fig. 1b. Furthermore, this also agrees with the observation that the location of the peak remains in the same density region for various ensemble sizes (see Fig. 3) and driving strengths $\Omega_{0}$ (see Supplement ).

For $N>6$, an exact numerical calculation of $\hat{\rho}_{s}$ becomes rapidly intractable with standard resources. State-of-the art approximation methods in the field, such as Holstein-Primakoff Jenkins12 ; Bettles16a ; Lee16 or extended mean-field Kraemer15 , are not valid in the strong driving ($\Omega_{0}\gg\Gamma$) and strong dipole-dipole interaction ($\bar{g}\gg\Gamma$) regime respectively. The separation of time scales between the spontaneous emission and collective dephasing $\gamma_{c}\gg\Gamma$ also makes trajectory methods  Clemens03 unsuitable. For small emitter separations up to $\bar{r}/\lambda_{0}\sim 10^{-1}$, we find that a numerical diagonalization of the hamiltonian $\hat{H}_{A}+\hat{H}_{I}$ followed by a secular approximation to convert the master equation to a rate equation in the dressed basis Dalibard85 allows us to go up to $N=10$ emitters. For $\bar{r}/\lambda_{0}\gtrsim 10^{-1}$ the weak dipole-dipole interactions do not appreciably lift the degeneracy between the Dicke states with the same angular momentum projection $m$ but different total angular momentum $S$, rendering the dressed approach invalid. In this regime we proceed as follows. The Liouvillian of the master equation can be written as $\mathcal{L}=\mathcal{L}_{0}+\mathcal{L}_{1}$, where $\mathcal{L}_{0}\hat{\rho}\equiv\mathcal{L}_{\gamma}\hat{\rho}+\text{i}\Delta_{0}[\hat{S}^{z},\hat{\rho}]/2$ and $\mathcal{L}_{1}\equiv\mathcal{L}_{\Gamma}\hat{\rho}+(\text{i}\hbar)^{-1}[\hat{H}_{I}+\hbar\Omega_{0}\hat{S}^{x}/2,\hat{\rho}]$. In most of the regime where $\eta>1$, one has that $\gamma_{c}$ and $\Delta_{0}$ are much larger than $\bar{\Gamma}$, $\bar{g}$, and $\Omega_{0}$. Indeed, note that $\Delta_{0}\approx\Omega_{0}\sqrt{\gamma_{c}/\Gamma}$ when $\Omega_{0}^{2}\gg\Gamma\gamma_{c}$, which according to Eq. (11), is required to have an appreciable value of $\eta-1>0$. Under these assumptions, $\mathcal{L}_{0}$ describes faster dynamics than $\mathcal{L}_{1}$, and hence, one can adiabatically eliminate Cirac92 ; Lesanovsky13 the fast dynamics. This leads to an effective master equation in the dark subspace of $\mathcal{L}_{0}$, namely for states $\hat{\mu}$ such that $\mathcal{L}_{0}\hat{\mu}=0$. As shown in Supplement , the effective master equation is given by

[TABLE]

Here $\kappa\equiv\Omega_{0}^{2}\gamma_{c}/\left(\gamma_{c}^{2}+\Delta_{0}^{2}\right)$ is of the order of $\Gamma$ in the assumed parameter regime. Eq. (13) can be conveniently solved numerically via trajectory unravellings Clemens03 .

In Fig. 3, the results for the averaged $\eta$ over multiple random configurations, $\langle\eta\rangle$, for $8$ and $10$ TLEs are presented. For $\bar{r}/\lambda_{0}<10^{-1}$, $\langle\eta\rangle$ was calculated by averaging over $400$ random configuration with steady states calculated using the rate equation method. For $\bar{r}/\lambda_{0}>10^{-1}$, the mean was calculated over $200$ configurations from the steady state solutions determined by averaging over 500 trajectories each Qutip13 .

The statistical distribution of $\eta$, and checks to ensure that the approximate methods used for $N=8,10$ in Fig. 3 agree with the exact results for $N\leq 6$ are presented in Supplement . From the inset in Fig. 3, we see that while $\langle\eta\rangle$ increases with $N$ in the CE region up to $N=8$, there is no appreciable gain for $N=10$. In Supplement , we also demonstrate CE for an equidistant circular arrangement of TLEs. Since dipole-dipole interactions are permutation symmetric in this case, we find a significant CE also at smaller separations than in the random arrangement. We remark that the choice of detuning $\Delta_{0}$ maximizes $\langle S^{x}\rangle$ for independent emitters and that hence, an optimized choice for the collective case can certainly lead to even larger $\eta$. A more systematic study of the CE as a function of different arrangements of TLEs Kraemer16 , as well as developing efficient methods to numerically and analytically address larger number of TLEs is left for future work.

We conclude with some remarks relating our findings to the recent studies of dipole force on colour centres embedded in nano-diamonds Juan16 ; Bradac16 . In the experiment Juan16 , the underlying mechanism for the large dephasing at room temperature is mediated via phonon interactions Fu2009 and consequently changes rapidly with the temperature of the lattice. Since we have demonstrated that the presence of large collective dephasing is crucial to the observation of the enhanced dipole force, this raises the prospect of repeating the experiment Juan16 , or even life-time measurements in Bradac16 , at lower temperatures. At lower temperatures, the dephasing will be reduced which should lead to a strong modification or even suppression of collective effects. This is counter-intuitive to the study of collective effects in atomic systems. In connection to the proposal for levitated optomechanics with nano-diamonds Juan16a (as well as other proposals Xuereb12 concerning collective effects in optomechanics), it is interesting to explore if collective effects in dense ensembles lead to polarizabilities comparable or greater than the bulk polarizability of the embedding medium. A promising direction for further research is to explore other scenarios where large collective dephasing restores cooperative effects. Remarkably, in systems such as superconducting qubits Nissen13 ; Mlynek14 where the collective dephasing can be externally controlled, this could allow to observe cooperative effects even in the presence of inhomogeneities and dipole shifts.

This work is supported by the Austrian Federal Ministry of Science, Research,and Economy (BMWFW). We acknowledge discussions with J. I. Cirac, J. J. Garcia-Ripoll, C. Gonzalez-Ballestero, G. Kirchmair, H. Ritsch, and B. Vermersch.

I Random configuration choice - Details

Following SDamanet16 , in order to generate a random configuration with a given average separation $\bar{r}$, we first pick $N$ points uniformly distributed over a spherical container of arbitrary radius $R$. Following this, we rescale all distances by the average separation and multiply all co-ordinates by the prescribed average $\bar{r}$. If any two emitters are separated by a distance less than a cutoff (chosen as $10^{-4}\lambda_{0}$), the configuration is dropped. We note that the inset plot in Fig. 1 (b) of the main paper showing the mean and confidence intervals of the distribution for $\bar{g}\equiv\sum_{n\neq 1}\left|{g_{1n}}\right|$ and $\overline{\Gamma}\equiv\sum_{n\neq 1}\Gamma_{1n}/(N-1)$ has no qualitative dependence on which emitter we choose to index by $m=1$ within the ensemble since we are considering random distributions.

II Analytical calculation of $\eta$ for $N=2$

Our method closely follows the treatment in SFicek83 that allows us to solve the $N=2$ emitter case exactly without any approximations i.e. $\overline{\Gamma}\neq\Gamma/2$ and $\overline{g}\neq 0$. The idea in SFicek83 is to write down the equations of motion for the operator averages $\langle\hat{\sigma}_{m}^{\alpha}\rangle$, $\langle\hat{\sigma}_{m}^{\alpha}\hat{\sigma}_{n}^{\beta}\rangle$. For the $N=2$ case there are $15$ unique averages, exactly equal to the number of unique matrix elements of the density matrix. Thus solving this set of equations gives an exact solution to the master equation. Moreover, the $15$ equations also split up into two groups of $9$ and $6$ equations - one for permutation symmetric averages, such as $\langle\hat{\sigma}_{1}^{x}\rangle+\langle\hat{\sigma}_{2}^{x}\rangle,\langle\hat{\sigma}_{1}^{x}\hat{\sigma}_{2}^{y}\rangle+\langle\hat{\sigma}_{2}^{x}\hat{\sigma}_{1}^{y}\rangle$, denoted by the column vector $\mathbf{X}$, and anti-symmetric ones, such as $\langle\hat{\sigma}_{1}^{x}\rangle-\langle\hat{\sigma}_{2}^{x}\rangle$, denoted by the column vector $\mathbf{Y}$. The equation of motions have the form SFicek83 :

[TABLE]

Since $M_{Y}$ turns out to be invertible, $\mathbf{Y}$ vanishes in the steady state. Thus in order to solve for the steady state we just need to invert the matrix $M_{X}$ in Eq. (S14). Interestingly, we find that the determinant of $M_{X}$ is proportional to $\Gamma/2-\overline{\Gamma}$. This means we have two distinct solutions depending on whether $\overline{\Gamma}=\Gamma/2$ (collective) or not SFicek83 ; SFicekTanasReview02 . This is so since for $\overline{\Gamma}=\Gamma/2$, $\hat{\mathbf{S}}^{2}$ is a constant of motion and one of the equations in Eq. (S14) becomes redundant. We do not give the exact analytical expressions for $\eta$ as they are cumbersome (but have plotted the same in Fig. 2 of the main text) and consider only limiting cases which exhibit CE. We presented the result for the case with $\overline{\Gamma}\neq\Gamma/2$ in Eqs. (11,12) of the main paper in the limits $\gamma_{c}/\Gamma\rightarrow\infty$. An additional comment in this case is in order. From Fig. 2 of the main paper it is apparent that at large but finite values of $\gamma_{c}/\Gamma$, in order to have $\eta>1$ we will need $\Omega_{0}$, even larger than the value suggested on the RHS of Eq. (11), in order to overcome the dipole interaction induced energy shifts. In additionm we also note that in the case with perfect collective interaction $\overline{\Gamma}=\Gamma/2$, the equivalent expressions to Eqs. (11,12) are:

[TABLE]

From the above equation, it is clear that, unlike in the case with $\overline{\Gamma}\neq\Gamma/2$, provided $\Omega_{0}>\sqrt{2}\Gamma$ and $\Omega_{0}>112\sqrt{16\overline{g}^{2}+15\Gamma^{2}}/225$, there is CE both in the case with and without $\gamma_{c}$.

III Derivation of Effective Master Equation (13)

In order to detail the procedure for adiabatic elimination SCirac92 ; SLesanovsky13 , it is useful to first choose a basis that commutes with the faster Liouvillian $\mathcal{L}_{0}$. The eigenstates of $\hat{S}_{z}$ defined as:

[TABLE]

with $-N/2\leq m\leq N/2$ and $1\leq\alpha_{m}\leq\binom{N}{N/2+m}$ provide the required basis. Note that these are not the Dicke states. The idea then is to derive an equation for $\hat{\mu}=\mathcal{P}\rho$, where $\mathcal{P}$ is the the projector to the dark subspace of $\mathcal{L}_{0}$ and can be formally written as $\lim_{t\rightarrow\infty}e^{\mathcal{L}_{0}t}$. From the formal expression for $\mathcal{P}$ we can see immediately that the projection essentially suppresses all the block off-diagonal elements of the density matrix in the basis Eq. (S18) i.e. $\langle{m}|{e^{\mathcal{L}_{0}t}}|{n}\rangle=e^{i\Delta_{0}(m-n)t-\gamma_{c}(m-n)^{2}t}$. To determine the evolution of $\hat{\mu}$, we rewrite the density matrix as $\hat{\rho}=\mu+\mathcal{Q}\hat{\rho}$, with $\mathcal{Q}=\mathcal{I}-\mathcal{P}$ and use the full master equation to write separate equations for $\hat{\mu}$ and $\mathcal{Q}\hat{\rho}$. At this point it is also convenient to split the slow Liouvillian as $\mathcal{L}_{1}=\mathcal{L}_{1c}+\mathcal{L}_{1d}$. Here $\mathcal{L}_{1c}\rho=\mathcal{L}_{\Gamma}\rho+(\text{i}\hbar)^{-1}[\hat{H}_{I},\hat{\rho}]$ commutes with $\mathcal{L}_{0}$, and we have $\mathcal{L}_{1d}\rho=(\text{i}\hbar)^{-1}[\hbar\Omega_{0}\hat{S}^{x},\hat{\rho}]$. Using the properties $\mathcal{P}\mathcal{L}_{0}\mathcal{P}=\mathcal{P}\mathcal{L}_{0}\mathcal{Q}=\mathcal{P}\mathcal{L}_{1d}\mathcal{P}\equiv 0$, and $\mathcal{P}\mathcal{L}_{1}\mu=\mathcal{L}_{1c}\mu$, we can derive

[TABLE]

The key step in adiabatic elimination is to constrain the dynamics of $\mathcal{Q}\hat{\rho}$ to $\hat{\mu}$ via $\mathcal{Q}\hat{\rho}\approx-(\mathcal{Q}\mathcal{L}_{1}+\mathcal{Q}\mathcal{L}_{0})^{-1}\mathcal{Q}\mathcal{L}_{1}\hat{\mu}$ SCirac92 ; SLesanovsky13 , and using properties of Laplace transform to write:

[TABLE]

to first order in $\mathcal{L}_{1d}$. Substituting Eq. (S20) in Eq. (S19) and using $\mathcal{P}\mathcal{L}_{1}\mathcal{Q}=\mathcal{P}\mathcal{L}_{1d}$, we obtain:

[TABLE]

In order to simplify the above equation we note that

[TABLE]

Using $\hat{S}^{x}=\hat{S}^{+}+\hat{S}^{-}$, and the property $\langle{m,\alpha_{m}}|{\hat{S}^{\pm}}|{n,\alpha_{n}}\rangle\propto\delta_{m,n\pm 1}$ we can show, for example, that:

[TABLE]

Constructing similar simplifications for the rest of the terms in Eq. (S22), we arrive at Eq. 13 of the main paper.

IV Approximate methods

Fig. S1 compares $\langle\eta\rangle$ obtained from approximate methods (rate equation averaged over $1000$ random configurations and trajectory simulations averaged over $200$ configurations with 500 trajectories each) for $N=4,6$ emitters to that from a direct steady state calculation for the full master equation (averaged over $1000$ configurations). The agreement is good and suggests that the approximate methods used in Fig. 3 of the main paper for larger $N$ can be trusted. The shaded region in Fig. S2 displays the interval where $60\%$ of $\eta$ values lie for the $N=8,10$ case. In Fig. 3 of the paper we plotted only the averaged value $\langle\eta\rangle$. Note that even when $\bar{r}/\lambda_{0}>10^{-1}$, there are some random configurations that lead to very large $\gamma_{c}>\overline{g}\gg\Gamma$ leading to very slow evolution of the trajectory simulations. For such isolated configurations we calculate the steady state via a rate equation approximation to the effective master equation.

V Additional Numerical Results

Fig. S3 depicts the $\langle\eta\rangle$ for $N=6$ randomly distributed emitters as a function of the drive strength $\Omega_{0}$. This demonstrates that in the random distribution case, increasing $\Omega_{0}$ does not lead to enhancement at smaller separations due to the strong local dephasing induced by dipole-dipole interactions. We note that the position of the CE peak is not affected by the drive strength. This is consistent with the condition obtained in the main paper by analyzing the optimal separation for two emitters within the ensemble giving $r_{12}\approx 0.2\lambda_{0}$, independent of the drive.

Fig. S4 presents the numerical results for the 1000 random configurations for $N=6$. In addition to the mean value for $\eta$, already presented in Fig. 1b, we show here the values for $\eta$ for each configuration. In the region $\bar{r}/\lambda_{0}\in\left[10^{-2},10^{-1}\right]$, $\langle\eta\rangle$ is not appreciably above $1$ although a vast majority of random configurations exhibit $\eta>1$. Some configuration even show large values of $\eta$, which are likely to be related to geometrical configurations with strong symmetries (see the discussion concerning the equidistant circular arrangement). The most remarkable aspect of this figure concerns the region of maximum CE. As already mentioned, this region is consistent with the optimal density parameter anticipated from considering just $2$ emitters, $r_{12}\approx 0.2\lambda_{0}$, and moreover we see that all the configurations provide CE. This suggests $\bar{r}\approx 0.2\lambda_{0}$ to be a fairly general criterion to maximize the CE.

Fig. S5 illustrates the possibility of controlling $\eta$ for a given random distribution of emitters by tuning the drive strength $\Omega_{0}$. As in the case of $N=2$ emitters, depicted in Fig. 2 of the main paper, at a given value of the collective dephasing there is enhancement only when the drive strength is larger than some critical value that will in general vary from realisation to realisation.

VI Equidistant circular arrangement of TLEs

In Fig. S6, we consider $\eta$ for an equispaced arrangement of emitters on a two dimensional circle (see inset in Fig. S6 for a schematic of this configuration) as a function of the average separation $\bar{r}$. For $N=4,6$, $\eta$ was calculated from the steady state of the full master equation without any approximation. For $N=8,10$ ($N=12$), similar to the random configuration case, we used the rate equation approximation to the full master equation (effective master equation) for small separations and trajectory simulations for larger separations. Since the effective master equation is valid only when $\overline{g}\ll\gamma_{c}$, we do not consider very small separations for $N=12$. We can see from the results that for such a permutation symmetric configuration SGross82 , we find stronger CE that also extends over a larger region in $\bar{r}$ compared to the random distribution results, e.g. Fig. 1b of the main paper.

VII Collective dephasing master equation

We consider two situations in which a collective dephasing Markovian master equation of the form in Eq. (4) of the main paper may be justified. In the first example the TLEs of interest are NV centers in diamonds and the dephasing arises from their interaction with phonons in the diamond lattice. The interaction with phonons arises due to the modification of the energy level structure of the TLE under mechanical strain SAlbrecht13 , which in turn can be expressed in terms of lattice phonons. Let us now specialise to large wavelength acoustic phonons, that have the dispersion relation $\omega(\mathbb{k})=\omega_{k}=c_{s}|\mathbb{k}|$, with $c_{s}=1.2\times 10^{4}$ m s*-1*, the speed of sound in diamond. For the sake of simplicity, let us also restrict to just the longitudinal polarization modes parallel to the propogation momentum $\mathbb{k}$. Then we obtain the following TLE-phonon interaction hamiltonian:

[TABLE]

with $\hat{a}^{\dagger}_{\mathbb{k}}$ denoting the creation for a phonon with momentum $\mathbb{k}$, $M$ is the mass of the crystal, $u_{\mathbb{k}}=i\frac{\zeta}{2}\sqrt{\frac{\hbar}{2M\omega_{k}}}(k_{x}^{2}+k_{y}^{2}+8\Xi^{2}k_{z}^{2})/k$, and $k=|\mathbb{k}|$. In the above, $\zeta\simeq 610$THz is an energy scale emerging from the strain coupling of the NV center and $\Xi$ is a factor introduced to differentiate between the Nitrogen and Carbon atoms that contribute to the coupling SAlbrecht13 . Since we are interested in a simple description capturing only the key physics, we will set $8\Xi^{2}\equiv 1$ to obtain an isotropic coupling $u_{\mathbb{k}}=u_{k}=i\frac{\zeta}{2}\sqrt{\frac{\hbar}{2M\omega_{k}}}k$. The above hamiltonian has the same form as the ones studied in the context of pure dephasing of qubits in quantum information SPalma96 ; SLidar01 and is also known as the spin-boson model. Let us now detail some of the steps towards eliminating the phononic bath and constructing a master equation SBreuerPet02 ; SCarmichael99 for the TLEs. First we note that for simplicity we will neglect the driving of the TLEs (we comment on this aspect in the end) and consider the following bare hamiltonian for the system:

[TABLE]

The interaction picture hamiltonian is then given by:

[TABLE]

We follow the usual Born-Markov approach to derive a master equation for the reduced density matrix of the TLEs denoted by $\hat{\rho}(t)$. To this end we assume that the total state of the system and the phonon reservoir is a product state of the form $\hat{\rho}_{tot}(t)\approx\hat{\rho}(t)\otimes\hat{\rho}_{R}(0)$. The reservoir is at its initial state, chosen as a thermal state $\hat{\rho}_{R}(0)=\exp\left(-\beta\hbar\sum_{\mathbb{k}}\omega_{k}\hat{a}^{\dagger}_{\mathbb{k}}\hat{a}_{\mathbb{k}}\right)$ with $\beta=1/(k_{B}T)$, and is essentially affected very little by the system at all times. This is the content of the Born approximation, which is valid here since the typical magnitude of $u_{k}\sim 2\pi\times 0.3$ MHz (for a cubic diamond sample of size $\lambda_{0}$) and is much less than the optical transition frequency of the NV center. In order to justify the validity of the Markov approximation, we have to examine the two-time correlations of the reservoir operators $\hat{R}_{l}(t)$. For this, it is essential to choose a more specific model for the phonon density of states. For the sake of simplicity, we consider a cubical region of finite volume $V=L^{3}$ with periodic boundary conditions along all three directions. This immediately leads to the following allowed values of the wavenumbers SKittel96 :

[TABLE]

We will assume the Debye model and introduce a cut-off frequency $\omega_{D}=c_{s}(6\pi^{2}n)^{1/3}$ with $n$ the number density of diamond SKittel96 . The density of states is then given by:

[TABLE]

It is then straight forward to evaluate the correlation function $R_{ml}=\langle\hat{R}_{m}(0)\hat{R}_{l}(\tau)\rangle$ and it can be shown that it decays over the time scale $1/\omega_{T}$, where $\omega_{T}=k_{B}T/\hbar$ SCarmichael99 . Thus the Markovian approximation is justified as long as the system dynamics is over time scales slower than $1/\omega_{T}$. The master equation within the Born-Markov approximation then reads:

[TABLE]

with $\gamma_{ml}=R_{ml}+R_{lm}^{*}$ representing a positive definite dephasing matrix $\gamma$ and $S_{ml}=\left(R_{ml}-R_{lm}^{*}\right)/(2i)$ representing the hermitian Lamb-shift matrix SBreuerPet02 . The explicit for the dephasing rate is given by:

[TABLE]

with $G_{0}=\hbar\zeta^{2}/(8Mc_{s}^{2})$, which gives a measure of the NV center phonon coupling, and $\tilde{r}_{ml}=r_{ml}/c_{s}$. If we naively apply the limit above it is clear that we get $\gamma_{ml}=0$. But, since we are in a finite crystal with periodic boundary condition and the allowed phonon mode of lowest frequency has a wavenumber $k_{f}=2\pi/L=\omega_{f}/c_{s}$, we impose an infra-red cut off to the limit and evaluate the limit Eq. (S31) at $\omega_{k}=\omega_{f}$ to get:

[TABLE]

If the average emitter separation is much smaller compared to the fundamental mode of the diamond, the factor $\sin(k_{f}r_{ml})/k_{f}r_{ml}\sim 1$ (for example this factor is about $0.93$ for $k_{f}r_{ml}\approx 0.1\times 2\pi$) and $\gamma_{ml}=\gamma_{c}$ and the dephasing becomes collective, with:

[TABLE]

This is intuitive as it is simply given by the frequency $G_{0}$ characterizing the NV center phonon coupling times the occupation number $k_{B}T/(\hbar\omega_{f})$ of the lowest phonon mode. It is also clear that the Markov approximation is valid since $\gamma_{c}/\omega_{T}=4\pi^{2}G_{0}/\omega_{f}$ and is of the order of $10^{-8}$ for a diamond sample with $L=\lambda_{0}$ (the optical wavelength). For this case, we also get a collective dephasing rate of $\gamma_{c}=2\pi\times 0.1$ MHz at room temperature. We note that this collective dephasing rate is small compared to the spontaneous emission rate of the NV center optical transition and not as large as the values we use in the main paper. This is a limitation of the rather simple two-level model of the quantum emitter we have used to describe phonon induced dephasing. In priniciple, the excited state of NV centers are multiplets and phonon interactions could mediate coupling between such excited state multiplets. This can then lead to large dephasing seen in experiments which also scales very differently with temperature (see SFu2009 ). Finally, we should add that the vanishing of the dephasing rate when we do not impose the infra-red cut-off is very specific to the dimensionality of space and the behaviour of the coupling $u_{k}\propto\sqrt{k}$. For instance in 1-D, we can easily check that the limit $\omega_{k}\rightarrow 0$ leads to the finite dephasing rate:

[TABLE]

For the sake of completeness, we note that the Lamb-shift term in Eq. (S30) is given by:

[TABLE]

with the Debye wavenumber $k_{D}=\omega_{D}/c_{s}$. For average emitter separations of $r_{ml}\sim\lambda_{0}/10$ in a sample of size $L\sim\lambda_{0}$, the magnitude of $S_{ml}\approx 2\pi\times 0.4$ MHz. Since this term also breaks the collective symmetry of the system, let us compare this to the dipole-dipole interaction strength $\overline{g}$. For $r_{ml}\sim\lambda_{0}/10$, $\overline{g}\sim\Gamma$, spontaneous emission rate of the NV optical transition which is at least one order of magnitude larger than $S_{ml}$. Nonetheless, it would be interesting to study systematically how this affects collective enhancement as part of future work.

Finally, we would like to point out another mechanism to generate collective dephasing. Consider a collection of TLEs that are susceptible to externally applied magnetic fields with a general hamiltonian of the form:

[TABLE]

As detailed in the appendix A of SDorner12 , if the applied field is of the white noise type, namely $B(t)=\sqrt{\frac{\gamma_{c}}{2}}\xi(t)$ with $\langle\xi(t)\rangle=0$ and $\langle\xi(t)\xi(t^{\prime})\rangle=\delta(t-t^{\prime})$ and has a large correlation length in space i.e. $p_{l}\equiv 1$, then it leads to a collective dephasing master equation of the form Eq. (4) in the main paper.