Transfer operator approach to ray-tracing in circular domains

The computation of wave-energy distributions in the mid-to-high frequency regime can be reduced to ray-tracing calculations. Solving the ray-tracing problem in terms of an operator equation for the energy density leads to an inhomogeneous equation which involves a Perron-Frobenius operator defined on a suitable Sobolev space. Even for fairly simple geometries, let alone realistic scenarios such as typical boundary value problems in room acoustics or for mechanical vibrations, numerical approximations are necessary. Here we study the convergence of approximation schemes by rigorous methods. For circular billiards we prove that convergence of finite-rank approximations using a Fourier basis follows a power law where the power depends on the smoothness of the source distribution driving the system. The relevance of our studies for more general geometries is illustrated by numerical examples.


Introduction
Ray-tracing methods serve as an important toolkit in finding approximate solutions of linear wave equations in the high frequency limit. This approximation is used in a variety of fields providing, for example, the connection between Maxwell's equations and geometric optics, as well as between quantum mechanics and classical Hamiltonian mechanics [13]. The ray-tracing limit has also been considered in detail in acoustics, seismology and mechanical vibrations [23]. In engineering applications, ray tracing is employed in handling electromagnetic problems, such as coverage estimates for 5G or WiFi communication [12], room acoustics simulations [22] as well as structure-borne sound propagation in mechanical structures [6]. Finding closed form, analytical solutions to such engineering problems of sufficient complexity is generally impossible, even using ray-tracing techniques, and one has to use numerical methods instead.
For solving linear wave problems such as those listed above, the numerical methods used have to be adapted to the relevant length and frequency scales involved. In the low frequency regime, finite element methods (FEM) are routinely employed for resolving the full wave dynamics. However, the number of degrees of freedom in an FEM model needs to scale with the wavelength and there is thus an upper limit in frequency above which the required computational resources become unfeasible. At very high frequencies, power balance approaches can often be used as long as certain assumptions on the ergodicity of the underlying ray dynamics are satisfied [24]. In the mid-to-high frequency range, raytracing becomes the method of choice; standard ray-tracing techniques track all possible rays from a source to a receiver point [22] -a method which becomes cumbersome if many reflections need to be taken into account. As an alternative Dynamical Energy Analysis (DEA) was proposed and has proven to be useful in particular for structureborne sound problems [24,14]. Instead of tracking individual rays carrying vibrations across the complex structure -which is extremely challenging -in DEA, the problem is reformulated in terms of densities of rays, which are then mapped across a mesh representing the structure [7,8]. This reduces the ray tracing problem from tracking rays on complicated and curved domains to mapping ray segments across small, plane patches of a simple shape forming the mesh, typically triangular or quadrilateral mesh cells. The ray densities are then mapped from one cell of the mesh to adjacent ones and the overall transport problem can be formulated in terms of an inhomogeneous equation of the form where f 0 is the initial ray density, L a Perron-Frobenius type operator describing the evolution of ray densities and f the required final ray density. Using DEA, the distribution of vibrational energy in mechanical structures, such as ships, cars and tractors [15,14] can be calculated successfully. For such realistic geometries, equation (1) above cannot be solved analytically, so recourse is made to numerical schemes based on heuristic finite-dimensional matrix approximations of the operator L. To date, very little is known about the convergence properties of these schemes and the dependence of the convergence rate on the ray dynamics, as well as the discretisation techniques [8]. The precise form of convergence is likely to be highly sensitive to both the basis functions used in approximating the inhomogeneous equation (1), as well as dynamical and damping properties of the system under investigation [15]. For our study, we will therefore be concerned with the approximation of L by operators of finite rank. There is a plethora of papers on numerical approximation of Perron-Frobenius operators, starting with Ulam's method of phase space discretisation, finite section or Galerkin methods, and data-driven methods, see for example [2,10,11,17,19] to mention but a few. Surprisingly, the application of DEA (which falls into the Galerkin category) to even fairly simple geometries has not been dealt with at a rigorous level. Here, we shall thus focus on one of the simplest cases, the billiard dynamics given by the ballistic motion within a circular disk. We shall establish rigorous error bounds of finite-dimensional approximations for the resulting energy distribution.
In order to set up the required notation, consider a particle moving inside a circular billiard table D being specularly reflected at its boundary ∂D. We parametrise ∂D by the polar angle x ∈ R/2πZ and we denote by y ∈ [−π/2, π/2] the angle of reflection that the postcollisional velocity vector has with the inward normal to ∂D. Initially the collision angle is defined on an interval. It is, however, technically simpler to deal with cyclic variables. Since both angles −π/2 and π/2 correspond to a particle which sticks on the boundary we identify both angles so that the collision angle becomes a cyclic variable as well. With these conventions, the collision map T on the domain Ω = (R/2πZ) × (R/πZ) can be written as with its inverse φ = T −1 given by It is not difficult to see that the collision map T preserves the normalised Haar measure on Ω. The long-term statistical behaviour of T can thus be studied by investigating the associated Perron-Frobenius operator (see, for example, [5]), which for invertible measure-preserving maps is given by the composition operator C φ defined as where f : Ω → C. In the current work we are interested in the properties of a weighted Perron-Frobenius operator, also known as a transfer operator. In order to define it, let us first introduce a multiplication operator M w acting on functions f : Ω → C by where w : Ω → [0, ∞) is a suitable weight function, which in the DEA framework accounts for dissipation caused either by collisions with the wall or by in-flight dissipation. The transfer operator, understood to be acting on a suitable space of functions detailed in the following section, is now given by In the present article, we are interested in approximations of the solution to the operator equation (1) with f 0 : Ω → [0, ∞) interpreted as the initial boundary density of particles induced by the first boundary collision of particles emitted by a source located in the interior of D (see [24]). In the DEA approach this quantity represents the energy source. The resulting energy distribution is captured by the solution, f : Ω → [0, ∞), which gives the stationary boundary density generated by the collision dynamics. Given a suitable Banach space and a sequence of finite-rank projections (P K ) K∈N , an approximation method for (1) can be constructed by considering the projected finite-dimensional problem The aim of this work is to present a Banach space for f 0 and (P K ) K∈N , so that problem (7) has solutions, which converge in a suitable topology to the solution of (1) as K tends to infinity, with the speed of convergence being of the order K −α . The exponent α depends on the smoothness of f 0 and the requirements imposed on the type of convergence.
In passing we note that transfer operators have their roots in statistical mechanics [20,21] and nowadays play an important role in the ergodic theory of smooth expanding, or more generally, hyperbolic dynamical systems (see, for example, [3,4]). The main reason for their popularity in this context derives from the fact that for expanding or hyperbolic dynamical systems the transfer operator, when considered on a suitable function space, can be shown to have discrete peripheral spectrum, from which long-term statistical properties of the underlying system can be derived. In the elliptic setting, however, such as for the circular billiard considered in this article, analogous results cannot be expected, and, as a consequence, transfer operator methods have received little attention in this context. It is perhaps worth noting that in our setting we do not require discreteness of the peripheral spectrum of the transfer operator. The main onus is to show that the resolvent of the transfer operator exists at the point 1 (see equation (1)) and can be effectively approximated by finite-rank operators (see equation (7)).
As we intend to keep our presentation accessible to non-specialists, we will occasionally elaborate on aspects covered in the specialised literature but which may not be well known to a general audience. The remaining parts are organised as follows. In Section 2 we introduce Sobolev spaces, on which the transfer operator and its finitedimensional approximations are bounded operators with spectral radii bounded away from 1. In Section 3 we shall prove the convergence results for the operator equations (1) and (7) stated as Theorem 3.4. In the final Section 4 we summarise the main findings, compare the formal results with numerical simulations and explore the relevance of the current study in a wider context.

