Rainbow reflection and broadband energy absorption of water waves by graded arrays of vertical barriers

Abstract The rainbow reflection effect describes the broadband spatial separation of wave spectral components caused by a spatially graded array of resonators. Although mainly studied in optics and acoustics, this phenomenon has recently been demonstrated both theoretically and experimentally for water waves travelling through an array of vertical cylinders. The scattering of two-dimensional linear water waves by an array of vertical, surface-piercing barriers is considered here, in which both the submergence and spacing between the barriers are spatially graded. The rainbow reflection effect arises naturally as wave energy temporarily becomes amplified at different locations depending on frequency. Band diagram calculations are used to demonstrate that this is a consequence of the wave gradually slowing down throughout the array. The wave/barriers scattering problem is then augmented by positioning heave-restricted, rectangular floating bodies equipped with a linear damping mechanism between each adjacent pair of barriers. A solution to the resulting boundary value problem is obtained using an integral equation/Galerkin method. Using constrained optimisation, passive rainbow absorbers are designed that achieve near-perfect absorption (i) over a discrete set of frequencies and (ii) over an octave. This suggests potential applications of rainbow absorbers in the design of smart coastal technologies.


Introduction
There is growing interest throughout the wave sciences in metamaterials which have a subwavelength structure that is spatially varying in some way.This spatial variation can give rise to so called rainbow phenomena, through which broadband wave signals are slowed and spatially separated into their spectral components.These phenomena can be classified into rainbow trapping, where the separated energy becomes permanently localised, and rainbow reflection, where the energy is temporarily amplified before being reflected [He et al., 2012, Chaplain et al., 2020].First demonstrated by Tsakmakidis et al. [2007] for optical waves, these concepts have since inspired the development of analogues in several fields, including elastic waves [Skelton et al., 2018, Arreola-Lucas et al., 2019], seismic waves [Colombi et al., 2016], water waves [Bennetts et al., 2018] and acoustics.In this latter discipline, rainbow reflection has been prominent in the design of acoustic filters and sensors, where one-dimensional waveguides with graded resonant side-cavities have been shown to spatially separate frequency components [Zhu et al., 2013, Ni et al., 2014, Zhao and Zhou, 2019].The local energy amplification of acoustic waves has also been studied in two-dimensional arrays of (i) rigid cylinders with chirped spacing [Romero-García et al., 2013, Cebrecos et al., 2014] and (ii) C-shaped resonators with graded radii [Bennetts et al., 2019].Jiménez et al. [2017] further extended the study of rainbow reflecting devices by including a loss-inducing mechanism within each resonator of a graded array.Using numerical optimisation, they created a structure that achieves near-perfect sound absorption over a prescribed frequency interval.Our primary goal is to investigate the feasibility of achieving this rainbow absorption for water waves.
Recently, several other remarkable capabilities of metamaterials, originally studied for optical or acoustic waves, have been predicted in the context of water waves, where research is motivated by applications that include the development of coastal protection technology or wave energy parks.For instance, two-dimensional periodic arrays of uniform rigid cylinders immersed in water exhibit negative refraction, ultra-refraction and bandgaps-the latter characterising frequency intervals in which waves cannot propagate without a change in amplitude [McIver, 2000, Hu et al., 2004, Farhat et al., 2010].The lowest frequency bandgap of such arrays is shifted to lower frequencies when the cylinders are replaced with C-shaped cylindrical resonators, which suggests that they could be used as effective breakwaters [Hu et al., 2011, Dupont et al., 2017].Bennetts et al. [2018] demonstrated rainbow reflection in arrays of C-shaped cylindrical resonators with graded radii.In these arrays, the group velocity of the waves gradually approaches zero through the array at a location that depends on frequency, which results in local energy amplification.Arrays of rigid cylinders with chirped spacing also give rise to the rainbow reflection effect and associated local energy amplifications, which was verified experimentally by Archer et al. [2020].Here, the absence of a resonant substructure means that this effect is limited to higher frequency intervals when compared with arrays of C-shaped cylinders.Although local energy amplification in graded and periodic arrays is now fairly well understood, energy absorption within such systems remains mostly unexplored.
In the last few decades, wave energy converters have been the subject of an intensive research effort [Falcão, 2010].Several devices have been proposed which are capable of absorbing 100% of wave energy at a particular frequency, but their efficiency is usually sensitive to a small change in frequency [Salter, 1974, Evans et al., 1979, Evans and Porter, 1995].Because real ocean waves occur over a broad range of frequencies, the limited band of effectiveness of single absorbers is typically overcome using optimal control theory [see Coe et al., 2017 and references therein].Recent research has also focused on the design of arrays of multiple absorbers, or wave energy parks, which optimally extract energy from realistic sea states [Babarit, 2013, Göteman et al., 2020].By using shallow water theory, Porter [2021] recently demonstrated perfect absorption at all frequencies by a semi-infinite array of damped buoys.The buoys were modelled using a homogenised boundary condition at the upper surface of the two-dimensional fluid domain.To achieve perfect absorption, the first segment of the absorbing surface must have a negative damping coefficient and thus add energy to the fluid.Using full linear water wave theory, a similar device of finite length was shown to achieve near perfect absorption for sufficiently large frequencies.
In §2, we explore rainbow reflection in a one-dimensional waveguide consisting of a graded array of parallel, surfacepiercing rigid vertical barriers.The limited complexity of this problem compared with two-dimensional arrays allows us to reveal a clear relationship between the structural geometry and the local energy amplifications.We build upon the results of Ursell [1947] and Porter and Evans [1995], who studied the scattering of water waves by a single vertical barrier, and those of Evans and Morris [1972], Newman [1974] and McIver [1985], who studied resonators consisting of two vertical barriers.In the present study, we demonstrate how these resonators can be combined into an array, and how local energy amplification within this array relates back to the local geometry, as well as to the band structure of the corresponding periodic infinite array.
The local energy amplification in graded arrays of vertical barriers suggests that highly efficient broadband absorption could be achieved by positioning suitable energy absorbing mechanisms within the array.In §3 we modify the problem so that the surface between each successive pair of barriers is occupied by a damped rectangular piston that is restricted to move in heave.Our solution to the corresponding linear boundary value problem is obtained by reformulating it as a coupled system of integral equations, which we solve using a Galerkin method.In §4, we briefly discuss absorption by a single piston oscillating between two barriers.Then, we design rainbow absorbers which achieve near-perfect absorption over (i) a discrete set of frequencies, and (ii) over an octave.Our results extend those from acoustics obtained by Jiménez et al. [2017], and highlight the importance of non-absorbing structures in transferring and localising energy.

