\section{Introduction.} 
One of the most influential theoretical results for solar-cell design has been the ray-optical Yablonovitch limit \cite{yablonovitch1982statistical,green2002lambertian,yu2010fundamental,yu2011angular,sheng2011optimization,wang2012absorption,callahan2012solar,ganapati2014light}, which provides a bound to how much surface texturing can enhance the performance of an absorbing film averaged over a broad bandwidth and angular range. In this paper, we obtain approximate broad-band/angle absorption limits for a case in which the traditional Yablonovitch result is not useful: dilute arrays of \textquotedblleft metaparticles\textquotedblright (synthetic absorbers/scatterers). Known limits bound the absorption at every wavelength \cite{wolgamot2012interaction,buddhiraju2017theory}, but they tend to be loose when considering large bandwidths since coherent effects average out \cite{yu2010fundamentalopt, sheng2011optimization}. Here, we find limits on the absorption for arrays of particles that can be described by the radiative-transfer equation (RTE) \cite{ishimaru1978wave, tsang2000scattering}. In particular, we show that an isotropic diffusive regime is \emph{optimal} for maximizing absorption. This allows us both to obtain analytical upper bounds (Eqs.~\ref{eq6}, \ref{eq1}) and identify the ideal operating regime of absorbing metaparticle arrays.

\begin{figure}
\includegraphics[width=0.5\textwidth]{Figure1_new.png}
\caption{Upper left: We bound absorption for very general arrays of ``particles'', including arrays of buoys that extract energy from ocean waves. Upper right: Ocean surface displacement $\eta$ for a cylindrical buoy array \cite{grgur2016optimal} where $A$ is the amplitude of waves incident from left (arrow). Lower: Sketch of RTE system.}
\end{figure}

In optics contexts, scattering particles can be used to enhance absorption in thin-film or dye-sensitized solar cells \cite{nagel2010enhanced,wang2010enhancing,rothenberger1999contribution,galvez2014dye}. Most previous work focused on numerical optimization using the full-wave equations \cite{nagel2010enhanced,wang2010enhancing} or, in the case of dye-sensitized solar cells, RTE for random arrays \cite{rothenberger1999contribution,galvez2014dye}. In \citeasnoun{mupparapu2015path}, approximate analytical estimations of absorption enhancement were given in cases of optically-thin/thick layers under assumptions of weak absorption, normal incidence and isotropic differential cross section. In this work, we were actually motivated by arrays of buoys designed to extract energy from \emph{ocean} waves \cite{falnes2007review,tollefson2014blue, stratigaki2014wave, penesis2016performance} depicted in Fig.~1. Previous numerical-optimization work \cite{cruz2009wave, child2010optimal,  grgur2016optimal, tokic2019hydrodynamics}, in particular a recent extensive computational study on large arrays \cite{grgur2016optimal, tokic2019hydrodynamics}, showed promising results through the design of buoy positions. The question we are trying to answer in this work is more general:  given the absorbing/scattering properties of an individual metaparticle, is there a limit on the total enhancement and how can it be reached? The Yablonovitch limit cannot be applied to all metaparticle arrays since it requires an effective-medium approximation, which is only accurate for either dilute weakly interacting dipolar particles \cite{choy2015effective} or for strongly interacting particles with sufficiently subwavelength separation \cite{smith2002determination}, neither of which is true of the ocean-power problem. Moreover, the Yablonovitch limit is independent of the precise nature of the scattering texture, whereas in our case the whole point is to extrapolate the array properties from the individual-scatterer properties. 
\\

In this paper, we define the \emph{interaction factor} $q(\theta)$ \cite{falnes1980radiation, evans1981power} as the ratio of the power extracted by the array to that of the equivalent number of isolated particles for a given incident angle $\theta$. We first point out that previously known limits in both solar cells and ocean buoys arise from reciprocity constraints on the full-wave equations (Section II). The use of reciprocity in the radiative-transfer equation leads to a general limit (Eq.~\ref{eq6}), valid for any geometrical configuration in RTE regime, that is reached through an isotropic distribution of intensity in the ideal case of small absorption (Section III). This optimal solution justifies the use of a \emph{corrected} radiative-diffusion model (Eq.~\ref{eq1}) that predicts the frequency-averaged performance of random arrays, but also the angle/frequency-averaged performance of the optimized \emph{periodic} array with better than 5\% accuracy. This corrected model can be used to estimate the upper bound on $q$ (which is proportional to the spatially-averaged intensity in RTE framework) even in regimes where the standard diffusion model is not expected to be accurate. This result allows us to quickly evaluate the performance benefits of different metaparticle designs and array configurations, and we show that substantial improvement is possible if the scattering cross-section is increased (relative to the absorption cross-section) and/or if partially reflecting strips are placed on either side of the array (Section IV). More specifically, we show that the use of bending membranes on the water's surface around the buoys significantly increases the interaction factor. We finally use the corrected radiative-diffusion model to find optimal parameters that maximize $q$. 

