# Mesoscopic Cooperative Emission from a Disordered System

## I Introduction

The study of cooperative phenomena in optics was initiated by the pioneering work of Dicke. The underlying physics of the cooperative emission can be readily understood using a classical approach. Suppose, that a large number $`N`$ of identical oscillators with frequency $`\omega _0`$ are confined within a small volume with characteristic size $`L2\pi c/\omega _0=\lambda _0`$, where $`\lambda _0`$ is the radiation wavelengh; this is reffered to as a “point” sample. If $`\tau `$ is the radiative lifetime of an isolated oscillator, then according to Dicke, the $`N`$ eigenmodes of the system of oscillators consist of one mode with a short lifetime $`\tau /N\tau `$, and $`N1`$ modes with lifetimes much longer than $`\tau `$ \[by a factor $`(L/\lambda _0)^2`$\]. Correspondingly, the emission spectrum of this system consists of superimposed broad (superradiant) and a narrow (subradiant) bands. The intensities ratio of these bands is determined by the details of the excitation. This type of lifetimes redistribution is caused by the interactions among the oscillators through their radiation fields.

Certainly, the classical picture does not describe all aspects of the cooperative emission. In fact, the original work of Dicke primarily addressed the time evolution of the radiation emission, provided that at the initial moment, $`t=0`$, all the oscillators are coherently excited. For this situation, the classical picture helps in understanding that the radiation is released during a short time, $`\tau /N`$; understanding of the initial stages of the emission process (the delay time statistics) requires, however, a quantum description. The original treatment in Ref.  also ignored the dipole–dipole interactions, which give rise to a spread in the oscillators frequencies (dephasing). The question whether or not this dephasing would completely destroy the cooperative emission is very non–trivial and was addressed in a number of later works.

In the previous considerations of cooperative emission, it was assumed that all $`N`$ oscillators (atoms, molecules or excitons) have the same frequencies. Such a restriction was adequate for the experimental situation in both gases and single crystals. To the best of our knowledge, the only account of disorder in the frequencies of the oscillators was given in Ref. , which addressed the transient behavior of the cooperative emission from molecular aggregates. The case of $`J`$–aggregates corresponds to a symmetrical arrangement of oscillators in a circle. The authors treated the disorder within the perturbation theory and averaged the second–order correction to the decay rates (the first–order correction vanishes upon averaging) with a Gaussian distribution. The advantage of the work in Ref. is that the nearest neighbors dipole–dipole interactions were taken into account exactly. The drawback is in the perturbative approach, which rules out certain $`\mathrm{q}\mathrm{u}\mathrm{a}\mathrm{l}\mathrm{i}\mathrm{t}\mathrm{a}\mathrm{t}\mathrm{i}\mathrm{v}\mathrm{e}`$ physical effects (see below).

Whereas Ref. addressed a rather particular situation, the following general questions might be asked. Suppose that the oscillators frequencies are randomly distributed with a characteristic width $`\mathrm{\Omega }`$. Obviously, as $`\mathrm{\Omega }`$ increases, it would eventually destroy the cooperative features in the emission spectrum. Then what is the critical magnitude of $`\mathrm{\Omega }`$? How does this magnitude depend on the parameters of the system $`N`$, $`L`$, and $`\lambda _0`$? What is the structure of the emission spectrum when disorder is smaller than critical?

These questions have become not purely academic due to the recent advances in the field of laser–action in $`\pi `$–conjugated polymers. Some experiments provide a strong evidence for cooperative emission from an ensemble of excitons in these materials, for excitation intensities exceeding a certain characteristic threshold. On the other hand, it is well known that the films of $`\pi `$-conjugated polymers are strongly disordered (in the absense of disorder, cooperative emission by polymer chain was considered in Ref. ). They contain impurities and defects which break the polymer chains into segments of relatively short conjugation length, with a distribution depending on the film quality. This has a direct effect on the exciton energy, $`\mathrm{}\omega `$, since it has been found that $`\mathrm{}\omega `$ directly depends on the chain conjugation length.

The questions formulated above are addressed in the present paper. We study here the effect of disorder on the cooperative emission spectrum of the system of classical oscillators. We consider the situation of incoherent excitation, which is most relevant to the experiment. In contrast to Ref. , we are interested in the nonaveraged (but universal) properties of the emission spectrum. In other words, our goal is to assess the $`\mathrm{m}\mathrm{e}\mathrm{s}\mathrm{o}\mathrm{s}\mathrm{c}\mathrm{o}\mathrm{p}\mathrm{i}\mathrm{c}`$ aspects of the cooperative emission. By mesoscopic we mean that, in the presence of a disorder, the emission spectrum of a large number of oscillators develops a fine structure. The actual shape of this spectral structure represents the fingerprints of the distribution of the oscillator frequencies and positions for a given realization. At the same time, the characteristic period and amplitude of the fine structure are determined by the net parameters of the system: $`N`$, $`L`$, and $`\mathrm{\Omega }`$.

The paper is organized as follows. In Section II we derive the expression for the emission spectrum of a system of classical oscillators coupled by their radiation fields. In Section III we study in detail a simplified model in which the coupling among the oscillators is independent of distance. The eigenmodes of a “point” sample in the presence of disorder are analyzed in Section IV. The universal properties of the mesoscopic stucture in the emission spectrum for small and large (but still smaller than $`\lambda _0`$) sizes $`L`$ are discussed in Sections V and VI, respectively. In Section VII the effect of the dipole–dipole interactions is addressed. The conclusions are given in Section VIII.

## II The Basic Equations

We consider a system of $`N`$ oscillators located at random points $`𝐫_i`$, with frequencies $`\omega _i`$ randomly distributed around a central frequency $`\omega _0`$ with a characteristic width $`\mathrm{\Omega }`$. Each oscillator is driven by the radiation field $`𝐄(𝐫,t)`$ produced by all oscillators. The equation of motion for the displacement $`u_i`$ of a given oscillator $`i`$ reads

$$\ddot{u}_i+\omega _i^2u_i=\frac{e}{m}𝐧_i𝐄(𝐫_i,t),$$
(1)   
where $`e`$ and $`m`$ are the dipole characteristics (effective charge and mass) and $`𝐧_i`$ is a unit vector in the direction of the dipole moment.

The current density, associated with the oscillators motion, can be written as

$$𝐉(𝐫,t)=e\underset{i}{}𝐧_i\dot{u}_i\delta (𝐫𝐫_i).$$
(2)   
The current $`𝐉`$ plays the role of a source, which generates the electric field $`𝐄(𝐫,t)`$ according to

$$\mathrm{\Delta }𝐄\frac{1}{c^2}\ddot{𝐄}=\frac{4\pi }{c^2}\dot{𝐉},$$
(3)   
where $`c`$ is the speed of light.

Within the classical approach, the emission spectrum of the system should be calculated as follows. We assume that at the initial moment, $`t=0`$, all oscillators are excited with different phases $`\varphi _i`$, and that the radiation field at the initial moment is zero, $`𝐄(𝐫,0)=0`$. The evolution of $`𝐄`$ with time can be then obtained by solving Eqs. (1)–(3). After taking the limit $`r\mathrm{}`$ and expanding the field into harmonics, the spectral intensity can be obtained as $`I(\omega )=|𝐄(\mathrm{},\omega )|^2`$.

