2.1 Reflection and transmission kernels for a single cylinder

Consider the subregion of the horizontal plane

where the limits $x_{\pm }$ are chosen so that it contains a single cylinder only (cylinder $n$ , say). Figure 3(a) shows a schematic of the subregion. The wave field incident on the cylinder in $\unicode[STIX]{x1D6FA}$ , $\unicode[STIX]{x1D719}_{in}$ , consists of the ambient incident wave and the scattered wave fields produced by all other cylinders in the array. Let the incident field be decomposed as $\unicode[STIX]{x1D719}_{in}=\unicode[STIX]{x1D719}_{in-}+\unicode[STIX]{x1D719}_{in+}$ , where $\unicode[STIX]{x1D719}_{in-}$ represents the part of the wave field incident from the left-hand boundary of $\unicode[STIX]{x1D6FA}$ , and $\unicode[STIX]{x1D719}_{in+}$ represents the part incident from the right-hand boundary, as illustrated in figure 3(a).

The incident fields are expressed in the plane-wave integral form

where the integration contours are defined by

as shown by the grey curve in figure 3(b), and $\unicode[STIX]{x1D6E4}_{+}=\unicode[STIX]{x1D6E4}_{-}+\unicode[STIX]{x03C0}$ . On the real branches of $\unicode[STIX]{x1D6E4}_{\pm }$ , $\unicode[STIX]{x1D711}$ defines plane waves propagating rightwards for $\unicode[STIX]{x1D6E4}_{-}$ and leftwards for $\unicode[STIX]{x1D6E4}_{+}$ . On the complex branches, it defines plane waves that decay rightwards for $\unicode[STIX]{x1D6E4}_{-}$ and leftwards for $\unicode[STIX]{x1D6E4}_{+}$ . Therefore, equation (2.2) expresses the incident wave field $\unicode[STIX]{x1D719}_{in-}$ as a superposition of plane waves propagating/decaying rightwards, and $\unicode[STIX]{x1D719}_{in+}$ as a superposition of plane waves propagating/decaying leftwards, weighted by amplitude functions $A_{\pm }$ , respectively.

Figure 3. (a) Schematic of region $\unicode[STIX]{x1D6FA}$ defined in (2.1), and corresponding incident and scattered fields, as defined in (2.2) and (2.7b ), respectively. (b) Integration contour $\unicode[STIX]{x1D6E4}_{-}$ (–) in the complex plane, as defined in (2.3), and sampled points used for numerical approximation (●, see § 2.2).

The individual plane-wave components of the incident field, $\unicode[STIX]{x1D711}$ , have the Fourier–Bessel series representations

where $r=\sqrt{(x-x_{n})^{2}+(y-y_{n})^{2}}$ and $\unicode[STIX]{x1D703}=\arg \{(x-x_{n})+\text{i}(y-y_{n})\}$ are, respectively, the radial and azimuthal coordinates associated with the cylinder in $\unicode[STIX]{x1D6FA}$ , and $\text{J}_{m}$ is the first-kind Bessel function of order $m$ . The scattered wave field produced by this incident component is

where $\text{H}_{m}$ is the first-kind Hankel function of order $m$ and $Z_{m}=-\text{J}_{m}^{\prime }(ka)/\text{H}_{m}^{\prime }(ka)$ (e.g. Martin Reference Martin2006). It has the plane-wave integral representation

for $\pm (x-x_{n})\geqslant 0$ and $r>a$ , derived using the Sommerfeld integral representation for Hankel functions (Sommerfeld Reference Sommerfeld1949, § 19).

For the full incident wave field, $\unicode[STIX]{x1D719}_{in}$ , defined by (2.2), the scattered wave field is

for $\pm (x-x_{n})\geqslant 0$ and $r>a$ , where

for $\unicode[STIX]{x1D712}\in \unicode[STIX]{x1D6E4}_{\mp }$ , are scattered amplitude functions. Equation (2.7b ) expresses the scattered field as a superposition of waves propagating/decaying leftwards on the left-hand side of the cylinder, and propagating/decaying rightwards on its right-hand side, as shown in figure 3(a).

The total wave field, $\unicode[STIX]{x1D719}$ , is decomposed into fields propagating/decaying rightwards and leftwards on the left- and right-hand sides of the source via the integral representations

for $\pm (x-x_{n})\geqslant 0$ and $r>a$ , where $B_{\pm }=A_{\mp }+\tilde{B}_{\pm }$ represent outgoing amplitude functions. The $A_{\pm }$ are, accordingly, referred to as incoming amplitude functions. The outgoing amplitude functions are expressed in terms of the incoming amplitude functions using the scattering relations

Here, $R$ and $T$ are, respectively, reflection and transmission kernels, which are defined by

They determine the outgoing amplitude response in direction $\unicode[STIX]{x1D712}$ due to an incoming amplitude forcing in direction $\unicode[STIX]{x1D713}$ .

