Superradiance with local phase-breaking effects

TL;DR

This paper analyzes how phase-breaking effects like non-radiative losses and dephasing influence superradiance in two-level systems, revealing that superfluorescent bursts can still occur under realistic solid-state conditions.

Contribution

It introduces a comprehensive model of superradiance considering phase-breaking mechanisms and applies it to dilute excitation regimes in solid-state systems.

Findings

Superfluorescent bursts are observable despite high dephasing rates.

Dephasing and losses affect inter-mode scattering and quasiparticle lifetimes.

Superradiance persists in solid-state systems with realistic phase-breaking effects.

Abstract

We study the superradiant evolution of a set of $N$ two-level systems spontaneously radiating under the effect of phase-breaking mechanisms. We investigate the dynamics generated by non-radiative losses and pure dephasing, and their interplay with spontaneous emission. Our results show that in the parameter region relevant to many solid-state cavity quantum electrodynamics experiments, even with a dephasing rate much faster than the radiative lifetime of a single two-level system, a sub-optimal collective superfluorescent burst is still observable. We also apply our theory to the dilute excitation regime, often used to describe optical excitations in solid-state systems. In this regime, excitations can be described in terms of bright and dark bosonic quasiparticles. We show how the effect of dephasing and losses in this regime translates into inter-mode scattering rates and quasiparticle lifetimes.

Tables

Table 1. Table 1: The table shows the derivatives of d d ​ t ​ m 𝑑 𝑑 𝑡 𝑚 \tfrac{d}{dt}{m} and d d ​ t ​ j 𝑑 𝑑 𝑡 𝑗 \tfrac{d}{dt}{j} in some characteristic points of the Dicke space (see Fig.  1 for a reference to the symbols) according to Eq. ( 11 ). In the Dicke triangle, shown in Fig.  1 , d d ​ t ​ m 𝑑 𝑑 𝑡 𝑚 \tfrac{d}{dt}{m} affects the vertical component of the system dynamics and d d ​ t ​ m 𝑑 𝑑 𝑡 𝑚 \tfrac{d}{dt}{m} the horizontal one. For each rate γ S subscript 𝛾 S \gamma_{\text{S}} , γ L subscript 𝛾 L \gamma_{\text{L}} , and γ D subscript 𝛾 D \gamma_{\text{D}} the change is given to the leading order in N 𝑁 N . Notice that only the spontaneous emission of the superradiant state | N 2 , 0 ⟩ ket 𝑁 2 0 \mathinner{|{\tfrac{N}{2},0}\rangle} is collectively enhanced ( ∝ N 2 ) proportional-to absent superscript 𝑁 2 (\propto N^{2}) .

	$○$	$\text{blacksquare}$	$\odot$	$\text{bigstar}$	$\square$
	$|\frac{N}{2},\frac{N}{2}⟩$	$|\frac{N}{4},\frac{N}{4}⟩$	$|\frac{N}{2},0⟩$	$|0,0⟩$	$|\frac{N}{4},-\frac{N}{4}⟩$
$\frac{d}{dt}m\updownarrow$	$-\left({\gamma }_{\text{S}}+{\gamma }_{\text{L}}\right)N$	$-\left(\frac{{\gamma }_{\text{S}}}{2}+\frac{3}{4}{\gamma }_{\text{L}}\right)N$	$-\frac{{\gamma }_{\text{S}}}{4}{N}^2-\frac{{\gamma }_{\text{L}}}{2}N$	$-\frac{{\gamma }_{\text{L}}}{2}N$	$-\frac{{\gamma }_{\text{L}}}{4}N$
$\frac{d}{dt}j\leftrightarrow$	$-{\gamma }_{\text{L}}N$	$\frac{{\gamma }_{\text{D}}}{2}-\frac{3}{4}{\gamma }_{\text{L}}N$	$-\left(\frac{{\gamma }_{\text{D}}}{4}+\frac{{\gamma }_{\text{L}}}{4}\right)N$	$\left(\frac{{\gamma }_{\text{D}}}{2}+{\gamma }_{\text{L}}\right)N$	$\frac{{\gamma }_{\text{D}}}{2}+\frac{{\gamma }_{\text{L}}}{4}N$

Equations

$$[J_{x,n},J_{y,n^{\prime }}]$$

$$H_0$$

$$H_{\text{int}}$$

$$J^2|j,m,\alpha ⟩=j(j+1)|j,m,\alpha ⟩,J_z|j,m,\alpha ⟩=m|j,m,\alpha ⟩,J_{\pm }|j,m,\alpha ⟩=\sqrt{(j\mp m)(j\pm m+1)}|j,m\pm 1,\alpha ⟩,$$

$$D_j$$

$$d_m$$

$$\frac{d}{dt}\rho$$

$$\frac{d}{dt}⟨J_z⟩$$

$$\frac{d}{dt}⟨J^2⟩$$

$$\frac{d}{dt}m$$

$$\frac{d}{dt}j$$

$$\frac{d}{dt}{m}_{\text{S}}$$

$$\frac{d}{dt}{j}_{\text{D}}$$

$$\frac{d}{dt}{m}_{\text{L}}$$

$$\frac{d}{dt}{j}_{\text{L}}$$

$$\frac{d}{dt}⟨J_z^2⟩$$

$$\frac{d}{dt}⟨J_{+}^p{J}_z^r{J}_{-}^q⟩$$

$$t_{\text{d}}$$

$$J_{-}$$

$$b_p^{†}$$

$$\sum _n{f}_n^p{\bar{f}}_n^q={\delta }_{p,q},f_n^0=\frac{1}{\sqrt{N}}\forall n.$$

$$J_z=\sum _p\left(b_p^{†}{b}_p-\frac{1}{2}\right),$$

$$[b_p,b_q^{†}]$$

$$H_{\text{int}}$$