To carry out this program, it is convenient to employ the Laplace transformation. The transformed functions $`\overline{u}_i(p)`$ and $`\overline{𝐄}(𝐫,p)`$ satisfy the following system of equations

$`(\omega _i^2+p^2)\overline{u}_i(p)`$ $`=`$ $`{\displaystyle \frac{e}{m}}𝐧_i\overline{𝐄}(𝐫_i,p)+u_0(p\mathrm{cos}\varphi _i\omega _i\mathrm{sin}\varphi _i),`$ (4)   
$`\mathrm{\Delta }\overline{𝐄}(𝐫,p){\displaystyle \frac{p^2}{c^2}}\overline{𝐄}(𝐫,p)`$ $`=`$ $`{\displaystyle \frac{4\pi e}{c^2}}{\displaystyle \underset{i}{}}𝐧_i\left[p^2\overline{u}_i(p)u_0(p\mathrm{cos}\varphi _i\omega _i\mathrm{sin}\varphi _i)\right]\delta (𝐫𝐫_i),`$ (5)   
where $`u_0\mathrm{sin}\varphi _i`$ and $`\omega _iu_0\mathrm{cos}\varphi _i`$ are the respective initial displacement and velocity of the $`i`$th oscillator. The solution of Eq. (5) for $`\overline{𝐄}(𝐫,p)`$ can be presented as a superposition of eigenmodes, $`𝐄_\nu (𝐫)`$, of the wave equation for the electromagnetic field,

$$\mathrm{\Delta }𝐄_\nu (𝐫)+\frac{\omega _\nu ^2}{c^2}𝐄_\nu (𝐫)=0,$$
(6)   
where $`\omega _\nu `$ is the eigenfrequency. Assuming that the modes are normalized, $`𝑑\mathrm{𝐫𝐄}_\nu ^2(𝐫)=1`$, we obtain the following expression for $`\overline{𝐄}(𝐫,p)`$

$$\overline{𝐄}(𝐫,p)=4\pi e\underset{i\nu }{}\left[p^2\overline{u}_i(p)u_0(p\mathrm{cos}\varphi _i\omega _i\mathrm{sin}\varphi _i)\right]\frac{𝐧_i𝐄_\nu (𝐫_i)}{\omega _\nu ^2+p^2}𝐄_\nu (𝐫).$$
(7)   
Substituting Eq. (7) into Eq. (4), we get a system of coupled equations for the amplitudes $`\overline{u}_i(p)`$

$`(\omega _i^2+p^2)\overline{u}_i(p)={\displaystyle \frac{4\pi e^2}{m}}{\displaystyle \underset{j\nu }{}}{\displaystyle \frac{\left[𝐧_i𝐄_\nu (𝐫_i)\right]\left[𝐧_j𝐄_\nu (𝐫_j)\right]}{\omega _\nu ^2+p^2}}`$ $`\left[p^2\overline{u}_j(p)u_0(p\mathrm{cos}\varphi _j\omega _j\mathrm{sin}\varphi _j)\right]`$ (9)   
$`+u_0(p\mathrm{cos}\varphi _i\omega _i\mathrm{sin}\varphi _i).`$   
To simplify Eq. (9), it is convenient to introduce new variables $`v_i(p)`$:

$$v_i(p)=\frac{p^2\overline{u}_i(p)}{u_0}(p\mathrm{cos}\varphi _i\omega _i\mathrm{sin}\varphi _i).$$
(10)   
Then Eq. (9) takes the form

$$(\omega _i^2+p^2)v_i+\underset{j}{}S_{ij}v_j=\omega _i^2(\omega _i\mathrm{sin}\varphi _ip\mathrm{cos}\varphi _i),$$
(11)   
where the coefficients

$$S_{ij}(p)=\frac{4\pi e^2p^2}{m}\underset{\nu }{}\frac{\left[𝐧_i𝐄_\nu (𝐫_i)\right]\left[𝐧_j𝐄_\nu (𝐫_j)\right]}{\omega _\nu ^2+p^2}$$
(12)   
describe the coupling between oscillators $`i`$ and $`j`$ via the radiation field.

Let us now express the intensity, $`I(\omega )`$, in terms of the variables $`v_i(p)`$. The expression for $`\overline{𝐄}(𝐫,p)`$ follows from Eqs. (7) and $`(\text{10})`$

$$\overline{𝐄}(𝐫,p)=4\pi e\underset{i\nu }{}v_i(p)\frac{𝐧_i𝐄_\nu (𝐫_i)}{\omega _\nu ^2+p^2}𝐄_\nu (𝐫).$$
(13)   
The Fourier transform of the electric field is obtained by replacing $`p`$ by $`i\omega `$ in Eq. (13). In the limit $`r\mathrm{}`$, only the pole $`\omega _\nu =\omega `$ contributes to the sum over $`\nu `$, so that

$$𝐄(𝐫,\omega )|_r\mathrm{}\underset{i\nu }{}v_i\left[𝐧_i𝐄_\nu (𝐫_i)\right]𝐄_\nu (𝐫)\delta (\omega _\nu ^2\omega ^2).$$
(14)   
This corresponds to taking the continuum limit for electromagnetic modes. The terms proportional to $`\left[𝐄_\nu (𝐫)𝐄_\mu (𝐫)\right]`$, which appear after calculating $`|𝐄(𝐫,\omega )|^2`$ from Eq. (14), oscillate rapidly if $`\mu \nu `$. Therefore, only the terms with $`\mu =\nu `$ survive at large $`r`$. These terms contain products of the form $`\left[𝐧_i𝐄_\nu (𝐫_i)\right]\left[𝐧_j𝐄_\nu (𝐫_j)\right]`$. Note that the same products enter into the coupling coefficients, $`S_{ij}`$, defined by Eq. (12). This allows us to present the final expression for the spectral intensity in a compact form

$$I(\omega )\underset{ij}{}v_i(\mathrm{Im}S_{ij})v_j^{},$$
(15)   
where $`v_i`$ and $`S_{ij}`$ are calculated at $`p=i\omega `$.

We assume that the spread of the oscillators frequencies due to the disorder is much smaller than the central frequency, $`\mathrm{\Omega }\omega _0`$. This means that the frequency dependence of the coupling constants is weak, so that $`S_{ij}(i\omega )`$ can be evaluated at $`\omega =\omega _0`$. The real part of $`S_{ij}`$, which comes from the principle value of the sum over modes in Eq. (12), diverges for $`i=j`$. This divergency is the manifestation of the Lamb shift, well–known in quantum electrodynamics, and can be absorbed into $`\omega _i`$. At the same time, the imaginary part of $`S_{ii}`$, which results from the pole $`\omega _\nu =\omega _0`$, is finite. It determines the radiative lifetime, $`\tau `$, of an individual oscillator via the relation $`\mathrm{Im}S_{ii}(i\omega _0)=2\omega _0/\tau `$. For a single oscillator in vacuum, the modes $`𝐄_\nu `$ are simply plane waves, and the summation over $`\nu `$ in Eq. (12) recovers the textbook result

