Cooperative spontaneous emission via renormalization approach: Classical versus semi-classical effects

TL;DR

This paper develops a renormalized perturbation theory to analyze cooperative spontaneous emission in large atomic ensembles, distinguishing classical from semi-classical effects and enabling efficient simulation of complex quantum dynamics.

Contribution

It introduces an analytical solution for coupled-dipole equations that captures classical and semi-classical emission differences and allows scalable computation of large atom systems.

Findings

Reduction of radiated power cooperativity at higher excitation levels

Accurate second-order emission predictions over multiple lifetimes

Efficient simulation of hundred-thousand-atom dynamics

Abstract

We address the many-atom emission of a dilute cloud of two-level atoms through a renormalized perturbation theory. An analytical solution for the truncated coupled-dipole equations is derived, which contains an effective spectrum associated to the initial conditions. Our solution is able to distinguish precisely classical from semi-classical predictions for large atomic ensembles. This manifests as a reduction of the cooperativity in the radiated power for higher atomic excitation, in disagreement with the fully classical prediction from linear optics. Moreover, the second-order cooperative emission appears accurate over several single atom lifetimes and for interacting regimes stronger than those permitted in conventional perturbation theory. We can compute the semiclassical dynamics of hundred thousand of interacting atoms with ordinary computational resources, which makes our formalism particularly promising to probe the nonlinear dynamics of quantum many-body systems that emerge from cumulant expansions.

Figures

Equations

$$P\left(t\right)\propto\int_{0}^{2\pi}\int_{0}^{\pi}\Bigl{\langle}\hat{E}^{\dagger}\left(\mathbf{r},t\right)\hat{E}\left(\mathbf{r},t\right)\Bigr{\rangle}r^{2}d\Omega,$$

$$\hat{E}\left(\hat{\mathbf{n}},t\right)\propto \frac{1}{r}\sum _{m=1}^N{e}^{-i{k}_0\hat{\mathbf{n}}\cdot {\mathbf{r}}_m}{\hat{\sigma }}_m^{-}\left(t\right),$$

$$P_{\mathrm{sc}}\propto\sum_{m=1}^{N}\biggl{[}\frac{1+\bigl{\langle}\hat{\sigma}_{m}^{z}\bigr{\rangle}}{2}+\sum_{n\neq m}\Re\left(u_{mn}\bigl{\langle}\hat{\sigma}_{m}^{+}\bigr{\rangle}\bigl{\langle}\hat{\sigma}_{n}^{-}\bigr{\rangle}\right)\biggr{]},$$

$$\displaystyle\overset{.}{\bigl{\langle}\hat{\sigma}_{m}^{-}\bigr{\rangle}}$$

$$\displaystyle\overset{.}{\bigl{\langle}\hat{\sigma}_{m}^{z}\bigr{\rangle}}$$

$$\left(\bigl{\langle}\hat{\sigma}_{m}^{-}\left(0\right)\bigr{\rangle},\bigl{\langle}\hat{\sigma}_{m}^{z}\left(0\right)\bigr{\rangle}\right)=\left(e^{-i\mathbf{k}_{0}\cdot\mathbf{r}_{m}}\beta,2p_{e}-1\right),$$

$${\dot{x}}_l=-V_l-M_{lm}{x}_m-T_{lmn}{x}_m{x}_n,$$

$$\begin{gathered} x_l(t)={\int }_0^{\infty }d{t}^{\prime }{G}_{lm}(t,t^{\prime })V_m-G_{lm}(t,0)x_m(0) \\ +{\int }_0^{\infty }d{t}^{\prime }{G}_{lm}\left(t,t^{\prime }\right)T_{mno}{x}_n\left(t^{\prime }\right)x_o\left(t^{\prime }\right), \end{gathered}$$

$$\frac{d{G}_{lm}}{dt}+M_{ln}G{}_{nm}=-\delta \left(t-t^{\prime }\right){\delta }_{lm}.$$

$$G_{lm}\left(t-t^{\prime }\right)=-\sum _{\alpha }{\psi }_l^{\left(\alpha \right)}{\varphi }_m^{\left(\alpha \right)}\theta \left(t-t^{\prime }\right)e^{-{\gamma }^{\left(\alpha \right)}\left(t-t^{\prime }\right)},$$

$${\varphi }_l^{\left(\alpha \right)}\equiv \sum _{\beta }{\left[{\mathbf{J}}^{-1}\right]}^{\left(\alpha \beta \right)}{\psi }_l^{\left(\beta \right)},$$

$$\begin{gathered} x_l(t)=x_l^{(0)}(t)-\sum _{\alpha }{\psi }_l^{\left(\alpha \right)}{\varphi }_m^{\left(\alpha \right)}{T}_{mno}{e}^{-{\gamma }^{\left(\alpha \right)}t} \\ \times {\int }_0^{t}d{t}^{\prime }{x}_n\left(t^{\prime }\right)x_o\left(t^{\prime }\right)e^{{\gamma }^{\left(\alpha \right)}{t}^{\prime }}, \end{gathered}$$

