Cooperative scattering of scalar waves by optimized configurations of point scatterers

TL;DR

This study explores how arranging point scatterers optimally can significantly enhance scattering or create long-lived resonances, with configurations mostly aligned linearly, revealing potential for quantum memory and mirror applications.

Contribution

It introduces numerical optimization of scatterer positions to maximize scattering or minimize decay, revealing linear arrangements as optimal configurations for scalar wave scattering.

Findings

Maximal scattering cross section scales quadratically with number of scatterers.

Minimal decay rate decreases exponentially with number of scatterers.

Linear configurations are optimal for both maximizing scattering and minimizing decay.

Abstract

We investigate multiple scattering of scalar waves by an ensemble of $N$ resonant point scatterers in three dimensions. For up to $N = 21$ scatterers, we numerically optimize the positions of the individual scatterers, such as to maximize the total scattering cross section for an incoming plane wave, on the one hand, and to minimize the decay rate associated to a long-lived scattering resonance, on the other hand. In both cases, the optimimum is achieved by configurations where all scatterers are placed on a line parallel to the direction of the incoming plane wave. The associated maximal scattering cross section increases quadratically with the number of scatterers for large $N$, whereas the minimal decay rate -- which is realized by configurations that are not the same as those that maximize the scattering cross section -- decreases exponentially as a function of $N$. Finally, we also analyze the stability of our optimized configurations with respect to small random displacements of the scatterers. These results demonstrate that optimized configurations of scatterers bear a considerable potential for applications such as quantum memories or mirrors consisting of only a few atoms.

Tables

Table 1. Table 1: Positions in units of k − 1 superscript 𝑘 1 k^{-1} , scattering cross sections σ 𝜎 \sigma in units of the single scatterer cross section σ max ( 1 ) = 4 ​ π / k 2 subscript superscript 𝜎 1 max 4 𝜋 superscript 𝑘 2 \sigma^{(1)}_{\rm max}=4\pi/k^{2} and the smallest decay rates γ m ​ i ​ n subscript 𝛾 𝑚 𝑖 𝑛 \gamma_{min} in units of the decay rate γ 𝛾 \gamma of a single scatterer of the optimized configurations for eight and nine scatterers. The three different types of configurations (narrow, wide and decay) correspond to the ones shown in Fig. 2 (narrow), Fig. 3 (wide) and Fig. 7 (decay).

$N$	type	$k{z}_1,..,k{z}_N$	$\sigma$/${\sigma }_{\max }^{\left(1\right)}$	${\gamma }_{min}/\gamma$
8	narrow	-8.341, -5.411, -2.908, -0.830, 0.830, 2.908, 5.411, 8.341	23.6	$6.8\cdot {10}^{-3}$
8	wide	-8.452, -5.667, -3.453, -2.004, 2.004, 3.453, 5.667, 8.452	23.9	$1.7\cdot {10}^{-2}$
8	decay	-3.32587, -1.92458, -0.95543, -0.25, 0.25, 0.95543, 1.92458, 3.32587	1.3	$6.2\cdot {10}^{-10}$
9	narrow	-8.340, -5.458, -2.960, -0.843, 0.843, 2.851, 5.273, 8.100, 11.302	26.1	$4.4\cdot {10}^{-3}$
9	wide	-8.423, -5.679, -3.444, -1.989, 1.989, 3.460, 5.599, 8.301, 11.479	26.2	$1.2\cdot {10}^{-2}$
9	decay	-3.31104, -1.91837, -0.95357, -0.25, 0.25, 0.95338, 1.91681, 3.30431, 7.26926	1.7	$4.6\cdot {10}^{-10}$

Equations

$${\vec{\nabla }}^2\Psi (\mathbf{r})+k^2[1+\chi (\mathbf{r})]\Psi (\mathbf{r})=0$$

$$\tilde{\Psi }(\mathbf{r},t)=\Psi (\mathbf{r})e^{-i\omega t}$$

$${\Psi }_0(\mathbf{r})={\Psi }_0{e}^{i{\mathbf{k}}_{\mathrm{in}}\cdot \mathbf{r}}$$

$$\Psi (\text{bi}r)$$

$${\sigma }^{(1)}=4\pi |f{|}^2$$

$$k|f{|}^2=\mathrm{Im}\left[f\right].$$

$$f=\frac{1}{k(\alpha -i)},\alpha \in \mathbb{R}.$$

$${\sigma }_{\max }^{(1)}=\frac{4\pi }{k^2}.$$

$$f=\frac{1}{k\left(\frac{\delta }{\gamma }-i\right)}$$

$${\Psi }_i={\Psi }_0({\mathbf{r}}_i)+k\sum _{j=1}^N{G}_{ij}{f}_j{\Psi }_j$$

$$G_{ij}=\cases{0&if $i=j$\\ \frac{e^{ik|{\bf r}_{i}-{\bf r}_{j}|}}{k|{\bf r}_{i}-{\bf r}_{j}|}&if $i\neq j$\\ }$$

$$\Psi (\mathbf{r})={\Psi }_0(\mathbf{r})+\sum _{i=1}^N\frac{e^{ik|\mathbf{r}-{\mathbf{r}}_i|}}{|\mathbf{r}-{\mathbf{r}}_i|}{f}_i{\Psi }_i$$