$$\tau =\frac{3mc^3}{e^2\omega _0^2}.$$
(16)   
For $`ij`$, the coupling $`S_{ij}`$ between two oscillators depends on the ratio $`r_{ij}/\lambda _0`$, where $`r_{ij}`$ is the distance between the oscillators, and $`\lambda _0=2\pi c/\omega _0`$ is the radiation wavelength. For $`r_{ij}\lambda _0`$, both real and imaginary parts of $`S_{ij}`$ oscillate rapidly with $`r_{ij}`$, and the effect of coupling is negligibly small for a large ensemble of oscillators. For $`r_{ij}\lambda _0`$, the real part of $`S_{ij}`$ represents the dipole–dipole interaction of the oscillators $`i`$ and $`j`$. It is convenient to present $`S_{ij}`$ in the form

$$S_{ij}=\frac{2\omega _0}{\tau }(\beta _{ij}+i\alpha _{ij}),$$
(17)   
where $`\beta _{ij}`$ and $`\alpha _{ij}`$ are the dimensionless matrices of coupling between the oscillators, defined as

$$\beta _{ij}=\left(\frac{\lambda _0}{2\pi r_{ij}}\right)^3\left[(𝐧_i𝐧_j)\frac{3(𝐧_i𝐫_{ij})(𝐧_j𝐫_{ij})}{r_{ij}^2}\right],$$
(18)   
and

$`\alpha _{ij}=𝐧_i𝐧_j{\displaystyle \frac{1}{5}}\left({\displaystyle \frac{2\pi r_{ij}}{\lambda _0}}\right)^2\left[(𝐧_i𝐧_j){\displaystyle \frac{(𝐧_i𝐫_{ij})(𝐧_j𝐫_{ij})}{2r_{ij}^2}}\right].`$ (19)   
Turning back to Eq. (11), we note that since the distribution of oscillators frequencies is relatively narrow, that is $`\mathrm{\Omega }\omega _0`$, we can make some simplifications. Namely, for $`p=i\omega `$, the factor $`(\omega _i^2+p^2)`$ in the lhs can be replaced by $`2\omega _0(\omega _i\omega )`$, and the rhs can be written as $`i\omega _0^3e^{i\varphi _i}`$. Finally, after rescaling $`v_i`$ by factor $`\omega _0^2`$, Eq. (11) takes the form

$$(\omega _i\omega )v_i+\frac{1}{\tau }\underset{j}{}(\beta _{ij}+i\alpha _{ij})v_j=\frac{i}{2}e^{i\varphi _i}.$$
(20)   
Equation (20) together with Eqs. (15) and (17)–(19) allow us to calculate the spectral intensity $`I(\omega )`$ for any set of initial oscillators phases. For $`\alpha _{ij}=\beta _{ij}=0`$ ($`ij`$), the eigenfrequencies of the system are simply the frequencies of individual oscillators, and the the emission spectum represents a superposition of Lorentzian peaks centered at $`\omega _i`$. In the presence of nondiagonal coupling, the eigenfrequencies are those of cooperative eigenmodes which, in turn, are determined by the imaginary part of the coupling, $`\alpha _{ij}`$. In the experiment, the measured spectrum represents the result of averaging over many excitation pulses. In order to simulate the experimental situation, we will assume the phases $`\varphi _i`$ to be uncorrelated random numbers and average the result for the spectral intensity over all $`\varphi _i`$.

## III A Simple Model

In this section we consider a simplified situation, in which Eq. (20) with random frequencies $`\omega _i`$ can be solved exactly and the expression for the spectral intensity can be obtained in a closed form. Following Dicke, we disregard the dipole–dipole interactions by setting $`\beta _{ij}=0`$. Although this approximation is rather common, later on we will discuss it in more detail. Turning to $`\alpha _{ij}`$, we note that since $`L^2/\lambda _0^21`$, the second term in Eq. (19) is a small correction to the first term. We therefore approximate the non–diagonal elements of $`\alpha _{ij}`$ by replacing $`r_{ij}^2/\lambda _0^2`$ with its average,

$$\alpha _{ij}=\alpha 𝐧_i𝐧_j,\alpha _{ii}=1,$$
(21)   
where the coupling constant $`\alpha `$, with a typical value $`(1\alpha )L^2/\lambda _0^21`$, is the same for all pairs. Note however, that the disorder coming from random orientations of $`𝐧_i`$ is still included. Later we will use this model for the analysis of the system (20) with realistic $`\alpha _{ij}`$.

### A General solution

For the model coupling (20), the system of equations (20) takes the form

$$\left[\omega _i\omega +\frac{i}{\tau }(1\alpha )\right]v_i+\frac{i}{\tau }\alpha 𝐧_i𝐬=\frac{i}{2}e^{i\varphi _i},$$
(22)   
with vector $`𝐬`$ defined as

$$𝐬=\underset{i=1}{\overset{N}{}}v_i𝐧_i.$$
(23)   
A closed equation for $`𝐬`$ can be obtained by multiplying $`v_i`$, found from Eq. (22), by $`𝐧_i`$ and taking the sum over $`i`$. This yields

$$𝐬+\frac{i\alpha }{\tau }\underset{i}{}\frac{𝐧_i(𝐧_i𝐬)}{\omega _i\omega +i(1\alpha )/\tau }=\frac{i}{2}\underset{i}{}\frac{𝐧_ie^{i\varphi _i}}{\omega _i\omega +i(1\alpha )/\tau }.$$
(24)   
Solving Eqs. (22) and (24) for $`v_i`$ and substituting the result into Eq. (15), we obtain for the spectral intensity after some algebra

$$I(\omega )\mathrm{Im}\left[f(\omega )\frac{i\alpha }{\tau }\underset{\mu \nu }{}g_\mu ^{}\left(1+\frac{i\alpha }{\tau }F\right)_{\mu \nu }^1g_\nu ^+\right],$$
(25)   
where we introduced a function

$$f(\omega )=\underset{i}{}\frac{1}{\omega _i\omega +i(1\alpha )/\tau },$$
(26)   
a vector

$$g_\mu ^\pm (\omega )=\underset{i}{}\frac{e^{\pm i\varphi _i}n_{i\mu }}{\omega _i\omega +i(1\alpha )/\tau },$$
(27)   
and a tensor

$$F_{\mu \nu }(\omega )=\underset{i}{}\frac{n_{i\mu }n_{i\nu }}{\omega _i\omega +i(1\alpha )/\tau },$$
(28)   
where $`n_{i\mu }`$ are the components of $`𝐧_i`$.

### B Identical oscillators

Let us first consider the case of $`N`$ identical oscillators having the same frequencies $`\omega _i=\omega _0`$, and dipole momenta all aligned in the same direction. Then we find from Eq. (28),

$$F_{\mu \nu }(\omega )=\delta _{\mu \nu }f(\omega )=\frac{\delta _{\mu \nu }N}{\omega _i\omega +i(1\alpha )/\tau },$$
(29)   
and after averaging over the initial phases $`\varphi _i`$, we obtain from Eq. (25)