$$x_l^{(0)}(t)\equiv \sum _{\alpha }{\psi }_l^{\left(\alpha \right)}\left(b^{\left(\alpha \right)}{e}^{-{\gamma }^{\left(\alpha \right)}t}-a^{\left(\alpha \right)}\right),$$

$$x_{l}^{\left(1\right)}\left(t\right)=-\sum_{\alpha}\psi_{l}^{\left(\alpha\right)}e^{-\gamma^{\left(\alpha\right)}t}\Biggl{[}A^{\left(\alpha\right)}\int_{0}^{t}dt^{\prime}e^{\gamma^{\left(\alpha\right)}t^{\prime}}\\ -\sum_{\beta}B^{\left(\alpha\beta\right)}\int_{0}^{t}dt^{\prime}e^{\left(\gamma^{\left(\alpha\right)}-\gamma^{\left(\beta\right)}\right)t^{\prime}}\\ +\sum_{\beta,\delta}C^{\left(\alpha\beta\delta\right)}\int_{0}^{t}dt^{\prime}e^{\left(\gamma^{\left(\alpha\right)}-\gamma^{\left(\beta\right)}-\gamma^{\left(\delta\right)}\right)t^{\prime}}\Biggr{]}.$$

$$A^{\left(\alpha \right)}$$

$$B^{\left(\alpha \beta \right)}$$

$$C^{\left(\alpha \beta \delta \right)}$$

$$x_{l}\left(t\right)\approx\sum_{\alpha}\psi_{l}^{\left(\alpha\right)}\Biggl{[}\left(b^{\left(\alpha\right)}+tS^{\left(\alpha\right)}\right)e^{-\gamma^{\left(\alpha\right)}t}-a^{\left(\alpha\right)}\\ -A^{\left(\alpha\right)}\frac{1-e^{-\gamma^{\left(\alpha\right)}t}}{\gamma^{\left(\alpha\right)}}-\overset{{\scriptscriptstyle\left(\gamma^{\left(\alpha\right)}\neq\gamma^{\left(\beta\right)}\right)}}{\sum_{\beta}}B^{\left(\alpha\beta\right)}\frac{e^{-\gamma^{\left(\alpha\right)}t}-e^{-\gamma^{\left(\beta\right)}t}}{\gamma^{\left(\alpha\right)}-\gamma^{\left(\beta\right)}}\\ +\overset{{\scriptscriptstyle\left(\gamma^{\left(\alpha\right)}\neq\gamma^{\left(\beta\right)}+\gamma^{\left(\delta\right)}\right)}}{\sum_{\beta,\delta}\quad C^{\left(\alpha\beta\delta\right)}}\frac{e^{-\gamma^{\left(\alpha\right)}t}-e^{-\left(\gamma^{\left(\beta\right)}+\gamma^{\left(\delta\right)}\right)t}}{\gamma^{\left(\alpha\right)}-\gamma^{\left(\beta\right)}-\gamma^{\left(\delta\right)}}\Biggr{]},$$

$$S^{\left(\alpha \right)}\equiv \stackrel{\left({\gamma }^{\left(\alpha \right)}={\gamma }^{\left(\beta \right)}\right)}{\sum _{\beta }}{B}^{\left(\alpha \beta \right)}-\stackrel{\left({\gamma }^{\left(\alpha \right)}={\gamma }^{\left(\beta \right)}+{\gamma }^{\left(\delta \right)}\right)}{\sum _{\beta ,\delta }{C}^{\left(\alpha \beta \delta \right)}}.$$

$$e^{-{\gamma }^{\left(\alpha \right)}t}=e^{-{\bar{\gamma }}^{\left(\alpha \right)}t}\left(1-\delta {\gamma }^{\left(\alpha \right)}t+\cdots \right),$$

$$x_{l}\left(t\right)\approx\sum_{\alpha}\psi_{l}^{\left(\alpha\right)}\Biggl{[}b^{\left(\alpha\right)}e^{-\bar{\gamma}^{\left(\alpha\right)}t}-a^{\left(\alpha\right)}-A^{\left(\alpha\right)}\frac{1-e^{-\bar{\gamma}^{\left(\alpha\right)}t}}{\gamma^{\left(\alpha\right)}}\\ -\overset{{\scriptscriptstyle\left(\gamma^{\left(\alpha\right)}\neq\gamma^{\left(\beta\right)}\right)}}{\sum_{\beta}}B^{\left(\alpha\beta\right)}\frac{e^{-\bar{\gamma}^{\left(\alpha\right)}t}-e^{-\bar{\gamma}^{\left(\beta\right)}t}}{\gamma^{\left(\alpha\right)}-\gamma^{\left(\beta\right)}}\\ +\overset{{\scriptscriptstyle\left(\gamma^{\left(\alpha\right)}\neq\gamma^{\left(\beta\right)}+\gamma^{\left(\delta\right)}\right)}}{\sum_{\beta,\delta}\quad C^{\left(\alpha\beta\delta\right)}}\frac{e^{-\bar{\gamma}^{\left(\alpha\right)}t}-e^{-\left(\bar{\gamma}^{\left(\beta\right)}+\bar{\gamma}^{\left(\delta\right)}\right)t}}{\gamma^{\left(\alpha\right)}-\gamma^{\left(\beta\right)}-\gamma^{\left(\delta\right)}}\Biggr{]}.$$