Preliminaries
We consider the scattering of linear water waves by N + 1 parallel surface-piercing vertical barriers in a fluid domain of infinite horizontal extent and constant, finite depth H.A two-dimensional Cartesian coordinate system (x, z) is used, where the z-axis points vertically upwards and the horizontal coordinate x corresponds to the direction of propagation of the waves.The fluid domain is bounded vertically by a flat sea-bed at z = −H and a free surface located at z = 0 when at equilibrium.Each vertical barrier is of infinitesimal thickness and is situated at x = x (n) , where x (n) ∈ R for all n ∈ {0, . . ., N }, and x (n) > x (n−1) for n ≥ 1.Each barrier is defined by the set of points Γ is the submergence depth of the n th barrier.A schematic of this geometry is given in figure 1.
In order to apply linear water wave theory to this problem, we assume that the steepness of the free surface elevation of the incident waves is small.Thus, we consider the fluid to be inviscid, incompressible and undergoing irrotational, timeharmonic motion with angular frequency ω.We assume that there is no normal flow across the seabed and the barriers.
The velocity field of the fluid can therefore be expressed as the spatial gradient of some potential The complex-valued spatial function φ then satisfies the following boundary value problem where (3c) is the linear free-surface condition with g denoting acceleration due to gravity.
Since φ is twice differentiable, it must be continuous and have a continuous first partial derivative with respect to x across the gaps beneath each barrier.This can be expressed in terms of one-sided limits as lim for all n = 0, . . ., N where z ∈ (−H, −d (n) ).We require that |φ| remains finite as x → ±∞.The solution must also obey an appropriate Sommerfeld radiation condition, which enforces that the scattered plane waves must propagate away from the array of barriers (see ( 10)).
The submerged edges of the barriers located at (x (n) , −d) create singularities of order −1/2 in the fluid velocity [Evans andMorris, 1972, Mosig, 2018].Thus we require that for all n ∈ {0, . . ., N } The general solution to the (3a)-(3c) is obtained using separation of variables The quantities k m are the solutions to the free surface dispersion relation k tanh kH = ω 2 /g.We use the convention that k 0 = 2π/λ ∈ R + , where λ is the wavelength of the propagating waves, and that −ik n ∈ ((n − 1)π/H, nπ/H) for all n ∈ N. The vertical eigenfunctions are defined as where the superscript (0) has been included for consistency with §3.The normalisation constants β m = sinh(2k m H)/(4k m H) + 1 2 are chosen so that the vertical eigenfunctions satisfy for all m, p ∈ {0} ∪ N, where δ mp is the Kronecker delta.
The fluid motion is forced by incident plane waves with angular frequency ω.This forcing determines the incident potential on the exterior regions of the fluid, which is of the form We will mostly consider problems where the forcing originates from the left.In this case, B = 0 for all m ≥ 1.The Sommerfeld radiation condition is thus for all z ∈ (−H, 0).

Overview of the solution technique
We invoke the continuity conditions (4a,b) and the no flow condition (3d) for each barrier in order to map the amplitudes of waves incident on the nth barrier to the amplitudes of the waves scattered by that barrier.This mapping is approximated by the scattering matrix relation for n ∈ {0, . . ., N }.Here, A (n) and B (n) are vectors with entries A (n) m and B (n) m , for 0 ≤ m ≤ M , where M is a truncation parameter.
The entries of each scattering matrix S (n,n+1) are determined as a function of the barrier submergence d (n) as follows: (i) express the scattered wave amplitudes in terms of the horizontal fluid velocity u(z) = ∂ x φ on the interval (−H, −d (n) ), (ii) invoke the continuity conditions to construct an integral equation satisfied by u, and (iii) solve the integral equation via a Galerkin method which involves expanding u over a basis of weighted Chebyshev polynomials, which encapsulate the singularities at (x (n) , −d (n) ) [see Porter and Evans, 1995 for more details].
The global scattering matrix is used to compute the amplitudes of the scattered waves in Ω (0) and Ω (N +1) .Finally, the wave amplitudes in the interior regions Ω (n) are computed recursively for n ∈ {1, . . ., N }.The algorithm used to construct the global scattering matrix and compute the unknown coefficients is based on the method developed by Ko and Sambles [1988] to study electromagnetic wave propagation in layered media.This method has been adapted for problems involving the multiple scattering of linear water waves, yielding numerically stable results [Bennetts andSquire, 2009, Montiel et al., 2015].