$$n_b$$

$$n_d$$

$$n_b$$

$$n_d$$

$$n_b$$

$$n_d$$

$$\frac{d}{dt}{n}_b$$

$$\frac{d}{dt}{n}_d$$

Full text

I Introduction

The time a two-level system initially prepared in the excited state $\mathinner{|{\uparrow}\rangle}$ takes to relax into its ground state $\mathinner{|{\downarrow}\rangle}$ emitting a photon is set by the spontaneous emission rate $\gamma_{\text{S}}$, which quantifies its coupling with the electromagnetic environment. When $N$ identical copies of such a system are prepared in their excited state, $\mathinner{|{\uparrow...\uparrow}\rangle}$, each two-level system emits independently, leading to a total emission rate $N\gamma_{\text{S}}$. Under the condition that the linear size of the ensemble is smaller than the resonance wavelength, subsequent photon emissions then cause a build-up of quantum correlations between different systems, leading after a delay time $t_{\text{d}}=\frac{\log{N}}{N\gamma_{\text{S}}}$ Dicke (1954); Gross and Haroche (1982) to the emission of light in a short and intense burst of peak intensity $\propto N^{2}$, with a narrow pulse duration $\tau_{\text{p}}\sim\frac{1}{N\gamma_{\text{S}}}$. Such phenomenon, referred to as superfluorescence or Dicke’s superradiance, is an example of spontaneous symmetry breaking triggered by quantum fluctuations and leading to the emergence of macroscopic quantum correlations in time Vrehen et al. (1980).

Superfluorescence is due to the fully-symmetric coupling with the electromagnetic field. If the system is initially prepared in the totally-excited, maximally-symmetric state $\mathinner{|{\uparrow...\uparrow}\rangle}$, photon emission will only cause transitions within the $N+1$ symmetric states, eventually reaching the ground state $\mathinner{|{\downarrow...\downarrow}\rangle}$. The evolution of the system will thus be bound to the $(N+1)$-dimensional, symmetric subspace of the full $2^{N}$-dimensional Hilbert space. Radiative transitions, allowed only between subsequent states in the discrete energy ladder, are usually characterised by dipole matrix elements of order $\sqrt{N}$, except for the states close to half-filling, where instead the matrix element becomes of order $N$. As initially noticed by Dicke Dicke (1954), the emission thus speeds up, leading to the occurrence of a superfluorescent peak midway through the radiative cascade.

In most physical systems there are also other mechanisms competing with the light-matter interaction. While their specific sources and microscopic nature may vary, their effect on the dynamics can generally be described in terms of either nonradiative losses, or energy-conserving pure dephasing. These mechanisms are detrimental to the development of a coherent, collective dynamics, posing limitations to the observation of cooperative light emission. Their study has been the subject of several works that have considered non-ideal conditions for superfluorescence, such as inhomogeneous distribution of the two-level systems resonances, nonradiative emission, radiative emission through a local photonic environment, and effects of the light propagation in large samples that are optically thick Lehmberg (1970); Agarwal (1971); Bonifacio et al. (1971a, b); Friedberg et al. (1972); Skribanowitz et al. (1973); Bonifacio and Lugiato (1975); Ressayre and Tallet (1973); Jodoin and Mandel (1974); Ressayre and Tallet (1974); Friedberg and Hartmann (1974); Lee (1976); MacGillivray and Feld (1976); Schuurmans and Polder (1979); Leonardi and Vaglica (1982a, b); Malcuit et al. (1987); Maki et al. (1989); Drummond and Raymer (1991); Temnov and Woggon (2005); Ishikawa et al. (2016); Lambert et al. (2016). Under the continuous pumping condition, several works have also investigated the effects of non-cooperative mechanisms on the cooperative dynamics, with regard to the properties of light emission Zheleznyakov et al. (1989); Belyanin et al. (1998); Vasil’ev (1999); Temnov and Woggon (2009); Meiser and Holland (2010a, b); Bohnet et al. (2012); Pleinert et al. (2017); Sáez-Blázquez et al. (2017), and to the related phenomena of optical bistability Bonifacio and Lugiato (1978); Gronchi and Lugiato (1978); Bonifacio (1982); Carmichael et al. (1986); Sarkar and Satchell (1987a, b); Amthor et al. (2015) and superradiant phase transitions Delanty et al. (2011); Gelhausen et al. (2017); Dalla Torre et al. (2016); Kirton and Keeling (2017); Gegg et al. (2017).

Here we will exploit a theory able to describe the incoherent interaction of $N$ two-level systems with the electromagnetic field for arbitrary losses and dephasing, in order to gain a microscopic understanding of the trajectory of the system in its evolution away from symmetric states, under the effect of nonradiative losses and pure dephasing. Our results shed light on the parameters required to observe superfluorescence in solid-state environments in which the nonradiative loss rate $\gamma_{\text{L}}$ and the pure dephasing rate $\gamma_{\text{D}}$ are usually larger than the spontaneous emission rate of the single two-level system $\gamma_{\text{S}}$. In order to accomplish this, we will investigate the different dynamics generated by both phase-breaking mechanisms, and their interplay with radiative emission. We begin by introducing a master-equation description of the quantum dynamics. After testing the reliability of the theory against exact simulations of the system dynamics for small $N$, we use it to investigate a wide parameter range for large $N$. In the regime relevant for solid-state cavity quantum electrodynamics, in which collective strong coupling is observable notwithstanding dephasing effects, $N\gamma_{\text{S}}\gg\gamma_{\text{D}}\gg\gamma_{\text{S}}$, we show that the main effect of dephasing and losses is that of making superfluorescence sub-optimal, with only a limited amount of the energy released superradiantly. Nonetheless, the delay time of pure superfluorescence, $t_{\text{d}}$, remains a good predictor of the occurrence of the superfluorescent peak. Subradiance dominates the later stage of the dynamics, occurring on a timescale set by the phase-breaking mechanisms.