Sobolev spaces and transfer operators
We will be interested in certain subspaces of L 2 (Ω) = L 2 (Ω, m) where dm = dxdy/(2π 2 ) is the normalised two-dimensional Lebesgue measure on Ω. The natural inner product is given by An orthonormal basis of L 2 (Ω) is given by One can rewrite this definition in terms of Fourier coefficients. Using the fact that (8) can be expressed as Using equation (9) we can define fractional Sobolev spaces H s (Ω) for s = (s 1 , s 2 ) ∈ R 2 + as which are Hilbert spaces when equipped with the inner product given in equation (9) with m replaced by s. We shall next investigate the properties of the composition operator C φ associated with the map φ in (3) on the fractional Sobolev space H s (Ω).
Proof. For any n ∈ N and (x, y) ∈ Ω we have φ n (x, y) = (x − nπ + 2ny, y), and thus for any k ∈ Z 2 . ‡ This choice of inner product is sometimes referred to as the modified inner product, in contrast with the classical one (see, for example, [18, Def 2.2]).
In order to show that the operator is bounded we will need the following general inequality. Let (x, y) ∈ [0, ∞) 2 and t ≥ 0, then Using equation (11) we obtain the bound Since (C n φ e k , C n φ e l ) H s = 0 for k = l, the operator norm of C n φ is bounded from above by (1 + (2n) 2s 2 max(1, 2 2s 2 −1 )) 1/2 , resulting in the upper bound for the spectral radius In order to see that the inequality above is an equality, observe that the operator norm of C n φ is bounded from below by 1 as C n φ e 0 H s = e 0 H 2 . Thus r(C φ ) = 1.
Before proceeding we note that by (10), the action of the composition operator on H s (Ω) can be represented by the action of the matrix A = 1 0 1 1 on Fourier coefficients. In particular, we have . Then for any n ∈ N 0 we can define a finite-rank operator Lemma 2.4. Let C φ and P Λ K be as above. Then Proof. This follows by checking the equality for any basis element e k and noting that A n is invertible.
Definition 2.5. Let M w denote the multiplication operator as defined in equation (5), considered as an operator on H s (Ω), with a smooth weight function w : Ω → [0, ∞). In addition, we assume that w has the following properties: (a) w ∞ = sup x∈Ω |w(x)| < 1; (b) w is bounded away from zero; (c) w(x, y) = w(x ′ , y) for any (x, y), (x ′ , y) ∈ Ω, that is, the weight w does not depend on the first argument.
Remark 2.6. The operator M w models the effect of damping on the motion of the billiard particle. Assumptions (a) and (b) imply that the damping is well-behaved, while assumption (c) is innocuous, given the circular symmetry of the billiard table.
The following two lemmas summarise basic properties of M w and C φ .
Lemma 2.7. Let M w , C φ and P Λ K be as above. Then we have the following.
Proof. Items (i) and (iii) follow from Definition 2.5(c); items (ii) and (iv) follow by direct computation using the map We write L K = P Λ K M w C φ P Λ K for the finite-rank approximation of L = M w C φ . Using Lemma 2.7 (i) and Lemma 2.4, we can write L n K for n ∈ N as In order to state the properties of L and L K we need to introduce the following multi-index notation: an n-dimensional multi-index is an n-tuple i n = (i 1 , i 2 , . . . , i n ) of non-negative integers of order |i n | = i 1 +i 2 +· · ·+i n = m; the corresponding multinomial coefficient is given by Lemma 2.8. Let M w , C φ and P Λ l K be as above. Then we have the following.
Proof. Item (i) follows by induction over m using Lemma 2.7(iv) for the base case m = 1. For item (ii), the additional induction over n follows by rewriting (i) We are now ready to prove the main result of this section. Keeping in mind that we assume that the billiard dynamics is dissipative, that is, the weight is chosen so that w ∞ < 1, the following lemma shows that, given f 0 ∈ H s (Ω), the problem (1) and the projected version (7) have unique solutions f ∈ H s (Ω) and f K ∈ H s (Ω), respectively.
Proof. We shall only prove statement (i), as the proof of statement (ii) follows by almost identical arguments. In the following, we shall assume that s 1 ≥ s 2 ≥ 1, as the case s 1 s 2 = 0 follows by identical arguments. For f ∈ H s (Ω) we have Let p, q ∈ N with p ≤ s 1 and q ≤ s 2 . It is not difficult to see that for any f ∈ H s (Ω) and K ∈ N 0 the following holds.
Here, statements (c) and (d) follow by writing the L 2 norm of D p x f and D q y f using Parseval's identity.
Writing L n K as in equation (14) and using (a) and (b) above iteratively we have where the last inequality follows from the fact that the operator norm of C φ on L 2 (Ω) equals 1.
As D x commutes with any of the operators involved (Lemma 2.7 (ii,iii,vi)) we have in the second term on the right-hand side of (15) that D s 1 x L n K = L n K D s 1 x . By the same argument as above we have In order to bound the last term in equation (15) we are using Lemma 2.8(iii) and Hölder's inequality in order to write We shall first obtain a bound for D j y C n φ f L 2 . Using Lemma 2.8(ii) and decomposing the sum in terms of powers of D x and D y we obtain where we have used the multinomial formula |in|=k k in = n k . Thus, for j ≤ m we obtain using Hölder's inequality, the multinomial formula and upper bounds for 2 j−i n+1 where the last inequality uses (c) and (d).
Next we shall obtain a bound on the operator norm of A j for j ≤ s 2 . First note that M D l y w = M w M (D l y w)/w is well-defined as w is bounded away from zero. By using (a) and (b) iteratively, for any i n = (i 1 , . . . , i n ) with |i n | = s 2 we have where C s 2 = max i 1 ,...,in |in|=s 2 n l=1 D i l y w/w 2 ∞ ≤ max 0≤l≤s 2 D l y w/w 2s 2 ∞ is a constant independent of n. Using arguments analogous to those used to obtain inequality (18), we obtain the bound Using the estimates (16), (17), (18) and (19) in equation (15) we arrive at the bound withC n,s 2 ≤ (s 2 + 1)s 2 (n + 1) 4s 2 2 2s 2 C s 2 + 1. AsC n,s 2 is independent of K, the family (L K ) K∈N is a uniformly bounded family of bounded operators on H s (Ω). Finally, taking the right hand side of equation (15) to the power of 1/n and observing thatC n,s 2 grows polynomially in n, the upper bound for the spectral radius of L K follows.