$${\Psi }_{\mathrm{sc}}(\mathbf{R})=\frac{e^{ikR}}{R}\sum _{i=1}^N{e}^{-i{\mathbf{k}}_{\mathrm{out}}\cdot {\mathbf{r}}_i}{f}_i{\Psi }_i$$

$$\sigma =\int \mathrm{d}\Omega R^2\frac{|{\Psi }_{\mathrm{sc}}(\mathbf{R}){|}^2}{|{\Psi }_0{|}^2}=\int \mathrm{d}\Omega {\left|\sum _{i=1}^N{e}^{-i{\mathbf{k}}_{\mathrm{out}}\cdot {\mathbf{r}}_i}{f}_i\frac{{\Psi }_i}{{\Psi }_0}\right|}^2$$

$$\int \mathrm{d}\Omega e^{-i{\mathbf{k}}_{\mathrm{out}}\cdot ({\mathbf{r}}_i-{\mathbf{r}}_j)}=4\pi \left(\frac{G_{ij}-G_{ij}^{\ast }}{2i}+{\delta }_{ij}\right)$$

$$\sigma$$

$$|\Psi ⟩=(\text{mathbbm}1-kfG{)}^{-1}|{\Psi }_0⟩$$

$$\sigma (\delta )=\frac{4\pi }{k^2|{\Psi }_0{|}^2}\mathrm{Im}\left[⟨{\mathbf{\Psi }}_0\left|{\left(\frac{\delta }{\gamma }-i-G\right)}^{-1}\right|{\mathbf{\Psi }}_0⟩\right]$$

$${\sigma }_n(\delta )=\frac{4\pi \gamma }{k^2|{\Psi }_0{|}^2}\mathrm{Im}\left[\frac{⟨{\mathbf{\Psi }}_0|{\lambda }_{\mathbf{n}}^{(\mathbf{R})}⟩⟨{\lambda }_{\mathbf{n}}^{(\mathbf{L})}|{\mathbf{\Psi }}_0⟩}{\delta -{\delta }_n-i{\gamma }_n}\right]$$

$${\delta }_n$$

$${\gamma }_n$$

$$\sigma (\delta )=\sum _{n=1}^N{\sigma }_n(\delta )$$

$${\gamma }_n$$

$$\frac{\sin (k{r}_{ij})}{k{r}_{ij}}=\frac{1}{2}{\int }_{-1}^1\mathrm{d}\cos (\theta )e^{i\cos (\theta )k{r}_{ij}}=\int \frac{\mathrm{d}\Omega }{4\pi }{e}^{i({\mathbf{r}}_i-{\mathbf{r}}_j)\cdot {\mathbf{k}}_{\Omega }}$$

$$g_n(\mathbf{k})=\sum _{i=1}^N{c}_i^{(n)}{e}^{-i{\mathbf{r}}_i\cdot \mathbf{k}}$$

$${\sigma }_n\le \frac{4\pi \gamma }{k^2{\gamma }_n|{\Psi }_0{|}^2}\left|⟨{\mathbf{\Psi }}_0|{\lambda }_{\mathbf{n}}^{(\mathbf{R})}⟩⟨{\lambda }_{\mathbf{n}}^{(\mathbf{L})}|{\mathbf{\Psi }}_0⟩\right|$$

$$\frac{⟨{\mathbf{\Psi }}_0|{\lambda }_{\mathbf{n}}^{(\mathbf{R})}⟩⟨{\lambda }_{\mathbf{n}}^{(\mathbf{L})}|{\mathbf{\Psi }}_0⟩}{|{\Psi }_0{|}^2}=\sum _{ij}{c}_i^{(n)}{c}_j^{(n)}{e}^{i{\mathbf{k}}_{\mathrm{in}}\cdot ({\mathbf{r}}_i-{\mathbf{r}}_j)}=g_n({\mathbf{k}}_{\mathrm{in}})g_n(-{\mathbf{k}}_{\mathrm{in}})$$

$$\frac{{\sigma }_n}{{\sigma }_{\max }^{(1)}}\le \left|\frac{g_n({\mathbf{k}}_{\mathrm{in}})g_n(-{\mathbf{k}}_{\mathrm{in}})}{{\sum }_i(c_i^{(n)}{)}^2}\right|\frac{{\sum }_i|c_i^{(n)}{|}^2}{\int \frac{\mathrm{d}\Omega }{4\pi }|g_n({\mathbf{k}}_{\Omega }){|}^2}$$

Full text

Keywords: multiple scattering theory, point scatterers, cooperative scattering, subradiance

1 Introduction

In general, a collection of $N$ point scatterers, with distances between the scatterers of the order of the wavelength, scatters an incoming coherent wave in a different way than $N$ independent scatterers which are placed far away from each other. An obvious example of such ‘cooperative scattering’ is observed in a periodic crystal of scatterers, where the fields emitted by each individual scatterer interfere contructively with each other in certain directions and destructively in other directions (Bragg scattering [1]). Recently, Bragg scattering has been exploited to realize strong reflection of light from only one or two thousand atoms aligned along an optical fiber [2, 3].

Apart from interference between singly scattered waves, cooperative effects are also induced by multiple scattering, where waves emitted by individual scatterers are again scattered by other scatterers, and so on. If the wave is scattered back and forth many times between a few scatterers, this may lead to a strong enhancement of the local field intensity in the vicinity of these scatterers. These local field enhancements are interesting for technological applications, since they can be used to increase the efficiency of solar cells or other optical devices [4, 5], and have been observed in specifically tailored nanostructured materials consisting of, e.g. metallic nanoantennas [6], nanospheres [7] or nanoparticles [8].