Finally, we will apply the developed theory to the dilute regime, in which only few excitations are present. Such a regime, very different from the half-filling relevant to study superfluorescence, is well described by a linear bosonic theory, extensively used in solid-state physics Holstein and Primakoff (1940); Ciuti et al. (2005); Agranovich (2009); De Liberato and Ciuti (2009); Kyriienko et al. (2013). Specialising our master-equation approach to this regime, we will recover a theory of interacting bright and dark bosonic quasiparticles De Liberato and Ciuti (2009); Emary (2013), in which dephasing and losses are described in terms of intermode scattering rates and particle lifetimes.

II General theory

In this section we will develop the general theory to study the evolution of a set of $N$ two-level systems evolving under the effect of different incoherent interactions with their environment. Following Dicke Dicke (1954) we begin by describing the $n$th two-level system with the algebra of a spin $\frac{1}{2}$. We thus define the operators over its two-dimensional Hilbert space $J_{\xi,n}$, $\xi\in\{x,y,z,+,-\}$, obeying the angular momentum algebra

[TABLE]

and cyclic permutations, with $J_{\pm,n}=J_{x,n}\pm iJ_{y,n}$. We also define the corresponding total spin operators $J_{\xi}=\sum_{n=1}^{N}J_{\xi,n}$. The free Hamiltonian describing $N$ identical and independent two-level systems of energy $\hbar\omega_{0}$ can then be written in terms of spins as

[TABLE]

and their interaction Hamiltonian with the electric field $\mathbf{E}$ in the dipolar approximation

[TABLE]

where $\mathbf{d}$ is the dipolar moment of each spin. From the theory of addition of angular momenta we can diagonalise the free Hamiltonian $H_{0}$ on the basis indexed by the quantum numbers $j$, $m$, and $\alpha$, where

[TABLE]

with $j=\frac{N}{2},\frac{N}{2}-1,...,j_{\text{min}}+1,j_{\text{min}}$, and $j_{\text{min}}=0,\frac{1}{2}$ for $N$ even or odd, respectively, and $m$ an integer or semi-integer, respectively, such that $|m|\leq j$. The index $\alpha$, referred to as a symmetry parameter, runs over subspaces with the same value $j$ Tavis and Cummings (1968). As clear from Eq. (4), the collective spin operators cannot change $\alpha$ and we will thus not mark it explicitly in the following. This choice introduces a degeneracy for the Dicke state $\mathinner{|{j,m}\rangle}$ that can be easily computed. Notice that there are $N$ general states with one excitation, but the ground state $\mathinner{|{\frac{N}{2},-\frac{N}{2}}\rangle}$ is connected only to a single one-excitation state, the maximally-symmetric state $\mathinner{|{\frac{N}{2},-\frac{N}{2}+1}\rangle}\propto J_{+}\mathinner{|{\frac{N}{2},-\frac{N}{2}}\rangle}$. This implies that the degeneracy of the state $\mathinner{|{\frac{N}{2}-1,-\frac{N}{2}+1}\rangle}$ is $(N-1)$. From Eq. (4) we see that each state $\mathinner{|{\frac{N}{2}-1,m}\rangle}$ has the same degeneracy, and is said to belong to the same Dicke ladder. Analogously, the degeneracy of all of the other inequivalent subspaces can be iteratively computed, leading to the formula for the degeneracy of each Dicke state $\mathinner{|{j,m}\rangle}$ to be

[TABLE]

The degeneracy for a given eigenvalue $\hbar\omega_{0}m$ of Eq. (2) instead is Mandel and Wolf (1995)

[TABLE]

The coupling with the electromagnetic field, and with the other reservoirs responsible for phase randomisation and energy loss, can be described by introducing the master equation for the density matrix $\rho$ of the system

[TABLE]

where the Lindblad operator $\mathcal{L}_{O}[\rho]=2O\rho O^{\dagger}-O^{\dagger}O\rho-\rho O^{\dagger}O,$ describes the different scattering rates quantified by the spontaneous emission rate $\gamma_{\text{S}}$, the nonradiative relaxation rate $\gamma_{\text{L}}$, and the pure dephasing rate $\gamma_{\text{D}}$. The crucial element to notice in Eq. (9) is that the electromagnetic field couples symmetrically with all the spins, and can thus be expressed in terms of the collective spin operators. This expression is not possible for the other two processes, that instead couple with each individual spin independently. In this sense dephasing and losses are both phase-breaking mechanisms, that randomise the phases between different spins. This implies that the spontaneous emission conserves $j$, leading to an evolution simply described by a radiative cascade through the $(2j+1)$ levels of one of the subspaces indexed by fixed values of $j$ and $\alpha$. This is a huge simplification, restricting the size of the relevant Hilbert space from $2^{N}$ to $(2j+1)$, allowing one to solve the problem through direct numerical diagonalisation. The other two terms instead, by randomising phases, can couple states with different $j$ and $\alpha$. The problem in this form cannot be solved exactly and brute-force numerical simulations of the dynamics cannot explore the large-$N$ limit, with a tight bound given by the growth of the Hilbert space as $2^{N}$. We will nevertheless exploit the randomisation of the phase-breaking processes by employing a numerical method that exploits the properties of permutational-invariant density matrices, which scale as ${O}(N^{3})$ in general and as ${O}(N^{2})$ for special initial states, including Dicke states Chase and Geremia (2008); Baragiola et al. (2010); Damanet et al. (2016). Similarly, the fact that the Lindblad operators are invariant under SU(4) symmetry can be employed to reduce the numerical resources required to solve the dynamics Hartmann (2016); Xu et al. (2013). Hereafter we will initially begin to gain a qualitative understanding of the effect of each one of the terms of Eq. (9), before developing solid approximations that will allow us to simulate the system dynamics.

