IMAGE CONDITIONS FOR POLAR AND CYLINDRICAL COORDINATE SEPARATION-OF-VARIABLES ACOUSTIC MULTIPLE SCATTERING MODELS WITH PERFECTLY REFLECTING FLAT BOUNDARIES

Summary This article addresses efﬁcient implementation of the method of images for acoustic multiple scattering models (MSM) with perfectly reﬂecting ﬂat boundaries. Time-harmonic problems are ﬁrst solved in the polar coordinate system for circular scatterers; then the models are extended to the cylindrical coordinate system with (semi-)inﬁnitely long circular cylinders. The MSM in this article is based on the method of separation of variables and Graf’s addition theorem. Derivations are provided for ‘image conditions’ which relate the unknown coefﬁcients of outgoing waves from image scatterers with those of real counterparts. The method of images is applied to wedge-shaped domains with apex angles of π/ n rad for a positive integer n . Image conditions make the MSM numerically more efﬁcient: the system of linear equations for unknown coefﬁcients is formulated 2 n times faster; its memory requirements are reduced by 4 n 2 times for direct solvers. The proposed model is applied to a benchmark wedge in ocean environment with n = 64. Good agreementisobservedbetweentheMSMwithimageconditionsandtheboundaryelementmethod.Furthermore,half-andquarter-spacemeasurementsinananechoicchamberareinaccordancewiththecorrectuseofimageconditions.Incorrectimageconditionsreportedelsewhereforpolarcoordinatesarealsodiscussed.


Introduction
Multiple scattering (1) of waves has been one of the classical problems which has captured the imagination of research communities: acoustics and beyond. There are several numerical modalities such as the boundary element method and finite element method for analysis. These are very versatile in terms of scatterer shapes, boundary types, material properties, and so on. However, they are usually computation-intensive. A large number of scatterers can be challenging to consider.
In contrast, when all the scatterers are circles in 2D domains (equivalent to infinitely long circular cylinders in 3D) or can be approximated as such, scattering problems can be addressed much more efficiently by using the method of separation of variables and addition theorems (2)(3)(4)(5)(6)(7). These combinations are the constituents of the multiple scattering model addressed in this article, which is equivalent to the T-matrix method for circular scatterers (1). However, flat interfaces such as the ground or free surface of water pose challenges of implementation alongside circles. Therefore, the effect of such flat boundaries has been achieved by the method of images as long as these boundaries are infinitely long and perfectly reflecting (8): either rigid or pressure release.
Since the inception of sonic crystals using horizontal circular cylinders as a potential tool for outdoor noise mitigation, multiple scattering works have been carried out in conjunction with the method of images due to the inevitable presence of outdoor ground (9)(10)(11). Underwater sound scattering underneath the free surface has also been analysed with similar approaches (12,13). Beyond these half-plane problems, understanding of sound fields near right-angled corners or acuteangled wedges are also of great interest due to an abundance of such geometries, even without scatterers (14,15). Many man-made structures above ground or underwater can be considered as the vertical side of a corner or slant side of wedges subject to the wavelength of interest. There are ample examples of rigid corners in practice: for example, hard noise barriers and textbook demonstration of sound fields near a corner. A corner or wedges having both rigid and pressure release sides may be encountered during underwater sound propagation near a calm free surface and an adjacent rigid wall, a large vessel or slant sea floor near a beach. Multiple scattering problems within these corners and limited cases of wedges can also be analysed by using the method of images, and are topics of this presentation.
The method of images for multiple scattering models can be performed in two different ways. The first approach (16) is to assign additional computational addresses for the coefficients of mirrored scattered field during the formulation of the system of linear equations. This is straightforward, but inefficient due to the enlarged size of a matrix and does not provide any advancement in multiple scattering modelling. The second procedure is to use 'image conditions' to merge these addresses for image scatterers with those of the corresponding real scatterers in the matrix. This terminology of 'image condition' is named such because it is a necessary condition to make one scatterer an image of the other (17). Therefore, such use of image conditions reduces memory requirements and results in faster computation compared to the former approach. This latter procedure has been used by others for half-plane problems with a single flat interface (9,12,13,18); however, has not been attempted for corners and wedge-shaped domains. A previous attempt (16) to use the method of images for a rigid corner did not adopt image conditions. Detailed derivations of these image conditions and their verifications are the main objectives of this article.
For the current research, a constituent interface of the method of images is an infinitely long axis (or plane), especially an inclined one with an arbitrary angle from, say, the x-axis (or zx-plane). Once image conditions for such are developed, successive applications thereof can assemble image conditions for corners and wedges with apex angles of integer division of π rad, which leads to a finite number of image domains. These image conditions for an inclined half-plane, corner and wedge-shaped domains are novelty of this article.
Therefore, the method of images can realise acoustic domains shown in Fig. 1, for example. For 3D problems, cylinders are placed parallel to one another and can be semi-infinitely long as shown in the diagram; they can also be infinitely long without the dark grey planes. The geometric configurations of 2D problems are equivalent to those of infinitely long cylinders; then, a cross-section of domain perpendicular to the cylinder axes is effectively considered. The domains in Fig. 1 can be realised regardless of using the aforementioned image conditions. However, the image conditions can be introduced across the light grey planes to achieve numerical efficiency.
The application of the method of images leads to the geometric symmetry between image and real entities. Symmetries can also be exploited in multiple scattering problems, without using the method of images, when (real) scatterers are placed in specific arrangements such as infinite number of periodic identical circular cylinders (2,19,20); two equal cylinders along the centreline of a long Cylinders are semi-infinitely long; their axes are parallel to the z-axis. Image conditions are defined across the light grey planes, but not across the dark grey planes channel (21); circular or linear periodic arrays (22). Then, the reduction in the number of unknowns is achieved for real scatterers, while the reduction in the current article is for image scatterers. Unlike these works, the real scatterers in our work can be randomly located. The article is organised as follows. First, conventions of the article are explained, followed by the description of incident and scattered fields. Next, image conditions are derived for an arbitrarilyrotated flat interface: separately for pressure release and rigid interfaces. Use of the method of images is explained in building wedge-shaped domains; subsequently image conditions are derived for perfectly reflecting wedges. This enables novel multiple scattering models efficient for acoustic domains in the shape of corners and wedges; savings in computational resources are highlighted. A spherical-wave incidence upon infinitely long circular cylinders is also formulated, followed by the oblique incidence of a plane wave in 3D. Then, image conditions are also discussed for radiation problems and the modified Helmholtz equation. Verifications are presented, first for measurements in an anechoic chamber for half-and quarter-space problems. Numerical simulations are carried out for 2D and 3D corners. Further numerical examples are presented for an ocean-environment benchmark wedge configuration and are validated by the boundary element method. Discussion is held about incorrect image conditions. Finally, conclusions are drawn.