\section{Reciprocity} The original intuition behind the ray-optical Yablonovitch limit is that the optimal enhancement is achieved through an \emph{isotropic} distribution of light inside the device \cite{yablonovitch1982statistical,green2002lambertian}. This can be thought of as a reciprocity condition. Reciprocity \cite{tsang2000scattering} implies that rays at a given position cannot emerge in the same direction from two different paths. In consequence, if a given point in the absorber is to be reached from as many ray bounces as possible, the rays must be entering/exiting that point from all angles. More formally, we show in \SII{A} that reciprocity can be applied to the \emph{full} Maxwell's equations in order to relate the enhancement to the density of states (accomplished in another way by \citeasnoun{buddhiraju2017theory}), leading to:
\begin{equation} \langle q \rangle = \int_{4\pi}q(\theta)f(\theta)d\Omega \leq \frac{4\pi}{n} \frac{\langle \rho_d \rangle}{\rho_v} \max_{\theta}f   \end{equation}
where $\langle q \rangle$ refers here to the absorption enhancement compared to the single pass, averaged over both polarizations and over a directional spectrum $f(\theta)$ with normalized flux ($\int_{4\pi}|\cos \theta|f(\theta)d\Omega = 1$), $\langle \rho_d \rangle$ is the average density of states in the device, $\rho_v$ the free space density of states and $n$ the index of the absorbing medium. The previous equation becomes an equality in the case of isotropic incidence and small absorption. Yablonovitch limit can then be recovered in bulk media ($\rho_d = n^3 \rho_v$) for an incident field confined to a cone of aperture $2\theta_i$ ($f = \frac{1}{\pi \sin \theta_i^2}\delta(\theta<\theta_i)$): $\langle q \rangle \leq \frac{4n^2}{\sin \theta_i^2}$.

A similar procedure can be followed in the ocean-buoy problem. By applying the appropriate reciprocity relation derived from the wave equation, the Haskind--Hanaoka formula \cite{chiang2005theory}, to the absorption of an \emph{optimal} array of buoys \cite{falnes1980radiation}, one can bound the interaction factor $\langle q \rangle$ for a given directional spectrum $f(\theta)$ ($\int_{2\pi}f(\theta)d\theta$=1) by \citeSII{wolgamot2012interaction}{B}:
\begin{equation} \langle q\rangle = \int_{2\pi}q(\theta)f(\theta)d\theta \leq \frac{M}{k\sigma_a }2\pi \max_\theta f \end{equation}
where $k$ is the wavenumber, $\sigma_a$ the single-buoy absorption and $M$ the number of degrees of freedom for the buoy motion (1--6, e.g. 1 for only heave motion). This result implies that for \emph{isotropic incidence}, we have $\langle q \rangle \leq 1$ at the resonance frequency (the frequency at which the single buoy reaches its maximum absorption $M/k$), while it can in principle be larger at other frequencies. Although this sets a general limit valid at any frequency for any structure, we show in the following that it is not tight when considering the frequency-averaged performance. 

\section{RTE limits} 

We consider a two-dimensional array of scattering/absorbing particles distributed inside a \emph{region} $S$ bounded with a curve $C$ (Fig.~1). 

In the case of dilute and non-structured arrays, coherent scattering effects average out. This allows one to use the radiative-transfer equation (RTE) that only involves specific intensity $I(\vec{r},\theta)$, and that is applicable to ensemble averages of random arrays at a single frequency \cite{ishimaru1978wave, tsang2000scattering}: \begin{equation}\label{eqrte} \vec{e}_\theta \cdot \nabla_rI = -\rho\sigma_eI + \rho\sigma_s \int d\theta'p(\theta,\theta')I+\epsilon \end{equation} where $\sigma_s$, $\sigma_a$ and $\sigma_e$ denote respectively the scattering, absorbtion and extinction cross sections of the individual particles ($\sigma_e=\sigma_s+\sigma_a$), $p$ the normalized differential cross section, $\rho$ the particles' density, $\vec{e}_\theta$ the unit vector with direction $\theta$ and $\epsilon$ internal sources.
  
We conjecture that a similar averaging of coherent effects arises from averaging over frequency and/or angle, and below we demonstrate numerically that this allows RTE to make accurate predictions even for a small number of random samples or for optimized periodic arrays. This is similar to optical light trapping where Yablonovitch model can predict the frequency/angle-average performance of textured solar cells even though it cannot reproduce the exact spectral or angular response \cite{yu2010fundamentalopt, sheng2011optimization}.