III Dynamics analysis in the Dicke triangle

We begin by writing the evolution equations for the expectation values of $J_{z}$ and $J^{2}$. Using the shortcuts $\tfrac{d}{dt}\langle{J_{z}}\rangle=\text{Tr}\left[J_{z}\,\tfrac{d}{dt}{\rho}\,\right]$, $\tfrac{d}{dt}\langle{J^{2}}\rangle=\text{Tr}\left[J^{2}\,\tfrac{d}{dt}{\rho}\,\right]$, and $J^{2}=J^{2}_{z}-J_{z}+J_{+}J_{-}$, after some algebra we obtain

[TABLE]

We are now able to analyse the effect of the three scattering channels on an arbitrary state $\mathinner{|{j,m}\rangle}$. Projecting Eq. (10) on such a pure state we obtain

[TABLE]

We will represent the Hilbert space as done in Fig. 1, in the form of an isosceles triangle, with $m$ on the y-axis and $j$ on the x-axis. Spontaneous emission, conserving $j$ and reducing $m$ by one, can only couple a state with the one immediately below it, and indeed in Eq. (11) only $\tfrac{d}{dt}{m}$ depends on $\gamma_{\text{S}}$. In particular we obtain that the photonic emission rate from a state $\mathinner{|{j,m}\rangle}$ is

[TABLE]

The presence of many spins translates into a state-dependent superradiant enhancement of the spontaneous emission rate. From Eq. (12) we see that when $m$ is of order $j$, the squares tend to cancel, leading to a contribution at most of order $j$. Only when $m^{2}\ll j^{2}$ the $j^{2}$ term tends instead to dominate, leading to a superfluorescent emission rate of order $j^{2}$. In the literature of Dicke superradiance, $j$ is usually named cooperation number or cooperativity: it indeed measures the effective number of spins that coherently participate in the superradiant emission process, leading symmetric states to a burst of peak intensity proportional to $N^{2}$. The effect of spontaneous emission is illustrated in Fig. 1 by red vertical arrows with sizes corresponding to the emission rate. This result validates the intuitive picture of superfluorescence: after the system is prepared in the totally-excited state in the top left corner of the triangle in Fig. 1, it starts emitting photons, moving down and accelerating. When approaching $m=0$ it quickly emits most of its energy in a superfluorescent burst, then slowing down until it reaches the ground state in the bottom left corner.

The pure dephasing term in Eq. (9) on the contrary, cannot change the value of $m$ because its jump operators $J_{z,n}$ commute with $J_{z}$. It can only modify the value of the cooperativity $j$ as

[TABLE]

Its elastic, phase-randomising effect is shown as blue horizontal arrows in Fig. 1. Contrary to spontaneous emission, the effect of pure dephasing cannot lead to a superfluorescent enhancement of order $N^{2}$ due to the denominator in Eq. (13), which limits $\tfrac{d}{dt}{j}_{\text{D}}$ to be at most of order $N$. Finally, the nonradiative losses play a role on both $m$ and $j$ Lee (1976)

[TABLE]

and they are represented by green diagonal arrows in Fig. 1. Both contributions are at most of order $N$, implying again the absence of a collective enhancement to the loss rate. The previous results are summarised in Table 1 for some relevant points of the Hilbert space. The table prompts a comment on the feasibility of state preparation beyond the case of initial full excitation. Preparing a large ensemble into the maximally-entangled superradiant state $\mathinner{|{\frac{N}{2},0}\rangle}$ requires local control on the two-level systems, yet a $\tfrac{\pi}{2}-$pulse on the ground state initialises the system in a separable state whose dynamical evolution can be similarly superradiant Bonifacio et al. (1971b); Shammah (2016). Regarding the subradiant states $\mathinner{|{j,-j}\rangle}$, theoretical investigations have recently focused on the effect of a photonic cavity Gegg et al. (2017) and on the one-excitation subspace, setting thus $j=\tfrac{N}{2}-1$ and generalizing the formalism of Dicke states to the large sample regime Scully et al. (2006); Scully (2015); Vetter et al. (2016). This is the one that can be experimentally investigated in large atomic clouds Bienaimé et al. (2012). Yet enhanced state-preparation protocols and locally-controlled dynamics can be obtained in artificial atoms in solid state, especially in circuit QED systems, which have demonstrated to be a promising platform to scale up the ensemble size van Loo et al. (2013); Mlynek et al. (2014); Kakuyanagi et al. (2016).

From Eq. (10) we can easily ascertain that $\tfrac{d}{dt}{m}\leq 0$, as expected, due to the dissipative nature of the system. On the contrary, $\tfrac{d}{dt}{j}$ can be positive or negative depending on the value of both $j$ and $m$, and the ratio $\gamma_{\text{L}}/\gamma_{\text{D}}$. Indeed, the possibility of having a dephasing-driven increase of $j$ is not surprising in a finite-dimensional Hilbert space, where quantum revivals can be observed Haroche and Raimond (2006). In particular we can understand it by noticing that the ground state corresponds to $j=\frac{N}{2}$ and thus dissipative phenomena that initially tend to reduce the value of the cooperativity have eventually to increase it. In Fig. 1 we mark the boundary that separates the region of $\tfrac{d}{dt}{j}>0$ (lower-right concave area) and that of $\tfrac{d}{dt}{j}<0$ (upper-left convex area), showing that for increasing values of $\gamma_{\text{L}}/\gamma_{\text{D}}$, the region with positive derivative increases. This highlights the fact that the interplay of the two mechanisms can effectively steer the dynamics and could be used to access different points of the Dicke triangle.