$$\displaystyle\bigl{\langle}\sigma_{l}^{-}\left(t\right)\bigr{\rangle}$$

$$\displaystyle\bigl{\langle}\sigma_{l}^{z}\left(t\right)\bigr{\rangle}$$

$${\bar{\gamma }}_{\pm }^{\left(l\right)}$$

$${\bar{\gamma }}_z^{\left(l\right)}$$

$${\eta }_l\equiv \sum _{m\ne l}{u}_{lm}{e}^{i{\mathbf{k}}_0\cdot ({\mathbf{r}}_l-{\mathbf{r}}_m)}.$$

$$\delta \left(t-t^{\prime }\right)\left({\delta }_{lm}-\sum _{\alpha }{\psi }_l^{\left(\alpha \right)}{\varphi }_m^{\left(\alpha \right)}{e}^{-{\gamma }^{\left(\alpha \right)}\left(t-t^{\prime }\right)}\right)=0.$$

$$\sum _{\alpha }{\psi }_l^{\left(\alpha \right)}{\varphi }_m^{\left(\alpha \right)}={\delta }_{lm}.$$

$$\sum _{\alpha }{\psi }_l^{\left(\alpha \right)}{\varphi }_m^{\left(\alpha \right)}{\psi }_m^{\left(\gamma \right)}={\delta }_{lm}{\psi }_m^{\left(\gamma \right)}.$$

$$\sum _{\alpha ,\beta }{\psi }_l^{\left(\alpha \right)}{\left[{\mathbf{J}}^{-1}\right]}^{\left(\alpha \beta \right)}{\psi }_m^{\left(\beta \right)}{\psi }_m^{\left(\gamma \right)}$$

$$\sum _{\alpha ,\beta }{\psi }_l^{\left(\alpha \right)}{\left[{\mathbf{J}}^{-1}\right]}^{\left(\alpha \beta \right)}{\left[\mathbf{J}\right]}^{\left(\beta \gamma \right)}$$

$$\sum _{\alpha }{\psi }_l^{\left(\alpha \right)}{\delta }_{\alpha \gamma }$$

$${\psi }_l^{\left(\gamma \right)}$$

$$C\equiv 2\frac{\left|\sum_{m=1}^{N}\sum_{n\neq m}\Re\left(u_{mn}\bigl{\langle}\hat{\sigma}_{m}^{+}\left(0\right)\bigr{\rangle}\bigl{\langle}\hat{\sigma}_{n}^{-}\left(0\right)\bigr{\rangle}\right)\right|}{\left|\sum_{m=1}^{N}\left(1+\bigl{\langle}\hat{\sigma}_{m}^{z}\left(0\right)\bigr{\rangle}\right)\right|}.$$

Full text

Non-equilibrium collective phenomena can manifest in classical and quantum many-body systems. A particularly challenging task is the segregation of those with purely non-classical nature. Indeed, the description of many-body systems through a fully quantum model suffers from the exponential growth of the number of quantum correlations with the number of particles. Hence, exact numerical methods, such as diagonalization, are limited to a few tens of particles, even for supercomputers. This is the reason why truncation methods were developed, based on cumulant expansions Krämer, Sebastian and Ritsch, Helmut (2015); Schachenmayer et al. (2015a, b); Pucci et al. (2016, 2017); Piñeiro Orioli et al. (2017), in which higher-order connected correlations are neglected to obtain a more favorable scaling of the complexity with the number of particles. As a consequence, nonlinearities emerge from the decomposition of higher-order correlations into products of lower-order ones, which turns even trickier predictions on collective quantum phenomena.

At the lowest order of the cumulant expansion one can find the “semi-classical model”, sometimes called “mean-field approximation”. It can be derived by neglecting two-body and higher-order connected correlations, keeping the contribution from “single-body” auto-correlations. As a result, a set of nonlinear coupled equations is obtained only for one-operator expected values, that scales linearly with the number of particles Zhu et al. (2016). Even in this simpler case, numerical solutions for real physical systems demand long running times, due to the nonparallelizable nature of the numerical algorithms. Thus, approximated analytical solutions which allow for independent calculations might be the unique possible alternative.

In this letter, we propose to treat analytically the many-atom spontaneous emission, in the semi-classical approximation, using a renormalization method in perturbation theory. This method was explored in the nineteenth century when it was discovered that naive perturbation theory was giving rise to unphysical behavior in celestial mechanics problems Lindstedt (1883). A proposed solution to this plight arose from the so-called Poincaré-Lindstedt method Bender and A. Orszag (1999), also known as renormalization method. It consists in eliminating unphysical secular terms from perturbative solutions, by renormalizing the characteristic eigenvalues of the differential equations. As a consequence, great advantages can be obtained over conventional power series solutions. More recently, this approach has been used with variational methods Pelster et al. (2003) and applied to dynamics of Bose-Einstein condensates Vidanović et al. (2011); Al-Jibbouri et al. (2013).