In the present article, we will investigate scalar point scatterers as a paradigmatic model system for multiple scattering. A single point scatterer exhibits a maximum scattering cross section of $\sigma_{\rm max}^{(1)}=\lambda^{2}/\pi$, where $\lambda$ denotes the wavelength of the scattered light [9]. The fundamental questions which we address in the present article are: what is the largest scattering cross section that can be achieved with an ensemble of $N$ point scatterers, and how must the scatterers be positioned in space in order to realize this maximum? In particular, we will show that, by choosing optimized configurations of point scatterers, it is possible to exploit effects of cooperative scattering such that the total scattering cross section exceeds the value $N\sigma_{\rm max}^{(1)}$ for $N$ independent scatterers. Expressing the total scattering cross section as a sum over scattering resonances, we furthermore demonstrate that a necessary (but not sufficient) condition to achieve a large scattering cross section is the existence of a narrow scattering resonance with correspondingly small decay rate. In quantum optics, such states are known as ‘subradiant states’ or ‘dark states’ [10], and have recently been observed in a large cloud of cold atoms [11].

The article is structured as follows: After briefly outlining the theoretical frame (see chapter 2), we will describe the numerical optimization algorithm and present our results for optimized configurations which maximize the total scattering cross section, on the one hand, and those which minimize the smallest decay rate, on the other hand (chapter 3). Furthermore, we will analyze the stability of our optimized configurations with respect to small random displacements of the scatterers (chapter 4). The last section concludes this article and gives a short outlook on future perspectives.

2 Theoretical frame

In this section, we review basic theoretical concepts of scattering from a collection of point scatterers. After introducing fundamental equations describing scattering of a monochromatic, scalar wave by a single point scatterer, we treat the case of many point scatterers based on the Foldy-Lax formalism for multiple scattering [12, 13]. We derive the total scattering cross section for $N$ scatterers and express it in terms of scattering resonances, where the widths of these resonances correspond to the decay rates of the corresponding resonance states.

2.1 Scalar wave equation

We consider the following time-independent equation for a scalar wave $\Psi({\bf r})$ in three dimensional real space:

[TABLE]

with wavenumber $k$ and a spatial inhomogeneity $\chi({\bf r})$ that induces scattering. For the case of a quantum particle obeying the Schrödinger equation, for example, the spatial inhomogeneity is given by $\chi({\bf r})=-2mV({\bf r})/(\hbar^{2}k^{2})$ in terms of the scattering potential $V({\bf r})$ and the mass $m$. Solutions of the corresponding time-dependent wave equation are given by

[TABLE]

where the frequency $\omega$ is determined by $k$ through the dispersion relation (e.g. $\omega=\frac{\hbar k^{2}}{2m}$ in case of the Schrödinger equation). In the following, however, we will restrict ourselves to the stationary regime described by Eq. (1).

Electromagnetic waves propagating in dielectric media with space-dependent electric susceptibility $\chi({\bf r})$ are also described by Eq. (1) in a scalar approximation (i.e. ignoring the vectorial character of the electric field) which is often used in the theory of multiple scattering of light (e.g. [14]). Nevertheless, scalar waves and vectorial waves may behave differently, especially if the distances between scatterers are smalll [15]. Therefore, the optimal configurations of point scatterers investigated in this article apply, in principle, only to scalar waves, although we expect the optimal configurations for vector waves to have similar qualitative properties. In the vectorial case, cooperative effects of light scattering in one-dimensional atomic arrays have been studied recently [16], but only for periodic arrays (i.e. without optimizing the individual positions of the scatterers).

2.2 Single point scatterer

A point scatterer at position ${\bf r}_{0}$ is defined by the condition that the inhomogeneity $\chi({\bf r})$ in Eq. (1) is different from zero only in a small volume containing the point ${\bf r}_{0}$, i.e. if $|{\bf r}-{\bf r}_{0}|\ll 1/k$. The solution of Eq. (1) with incoming plane wave

[TABLE]

where ${\bf k}_{\rm in}$ with $|{\bf k}_{\rm in}|=k$ denotes the incoming wave vector, is then given by [18]:

[TABLE]

where $e^{ik|\bi{r}-\bi{r}_{0}|}/(|\bi{r}-\bi{r}_{0}|)$ represents a spherical wave at $\bi{r}$ originating from point $\bi{r}_{0}$, and $f$ denotes the scattering amplitude. The latter is isotropically distributed (i.e. exhibits no angular dependence) due to the negligible extent of the scatterer. In principle, the scattering amplitude $f$ can be calculated for a given $\chi({\bf r})$ through the solution of the Lippmann-Schwinger equation [17]. In the following, however, we will not be concerned with the actual physical realization of the point scatterer, but rather assume that its properties are given to us in terms of its scattering amplitude $f$.

The scattering cross section of a single point scatterer is given by:

[TABLE]

Moreover, the scattering amplitude obeys the optical theorem which guarantees flux conservation [19]:

[TABLE]

The general solution of equation (6) (for $f\neq 0$) can be written in the following form:

[TABLE]

