On the eigenmodes of periodic orbits for multiple scattering problems in 2D

TL;DR

This paper develops an asymptotic approximation for the phases of densities in periodic orbits of multiple scattering problems in 2D, enabling faster ray tracing independent of frequency.

Contribution

It introduces a Taylor series-based asymptotic approximation for phases in periodic orbits, applicable to multiple obstacles, improving computational efficiency in high-frequency wave problems.

Findings

Phase coefficients are independent of wavenumber and incident wave.

The method accelerates ray tracing schemes after initial iterations.

Explicit algorithms are provided for arbitrary 2D obstacles.

Abstract

Wave propagation and acoustic scattering problems require vast computational resources to be solved accurately at high frequencies. Asymptotic methods can make this cost potentially frequency independent by explicitly extracting the oscillatory properties of the solution. However, the high-frequency wave pattern becomes very complicated in the presence of multiple scattering obstacles. We consider a boundary integral equation formulation of the Helmholtz equation in two dimensions involving several obstacles, for which ray tracing schemes have been previously proposed. The existing analysis of ray tracing schemes focuses on periodic orbits between a subset of the obstacles. One observes that the densities on each of the obstacles converge to an equilibrium after a few iterations. In this paper we present an asymptotic approximation of the phases of those densities in equilibrium, in the form of a Taylor series. The densities represent a full cycle of reflections in a periodic orbit. We initially exploit symmetry in the case of two circular scatterers, but also provide an explicit algorithm for an arbitrary number of general 2D obstacles. The coefficients, as well as the time to compute them, are independent of the wavenumber and of the incident wave. The results may be used to accelerate ray tracing schemes after a small number of initial iterations.

Equations

$$u^{\text{s}}(\mathbf{x})={\int }_{\Gamma }K(\mathbf{x},\mathbf{y})v(\mathbf{y})\mathrm{d}s(\mathbf{y}).$$

$$v(\mathbf{y})=\sum _{m=1}^M{f}_m(\mathbf{y},k)e^{ik{g}_m(\mathbf{y})}.$$

$$V_j(\tau )=v_j^{\text{smooth}}(\tau )e^{ik{\varphi }_j(\tau )},j=1,\dots ,J,$$

$$L(\mathcal{T})=\Vert {\Gamma }_J({\tau }_J^{\ast })-{\Gamma }_1({\tau }_1^{\ast })\Vert +\sum _{j=1}^{J-1}\Vert {\Gamma }_j({\tau }_j^{\ast })-{\Gamma }_{j+1}({\tau }_{j+1}^{\ast })\Vert =\sum _{j=1}^J\Vert {\Gamma }_j({\tau }_j^{\ast })-{\Gamma }_{j+1}({\tau }_{j+1}^{\ast })\Vert .$$

$${\Delta }_{i,j}({\tau }_i,{\tau }_j)=\Vert {\Gamma }_i({\tau }_i)-{\Gamma }_j({\tau }_j)\Vert ,1\le i,j\le J.$$

$$f(x;k)\sim u(x;k)e^{ikg(x)},k\gg 1,$$

$$(S_{ji}{V}_i)(\mathbf{x})={\int }_{{\Gamma }_i}K(\mathbf{x},\mathbf{y})V_i(\mathbf{y})\mathrm{d}s(\mathbf{y}),\mathbf{x}\in {\mathbf{\Gamma }}_{\mathbf{j}}.$$