$`I(\omega )\left[{\displaystyle \frac{(N1)(1\alpha )/\tau }{(\omega _0\omega )^2+(1\alpha )^2/\tau ^2}}+{\displaystyle \frac{(1\alpha +\alpha N)/\tau }{(\omega _0\omega )^2+(1\alpha +\alpha N)^2/\tau ^2}}\right].`$ (30)   
The emission spectrum is a superposition of a wide and a narrow Lorentzians with spectral widths $`\mathrm{\Gamma }N/\tau `$ and $`\gamma =(1\alpha )/\tau `$, respectively. In accordance to the classical result, the eigenmodes of the system of $`N`$ identical oscillators coupled via their radiation field represent a single superradiant mode with short radiation time $`\tau /N`$, and $`N1`$ subradiant modes with radiation time much longer than that for an isolated oscillator, $`\tau /(1\alpha )\tau `$. The superradiant mode is a symmetric superposition of oscillator states and is strongly coupled to the radiation field, whereas the coupling of the subradiant modes to the radiation field is suppressed. In this case, the frequencies of all $`N1`$ subradiant modes are degenerate, and the spectrum consists of a single narrow peak of width $`\gamma `$ on top of much broader band of width $`\mathrm{\Gamma }`$, as shown in Fig. 1. As can be seen, with decreasing $`1\alpha `$, the hight of the subradiant peak increases, whereas the amplitude of the superradiand band diminishes.

### C Random frequencies

Consider now the case when the oscillators frequencies are random, but orientational disorder is still absent, i.e. all dipoles are aligned in one direction. Again we have $`F_{\mu \nu }(\omega )=\delta _{\mu \nu }f(\omega )`$, with $`f(\omega )=f^{}(\omega )+if^{\prime \prime }(\omega )`$ given by Eq. (26). Then a straightforward evaluation of Eq. (25) yields (after averaging over the phases)

$`I(\omega )\left[{\displaystyle \frac{\frac{\alpha }{\tau }f_1^{}\left(1\frac{\alpha }{\tau }f^{\prime \prime }\right)+\left(\frac{\alpha }{\tau }\right)^2f_1^{\prime \prime }f^{}}{\left(1\frac{\alpha }{\tau }f^{\prime \prime }\right)^2+\left(\frac{\alpha }{\tau }f^{}\right)^2}}f^{\prime \prime }\right].`$ (31)   
where the function $`f_1(\omega )=f_1^{}(\omega )+if_1^{\prime \prime }(\omega )`$ is defined as

$$f_1(\omega )=\underset{i}{}\frac{1}{\left[\omega _i\omega +i(1\alpha )/\tau \right]^2}.$$
(32)   
In order to clarify the underlying physics, it is useful to express the spectral intensity in terms of the system eigenmodes. The eigenfrequencies $`\stackrel{~}{\omega }_k`$ are determined by the equation:

$$1+\frac{i\alpha }{\tau }f(\stackrel{~}{\omega }_k)=0.$$
(33)   
Then the intensity Eq. (31) can be simply rewritten as

$$I(\omega )\underset{k}{}\frac{\stackrel{~}{\omega }_k^{\prime \prime }}{(\omega \stackrel{~}{\omega }_k^{})^2+\stackrel{~}{\omega }_k^{\prime \prime 2}},$$
(34)   
where $`\stackrel{~}{\omega }_k^{}=\mathrm{Re}\stackrel{~}{\omega }_k`$ is the eigenmode frequency and $`\stackrel{~}{\omega }_k^{\prime \prime }=\mathrm{Im}\stackrel{~}{\omega }_k`$ characterizes its width. Note that for $`\omega _i=\omega _0`$, we have $`N1`$ degenerate eigenmodes with $`\stackrel{~}{\omega }_k^{}=\omega _0`$, and Eq. (34) turns into Eq. (30).

### D Disorder in orientations

In the presence of the orientational disorder, the spectral intensity (25) depends, in principle, on the direction of each $`𝐧_i`$. However, for large $`N`$, one can replace the product $`n_{i\mu }n_{i\nu }`$ in Eq. (28) for $`F_{\mu \nu }`$, with its average,

$$n_{i\mu }n_{i\nu }=\frac{1}{3}\delta _{\mu \nu }.$$
(35)   
Thus, we have $`F_{\mu \nu }(\omega )=\frac{1}{3}\delta _{\mu \nu }f(\omega )`$, so that the expression for the spectral intensity is similar to Eq. (31) with the only difference that in the first term, the functions $`f(\omega )`$ and $`f_1(\omega )`$ are now multiplied by $`1/3`$. This results in a shrinkage of the superradiant emission band by the same factor. At the same time, the width of subradiant peak increases by a factor of 3. Thus, the orientational disorder has no qualitative effect on the cooperative emission spectrum. The reason is that the coupling (20) is separable, that is it depends on orientations via the product $`𝐧_i𝐧_j`$. Furthermore, for realistic $`\alpha _{ij}`$ given by Eq. (19), the main (first) term has the same separable form; therefore, the orientational disorder does not qualitatively affect the cooperative emission spectrum and will be disregarded in the rest of the paper.

### E Numerical results

In Fig. 2 we plot the normalized spectral intensity in the absence of coupling, i.e. $`\alpha =0`$, with increasing number of oscillators. Each spectrum is calculated with a computer generated set of $`N`$ random frequencies $`\omega _i`$, which we have chosen, for simplicity, to be uniformly distributed in the interval $`(\omega _0\mathrm{\Omega },\omega _0+\mathrm{\Omega })`$. For convenience, the spectra corresponding to different $`N`$ are normalized and shifted in the vertical direction. It can be seen that the peaks are resolved in the spectrum as long as the disorder, $`\mathrm{\Omega }`$, is larger than $`N/\tau `$. We also see that for sufficiently large $`N`$, the intensity peaks are washed out from the spectrum.

In Figs. 3–6 we present the results for $`I(\omega )`$ calculated using Eq. (31) for several values of $`\alpha `$ close to 1. The striking feature of the emission spectrum is its mesoscopic character. In the presence of disorder, the narrow subradiant peak of Eq. (30) (see Fig. 1) is not smeared out due to a large spread in $`\omega _i`$, as in the case of uncoupled oscillators (see Fig. 2), but rather splits into a multitude of narrow peaks corresponding to the eigenmodes of the disordered system. Furthermore, although the curves are calculated with different random sets of frequencies, the overall pattern of the emission spectrum exhibits certain universal features. In particular, it can be seen by comparing Figs. 3–6 that with increasing $`N`$, the random structure of the spectrum undergoes several transformations, and that the characteristic $`N`$, at which the changes in the pattern occur, is sensitive to the proximity of $`\alpha `$ to 1. This indicates a rather non–trivial structure of the eigenmodes, which we address in the next section.

## IV Structure of eigenmodes

The eigenmodes of a system of $`N`$ oscillators coupled through their radiation field are determined by the homogeneous part of Eq. (20) (we set $`\beta _{ij}=0`$ in this section)

$$(\omega _i\omega )v_i+\frac{i}{\tau }\underset{j}{}\alpha _{ij}v_j=0.$$
(36)   
Since the typical values of $`(1\alpha _{ij})r_{ij}^2/\lambda _0^2`$ are small, we split the second term in Eq. (36) into a sum of the main contribution, with $`\alpha _{ij}=1`$, and a correction proportional to $`(\alpha _{ij}1)`$. Analogously to the consideration in the previous section, we rewrite Eq. (36) as

$$(\omega _i\omega )v_i+\frac{i}{\tau }s(1+\sigma _i)=0,$$
(37)   
with

$$s=\underset{j}{}v_j,\sigma _i=\frac{1}{s}\underset{j}{}(\alpha _{ij}1)v_j.$$
(38)   
Expressing $`v_i`$ from Eq. (37) and taking the sum over $`i`$, we obtain

$$1+\frac{i}{\tau }\underset{j}{}\frac{1+\sigma _j}{\omega _j\omega }=0.$$
(39)   
The equation for $`\sigma _i`$ follows from substituting of $`v_j`$, found from Eq. (37), into the definition of $`\sigma _i`$, Eq. (38),