The optical theorem thus yields a fundamental upper limit for the scattering cross section of a point scatterer, which is independent of its physical realization. This maximum is achieved by choosing $\alpha=0$ in (7), leading to:

[TABLE]

In general, the real parameter $\alpha$ determining the scattering amplitude (7) depends on the wavenumber $k$ or, equivalently (due to the dispersion relation), on the frequency $\omega$ of the scattered wave. The condition $\alpha(\omega_{0})=0$ then corresponds to a resonance condition with resonance frequency $\omega_{0}$. For small deviations $\delta=\omega-\omega_{0}$ from the resonance, we can expand $\alpha$ in first order around $\omega_{0}$, such that

[TABLE]

where $\gamma$ defines the width of the resonance.

2.3 Many point scatterers

For $N$ point scatterers placed at positions ${\bf r}_{1},\dots,{\bf r}_{N}$, the following set of equations for the fields $\Psi_{i}$ incident on scatterer $i$ have been derived [12, 13]:

[TABLE]

where $f_{i}$ denotes the scattering amplitude of scatterer $i$, and the Green matrix $G$ is defined as

[TABLE]

According to (10) and (11), the field at scatterer $i$ is obtained as the incoming wave $\Psi_{0}({\bf r}_{i})$ plus the fields scattered from all other scatterers $j\neq i$. Similarly, the total wave function $\Psi({\bf r})$ (at any position ${\bf r}$ different from ${\bf r}_{i}$, $i=1,\dots,N$) is obtained as:

[TABLE]

In order to determine the scattering cross section $\sigma$, we calculate the flux of the scattered part $\Psi_{\rm sc}({\bf r}):=\Psi({\bf R})-\Psi_{0}({\bf r})$ of the wave through a sphere with large radius $R$. By definition, this flux is equal to the flux of the incoming wave through the area $\sigma$. Applying a far-field approximation ($|{\bf R}|\gg|{\bf r}_{i}|$, $i=1,\dots,N$), we obtain:

[TABLE]

where ${\bf k}_{\rm out}=k{\bf R}/|{\bf R}|$ denotes the wavevector of the outgoing wave. Since, in the far field, the flux density ${\bf J}_{\rm sc}=-\frac{1}{k}\frac{{\rm d}\omega}{{\rm d}k}{\rm Im}[\Psi_{\rm sc}\vec{\nabla}\Psi_{\rm sc}^{*}]$ points in the radial direction, the total scattering cross section is obtained by integrating the wave intensity (normalized by the incoming intensity) over the surface of a sphere with radius $R$:

[TABLE]

where $\Omega$ represents the angular variables defining the direction of ${\bf R}$ (and of ${\bf k}_{\rm out}$). Using

[TABLE]

see equation (24), we obtain:

[TABLE]

where we used (6) and (10) in the second line. Equation (16) expresses the optical theorem for the $N$-scatterer system: the total scattering cross section is associated with destructive interference of the waves emitted by scatterers $i=1,\dots,N$ with the incoming wave, or, in simpler words, due to flux conservation, the scattered flux is taken away from the incoming flux.

2.4 Scattering resonances

From now on, we assume that all $N$ scatterers are identical, i.e. $f_{i}=f$ for all $i=1,\dots,N$. Defining $N$-dimensional vectors $|{\bf\Psi}\rangle$ and $|{\bf\Psi_{0}}\rangle$ with components $\langle i|{\bf\Psi}\rangle:=\Psi_{i}$ and $\langle i|{\bf\Psi_{0}}\rangle:=\Psi_{0}({\bf r}_{i})$, a formal solution of (10) is obtained as follows:

[TABLE]

Using the explicit form (9) of the scattering amplitude, we can rewrite the scattering cross section (16) as a function of the detuning $\delta$ as follows:

[TABLE]

Since $\omega=\omega_{0}+\delta$ is related to $k$ through the dispersion relation, we must be aware of the fact that $k$ (and, consequently also $G$) depends on $\delta$. These dependencies can be neglected, however, if we assume $|\delta|\ll k\left|\frac{{\rm d}\omega}{{\rm d}k}\right|$ and $|\delta|\ll\left|\frac{{\rm d}\omega}{{\rm d}k}\right|/r_{ij}$ for all distances $r_{ij}$. This, in turn, is justified in the case of a narrow single-scatterer resonance (i.e. if $\gamma\ll k\left|\frac{{\rm d}\omega}{{\rm d}k}\right|$ and $\gamma\ll\left|\frac{{\rm d}\omega}{{\rm d}k}\right|/r_{ij}$), since the scattering cross section is negligibly small for $|\delta|\gg\gamma$.

From (18), we see that each eigenvalue $\lambda_{n}$ of $G$, with corresponding right-eigenvector $|{\bf\lambda_{n}^{(R)}}\rangle$ and left-eigenvector $\langle{\bf\lambda_{n}^{(L)}}|$ normalized such that $\langle{\bf\lambda_{n}^{(L)}}|{\bf\lambda_{n}^{(R)}}\rangle=1$, gives rise to a scattering resonance:

[TABLE]

where the position $\delta_{n}$ and the width $\gamma_{n}$ are determined by the real and imaginary part of the eigenvalue $\lambda_{n}$ as follows:

[TABLE]

If the matrix $G$ is diagonalisable – which is the case for all the optimized configurations discussed below (see also the last paragraph of A) – the total scattering cross section results as the sum over all $N$ resonances:

[TABLE]

Since, in general, the left- and right-eigenvectors differ from each other (see A), the factor $\langle{\bf\Psi_{0}}|{\bf\lambda_{n}^{(R)}}\rangle\langle{\bf\lambda_{n}^{(L)}}|{\bf\Psi_{0}}\rangle$ appearing in Eq. (19) may exhibit a non-vanishing imaginary part leading to an asymmetric, Fano-like profile of the scattering resonance $\sigma_{n}(\delta)$ [16, 20].

2.5 Decay rates

Setting $\Psi_{i}=\langle i|{\bf\lambda_{n}^{(R)}}\rangle$ and choosing a complex frequency $\omega=\omega_{0}+\delta_{n}+i\gamma_{n}$ [i.e. $\delta=\delta_{n}+i\gamma_{n}$ in equation (9), while neglecting the frequency dependence of $k$ and $G$, as mentioned above], we obtain a solution of the multiple scattering equation (10) without incoming wave, i.e. for $\Psi_{0}({\bf r}_{i})=0$. Since an imaginary part of $\omega$ leads to an exponential decay in time, see equation (2), the resonance width $\gamma_{n}$ yields the decay rate of the corresponding eigenstate $|{\bf\lambda_{n}^{(R)}}\rangle$. Of special interest are states which exhibit very small decay rates: if such a state is prepared at time $t=0$, it will be stored in the system for a very long time. The physical reason for the long lifetime is that the spherical waves emitted by the individual scatterers interfere destructively with each other.

Such long-lived, subradiant states arise naturally if two scatterers are placed very close to each other [21]. From a practical point of view, however, it is not possible to choose the distance arbitrarily small, since every physical scatterer has a finite size. We will therefore show below that strongly subradiant states can also be constructed with scatterers that are not very close together.

Choosing $\delta=\delta_{n}$ in (19), where $\gamma_{n}$ appears in the denominator, one could expect that states with very small decay rate $\gamma_{n}$ also lead to a very large scattering cross section. As it turns out, however, the enumerator $\langle{\bf\Psi_{0}}|{\bf\lambda_{n}^{(R)}}\rangle\langle{\bf\lambda_{n}^{(L)}}|{\bf\Psi_{0}}\rangle$, which describes the overlap of the resonance state with the incoming wave, also tends to zero if $\gamma_{n}\to 0$ (see A). This can be understood in terms of a reciprocity argument: states which do not emit outgoing waves are effectively isolated from the rest of the system and cannot be excited by any incoming wave. Maximizing the scattering cross section therefore amounts to a compromise between minimizing the decay rate and maximizing the overlap with the incoming wave. As shown in A, this compromise amounts to concentrating the angular profile of the scattered wave into the forward and the backward direction, whereas keeping scattering into all other directions as small as possible.

3 Optimized configurations

In this section, we present our numerical results for the maximization of the scattering cross section, as well as for the minimization of the smallest decay rate among the $N$ eigenstates of the Green matrix. In both cases, we look for optimal configurations by varying the positions of the $N$ scatterers.

3.1 Numerical algorithm

Our algorithm for the optimization of the positions of the scatterers in order to maximize the scattering cross section (18) is constructed as follows [23]: The scatterers are randomly distributed in a sphere with a radius $R$. At least 10 000 random configurations are chosen, and the scattering cross section is evaluated for each of them. For the best of these random positions, a downhill simplex-algorithm [22] searches for the local maximum. These steps are repeated at least 1000 times for each value of the radius $R$, and, finally, the whole procedure is repeated again for different radii ($kR=1,2,\dots,12$). At the end, we save the ten best configurations found in this way with the associated optimized scattering cross sections. For simplicity, we restrict the optimization to the case $\delta=0$, where the frequency of the incoming wave exactly coincides with the resonance frequency of the individual scatterers. As we have checked, an additional optimization of the incident frequency leads only to a small improvement of the scattering cross section (of the order of a few percent for odd numbers of scatterers, and even less for even numbers).

The program for the minimization of the decay rate [24] is designed in a similar manner, with the exception that the minimization is performed under the constraint of a minimal exclusion radius $r_{\rm excl}$ around each scatterer to ensure that the scatterers are not at the same place, see the discussion in chapter 2.5.

3.2 Maximization of the scattering cross section

For $N\leq 6$ scatterers, the optimal configurations found by the above algorithm exhibit the property that all scatterers are placed on a line parallel to the direction of the incoming plane wave. For $N>6$, it becomes increasingly difficult and time-consuming to find the optimal configuration in the $3N$-dimensional space of position variables. In this case, better results are obtained if the optimization is restricted to one-dimensional configurations in accordance with the finding for $N\leq 6$. For even $N$, it turns out to be advantageous to further restrict the search to symmetric configurations (see below). For large numbers of scatterers ($N>18$ for even $N$, and $N>10$ for odd $N$), we obtain the largest scattering cross sections by starting with optimized configuration for $N-2$, placing two additional scatterers at both ends of the chain, and applying the downhill simplex algorithm to this configuration.

Since an analytical upper bound of the scattering cross section is not available (see also A), there is no strict proof that our numerically optimized configurations are truly the optimal ones, although we suspect that this is the case. Apart from numerical evidence, another indication that our configurations are indeed the optimal ones is the fact that they behave in a regular way as a function of $N$ (see figures 2 and 3).