Conventions in this article
In this article, Cartesian (x, y) and polar coordinates (r, θ ) are related by x = r cos θ and y = r sin θ . The following conventions are adopted to define position vectors: see Fig. 2. A double-subscript A general field point is indicated by (x, y); b is a field point on the x -axis. For 3D problems with infinitely long cylinders, the diagram can be interpreted as a projection to the xy-plane. Note that the coordinate system 'local' to scatterers are not rotated vector r co is the position vector of an origin O c with respect to another origin O o . Its radial magnitude and polar angle are denoted respectively by r co and θ co . A single-subscript vector r c (with r c and θ c ) is the position vector of a general field point (x, y) with respect to a local origin O c . A non-subscript position vector r (with r and θ) represents the general field point (x, y), referenced to the global origin O o . These three vectors are related by r = r co + r c . A tilde ( ) indicates the image counterpart of a real entity: for example, O c represents the image version of O c . Therefore, a position vector r can also be represented by r = r co + r c . Furthermore, this subscript convention is also applicable to Cartesian coordinates: for example, (x co , y co ) corresponds to r co .
The following notations are used for outgoing (ψ n ) and regular (φ n ) wavefunctions in polar coordinates (1, p. 125): These are equally applicable to position vectors with subscripts. The Hankel and Bessel functions are denoted by H (1) n (·) and J n (·), respectively: both are of the first kind and order n. Use of H (1) n (·) for modelling an outgoing wave implies the time-harmonic convention of exp(−iωt). As widely used, i, ω and t denote the imaginary unit, angular frequency and time, respectively. Note that a wavenumber k, for a medium outside scatterers, is implicit at ψ n (r) and φ n (r).
For cylindrical coordinates (r, θ, z), a 3D position vector is represented by R whose radial variable is R = r 2 + z 2 . The aforementioned subscript conventions are also applied to R and its (r, θ, z). Furthermore, a wavenumber k is assumed for polar coordinates and 2D Cartesian; its 3D counterpart is represented with a bar likek. Wavenumbers also have components whose subscripts denote directions (or coordinate dimension) such as k = k 2 x + k 2 Here, k r is equivalent to k, but used exclusively in 3D problems.

Incident fields
Imagine the two-dimensional acoustic field, say, on the xy-plane with an infinitely long flat boundary, for example, the x -axis in Fig. 2. For cylindrical coordinates in three dimensions, Fig. 2 can be interpreted as a projection to the xy-plane. This flat boundary is perfectly reflecting: either rigid or pressure release. A finite number of circular objects or scatterers are placed at any location in the real domain (for example, y < 0). These scatterers can be of any size and property, as long as they are circles and do not overlap one another; they can also make contact with the flat boundary, but do not penetrate it. Then, a small-amplitude time-harmonic sound is incident from either a circularwavefront line source (or a spherical-wavefront point source in 3D) at any real location outside the scatterers or a plane wave from any direction (preferably away from the boundary) through a homogeneous fluid medium.
For an example of a line source in two dimensions, the incident pressure field at r = (x, y) may be represented by where B s is the strength of a source s located at O s in Fig. 2.
The presence of perfectly reflecting boundaries permits the use of the method of images (8). By analogy with optics, mirrored sources and scatterers are created in the image domain: for example, y > 0 in Fig. 2. The image counterpart of (3.1) may be expressed by where B s = B s across a rigid interface; B s = −B s across pressure release. It is possible to extend this 2D problem to the 3D equivalent for infinitely long cylinders (4, 18, 23), parallel to one another, whose circular cross-sections are invariant along their axes (for example, z-axis). Then, the incident pressure field excited by a point source located at O s in 3D may be evaluated by (24, p. 710) An image counterpart can be similarly expressed by replacing the subscript s by s. As the integrand of (3.3) suggests, even for this spherical-wave incidence, the underlying configuration can be addressed in polar coordinates with height (that is cylindrical coordinates); hence, the image conditions of the current article are equally applicable. The 3D topic is addressed in section 5.4; until then, this article considers only 2D configurations.