$$\sigma _i+\frac{i}{\tau }\underset{j}{}\frac{\alpha _{ij}1}{\omega _j\omega }(1+\sigma _j)=0.$$
(40)   
The solutions of Eqs. (39) and (40) determine the complex frequencies of the eigenmodes, $`\stackrel{~}{\omega }_k\stackrel{~}{\omega }_k^{}+i\stackrel{~}{\omega }_k^{\prime \prime }`$.

### A “Point” sample

Let us first analyze the effect of disorder on a system with all $`\alpha _{ij}=1`$, corresponding to the limit of a “point” sample, i.e. $`(L/\lambda _0)^21`$. With $`\sigma _i=0`$, the real and imaginary parts of Eq. (39) read

$`{\displaystyle \frac{1}{\tau }}{\displaystyle \underset{j}{}}{\displaystyle \frac{\omega _j\omega ^{}}{(\omega _j\omega ^{})^2+\omega ^{\prime \prime 2}}}=0,`$ (41)   
$`{\displaystyle \frac{1}{\tau }}{\displaystyle \underset{j}{}}{\displaystyle \frac{\omega ^{\prime \prime }}{(\omega _j\omega ^{})^2+\omega ^{\prime \prime 2}}}=1.`$ (42)   
This system of equations has two different solutions with a crossover between them governed by the parameter $`\mathrm{\Omega }\tau /N`$. For large disorder, $`\mathrm{\Omega }N/\tau `$, it can be readily seen that only one term in each of Eqs. (41) and (42) contributes to the sum. In this case, the solutions are simply $`\omega =\omega _j+i/\tau `$, as if the oscillators were uncoupled. In fact, this conclusion could be anticipated. The above parameter represents the ratio of the mean frequency spacing (MFS) of oscillators, $`\mathrm{\Omega }/N`$, and the inverse lifetime of an individual oscillator, $`1/\tau `$; when the former is much larger than the latter, $`\mathrm{\Omega }/N1/\tau `$, the oscillators do not “feel” each other.

In the opposite case of large $`N`$ (or weak disorder), $`N\mathrm{\Omega }\tau `$, the analysis of Eqs. (41) and (42) is carried out as follows. First note that in Eq. (41), which determines the real parts of the eigenfrequencies, $`\stackrel{~}{\omega }_k^{}`$, all the terms in the sum contribute now. Let us drop $`\omega ^{\prime \prime 2}`$ in the denominator of Eq. (41) (this step will be justified below). Then we obtain that the solutions $`\stackrel{~}{\omega }_k^{}`$ are given by the extrema of the polynomial $`P(\omega )=_j(\omega _j\omega )`$. These determine the frequencies of the $`N1`$ subradiant modes. At the same time, in Eq. (42), which determines the imaginary parts of the eigenfrequencies, $`\stackrel{~}{\omega }_k^{\prime \prime }`$, all the terms in the sum are positive, so that one should keep only the term with $`\omega _j`$ closest to $`\stackrel{~}{\omega }_k^{}`$. Since $`(\stackrel{~}{\omega }_k^{}\omega _j)\mathrm{\Omega }/N`$ for this term, we obtain the following estimate for the width of the subradiant mode: $`\stackrel{~}{\omega }_k^{\prime \prime }\gamma `$, where

$$\gamma \tau \mathrm{\Omega }^2/N^2.$$
(43)   
It can be seen that $`\gamma `$ is much smaller than the MFS (by the factor $`\mathrm{\Omega }\tau /N1`$). This justifies neglecting $`\omega ^{\prime \prime 2}`$ in the denominators of Eqs. (41) and (42).

The superradiant solution of Eqs. (41) and (42) corresponds to the case $`\omega ^{\prime \prime }\mathrm{\Omega }`$. Then we readily obtain $`\stackrel{~}{\omega }^{}=N^1_j\omega _j`$ and $`\stackrel{~}{\omega }^{\prime \prime }=\mathrm{\Gamma }N/\tau `$. We see that, indeed, $`\mathrm{\Gamma }/\mathrm{\Omega }N/\mathrm{\Omega }\tau 1`$. Therefore, the superradiant band in the spectral intensity is not affected by the disorder.

We therefore conclude that cooperative emission is not destroyed by disorder. The spectrum of the system consists of a single superradiant and $`N1`$ subradiant eigenmodes. For large $`N/\mathrm{\Omega }\tau `$, the subradiant modes are well defined, since their spectral widths are much smaller than the MFS.

### B Limit of weak disorder

In this subsection, we address a nontrivial question about the fate of cooperative eigenmodes when the disorder in frequencies vanishes. In this limit, $`\mathrm{\Omega }\tau /N0`$, all oscillator frequencies become equal, i.e. $`\omega _i\omega _0`$. In the absence of cooperative coupling, $`\alpha _{ij}=0`$ ($`ij`$), the eigenfrequencies of the system are those of individual oscillators with the energy width much larger than the MFS, $`1/\tau \mathrm{\Omega }/N`$, so that the spectrum of the system is degenerate.

However, the situation is more complicated in the presence of cooperative coupling, $`\alpha _{ij}0`$. Consider the case of a “point sample”, $`\alpha _{ij}=1`$. In this case, the width of subradiant modes is given by Eq. (43). Important is that although the MFS diminishes with decreasing $`\mathrm{\Omega }`$, the width $`\gamma `$ decreases even faster: $`\gamma /(\mathrm{\Omega }/N)\mathrm{\Omega }\tau /N0`$. In orher words, in the presence of even a very weak disorder, the narrow subradiant peaks do not overlap. Therefore, the cooperative modes remain distinct even though the “bare” oscillator modes were already degenerate.

In the case of general coupling, the width $`\gamma `$ of subradiant modes for small values of $`\mathrm{\Omega }\tau /N`$ will be determined by the fluctuations of $`\alpha _{ij}`$, as we will see below.

### C Fluctuations of $`\alpha _{ij}`$

Let us turn to the case with realistic coupling $`\alpha _{ij}`$. The eigenfrequencies $`\stackrel{~}{\omega }_k`$ should now be determined from Eq. (36), which in component form reads

$`{\displaystyle \frac{1}{\tau }}{\displaystyle \underset{j}{}}{\displaystyle \frac{(\omega _j\omega ^{})(1+\sigma _j^{})\omega ^{\prime \prime }\sigma _j^{\prime \prime }}{(\omega _j\omega ^{})^2+\omega ^{\prime \prime 2}}}=0,`$ (44)   
$`{\displaystyle \frac{1}{\tau }}{\displaystyle \underset{j}{}}{\displaystyle \frac{\omega ^{\prime \prime }(1+\sigma _j^{})+(\omega _j\omega ^{})\sigma _j^{\prime \prime }}{(\omega _j\omega ^{})^2+\omega ^{\prime \prime 2}}}=1,`$ (45)   
with $`\sigma _i(\omega )=\sigma _i^{}(\omega )+i\sigma _i^{\prime \prime }(\omega )`$ satisfying Eq. (40), or in component form