The corresponding optimized scattering cross sections are shown in figure 1 as a function of $N$. For each value of $N\geq 4$, we find two different optimized configurations (a ‘narrow’ one and a ‘wide’ one, see figures 2 and 3) which exhibit almost the same scattering cross section. For $N=8$ and $N=10$ the scattering cross section of the wide configuration (circles) is slightly larger than the scattering cross section of the narrow configuration (crosses), and vice versa for all other values of $N$ up to $N=21$. The fit function $\sigma(N)=aN^{2}+bN$ (solid line) indicates that the optimized scattering cross sections grow quadratically as a function of $N$ for large $N$. They hence exceed the linearly increasing result obtained for independent scatterers that are placed far away from each other (dotted line). This is remarkable since, as we have verified, the scattering cross sections of non-optimized, random configurations of scatterers (with distances between the scatterers of the order of the wavelength or smaller) are typically much smaller than $N\sigma^{(1)}_{\rm max}$. For one-dimensional, regular arrays (with lattice spacing $a$ between neighbouring scatterers), we have checked that the scattering cross section $\sigma(N)$ grows linearly with $N$ for large $N$ (and only slightly exceeds the value $N\sigma^{(1)}_{\rm max}$ for independent scatterers), even if the lattice spacing $a$ is optimized for each value of $N$. Therefore, only by choosing the positions of all scatterers in an optimal way, it is possible to induce cooperative effects of multiple scattering that increase the scattering cross section significantly beyond the value of $N$ independent scatterers.

The corresponding optimal positions are shown in figures 2 and 3 for the narrow and wide configurations, respectively. As already mentioned above, all scatterers are arranged on a line parallel to the direction ${\bf k}_{\rm in}$ of the incoming wave (here chosen as the $z$-direction). For even numbers of scatterers, we find that the optimized positions (in both configurations) are mirror-symmetric with respect to a certain symmetry point, which we choose as $z=0$ (since the cross section is invariant under constant translations of all scatterers). For odd numbers of scatterers, the configurations are similar to the symmetric ones for $N-1$ scatterers, with one additional scatterer placed at one side of the chain. (We choose the right-hand side in figures 2 and 3, but, due to the symmetry of the scattering cross section with respect to inverting the direction of the incoming wave, ${\bf k}_{\rm in}\to-{\bf k}_{\rm in}$, see A, the same scattering cross section is obtained when placing the additional scatterer at the left-hand side. Moreover, the point $z=0$ is choosen such that $z_{(N-1)/2}=-z_{(N+1)/2}$, see also table 1.)

As shown in A, the maximization of the scattering cross section amounts to concentrating the angular profile of the scattered wave into the forward and backward directions (${\bf k}_{\rm out}=\pm{\bf k}_{\rm in}$), while keeping the emission into other directions as small as possible. In figure 4, we therefore plot the angular profile $\sigma(\cos\theta)$ of the scattering cross section, defined by equation (14) without integral over $\Omega$, i.e. $\sigma(\cos\theta)=2\pi\left|\sum_{i}e^{ikz_{i}\cos\theta}f_{i}\Psi_{i}/\Psi_{0}\right|^{2}$ (where we used the fact that the optimized configurations are one-dimensional, and therefore the differential cross section is independent of the azimuthal angle $\phi$). Indeed, we see two pronounced peaks in the forward and backward directions, $\cos\theta=+1$ and $\cos\theta=-1$, respectively. According to this result, it is not surprising that the optimal configurations are one-dimensional (with scatterers placed on a line parallel to ${\bf k}_{\rm in}$), since, for such a configuration, it is intuitively plausible that a large asymmetry between the parallel and perpendicular directions can be achieved.

3.3 Scattering resonances of optimized configurations

In this section, we will further analyze the properties of the optimized configurations in terms of their scattering resonances. As an example, figure 5 shows the scattering cross sections $\sigma_{n}(\delta)$ associated to the individual resonances, as well as the total scattering cross section $\sigma(\delta)$, see Eq. (22), for the optimal wide configuration of $N=8$ scatterers as a function of the detuning $\delta$ of the incoming wave’s frequency from the scatterers’ resonance frequency. Since the configuration is optimized for the case $\delta=0$, it is not surprising that the total scattering cross section exhibits a sharp maximum around $\delta=0$. This maximum arises from a single, narrow resonance that gives the main contribution to the total scattering cross section at $\delta=0$. As we have checked, the same is true also for the optimal narrow configuration (where the scattering cross sections are similar to the ones shown in Fig. 5), and also for other numbers $N$ of scatterers. We therefore conclude that the maximization of the scattering cross section essentially relies on the formation of one single narrow scattering resonance, rather than on finetuning the interplay between several resonances.

3.4 Minimization of the smallest decay rate