\subsection{General limit} 

Similarly to our previous discussion of reciprocity-based limits from the wave equation, we now use reciprocity constraints on RTE to obtain general limits on the interaction factor $q$. 

One can first define a surface Green's function $G_s(\vec{r},\theta;\vec{r'},\theta')$ \cite{case1967linear} giving $I(\vec{r},\theta)$ for an incident field $I_i(\vec{r_i},\theta_i)=\delta (\theta_i-\theta') \delta (\vec{r_i}-\vec{r'})$ and no internal sources $\epsilon = 0$. Similarly, a volume Green's function $G_p(\vec{r},\theta;\vec{r'},\theta')$ can be defined as the intensity $I(\vec{r},\theta)$ obtained with no incident field $I_i=0$ and a point source $\epsilon(\vec{r_i},\theta_i) = \delta (\theta_i-\theta') \delta (\vec{r_i}-\vec{r'})$.

We recall that the flux density $\vec{F}$ is defined as $\int_{2\pi}I\vec{e}_\theta d\theta$. Conservation of energy \cite{ishimaru1978wave} then leads to $\int_C \vec{F}\cdot \vec{n}_{out}d\vec{r}=P_e-P_a$ where $P_e$ and $P_a$ are the generated and absorbed power respectively. For a unit source, we have $P_e=\int_S\epsilon(\vec{r},\theta) d\vec{r}d\theta=1$ so that: \begin{equation}\label{p-bound} \int_C\int_{\vec{e_{\theta}}\cdot \vec{n_{out}>0}}G_p(\vec{r},\theta;\vec{r'},\theta')(\vec{e}_\theta\cdot \vec{n}_{out})d\vec{r}d\theta = 1-P_a \end{equation}

To bound this last expression, we need a \textit{lower}-bound for $P_a$. By noting that the intensity at any point is larger than the \emph{single pass} value (obtained after extinction without multiple scattering), we have: \begin{equation} \label{abs-bound} \begin{split} P_a &= \rho \sigma_a \int_S\int_{2\pi}G_p(\vec{r},\theta;\vec{r'},\theta')d\vec{r}d\theta \\& \geq \rho \sigma_a\int_S \frac{e^{-\rho \sigma_e|\vec{r}-\vec{r'}|}}{|\vec{r}-\vec{r'}|}\delta[angle(\vec{r}-\vec{r'})-\theta']d\vec{r}\\& = \frac{\sigma_a}{\sigma_e} H_{\rho\sigma_e}(\vec{r'},\theta') \end{split} \end{equation} 
where $H_{\rho \sigma_e}(\vec{r'},\theta')$ defined in the previous equation can be interpreted as the power absorbed by a medium \emph{without} scattering and with an absorption coefficient $\rho\sigma_e$ in the presence of a unit source at the point $\vec{r'}$ emitting in direction $\theta'$.

Finally, we relate $G_s$ to $G_p$ through reciprocity using $G_p(\vec{r},\theta;\vec{r'},\theta')|\vec{e}_\theta\cdot \vec{n}_{out}| = G_s(\vec{r'},\pi-\theta';\vec{r},\pi-\theta)$ \cite{case1967linear}. We conclude from Eq. (\ref{p-bound}) and Eq. (\ref{abs-bound}) after a simple change of variable that: \begin{equation} \int_C\int_{\vec{e_{\theta}}\cdot \vec{n_{out}<0}}G_s(\vec{r'},\theta';\vec{r},\theta)d\vec{r}d\theta \leq 1-\frac{\sigma_a}{\sigma_e} H_{\rho\sigma_e}(\vec{r'},\pi-\theta') \end{equation} 
with equality \emph{always} realized in the absence of absorption.

\begin{figure}
\includegraphics[width=0.5\textwidth]{Figure2_new.pdf}
\caption{Upper: Frequency-averaged interaction factor $q_s$ vs incident angle $\theta$ for $N_x \times N_y = 3 \times 30$ arrays of buoys from exact solution \cite{grgur2016optimal} (solid lines), compared to standard-diffusion (black dashed lines), corrected-diffusion (red dashed lines) and radiative-transfer (RTE with Monte Carlo simulation, dots) models. ($q$ = array absorption / isolated-buoys absorption.)  The average buoy spacings (randomly chosen via a Gamma distribution) are $d_x/h = 1.73$, $d_y/h = 3.63$, with $h$ = ocean depth (the density is $\rho=1/d_yd_x$). Numbers in legend are $q_s$ averaged over $\theta$ for a typical ocean-wave directional spectrum $\cos^{2s}\theta$ with $s=4$ \cite{mitsuyasu1975observations}.  Inset: $q$ vs. wavelength at $\theta=0$, where shaded regions is one standard deviation from mean value (blue line) for 100 random structures. Lower: $\langle q \rangle$ for over isotropic incidence. Results compared to limit in Eq.~(\ref{eq6}).}
\end{figure}

Since the interaction factor in RTE is given by $q = \langle \int_0^{2\pi} I(\vec{r'},\theta')d\theta'\rangle_{\vec{r'}}/I_i$ where $I_i$ is the incident intensity and $\langle . \rangle_\vec{r'}$ is the average over $\vec{r'}$ in S, we can therefore bound the interaction factor $q$ for a given directional spectrum $f(\theta)$ [fraction of power incident from angle $\theta$]:
\begin{equation} \label{eq6} \begin{split} \langle q \rangle & = \int_{C}\int_{2\pi}\int_{\vec{e_{\theta}}\cdot \vec{n_{out}<0}} f(\theta) \langle G_s(\vec{r'},\theta';\vec{r},\theta)\rangle_{\vec{r'}} d\vec{r}d\theta' d\theta \\& \leq 2\pi \left[ 1- \frac{\sigma_a}{\sigma_e} h(\rho\sigma_e)\right] \max_{\theta} f \end{split} \end{equation} 
where $h(\rho\sigma_e) = \langle H_{\rho\sigma_e}(\vec{r'},\theta')\rangle_{\vec{r'},\theta'}\geq 0$. In the case of a ``slab" of thickness $d$, we can show that [\SII{C}]: \begin{equation} h(x) = 1-\frac 2 \pi \frac{1-e_1(xd)}{xd}, e_i(x)=\int_0^{\pi/2}e^{-x\sec \alpha}\cos^i\alpha d\alpha\end{equation}

Note that the bound in Eq.~(\ref{eq6}) reaches its maximal value $2\pi \max f$ in the limit of small absorption. This maximal value, which does not assume optimal single-buoy absorption, generalizes then the previous ocean-buoy bound, giving $\langle q\rangle \leq 1$ for isotropic incidence $f = 1/2\pi$ at any wavelength in RTE regime. In addition, $\langle q\rangle = 1$ is always realized in the small absorption limit. This special case is sometimes referred to as Aronson's theorem \cite{aronson1971theorem}. 

The equality in Eq.~(\ref{eq6}) is reached for: \begin{equation}\int\limits_{C_{\{\vec{e_{\theta}}\cdot \vec{n_{out}}<0\}}}\hspace{-1.5em}G_s(\vec{r'},\theta';\vec{r},\theta)d\vec{r} = \left[1-\frac{\sigma_a}{\sigma_e}H_{\rho\sigma_e}(\vec{r}',\theta')\right]\delta(\theta-\theta_m)\end{equation} where $\theta_m = \argmax f$. This means that the interaction factor should be equal to zero for any incident angle different from $\theta_m$. In the ideal case of small absorption, the optimal $G_s$ becomes independent of $\theta'$, which corresponds to \emph{isotropic} interior intensity, similar to the Yablonovitch model. Therefore, in order to explore \emph{optimal} solutions of RTE, we solve it under the assumption of nearly isotropic intensity, which is well known to lead to a diffusion model \cite{case1967linear, ishimaru1978wave,tsang2000scattering}. We emphasize that not all RTE systems are diffusive, but our result above shows that the optimal $\langle q \rangle$ is attained in an isotropic diffusive regime.

\subsection{Radiative-diffusion model}

Unless otherwise stated, we restrict ourselves to scatterers distributed inside a slab of thickness $d$ (Fig.~1). 

In addition to RTE parameters and reflection coefficients at the boundaries ($R_i$), the radiative-diffusion solution uses an asymmetry factor ($\mu$) \citeSII{joseph1976delta}{F} of the single particle (Fig.~3). The intensity is then given by $ I = I_{ri} + I_d $: $I_{ri}$ is the reduced intensity, solution of $\cos\theta \partial_xI_{ri} = -\rho\sigma_eI_{ri}$, and $I_d$ is the diffuse intensity approximated by $U(x)+\frac 1 \pi \vec{F}(x)\cdot \vec{e}_\theta$ where $U$ verifies a diffusion equation with flux-matching boundary conditions [\SII{D}]. By defining the cross sections per unit of length as $\upsilon_{s,a,e} = \rho d\sigma_{s,a,e}$, the model predicts an interaction factor $q$ of:  \begin{equation}\label{eq1} q (\theta)= q_0(\theta)\left(\eta \left[ D\frac{\xi(\upsilon_d)}{\xi(\upsilon_e\sec\theta)}+C\right]+1\right)\end{equation}  where $\upsilon_d^2=\gamma \upsilon_a(\upsilon_e-\upsilon_s\mu)$ is the diffusion coefficient [$\gamma = 2$ (resp. $=3$) in 2D (resp. 3D)], $\xi(x)$ is the function $(1-e^{-x})/x$, $C =\gamma [\upsilon_s(\upsilon_e+\mu\upsilon_a)]/[\upsilon_d^2-(\upsilon_e\sec \theta)^2]$, $D$ is given by the boundary conditions,  $q_0(\theta)$ is the reduced factor and $\eta$ is an additional correction term that we discuss later. General formulas for $q_0(\theta)$ and $D$ are given in \SII{E}, but in the absence of reflecting walls ($R_i=0$) they simplify to $q_0(\theta) = \xi(\upsilon_e\sec\theta)$ and: \begin{equation} D = -\frac{C(1+e^{-\upsilon_e\sec\theta})+\beta \frac{ (C+\gamma p_1\cos^2\theta)}{(1-p_1)\cos\theta}(1-e^{-\upsilon_e\sec\theta})}{(1+e^{-\upsilon_d})+\beta \frac{\upsilon_d}{\upsilon_e(1-p_1)}(1-e^{-\upsilon_d})},\end{equation}where $p_1=\sigma_s\mu/\sigma_e$ and $\beta = \pi/4$ (resp. $=1$) in 2D (resp. 3D). 

Equation (\ref{eq1}) with $\eta=1$ is obtained from the standard diffusion model. However, it is also known that the diffusion solution is inaccurate for small thicknesses \cite{kim2011correcting,chen2015extending,tricoli2018optimized}. A major problem is that it does \emph{not} guarantee $\langle q\rangle=1$ for isotropic incidence and negligible absorption, even though we previously mentioned that this is the case for any solution of RTE. The reason behind this problem is that the term $I_{ri}$ is not isotropic even for an isotropic incidence. For large thicknesses, however, the contribution of the term $I_{ri}$ becomes negligible and the diffuse term $I_d$ can ensure an isotropic solution. This simply means that the higher order terms in the expression of $I_d$ cannot be neglected for small thicknesses. In order to keep the simplicity of the diffusion solution, we suppose that the effects of higher order terms can be incorporated by the introduction of a scalar term in the diffuse intensity $\eta I_d$ instead of $I_d$. $\eta$ is then defined so as ensure the condition $\langle q\rangle=1$ for isotropic incidence and zero absorption. This procedure is somewhat similar to the approach in \citeasnoun{chen2015extending} except that we use a constant factor $\eta$ since we are interested in the total $q$ and not the spatially resolved $I$. In order to define $\eta$, we study the limit of negligible absorption for which $\upsilon_d \to 0$, $C \to -2\cos^2\theta$ and $D\to \cos^2\theta(1-e^{-\upsilon_e\sec\theta})+\frac \pi 4 \cos\theta (1-e^{-\upsilon_e\sec\theta})$. After simplification, the condition $\langle q\rangle=1$ allows to define $\eta$ as: \begin{equation}\label{eq2} \eta = \frac{\frac \pi 2-\frac{1}{\upsilon_e}[1-e_{1}(\upsilon_e)]}{\beta+\frac \pi 8 \gamma-\frac{2\gamma}{3\upsilon_e}-\beta e_1(\upsilon_e)+\frac \gamma 2 e_2(\upsilon_e)+\frac{\gamma}{\upsilon_e}e_3(\upsilon_e)}.\end{equation}  We note that, as expected, $\eta \to 1$ for an absorber that is thick compared to the extinction length. From our discussion above, this corrected radiative-diffusion model can now be used to estimate the upper bound on the interaction factor even in regimes where the standard diffusion model is not expected to be accurate (optically thin or large absorption). 


\section{Ocean-buoy arrays} 

\subsection{Example} We now present a validation of the accuracy of Eq.~(\ref{eq1}) in a model of ocean-wave energy converter (WEC) consisting of a truncated cylinder in heave motion (Fig.~1). The isolated-buoy properties can be obtained analytically \cite{garrett1971wave, yeung1981added, bhatta2003scattering} and are depicted in Fig.~3: they are designed \cite{grgur2016optimal} to have an absorption resonance that matches the typical Bretschneider spectrum \cite{bretschneider1959wave} of ocean waves. We choose the array density based on an earlier optimized periodic 3-row WEC arrangement \cite{grgur2016optimal}.  For this density, we then compare the exact numerical scattering solution calculated for both random and optimized-periodic arrays (using the method from \citeasnoun{grgur2016optimal}) to both the analytical radiation-diffusion $q$ from Eq.~(\ref{eq1}), with and without the correction $\eta$, and the numerical solution of RTE model by a Monte Carlo method \cite{marchuk2013monte}.

In Fig.~2 (upper plot), our corrected model agrees to $< 2\%$ accuracy with exact solutions for random arrays at $\theta < 80^\circ$, as long as the results are frequency-averaged. The importance of frequency averaging is shown by the $q$ frequency spectrum shown in the inset for $\theta=0^\circ$. For an ensemble of random structures, this spectrum exhibits a large standard deviation (gray shaded region), due to the many resonance peaks that are typical of absorption by randomized thin films \cite{yu2010fundamental, sheng2011optimization}, but the \emph{frequency average} mostly eliminates this variance and matches our predicted $q(\theta)$. Precisely such an average over many resonances is what allows the Yablonovitch model to accurately predict the performance of textured solar cells even though it cannot reproduce the detailed spectrum \cite{yu2010fundamentalopt, sheng2011optimization}. 

At first glance, our model does \emph{not} agree in Fig.~2 with the performance of the optimized periodic array from \citeasnoun{grgur2016optimal}: the periodic array, which was optimized for waves near normal incidence, is better at $\theta$ near $0^\circ$ and worse elsewhere. However, when we \emph{also} average over $\theta$ (from a typical ocean-wave directional spectrum \cite{mitsuyasu1975observations}), the result (shown as a parenthesized number in the legend of Fig.~2) matches Eq.~(\ref{eq1}) within 5\%.  If we average over all angles assuming an \emph{isotropic} distribution of incident waves, the results match within 1\%.  Similar results have been observed for thin-film solar cells, in which an optimized structure can easily exceed the $4n^2$ Yablonovitch limit for particular incident angles, but the Yablonovitch result is recovered upon angle-averaging \cite{yu2010fundamentalopt, yu2011angular, sheng2011optimization, ganapati2014light}.

Finally, we note in Fig.~2 (lower plot) that RTE results respect indeed the bound in Eq.~(\ref{eq6}) for isotropic incidence. In particular, we confirm that random arrays achieve $\langle q \rangle = 1$ for small absorption (i.e. small wavelength in our case). The periodic array, on the other hand, doesn't satisfy this relation unless it is frequency averaged. We also mention that the limit Eq.~(\ref{eq6}) is very loose for anisotropic incidence and cannot be reached without using external reflectors as discussed in Section IV-B below.

\subsection{Larger interaction factor} 

\begin{figure}
\includegraphics[width=0.5\textwidth]{Figure3_new.pdf}
\caption{Properties of a single truncated-cylinder wave energy converter (WEC) in heave (vertical) motion, with radius $a=0.3h$ and draft $H=0.2h$ where $h$ is the ocean depth. The WEC has an isotropic response with respect to the direction of the incident field. Left: Scattering and absorption cross sections of a single buoy normalized to the cylinder diameter ($\sigma/2a$). The ocean spectral energy density (energy per horizontal surface) is chosen as Bretschneider \cite{bretschneider1959wave} with resonant frequency matching that of the body and is shown in units of $\rho g H^2T$ ($\rho$ is the water density, $g$ the acceleration of gravity, $T$ the mean wave period and $H$ the significant wave height). Right: Asymmetry factors, defined as the average of $\cos\phi$ and $\cos2\phi$ for the two-dimensional differential scattering cross section. These parameters enter into the diffusion equation as $\mu  = (\mu_1-\mu_2)/(1-\mu_2)$ and with $\sigma_s$ replaced by $\sigma_s (1-\mu_2)$ \citeSII{joseph1976delta}{F}.}
\end{figure}

Given this model, we can now explore ways to increase the interaction factor $q$. By examining the dependence of $q$ in Eq.~(\ref{eq1}) on the parameters (Fig.~4), we find that for a fixed scattering-to-absorption ratio $\sigma_s/\sigma_a$, $q$ reaches a maximum $q_{max}$ for an intermediate value of scattering per unit length $ \rho d\sigma_s$, whereas it increases monotonically with $\mu$. A maximum $\mu$ is achieved by increasing $\mu_1$ (forward scattering) and decreasing $\mu_2$ (lateral scattering). The optimal value of $\rho d\sigma_s $ and $q_{max}$ both increase with $\sigma_s/\sigma_a$; as the single particle absorbs more, the interaction factor decreases and the optimal configuration requires a larger spacing between the particles. The maximum $q$ is then achieved in the limit of small absorption ($ \rho d\sigma_a\ll 1$) and large scattering ($ \rho d\sigma_s\gg 1$) for which we obtain a perfect isotropic diffuse intensity. 

From Fig.~3, we see that we have $\sigma_a/\sigma_s \approx 1$ at the resonance of the WEC. In this case, the enhancement is expected to be smaller than 1 around the resonance and optimal structures will tend to have a large spacing $d_y$. (If the array were optimized for small wavelengths $\lambda$, where $\sigma_s \gg \sigma_a$, then a larger $q$ could be obtained at those wavelengths, but $q_s$ would be worse because the optimal spacing in this case is too small for good performance at the resonance.) Still, multiple scattering significantly improves the broadband performance of our array: our $\langle q \rangle \approx 0.99$ is larger than the $\langle q_0 \rangle \approx 0.78$ that is obtained from RTE in absence of multiple scattering (reduced factor $q_0$). The performance is still lower than the 1.65 that would be obtained for $\sigma_s \gg \sigma_a$ in the ideal isotropic regime discussed below, essentially because $\sigma_a/\sigma_s$ is too small and the structure is too thin (as for example quantified by the transport mean-free path $d/l_{tr}={\upsilon_s(1-\mu)} \approx 0.5$ for $\frac{2a}{\lambda} \gtrsim 0.3$) to practically achieve an isotropic diffuse intensity.

Alternatively, we show that $q$ can be enhanced by putting partially reflecting strips around the array. Similar to light-trapping by total internal reflection \cite{yablonovitch1982statistical,green2002lambertian}, one possibility is to use a strip of a lower-``index'' \cite{chiang2005theory} medium (compared to the array's ambient medium) on either side of the array. In the ocean-buoy problem, this can for example be achieved by either a depth change or the use of a tension/bending surface membrane which can lead to near-zero index \cite{zhang2014broadband, bobinski2015experimental}. This modifies equations (2--4) with additional reflection coefficients $R_i$, as given in \SII{E}.

In Fig.~4, we show the effect of an increase in the scattering cross section and/or the index contrast for the same array studied before. By combining both effects, a large ($> 3$) spectral interaction factor can be achieved at normal incidence. At the same time, waves incident at large angles will be reflected out, so that the interaction factor integrated over isotropic incidence is still smaller than 1. For a given directional spectrum and scattering cross section of a single buoy, the optimal interaction factor is achieved for a specific value of the index contrast as can be seen in Fig.~4 (right). 

\begin{figure}[t!]
\includegraphics[width=0.5\textwidth]{Figure4.pdf}
\caption{Upper: Dependence of $q(0^\circ)$ on parameters in absence of reflecting boundaries. In the left plot, we take $\sigma_s/\sigma_a = 5$. In the right plot, we show the optimal $\rho d\sigma_s$ and $q_{max}$ for different values of $\sigma_s/\sigma_a$ and $\mu$. Lower: Effect of a change in the index contrast and scattering cross section on the bandwidth-averaged factor $q_s$ for the same array in Fig.~2. We tune the index $n_1$ along a strip surrounding the array, with $n_0$ being the index of the array's ambient medium. We suppose that the WEC has new scattering cross section $\tilde{\sigma}_s$, but keep the same absorption cross section. Left: $q_s$ at normal incidence. Right: $q_s$ averaged over $\theta$ with a directional spectrum of $\cos^{2s}\theta$ and s = 4.}
\end{figure}

Finally, it is instructive to look at the ideal case of small absorption and large scattering, for which Eq.~(\ref{eq1}) simplifies to: \begin{equation} \label{eq7}q(\theta) =  \left[1-R_1(\theta)\right] \left(\frac{\pi}{4\alpha}+\cos \theta\right)\cos\theta \end{equation} where $R_1$ is the reflection coefficient of the front-surface and $\alpha = (1-r_1)/(1+r_2)$ with $r_i = \int_{-\pi/2}^{\pi/2}R_1(\theta)\cos^i(\theta)\text{d}\theta / \int_{-\pi/2}^{\pi/2}\cos^i(\theta)\text{d}\theta$. Equation (\ref{eq7}) still gives 1 when averaged over isotropic incidence, but the interaction factor is larger near normal incidence. Without reflectors ($R_1=0$), the maximum value of $q$ at normal incidence is $1+\frac \pi 4$, and the previous directional spectrum gives $\langle q \rangle \approx 1.65$. This maximum value of $q(0)$ does not reach the arbitrarily large enhancement allowed by Eq.~(\ref{eq6}). However, $q(0)$ can still be made sufficiently large by including a reflector designed for transmission near normal incidence and reflection elsewhere (since $\alpha \rightarrow 0$).

\subsection{Surface membrane}

We now use a specific example to demonstrate a larger interaction factor $q$ using surface membranes surrounding the WEC array. For large scale applications, such membranes could be designed to have the desired properties by connecting floating pontoons with elastic elements of appropriate stiffness.

A thin bending membrane on the water surface changes the ``refractive index'' ($\sim k/\omega$) through the following dispersion relation (e.g. \citeasnoun{fox1994oblique}):
 \begin{equation}\label{eq-dispersion} \omega^2 = gk\tanh (kh) \frac{1+C_b(kh)^4}{1+m\cdot kh\tanh(kh)} \end{equation} where $\omega$ is the frequency, $g$ the acceleration of gravity, $k$ the wavenumber, $C_b$ is a dimensionless bending coefficient, $m$ is the mass of the membrane relative to the mass of the water beneath it and $h$ is the water depth. We simply assume $m=0$ in the following. 

At a fixed $\omega$, the membrane decreases $k$ (decreases the ``index'') compared to the surrounding medium. This change of index leads to a reflection off the membrane's edges. In particular, total internal reflection traps the water waves similarly to light trapping in solar cells, which increases the interaction factor $q$. The reflection coefficient, which depends on $\omega$, $C_b$, the incident angle and the membrane's width $w$, can be computed by applying appropriate boundary conditions on either side of the membrane and using a transfer-matrix method as reviewed in \SII{G}. We note that evanescent modes need to be included because of the change in dispersion relations.

The index contrast increases with $C_b$ (increasing stiffness), which increases the range of angles undergoing total internal reflection, making a more effective mirror. Since no waves are coming from the rear of the array, the optimal membrane behind the array should be a perfect reflector ($C_b \to \infty$, limited only by the attainable practical $C_b$). 

\begin{figure}[t]
\includegraphics[width=0.5\textwidth]{Figure5_new.pdf}
\caption{Upper: $\langle q_s \rangle$ with a directional spectrum of $\cos^{2s}\theta$ and s = 4 for different values of $C_{b1}$ and $C_{b2}$ corresponding to the front and back membranes respectively. Each point is obtained after optimizing over the membranes' thicknesses. Lower: Frequency-averaged interaction factor $q_s$ vs incident angle $\theta$ for the previously studied array using additional membranes with parameters $(C_{b1},C_{b2})=(0.048,2)$ and $w_1=w_2=1.6h$.} 
\end{figure}

We can now use our corrected diffusion model to predict the upper-bound for the previously studied array as we change $C_{b}$. For each value of $C_{b1}$ and $C_{b2}$ representing the front and rear membranes, respectively, we find the optimal membrane widths that maximize the radiative-diffusion bound. The resulting optimized $\langle q_s\rangle$ values are shown are shown in Fig.~5 (upper plot). We first note that the frequency/angle-averaged interaction factor $\langle q_s\rangle$ increases significantly ($>1.8$) compared to the $\langle q_s\rangle = 1.00$ without the membranes. We also confirm that $\langle q_s\rangle$ increases with $C_{b2}$ (rear membrane) as expected. On the other hand, there is an optimal value for $C_{b1}$ depending on the directional spectrum $f(\theta)$. For a \emph{focused} incident field, only angles near normal incidence matter so that $C_{b1}$ can be increased allowing more of the waves scattered by the WECs to be trapped. On the other hand, for a broad directional spectrum, a large value of $C_{b1}$ prevents waves incident from wide angles from reaching the WECs. 

For our array, supposing for example that the maximal attainable value of $C_{b2}$ is equal to $2$, the optimal value for $C_{b1}$ is $0.048$ with optimal widths equal to $1.6h$ for both the front and rear membranes. The frequency-averaged interaction factor $q_s$ for the optimal parameters is shown in Fig.~5 (lower plot). Our predicted bound (red dashed line = corrected diffusion) is indeed larger than the actual performance of the array as modeled by RTE (orange dots). That is mainly due to the relatively small scattering cross section compared to the absorption cross section. As illustrated in the inset of Fig.~5 at small wavelengths where $\sigma_s$ is large (Fig.~3), we see that an increase in the scattering cross section leads to arrays with performance closer to the radiative-diffusion bound. 

We finally mention that in the case of using a perfect back-reflector, $\langle q_s\rangle$ can reach a value of 2.26 for $C_{b1}=0.06$ and $w_1=1.65h$.

\section{Conclusion.} 
We believe that the angle/frequency-averaged limits presented in this paper provide guidelines for future designs to achieve a large $q$ factor which may open the path for the realization of large arrays of buoys for efficient ocean energy harvesting. In particular, the use of external reflecting elements such as surface membranes seems a promising approach. The results are also applicable to other problems where multiple scattering effects are used to achieve enhancement, including scattering particles inside an absorbing layer. One can, for example, recover the standard Yablonovitch-$4n^2$ result from our approach in an appropriate limit [\SII{H}], but the real power of our result is that it allows to study the effect of single-metaparticle properties, angle of incidence and reflecting boundaries. 

\begin{acknowledgments}
This work was supported in part by the Army Research Office under Cooperative Agreement Number W911NF-18-2-0048. 
\end{acknowledgments}