Convergence properties
In the previous section we established (see Lemma 2.9) that given f 0 ∈ H s (Ω), the problem (1) and the projected version (7) have unique solutions f ∈ H s (Ω) and f K ∈ H s (Ω), respectively. We shall now turn to establishing the convergence of f K to f . This would be straightforward if we knew that L K → L as K → ∞ in the operator norm on H s (Ω), since then, using the so-called second resolvent identity we would have from which convergence of f K → f in H s (Ω) could be readily obtained. This, however, cannot be the case, as if L K → L as K → ∞ in the operator norm on H s (Ω), then L, as a uniform limit of finite-rank operators, would be compact on H s (Ω). However, as L has a bounded inverse on H s (Ω), it cannot be compact.
We thus need to resort to a slightly weaker notion of convergence, that is, we shall consider the transfer operator as an operator between Sobolev spaces of different order. In passing, we remark that this idea is also at the heart of one of the most successful techniques to obtain spectral approximation results of transfer operators, where perturbation sizes are measured in 'triple' norms (see, for example, [16]).
In the following we shall explain this idea in more detail. We start with the following important observation. For t, s ∈ [0, ∞) 2 with s 1 ≥ s 2 > t 1 ≥ t 2 ≥ 0, functions in H s (Ω) can be identified with functions in H t (Ω) using the embedding operator J : H s (Ω) ֒→ H t (Ω) given by J f = f . This operator is not just continuous, but also compact, as the following lemma shows. Proof. Let f ∈ H s (Ω). Using the notation a t (k) = 1 + |k 1 | 2t 1 + |2k 2 | 2t 2 we have For this, first observe that which follows by Hölder's inequality and s 1 ≥ s 2 . Then, Now, by bounding from above each (1 + |k 1 | 2 + |2k 2 | 2 ) −α with its maximal value in each of the sums, we obtain We are now able to show that L can be approximated by finite-rank operators when considered as operators from H s to H t . Proposition 3.2. Let L K = P K LP K be the finite-rank approximation of L on H s (Ω) with s ∈ N 2 and s 1 ≥ s 2 . Let J be as above and t ∈ N 2 0 with s 2 > t 1 ≥ t 2 . Then for some C > 0 and α = s 2 − t 1 .
Proof. Let L ′ denote the transfer operator when considered on the larger space H t (Ω). Then using the property J L = L ′ J , we have Thus, where we have used Lemma 2.9, Lemma 3.1 and P K H s →H s ≤ 1. Proposition 3.3. Let L and the family (L K ) K be as above, considered as operators on H s (Ω) where s ∈ N 2 with s 1 ≥ s 2 . Then, for t ∈ N 2 0 with s 2 > t 1 ≥ t 2 and for all K ∈ N we have Proof. As r(L) ≤ w ∞ < 1 by Lemma 2.9, the operator (I−L) −1 exists and is bounded. Let (L ′ K ) K denote the family of transfer operators when considered on the larger space H t (Ω). Similarly, as ρ(L ′ K ) ≤ w ∞ < 1 and the norms of (L ′ K ) n are bounded uniformly in K by Lemma 2.9, the sums ∞ n=0 L K n H t →H t are bounded by a constant independent of K and therefore (I − L K ) −1 H t →H t is uniformly bounded in K. Using the property J (I − L K ) = (I − L ′ K )J and the second resolvent identity (see equation (20)) we have Using Proposition 3.2 for the bound on J (L − L K ) H s →H t finishes the proof.
We are finally able to state and prove our main convergence result.
Theorem 3.4. Let L and the family (L K ) K be as above, considered as operators on H s (Ω) with s 1 , s 2 ∈ N and s 1 ≥ s 2 > t 1 ≥ t 2 ≥ 0. Then for f 0 ∈ H s (Ω) the operator equations (1) and (7) have unique solutions f ∈ H s (Ω) and f K ∈ H s (Ω), respectively. Moreover there exist a constant C > 0 such that for all K ∈ N we have Proof. The statement follows by writing f = (I − L) −1 f 0 , f K = (I − L K ) −1 f 0 and using Proposition 3.3.