In the previous section, we have investigated the total scattering cross section resulting from the sum of all $N$ scattering resonances, which are associated to the $N$ eigenvalues of the Green matrix, see chapter 2.4. In this section, we are interested in minimizing the decay rate resulting from the smallest imaginary part among these $N$ eigenvalues, see equation (21). The result of our numerical optimization is shown in figure 6, for different values of the exclusion radius $r_{\rm excl}$ defining the smallest possible distance between two scatterers. We see that the minimal decay rate decreases exponentially with $N$ for all three values of $r_{\rm excl}$, all of which fulfill the condition $r_{\rm excl}<\pi k^{-1}$ (see also the discussion at the end of this section). This decrease is faster for smaller values of $r_{\rm excl}$. For the smallest value $r_{\rm excl}=0.5~{}k^{-1}$ of the exclusion radius (crosses), the decay rate is suppressed by 12 orders of magnitude using only $N=10$ scatterers. In contrast, for equally spaced scatterers with distance smaller than $\pi k^{-1}$ (e.g., $r_{\rm excl}=0.5~{}k^{-1}$, see the solid line in Fig. 6), the minimal decay rate decreases only algebraically as a function of $N$ (with exponent $-3$ [24, 25]).

The corresponding optimized configurations, see figure 7, display similar properties as those that maximize the scattering cross section: The configurations that achieve the lowest decay rates are always one-dimensional configurations with an increasing spacing between the scatterers from the center towards the outer edge. While the scatterers are arranged symmetrically for even $N$, the optimal configuration for odd numbers of scatterers is realised by placing an additional scatterer well outside the exclusion radius of the others. As a result, we obtain in figure 6 (crosses) again a pattern originating from the change between even and odd numbers of scatterers, which is more pronounced than for the scattering cross section, compare with figure 1. In spite of the extremely long lifetimes of the optimized scattering resonances, the associated scattering cross sections are rather small, see, e.g., the values indicated in table 1 in B for $N=8$ and $N=9$. As already discussed at the end of section 2.5, this is due to the fact that it is difficult to excite the long-lived resonance by an incoming plane wave.

Finally, let us shortly comment on the case where the distances between the atoms exceed the value $\pi k^{-1}$. This case is not included in Fig. 6, since, as we have found, the suppression of the minimal decay rate due to collective effects then becomes considerably less efficient: both for equally spaced scatterers [24] and for optimized configurations, our numerical investigations indicate that, in this case, the minimal decay rate is always bounded between $\gamma/2$ and $\gamma$ for all values of $N$, and, in particular, does not decrease to zero in the limit $N\to\infty$.

4 Stability with respect to small deviations from the optimized positions

In chapters 3.2 and 3.4, we have seen that the maximal scattering cross section grows quadratically and the minimal decay rate decreases exponentially with the number $N$ of scatterers – provided that, in each case, the positions of all scatterers are chosen in an optimal way. From a practical point of view, it is important to know how accurately the optimized configurations must be realized, or, in other words, how stable the results are if the positions of the scatterers are allowed to fluctuate around their optimized values.

To answer this question, we start from the corresponding optimized configurations displayed in figures 2, 3 and 7, and allow the position of each scatterer to deviate from its optimal position by at most $\delta r=0.005k^{-1}$ (filled symbols) and $\delta r=0.05k^{-1}$ (open symbols). More precisely, each position is chosen randomly according to a uniform distribution inside a three-dimensional sphere with radius $\delta r$ surrounding the optimal position. Then, we determine the average scattering cross section or the average smallest decay rate, respectively, where the average is taken over $10^{5}$ random configurations.

In figure 8, we see that, with increasing number $N$ of scatterers, the optimized configurations are increasingly sensitive against small random fluctuations. This applies especially to the narrow configurations (circles), which are less robust than the wide ones (squares). Furthermore, it is evident that, for larger values of $N$, the configurations which have been optimized under the assumption that the positions of the scatterers can be precisely controlled, are no longer optimal when taking into account small fluctuations: after reaching a maximum, the average scattering cross section decreases as a function of $N$ and even drops below the value reached by independent scatterers (dotted line).

Fig. 9 shows the corresponding results for the smallest decay rate, for fluctuations around the configurations depicted in figure 7 (with exclusion radius $r_{\rm excl}=0.5k^{-1}$). The average smallest decay rate saturates as a function of $N$ at a value which is determined by the size of the fluctuations. Remarkably, even in the presence of fluctuations with size $\delta r=0.05k^{-1}$ ($\delta r=0.005k^{-1}$), it is possible to suppress the decay rate by more than three (five) orders of magnitude below the decay rate $\gamma$ of a single scatterer.

5 Conclusion & Outlook

We investigated cooperative effects of multiple scattering of scalar waves by finitely many resonant point scatterers in 3D. In particular, we found the maximal scattering cross section as well as the minimal decay rate by numerically optimizing the positions of the individual scatterers. In both cases, the optimum is realized by configurations where the scatterers are positioned on a line parallel to the direction of the incoming light. Exploiting multiple scattering, we consequently could achieve a quadratic increase of the scattering cross section as well as an exponential decrease of the minimal decay rate (in the presence of an exclusion radius $r_{\rm excl}<\pi k^{-1}$ limiting the smallest distance between two atoms) with increasing number of scatterers. We showed that the existence of a narrow scattering resonance is a necessary, but not sufficient condition for achieving a large scattering cross section. In particular, the maximization of the total scattering cross section goes along with a concentration of its angular profile into the backward and the forward direction, whereas scattering into all other directions is kept as small as possible. Allowing for experimental imperfections in the precise realization of the optimized configurations, we also investigated their stability with respect to small fluctuations of the scatterers’ positions.