Scattered fields
A total external acoustic field in scattering problems is therefore composed of both incident and scattered fields. Under the method of images, contributions from both real and image sources are known to satisfy perfectly reflecting boundary conditions; in fact, it is elementary to prove that (3.1) and (3.2) satisfy the Neumann condition for B s = B s and the Dirichlet condition for B s = −B s . Therefore, it suffices to consider only the fields scattered by scatterers to derive image conditions due to linearity. Furthermore, it is also sufficient to examine a single pair of real and image scatterers, that is cth pair for example. The rest pairs should be subject to the same image conditions. In this article, the multiple scattering problems are addressed by solving the Helmholtz equation. It is well known that the Helmholtz equation separates in polar and cylindrical coordinates. Separation of variables in polar coordinates is mainly addressed until section 5.4; then its findings can be easily extended to the cylindrical coordinate system since the z-dimension is also separated. Therefore, in polar coordinates, the pressure field scattered from a real object centred at O c can be expressed as (1, p. 125) which satisfies the Sommerfeld radiation condition and (∇ 2 + k 2 )p c = 0 with ∇ 2 = ∂ 2 /∂r 2 + (∂/∂r)/r + (∂ 2 /∂θ 2 )/r 2 . The scattered field from the corresponding image object at O c may be represented by One of the aims of this article is to show that the unknown coefficient C c m can be expressed in terms of C c m as long as a flat boundary is perfectly reflecting; that is C c m will be the only unknown to be solved. This relationship between C c m and C c m is defined as an 'image condition'. Image conditions for a half-plane with its flat interface being either x-or y-axis were presented by others (9, 12, 13, 18). To analyse a corner or wedges, image conditions for an infinitely long flat boundary oblique from, say, the x-axis have to be established. First, an image condition for a pressure release flat boundary is derived, and then followed by that for a rigid straight interface. Image conditions for single boundaries can be combined to construct those for corners and wedges.

Image condition: pressure release interface
A pressure release boundary is characterised by a vanishing acoustic pressure on the boundary: that is the Dirichlet condition. An approximate example is underwater acoustic incidence to a calm free surface (13). A pressure release boundary is also known as 'sound-soft' in literature (1, p. 13) and has the plane-wave pressure reflection coefficient of −1.
Imagine an inclined flat interface of y = 0, in Fig. 2, with an anticlockwise rotation angle β: note that β ± π leads to the same interface since the interface stretches to infinity. At a receiver b on the flat boundary, the acoustic fields scattered from a pair of real and image objects may be observed as follows from (3.4) and (3.5): and Then, p c (r bc ) + p c (r b c ) = 0 is expected for the Dirichlet condition. Under the method of images, the location of an image scatterer is restricted with regard to that of the corresponding real scatterer: r b c = r bc , θ b c + θ bc = 2nπ + 2β with integer n. Therefore, the sum of pressures at the point b becomes To have e imθ bc common in (3.8), only the terms associated with C c m have their index m changed in sign under the framework of ∞ m=−∞ : that is m with −m. Then, H −m (·) is converted back to H (1) m (·) by using the following property (25, p. 358) Therefore, to satisfy the Dirichlet condition under general circumstances, the expression inside the square bracket should vanish. This yields an image condition between the paired coefficients: in which the sign of m has been changed just for the purpose of display. Equation (3.10) do not appear to be available elsewhere. Two special cases of (3.10) can be discussed. For a x-axis boundary, setting β = 0 or π results in This condition was adopted by Ye and Chen (13): see their Eq. (5). Although these special interfaces are referred to as x-and y-axes for convenience of addressing them, what matters more is the relation between the polar angle and the interface. For an example of (x, y) ≡ r(sin θ, cos θ ), the roles of x-and y-axes are exchanged.
The difference between (3.11) and (3.12), or (3.10) in general, demonstrates that image conditions depend on the orientation of a flat interface relative to the 'global' coordinate system. This is because the coordinate system 'local' to scatterers are not rotated. Such rotation is not necessary since scatterers are circular.

Image condition: rigid interface
A rigid boundary can be characterised by a disappearing particle velocity normal to the boundary, which is equivalent to a vanishing pressure gradient, that is satisfying the Neumann condition. An acceptable example can be the ground during outdoor sound propagation (10,11). A rigid boundary is also known as 'sound-hard' (1, p. 13) and has the plane-wave pressure reflection coefficient of +1.
Therefore, image conditions for rigid flat boundaries can be derived by enforcing the Neumann condition, that is applying the derivatives of (3.6) and (3.7) at an observation point b with respect to the normal to the interface. However, image conditions for a rigid boundary can be obtained more easily by exploiting the phenomenon that the acoustic fields in real and image domains are perfectly replicated with symmetry. Therefore, one can expect p c (r bc ) = p c (r b c ), which means that image conditions of rigid flat interfaces will be of opposite sign to those of corresponding pressure release boundaries.
For a rigid flat interface of y = 0, in Fig. 2, with an anticlockwise rotation angle β, the following image condition is obtained: by taking an opposite sign to (3.10). Equation (3.13) may not have been documented elsewhere. Two special cases can be obtained: a flat interface being either x-or y-axis. For a x-axis boundary, β = 0 or π makes the following image condition: (3.14) This condition was known to Romero-García et al. (9): see their Eq. (27). For a y-axis boundary, which is the same as Eq. (9) by Lui and Li (18).