The fundamental resonance of a pair of barriers
We first consider the response of an isolated resonator consisting of a pair of barriers, i.e. when N = 1.This problem has been extensively studied when d (0) = d (1) .Such a pair of identical barriers is associated with an infinite collection of resonant frequencies [Evans and Morris, 1972].Here, we are exclusively concerned with the fundamental resonance (i.e. the lowest frequency resonance), which is analogous to the Helmholtz resonance in acoustic resonators.The corresponding fundamental frequency, which we denote ω fund , can be approximated by g/d (0) .This approximation is valid when the barriers are close to each other and submerged in deep water, i.e. k 0 l (1) << 1 and d (0) /H << 1.Under these assumptions, the mass of fluid between the barriers can be treated as a heaving rectangular body oscillating about its hydrostatic equilibrium [Newman, 1974].As l (1) increases, this approximation loses validity and the fundamental frequency decreases [McIver, 1985].When d (0) = d (1) , the fundamental frequency of the corresponding resonator is similar to that of a pair of identical barriers with submergence (d (0) + d (1) )/2, provided |d (1) /d (0) − 1| << 1.This result is shown in figure 2a.Throughout this paper, we fix the depth of the fluid is to be H = 20 m.
The bandwidth of the fundamental resonance of a pair of identical barriers is governed by the ratio l (1) /d (0) .To quantify this relationship, we first introduce the following analogue for the average power of the oscillations between two barriers where φ is the velocity potential of the fluid when excited by a plane wave of amplitude 0.1 m incident from the left.We define ∆ω fund to be the full width at half the maximum of the peak of P associated with the fundamental resonance.
The Q-factor is then Q = ω fund /∆ω fund .As shown in 2b, Q and l (1) /d (0) are inversely proportional.This extends the assertion of McIver [1985] that the fundamental resonance 'detunes' as the spacing between the barriers increases, which was quantified by the distance between the zeros of the reflection and transmission coefficients associated with the resonance.To the best of our knowledge, the inverse relationship between Q and l (1) /d (0) has not been reported before.