In order to calculate wave interactions between adjacent subregions, it is convenient to normalise the wave phases to the left- or right-hand side of the region, using the phase-shifted amplitude functions

so that the potential can be expressed as

for $\pm (x-x_{n})\geqslant 0$ and $r>a$ , where $\hat{\unicode[STIX]{x1D711}}_{\pm }=\text{e}^{\text{i}k\{(x-x_{\pm })\text{cos}\,\unicode[STIX]{x1D713}+y\,\text{sin}\,\unicode[STIX]{x1D713}\}}$ . The phase-shifted amplitude functions satisfy scattering relations analogous to (2.10a–b ), but with reflection and transmission kernels, respectively,

2.2 Numerical approximation

Following Montiel et al. (Reference Montiel, Squire and Bennetts2015, Reference Montiel, Squire and Bennetts2016), for numerical calculations the complex branches of the contour $\unicode[STIX]{x1D6E4}_{-}$ (and, hence, $\unicode[STIX]{x1D6E4}_{+}$ ) are truncated, with the magnitude of the imaginary parts restricted to being less than or equal to a prescribed value $\unicode[STIX]{x1D712}_{mx-\text{i}}\in \mathbb{R}_{+}$ . Further, the truncated contours are sampled, with the real branch sampled at $J_{\text{r}}$ regular points, and the complex branches sampled at $J_{\text{i}}+1$ regular points. The sampled points at the ends of the real branch overlap with the sampled points at the real-valued ends of the complex branches, so that the composite contours $\unicode[STIX]{x1D6E4}_{\pm }$ each contain $2J_{\text{i}}+J_{\text{r}}$ sampled points. The black bullets on figure 3(b) are an example of the truncated and sampled version of the contour $\unicode[STIX]{x1D6E4}_{-}$ . The infinite series stemming from the Fourier–Bessel representations are also truncated, with the truncation limits denoted $\pm M$ . Results presented in §§ 3 and 4 use the truncations $M=5$ and $\unicode[STIX]{x1D712}_{mx-\text{i}}=2$ , and discretisations $J_{\text{r}}=J_{\text{i}}=100$ , which were found to provide sufficient accuracy, noting Montiel et al. (Reference Montiel, Squire and Bennetts2016) present relevant convergence analyses for a related problem.

The numerical approximations convert the operator form of the scattering relations satisfied by the phase-shifted amplitude functions to the array relations

where the amplitude vectors $\hat{\boldsymbol{a}}_{\pm }$ and $\hat{\boldsymbol{b}}_{\pm }$ have entries

for $j=1,\ldots ,2J_{\text{i}}+J_{\text{r}}$ , and the matrices $\hat{\unicode[STIX]{x1D64D}}_{\pm }$ and $\hat{\unicode[STIX]{x1D64F}}_{\pm }$ have entries

for $i,j=1,\ldots ,2J_{\text{i}}+J_{\text{r}}$ . The matrices combine the reflection/transmission kernels and numerical integration, where $w_{i}$ ( $i=1,\ldots ,2J_{\text{i}}+J_{\text{r}}$ ) denote the quadrature weights – a composite trapezoidal rule is used for calculations, as in Montiel et al. (Reference Montiel, Squire and Bennetts2015, Reference Montiel, Squire and Bennetts2016).

2.3 Spectrum for an infinite array of cylinders

Let the transfer matrix, $\unicode[STIX]{x1D64B}$ , which maps the amplitude functions on the left-hand side of subregion $\unicode[STIX]{x1D6FA}$ to the amplitude functions on the right-hand side, be defined by

Centring the subregion around the cylinder it contains and setting its width to be $d$ , i.e. $x_{\pm }=x_{n}\pm d/2$ are the midpoints between the cylinder in $\unicode[STIX]{x1D6FA}$ and the adjacent cylinders, the spectrum of the transfer matrix determines the modes supported by the unperturbed array. Figure 4(a) shows an example of the eigenvalues $\unicode[STIX]{x1D706}=\unicode[STIX]{x1D706}_{i}\;(i=1,\ldots ,4J_{\text{i}}+2J_{\text{r}})$ of the transfer matrix in the complex plane. The example problem considered is the resonant case Maniar & Newman (Reference Maniar and Newman1997) identified for $N=100$ cylinders, at frequency $kd=2.78142\;(kd/\unicode[STIX]{x03C0}\approx 0.88535)$ . The eigenvalues appear in reciprocal pairs, i.e. for $|\unicode[STIX]{x1D706}|=1$ if $\exp \{\text{i}\arg (\unicode[STIX]{x1D706})\}$ is an eigenvalue then so is $\exp \{-\text{i}\arg (\unicode[STIX]{x1D706})\}$ , and for $\arg (\unicode[STIX]{x1D706})=0$ if $|\unicode[STIX]{x1D706}|$ is an eigenvalue then so is $1/|\unicode[STIX]{x1D706}|$ . Repeated eigenvalues relate to symmetric and antisymmetric eigenfunctions with respect to $\unicode[STIX]{x1D712}$ .