Since many experiments on multiple scattering are performed with light waves, an obvious extension of our work will be to repeat our analysis for the case of vectorial waves. We expect that the optimal configurations for vectorial scatterers display similar properties as for scalar waves and, in particular, that the scattering cross section (the smallest decay rate) can be significantly enhanced (reduced) by choosing optimized configurations as compared to regular arrangements with equal spacing between the scatterers. With a view at the experiments on Bragg reflection of light by atoms mentioned in the introduction [2, 3], such optimized configurations could be used to achieve even stronger reflection with even smaller numbers of atoms, provided that the positions of the atoms can be controlled with the necessary precision.

Another potential application of our work is the realization of strongly subradiant states, which, due to their long lifetimes, could be good candidates for new types of quantum memories. As we have shown, the decay rate can be suppressed by several orders of magnitude by choosing suitable configurations of point-like emitters. However, we also pointed out that it is difficult to excite these long-lived states by an incident wave. Possible remedies are to use a non-optical way of excitation, or to shift the positions of the emitters after they have been excited. Finally, it will be an interesting, but also very challenging task, to extend our studies to the quantum mechanical regime of multiple excitations.

Appendix A An attempt to derive an upper bound of the scattering cross section

We first calculate the decay rate of a scattering resonance $|{\bf\lambda_{n}^{(R)}}\rangle$, which is defined as right-eigenvector of $G$, i.e. $G|{\bf\lambda_{n}^{(R)}}\rangle=\lambda_{n}|{\bf\lambda_{n}^{(R)}}\rangle$ Writing $|{\bf\lambda_{n}^{(R)}}\rangle=\sum_{i}c^{(n)}_{i}|i\rangle$ with $\sum_{i}|c^{(n)}_{i}|^{2}=1$ (i.e. $\langle{\bf\lambda^{(R)}_{n}}|{\bf\lambda^{(R)}_{n}}\rangle=1$), we obtain, see Eq. (21):

[TABLE]

Here, $\int{\rm d}\Omega$ denotes an integral over the angular variables of ${\bf k}_{\Omega}$ (with $|{\bf k}_{\Omega}|=k$). We used:

[TABLE]

to arrive at the third line of (23). We define

[TABLE]

We note that $g_{n}({\bf k})$ (with $|{\bf k}|=k$) is proportional to the wave amplitude emitted from the $n$-th resonance into the direction ${\bf k}$, see also (13). According to equation (23), the decay rate $\gamma_{n}$ is determined by the total flux emitted into all directions.

In order to derive an upper bound for $\sigma_{n}(\delta)$, see equation (19), we replace the imaginary part by the absolute value. Then, the maximum achieved at $\delta=\delta_{n}$ reads:

[TABLE]

Since the matrix $G$ is complex symmetric (i.e. $G_{ij}=G_{ji}$), we can choose the left eigenvector as the complex conjugate of the right eigenvector, i.e. $\langle{\bf\lambda_{n}^{(L)}}|=\sum_{i}c^{(n)}_{i}\langle i|$. Thereby, we obtain:

[TABLE]

This shows that the scattering cross section is invariant with respect to changing the direction of the incoming wave (${\bf k}_{\rm in}\to-{\bf k}_{\rm in}$). The normalization condition $\langle{\bf\lambda_{n}^{(L)}}|{\bf\lambda_{n}^{(R)}}\rangle=1$ reads $\sum_{i}(c^{(n)}_{i})^{2}=1$, which differs from the above condition $\sum_{i}|c^{(n)}_{i}|^{2}=1$ for $\langle{\bf\lambda_{n}^{(R)}}|{\bf\lambda_{n}^{(R)}}\rangle=1$. Taking this into account, we finally obtain:

[TABLE]

The upper bound is saturated if and only if the expression inside the absolute value is real and positive, i.e. if $g_{n}({\bf k}_{\rm in})g_{n}(-{\bf k}_{\rm in})/\sum_{i}(c^{(n)}_{i})^{2}>0$.

This result clarifies the qualitative discussion in chapter 2.5. In order to maximize the scattering cross section, one has to achieve large emissions $|g_{n}({\bf k}_{\rm in})|$ and $|g_{n}(-{\bf k}_{\rm in})|$ into the forward and backward directions $\pm{\bf k}_{\rm in}$ while keeping the emission into all other directions as small as possible. At the same time, the amplitudes $c_{i}$ of the resonant state should be optimized such that $\left|\sum_{i}(c^{(n)}_{i})^{2}\right|$ is minimized under the constraint $\sum_{i}|c^{(n)}_{i}|^{2}=1$.

The latter expression is, in principle, not bounded from below. This causes the main difficulty in deriving an upper bound of the scattering cross section. Indeed, it may even happen that $\left|\sum_{i}(c^{(n)}_{i})^{2}\right|=0$ for certain configurations. This is the case if the Green matrix $G$ is not diagonalisable. An example of such a case for $N=3$ is obtained if the distances between the three scatterers are chosen as $kr_{12}=kr_{13}=3\pi(4+\sqrt{2})/7$ and $kr_{23}=3\pi(1+2\sqrt{2})/14$ (which corresponds to a two-dimensional configuration). In this case, contributions from invidual scattering resonances give a diverging result, although, as we have checked, the total scattering cross section defined by (18) remains finite (and smaller than the scattering cross section of the optimized configuration).

Appendix B Optimized configurations for $N=8$ and $N=9$

In Table 1, we give the positions of all scatterers in the optimized configurations for $N=8$ and $N=9$, together with their total cross sections and minimal decay rates.