Local energy amplifications in graded arrays
The observations from the previous section motivate the design of a rainbow trapping device consisting of multiple vertical barriers, where the submergence and spacing of the barriers increases throughout the array.We aim to adjust the parameters d (n) and l (n) so that the fundamental resonances of each cavity are evenly spaced over the octave ω ∈ [0.8, 1.6] s −1 .The resonant frequencies of the cavities must decrease in the direction of wave propagation, in order to prevent the deep barriers of the low frequency cavities from sheltering the high frequency cavities.This is because the transmission coefficient of the nth barrier in isolation is a decreasing function of ω 2 d (n) [Ursell, 1947].We design the array so as to excite local energy amplifications when the waves are incident from the left.Thus for N = 9 the fundamental resonance of the nth cavity is given by ω  0) .There is a very strong agreement between the data used to generate this plot and the model 1/Q = a l (1) /d (0) , where a ∈ R is a coefficient.This demonstrates that Q and l (1) /d (0) are inversely proportional.In both subfigures we use d (0) = 10 m, and the depth of the fluid is fixed to be H = 20 m here and throughout this paper.1: The submergences d(n) of pairs of identical barriers whose fundamental resonance occurs at ω (n) , which were obtained by maximising the time-average power of free surface oscillations.The optimised submergence values occur at approximately 0.9 times the estimate g/(ω (n) ) 2 .
First, we find the submergence d(n) for pairs of identical barriers so that the fundamental resonance is ω (n) , while fixing the separation of the barriers to be l(n) = d(n) /5.This is done by maximising P (ω (n) ) as a function of d(n) ∈ (0, H), which yields the optimised submergences of the barriers presented in table 1.Second, we use these optima to choose parameters for the array of barriers.We want the geometry of the nth cavity to be similar to the pair of identical barriers which produces a resonance at ω (n) .Thus, we select l (n) = l(n) .In addition, we require that the average submergence of the two barriers forming each cavity is equal to the depth of the optimised pair of barriers, that is for all n ∈ {1, . . ., N }.This is an underdetermined system for the N + 1 unknown parameters d (n) .We choose the solution which satisfies the constraints 0 < d (1) < • • • < d (N ) < H and minimises the grading of the array (see appendix A for further details).
Colour plots of |φ(x, z)| for the array parametrised by the submergences reported in table 1 are given in figure 3 at all the frequencies ω (n) .As expected, the strongest local energy amplification in the graded array for a given frequency occurs in the cavity associated with that fundamental frequency.In results not shown here, this does not occur when waves are incident from the right.This is a consequence of the low transmission of high frequency waves discussed earlier.
In order to analyse this phenomenon further, we consider the propagation of water waves through an infinite periodic array of uniform vertical barriers.Using the present notation, this is the case where x (n) = nl and d (n) = d for all n ∈ Z, where l > 0 and d ∈ (0, H) are fixed.Our intention is to locally approximate the wave field within a finite graded array of vertical barriers using solutions to the infinitely periodic case, which is reasonable provided that the grading is weak [Bennetts et al., 2019].Infinite periodic arrays of uniform vertical barriers can be studied using Bloch's theorem, which motivates seeking solutions of the form In the above equation we have introduced the unknown Bloch wavenumber q ∈ C. To compute the corresponding Bloch modes, ( 13) is applied to the horizontal boundary of the unit cell {(x, z)|x ∈ (−l/2, l/2) and z ∈ (−H, 0)}, which yields These quasi-periodic boundary conditions are then solved in conjunction with (3a-3d), where suitable modifications are made so that x ∈ (−l/2, l/2).
A numerical solution is obtained by solving a generalised eigenvalue problem, which incorporates the scattering matrix for a single vertical barrier [see e.g.Peter and Meylan, 2009].For each frequency ω, propagating Bloch modes (i.e.q ∈ R) are sought in the irreducible Brillouin zone q ∈ [0, π/l].We obtain dispersion curves in the q − ω plane, which form so-called band diagrams.We remark that frequency intervals occupied by the dispersion curves are called passbands, and these frequencies are associated with a propagating mode which can travel through the infinite array without a change in amplitude.The slope of the dispersion curve dω/dq determines the group velocity of this Here, the values of d (n) and l (n) are those of the graded array shown in figure 3.In each subfigure, the dotted line is the dispersion curve for unobstructed waves.Notably, the bandgap begins above the fundamental frequency of a pair of barriers with separation l and submergence d.These fundamental frequencies, which are represented by dashed lines, are (a) propagating mode.The frequency intervals in which there are no real-valued Bloch wavenumbers are called bandgaps.At these frequencies, waves cannot propagate without a change in amplitude.Some examples of band diagrams with different barrier submergences and spacings are shown in figure 4. In each case, there is only one passband and one bandgap.We observe that the fundamental frequency of a pair of barriers with separation l and submergence d is a lower bound of the bandgap.The fundamental frequency decreases as the barriers become deeper, which coincides with the dispersion curves becoming more compressed.This mirrors a similar result found by Bennetts et al. [2018] for infinite two-dimensional arrays of C-shaped cylindrical resonators.In that setting, the resonant frequency of the isolated cylinder is an upper bound of the passband.In either cases, the bandgap can be shifted by modifying the resonant substructure of the array.
In the finite array of barriers, the potential surrounding each successive pair of barriers can be understood in terms of the corresponding Bloch problem.In particular, the dispersion curve of the infinite array of barriers with submergence (d (n−1) + d (n) )/2 and spacing l (n) approximately determines the behavior in the nth cavity for a weakly graded structure, and in this setting it is referred to as the local dispersion curve.Since the fundamental frequencies of the cavities decrease from left to right, the local dispersion curves become more compressed and the frequency interval which supports propagating modes becomes narrower.Thus low frequency waves propagate further towards the right than high frequency waves.
Finite difference estimates of the group velocity of the propagating mode are tabulated in table 2 for all cavities and for a range of frequencies.Waves slow down as they propagate through the finite, graded array because the group velocity is a decreasing function of n, which attains a minimum before becoming undefined.This minimum typically occurs in the cavity whose fundamental frequency is the frequency of the wave.The wave can only attenuate beyond this point [which is called the turning point Romero-García et al., 2013], since propagation at this frequency is no longer supported by the array.A local energy amplification occurs near the turning point because the propagating mode has a relatively small group velocity, so its energy is localised for a relatively large duration.Moreover, the group velocity is a decreasing function of ω where it is defined, which means that higher frequency waves slow down more rapidly than lower frequency waves and amplify further towards the left, resulting in spatial separation.Thus, graded array of vertical barriers gives rise to the rainbow reflection effect.
However, the graded array does not give rise to the rainbow trapping effect.Because the unit cell of the Bloch problem possesses axial symmetry, the infinite array supports both left and right propagating modes at passband frequencies, which have equal and opposite group velocity.Consequently, the strong coupling between the right and left propagating modes near the turning point results in the wave changing direction [He et al., 2012, Chaplain et al., 2020].The left-travelling wave then accelerates as it propagates in reverse, carrying energy that is ultimately reflected by the array.In the following section, we will modify the system so that the localised energy can be absorbed within the array rather than being reflected.2: Finite difference estimates of the group velocity (m s −1 ) of waves travelling through the infinite array of barriers with d = d (n−1) and l = l (n) , with frequency ω.Dashes (-) are used to denote the absence of a propagating Bloch mode.The group velocity is a decreasing function of both ω and n. 3 Solution of multiple piston cavity problem