Figure 4. Spectrum of transfer matrix $\unicode[STIX]{x1D64B}$ , defined in (2.18), for $a/d=0.25$ at the resonant frequency for the 100-cylinder array, $kd/\unicode[STIX]{x03C0}\approx 0.88535$ . (a) Eigenvalues $\unicode[STIX]{x1D706}$ in the complex plane, with eigenvalues supporting head-on incidence in the positive $x$ -direction,  $\exp (\text{i}kd)$ , and the Rayleigh–Bloch mode propagating in the positive $x$ -direction,  $\exp (\text{i}\unicode[STIX]{x1D6FD}_{0}d)$ , indicated. (b) The moduli of the corresponding eigenfunctions associated with rightward propagating/decaying motions, $v_{kd+}$ and $v_{\unicode[STIX]{x1D6FD}_{0}d+}$ , respectively, as functions of the spectral parameter $\unicode[STIX]{x1D712}$ .

The eigenvalues forming the backwards

-shape on the unit circle are the discrete approximation of the continuous spectrum, corresponding to solutions in which wave energy propagates to the far field. In the infinite-array setting, these are solutions forced by a plane ambient incident wave. The angle of the ambient incident wave with respect to the $x$ -axis, $\unicode[STIX]{x1D713}_{am}$ , is related to the argument of the eigenvalue via $\arg (\unicode[STIX]{x1D706})=kd\cos (\unicode[STIX]{x1D713}_{am})$ . Therefore, the eigenvalues at $\unicode[STIX]{x1D706}=1$ correspond to incident waves normal to the array, propagating parallel to the $y$ -axis, and the eigenvalues at the tips of the backwards -shape correspond to forcing at grazing incidence, i.e. parallel to the $x$ -axis. The eigenvalue at the upper tip, which is labelled $\exp (\text{i}kd)$ , corresponds to a rightward propagating incident wave, and the eigenvalue at the lower tip, $\exp (-\text{i}kd)$ (unlabelled), corresponds to a leftward propagating incident wave.

The subset of eigenvalues on the positive real axis form the discrete approximation to the continuous spectrum, generalised to forcing waves that decay exponentially. Eigenvalues smaller than unity decay rightwards, and eigenvalues greater than unity decay leftwards. The set of eigenvalues greater than unity extends beyond the figure limits.

The pair of eigenvalues on the unit circle closest to negative unity correspond to Rayleigh–Bloch wave modes. These eigenvalues are denoted $\exp (\pm \text{i}\unicode[STIX]{x1D6FD}_{0}d)$ , where $\unicode[STIX]{x1D6FD}_{0}\in \mathbb{R}_{+}$ is the Rayleigh–Bloch wavenumber. They form the discrete spectrum, i.e. they correspond to solutions for which wave energy decays with distance away from the array. The eigenvalue in the upper half of the complex plane, which is labelled $\exp (\text{i}\unicode[STIX]{x1D6FD}_{0}d)$ , relates to the mode that propagates rightwards along the array, and the eigenvalue in the lower half of the complex plane, $\exp (-\text{i}\unicode[STIX]{x1D6FD}_{0}d)$ (unlabelled), to the mode that propagates leftwards. As these eigenvalues are separated from the

-shape $\unicode[STIX]{x1D6FD}_{0}>k$ , confirming that the Rayleigh–Bloch waves are shorter than the waves in the surrounding open water.

The definitions of continuous and discrete spectra used here are in terms of quasi-periodicities of possible solutions for a given frequency. This differs from, for example, Porter & Evans (Reference Porter and Evans1999), who define the spectra in terms of frequencies at which solutions exist for a given quasi-periodicity. In simple terms, the spectral parameter used here relates to the directional spectrum, whereas in Porter & Evans (Reference Porter and Evans1999) and others it relates to the frequency spectrum.

Figure 4(b) shows the moduli of the rightward-propagating/decaying components of the eigenfunctions, $v_{kd+}(\unicode[STIX]{x1D712})$ and $v_{\unicode[STIX]{x1D6FD}_{0}d+}(\unicode[STIX]{x1D712})$ , associated with the labelled eigenvalues $\exp (\text{i}kd)$ and $\exp (\text{i}\unicode[STIX]{x1D6FD}_{0}d)$ , respectively. The eigenfunctions are approximated using the eigenvectors associated with the eigenvalues, and are normalised so that their $L_{2}$ -norm is unity. The eigenfunction $v_{kd+}$ displays a sharp spike at $\unicode[STIX]{x1D712}=0$ . This is an approximation of a Dirac delta function, corresponding to the plane ambient incident wave in the infinite-array problem. In contrast, the eigenfunction $v_{\unicode[STIX]{x1D6FD}_{0}d+}$ is smooth, indicating that it is an unforced solution of the infinite-array problem. Its dominant energy is spread around $\unicode[STIX]{x1D712}=0$ , and it has sidebands extending onto the complex branches ( $|\text{Re}(\unicode[STIX]{x1D712})/\unicode[STIX]{x03C0}-\text{Im}(\unicode[STIX]{x1D712})|>0.5$ ).