From the previous arguments we can qualitatively predict how the system will behave upon initial excitation in the totally-excited state $\mathinner{|{\frac{N}{2},\frac{N}{2}}\rangle}$. If the radiative decay is the faster process, $\gamma_{\text{S}}\gg\gamma_{\text{L}},\gamma_{\text{D}}$, the system will move mainly downward, resulting in the emission of a superfluorescent burst, although at a reduced intensity with respect to pure fluorescence, as the phase-breaking processes reduce the cooperation number $j$ from its maximum value of $\frac{N}{2}$. If the opposite condition $N\gamma_{\text{S}}\ll\gamma_{\text{L}},\gamma_{\text{D}}$ is valid, then most excitations will be lost through nonradiative decay and the system will thus arrive to the ground state $\mathinner{|{\frac{N}{2},-\frac{N}{2}}\rangle}$ hopping down states lying on, or close to, the two diagonal sides of the triangle in Fig. 1. In the intermediate regime, $\gamma_{\text{S}}\ll\gamma_{\text{L}},\gamma_{\text{D}}$ but $N\gamma_{\text{S}}\gg\gamma_{\text{L}},\gamma_{\text{D}}$, the cooperation number will initially decrease as excitations are nonradiatively lost, but as the system approaches the middle of the Dicke triangle, cooperative light emission eventually kicks in leading to a superfluorescent burst.

IV Quantum dynamics through the truncated hierarchy

The system given by Eq. (10) is not closed, as it contains the average of $J_{z}^{2}$. While we could write the equation of motion for its expectation value, it will in turn depend on the evolution of products of three operators and so on, leading to a hierarchy of $O(N^{2})$ coupled equations Carmichael (1999, 2008). One usual way to deal with this kind of problem, ubiquitous in the study of nonlinear quantum systems, is to truncate the hierarchy at a certain order, factorising all the products of more than a given number of operators. The system in Eq. (11), in which we calculated the scattering rates from a specific $\mathinner{|{j,m}\rangle}$ eigenvector, effectively corresponds to a first-order approximation, $\langle J_{z}^{2}\rangle\simeq\langle J_{z}\rangle\langle J_{z}\rangle$.

With the final objective to calculate the quantum dynamics of the system under losses and dephasing, here we will check the solidity of such an approximation. In order to do this we will compare the first-order approximation of Eq. (10) both to an exact simulation for up to $N=50$, obtained by a permutational invariant method Chase and Geremia (2008) requiring only ${O}(N^{2})$ resources to model the dynamics, and to a second-order approximation for large $N$. To go to the second order, we complete the system in Eq. (10) with the equation for the expectation value of $J_{z}^{2}$

[TABLE]

and then factorise the products of three operators assuming $\langle J_{z}^{3}\rangle\simeq\langle J_{z}\rangle\langle J_{z}^{2}\rangle$ and $\langle J_{z}{J}^{2}\rangle\simeq\langle J_{z}\rangle\langle{J}^{2}\rangle$, leading to

[TABLE]

Thanks to the factorisation now Eq. (10) and Eq. (16) form a closed, nonlinear system that can be numerically solved with the proper initial conditions. With the aim of allowing an easy extension of the factorisation scheme to higher orders, here we report also the general evolution equation for any product of collective operators $J_{+}^{p}J_{z}^{r}J_{-}^{q}$ Carmichael (1999) with $p,q,r\in\mathbb{N}+\{0\}$, generated by Eq. (9)

[TABLE]

V Numerical study

We begin by assessing the validity of the approximations made in the previous Section for small $N$, as in that limit it is possible to solve the dynamics with an exact numerical diagonalisation Johansson et al. (2012, 2013). We observe that the Lindblad dynamics of Eq. (9) is permutational-invariant upon exchange of any given couple of two-level systems. The time evolution of a Dicke state (or, more generally, any permutational invariant state) will thus be restricted to the subset of other permutational-invariant density matrices, which can be represented with only ${O}(N^{3})$ terms, thanks to the fact that matrix elements corresponding to the $D_{j}$ degenerate subspaces indexed by $\alpha$ are identical. Moreover the density matrix is a block-diagonal matrix in which the coherences $\mathinner{\langle{j,m}|}\rho\mathinner{|{j^{\prime},m^{\prime}}\rangle}$ are non-zero only for $j=j^{\prime}$ Chase and Geremia (2008). Thanks to these facts, the dynamics of Eq. (9) decouples the evolution of coherences $\mathinner{\langle{j,m}|}\rho\mathinner{|{j,m^{\prime}}\rangle}$ to that of the populations $\mathinner{\langle{j,m}|}\rho\mathinner{|{j,m}\rangle}$, so that a set of rate equations for the evolution of populations of the Dicke states requires only ${O}(N^{2})$ terms Chase and Geremia (2008); Baragiola et al. (2010); Damanet et al. (2016). This is a considerable simplification, because a brute-force numerical simulation of the Liouvillian superoperator requires $(2^{2N}\times 2^{2N})$-dimensional matrices, which would become numerically challenging already for $N$ of the order of 10. By exploiting permutational invariance we can easily explore ensembles of up to $N=50$ two-level systems. We point out that there are several approaches to reduce the computational resources related to the dimension of the Lindbladian operator matrices. Permutational-invariant methods have been applied to the study of the dynamics of spin-squeezing of $N$ spin-$\frac{1}{2}$ particles Chase and Geremia (2008); Baragiola et al. (2010) and for the purpose of entanglement and state estimation Moroder et al. (2012); Novo et al. (2013). Permutational invariance has also been used to treat the more general case of multi-level systems dynamics Gegg and Richter (2016, 2017). The dimension of the Liouvillian space can be reduced to ${O}(N^{3})$ by exploiting the SU(4) symmetry of the Lindblad superoperators Hartmann (2016); Xu et al. (2013); indeed in Ref. Xu and Holland (2015) a master equation that includes the terms of Eq. (9) has been treated for the study of Ramsey spectroscopy of atomic ensembles. Recently these methods have been further developed and they have been employed to study the cooperative effects arising from the interaction with light of an ensemble of $N$ two-level systems both in the context of light emission Damanet et al. (2016); Gong et al. (2016) and that of superradiant phase transition Kirton and Keeling (2017); Gegg et al. (2017). With regard to the study of driven-dissipative steady states in open systems other tools from condensed matter theory have also been applied Sieberer et al. (2016); Dalla Torre et al. (2016); Scarlatella and Schiró (2016).