Preliminaries
We adapt the problem considered in §2 by introducing N rigid rectangular floating bodies, or 'pistons', which completely occupy the horizontal gap between each adjacent pair of barriers.The pistons are restricted to move vertically in heave and are less dense than the fluid.Further, the density of each piston is ρ A is the area density of the piston with respect to the wetted surface, and δ (n) is the thickness.The motion of each piston is governed by a damping term proportional to the piston's vertical velocity, with the damping coefficient denoted µ (n) .A schematic of the geometry of this problem is given in figure 5.
We assume that pistons are undergoing small, time-harmonic oscillations about their equilibria.This implies that the submergence of the bottom face of each piston is where s (n) and A /ρ are the complex amplitude and the equilibrium submergence of the n th piston respectively.Here, ρ is the density of the fluid.The relevant forces which govern the dynamics of the pistons are due to hydrostatic pressure, hydrodynamic pressure, gravity and viscous damping.After some algebra, the linearised equations of motion can be stated as where , 1971].
We redefine the subregions Ω (n) to be the sets of points beneath the n th piston when at equilibrium, that is )} for all n ∈ {1, . . ., N }.The definitions of Ω (0) and Ω (N +1) are unchanged.The velocity potential of the fluid φ satisfies the boundary value problem given by (3a), (3b), (3d) and the additional conditions The potential is also subject to (4a,b), ( 5), ( 9) and (10).

Derivation of a system of integral equations
The expansion for the solution given in (6) still applies in subregions Ω (0) and Ω (N +1) , since there is no piston at the surface of the fluid.In the remaining subregions Ω (n) , we express the general solution as the sum of a homogeneous solution φ (n) h and a particular solution φ p , for all n ∈ {1, . . ., N }.The homogeneous solution satisfies (3a), (3b) and for x ∈ (x (n−1) , x (n) ).These three equations constitute a Sturm-Liouville problem in the z-direction.By using separation of variables and choosing an expansion centered at the midpoints of each piston m where κ n) , and the vertical eigenfunctions have been defined as for all m ∈ N.These eigenfunctions satisfy the following orthogonality relation for all m, p ∈ {0} ∪ N and for all n ∈ {1, . . ., N }.
The particular solution φ (n) p satisfies (3a), (3b) and (17b).We find a particular solution with the ansatz that φ p is quadratic in both x and z, which gives where we have imposed symmetry about the horizontal midpoints x = m (n) .This treatment of the inhomogeneous boundary condition is in complete analogy to the solution obtained by Yeung [1981] for the problem of an oscillating vertical circular cylinder in a fluid of finite depth.

Numerical solution using a Galerkin method
To find an approximate solution to the system of integral equations, we first suppose that for each n = 0, . . ., N , the auxiliary function u (n) can be expanded in terms of an appropriately chosen basis {v n) ).We then approximate each auxiliary function u (n) by where we have truncated each expansion at some sufficiently large M aux ∈ N.
After substitution of (32) into the system ( 29)-( 31), we multiply each expression by v (n) p (z), for some p ∈ {1, . . ., M aux }.Here, n ∈ {0, . . ., N } is the index of the barrier associated with each integral equation.For instance, (29) was derived by applying (4a) at x (0) , so this expression is multiplied by v (0) p (z) following substitution of (32).Integrating each of these expressions with respect to z over (−H, −d (n) ) subsequently yields the following linear system of equations where ( 34) holds for all n ∈ {1, . . ., N − 1} and we have defined Additionally, we have defined The singularities at (x (n) , −d (n) ) motivate the following choice of basis functions in which T n is the Chebyshev polynomial of the first kind of degree n.This basis was first introduced by Porter and Evans [1995] to solve the scattering of linear waves by a single surface-piercing barrier.By using equation (1.10.2) in Bateman [1954], the integrals in (36a) and (36b) can be evaluated as )), where J n is the Bessel function of the first kind of order n, and I n is the corresponding modified Bessel function.Furthermore, by applying the symmetry and orthogonality of the even Chebyshev polynomials, the integrals in (36c) can be evaluated as Equations ( 33), ( 34) and ( 35) form a linear system of (N + 1)M aux equations with (N + 1)M aux + 3N unknowns.To obtain a further 2N equations, we once again apply (4b) as x → x (n−1)+ and x → x (n)− .However, in contrast to the derivation of ( 27) and ( 28), we use the fact that ψ m is orthogonal to the constant function 1 for all m ∈ N.After referring to (36c), we have for n ∈ {1, . . ., N }.
The remaining N equations are derived from the equations of motion of the pistons ( 16).The integral of the hydrodynamic pressure is evaluated by invoking symmetry about x = m (n) , then expressed in terms of the coefficients of the auxiliary functions by referring to ( 27), ( 32) and (36b) After some rearrangement, substitution of ( 15) and ( 41) into (16) yields where we have defined Finally, we truncate the infinite sums in ( 33), ( 34), ( 35) after M ker ∼ 1000 terms, then combine them with (39), ( 40) and ( 43) to obtain a matrix equation of dimension M aux (N + 1) + 3N .This allows us to simultaneously solve for 3.41 3.88 4.46 5.17 6.06 7.22 8.74 10.78 13.58 μ(n) (kg m −1 s −1 ) 107 130 159 199 253 332 460 692 1170 Table 3: The optimised parameters of a SCA which maximise absorption at ω (n) , i.e. α(ω (n) ) = 0.5 is attained.the coefficients of the auxiliary functions c (n) j for all n ∈ {0, . . ., N } and j ∈ {1, . . ., M aux }, as well as the quantities and s (n) for all n ∈ {1, . . ., N }.The size of M aux must be chosen in accordance with M ker ∼ 1000.By choosing M ker = 1000 and M aux = 50, we obtain 4-digit accuracy in the unknown coefficients R, T and s (n) .The convergence properties of this class of Galerkin method is discussed in more detail by Linton and McIver [2001].