3.1 Solution method for finite, perturbed arrays

Positional disorder is introduced into the problem via random perturbations in the locations of the cylinders. The perturbed locations of the cylinder centres are denoted

A parameter $\unicode[STIX]{x1D700}$ is used to control the perturbation strength, with the perturbation vectors expressed as

where $\unicode[STIX]{x1D707}_{n}$ and $\unicode[STIX]{x1D708}_{n}$ are randomly selected from uniform distributions over the interval $[-0.5,0.5]$ . The perturbation strength is bounded by $\unicode[STIX]{x1D700}<1-2a/d$ so that cylinders cannot overlap. Figure 5 shows a schematic of the perturbed geometry.

Figure 5. Schematic plan view of the perturbed geometry.

A contiguous sequence of subregions $\unicode[STIX]{x1D6FA}_{n}\;(n=1,\ldots ,N)$ is chosen, as shown in figure 5, with each subregion containing the cylinder of the corresponding index. (As indicated in figure 5, it is unnecessary for the subregions to contain the corresponding unperturbed cylinders.) The amplitude functions and reflection/transmission matrices associated with subregion $\unicode[STIX]{x1D6FA}_{n}$ are assigned superscripts $(n)$ , with the reflection/transmission matrices calculated using the method outlined in §§ 2.1–2.2.

For $1\leqslant p\leqslant q\leqslant N$ , $\unicode[STIX]{x1D6FA}_{p,q}$ denotes the region contained between the left-hand boundary of $\unicode[STIX]{x1D6FA}_{p}$ and the right-hand boundary of $\unicode[STIX]{x1D6FA}_{q}$ , with figure 5 showing $\unicode[STIX]{x1D6FA}_{1,N}$ as an example. Let $\hat{\unicode[STIX]{x1D64D}}_{\pm }^{(p,q)}$ and $\hat{\unicode[STIX]{x1D64F}}_{\pm }^{(p,q)}$ denote the reflection and transmission matrices for this region, i.e.

The reflection and transmission matrices for $\unicode[STIX]{x1D6FA}_{p,q}$ with $\unicode[STIX]{x1D6FA}_{q+1}$ appended to its right-hand side, $\hat{\unicode[STIX]{x1D64D}}_{\pm }^{(p,q+1)}$ and $\hat{\unicode[STIX]{x1D64F}}_{\pm }^{(p,q+1)}$ , can be calculated from $\hat{\unicode[STIX]{x1D64D}}_{\pm }^{(p,q)}$ , $\hat{\unicode[STIX]{x1D64F}}_{\pm }^{(p,q)}$ , $\hat{\unicode[STIX]{x1D64D}}_{\pm }^{(q+1)}$ and $\hat{\unicode[STIX]{x1D64F}}_{\pm }^{(q+1)}$ via

where $\unicode[STIX]{x1D644}$ is the identity matrix of dimension $2J_{\text{i}}+J_{\text{r}}$ . These relations are derived from the scattering relations for an individual region (2.15), on the basis that (i) the outgoing wave field on the right-hand side of cylinder $q$ is the incoming field on the left-hand side of cylinder $q+1$ , and (ii) the outgoing field on the left-hand side of cylinder $q+1$ is the incoming field on the right-hand side of cylinder $q$ , i.e.

Versions of relations (3.4) have been used by Botten et al. (Reference Botten, White, Asatryan, Langtry, de Sterke and McPhedran2004), Peter & Meylan (Reference Peter and Meylan2009) and Bennetts (Reference Bennetts2011), among others, for electromagnetic, water wave and acoustic problems, respectively.

The reflection and transmission matrices $\hat{\unicode[STIX]{x1D64D}}_{\pm }^{(1,n)}$ and $\hat{\unicode[STIX]{x1D64F}}_{\pm }^{(1,n)}$ for $n=1,\ldots ,N$ are obtained by beginning with the leftmost region, $\unicode[STIX]{x1D6FA}_{1}$ , recursively appending the next region to the right and applying (3.4). The process is initialised with the identities

Similarly, the reflection and transmission matrices $\hat{\unicode[STIX]{x1D64D}}_{\pm }^{(n,N)}$ and $\hat{\unicode[STIX]{x1D64F}}_{\pm }^{(n,N)}$ for $n=1,\ldots ,N$ are obtained by beginning at the rightmost region, $\unicode[STIX]{x1D6FA}_{N}$ , then recursively appending the next region to the left and applying similar equations to (3.4) for region $\unicode[STIX]{x1D6FA}_{p,q}$ with $\unicode[STIX]{x1D6FA}_{p-1}$ appended to its left-hand side.