In Fig. 2(a)-(c) we show, for $N=50$, the time evolution of the normalised values of $\langle{J}_{z}(t)\rangle$, $\langle{J}^{2}(t)\rangle$, and $\langle{J}_{+}{J}_{-}(t)\rangle$ calculated exactly (blue solid curves), with the first-order approximation of Eq. (10), setting $\langle{J}_{z}^{2}\rangle=\langle{J}_{z}\rangle^{2}$ (red short-dash curves), and with the second-order approximation of Eq. (10) and Eq. (16) (black long-dash curves). We set $\gamma_{\text{S}}=1$, $\gamma_{\text{L}}=0.1$, and tune $\gamma_{\text{D}}=1$ (thin curves), 10 (medium curves), and 100 (thick curves). One can observe that qualitatively the first- and second-order approximations match well the exact dynamics, with the second-order faring better than the first-order. For $\gamma_{\text{D}}=1$, both dephasing and losses are small compared to the superradiant rate $N\gamma_{\text{S}}$. This is the condition closest to the standard case of pure superfluorescence, $\gamma_{\text{D}}=\gamma_{\text{L}}=0$, and we report it as a reference point or benchmark. In absence of dephasing and losses the peak of the superradiant pulse would occur after the delay time

[TABLE]

which is calculated by linearizing Eq. (12) Gross and Haroche (1982). As seen in Fig. 2(a) (thin curves), for $\gamma_{\text{D}}=1$ the ensemble’s energy is lost on the scale of the superradiance delay time $t_{\text{d}}$ and, as shown in Fig. 2(b), the cooperativity of the system is almost constant, marking the fact that the atomic correlations enhance the light emission intensity, which is $\propto N^{2}$, as seen in Fig. 2(c). For $\gamma_{\text{D}}=10$ (medium curves), still most of the dynamical changes are occurring on the scale of $t_{\text{d}}$, yet the cooperativity, shown in Fig. 2(b), is not constant. This is an example of sub-optimal superfluorescence, for which radiative and phase-breaking mechanisms compete ($\gamma_{\text{D}}=\gamma_{\text{S}}$), yet the collective enhancement is greater than the other processes $N\gamma_{\text{S}}\gg\gamma_{\text{D}}$. For $\gamma_{\text{D}}=100$ (thick curves) the system’s evolution is dominated by the phase-breaking processes and one retrieves that the dynamics is simply set by the characteristic time $t_{0}=(\gamma_{\text{S}}+\gamma_{\text{L}})^{-1}$ of exponential decay, as shown in Fig. 2(c).

In Fig. 2(d)-(f) we compare the results for the time evolutions for large $N$, with $N=10^{2}$ (thick curves), $N=10^{3}$ (medium curves), and $N=10^{4}$ (thin curves). We can access the large-$N$ limit thanks to the first-order (red short-dash curves) and second-order (black long-dash curves) approximations deriving from Eq. (10), which become almost perfectly overlapped. This is not surprising considering that the factorisation schemes, as mean field approximations, are expected to improve with increasing $N$. Since in the solid state, dephasing is generally the fastest mechanism, we have set now radiative decay as the weakest one, $\gamma_{\text{S}}=1$, $\gamma_{\text{L}}=10$, and $\gamma_{\text{D}}=100$. Notice that for $N=10^{2}$ (thin curves) we have $\gamma_{\text{D}}=N\gamma_{\text{S}}$, the dynamics is completely incoherent and the decay is exponential with characteristic time $t_{0}$. As $N$ is increased to $N=10^{3}$ and $N=10^{4}$ (medium and thick curves), the dominant dynamic drive is that of collective decay, since $N\gamma_{\text{S}}\gg\gamma_{\text{D}},\gamma_{\text{L}}$ even if the system is still in the limit of strong dephasing, $\gamma_{\text{D}}\gg\gamma_{\text{L}}\gg\gamma_{\text{S}}$.

We are now in a position to study numerically the dynamics for large $N$ and investigating how superfluorescence is affected by losses and dephasing. We can define operatively a superfluorescence effective delay time $t_{\text{d}}^{\text{eff}}$ in the presence of phase-breaking mechanisms as the time at which the system reaches half-filling, $\langle J_{z}(t_{\text{d}}^{\text{eff}})\rangle=0$. In the limiting case in which $\gamma_{\text{D}}=\gamma_{\text{L}}=0$, $t_{\text{d}}^{\text{eff}}$ corresponds to the pure superfluorescence delay time $t_{\text{d}}$ of Eq. (18), while when the incoherent mechanisms dominate the dynamics, this is simply the characteristic time $t_{0}$ of the exponential decay.