Remark 3.5.
Note that for f 0 ∈ Λ K = P K (H s (Ω)), the unique solution f K to (7) also lies in the finite-dimensional space Λ K , so that (7) can be solved as a truly finitedimensional problem.

Discussion and numerical experiments
Let us first summarise and rephrase our results in intuitive terms. Since the linear operator in equation (1) fails to be compact, any finite-dimensional matrix representation would not reflect properties of the operator at all. Nevertheless the finite-dimensional representation in (7) provides a meaningful approximation for the solution of the inhomogeneous equation. For smooth periodic functions in location and angle of reflection, the solution of the approximated problem (7) converges to the solution of (1) in the Sobolev norm. The approximation error depends on the degree of smoothness of the inhomogeneous part. In addition, the approximation error is measured in a weaker norm, for instance the frequently used L 2 norm for the choice t = (0, 0). The properties of this weaker norm also determine the speed of convergence. Broadly speaking, the convergence rate obeys a power law with the exponent being determined by the smoothness of the energy source and the norm used to measure the approximation error. A finite amount of dissipation is a crucial ingredient in the entire approach, that is, the weight w has to satisfy w ∞ < 1. The simplest choice of a constant weight, w(x, y) = µ < 1, corresponds to a dissipation which occurs at each collision at the boundary, for example, an attenuation of the sound wave caused by an inelastic reflection at the boundary of the cavity. Proper modelling of the damping parameters involved is a crucial aspect of the method and is necessary to describe realistic problems accurately [14]. For example, a linear attenuation in the medium would result in a path-length dependent weight w(x, y) = exp(−2µ cos(y)). This choice, however, does not obey the stipulated bound as orbits with angles close to y = ±π/2 have arbitrarily small path length, and hence small dissipation between subsequent collisions. We could overcome this particular problem by restricting the angle of reflection to nontangential collisions, that is, y ∈ (−(1 − ǫ)π/2, (1 − ǫ)π/2) for a small ǫ > 0, effectively constraining the permitted type of energy source. This however requires changing the Hilbert space and the projection operators, as the validity of c k (D m y f ) = (i2k) m c k (f ) and D y P k = P k D y is no longer given for a smooth function f on an interval instead of on a circle. One suitable choice could be the space of functions in H s (Ω ǫ ) with vanishing weak derivatives D ν f on the boundary. A suitable basis is then the basis of Daubechies wavelets [26].
To illustrate the impact of Theorem 3.4, we perform numerical simulations of circular billiards with constant damping w(x, y) = µ. As a proxy for the error estimate we use the distance between approximations of subsequent order f K+1 − f K H t , which obeys essentially the same upper bound Strictly speaking we have established this bound for integer vales of t ℓ and s ℓ only. With a little more effort this could be remedied by appealing to interpolation theory [25]. For simplicity of exposition we shall not pursue this here. For our numerical considerations we take the liberty to apply the bound above for non-integer values. For the norm · H t , which estimates the truncation error, we use the choices t = (0, 0), that is, the L 2 norm, and t = (1, 1), a norm which is just outside the set of exponents guaranteeing pointwise convergence.
The transfer operator's action on Fourier modes is given in equation (12). In order to use it for a numerical test, we have to use a representation for all Fourier modes, see equation (A.2). We show results for three different choices of the initial boundary density f 0 . They have in common that their support is given by supp(f 0 ) = {(x, y) : x ∈ [π/6, π/6 + 4π/3] , y ∈ [−0.8, 1.2]} .
• Case W 2 : a smooth function given by f 0 (x, y) = x(1 −x) ỹ(1 −ỹ) The data shown in Figure 1 confirm the upper bound in Theorem 3.4. For the L 2 norm, t 1 = t 2 = 0, we observe, in each case, convergence at a rate which is slightly faster than the theoretical prediction α = s 2 −t 1 . The power law decay of the truncation error shows up for large values of K and the onset of this scaling region shifts towards larger values if the initial boundary density becomes smooth. This should not come as a surprise, since the resolution of higher order derivatives requires higher order Fourier modes. For the parameter at the boundary of point-wise convergence t = (1, 1), we see that the discontinuous boundary density fails to converge in line with our theoretical predictions. While Theorem 3.4 does not guarantee convergence in case W 1 either, the numerical data suggest an extremely slow convergence which is still consistent with the upper bound estimate α = s 2 −t 1 = 1 −1 = 0. Finally, for the smooth boundary density (case W 2 ) we observe a convergence rate slightly faster than the theoretical prediction.
From a dynamical perspective, circular billiards are trivial since the billiard map (2) is an integrable twist map. In order to get an idea of how dynamical properties impact on convergence properties we show numerical results for a deformed circle billiard which displays mixed regular and chaotic dynamics. For the deformation we choose the radius to depend on the polar angle x according to where we choose m = 3 in the following. Deformations of this kind are known in the literature as Limaçon billiards [1]. We will cover the cases δ = 0.01 and δ = 0.1. For larger values, the billiard fails to be convex. In order to demonstrate the change in dynamical behaviour, Figure 2 shows the Poincare plot of the collision map T . For a small value of the deformation, δ = 0.01, one still observes a fairly large number of invariant tori in accordance with general KAM folklore. The larger perturbation shown in Figure 2, δ = 0.1, destroys most of the regular motion and renders the system chaotic with a few exceptions, for example, the highlighted period-3 island. In order to calculate the convergence of the energy distribution we have to evaluate the matrix elements of the transfer operator. For the circular billiard, the only non-zero entries take the value ±µ and follow the structure given by equation (A.2). Once the circle has been deformed, the analytic calculation of the matrix elements is no longer possible. Even worse, the collision map is not given in closed analytic form either, so that an efficient numerical calculation becomes a nontrivial task (see the appendix for details). However, we are able to reduce the calculation of the matrix elements to double integrals with the kernel being given in closed analytic form, see equation (A.3). Nevertheless, the numerical evaluation is still time consuming, in particular, since the matrix is no longer sparse. Hence, we can only calculate finite approximations up to K = 30. In order to reach the scaling regime (see Figure 1 for comparison) we employ a stronger damping of µ = 0.1. The results for the error measured in L 2 norm, that is, for the choice t 1 = t 2 = 0, are shown in Figure 3.
It is quite remarkable that the decay of the error is apparently almost unaffected by  the degree of chaoticity. Hence the rigorous error estimate of Theorem 3.4 which covers the case δ = 0 seems to have a wider range of applicability. While intuitively such an observation would not be surprising for nearly integrable cases it is quite counterintuitive that the same error estimate may hold as well in strongly chaotic situations. However, our proof does not cover any of the deformed billiards and there does not seem to be an obvious way how the methodology can be generalised to these complicated cases. Nevertheless, it is reaffirming that our study of a simple dynamical system like the circular billiard has relevance for more complex dynamical behaviour.
matrix elements M l,k of the transfer operator read M l,k = 1 2π 2 2π 0 π/2 −π/2 (C φ e k )(x, y) e l (x, y) dydx = 1 2π 2 2π 0 π/2 −π/2 e ik 1 φx(x,y)−il 1 x+2ik 2 φy(x,y)−2il 2 y dydx with k = (k 1 , k 2 ) and l = (l 1 , l 2 ). In case of the perfect circle we get a representation which is given by a sparse matrix with only a few non-zero elements, close to the main diagonal, namely (C φ e l )(x, y) = This is the extension of equation (12) to all Fourier modes and it was used to calculate the values for Figure 1.
In order to eliminate the implicitly defined collision map we change integration variables from (x, y) to (x, x ′ ). Using y 1 (x, x ′ ) = y and y 2 (x, x ′ ) = y ′ for the two scattering angles the matrix elements become M l,k = 1 2π 2 2π 0 2π 0 ∂y 2 (x, x ′ ) ∂x e i(k 1 x−l 1 x ′ ) e 2i(k 2 y 1 (x,x ′ )−l 2 y 2 (x,x ′ )) dx ′ dx, (A. 3) where the additional factor is the Jacobian of the coordinate transformation. In contrast to the collision map T , the expressions y 1 (x, x ′ ) and y 2 (x, x ′ ) can be obtained in closed analytic form so that equation (A.3) is easier to implement numerically. Figure A1 shows a sketch of two subsequent collisions. The first scattering angle y 1 is given in terms of an inner product sin(y 1 ) = d · t/(|d||t|) .
Since the position vector of the initial point is given by r(x)e r the tangent is easily obtained as t = r ′ (x)e r + r(x)e ϕ . The vector separating the two points of collision is given in terms of the local basis vectors by d = (r(x ′ ) cos(x ′ − x) − r(x))e r + r(x ′ ) sin(x ′ − x)e ϕ .
Hence the closed form expression for the first scattering angle reads t y e r e ϕ d Figure A1. Geometric configuration of two subsequent collisions in a convex billiard with a particle moving from point 1 (with parameter value x) to point 2 (with parameter value x ′ ). We also depict the ray vector d, the tangent vector t, and the unit vectors e r and e ϕ in polar coordinates.