Figure 6. (a) As in figure 2 but for a Rayleigh–Bloch ambient incident wave. (b) Corresponding indices of the cylinders experiencing the maximum loads.

The sampled amplitude functions throughout the array are obtained from the reflection and transmission matrices via

and

for $n=1,\ldots ,N-1$ . The vectors $\hat{\boldsymbol{a}}_{-}^{(1)}$ and $\hat{\boldsymbol{a}}_{+}^{(N)}$ are, respectively, the sampled directional spectra for the prescribed rightward and leftward propagating/decaying components of the ambient incident wave field, $\hat{A}_{am-}=\hat{A}_{-}^{(1)}$ and $\hat{A}_{am+}=\hat{A}_{+}^{(N)}$ , defined via

where $\hat{\unicode[STIX]{x1D711}}_{-}^{(1)}$ and $\hat{\unicode[STIX]{x1D711}}_{+}^{(N)}$ are the phase-shifted plane-wave potentials defined below (2.13) for subregions $\unicode[STIX]{x1D6FA}_{1}$ and $\unicode[STIX]{x1D6FA}_{N}$ , respectively. Bennetts & Squire (Reference Bennetts and Squire2009) originally presented equivalent expressions to (3.7) for a problem involving multiple rows of floating bodies. Calculation of the amplitude functions completes the solution process.

Figure 7. Example realisations of load profiles on 100-cylinder arrays for $a/d=0.25$ , $kd/\unicode[STIX]{x03C0}\approx 0.88535$ , and (a) no perturbation, (b) perturbation strength $\unicode[STIX]{x1D700}=0.01$ , and (c) $\unicode[STIX]{x1D700}=0.2$ .

3.2 Numerical results

To simplify the analysis, the ambient incident wave field is set to be the rightward-propagating Rayleigh–Bloch mode, i.e. $\hat{A}_{am-}=v_{\unicode[STIX]{x1D6FD}_{0}d+}$ and $\hat{A}_{am+}=0$ . Figure 6(a) shows the maximum (normalised) load on a 100-cylinder array computed from (1.1), as a function of frequency, $kd/\unicode[STIX]{x03C0}$ , for different perturbation strengths, $\unicode[STIX]{x1D700}$ . It is equivalent to figure 2 for the Rayleigh–Bloch incident wave rather than the plane incident wave, with the frequency range truncated due to the cutoff at $kd/\unicode[STIX]{x03C0}\approx 0.88574$ . For frequencies immediately beyond the cutoff, the Rayleigh–Bloch wavenumbers become complex and the Rayleigh–Bloch wave decays along the array (as identified by Thompson et al. Reference Thompson, Linton and Porter2008).

For the unperturbed problem, the Rayleigh–Bloch incident wave excites a resonance at the same frequency as the plane incident wave, $kd/\unicode[STIX]{x03C0}\approx 0.88535$ , preceded by a series of maxima and minima. This is expected, as these features are caused by interactions between rightward- and leftward-propagating Rayleigh–Bloch waves (as explained in § 1, cf. Thompson et al. Reference Thompson, Linton and Porter2008). Moreover, in an identical fashion to the plane-wave case, perturbations damp the resonance for the Rayleigh–Bloch incident wave field.

Figure 6(b) shows the index of the cylinder at which the maximum load is attained, as a function of frequency and for the different perturbation strengths. For low frequencies, the maximum is attained at the farthest cylinder from the source of the incident wave field (the right-hand end of the array, $n=100$ ), due to cylinders refracting the rightward-propagating incident wave energy in towards the array (Maniar & Newman Reference Maniar and Newman1997). This low-frequency behaviour is consistent for the different perturbation strengths, as the perturbations are small in comparison to the wavelength, $2\unicode[STIX]{x03C0}/k$ , in this regime. For the unperturbed array, the maximum load is attained at the right-hand end of the array up to $kd/\unicode[STIX]{x03C0}\approx 0.4{-}0.5$ . As the frequency increases beyond this limit, the location of the maximum begins to jump around erratically, as it moves between different local maxima of comparable magnitudes along the array. For the perturbed arrays, the maximum load departs the right-hand end of the array for lower frequencies than the unperturbed array, with the departure frequency decreasing as the perturbation strength increases. The location of the maximum then transitions relatively smoothly towards the front cylinder (left-hand end of the array, $n=1$ ), as frequency increases.