Method of images
In two dimensions, the method of images can be applied to construct acoustic domains other than a half-plane, regardless of using the image conditions. Acoustic fields in the concave part of a corner or acute-angled wedges can be analysed. However, there are certain restrictions in realising a wedgeshaped acoustic domain through the method of images (17,26). Having a finite number of image scatterers requires a wedge angle of β = π/Q where Q is the number of total wedges in a given halfplane; hence, 2Q wedges in total. When the wedge apex is placed at the global origin O o , images of a cth scatterer are constructed according to r q co = r 1 co for integer q (2 ≤ q ≤ 2Q); θ q co = (q − 1)β + θ 1 co for odd q; θ q co = qβ − θ 1 co for even q. Here, superscripts represent wedge numbers. For wedges with both rigid sides, Q can be any positive integer: Q = 1, 2 are equivalent to halfplane and corner, respectively. All component interfaces are of rigid as shown for an example of Q = 5 in Fig. 3a, where wedges are identified by numbers with the real wedge indicated by 1. Signs of all image wedges are the same as that of the real wedge. A rule of signs for wedges is that signs are toggled when broken lines (for pressure release interfaces) are crossed both between adjacent wedges and across any half-plane boundary, but are kept unchanged across continuous lines (for rigid boundaries).
For wedges with both pressure release sides, Q can also be any positive integer. All component flat boundaries are of pressure release; thus, signs are alternating as demonstrated in Fig. 3b: e.g. + for odd q and − for even q when +sign is assigned to the real first wedge.
However, there appears to be further restrictions for wedges mixed with rigid and pressure release sides. First, Q needs to be even. Then, the pressure release and rigid interfaces are placed as illustrated in Fig. 3d: that is (2j − 1)th and 2jth wedges have +sign for odd j; −sign for even j with 1 ≤ j ≤ Q, when +sign is assigned to the real first wedge. There is no contradiction in signs assigned to all  Fig. 3c, signs of bottom-half wedges cannot be determined uniquely. For example, for the sixth wedge to be an image of the real first wedge, its sign needs to be +; but to be an image of the fifth wedge, its sign must be −. For wedges and corners, image conditions derived for single flat interfaces can be combined to construct image conditions customised for them. Single-boundary image conditions in (3.10) and (3.13) are applied whenever their corresponding interfaces are crossed. Alternatively, the rigid boundary image condition of (3.13) is invoked across each single interface; then, its sign is toggled if the boundary is of pressure release.

Wedge image conditions
Without loss of generality, imagine a wedge-shaped domain, one of which side is always the positive x-axis and its real domain lies in y > 0. Further, confine the discussion to wedges with all rigid sides for the time being. Then the image condition of (3.13) can be applied successively, first for an interface with an angle of β(= π/Q), then 2β, 3β till (2Q − 1)β. For instance, up to q = 4: Then, partial sums of alternating series are employed: for example, This leads to image conditions for wedges with arbitrary Q, collectively expressed by where q identifies various image wedges from 2 to 2Q. For the last image wedge of q = 2Q, (4.2) becomes (3.14). In addition, for q = Q and even Q, the corresponding image condition is equal to (3.15). Equation (4.2) appears to be novel. For the other two types of wedges, the sign of (4.2) needs adjusting as discussed in section 4.1 and Fig. 3. These signs are also applied to the strength (B s ) of image sources.

Corner image conditions.
A right-angled corner is certainly an interesting subset of wedges. Corner image conditions for three image quadrants may be expressed by setting β = π/2 in (4.1) or (4.2): The ± sign is determined according to: quadrants (2 nd , 3 rd , 4 th ) = (+, +, +) for all rigid sides; (−, +, −) for all pressure release sides; (+, −, −) for +x-axis pressure release and +y-axis rigid sides. Equation (4.3) does not appear to be reported elsewhere. Hasheminejad and Alibakhshi (16) evaluated scattering of a single circular cylinder near a rigid 2D corner. They did not however make use of image conditions.

Multiple scattering model
Formulations of multiple scattering models, except the use of image conditions, are identical between free-space models and the ones with the method of images. These free-space formulations are well established, especially for homogeneous circular scatterers; therefore, derivations are not provided here. Detailed derivations can be found in other works (1, 5, 6, 9).

System of linear equations in two dimensions
When image entities are explicitly expressed for the method of images, the system of linear equations to be solved in two dimensions may be written as follows: s : see Fig. 3 for example. The propagation angle α of plane-wave incidence is defined by k α x = k cos α and k α y = k sin α, where k α x and k α y are the x-and y-components of a wavenumber k for the incidence. Since the origin of a plane-wave source cannot be known, I α c ≡ e ikr co cos(θ co −α) is required to take into account the position of a scatterer relative to others. This term is often known as 'phase factor' in literature (5). Furthermore, α q is the qth-image wedge counterpart of α 1 : α q = (q − 1)β + α 1 for odd q; α q = qβ − α 1 for even q.
The characteristics of fluid scatterers are embodied in (5) .
Here, J n (z) = (d/dz)J n (z) and H n (ka c )/J n (ka c ) (5). The inverse of W c n can be seen as the diagonal components of T-matrix for a circular scatterer c in isolation (1).
The particular form of the simultaneous equations (5.1) is the outcome of dividing the original set of equations by W c n,D , which is the denominator of W c n , to present (5.1) concisely. Numerically, however, the system can become ill-conditioned; therefore, it is advisable to return to the original equations when it comes to numerical implementation. Especially for either rigid or pressure release scatterers, W c n,D can vanish.