In general, cooperative spontaneous emission presents all features of a truly quantum many-body problem. Yet the accelerated emission of radiation, known as superradiance, was understood to be essentially a semi-classical cooperative effect, as the radiative decay of a nearly fully inverted system is well described by the Bloch-Maxwell equations Arecchi et al. (1972); MacGillivray and Feld (1976). Only almost four decades later, a similar phenomenon was revealed in the low excitation regime Scully et al. (2006); Mazets and Kurizki (2007), described by a fully-classical Coupled-Dipole Model (CDM) Svidzinsky et al. (2008); Bux et al. (2010); Svidzinsky et al. (2010); Bienaimé et al. (2011, 2013); Svidzinsky et al. (2015); Cottier et al. (2018), with experimental verifications provided in the linear-optics regime Bux et al. (2010); de Oliveira et al. (2014); Araújo et al. (2016); Roof et al. (2016). Long-lived radiation modes, know as subradiance, were equally well captured by the CDM Eberly (2006); Bienaimé et al. (2012); Guerin et al. (2016); Cottier et al. (2018). Thus, the divergence between classical and semi-classical cooperative effects represents an opportunity to demonstrate the efficiency of the renormalization method.

A natural consequence of eigenvalues renormalization is the rising of an analytical effective eigenspectrum, with strong dependence on initial conditions. We demonstrate that the average many-atom decay rate of this spectrum scales with an excitation-dependent cooperativity parameter, which is a modification of the linear-optics optical thickness. As a result, super- and subradiance are substantially reduced when the excited state population of the atoms is increased. While such cooperativity suppression cannot be predicted by the classical model, it is accurately captured by a first-order renormalized solution for the semi-classical equations.

Let us first introduce classical and semi-classical models which describe the spontaneous emission from $N$ two-level atoms . This emission process can be monitored by the total power radiated throughout a sphere of radius $r$,

[TABLE]

with $d\Omega$ an infinitesimal solid angle. In the far-field limit, the scalar electric field operator $\hat{E}\left(\hat{\mathbf{n}},t\right)$, at the observation point $\hat{\mathbf{n}}$, relates to the radiating dipoles as

[TABLE]

where $\mathbf{r}_{m}$ are the constant center-of-mass position of the atoms, randomly distributed in the three-dimensional space. \hat{\sigma}_{m}^{-}=\bigl{|}g_{m}\bigr{\rangle}\bigl{\langle}e_{m}\bigr{|} and $\hat{\sigma}_{m}^{+}=\left(\hat{\sigma}_{m}^{-}\right)^{\dagger}$ correspond, respectively, to the lowering and rising dipole operators, with \bigl{|}g_{m}\bigr{\rangle} ( \bigl{|}e_{m}\bigr{\rangle}) the ground (excited) state. For the derivation of Eq.(2), one considers the standard atom-light interaction, with a continuous integration over all vacuum modes of the electromagnetic field Bienaimé et al. (2011, 2013), together with the Markov and rotating wave approximations. We assume that the atomic cloud is previously prepared by a laser field, with wave vector $\mathbf{k}_{0}=k_{0}\hat{\mathbf{z}}$ and wavelength $\lambda=2\pi/k_{0}$, detuned by $\Delta$ from the atomic resonance. We then observe the atomic dipoles spontaneously decaying after switching off the pump.

By substituting Eq.(2) in Eq.(1), the radiated power can be expressed as a sum over two-atom quantum correlators. Therefore, a solution for the set of $2^{N}$ coupled differential equations, that contains $N$-body correlations between atoms, is required. In order to prevent such exponential growth, one typically truncates the system of equations up to some order in connected correlations Krämer and Ritsch (2015); Young et al. (2018). In this work, we neglect the correlations between the quantum fluctuations of different atomic dipoles (\bigl{\langle}\hat{\sigma}_{m}^{\alpha}\hat{\sigma}_{n}^{\beta}\bigr{\rangle}\approx\bigl{\langle}\hat{\sigma}_{m}^{\alpha}\bigr{\rangle}\bigl{\langle}\hat{\sigma}_{n}^{\beta}\bigr{\rangle}, for $m\neq n$), keeping the single-body correlations of equal indexes (\bigl{\langle}\hat{\sigma}_{m}^{+}\hat{\sigma}_{m}^{-}\bigr{\rangle}=\left(1+\bigl{\langle}\hat{\sigma}_{m}^{z}\bigr{\rangle}\right)/2), with \hat{\sigma}_{m}^{z}=\bigl{|}e_{m}\bigr{\rangle}\bigl{\langle}e_{m}\bigr{|}-\bigl{|}g_{m}\bigr{\rangle}\bigl{\langle}g_{m}\bigr{|} the inversion population operator. This is the origin of the name “semi-classical”.

Using this approximation, the total power (1) reduces to

[TABLE]

where the single-atom contributions have been explicitly separated from the second sum. The time evolution of the total power is then completely determined by solving the following set of nonlinear equations:

[TABLE]