Figure 7 shows example profiles of the loads on the array at the resonant frequency $kd/\unicode[STIX]{x03C0}\approx 0.88535$ . It shows the profile for the unperturbed array and for two random realisations of the perturbed array, for perturbation strengths $\unicode[STIX]{x1D700}=0.01$ and 0.2. The profile for the unperturbed array is qualitatively identical to that presented by Maniar & Newman (Reference Maniar and Newman1997) for a plane incident wave – it is an approximately symmetric hump, with the maximum load attained at the centre ( $n=51$ ), due to coherence between strongly excited rightward- and leftward-propagating Rayleigh–Bloch modes (Thompson et al. Reference Thompson, Linton and Porter2008). The weak perturbation ( $\unicode[STIX]{x1D700}=0.01$ ) preserves the near symmetry but slightly shortens the length of the resonant hump and damps the maximum load. The strong perturbation ( $\unicode[STIX]{x1D700}=0.2$ ) fundamentally changes the qualitative properties of the profile, eliminating symmetry and resonance altogether. The profile attenuates from the front of the array to the rear, with accompanying oscillations, indicating a localised state has been reached.

4.1 Perturbed Rayleigh–Bloch wavenumber calculation

For the unperturbed array forced by a Rayleigh–Bloch ambient incident wave, the wave field along the array is dominated by rightward- and leftward-propagating Rayleigh–Bloch modes. For a perturbed array, it is hypothesised that, similarly, the wave field is dominated by perturbed Rayleigh–Bloch modes, with complex-valued wavenumbers $\unicode[STIX]{x1D6FD}_{+}(\unicode[STIX]{x1D700})\in \mathbb{R}_{+}+\text{i}\mathbb{R}_{+}$ for the rightward-propagating mode and $\unicode[STIX]{x1D6FD}_{-}(\unicode[STIX]{x1D700})\in \mathbb{R}_{-}+\text{i}\mathbb{R}_{-}$ for the leftward-propagating mode. It is expected that $\unicode[STIX]{x1D6FD}_{-}\approx -\unicode[STIX]{x1D6FD}_{+}$ , on average at least, due to symmetry of the perturbations, and for the unperturbed array $\unicode[STIX]{x1D6FD}_{\pm }(0)=\pm \unicode[STIX]{x1D6FD}_{0}$ . The change in the real parts of the wavenumbers, with respect to the unperturbed Rayleigh–Bloch wavenumbers, induces phase shifts in the modes, and the introduction of imaginary components produces attenuation of the modes, resulting in localisation for sufficiently large attenuation rates.

The perturbed wavenumbers, $\unicode[STIX]{x1D6FD}_{\pm }$ , are extracted from wave fields for the perturbed arrays, using a version of an approach proposed by Bennetts & Peter (Reference Bennetts and Peter2013). For each realisation of the perturbed array, the wave fields are calculated for

1. (i) a rightward-propagating Rayleigh–Bloch incident wave ( $\hat{A}_{am-}=v_{\unicode[STIX]{x1D6FD}_{0}d+}$ and $\hat{A}_{am+}=0$ , as in § 3),

2. (ii) and a leftward-propagating Rayleigh–Bloch incident wave ( $\hat{A}_{am-}=0$ and $\hat{A}_{am+}=v_{-\unicode[STIX]{x1D6FD}_{0}d-}=v_{\unicode[STIX]{x1D6FD}_{0}d+}$ ),

where $v_{-\unicode[STIX]{x1D6FD}_{0}d-}$ is the component of the eigenfunction associated with $-\unicode[STIX]{x1D6FD}_{0}$ that propagates/decays leftwards, and $v_{\unicode[STIX]{x1D6FD}_{0}d+}$ is defined in § 2.3.

The rightward- and leftward-propagating directional amplitude functions along the array for these wave fields are denoted

for $n=1,\ldots ,N+1$ , respectively, where the undefined amplitude functions at the boundaries of the array are replaced by the equivalent functions of the ambient incident wave fields,

The hypothesis that the wave fields are dominated by perturbed Rayleigh–Bloch modes is translated into the ansatzes

where $a_{\pm }^{(\bullet )}(\unicode[STIX]{x1D700})$ and $b_{\pm }^{(\bullet )}(\unicode[STIX]{x1D700})$ are coefficients/amplitudes, $v_{\unicode[STIX]{x1D6FD}_{\pm }d+}$ is the component of the eigenfunction associated with $\unicode[STIX]{x1D6FD}_{\pm }$ that propagates/decays rightwards, and $v_{\unicode[STIX]{x1D6FD}_{\pm }d-}$ is the component propagating/decaying leftwards.

Ansatzes (4.3a–b ) are combined, as in Bennetts & Peter (Reference Bennetts and Peter2013), to produce the expression

where

and $\text{eig}(\bullet )$ denotes the set of eigenvalues of the included matrix. The values of the amplitude functions corresponding to wave components propagating parallel to the $x$ -axis only are used because the ansatzes essentially reduce the problem to a single dimension.

Figure 8. Eigenvalues of matrix $\unicode[STIX]{x1D648}^{(n)}$ , given in (4.4), corresponding to rightward-propagating (a,b) and leftward-propagating (c,d) Rayleigh–Bloch modes for $a/d=0.25$ , $kd/\unicode[STIX]{x03C0}\approx 0.88535$ and $\unicode[STIX]{x1D700}=0.2$ . (a,c) Real parts of the natural logarithms of the modes (phases) and (b,d) imaginary parts (attenuation), calculated as means of an ensemble of 500 simulations (○). Straight-line fits (–) used to extract the Rayleigh–Bloch wavenumbers are overlaid.

