Solving Dicke superradiance analytically: A compendium of methods

1. Problem Statement

The paper addresses the Dicke superradiance problem, which describes the collective spontaneous emission of an ensemble of $N$ identical two-level systems (TLS) initially in a fully inverted state.

• Physical Context: Unlike independent emitters where radiation scales linearly with $N$ ($N\Gamma$), a collectively decaying system exhibits a burst of radiation scaling as $N^2$. This occurs due to particle-particle correlations, peaking when the system reaches a half-inverted state.
• Mathematical Challenge: The goal is to find the exact time evolution of the density operator $\rho(t)$ governed by the Lindblad master equation:
  $$\dot{\rho}(t) = \Gamma \left( S\rho S^\dagger - \frac{1}{2} \{S^\dagger S, \rho\} \right)$$
  where $S = \sum_{j=1}^N \sigma_j$ is the collective collapse operator.
• Symmetry: Due to the permutation symmetry of $S$, the dynamics are restricted to the symmetric subspace spanned by $N+1$ Dicke states $|m\rangle$ (where $m$ is the number of excitations, $0 \le m \le N$). The density matrix remains diagonal in this basis, reducing the problem to finding the populations $\rho_m(t)$.
• Historical Context: While exact solutions exist (e.g., Lee, 1977), they are often considered computationally cumbersome or "convoluted," leading recent research to focus on large-$N$ approximations. This paper aims to provide a unified, accessible, and fully analytical framework valid for any $N$ and any time $t$.

2. Methodology

The authors present five distinct analytical approaches to solve the master equation, all converging on the same result. They utilize the formal time-series expansion $\rho(t) = \sum \frac{t^n}{n!} \mathcal{L}^n[\rho_0]$ and various mathematical techniques to sum this series.

A. Quantum Jump Unraveling (Primary Approach)

This is highlighted as the most elegant method.

• Concept: The master equation is unraveled into quantum trajectories. Each trajectory involves $N$ discrete jumps from state $|N\rangle$ to $|0\rangle$ at random times determined by uniform random variables.
• Derivation: The population $\rho_m(t)$ is calculated as the probability that a system remains in state $|m\rangle$ at time $t$ averaged over all possible trajectories.
• Result: By transforming the integration variables from random probabilities to time intervals and utilizing the integral representation of the Heaviside step function, the solution is derived as a closed contour integral in the complex plane.

B. Combinatorics Approach

• Concept: The superoperator $\mathcal{L}$ is split into a decay part ($B$) and a non-decay part ($C$). The time evolution is expanded as a binomial series of non-commuting operators $(B+C)^n$.
• Derivation: The authors sum over all possible "paths" (sequences of $B$ and $C$ operators) that lead from the initial state to a target state $m$.
• Result: The constrained sums over operator sequences are identified as coefficients in a generating function, which are then evaluated using residue calculus.

C. Probabilistic Approach

• Concept: The time evolution is discretized into small steps $\Delta t$. At each step, the system either jumps (probability $d_k$) or stays (probability $s_k$).
• Derivation: The population is expressed as a sum over all possible sequences of jumps and stays. This is mapped to a generating function problem similar to the combinatorics approach.
• Result: In the limit $\Delta t \to 0$, the discrete sum converges to the same contour integral solution.

D. Rate Equations (Recursive) Approach

• Concept: The master equation is converted into a set of coupled recursive rate equations for the coefficients of the time series expansion.
• Derivation: A general recursive formula is derived connecting the $n$-th order term of state $m-1$ to the $n$-th and $(n+1)$-th order terms of state $m$.
• Result: By solving the recursion and summing the series, the authors recover the residue sum solution, noting that the form naturally suggests a generalization to all states via the residue theorem.

E. Non-Hermitian Hamiltonian Approach

• Concept: The rate equations are cast into a matrix form $\dot{\rho} = \Gamma H \rho$, where $H$ is a non-Hermitian matrix.
• Key Feature: Due to the symmetry $h_m = h_{N-m}$, the matrix $H$ possesses Exceptional Points (EPs) (specifically EP2s) where eigenvalues and eigenvectors coalesce.
• Derivation: The solution requires a Jordan decomposition ($H = T J T^{-1}$). The presence of Jordan blocks of size 2 leads to decay terms of the form $t e^{-\lambda t}$ rather than pure exponentials.
• Result: The solution is obtained via matrix exponentiation or Laplace transforms (resolvent method), yielding the same contour integral form.

3. Key Contributions and Results

The Unified Analytical Solution

The paper culminates in a compact, exact solution for the population of any Dicke state $|m\rangle$ at time $t$:
$$\rho_m(t) = \frac{1}{2\pi i} \oint_C dz \, f_m(z, t) = \sum_{p \in \text{poles}} \text{Res}\left[ f_m(z, t) \right]_{z=p}$$
where the integrand is:
$$f_m(z, t) = \frac{(-1)^{N-m} h_N \cdots h_{m+1}}{(z-h_N)\cdots(z-h_m)} e^{-z\Gamma t}$$
Here, $h_m = m(N+1-m)$ are the coupling coefficients.

Key features of the solution:

1. Residue Sum: The time evolution is expressed entirely as a sum of residues of a meromorphic function in the complex plane.
2. Pole Structure:
   • For states above the "equator" ($m \ge (N+1)/2$), the poles are simple.
   • For states below the equator, degenerate pairs of poles ($h_m = h_{N-m}$) create double poles, requiring the derivative of the function for the residue calculation. This mathematically encodes the non-exponential decay behavior associated with Exceptional Points.
3. Generality: The solution is generalized to arbitrary initial Dicke states $|m_0\rangle$, not just the fully inverted state.

Computational Efficiency

The authors demonstrate that while the solution involves a sum over poles, it is computationally efficient compared to direct numerical integration of the master equation, especially for large $N$, as it avoids the $2^N$ Hilbert space dimension by restricting to the symmetric subspace and using analytical summation.

4. Significance and Outlook

• Revival of Exact Solutions: The paper revitalizes the exact analytical solutions proposed by Lee (1977), which had been largely overlooked due to their perceived complexity. The authors reframe these solutions using modern tools (contour integration, quantum trajectories) to make them accessible.
• Beyond Dicke Superradiance: The authors argue that the existence of such analytical solutions is not unique to this problem. The methods (particularly the quantum jump unraveling leading to contour integrals) could be applied to more complex open quantum systems, such as:
  • Coherently driven Dicke systems.
  • Incoherently pumped systems (relevant for superradiant lasers).
• Physical Insight: The derivation explicitly links the mathematical structure of the solution (double poles/Jordan blocks) to the physical phenomenon of Exceptional Points in non-Hermitian physics, providing a clear theoretical foundation for the non-exponential decay dynamics observed in superradiance.
• Experimental Relevance: The paper notes that while the original Dicke model is hard to realize, modern platforms like cavity QED, waveguide QED, and superconducting circuits allow for the realization of these dynamics, making the analytical solution directly applicable to current experimental data analysis.

In summary, this paper provides a comprehensive "compendium" that unifies various mathematical techniques to solve a fundamental problem in open quantum systems, offering a powerful, exact, and elegant analytical tool for studying collective quantum phenomena.