Computational resources
Without use of image conditions, the system of linear equations in (5.1) has 2QN R (2M + 1) unknowns when a constant M is assumed across scatterers for simplicity of presentation. With use of image conditions, however, the number of unknowns is reduced to N R (2M + 1) because C jq m is expressed by C j1 m . This follows that the memory requirement for the matrix in (5.1) is reduced by 4Q 2 -fold for direct solvers; the computation time for the matrix speeds up by 2Q times.

Field evaluation
Once C j1 m is resolved from (5.1), the total exterior sound field p E at r = (x, y) in the real first wedge can be calculated by for r c a c regardless of the type of sources. Note that, except C c1 m , the superscript 1 is implicit or unnecessary. The appearance of 2i/πka c is a result of applying the Wronskian for Bessel and Hankel functions. For rigid and pressure release scatterers, (5.4) is not applicable.

Spherical-wave incidence.
So far, the current article has addressed 2D problems where the scatterers are circles and hence a concentrated source has a circular wavefront. The formulations can be equally applicable to 3D problems with infinitely long circular cylinders with uniform cross sections and line sources in which both axes from scatterers and a source are parallel. Such arrangement can be imagined for horizontal sonic crystals alongside the road with continuous presence of traffic. More practically, however, the incidence of a spherical wave can also be simulated by extending the multiple scattering formulations addressed in this article. This has been presented elsewhere for a single cylinder (4,18); also for an array of cylinders (23). Under the method of images, the same image conditions derived for 2D configurations can also be used as reported for a single cylinder parallel to the rigid half-space (18).
Equation (3.3) suggests that a spherical wave can be constructed through the inverse Fourier transform of cylindrical waves. Therefore, one can also use the pair of Fourier transform (23) for multiple scattering. Then, equations from (5.1) to (5.4) can be treated as if they were the results of the forward Fourier transform along the z-dimension. To that effect, every single k in these equations would be replaced by k r = (k 2 − k 2 z ) 1/2 ; every single κ c by (κ 2 c − k 2 z ) 1/2 withk = ω/v andκ c = ω/v c . Then, the inverse transform can be performed to (5.3) and (5.4) to obtain the response under the spherical-wave incidence. The exterior fields, for example, can be evaluated as follows, with the upper alternative of (5.3): where a superscript 3D is used to represent three dimensions. Image conditions are still applicable for the procedure leading to p E (r; |k z |), in which |k z | highlights it is an even function of k z . Within (5.5), the source-field contribution can be calculated analytically by using (3.3). For the remaining scattered contributions, (5.5) can be evaluated by a direct numerical integration. Note, below, that I n (·) and K n (·) are the modified Bessel functions of the first and second kinds, respectively, with order n. For |k z | >k, H (1) n (k r d) in (5.1), (5.2) and (5.3) can be expressed in terms of K n (d(k 2 z − k 2 ) 1/2 ), where d = a c , r cs , r c , r j , etc. Since K n (x) decreases monotonically as x → ∞ and it is a c < r cs , r cj , r c for exterior problems, the combination of H (1) n (k r d) for various d yields overall attenuation as k z → ±∞ due to the specific position of their respective H (1) n (·). At the same time, J n (k r a c ) in (5.2) behaves like I n (a c (k 2 z −k 2 ) 1/2 ) which is diverging as k z → ±∞; however, their effect is suppressed by other H (1) n (k r d). Therefore, the integration limit can be truncated with good accuracy. Standalone, H (1) n (x) becomes singular as x → 0, which occurs while |k z | →k. This can pose a difficulty in principle; however, Li and Ueda (4) showed, for a single circular cylinder, the integrand of (5.5) reaches a limiting value in the vicinity of k z = ±k for n = 0 and has a logarithmic singularity for n = 0. We assume this can also be the case for multiple cylinders; hence, the direct numerical integration can be made.
For infinitely long cylinders, an additional perfectly reflecting plane perpendicular to their axes can be introduced under the method of images. Without loss of generality, the plane can be set to have z = 0: that is the xy-plane. Image sources are introduced below z = 0 when the real domain belongs to z > 0. There is no need to add image cylinders across the new plane, since the existing ones are of infinite length already covering the image domain; the effect will be equivalent to the length of real cylinders being cut half: that is semi-infinite. Possible configurations are shown in Fig. 1, among which Fig. 1a was demonstrated elsewhere (23) for a rigid ground without vertical flat boundaries. Then, the acoustic response outside the cylinders may be obtained by where z s = z − z so , z s = z − z so and z so = −z so . The +sign in ± indicates the xy-plane to be rigid which corresponds to the alternative with two cos(·) functions; the −sign, to be pressure release, with two sin(·) functions. It is noted that image conditions are not defined for this additional plane perpendicular to the axes of cylinders. For example, Fig. 1a has no image condition defined. Across the light-grey planes, Fig. 1b has half-plane image conditions such as (3.14) and (3.15) among others; Fig. 1c has 2D corner image conditions in (4.3); Fig. 1d has wedge image conditions in (4.2) with adjustment of wedge signs.