Figure 8 shows an example of natural logarithms of eigenvalues of the matrix $\unicode[STIX]{x1D648}^{(n)}$ , as a function of cylinder index $n$ , for a 100-cylinder array with perturbation strength $\unicode[STIX]{x1D700}=0.2$ at the resonant frequency for the unperturbed array, $kd/\unicode[STIX]{x03C0}\approx 0.88535$ . The eigenvalues are separated into real and imaginary parts (a,c and b,d, respectively), and those corresponding to $\unicode[STIX]{x1D6FD}_{+}$ and $\unicode[STIX]{x1D6FD}_{-}$ (a,b and c,d, respectively). The results are the means of ensembles of 500 randomly generated realisations of the perturbations.

The log-eigenvalue profiles are approximately linear – noise is visible in the imaginary parts but not in the real parts, which are dominated by the eigenvalues of the underlying unperturbed array. Straight-line fits are overlaid on the profiles, from which the non-dimensional Rayleigh–Bloch wavenumbers, $\unicode[STIX]{x1D6FD}_{\pm }d$ , are extracted as the slopes. The fits are made to the interior of the profiles only to avoid contamination by local effects of the array ends, just visible in the imaginary parts. The extracted Rayleigh–Bloch wavenumbers display the expected symmetry (to three decimal places), with $\unicode[STIX]{x1D6FD}_{+}d\approx 3.02058+0.01536\text{i}$ and $\unicode[STIX]{x1D6FD}_{-}d\approx -3.02055-0.01548\text{i}$ .

Extracting array wavenumbers (in this case Rayleigh–Bloch wave numbers) from mean log-eigenvalue profiles is more stable than Bennetts & Peter (Reference Bennetts and Peter2013)’s method of extracting array wavenumbers from log-eigenvalue profiles of individual realisations and then averaging. The Rayleigh–Bloch wavenumbers are not those of the effective (mean) wave field, as the amplitudes and phases are separated prior to averaging.

Figure 9. (a,c) Phase changes of perturbed Rayleigh–Bloch wavenumbers relative to the unperturbed wavenumbers, as functions of disorder strength for $a/d=0.25$ and $kd/\unicode[STIX]{x03C0}\approx 0.88535$ . (b,d) Corresponding attenuation rates. (a,b) Results for the rightward-propagating/attenuating Rayleigh–Bloch mode, and (c,d) results for the leftward-propagating/attenuating mode. Results are calculated using a 100-cylinder array (▫) and a 200-cylinder array (▪).

4.2 Numerical results

Figure 9 shows the extracted wavenumbers as functions of perturbation strength for the resonant frequency, and using ensembles of 100 realisations for each perturbation strength. The real parts (phases) are shown as relative deviations from the Rayleigh–Bloch wavenumbers of the unperturbed array,

in (a,c), with the imaginary parts (attenuation rates) shown in (b,d). Wavenumbers supporting rightward-propagating/attenuating modes are shown in (a,b), with those for the leftward-propagating/attenuating modes in (c,d). Results are shown for 100- and 200-cylinder arrays.

The results again show the expected symmetry $\unicode[STIX]{x1D6FD}_{-}\approx -\unicode[STIX]{x1D6FD}_{+}$ . The phase changes are similar for the 100- and 200-cylinder arrays, indicating that the real parts of the Rayleigh–Bloch wavenumbers do not depend on array length (beyond $N=100$ , at least). Deviations between the phases for the two array lengths are evident for weak perturbations, approximately $\log _{10}\unicode[STIX]{x1D700}<-1.5\;(\unicode[STIX]{x1D700}<0.03162)$ , with the phase changes approximately zero in this regime for the 200-cylinder arrays, i.e. $\unicode[STIX]{x1D6FD}_{\pm }(\unicode[STIX]{x1D700})\approx \pm \unicode[STIX]{x1D6FD}_{0}$ for $\unicode[STIX]{x1D700}\ll 1$ as expected, but being a small finite value for the 100-cylinder array. The errors for the 100-cylinder array are due to 100 cylinders not providing an interval long enough to average out the presence of oscillations caused by the circular wave fields generated at the array ends, which compromise the validity of ansatzes (4.3). Stronger perturbations, $\log _{10}\unicode[STIX]{x1D700}>-1.5$ , clearly reduce the modulus of the phase of the Rayleigh–Bloch wavenumbers, with the reduction steadily increasing as the perturbation strength increases.