whose number of variables increases linearly with $N$. For the derivation of Eqs.(3-4), the same conditions used to derive Eq.(2) are assumed Zhu et al. (2016); Pucci et al. (2017). In particular, the single atom decay rate $\Gamma$ and the light-mediated long-range interactions u_{mn}\equiv\exp\left(ik_{0}\bigl{|}\mathbf{r}{}_{m}-\mathbf{r}{}_{n}\bigr{|}\right)/ik_{0}\bigl{|}\mathbf{r}{}_{m}-\mathbf{r}{}_{n}\bigr{|} arise, by first principles, from the continuous integration over all vacuum modes of the electromagnetic field. The effective interactions $u_{mn}$ are valid in the optically dilute regime, where near-field terms do not affect the scattering.

In the low excitation regime, all the atoms evolve close to the ground state (\bigl{\langle}\hat{\sigma}_{m}^{z}\bigr{\rangle}\approx-1), and even the single-atom correlations separate: \bigl{\langle}\hat{\sigma}_{m}^{+}\hat{\sigma}_{m}^{-}\bigr{\rangle}\approx\bigl{\langle}\hat{\sigma}_{m}^{+}\bigr{\rangle}\bigl{\langle}\hat{\sigma}_{m}^{-}\bigr{\rangle}. In such case, Eqs.(3-4) reduce even further, giving rise to the linear-optics CDM, that can be derived from the classical Maxwell equations Svidzinsky et al. (2010); Cottier et al. (2018).

Regarding the initial conditions, we hereafter assume a mean-field-like initial state,

[TABLE]

where all atoms have an equal probability $p_{e}$ of being excited by the laser. However, different phases $\mathbf{k}_{0}\cdot\mathbf{r}_{m}$ are kept for the coherences of the atoms, which are proportional to the real amplitude $\beta=\sqrt{p_{e}\left(1-p_{e}\right)}$. In the low excitation regime ($p_{e}\rightarrow 0$), the classical case is recovered ($\beta\approx\sqrt{p_{e}}$). Oppositely, in the strong excitation limit ($p_{e}\rightarrow 1$), the coherence vanishes ($\beta\rightarrow 0$).

In order to extract an approximate solution from Eq.(4), let us consider a broader class of inhomogeneous nonlinear differential equations,

[TABLE]

with V a complex vector, M a non-Hermitian diagonalizable matrix and T a complex tensor, all of them constant in time. Hereafter, the Einstein summation rule is adopted for Latin indices but not for the Greek ones. We then transform Eq.(6) into a integral equation:

[TABLE]

where the propagator $G_{lm}(t,t^{\prime})$ satisfies Green’s linear equation

[TABLE]

In the Supplementary Material (SM), we prove that the solution,

[TABLE]

satisfies Eq.(8), where $\gamma^{\left(\alpha\right)}$ and $\psi_{m}^{\left(\alpha\right)}$ are, respectively, the eigenvalues and eigenvectors of $M_{lm}$, with $\theta\left(t-t^{\prime}\right)$ the Heaviside step function. We have also defined the vector

[TABLE]

with $\left[\mathbf{J}\right]^{\left(\alpha\beta\right)}\equiv\psi_{l}^{\left(\alpha\right)}\psi_{l}^{\left(\beta\right)}$. $\psi_{l}^{\left(\alpha\right)}$ and $\phi_{m}^{\left(\alpha\right)}$ satisfy the closure relation $\sum_{\alpha}\psi_{l}^{\left(\alpha\right)}\phi_{m}^{\left(\alpha\right)}=\delta_{lm}$, and no orthogonality condition is assumed between them (see SM).

Substituting (9) in the integral equation (7), we obtain the more explicit expression:

[TABLE]

where we have found analytically the zeroth-order solution

[TABLE]

with $a^{\left(\alpha\right)}\equiv\phi_{m}^{\left(\alpha\right)}V_{m}/\lambda^{\left(\alpha\right)}$ the steady-state amplitude and $b^{\left(\alpha\right)}\left(0\right)\equiv a^{\left(\alpha\right)}+\phi_{m}^{\left(\alpha\right)}x_{m}\left(0\right)$ the transient amplitude.

The next step is to iterate solution (11), keeping up to first-order terms in $T_{mno}$. As a result, the variables $x_{n}\left(t^{\prime}\right)$ and $x_{o}\left(t^{\prime}\right)$, in Eq.(11), are replaced by the zeroth-order expression (12). The total solution can approximately be written as $x_{l}(t)\approx x_{l}^{(0)}(t)+x_{l}^{(1)}(t)$, where we define the first-order temporal solution

[TABLE]

The terms proportional to the tensor components $T_{mno}$ have been absorbed into the following coefficients:

[TABLE]

Note that no assumption has been made regarding the nature of the generalized interactions $T_{mno}$.

So far we have performed nothing more than what one could call naive perturbation theory. The main problem of simply applying the iterative procedure is the appearance of unphysical secular terms in the perturbative solution. This becomes clear when the time integrals in Eq.(13) are analytically solved, after separating the cases $\gamma^{\left(\alpha\right)}=\gamma^{\left(\beta\right)}$ and $\gamma^{\left(\alpha\right)}=\gamma^{\left(\beta\right)}+\gamma^{\left(\delta\right)}$ from $\gamma^{\left(\alpha\right)}\neq\gamma^{\left(\beta\right)}$ and $\gamma^{\left(\alpha\right)}\neq\gamma^{\left(\beta\right)}+\gamma^{\left(\delta\right)}$. The result at first order reads:

[TABLE]

where a linear growth in time appears in the first line of Eq.(15), multiplying the secular coefficient

[TABLE]

The temporal behavior $te^{-\gamma^{\left(\alpha\right)}t}$ is problematic for Markovian dissipative systems since it represents an energy gain for short times. Such type of growth might even turn first-order terms larger than zeroth-order ones.

In order to eliminate secular terms from the dynamics, we resort to the Poincaré-Lindstedt method  Bender and A. Orszag (1999). It consists in rewriting the eigenvalues in the form $\gamma^{\left(\alpha\right)}=\bar{\gamma}^{\left(\alpha\right)}+\delta\gamma^{\left(\alpha\right)}$, where we define the renormalized eigenspectrum $\bar{\gamma}^{\left(\alpha\right)}\equiv\gamma^{\left(\alpha\right)}-\delta\gamma^{\left(\alpha\right)}$. The frequency $\delta\gamma^{\left(\alpha\right)}$ is a mathematical artifice to be calculated later. We now replace all exponentials in Eq.(15) by the series expansion

[TABLE]

taking care of neglecting second-order terms and higher in $\delta\gamma^{\left(\alpha\right)}$. It is straightforward to realize that the choice $\delta\gamma^{\left(\alpha\right)}=S^{\left(\alpha\right)}/b^{\left(\alpha\right)}$ cancels out the linear growth, and the following expression is finally obtained:

[TABLE]

In Eq.(18), each rectified exponential $e^{-\bar{\gamma}^{\left(\alpha\right)}t}$ is equivalent to the sum over the contribution of an infinite number of terms, coming from all perturbative orders, representing thus an immense improvement over naive perturbative solutions. It is worth mentioning that the new analytical eigenspectrum $\bar{\gamma}^{\left(\alpha\right)}$ does depend on the initial conditions through amplitudes $b^{\left(\alpha\right)}\left(0\right)$, a characteristic absent in linear algebra approaches.

In order to clearly understand the implications of our general formalism, let us consider the semi-classical Eqs.(4) for cooperative spontaneous emission as a particular case of Eq.(6). Replacing the general components $x_{l}$ by the expectation values \bigl{\langle}\sigma_{l}^{\alpha}\bigr{\rangle}, for $\alpha=\pm,z$, M becomes a $3N\times 3N$ diagonal matrix, where the diagonal block of three single-atom eigenvalues, $\left(\Gamma/2-\Delta i,\Gamma/2+\Delta i,\Gamma\right)$, is the same for all atoms. Hence the associated eigenvectors correspond to the well-known canonical basis, which reduces the general expression (18) to

[TABLE]

where we have introduced the semi-classical eigenspectrum at first-order

[TABLE]

and the structure factor

[TABLE]

Note that the equations above incorporate the initial conditions  (5). Moreover, as a consequence of the renormalization protocol, $\eta_{l}$ and the excited population $p_{e}$ bring cooperative corrections to the independent atoms eigenvalues.

Let us now turn our attention to the linear portion of the analytical eigenspectrum, in Eq.(20a). Note that the dimensionless decay rates $1+\Re\left(\eta_{l}\right)$ and energy shifts $\Im\left(\eta_{l}\right)$ contain information on the initial state (5), where a laser induces phases in the dipoles (see Eq.(21)). Thus, even in the classical regime ($p_{e}\ll 1$), the effective spectrum (20a) differs substantially from the one extracted by diagonalizing the linear optics equations. We compare both in Fig.1 for a Gaussian atomic cloud of standard deviation $R$, with a degree of excitation $p_{e}=0.01$. For the exact linear-optics eigenvalues, the average of the decay rates is the single-atom linewidth for any initial condition. When the same average is performed in Eq.(20a), the result is $1+b_{0}/8$, with $b_{0}=2N/k_{0}^{2}R^{2}$ the cloud resonant optical thickness (see SM), i.e., the superradiant modes are largely populated. In the SM, we show that uncorrelated phases from the atomic positions bring the average decay rate back to the single atom one, which confirms that the superradiance stems from the phase coherence Scully et al. (2006).

Before we proceed to studying the radiated power, let us remind that the linear scaling with $b_{0}$ for super- Bienaimé et al. (2011); Araújo et al. (2016) and subradiance Guerin et al. (2016, 2017) has been observed in the linear optics regime, and this parameter has been identified as the control parameter for cooperative phenomena in dilute clouds. For the nonlinear regime, the first-order rate $\bar{\gamma}_{z}^{\left(l\right)}$ indicates that $b_{0}$ has to be corrected by a factor $1-p_{e}$, so we define $\mathcal{C}\equiv\left(1-p_{e}\right)b_{0}/8$ as the semi-classical cooperative parameter. Accordingly, solution (19) must be valid in the regime $\mathcal{C}\ll 1$.