5.4.2
Plane-wave incidence in three dimensions. The total acoustic field under plane-wave incidence in 3D can also be evaluated by using (5.3) and its associated equations. Physically, a plane-wave incidence can be emulated by a spherical-wave incidence from a faraway source; however, such adaptation is mathematically rather inefficient. In fact, oblique incidence of plane waves can be considered without the integration encountered in (5.5) and (5.6). Suppose that a plane wave propagates with angles of γ and α 1 where k x =k sin γ cos α 1 , k y =k sin γ sin α 1 and k z =k cos γ . The angle α 1 has the same role as in (5.1) and (5.3); γ is independent of wedge locations as long as the xy-plane is not introduced. First, one can replace every single k in (5.3) and others by k r = (k 2 − k 2 z ) 1/2 =k sin γ ; every single κ c by (κ 2 c − k 2 z ) 1/2 . Secondly, one needs to multiply exp(ik z z) to (5.3) (3): p 3D E (R) = p E (r; |k z |)e ik z z ; therefore, extra computational cost incurred is marginal, in comparison with 2D solutions for circles. If a perfectly reflecting xy-plane is introduced, the exterior response of semi-infinitely long cylinders can be evaluated by 2i sin(kz cos γ ) , (5.8) where γ represents image plane waves from z < 0 and γ + γ = π has been used. The convention for the use of ± is as before; the alternative with cos(·) is for a rigid xy-plane, the one with sin(·) is for pressure release.

Other separation-of-variables models
Equations (3.4) and (3.5) are not the only ways to model outgoing or scattered acoustic waves in two dimensions. There are other variants which may have implications for image conditions. For example, other works (5, 10, 11) have rescaled the unknown coefficients as C c m = A c m Z c m by using (5), the image conditions for A c m are identical to those for C c m . Furthermore, the radiation problems from directional finite-sized source in a circular shape are also subject to the same image conditions. 5.5.1 Type of scatterers. Equations (5.1) and (5.2) are mainly concerned about homogeneous fluid-like scatterers. There is a free-space formulation addressing thin-walled circular scatterers with slits on the circumference (27). Since their outgoing scattered fields obey (3.4), their model is also expected to follow the same image conditions. An advantage of using image conditions, in addition to the savings in resources, is that image scatterers do not need rotation to mirror the position of slits; such rotation is however required if image conditions are not used.
In fact, image conditions can be applied regardless of scatterer types as long as scatterers are circular-shaped and outgoing acoustic fields are modelled by (3.4). This is because, during the derivation of image conditions, the boundary conditions of individual scatterers were not required. The only boundary conditions necessary during the derivation were the ones for perfectly reflecting flat interfaces. For fluid-like scatterers, their boundary conditions are taken into account through (5.2) when the system of linear equations is formed in (5.1). where E c m is an unknown coefficient. Unlike (3.9), the modified Bessel function of the second kind is subject to K −m (z) = K m (z) for an integer order m (25, p. 375). Therefore, their image conditions are derived as

Modified Helmholtz equation.
with the aforementioned usage for ± sign. Furthermore, corresponding wedge image conditions become for rigid wedges. For mixed-boundary and pressure release wedges, the sign of wedges needs adjusting.

Verification
It is straightforward to numerically verify the image conditions. Two runs of multiple scattering calculations for the same configuration with and without using them will confirm the veracity of the image conditions. Presenting such comparisons, however, does not have merit since two outputs will be numerically identical around the order of a machine epsilon and will not be distinguishable in a printed material. It is simply confirmed that was the case for our test runs. Instead, the multiple scattering model (MSM) is compared with measurements for half-and quarter-space problems. Then, numerical comparisons between the MSM and boundary element method (BEM) are presented for 2D and 3D corners, and finally for a benchmark wedge geometry established by the Acoustical Society of America.
For display of results, a measure of 'pressure loss' (PL) is adopted: where p emp represents the field of empty flat interfaces; p tot , the total pressure level with scatterers present.