In Fig. 3(a) we set $\gamma_{\text{S}}=1$, $\gamma_{\text{L}}=10$ and report the value of $t_{\text{d}}^{\text{eff}}$ calculated from the closed set of Eq. (10) and Eq. (16) for varying $N$ and $\gamma_{\text{D}}$. The threshold between the exponential decay (white background) and the cooperative one (dark shading) is clear, occurring when $N\gamma_{\text{S}}\gg\gamma_{\text{D}}$. This is consistent with previous analyses that estimated the threshold as $T_{2}^{*}=\sqrt{\tau_{\text{p}}t_{\text{d}}}$, where $T_{2}^{*}$ is the critical pure dephasing time, $\tau_{\text{p}}$ the pulse duration Leonardi and Vaglica (1982b); Malcuit et al. (1987); Maki et al. (1989); Ishikawa et al. (2016), i.e. $\gamma_{\text{D}}^{*}\simeq\frac{\gamma_{\text{S}}N}{\sqrt{\ln{N}}}$, plotted as a black dashed curve in Fig. 3(a). In order to better understand the global dynamics of the system, we can now perform a qualitative analysis of a specific point in the parameter space of Fig. 3(a), choosing one for which superfluorescence occurs but it is relatively near to the threshold. We consider $N=10^{3}$ and $\gamma_{\text{D}}=100$ (red dot in panel (a)), and recall that $\gamma_{\text{S}}=1$ and $\gamma_{\text{L}}=10$. This means that this point is in the regime $N\gamma_{\text{S}}>\gamma_{\text{D}}\gg\gamma_{\text{L}}>\gamma_{\text{S}}$, where we expect superfluorescence to occur. In Fig. 3(b) we show for this choice of parameters the trajectory of the system in the $(j,m)$ Dicke space (arrows guide the eye to indicate the time evolution). The system does not evolve along the left side of the triangular phase space down the longest Dicke ladder for which $j=\tfrac{N}{2}$, as in pure superfluorescence; instead it explores the inner area of the Dicke space until it reaches a sub-optimal dark state $\mathinner{|{j,-j}\rangle}$ on the lower right side of the triangle. From then on, as already analysed qualitatively in the previous section, see Table 1 for the state $\mathinner{|{\tfrac{N}{4},-\tfrac{N}{4}}\rangle}$, the dynamics is governed by nonradiative energy loss. The inset of Fig. 3(b) further shows that the two different stages of the dynamics are a first fast one, and a second one, given by the exponential nonradiative decay processes, of order of $t_{0}$. The effect of this two-stage dynamics for the controlled generation of subradiant states has also been analysed in Ref. Gegg et al. (2017) in the regime of small dephasing. A color plot of the Dicke space shows the intensity of light emission enhancement in a given point ($j,m$), through $\tfrac{d}{dt}{m}_{\text{S}}$ from Eq. (12) expressed in units of $N^{2}\gamma_{\text{S}}$, from which it is clear that although the system does not undergo pure superfluorescence, it crosses the inner superfluorescent area and light emission becomes superradiantly enhanced.

VI Bosonisation in the dilute regime

Having elucidated the impact of losses and dephasing on superfluorescence, here we apply the theory we developed to the lower left corner of the Dicke triangle, characterised by average numbers of excited spins $(m+\frac{N}{2})$ much smaller than their total number $N$. In this regime the system’s excitations can be treated as approximately bosonic. Such an approach, largely applied in solid-state physics Ciuti et al. (2005); Agranovich (2009); De Liberato and Ciuti (2009); Kyriienko et al. (2013) is based on the general idea that a system far from saturation is essentially harmonic. It can be formalised through the Holstein-Primakoff transformation that exactly maps the algebra of a spin $j$ into the one of an harmonic oscillator Holstein and Primakoff (1940); Ressayre and Tallet (1975); Lambert et al. (2004); Kessler et al. (2012)

[TABLE]

with $\left[b,b^{\dagger}\right]=1$. The number of bosonic excitations in a generic state $\mathinner{\langle{j,m}|}b^{\dagger}b\mathinner{|{j,m}\rangle}$ is thus equal to $(j+m)$, that is the number of excitations over the lowest lying state for a fixed value of $j$. If such a number is much smaller than $j$ (dilute regime) we can then do a lowest-order expansion of the square roots in the small parameter $b^{\dagger}b/2j$, recovering standard quadratic, bosonic Hamiltonians. In our case this formal approach is not directly applicable because the mapping explicitly depends on $j$ and phase-breaking effects would thus have an impact on the very definition of the bosonic operators Dalla Torre et al. (2016). In order to understand how the bosonic approach maps into the full picture of a spin Hamiltonian, we will instead proceed more phenomenologically, studying the structure of the lower left corner of the Hilbert space in Fig. 4, corresponding to the dilute regime. We begin by formally defining the operators

[TABLE]

where the $f$s form an orthonormal basis

[TABLE]

Using Eq. (20) we can then calculate

[TABLE]

which stands as a generalisation of the third term of Eq. (19). The operators of Eq. (20) are defined such that the $N$ states $b^{\dagger}_{p}\mathinner{|{\frac{N}{2},-\frac{N}{2}}\rangle}$ span the one excitation subspace. The bosonic approximation consists in assuming that those modes remain orthogonal and bosonic also in the excited manifolds, that is

[TABLE]

Using Eq. (21), we can rewrite Eq. (3) in the form

[TABLE]

from which we see that only the $p=0$ (bright) mode couples with light, with a dipole $\sqrt{N}$ times the bare one Kakuyanagi et al. (2016), while the other $N-1$ modes are dark. Notice that the collective $\sqrt{N}$ enhancement to the dipole has an impact on the system’s spectrum, allowing us to explore non-perturbative coupling regimes (strong Kavokin et al. (2007), ultra Ciuti et al. (2005), or deep De Liberato (2014)), but it does not affect the emission rate. Since the emission rate is proportional to the square of the dipole, we recover a total emission that scales as $N$, as if each dipole was emitting independently. This is as expected from the discussion of the previous section, in which we saw that an enhancement of the emission rate can only be obtained for a system close to half-filling $\lvert m\lvert\ll j$, where the bosonic approximation does not hold.