The first-order solution (19) allows us to obtain a second-order approximation for the total power radiated (Eq.(3)). In the regime of vanishing interactions ($u_{lm}\rightarrow 0$, $\forall m\neq l$), the decay of independent atoms ($Np_{e}e^{-\Gamma t}$) is recovered, as expected. Interestingly, the same free decay is equally obtained when all dipoles tend to become fully-inverted $p_{e}\rightarrow 1$, an evidence that cooperation might be suppressed by atomic population ($\mathcal{C}\rightarrow 0$). Yet the complete absence of cooperation for $p_{e}=1$, which is observed both in numerical simulations and from the renormalization approach (see SM), is an incorrect prediction of the semi-classical truncation, as the emission of the first photons must be treated fully quantum-mechanically Arecchi et al. (1972); MacGillivray and Feld (1976).

In Fig.4(a), exact classical (linear optics) and semi-classical total powers are compared with the one predicted by the analytical solution (19), for half excited population and a cooperativity parameter $\mathcal{C}\approx 0.04$. The linear-optics decay appears farther from the single-atom case than that calculated from the semi-classical model and fails at giving an accurate prediction for the semi-classical emission. However, the renormalization approach provides a very good approximation for the exact semi-classical emission over several $\Gamma t$, including the whole superradiant stage and a long period of the subradiant stage. Yet the superradiance effect is to small to be observed.

In order to clearly see superradiance, in Fig.4(b) we have increased the cooperation to $\mathcal{C}\approx 0.16$. Note that the analytical solution is still capable to describe the semi-classical emission, which is expressively less accelerated than the classical case. However, the agreement with the numerical prediction ends up lasting for a smaller temporal range. Therefore, the parameter $\mathcal{C}$ also controls the temporal validity of solution (19): an accurate description of the dynamics is achieved for stronger cooperation, at the expense of a reduction of the trustful temporal range. It happens because, as the cooperation increases, the most subradiant rates becomes negative, predicting an unphysical energy gain in the emitted power.

To conclude, in this work we have proposed a renormalization formalism to solve truncated nonlinear equations of quantum many-body systems. The feasibility of our method was demonstrated by exploring cooperative spontaneous emission in the semi-classical regime, where a quadratic nonlinearity was addressed. However, the procedure can be straightforwardly extended to other higher-order cumulant expansions. The renormalization technique presents the advantage of solving these nonlinear problems for larger physical systems, with real possibilities of describing properly the thermodynamic limit (see SM). Different from numerical methods, we obtain analytical expressions for many-body dissipative rates and energy shifts, allowing for easy computational parallelization and more accurate physical interpretations.

As shown here at the lowest-order in the cumulant expansion, the interaction for renormalized solutions can be strong enough to distinguish higher-order quantum effects in connected correlations from lower-order ones. Hence collective effects of a truly quantum nature can be associated with a given order of the cumulant expansion for larger systems. The approach is particularly promising to investigate non-equilibrium many-body physics in two and three dimensions, for short and long-range interacting systems.

I Supplementary Material

I.1 Green’s function proof

In this section, we prove that solution (9) satisfies the Green’s differential equation (8). After substituting expression (9) in Eq.(8), we can use the eigenvalue equation $M_{ln}\psi_{n}^{\left(\alpha\right)}=\gamma^{\left(\alpha\right)}\psi_{l}^{\left(\alpha\right)}$ to obtain the condition

[TABLE]

In the step above, it was assumed that the derivative of the Heaveside function results in the Dirac delta function. Note that condition Eq.(S-1) is always true for $t\neq t^{\prime}$. For $t=t^{\prime}$, we obtain the following generalized closure relation:

[TABLE]

Multiplying both sides of Eq.(S-2) by $\psi_{m}^{\left(\gamma\right)}$, we end up with

[TABLE]

Using the definition (10) in (S-3), we finally prove our goal:

[TABLE]

I.2 Alternative way to derive the parameter $\mathcal{C}$ in the large $N$ limit

It is possible to derive $\mathcal{C}$ from the total power radiated Eq.(3) at $t=0$, by evaluating the magnitude of the interacting term of over the corresponding single atom term:

[TABLE]

Therefore, the light-mediated cooperation become substantial when $\mathcal{C}\gtrsim 1$. Considering the initial conditions (5), we obtain $\mathcal{C}=\left(1-p_{e}\right)\left|\Re\left(\left\langle\eta\right\rangle\right)\right|$, with $\left\langle\eta\right\rangle=\sum_{m=1}^{N}\eta_{m}/N$ the average value of the microscopic factor $\eta_{m}$ and $1-p_{e}$ the ground state population.

For macroscopic atomic clouds, with large number of atoms ($N\gg 1$), we can calculate $\left\langle\eta\right\rangle$ substituting sums by integrals, namely,

[TABLE]

with $\rho\left(\mathbf{r}\right)$ the density distribution of the cloud. The change of variables, $\mathbf{r}=\mathbf{r}_{1}-\mathbf{r}_{2}$ and $\mathbf{s}=\left(\mathbf{r}_{1}+\mathbf{r}_{2}\right)/2$, yields

[TABLE]