Energy absorption
Due to the presence of damping, the class of devices introduced in §3 is capable of absorbing energy from the incident waves.Here, absorption is defined as α(ω) = 1 − |T | 2 − |R| 2 , where the transmission and reflection coefficients T and R are calculated when the system is forced from the left by a plane wave of frequency ω.By considering the energy flux of the incident and scattered plane waves [Linton and McIver, 2001] and the time-average work done by the damping force, we can write This second form allows us to consider the contribution of each piston to the absorption by the entire device.It is clear that the minimum value of absorption is α = 0, which occurs, for example, in the absence of any damping.Further, the maximum value of absorption is α = 1, which occurs when no energy is scattered by the device (i.e.R = T = 0).

Absorption by a single cavity
We first consider absorption by a single cavity absorber (SCA), here defined as a single piston surrounded by a pair of barriers.When the two barriers are identical, an established theoretical result implies that the maximum value of absorption is α = 0.5 [Budal and Falnes, 1975, Evans, 1976, Newman, 1977, Mei, 1976].For any given frequency, this maximum can be attained by tuning the damping coefficient of the piston and the submergence of the barriers.This tuning is simplified by fixing l (1) = d (0) /5 and ρ (1) A = 500 kg m −2 .The parameter which maximise absorption for the frequencies considered in §2.3 are presented in table 3, where we use a tilde to distinguish parameters that relate to SCAs.We observe that the submergence depths are nearly identical to those in table 1, which were obtained by maximising the time-average power of the free surface oscillations.This is because the fundamental resonance considered in §2.3 is essentially a vertical fluid motion between the barriers, which is the fluid motion that maximises energy absorption by a heaving piston.
We remark that larger values of α can be attained when d (1) > d (0) .Perfect absorption (i.e.α = 1) occurs when the wave radiated by the piston destructively interferes with the scattered wave, which is impossible for a symmetric  4: Parameters of the device obtained using the modified algorithm of Jiménez et al. [2017].The absorption spectrum of this device is shown in figure 6i.device oscillating in one mode of motion [see Falnes, 1998 for more details].However, in the limiting case where d (1) = H, the scattered wave propagates in one direction only as no transmission is allowed.This implies that complete destructive interference of the scattered wave by the radiated wave produced by a heave-restricted piston is possible.This limiting case is analogous to the oscillating water columns positioned against a wall considered by Evans and Porter [1995], where perfect absorption at a given frequency was achieved.

Optimised absorption at N frequencies
Consider now an array of absorbing cavities whose absorption spectrum has maxima at the N = 9 prescribed frequencies considered in §4.1.We construct this array by modifying the algorithm which was used by Jiménez et al. [2017] to design an acoustic rainbow absorber.The algorithm is a sequence of optimisation problems, whereby the size of the device and the set of frequencies at which absorption is maximised both increase with each iteration.The details of the algorithm are provided in the following paragraph.
For n iteratively incrementing from 1 to N , we maximise n j=1 α(ω (N −j+1) ) for a device consisting of n cavities, by tuning d (j) n for j ∈ {0, . . ., n}, and µ (j) n for j ∈ {1, . . ., n}.Here, the subscript n refers to the stage of the iteration.The optimisation is performed while fixing l n )/10 for j ∈ {1, . . ., n}, and with the constraints The absorption spectra of all of the intermediate devices obtained by this process are given in figure 6, and the parameters of the final device with n = 9 cavities are given in table 4. As desired, we find that the peaks of the spectra of all intermediate devices coincide with the prescribed frequencies.In analysis not presented here, we find that at the jth peak of each spectrum, the largest contribution to absorption is by the (n − j + 1)th piston.For the final device, we obtain α(ω (j) ) > 0.98 for j ∈ {1, . . ., N − 1}.The constraint that d (N ) < 0.8H (which was imposed to ensure the solution converges rapidly) hinders absorption at the lowest frequency, where we obtain α(ω (N ) ) ≈ 0.93.In contrast, Jiménez et al. [2017] was able to obtain an acoustic device for which α = 1 at all N frequencies as the result of using a less-restrictive absorption mechanism.We also observe that the peaks of the spectrum in figure 6i become increasingly sharp as ω increases.As a consequence, absorption becomes increasingly lower in the intervals between the prescribed frequencies.Our next task is to rectify this by using a broadband absorption metric to design a rainbow absorber (RA), whose absorption spectrum remains close to 1 over the entire interval [ω (N ) , ω (1) ].