In order to investigate the role of losses and dephasing on the dynamics of bright and dark modes in the bosonic approximation, in Fig. 4 we plot a zoom of the lowest left corner of the Hilbert space, in which the bosonic approximation is relevant, explicitly showing the degeneracy of each state, $D_{j}$ Mandel and Wolf (1995). In the one excitation manifold, which is exactly spanned by the $N$ vectors $b^{\dagger}_{p}\mathinner{|{\frac{N}{2},-\frac{N}{2}}\rangle}$, the effect of both spontaneous emission and of nonradiative losses is trivial: as both need to lower the number of excitations and there is a single ground state, they will scatter each state into the ground state $\mathinner{|{\frac{N}{2},-\frac{N}{2}}\rangle}$. From Eq. (13) we can see that the pure dephasing term will lead, to the dominant order in $N$, to $\tfrac{d}{dt}{j}_{\text{D}}=-\gamma_{\text{D}}$ for the bright state and to $\tfrac{d}{dt}{j}_{\text{D}}=\frac{1}{N}\gamma_{\text{D}}$ for the dark ones. This can be easily interpreted recalling that the dephasing randomises the phases between the different spins, thus transforming one mode into the other. Given that there is a single bright mode and $(N-1)$ dark ones, any phase change will transform a bright mode into a dark one, decreasing $j$, but the majority of the phase changes will not influence the population of the dark modes. This describes well, for example, intersubband transitions in doped quantum wells, where the line width of the bright mode coupled with the electromagnetic field is determined by its dephasing rate, that transforms it into a dark, uncoupled excitation that eventually relaxes nonradiatively Liu and Capasso (1999); Shammah and De Liberato (2015).

This picture can be extended to the higher excited manifolds taking into consideration that the addition of a dark excitation in the dilute regime effectively translates only into an upward-right diagonal step, as shown in Fig. 4, reducing the value of $j$ by 1 while increasing $m$, due to the greater weight of the degeneracy $D_{j}$ for the Dicke ladders of lower cooperation number $j$, for $j=\tfrac{N}{2}$. This can be verified by counting the available modes: there are ${O}(N)$ dark modes and a single bright one, so that there are ${O}(N^{2})$ states with two dark excitations and ${O}(N)$ states with a dark and a bright excitation and analogously for higher excitation numbers. We are using orders of magnitude here instead of the exact mode counting because considering $N$ independent bosonic modes leads to over-counting the number of states, with an error of the order of the total number of excitations divided $N$. Once again we see that the bosonic approximation breaks down when the number of excitations becomes comparable with the total number of spins.

In order to describe the effect of phase-breaking mechanisms on the multi-excitation dynamics in the bosonic approximation, we start by introducing the bright and dark excitation populations

[TABLE]

Using Eq. (20) and Eq. (21) in Eq. (25) we obtain

[TABLE]

that in the dilute regime, to the dominant order in $N$ leads to

[TABLE]

The physical content of Eq. (27) can be immediately verified, as the total number of excitations is clearly given by $n_{b}+n_{d}=m+\tfrac{N}{2}$, whereas the number of dark excitations increases by 1 when reducing $j$ by 1.

We can now capture the dynamics of the populations of bright and dark states in Eq. (25) either directly from Eq. (26) using Eq. (10) or more formally by a derivation in time of Eq. (25) using Eq. (20) and Eq. (22), as both give to the dominant order in $N$

[TABLE]

Such a system is consistent with our initial interpretation: All the scattering processes destroy the bright mode, while the dark modes can decay only through nonradiative losses and they are replenished by dephasing that transforms a bright mode into a dark one. Note that in the second line of Eq. (28), the absence of a term proportional to $\gamma_{\text{L}}n_{b}$ that one would expect from the previous discussion, as the phase-randomising effect of the nonradiative losses, should transform bright into dark modes. This effect is absent because from Eq. (14) we see that in the symmetric-states subspace $j=\frac{N}{2}$, $n_{b}=-\frac{N}{2}+m\ll N$, $\tfrac{d}{dt}{j}_{\text{L}}$ is of order $\frac{1}{N}$. That is, in the dilute excitation regime nonradiative losses do not cause dephasing. An equivalent result could be obtained also by neglecting the jump term $\tfrac{\gamma_{\text{L}}}{2}\sum_{n}J_{-,n}\rho J_{+,n}$ in the Lindblad superoperator $\mathcal{L}_{J_{-,n}}[\rho]$ in Eq. (9), a procedure that is justified in the case of weak pumping, which falls into the dilute regime treated here Visser and Nienhuis (1995); Brecha et al. (1999); Sáez-Blázquez et al. (2017).

VII Conclusions

We have performed a study of the interplay of superradiant light emission and local phase-breaking mechanisms. Describing dephasing and nonradiative decay with a master equation formalism, solved through a mix of approximate and exact approaches, we investigated how the interplay of the different scattering channels influences the characteristics, and even the occurrence of the superfluorescent burst. Our results set clear requirements to observe superfluorescence in different condensed-matter systems. We also show how our treatment, when applied to the dilute excitation regime, allows one to describe phase-breaking mechanisms in term of scattering of bright and dark modes.

VIII Acknowledgements

We thank Adam Miranowicz, Anton Frisk Kockum, Carlos Sánchez Muñoz, François Damanet, and Peter Kirton for useful discussions and comments. S.D.L. acknowledges support from EPSRC Grant EP/M003183/1. S.D.L. is a Royal Society Research Fellow. F.N. was partially supported by the RIKEN iTHES Project, MURI Center for Dynamic Magneto-Optics via the AFOSR Award No. FA9550-14-1-0040, the Japan Society for the Promotion of Science (KAKENHI), the IMPACT program of JST, and CREST grant No. JPMJCR1676. F. N. and N.L. acknowledge support from the RIKEN-AIST Joint Research Fund. F. N., N.S., and N.L. acknowledge support from the Sir John Templeton Foundation. N.S. acknowledges support from EPSRC Grant EP/L020335/1.