In the weak-perturbation regime, $\log _{10}\unicode[STIX]{x1D700}<-1.5$ , deviations between the two array lengths are exacerbated for the attenuation rates, as the attenuation is negligible so that the oscillations caused by the circular wave fields dominate the imaginary-component profiles. The 200-cylinder array provides attenuation rates that are approximately zero for very weak perturbations, approximately $\log _{10}\unicode[STIX]{x1D700}<-2\;(\unicode[STIX]{x1D700}<0.01)$ , as expected, whereas attenuation rates for the 100-cylinder array are visibly non-zero in this regime. For stronger perturbations, $\log _{10}\unicode[STIX]{x1D700}>-1.5$ , the moduli of the attenuation rates increase with increasing perturbation strength and the attenuation rates for the 100-cylinder array rapidly transition to agreement with those of the 200-cylinder array. This marks the regime in which attenuation dominates the profiles of the imaginary components of the eigenvalues of $\unicode[STIX]{x1D648}^{(n)}$ .

Figure 10. (a) Maximum load (–) and load on the first cylinder (– - –), caused by a rightward-propagating Rayleigh–Bloch ambient incident wave on a 200-cylinder array for $a/d=0.25$ and $kd/\unicode[STIX]{x03C0}\approx 0.88535$ , as functions of disorder. (b) A corresponding load-profile realisation for $\log _{10}\unicode[STIX]{x1D700}=-2$ (●), and the profile for the unperturbed problem with $kd/\unicode[STIX]{x03C0}\approx 0.88532$ (○). (c) As in (b) but for $\log _{10}\unicode[STIX]{x1D700}=-1.25$ (●) and $kd/\unicode[STIX]{x03C0}\approx 0.88487$ (○).

Figure 10(a) shows the corresponding maximum load on the 200-cylinder array produced by a rightward-propagating Rayleigh–Bloch ambient incident wave field and the load on the front cylinder ( $n=1$ ), as functions of perturbation strength. The resonance in the maximum load for the unperturbed 200-cylinder array at $kd/\unicode[STIX]{x03C0}\approx 0.88535$ manifests as two symmetric humps in the load profile, where each hump is almost identical to that shown in figure 7(a) for the 100-cylinder array. (A stronger, single-hump resonance for the 200-cylinder array, analogous to that shown by figure 7(a) for the 100-cylinder array, is obtained at a slightly higher frequency.)

For weak perturbations, $\log _{10}\unicode[STIX]{x1D700}<-1.5$ , resonance in the maximum load is damped, as indicated in figure 7(b) for the 100-cylinder array, with the damping gradually increasing as the perturbation strength increases. Figure 10(b) shows an individual load profile in the weak-perturbation regime, with $\log _{10}\unicode[STIX]{x1D700}=-2$ , for which $\unicode[STIX]{x1D6FD}_{\pm }d\approx \pm 3.109\pm 0.000\text{i}$ . The profile is similar to the two-hump resonant profile occurring for the unperturbed array at this frequency, but with the resonances damped and the hump lengths reduced. The load profile for the unperturbed array with frequency $kd/\unicode[STIX]{x03C0}\approx 0.88532$ , for which $\unicode[STIX]{x1D6FD}_{0}d\approx 3.109$ , is superimposed on the panel. It is almost identical to the profile for the perturbed array, indicating that a weakly perturbed array acts like an unperturbed array with a slightly shortened cylinder spacing and/or to slightly shifted forcing. (A similar profile can be obtained at the resonant frequency by increasing the number of cylinders in the unperturbed array to $N=204$ . Alternatively, reducing the perturbed array to $N=196$ cylinders produces a profile approximating the two-hump resonance for the unperturbed array with $N=200$ .)

For stronger perturbations, $\log _{10}\unicode[STIX]{x1D700}>-1.5$ , the maximum load is approximately equal to the load on the front cylinder, as the attenuation rates are, in general, sufficiently large to produce localised states, as shown in figure 7(c) for the 100-cylinder array. A weak resonance in the maximum load around $\log _{10}\unicode[STIX]{x1D700}=-1.25\;(\unicode[STIX]{x1D700}\approx 0.05623)$ causes it to deviate away from the load on the front cylinder. Figure 10(c) shows an individual load profile for $\log _{10}\unicode[STIX]{x1D700}=-1.25$ , for which $\unicode[STIX]{x1D6FD}_{\pm }\approx \pm 3.094\pm 0.003\text{i}$ , with the load profile for the unperturbed array at frequency $kd/\unicode[STIX]{x03C0}\approx 0.88487$ superimposed, for which $\unicode[STIX]{x1D6FD}_{0}\approx 3.094$ . The profile for the unperturbed array is resonant, as indicated by the three humps and the symmetry of the profile, although it is a weaker resonance than the two-hump resonance at $kd/\unicode[STIX]{x03C0}\approx 0.88535$ . The three humps are visible in the profile for the perturbed array, driving the resonance in the maximum load around $\log _{10}\unicode[STIX]{x1D700}=-1.25$ . However, attenuation produced by the imaginary component of the Rayleigh–Bloch wavenumber significantly reduces the hump peaks, so that the resonance is weak.