$`\sigma _i^{\prime \prime }+{\displaystyle \frac{1}{\tau }}{\displaystyle \underset{j}{}}(\alpha _{ij}1){\displaystyle \frac{(\omega _j\omega ^{})(1+\sigma _j^{})\omega ^{\prime \prime }\sigma _j^{\prime \prime }}{(\omega _j\omega ^{})^2+\omega ^{\prime \prime 2}}}=0,`$ (46)   
$`\sigma _i^{}{\displaystyle \frac{1}{\tau }}{\displaystyle \underset{j}{}}(\alpha _{ij}1){\displaystyle \frac{\omega ^{\prime \prime }(1+\sigma _j^{})+(\omega _j\omega ^{})\sigma _j^{\prime \prime }}{(\omega _j\omega ^{})^2+\omega ^{\prime \prime 2}}}=0.`$ (47)   
For $`\omega ^{\prime \prime }\mathrm{\Omega }/N`$, the system (44)–(47) can be approximately solved in the same way as for a “point” sample. The corresponding condition will be derived in Section V.

When evaluating the contribution to the lhs of Eq. (45) coming from the first term in the numerator, one should keep only one term in the sum with $`\omega _j`$ closest to $`\stackrel{~}{\omega }_k^{}`$: $`(\omega _j\stackrel{~}{\omega }_k^{})\mathrm{\Omega }/N`$. Then we obtain

$$\stackrel{~}{\omega }_k^{\prime \prime }\frac{\tau \mathrm{\Omega }^2}{N^2}\left(1\frac{1}{\tau }\underset{j}{}\frac{\sigma _j^{\prime \prime }}{\omega _j\stackrel{~}{\omega }_k^{}}\right),$$
(48)   
where we again dropped $`\omega ^{\prime \prime 2}`$ in the denominator. Since $`\sigma _i^{}\sigma _i^{\prime \prime }1`$ (see Section V), the frequencies $`\stackrel{~}{\omega }_k^{}`$ in Eq. (48) are the same as for the case $`\alpha _{ij}=1`$. Finding $`\sigma _i^{\prime \prime }`$ in the first order from Eq. (46), and substituting the result into Eq. (48), we obtain

$$\stackrel{~}{\omega }_k^{\prime \prime }\frac{\tau \mathrm{\Omega }^2}{N^2}\left[1+\frac{1}{\tau ^2}\underset{ij}{}\frac{\alpha _{ij}1}{(\omega _i\stackrel{~}{\omega }_k^{})(\omega _j\stackrel{~}{\omega }_k^{})}\right].$$
(49)   
The second term is the sought correction to the width of the subradiant modes. Remarkably, this term turns to zero if the matrix elements $`\alpha _{ij}`$ are replaced by their average $`\overline{\alpha }`$. Indeed, in this case the double sum in Eq. (49) would factorize into a product of two sums, each vanishing due to the fact that $`\stackrel{~}{\omega }_k^{}`$ are the solutions of Eq. (41) (corresponding to $`\alpha _{ij}=1`$). Therefore, the widths of the subradiant modes are determined by the fluctuations, $`\delta \alpha _{ij}`$, of the coupling parameters $`\alpha _{ij}`$ rather than the deviation of their average, $`\overline{\alpha }`$, from unity. It should be noted that this property is general: one can easily see by comparing Eqs. (44) and (45) to Eqs. (46) and (47) that for $`\alpha _{ij}=const`$, we have $`\sigma _i^{\prime \prime }=0`$ and $`\sigma _i^{}1`$, so that the eigenfrequencies $`\stackrel{~}{\omega }_k`$ are unaffected.

### D Discussion of the numerical results

We are now in the position to explain the spectra shown in Figs. 3–6. For the model coupling: $`\alpha _{ij}=\alpha `$, $`\alpha _{ii}=1`$, fluctuations only in the diagonal elements are finite: $`\delta \alpha _{ij}=(1\alpha )\delta _{ij}`$. Substituting this $`\delta \alpha _{ij}`$ into Eq. (49) \[instead of $`(\alpha _{ij}1)`$\] and keeping only the term with $`(\stackrel{~}{\omega }_k^{}\omega _j)\mathrm{\Omega }/N`$ in the remaining sum, we obtain

$$\gamma \left[\frac{\tau \mathrm{\Omega }^2}{N^2}+\frac{1}{\tau }(1\alpha )\right].$$
(50)   
The above expression indicates that after the cooperative modes have been formed (at $`N\mathrm{\Omega }\tau `$), the system can be found in two different regimes characterized by the relative magnitude of the first and second terms in the rhs. For intermediate number of oscillators, $`\mathrm{\Omega }\tau N\mathrm{\Omega }\tau (1\alpha )^{1/2}`$, the width decreaeses with increasing $`N`$, as can be seen by comparing the bottom second and third curves in each of Figs. 3–6 (note that the lowest curves with $`N=2`$ show no sign of cooperative emission). In this regime, the system behaves in the same way as a “point” sample. With increasing number of oscillators, the dependence on $`N`$ saturates, and the width is dominated by the fluctuations of $`\alpha _{ij}`$. Correspondingly, the change in the pattern of the peaks in Figs. 3–6, calculated for different values of $`(1\alpha )`$, occurs at different $`N`$, as can be seen by comparing the next two curves in each figure. Note, however, that with further increase in $`N`$, the curves exhibit yet another change in pattern. Namely, the peaks get smeared out (the top two curves in each figure). This occurs when the value of $`(1\alpha )/\tau `$ exceeds the MFS, $`\mathrm{\Omega }/N`$, which is inconsistent with the above analysis. The reason for such a discrepancy is that for large $`N`$, the model with coupling $`\alpha _{ij}`$ independent of separation between oscillators, becomes inadequate, as we mentioned above. For the correct description of the peaks smearing at large $`N`$, the spatial dependence of $`\alpha _{ij}`$ is crucial; this question is addressed in Section VI. Nevertheless, for $`N\mathrm{\Omega }\tau /(1\alpha )`$, this model describes accurately the mesoscopic features of the spectral intensity, as shown in the next section.

## V Strong Mesoscopics Regime

Let us now estimate the typical width of the radiation eigenmodes due to the fluctuations in $`\alpha _{ij}`$. Since the configurational average of the second term in Eq. (49) \[with $`\delta \alpha _{ij}`$ instead of $`(\alpha _{ij}1)`$\] vanishes, we need to evaluate $`(\stackrel{~}{\omega }_k^{\prime \prime })^2`$. Using the fact that only diagonal terms in the average $`\delta \alpha _{ij}\delta \alpha _{i^{}j^{}}`$ survive and omitting the first term in Eq. (49), we write

$$(\stackrel{~}{\omega }_k^{\prime \prime })^2\left(\frac{\mathrm{\Omega }^2}{\tau N^2}\right)^2\underset{ij}{}\frac{(\delta \alpha _{ij})^2}{(\omega _i\stackrel{~}{\omega }_k^{})^2(\omega _j\stackrel{~}{\omega }_k^{})^2}.$$
(51)   
The sum is dominated by the terms with $`(\omega _i\stackrel{~}{\omega }_k^{})(\omega _j\stackrel{~}{\omega }_k^{})\mathrm{\Omega }/N`$. Since the typical spatial separation between two oscillators with close frequencies is $`L`$, the separation fluctuations are of the same order. Thus, the typical fluctuation of $`\alpha _{ij}`$ is $`\delta \alpha \sqrt{(\delta \alpha _{ij})^2}\left(L/\lambda _0\right)^2`$, and we finally obtain the typical width of a subradiant mode, $`\gamma \sqrt{(\stackrel{~}{\omega }_k^{\prime \prime })^2}`$, as