Optimised broadband absorption
Using a similar metric to Jiménez et al. [2017], we quantify the broadband performance of the RA as the average value of α over the interval [ω (N ) , ω (1) ], i.e.
It is clear that E ∈ [0, 1].Naively, we can increase E by increasing the length of the cavities, since this decreases the Q-factors of the resonant peaks and allows them to overlap more.Broadband performance can also be increased by adding more appropriately-tuned cavities to the RA (i.e.increasing N ), since this increases the number of absorption peaks in a given interval.Both of these avenues to increasing E also increase the total length of the RA given by L = N n=1 l (n) , which is undesirable due to practical considerations.This motivates the present problem of maximising E while constraining L, where we will fix N = 9. : Absorption spectra of the devices obtained using the steps outlined in §4.2, which is a modification of that used by Jiménez et al. [2017].The frequencies which absorption is maximised at each stage of the algorithm are denoted by vertical dotted lines.As desired, the peaks of the spectra coincide with the desired peaks of absorption.To achieve this, we simultaneously tune the parameters d (n) , µ (n) and l (n) in order to maximise E, subject to the constraint that L < λ 0.8 /5, where λ 0.8 ≈ 86.4 m is the wavelength of a plane wave with frequency ω (N ) = 0.8 s −1 .We also require that l (n) and µ (n) are increasing, and that 2D (1) < d (0) < • • • < d (N ) < 0.8H.The initial guess is the device with 9 cavities obtained in §4.2, which has an absorption peak at 9 discrete frequencies and for which E ≈ 0.787.Our is conducted using the MATLAB routine fmincon, and E is approximated using trapezoidal quadrature.
The parameters of the optimised RA, which attains a value of E ≈ 0.982, are given in table 5, and its absorption spectrum is shown in figure 7. The spectrum is predictably flat over the interval [0.8, 1.6] s −1 , and we observe two additional absorption peaks at ω ≈ 1.84 s −1 and ω ≈ 2.07 s −1 .By using (44) to partition the area under the absorption spectrum in figure 7, we find that the two additional peaks are the resonances of the third and second cavity respectively.Additionally, the first three pistons have a minor contribution to absorption over the target interval, which can be demonstrated by computing However, the RA obtained by removing the first three cavities, with absorption spectrum given by the dotted line in figure 8a, only attains a value of E ≈ 0.909.Therefore the purpose of first three cavities cannot be understood in terms of absorption alone.Indeed, the primary role of cavities 1-3 is to intensify the local energy amplification in cavities 4-9, which results in greater absorption by the pistons in these cavities.We interpret that this is because the coupling between the incident plane wave and the resonant mode in the cavity where the amplification occurs is strengthened by the inclusion of cavities 1-3.To illustrate this, we consider (i) a SCA with the same parameters as cavity 4 of the optimised RA, and (ii) an array of 4 cavities consisting of the first four cavities of the optimised RA, but where µ (n) = 0 for n ∈ {1, 2, 3}.Hence both devices contain a single absorbing cavity whose parameters are identical, but device (ii) additionally and where the damping coefficients of the first three pistons are changed to zero (solid line).The absorption spectrum of the unaltered RA also shown in figure 7 is included for comparison (dashed line).(b) The absorption spectra of devices (i) (dotted line) and (ii) (solid line) in the text.As was the case in (a), device (ii) differs from device (i) by the inclusion of 3 non-absorbing cavities.contains 3 non-absorbing cavities positioned to the left.Cavities 5-9 have been omitted from both devices in order to simplify their absorption spectra, which are contrasted in figure 8b.Importantly, the absorption peak of device (ii) is significantly higher than that of device (i), which suggests that the resonant mode in the absorbing cavity of device (ii) is more strongly coupled to the incident plane wave.This coupling occurs sequentially via the propagating Bloch modes of the three non-absorbing cavities, through which the group velocity gradually decreases (recall from §2.4 that the group velocity of the resonant mode is close to zero).In device (i), the coupling between the plane wave and the resonant mode must occur directly, which evidently results in a less-intense energy amplification.In figure 8a, this experiment is also conducted without omitting cavities 5-9 (see solid line), so that both of the devices being compared have 6 identical absorbing cavities.We observe that absorption over [0.8, 1.6] s −1 is reduced much less by setting the damping coefficient of the first three pistons to zero than by removing the first three cavities.This further demonstrates the important role that cavities 1-3 have in facilitating the transfer of energy from the incident plane wave to pistons 4-9 in the optimised RA.We conclude that the two additional absorption peaks in figure 7 are a by-product of the geometry that cavities 2 and 3 must have in order to optimise this energy transfer.

Concluding remarks
Graded arrays of surface piercing vertical barriers have been proposed as a model for studying the rainbow reflection of linear water waves.An example consisting of 10 barriers was studied, where the parameters were selected so that the fundamental frequency of each successive pair of barriers gradually decreases throughout the array.An incident plane wave causes a local energy amplification in this array in the pair of barriers with that fundamental frequency.Using band diagram calculations for the cognate infinite array problem, we demonstrated that the local energy amplification originates from the rainbow reflection effect, since broadband wave signals slow down and become spatially separated based on frequency.
The graded array was modified to incorporate absorption by adding damped, heave-restricted pistons between each adjacent pair of barriers.We solved the resulting boundary value problem using an integral equation/Galerkin technique, giving accurate numerical results.We found parameters for devices which maximise absorption over a discrete set of frequencies by modifying the algorithm of Jiménez et al. [2017].Further tuning using a broadband absorption metric produced a device that achieves near-perfect absorption over a prescribed frequency interval.While several of the cavities in the optimised rainbow absorber absorb a negligable proportion of energy over this interval, they were shown to be indirectly important for transferring energy from the incident wave to other absorbing cavities.
The study of rainbow absorption presented here is a proof of concept that should motivate further exploration of the applications of rainbow phenomena to water waves.For instance, a more realistic three dimensional fluid domain containing a two-dimensional arrays of resonant absorbers could be considered.Future work also seek to validate the local energy amplification predicted using linear water wave theory, which could involve the use of computational fluid dynamics software or wave-tank experiments.Finally, more realistic energy absorption models, including models of wave-to-wire energy conversion, could be explored in the context of rainbow absorption for applications to the design of wave energy converters.