Measurements
Two sets of pressure-loss measurements were conducted in an anechoic chamber: half-and quarterspace problems. Five identical plastic circular pipes with 55 mm in nominal diameter were placed parallel to one another in contact with a horizontal board; for a corner, a vertical board was also erected. An omni-directional source and a microphone was set up above the horizontal board: Brüel & Kjaer models 4295 and 4189 respectively. The source-receiver axis was aligned to cross perpendicularly the midpoint of 2-m long cylinders. To measure the pressure loss, the acoustic responses of the bare boards were recorded separately from those of the boards populated with the cylinders by using the maximum-length sequence signals which were generated and acquired by using a National Instruments data acquisition board controlled by a Matlab ® toolbox for data acquisition. The source signal was amplified through an ordinary audio amplifier; the receiver signal was conditioned via Brüel & Kjaer model 2636.
For the MSM, the boards and cylinders were all treated as rigid. Although the circular cylinders were of finite length, they are assumed to be infinitely long under the spherical-wave incidence during the MSM. Direct numerical integration of (5.5) has been performed, during which the image conditions have been adopted also. 6.1.1 Half-space measurement. For the MSM, one must choose a coordinate orientation relative to the flat interface. When the horizontal board is assumed to be along the x-axis, the source and receiver were placed at (x, y) = (0, 8) cm and (95, 5.95) cm respectively, with the same z-coordinate. The x-coordinates of the cylinder centres were from 12.75 to 92.75 cm in 20-cm spacing. An estimated speed of sound was 343.3 m s −1 based on temperature reading.
For this configuration of a x-axis (that is zx-plane) boundary, a correct image condition for the MSM is (3.14) which results in the black line in Fig. 4. The measured PL is presented by the grey curve. The level of the agreement is convincing although the measurement was carried out for finite-length circular cylinders in contrast to infinitely long ones assumed for the MSM. The measured pressure loss for the corner is drawn in the grey line in Fig. 5. The MSM result using image conditions (4.3) for the corner is shown in the black line. It is demonstrated that the use of image conditions for a corner can lead to reasonable agreement with the measured pressure loss for a wide range of frequencies.
6.2 Simulations 6.2.1 2D corner. Numerical investigation has been carried out in the presence of two real circular scatterers of the same radius a near a 2D corner. A plane wave propagates towards the corner with 45 • incidence angle. The outer medium is water, with wavenumber k w , having the density and sound speed of 999 kg m −3 and 1482 m s −1 , respectively. A scatterer whose centre is located at (2a, 4a) is chosen as air (1.2 kg m −3 , 343 m s −1 ) with wavenumber k a . Another scatterer is placed at (4a, 2a) with the properties of mercury (13534 kg m −3 , 1451 m s −1 ) as an example of a harder scatterer. The positive x-axis side of the corner is set to be pressure release; the positive y-axis is rigid.
The MSM was run by using image conditions for a mixed-boundary corner shown in (4.3). For the BEM runs, between the outer medium and scatterers, multi-domain transmission problems were established. Both sides of the corner were also discretised. For a reference field for the pressure loss, an empty corner was also calculated by the BEM. The rigid side of the corner had the admittance of zero, while the pressure release side had the admittance of from-water-to-air incidence. Therefore, this pressure release interface is not exact, but a good approximation: for example, under normal incidence, the plane-wave pressure reflection coefficient reaches −0.9994.  Fig. 6, the comparison of pressure loss between the MSM and BEM is shown for a field point (3.25a, 2a) inside the mercury scatterer: the agreement is good. There are several abrupt changes. The one at the lowest frequency is related to a pulsating mode of the air scatterer. The rest, starting from k a a ≈ 1.84, occur because inside the air pocket its air-to-water boundary is near rigid. These resonant frequencies closely match the k a a values for the extrema of various orders of Bessel function of the first kind, although some resonances are not detected in this field point and graphical resolution. 6.2.2 3D corner. Two circular cylinders (with height) have been placed in a 3D corner, as shown schematically in Fig. 1c. The origin of the corner is (0, 0, 0) in Cartesian coordinates; the real domain is the first octant. The centres of two cylinders are at (x, y) = (40, 40) and (50, 50) cm. For the MSM, the axes of two cylinders are semi-infinite along the z-axis; for the BEM, the heights are only 2 m. The diameters are 5.5 cm for both modalities. A point monopole source is located at (75, 75, 50) cm. The medium is air. The cylinders and corner are treated as rigid.
For the MSM, 2D-corner image conditions in (4.3) have been adopted for the light grey sides in Fig. 1c. For the dark grey side perpendicular to cylinders, (5.6) with its upper alternative has been used. For the BEM, a 3D model has been implemented with the surface of cylinders discretised by using a third-party software Gmsh © (28). The surface of the corner was not meshed; therefore, the method of images has been explored with a 3D-corner Green's function. Figure 7 shows the comparison between the MSM and 3D-BEM. Frequency spectra of pressure loss are presented in Fig. 7a  two cylinders: the agreement is good overall. The number of triangular meshes was adjusted over frequency: from 2120 at 1 kHz to 21452 at 5 kHz per cylinder. Figure 7b demonstrates the variation of pressure loss along the height at (x, y) = (45, 45) cm for 2.1 kHz with 21452 elements per cylinder. Both modalities match acceptably below the 2-m height, above which discrepancy is observed. Such deviation is anticipated since the BEM employed 2-m high cylinders while the height is infinite in the MSM.

Wedge
The Acoustical Society of America identified benchmark configurations for ocean acoustic problems (14,15). A wedge-shaped waveguide with homogeneous water (1000 kg m −3 , 1500 m s −1 ) was established with a pressure release flat sea surface on top and three different boundary conditions for a sloping sea floor. A pressure release bottom was considered, but a rigid sea floor was not part of the benchmark problem. The wedge has 4-km horizontal range from its apex and 200-m depth from the sea surface with the wedge angle of approximately 2.863 • .
A multiple scattering problem within this benchmark wedge is demonstrated by using image conditions. However, the geometry of the wedge is modelled approximately. It is approximate only because the method presented in this article cannot address exactly the wedge angle of arctan(1/20). A rigid condition is chosen instead for the sloping boundary. For a wedge with both pressure release and rigid sides, an acceptable wedge angle is either β ≈ 2.903 • by Q = 62 or β ≈ 2.813 • by Q = 64, since Q = 63 does not meet the criterion addressed earlier in section 4.1: see Fig. 3c for odd Q. Therefore, Q = 64 is chosen for the ensuing calculations, which results in the relative error of 1.74% for β. The water depth at the 4-km range is now 196.5 m. It is emphasised that, despite the limited choice of a wedge angle, our model is capable of including circular objects in the wedge interior, unlike other arbitrarily-angled wedge formulations (14,15).
As a scatterer, a rigid circle with 10 m in diameter was chosen to approximate the cross section of, say, a large submarine. Note that, despite a single real scatterer, it is a multiple scattering problem involving 128 scatterers due to the method of images. Figure 8 shows the comparison of pressure loss between the MSM and BEM due to the presence of a scatterer centred at (range, depth) = (3.8, 0.1) km and a line source at (4, 0.1) km with the wedge apex at (0, 0). Range dependence in every 4 m is demonstrated in Fig. 8a, for which each side of the wedge was set to 5 km in length for the BEM and discretised. The BEM with a 2D free-space Green's function was run twice: with and without a circle.
Frequency dependence in every 5 Hz is illustrated in Fig. 8b, for which the BEM was based on the method of images by using a wedge Green's function. The truncation order M was increased gradually from 4 at 25 Hz to 20 at 500 Hz. Good agreement is found in Fig. 8. Non-zero pressure loss in dB indicates that the acoustic field is disturbed by the introduction of a scatterer.
The data for Fig. 8 can be obtained without using wedge image conditions. However, a penalty can be severe especially for higher frequencies and a large number of scatterers. Due to 2Q = 128, the size of the matrix in (5.1) is 16384-times larger. Building the system of equations requires 128-times more computations. For example, by using image conditions, the 500-Hz data point in Fig. 8b for a single (real) scatterer requires the memory space of about 25 kB for a double-precision complexvalued matrix, but would take as much as 400 MB without the image conditions. When the size of the linear system becomes too large, direct solvers may not be efficient; hence iterative solvers can be more suitable (6, 7). It is noted however that the use of our image conditions will aid the computation regardless of the type of solvers since the size of the problem is reduced in the first place.