$$\gamma \frac{\delta \alpha }{\tau }\frac{1}{\tau }\left(\frac{L}{\lambda _0}\right)^2.$$
(52)   
Comparing Eq. (52) to Eq. (49), we see that fluctuations in $`\alpha _{ij}`$ dominate the width $`\gamma `$ for $`N\mathrm{\Omega }\tau (\lambda _0/L)`$.

In order to characterize the fine structure in the emission spectrum, it is convenient to introduce the dimensionless parameter

$$\kappa =\frac{\mathrm{\Omega }\tau }{N}\left(\frac{\lambda _0}{L}\right)^2.$$
(53)   
It represents the product of a small, $`\mathrm{\Omega }\tau /N`$, and a large, $`(\lambda _0/L)^2`$, factors, which characterize the disorder and the system size, respectively. In terms of $`\kappa `$, the condition for the formation of the cooperative modes, $`N/\mathrm{\Omega }\tau 1`$, can be presented as $`\kappa \left(\lambda _0/L\right)^2`$.

Using Eq. (53), the width (52) can be expressed in terms of the MFS as

$$\gamma \frac{1}{\kappa }\left(\frac{\mathrm{\Omega }}{N}\right).$$
(54)   
This result applies when $`\kappa \lambda _0/L`$. On the other hand, it was implicit in the above derivation (Section IV) that typical $`\sigma _i^{}`$ and $`\sigma _i^{\prime \prime }`$ are smaller than unity . The latter parameters can be estimated in a similar way from Eqs. (46) and (47) with the result: $`\sigma ^{\prime \prime }\kappa ^1`$ and $`\sigma ^{}\sigma ^{\prime \prime 2}\kappa ^2`$. Thus, the lower boundary for $`\kappa `$, at which Eq. (54) applies, is $`\kappa 1`$. For $`\kappa 1`$, all terms in Eqs. (44)–(47) become of the same order of magnitude, and for smaller $`\kappa `$ this system has no subradiant solutions, as discussed above.

Since the MFS exceeds the width $`\gamma `$ within the entire domain $`1\kappa \lambda _0/L`$, the fine structure in the spectral intensity $`I(\omega )`$ is well pronounced. In other words, this domain corresponds to the strong mesoscopics regime. The opposite case $`\kappa 1`$ is considered in the next section.

## VI Weak Mesoscopics Regime

In the domain $`\kappa 1`$, the system cannot sustain eigenmodes that involve all $`N`$ oscillators. As a result, the eigenmodes become localized, in the sense that each eigenmode would comprise some $`N_cN`$ oscillators and occupy the volume with characteristic size $`L_cL`$. The magnitude of $`L_c`$ and $`N_c`$ can be estimated from the following argument. Let us divide the system of oscillators into subsystems of increasingly smaller size. When the size of the subsystem becomes $`L_c`$, the system of equations (44) and (45), applied to a subsystem, first acquires a solution. This happens when the width $`\gamma _c\tau ^1(L_c/\lambda _0)^2`$, determined from Eq. (52) for a subsystem, becomes of the order of MFS within a subsystem, i.e.

$$\frac{1}{\tau }\left(\frac{L_c}{\lambda _0}\right)^2\frac{\mathrm{\Omega }}{N_c}.$$
(55)   
Taking into account that $`N_c=N(L_c/L)^3`$, we find

$$L_c\kappa ^{1/5}L,N_c\kappa ^{3/5}N.$$
(56)   
Substituting these results back into Eq. (55), we find for the eigenmodes width

$$\gamma =\gamma _c\frac{1}{\kappa ^{3/5}}\left(\frac{\mathrm{\Omega }}{N}\right)\frac{\mathrm{\Omega }}{N}.$$
(57)   
From Eq. (56), we can also estimate how the relative amplitude of mesoscopic fluctuations in the spectral intensity $`I(\omega )`$ falls off with decreasing $`\kappa `$:

$$\frac{\delta I}{I}\left(\frac{N_c}{N}\right)^{1/2}=\kappa ^{3/10}.$$
(58)   
It is apparent that smearing of the fine structure in the cooperative emission spectrum with decreasing $`\kappa `$ occurs rather slowly.

## VII Dipole-dipole interactions

In this section we study the effect of dipole-dipole interactions on the cooperative emission from a disordered system. Note that for a “point” sample with $`L\lambda _0`$, the typical magnitude of the dipole–dipole interaction between two oscillators is much larger than their superradiant coupling, $`\beta _{ij}/\alpha _{ij}(\lambda _0/L)^31`$. The structure of the eigenmodes in the absence of superradiant coupling, given by Eq. (20) with $`\alpha _{ij}=0`$, was considered in several papers. Renormalization–group arguments of Ref.  (see also Ref. ) suggest that all eigenmodes are delocalized. Numerical studies indicate a wide range of spatial scales in eigenmodes and, thus, seem to support this conclusion. In Ref. , the role of general random-matrix perturbation in the spectrum of multilevel system was studied analytically; he ensemble-averaged renormalization of the spectrum of the system was derived which does not capture, however, the mesoscopic effects. Below we argue that finite disorder in combination with superradiant coupling lead to a certain “resistance” of the system to large, but zero on average, dipole–dipole terms because of the formation of cooperative modes.

In the absence of superradiant coupling ($`\alpha _{ij}=0`$) the dipole–dipole interactions lead to the shifts in the frequencies of individual oscillators. The resulting additional spread in $`\omega _i`$ is, in general, much larger than the “bare” spread $`\mathrm{\Omega }`$. This can be readily seen from the lowest–order correction to the frequency, $`\delta \omega _i`$, which has the form

$$\delta \omega _i=\frac{1}{\tau ^2}\underset{ji}{}\frac{\beta _{ij}^2}{\omega _i\omega _j},$$
(59)   
(since $`\beta _{ii}=0`$, the lowest–order correction to $`\omega _i`$ is quadratic). The main contribution to the sum comes from pairs of oscillators located closely in space, with $`r_{ij}LN^{1/3}`$ (nearest neighbor interaction), so that

$$\beta _{ij}N\left(\frac{\lambda _0}{L}\right)^3.$$
(60)   
Since the typical frequency difference for such pairs is $`\mathrm{\Omega }`$, we obtain

$$\delta \omega _i\frac{N^2}{\mathrm{\Omega }\tau ^2}\left(\frac{\lambda _0}{L}\right)^6.$$
(61)   
On the other hand, the results obtained in the previous sections apply only if the additional disorder, caused by dipole–dipole interactions, does not affect the MFS. This requires the condition $`\delta \omega _i\mathrm{\Omega }`$ to be met. Using Eq. (61), this condition could be rewritten as $`\lambda _0/L(\mathrm{\Omega }\tau /N)^{1/3}`$. Since the formation of cooperative modes occurs only if $`\mathrm{\Omega }\tau /N1`$, one could draw the conclusion that neglecting the dipole–dipole interactions would be inconsistent with our basic assumption $`L\lambda _0`$.