For a Gaussian distribution $\rho\left(\mathbf{r}\right)=\rho_{0}\exp{\left(-r^{2}/2R^{2}\right)}$, with standard deviation $R$ and peak $\rho_{0}\equiv N/R^{3}\left(2\pi\right)^{3/2}$, note that the following property holds:

[TABLE]

so the radial integrals completely decouple from each other.

Considering that the incident plane wave propagates over the $z$ direction, we can express the phase of the exponential as $\mathbf{k}_{0}\cdot\mathbf{r}=k_{0}r\cos\theta$ to trivially solve the angular integral. In the limit of large atomic clouds ($k_{0}R\gg 2\pi$), the remaining radial integrals

[TABLE]

results in $\left\langle\eta\right\rangle\approx b_{0}/8$ and $\mathcal{C}\approx(1-p_{e})b_{0}/8$, with $b_{0}=2N/k_{0}^{2}R^{2}$ the cloud resonant optical thickness.

I.3 Further discussions

In the following, we demonstrate numerically that the perturbation parameter $\mathcal{C}$ is consistent with the cooperation reduction also for the exact models. We highlight that, we have detected no relevant differences between random realizations of the total power radiated for $N\geqslant 20000$. Also, in order to avoid significant contributions of atomic pairs, an exclusion volume of $0.3\lambda$ was assumed for all atoms in the entire work.

In Fig.3(a), we present the normalized total power radiated, extracted from exact classical (linear optics) and semi-classical (nonlinear) models, where a superradiant stage is followed by a subradiant one. When the excited population is very low ($p_{e}\approx 10^{-3}$), note that $\mathcal{C}\approx b_{0}/8$ and the nonlinear model presents a behavior in excellent agreement with the linear optics model. However, as the excited population increases, a cooperativity reduction of $\mathcal{C}<b_{0}/8$ is expected for the nonlinear model as compared to the linear case. Indeed, the classical curve does not exhibit any change for $p_{e}=0.5$ (not shown since all curves overlap exactly with the one at low $p_{e}$), whereas the semi-classical radiation become closer to the single-atom curve. In this regime, the excited population plays an important role on the coherence dynamics, which is not accounted by the linear optics model.

In the extreme case of a fully inverted system ($p_{e}=1$), a complete absence of cooperation is expected ($\mathcal{C}=0$), as a consequence of the coherence cancelling in the Bloch state (5). In agreement with this prediction, we verify in Fig.3(a) that the numerical curve for the nonlinear model coincides with the single atom emission. In the superradiant cascade, the necessity of treating the emission of the first photons totally quantum-mechanically, until the vacuum modes are significantly populated and the semi-classical modelling can take over, has been pointed out in Ref.Arecchi et al. (1972); MacGillivray and Feld (1976). The present analysis thus confirms that the semi-classical approach does not describe properly that limit. Our analytical solution is naturally loyal to such mathematical prediction, in agreement with $\mathcal{C}=0$. These results show that the renormalization approach predicts qualitatively the many-particle effects of the exact semi-classical model for stronger regimes.

In order to demonstrate the efficiency of the analytical method in reaching the thermodynamic limit, we present in Fig.3(b) a single realization of the averaged dipole dynamics for $N=10^{5}$ atoms. We have used a common desktop computer, equipped with a 6th Generation Intel® Core™ i7 processor, which is capable to run 8 independent tasks in parallel. By calculating 8 temporal points of the dynamics in parallel, from a total of 72 (see Fig.3(b)), we have got a running time of about 52 hours. For direct numerical integration of Eq.(4), typically done with Runge Kutta algorithms, such simulations would last at least 7 times more on the same device. Obviously, such discrepancy increases for workstations and clusters. Additionally to time steps parallelization, our analytical solution also allows for calculating the dynamics of each particle separately, even for long-range interacting systems.

For a different initial state, the eigenspectrum selected by the renormalization approach can be qualitatively different. In Fig.4(a), for example, the case of atoms prepared with a random phase is presented. The absence of correlations between the positions and phases of the atomic dipoles makes the radiation much closer to the single-atom one. As a consequence, the eigenspectrum is centered around the single-atom decay rate. Accordingly, the decay of the total power radiated shows very limited superradiance, which is even weaker in the case of a large excited population ($p_{e}=0.5$). It happens because cooperative effects in large atomic clouds are a consequence of phase profiles strongly correlated with atomic positions Scully (2007).

Another interesting result to point out is the fact that the renormalized decay rates (20) do not depend explicitely on the detuning $\Delta$. This is known feature in linear optics, where the eigenvalues are derived by diagonalization of the coupling matrix underlying the linear dynamics; indeed, the detuning appears there as a diagonal term, so it does not affect directly the eigenspectrum and the associated decay rates. In the renormalization approach, the detuning would appears during the preparation of the system, i.e., in the initial conditions. In this context, a pump largely detuned from the atomic resonance ($(\Delta/\Gamma)^{2}\gg b_{0}$) allows one to avoid multiple scattering and reach a state close to (5). Then, despite being in a single-scattering regime, cooperative effects can be observed Araújo et al. (2016); Roof et al. (2016); Guerin et al. (2016). Our work generalizes this result to the semi-classical dynamics.