1. Introduction
Ocean waves offer an abundant and consistent source of renewable energy (Pecher & Kofoed Reference Pecher and Kofoed2017). For over half a century, researchers have sought to harness this power using wave energy converters (WECs) (Cruz Reference Cruz2008), although they are yet to reach economic viability (Gallutia et al. Reference Gallutia, Fard, Soto and He2022). Heaving buoys with attached power take-off mechanisms (PTOs) are one of the most popular proposed/tested WEC concepts for commercial scale ocean-wave energy conversion (IRENA 2020). Like most WECs, heaving buoys are tuned to resonate, such that maximum absorption is achieved through optimum wave interference at the resonant frequency (Falnes Reference Falnes2005). But a naïve deployment of WECs in the inherently broadband, seasonal, ocean-wave climate (Pecher & Kofoed Reference Pecher and Kofoed2017) will not provide the broadband absorption needed to make them a competitive energy source, as they are efficient only in a narrow frequency range due to their physical characteristics (Falnes & Hals Reference Falnes and Hals2012).
WECs are deployed in arrays (or ‘farms’) to help meet energy demands and improve cost-effectiveness (Pecher & Kofoed Reference Pecher and Kofoed2017). Wave interactions between the WECs in an array can either enhance or degrade the performance of the array in comparison with the same number of WECs operating independently (Göteman et al. Reference Göteman, Giassi, Engström and Isberg2020; Gallutia et al. Reference Gallutia, Fard, Soto and He2022). Designing WEC-arrays to enhance power capture is an active research topic, with most studies based on numerical and mathematical models, of which the majority use linear potential-flow theory (Folley Reference Folley2016). Approaches include optimising WEC geometries/layouts (Babarit Reference Babarit2013; Göteman et al. Reference Göteman, Giassi, Engström and Isberg2020; Edwards & Yue Reference Edwards and Yue2022), employing control strategies via the PTO mechanism to maintain/create resonant conditions over broad frequency bands (Bacelli & Ringwood Reference Bacelli and Ringwood2013; Pecher & Kofoed Reference Pecher and Kofoed2017), or a combination of optimising the WEC layout and a control strategy (Garcia-Rosa, Bacelli & Ringwood Reference Garcia-Rosa, Bacelli and Ringwood2015; Golbaz et al. Reference Golbaz, Asadi, Amini, Mehdipour, Nasiri, Etaati, Naeeni, Neshat, Mirjalili and Gandomi2022).
Theoretical progress has been made for uniform WEC-arrays, i.e. identical and equally spaced WECs with identical PTO parameters. Without absorption, the so-called Bloch waves supported by the array separate into passbands (where the Bloch waves propagate through the array) and bandgaps (where the Bloch waves are unable to propagate). Bandgaps are related to destructive wave interference between WECs in an array, which considerably reduces power capture, whereas frequencies on passbands just below the bandgap are related to constructive wave interference over the array that enhances power capture (Garnaud & Mei Reference Garnaud and Mei2010; Tokić & Yue Reference Tokić and Yue2019). In general, bandgaps can be created by Bragg resonance, which is controlled by the spacing of the WECs in the array (Tokić & Yue Reference Tokić and Yue2019), or local resonances of individual WECs in the array. In the latter case, the bandgap structure has been shown experimentally to be robust to nonlinear wave breaking caused by the local resonances (motivated by coastal protection and without absorption (Dupont et al. Reference Dupont, Remy, Kimmoun, Molin, Guenneau and Enoch2017)).
In order to be considered broadband, a WEC-array would need to achieve high efficiency over a broad power capture interval, most likely lying within wave periods between 
$5$ and 20 s, which are feasible for ocean wave energy capture (Coe, Bacelli & Forbush Reference Coe, Bacelli and Forbush2021). The PTO parameters can be adjusted in time according to predicted sea states to achieve broadband absorption, but this approach requires accurate wave prediction capability (Fusco & Ringwood Reference Fusco and Ringwood2010; Gallutia et al. Reference Gallutia, Fard, Soto and He2022) and is sensitive to model errors (Cruz Reference Cruz2008; Folley Reference Folley2016), while increasing cost, system complexity and maintenance requirements (Garnaud & Mei Reference Garnaud and Mei2009; Pecher & Kofoed Reference Pecher and Kofoed2017). Designing WEC-arrays to be capable of broadband absorption with PTO parameters that are constant in time would improve economic viability (Pecher & Kofoed Reference Pecher and Kofoed2017; Sergiienko et al. Reference Sergiienko, Cazzolato, Ding, Hardy and Arjomandi2017).
Spatially graded arrays are a promising design strategy (Bennetts, Peter & Craster Reference Bennetts, Peter and Craster2018) to achieve broadband absorption (Porter Reference Porter2021; Wilks, Montiel & Wakes Reference Wilks, Montiel and Wakes2022). Grading (without absorption) is used for spatially controlled amplification of wave energy within an array according to frequency, as wave energy at a certain frequency is gradually slowed and accumulates as it reaches the location in the array at which the local periodicity defines the transition from a passband to a bandgap (Bennetts et al. Reference Bennetts, Peter and Craster2018). Predictions of this phenomenon from linear theory have been confirmed experimentally, in spite of nonlinear wave breaking where the waves are amplified, which merely reduces the amplification factors (Archer et al. Reference Archer, Wolgamot, Orszaghova, Bennetts, Peter and Craster2020). Without absorption (e.g. Bennetts et al. Reference Bennetts, Peter and Craster2018; Archer et al. Reference Archer, Wolgamot, Orszaghova, Bennetts, Peter and Craster2020; Xu, Ning & Chen Reference Xu, Ning and Chen2024), the graded array ultimately reflects the amplified wave energy (hence, the phenomenon is referred to as rainbow reflection), but with absorption the energy can potentially be captured (e.g. Chaplain et al. Reference Chaplain, Pajer, De Ponti and Craster2020). The resulting rainbow absorption has been demonstrated in acoustics (e.g. Romero-García et al. Reference Romero-García, Theocharis, Richoux and Pagneux2016; Jiménez et al. Reference Jiménez, Romero-García, Pagneux and Groby2017), and elasticity (e.g. Chaplain et al. Reference Chaplain, Pajer, De Ponti and Craster2020).
Of particular relevance to the present study, Wilks et al. (Reference Wilks, Montiel and Wakes2022) introduced the concept of rainbow absorption to ocean wave energy harvesting, using a two-dimensional (2-D), linear model of wave–structure interactions. Closely following the approach of Jiménez et al. (Reference Jiménez, Romero-García, Pagneux and Groby2017) in acoustics, Wilks et al. (Reference Wilks, Montiel and Wakes2022) designed a water-wave rainbow reflection structure from an array of surface-piercing, rigid, vertical barriers by grading the width between barriers and their submergence. They then added heaving pistons with attached linear damping to create rainbow absorption, and optimised the graded geometry to achieve broadband absorption of 98.2 % on an angular frequency interval of 0.8–1.6 rad s
$^{-1}$ (wave periods 
$\approx$4–8 s). The barrier submergence increases with distance along the array, reaching up to 80 % of the water depth, such that the isolated barrier allows virtually no transmission for the wave frequencies of interest (B. Wilks, personal communication). The corresponding periodic arrays support a single passband at low frequencies, followed by a bandgap that extends to infinity in frequency (although interspersed with very narrow high-frequency passbands (Wilks et al. Reference Wilks, Montiel, Bennetts and Wakes2024)). Their approach is broadly similar to positioning a WEC-array near a rigid wall or breakwater (e.g. Konispoliatis & Mavrakos Reference Konispoliatis and Mavrakos2020) to increase absorption by harnessing reflections (Falnes & Hals Reference Falnes and Hals2012).
In this study, we show that it is possible to create rainbow absorption with a generic array of heaving buoy type WECs by grading the WEC-resonances using the PTOs while maintaining uniform WEC geometries and spacing, and without the need for a surrounding structure. We achieve near-perfect absorption (99 % efficiency) over a broad, targeted frequency range of practical interest (angular frequencies 0.3–0.65 rad s
$^{-1}$ or wave periods 
$\approx$10–20 s) using a small number of WECs (typically five) with time-independent PTO parameters. The WECs do not act like rigid barriers, and destructive wave interactions within the array are necessary to generate bandgaps. Over the frequency range considered, multiple passbands and bandgaps exist, which results in interplay between local resonance and Bragg resonance (Guo et al. Reference Guo, Cao, Xiao, Shen and Wen2020) not seen in previous cognate studies.
The governing equations and solutions methods used in the study are set out in §§ 2 and 3, respectively. The model (based on 2-D linear potential-flow theory) and methods (eigenfunction matching for individual WECs and transfer matrices for WEC interactions) are standard, but are presented for completeness and to establish key notations. The primary contribution of the article is to outline an approach to design an array that achieves near-perfect broadband absorption using bandgaps, zeros of reflection for complex-valued frequencies, and efficient optimisation algorithms. The approach is introduced in stages: by introducing Bloch waves for a uniform non-absorbing array and connecting the complex resonances of the finite array with the band structure (passbands and bandgaps) of the corresponding infinite array (§ 4); grading the PTO stiffness to achieve rainbow reflection, and then tuning the PTO damping to achieve complex zeros of reflection and, hence, rainbow absorption (§ 5); and using the knowledge gained to inform optimisation algorithms that generate near-perfect broadband absorption, which is subsequently extended to multiple degrees of freedom (surge and pitch released) and applied to irregular sea states (§ 6). The performance of one of the proposed arrays is demonstrated in the time domain for incident wave packets, for which a novel decomposition into the time-domain Bloch wave modes is used to connect with the underlying theory (§ 7). The relevance of the outlined approach to an equivalent three-dimensional (3-D) model is discussed, along with other practical next steps (§ 8).
2. Preliminaries
Consider a 2-D water domain, 
$\varOmega$, in which locations are defined by a Cartesian coordinate system 
$(x,z)$, where the 
$x$-axis is chosen to coincide with the undisturbed free surface, and the 
$z$-axis is directed out of the water. The water domain is bounded below by an impermeable seabed at 
$z=-h$, and above by a free surface, 
$z=\eta (x,t)$, where 
$t$ is time. An array of 
$N$ geometrically identical WECs (WEC 1, WEC 2, …, WEC 
$N$) occupies the water domain. Each WEC involves a heaving buoy and a PTO mechanism, broadly representative of the bottom-referenced heaving buoy developed by CorPower Ocean (2024). For ease of computation, the buoys are square with side lengths 
$2L$, but should have rounded vertices in practice to prevent flow separation and vortices (Yeung & Jiang Reference Yeung and Jiang2014). The buoys are evenly spaced along the free surface with separation distance 
$d$, and the first buoy is centred at 
$x=0$ (figure 1).
Schematic of the problem involving an array of 
$N$ WECs forced by an incident wave from 
$x\to -\infty$ (blue arrow).
Water motions are modelled using linear potential-flow theory (inviscid, incompressible and irrotational fluid with small amplitude relative to wavelength 
$\lambda$), such that the velocity field is the gradient of a scalar velocity potential, 
$\varPhi (x,z,t)$ (Linton & McIver Reference Linton and McIver2001). The free surface is related to the velocity potential via the linearised dynamic condition
\begin{equation} [\partial_{t}\varPhi]_{z=0} ={-}g\eta, \end{equation}
where 
$g \approx 9.81\ \text {m}\ \text {s}^{-2}$ is gravitational acceleration and 
$\partial _{\bullet }$ denotes the partial derivative with respect to the subscript (Linton & McIver Reference Linton and McIver2001). Under the assumption of time-harmonic motions at a prescribed angular frequency 
$\omega \in \mathbb {R}>0$, the velocity potential and surface elevation can be written as
\begin{align} \varPhi(x,z,t) = \operatorname{Re} \left \{ \frac{gA}{\mathrm{i} \omega}\phi(x,z)\exp({-\mathrm{i} \omega t})\right\}\quad\text{and}\quad \eta(x,t)= \operatorname{Re}\{A\zeta(x)\exp({-\mathrm{i} \omega t})\}, \end{align}
where 
$\phi \in \mathbb {C}$ and 
$\zeta \in \mathbb {C}$ are a reduced potential and surface elevation, respectively, 
$A$ is the incident wave amplitude and 
$\mathrm {i} = \sqrt {-1}$ is the imaginary unit. The potential, 
$\phi$, satisfies Laplace's equation
\begin{equation} \nabla^2\phi = 0 \quad\text{throughout } \varOmega, \end{equation}and a no-normal flow condition at the seabed
\begin{equation} \partial_z \phi = 0 \quad \text{on } z ={-} h. \end{equation}
At the free surface, the potential is related to the surface elevation, 
$\zeta$, by dynamic and kinematic conditions, respectively,
\begin{equation} \zeta=[\phi]_{z=0}\quad\text{and}\quad \zeta \omega^2 / g = [\partial_{z}\phi]_{z=0}, \end{equation}which can be combined into the single boundary condition
\begin{equation} \partial_{z} \phi =\frac{\omega^2}{g} \phi \quad\text{at } z = 0, \end{equation}for the potential only. Radiation conditions are applied in the far-field (Linton & McIver Reference Linton and McIver2001).
The PTOs are modelled as linear spring–damper mechanisms, defined by stiffness and damping coefficients, respectively,
\begin{equation} c_{(n)}^{PTO}\quad\text{and}\quad b_{(n)}^{PTO}\quad\text{for}\ n = 1, 2, \ldots N. \end{equation}
The damping coefficients are non-negative but the stiffness coefficients can be positive or negative (Todalshaug et al. Reference Todalshaug, Ásgeirsson, Hjálmarsson, Maillet, Möller, Pires, Guérinel and Lopes2016; Kurniawan & Zhang Reference Kurniawan and Zhang2018). The complex amplitude of the vertical (heave) oscillations of WEC 
$n$ is denoted by 
$\xi _{(n)}^h\exp ({-\mathrm {i} \omega t})$ (
$n=1,2,\ldots,N$). Linearised boundary conditions on the wetted surfaces of each WEC are given by
\begin{equation} \partial_z \phi = \frac{\omega^2}{g}\xi_{(n)}^{h} \quad\text{for}\ x_L^{(n)} \leq x \leq x_R^{(n)}\quad\text{and}\quad z ={-}L \end{equation}and
\begin{equation} \partial_x \phi = 0\quad\text{for}\ x = x_L^{(n)},\ x_R^{(n)}\quad\text{and}\quad -L < z < 0, \end{equation}
for 
$n=1,2,\ldots,N$, where 
$x_L^{(n)}$ and 
$x_R^{(n)}$ denote the left- and right-hand edges of WEC 
$n$, respectively (figure 1).
For an incident plane wave with angular frequency 
$\omega$, the heave amplitude 
$\xi _{(n)}^h$ of WEC 
$n$, where 
$n = 1,2, \ldots, N$, is obtained from the equations of motion
\begin{equation} [-\omega^2(\mathcal{M}+a(\omega))-\mathrm{i}\omega (b(\omega)+b_{(n)}^{PTO}) + (c+c_{(n)}^{PTO})]\xi_{(n)}^h(\omega) = F_{(n)}(\omega), \end{equation}
given the excitation force 
$F_{(n)}(\omega )$ on WEC 
$n$, which depends on the wave amplitudes within an array (i.e. it accounts for WEC interactions). The heave excitation force arises from the hydrodynamic pressure exerted on the underside of a fixed WEC, while the added mass, 
$a(\omega )$, and radiation damping, 
$b(\omega )$, result from forced oscillations in heave, and are proportional to the acceleration and velocity of the body, respectively (Linton & McIver Reference Linton and McIver2001; Mei, Stiassnie & Yue Reference Mei, Stiassnie and Yue2005). Assuming a suitable static force component of the PTO, the mass per unit breadth of a single WEC is 
$\mathcal {M} = 4100 L^2\ \text {kg}\ \text {m}^{-3}$, and 
$c$ is the hydrostatic stiffness (Mei et al. Reference Mei, Stiassnie and Yue2005).
The proportion of incident wave energy captured by the WEC-array, 
$\alpha (\omega )$, can be determined from the scattering coefficients of the array (reflection, 
$R\in \mathbb {C}$, and transmission, 
$T\in \mathbb {C}$, for the array of heaving WECs), using the relation
\begin{equation} \alpha(\omega) = 1-|R(\omega)|^2-|T(\omega)|^2. \end{equation}
To achieve perfect absorption (
$\alpha =1$) at some 
$\omega _{0} \in \mathbb {R}^+$, the WEC-array must prevent reflection and transmission of waves at 
$\omega _{0}$, implying
\begin{equation} |R(\omega_{0})|^2=0\quad\text{and}\quad |T(\omega_{0})|^2 = 0. \end{equation}
For a single WEC in isolation (i.e. an array with 
$N=1$ WEC), the optimal PTO parameters for maximum power capture at a specified resonance frequency, 
$\omega =\omega _{0}$, can be derived as
\begin{equation} c^{PTO} = \omega_0^2(\mathcal{M} + a(\omega_0))- c \quad \text{and}\quad b^{PTO} = b(\omega_0), \end{equation}
which gives 
$\alpha (\omega _{0})=0.5$ for a rigid, axisymmetric, heaving buoy (Evans Reference Evans1981; Falnes Reference Falnes2005). This can be interpreted as setting the stiffness coefficient to tune the natural WEC resonance to the chosen frequency and the damping coefficient such that the radiation cancels the maximum possible amount of diffracted energy.
Broadband absorption over the target frequency interval
\begin{equation} \omega_{\alpha}\equiv[\omega_{lb},\omega_{ub}]=[0.3,0.65]\ \text{rad}\ \text{s}^{{-}1} \end{equation}
is defined in terms of the mean absorption over the interval, 
$\hat {\alpha }$, such that
\begin{equation} \hat{\alpha}= \frac{1}{\omega_{ub}-\omega_{lb}} \int_{\omega_{lb}}^{\omega_{ub}}\underbrace{(1-|R(\omega)|^2-|T(\omega)|^2)}_{\alpha(\omega)}\ \mathrm{d}\omega. \end{equation}
For the purposes of this study, near-perfect absorption is defined as 
$\alpha \geq 0.990$ at a single frequency and 
$\hat {\alpha } \geq 0.990$ over the target interval, which covers two thirds of usable ocean wave frequencies (Coe et al. Reference Coe, Bacelli and Forbush2021).
Parameter values broadly representative of CorPower Ocean's C4 device are applied to model the WEC-array (Alday, Raghavan & Lavidas Reference Alday, Raghavan and Lavidas2023). Each WEC has width 
$2L=10$ m (Babarit Reference Babarit2013), and the water depth is 
$h=50$ m (Liu et al. Reference Liu, Meucci, Liu, Babanin, Ierodiaconou, Xu and Young2023), as is typical of heaving buoy arrays. However, different dimensions and wavelengths could equivalently be employed. Negative stiffness coefficients are used to decrease the natural WEC-resonances in order to achieve resonance on the target interval. A negative stiffness mechanism can broaden the resonance bandwidth of the WECs, and has been successfully applied to tune the C4 devices for power capture over the target interval 
$\omega _\alpha$ (Todalshaug et al. Reference Todalshaug, Ásgeirsson, Hjálmarsson, Maillet, Möller, Pires, Guérinel and Lopes2016; Kurniawan & Zhang Reference Kurniawan and Zhang2018; Satymov et al. Reference Satymov, Bogdanov, Dadashi, Lavidas and Breyer2024).
Without loss of generality, the array is forced by an incident plane wave of amplitude 
$A=1$ m from the left of the domain. Based on the corresponding wavelengths (
$\lambda \in [142,429]$ m), the steepness of the incident waves, 
$\epsilon$, satisfies 
$0.0149 < \epsilon < 0.042$, where 
$\epsilon = kA \ll 1$ is required to operate in a linear regime (Ding & Ning Reference Ding and Ning2022).
3. Solution method
3.1. Single cell problem
Unit cell 
$n$ is defined as the subregion 
$\varOmega _{(n)}$ of length 
$W=2L+d$, containing WEC 
$n$ at its centre. The unit cell itself is partitioned into three regions, with interfaces defined by the left (
$x_L^{(n)}$) and right (
$x_R^{(n)}$) edges of WEC 
$n$ (figure 1). Without loss of generality, suppose WEC 
$n$ is centred at 
$x=0$ so that the unit cell is defined on 
$|x|\leq W/2$, with the left- and right-hand edges of the WEC given by 
$x_L^{(n)} = -L$ and 
$x_R^{(n)} =L$, respectively.
Separation of variables is used to solve for 
$\phi \equiv \phi _{(n)}$ on 
$\varOmega _{(n)}$ (using the code of Yiew et al. (Reference Yiew, Bennetts, Meylan, French and Thomas2016)), which is expressed as an eigenfunction expansion in each region (Linton & McIver Reference Linton and McIver2001). Dynamic and kinematic boundary conditions (ensuring continuity of pressure and horizontal velocity) are applied at interfaces between the regions (
$x={\pm }L$) to determine the unknown coefficients of the expansions (Linton & McIver Reference Linton and McIver2001), which leads to a system of linear equations involving the amplitudes, denoted 
$a^{\pm }_{m}$ in Region 1 and 
$b^{\pm }_{m}$ in Region 3 (
$m=0,1,\ldots$), as well as amplitudes in Region 2 that are not required for subsequent analysis (see Appendix A). Truncating each of the eigenfunction expansions at 
$m=M$ for a sufficiently large 
$M$ (
$M=25$ is used in computations presented, see Appendix B), the wave fields on either side of WEC 
$n$ are expressed as
\begin{align} & \phi_{(n)}(x,z) \nonumber\\ &\quad \approx \begin{cases} \underset{m=0}{\overset{M}{\sum}}(a_m^+\exp({\mathrm{i} k_m (x+L)})+a_m^-\exp({-\mathrm{i} k_m (x+L)})) \dfrac{\cosh(k_m(z+h))}{\cosh(k_m h)} & x <{-}L, \\ \underset{m=0}{\overset{M}{\sum}}(b_m^+\exp({\mathrm{i} k_m (x-L)})+b_m^-\exp({-\mathrm{i} k_m (x-L)})) \dfrac{\cosh(k_m(z+h))}{\cosh(k_m h)} & x > L. \end{cases} \end{align}
The wavenumbers 
$k_m$ are the roots 
$k$ of the dispersion equation
\begin{equation} g k \tanh(kh) = \omega^2, \end{equation}
where 
$k_{0}$ is the positive, real root and 
$k_{m}$ (
$m \geq 1$) are purely imaginary, in the upper half of the complex plane and ordered in increasing magnitude. Thus, amplitudes with superscripts 
$+$/
$-$ are related to wave modes that propagate (
$m=0$) or decay (
$m>0$) rightwards/leftwards.
The scattering properties of WEC 
$n$ (the heaving buoy plus PTO, see Appendix A) are defined by 
$(M+1)\times (M+1)$ reflection and transmission matrices, 
${{\boldsymbol{\mathsf{R}}}}_{(n)}$ and 
${{\boldsymbol{\mathsf{T}}}}_{(n)}$, respectively. These map the incoming waves on WEC 
$n$ (
$a^{+}_{m}$ and 
$b^{-}_{m}$) to the outgoing waves (
$a^{-}_{m}$ and 
$b^{+}_{m}$). The relationships can be expressed in terms of a scattering matrix, 
${{\boldsymbol{\mathsf{S}}}}_{(n)}$, such that
\begin{equation} \begin{bmatrix} \boldsymbol{a^{-}}_{(n)} \\ \boldsymbol{b^{+}}_{(n)} \end{bmatrix} = \begin{bmatrix} {{\boldsymbol{\mathsf{R}}}}_{(n)} & {{\boldsymbol{\mathsf{T}}}}_{(n)} \\ {{\boldsymbol{\mathsf{T}}}}_{(n)} & {{\boldsymbol{\mathsf{R}}}}_{(n)} \end{bmatrix} \begin{bmatrix} \boldsymbol{a^{+}}_{(n)} \\ \boldsymbol{b^{-}}_{(n)} \end{bmatrix},\quad\text{where}\ {{\boldsymbol{\mathsf{S}}}}_{(n)}\equiv \begin{bmatrix} {{\boldsymbol{\mathsf{R}}}}_{(n)} & {{\boldsymbol{\mathsf{T}}}}_{(n)} \\ {{\boldsymbol{\mathsf{T}}}}_{(n)} & {{\boldsymbol{\mathsf{R}}}}_{(n)} \end{bmatrix} \end{equation}
and the 
$M+1$ column vectors
\begin{equation} \boldsymbol{a^{{\pm}}}_{(n)}\equiv \begin{bmatrix} a^{{\pm}}_{n,0}\\ \vdots \\ a^{{\pm}}_{n,M} \end{bmatrix} \quad\text{and}\quad \boldsymbol{b^{{\pm}}}_{(n)}\equiv \begin{bmatrix} b^{{\pm}}_{n,0}\\ \vdots \\ b^{{\pm}}_{n,M} \end{bmatrix}, \end{equation}contain the amplitudes, where additional subscripts have been included to specify that the amplitudes belong to unit cell n. In general, the PTO parameters differ between the WECs in the array, so that each of the WECs has a different reflection and transmission matrix.
3.2. WEC array
Let the scattering matrix for WEC 
$p$ to WEC 
$q$ be
\begin{equation} {{\boldsymbol{\mathsf{S}}}}_{(p,q)} \equiv \begin{bmatrix} {{\boldsymbol{\mathsf{R}}}}_{(p,q)}^- & {{\boldsymbol{\mathsf{T}}}}_{(p,q)}^- \\ {{\boldsymbol{\mathsf{T}}}}_{(p,q)}^+ & {{\boldsymbol{\mathsf{R}}}}_{(p,q)}^+ \end{bmatrix}, \end{equation}which includes multiple scattering effects between the WECs. Its reflection and transmission matrices are given by (Bennetts & Squire Reference Bennetts and Squire2009)
$$\begin{gather} {{\boldsymbol{\mathsf{R}}}}_{(p,q)}^-= {{\boldsymbol{\mathsf{R}}}}_{(p,q-1)}^{-}+ {{\boldsymbol{\mathsf{T}}}}_{(p,q-1)}^+[{{\boldsymbol{\mathsf{I}}}} - {{\boldsymbol{\mathsf{R}}}}_{(q)}^- {{\boldsymbol{\mathsf{R}}}}_{(p,q-1)}^+]^{{-}1} {{\boldsymbol{\mathsf{R}}}}_{(q)}^- {{\boldsymbol{\mathsf{T}}}}_{(p,q-1)}^-, \end{gather}$$
$$\begin{gather}{{\boldsymbol{\mathsf{T}}}}_{(p,q)}^-= {{\boldsymbol{\mathsf{T}}}}_{(p,q-1)}^- [{{\boldsymbol{\mathsf{I}}}} - {{\boldsymbol{\mathsf{R}}}}_{(q)}^- {{\boldsymbol{\mathsf{R}}}}_{(p,q-1)}^+]^{{-}1} {{\boldsymbol{\mathsf{T}}}}_{(q)}^-, \end{gather}$$
$$\begin{gather}{{\boldsymbol{\mathsf{R}}}}_{(p,q)}^+= {{\boldsymbol{\mathsf{R}}}}_{(q)}^+{+} {{\boldsymbol{\mathsf{T}}}}_{(q)}^- [{{\boldsymbol{\mathsf{I}}}} - {{\boldsymbol{\mathsf{R}}}}_{(q)}^- {{\boldsymbol{\mathsf{R}}}}_{(p,q-1)}^+]^{{-}1} {{\boldsymbol{\mathsf{R}}}}_{(p,q-1)}^+ {{\boldsymbol{\mathsf{T}}}}_{(q)}^+, \end{gather}$$
$$\begin{gather}\text{and}\quad {{\boldsymbol{\mathsf{T}}}}_{(p,q)}^+= {{\boldsymbol{\mathsf{T}}}}_{(q)}^+ [{{\boldsymbol{\mathsf{I}}}} - {{\boldsymbol{\mathsf{R}}}}_{(q)}^- {{\boldsymbol{\mathsf{R}}}}_{(p,q-1)}^+]^{{-}1} {{\boldsymbol{\mathsf{T}}}}_{(p,q-1)}^+, \end{gather}$$
where 
${{\boldsymbol{\mathsf{I}}}}$ is the 
$(M+1)\times (M+1)$ identity matrix. A recursive algorithm (Bennetts & Squire Reference Bennetts and Squire2009) is used to calculate the reflection and transmission matrices, 
${{\boldsymbol{\mathsf{R}}}}_{(1,n)}^{\pm }$ and 
${{\boldsymbol{\mathsf{T}}}}_{(1,n)}^{\pm }$, for 
$n=2,3,\ldots, N$, and initialised with 
${{\boldsymbol{\mathsf{S}}}}_{(1,1)}={{\boldsymbol{\mathsf{S}}}}_{(1)}$ (3.3). Wave amplitudes within the array (to the right of WEC 
$n$ for 
$n =1,2,\ldots, N-1$) are determined, based on the incident wave amplitude, using a combination of the left-to-right and right-to-left scattering matrices (Bennetts & Squire Reference Bennetts and Squire2009), such that
$$\begin{gather} \boldsymbol{a^+}_{(n+1)} = [{{\boldsymbol{\mathsf{I}}}} - {{\boldsymbol{\mathsf{R}}}}_{(1,n-1)}^+{{\boldsymbol{\mathsf{R}}}}_{(N,n)}^-]^{{-}1} [{{\boldsymbol{\mathsf{T}}}}_{(1,n-1)}^+ \boldsymbol{a^+}_{(1)} + {{\boldsymbol{\mathsf{R}}}}_{(1,n-1)}^+{{\boldsymbol{\mathsf{T}}}}_{(N,n)}^-\boldsymbol{b^-}_{(N)}], \end{gather}$$
$$\begin{gather}\text{and}\quad \boldsymbol{b^-}_{(n)} = [{{\boldsymbol{\mathsf{I}}}} - {{\boldsymbol{\mathsf{R}}}}_{(N,n)}^- {{\boldsymbol{\mathsf{R}}}}_{(1,n-1)}^+]^{{-}1} [{{\boldsymbol{\mathsf{R}}}}_{(N,n)}^-{{\boldsymbol{\mathsf{T}}}}_{(1,n-1)}^+\boldsymbol{a^+}_{(1)} + {{\boldsymbol{\mathsf{T}}}}_{(N,n)}^-\boldsymbol{b^-}_{(N)}]. \end{gather}$$3.3. Wide-spacing approximation
The wide-spacing approximation is used, in which evanescent modes (
$m \geq 1$) are neglected in interdevice interactions (Linton & McIver Reference Linton and McIver2001). Consequently, only the reflection and transmission coefficients associated with propagating modes (
$m=0$) are retained in scatterings relations (
${{\boldsymbol{\mathsf{R}}}}$ and 
${{\boldsymbol{\mathsf{T}}}}$ are reduced to scalars in (3.3)–(3.11)), i.e. it is assumed that
\begin{align} \phi_{(n)} &\sim \frac{\cosh(k_0(z+h))}{\cosh(k_0 h)} \nonumber\\ &\quad \times\begin{cases} a_{(n)}^+\exp({\mathrm{i} k_0 (x-\varOmega^L_{(n)})})+a_{(n)}^- \exp({-\mathrm{i} k_0(x-\varOmega^L_{(n)})}) & x \rightarrow \varOmega^L_{(n)}, \\ b_{(n)}^+\exp({\mathrm{i} k_0 (x-\varOmega^R_{(n)})})+b_{(n)}^- \exp({-\mathrm{i} k_0 (x-\varOmega^R_{(n)})}) & x \rightarrow \varOmega^R_{(n)}, \end{cases} \end{align}
where 
$\varOmega ^L_{(n)}= (2n-3)W/2$ and 
$\varOmega ^R_{(n)}= (2n-1)W/2$ specify the left- and right-hand edges of 
$\varOmega _{(n)}$, respectively. Consequently, the far field (
$x \rightarrow \pm \infty$) can be obtained using the 
$2\times 2$ global scattering matrix, as
\begin{equation} \phi \sim \frac{\cosh(k_0(z+h))}{\cosh(k_0 h)} \times \begin{cases} \exp({\mathrm{i}k_0x})+R\exp({-\mathrm{i}k_0x}) & x \rightarrow -\infty \\ T\exp({\mathrm{i}k_0x}) & x \rightarrow \infty, \end{cases} \end{equation}
where 
$R$ and 
$T$ are the reflection and transmission coefficients for the array of heaving WECs, which can be used to calculate the absorption of the array (2.11). In the absence of PTO damping, 
$R$ and 
$T$ satisfy the energy conservation identity 
$| R|^{2}+| T|^{2}=1$.
Formally, the wide-spacing assumption requires 
$kd \gg 1$ and 
$d/2L \gg 1$ (Linton & McIver Reference Linton and McIver2001). For the problem considered, accurate solutions are obtained even when these conditions are violated, as demonstrated using the WEC amplitudes,
\begin{equation} \xi_n = (a^+_{(n)}\exp({\mathrm{i}kd/2}) + b^-_{(n+1)}\exp({\mathrm{i}kd/2}))\xi_{(n)}^h, \end{equation}
in figure 2, where a spacing of 
$d=4$ m is applied to a uniform array of three WECs, tuned to resonate at 
$\omega \approx 0.44\ \mathrm {rad}\ \mathrm {s}^{-1}$ (2.13a). (Additional justification is provided in Appendix B.) This spacing is selected to separate resonances of individual WECs from resonances driven by the spacing of WECs in the array (§§ 4–5).
The amplitudes of three non-absorbing WECs (
$\omega _0=0.44\ \mathrm {rad}\ \mathrm {s}^{-1}$) in a uniform array obtained using the wide-spacing approximation (—) when 
$kd \in [0.0372, 0.2024]$ and 
$d/2L = 0.4$, compared with the amplitudes of (a) WEC 1, (b) WEC 2 and (c) WEC 3 obtained when including evanescent modes (
$\cdots$) in WEC-interaction calculations.
4. Uniform arrays, Bloch waves and complex resonances
4.1. A uniform array of non-absorbing WECs
Figure 3(a) shows the surface elevation over a range of frequencies (broader than the target power capture interval) for a uniform array of five non-absorbing (
$b^{PTO}=0$), identically tuned WECs (
$\omega _0=0.44\ \mathrm {rad}\ \mathrm {s}^{-1}>\omega _{lb}$, as described in § 3.3), with spacing 
$d=4$ m (
$W=14$ m). For frequencies below approximately 
$\omega =\omega _{0}$, incident waves propagate through the array (
$|\zeta | \in [0.88,1]$ m beyond WEC 5). Around the resonant frequency, 
$\omega =\omega _{0}$, there is a series of large resonant WEC responses (
$|\xi _n|$ up to 36 m). The resonances mark the transition to the array prohibiting wave propagation, up to approximately 
$\omega =0.65$ rad s
$^{-1}$ (
$|\zeta |\approx 0$ m and 
$|\xi _n| \approx 0$ m beyond WEC 2). Between 
$\omega =0.65$ rad s
$^{-1}\equiv \omega _{ub}$ and 
$0.82$ rad s
$^{-1}$, there is an alternating pattern of narrow bands of moderately large responses of the array (
$|\xi _n|$ up to 7.5 m) and small responses (near zero) to the incident waves. The array then transitions back to prohibiting wave propagation up to the highest frequency considered.
The modulus of the surface elevation 
$| \zeta |$ for a uniform array of five WECs (
$\omega _0=0.44\ \mathrm {rad}\ \mathrm {s}^{-1}$) is shown on the (a) 
$\omega$–
$x$ axis, with the corresponding WEC amplitudes overlaid on 
$\omega$–
$z$ axes. The local wave field closely resembles the Bloch waves on the corresponding unit cell in (b) the dispersion diagram, with the array supporting propagation in passbands and preventing transmission in bandgaps.
4.2. Bloch waves
The wave field in a given cell can be decomposed into rightward and leftward Bloch wave modes, 
$\psi _{(n)}^+$ and 
$\psi _{(n)}^-$ (
$n=1,2,\ldots,N$), respectively. They are such that
\begin{equation} \psi_{(n)}^\pm(x,0) = \begin{cases} v_{(n)}^\pm \exp({\mathrm{i}k_0(x-\varOmega^L_{(n)})})\\ \quad + v_{(n)}^\mp \exp({-\mathrm{i}k_0(x-\varOmega^L_{(n)})}), & \text{for } x \leq x_L^{(n)} \\ \exp({\pm \mathrm{i}\beta W})(v_{(n)}^\pm \exp({\mathrm{i}k_0(x-\varOmega^R_{(n)})})\\ \quad + v_{(n)}^\mp \exp({-\mathrm{i}k_0(x-\varOmega^R_{(n)})})) & \text{for } x \geq x_R^{(n)}. \end{cases} \end{equation}
The Bloch wavenumber, 
$\beta _{(n)}$, and amplitudes, 
$v_{(n)}^{\pm }$, are calculated from the 
$2 \times 2$ transfer matrix, 
${{\boldsymbol{\mathsf{P}}}}_{(n)}$, which is defined such that (Porter & Porter Reference Porter and Porter2003)
\begin{equation} \begin{bmatrix} b^{+}_{(n)} \\ b^{-}_{(n)} \end{bmatrix} = {{\boldsymbol{\mathsf{P}}}}_{(n)} \begin{bmatrix} a^{+}_{(n)} \\ a^{-}_{(n)} \end{bmatrix}\quad \text{where}\ {{\boldsymbol{\mathsf{P}}}}_{(n)} = \frac{1}{T_{(n)}} \begin{bmatrix} T_{(n)}^2 - R_{(n)}^2 & R_{(n)} \\ -R_{(n)} & 1 \end{bmatrix}. \end{equation}The transfer matrix is diagonalised as
\begin{equation} {{\boldsymbol{\mathsf{P}}}}_{(n)} = \begin{bmatrix} v^+_{(n)} & v^-_{(n)} \\ v^-_{(n)} & v^+_{(n)} \end{bmatrix} \begin{bmatrix} \mu_{(n)} & 0 \\ 0 & \mu_{(n)}^{{-}1} \end{bmatrix} \begin{bmatrix} v^+_{(n)} & v^-_{(n)} \\ v^-_{(n)} & v^+_{(n)} \end{bmatrix}^{{-}1}, \end{equation}
so that the Bloch wavenumber is calculated from the eigenvalue of the transfer matrix, 
$\beta _{(n)} = -\mathrm {i}\ln (\mu _{(n)})/W$, and the amplitudes are the entries of the eigenvectors. Note the symmetry of the rightward and leftward Bloch waves, which is due to the symmetry of the unit cells.
The band structure of the unit cell is visualised by plotting the real parts of the Bloch wavenumbers against frequency. For the uniform array, all unit cells are identical, and so the band structure of any unit cell is representative of the array. Passbands are defined by real-valued Bloch wavenumbers, and indicate frequency ranges over which Bloch waves propagate across the unit cell. Bandgaps are defined by complex-valued Bloch wavenumbers (
${\rm Re}(\beta _{(n)})=0$ or 
${\rm \pi}$), and indicate frequency intervals over which Bloch waves decay across the unit cell.
For the uniform array and frequency range considered in § 4.1, there are two passbands and two bandgaps (figure 3b). The first passband corresponds to the low-frequency interval over which the incident waves propagate through the array, plus the narrow interval of resonances around 
$\omega _{0}$. The first bandgap is connected with the resonances of the individual WECs and is often referred to as the local resonance bandgap. It starts just above the resonances of the individual WECs, and corresponds to the lowest-frequency interval over which there is no propagation through the array, i.e. its upper bound is 
${\approx }\omega _{ub}$. The second passband covers the interval of alternating narrow bands of large and small responses. Over this interval, waves propagate along the array but partially reflect at the ends of the array and constructively or destructively interfere following rereflections. The second bandgap is connected with the WEC spacing along the array and is referred to as the Bragg bandgap, as it is created by Bragg resonance. It extends to the highest frequencies considered, and also corresponds to an interval over which there is no propagation through the array.
4.3. Complex resonances
The modulus of the transmission coefficient squared for the array, 
$|{}T|^{2}$, i.e. the proportion of transmitted energy, also shows the band structure (figure 4b). Almost full transmission (
$|{}T|^{2}\approx 1$) occurs in the first passband at frequencies up to 
$\omega \approx {0.3}$ rad s
$^{-1}\equiv \omega _{lb}$. There is a sequence of increasingly sharp, narrow peaks and troughs at the high-frequency end of the first passband, just below the resonant frequency, 
$\omega =\omega _{0}$. Transmission is zero in both bandgaps and oscillates in the second passband that separates the bandgaps. The modulus of the reflection coefficient squared, 
$|{}R|^{2}$, shows the band structure in a complementary manner, i.e. near zero reflection or oscillations in the passbands and full reflection in the bandgaps, due to the energy conservation identity (2.11) for the non-absorbing array (
${\alpha \equiv 0}$).
(a) The structure of transmission coefficient (
$T\in \mathbb {C}$) as a function of 
$\omega \in \mathbb {C}$ for a uniform array of five WECs with 
$b^{PTO}=0$, visualised by colour-coding the phase (
$\arg (T)$) to form a phase portrait (Wegert Reference Wegert2012), where the magnitude of the phase is represented by hue. Poles (WEC-resonances; marked with an X) and zeros are identified by rapid phase changes and distinguished by the ordering of colours in the anticlockwise direction (Wegert Reference Wegert2012). A white 
$\bullet$ marks the location where the complex zeros coalesce on 
$\operatorname {Im}(\omega )=0$ at 
$\omega \approx 0.45\ \text {rad}\ \text {s}^{-1}$, resulting in (b) 
$|T|^2=0$ and (c) 
$|R|^2=1$ for 
$\omega \in \mathbb {R}^+$, and the first bandgap in figure 3(b). Resonances in the complex plane produce a local maximum of 
$|T|^2$ and a local minimum of 
$|R|^2$ for 
$\omega \in \mathbb {R}^+$. The real-valued WEC-resonance is denoted 
$\omega _0$.
Transmission oscillations in the passbands are governed by the structure of the transmission coefficient in the complex frequency plane (
$\omega \in \mathbb {C}$; figure 4a). Each peak in transmission on the real frequency axis is associated with a pole in the transmission coefficient in the lower half of the complex plane, which are known as complex resonances (Romero-García et al. Reference Romero-García, Theocharis, Richoux and Pagneux2016; Meylan & Fitzgerald Reference Meylan and Fitzgerald2018; De Chowdhury, Bennetts & Manasseh Reference De Chowdhury, Bennetts and Manasseh2023). The closer the complex resonances are to the real frequency axis, the sharper the transmission peaks on the real axis are. Each complex resonance corresponding to a peak in the first passband is associated with a specific WEC in the array and labelled accordingly. The association indicates that the location of the complex resonance is sensitive to variations in the properties of the WEC (not shown), although all of complex resonances move to a certain degree when the properties of any WEC in the array or the array spacing vary. The resonance associated with WEC 5 is outside the axes limits.
The complex resonances associated with transmission peaks in the second passband are created by wave interference along the array (they have no analogue in the single body problem). Thus, their locations depend predominantly on the WEC spacing – for example, a larger spacing would push them to lower frequencies. There are four overlapping zeros in transmission on the real-frequency axis very close to the resonant frequency, 
$\omega =\omega _{0}$, and at the point where the Bloch wavenumber jumps between branches (shown by the discontinuity in Re
$(\beta )$, figure 3b).
5. Grading the resonant properties of WECs
5.1. Rainbow reflection
Consider the uniform array of five non-absorbing WECs (as in § 4) but in which the stiffness coefficients are graded with distance along the array, so that the WEC-resonant frequencies increase from right to left (from WEC 5 to WEC 1). The grading has the effect of frequency upshifting the first bandgaps (in the unit cells, from right to left), as well as narrowing their frequency ranges, as they are bounded above by the second (Bragg) bandgaps, which are insensitive to the WEC-resonances (figure 5a–c). The grading is sufficiently gradual that the first bandgaps for adjacent unit cells overlap (the resonant frequency for WEC 
$n$ lies in the bandgap associated with WEC (
$n+1$)). This creates a wide effective bandgap for the array that covers the target interval, i.e. array transmission is approximately zero (
$| T|^{2}\approx 0$ and 
$| R|^{2}\approx 1$) for 
$\omega \in [0.3,0.8]\ \text {rad}\ \text {s}^{-1} \supset \omega _{\alpha }$. The effective bandgap manifests as separation of the complex zeros in the phase portrait of 
$T(\omega )$ along the real frequency axis over the target frequency range (figure 6a).
(a–c) Band diagrams of the Bloch waves for an array with 
$W=14$ m (
$d = 4$ m and 
$L = 5$ m), for increasing 
$c^{PTO}$ values and 
$b^{PTO}=0$. Real-valued WEC-resonances 
$\omega _0 \in \mathbb {R}$ are marked by a red x, where (a) 
$\omega _0 = 0.31\ \text {rad}\ \text {s}^{-1}$, (b) 
$\omega _0 = 0.48\ \text {rad}\ \text {s}^{-1}$, (c) 
$\omega _0 =0.72\ \text {rad}\ \text {s}^{-1}$. Increasing 
$c^{PTO}$ shifts the first bandgap (shaded green) to higher frequencies, and reduces the interval of the second passband. The second bandgap (grey) is caused by Bragg resonance. Grading the WEC-resonances in a finite array of five WECs (first, third and fifth WECs correspond to panels (c–a), respectively) forms (d) an effective bandgap on 
$\omega \in [0.3,0.8]\ \text {rad}\ \text {s}^{-1}$, where 
$|R|^2\approx 1$ and 
$|T|^2\approx 0$.
Phase portraits of (a) 
$T(\omega )$ and (b) 
$R(\omega )$ are shown as a function of 
$\omega \in \mathbb {C}$ for the graded array in figure 5. Complex WEC-resonances (X) can be identified from (a) or (b). White circles (o) denote the complex zeros in 
$R$ and 
$T$, and the corresponding 
$|R(\omega )|^2$ and 
$|T(\omega )|^2$ for 
$\omega \in \mathbb {R}$ are superimposed. Vertical dotted lines demarcate the effective bandgap induced by the grading (figure 5), and red xs correspond to the real-valued WEC-resonance of each WEC. The complex resonances preceding the Bragg bandgap are marked by 
$\blacklozenge$ in (a), with the corresponding pole–zero pairs in (b) not visible at the current scale.
The WEC-resonances are graded from high-to-low from left-to-right to create the rainbow reflection effect, such that incident waves at frequencies in the target range penetrate into the array, with penetration distances that depend on the wave frequency. (Grading low-to-high would prevent propagation into the array.) For a given 
$\omega \in \omega _{\alpha }$, an incident wave will propagate into the graded array until reaching a cutoff point, which is approximately where the local Bloch wave has zero group velocity, and occurs in the unit cell at which the frequency lies in the corresponding bandgap. Put simply, the Bloch wave propagates until it reaches a WEC with a comparable resonant frequency. The wave is amplified just before the cutoff point due to the slow-down of wave energy transport, which creates a near-resonant (i.e. finite) response.
The near-resonances are connected with complex resonances in the reflection coefficient in the lower-half complex frequency plane (black crosses in figure 6b), which have real parts approximately covering the target interval. The near full reflection over the real frequency interval creates near symmetry in the real frequency axis of the reflection phase and reciprocity of the modulus, so that complex zeros of reflection occur at approximately the complex conjugates of the complex resonant frequencies (figure 7a; Bennetts & Meylan (Reference Bennetts and Meylan2021)). (Note the symmetry is weakest for the WEC 4 pole–zero pair, where full reflection is not achieved.) Manipulating the properties of the pole–zero pairs in complex frequency space provides a novel means to analyse and control the array response based on the resonant properties of individual WECs. The PTO stiffness is used to influence the locations of real parts of the pole–zero pairs (figure 7b). However, interactions between the phase structures supported by the pole–zero pairs (i.e. their bandwidths) prevent full control. The location of the pole–zero pair associated with WEC 1 is restricted towards high frequencies by pole–zero pairs related to complex resonances in the second passband.
Phase portrait of 
$R(\omega )$, when 
$\omega \in \mathbb {C}$, for a single WEC in a graded array (of five WECs). When 
$b^{PTO} =0$, (a) the complex zero (white circle) is located above 
$\operatorname {Im}(\omega )=0$, and the complex WEC-resonance (X) below 
$\operatorname {Im}(\omega )=0$. Increasing 
$c^{PTO}$ (b) moves the WEC-resonance to the right, tuning the WEC to a higher frequency. Positive 
$b^{PTO}$ (c) moves the complex zero towards 
$\operatorname {Im}(\omega )=0$. Perfect absorption (d) is achieved when 
$b_*^{PTO}>0$ places the complex zero on 
$\operatorname {Im}(\omega )=0$. Simultaneously, a near-zero minimum of 
$|R|^2$ is obtained for 
$\omega \in \mathbb {R}$ (overlaid).
5.2. Rainbow absorption
In complex frequency space, adding PTO damping to WEC 
$n$ moves the corresponding complex zero of reflection towards the real axis and the pole away from it (figure 7c). When the pole–zero pairs are sufficiently separated, a value of the PTO damping can be found such that the complex zero of reflection lies on the real frequency axis (figure 7d). This approach is used as a basis to create rainbow absorption, i.e. high absorption over a wide frequency range, where the absorption peaks are spatially separated according to frequency.
The rainbow absorption approach involves adding the WECs to the array one at a time, starting with the rightmost WEC (WEC 5 in this example) and moving leftwards, so that the number of WECs in the array increases with each iteration. The rightmost WEC is used only to create a zero in transmission close to the low-frequency end of the target frequency interval, rather than as an absorbing device (figure 8a). As each WEC is added to the array, its PTO stiffness is tuned to a chosen frequency (noting that pole–zero pair interactions mean the chosen frequency is not arbitrary) and then its damping is tuned to create a zero in reflection at a nearby (real) frequency. The tuning of the additional WEC PTO tends to slightly frequency-downshift the existing pole–zero pairs and displace the complex zeros from the real frequency axis. Thus, the PTOs of the existing WECs are simultaneously retuned increasingly far apart in the complex plane as the WEC-resonance frequency and corresponding resonance bandwidth increase (Romero-García et al. Reference Romero-García, Theocharis, Richoux and Pagneux2016), which forces the reflection zeros to become spaced increasingly far apart.
The transmission (–, blue) and reflection (–, red) of the array versus frequency as (a) WEC 
$5$, (b) WEC 
$4$, (c) WEC 
$3$, (d) WEC 
$2$ and (e) WEC 
$1$ are added to the array (from right to left, with 
$W=14$ m) and tuned. The location of 
$|R|^2\approx 0$ associated with each WEC is denoted 
$\omega _n$, and corresponds to the location of complex zeros (white circles) of the absorbing WECs in the phase portrait of 
$R(\omega )$ for 
$\omega \in \mathbb {C}$ for the graded array shown in ( f). The non-absorbing WEC 5 is denoted 
$\omega _{{low}}$. Complex WEC-resonances are denoted X and open circles denote the zeros when 
$b_{(n)}^{PTO}=0$ (
$n=1,2,\ldots,N$).
The process is accomplished manually for WEC 4 down to WEC 2 (figure 8b–d). However, it becomes increasingly challenging to separate the pole–zero pairs on the restricted interval (bounded above by the complex resonances associated with the second passband). The increasing distance between each pole–zero pair as frequency increases is responsible for increasing bandwidths around the zeros and the need for increasing separation between the zeros (e.g. figure 8e). Only four near zeros of reflection (
$|R|^{2}\approx 0$) are found manually after WEC 1 is added. This solution is used as an initial guess in an optimisation algorithm (the sequential quadratic programming algorithm in the MATLAB function fmincon) with the objective
\begin{equation} \text{minimise } \int_{\omega_{lb}}^{\omega_{ub}} |R(\omega)|^2 \ \mathrm{d}\omega \quad\text{with respect to}\ c^{PTO}_{(n)}\ {\rm and}\ b^{PTO}_{(n)}\quad (n=1,2,\ldots,N), \end{equation}to obtain four zeros in reflection (figure 8e).
The above approach yields reflection zeros at 
$N-1=4$ real frequencies, which corresponds to four complex zeros of reflection being translated to the real frequency axis by adding PTO damping (figure 8f). The reflected energy remains small between the zeros (
$|R|^{2}<0.003$), such that near-perfect broadband absorption of the incident wave energy is almost achieved over the target interval (
$\hat {\alpha }=0.984$). A small amount of transmitted energy away from its low-frequency zero prevents near-perfect absorption (
${\int } |T(\omega )|^2 \ \mathrm {d}\omega = 0.014$ on 
$\omega _{\alpha }$).
6. Near-perfect broadband absorption
Both 
$|R|^2$ and 
$|T|^2$ should be minimised simultaneously over the target frequency interval (
$\hat {\alpha }$ maximised) to achieve near-perfect, broadband power capture. Thus, optimisation is performed with the objective to
\begin{equation} \text{minimise } \int_{\omega_{lb}}^{\omega_{ub}} \{|R(\omega)|^2 + |T(\omega)|^2\}\ \mathrm{d}\omega \quad\text{with respect to}\ c^{PTO}_{(n)}\ {\rm and}\ b^{PTO}_{(n)} \end{equation}
for 
$n=1,\ldots,N$. The aim is to use the theory for rainbow absorption (§§ 4–5) to determine an initialisation strategy and constraints on the optimisation parameters that gives an efficient automated algorithm for near-perfect absorption.
6.1. Generic algorithm
To cover the target power capture interval, 
$\omega _{\alpha }=[\omega _{lb},\omega _{ub}]$, the stiffness coefficients of the first and last WECs in the array (
$c^{PTO}_{(1)}$ and 
$c^{PTO}_{(N)}$, respectively) are initialised so that the corresponding resonances (
$\omega ^{(1)}_0$ and 
$\omega ^{(N)}_0$) coincide with the upper and lower bounds of the target interval (
$\omega _{ub}$ and 
$\omega _{lb}$). The WEC 
$N$ resonance induces a zero in transmission at a frequency just above the lower bound (
$\omega _{lb}$, as in figure 8a), and the stiffness coefficient of the penultimate WEC (
$c^{PTO}_{(N-1)}$) is initialised to place its resonance (
$\omega ^{(N-1)}_0$) at the frequency of the transmission zero. The remaining stiffness coefficients are initialised to space the resonances of WEC 2 to WEC 
$(N-2)$ evenly between 
$\omega _0^{(N-1)}$ and 
$\omega _0^{(1)}$, such that
\begin{equation} \omega_0^{(N)} < \omega_0^{(N-1)}< \cdots < \omega_0^{(1)}. \end{equation}
The damping coefficients 
$b^{PTO}_{(1)}$, 
$b^{PTO}_{(2)}, \ldots, b^{PTO}_{(N-1)}$ are initialised at the optimal values for the corresponding WECs in isolation (2.13b). The final WEC is non-absorbing (
$b^{PTO}_{(N)}=0$).
To account for the leftward movement of WEC-resonances due to damping and WEC interactions on the restricted interval, the stiffness coefficients 
$c^{PTO}_{(n)}$ (
$n=2,3,\ldots, (N-1)$) are bounded below by their initial values and above by the initial value of the preceding WEC in the array (i.e. the initial value of 
$c^{PTO}_{(n-1)}$). The upper bound for the stiffness coefficient of the first WEC (
$c^{PTO}_{(1)}$) is above the target interval and not an initial value. For arrays of 
$N < 7$ WECs, the upper bound is set to constrain 
$c^{PTO}_{(1)}$ to a frequency interval length comparable to the constraints on 
$c^{PTO}_{(2)}$, provided 
$\omega _0^{(1)} < 0.79\ \text {rad}\ \text {s}^{-1}$ (i.e. the pole–zero pair remains on 
$\omega _{\alpha }$). When 
$N \geq 7$, the upper bound corresponds to 
$\omega _0^{(1)} = 0.72\ \text {rad}\ \text {s}^{-1}$, which is approximately the minimum WEC-resonance required to maintain a bandgap at 
$\omega _{ub}$ (depending on pole–zero pair interactions). The damping coefficients 
$b^{PTO}_{(n)}$ (
$n=1,2,\ldots, (N-1)$) are constrained below such that they do not become negative (i.e. do not add energy to the system) and above such that they do not exceed twice the initial value, which does not limit the optimisation in practice (in tests conducted).
The generic algorithm creates near-perfect broadband absorption from an array of five WECs (
$\hat {\alpha }=0.990$), which slightly improves on the absorption given by the graded array of five WECs with zeros in reflection (
$\hat {\alpha }=0.984$; § 5.2). In comparison with the array with zeros in reflection (table 1), the WEC-resonances are graded more gradually, over a narrower frequency range (table 2). This slightly increases the absorption over the majority of the target power capture interval (figure 9a). In complex-frequency space, the zeros in reflection are slightly displaced from the real frequency axis (figure 9b) to balance reduced transmission. However, the absorption is lower than the array with zeros in reflection towards the ends of the target frequency interval, particularly at the high frequency end, which is influenced by the complex zeros corresponding to the second passband.
The PTO parameters and resulting WEC-resonance of each WEC in the graded array of five WECs shown in figure 8.
Optimised PTO parameters for broadband absorption by the (generic) graded array of five WECs in figure 9(b).
The absorption of the optimised, graded array (—) of five WECs is shown in (a) with the corresponding phase portrait of 
$R(\omega )$ for 
$\omega \in \mathbb {C}$ in (b). The absorption of the array in figure 8( f) is overlaid (
$\cdots$) in (a), with the associated complex zeros in reflection marked by white os in (b). Absorption improves over the majority of 
$\omega _\alpha$ as a result of a more gradual grading in the generic solution (reduces 
$| T |^2$), which is facilitated by increasing the number of WECs to 
$N=10$ in (a), using the array of 
$N=10$ WECs in figure 10(b).
The average absorption on 
$\omega _\alpha$ (a) increases with the number of WECs in both algorithms. Constraining the PTO parameters based on problem knowledge in the hybrid algorithm reduces run time and produces almost identical absorption (
$\cdots$) to the generic algorithm (—). However, the array properties differ, as demonstrated for an array of 10 WECs using the phase portraits of 
$R(\omega )$ for 
$\omega \in \mathbb {C}$ corresponding to the solution of the generic (b), and hybrid (c) algorithms, respectively. Pole–zero pairs are more separated in (c) allowing for more near-zeros in reflection to be obtained compared with (b). In both algorithms, the pole–zero pair of an absorbing WEC is pushed beyond axes limits (WEC 
$10$ lies below the axes limits).
The grading of WEC-resonances partially controls the lower bound of the second passband. By grading the array more gradually to reduce transmission, this lower bound is downshifted, causing reflection and transmission to rise near 
$\omega _{ub}$ which decreases absorption. Increasing the number of WECs (e.g. to 
$N=10$) improves the absorption at the high-frequency end of the target interval (and more generally across the interval, figure 9a), by facilitating a gradual grading across the entire target interval, which simultaneously reduces the width of the second passband, and the transmission over 
$\omega _\alpha$.
The generic algorithm becomes prohibitively expensive as the number of WECs increases. For 
$N=6$, the algorithm takes 8 h on an Intel(R) Xeon(R) CPU E5-2699 v3 @ 2.30 GHz with 251 GB RAM. Moreover, MATLAB exceeds the default tolerated number of function evaluations for 
$N\geq 7$. Consequently, an approach was developed for 
$N\geq {}7$, in which approximate solutions were obtained using a coarse frequency resolution (
$\omega = 0.01257\ \text {rad}\ \text {s}^{-1}$) after the tolerated number of function evaluations, and then used as initial guesses in a second application of the algorithm with a higher frequency resolution (
$\omega = 0.006283\ \text {rad}\ \text {s}^{-1}$). In the second application, the stiffness coefficient constraints are updated to be bounded below by the initial guess, and above by the preceding WEC resonance. Here, WEC 1 is constrained such that 
$\omega _0^{(1)} \in [0.72,0.78]\ \text {rad}\ \text {s}^{-1}$. The array of 
$N = 10$ WECs requires three applications of the algorithm.
6.2. Hybrid algorithm
The generic algorithm is adapted into a hybrid algorithm that incorporates more knowledge gained from the analysis of the graded array (§ 5.2) to accelerate the optimisation process. To cater for the leftward shift of the WEC-resonances and ensure absorption at 
$\omega _{ub}$, WEC 
$1$ is initialised and constrained such that its real-valued resonance lies above the power capture interval, with 
$c^{PTO}_{(1)}$ constrained around 
$\omega \approx 0.75\ \text {rad}\ \text {s}^{-1}$. Similarly, 
$c^{PTO}_{(N-1)}$ is constrained to a narrow interval above 
$\omega _{lb}$, around 
$\omega \approx 0.33\ \text {rad}\ \text {s}^{-1}$. Further, the initial damping constants for the first and second last WECs in the array (WEC 1 and WEC 
$(N-1)$) are set lower than optimal values for the WECs in isolation to provide greater control over the location of pole–zero pairs of reflection and transmission in complex frequency space, i.e. consistent with observations in table 1. The stiffness and damping coefficients, 
$c^{PTO}_{(n)}$ (
$n=2,\ldots, (N-2)$) and 
$b^{PTO}_{(n)}$ (
$n=1,2,\ldots, N$), respectively, are constrained analogously to the generic algorithm, with the exception of the upper bound for WEC 
$2$, which is constrained above by the lower constraint bound of WEC 1.
The hybrid algorithm gives almost identical broadband absorption to the generic algorithm for 
$N\geq 5$ (figure 10a). It does so at approximately half the runtime of the generic algorithm when 
$N=6$. The algorithm was reapplied (as in § 6.1) for 
$N=10$, with the constraints held fixed. The differences between the algorithms are manifest in the structure of the reflection coefficient in complex frequency space, for which the generic algorithm stacks the pole–zero pairs at low frequencies (figure 10b), whereas the hybrid algorithm keeps them separated (figure 10c). The increased separation allows for greater control over the location of pole–zero pairs, as evidenced by additional near-zeros in reflection in the hybrid solution. A pole–zero pair associated with an absorbing WEC is pushed beyond the axes limits in both algorithms as a result of interactions on the restricted interval.
6.3. Multiple degrees of freedom
Allowing the buoys to surge and pitch (as non-absorbing degrees of freedom), as well as heave (the absorbing degree of freedom), has little impact on the near-perfect, broadband absorption achieved by the optimised arrays. For example, the average absorption of the array of five WECs produced by the generic algorithm (table 2) decreases only from 
$\hat {\alpha } = 0.990$ to 
$\hat {\alpha } = 0.988$ when the surge and pitch motions are released (figure 11). Reoptimising the PTO parameters for heave motion for the problem including surge and pitch (starting from the optimised solution for heave only) brings the average absorption back to 
$\hat {\alpha } = 0.990$. The optimal PTO parameters for the multiple-degrees-of-freedom case (table 3) differ from the single-degree-of-freedom case (table 2) by only 5 %, on average. Small differences for the absorption versus frequency are visible for the two optimised arrays (figure 11).
Releasing surge and pitch (
$\cdots$) in the optimal generic solution (–) for the array of 
$N = 5$ WECs in table 2 reduces the average absorption by 0.0017 (
$\hat {\alpha } = 0.988$). Reoptimising the solution to account for the uncontrolled surge and pitch motions (
$\textbf{-}{\bullet}\textbf{-}$) over the target interval results in an average absorption of 
$\hat {\alpha } = 0.990$.
Optimised PTO parameters for the graded array of five WECs in figure 11 when surge and pitch are released as uncontrolled degrees of freedom.
6.4. Ocean wave spectra
The JONSWAP ( Joint North Sea Wave Project (Hasselmann et al. Reference Hasselmann1973)) spectrum is commonly used to model realistic (irregular, random) sea states (Chakrabarti Reference Chakrabarti2005). The JONSWAP spectrum is defined by the spectral density function
\begin{equation} S(\omega) = \frac{0.0081 g^2}{\omega^5} \exp({-1.25({\omega_p}/{\omega})^4}) \gamma^r\quad\text{where}\ r = \exp({-{(\omega-\omega_p)^2}/{2\sigma^2\omega_p^2}}) \end{equation}and the spectral width is
\begin{equation} \sigma =\begin{cases} 0.007 & \omega < \omega_p, \\ 0.009 & \omega \geq \omega_p. \end{cases} \end{equation}
The peak enhancement factor 
$\gamma$ characterises the sharpness of the spectrum around the peak frequency 
$\omega _p = 2{\rm \pi} /T_p$ (peak period 
$T_p$, with maximum energy in the spectrum). Smaller values are indicative of broader banded sea states, with 
$\gamma \leq 2.4$ generally suited to coastal regions (Mazzaretto, Menéndez & Lobeto Reference Mazzaretto, Menéndez and Lobeto2022).
Let the JONSWAP spectrum be approximated by a finite sum of monochromatic waves (Cruz Reference Cruz2008), at frequencies 
$\omega _{0}=0.22$ rad s
$^{-1}$, 
$\omega _{1}=\omega _{0}+\Delta \omega, \ldots, \omega _{200} =\omega _{0}+200\Delta \omega =1.26$ rad s
$^{-1}$ (
$\Delta \omega = 0.0052$), then the corresponding absorption of the spectrum is defined as
\begin{equation} \alpha_s(\omega_i) = \sum_{i = 0}^{200}\alpha(\omega_i) | A(\omega_i) |^2 \quad\text{where}\ A(\omega_i) = \sqrt{2S(\omega_i) \Delta \omega}. \end{equation}The proportion of the spectrum absorbed is
\begin{equation} \hat{\alpha}_{s} = \frac{\int\alpha_s(\omega)\mathrm{d}\omega}{\int |A(\omega)|^2 \ \mathrm{d}\omega}. \end{equation}
The optimised array of five WECs (generic algorithm, table 2) captures 
$\hat {\alpha }_s=0.95$ of the incident energy in a JONSWAP spectrum with a narrow-banded sea state (
$\gamma =3.3$) and a peak period 
$T_p = 17$ s, as the spectral peak lies within the target power capture interval (figure 12a). The array still captures over 
$\hat {\alpha }_{s}>0.75$ of the incident energy when the spectral peak is moved close to the upper boundary of the target interval (
$T_{p}=10$ s, figure 12a). The array captures over 85 % of the incident energy for peak periods of approximately 11–20 s. For a broader sea state (
$\gamma =1.54$) this 
$T_{p}$-interval is shifted to 12–21 s, and the peak absorption is reduced slightly to 
$\hat {\alpha }_s=0.936$ at 
$T_{p}=17$ s.
The absorption of spectra (
$\gamma = 3.3$) with peak periods located in 
$\omega _\alpha$ is shown in (a) for the graded array of five WECs in table 2. The proportion of spectra captured by the array is shown as a function of peak period in (b) for a constant significant wave height of 1 m. Absorption decreases as the peak period shifts farther from the targeted interval and the energy of the spectrum is located outside the designed frequency range of the array. On average, 
$\widehat {\alpha _s} > 0.85$ (
$\cdots$) on 
$\omega _\alpha$ for both 
$\gamma = 1.54$ and 
$\gamma = 3.3$.
7. Absorption of Bloch modes in time domain
An incident wave packet of unit amplitude in the time domain is specified via the Fourier transform of its ocean surface displacement, 
$\eta (x,t)$. The Fourier transform maps between the time and frequency domains and is expressed in terms of wavenumber, 
$k(\omega )$, as
\begin{equation} {F}\lbrace \eta(x,0)\rbrace\equiv \hat{f}(k) = \frac{1}{\sigma\sqrt{2{\rm \pi}}}\exp({-(k-k_0)^2/{2\sigma^2}}), \end{equation}
where 
$k_{0}$ is the prescribed central wavenumber and 
$\sigma$ is the prescribed packet width (in wavenumber space). The response of the WECs and the free surface elevation are calculated from the frequency domain solutions, as, respectively,
\begin{equation} \varXi_{n}(t)=\operatorname{Re}\left\lbrace \frac{1}{\rm \pi} \int^\infty_0 \hat{f}(k)\xi_{n}(k) \exp({-\mathrm{i}\omega(k) t})\ \mathrm{d}k\right\rbrace \quad\text{for}\ n=1,2,\ldots,N, \end{equation}and
\begin{equation} \eta(x,t) = \operatorname{Re}\left\lbrace \frac{1}{\rm \pi} \int^\infty_0 \hat{f}(k)\zeta(x:k) \exp({-\mathrm{i}\omega(k) t})\ \mathrm{d}k\right\rbrace, \end{equation}where the dependencies of frequency domain quantities on wavenumber are made explicit.
In the frequency domain, the wave field in each unit cell is decomposed into the rightward and leftward Bloch modes, by writing
\begin{equation} \phi_{(n)}(x,z:k) \approx w^{+}_{(n)}(k)\psi_{(n)}^{+}(x,z:k) +w^{-}_{(n)}(k)\psi_{(n)}^{-}(x,z:k) \quad\text{for}\ n=1,\ldots,N, \end{equation}
where the weights 
$w^{\pm }_{(n)}$ (
$n=1,\ldots,N$) are calculated as
\begin{equation} \begin{bmatrix} w^{+}_{(n)} \\ w^{-}_{(n)} \end{bmatrix}=\begin{bmatrix} v^+_{(n)} & v^-_{(n)} \\ v^-_{(n)} & v^+_{(n)} \end{bmatrix}^{{-}1} \begin{bmatrix} a^{+}_{(n)} \\ a^{-}_{(n)} \end{bmatrix}. \end{equation}Similarly, the free surface elevation and heave amplitudes are decomposed as
\begin{equation} \zeta(x:k) = w^{+}_{(n)}(k)\psi_{(n)}^{+}(x,0:k) + w^{-}_{(n)} (k)\psi_{(n)}^{-}(x,0:k)\quad \text{for } \varOmega^{L}_{(1)} \leq x \leq \varOmega^{R}_{(N)}, \end{equation}and
\begin{equation} \xi_{n}(k) = w^{+}_{(n)}(k) \gamma^{+}_{(n)}(k) + w^{-}_{(n)}(k) \gamma^{-}_{(n)}(k)\quad\text{for}\ n=1,\ldots,N, \end{equation}where
\begin{equation} \gamma^{{\pm}}_{(n)} = (v^{{\pm}}_{(n)}\exp({\mathrm{i} k_0 (x-\varOmega^L_{(n)})}) + \exp({\pm \mathrm{i}\beta W}) v^{{\mp}}_{(n)}\exp({-\mathrm{i} k_0 (x-\varOmega^R_{(n)})})) \xi_{(n)}^h. \end{equation}Equivalent decompositions of the surface elevation and heave displacements are made in the time domain, with
\begin{equation} \eta(x,t) =\varPsi_{n}^{+}(x,t)+\varPsi_{n}^{-}(x,t)\quad \text{for}\ \varOmega^{L}_{(n)} \leq x \leq \varOmega^{R}_{(n)}, \end{equation}and
\begin{equation} \varXi_n(t)=\varXi^{+}_n(t)+\varXi^{-}_n(t), \end{equation}where
\begin{equation} \varPsi^{{\pm}}_n(x,t) = \operatorname{Re}\left\lbrace \frac{1}{\rm \pi} \int^{\infty}_{0} \hat{f}(k) w^\pm_{(n)}(k)\psi^\pm_{(n)}(x,0:k)\exp({-\mathrm{i}\omega(k) t})\ \mathrm{d}k\right\rbrace \end{equation}and
\begin{equation} \varXi^{{\pm}}_n(t) = \operatorname{Re}\left\lbrace \frac{1}{\rm \pi} \int^{\infty}_{0} \hat{f}(k)w^\pm_{(n)}(k)\gamma^{{\pm}}_{(n)} \exp({-\mathrm{i}\omega(k) t})\ \mathrm{d}k\right\rbrace \end{equation}
for 
$n=1,2,\ldots,N$. Note that the rightward and leftward Bloch modes are not continuous at the interfaces between the unit cells (only their sum is continuous).
As an example, consider the case of the optimised array of 
$N=5$ WECs (table 2), which gives near-perfect broadband absorption over the target interval (
$\hat {\alpha }=0.990$). The spatiotemporal responses to two separate incident wave packets are shown, where the packets are centred around a wavenumber at the resonant frequency of WEC 2 (
$k_{0}=0.03441$, figure 13) and WEC 4 (
$k_{0}=0.01703$, figure 14). Both have a packet width of 
$\sigma = 0.002673\ \text {m}^{-1}$, so that the packet centred at WEC 2 covers the resonant frequencies of WECs 1–3, and the packet centred at WEC 4 covers WECs 2–5. The total responses are shown (figures 13a,b and 14a,b), as well as their decompositions into rightward and leftward Bloch modes (figures 13c,d and 14c,d, and figures 13e, f and 14e, f, respectively). The responses are shown for the non-absorbing array (
$b^{PTO}_{(n)}=0$ for 
$n=1,2,\ldots,5$, figures 13a,c,e and 14a,c,e), as well as the absorbing array (figures 13b,d, f and 14b,d, f). (The artificial periodicity introduced via the discrete implementation of the Fourier transform is 
${\approx }40$ min, which is much greater than the 
${\approx }3$ min duration of the responses shown.)
The spatiotemporal behaviour of the non-absorbing (a,c,e) and absorbing arrays (b,d, f), are shown in (a,b), respectively, when forced by a wave packet centred at 
$k_0=0.03441$. The total wave fields are decomposed into rightward Bloch modes in (c,d), and leftward Bloch modes in (e, f), respectively. The WEC displacements are overlaid on 
$z$–
$t$ axes at the 
$x$-location of WECs in the 
$x$–
$t$ domain. In the non-absorbing array, the rightward Bloch mode cuts off at WEC 
$2$ (c). The leftward mode (e) then drives almost total reflection of frequencies above the cutoff. In the absorbing array, the total wave field (b) and WEC displacements are dominated by the rightward Bloch mode (d), with little excitement of the leftward Bloch mode ( f) through the absorption of incident energy.
The spatiotemporal behaviour of the non-absorbing (a,c,e) and absorbing arrays (b,d, f) when forced by a wave packet centred at 
$k_0 = 0.01703$. The total wave fields are decomposed into rightward Bloch modes in (c,d), and leftward Bloch modes in (e, f), respectively. The WEC displacements are overlaid on 
$z$–
$t$ axes at the 
$x$-location of WECs in the 
$x$–
$t$ domain. The total wave field (a) is dominated by the rightward Bloch mode (c) in the non-absorbing array before the cut off at WEC 
$4$ is reached. The leftward mode (e) drives high reflection above the cutoff. With PTO damping, the total wave field (b) and WEC displacements in the absorbing array are governed predominantly by the rightward Bloch mode (d), with little excitement of the leftward Bloch mode ( f), and thus, near-zero reflection.
For the non-absorbing array, the incident wave packets gradually reduce in amplitude as they propagate through the array due to a series of partial reflections by the WECs (figures 13a and 14a). There is a more abrupt cutoff after the WECs corresponding to the central wavenumbers of the incident packets (WEC 2 in figure 13a and WEC 4 in figure 14a), such that the amplitudes reduce to 
$\approx$20 % of their incident values. The overall reflection is indicated by the durations of the wave signals decreasing with distance into the array (up to the cutoffs). The incident and reflected packets are isolated from one another by the decomposition into Bloch modes, such that the incident wave packets propagate through the array in the form of rightward Bloch modes (figures 13c and 14c) and are reflected back out of the array in the form of leftward Bloch modes (figures 13e and 14e).
The incident wave packets propagate similar distances into the absorbing arrays (figures 13b and 14b), as their non-absorbing counterparts. However, there is almost no reflected wave packet generated, as indicated by similar durations of the wave signal along the array (up to the cutoffs), and the near-resonant WEC responses are reduced by 30 %–40 % compared with the non-absorbing arrays. Thus, absorption is the cause of reductions in the amplitude of the incident packet. The decomposition into Bloch modes gives further evidence of the absence of reflection and dominance of absorption, with the rightward Bloch modes indistinguishable from the total wave fields (figures 13d and 14d) and virtually no leftward Bloch modes excited (figures 13f and 14f).
8. Conclusions and discussion
Arrays of heaving point-absorber WECs have been designed to achieve near-perfect, broadband wave energy absorption, using 2-D linear potential-flow theory. The WECs consisted of floating buoys plus spring–damper PTOs, and used parameters inspired by CorPower Ocean's C4 device. Broadband absorption was demonstrated in both the frequency and time domains over a typical frequency interval for power capture (equivalent to a wave period interval 
$\approx$10–20 s) using time-independent PTO parameters. In particular, an array of only five heaving buoy type WECs was shown to capture 99 % of the incident wave energy over the target interval. The array performance was analysed for irregular sea states, and near-perfect broadband absorption was also demonstrated when surge and pitch motions were released.
The approach to designing near-perfect broadband absorption was based on theories for graded arrays. From an underlying uniform array of non-absorbing WECs, the rainbow reflection effect was generated by grading the PTO stiffness coefficients to prohibit transmission of wave frequency bands at spatially controlled locations within the array, thus producing a wide effective bandgap over the targeted frequency range. The approach to this first stage was based on decomposing the local wave fields into Bloch wave modes, and requires the WEC-resonances to decrease in the direction of the incident wave. Rainbow absorption was then created by adding PTO damping to capture the reflected energy by manipulating the associated complex zeros over the targeted interval. This second stage used analysis of the reflection coefficient in complex-frequency space, and control of the pole–zero pair locations (via the PTO parameters) for zero reflection at discrete frequencies. A rainbow absorbing array of five WECs was presented, which gave high broadband absorption (>98 % over the target frequency interval), although near-perfect absorption was prohibited by small amounts of transmission. The rainbow absorbing theory was used to inform algorithms that minimise the sum of the reflected and transmitted energies to generate near-perfect broadband absorption. An algorithm in which the initialisation incorporates knowledge about the manual process to move complex zeros to the real-frequency axis was found to create a far more efficient optimisation algorithm.
The high-frequency end of the target power capture interval is the most difficult to control as the cluster of pole–zero pairs close to the real frequency axis just above the target frequency interval (related to the second passband) force spikes in transmission and reflection, hence, lowering absorption. Near-perfect absorption at the high-frequency end of the target power capture requires tuning a sufficient number of WECs to high frequencies, although this increases the difficulty of separating the complex zeros close to the real frequency axis without them interfering with one another. Pole–zero pairs had to be forced towards the upper bound by adjusting the constraints and initialisations in the optimisation algorithms. This issue does not seem to have been encountered in previous designs of rainbow absorbing arrays (Jiménez et al. Reference Jiménez, Romero-García, Pagneux and Groby2017; Wilks et al. Reference Wilks, Montiel and Wakes2022), due to near total reflection above the target intervals, and more weakly coupled systems (Romero-García et al. Reference Romero-García, Theocharis, Richoux and Pagneux2016).
The distance between WECs (i.e. the widths of the unit cells) can be used to control the separation of the pole–zero pairs corresponding to the first and second passbands, such that shorter WEC spacing gives greater separation. A spacing similar to the WEC dimensions (
$d/L=0.8$) was chosen to give a reasonable separation between the passbands/pole–zero pairs. A shorter spacing is unlikely to be feasible in applications. From a different perspective, using a longer WEC separation would have the benefit of compressing the second passband and bringing the second bandgap (due to Bragg resonance) to lower frequencies, such that it could be used to extend the target power capture interval to higher frequencies. Thus, the target interval would be covered by the first bandgap (due to local/WEC resonances) and the second bandgap (due to Bragg resonance). The combined bandgap (in the non-absorbing case) is similar to the super bandgap idea being developed in vibroacoustic systems (Guo et al. Reference Guo, Cao, Xiao, Shen and Wen2020; Cleante et al. Reference Cleante, Brennan, Gonçalves and Carneiro2022). However, the second bandgap/pole–zero pairs will occupy the target interval and cause absorption to drop over a frequency interval, such that near-perfect absorption would be more challenging to achieve.
The broadband extraction of wave energy demonstrated by WEC-arrays in this study could help protect coastlines against extreme events (e.g. Ozkan, Mayo & Passeri Reference Ozkan, Mayo and Passeri2022), by countering erosion at vulnerable coastlines as a by-product of power capture (Abanades et al. Reference Abanades, Flor-Blanco, Flor and Iglesias2018), or in the design of WEC-arrays for coastal protection (e.g. Cui et al. Reference Cui2024). In this guise, graded WEC-arrays would provide a broadband extension for coastal protection to the proposed uniform arrays of non-absorbing resonators that reflect (Dupont et al. Reference Dupont, Remy, Kimmoun, Molin, Guenneau and Enoch2017) or attenuate (Zhang et al. Reference Zhang, Jin, Zheng and Xu2024) low-frequency bands of damaging wavelengths, or utilise local and Bragg resonant effects for attenuation over wider frequency ranges (Lorenzo et al. Reference Lorenzo, Pezzutto, De Lillo, Ventrella, De Vita, Bosia and Onorato2023).
The design for broadband absorption can be extended to different modes of motion, devices and operating principles, including multiresonant devices which, as individual devices, are efficient over a relatively broad frequency range for fixed PTO parameters (e.g. Evans & Porter Reference Evans and Porter2012; Crowley, Porter & Evans Reference Crowley, Porter and Evans2013), and oscillating wave surge converters which can achieve high capture factors (e.g. Renzi & Dias Reference Renzi and Dias2012; Noad & Porter Reference Noad and Porter2015; Huang & Porter Reference Huang and Porter2024). In some instances, grading geometry may be more appropriate, as in Wilks et al. (Reference Wilks, Montiel and Wakes2022), which can be applied to oscillating water columns with multiple chambers (e.g. Zhao et al. Reference Zhao, Li, Zhou, Geng, Zou and Qin2023), or by altering flap length in arrays of oscillating wave surge converters (Noad & Porter Reference Noad and Porter2015).
The computational efficiency gained by using a 2-D model facilitated the in-depth analysis of, in particular, complex-frequency space and the optimisation algorithms. However, the key components of the proposed approach to design highly efficient broadband arrays of WECs, namely Bloch waves and complex resonances, also exist in 3-D models. Therefore, we envisage 3-D models of WECs (e.g. Tokić & Yue Reference Tokić and Yue2019, Reference Tokić and Yue2023), arranged into rows (otherwise known as gratings with associated Bloch waves (Bennetts & Squire Reference Bennetts and Squire2009; Peter & Meylan Reference Peter and Meylan2010)), graded such that they produce rainbow reflection (similar to Bennetts, Peter & Craster (Reference Bennetts, Peter and Craster2019)), and with PTOs tuned to achieve rainbow absorption.
Future research directions include incorporating factors such as physical constraints on the WEC operation and motion to develop practical designs for WEC-arrays (Bacelli & Ringwood Reference Bacelli and Ringwood2013; Wang et al. Reference Wang, Engström, Göteman and Isberg2015), and nonlinearities particularly relevant to accurately modelling heaving buoys under operational conditions, for example, through viscous drag and nonlinear Froude–Krylov forces, in partially nonlinear models (Giorgi & Ringwood Reference Giorgi and Ringwood2017; Penalba, Giorgi & Ringwood Reference Penalba, Giorgi and Ringwood2017), or weakly nonlinear approaches (e.g. Michele, Sammarco & d'Errico Reference Michele, Sammarco and d'Errico2018). Based on experimental results relating to band structures (Dupont et al. Reference Dupont, Remy, Kimmoun, Molin, Guenneau and Enoch2017; Lorenzo et al. Reference Lorenzo, Pezzutto, De Lillo, Ventrella, De Vita, Bosia and Onorato2023) and rainbow trapping (Archer et al. Reference Archer, Wolgamot, Orszaghova, Bennetts, Peter and Craster2020) in arrays of non-absorbing resonators, it is likely that the demonstrated rainbow absorption will be robust to nonlinear phenomena, but that the efficiency will be reduced.
In conclusion, we have shown that a handful of heaving point-absorber WECs can achieve near-perfect, broadband absorption, by grading their resonant properties using time-independent PTO parameters. The novel approach for tuning arrays of heaving buoys revolved around manipulating band structures local to each WEC in the array and taking an analytic continuation of the problem into complex-frequency space to connect the near-resonances generated as part of the rainbow reflection phenomenon to pole–zero pairs of the reflection coefficient. Near-perfect, broadband absorption was generated by choosing the PTO damping to move the complex zeros in reflection to the real frequency axis.
Acknowledgements
The authors thank L.J. Yiew for code to solve the single cell problem. A.-R.W. and L.G.B. thank the Isaac Newton Institute for Mathematical Sciences, Cambridge, for support and hospitality during the MWS programme where work on this paper was undertaken (supported by EPSRC grant no EP/R014604/1).
Funding
This work was funded by the Australian Research Council (grant number LP180101109). A.-R.W. is funded by a University of Adelaide PhD scholarship. L.G.B. is funded by an Australian Research Council Future Fellowship (FT190100404). N.S. is the recipient of an Australian Research Council Early Career Industry Fellowship (IE230100545) funded by the Australian Government.
Declaration of interests
The authors report no conflict of interest.
Appendix A. Scattering matrix for the single cell problem
The wave field in cell 
$n$ (subscript omitted for ease of notation) is obtained by solving separate diffraction (
$^d$) and radiation (
$^r$) subproblems, which are combined to return the full wave field 
$\phi = \phi ^d + \xi _h \phi ^r$. The unknown coefficients in the diffraction problem are determined by truncating the eigenfunction expansions at 
$M$ modes and applying kinematic and dynamic boundary conditions at the interfaces between regions, namely
\begin{equation} \phi_i^d = \phi^d_{i+1} \quad \text{and}\quad \frac{\partial \phi^d_{i}}{\partial x} = \frac{\partial \phi^d_{i+1}}{\partial x} \quad \text{at } x ={\pm} L \quad \text{for } i=1,2. \end{equation}Multiplying by an appropriate test function and integrating over depth produces a system of equations for the unknown diffraction coefficients in terms of prescribed incident waves
\begin{equation} \begin{bmatrix} \boldsymbol{a^-}_d \\ \boldsymbol{b^+}_d \\ \boldsymbol{\alpha^+}_d \\ \boldsymbol{\alpha^-}_d \end{bmatrix} = \begin{bmatrix} {{\boldsymbol{\mathsf{R}}}}^{{{\boldsymbol{\mathsf{d}}}}} & {{\boldsymbol{\mathsf{T}}}}^{{{\boldsymbol{\mathsf{d}}}}} \\ {{\boldsymbol{\mathsf{T}}}}^{{{\boldsymbol{\mathsf{d}}}}} & {{\boldsymbol{\mathsf{R}}}}^{{{\boldsymbol{\mathsf{d}}}}} \\ {{{\boldsymbol{\mathsf{G}}}}}^{{{\boldsymbol{\mathsf{d}}}}} & {{{\boldsymbol{\mathsf{H}}}}}^{{{\boldsymbol{\mathsf{d}}}}} \\ {{{\boldsymbol{\mathsf{H}}}}}^{{{\boldsymbol{\mathsf{d}}}}} & {{{\boldsymbol{\mathsf{G}}}}}^{{{\boldsymbol{\mathsf{d}}}}} \end{bmatrix} \begin{bmatrix} \boldsymbol{a^+}\\ \boldsymbol{b^-} \end{bmatrix}. \end{equation}
The vectors 
$\boldsymbol {\alpha ^\pm }_d$ contain the coefficients corresponding to the potential underneath the WEC (Region 2),
\begin{align} \phi_2^d &= \alpha_0^{d-} + \alpha_0^{d+}x + \sum_{m = 1}^{\infty} (\alpha_m^{d+} \exp({{\rm i} \kappa_m (x +L)}) \nonumber\\ &\quad + \alpha_m^{d-} \exp({-{\rm i} \kappa_m (x - L)})) \frac{\cosh(\kappa_m(z+h))}{\cosh(\kappa_m(h-d))}, \end{align}
where 
$\kappa _m = \mathrm {i}m{\rm \pi} /(h-d)$ are purely imaginary wavenumbers. Assuming unit heave amplitude (
$\xi _h = 1$), the radiation potentials in Regions 1 and 3,
$$\begin{gather} \phi^r_1 = \sum_{m =0}^{\infty}(a_m^{r-}\exp({-{\rm i} k_m(x+L)})) \frac{\cosh(k_m(z+h))}{\cosh(k_mh)} \end{gather}$$
$$\begin{gather}\phi^r_3 = \sum_{m =0}^{\infty}(b_m^{r+}\exp({{\rm i} k_n(x-L)})) \frac{\cosh(k_m(z+h))}{\cosh(k_mh)}, \end{gather}$$
are obtained analogously to the diffraction potential. Underneath the WEC, 
$\phi ^r_2$ is composed of a homogeneous and particular solution which satisfies the boundary condition 
${\partial _z \phi } = { \xi _h\omega ^2/g}$,
\begin{align} \phi_2^r &= \alpha_0^{r-} + \alpha_0^{r+}x + \sum_{m = 1}^{\infty} (\alpha_m^{r+} \exp({{\rm i} \kappa_m (x +L)}) \nonumber\\ &\quad + \alpha_m^{r-} \exp({-{\rm i} \kappa_m (x - L)}))\zeta_m + \frac{\omega^2}{2g(h-d)}[(z+h)^2-x^2]. \end{align} To determine the unknown heave amplitude (2.10), the hydrodynamic forcing is found by integrating the pressure over the wetted surface of the body 
$S_B$, with respect to the outward pointing unit normal to the surface of the body,
\begin{equation} \rho g \int_{S_B} \phi n \ \mathrm{d}S = \underbrace{\rho g \int_{S_B} \phi^d n \ \mathrm{d}S}_{{excitation\ force}} + \underbrace{\rho g \int_{S_B} \xi_h \phi^r n \ \mathrm{d}S}_{{radiation\ component}}. \end{equation}The radiation component resulting from forced oscillations of the body is separated into real and imaginary parts in phase with the acceleration and velocity of the body, respectively, to give the added mass and radiation damping (Linton & McIver Reference Linton and McIver2001; Mei et al. Reference Mei, Stiassnie and Yue2005)
\begin{equation} [\omega^2 a(\omega) + \mathrm{i}{\omega} b(\omega)]\xi_h = \rho g \left(\int_{S_B} \xi_h \phi^r n \ \mathrm{d}S\right). \end{equation}
The excitation force 
$F$ produces a vertical force at the underside of the WEC which depends on the coefficients of the incident waves, and therefore, interactions within the array. This dependence is incorporated in the initial scattering matrix of a WEC by expressing the excitation force in the radiation problem in terms of the incident amplitudes,
\begin{align} \int_{S_B} \phi^d n \ \mathrm{d}S &= \boldsymbol{V^+}\boldsymbol{\alpha^{+}_d} + \boldsymbol{V^-}\boldsymbol{\alpha^{-}_d} \end{align}
\begin{align} &= [\boldsymbol{V^+}{{{\boldsymbol{\mathsf{G}}}}}^{{{\boldsymbol{\mathsf{d}}}}} + \boldsymbol{V^-}{{{\boldsymbol{\mathsf{H}}}}}^{{{\boldsymbol{\mathsf{d}}}}}] \boldsymbol{a^+} +[\boldsymbol{V^+}{{{\boldsymbol{\mathsf{H}}}}}^{{{\boldsymbol{\mathsf{d}}}}} + \boldsymbol{V^-}{{{\boldsymbol{\mathsf{G}}}}}^{{{\boldsymbol{\mathsf{d}}}}}] \boldsymbol{b^-}, \end{align}
where 
$\boldsymbol {V^\pm }$ are vectors with 
$M+1$ elements. From (2.10), the heave amplitude is
\begin{align} \xi^h(\omega) &= [-\omega^2(\mathcal{M}+a(\omega))-\mathrm{i}\omega (b(\omega)+b^{PTO}) + (c+c^{PTO})]^{{-}1} F(\omega) \end{align}
\begin{align} &= \boldsymbol{\hat{V}^+}\boldsymbol{a^+} + \boldsymbol{\hat{V}^-}\boldsymbol{b^-}. \end{align}
The scattering matrix 
${{\boldsymbol{\mathsf{S}}}}$ in (3.3) for the heaving buoy is
\begin{equation} {{\boldsymbol{\mathsf{S}}}} = {{\boldsymbol{\mathsf{S}}}}^{{{\boldsymbol{\mathsf{d}}}}}+ {{\boldsymbol{\mathsf{S}}}}^{{{\boldsymbol{\mathsf{r}}}}} = \begin{bmatrix} {{\boldsymbol{\mathsf{R}}}}^{{{\boldsymbol{\mathsf{d}}}}} + \boldsymbol{a^-}_r \boldsymbol{\hat{V}^+} & {{\boldsymbol{\mathsf{T}}}}^{{{\boldsymbol{\mathsf{d}}}}} + \boldsymbol{a^-}_r \boldsymbol{\hat{V}^-} \\ {{\boldsymbol{\mathsf{T}}}}^{{{\boldsymbol{\mathsf{d}}}}} + \boldsymbol{b^+}_r \boldsymbol{\hat{V}^+} & {{\boldsymbol{\mathsf{R}}}}^{{{\boldsymbol{\mathsf{d}}}}} + \boldsymbol{b^+}_r \boldsymbol{\hat{V}^-} \end{bmatrix}, \end{equation}
with the corresponding coefficients in (3.1) given by 
$a_m^- = a_m^{d-}+ \xi _h a_m^{r-}$ and 
$b_m^+ = b_m^{d+}+ \xi _h b_m^{r+}$. The multiple scattering routine resolves the coupling between WECs, and the forces on a particular WEC can be recovered from (A.10).
Appendix B. Convergence results
The truncation of 
$M = 25$ was determined based on the convergence of the WEC amplitudes and scattering coefficients in the single cell problem. Higher truncation values produce indistinguishable solutions (converged to four decimal places), as illustrated in figure 15 where the absorption for the array of 
$N=5$ WECs in table 2 with 
$M = 25$ evanescent modes is compared with the absorption when 
$M = 100$ evanescent modes are included. Here, 
$N$ refers to the number of evanescent modes included in WEC interactions, with 
$N = 0$ indicating the wide-spacing approximation is applied to the array. The solution obtained when including evanescent modes (
$N>0$) in WEC interactions is overlaid to further justify the use of the wide-spacing approximation. The average absorption of the array increases from 
$\hat {\alpha } = 0.9900$ to 
$\hat {\alpha } = 0.9901$ with the inclusion of evanescent modes. Only 
$N = 3$ evanescent modes are required for convergence (solution is unchanged for larger truncation values, e.g. 
$N = 10$ with 
$M = 100$ in figure 15). The error incurred was deemed to be acceptably small so as to motivate the reduced computational effort associated with the wide-spacing approximation.
The absorption of the graded array of 
$N=5$ WECs in table 2 is shown for combinations of the truncation 
$M$ in the single body problem, and the number of evanescent modes 
$N$ included in WEC interactions. While sufficient evanescent modes are included in the single body problem for convergence of the final solution, the wide-spacing approximation (
$N = 0$) incurs a small error in the average absorption of 0.0001.
Open access