Declaration of Interests
The authors report no conflict of interest.

A
This appendix contains a discussion of how the solution to the underdetermined system of equations given in ( 12) is uniquely specified in the paper.This system can be written in matrix form as where d ∈ R N +1 has entries d (n+1) and d ∈ R N has entries d(n) .The N × (N + 1) matrix C is bidiagonal, where each of the diagonal and superdiagonal entries are 0.5.The null space of C is the span of a single vector v, where (v) j = (−1) (j−1) .The set of solutions to ( 46 ) 2 , so that the grading of the array is as weak as possible.

Figure 1 :
Figure 1: Schematic of the multiple vertical barrier problem.
. Since |φ| remains finite as x → ±∞, we must have A

Figure 2 :
Figure 2: (a) A comparison of the fundamental frequency of two unequal barriers (solid curve) with the fundamental frequency of two equal barriers with submergence (d (0) + d (1) )/2 (dotted curve).Here, l (1) = 1 is fixed.(b) A plot of 1/Q as a function of l (1) /d(0) .There is a very strong agreement between the data used to generate this plot and the model 1/Q = a l (1) /d(0) , where a ∈ R is a coefficient.This demonstrates that Q and l (1) /d (0) are inversely proportional.In both subfigures we use d (0) = 10 m, and the depth of the fluid is fixed to be H = 20 m here and throughout this paper.

Figure 3 :
Figure 3: Colour plots of the velocity potential |φ(x, z)| when forced by 0.1 m waves which are incident from the left of the array of barriers.The frequencies of the incident waves shown here are the desired resonant frequencies, which were used to determine the parameters of the array.Specifically, these frequencies are (a) ω = 1.6 s −1 , (b) 1.5 s −1 , (c) 1.4 s −1 , (d) 1.3 s −1 , (e) 1.2 s −1 , (f) 1.1 s −1 , (g) 1.0 s −1 , (h) 0.9 s −1 (i) 0.8 s −1 .The position of the strongest local energy amplification moves towards the right as the frequency decreases.

Figure 4 :
Figure 4: (a-d) Dispersion curves for infinite periodic arrays of uniform vertical barriers (solid lines) in which d = (d (n−1) + d (n) )/2 and l = l (n) for (a) n = 2, (b) 4, (c) 6 and (d) 8.Here, the values of d(n) and l(n) are those of the graded array shown in figure3.In each subfigure, the dotted line is the dispersion curve for unobstructed waves.Notably, the bandgap begins above the fundamental frequency of a pair of barriers with separation l and submergence d.These fundamental frequencies, which are represented by dashed lines, are (a) ω (n) = 1.5 s −1 , (b) 1.3 s −1 , (c) 1.1 s −1 and (d) 0.9 s −1 .

Figure 5 :
Figure 5: Schematic of the multiple piston cavity problem.Viscous damping is symbolised by dashpots.
For n = 1 the initial guesses of the parameters are d (N) .At each subsequent stage n > 1, we initialise the optimisation with d for j ∈ {2, . . ., n}.In other words, a new cavity which maximises absorption near ω (N −n+1) is added to the left of the array constructed during step n − 1.
Figure6: Absorption spectra of the devices obtained using the steps outlined in §4.2, which is a modification of that used byJiménez et al. [2017].The frequencies which absorption is maximised at each stage of the algorithm are denoted by vertical dotted lines.As desired, the peaks of the spectra coincide with the desired peaks of absorption.

Figure 7 :
Figure 7: (Solid line) the absorption spectrum of the optimised RA.The coloured areas under the curve show the relative contribution of each piston to the total absorption by the device.

Figure 8 :
Figure8: (a) The absorption spectrum of the optimised RA where the first three cavities have been removed (dotted line) and where the damping coefficients of the first three pistons are changed to zero (solid line).The absorption spectrum of the unaltered RA also shown in figure7is included for comparison (dashed line).(b) The absorption spectra of devices (i) (dotted line) and (ii) (solid line) in the text.As was the case in (a), device (ii) differs from device (i) by the inclusion of 3 non-absorbing cavities.
) is the line {C + d + tv|t ∈ R}, where C + denotes the Moore-Penrose inverse of C. The set of points satisfying the constraints 0< d (0) < d (1) < • • • < d (N ) < H is the intersection of N + 1 half-spaces,which is a bounded convex set.The set of feasible solutions to (46) is the line segment F = {C + d + tv|t ∈ [t min , t max ]}, where t min , t max ∈ R are determined by the constraints.We select d ∈ F which minimises N n=1 (d (n) − d(n−1)

Table 5 :
Graded arrays of vertical barriers: rainbow reflection and broadband energy absorption A PREPRINT Parameters of the device designed by maximising the broadband performance metric E. The absorption spectrum of this device is shown in figure7.