The resolution of this apparent contradiction lies in the observation that, in the presence of superradiant coupling, i.e. $`\alpha _{ij}0`$, the true eigenmodes of the system are cooperative modes comprised from a large number of oscillators. Therefore, the relevant condition should involve the shifts, $`\delta \stackrel{~}{\omega }_k^{}`$, of the eigenmodes frequencies, rather than $`\delta \omega _i`$. In the first order, $`\delta \stackrel{~}{\omega }_k^{}`$ is given by an expression similar to the second term in the rhs of Eq. (49) \[with $`\beta _{ij}`$ instead of $`(\alpha _{ij}1)`$\]. Since this term vanishes on average, as discussed above, the typical shift, $`\delta \stackrel{~}{\omega }^{}\sqrt{(\delta \stackrel{~}{\omega }_k^{})^2}`$, can be estimated from \[compare with Eq. (51)\]

$$(\delta \stackrel{~}{\omega }_k^{})^2\left(\frac{\mathrm{\Omega }^2}{\tau N^2}\right)^2\underset{ij}{}\frac{(\beta _{ij})^2}{(\omega _i\stackrel{~}{\omega }_k^{})^2(\omega _j\stackrel{~}{\omega }_k^{})^2}.$$
(62)   
There are two main contributions to the sum in the rhs. The first comes from the nearest–neighbor interaction with $`\beta _{ij}`$ given by Eq. (60). The second contribution originates from the pairs ($`ij`$) which are close in frequency; for such pairs, $`\beta _{ij}(\lambda _0/L)^3`$. Both contributions turn out to be of the same order of magnitude, resulting in

$$\delta \stackrel{~}{\omega }^{}\frac{1}{\tau }\left(\frac{\lambda _0}{L}\right)^3.$$
(63)   
This result is smaller than $`\delta \omega _i`$ in Eq. (61) by the factor $`N(N/\mathrm{\Omega }\tau )(\lambda _0/L)^31`$. Such a dramatic difference illustrates the “resistance” of a coupled system of oscillators with disorder in frequencies to dipole–dipole interactions, as mentioned above. This property can also be qualitatively explained as follows. The dipole–dipole interaction between two subradiant modes can be viewed as an interaction between a mode and the electric field, $`\stackrel{~}{𝐄}(𝐫)`$, created by the dipole moments of oscillators making up the other mode. Since the number of oscillators in a mode is large, their electric fields effectively cancel each other, so that the resulting net field, $`\stackrel{~}{𝐄}(𝐫)`$, varies in space much slower than those of the individual oscillators. Note now that a slowly varying electric field couples only weakly to a subradiant mode. In fact, the suppression of the dipole–dipole interaction between subradiant modes has the same physical origin as their decoupling from the radiation field: had the electric field $`\stackrel{~}{𝐄}`$ been uniform, the cooperative modes would not interact at all with each other. This is the reason why the corrections to $`\stackrel{~}{\omega }_k^{}`$ and $`\stackrel{~}{\omega }_k^{\prime \prime }`$ vanish on average, and consequently the typical $`\delta \stackrel{~}{\omega }^{}`$ and $`\gamma `$ are determined by the fluctuations of $`S_{ij}`$. In contrast, the frequency shifts of individual oscillators are due to their interactions with the nearest neighbors, so that no cancellations occur.

Thus we arrive at the condition $`\lambda _0/L\left(\mathrm{\Omega }\tau \right)^{1/3}`$, or, in terms of the parameter $`\kappa `$,

$$\kappa \frac{1}{N}\left(\frac{\lambda _0}{L}\right)^5.$$
(64)   
This condition should be consistent with the condition for the formation of the cooperative modes, $`\kappa (\lambda _0/L)^2`$. We see that both conditions are satisfied for sufficiently large $`N`$, i.e. $`N(\lambda _0/L)^3`$. To account for different mesoscopics regimes, it is convenient to present Eq. (64) in the form

$$N\left(\frac{\lambda _0}{L}\right)^n.$$
(65)   
Then $`n=3`$, 4, and 5 correspond to the “point” sample, strong mesoscopics, and weak mesoscopics regimes, respectively.

## VIII Conclusions

The main result of the present paper is that disorder in oscillators frequencies does not destroy the cooperative character of the emission from a “point” sample, as long as the MFS, $`\mathrm{\Omega }/N`$, is smaller than the linewidth of an individual oscillator, $`\tau ^1`$. In the opposite case, when $`\mathrm{\Omega }/N\tau ^1`$, the spectrum represents a system of non–overlapping Lorentzians with the width $`\tau ^1`$.

It is convenient to characterize the disorder in terms of the dimensionless parameter $`\kappa =\mathrm{\Omega }\tau \lambda _0^2/NL^2`$. Below we summarize our results for the characteristic width, $`\gamma `$, of the subradiant peaks (in units of $`\mathrm{\Omega }/N`$) for different domains of $`\kappa `$:

$$\gamma =\frac{\mathrm{\Omega }}{N}\mathrm{\Phi }(\kappa ,L/\lambda _0),$$
(66)   
where the dimensionless function $`\mathrm{\Phi }`$ has the following asymptotes

$`\mathrm{\Phi }`$ $`=`$ $`\kappa \left({\displaystyle \frac{L}{\lambda _0}}\right)^2,\mathrm{for}{\displaystyle \frac{\lambda _0}{L}}\kappa \left({\displaystyle \frac{\lambda _0}{L}}\right)^2,\text{“point” sample},`$ (67)   
$`\mathrm{\Phi }`$ $`=`$ $`\kappa ^1,\mathrm{for}\mathrm{\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}1}\kappa {\displaystyle \frac{\lambda _0}{L}},\text{strong mesoscopics},`$ (68)   
$`\mathrm{\Phi }`$ $`=`$ $`\kappa ^{3/5},\mathrm{for}\kappa 1,\text{weak mesoscopics}.`$ (69)   
For $`\kappa (\lambda _0/L)^2`$, the spectrum corresponds to uncoupled oscillators.

In Section III, we presented an exact solution of a model with simplified (separable) coupling Eq. (21). This model describes accurately the first two (“point” sample and strong mesoscopics) regimes in Eq. (67). It becomes, however, inadequate in the third (weak mesoscopics) regime, giving a $`\kappa ^1`$ instead of the correct $`\kappa ^{3/5}`$ dependence for the period of mesoscopic structure in the cooperative emission spectrum.

Throughout the paper we have considered a three-dimensional system of oscillators. When the oscillators are confined to a plane, only the results for $`\kappa 1`$ should be modified. In this case, repeating the consideration of Section VI, we obtain $`\mathrm{\Phi }=\kappa ^{1/2}`$. Also for the relative magnitude of mesoscopic fluctuations, $`\left(\delta I/I\right)`$, instead of Eq. (58) we obtain $`\left(\delta I/I\right)\kappa ^{1/4}`$.

Note finally that in experiments, such as photoexcited excitons in polymer films, the number of oscillators $`N`$ is governed by the excitation intensity. Thus, for a given disorder, the crossover from the strong mesoscopics regime ($`\kappa >1`$) to the weak mesoscopics regime ($`\kappa <1`$) can be simply achieved by increasing the excitation intensity level.

###### Acknowledgements.

The authors are grateful to S. Mukamel and M. I. Stockman for helpful discussions. The work at Vanderbilt University was supported by NSF grant ECS-9703453, and the work at the University of Utah was supported by NSF grant DMR-9732820. M.E.R. was also supported by the Petroleum Research Fund under grant ACS-PRF# 34302-AC6.    

FIG. 1

FIG. 2

FIG. 3

FIG. 4

FIG. 5

FIG. 6