Incorrect image conditions
It may be informative to discuss the incorrect image conditions reported elsewhere. van der Aa and Forssén (10) investigated a large number of scatterers up to 1128 cylinders above a rigid ground, and hence expressed need for more efficient algorithm. However, they cited (3.15), while the orientation of their coordinate system was compatible only with (3.14). In their follow-up publication (11), an image condition was not mentioned at all, despite a high workload that 15 design parameters of sonic crystals were optimised by using a multi-objective genetic algorithm for arrays of up to 141 cylinders above a rigid ground. This may suggest that an image condition correct for their configuration was not known to them. Hasheminejad and Alibakhshi (29) evaluated the scattering of a poroelastic circle above a flat interface. Hasheminejad and Azarpeyvand (30) studied the vibrations of a circular radiator. Qiao et al. (31) addressed the radiation force on a circular cylinder near a flat boundary. However, their image condition was different from those addressed in the current article and elsewhere (9, 12, 13, 18). Instead, their image condition appears to have been adapted from the axisymmetric case of a spherical-coordinate image condition (for a sphere) in an earlier publication by Hasheminejad (32): across a rigid flat interface. Otherwise, no derivation was given in (29)(30)(31). In fact, their configuration conforms to (3.12) and (3.15) for pressure release and rigid boundaries, respectively. Results of using image conditions for five circular cylinders on a rigid plate. Continuous line, using a correct image condition; dashed line, using an image condition wrong for the coordinate orientation; dotted line, using a spherical coordinate image condition. The frequency resolution is 20 Hz Figure 9 shows the comparison of using a correct image condition and two wrong ones for a given choice of the coordinate orientation. The geometric configuration and material properties identical to those for Fig. 4 are revisited. When a rigid flat boundary is chosen to be the x-axis (that is zx-plane) as for Fig. 4, the MSM result using a correct image condition of (3.14) is shown in the continuous line in Fig. 9: this curve is identical to the simulation in Fig. 4. The outcome by an image condition of (3.15) for a y-axis (that is yz-plane) interface is drawn in the dashed line. Use of (7.1) from an axisymmetric spherical coordinate image condition results in the dotted line. Stark discrepancy is observed between the continuous curve and the other two discontinuous. Figure 4 has demonstrated that the continuous curve in Fig. 9 is in good agreement with the measurement. Therefore, caution is warranted to avoid wrong choice of image conditions. Especially, it is not surprising that spherical coordinate image conditions should not be used in place of polar coordinate image conditions due to the difference in their respective wavefunctions. More details of spherical coordinate image conditions can be found in (17).

Conclusions
Image conditions proposed in this article enable efficient implementation of the method of images for polar and cylindrical coordinate multiple scattering models (MSM), based on the method of separation of variables approach and Graf's addition theorem. Scatterers can be circles in polar coordinates and (semi-)infinitely long circular cylinders in cylindrical coordinates. First, image conditions were derived for a single flat boundary of either pressure release or rigid nature. Then, these were combined to produce image conditions for wedges and corners. For wedge-shaped acoustic domains with a wedge angle of π/Q rad with an integer Q, using wedge image conditions reduces computation of the system of linear equations by 2Q times and memory requirements for a matrix by 4Q 2 times compared to the approach without using the image conditions.
The MSM with image conditions was validated through measurements for the configurations of half-space and a 2D corner in an anechoic chamber. Numerical simulations were conducted for a pressure release and rigid 2D corner and a rigid 3D corner, both of which were confirmed by 2D and 3D boundary element method (BEM), respectively. A benchmark configuration for a wedgeshaped ocean environment with Q = 64 was also evaluated by using wedge image conditions and subsequently verified by the BEM.
Awareness was raised of the danger of using wrong image conditions reported elsewhere. Image conditions are sensitive to the orientation of the coordinate system relative to flat interfaces of interest. It is wrong to use a spherical coordinate image condition (for spheres) to polar coordinate multiple scattering problems (for circles).
The MSM was initially confined to free-space problems. Use of the method of images extends its appeal to wider applications. Image conditions from this article further improve its numerical efficiency.