$$\{\begin{array}{ll}{S}_{11}{m}_1+S_{12}{m}_2& =-u_1^{\text{inc}},\\ S_{21}{m}_1+S_{22}{m}_2& =-u_2^{\text{inc}}.\end{array}$$

$$T=S_{11}^{-1}{S}_{12}{S}_{22}^{-1}{S}_{21},$$

$$T{V}_1=\lambda V_1.$$

$$S_{21}{V}_1$$

$$H_0^{(1)}(z)\sim \frac{\sqrt{2}}{\sqrt{\pi z}}{e}^{i(z-\pi /4)},z\to \infty .$$

$$S_{21}{V}_1$$

$${\Delta }_{1,2}({\tau }_1,{\tau }_2)=\Vert {\Gamma }_2({\tau }_2)-{\Gamma }_1({\tau }_1)\Vert .$$

$$g_{2,1}({\tau }_2,{\tau }_1)={\Delta }_{1,2}({\tau }_1,{\tau }_2)+{\varphi }_1({\tau }_1).$$

$$\frac{\partial g_{2,1}}{\partial {\tau }_1}({\tau }_2,{\chi }_1({\tau }_2))=0.$$

$$g_{2,1}({\tau }_2,{\chi }_1({\tau }_2))={\Delta }_{1,2}({\chi }_1({\tau }_2),{\tau }_2)+{\varphi }_1({\chi }_1({\tau }_2)).$$

$${\chi }_1({\tau }_2^{\ast })={\tau }_1^{\ast }\text{mbox}and{\chi }_2({\tau }_1^{\ast })={\tau }_2^{\ast }.$$

$${\varphi }_1({\tau }_1)=g_{1,2}({\tau }_1,{\chi }_2({\tau }_1))={\Delta }_{2,1}({\chi }_2({\tau }_1),{\tau }_1)+{\varphi }_2({\chi }_2({\tau }_1)).$$

$$S_{11}{V}_1({\tau }_1)=u({\tau }_1)={\int }_{{\Gamma }_1}\frac{i}{4}{H}_0^{(1)}\left(k\Vert {\Gamma }_1({\tau }_1)-{\Gamma }_1(\tau )\Vert \right)\left(v_1^{\text{smooth}}(\tau )e^{ik{\varphi }_1(\tau )}\right)\Vert \nabla {\Gamma }_1(\tau )\Vert d\tau .$$

$$H_0^{(1)}(z)=\frac{1}{\pi i}{\int }_{-\infty -\pi i/2}^{+\infty +\pi i/2}{e}^{iz\cosh (t)}\mathrm{d}t.$$

$$\{\begin{array}{l}\frac{\partial g_{2,1}}{\partial {\tau }_1}({\tau }_2,{\chi }_1({\tau }_2))=0\\ {\varphi }_2({\tau }_2)={\Delta }_{1,2}({\chi }_1({\tau }_2),{\tau }_2)+{\varphi }_1({\chi }_1({\tau }_2))\\ \frac{\partial g_{1,2}}{\partial {\tau }_2}({\tau }_1,{\chi }_2({\tau }_1))=0\\ {\varphi }_1({\tau }_1)={\Delta }_{2,1}({\chi }_2({\tau }_1),{\tau }_1)+{\varphi }_2({\chi }_2({\tau }_1))-\mu .\end{array}$$

$$g_{1,2}({\tau }_1,{\tau }_2)={\Delta }_{2,1}({\tau }_2,{\tau }_1)+{\varphi }_2({\tau }_2).$$

$$\mu =L(\mathcal{T}),$$

$$T=S_{11}^{-1}{S}_{1J}\dots S_{33}^{-1}{S}_{32}{S}_{22}^{-1}{S}_{21}.$$

$$g_{j,i}({\tau }_j,{\tau }_i)={\Delta }_{i,j}({\tau }_i,{\tau }_j)+{\varphi }_i({\tau }_i),$$

$$\{\begin{array}{l}\frac{\partial g_{j+1,j}}{\partial {\tau }_j}({\tau }_{j+1},{\chi }_j({\tau }_{j+1}))=0,\\ {\varphi }_{j+1}({\tau }_{j+1})={\Delta }_{j,j+1}({\chi }_j({\tau }_{j+1}),{\tau }_{j+1})+{\varphi }_j({\chi }_j({\tau }_{j+1})),\\ \frac{\partial g_{1,J}}{\partial {\tau }_J}({\tau }_1,{\chi }_J({\tau }_1))=0,\\ {\varphi }_1({\tau }_1)={\Delta }_{J,1}({\chi }_J({\tau }_1),{\tau }_1)+{\varphi }_J({\chi }_J({\tau }_1))-\mu .\end{array}j=1,\dots ,J-1.$$

$$\mu =L(\mathcal{T}).$$

$${\chi }_j({\tau }_{j+1}^{\ast })={\tau }_j^{\ast },$$

$${\Gamma }_1(\tau )$$

Full text

Keywords: Boundary Element Method, Oscillatory integration, High-frequency scattering, Multiple obstacles, Periodic orbits, Phase extraction.

MSC: Primary: 65N38, Secondary : 45A05, 45M05, 65R20.

1 Introduction

Numerical simulations in acoustics are often based on a boundary integral equation reformulation of the Helmholtz equation, see for example [12, 32]. An incoming wave that is scattered by an obstacle $\Omega$ with boundary $\Gamma$ results in a scattered field $u^{\text{s}}(\mathbf{x})$, that can be represented by the so-called single layer potential,

[TABLE]

Here, $K(\mathbf{x},\mathbf{y})$ is the Green’s function of the Helmholtz equation with wavenumber $k$ and $v(\mathbf{y})$ is the unknown density function, defined on $\Gamma$. Sound-soft and sound-hard obstacles give rise to Dirichlet or Neumann boundary conditions respectively, with zero pressure or zero normal velocity on $\Gamma$. These lead to integral equations of the first or second kinds. Though other representations of the scattered field exist, which give rise to different integral equations, in this paper we present our results using the representation above. The analytical results depend solely on the parameters of the physical problem, such as the geometry and relative locations of the scattering obstacles, and not on the particular integral equation that is used.

The numerical discretisation of, for example, the Dirichlet problem $u^{\text{s}}(\mathbf{x})=f(\mathbf{x})$, $\mathbf{x}\in\Gamma$, leads to a large and dense linear system [5]. A family of hybrid numerical-asymptotic methods aims to significantly reduce the size of the linear system by incorporating information about the solution from asymptotic analysis (see e.g. [1, 8, 3, 4, 19], as well as the review [10] and references therein). In particular, phase-extraction methods use information about the phase $g$ (or multiple phases $g_{m}$) of the solution $v(\mathbf{y})$, in order to discretize only the remaining non-oscillatory parts $f_{m}$ in the following factorization:

[TABLE]

The phases $g_{m}$ can be determined from asymptotic analysis, for example via ray tracing. Ray tracing methods for scattering configurations with multiple obstacles at high frequencies lead to a sequence of single-scattering problems, mimicking exactly the tracing of wavefronts, in such a way that the phase can be extracted from a single density on a single obstacle in each case [17, 16, 2, 11]. The methodology is illustrated in Figure 1 for the case of two circular obstacles. The first computation is the scattering of the incoming wave by just one of the obstacles, disregarding the other one. In this case, the density mostly inherits the phase of the incoming wave, evaluated along the boundary. Next, the reflected field is considered as an incoming wave for the other obstacle, and so on. This process is repeated and care has to be taken to include all reflections by all obstacles until convergence.

It was shown by Ecevit et al. in a series of papers that for a set of convex obstacles the sequence of single-scattering problems can be organised according to certain periodic orbits [16, 2, 14, 15, 7]. Examples of these are illustrated in Figure 2. Each periodic orbit consists of the shortest path between a subset of obstacles, and one can intuitively understand them as follows. Rays originate in a source and reflect off obstacles in the scene. Rays that do not reflect onto another obstacle leave the scene forever. For rays that do reflect again, one can repeat the reasoning on the next iteration of reflections. Fewer and fewer rays will remain, and those that do can only be close to the aforementioned periodic orbits. Indeed, a ray that travels exactly on the shortest path between a collection of obstacles remains trapped in the scene forever, while any other ray that deviates from these particular trajectories ultimately leaves the scene.

The rays in periodic orbits settle down to an asymptotic limit – asymptotic both in the number of iterations and in the frequency – which we call a mode of the configuration of scatterers. From experiments, one observes that the phases of the densities involved converge to some equilibrium after a number of iterations. Indeed, this can be seen in for example [17, Fig. 1 & 2], although this feature was not further investigated there.

Say a periodic orbit involves $J$ obstacles $\Omega_{j}$, $j=1,\ldots,J$. Rays that travel along the orbit induce a density $V_{j}$ on each of the participating obstacles. Starting with density $V_{1}$ on the first obstacle, the observation is that after a full cycle of reflections the induced density on the first obstacle is (asymptotically) proportional to $V_{1}$ again. We will assume the densities to have the form

[TABLE]

where the phase $\phi_{j}(\tau)$ corresponds to the conjectured limiting phase described above. We expect that this phase is only well-defined for $\tau$ corresponding to the portions visible to the previous scatterer in the periodic orbit, but $v_{j}^{\text{smooth}}(\tau)$ is very small for other $\tau$ anyway.

The ‘modes’ can be seen as eigenfunctions of the multiple scattering problem. Indeed, the existence of a limit indicates that a wave with the density $V_{1}$ on the first obstacle induces a density $V_{2}$ on the second obstacle, $V_{3}$ on the third obstacle and so on, such that the last density $V_{J}$ induces precisely the density $V_{1}$ again on the first obstacle, up to a constant factor. This constant factor can be seen as an eigenvalue, and the collection of densities on each of the obstacles is an eigenfunction of the sequence of single scattering problems. It only depends on the geometry of the scatterers and not on the incident wave.

Similar periodic or short orbits have appeared in other settings, in connection to Random Matrix Theory [25, 21]. Certain quantum graphs are employed to study disordered systems. In principle, rays can keep bouncing from an obstacle to another arbitrary obstacle and this wave chaos leads to the application of Random Matrix Theory to obtain spectral statistics of those graphs. Although wave chaos theory has been developed primarily in the context of quantum mechanics, the similarities with statistical methods such as Statistical Energy Analysis are explained in [30]. The case of an array of cylinders has been studied in [26] and can be viewed in analogy with the case of two disks studied here.

The goal of the current paper is to compute those limiting phases, assuming their existence. We formulate a highly nonlinear system of differential equations that describes the phases. We show that our goal can be achieved if we restrict ourselves to the computation of Taylor series expansions of the phases around the points where the periodic orbits touch the obstacles. This computation is independent of any incident wave, and it is also computationally independent of the wavenumber $k$ of the overall Helmholtz problem and allows solving the system numerically in a simple setting. Still, having an approximation close to the periodic orbit may already be sufficient for practical puposes, as the modulus of $V_{j}(\tau)$ is maximal for $\tau$ corresponding to rays near the periodic orbit. Physically, this is due to rays rapidly leaving the scene the further away they are from the periodic orbit.

Our implementation of the computation of Taylor series and other results in this paper is publicly available on GitHub [29] and a copy of the code can always be obtained from the authors.

Firstly, we exploit symmetry for the simple case of two circular scatterers, as shown in Figure 3. In this case, which was also presented in [22], there is a single periodic orbit, namely the shortest path between the two disks. Using symmetry, it is sufficient to study a single phase function. Secondly, a general procedure is presented to approximate the phase of densities associated with periodic orbits with arbitrary numbers of obstacles including quite general geometries. We analyze the limiting phase via its Taylor series around the point where the shortest path hits the boundary. The advantage of the Taylor series approximation is that it leads to fully explicit expressions.

We formulate the problem as an eigenvalue problem for an oscillatory integral operator that represents scattering. Looking for eigenfunctions of the form 2 leads to oscillatory integrals. In order to analyze these asymptotically for a large wavenumber $k$, we use the well-known method of steepest descent [31, 6]. It deforms the path of integration into paths of steepest descent in the complex plane, such that oscillating integrands are transformed to rapidly decaying integrands. A considerable simplification is that we do not need to determine these paths exactly, as we are only interested in the phase of these oscillatory integrals. We do remark that a more extensive analysis of the integrals could lead to an asymptotic description of the amplitudes $v_{j}^{\text{smooth}}(\tau)$ in 2 as well as the phase.

A similar methodology was employed in [20] in order to obtain the asymptotic expansion of the solution of a scattering problem, possibly involving multiple obstacles as well. In particular, the method of steepest descent was used in [20] to track the forward propagation of rays, given an initial incoming wave. The main difference is that in our current setting, we are concerned with an eigenvalue problem instead.

The paper is structured as follows. First, we detail the formulation of a non-linear system of equations that describes the phases on two obstacles, later generalised to multiple obstacles. In § 3, we derive explicit Taylor expansions for the phases of the mode of a periodic orbit with two disks. In the next section this is again generalized to a larger number of obstacles with general shapes. A geometrical interpretation of the phase is given in § 5. The next section details the numerical approximation of the phase that we use to validate our results in § 7. Finally, we present some concluding remarks.

Remark 1

In this paper we compute the asymptotic phases of the densities in periodic orbits, under the assumption that these orbits exist. This is the case at least when the configuration of obstacles satisfies the conditions stated by Ecevit et al. in [16], where 2D multiple scattering problems involving multiple obstacles were first analyzed. Thus, following [16], we assume a set of $J$ disjoint, smooth and convex obstacles, where in addition rays between two obstacles never intersect a third obstacle (the no-occlusion condition). Yet, we show that the expressions do seem to be valid for some more general configurations in § 7.

2 Problem statement

In this section, we arrive at the formulation of a non-linear system of equations that describes the phases on each obstacle. We focus on $J=2$ scatterers and generalise this to an arbitrary number at the end of this section.

2.1 Notation

As in 2, we will refer to functions on obstacle $\Omega_{j}$ with subscript $j$, where $j=1,\ldots,J$. Each obstacle $\Omega_{j}$ has a boundary $\Gamma_{j}$, and we assume a parameterization which we denote by $\mathbf{y}=\Gamma_{j}(\tau)$. Here, $\mathbf{y}$ is a point on the obstacle corresponding to a unique value $\tau\in[0,1]$ in the parameter domain, which we assume to be the interval $[0,1]$ for all obstacles. We will use $\mathbf{y}$ and $\tau$ interchangeably to denote a point. For clarity, we will often use an indexed parameter variable $\tau_{j}$ to denote a point on $\Gamma_{j}$.

We assume that the periodic orbit with $J$ obstacles is known: i.e., the obstacles, their order and the corresponding reflection points $\Gamma_{j}(\tau_{j}^{*})$ on them. Given an ordered collection of obstacles $\{\Omega_{j}\}_{j=1}^{J}$, the points ${\mathcal{T}}=\{\tau_{j}^{*}\}$ in the corresponding periodic orbit minimise the total length

[TABLE]

In the latter sum as well as further on, we define $\Gamma_{J+1}=\Gamma_{1}$ in order to shorten the formulae.

We assume that there is a unique solution to this minimisation problem, which is clear if all the obstacles involved are convex [16]. However, we do present an example later on with numerical results where this condition is not met. It appears, for example, that the formulae continue to hold if all objects are locally convex and a local minimizer of the form above can be found.

In the expressions above, we have assumed for simplicity of presentation that the periodic orbit contains exactly the obstacles $\Omega_{j}$, $j=1,\ldots,J$, in that order. Of course, one may consider a periodic orbit between any subset of obstacles, say, $\Omega_{3}$ and $\Omega_{5}$. The formulae that follow remain valid simply by considering the multiple scattering problem involving only $\Omega_{3}$ and $\Omega_{5}$, i.e., with $J=2$.

The mode that corresponds to the orbit ${\mathcal{T}}$ is the collection of densities $\{V_{j}(\tau)\}_{j=1}^{J}$ on each of the obstacles in ${\mathcal{T}}$. The distance function between points on two consecutive obstacles is briefly written as

[TABLE]

Finally, we call the function $g(x)$ in an asymptotic expression of the form

[TABLE]

the phase of $f$. The notation above indicates that $f$ and $u$ depend on a wavenumber $k$. We note right away that the phase is not uniquely defined. In the decomposition we assume that $u$ is non-oscillatory, in the sense that its derivatives do not grow with increasing $k$. Even then, any constant $c$ added to $g$ can readily be absorbed into $u$. Hence, $g$ can only be defined up to a constant shift.

2.2 Coupling operators and the eigenvalue problem

In order to formulate the eigenvalue problem, recall the single layer potential 1. We shall use the notation

[TABLE]

for its application to a density $V_{i}$ defined on $\Gamma_{i}$, and with the point $\mathbf{x}$ on $\Gamma_{j}$. An incident wave on $\Gamma_{j}$ is written as $u_{j}^{\textrm{inc}}$. Disregarding the possibility of resonances, for the time being, the single layer potential leads to coupled integral equations of the first kind for a multiple scattering sound-soft Dirichlet problem. For a set of two obstacles, the coupled system is given by

[TABLE]

Here, $m_{1}$ and $m_{2}$ are the densities on obstacles $\Gamma_{1}$ and $\Gamma_{2}$ and the operators $S_{ji}$ are the single layer potentials 5. One way to solve this problem, alluded to as ray tracing in the introduction, is to first ignore the other obstacle by leaving out the cross-terms and then to trace consecutive reflections. We numerically illustrate that this converges at the start of § 7.1 and this technique is applied analytically in [20] as well.

The operator under investigation in this paper for $J=2$ is given succinctly by

[TABLE]

which represents a full cycle of reflections over two obstacles, independently of the incident wave. With this formulation, the goal of the paper is to find an eigenfunction of $T$ of the form 2, such that

[TABLE]

We proceed by analyzing the operator $T$ asymptotically for large $k$. Expression 7 consists of a number of steps. It is sufficient to analyze only $S_{11}$ and $S_{21}$, since $S_{22}$ and $S_{12}$ are entirely analogous in the case of just two obstacles. We proceed with the latter operator $S_{21}$ first.

2.3 Field scattered by the first obstacle onto the second

When we calculate the field that $\Gamma_{1}$ scatters onto $\Gamma_{2}$, we obtain

[TABLE]

where we have substituted the Green’s function $K(x,y)$ of 1 by its explicit expression for the Dirichlet 2D Helmholtz problem. In order to analyze the integral asymptotically for large $k$, we want to explicitly characterize its total phase.

Since the points $\Gamma_{1}(\tau_{1})$ and $\Gamma_{2}(\tau)$ are well separated for all $\tau\in[0,1]$, and since $k$ is assumed to be large, we can replace the Hankel function by the leading order term of its expansion for large arguments [27, (10.17.5)],

[TABLE]

Furthermore, we assume $V_{1}$ to have the form given by 2. Finally, we parameterize the boundary $\Gamma_{1}$ using the parameterization $\Gamma_{1}(\tau)$ with the same name and with Jacobian $\|\nabla\Gamma_{1}(\tau)\|$. This results in the asymptotic approximation

[TABLE]

where $\Delta_{1,2}$ is the Euclidean distance between two points on the two obstacles, recall 4,

[TABLE]

We are mainly interested in the leading order behaviour of the integral. Higher order terms can be found by substituting the full expansion of the Hankel function, a process detailed in [20]. Yet, this entails also a full asymptotic expansion of $v_{1}^{\text{smooth}}(\tau_{1})$, which as of yet is unknown. Importantly, the total phase of the integrand 10 is fully explicit,

[TABLE]

As is standard in asymptotic expansions of integrals [31, 23, 6], one notes that the leading order contribution to an oscillatory integral originates in a region around the stationary points of the phase of the integrand. A stationary point $\chi$ of a phase function $g(x)$ is such that $g^{\prime}(\chi)=0$, i.e., it is a root of the derivative of $g$. In the context of scattering, such stationary points correspond precisely to points of reflection [20].

In order to further approximate 10 asymptotically, we look for a stationary point $\chi_{1}(\tau_{2})$ for each value of $\tau_{2}$. Physically, this means that a ray hitting obstacle $\Gamma_{2}$ at the point $\tau_{2}$ originated in the point $\tau_{1}=\chi_{1}(\tau_{2})$ on $\Gamma_{1}$: it is like inverse ray tracing. The condition for $\chi_{1}(\tau_{2})$ is that the derivative of the total phase vanishes with respect to the integration variable $\tau$, hence

[TABLE]

This is an implicit definition of the function $\chi_{1}$. The phase of the resulting function $S_{12}V_{1}$ is given by

[TABLE]

It equals the phase $\phi_{1}$ at the reflection point, plus the distance from that point to $\Gamma_{2}(\tau_{2})$. This agrees with the principle that the phase along a ray equals the distance travelled.

Recall that $\tau_{j}^{*}$ are the points of the periodic orbit. We know by construction that $\chi_{1}(\tau_{2}^{*})=\tau_{1}^{*}$: this expresses that the periodic orbit contains the ray from $\tau_{1}^{*}$ to $\tau_{2}^{*}$. Thus, though $\chi_{1}$ is characterized only implicitly above, we know at least its value in one point. Furthermore, due to the conditions on the obstacles, the function is necessarily at least locally single-valued.

Finally, we can also define the function $\chi_{2}(\tau_{1})$ to find the point on $\Gamma_{2}$ from which a ray reflects onto $\Gamma_{1}(\tau_{1})$. Note that $\chi_{2}$ is not the inverse of $\chi_{1}$, but we do have in this scenario of just two obstacles in the orbit that both

[TABLE]

2.4 Single scattering problem on the first obstacle

Next, we consider $S_{11}^{-1}$. A full computation of the asymptotic expansion of the solution to a single scattering problem, even with phase extraction, is rather involved. The computation is carried out in [20]. The problem of computing the phase of the solution is simpler since, in fact, it is known that the phase of the solution equals the phase of the incoming wave from $\Gamma_{2}$ under the conditions in [16]. As such,

[TABLE]

For completeness, the process of obtaining full asymptotics is as follows. Assuming a density of the form 2, the application of the integral operator on $\Gamma_{1}$ yields

[TABLE]

A technical complication is the appearance of the logarithmically singular Hankel function, which in this case can not be replaced by its expansion for large arguments as before in 10, where it led to a simpler oscillatory exponential. An alternative is to substitute the Hankel function by its Mehler-Sonine integral representation, which does have a complex exponential in the integrand [27, (10.9.10)]

[TABLE]

After splitting the integration interval around $\tau_{1}^{*}$, inverting the distance function and interchanging integration variables, a regular steepest descent analysis can be applied on the resulting bivariate integral.

The main observation to make here is that the only relevant asymptotic contribution to the singular oscillatory integral originates in the singularity at $\tau_{1}=\tau$. As one can apply a path of steepest descent such that the Mehler-Sonine integral representation of $H_{0}^{(1)}$ does not contribute to the phase of the integrand (because the singularity is precisely where the argument of $H_{0}^{(1)}$ is zero), the phase of $S_{11}V_{1}$ indeed equals the phase of $V_{1}$.

2.5 System of nonlinear equations for two obstacles

We have sufficient information to describe the action of the operator $T$ given by 7 on functions of the form 2. We look for a function $V_{1}$ with phase $\phi_{1}(\tau)$, such that the density after a full cycle of reflections has the phase $\phi_{1}(\tau)+\mu$, where $\mu$ represents a constant phase shift. Using 12–15, the nonlinear system of differential-algebraic equations that $\phi_{1}(\tau)$, $\phi_{2}(\tau)$, $\chi_{1}(\tau)$ and $\chi_{2}(\tau)$ should satisfy is:

[TABLE]

Analogously to (11), we have introduced the phase of the integral representation of $S_{12}V_{2}$,

[TABLE]

The first three equations in (18) were introduced before. The final equation introduces the constant phase shift $\mu$, and ensures that $\phi_{1}$ is recovered after a full orbit up to $\mu$. Since the phase equals the distance travelled, this implies that

[TABLE]

where $L({\mathcal{T}})$ is the length of the periodic orbit 3.

We have thus derived at a system of $4$ equations, which is highly nonlinear in the four unknown functions $\phi_{1}(\tau)$, $\phi_{2}(\tau)$, $\chi_{1}(\tau)$ and $\chi_{2}(\tau)$. Even in the case of two circular obstacles with equal radius, the system does not seem to allow an explicit analytic solution.

It is not obvious whether the system has a solution at all and, if so, whether the solution is unique.

Remark 2

We make no rigorous claims regarding the (unique) solvability of 18 globally for $\tau_{1},\tau_{2}\in[0,1]$. Indeed, since our extraction of phases is based on geometrical optics, it is not valid in the shadow region of the obstacles (the shadow relative to the previous obstacle in the orbit). By the no-occlusion condition of [16] and formulated in Remark 1, the periodic orbit is such that rays close to it are amenable to geometrical optics. We assert that a solution to 18 exists in a neighbourhood of the points $\tau_{j}^{*}$ on each $\Gamma_{j}$, $j=1,\ldots,J$.

In order to find a local solution to 18, we augment the system with the initial conditions 14. Furthermore, it is clear that the phase $\phi_{1}$ is only determined up to a constant, which we may choose, e.g., by prescribing a value for $\phi_{1}(\tau_{1}^{*})$.

In the following section, we will resort to Taylor series expansions of the unknown functions to approximate the local solution of 18. This approach does lead to fully explicit expressions for the coefficients. The series is computed around the periodic orbit, where the mode $V_{j}(\tau)$ has a (nearby) local maximum. Though our derivation is not presented in a fully mathematically rigorous way in this paper, the fact that the computation succeeds is a strong indication that the system 18, along with the initial conditions for $\chi$ and the function value $\phi_{1}(\tau^{*})$, is indeed uniquely solvable in a neighourhood of the periodic orbit.

2.6 System for a general number of obstacles

We can readily generalize the arguments of this section to an arbitrary number of obstacles $J$. The orbit operator $T$ becomes

[TABLE]

We can consider each pair of obstacles involved in the reflection of rays, the obstacles $(j+1,j)$, for $j=1,\ldots,J-1$, and the final pair $(J,1)$. An increasing $j$ means the rays progress forward in time through the obstacles and hence diverge away from the periodic orbit. The function $\chi_{j}(\tau_{j+1})$ instead contracts around the periodic orbit, i.e., it traces the rays backwards in time from obstacle $j+1$ back to $j$.

We define the phase function of the integral representation of $S_{ji}V_{i}$ (recall (16)) as

[TABLE]

where the distance $\Delta_{i,j}$ between points on obstacles $\Gamma_{i}$ and $\Gamma_{j}$ is given by 4.

With that notation, we can formulate the following $2J\times 2J$ nonlinear system:

[TABLE]

Here, as before, the phase shift after a full cycle equals the length of the periodic orbit,

[TABLE]

This means that the phase of the ray that follows the periodic orbit increases exactly by the length of the path that is travelled.

In order to characterize a local solution to 21, we augment the system with the initial conditions

[TABLE]

and we choose a value for $\phi_{1}(\tau_{1}^{*})$.

3 The special case of two disks

In this section, we derive explicit Taylor expansions for the phases $\phi_{1}$ and $\phi_{2}$ of the mode $\{V_{1},V_{2}\}$ of a periodic orbit between two disks.

3.1 Parameterization of the boundary

The non-overlapping disks in Figure 3 are separated by a distance $d$ and the corresponding shortest path is the periodic orbit along which rays will be trapped forever. We choose to parameterize the circular boundaries with radius $r$ specifically as, with $\tau\in[0,1]$,

[TABLE]

Since the length of the periodic orbit is twice the distance, we have $\mu=L({\mathcal{T}})=2d$. The points in the periodic orbit are $\tau_{1}^{*}=\tau_{2}^{*}=0$. When we show numerical values further on, we choose the parameters

[TABLE]

The unknown phase function for the density on $\Gamma_{1}$ was denoted by $\phi_{1}(\tau)$ before. Due to the geometrical symmetry, and because the parameterization is defined in opposite directions on the two obstacles, we have that $\phi_{2}(\tau)=\phi_{1}(\tau)+\nu$, where $\nu$ is a constant phase factor. It follows from 18 at $\tau^{*}=0$ that $\nu=d$, hence

[TABLE]

In the following we omit the subscript, i.e. $\phi_{1}(\tau)=\phi(\tau)$. Similarly, we write $\chi_{1}(\tau)=\chi_{2}(\tau)=\chi(\tau)$.

3.2 Procedure

We look for a Taylor series around $\tau^{*}=0$ of the form

[TABLE]

We are free to choose the constant $c_{0}$, since it corresponds purely to a constant phase shift and eigenfunctions are only determined up to a constant factor. We arbitrarily set $c_{0}=d$.

We also expand $\chi(\tau)$ into a Taylor series,

[TABLE]

as well as the residuals of the first two equations in (18):

[TABLE]

Note that the third and fourth equation of (18) lead to the exact same expansion, recalling that $\phi_{2}(\tau)=\phi_{1}(\tau)+d$ and $\mu=2d$.

We already know that $a_{0}=0$ since $\chi(\tau_{1}^{*})=\tau_{2}^{*}=0$. Substituting $\tau=0$ into (27) also shows that $\omega_{0}=0$. All other coefficients are computed in a recursive procedure in the next subsection. For each index $i$, starting from $i=0$, we solve

[TABLE]

for the unknowns $a_{i}$ and $c_{i+1}$. These conditions express that the residuals of the nonlinear equations vanish to all orders. The newly found coefficients for index $i$ can be used in the next iteration for index $i+1$.

3.3 Explicit Taylor series coefficients

In the first iteration of (28) for $i=0$, starting from $\omega_{1}=0=\psi_{0}$ we find that

[TABLE]

This corresponds to $\phi^{\prime}(0)=0$, which indicates that $\phi$ has a local extremum at the point $\tau^{*}=0$ as we expect.

For $i=1$, we need the series expansion of the left hand side of 27 up to order two, and of the left hand side of 26 up to order $1$. Substituting the expansions for $\chi$ and $\phi$ leads to

[TABLE]

where the latter is evaluated at $\tau_{\alpha}=a_{0}+a_{1}\tau+a_{2}\tau^{2}$. We also expand the distance function $\Delta_{1,2}$ into a Taylor series. This requires the expansion of the sines, cosines and the square root function in the expressions above, which can all be performed explicitly. Inserting the newly found value $c_{1}=0$, the equations (28) become the system of quadratic equations

[TABLE]

from which we deduce

[TABLE]

Because (33) is quadratic in $a_{1}$, there is a second solution, which is given by $a_{1}=3+2\sqrt{2},c_{2}=-\sqrt{2}\pi^{2}$. This solution is not physically relevant, as the left side of Figure 6 indicates that $c_{2}>0$ and, furthermore, $a_{1}>1$ would imply that the rays move away from the periodic orbit shown in the left part of Figure 3. Instead, we want $\chi$ to be contracting.

For higher $i$, the systems of equations that arise are linear in all coefficients and the solution is unambiguous. Due to the symmetry in the current problem of two disks, we find that $a_{i}=c_{i+1}=0$ for even values of $i$. Further computation leads to the next nonzero coefficients

[TABLE]

and

[TABLE]

Recall that these values are specific to the choice $d=1$ and $r=1/2$, though they can easily be generalized by repeating the derivation.

The accuracy of these formulae is depicted graphically in the right panel of Figure 3. We have computed a numerical approximation to the eigenfunction of the operator $T$ (described further on in § 6), and extracted the Taylor series approximation of the phase. The mode is non-oscillatory in a neighbourhood of the critical point $\tau=0$, it is supported in the region where the disks can ‘see’ each other, and the oscillatoriness decreases with increasing number of terms in the expansion.

4 The case of multiple obstacles with general convex shapes

We derive explicit Taylor expansions for the phases of the mode $\{V_{j}\}_{j=1}^{J}$ of a periodic orbit with a general number $J$ of obstacles with a general shape. The process is largely comparable to the special case of two disks, yet slightly more involved. One issue that complicates the algebra is the expansion of the distance function $\Delta_{j,j+1}(\tau_{j},\tau_{j+1})$, which is a function of two arguments. It could be done explicitly for the case of two disks, but in general it can only be expanded in terms of the Taylor series expansions of the parameterization of the boundaries at the points $\tau_{j}^{*}$.

4.1 Computation of the periodic orbit

In the analytical derivation that follows, we assume that the periodic orbit ${\mathcal{T}}=\{\tau_{j}^{*}\}_{j=1}^{J}$ is known, given the ordered collection of $J$ obstacles. In practice, of course, one has to compute the reflection points $\Gamma_{j}(\tau_{j}^{*})$ on each of the obstacles $\Gamma_{j}$. Recall that the periodic orbit minimizes the sum of the distances between each pair of consecutive obstacles. To that end, one has to solve the minimization problem 3.

In our implementation, we find the corresponding parameters $\tau_{j}^{*}$ on each obstacle $\Gamma_{j}$ by a nonlinear optimisation, more specifically the quasi-Newton algorithm in fminunc in Matlab based on supplied gradient and Hessian, and with all tolerances set to machine precision. The initialisation is via the parameters $\tau_{j}$ giving the point approximately closest to the mean of $100$ approximately equidistant points on each obstacle. Though the solution of the minimisation problem is unique for a collection of convex obstacles, the numerical procedure typically also converges to a periodic orbit for non-convex obstacles. We present an example further on in § 7.3.

4.2 Series expansion of the bivariate distance function

The distance function $\Delta_{j,j+1}(\tau_{j},\tau_{j+1})$ appears in the system 21 and has to be expanded in both of its arguments. For this purpose, we determine a double series expansion of the distance between two points $\tau_{j}$ and $\tau_{j+1}$ on the obstacles $\Gamma_{j}$ and $\Gamma_{j+1}$, where we define $\Gamma_{J+1}=\Gamma_{1}$ for convenience of notation.

Since we work in two dimensions, the parameterization $\Gamma_{j}$ is vector-valued. We define its vector-valued Taylor series in terms of the Taylor series of the $x$ and $y$-components,

[TABLE]

The coefficients $\Gamma_{j,x,n}$ and $\Gamma_{j,y,n}$ can be obtained directly from the parameterization $\Gamma_{j}(\tau_{j})$ by differentiation in the point $\tau_{j}^{*}$.111Note that the index $n$ starts at $1$ here, chosen such that there is an exact match between the expressions in this paper and the formulae in our implementation in Matlab [29] (in which vector indices start at $1$).

Substituting 39 into expression 4 for the distance function leads, after some algebraic manipulation further detailed below, to

[TABLE]

Here, we have introduced several intermediate variables. The values $\Lambda_{j,x,i,l}$ and $\Lambda_{j,y,i,l}$ arise from the full expansion of the square of the difference of two expansions of the form 39 for $j$ and $j+1$. For the $x$-component, the relation is

[TABLE]

from which we derive the expression for $\Lambda_{j,x,m,l}$ after the calculation of two convolutions:

[TABLE]

The expression for $\Lambda_{j,y,i,l}$ is entirely analogous. We have used the Kronecker delta ($\delta_{n,k}=1$ if $n=k$ and zero elsewhere) to describe the constant terms when $m=1$ or $l=1$.

The coefficients $z_{j,m,i,k}$ arise from the expansion of the square root, and they are given by

[TABLE]

The coefficient $z_{j,1,1,1}$ is not used. As an implementation note, when the two last indices in $z_{j,m,i,k}$ are smaller than $a$ the range of the subscript $m$ relevant for the computation of the expansion is limited to $1\leq m\leq 2a-2$.

The final expression for the distance function has the form

[TABLE]

where $f_{j,1,1}$ is the distance between two consecutive points in the orbit,

[TABLE]

This is the value around which the square root has been expanded above.

Finally, the coefficients that we are looking for in this section are

[TABLE]

Since $\tau_{j}^{*}$ is part of the periodic orbit which minimizes the total distance, we have that

[TABLE]

This leads to the property

[TABLE]

4.3 Procedure

More general than in § 3.2, the points $\tau_{j}^{*}$ may be nonzero. We look for a Taylor series of the form

[TABLE]

As before, we are free to choose the constant $c_{0}$ and we set $c_{j,0}=d_{j}$, the distance to the next obstacle on the periodic orbit. For notational convenience we again define $\Gamma_{J+1}=\Gamma_{1}$, $\phi_{J+1}=\phi_{1}$ and so on.

We expand $\chi$ and the residuals of the equations in 21 into Taylor series as well:

[TABLE]

We already know that $a_{j,1,0}=\tau_{j}^{*}$ since $\chi_{j}(\tau_{j+1}^{*})=\tau_{j}^{*}$. The coefficients $\omega_{j,i}$ and $\psi_{j,i}$ can be computed explicitly in terms of $a_{j,1,i}$ and $c_{j,i}$, using the bivariate Taylor series expansion of the distance function derived above in § 4.2. All other coefficients are computed in a recursive procedure explained in the next subsection. As before, we want the coefficients in the expansion of the residuals to be zero, with the exception of $\omega_{j,0}$ which may be nonzero because of 15.

4.4 Formulae for general 2D obstacles

We explicitly expand the power of an expansion using a convolution, leading to

[TABLE]

with

[TABLE]

and with $a_{j,i,i}=(a_{j,1,1})^{i}$.

Then, 49 and 50 respectively become

[TABLE]

This gives a set of $2J$ equations in $2J$ unknowns $a_{j,1,i}$ and $c_{j,i}$. For $i=1$, this becomes

[TABLE]

Note that, when considering the analogous equations for the obstacles $j+1$ and $j+2$, the first equation above becomes $c_{j+1,1}=-f_{j+1,2,1}$. Combined with the second equation above, this is consistent with 46.

From the Taylor series 47 it is clear that the value $c_{j,1}$ is the derivative of the phase at the point $\tau_{j}^{*}$ on the periodic orbit. If it is zero, the density has an extremum at the critical point, which is related to that the periodic orbit is exactly the shortest distance between the obstacles. Since the signs of $c_{j,1}$ and of $f_{j,2,1}$ are opposite, the phase $\phi_{j}(\tau_{j})$ decreases when moving towards the minimum of the distance $\Delta_{j,j+1}(\tau_{j},\tau_{j+1}^{*})$. This minimum does not correspond to the minimum of the phase $\phi_{j}^{\prime}(\tau_{j}^{M})=0$. Physically, the latter should be the point where a ray leaves orthogonal to $\Gamma_{j}$ and reaches the periodic orbit after an infinite number of reflections, where $\chi_{j}(\tau_{j+1}^{*})=\tau_{j}^{*}\neq\tau_{J}^{M}$ if $c_{j,1}\neq 0$, see Figure 4. Note that the phase depends on the direction of the periodic orbit as well now, in contrast to a periodic orbit with only two obstacles.

The coefficient $a_{j,1,1}$ is still unknown. So in the next steps over increasing $i$, all $\omega_{j,i}$ are combined with $\psi_{j,i-1}$ to obtain $c_{j,i}$ and $a_{j,1,i-1}$ simultaneously for all $j$. For $i=2$, 53 and 51 become

[TABLE]

The latter simplification was possible due to the equation for $\psi_{j,1}$. This gives a system of $2J$ quadratic equations which can be solved with a nonlinear solver with an initial guess $c_{j+1,2}=f_{j,1,3}$ and $a_{j,1,1}=0$. A motivation for the latter is that $\chi_{j}(\tau_{j+1})$ should be a function that contracts around the periodic orbit, so $a_{j,1,1}$ should be small: more specifically, $|a_{j,1,1}|<|f_{j,1,2}/f_{j,2,1}|$. We add this check to our code as well as $c_{j,2}>0$, because we expect $\tau_{j}^{*}$ to lie close to a (local) minimum of $\phi_{j}$ in analogy to the case $J=2$. The initial guess $c_{j+1,2}=f_{j,1,3}$ makes sure we already satisfy $\omega_{j,2}$, one of both equations. In the case of the two disks, $f_{1,2,2}=-\pi^{2},f_{1,3,1}=3\pi^{2}/2=f_{1,1,3}$ such that the above equations correspond to 33 and 34. We expect that if $f_{j,i,m}=0$ for all $i+m\leq n+1$ and $j\in[1,J]$, a polynomial system of $2J$ equations, each of which has a total degree $\leq n$, will give the solution $a_{j,1,1}$ and $c_{j,n}$ for $1\leq j\leq J$. Physically, an increasing $n$ means that both obstacles become flatter at the periodic orbit, resulting in less oscillatory modes.

For $i\geq 3$ in 51 and 53, the following $2J\times 2J$ linear system was implemented

[TABLE]

The validity of our approach is confirmed by the numerical experiments further on involving three obstacles.

5 A geometrical interpretation of $\chi$ and $\phi$

We can assign a geometrical meaning to the stationary point function $\chi(\tau)$ and the phase $\phi(\tau)$ in the case of two disks, by expanding the analogy to ray tracing. Moreover, the reasoning extends to general periodic orbits.

Recall that $\chi(\tau_{1})$ is the point on $\Gamma_{2}$ from which a ray hits the point $\Gamma_{1}(\tau_{1})$. By symmetry, the converse also holds: $\chi(\chi(\tau_{1}))$ yields a point back on $\Gamma_{1}$. In the sequence $\chi(\chi(\tau_{1}))\to\chi(\tau_{1})\to\tau_{1}$, it should be true that the incoming and reflected ray at the point $\chi(\tau_{1})$ on $\Gamma_{2}$ are at equal angles with the normal direction. In our setting of two disks or radius $r$ and at a distance $d$ of each other, this leads to the equation

[TABLE]

This equation is again highly nonlinear, and it is hard to obtain a symbolic solution on account of the application of $\chi$ on itself. However, one can compute the Taylor series of this expression for $\tau_{1}$ near zero, and equate all its coefficients to zero, and we checked this to result in the same $a_{i}$ as in the procedure described before.

One can also numerically approximate $\chi$ using 62. In turn, we can obtain the phase in a limited range of $\tau$ from (11)–(12). Indeed, differentiating $\eqref{eq:g21}$ with respect to $\tau_{1}$, and substituting into (12), leads to

[TABLE]

Multiplying both sides by $\chi_{1}^{\prime}(\tau_{2})$ and integrating with respect to $\tau_{2}$ yields the explicit expression

[TABLE]

This determines $\phi_{1}$, up to inversion of $\chi(\tau)$ in the argument of $\phi_{1}$ above. The accuracy of the approximation given by 63 is included in the right panel of Figure 6, which shows uniform accuracy in a range of $\tau$. This could be extended to a collection of $J$ general obstacles, but that would require repeated integration and differentiation.

We may characterize the phase further in geometrical terms. As in ray tracing, one can track the phase of a ray in a homogeneous medium along straight lines. The phase of a ray is equal to its initial phase plus the travelled distance. Hence, the phase is the sum of the lengths of all the reflecting rays. We define the $r$-times iterated application of $\chi$ by $\chi^{[r]}(\tau)$, such that $\chi^{[0]}(\tau)=\tau$, $\chi^{[1]}(\tau)=\chi(\tau)$, $\chi^{[2]}(\tau)=\chi(\chi(\tau))$ and so on. We say that a point $\Gamma_{1}(\tau)$ is in the range of the critical point $\tau^{*}$ if the repeated application of $\chi$ leads to the critical point,

[TABLE]

This amounts to inverse ray tracing.

Summing the lengths of the reflected rays leads to, for $\tau$ in the range of $\tau^{*}$,

[TABLE]

This formula holds for the case of two disks. Each term in the sum corresponds to a cycle from $\Gamma_{1}$ to $\Gamma_{2}$ and back. Note that we subtract the length of the periodic orbit $L({\mathcal{T}})=2d$ in each cycle, otherwise the sum would diverge. The formula agrees with the repeated application of the fourth and second equation in 18, each time recursively replacing $\phi_{1}$ or $\phi_{2}$ in the right hand sides by their expression in terms of the other function.

For a collection of $J$ obstacles and a periodic orbit ${\mathcal{T}}$, the reasoning can be extended by tracing a point $\tau_{1}$ on $\Gamma_{1}$ back to the critical point $\tau_{1}^{*}$ on the periodic orbit. One sums the lengths of each cycle of reflections, following $\chi_{J}$ down to $\chi_{1}$, and subtracts in each cycle the total length $L({\mathcal{T}})$ of the orbit. As above, this corresponds to the recursive application of the second and fourth equation of 21.

The phase at $\tau_{1}$ thus computed is the distance of the shortest path traveled by the ray to the periodic orbit, in agreement with the principle of Fermat. At each obstacle, there is reflection at equal angles, and in each cycle the distance between the obstacles along the periodic orbit is subtracted.

6 Numerical approximation

We briefly sketch a numerical method to validate the expressions obtained in this paper, based on an implementation of the standard boundary element method for the coupled system 6. Each operator $S_{ji}$ results in a matrix $A_{j,i}$, and we can realize the operator $T$ given by 7 (for two obstacles) or 19 (for multiple obstacles) numerically.

We assume for each boundary $\Gamma_{j}$ a set of functions $\psi_{j,n}(\tau)$ to discretize the densities with $N_{j}$ unknowns, see 67 further on. The discretization of a coupled problem with $J$ obstacles results in a dense BEM matrix with block structure form

[TABLE]

Here, $A_{j,i}\in\mathbb{C}^{N_{j}\times N_{i}}$ represents the discretization of $S_{ji}$. Thus, the blocks $A_{j,j}$ on the diagonal are the BEM matrices of the single scattering problems involving only $\Gamma_{j}$. The off-diagonal blocks $A_{j,i}$ represent the coupling from $\Gamma_{i}$ to $\Gamma_{j}$. Our experiments are based on a collocation approach, in which we have used piecewise linear basis functions and $N_{j}=10k_{j}$ evenly spaced collocation points on each boundary $\Gamma_{j}$ [24].

In the ray-tracing scheme for multiple scattering problems, one only uses the inverse of the diagonal blocks, as these correspond to a single scattering problem [16]. It is shown in [16] that, under certain conditions, the norms of the off-diagonal blocks are sufficiently small to allow the iterative process to converge. In addition to the convexity and non-occlusion condition already mentioned, the obstacles should be separated by more than the inverse of the wavenumber. More specifically, regularity estimates are available if $\Delta_{j,i}(\tau_{j},\tau_{i})\geq 1/k$ for all $1\leq i,j\leq J$ and $\tau_{i},\tau_{j}\in[0,1]$, see [18, §4.3] and [9].

The numerical analogue of $T$, a full cycle of reflections in a periodic orbit given by 19, is obtained simply by replacing any operator $S_{ji}$ with its discretization $A_{j,i}$,

[TABLE]

We can numerically compute the eigenvalue decomposition $M=EDE^{-1}$, ordered such that its largest (in modulus) eigenvalue is $D_{1,1}$. Experiments confirm that the complex argument of $D_{1,1}$ equals $L({\mathcal{T}})$, the length of the periodic orbit, up to a multiple of $2\pi/k$. We proceed with only the first eigenvalue and eigenvector in what follows, noting only that the further eigenvalues seemingly decay rapidly, while the phases of the corresponding eigenvectors are similar to that of the first eigenvector. The modulus of the eigenvalues seems to only depend on the geometry and not on (high) $k$ nor the number of points per wavelength, where increasing the latter two only adds negligible eigenvalues.

We obtain an approximation to the phase of the eigenmode from the phase of 67, where the coefficients in the expansion are given by the first eigenvector. Recall the factorization of 2 of $V_{j}$ into an oscillatory and non-oscillatory part, which is approximated numerically using the first eigenvector in $E$,

[TABLE]

To obtain the first eigenvectors on other obstacles than the first, one can reorder 66 and recompute the eigenvalue decomposition $M=EDE^{-1}$, but in our implementation we simply multiply $E_{n,1}$ on the first obstacle with the corresponding pairs of submatrices in 66. We can not merely consider the complex phase of point evaluations of $V_{j}$ to approximate $\phi_{j}$, since $v_{j}^{\text{smooth}}$ may also be complex-valued. However, assuming that the complex argument of $v_{j}^{\text{smooth}}$ varies much more slowly than Arg$(\exp[ik\phi_{j}(\tau)])$ for moderately large $k$, we can estimate the phase $\phi_{j}$ at a point $\tau$ from a finite difference with small positive or negative $\delta$,

[TABLE]

Note that the final term simply removes any jump of size $2\pi$ in the argument that may have occurred from one point to the next. In our implementation, $|\delta|$ is the distance between collocation points in the $\tau$-domain and because we use piecewise linear basis functions on each boundary, $\tilde{V}_{j}(\tau-\delta)$ simply equals $E_{n,1}$ at the collocation point $\tau_{1,n}=\tau-\delta$. The numerical phase function was already illustrated in the right panel of Figure 3 and used to validate our Taylor series expressions for $\phi_{1}(\tau)=\phi_{2}(\tau)$ in the case of two disks. Several other examples are given in the next section.

Though in this paper we present no analysis of the amplitude function $v_{j}^{\text{smooth}}$ of the eigenmode, it seems concentrated in a region around the critical point $\tau_{j}^{*}$. As one would expect, the support is effectively limited to that part of the obstacles that ‘see’ each other. The amplitude decays by a factor equal to the modulus of the largest eigenvalue after each cycle of reflections. This is also illustrated in the next section.

7 Numerical results

7.1 Two disks

We have already illustrated the accuracy of the Taylor series approximation to the limiting phase function $\phi$ in the right panel of Figure 3. There, the eigenmode was computed numerically using the numerical method outlined in the previous section, i.e., as the eigenvector of $M$ given by 66. The Taylor series of the phase is computed using the explicit expressions for the coefficients given in § 3.3. The figure confirms, at least visually, that the remainder after factoring out the Taylor series of the phase is non-oscillatory, in a neighbourhood of the critical point $\tau^{*}$.

We elaborate on the case of two disks with further experiments that are more quantitative.

First, we illustrate in Figure 5 the convergence of the phase to a limiting function in consecutive iterations of the multiple scattering problem. We consider two circular obstacles again, which is similar to for example [17, Fig. 1 & 2] mentioned in § 1. The incident wave is a plane wave from the left (left panel) or a point source incidence (right panel). The density that is induced by the incident wave on the second disk is represented as a blue dotted line. In case of the plane wave, the density has a maximum and is least oscillatory near $\tau_{2}=-1/4$, which corresponds to the leftmost point on the disk. Already after the first reflection, the peak of the density is quite close to $\tau^{*}=0$. After two reflections, the peak is visually indistinguishable from the critical point $\tau^{*}$. Furthermore, the amplitude decreases by a certain factor after each reflection due to rays leaving the scene. In the right panel of the same figure, the results are repeated for a point source incident wave, and the conclusions are similar. These experiments were performed using submatrices from the numerical method outlined in the previous section.

Next, we turn to the computation of the limiting phase function in a range of $\tau$, rather than as a Taylor series. As outlined in § 5, we can compute $\chi$ from 62. We have applied six Newton-type iterations on 62 with an initial guess given by $\chi(\tau)\approx(3-2\sqrt{2})\tau$, using the Chebfun package [13]. The resulting function satisfied 62 with a residual close to machine precision. Furthermore, its Taylor series coefficients closely match our coefficients $a_{i}$. We have added a few iterations of $\chi$ as dashed purple lines in the left panel of Figure 3. Recall that applying $\chi$ corresponds to tracing a ray back in time, and as such it contracts to the periodic orbit, as the figure illustrates.

The right panel of Figure 6 quantifies the convergence rate for a diverse set of approximations to the phase. Here, we have chosen a wavenumber $k=2^{9}$. We compute the phase using the geometric sum of distances given by 64. We end the summation after $r=10$ reflections, as the next term is close to machine precision in absolute value. The Taylor series with series expansion coefficients $c_{i}$ computed in § 3.3 converges to this phase at a rate $\tau^{T}$ when taking $T$ terms. Note that the accuracy of $\tilde{\phi}$, recall 68, is limited to $1e-6$ in this figure. This is due to the discretization error of the boundary element method, the composition of BEM submatrices, the computation of an eigenvector, approximating its argument and trying to undo the modulo $2\pi$ in it. The phase can indeed also be computed as 63, although in a smaller domain due to the inverse of $\chi$.

In the left panel of Figure 6 we illustrate that the limiting phase seems well approximated by two physical distances $\zeta$ and $\xi$ that are easily computed and are illustrated in the left panel of Figure 3. The value $\zeta$ is the distance from a point on $\Gamma_{1}$ to the point $\Gamma_{2}(0)$, i.e., the distance to the critical point on the other obstacle. The value $\xi$ is the distance from a point on $\Gamma_{1}$ to its closest point on $\Gamma_{2}$. Neither of these values is an exact expression of the phase, since they do not take into account reflections. However, we merely want to point out that both values are reasonable approximations, and they may prove to be useful as starting values for the solution of the non-linear system of equations 12–15.

7.2 Two general obstacles

The scattering configuration shown in Figure 7 does not exhibit the symmetry of the two disks. Moreover, the first obstacle is not convex. Still, we can consider the eigenmode corresponding to a periodic orbit connecting the two closest points on the two obstacles. In the right panel, we see that the mode is again locally nonoscillatory at $\tau_{1}^{*}\approx 0.62$ and its modulus still has a local maximum there. Our expressions for the Taylor series coefficients remain applicable. Factoring out a cubic or quintic approximation of the phase around $\tau_{1}^{*}$ does result in non-oscillatory behaviour of the eigenfunction near the critical point, as shown in the right panel of Figure 7. However, the quintic approximation does not seem to be much better than the cubic one.

The accuracy is quantified in Figure 8. The numerical results for the two obstacles are shown in the left panel. Higher degree Taylor series approximations of the phase do lead to higher accuracy in a larger neighbourhood of the critical point. The accuracy in this figure is limited to $4e-8$ because of similar reasons as with the two disks. For relative errors well above this limit, we also see that the slope increases when adding more terms, as expected. So although less clear than the right part of Figure 6, Figure 8 suggests we obtain the correct Taylor approximation of the phase.

We see a decline in Figure 7 of the absolute value of the mode at the right of $\tau_{1}^{*}\approx 0.62$, which corresponds to the convex part at the bottom of $\Gamma_{1}$. More reflections are possible in the nonconvex part at the left of $\Gamma_{1}$, resulting in a slower decrease of $|V_{1}(\tau_{1})|$ for $\tau_{1}$ decreasing from $0.62$. Consequently, one can now distinguish two different curves for each number of terms in Figure 8, where the one with the highest error corresponds to $\tau_{1}<0.62$.

7.3 Three general obstacles

Finally, we consider a configuration involving three obstacles ($J=3$), shown in Figure 9. Here, the critical point $\tau_{3}^{*}$ on the periodic orbit is in a convex part. However, the third obstacle is strongly non-convex, and has a near-inclusion with extremely complicated high-frequency behaviour in the range $0.33<\tau_{3}<0.67$, $\tau_{3}^{*}$ is just outside this range and the periodic orbit is quite flat. Moreover, there is an additional $J=3$ periodic orbit near $\tau_{3}=0.67$, which may influence the one under investigation.

Still, even though the phase of the mode on the ellipse is influenced by the near-inclusion obstacle, the right part of Figure 9 shows we obtain a good approximation of it. The error in the phase is quantified in the right panel of Figure 8, and it is clearly worse than in the previous example. Yet, somewhat surprisingly, even at the relatively small wavenumber $k=128$, high accuracy is reached near the critical point. The plots as a function of $|\tau_{3}-\tau_{3}^{*}|$ now again have two clearly visible distinguished parts: the top one corresponds to $\tau_{3}>\tau_{3}^{*}$ which is the inclusion part, while the bottom one corresponds to $\tau_{3}<\tau_{3}^{*}$, the convex part. The accuracy is clearly better in the convex part of the domain.

8 Concluding remarks and future research

The asymptotic expressions in this paper are independent of the wavenumber and the incident wave. As a result, one can approximate the limiting phases of the densities on each obstacle in a periodic orbit at a cost that is independent of the wavenumber. That is a significant advantage. However, in any set of $J$ obstacles, there is potentially a very large number of periodic orbits. Indeed, there is a periodic orbit for any subset of the $J$ obstacles. This combinatorial cost of multiple scattering problems is inherent to the ray-tracing methodology: the cost of ray-tracing is independent of the wavenumber, but increases rapidly with an increasing number of obstacles.

Still, the results of this paper may be used to accelerate an iterative hybrid numerical-asymptotic method (such as the methods of Ecevit et al in [16, 2, 14, 15, 7]) as follows. A large cost of the ray-tracing is the long tail, i.e., a large number of iterations has to be performed before the difference between two consecutive iterations is small. Previous efforts have focused on accelerating the computation of the tail. However, the phases settle down quickly after only a few iterations. By extracting the asymptotic phase simultaneously for all densities on all obstacles at once, it may be possible to solve the tail at once as a coupled non-oscillatory problem. That is, we propose to perform a few iterations using phase-extraction in every step, followed by the solution of a coupled problem with the asymptotic phase extracted from each density. Since the consecutive iterations lead to the same phase in every cycle, the solution of the coupled problem is expected to be non-oscillatory. This is a topic of possible future research.

The cost of phase-extraction in the first few iterations can be further reduced if it is known how the phase of the incident wave evolves to the limiting phase we computed as a function of the number of reflections. A series expansion of $v_{j}^{\text{smooth}}(\tau)$ for the limiting density would even allow this without needing a mesh for $V_{j}(\tau)$ at all, from which also the magnitude of the associated eigenvalues can be derived.

Acknowledgments

The authors would like to thank Samuel Groth and Niels Søndergaard for interesting and useful discussions on topics related to this paper. The authors were supported by FWO Flanders [projects G.0617.10, G.0641.11 and G.A004.14].

Список литературы

• [1]

T. Abboud, J.-C. Nédélec, and B. Zhou.

Méthode des équations intégrales pour les hautes fréquences.

CR Acad. Sci. I Math., 318:165–170, 1994.

• [2]

A. Anand, Y. Boubendir, F. Ecevit, and F. Reitich.

Analysis of multiple scattering iterations for high-frequency scattering problems. II: The three-dimensional scalar case.

Numer. Math., 114:373–427, 2010.

• [3]

A. Asheim and D. Huybrechs.

Local solutions to high-frequency 2D scattering problems.

J. Comput. Phys., 229:5357–5372, 2010.

• [4]

A. Asheim and D. Huybrechs.

Extraction of uniformly accurate phase functions across smooth shadow boundaries in high frequency scattering problems.

SIAM J. Appl. Math., 74(2):454–476, 2014.

• [5]

V. M. Babich and V. S. Buldyrev.

Short-wavelength diffraction theory.

Springer-Verlag, Berlin, 1991.

• [6]

N. Bleistein and R. A. Handelsman.

Asymptotic expansions of integrals.

Dover Publications, Inc., 31 East 2nd Street, Mineola, N.Y., 1986.

• [7]

Y. Boubendir, F. Ecevit, and F. Reitich.

Acceleration of an iterative method for the evaluation of high-frequency multiple scattering effects.

SIAM Journal on Scientific Computing, 39(6):B1130–B1155, 2017.

• [8]

O. Bruno, C. Geuzaine, J. J. Monro, and F. Reitich.

Prescribed error tolerances within fixed computational times for scattering problems of arbitrarily high frequency: the convex case.

Phil. Trans. R. Soc. Lond. A, 362:629–645, 2004.

• [9]

S. N. Chandler-Wilde, A. Gibbs, A. Moiola, and S. Langdon.

Hybrid numerical asymptotic approximation for multiple scattering problems.

In Proceedings of the 12th International Conference on Mathematical and Numerical Aspects of Wave Propagation, pages 130–131, Karlsruhe, Germany, July 2015.

Book of abstracts available at: http://waves2015.math.kit.edu/boa/book-of-abstracts.pdf.

• [10]

S. N. Chandler-Wilde, I. Graham, S. Langdon, and E. Spence.

Numerical-asymptotic boundary integral methods in high-frequency acoustic scattering.

Acta Numer., pages 89–305, 2012.

• [11]

S. N. Chandler-Wilde, D. P. Hewett, S. Langdon, and A. Twigger.

A high frequency boundary element method for scattering by a class of nonconvex obstacles.

Numerische Mathematik, 129(4):647–689, 2015.

• [12]

D. L. Colton and R. Kress.

Integral equation methods in scattering theory.

John Wiley & Sons, New York (N.Y.), 1983.

• [13]

T. A. Driscoll, N. Hale, and L. N. Trefethen.

Chebfun Guide.

Pafnuty Publications, 2014.

http://www.chebfun.org.

• [14]

F. Ecevit and H. c. Özen.

Frequency-adapted galerkin boundary element methods for convex scattering problems.

Numer. Math., 135:27–71, 2017.

• [15]

F. Ecevit and F. Reitich.

A high-frequency integral equation method for electromagnetic and acoustic scattering simulations: rate of convergence of multiple scattering iterations.

In Proceedings of the Seventh International Conference on Mathematical and Numerical Aspects of Wave Propagation, pages 145–147, 2005.

• [16]

F. Ecevit and F. Reitich.

Analysis of multiple scattering iterations for high-frequency scattering problems. I: The two-dimensional case.

Numer. Math., 114:271–354, 2009.

• [17]

C. Geuzaine, O. Bruno, and F. Reitich.

On the O(1) solution of multiple-scattering problems.

IEEE Trans. Magn., 41(5):1488–1491, 2005.

• [18]

A. Gibbs.

Numerical methods for high frequency scattering by multiple obstacles.

PhD thesis, University of Reading, 2017.

www.reading.ac.uk/web/files/maths/THESIS_Gibbs.pdf.

• [19]

S. P. Groth, D. P. Hewett, and S. Langdon.

Hybrid numerical-asymptotic approximation for high-frequency scattering by penetrable convex polygons.

IMA J. Appl. Math., 80:324–353, 2015.

• [20]

S. P. Groth, D. Huybrechs, and P. Opsomer.

High-order terms in the ray expansion for high-frequency scattering by single and multiple obstacles.

In preparation.

• [21]

J. Hart, T. Anthonsen, and E. Ott.

Effect of short ray trajectories on the scattering statistics of wave chaotic systems.

Phys. Rev., 80:041109–1/15, 2009.

• [22]

D. Huybrechs and P. Opsomer.

Coupling modes in high-frequency multiple scattering problems: the case of two circles.

In D. J. Chappell, editor, Proceedings of the Eleventh UK Conference on Boundary Integral Methods (UKBIM 11), pages 99–108. Nottingham Trent University: Publications, 2017.

• [23]

D. Huybrechs and S. Vandewalle.

On the evaluation of highy oscillatory integrals by analytic continuation.

SIAM J. Numer. Anal., 44(3):1026–1048, 2006.

• [24]

D. Huybrechs and S. Vandewalle.

An efficient implementation of boundary element methods for computationally expensive Green’s functions.

Eng. Anal. Boundary Elem., 32(8):621–632, August 2008.

• [25]

T. Kottos and U. Smilansky.

Periodic orbit theory and spectral statistics for quantum graphs.

Ann. Phys., 274:76–124, 1999.

• [26]

C. Linton and D. Evans.

The interaction of waves with arrays of vertical circular cylinders.

J. Fluid Mech., 215:549–569, 1990.

• [27]

NIST.

NIST Digital Library of Mathematical Functions.

http://dlmf.nist.gov/, 2016.

Online companion to [28].

• [28]

F. W. J. Olver, D. W. Lozier, R. F. Boisvert, and C. W. Clark, editors.

NIST Handbook of Mathematical Functions.

Cambridge University Press, New York, NY, 2010.

Print companion to [27].

• [29]

P. Opsomer.

Eigenmodes of periodic orbits for 2d multiple scattering problems, 2018.

asyBEM release 2 https://github.com/popsomer/asyBEM.

• [30]

G. Tanner and N. Søndergaard.

Wave chaos in acoustics and elasticity.

J. Phys. A: Math. Theor., 40:R443–R509, 2007.

• [31]

R. S. Wong.

Asymptotic approximations of integrals.

SIAM, 3600 Market Street, 6th floor, Philadelphia (Pa.) 19104-2688 USA, 2001.

Republication of 1944.

• [32]

T. Wu.

Boundary element acoustics.

WIT Press, 2000.

Reprint 2005.